# Implicit kinetic schemes for the Saint-Venant system

Chourouk El Hassanieh, Mathieu Rigal, Jacques Sainte-Marie

► To cite this version:

Chourouk El Hassanieh, Mathieu Rigal, Jacques Sainte-Marie. Implicit kinetic schemes for the Saint-Venant system. 2024. hal-04048832v2

HAL Id: hal-04048832

https://hal.science/hal-04048832v2

Preprint submitted on 28 Mar 2024

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

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

<figure><img src="image_1.png" /></figure>
# Implicit kinetic schemes for the Saint-Venant system

Chourouk El Hassanieh, Mathieu Rigal, Jacques Sainte-Marie

## ► To cite this version:

Chourouk El Hassanieh, Mathieu Rigal, Jacques Sainte-Marie. Implicit kinetic schemes for the Saint-Venant system. ESAIM: Mathematical Modelling and Numerical Analysis, 2025, 59, pp.1863-1908. 10.1051/m2an/2025038. hal-04048832v3

HAL Id: hal-04048832

https://hal.science/hal-04048832v3

Submitted on 22 Sep 2025

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

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

<figure><img src="image_1.png" /></figure>
Implicit kinetic schemes for the Saint-Venant system

Chourouk El Hassanieh¹, Mathieu Rigal², and Jacques Sainte-Marie¹

¹Centre Inria de Sorbonne université, 2 rue Simone Iff, CS 42112, 75589 Paris Cedex 12 and Sorbonne Université, Univ. Paris Diderot, SPC, CNRS, LJLL, F-75005 Paris

²Institut de Mathématiques de Bordeaux, équipe EDP et physique-mathématique, Université de Bordeaux, 351 cours de la Libération, F-33405 Talence, France

Emails: (chourouk.el-hassanieh, jacques.sainte-marie)@inria.fr,
mathieu.rigal@math.u-bordeaux.fr

March 28, 2024

**Abstract**

Explicit (in time) kinetic schemes applied to the nonlinear shallow water equations have been extensively studied in the past. The novelty of this paper is to investigate an implicit version of such methods in order to improve their stability properties. In the case of a flat bathymetry we obtain a fully implicit kinetic solver satisfying a discrete entropy inequality and keeping the water height non negative without any restriction on the time step. Remarkably, a simplified version of this nonlinear implicit scheme allows to express the update explicitly which we implement in practice. The case of varying bottoms is then dealt with through an iterative solver combined with the hydrostatic reconstruction technique. We show that this scheme preserves the water height non-negativity under a CFL condition and satisfies a discrete entropy inequality without error term, which is an improvement over its explicit version. An extension of the implicit and iterative methods to the two dimensional case is also discussed. Finally we perform some numerical validations underlining the advantages and the computational cost of our strategy.

This work has been supported by a DIM MathInnov funding from Région Île-de-France.

*AMS classification.* 65M12, 74S10, 76M12, 35L65.

*Keywords.* Shallow water equations, Implicit and iterative kinetic solvers, Fully discrete entropy inequality, Well-balanced schemes, Hydrostatic reconstruction.
# IMPLICIT KINETIC SCHEMES FOR THE SAINT-VENANT SYSTEM

CHOUROUK EL HASSANIEH<fnref n="1" />,<fnref n="2" /> ID, MATHIEU RIGAL<fnref n="3" />,* ID AND JACQUES SAINTE-MARIE<fnref n="1" />,<fnref n="2" /> ID

**Abstract.** Explicit (in time) kinetic schemes applied to the nonlinear shallow water equations have been extensively studied in the past. The novelty of this paper is to investigate an implicit version of such methods in order to improve their stability properties. In the case of a flat bathymetry we obtain a fully implicit kinetic solver satisfying a discrete entropy inequality and keeping the water height non negative without any restriction on the time step. Remarkably, a simplified version of this nonlinear implicit scheme allows to express the update explicitly which we implement in practice. The case of varying bottoms is then dealt with through an iterative solver combined with the hydrostatic reconstruction technique. We show that this scheme preserves the water height non-negativity under a CFL condition and satisfies a discrete entropy inequality without error term, which is an improvement over its explicit version. An extension of the implicit and iterative methods to the two dimensional case is also discussed. Finally we perform some numerical validations underlining the advantages and the computational cost of our strategy.

**Mathematics Subject Classification.** 65M12, 74S10, 76M12, 35L65.

Received May 13, 2024. Accepted May 16, 2025.

## 1. INTRODUCTION

### 1.1. General motivation and goals

Free surface flows are a subject of great importance whose understanding is at the core of many present challenges. These include the issue of safety related to submersion and tsunami waves affecting urbanized areas, rogue waves that can be a hazard for sailing, or the transport in rivers or lakes. Another relevant topic is that of renewable energy generated from hydraulic dams or sea buoys, and the impact of these systems on the aquatic fauna to give only a few examples. To describe such complex fluid flows, there exists numerous mathematical models, and when restricting to shallow areas where it is reasonable to neglect dispersive effects but not nonlinear ones the Saint-Venant equations [18, 28], also known as the shallow water equations, are often preferred to the incompressible Navier-Stokes system with free surface for practical applications. In fact, although the latter model is more accurate, it is also significantly more complex to study theoretically and numerically. On the

<footnote>*Keywords and phrases.* Shallow water equations, implicit and iterative kinetic solvers, fully discrete entropy inequality, well-balanced schemes, hydrostatic reconstruction.</footnote>

<footnote n="1">¹ Centre Inria de Sorbonne université, 2 rue Simone Iff, CS 42112, 75589 Paris Cedex 12, France.</footnote>

<footnote n="2">² Sorbonne Université, Univ. Paris Diderot, SPC, CNRS, LJLL, F-75005 Paris, France.</footnote>

<footnote n="3">³ Institut de Mathématiques de Bordeaux, équipe EDP et physique-mathématique, Université de Bordeaux, 351 cours de la Libération, F-33405 Talence, France.</footnote>

<footnote>*Corresponding author: mathieu.rigal@math.u-bordeaux.fr</footnote>
# Contents

<table><tbody><tr><td><strong>1 Introduction</strong></td><td><strong>2</strong></td></tr><tr><td>1.1 General motivation and goals</td><td>2</td></tr><tr><td>1.2 State of the art and outline of the paper</td><td>3</td></tr><tr><td>1.3 The Saint-Venant system and its kinetic interpretation</td><td>4</td></tr><tr><td>1.4 Explicit kinetic schemes for the Saint-Venant system</td><td>6</td></tr><tr><td><strong>2 A fully implicit kinetic scheme over a flat bathymetry</strong></td><td><strong>8</strong></td></tr><tr><td>2.1 Kinetic discretization and properties</td><td>8</td></tr><tr><td>2.2 Practical computation of the implicit kinetic update</td><td>11</td></tr><tr><td>2.3 Macroscopic implicit scheme</td><td>12</td></tr><tr><td>2.4 Boundary conditions</td><td>13</td></tr><tr><td>2.5 Implementation and computational costs</td><td>15</td></tr><tr><td><strong>3 An iterative resolution scheme</strong></td><td><strong>18</strong></td></tr><tr><td>3.1 Case of a flat topography</td><td>18</td></tr><tr><td>3.2 Case of a non flat topography</td><td>24</td></tr><tr><td><strong>4 Towards the two dimensional Saint-Venant system</strong></td><td><strong>28</strong></td></tr><tr><td>4.1 The two dimensional system and its kinetic representation</td><td>29</td></tr><tr><td>4.2 Fully implicit kinetic scheme over a flat bathymetry</td><td>29</td></tr><tr><td>4.3 Iterative kinetic scheme with hydrostatic reconstruction</td><td>31</td></tr><tr><td><strong>5 Numerical simulations</strong></td><td><strong>32</strong></td></tr><tr><td>5.1 Fully implicit kinetic scheme in the flat bottom case</td><td>32</td></tr><tr><td>5.2 Iterative kinetic scheme with hydrostatic reconstruction</td><td>35</td></tr><tr><td>5.3 Two dimensional simulation</td><td>36</td></tr><tr><td><strong>A Expression of the numerical fluxes</strong></td><td><strong>39</strong></td></tr><tr><td><strong>B Computations of the fluxes involving the boundary conditions</strong></td><td><strong>42</strong></td></tr></tbody></table>

# 1 Introduction

## 1.1 General motivation and goals

Free surface flows are a subject of great importance whose understanding is at the core of many present challenges. These include the issue of safety related to submersion and tsunami waves affecting urbanized areas, rogue waves that can be a hazard for sailing, or the transport in rivers or lakes. Another relevant topic is that of renewable energy generated from hydraulic dams or sea buoys, and the impact of these systems on the aquatic fauna to give only a few examples. To describe such complex fluid flows, there exists numerous mathematical models, and when restricting to shallow areas where it is reasonable to neglect dispersive effects but not nonlinear ones the Saint-Venant equations [27, 18], also known as the shallow water equations, are often preferred to the incompressible Navier-Stokes system with free surface for practical applications. In fact, although the latter model is more accurate, it is also significantly more complex to study theoretically and numerically. On the other hand the Saint-Venant system has a reduced complexity since it is a vertically averaged nonlinear model belonging to the class of hyperbolic systems of balance laws. Despite being simpler, it can be successfully used to approximate various geophysical flows such as rivers, coastal areas, and oceans when completed with a Coriolis term, and granular flows when completed with friction terms.

The analytical solutions of the Saint-Venant equations satisfy important properties, among which the conservation of the total water and of the total momentum when the bathymetry is flat, the positivity of the water height, or the existence of the lake at rest steady state that becomes a nontrivial equilibrium in presence of a varying bathymetry, which is accounted for through a source term in the
other hand the Saint-Venant system has a reduced complexity since it is a vertically averaged nonlinear model belonging to the class of hyperbolic systems of balance laws. Despite being simpler, it can be successfully used to approximate various geophysical flows such as rivers, coastal areas, and oceans when completed with a Coriolis term, and granular flows when completed with friction terms.

The analytical solutions of the Saint-Venant equations satisfy important properties, among which the conservation of the total water and of the total momentum when the bathymetry is flat, the positivity of the water height, or the existence of the lake at rest steady state that becomes a nontrivial equilibrium in presence of a varying bathymetry, which is accounted for through a source term in the momentum equation. Because the solutions of the Saint-Venant system are non-unique, it is usual to select the one obtained through a viscosity perturbation, and which satisfies additional entropy inequalities that can be interpreted as the dissipation of some energy. For the derivation of a numerical scheme approximating the Saint-Venant system, the main difficulty is to construct a structure preserving numerical method, that is to say a method able to preserve the aforementioned properties at the discrete level so that the numerical approximation is qualitatively relevant. The present paper is dedicated to this issue. More specifically we focus on a kinetic solver, which provides a good framework regarding positivity and entropy stability, and investigate the advantages of combining it with an implicit time integrator compared to the more common explicit approach.

There are several contributions in this work. In the case of the one dimensional Saint-Venant system with flat bathymetry, we propose a new, fully implicit kinetic scheme which is structure preserving without constraint on the time step. The case of varying bathymetries is dealt with by the mean of the hydrostatic reconstruction technique together with a fixed point method approximating the implicit update. Unlike its explicit version which does not always dissipate the entropy, our iterative scheme is shown to be structure preserving under a CFL condition. We then point to the possibility of extending these two numerical methods to the two dimensional Saint-Venant equations, and numerical simulations are performed to evaluate the interest of the proposed approaches.

## 1.2. State of the art and outline of the paper

The derivation of an efficient, robust and stable numerical scheme for the Saint-Venant system has received an extensive coverage, we refer the reader to [11, 20, 23, 30] and references therein. In particular, kinetic schemes approximating macroscopic systems of conservation laws have become popular over the last decades. Such methods consist to substitute the macroscopic system, in our case the Saint-Venant equations, by a scalar kinetic equation featuring a transport term with a collision operator, and which is easier to discretize. In [9], Bouchut developed a general theory to construct kinetic representations involving a BGK collision operator and admitting a family of kinetic entropies. As we will recall later, kinetic entropies are an important tool allowing one to recover the dissipation mechanism of the macroscopic entropy at the kinetic level. Once a kinetic representation with a kinetic entropy is known, it is then possible to design a numerical scheme satisfying a discrete entropy inequality as explained in [10].

In presence of a source term, such as the one induced by a varying bathymetry in the Saint-Venant equations, one of the challenges involves the design of a well-balanced scheme, that is to say a discretization able to preserve some characteristic stationary solutions. In [27], Perthame and Simeoni proposed a well-balanced kinetic scheme where the bathymetry is seen as a stair-shaped parametrization against which fluid particles can be reflected. A drawback of this is that the associated macroscopic update cannot be written explicitly, and instead some integrals have to be approximated by a quadrature. They also provided a Maxwellian with compact support that differs from the usual Gaussian distribution, and that has been used later in [6, 14, 21]. This allows the scheme to be positive and entropy stable under a CFL condition depending on the size of the support. In a recent work [6], some of the authors have proposed a kinetic scheme coupled with the hydrostatic reconstruction technique introduced in [4, 5] for the numerical treatment of the topography source term. With this approach an explicit writing of the macroscopic update is available, however an error term with positive sign is present in the discrete entropy inequality; as a result the entropy can increase in some test cases. Despite this, Bouchut and Lhebrard have proved the convergence of the kinetic hydrostatic reconstruction scheme towards entropy
momentum equation. Because the solutions of the Saint-Venant system are non-unique, it is usual
to select the one obtained through a viscosity perturbation, and which satisfies additional entropy
inequalities that can be interpreted as the dissipation of some energy. For the derivation of a numer-
ical scheme approximating the Saint-Venant system, the main difficulty is to construct a structure
preserving numerical method, that is to say a method able to preserve the aforementioned properties
at the discrete level so that the numerical approximation is qualitatively relevant. The present paper
is dedicated to this issue. More specifically we focus on a kinetic solver, which provides a good frame-
work regarding positivity and entropy stability, and investigate the advantages of combining it with
an implicit time integrator compared to the more common explicit approach.

There are several contributions in this work. In the case of the one dimensional Saint-Venant system with flat bathymetry, we propose a new, fully implicit kinetic scheme which is structure preserving without constraint on the time step. The case of varying bathymetries is dealt with by the mean of the hydrostatic reconstruction technique together with a fixed point method approximating the implicit update. Unlike its explicit version which does not always dissipate the entropy, our iterative scheme is shown to be structure preserving under a CFL condition. We then point to the possibility of extending these two numerical methods to the two dimensional Saint-Venant equations, and numerical simulations are performed to evaluate the interest of the proposed approaches.

1.2 State of the art and outline of the paper

The derivation of an efficient, robust and stable numerical scheme for the Saint-Venant system has re-
ceived an extensive coverage, we refer the reader to [10, 22, 19, 29] and references therein. In particular,
kinetic schemes approximating macroscopic systems of conservation laws have become popular over the
last decades. Such methods consist to substitute the macroscopic system, in our case the Saint-Venant
equations, by a scalar kinetic equation featuring a transport term with a collision operator, and which
is easier to discretize. In [9], Bouchut developed a general theory to construct kinetic representations
involving a BGK collision operator and admitting a family of kinetic entropies. As we will recall later,
kinetic entropies are an important tool allowing one to recover the dissipation mechanism of the macro-
scopic entropy at the kinetic level. Once a kinetic representation with a kinetic entropy is known, it is
then possible to design a numerical scheme satisfying a discrete entropy inequality as explained in [12].

In presence of a source term, such as the one induced by a varying bathymetry in the Saint-Venant
equations, one of the challenges involves the design of a well-balanced scheme, that is to say a dis-
cretization able to preserve some characteristic stationary solutions. In [24], Perthame and Simeoni
proposed a well-balanced kinetic scheme where the bathymetry is seen as a stair-shaped parametriza-
tion against which fluid particles can be reflected. A drawback of this is that the associated macroscopic
update cannot be written explicitly, and instead some integrals have to be approximated by a quadra-
ture. They also provided a Maxwellian with compact support that differs from the usual Gaussian
distribution, and that has been used later in [20, 14, 4]. This allows the scheme to be positive and
entropy stable under a CFL condition depending on the size of the support. In a recent work [4], some
of the authors have proposed a kinetic scheme coupled with the hydrostatic reconstruction technique
introduced in [3, 6] for the numerical treatment of the topography source term. With this approach
an explicit writing of the macroscopic update is available, however an error term with positive sign
is present in the discrete entropy inequality; as a result the entropy can increase in some test cases.
Despite this, Bouchut and Lhebrard have proved the convergence of the kinetic hydrostatic reconstruc-
tion scheme towards entropy solutions of the Saint-Venant system, see [11]. The implicit version of the
scheme proposed in [4] gives a stronger stability result where the discrete entropy inequality is free of
the positive error term and the entropy is always dissipated.

The capability of implicit time integrators to improve the stability properties of numerical schemes
has been studied in the last years. Explicit finite volume approaches for the approximation of conser-
vation laws have to deal with a CFL constraint that can be very restrictive for some applications where
large time scales and significant wave velocities have to be considered. This is for instance the case
in the low Froude regime, where the surface gravity waves travel at a much larger velocity than the
fluid particles. A solution consists in implicit-explicit schemes [21, 8, 28, 2] where the fast dynamics
are approximated implicitly, allowing the use of larger time steps while keeping the scheme stable.
solutions of the Saint-Venant system, see [12]. The implicit version of the scheme proposed in [6] gives a stronger stability result where the discrete entropy inequality is free of the positive error term and the entropy is always dissipated.

The capability of implicit time integrators to improve the stability properties of numerical schemes has been studied in the last years. Explicit finite volume approaches for the approximation of conservation laws have to deal with a CFL constraint that can be very restrictive for some applications where large time scales and significant wave velocities have to be considered. This is for instance the case in the low Froude regime, where the surface gravity waves travel at a much larger velocity than the fluid particles. A solution consists in implicit-explicit schemes [2, 8, 22, 29] where the fast dynamics are approximated implicitly, allowing the use of larger time steps while keeping the scheme stable. Another important question is the computational costs, in the sense that for explicit schemes, the CFL constraint implies small time steps whereas for implicit schemes, the computation of the numerical fluxes can be costly. Indeed in the explicit setting, the numerical fluxes at an interface depend on the value of the variables at the two neighboring cells whereas in the implicit context, the numerical fluxes depend on the value of the variables of all the cells. In [16], a fully implicit lattice Boltzmann scheme was proposed, with the advantage of having the computational cost of an explicit approach thanks to a sweeping algorithm, however to our knowledge this method lacks a fully discrete entropy inequality. The novelties of the paper are:

- a fully implicit numerical scheme for the Saint-Venant system is detailed and analyzed;

- the proposed numerical scheme satisfies a fully discrete entropy inequality;

- whereas, in general, an implicit scheme often requires to invert an operator – typically a matrix – at each time step, our fully implicit kinetic scheme offers a very favorable context since we have an explicit expression of the inverse of the operator. Hence, one can hardly imagine a truly implicit entropy stable scheme for the Saint-Venant system with a lower computational cost than the kinetic solver proposed here.

This paper is organized as follows. In the remainder of this section we recall the formulation of the Saint-Venant system, its kinetic description and the framework of its numerical approximation in the context of a kinetic solver. Especially, the notion of kinetic entropy is detailed. Then in Section 2, the implicit kinetic scheme for the Saint-Venant system with flat topography is proposed and studied in the one dimensional case. An iterative kinetic scheme is proposed in Section 3 where the topography can be taken into account through the hydrostatic reconstruction technique. In Section 4, we briefly discuss the extension of the previous approaches to the two dimensional Saint-Venant system. Finally in Section 5, numerical examples are given to evaluate the interest of the proposed schemes.

## 1.3. The Saint-Venant system and its kinetic interpretation

The classical Saint-Venant system for shallow water flows describes the evolution of the height of water $h(t, x) \ge 0$, and the water velocity $u(t, x) \in \mathbb{R}$ ($x$ denotes a coordinate in the horizontal direction) over a slowly varying topography $z(x)$, and reads

$$ 
\begin{aligned}
\partial_t h + \partial_x(hu) &= 0, \\
\partial_t(hu) + \partial_x(hu^2 + g\frac{h^2}{2}) + gh\partial_x z &= 0,
\end{aligned} 
\quad (1.1) $$

where $g > 0$ is the gravity constant. This system is completed with an entropy (energy) inequality

$$ \partial_t \left( h \frac{u^2}{2} + g \frac{h^2}{2} + ghz \right) + \partial_x \left( \left( h \frac{u^2}{2} + gh^2 + ghz \right) u \right) \le 0. \quad (1.2) $$

We shall denote $U = (h, hu)^T$ and

$$ \eta(U) = h \frac{u^2}{2} + g \frac{h^2}{2}, \quad G(U) = \left( h \frac{u^2}{2} + gh^2 \right) u, \quad (1.3) $$
Another important question is the computational costs, in the sense that for explicit schemes, the CFL constraint implies small time steps whereas for implicit schemes, the computation of the numerical fluxes can be costly. Indeed in the explicit setting, the numerical fluxes at an interface depend on the value of the variables at the two neighboring cells whereas in the implicit context, the numerical fluxes depend on the value of the variables of all the cells. In [16], a fully implicit lattice Boltzmann scheme was proposed, with the advantage of having the computational cost of an explicit approach thanks to a sweeping algorithm. However this method requires the boundary conditions to be non periodic, and no discrete entropy inequality was provided. The novelties of the paper are:

* a fully implicit numerical scheme for the Saint-Venant system is detailed and analysed;

* the proposed numerical scheme satisfies a fully discrete entropy inequality;

* whereas, in general, an implicit scheme often requires to invert an operator – typically a matrix – at each time step, our fully implicit kinetic scheme offers a very favorable context since we have an explicit expression of the inverse of the operator. Hence, one can hardly imagine a truly implicit entropy stable scheme for the Saint-Venant system with a lower computational cost than a kinetic solver proposed here.

This paper is organized as follows. In the remainder of this section we recall the formulation of the Saint-Venant system, its kinetic description and the framework of its numerical approximation in the context of a kinetic solver. Especially, the notion of kinetic entropy is detailed. Then in Section 2, the implicit kinetic scheme for the Saint-Venant system with flat topography is proposed and studied in the one dimensional case. An iterative kinetic scheme is proposed in Section 3 where the topography can be taken into account through the hydrostatic reconstruction technique. In Section 4, we briefly discuss the extension of the previous approaches to the two dimensional Saint-Venant system. Finally in Section 5, numerical examples are given to evaluate the interest of the proposed schemes.

## 1.3 The Saint-Venant system and its kinetic interpretation

The classical Saint-Venant system for shallow water flows describes the evolution of the height of water $h(t, x) \ge 0$, and the water velocity $u(t, x) \in \mathbb{R}$ ($x$ denotes a coordinate in the horizontal direction) over a slowly varying topography $z(x)$, and reads

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

where $g > 0$ is the gravity constant. This system is completed with an entropy (energy) inequality

$$ (1.2) \qquad \partial_t \left( h \frac{u^2}{2} + g \frac{h^2}{2} + ghz \right) + \partial_x \left( \left( h \frac{u^2}{2} + gh^2 + ghz \right) u \right) \le 0. $$

We shall denote $U = (h, hu)^T$ and

$$ (1.3) \qquad \eta(U) = h \frac{u^2}{2} + g \frac{h^2}{2}, \qquad G(U) = \left(h \frac{u^2}{2} + gh^2\right) u, $$

the entropy and entropy fluxes without topography. The reader can refer to [4] and references therein for a complete presentation of the description of the Saint-Venant system.

The classical kinetic Maxwellian (see e.g. [24]) is given by

$$ (1.4) \qquad M(U, \xi) = \frac{1}{g\pi} \left( 2gh - (\xi - u)^2 \right)_+^{1/2}, $$

where $\xi \in \mathbb{R}$ and $x_+ = \max(0, x)$ for any $x \in \mathbb{R}$. It satisfies the following moment relations,

$$ (1.5) \qquad \begin{aligned} \int_{\mathbb{R}} M(U, \xi) \, d\xi &= h, & \int_{\mathbb{R}} \xi M(U, \xi) \, d\xi &= hu, \\ \int_{\mathbb{R}} \xi^2 M(U, \xi) \, d\xi &= hu^2 + g \frac{h^2}{2}. \end{aligned} $$
the entropy and entropy fluxes without topography. The reader can refer to [6] and references therein for a complete presentation of the description of the Saint-Venant system.

The classical kinetic Maxwellian (see e.g. [27]) is given by

$$M(U, \xi) = \frac{1}{g\pi} (2gh - (\xi - u)^2)_+^{1/2}, \qquad (1.4)$$

where $\xi \in \mathbb{R}$ and $x_+ \equiv \max(0, x)$ for any $x \in \mathbb{R}$. It satisfies the following moment relations,

$$\int_{\mathbb{R}} M(U, \xi) \, d\xi = h, \qquad \int_{\mathbb{R}} \xi M(U, \xi) \, d\xi = hu, \\ \int_{\mathbb{R}} \xi^2 M(U, \xi) \, d\xi = hu^2 + g \frac{h^2}{2}. \qquad (1.5)$$

These definitions allow us to obtain a kinetic representation of the Saint-Venant system.

**Lemma 1.1.** If the topography $z(x)$ is Lipschitz continuous, the pair of functions $(h, hu)$ is a weak solution to the Saint-Venant system (1.1) if and only if $M(U, \xi)$ satisfies the kinetic equation

$$\partial_t M + \xi \partial_x M - g(\partial_x z) \partial_\xi M = Q, \qquad (1.6)$$

for some “collision term” $Q(t, x, \xi)$ that satisfies, for a.e. $(t, x)$,

$$\int_{\mathbb{R}} Q \, d\xi = \int_{\mathbb{R}} \xi Q \, d\xi = 0. \qquad (1.7)$$

*Proof.* If (1.6) and (1.7) are satisfied, we can multiply (1.6) by $(1, \xi)^T$, and integrate with respect to $\xi$. Using (1.5) and (1.7) and integrating by parts the term in $\partial_\xi M$, we obtain (1.1). Conversely, if $(h, hu)$ is a weak solution to (1.1), just define $Q$ by (1.6); it will satisfy (1.7) according to the same computations. $\square$

The standard way to use Lemma 1.1 is to write a kinetic relaxation equation [9, 10, 15, 24, 25], like

$$\partial_t f + \xi \partial_x f - g(\partial_x z) \partial_\xi f = \frac{M - f}{\epsilon}, \qquad (1.8)$$

where the distribution $f(t, x, \xi)$ is positive, where $M = M(U, \xi)$ with $U(t, x) = \int (1, \xi)^T f(t, x, \xi) d\xi$, and where $\epsilon > 0$ is a relaxation time. In the limit $\epsilon \to 0$ we recover formally the formulation (1.6), (1.7). We refer to [9] for general considerations on such kinetic relaxation models without topography, the case with topography being introduced in [27]. Note that the notion of kinetic representation as (1.6), (1.7) differs from the so called kinetic formulations where a large set of entropies is involved, see [26]. For systems of conservation laws, these kinetic formulations include non-adptiveive terms that prevent from writing down simple approximations. In general, kinetic relaxation approximations can be compatible with just a single entropy. Nevertheless this is enough for proving the convergence as $\epsilon \to 0$, see [7].

The interest of the particular form (1.4) lies in its link with a kinetic entropy. Consider the kinetic entropy,

$$H(f, \xi, z) = \frac{\xi^2}{2} f + \frac{g^2 \pi^2}{6} f^3 + gz f, \qquad (1.9)$$

where $f \ge 0$, $\xi \in \mathbb{R}$ and $z \in \mathbb{R}$, and its version without topography

$$H_0(f, \xi) = \frac{\xi^2}{2} f + \frac{g^2 \pi^2}{6} f^3. \qquad (1.10)$$
These definitions allow us to obtain a kinetic representation of the Saint-Venant system.

**Lemma 1.1** If the topography $z(x)$ is Lipschitz continuous, the pair of functions $(h, hu)$ is a weak solution to the Saint-Venant system (1.1) if and only if $M(U, \xi)$ satisfies the kinetic equation

$$ (1.6) \qquad \partial_t M + \xi \partial_x M - g(\partial_x z) \partial_\xi M = Q, $$

for some “collision term” $Q(t, x, \xi)$ that satisfies, for a.e. $(t, x)$,

$$ (1.7) \qquad \int_{\mathbb{R}} Q d\xi = \int_{\mathbb{R}} \xi Q d\xi = 0. $$

**Proof.** If (1.6) and (1.7) are satisfied, we can multiply (1.6) by $(1, \xi)^T$, and integrate with respect to $\xi$. Using (1.5) and (1.7) and integrating by parts the term in $\partial_\xi M$, we obtain (1.1). Conversely, if $(h, hu)$ is a weak solution to (1.1), just define $Q$ by (1.6); it will satisfy (1.7) according to the same computations. $\square$

The standard way to use Lemma 1.1 is to write a kinetic relaxation equation [23, 26, 15, 9, 12], like

$$ (1.8) \qquad \partial_t f + \xi \partial_x f - g(\partial_x z) \partial_\xi f = \frac{M - f}{\epsilon}, $$

where the distribution $f(t, x, \xi)$ is positive, where $M = M(U, \xi)$ with $U(t, x) = \int (1, \xi)^T f(t, x, \xi) d\xi$, and where $\epsilon > 0$ is a relaxation time. In the limit $\epsilon \to 0$ we recover formally the formulation (1.6), (1.7). We refer to [9] for general considerations on such kinetic relaxation models without topography, the case with topography being introduced in [24]. Note that the notion of kinetic representation as (1.6), (1.7) differs from the so called kinetic formulations where a large set of entropies is involved, see [25]. For systems of conservation laws, these kinetic formulations include non-advective terms that prevent from writing down simple approximations. In general, kinetic relaxation approximations can be compatible with just a single entropy. Nevertheless this is enough for proving the convergence as $\epsilon \to 0$, see [7].

The interest of the particular form (1.4) lies in its link with a kinetic entropy. Consider the kinetic entropy,

$$ (1.9) \qquad H(f, \xi, z) = \frac{\xi^2}{2} f + \frac{g^2 \pi^2}{6} f^3 + gz f, $$

where $f \ge 0$, $\xi \in \mathbb{R}$ and $z \in \mathbb{R}$, and its version without topography

$$ (1.10) \qquad H_0(f, \xi) = \frac{\xi^2}{2} f + \frac{g^2 \pi^2}{6} f^3. $$

Then one can check the relations

$$ (1.11) \qquad \int_{\mathbb{R}} H(M(U, \xi), \xi, z) d\xi = \eta(U) + ghz, $$

$$ (1.12) \qquad \int_{\mathbb{R}} \xi H(M(U, \xi), \xi, z) d\xi = G(U) + ghzu. $$

One has the following entropy relations.

**Lemma 1.2** Let $f(\xi) \ge 0$ satisfy $\int f(\xi) d\xi = h$ and $\int \xi f(\xi) d\xi = hu$ (assumed finite). The half-disk Maxwellian (1.4) satisfies the three properties below.

(i) For any $\xi \in \mathbb{R}$ one has

$$ (1.13) \qquad H_0(f, \xi) \ge H_0(M(U, \xi), \xi) + \eta'(U) \cdot \binom{1}{\xi} (f - M(U, \xi)). $$
Then one can check the relations

$$ \int_{\mathbb{R}} H(M(U, \xi), \xi, z) \, d\xi = \eta(U) + ghz, \qquad (1.11) $$

$$ \int_{\mathbb{R}} \xi H(M(U, \xi), \xi, z) \, d\xi = G(U) + ghzu. \qquad (1.12) $$

One has the following entropy relations.

**Lemma 1.2.** Let $f(\xi) \ge 0$ satisfy $\int f(\xi) \, d\xi = h$ and $\int \xi f(\xi) \, d\xi = hu$ (assumed finite). The half-disk Maxwellian (1.4) satisfies the two properties below.

(i) For any $\xi \in \mathbb{R}$ one has

$$ H_0(f, \xi) \ge H_0(M(U, \xi), \xi) + \eta'(U) \cdot \binom{1}{\xi} (f - M(U, \xi)). \qquad (1.13) $$

(ii) For $U = (h, hu)^T$ there holds the inequality

$$ \eta(U) = \int_{\mathbb{R}} H_0(M(U, \xi), \xi) \, d\xi \le \int_{\mathbb{R}} H_0(f(\xi), \xi) \, d\xi. \qquad (1.14) $$

*Proof of Lemma 1.2.* (i) A proof of inequality (1.13) is given in [6] and relies on the following relation

$$ \partial_f H_0(M(U, \xi), \xi) = \begin{cases} \eta'(U) \cdot \binom{1}{\xi} & \text{if } \xi \in \text{supp } M(U, \cdot) \\ \ge \eta'(U) \cdot \binom{1}{\xi} & \text{otherwise,} \end{cases} \qquad (1.15) $$

which we will use in Section 3 to obtain discrete entropy inequalities. The property (ii) was proved by Perthame and Simeoni in [27]. It is simply recovered from (i) by integrating (1.13) over $\xi \in \mathbb{R}$ and using the fact that $f - M(U, \xi)$ is a collision term. $\square$

One can recover the macroscopic entropy inequality (1.2) thanks to Lemma 1.2. For this, assuming $H$ to be at least $C^1$, we multiply the kinetic relaxation equation (1.8) by $\partial_f H(f, \xi, z)$ to get

$$ \partial_t H(f, \xi, z) + \xi(\partial_x H(f, \xi, z) - \partial_z H(f, \xi, z)\partial_x z) - g(\partial_x z)(\partial_\xi f)\partial_f H(f, \xi, z) = \partial_f H(f, \xi, z) \frac{M-f}{\epsilon}. $$

By convexity of $H$ we have the following inequality

$$ \partial_t H(f, \xi, z) + \xi \partial_x H(f, \xi, z) - \partial_x z \left( \xi \partial_z H(f, \xi, z) + g(\partial_\xi f) \partial_f H(f, \xi, z) \right) \le \frac{H(M, \xi, z) - H(f, \xi, z)}{\epsilon}, $$

and using $\partial_z H(f, \xi, z) = gf$ and $\partial_f H(f, \xi, z) = \frac{1}{2}\xi^2 + \frac{1}{2}(g\pi f)^2 + gz$, this yields

$$ \partial_t H(f, \xi, z) + \xi \partial_x H(f, \xi, z) - g \partial_x z \left( \xi f + \left( \frac{\xi^2}{2} + gz \right) \partial_\xi f + \frac{g^2 \pi^2}{6} \partial_\xi f^3 \right) \le \frac{H(M, \xi, z) - H(f, \xi, z)}{\epsilon}. \qquad (1.16) $$

In the general setting of weak solutions, one cannot assume $H$ to be smooth, in which case the above steps cannot be undertaken. Rather, the inequality (1.16) should be seen as an entropy condition. Notice that in the
(ii) For $U = (h, hu)^T$ there holds the inequality

$$ (1.14) \qquad \eta(U) = \int_{\mathbb{R}} H_0(M(U, \xi), \xi) d\xi \le \int_{\mathbb{R}} H_0(f(\xi), \xi) d\xi. $$

(iii) There holds the kinetic entropy equality

$$ \partial_t H(M, \xi, z) + \partial_x(\xi H(M, \xi, z)) - g(\partial_x z)\partial_\xi H(M, \xi, z) = \partial_f H(M, \xi, z)Q, $$

and under the additional assumption $Q(t, x, \xi) \le 0$ for all $\xi \notin \text{supp } M(U(t, x), \cdot)$ we have the macroscopic inequality

$$ \partial_t \int_{\mathbb{R}} H(M, \xi, z) d\xi + \partial_x \int_{\mathbb{R}} \xi H(M, \xi, z) d\xi \le 0, $$

which is exactly (1.2) owing to the relations (1.11)-(1.12).

**Proof of Lemma 1.2.** (i) A proof of inequality (1.13) is given in [4] and relies on the following relation

$$ (1.15) \qquad \partial_f H_0(M(U, \xi), \xi) = \begin{cases} \eta'(U) \cdot \begin{pmatrix} 1 \\ \xi \end{pmatrix} & \text{if } \xi \in \text{supp } M(U, \cdot) \\ \ge \eta'(U) \cdot \begin{pmatrix} 1 \\ \xi \end{pmatrix} & \text{otherwise} \end{cases}, $$

which we will use in Section 3 to obtain discrete entropy inequalities. The property (ii) was proved by Perthame and Simeoni in [24]. It is simply recovered from (i) by integrating (1.13) over $\xi \in \mathbb{R}$ and using the fact that $f - M(U, \xi)$ is a collision term. In order to prove (iii), we multiply (1.6) by $\partial_f H(M, \xi, z)$, we bound the right hand side using (1.15) yielding

$$ \partial_f H(M, \xi)Q = (\partial_f H_0(M, \xi) + gz)Q \le (\eta'(U) \cdot \begin{pmatrix} 1 \\ \xi \end{pmatrix} + gz)Q, $$

and owing to (1.7) an integration in $\xi$ of the obtained inequality gives the result. $\square$

## 1.4 Explicit kinetic schemes for the Saint-Venant system

We intend to approximate the solution $U(t, x)$ of the system (1.1) for $x \in \mathbb{R}$ and $t \ge 0$ by discrete values $U_i^n$, with $i \in \mathbb{Z}$ the spatial index and $n \in \mathbb{N}$ the time index. For this, we consider a grid of real points $(x_{i+1/2})_{i \in \mathbb{Z}}$ and discrete times $(t^n)_{n \in \mathbb{N}}$ satisfying

$$ x_{i+1/2} - x_{i-1/2} = \Delta x_i, \qquad t^{n+1} - t^n = \Delta t^n, $$

with $\Delta x_i > 0$ the spatial step and $\Delta t^n > 0$ the time step. Defining the cell $C_i = (x_{i-1/2}, x_{i+1/2})$ of length $\Delta x_i$, we shall consider $U^n(x)$ the piecewise constant function approximating the solution at time $t^n$ and corresponding to

$$ (1.16) \qquad U^n(x) = U_i^n, \quad \text{for } x \in C_i. $$

A finite volume scheme for (1.1) is a formula providing a relation between $U_i^{n+1}$ and the set $\{U_j^n\}_{j \in \mathbb{Z}}$ of the form

$$ (1.17) \qquad U_i^{n+1} = U_i^n - \sigma_i(F_{i+1/2-} - F_{i-1/2+}), $$

where $\sigma_i = \Delta t^n / \Delta x_i$, and where $F_{i\pm 1/2\mp}$ are two-points numerical fluxes given by

$$ (1.18) \qquad F_{i+1/2-} = \mathcal{F}_l(U_i^{n+p}, U_{i+1}^{n+p}, z_{i+1} - z_i), \qquad F_{i+1/2+} = \mathcal{F}_r(U_i^{n+p}, U_{i+1}^{n+p}, z_{i+1} - z_i), $$
left hand side of (1.16), the term in factor of $g\partial_x z$ is macroscopically zero, therefore integrating this inequality w.r.t. $\xi$ results in

$$ \partial_t \int_{\mathbb{R}} H(f, \xi, z) \, d\xi + \partial_x \int_{\mathbb{R}} \xi H(f, \xi, z) \, d\xi \le \frac{1}{\epsilon} \int_{\mathbb{R}} (H_0(M, \xi) + gzM - H_0(f, \xi) - gzf) \, d\xi \le 0, $$

where we used the second point of Lemma 1.2 to get the second inequality. Finally, considering that the solution $f$ of (1.8) formally converges to $M$ as $\epsilon \to 0$, together with the relations (1.11) and (1.12), one obtains the desired macroscopic entropy inequality (1.2).

## 1.4. Explicit kinetic schemes for the Saint-Venant system

We intend to approximate the solution $U(t, x)$ of the system (1.1) for $x \in \mathbb{R}$ and $t \ge 0$ by discrete values $U_i^n$, with $i \in \mathbb{Z}$ the spatial index and $n \in \mathbb{N}$ the time index. For this, we consider a grid of real points $(x_{i+1/2})_{i \in \mathbb{Z}}$ and discrete times $(t^n)_{n \in \mathbb{N}}$ satisfying

$$ x_{i+1/2} - x_{i-1/2} = \Delta x_i, \quad t^{n+1} - t^n = \Delta t^n, $$

with $\Delta x_i > 0$ the spatial step and $\Delta t^n > 0$ the time step. Defining the cell $C_i = (x_{i-1/2}, x_{i+1/2})$ of length $\Delta x_i$, we shall consider $U^n(x)$ the piecewise constant function approximating the solution at time $t^n$ and corresponding to

$$ U^n(x) = U_i^n, \quad \text{for } x \in C_i. \tag{1.17} $$

A finite volume scheme for (1.1) is a formula providing a relation between $U_i^{n+1}$ and the set $\{U_j^n\}_{j \in \mathbb{Z}}$ of the form

$$ U_i^{n+1} = U_i^n - \sigma_i (F_{i+1/2-} - F_{i-1/2+}), \tag{1.18} $$

where $\sigma_i = \Delta t^n / \Delta x_i$, and where $F_{i\pm 1/2\mp}$ are two-points numerical fluxes given by

$$ F_{i+1/2-} = \mathcal{F}_l(U_i^{n+p}, U_{i+1}^{n+p}, z_{i+1} - z_i), \quad F_{i+1/2+} = \mathcal{F}_r(U_i^{n+p}, U_{i+1}^{n+p}, z_{i+1} - z_i), \tag{1.19} $$

with $p \in \{0, 1\}$ and with $\mathcal{F}_l, \mathcal{F}_r$ some functions valued in $\mathbb{R}^2$ consistent with the analytic flux of the Saint-Venant system, see [11]. The value $p = 0$ classically corresponds to a first order explicit time scheme for solving (1.1) whereas $p = 1$ means an implicit time scheme. In this paper, we focus on the case $p = 1$, but we recall how to construct kinetic schemes in the classical explicit framework.

A kinetic scheme is a numerical discretization of the BGK-type kinetic equation (1.8) which then provides a numerical scheme for the Saint-Venant system using to the moment relations (1.5) after integration against $(1, \xi)$. In [27], a splitting between the transport term and the collision term is considered, requiring to solve the following problem at every time step

$$ \left\{ \begin{array}{l} \partial_t f + \xi \partial_x f - g(\partial_x z) \partial_\xi f = 0 \\ f(t^n, x, \xi) = M(U^n(x), \xi), \end{array} \right. \quad t \in (t^n, t^{n+1}). \tag{1.20} $$

Especially, we see that the initial kinetic density is projected onto the space of Maxwellians. Denoting $f_i^{n+1-}$ an approximation of the solution of (1.20) at time $t^{n+1}$ in cell $C_i$, the macroscopic update is given by

$$ U_i^{n+1} = \int_{\mathbb{R}} \left( \frac{1}{\xi} \right) f_i^{n+1-}(\xi) \, d\xi. \tag{1.21} $$

It remains to explain how to discretize (1.20). This is achieved by the mean of an upwind scheme accounting for the sign of the velocity $\xi$, and which for a flat bottom reads

$$ f_i^{n+1-} = M_i - \sigma_i \xi (1_{\xi>0} M_i + 1_{\xi<0} M_{i+1} - 1_{\xi<0} M_i - 1_{\xi>0} M_{i-1}). \tag{1.22} $$
with $p \in \{0, 1\}$ and with $\mathcal{F}_l, \mathcal{F}_r$ some functions valued in $\mathbb{R}^2$ consistent with the analytic flux of the Saint-Venant system, see [10]. The value $p = 0$ classically corresponds to a first order explicit time scheme for solving (1.1) whereas $p = 1$ means an implicit time scheme. In this paper, we focus on the case $p = 1$, but we recall how to construct kinetic schemes in the classical explicit framework.

A kinetic scheme is a numerical discretization of the BGK-type kinetic equation (1.8) which then provides a numerical scheme for the Saint-Venant system using to the moment relations (1.5) after integration against $(1, \xi)$. In [24], a splitting between the transport term and the collision term is considered, requiring to solve the following problem at every time step

$$ (1.19) \qquad \begin{cases} \partial_t f + \xi \partial_x f - g(\partial_x z) \partial_\xi f = 0 \\ f(t^n, x, \xi) = M(U^n(x), \xi) \end{cases}, \quad t \in (t^n, t^{n+1}). $$

Especially, we see that the initial kinetic density is projected onto the space of Maxwellians. Denoting $f_i^{n+1-}$ an approximation of the solution of (1.19) at time $t^{n+1}$ in cell $C_i$, the macroscopic update is given by

$$ (1.20) \qquad U_i^{n+1} = \int_{\mathbb{R}} \binom{1}{\xi} f_i^{n+1-}(\xi) \, d\xi. $$

It remains to explain how to discretize (1.19). This is achieved by the mean of an upwind scheme accounting for the sign of the velocity $\xi$, and which for a flat bottom reads

$$ (1.21) \qquad f_i^{n+1-} = M_i - \sigma_i \xi (1_{\xi>0} M_i + 1_{\xi<0} M_{i+1} - 1_{\xi<0} M_i - 1_{\xi>0} M_{i-1}). $$

In particular we see that (1.20)-(1.21) corresponds to the macroscopic update (1.17) with numerical fluxes (1.18) defined as

$$ \mathcal{F}_{l/r}(U_l, U_r, \Delta z) = \int_{\mathbb{R}} \xi \binom{1}{\xi} (1_{\xi>0} M(U_l, \xi) + 1_{\xi<0} M(U_r, \xi)) \, d\xi. $$

In [4], the scheme (1.21) was extended to varying bathymetries thanks to the hydrostatic reconstruction technique. The resulting discretization reads

$$ (1.22) \qquad f_i^{n+1-} = M_i - \sigma_i \xi (1_{\xi<0}(M_{i+1/2-} - M_{i-1/2+}) + 1_{\xi>0}(M_{i+1/2-} - M_{i-1/2-})) \\ \qquad \qquad + \sigma_i(\xi - u_i)(M_{i+1/2-} - M_{i-1/2+}), $$

where indices $i + 1/2 \pm$ denote reconstructed values at the interface $i + 1/2$ and whose expression will be detailed later (see Section 3.2). In [4], some of the authors have proved the following properties for the schemes (1.21) and (1.22).

**Proposition 1.3** (Explicit kinetic schemes [4]). Under the CFL condition $\sigma\xi_{\max} \le 1$, where $\xi_{\max} = \sup\{|\xi|, \xi \in \cup_{i \in \mathbb{Z}} \text{supp } M(U_i^n, \cdot)\}$ is finite if $(U_i^n)_i \in \ell^\infty(\mathbb{Z}; \mathbb{R}^2)$, we have that

1. the scheme (1.21) keeps $f_i^{n+1-}$ non-negative and satisfies a fully discrete entropy inequality of the form

$$ H_0(f_i^{n+1-}, \xi) \le H_0(M_i^n, \xi) - \sigma_i \xi (H_{0,i+1/2}^n(\xi) - H_{0,i-1/2}^n(\xi)), $$

2. the scheme (1.22) keeps $f_i^{n+1-}$ non-negative and satisfies a fully discrete entropy inequality of the form

$$ H(f_i^{n+1-}, z_i, \xi) \le H(M_i^n, z_i, \xi) - \sigma_i \xi (\tilde{H}_{i+1/2}^n(\xi) - \tilde{H}_{i-1/2}^n(\xi)) + D_i^n(\xi) + E_i^n(\xi), $$

where $D_i^n \le 0$ is a dissipation and $E_i^n \ge 0$ an error term that can dominate $D_i^n$ in many cases.
In particular we see that (1.21) and (1.22) corresponds to the macroscopic update (1.18) with numerical fluxes (1.19) defined as

$$ \mathcal{F}_{l/r}(U_l, U_r, \Delta z) = \int_{\mathbb{R}} \xi \left(\frac{1}{\xi}\right) \left(\mathbb{1}_{\xi>0} M(U_l, \xi) + \mathbb{1}_{\xi<0} M(U_r, \xi)\right) d\xi. $$

In [6], the scheme (1.22) was extended to varying bathymetries thanks to the hydrostatic reconstruction technique. The resulting discretization reads

$$ f_i^{n+1-} = M_i - \sigma_i \xi (\mathbb{1}_{\xi<0} (M_{i+1/2+} - M_{i-1/2+}) + \mathbb{1}_{\xi>0} (M_{i+1/2-} - M_{i-1/2-})) \\ + \sigma_i (\xi - u_i) (M_{i+1/2-} - M_{i-1/2+}), \qquad (1.23) $$

where indices $i + 1/2 \pm$ denote reconstructed values at the interface $i + 1/2$ and whose expression will be detailed later (see Sect. 3.2). In [6], some of the authors have proved the following properties for the schemes (1.22) and (1.23).

**Proposition 1.3** (Explicit kinetic schemes [6]). Under the CFL condition $\sigma\xi_{\max} \le 1$, where $\xi_{\max} = \sup\{|\xi|, \xi \in \cup_{i \in \mathbb{Z}} \text{supp } M(U_i^n, \cdot)\}$ is finite if $(U_i^n)_i \in \ell^\infty(\mathbb{Z}; \mathbb{R}^2)$, we have that

(1) the scheme (1.22) keeps $f_i^{n+1-}$ non-negative and satisfies a fully discrete entropy inequality of the form

$$ H_0(f_i^{n+1-}, \xi) \le H_0(M_i^n, \xi) - \sigma_i \xi (H_{0,i+1/2}^n(\xi) - H_{0,i-1/2}^n(\xi)), $$

(2) the scheme (1.23) keeps $f_i^{n+1-}$ non-negative and satisfies a fully discrete entropy inequality of the form

$$ H(f_i^{n+1-}, z_i, \xi) \le H(M_i^n, z_i, \xi) - \sigma_i \xi (\tilde{H}_{i+1/2}^n(\xi) - \tilde{H}_{i-1/2}^n(\xi)) + D_i^n(\xi) + E_i^n(\xi), $$

where $D_i^n \le 0$ is a dissipation and $E_i^n \ge 0$ an error term that can dominate $D_i^n$ in many cases.

## 2. A FULLY IMPLICIT KINETIC SCHEME OVER A FLAT BATHYMETRY

In this section we consider the model (1.1) with a flat topography, and propose an implicit kinetic scheme based on (1.20) and (1.21) i.e. an implicit version of (1.22). Since in practice it is not possible to approximate solutions over the whole space $x \in \mathbb{R}$, we shall work in a bounded domain and will need to enforce boundary conditions. First in Section 2.1 we present and study the implicit kinetic scheme; we show that it is unconditionally positive and entropy stable by using the kinetic entropy (1.10). Next in Section 2.2 we detail how to compute explicitly the update at the kinetic level, which is possible thanks to the particular structure of the underlying matrix. Then we consider the associated macroscopic scheme in Section 2.3. Unfortunately, it is not possible to compute the integral of the kinetic update against $(1, \xi)$, instead we propose a compromise consisting to replace the standard Maxwellian by a simpler one. The question of boundary conditions is then treated in Section 2.4, and finally we discuss the practical implementation and the computational costs in Section 2.5.

### 2.1. Kinetic discretization and properties

In the case of a flat topography, the kinetic scheme is a flux vector splitting scheme [10]. We denote $P \in \mathbb{N}$ the number of cells, and the update (1.22) approaching the solution of (1.20) is replaced for all $1 \le i \le P$ by the implicit kinetic scheme

$$ f_i^{n+1-} = M_i - \sigma \xi (\mathbb{1}_{\xi<0} f_{i+1}^{n+1-} + \mathbb{1}_{\xi>0} f_i^{n+1-} - \mathbb{1}_{\xi<0} f_i^{n+1-} - \mathbb{1}_{\xi>0} f_{i-1}^{n+1-}), \qquad (2.1) $$
## 2 A fully implicit kinetic scheme over a flat bathymetry

In this section we consider the model (1.1) with a flat topography, and propose an implicit kinetic scheme based on (1.19)-(1.20) i.e. an implicit version of (1.21). Since in practice it is not possible to approximate solutions over the whole space $x \in \mathbb{R}$, we shall work in a bounded domain and will need to enforce boundary conditions. First in Section 2.1 we present and study the implicit kinetic scheme; we show that it is unconditionally positive and entropy stable by using the kinetic entropy (1.10). Next in Section 2.2 we detail how to compute explicitly the update at the kinetic level, which is possible thanks to the particular structure of the underlying matrix. Then we consider the associated macroscopic scheme in Section 2.3. Unfortunately, it is not possible to compute the integral of the kinetic update against $(1, \xi)$, instead we propose a compromise consisting to replace the standard Maxwellian by a simpler one. The question of boundary conditions is then treated in Section 2.4, and finally we discuss the practical implementation and the computational costs in Section 2.5.

### 2.1 Kinetic discretization and properties

In the case of a flat topography, the kinetic scheme is a flux vector splitting scheme [12]. We denote $P \in \mathbb{N}$ the number of cells, and the update (1.21) approaching the solution of (1.19) is replaced for all $1 \le i \le P$ by the implicit kinetic scheme

$$ (2.1) \qquad f_i^{n+1-} = M_i - \sigma \xi (1_{\xi<0} f_{i+1}^{n+1-} + 1_{\xi>0} f_i^{n+1-} - 1_{\xi<0} f_i^{n+1-} - 1_{\xi>0} f_{i-1}^{n+1-}), $$

with $\sigma = \Delta t^n / \Delta x$. We rewrite the previous equations under the form

$$ (2.2) \qquad \begin{cases} -\sigma 1_{\xi>0} \xi f_{i-1}^{n+1-} + (1+\sigma|\xi|)f_i^{n+1-} + \sigma 1_{\xi<0} \xi f_{i+1}^{n+1-} = M_i & 2 \le i \le P-1 \\ (1+\sigma|\xi|)f_1^{n+1-} + \sigma 1_{\xi<0} \xi f_2^{n+1-} = M_1 + \sigma 1_{\xi>0} \xi M_0^{n+1} \\ -\sigma 1_{\xi>0} \xi f_{P-1}^{n+1-} + (1+\sigma|\xi|)f_P^{n+1-} = M_P - \sigma 1_{\xi<0} \xi M_{P+1}^{n+1} \end{cases} $$

The quantities $M_0^{n+1} = M(U_0^{n+1}, \xi)$ and $M_{P+1}^{n+1} = M(U_{P+1}^{n+1}, \xi)$ appearing in the last two lines of (2.2) account for the imposed boundary conditions. In a first step, we assume that $M_0^{n+1}$ and $M_{P+1}^{n+1}$ are two known kinetic Maxwellian, their expressions will be discussed in more details in the paragraph devoted to the practical computation of the implicit variables, see paragraph 2.4.

With obvious notations, the system (2.2) consists in finding $f^{n+1} = \{f_i^{n+1-}\}_{i \in \{1, \dots, P\}}$ satisfying

$$ (2.3) \qquad (\mathbf{I} + \sigma \mathbf{L}) f^{n+1} = M + \sigma B^{n+1}, $$

where $\mathbf{I}$ is the identity matrix of length $P$ and $\mathbf{L} \in \mathbb{R}^{P \times P}$ is given by

$$ (2.4) \qquad \mathbf{L} = \begin{pmatrix} |\xi| & \xi 1_{\xi<0} & 0 & \dots & 0 \\ -\xi 1_{\xi>0} & |\xi| & \xi 1_{\xi<0} & \dots & \vdots \\ 0 & \dots & \dots & \dots & 0 \\ \vdots & \dots & -\xi 1_{\xi>0} & |\xi| & \xi 1_{\xi<0} \\ 0 & \dots & 0 & -\xi 1_{\xi>0} & |\xi| \end{pmatrix}. $$

The three vectors $f^{n+1}$, $M$ and $B^{n+1}$ of $\mathbb{R}^P$ are defined by

$$ (2.5) \qquad f^{n+1} = \begin{pmatrix} f_1^{n+1-} \\ \vdots \\ f_i^{n+1-} \\ \vdots \\ f_P^{n+1-} \end{pmatrix}, \quad M = \begin{pmatrix} M_1 \\ \vdots \\ M_i \\ \vdots \\ M_P \end{pmatrix} \quad \text{and} \quad B^{n+1} = \begin{pmatrix} 1_{\xi>0} \xi M_0^{n+1} \\ 0 \\ \vdots \\ 0 \\ -1_{\xi<0} \xi M_{P+1}^{n+1} \end{pmatrix}. $$

The practical computation of the densities vector $f^{n+1}$ will be discussed in paragraph 2.2. Hereafter, we focus on the properties of the numerical scheme (2.2) and the two following results hold.
with $\sigma = \Delta t^n/\Delta x$. We rewrite the previous equations under the form

$$ \begin{cases} -\sigma \mathbb{1}_{\xi>0} \xi f_{i-1}^{n+1-} + (1+\sigma|\xi|) f_i^{n+1-} + \sigma \mathbb{1}_{\xi<0} \xi f_{i+1}^{n+1-} = M_i, & 2 \le i \le P-1 \\ (1+\sigma|\xi|) f_1^{n+1-} + \sigma \mathbb{1}_{\xi<0} \xi f_2^{n+1-} = M_1 + \sigma \mathbb{1}_{\xi>0} \xi M_0^{n+1} \\ -\sigma \mathbb{1}_{\xi>0} \xi f_{P-1}^{n+1-} + (1+\sigma|\xi|) f_P^{n+1-} = M_P - \sigma \mathbb{1}_{\xi<0} \xi M_{P+1}^{n+1}. \end{cases} \quad (2.2) $$

The quantities $M_0^{n+1} = M(U_0^{n+1}, \xi)$ and $M_{P+1}^{n+1} = M(U_{P+1}^{n+1}, \xi)$ appearing in the last two lines of (2.2) account for the imposed boundary conditions. In a first step, we assume that $M_0^{n+1}$ and $M_{P+1}^{n+1}$ are two known kinetic Maxwellian, their expressions will be discussed in more details in the paragraph devoted to the practical computation of the implicit variables, see paragraph 2.4.

With obvious notations, the system (2.2) consists in finding $f^{n+1} = \{f_i^{n+1-}\}_{i \in \{1, \dots, P\}}$ satisfying

$$ (\mathbf{I} + \sigma \mathbf{L}) f^{n+1} = M + \sigma B^{n+1}, \quad (2.3) $$

where $\mathbf{I}$ is the identity matrix of length $P$ and $\mathbf{L} \in \mathbb{R}^{P \times P}$ is given by

$$ \mathbf{L} = \begin{pmatrix} |\xi| & \xi \mathbb{1}_{\xi<0} & 0 & \dots & 0 \\ -\xi \mathbb{1}_{\xi>0} & |\xi| & \xi \mathbb{1}_{\xi<0} & \dots & \vdots \\ 0 & \dots & \dots & \dots & 0 \\ \vdots & \dots & -\xi \mathbb{1}_{\xi>0} & |\xi| & \xi \mathbb{1}_{\xi<0} \\ 0 & \dots & 0 & -\xi \mathbb{1}_{\xi>0} & |\xi| \end{pmatrix}. \quad (2.4) $$

The three vectors $f^{n+1}$, $M$ and $B^{n+1}$ of $\mathbb{R}^P$ are defined by

$$ f^{n+1} = \begin{pmatrix} f_1^{n+1-} \\ \vdots \\ f_i^{n+1-} \\ \vdots \\ f_P^{n+1-} \end{pmatrix}, \quad M = \begin{pmatrix} M_1 \\ \vdots \\ M_i \\ \vdots \\ M_P \end{pmatrix} \quad \text{and} \quad B^{n+1} = \begin{pmatrix} \mathbb{1}_{\xi>0} \xi M_0^{n+1}, \\ 0 \\ \vdots \\ 0 \\ -\mathbb{1}_{\xi<0} \xi M_{P+1}^{n+1} \end{pmatrix}. \quad (2.5) $$

The practical computation of the densities vector $f^{n+1}$ will be discussed in paragraph 2.2. Hereafter, we focus on the properties of the numerical scheme (2.2) and the two following results hold.

**Lemma 2.1.** The matrix $\mathbf{I} + \sigma \mathbf{L}$ defined by through (2.4)

(i) is invertible for any $\sigma$ and $\xi$,

(ii) its inverse $(\mathbf{I} + \sigma \mathbf{L})^{-1}$ has only positive coefficients.

**Proposition 2.2.** The numerical scheme (2.2) satisfies the following properties

(i) the discretization (2.2) is consistent with (1.1),

(ii) the system (2.2) – or equivalently the system (2.3) – admits an unique solution and the solution satisfies

$$ f_i^{n+1-} = f_i^{n+1-}(\xi) \ge 0, \quad \forall 1 \le i \le P, \forall \xi \in \mathbb{R}. $$
**Lemma 2.1** The matrix $I + \sigma L$ defined by through (2.4)

(i) is invertible for any $\sigma$ and $\xi$,

(ii) its inverse $(I + \sigma L)^{-1}$ has only positive coefficients.

**Proposition 2.2** The numerical scheme (2.2) satisfies the following properties

(i) the discretization (2.2) is consistent with (1.1),

(ii) the system (2.2) – or equivalently the system (2.3) – admits a unique solution and the solution satisfies

$$f_i^{n+1-} = f_i^{n+1-}(\xi) \ge 0, \quad \forall 1 \le i \le P, \forall \xi \in \mathbb{R}.$$

To state the main result regarding the entropy stability of the fully implicit kinetic scheme we need to introduce the function $\Psi$ defined by

$$(2.6) \qquad \Psi: \mathbb{R}^2 \ni (a, b) \mapsto \frac{g^2 \pi^2}{6} (b + 2a)(b - a)^2,$$

which is positive on $\mathbb{R}_+^2$.

**Proposition 2.3** The system (2.3) admits a unique solution of positive quantities and defines an implicit kinetic scheme. Moreover, this scheme satisfies the fully discrete entropy equality

$$(2.7) \qquad H_0(f_i^{n+1-}) = H_0(M_i) - \sigma (H_{0,i+1/2}^{n+1-} - H_{0,i-1/2}^{n+1-}) - \Psi(f_i^{n+1-}, M_i) \\ \qquad \qquad + \sigma \xi (1_{\xi<0} \Psi(f_i^{n+1-}, f_{i+1}^{n+1-}) - 1_{\xi>0} \Psi(f_i^{n+1-}, f_{i-1}^{n+1-}))$$

where $H_{0,i+1/2}^{n+1-}$, $H_{0,i-1/2}^{n+1-}$ are given by

$$(2.8) \qquad H_{0,i+1/2}^{n+1-} = \xi 1_{\xi<0} H_0(f_{i+1}^{n+1-}) + \xi 1_{\xi>0} H_0(f_i^{n+1-}),$$

$$(2.9) \qquad H_{0,i-1/2}^{n+1-} = \xi 1_{\xi<0} H_0(f_i^{n+1-}) + \xi 1_{\xi>0} H_0(f_{i-1}^{n+1-}).$$

Since $\Psi$ given in (2.6) is positive on $\mathbb{R}_+^2$, the last three terms of equality (2.7) define a nonpositive dissipative term.

Notice that the results obtained in the two Propositions 2.2 and 2.3 do not require any CFL condition. A consequence of Proposition 2.3 is that, when using the classical Maxwellian (1.4), the macroscopic scheme associated to (2.1) will satisfy a discrete entropy inequality that always dissipates the energy. In fact since the Maxwellian (1.4) minimizes the functional (1.14) we have the following upper bound on the macroscopic entropy $\eta(U_i^{n+1})$

$$\eta(U_i^{n+1}) = \int_{\mathbb{R}} H_0(M(U_i^{n+1}, \xi), \xi) d\xi \le \int_{\mathbb{R}} H_0(f_i^{n+1-}(\xi), \xi) d\xi.$$

We then use equality (2.7) yielding

$$(2.10) \qquad \eta(U_i^{n+1}) \le \eta(U_i^n) - \sigma \left( \int_{\mathbb{R}} H_{0,i+1/2}^{n+1-}(\xi) d\xi - \int_{\mathbb{R}} H_{0,i-1/2}^{n+1-}(\xi) d\xi \right) + \int_{\mathbb{R}} D_i(\xi) d\xi,$$

where $D_i = \mathcal{O}(\Delta t)$ is the negative dissipation term corresponding to the terms of (2.7) involving $\Psi$, making the inequality (2.10) consistent with the entropy inequality (1.2).
To state the main result regarding the entropy stability of the fully implicit kinetic scheme we need to introduce the function $\Psi$ defined by

$$ \Psi : \mathbb{R}^2 \ni (a, b) \mapsto \frac{g^2 \pi^2}{6} (b + 2a)(b - a)^2, \qquad (2.6) $$

which is positive on $\mathbb{R}_+^2$.

**Proposition 2.3.** The system (2.3) admits a unique solution of positive quantities and defines an implicit kinetic scheme. Moreover, this scheme satisfies the fully discrete entropy equality

$$ H_0(f_i^{n+1-}) = H_0(M_i) - \sigma(H_{0,i+1/2}^{n+1-} - H_{0,i-1/2}^{n+1-}) - \Psi(f_i^{n+1-}, M_i) \\ + \sigma\xi(1_{\xi<0}\Psi(f_i^{n+1-}, f_{i+1}^{n+1-}) - 1_{\xi>0}\Psi(f_i^{n+1-}, f_{i-1}^{n+1-})) \qquad (2.7) $$

where $H_{0,i+1/2}^{n+1-}$, $H_{0,i-1/2}^{n+1-}$ are given by

$$ H_{0,i+1/2}^{n+1-} = \xi 1_{\xi<0} H_0(f_{i+1}^{n+1-}) + \xi 1_{\xi>0} H_0(f_i^{n+1-}), \qquad (2.8) $$

$$ H_{0,i-1/2}^{n+1-} = \xi 1_{\xi<0} H_0(f_i^{n+1-}) + \xi 1_{\xi>0} H_0(f_{i-1}^{n+1-}). \qquad (2.9) $$

Since $\Psi$ given in (2.6) is positive on $\mathbb{R}_+^2$, the last three terms of equality (2.7) define a nonpositive dissipative term.

Notice that the results obtained in the two Propositions 2.2 and 2.3 do not require any CFL condition. A consequence of Proposition 2.3 is that, when using the classical Maxwellian (1.4), the macroscopic scheme associated to (2.1) will satisfy a discrete entropy inequality that always dissipates the energy. In fact since the Maxwellian (1.4) minimizes the functional (1.14) we have the following upper bound on the macroscopic entropy $\eta(U_i^{n+1})$

$$ \eta(U_i^{n+1}) = \int_{\mathbb{R}} H_0(M(U_i^{n+1}, \xi), \xi) d\xi \le \int_{\mathbb{R}} H_0(f_i^{n+1-}(\xi), \xi) d\xi. $$

We then use equality (2.7) yielding

$$ \eta(U_i^{n+1}) \le \eta(U_i^n) - \sigma \left( \int_{\mathbb{R}} H_{0,i+1/2}^{n+1-}(\xi) d\xi - \int_{\mathbb{R}} H_{0,i-1/2}^{n+1-}(\xi) d\xi \right) + \int_{\mathbb{R}} D_i(\xi) d\xi, \qquad (2.10) $$

where $D_i = \mathcal{O}(\Delta t)$ is the negative dissipation term corresponding to the terms of (2.7) involving $\Psi$, making the inequality (2.10) consistent with the entropy inequality (1.2).

**Proof of Lemma 2.1.** The matrix $\mathbf{I} + \sigma\mathbf{L}$ writes

$$ \mathbf{I} + \sigma \mathbf{L} = \begin{pmatrix} 1 + \sigma |\xi| & \sigma \xi 1_{\xi<0} & 0 & \dots & 0 \\ -\sigma \xi 1_{\xi>0} & 1 + \sigma |\xi| & \sigma \xi 1_{\xi<0} & \dots & \vdots \\ 0 & \dots & \dots & \dots & 0 \\ \vdots & \dots & -\sigma \xi 1_{\xi>0} & 1 + \sigma |\xi| & \sigma \xi 1_{\xi<0} \\ 0 & \dots & 0 & -\sigma \xi 1_{\xi>0} & 1 + \sigma |\xi| \end{pmatrix}, $$

and it is easy to see that the matrix $\mathbf{I} + \sigma\mathbf{L}$ is strictly diagonally dominant and hence invertible. Moreover the matrix $\mathbf{\Lambda} = \mathbf{I} + \sigma\mathbf{L}$ is such that

$$ \mathbf{\Lambda}_{i,i} > 0, \quad \text{and} \quad \mathbf{\Lambda}_{i,j} \le 0, \text{ when } i \ne j, $$
**Proof of Lemma 2.1.** The matrix $\mathbf{I} + \sigma\mathbf{L}$ writes

$$ \mathbf{I} + \sigma\mathbf{L} = \begin{pmatrix} 1 + \sigma|\xi| & \sigma\xi 1_{\xi<0} & 0 & \dots & 0 \\ -\sigma\xi 1_{\xi>0} & 1 + \sigma|\xi| & \sigma\xi 1_{\xi<0} & \dots & \vdots \\ 0 & \dots & \dots & \dots & 0 \\ \vdots & \dots & -\sigma\xi 1_{\xi>0} & 1 + \sigma|\xi| & \sigma\xi 1_{\xi<0} \\ 0 & \dots & 0 & -\sigma\xi 1_{\xi>0} & 1 + \sigma|\xi| \end{pmatrix}, $$

and it is easy to see that the matrix $\mathbf{I} + \sigma\mathbf{L}$ is strictly diagonally dominant and hence invertible. Moreover the matrix $\mathbf{\Lambda} = \mathbf{I} + \sigma\mathbf{L}$ is such that

$$ \mathbf{\Lambda}_{i,i} > 0, \quad \text{and} \quad \mathbf{\Lambda}_{i,j} \le 0, \text{ when } i \neq j, $$

meaning $\mathbf{I} + \sigma\mathbf{L}$ is a monotone matrix and hence the solution of (2.3) satisfies

$$ f_i^{n+1-} = ((\mathbf{I} + \sigma\mathbf{L})^{-1}(M + \sigma B^{n+1,k}))_i \ge 0, \quad \forall i, $$

proving the result.

Denoting $\mathbf{L}^d$ (resp. $\mathbf{L}^{nd}$) the diagonal (resp. non diagonal) part of $\mathbf{L}$ we can write

$$ \mathbf{I} + \sigma\mathbf{L} = (\mathbf{I} + \sigma\mathbf{L}^d)(\mathbf{I} - (\mathbf{I} + \sigma\mathbf{L}^d)^{-1}(-\sigma\mathbf{L}^{nd})), $$

where all the entries of the matrix

$$ \mathbf{J} = (\mathbf{I} + \sigma\mathbf{L}^d)^{-1}(-\sigma\mathbf{L}^{nd}), $$

are non negative and less than 1. And hence, we can write

$$ (\mathbf{I} + \sigma\mathbf{L})^{-1} = (\mathbf{I} - \mathbf{J})^{-1} (\mathbf{I} + \sigma\mathbf{L}^d)^{-1} = \sum_{k=0}^{\infty} \mathbf{J}^k (\mathbf{I} + \sigma\mathbf{L}^d)^{-1}, $$

proving all the entries of $(\mathbf{I} + \sigma\mathbf{L})^{-1}$ are non negative. $\square$

**Proof of prop. 2.2.** (i) The four terms in parentheses in (2.1) are conservative, and are classically consistent with $\xi\partial_x f$ in (1.19).

(ii) This is a direct consequence of Lemma 2.1. $\square$

The proof of Proposition 2.3 makes use of the following Lemma which will also be useful later.

**Lemma 2.4** The following identity holds for any real pair $(a, b)$ and for any $\xi \in \mathbb{R}$

$$ (2.11) \qquad H_0(b, \xi) = H_0(a, \xi) + \partial_f H_0(a, \xi)(b-a) + \Psi(a, b), $$

with the function $\Psi$ defined in (2.6). Especially, we recover the convexity of $H_0(\cdot, \xi)$ on $\mathbb{R}_+$ thanks to the positivity of $\Psi$ on $\mathbb{R}_+^2$. Equality (2.11) remains satisfied if we replace $H_0$ by $H$.

**Proof of Lemma 2.4.** For any $(a, b)$ in $\mathbb{R}^2$ there holds

$$ \begin{align*} \partial_f H_0(a)(b-a) &= \frac{\xi^2}{2}b + \frac{g^2\pi^2}{2}a^2b - \frac{\xi^2}{2}a - \frac{g^2\pi^2}{2}a^3 \\ &= H_0(b) + \frac{g^2\pi^2}{2}a^2b - \frac{g^2\pi^2}{6}b^3 - H_0(a) - \frac{g^2\pi^2}{2}a^3 + \frac{g^2\pi^2}{6}a^3 \\ &= H_0(b) - H_0(a) - \frac{g^2\pi^2}{6}(b^3 - a^3 - 3a^2(b-a)), \end{align*} $$

and equality (2.11) is recovered using the formula

$$ b^3 - a^3 - 3a^2(b - a) = (b + 2a)(b - a)^2. $$

This result is extended to the kinetic entropy $H$ owing to the relation $H(f, \xi) = H_0(f, \xi) + gzf$. $\square$
meaning $\mathbf{I} + \sigma\mathbf{L}$ is a monotone matrix and hence the solution of (2.3) satisfies

$$f_i^{n+1-} = ((\mathbf{I} + \sigma\mathbf{L})^{-1}(M + \sigma B^{n+1,k}))_i \geq 0, \quad \forall i,$$

proving the result.

Denoting $\mathbf{L}^d$ (resp. $\mathbf{L}^{nd}$) the diagonal (resp. non diagonal) part of $\mathbf{L}$ we can write

$$\mathbf{I} + \sigma\mathbf{L} = (\mathbf{I} + \sigma\mathbf{L}^d)(\mathbf{I} - (\mathbf{I} + \sigma\mathbf{L}^d)^{-1}(-\sigma\mathbf{L}^{nd})),$$

where all the entries of the matrix

$$\mathbf{J} = (\mathbf{I} + \sigma\mathbf{L}^d)^{-1}(-\sigma\mathbf{L}^{nd}),$$

are non negative and less than 1. And hence, we can write

$$(\mathbf{I} + \sigma\mathbf{L})^{-1} = (\mathbf{I} - \mathbf{J})^{-1}(\mathbf{I} + \sigma\mathbf{L}^d)^{-1} = \sum_{k=0}^{\infty} \mathbf{J}^k (\mathbf{I} + \sigma\mathbf{L}^d)^{-1},$$

proving all the entries of $(\mathbf{I} + \sigma\mathbf{L})^{-1}$ are non negative.

*Proof of Proposition 2.2.* (i) The four terms in parentheses in (2.1) are conservative, and are classically consistent with $\xi\partial_x f$ in (1.20).

(ii) This is a direct consequence of Lemma 2.1.

The proof of Proposition 2.3 makes use of the following lemma which will also be useful later.

**Lemma 2.4.** The following identity holds for any real pair $(a, b)$ and for any $\xi \in \mathbb{R}$

$$H_0(b, \xi) = H_0(a, \xi) + \partial_f H_0(a, \xi)(b-a) + \Psi(a, b), \qquad (2.11)$$

with the function $\Psi$ defined in (2.6). Especially, we recover the convexity of $H_0(\cdot, \xi)$ on $\mathbb{R}_+$ thanks to the positivity of $\Psi$ on $\mathbb{R}_+^2$. Equality (2.11) remains satisfied if we replace $H_0$ by $H$.

*Proof of Lemma 2.4.* For any $(a, b)$ in $\mathbb{R}^2$ there holds

$$\begin{aligned} \partial_f H_0(a)(b-a) &= \frac{\xi^2}{2}b + \frac{g^2\pi^2}{2}a^2b - \frac{\xi^2}{2}a - \frac{g^2\pi^2}{2}a^3 \\ &= H_0(b) + \frac{g^2\pi^2}{2}a^2b - \frac{g^2\pi^2}{6}b^3 - H_0(a) - \frac{g^2\pi^2}{2}a^3 + \frac{g^2\pi^2}{6}a^3 \\ &= H_0(b) - H_0(a) - \frac{g^2\pi^2}{6}(b^3 - a^3 - 3a^2(b-a)), \end{aligned}$$

and equality (2.11) is recovered using the formula

$$b^3 - a^3 - 3a^2(b - a) = (b + 2a)(b - a)^2.$$

This result is extended to the kinetic entropy $H$ owing to the relation $H(f, \xi) = H_0(f, \xi) + gzf$.

*Proof of Proposition 2.3.* The proof follows similar lines as what was done in the case of the fully explicit version of the kinetic scheme in [6]. Instead of multiplying equation (2.1) by $\partial_f H_0(f_i^n)$, we multiply it by $\partial_f H_0(f_i^{n+1-})$, which leads to
**Proof of prop. 2.3.** The proof follows similar lines as what was done in the case of the fully explicit version of the kinetic scheme in [4]. Instead of multiplying Equation (2.1) by $\partial_f H_0(f_i^n)$, we multiply it by $\partial_f H_0(f_i^{n+1-})$, which leads to

$$ (2.12) \qquad \partial_f H_0(f_i^{n+1-})(f_i^{n+1-} - M_i) = -\sigma\xi 1_{\xi<0} \partial_f H_0(f_i^{n+1-})(f_{i+1}^{n+1-} - f_i^{n+1-}) \\ \qquad \qquad + \sigma\xi 1_{\xi>0} \partial_f H_0(f_i^{n+1-})(f_{i-1}^{n+1-} - f_i^{n+1-}). $$

In (2.12) we have three terms of the form $\partial_f H(a)(b-a)$ with $a = f_i^{n+1-}$ and $b \in \{f_{i-1}^{n+1-}, M_i, f_{i+1}^{n+1-}\}$. Taking advantage of Lemma 2.4 we can write

$$ H_0(f_i^{n+1-}) - H_0(M_i) + \Psi(f_i^{n+1-}, M_i) = \\ - \sigma\xi 1_{\xi<0} (H_0(f_{i+1}^{n+1-}) - H_0(f_i^{n+1-}) - \Psi(f_i^{n+1-}, f_{i+1}^{n+1-})) \\ + \sigma\xi 1_{\xi>0} (H_0(f_{i-1}^{n+1-}) - H_0(f_i^{n+1-}) - \Psi(f_i^{n+1-}, f_{i-1}^{n+1-})), $$

and we conclude by grouping the expressions. $\square$

## 2.2 Practical computation of the implicit kinetic update

When dealing with implicit schemes, one has often to invert an operator and the key point of the numerical scheme (2.3) is the computation of the inverse of the matrix $I + \sigma L$, where $L$ has been defined in (2.4). In our case, it will be possible to compute analytically this inverse thanks to the triangular structure of the matrix, which is due to the upwinding of the fluxes in (2.1). More precisely we decompose $(I + \sigma L)^{-1}$ as the contributions of the left- and right-going information, which gives

$$ (2.13) \qquad (I + \sigma L)^{-1} = (I + \sigma L^{+})^{-1} 1_{\xi<0} + (I + \sigma L^{-})^{-1} 1_{\xi>0}, $$

with the upwinding matrices $L^+$ and $L^-$ corresponding to

$$ (2.14) \qquad L^+ = \begin{pmatrix} -\xi & \xi & 0 & \dots & 0 \\ 0 & -\xi & \ddots & \vdots & \vdots \\ \vdots & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \xi \\ 0 & \dots & 0 & -\xi & \end{pmatrix}, \qquad L^- = \begin{pmatrix} \xi & 0 & \dots & 0 \\ -\xi & \ddots & \ddots & \vdots \\ 0 & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \dots & 0 & -\xi \end{pmatrix}. $$

Introducing $J^+$ and $J^-$ the matrices of $\mathbb{R}^{P \times P}$ defined as

$$ (J^+)^{i,j} = \begin{cases} 1 & \text{if } i = j - 1 \\ 0 & \text{otherwise} \end{cases}, \qquad (J^-)^{i,j} = \begin{cases} 1 & \text{if } i = j + 1 \\ 0 & \text{otherwise} \end{cases}, $$

we can write

$$ \left\{ \begin{aligned} (I + \sigma L^{-})^{-1} &= ((1 + \sigma\xi)I - \sigma\xi J^{-})^{-1} = \frac{1}{1 + \sigma\xi} \left(I - \frac{\sigma\xi}{1 + \sigma\xi} J^{-}\right)^{-1} \\ (I + \sigma L^{+})^{-1} &= ((1 - \sigma\xi)I + \sigma\xi J^{+})^{-1} = \frac{1}{1 - \sigma\xi} \left(I + \frac{\sigma\xi}{1 - \sigma\xi} J^{+}\right)^{-1}. \end{aligned} \right. $$

The above inverses can be computed through geometric sums since $J_P^+$ and $J_P^-$ have a spectral radius equal to zero. More specifically these two matrices are nilpotent, which implies that the geometric sums in question contain a finite number of nonzero terms and are given below

$$ (I + \sigma L^{-})^{-1} = \sum_{k=0}^{P} \frac{(\sigma\xi)^k}{(1 + \sigma\xi)^{k+1}} (J^{-})^k, \qquad (I + \sigma L^{+})^{-1} = \sum_{k=0}^{P} \frac{(-\sigma\xi)^k}{(1 - \sigma\xi)^{k+1}} (J^{+})^k. $$
$$ \partial_f H_0(f_i^{n+1-}) (f_i^{n+1-} - M_i) = -\sigma \xi \mathbb{1}_{\xi<0} \partial_f H_0(f_i^{n+1-}) (f_{i+1}^{n+1-} - f_i^{n+1-}) + \sigma \xi \mathbb{1}_{\xi>0} \partial_f H_0(f_i^{n+1-}) (f_{i-1}^{n+1-} - f_i^{n+1-}). \quad (2.12) $$

In (2.12) we have three terms of the form $\partial_f H(a)(b-a)$ with $a = f_i^{n+1-}$ and $b \in \{f_{i-1}^{n+1-}, M_i, f_{i+1}^{n+1-}\}$. Taking advantage of Lemma 2.4 we can write

$$ H_0(f_i^{n+1-}) - H_0(M_i) + \Psi(f_i^{n+1-}, M_i) = -\sigma \xi \mathbb{1}_{\xi<0} (H_0(f_{i+1}^{n+1-}) - H_0(f_i^{n+1-}) - \Psi(f_i^{n+1-}, f_{i+1}^{n+1-})) \\ + \sigma \xi \mathbb{1}_{\xi>0} (H_0(f_{i-1}^{n+1-}) - H_0(f_i^{n+1-}) - \Psi(f_i^{n+1-}, f_{i-1}^{n+1-})), $$

and we conclude by grouping the expressions. $\square$

## 2.2. Practical computation of the implicit kinetic update

When dealing with implicit schemes, one has often to invert an operator and the key point of the numerical scheme (2.3) is the computation of the inverse of the matrix $\mathbf{I} + \sigma \mathbf{L}$, where $\mathbf{L}$ has been defined in (2.4). In our case, it will be possible to compute analytically this inverse thanks to the triangular structure of the matrix, which is due to the upwinding of the fluxes in (2.1). More precisely we decompose $(\mathbf{I} + \sigma \mathbf{L})^{-1}$ as the contributions of the left- and right-going information, which gives

$$ (\mathbf{I} + \sigma \mathbf{L})^{-1} = (\mathbf{I} + \sigma \mathbf{L}^{+})^{-1} \mathbb{1}_{\xi<0} + (\mathbf{I} + \sigma \mathbf{L}^{-})^{-1} \mathbb{1}_{\xi>0}, \quad (2.13) $$

with the upwinding matrices $\mathbf{L}^+$ and $\mathbf{L}^-$ corresponding to

$$ \mathbf{L}^+ = \begin{pmatrix} -\xi & \xi & 0 & \dots & 0 \\ 0 & -\xi & \ddots & & \vdots \\ \vdots & \ddots & \ddots & \ddots & 0 \\ \vdots & & \ddots & \ddots & \xi \\ 0 & \dots & 0 & -\xi \end{pmatrix}, \quad \mathbf{L}^- = \begin{pmatrix} \xi & 0 & \dots & 0 \\ -\xi & \ddots & \ddots & \vdots \\ 0 & \ddots & \ddots & \vdots \\ \vdots & & \ddots & 0 \\ 0 & \dots & 0 & -\xi \end{pmatrix}. \quad (2.14) $$

Introducing $\mathbf{J}^+$ and $\mathbf{J}^-$ the matrices of $\mathbb{R}^{P \times P}$ defined as

$$ (\mathbf{J}^+)^{i,j} = \begin{cases} 1 & \text{if } i = j - 1 \\ 0 & \text{otherwise,} \end{cases} \quad (\mathbf{J}^-)^{i,j} = \begin{cases} 1 & \text{if } i = j + 1 \\ 0 & \text{otherwise,} \end{cases} $$

we can write

$$ \left\{ \begin{array}{l} (\mathbf{I} + \sigma \mathbf{L}^{-})^{-1} = ((1 + \sigma \xi)\mathbf{I} - \sigma \xi \mathbf{J}^{-})^{-1} = \frac{1}{1 + \sigma \xi} \left( \mathbf{I} - \frac{\sigma \xi}{1 + \sigma \xi} \mathbf{J}^{-} \right)^{-1} \\ (\mathbf{I} + \sigma \mathbf{L}^{+})^{-1} = ((1 - \sigma \xi)\mathbf{I} + \sigma \xi \mathbf{J}^{+})^{-1} = \frac{1}{1 - \sigma \xi} \left( \mathbf{I} + \frac{\sigma \xi}{1 - \sigma \xi} \mathbf{J}^{+} \right)^{-1}. \end{array} \right. $$

The above inverses can be computed through geometric sums since $\mathbf{J}_P^+$ and $\mathbf{J}_P^-$ have a spectral radius equal to zero. More specifically these two matrices are nilpotent, which implies that the geometric sums in question contain a finite number of nonzero terms and are given below

$$ (\mathbf{I} + \sigma \mathbf{L}^{-})^{-1} = \sum_{k=0}^{P} \frac{(\sigma \xi)^k}{(1 + \sigma \xi)^{k+1}} (\mathbf{J}^{-})^k, \quad (\mathbf{I} + \sigma \mathbf{L}^{+})^{-1} = \sum_{k=0}^{P} \frac{(-\sigma \xi)^k}{(1 - \sigma \xi)^{k+1}} (\mathbf{J}^{+})^k. $$
To conclude we give the analytic expression of the inverse:

$$ (2.15) \qquad (\mathbf{I} + \sigma \mathbf{L}^-)^{-1}_{i,j} = \begin{cases} \frac{(\sigma \xi)^{i-j}}{(1 + \sigma \xi)^{i-j+1}} & \text{if } i \ge j, \\ 0 & \text{else} \end{cases}, $$

$$ (2.16) \qquad (\mathbf{I} + \sigma \mathbf{L}^+)^{-1}_{i,j} = \begin{cases} \frac{(-\sigma \xi)^{j-i}}{(1 - \sigma \xi)^{j-i+1}} & \text{if } i \le j, \\ 0 & \text{else} \end{cases}. $$

Especially we recover the properties enumerated in Lemma 2.1, since we see that all the coefficients of the inverse (2.13) are comprised between zero and one respectively when $\xi \ge 0$ and $\xi \le 0$.

## 2.3 Macroscopic implicit scheme

We now focus on obtaining an explicit writing of the macroscopic update (1.19) associated to (2.1). Using the expressions (2.15)-(2.16) of the matrix inverse, we need to compute the integral of

$$ \mathbb{1}_{\pm\xi>0} \frac{(\pm\sigma\xi)^k}{(1 \pm \sigma\xi)^{k+1}} M(U, \xi) $$

against 1, $\xi$ and $\xi^2$ for all $0 \le k \le P-1$. This seems hardly possible with the classical Maxwellian (1.4) and we use the result presented in the following remark.

**Remark 2.5** Let us recall that the half-disk Maxwellian $M(U, \xi)$ defined by (1.4) has some optimal properties presented in Lemma 1.2, which allow to obtain the discrete entropy inequality (2.10) at the macroscopic scale. Other choices of Maxwellians are possible but the discrete entropy inequality (2.7) is not granted to hold anymore. A general possibility is to choose $M(U, \xi)$ of the form

$$ M(U, \xi) = \frac{h}{c} \chi\left(\frac{\xi - u}{c}\right), $$

with

$$ c = \sqrt{\frac{gh}{2}}, $$

where $\chi$ is a non-negative, compactly supported even function satisfying

$$ \int_{\mathbb{R}} \chi(z) dz = \int_{\mathbb{R}} z^2 \chi(z) dz = 1. $$

There are a lot of possible choices for $\chi$ e.g.

$$ (2.17) \qquad \chi_1(z) = \frac{1}{2\sqrt{3}} \mathbb{1}_{|z|\le\sqrt{3}}, \quad \text{or} \quad \chi_2(z) = \frac{3}{20\sqrt{5}} z^2 + \frac{3}{4\sqrt{5}} \mathbb{1}_{|z|\le\sqrt{5}}. $$

Notice that the definition (1.4) of the half-disk Maxwellian corresponds to the choice

$$ \chi_3(z) = \frac{1}{\pi} \sqrt{1 - \frac{z^2}{4}} \mathbb{1}_{|z| \le 2}. $$

In order to compute explicitly the macroscopic scheme i.e. without approximate quadrature formula, we now use in this section the following expression for the Maxwellian

$$ (2.18) \qquad M(U, \xi) = \frac{h}{c} \chi_1\left(\frac{\xi - u}{c}\right) = \frac{h}{2\sqrt{3c}} \mathbb{1}_{|\xi - u| \le \sqrt{3c}} $$

and referred to as the index Maxwellian. This is the simplest choice we can make, and it will enable us to obtain analytic expressions for the aforementioned integrals. Furthermore note that it satisfies all the moment relations (1.5), which is important for consistency with the macroscopic Saint-Venant system.
To conclude we give the analytic expression of the inverse:

$$ (\mathbf{I} + \sigma \mathbf{L}^{-})^{-1}_{i,j} = \begin{cases} \frac{(\sigma \xi)^{i-j}}{(1 + \sigma \xi)^{i-j+1}} & \text{if } i \geq j \\ 0 & \text{else,} \end{cases} \qquad (2.15) $$

$$ (\mathbf{I} + \sigma \mathbf{L}^{+})^{-1}_{i,j} = \begin{cases} \frac{(-\sigma \xi)^{j-i}}{(1 - \sigma \xi)^{j-i+1}} & \text{if } i \leq j \\ 0 & \text{else.} \end{cases} \qquad (2.16) $$

Especially we recover the properties enumerated in Lemma 2.1, since we see that all the coefficients of the inverse (2.13) are comprised between zero and one respectively when $\xi \ge 0$ and $\xi \le 0$.

## 2.3. Macroscopic implicit scheme

We now focus on obtaining an explicit writing of the macroscopic update (1.20) associated to (2.1). Using the expressions (2.15) and (2.16) of the matrix inverse, we need to compute the integral of

$$ \mathbb{1}_{\pm\xi>0} \frac{(\pm\sigma\xi)^k}{(1 \pm \sigma\xi)^{k+1}} M(U, \xi) $$

against 1, $\xi$ and $\xi^2$ for all $0 \le k \le P - 1$. This seems hardly possible with the classical Maxwellian (1.4) and we use the result presented in the following remark.

**Remark 2.5.** Let us recall that the half-disk Maxwellian $M(U, \xi)$ defined by (1.4) has some optimal properties presented in Lemma 1.2, which allow to obtain the discrete entropy inequality (2.10) at the macroscopic scale. Other choices of Maxwellians are possible but the discrete entropy inequality (2.7) is not granted to hold anymore. A general possibility is to choose $M(U, \xi)$ of the form

$$ M(U, \xi) = \frac{h}{c} \chi\left(\frac{\xi - u}{c}\right), $$

with

$$ c = \sqrt{\frac{gh}{2}}, $$

where $\chi$ is a non-negative, compactly supported even function satisfying

$$ \int_{\mathbb{R}} \chi(z) dz = \int_{\mathbb{R}} z^2 \chi(z) dz = 1. $$

There are a lot of possible choices for $\chi$ e.g.

$$ \chi_1(z) = \frac{1}{2\sqrt{3}} \mathbb{1}_{|z| \le \sqrt{3}}, \quad \text{or} \quad \chi_2(z) = \frac{3}{20\sqrt{5}} z^2 + \frac{3}{4\sqrt{5}} \mathbb{1}_{|z| \le \sqrt{5}}. \qquad (2.17) $$

Notice that the definition (1.4) of the half-disk Maxwellian corresponds to the choice

$$ \chi_3(z) = \frac{1}{\pi} \sqrt{1 - \frac{z^2}{4}} \mathbb{1}_{|z| \le 2}. $$

In order to compute explicitly the macroscopic scheme i.e. without approximate quadrature formula, we now use in this section the following expression for the Maxwellian

$$ M(U, \xi) = \frac{h}{c} \chi_1 \left( \frac{\xi - u}{c} \right) = \frac{h}{2\sqrt{3c}} \mathbb{1}_{|\xi - u| \le \sqrt{3c}} \qquad (2.18) $$
In order to compute the macroscopic version of the scheme (2.2), we proceed as follows. The strategy consists to dissociate the contribution of the information coming from the interior of the computational domain, and the one coming from the boundaries

$$ (2.19) \qquad \left\{ \begin{array}{l} U^{\text{int}} = \int_{\mathbb{R}} \left(\frac{1}{\xi}\right) (\mathbf{I} + \sigma \mathbf{L})^{-1} M \, d\xi \\ U^{\text{ext}} = \int_{\mathbb{R}} \left(\frac{1}{\xi}\right) (\mathbf{I} + \sigma \mathbf{L})^{-1} \sigma B^{n+1} \, d\xi \end{array} \right. . $$

The final update is then set as $U^{n+1} = U^{\text{int}} + U^{\text{ext}}$, which coincides with definition (1.20). We postpone the details about the computation of $B^{n+1}$ to the next section, and assume that it is known for now. First for $U^{\text{int}}$, we have

$$ U^{\text{int}} = \int_{\mathbb{R}} \mathbb{1}_{\xi \le 0} \left(\frac{1}{\xi}\right) (\mathbf{I} + \sigma \mathbf{L}^{+})^{-1} M d\xi + \int_{\mathbb{R}} \mathbb{1}_{\xi \ge 0} \left(\frac{1}{\xi}\right) (\mathbf{I} + \sigma \mathbf{L}^{-})^{-1} M d\xi. $$

Plugging the analytic expressions (2.15) and (2.16) in the above integrals, we can express the i-th component of $U^{\text{int}}$ as

$$ (2.20) \qquad \begin{aligned} U_i^{\text{int}} &= \int_{\xi<0} \sum_{j=1}^{P} \left(\frac{1}{\xi}\right) (\mathbf{I} + \sigma \mathbf{L}^{+})^{-1}_{i,j} M_j d\xi + \int_{\xi>0} \sum_{j=1}^{P} \left(\frac{1}{\xi}\right) (\mathbf{I} + \sigma \mathbf{L}^{-})^{-1}_{i,j} M_j d\xi \\ &= \int_{\xi<0} \sum_{j=i}^{P} \left(\frac{1}{\xi}\right) \frac{(-\sigma\xi)^{j-i}}{(1-\sigma\xi)^{j-i+1}} M_j d\xi + \int_{\xi>0} \sum_{j=1}^{i} \left(\frac{1}{\xi}\right) \frac{(\sigma\xi)^{i-j}}{(1+\sigma\xi)^{i-j+1}} M_j d\xi. \end{aligned} $$

A detailed expression of the quantities appearing in relation (2.20) is given in Appendix A. Similarly for the exterior contribution we have

$$ U^{\text{ext}} = \sigma \int_{\xi<0} \left(\frac{1}{\xi}\right) (\mathbf{I} + \sigma \mathbf{L}^{+})^{-1} B^{n+1} d\xi + \sigma \int_{\xi>0} \left(\frac{1}{\xi}\right) (\mathbf{I} + \sigma \mathbf{L}^{-})^{-1} B^{n+1} d\xi. $$

Using definition (2.5) and equalities (2.15)–(2.16), the i-th component of $U^{\text{ext}}$ is

$$ (2.21) \qquad U_i^{\text{ext}} = \int_{\xi<0} \left(\frac{1}{\xi}\right) \frac{(-\sigma\xi)^{P-i+1}}{(1-\sigma\xi)^{P-i+1}} M_{P+1}^{n+1} d\xi + \int_{\xi>0} \left(\frac{1}{\xi}\right) \frac{(\sigma\xi)^i}{(1+\sigma\xi)^i} M_0^{n+1} d\xi. $$

The expressions necessary to compute (2.20) and (2.21) are detailed in Appendix A and Appendix B.

## 2.4 Boundary conditions

In this paragraph we discuss how to enforce the boundary conditions associated with the contribution $U^{\text{ext}}$ in (2.19). To achieve this, it will be necessary to introduce two ghost cells numbered 0 and $P+1$ and neighbours of the cells $C_1$ and $C_P$ and we have to define $U_0^{n+1}$ and $U_{P+1}^{n+1}$. The problem we are facing is that these ghost quantities depend on the neighboring values in cells $C_1$ and $C_P$ at time $t^{n+1}$, and which are themselves unknown. Hence we have an implicit problem where the relation between the ghost and border terms can be nonlinear depending on the type of boundary conditions. In practice we will avoid this issue by substituting $B^{n+1}$ with $B^n$ in the definition (2.19) of $U^{\text{ext}}$. Doing so can be interpreted as a first order approximation in time since we have

$$ U_0^{n+1} = U_0^n + O(\Delta t), \quad U_{P+1}^{n+1} = U_{P+1}^n + O(\Delta t). $$

Note that explicating the ghost values does not prevent the scheme from satisfying a discrete entropy as the one in Proposition 2.3. In fact, equality (2.7) still holds in the interior cells $2 \le i \le P-1$, and in border cells $i=1, N$ one obtains

$$ (2.22) \qquad H_0(f_1^{n+1-}) = H_0(M_1) - \sigma (H_{0,3/2}^{n+1-} - H_{0,1/2}^{n+1-}) - \Psi(f_1^{n+1-}, M_1) $$
and referred to as the index Maxwellian. This is the simplest choice we can make, and it will enable us to obtain analytic expressions for the aforementioned integrals. Furthermore note that it satisfies all the moment relations (1.5), which is important for consistency with the macroscopic Saint-Venant system.

In order to compute the macroscopic version of the scheme (2.2), we proceed as follows. The strategy consists to dissociate the contribution of the information coming from the interior of the computational domain, and the one coming from the boundaries

$$ \left\{ \begin{aligned} U^{\text{int}} &= \int_{\mathbb{R}} \binom{1}{\xi} (\mathbf{I} + \sigma \mathbf{L})^{-1} M \, d\xi \\ U^{\text{ext}} &= \int_{\mathbb{R}} \binom{1}{\xi} (\mathbf{I} + \sigma \mathbf{L})^{-1} \sigma B^{n+1} \, d\xi. \end{aligned} \right. \qquad (2.19) $$

The final update is then set as $U^{n+1} = U^{\text{int}} + U^{\text{ext}}$, which coincides with definition (1.21). We postpone the details about the computation of $B^{n+1}$ to the next section, and assume that it is known for now. First for $U^{\text{int}}$, we have

$$ U^{\text{int}} = \int_{\mathbb{R}} \mathbb{1}_{\xi \le 0} \binom{1}{\xi} (\mathbf{I} + \sigma \mathbf{L}^{+})^{-1} M \, d\xi + \int_{\mathbb{R}} \mathbb{1}_{\xi \ge 0} \binom{1}{\xi} (\mathbf{I} + \sigma \mathbf{L}^{-})^{-1} M \, d\xi. $$

Plugging the analytic expressions (2.15) and (2.16) in the above integrals, we can express the i-th component of $U^{\text{int}}$ as

$$ \begin{aligned} U_i^{\text{int}} &= \int_{\xi<0} \sum_{j=1}^{P} \binom{1}{\xi} (\mathbf{I} + \sigma \mathbf{L}^{+})^{-1}_{i,j} M_j \, d\xi + \int_{\xi>0} \sum_{j=1}^{P} \binom{1}{\xi} (\mathbf{I} + \sigma \mathbf{L}^{-})^{-1}_{i,j} M_j \, d\xi \\ &= \int_{\xi<0} \sum_{j=i}^{P} \binom{1}{\xi} \frac{(-\sigma\xi)^{j-i}}{(1-\sigma\xi)^{j-i+1}} M_j \, d\xi + \int_{\xi>0} \sum_{j=1}^{i} \binom{1}{\xi} \frac{(\sigma\xi)^{i-j}}{(1+\sigma\xi)^{i-j+1}} M_j \, d\xi. \end{aligned} \qquad (2.20) $$

A detailed expression of the quantities appearing in relation (2.20) is given in Appendix A. Similarly for the exterior contribution we have

$$ U^{\text{ext}} = \sigma \int_{\xi<0} \binom{1}{\xi} (\mathbf{I} + \sigma \mathbf{L}^{+})^{-1} B^{n+1} \, d\xi + \sigma \int_{\xi>0} \binom{1}{\xi} (\mathbf{I} + \sigma \mathbf{L}^{-})^{-1} B^{n+1} \, d\xi. $$

Using definition (2.5) and equalities (2.15)-(2.16), the i-th component of $U^{\text{ext}}$ is

$$ U_i^{\text{ext}} = \int_{\xi<0} \binom{1}{\xi} \frac{(-\sigma\xi)^{P-i+1}}{(1-\sigma\xi)^{P-i+1}} M_{P+1}^{n+1} \, d\xi + \int_{\xi>0} \binom{1}{\xi} \frac{(\sigma\xi)^i}{(1+\sigma\xi)^i} M_0^{n+1} \, d\xi. \qquad (2.21) $$

The expressions necessary to compute (2.20) and (2.21) are detailed in Appendices A and B.

## 2.4. Boundary conditions

In this paragraph we discuss how to enforce the boundary conditions associated with the contribution $U^{\text{ext}}$ in (2.19). To achieve this, it will be necessary to introduce two ghost cells numbered 0 and $P+1$ and neighbors of the cells $C_1$ and $C_P$ and we have to define $U_0^{n+1}$ and $U_{P+1}^{n+1}$. The problem we are facing is that these ghost quantities depend on the neighboring values in cells $C_1$ and $C_P$ at time $t^{n+1}$, and which are themselves unknown. Hence we have an implicit problem where the relation between the ghost and border terms can be nonlinear depending on the type of boundary conditions. In practice we will avoid this issue by substituting $B^{n+1}$ with $B^n$ in the definition (2.19) of $U^{\text{ext}}$. Doing so can be interpreted as a first order approximation in time since we have

$$ U_0^{n+1} = U_0^n + O(\Delta t), \quad U_{P+1}^{n+1} = U_{P+1}^n + O(\Delta t). $$
$$ (2.23) \qquad \begin{aligned} & + \sigma \xi (1_{\xi<0} \Psi(f_1^{n+1-}, f_2^{n+1-}) - 1_{\xi>0} \Psi(f_1^{n+1-}, M_0^n)), \\ & H_0(f_P^{n+1-}) = H_0(M_P) - \sigma (H_{0,P+1/2}^{n+1-} - H_{0,P-1/2}^{n+1-}) - \Psi(f_P^{n+1-}, M_P) \\ & + \sigma \xi (1_{\xi<0} \Psi(f_P^{n+1-}, f_{P-1}^{n+1-}) - 1_{\xi>0} \Psi(f_P^{n+1-}, M_{P+1}^n)), \end{aligned} $$

where the border entropy fluxes $H_{0,1/2}^{n+1-}$, $H_{0,P+1/2}^{n+1-}$ are given by

$$ H_{0,1/2}^{n+1-} = \xi 1_{\xi<0} H_0(f_1^{n+1-}) + \xi 1_{\xi>0} H_0(M_0^n), \quad H_{0,P+1/2}^{n+1-} = \xi 1_{\xi<0} H_0(M_P^n) + \xi 1_{\xi>0} H_0(f_P^{n+1-}). $$

The benefit of expliciting the ghost Maxwellian $M_0, M_{P+1}$ is that we can more easily determine the macroscopic ghost quantities $U_0^n, U_{P+1}^n$ at time $t^n$ based on $U_1^n, U_P^n$ following the procedure described hereafter and similar to that of Bristeau and Coussin in [13]. We will focus on fluvial flows where the material velocity of particles $|u|$ is smaller than the celerity of surface gravity waves $\sqrt{gh}$; in particular low Froude flows enter this regime. Since in this case the eigenvalues $u - \sqrt{gh}$ and $u + \sqrt{gh}$ have opposite sign, at each boundary we have exactly one wave entering the domain and one wave leaving it. Hence we dispose of a single degree of freedom to set the ghost values, which generally consists in enforcing either a given water height or a discharge. The ghost state is then fully determined by asking the outgoing Riemann invariant to remain constant through the interface.

In the following we detail the handling of the left boundary condition when enforcing either the water height or the discharge; we put aside the case of the right boundary condition as it is symmetric.

**Given water height.** First we consider the case where the water height is enforced at the left boundary of the domain, and we denote by $h_{g,l}$ the value attributed to the ghost cell. Together with the condition on the outgoing Riemann invariant, we get the following nonlinear system which is solved analytically

$$ \left\{ \begin{array}{l} h_0^n = h_{g,l} \\ u_0^n - 2\sqrt{gh_0^n} = u_1^n - 2\sqrt{gh_1^n} \end{array} \right. \quad \Rightarrow \quad U_0^n = h_{g,l} \left( u_1^n - 2(\sqrt{gh_1^n} - \sqrt{gh_{g,l}}) \right). $$

**Given flux.** Another possibility is to enforce the discharge at the boundary, and we denote by $q_{g,l}$ the left ghost value. This time, the constraint on the Riemann invariant will enable to determine the ghost water height. Indeed we have the system

$$ (2.24) \qquad \left\{ \begin{array}{l} q_0^n = q_{g,l} \\ u_0^n - 2\sqrt{gh_0^n} = u_1^n - 2\sqrt{gh_1^n} \end{array} \right. , $$

and the second equality involving the outgoing Riemann invariant requires to find the real roots of the third order polynomial in $X = \sqrt{h_0^n}$ given by

$$ (2.25) \qquad \mathcal{P}(X) := X^3 + \frac{u_1^n - 2\sqrt{gh_1^n}}{2\sqrt{g}} X^2 - \frac{q_{g,l}}{2\sqrt{g}} = 0. $$

Its derivative $\mathcal{P}'$ admits two real distinct roots

$$ X_0 = 0 \quad \text{and} \quad X_1 = -\frac{u_1^n - 2\sqrt{gh_1^n}}{3\sqrt{g}}, \quad \text{with} \quad \mathcal{P}(X_1) = -\frac{1}{2}X_1^3 - \frac{q_{g,l}}{2\sqrt{g}}, $$

and we recall that we make the assumption that $u_1^n - \sqrt{gh_1^n} \le 0$ so that $X_1 \ge 0$. If $q_{g,l} \ge -\sqrt{g}X_1^3$, one has $\mathcal{P}(X_1) \le 0$; together with the fact that $\mathcal{P}(X) \to +\infty$ as $X$ goes to infinity, there exists one real root larger than $X_1 \ge 0$. If additionally $q_{g,l} > 0$, this root is the only real one since $\mathcal{P}(X_0) < 0$; if $q_{g,l} \le 0$ then $\mathcal{P}(X_0) \ge 0$ and there exists a second positive root between $X_0$ and $X_1$. In this case we take the convention to retain the smaller of the two. On the other hand, when $q_{g,l} < -\sqrt{g}X_1^3$ we have that $\mathcal{P}(X_0), \mathcal{P}(X_1)$ are both positive, and there doesn't exist a positive root. In this case one can choose to enforce the water height instead.
Note that explicating the ghost values does not prevent the scheme from satisfying a discrete entropy as the one in Proposition 2.3. In fact, equality (2.7) still holds in the interior cells $2 \le i \le P - 1$, and in border cells $i = 1, N$ one obtains

$$ H_0(f_1^{n+1-}) = H_0(M_1) - \sigma(H_{0,3/2}^{n+1-} - H_{0,1/2}^{n+1-}) - \Psi(f_1^{n+1-}, M_1) \quad (2.22) \\ + \sigma\xi(1_{\xi<0}\Psi(f_1^{n+1-}, f_2^{n+1-}) - 1_{\xi>0}\Psi(f_1^{n+1-}, M_0^n)), $$

$$ H_0(f_P^{n+1-}) = H_0(M_P) - \sigma(H_{0,P+1/2}^{n+1-} - H_{0,P-1/2}^{n+1-}) - \Psi(f_P^{n+1-}, M_P) \quad (2.23) \\ + \sigma\xi(1_{\xi<0}\Psi(f_P^{n+1-}, f_{P-1}^{n+1-}) - 1_{\xi>0}\Psi(f_P^{n+1-}, M_{P+1}^n)), $$

where the border entropy fluxes $H_{0,1/2}^{n+1-}$, $H_{0,P+1/2}^{n+1-}$ are given by

$$ H_{0,1/2}^{n+1-} = \xi 1_{\xi<0} H_0(f_1^{n+1-}) + \xi 1_{\xi>0} H_0(M_0^n), \quad H_{0,P+1/2}^{n+1-} = \xi 1_{\xi<0} H_0(M_{P+1}^n) + \xi 1_{\xi>0} H_0(f_P^{n+1-}). $$

The benefit of explicating the ghost Maxwellian $M_0, M_{P+1}$ is that we can more easily determine the macroscopic ghost quantities $U_0^n, U_{P+1}^n$ at time $t^n$ based on $U_1^n, U_P^n$ following the procedure described hereafter and similar to that of Bristeau and Coussin in [13]. We will focus on fluvial flows where the material velocity of particles $|u|$ is smaller than the celerity of surface gravity waves $\sqrt{gh}$; in particular low Froude flows enter this regime. Since in this case the eigenvalues $u - \sqrt{gh}$ and $u + \sqrt{gh}$ have opposite sign, at each boundary we have exactly one wave entering the domain and one wave leaving it. Hence we dispose of a single degree of freedom to set the ghost values, which generally consists in enforcing either a given water height or a discharge. The ghost state is then fully determined by asking the outgoing Riemann invariant to remain constant through the interface.

In the following we detail the handling of the left boundary condition when enforcing either the water height or the discharge; we put aside the case of the right boundary condition as it is symmetric.

**Given water height.** First we consider the case where the water height is enforced at the left boundary of the domain, and we denote by $h_{g,l}$ the value attributed to the ghost cell. Together with the condition on the outgoing Riemann invariant, we get the following nonlinear system which is solved analytically

$$ \left\{ \begin{array}{l} h_0^n = h_{g,l} \\ u_0^n - 2\sqrt{gh_0^n} = u_1^n - 2\sqrt{gh_1^n} \end{array} \right. \implies U_0^n = h_{g,l} \left( \begin{array}{l} 1 \\ u_1^n - 2(\sqrt{gh_1^n} - \sqrt{gh_{g,l}}) \end{array} \right). $$

**Given flux.** Another possibility is to enforce the discharge at the boundary, and we denote by $q_{g,l}$ the left ghost value. This time, the constraint on the Riemann invariant will enable to determine the ghost water height. Indeed we have the system

$$ \left\{ \begin{array}{l} q_0^n = q_{g,l} \\ u_0^n - 2\sqrt{gh_0^n} = u_1^n - 2\sqrt{gh_1^n}, \end{array} \right. \quad (2.24) $$

and the second equality involving the outgoing Riemann invariant requires to find the real roots of the third order polynomial in $X = \sqrt{h_0^n}$ given by

$$ \mathcal{P}(X) := X^3 + \frac{u_1^n - 2\sqrt{gh_1^n}}{2\sqrt{g}} X^2 - \frac{q_{g,l}}{2\sqrt{g}} = 0. \quad (2.25) $$

Its derivative $\mathcal{P}'$ admits two real distinct roots

$$ X_0 = 0 \quad \text{and} \quad X_1 = -\frac{u_1^n - 2\sqrt{gh_1^n}}{3\sqrt{g}}, \quad \text{with} \quad \mathcal{P}(X_1) = -\frac{1}{2}X_1^3 - \frac{q_{g,l}}{2\sqrt{g}}, $$
Note that when enforcing the flux, our approach differs from that of Bristeau and Coussin in [13], where the ghost value is chosen such that the resulting numerical flux at the interface coincides with the boundary discharge. Instead we do not enforce any value at the interface but directly in the ghost cell, which can be seen as a first order simplification in space. A common practice for channel flows is to enforce the water height at the inlet and the flux at the outlet.

**Remark 2.6** When substituting $B^{n+1}$ with $B^n$ in the implicit kinetic scheme (2.3), the corresponding update can be reformulated as $(\overline{\mathbf{I}} + \sigma \overline{\mathbf{L}})\overline{f}^{n+1} = \overline{M}^n$ with

$$ \overline{\mathbf{I}} + \sigma \overline{\mathbf{L}} = \begin{pmatrix} 1 & & \\ & -\sigma \xi \mathbb{1}_{\xi>0} & \\ 0 & \mathbf{I} + \sigma \mathbf{L} & \vdots \\ & \vdots & \sigma \xi \mathbb{1}_{\xi<0} \\ & & \frac{1}{1} \end{pmatrix}, \quad \overline{f}^{n+1} = \begin{pmatrix} f_0^{n+1} \\ f_1^{n+1} \\ f_{P+1}^{n+1} \end{pmatrix}, \quad \overline{M}^n = \begin{pmatrix} M_0^n \\ M_1^n \\ M_{P+1}^n \end{pmatrix} $$

As a consequence the maximum principle $\|\overline{f}^{n+1}(\xi)\|_{\infty} \le \|\overline{M}^n(\xi)\|_{\infty}$ holds for any $\xi$ in $\mathbb{R}$ during the transport step. In fact we can verify that matrix $(\overline{\mathbf{I}} + \sigma \overline{\mathbf{L}})$ is monotone, and following the argument involved in Lemma 5.1 from [1] we can write

$$ 0 \le \overline{f}^{n+1} = (\overline{\mathbf{I}} + \sigma \overline{\mathbf{L}})^{-1} \overline{M}^n \le (\overline{\mathbf{I}} + \sigma \overline{\mathbf{L}})^{-1} (\|\overline{M}^n\|_{\infty} \mathbf{1}), $$

with $\mathbf{1}$ the vector from $\mathbb{R}^{P+2}$ whose entries are all equal to one. Using equality $(\overline{\mathbf{I}} + \sigma \overline{\mathbf{L}})^{-1} \mathbf{1} = \mathbf{1}$ allows to conclude. Note however that there is no such principle at the macroscopic scale, similarly to the continuous Saint-Venant system.

## 2.5 Implementation and computational costs

It is important to keep a reasonable algorithmic complexity so that the implicit method presented in the previous lines and especially the formulae (2.20),(2.21) can be used in practice. We discuss here how to improve its computational cost by a substantial margin. In Appendix A, we show that the $i$-th component of vectors $h^{\text{int}}$ and $(hu)^{\text{int}}$ have the form

$$ (2.26) \qquad \left\{ \begin{aligned} h_i^{\text{int}} &= \frac{1}{2\sigma\sqrt{3}} \left( \sum_{j=i}^{P} \sqrt{\frac{2h_j^n}{g}} (Ah)_{i,j} + \sum_{j=1}^{i} \sqrt{\frac{2h_j^n}{g}} (Bh)_{i,j} \right) \\ (hu)_i^{\text{int}} &= \frac{1}{2\sigma^2\sqrt{3}} \left( -\sum_{j=i}^{P} \sqrt{\frac{2h_j^n}{g}} (Ahu)_{i,j} + \sum_{j=1}^{i} \sqrt{\frac{2h_j^n}{g}} (Bhu)_{i,j} \right), \end{aligned} \right. $$

where $Ah, Ahu$ are dense upper triangular matrices, and $Bh, Bhu$ are dense lower triangular matrices. Therefore computing $h^{\text{int}}$ and $(hu)^{\text{int}}$ through (2.26) is analog to performing a matrix-vector product which has a quadratic complexity $O(P^2)$, and we cannot hope to do better than that. However the coefficients (A.3)-(A.6) of the above matrices involve a summation, and at a first glance the cost to assemble them is seemingly cubic. This is quite expensive and can render the method pretty much inefficient. However this complexity can be reduced to a quadratic cost by computing the coefficients in the correct order. More specifically we show that all the matrices above can be defined through a recurrence relation allowing to compute each coefficient from a previous one in $O(1)$ operation. In fact, denoting $y = x/(1+x)$ and $z = \ln|1+x|$, the matrix $Ah$ is given by

$$ \left( \begin{array}{c|c|c|c} [z]^{-\min(0,a_1)\sigma}_{-\min(0,b_1)\sigma} & [z-y]^{-\min(0,a_2)\sigma}_{-\min(0,b_2)\sigma} & \cdots & [z - \sum_{l=1}^{P-1} y^l/l]^{-\min(0,a_P)\sigma}_{-\min(0,b_P)\sigma} \\ 0 & [z]^{-\min(0,a_2)\sigma}_{-\min(0,b_2)\sigma} & \cdots & \vdots \\ \vdots & \vdots & \ddots & \vdots \\ \vdots & \vdots & \ddots & [z-y]^{-\min(0,a_P)\sigma}_{-\min(0,b_P)\sigma} \\ 0 & \cdots & \cdots & [z]^{-\min(0,a_P)\sigma}_{-\min(0,b_P)\sigma} \end{array} \right), $$
and we recall that we make the assumption that $u_1^n - \sqrt{gh_1^n} \le 0$ so that $X_1 \ge 0$. If $q_{g,l} \ge -\sqrt{g}X_1^3$, one has $\mathcal{P}(X_1) \le 0$; together with the fact that $\mathcal{P}(X) \to +\infty$ as $X$ goes to infinity, there exists one real root larger than $X_1 \ge 0$. If additionally $q_{g,l} > 0$, this root is the only real one since $\mathcal{P}(X_0) < 0$; if $q_{g,l} \le 0$ then $\mathcal{P}(X_0) \ge 0$ and there exists a second positive root between $X_0$ and $X_1$. In this case we take the convention to retain the smaller of the two. On the other hand, when $q_{g,l} < -\sqrt{g}X_1^3$ we have that $\mathcal{P}(X_0), \mathcal{P}(X_1)$ are both positive, and there doesn't exist a positive root. In this case one can choose to enforce the water height instead.

Note that when enforcing the flux, our approach differs from that of Bristeau and Coussin in [13], where the ghost value is chosen such that the resulting numerical flux at the interface coincides with the boundary discharge. Instead we do not enforce any value at the interface but directly in the ghost cell, which can be seen as a first order simplification in space. A common practice for channel flows is to enforce the water height at the inlet and the flux at the outlet.

**Remark 2.6.** When substituting $B^{n+1}$ with $B^n$ in the implicit kinetic scheme (2.3), the corresponding update can be reformulated as $(\overline{\mathbf{I} + \sigma \mathbf{L}}) \bar{f}^{n+1} = \overline{M}^n$ with

$$ \overline{\mathbf{I} + \sigma \mathbf{L}} = \begin{pmatrix} 1 & & \\ & -\sigma \xi \mathbf{1}_{\xi>0} & \vdots \\ 0 & \mathbf{I} + \sigma \mathbf{L} & 0 \\ \vdots & & \sigma \xi \mathbf{1}_{\xi<0} \\ & & \frac{1}{1} \end{pmatrix}, \quad \bar{f}^{n+1} = \begin{pmatrix} f_0^{n+1} \\ f_1^{n+1} \\ f_{P+1}^{n+1} \end{pmatrix}, \quad \overline{M}^n = \begin{pmatrix} M_0^n \\ M_1^n \\ M_{P+1}^n \end{pmatrix}. $$

As a consequence the maximum principle $\|\bar{f}^{n+1}(\xi)\|_{\infty} \le \|\overline{M}^n(\xi)\|_{\infty}$ holds for any $\xi$ in $\mathbb{R}$ during the transport step. In fact we can verify that matrix $(\overline{\mathbf{I} + \sigma \mathbf{L}})$ is monotone, and following the argument involved in Lemma 5.1 from [1] we can write

$$ 0 \le \bar{f}^{n+1} = (\overline{\mathbf{I} + \sigma \mathbf{L}})^{-1} \overline{M}^n \le (\overline{\mathbf{I} + \sigma \mathbf{L}})^{-1} (\|\overline{M}^n\|_{\infty} \mathbf{1}), $$

with $\mathbf{1}$ the vector from $\mathbb{R}^{P+2}$ whose entries are all equal to one. Using equality $(\overline{\mathbf{I} + \sigma \mathbf{L}})^{-1} \mathbf{1} = \mathbf{1}$ allows to conclude. Note however that there is no such principle at the macroscopic scale, similarly to the continuous Saint-Venant system.

## 2.5. Implementation and computational costs

It is important to keep a reasonable algorithmic complexity so that the implicit method presented in the previous lines and especially the formulae (2.20), (2.21) can be used in practice. We discuss here how to improve its computational cost by a substantial margin. In Appendix A, we show that the $i$-th component of vectors $h^{\text{int}}$ and $(hu)^{\text{int}}$ have the form

$$ \left\{ \begin{aligned} h_i^{\text{int}} &= \frac{1}{2\sigma\sqrt{3}} \left( \sum_{j=i}^{P} \sqrt{\frac{2h_j^n}{g}} (Ah)_{i,j} + \sum_{j=1}^{i} \sqrt{\frac{2h_j^n}{g}} (Bh)_{i,j} \right) \\ (hu)_i^{\text{int}} &= \frac{1}{2\sigma^2\sqrt{3}} \left( -\sum_{j=i}^{P} \sqrt{\frac{2h_j^n}{g}} (Ahu)_{i,j} + \sum_{j=1}^{i} \sqrt{\frac{2h_j^n}{g}} (Bhu)_{i,j} \right), \end{aligned} \right. \qquad (2.26) $$

where $Ah, Ahu$ are dense upper triangular matrices, and $Bh, Bhu$ are dense lower triangular matrices. Therefore computing $h^{\text{int}}$ and $(hu)^{\text{int}}$ through (2.26) is analog to performing a matrix-vector product which has a quadratic complexity $O(P^2)$, and we cannot hope to do better than that. However the coefficients (A.3)–(A.6) of the above matrices involve a summation, and at a first glance the cost to assemble them is seemingly cubic. This is quite expensive and can render the method pretty much inefficient. However this complexity can be reduced to a
where $a_j^n = u_j^n - \sqrt{3} c_j^n$ and $b_j^n = u_j^n + \sqrt{3} c_j^n$. This corresponds to the recursive definition below

$$ (2.27) \qquad (Ah)_{i,j} = \begin{cases} 0 & \text{if } j < i \\ \left[ \ln(|1+x|) \right]_{-\min(0,a_j)\sigma}^{-\min(0,a_j)\sigma} & \text{if } i=j \\ (Ah)_{i+1,j} - \frac{1}{j-i} \left[ y^{j-i} \right]_{-\min(0,b_j)\sigma}^{-\min(0,a_j)\sigma} & \text{if } j > i \end{cases}. $$

Likewise, the lower triangular matrix $Bh$ is given by

$$ \left( \begin{array}{cccc} [z]_{\max(0, a_1^n)\sigma} & 0 & \cdots & 0 \\ [z-y]_{\max(0, a_1^n)\sigma} & \ddots & & \\ \vdots & & \ddots & \\ [z - \sum_{l=1}^{P-1} y^l/l]_{\max(0, a_1^n)\sigma} & \cdots & [z-y]_{\max(0, a_{P-1}^n)\sigma} & 0 \\ & & [z]_{\max(0, b_{P-1}^n)\sigma} & 0 \\ & & [z-y]_{\max(0, b_{P-1}^n)\sigma} & [z]_{\max(0, a_P^n)\sigma} \end{array} \right), $$

and can be defined by the following recurrence formula

$$ (2.28) \qquad (Bh)_{i,j} = \begin{cases} 0 & \text{if } i < j \\ \left[ \ln(|1+x|) \right]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma} & \text{if } i = j \\ (Bh)_{i-1,j} - \frac{1}{i-j} \left[ y^{i-j} \right]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma} & \text{if } i > j \end{cases}. $$

Hence it is more efficient to assemble matrices $Ah$ and $Bh$ column wise, starting from the diagonal coefficient and moving towards the first or last row. This way we only have to subtract one term to the previous coefficient so as to get the next one, and the cost of this operation is in $O(1)$. Since there are $P(P+1)/2$ coefficients to compute in total, the assembly of $Ah$ and $Bh$ following this strategy requires $O(P^2)$ steps.

A similar conclusion is achieved for $Ahu$ and $Bhu$, although the recurrence relation is less straightforward to obtain. We first remark that, introducing $(l)_{i,j} = i - j + 1$ the relations (A.5) and (A.6) become

$$ (Ahu)_{i,j} = \mathbf{1}_{j \ge i} \left[ -(l)_{j,i} \ln|1+x| + x + \sum_{k=1}^{j-i} k \frac{y^{(l)_{j,i}-k}}{(l)_{j,i}-k} \right]_{-\min(0,a_j)\sigma}^{-\min(0,b_j)\sigma}, $$

$$ (Bhu)_{i,j} = \mathbf{1}_{i \ge j} \left[ -(l)_{i,j} \ln|1+x| + x + \sum_{k=1}^{i-j} k \frac{y^{(l)_{i,j}-k}}{(l)_{i,j}-k} \right]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma}. $$

Performing the change of index $r = (l)_{j,i} - k$ for matrix $Ahu$ and $s = (l)_{i,j} - k$ for matrix $Bhu$ we find

$$ (Ahu)_{i,j} = \left[ -(l)_{i,j} \ln|1+x| + x + (l)_{i,j} \sum_{r=1}^{j-i} \frac{y^r}{r} - \sum_{r=1}^{j-i} y^r \right]_{-\min(0,a_j)\sigma}^{\min(0,b_j)\sigma}, $$

$$ (Bhu)_{i,j} = \left[ -(l)_{i,j} \ln|1+x| + x + (l)_{i,j} \sum_{s=1}^{i-j} \frac{y^s}{s} - \sum_{s=1}^{i-j} y^s \right]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma}. $$

Next we introduce the matrices defined column wise in a recursive manner

$$ (UA)_{i,j} = \begin{cases} 0 & \text{if } j \le i \\ (UA)_{i+1,j} + \left[ y^{j-i} \right]_{-\min(0,a_j)\sigma}^{-\min(0,b_j)\sigma} & \text{if } j > i \end{cases}, $$
quadratic cost by computing the coefficients in the correct order. More specifically we show that all the matrices above can be defined through a recurrence relation allowing to compute each coefficient from a previous one in $O(1)$ operation. In fact, denoting $y = x/(1+x)$ and $z = \ln|1+x|$, the matrix $Ah$ is given by

$$ \begin{pmatrix} [z]_{-\min(0,b_1)\sigma}^{-\min(0,a_1)\sigma} & [z-y]_{-\min(0,b_2)\sigma}^{-\min(0,a_2)\sigma} & \cdots & [z - \sum_{l=1}^{P-1} y^l/l]_{-\min(0,b_P)\sigma}^{-\min(0,a_P)\sigma} \\ 0 & [z]_{-\min(0,b_2)\sigma}^{-\min(0,a_2)\sigma} & & \vdots \\ \vdots & & \ddots & \vdots \\ & & & [z-y]_{-\min(0,b_P)\sigma}^{-\min(0,a_P)\sigma} \\ 0 & \cdots & \cdots & 0 \\ [z]_{-\min(0,b_P)\sigma}^{-\min(0,a_P)\sigma} & [z]_{-\min(0,b_P)\sigma}^{-\min(0,a_P)\sigma} & \cdots & 0 \end{pmatrix}, $$

where $a_j^n = u_j^n - \sqrt{3} c_j^n$ and $b_j^n = u_j^n + \sqrt{3} c_j^n$. This corresponds to the recursive definition below

$$ (Ah)_{i,j} = \begin{cases} 0 & \text{if } j < i \\ \left[ \ln(|1+x|) \right]_{-\min(0,b_j)\sigma}^{-\min(0,a_j)\sigma} & \text{if } i = j \\ (Ah)_{i+1,j} - \frac{1}{j-i} [y^{j-i}]_{-\min(0,b_j)\sigma}^{-\min(0,a_j)\sigma} & \text{if } j > i. \end{cases} \quad (2.27) $$

Likewise, the lower triangular matrix $Bh$ is given by

$$ \begin{pmatrix} [z]_{\max(0,a_1^n)\sigma} & 0 & \cdots & \cdots & 0 \\ [z-y]_{\max(0,a_1^n)\sigma} & \ddots & & & \vdots \\ \vdots & & \ddots & & \vdots \\ [z]_{\max(0,a_{P-1}^n)\sigma} & 0 & & & \vdots \\ [z - \sum_{l=1}^{P-1} y^l/l]_{\max(0,a_1^n)\sigma} & \cdots & \cdots & [z-y]_{\max(0,a_{P-1}^n)\sigma} & \cdots \\ [z]_{\max(0,a_{P-1}^n)\sigma} & [z-y]_{\max(0,a_{P-1}^n)\sigma} & \cdots & [z]_{\max(0,a_P^n)\sigma} & \cdots \end{pmatrix}, $$

and can be defined by the following recurrence formula

$$ (Bh)_{i,j} = \begin{cases} 0 & \text{if } i < j \\ \left[ \ln(|1+x|) \right]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma} & \text{if } i = j \\ (Bh)_{i-1,j} - \frac{1}{i-j} [y^{i-j}]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma} & \text{if } i > j. \end{cases} \quad (2.28) $$

Hence it is more efficient to assemble matrices $Ah$ and $Bh$ column wise, starting from the diagonal coefficient and moving towards the first or last row. This way we only have to subtract one term to the previous coefficient so as to get the next one, and the cost of this operation is in $O(1)$. Since there are $P(P+1)/2$ coefficients to compute in total, the assembly of $Ah$ and $Bh$ following this strategy requires $O(P^2)$ steps.

A similar conclusion is achieved for $Ahu$ and $Bhu$, although the recurrence relation is less straightforward to obtain. We first remark that, introducing $(l)_{i,j} = i - j + 1$ the relations (A.5) and (A.6) become

$$ (Ahu)_{i,j} = \mathbb{1}_{j \ge i} \left[ -(l)_{j,i} \ln|1+x| + x + \sum_{k=1}^{j-i} k \frac{y^{(l)_{j,i}-k}}{(l)_{j,i}-k} \right]_{-\min(0,b_j)\sigma}^{-\min(0,a_j)\sigma}, $$
$$ (VA)_{i,j} = \begin{cases} 0 & \text{if } j \le i \\ (VA)_{i+1,j} + \left[ \frac{y^{j-i}}{j-i} \right]^{-\min(0,a_j)\sigma} & \text{if } j > i \end{cases}. $$

Then we can write that

$$ (2.29) \qquad (Ahu)_{i,j} = \begin{cases} 0 & j < i \\ [x - \ln|1 + x|]^{-\min(0,a_j)\sigma} & j = i \\ (l)_{j,i}(VA)_{i,j} - (UA)_{i,j} + [x - (l)_{j,i} \ln|1 + x|]^{-\min(0,a_j)\sigma} & j > i \end{cases}. $$

Similarly we introduce

$$ (UB)_{i,j} = \begin{cases} 0 & \text{if } i \le j \\ (UB)_{i-1,j} + [y^{i-j}]_{\max(0,a_j)\sigma} & \text{if } i > j \end{cases}, $$

$$ (VB)_{i,j} = \begin{cases} 0 & \text{if } i \le j \\ (VB)_{i-1,j} + \left[ \frac{y^{i-j}}{i-j} \right]_{\max(0,a_j)\sigma} & \text{if } i > j \end{cases}, $$

so that we have

$$ (2.30) \qquad (Bhu)_{i,j} = \begin{cases} 0 & i < j \\ [x - \ln|1 + x|]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma} & i = j \\ (l)_{i,j}(VB)_{i,j} - (UB)_{i,j} + [x - (l)_{i,j} \ln|1 + x|]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma} & i > j \end{cases}. $$

To conclude, through relations (2.29) and (2.30) we are also able to assemble matrices *Ahu* and *Bhu* with a quadratic cost with respect to the number of cells, which means that the overall method has a $O(P^2)$ complexity.

**Remark 2.7** We have considered here the specific case of a kinetic solver and one can imagine that an implicit scheme for another finite volume solver can lead to reduced numerical costs. But it is worth noticing that since the explicit expression of the inverse of the matrix $I + \sigma L$ is accessible in the kinetic context, one can hardly find a more efficient implicit technique.

Obviously the proposed implicit scheme is not constrained by any CFL condition associated with an explicit scheme, nevertheless it is important to compare the computational costs of the explicit and implicit strategies in the context of a kinetic solver. This comparison is performed below in the case of a flat topography.

**Explicit scheme.** Let $\Delta t^n$ be the time step allowing to satisfy the CFL constraint. In order to obtain the expression of $U^{n+1}$ from $U^n$, approximately $4P$ numerical fluxes have to be computed (2 numerical fluxes at each interface for each variable $h$ and $hu$). The explicit kinetic scheme is fully detailed in [6, 4].

**Implicit scheme.** The CFL constraint being relaxed, we can consider a time step $\Delta t_{imp}^n \gg \Delta t^n$. The results obtained in this paragraph shows that the update from $U^{n+1}$ from $U^n$ requires approximately $P^2$ numerical fluxes to compute.

We conclude that the implicit strategy is less expensive when

$$ (2.31) \qquad \frac{\Delta t_{imp}^n}{\Delta t^n} \gg \frac{P^2}{4P} = \frac{P}{4}. $$

Note however that the computational cost is not the only factor to account for, and one should also consider the efficiency of the scheme, that is to say the relation between the error and the computational time. Generally, taking a very coarse resolution in time results in poorly accurate results, in which case it is not desirable to have (2.31). Nevertheless there are some cases where the fast dynamics
$$ (Bhu)_{i,j} = \mathbb{1}_{i \ge j} \left[ -(l)_{i,j} \ln|1+x| + x + \sum_{k=1}^{i-j} k \frac{y^{(l)_{i,j}-k}}{(l)_{i,j}-k} \right]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma} . $$

Performing the change of index $r = (l)_{j,i} - k$ for matrix *Ahu* and $s = (l)_{i,j} - k$ for matrix *Bhu* we find

$$ (Ahu)_{i,j} = \left[ -(l)_{i,j} \ln|1+x| + x + (l)_{i,j} \sum_{r=1}^{j-i} \frac{y^r}{r} - \sum_{r=1}^{j-i} y^r \right]_{-\min(0,b_j)\sigma}^{-\min(0,a_j)\sigma}, $$

$$ (Bhu)_{i,j} = \left[ -(l)_{i,j} \ln|1+x| + x + (l)_{i,j} \sum_{s=1}^{i-j} \frac{y^s}{s} - \sum_{s=1}^{i-j} y^s \right]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma}. $$

Next we introduce the matrices defined column wise in a recursive manner

$$ (UA)_{i,j} = \begin{cases} 0 & \text{if } j \le i \\ (UA)_{i+1,j} + [y^{j-i}]_{-\min(0,b_j)\sigma}^{-\min(0,a_j)\sigma} & \text{if } j > i, \end{cases} $$

$$ (VA)_{i,j} = \begin{cases} 0 & \text{if } j \le i \\ (VA)_{i+1,j} + \left[ \frac{y^{j-i}}{j-i} \right]_{-\min(0,b_j)\sigma}^{-\min(0,a_j)\sigma} & \text{if } j > i. \end{cases} $$

Then we can write that

$$ (Ahu)_{i,j} = \begin{cases} 0 & j < i \\ [x - \ln|1 + x|]_{-\min(0,b_j)\sigma}^{-\min(0,a_j)\sigma} & j = i \\ (l)_{j,i}(VA)_{i,j} - (UA)_{i,j} + [x - (l)_{j,i} \ln|1 + x|]_{-\min(0,b_j)\sigma}^{-\min(0,a_j)\sigma} & j > i. \end{cases} \quad (2.29) $$

Similarly we introduce

$$ (UB)_{i,j} = \begin{cases} 0 & \text{if } i \le j \\ (UB)_{i-1,j} + [y^{i-j}]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma} & \text{if } i > j, \end{cases} $$

$$ (VB)_{i,j} = \begin{cases} 0 & \text{if } i \le j \\ (VB)_{i-1,j} + \left[ \frac{y^{i-j}}{i-j} \right]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma} & \text{if } i > j, \end{cases} $$

so that we have

$$ (Bhu)_{i,j} = \begin{cases} 0 & i < j \\ [x - \ln|1 + x|]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma} & i = j \\ (l)_{i,j}(VB)_{i,j} - (UB)_{i,j} + [x - (l)_{i,j} \ln|1 + x|]_{\max(0,a_j)\sigma}^{\max(0,b_j)\sigma} & i > j. \end{cases} \quad (2.30) $$

To conclude, through relations (2.29) and (2.30) we are also able to assemble matrices *Ahu* and *Bhu* with a quadratic cost with respect to the number of cells, which means that the overall method has a $O(P^2)$ complexity.

**Remark 2.7.** We have considered here the specific case of a kinetic solver and one can imagine that an implicit scheme for another finite volume solver can lead to reduced numerical costs. But it is worth noticing that since the explicit expression of the inverse of the matrix $I + \sigma L$ is accessible in the kinetic context, one can hardly find a more efficient implicit technique.
do not play an important role such as in the low Froude regime. Then it might be advantageous to consider large time steps. We will see through the upcoming numerical results from Section 5.1 that the interest of the implicit kinetic scheme is rather limited when it comes to efficiency, at least for the considered test cases. Hence the explicit strategy is preferable to the implicit one, unless we account for the greater stability offered by the latter in terms of discrete entropy inequality.

## 3 An iterative resolution scheme

The kinetic scheme (2.3) requires to solve a linear system and in the previous section, we have seen that it was possible to have an analytic expression for the inverse of the matrix $I + \sigma L$ with $L$ given in (2.4). For the numerical approximation of PDEs e.g. in finite elements methods when the linear system to solve is large an iterative strategy is singled out compared to a direct inversion of the matrix. We propose to follow the same idea here, with mainly two benefits. First it will allow us to use the half disk Maxwellian (1.4), for which we recall the integrals (2.19) could not be computed analytically in the case of the fully implicit kinetic scheme. This is important as it will enable to prove some discrete entropy inequality at the macroscopic scale thanks to the existence of a kinetic entropy $H_0$ given by (1.10) that satisfies (1.13)-(1.14), while having an explicit writing of the update. The second advantage lies in the possibility to couple the iterative strategy with the hydrostatic reconstruction to obtain a well balanced treatment for varying bottoms, which we will discuss in Section 3.2. Before this, the simpler case of flat bathymetries is investigated as a toy problem in Section 3.1, and will give some insight.

More precisely we shall use a Gauss-Jacobi type decomposition, which consists to write the matrix $I + \sigma L$ as $D - N$, where $D$ and $N$ are two matrices from $\mathbb{R}^{P \times P}$ with $D$ is invertible. Then the scheme (2.3) also writes

$$Df^{n+1} = Nf^{n+1} + M + \sigma B^{n+1},$$

and one can propose an iterative resolution of the previous equation under the form

$$f^{n+1,k+1} = D^{-1}Nf^{n+1,k} + D^{-1}(M + \sigma^k B^{n+1}), \quad \sigma^k = \frac{\Delta t^k}{\Delta x},$$

where the timestep $\Delta t^k$ is allowed to vary at each subiterations. If it converges, the sequence $(f^{n+1,k})_{k \in \mathbb{N}}$ converges towards the solution of (2.3).

### 3.1 Case of a flat topography

In this section, we study iterative strategies with the particular choice

$$D = (1 + \alpha)I, \quad \text{and} \quad N = \alpha I - \sigma^k L,$$

when the bathymetry is flat and where $\alpha \in (0, \infty)$ is a relaxation parameter. The main goal is to put forward the differences between such iterative methods and the fully implicit kinetic scheme proposed in Section 2; the first one being that it will no longer be possible to use arbitrary large time steps since a CFL condition will be needed to ensure the convergence of the sub-iterations and the positivity of the water height; the second one being that the iterative scheme enables one to have an explicit writing of the update while using the half-disk Maxwellian. As long as we are restricted to flat bathymetries, it should also be noticed that an iterative strategy proves less advantageous than the forward Euler kinetic scheme, which is more efficient while preserving the water height positivity and being entropy stable. The iterative kinetic scheme presents however a clear advantage over the explicit scheme in terms of entropy stability in the case of varying bathymetries, which is treated in Section 3.2.

We consider two iterative schemes below; they differ only in wether or not the vector of kinetic densities $f^{n+1,k}(\xi)$ is projected onto the space of Maxwellians at every sub-iteration or just at $k = 0$. In the latter and easier setting, we shall obtain a CFL condition to ensure the convergence of sub-iterations. In the former case, this question is also addressed, although more restrictive assumptions are needed. The question of entropy stability is only investigated when the projection step occurs at every sub-iteration.
Obviously the proposed implicit scheme is not constrained by any CFL condition associated with an explicit scheme, nevertheless it is important to compare the computational costs of the explicit and implicit strategies in the context of a kinetic solver. This comparison is performed below in the case of a flat topography.

**Explicit scheme.** Let $\Delta t^n$ be the time step allowing to satisfy the CFL constraint. In order to obtain the expression of $U^{n+1}$ from $U^n$, approximately $4P$ numerical fluxes have to be computed (2 numerical fluxes at each interface for each variable $h$ and $hu$). The explicit kinetic scheme is fully detailed in [4, 6].

**Implicit scheme.** The CFL constraint being relaxed, we can consider a time step $\Delta t_{imp}^n \gg \Delta t^n$. The results obtained in this paragraph shows that the update from $U^{n+1}$ from $U^n$ requires approximately $P^2$ numerical fluxes to compute.

We conclude that the implicit strategy is less expensive when

$$ \frac{\Delta t_{imp}^n}{\Delta t^n} \gg \frac{P^2}{4P} = \frac{P}{4}. \qquad (2.31) $$

Note however that the computational cost is not the only factor to account for, and one should also consider the efficiency of the scheme, that is to say the relation between the error and the computational time. Generally, taking a very coarse resolution in time results in poorly accurate results, in which case it is not desirable to have (2.31). Nevertheless there are some cases where the fast dynamics do not play an important role such as in the low Froude regime. Then it might be advantageous to consider large time steps. We will see through the upcoming numerical results from Section 5.1 that the interest of the implicit kinetic scheme is rather limited when it comes to efficiency, at least for the considered test cases. Hence the explicit strategy is preferable to the implicit one, unless we account for the greater stability offered by the latter in terms of discrete entropy inequality. In this regard, a semi-implicit in time strategy could ensure the stability of the scheme with a better efficiency (see for instance [19]).

### 3. AN ITERATIVE RESOLUTION SCHEME

The kinetic scheme (2.3) requires to solve a linear system and in the previous section, we have seen that it was possible to have an analytic expression for the inverse of the matrix $I + \sigma L$ with $L$ given in (2.4). For the numerical approximation of PDEs e.g. in finite elements methods when the linear system to solve is large an iterative strategy is singled out compared to a direct inversion of the matrix. We propose to follow the same idea here, with mainly two benefits. First it will allow us to use the half disk Maxwellian (1.4), for which we recall the integrals (2.19) could not be computed analytically in the case of the fully implicit kinetic scheme. This is important as it will enable to prove some discrete entropy inequality at the macroscopic scale thanks to the existence of a kinetic entropy $H_0$ given by (1.10) that satisfies (1.13) and (1.14), while having an explicit writing of the update. The second advantage lies in the possibility to couple the iterative strategy with the hydrostatic reconstruction to obtain a well balanced treatment for varying bottoms, which we will discuss in Section 3.2. Before this, the simpler case of flat bathymetries is investigated as a toy problem in Section 3.1, and will give some insight.

More precisely we shall use a Gauss-Jacobi type decomposition, which consists to write the matrix $I + \sigma L$ as $D - N$, where $D$ and $N$ are two matrices from $\mathbb{R}^{P \times P}$ with $D$ is invertible. Then the scheme (2.3) also writes

$$ Df^{n+1} = Nf^{n+1} + M + \sigma B^{n+1}, $$

and one can propose an iterative resolution of the previous equation under the form

$$ f^{n+1,k+1} = D^{-1}Nf^{n+1,k} + D^{-1}(M + \sigma^k B^{n+1}), \quad \sigma^k = \frac{\Delta t^k}{\Delta x}, $$

where the time step $\Delta t^k$ is allowed to vary at each sub-iteration in order to satisfy some CFL condition. If it converges, the sequence $(f^{n+1,k})_{k \in \mathbb{N}}$ converges towards the solution of (2.3).
When developed, the iterative process with a single projection at $k = 0$ reads:

$$ (3.1) \qquad \left\{ \begin{array}{l} f^{n+1,0} = M \\ (1+\alpha)f^{n+1,k+1} = (\alpha\mathbf{I} - \sigma^k\mathbf{L})f^{n+1,k} + M + \sigma^k B^{n+1,k}, \\ \forall 1 \le i \le P, \quad U_i^{n+1,k} = \int_{\mathbb{R}} \left(\frac{1}{\xi}\right) f_i^{n+1,k}(\xi) \, d\xi \end{array} \right. , $$

with $B^{n+1,k}$ the boundary condition associated with the macroscopic state $U^{n+1,k}$ as explained in Section 2.4, see (2.5). The following Proposition highlights the main compromise linked with such an iterative approach, which is the requirement for a CFL condition in order for the method to converge.

**Proposition 3.1** Assume that $B^{n+1,k}$ remains constant equal to $B^n$ for any $k$ in $\mathbb{N}$. Then (3.1) defines an arithmetico-geometric sequence which converges if the CFL condition $\sigma|\xi| < 1 + 2\alpha$ holds for all $\xi$ belonging to $\text{supp } M \cup \text{supp } B^n$.

**Proof.** By recurrence, we can show that for any $k \in \mathbb{N}$ the support of $f^{n+1,k}$ is included in $\text{supp } M \cup \text{supp } B^n$, which is why we restrict to velocities $\xi$ belonging to this set. Consider $f$ the solution of

$$ f = \mathbf{D}^{-1}\mathbf{N}f + \mathbf{D}^{-1}(M + \sigma^k B^n). $$

The sequence $(g^k)_k$ defined by $g^k = f^{n+1,k} - f$ satisfies $g^{k+1} = \mathbf{D}^{-1}\mathbf{N}g^k$ and converges to zero as soon as the spectral radius of $\mathbf{D}^{-1}\mathbf{N}$ is strictly less than one. Since $\mathbf{D}^{-1}\mathbf{N}$ is a triangular matrix, its eigenvalues are given by its diagonal coefficients, all equal to $(1+\alpha)^{-1}(\alpha - \sigma^k|\xi|)$. Under the assumption $\sigma^k|\xi| < 1 + 2\alpha$, this quantity is strictly less than one in absolute value, which concludes the proof. $\square$

**Remark 3.2** As we did in Section 2.4 for the fully implicit scheme, we can replace $B^{n+1,k}$ by $B^n$ in the iterative process (3.1). In fact this constitutes a first order approximation in time since we have $f^{n+1,k} = M + \mathcal{O}(\Delta t)$. Under this simplification, the assumption from Proposition 3.1 that $B^{n+1,k}$ does not depend on $k$ becomes automatically satisfied.

In practice, we wish to apply an iterative method directly at the macroscopic level. An issue with (3.1) is that the distribution involved in the kinetic flux (i.e. the term in factor of $\sigma^k\mathbf{L}$) is not a vector of Maxwellians, which prevents us to write the recurrence relation at the macroscopic level since there is no general expression for the numerical flux. To bypass this issue, we propose the following modification of (3.1), where we replace all occurrences of $f^{n+1,k}$ on the right hand side by a vector of Maxwellians $M^{n+1,k}$, which yields

$$ (3.2) \qquad \left\{ \begin{array}{l} \bar{f}^{n+1,0}(\xi) = M \\ (1+\alpha)\bar{f}^{n+1,k+1}(\xi) = (\alpha\mathbf{I} - \sigma^k\mathbf{L})\bar{M}^{n+1,k} + M + \sigma^k \bar{B}^{n+1,k}. \\ \bar{M}^{n+1,k+1} = \bar{f}^{n+1,k+1} + \Delta t^k \bar{Q}^{n+1,k+1} \end{array} \right. . $$

This new iterative process is alternating two stages, the first one being the usual transport step, while the second one is a projection step onto the set of Maxwellians yielding $\bar{M}^{n+1,k+1}$. In this sense (3.2) is an iterative BGK splitting approach where the projection step doesn't modify the macroscopic quantities of interest since the term $\bar{Q}^{n+1,k}$ is a vector of collision operators satisfying the conservation constraints (1.7). Note that the time stepping $\Delta t^k$ is made dependent on $k$ as the support of $\bar{M}^{n+1,k}$ can change from iteration to iteration.

**Remark 3.3** For each value of $k$, each row of the density $\bar{f}^{n+1,k+1}(\xi)$ defined by (3.2) is a combination of half disks. Hence $\bar{f}_i^{n+1,k+1}(\xi)$ the that is the row $i$ of $\bar{f}^{n+1,k+1}(\xi)$ is a combination of $M_i^{n+1,k}$, $M_i$ and possibly $B_i^{n+1,k}$ if $i = 1$ or $i = P$, see Fig. 1. This illustrates the role of the collision term $\bar{Q}^{n+1,k+1}$.

The practical implementation of scheme (3.2) is based on its macroscopic version given by

$$ (3.3) \qquad (1+\alpha)U_i^{n+1,k+1} = \alpha U_i^{n+1,k} + U_i - \sigma \left( \mathcal{F}(U_i^{n+1,k}, U_{i+1}^{n+1,k}) - \mathcal{F}(U_{i-1}^{n+1,k}, U_i^{n+1,k}) \right), $$
### 3.1. Case of a flat topography

In this section, we study iterative strategies with the particular choice

$$ \mathbf{D} = (1 + \alpha)\mathbf{I}, \quad \text{and} \quad \mathbf{N} = \alpha\mathbf{I} - \sigma^k\mathbf{L}, $$

when the bathymetry is flat and where $\alpha \in (0, \infty)$ is a relaxation parameter. The main goal is to put forward the differences between such iterative methods and the fully implicit kinetic scheme proposed in Section 2; the first one being that it will no longer be possible to use arbitrary large time steps since a CFL condition will be needed to ensure the convergence of the sub-iterations and the positivity of the water height; the second one being that the iterative scheme enables one to have an explicit writing of the update while using the half-disk Maxwellian. As long as we are restricted to flat bathymetries, it should also be noticed that an iterative strategy proves less advantageous than the forward Euler kinetic scheme, which is more efficient while preserving the water height positivity and being entropy stable. The iterative kinetic scheme presents however a clear advantage over the explicit scheme in terms of entropy stability in the case of varying bathymetries, which is treated in Section 3.2.

We consider two iterative schemes below; they differ only in whether or not the vector of kinetic densities $f^{n+1,k}(\xi)$ is projected onto the space of Maxwellians at every sub-iteration or just at $k = 0$. In the latter and easier setting, we shall obtain a CFL condition to ensure the convergence of sub-iterations. In the former case, this question is also addressed, although more restrictive assumptions are needed. The question of entropy stability is only investigated when the projection step occurs at every sub-iteration.

When developed, the iterative process with a single projection at $k = 0$ reads:

$$ \left\{ \begin{array}{l} f^{n+1,0} = M \\ (1 + \alpha) f^{n+1,k+1} = (\alpha \mathbf{I} - \sigma^k \mathbf{L}) f^{n+1,k} + M + \sigma^k B^{n+1,k}, \\ \forall 1 \le i \le P, \quad U_i^{n+1,k} = \int_{\mathbb{R}} \left(\frac{1}{\xi}\right) f_i^{n+1,k}(\xi) \, d\xi \end{array} \right. \qquad (3.1) $$

with $B^{n+1,k}$ the boundary condition associated with the macroscopic state $U^{n+1,k}$ as explained in Section 2.4, see (2.5). The following proposition highlights the main compromise linked with such an iterative approach, which is the requirement for a CFL condition in order for the method to converge.

**Proposition 3.1.** Assume that $B^{n+1,k}$ remains constant equal to $B^n$ for any $k$ in $\mathbb{N}$. Then (3.1) defines an arithmetico-geometric sequence which converges if the CFL condition $\sigma|\xi| < 1 + 2\alpha$ holds for all $\xi$ belonging to $\text{supp } M \cup \text{supp } B^n$.

*Proof.* By recurrence, we can show that for any $k \in \mathbb{N}$ the support of $f^{n+1,k}$ is included in $\text{supp } M \cup \text{supp } B^n$, which is why we restrict to velocities $\xi$ belonging to this set. Consider $f$ the solution of

$$ f = \mathbf{D}^{-1}\mathbf{N}f + \mathbf{D}^{-1}(M + \sigma B^n). $$

The sequence $(g^k)_k$ defined by $g^k = f^{n+1,k} - f$ satisfies $g^{k+1} = \mathbf{D}^{-1}\mathbf{N}g^k$ and converges to zero as soon as the spectral radius of $\mathbf{D}^{-1}\mathbf{N}$ is strictly less than one. Since $\mathbf{D}^{-1}\mathbf{N}$ is a triangular matrix, its eigenvalues are given by its diagonal coefficients, all equal to $(1 + \alpha)^{-1}(\alpha - \sigma^k|\xi|)$. Under the assumption $\sigma^k|\xi| < 1 + 2\alpha$, this quantity is strictly less than one in absolute value, which concludes the proof. $\square$

**Remark 3.2.** As we did in Section 2.4 for the fully implicit scheme, we can replace $B^{n+1,k}$ by $B^n$ in the iterative process (3.1). In fact this constitutes a first order approximation in time since we have $f^{n+1,k} = M + \mathcal{O}(\Delta t)$. Under this simplification, the assumption from Proposition 3.1 that $B^{n+1,k}$ does not depend on $k$ becomes automatically satisfied.

In practice, we wish to apply an iterative method directly at the macroscopic level. An issue with (3.1) is that the distribution involved in the kinetic flux (i.e. the term in factor of $\sigma^k\mathbf{L}$) is not a vector of Maxwellians,
<figure><img src="image_2.png" /><figcaption>Figure 1: The maxwellians $M_i$, $M_{i-1}$ and $M_{i+1}$. For each quantities appearing over the figure, the superscripts have been omitted.</figcaption></figure>

for all $1 \le i \le P$, where the numerical flux $\mathcal{F}$ is defined as

$$ (3.4) \qquad \mathcal{F}(U_L, U_R) = \int_{\mathbb{R}} \xi \left(\frac{1}{\xi}\right) \left(1_{\xi>0} M(U_L, \xi) + 1_{\xi<0} M(U_R, \xi)\right) d\xi, $$

and where the vectors $U_0^{n+1,k}$, $U_{P+1}^{n+1,k}$ appearing in the border cells $i \in \{1, P\}$ are respectively functions of $U_1^{n+1,k}$ and $U_P^{n+1,k}$ since the boundary conditions are imposed though a ghost cell strategy fully described in [13, 1]. As previously, one can choose to make these ghost values independent from the sub-iteration index $k$. Notice that if the sequence $(U^{n+1,k})_{k \in \mathbb{N}} \subset (\mathbb{R}^2)^P$ from (3.3) converges in $(\mathbb{R}^2)^P$, its limit $U^{n+1}$ then satisfies

$$ (3.5) \qquad \forall 1 \le i \le P, \quad U_i^{n+1} = U_i^n - \sigma \left( \mathcal{F}(U_i^{n+1}, U_{i+1}^{n+1}) - \mathcal{F}(U_{i-1}^{n+1}, U_i^{n+1}) \right) $$

by continuity of the numerical flux (3.4).

Because of the projection step, the convergence of the iterative method (3.2) is more complex to prove than in Prop. 3.1. But the following result holds.

**Proposition 3.4** Let the sequence $\{\bar{f}^{n+1,k}(\cdot)\}_{k \in \mathbb{N}}$ be defined through (3.2), and let $\{\bar{Q}^{n+1,k}(\cdot)\}_{k \in \mathbb{N}}$ be the associated sequence of collision terms assumed to be bounded. Then under the CFL condition $2\sigma^k|\xi| < 1$ the sequence $\{\bar{f}^{n+1,k}(\cdot)\}_{k \in \mathbb{N}}$ converges.

**Proof of Prop. 3.4.** First, we see that under the CFL condition $\sigma^k|\xi| \le \alpha$, the quantity $\bar{f}^{n+1,k}$ defined by (3.2) is non-negative for any $\xi$.

Note that in (3.2) the convergence of $\{\bar{f}^{n+1,k}\}_{k \in \mathbb{N}}$ implies that of $\{\bar{Q}^{n+1,k}\}_{k \in \mathbb{N}}$, therefore we focus on the former sequence. The recurrence formula (3.2) also writes

$$ (3.6) \qquad (1 + \alpha) M^{n+1,k+1} = (\alpha \mathbf{I} - \sigma^k \mathbf{L}) M^{n+1,k} + M + \sigma^k \bar{B}^{n+1,k} + (1 + \alpha) \Delta t^k \bar{Q}^{n+1,k+1}. $$

The previous formula is nothing else than an arithmetico-geometric-like sequence having the form

$$ M^{n+1,k+1} = \mathbf{A}^{k,\alpha} M^{n+1,k} + C^{n,k}, $$

with

$$ (3.7) \qquad \begin{aligned} \mathbf{A}^{k,\alpha} &= \frac{\alpha \mathbf{I} - \sigma^k \mathbf{L}}{1 + \alpha}, \\ C^{n,k} &= \frac{M + \sigma^k \bar{B}^{n+1,k}}{1 + \alpha} + \Delta t^k \bar{Q}^{n+1,k+1}, \end{aligned} $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Example of Maxwellians $M_i$, $M_{i-1}$ and $M_{i+1}$. For each quantity appearing over the figure, the superscripts have been omitted.</figcaption></figure>

which prevents us to write the recurrence relation at the macroscopic level since there is no general expression for the numerical flux. To bypass this issue, we propose the following modification of (3.1), where we replace all occurrences of $f^{n+1,k}$ on the right hand side by a vector of Maxwellians $M^{n+1,k}$, which yields

$$ \left\{ \begin{aligned} \bar{f}^{n+1,0}(\xi) &= M \\ (1+\alpha)\bar{f}^{n+1,k+1}(\xi) &= (\alpha\mathbf{I} - \sigma^k\mathbf{L})M^{n+1,k} + M + \sigma^k \bar{B}^{n+1,k} \\ M^{n+1,k+1} &= \bar{f}^{n+1,k+1} + \Delta t^k \bar{Q}^{n+1,k+1}. \end{aligned} \right. \qquad (3.2) $$

This new iterative process is alternating two stages, the first one being the usual transport step, while the second one is a projection step onto the set of Maxwellians yielding $M^{n+1,k+1}$. In this sense (3.2) is an iterative BGK splitting approach where the projection step doesn't modify the macroscopic quantities of interest since the term $\bar{Q}^{n+1,k}$ is a vector of collision operators satisfying the conservation constraints (1.7). Note that the time stepping $\Delta t^k$ is made dependent on $k$ as the support of $M^{n+1,k}$ can change from iteration to iteration.

**Remark 3.3.** For each value of $k$, each row of the density $\bar{f}^{n+1,k+1}(\xi)$ defined by (3.2) is a combination of half disks. Hence $\bar{f}_i^{n+1,k+1}(\xi)$ that is the row $i$ of $\bar{f}^{n+1,k+1}(\xi)$ is a combination of $M_i^{n+1,k}$, $M_i$ and possibly $B_i^{n+1,k}$ if $i=1$ or $i=P$, see Figure 1. This illustrates the role of the collision term $\bar{Q}^{n+1,k+1}$.

The practical implementation of scheme (3.2) is based on its macroscopic version given by

$$ (1 + \alpha) U_i^{n+1,k+1} = \alpha U_i^{n+1,k} + U_i - \sigma^k \left( \mathcal{F}(U_i^{n+1,k}, U_{i+1}^{n+1,k}) - \mathcal{F}(U_{i-1}^{n+1,k}, U_i^{n+1,k}) \right), \qquad (3.3) $$

for all $1 \le i \le P$, where the numerical flux $\mathcal{F}$ is defined as

$$ \mathcal{F}(U_L, U_R) = \int_{\mathbb{R}} \xi \left( \frac{1}{\xi} \right) \left( 1_{\xi>0} M(U_L, \xi) + 1_{\xi<0} M(U_R, \xi) \right) d\xi, \qquad (3.4) $$

and where the vectors $U_0^{n+1,k}$, $U_{P+1}^{n+1,k}$ appearing in the border cells $i \in \{1, P\}$ are respectively functions of $U_1^{n+1,k}$ and $U_P^{n+1,k}$ since the boundary conditions are imposed though a ghost cell strategy fully described in [1, 13]. As previously, one can choose to make these ghost values independent from the sub-iteration index $k$. Notice that if the sequence $(U^{n+1,k})_{k \in \mathbb{N}} \subset (\mathbb{R}^2)^P$ from (3.3) converges in $(\mathbb{R}^2)^P$, its limit $U^{n+1}$ then satisfies

$$ \forall 1 \le i \le P, \quad U_i^{n+1} = U_i^n - \sigma^k \left( \mathcal{F}(U_i^{n+1}, U_{i+1}^{n+1}) - \mathcal{F}(U_{i-1}^{n+1}, U_i^{n+1}) \right) \qquad (3.5) $$

by continuity of the numerical flux (3.4).

Because of the projection step, the convergence of the iterative method (3.2) is more complex to prove than in Proposition 3.1. But the following result holds.
whose convergence is ensured as soon as the sequence

$$ \sum_{j=0}^{k} (\mathbf{A}^{j,\alpha})^j C^{n,k-j}, $$

converges when $k$ grows. Since the collision terms $\bar{Q}^{n+1,k+1}$ are assumed to be bounded independently of $k$, we have $|C^{n,k}| \le \bar{C}$ and the CFL condition ensures that

$$ 0 \le \|\mathbf{A}^{k,\alpha}\|_{\infty,\infty} = \frac{\alpha + 2\sigma^k |\xi|}{1 + \alpha} < 1, $$

where $\|\cdot\|_{\infty,\infty}$ is the maximum subordinate matrix norm. This allows to write

$$ (3.8) \qquad \sum_{j=0}^{k} |(\mathbf{A}^{j,\alpha})^j C^{n,k-j}| \le \sum_{j=0}^{k} \|\mathbf{A}^{k,\alpha}\|_{\infty,\infty}^j |\bar{Q}^{n+1,k-j}| \le \frac{\bar{C}}{1 - \|\mathbf{A}^{k,\alpha}\|_{\infty,\infty}}, $$

giving the absolute convergence of the series and thus proving the result. $\square$

Next we state the following result regarding the positivity of the iterative kinetic scheme (3.2).

**Proposition 3.5** Assume that the water height vectors $h^n$ and $h^{n+1,k}$ are positive. Then the update $f^{n+1,k+1}$ defined in the iterative scheme (3.2) is positive if for all $1 \le i \le P$ the CFL condition $\sigma^k|\xi| \le \alpha + M_i/M_i^{n+1,k}$ holds for any $\xi$ belonging to supp $M_i^{n+1,k}$. As a direct consequence, the water height vector $h^{n+1,k+1}$ from scheme (3.3) is positive under these assumptions.

We postpone the proof of Proposition 3.5 to the next section, where it is generalized to the case with varying bottom in Proposition 3.9. Finally we investigate the entropy stability of the iterative kinetic scheme (3.2). Note that when restricting to a flat bathymetry, the explicit kinetic scheme [4] already verifies a discrete entropy inequality, while being more efficient since it avoids the use of sub-iterations. The benefit of the iterative kinetic scheme compared to the fully explicit one will become clearer in the next section, where it is combined with the hydrostatic reconstruction. Nevertheless, the simpler case of a flat bottom is interesting in order to understand why the iterative kinetic scheme involving the half-disk Maxwellian (1.4) offers a favorable setting to get a discrete entropy inequality.

**Proposition 3.6** The kinetic entropy of the iterative scheme (3.2) with half-disk Maxwellian (1.4) satisfies the following inequality

$$ (3.9) \qquad H_0(M_i^{n+1,k+1}) \le \frac{H_0(M_i) + \alpha H_0(M_i^{n+1,k})}{1 + \alpha} - \frac{\sigma^k \xi}{1 + \alpha} (H_{0,i+1/2}^{n+1,k} - H_{0,i-1/2}^{n+1,k}) \\ \qquad \qquad + \eta'(U_i^{n+1,k+1}) \cdot \left(\frac{1}{\xi}\right) Q_i^{n+1,k+1} + D_i^{n+1,k+1}, $$

with $Q_i^{n+1,k+1} = M_i^{n+1,k+1} - f_i^{n+1,k+1}$ a collision operator verifying the conservation constraints (1.7), and with $\eta$ the entropy given in (1.3). The interfacial kinetic entropies $H_{0,i\pm1/2}^{n+1,k}$ are

$$ H_{0,i-1/2}^{n+1,k} = \mathbb{1}_{\xi>0} H_0(M_{i-1}^{n+1,k}, \xi) + \mathbb{1}_{\xi<0} H_0(M_i^{n+1,k}, \xi), $$

$$ H_{0,i+1/2}^{n+1,k} = \mathbb{1}_{\xi>0} H_0(M_i^{n+1,k}, \xi) + \mathbb{1}_{\xi<0} H_0(M_{i+1}^{n+1,k}, \xi), $$

and the term $D_i^{n+1,k+1}$ is given by

$$ D_i^{n+1,k+1} = - \frac{1}{1+\alpha} \Psi(M_i^{n+1,k+1}, M_i) - \frac{\alpha - \sigma^k |\xi| \mathbb{1}_{\Xi}}{1+\alpha} \Psi(M_i^{n+1,k+1}, M_i^{n+1,k}) \\ - \frac{\sigma^k |\xi| \mathbb{1}_{\Xi}}{1+\alpha} \Psi(M_i^{n+1,k+1}, M_{i\pm1}^{n+1,k}), $$
**Proposition 3.4.** Let the sequence ${\bar{f}^{n+1,k}(\cdot)}_{k \in \mathbb{N}}$ be defined through (3.2), and let ${\bar{Q}^{n+1,k}(\cdot)}_{k \in \mathbb{N}}$ be the associated sequence of collision terms assumed to be bounded. Then under the CFL condition $\sigma^k|\xi| < \min(1, \alpha)$ the sequence ${\bar{f}^{n+1,k}(\cdot)}_{k \in \mathbb{N}}$ is non-negative and converges.

*Proof of Proposition 3.4.* First, we see that under the CFL condition $\sigma^k|\xi| \le \alpha$, the quantity $\bar{f}^{n+1,k}$ defined by (3.2) is non-negative for any $\xi$. In fact there holds

$$ (1 + \alpha) f_i^{n+1,k+1} = M_i + \alpha M_i^{n+1,k} - \sigma^k (\xi 1_{\xi \le 0} (M_{i+1}^{n+1,k} - M_i^{n+1,k}) - \xi 1_{\xi \ge 0} (M_{i-1}^{n+1,k} - M_i^{n+1,k})) \\ \ge M_i + (\alpha - \sigma^k |\xi|) M_i^{n+1,k} - \sigma^k \xi 1_{\xi \le 0} M_{i+1}^{n+1,k} + \sigma^k \xi 1_{\xi \ge 0} M_{i-1}^{n+1,k}, $$

with the last two terms of the lower bound non-negative. To obtain the convergence we write the recurrence formula (3.2) as

$$ (1 + \alpha) M^{n+1,k+1} = (\alpha \mathbf{I} - \sigma^k \mathbf{L}) M^{n+1,k} + M + \sigma^k \bar{B}^{n+1,k} + (1 + \alpha) \Delta t^k \bar{Q}^{n+1,k+1}. \quad (3.6) $$

The previous formula is nothing else than an arithmetico-geometric-like sequence having the form

$$ M^{n+1,k+1} = \mathbf{A}^{k,\alpha} M^{n+1,k} + C^{n,k} = \mathbf{A}^{k,\alpha} \dots \mathbf{A}^{0,\alpha} M + \sum_{j=0}^{k} \left( \prod_{l=0}^{k-j-1} \mathbf{A}^{k-l,\alpha} \right) C^{n,j}, $$

with

$$ \mathbf{A}^{k,\alpha} = \frac{\alpha \mathbf{I} - \sigma^k \mathbf{L}}{1 + \alpha}, \qquad C^{n,k} = \frac{M + \sigma^k \bar{B}^{n+1,k}}{1 + \alpha} + \Delta t^k \bar{Q}^{n+1,k+1}, \quad (3.7) $$

whose convergence is ensured as soon as the sequence

$$ \sum_{j=0}^{k} \left( \prod_{l=0}^{k-j-1} \mathbf{A}^{k-l, \alpha} \right) C^{n,j}, $$

converges when $k$ grows. Since the collision terms $\bar{Q}^{n+1,k+1}$ are assumed to be bounded independently of $k$, we have $|C^{n,k}|_{\infty} \le \bar{C}$ and the CFL condition ensures the existence of $0 \le \beta < 1$ such that

$$ \| \mathbf{A}^{k, \alpha} \|_{\infty, \infty} = \frac{\alpha + 2\sigma^k |\xi|}{1 + \alpha} \le \frac{\alpha + \beta}{1 + \alpha} < 1, $$

where $\|\cdot\|_{\infty,\infty}$ is the maximum subordinate matrix norm associated to $|\cdot|_{\infty}$. This allows to write

$$ \sum_{j=0}^{k} \left| \left( \prod_{l=0}^{k-j-1} \mathbf{A}^{k-l, \alpha} \right) C^{n,j} \right|_{\infty} \le \sum_{j=0}^{k} \left( \frac{\alpha + \beta}{1 + \alpha} \right)^{k-j-1} |C^{n,j}|_{\infty} \le \frac{\bar{C}}{1 - \frac{\alpha + \beta}{1 + \alpha}} \quad (3.8) $$

giving the absolute convergence of the series and thus proving the result. $\square$

Finally we investigate the entropy stability of the iterative kinetic scheme (3.2). Note that when restricting to a flat bathymetry, the explicit kinetic scheme [6] already verifies a discrete entropy inequality, while being more efficient since it avoids the use of sub-iterations. The benefit of the iterative kinetic scheme compared to the fully explicit one will become clearer in the next section, where it is combined with the hydrostatic reconstruction. Nevertheless, the simpler case of a flat bottom is interesting in order to understand why the iterative kinetic scheme involving the half-disk Maxwellian (1.4) offers a favorable setting to get a discrete entropy inequality.
where $\Xi = \text{supp } M^{n+1,k}$, where we recall that the function $\Psi$ defined in (2.6) is positive on $\mathbb{R}_+^2$, and with $i \pm 1 = i - \text{sgn } \xi$. As a consequence, if for any integer $k$ the CFL condition

$$ (3.10) \qquad \forall \xi \in \text{supp } M^{n+1,k}, \sigma^k |\xi| \le \alpha $$

holds, then $D_i^{n+1,k+1}$ is a dissipation term with negative sign and at each iteration the kinetic entropy is dissipated up to terms that are macroscopically zero, that is to say there exists a kinetic entropy flux $\tilde{H}_{0,i+1/2}^{n+1,k}$, a negative dissipation $\tilde{D}_i^{n+1,k+1}$ and a term $\tilde{Z}_i^{n+1,k+1}(\xi)$ whose integral over $\xi \in \mathbb{R}$ is zero such that

$$ (3.11) \qquad H_0(M_i^{n+1,k+1}, \xi) \le H_0(M_i, \xi) - \sigma^k \xi (\tilde{H}_{0,i+1/2}^{n+1,k} - \tilde{H}_{0,i-1/2}^{n+1,k}) + \tilde{D}_i^{n+1,k+1} + \tilde{Z}_i^{n+1,k+1}. $$

Before giving the proof we have the remark below.

**Remark 3.7** Even when the CFL condition (3.10) is not satisfied, we can ensure that the scheme (3.2) satisfies a discrete entropy inequality from some rank $k$ assuming the convergence of the method, which holds under the assumptions of Proposition 3.4. In fact, multiplying inequality (3.9) by $1 + \alpha$ it is possible to write

$$ (3.12) \qquad H_0(M_i^{n+1,k+1}, \xi) \le H_0(M_i, \xi) - \sigma^k \xi (H_{0,i+1/2}^{n+1,k} - H_{0,i-1/2}^{n+1,k}) + (1+\alpha)\eta'(U_i^{n+1,k+1}) \cdot \begin{pmatrix} 1 \\ \xi \end{pmatrix} Q_i^{n+1,k+1} \\ - \Psi(M_i^{n+1,k+1}, M_i) - \sigma^k |\xi| \mathbb{1}_{\Xi} \Psi(M_i^{n+1,k+1}, M_{i\pm 1}^{n+1,k}) \\ + \alpha (H_0(M_i^{n+1,k}, \xi) - H_0(M_i^{n+1,k+1}, \xi)) - (\alpha - \sigma^k |\xi| \mathbb{1}_{\Xi}) \Psi(M_i^{n+1,k+1}, M_i^{n+1,k}). $$

In the right hand side of (3.12), the quantity

$$ (1 + \alpha)\eta'(U_i^{n+1,k+1}) \cdot \begin{pmatrix} 1 \\ \xi \end{pmatrix} Q_i^{n+1,k+1} $$

does not cause any issue as it vanishes upon integration over $\xi \in \mathbb{R}$. This is due to the collision term $Q_i^{n+1,k+1}$ satisfying the conservation constraints (1.7), meaning that its integral against $(1, \xi)^T$ vanishes. Therefore in (3.12) the only problematic terms are contained in the last line, as their sign can be positive since we do not assume $\sigma^k |\xi| \mathbb{1}_{\Xi} \le \alpha$ anymore. Nevertheless, by regularity of $H_0(\cdot, \xi)$ and by definition (2.6) of $\Psi$, these terms write as a $\mathcal{O}(M_i^{n+1,k+1} - M_i^{n+1,k})$ and vanish as $k \to \infty$ owing to the convergence of the method. As a consequence, from some rank $k$ these two terms become negligible compared to $-\Psi(M_i^{n+1,k+1}, M_i) < 0$ which remains bounded away from zero, and we recover a dissipation with negative sign. (Note that if $M_i^{n+1,k}$ was converging to $M$ as $k \to \infty$, it would imply that $M$ solves the fixed point problem and thus $M_i^{n+1,k} = M$ for all $k$; putting aside this trivial case, this is why $\Psi(M_i^{n+1,k+1}, M_i)$ remains bounded away from zero).

**Proof of prop. 3.6.** First, we remark that for any $\xi$ there holds

$$ (3.13) \qquad \partial_f H_0(M_i^{n+1,k+1}, \xi) Q_i^{n+1,k+1} \le \eta'(U_i^{n+1,k+1}) \cdot \begin{pmatrix} 1 \\ \xi \end{pmatrix} Q_i^{n+1,k+1}, $$

which is implied by the relation (1.15) and by the fact that $Q_i^{n+1,k+1} = M_i^{n+1,k+1} - f_i^{n+1,k+1}$ is negative for any $\xi \notin \text{supp } M_i^{n+1,k+1}$. As a consequence of (3.13), inequality (3.9) (and equivalently (3.12)) is verified if there holds

$$ (3.14) \qquad H_0(M_i^{n+1,k+1}) = \frac{H_0(M_i) + \alpha H_0(M_i^{n+1,k})}{1 + \alpha} - \frac{\sigma^k \xi}{1 + \alpha} (H_{0,i+1/2}^{n+1,k} - H_{0,i-1/2}^{n+1,k}) \\ + \partial_f H_0(M_i^{n+1,k+1}) Q_i^{n+1,k+1} + D_i^{n+1,k+1}, $$
**Proposition 3.5.** The kinetic entropy of the iterative scheme (3.2) with half-disk Maxwellian (1.4) satisfies the following inequality

$$H_0(M_i^{n+1,k+1}) \le \frac{H_0(M_i) + \alpha H_0(M_i^{n+1,k})}{1+\alpha} - \frac{\sigma^k \xi}{1+\alpha} (H_{0,i+1/2}^{n+1,k} - H_{0,i-1/2}^{n+1,k}) + \eta'(U_i^{n+1,k+1}) \cdot \frac{1}{\xi} Q_i^{n+1,k+1} + D_i^{n+1,k+1}, \quad (3.9)$$

with $Q_i^{n+1,k+1} = M_i^{n+1,k+1} - f_i^{n+1,k+1}$ a collision operator verifying the conservation constraints (1.7), and with $\eta$ the entropy given in (1.3). The interfacial kinetic entropies $H_{0,i\pm 1/2}^{n+1,k}$ are

$$H_{0,i-1/2}^{n+1,k} = \mathbb{1}_{\xi>0} H_0(M_{i-1}^{n+1,k}, \xi) + \mathbb{1}_{\xi<0} H_0(M_i^{n+1,k}, \xi),$$

$$H_{0,i+1/2}^{n+1,k} = \mathbb{1}_{\xi>0} H_0(M_i^{n+1,k}, \xi) + \mathbb{1}_{\xi<0} H_0(M_{i+1}^{n+1,k}, \xi),$$

and the term $D_i^{n+1,k+1}$ is given by

$$D_i^{n+1,k+1} = -\frac{1}{1+\alpha}\Psi(M_i^{n+1,k+1}, M_i) - \frac{\alpha - \sigma^k |\xi| \mathbb{1}_{\Xi}}{1+\alpha} \Psi(M_i^{n+1,k+1}, M_i^{n+1,k}) \\ - \frac{\sigma^k |\xi| \mathbb{1}_{\Xi}}{1+\alpha} \Psi(M_i^{n+1,k+1}, M_{i\pm 1}^{n+1,k}),$$

where $\Xi = \text{supp } M^{n+1,k}$, where we recall that the function $\Psi$ defined in (2.6) is positive on $\mathbb{R}_+^2$, and with $i \pm 1 = i - \text{sign } \xi$. As a consequence, if for any integer $k$ the CFL condition

$$\forall \xi \in \text{supp } M^{n+1,k}, \sigma^k |\xi| \le \alpha \quad (3.10)$$

holds, then $D_i^{n+1,k+1}$ is a dissipation term with negative sign and at each iteration the kinetic entropy is dissipated up to terms that are macroscopically zero, that is to say there exists a kinetic entropy flux $\tilde{H}_{0,i+1/2}^{n+1,k}$, a negative dissipation $\tilde{D}_i^{n+1,k+1}$ and a term $\tilde{Z}_i^{n+1,k+1}(\xi)$ whose integral over $\xi \in \mathbb{R}$ is zero such that

$$H_0(M_i^{n+1,k+1}, \xi) \le H_0(M_i, \xi) - \sigma^k \xi (\tilde{H}_{0,i+1/2}^{n+1,k} - \tilde{H}_{0,i-1/2}^{n+1,k}) + \tilde{D}_i^{n+1,k+1} + \tilde{Z}_i^{n+1,k+1}. \quad (3.11)$$

Before giving the proof we have the remark below.

**Remark 3.6.** Even when the CFL condition (3.10) is not satisfied, we can ensure that the scheme (3.2) satisfies a discrete entropy inequality from some rank $k$ assuming the convergence of the method, which holds under the assumptions of Proposition 3.4. In fact, multiplying inequality (3.9) by $1 + \alpha$ it is possible to write

$$H_0(M_i^{n+1,k+1}, \xi) \le H_0(M_i, \xi) - \sigma^k \xi (H_{0,i+1/2}^{n+1,k} - H_{0,i-1/2}^{n+1,k}) + (1+\alpha) \eta'(U_i^{n+1,k+1}) \cdot \frac{1}{\xi} Q_i^{n+1,k+1} \\ - \Psi(M_i^{n+1,k+1}, M_i) - \sigma^k |\xi| \mathbb{1}_{\Xi} \Psi(M_i^{n+1,k+1}, M_{i\pm 1}^{n+1,k}) \\ + \alpha (H_0(M_i^{n+1,k}, \xi) - H_0(M_i^{n+1,k+1}, \xi)) - (\alpha - \sigma^k |\xi| \mathbb{1}_{\Xi}) \Psi(M_i^{n+1,k+1}, M_i^{n+1,k}). \quad (3.12)$$

In the right hand side of (3.12), the quantity

$$(1+\alpha)\eta'(U_i^{n+1,k+1}) \cdot \frac{1}{\xi} Q_i^{n+1,k+1}$$
To prove equality (3.14) we write the sub-iteration (3.2) as

$$ (3.15) \qquad M_i^{n+1,k+1} = \frac{1}{1+\alpha} \left( M_i + (\alpha - \sigma^k |\xi| \mathbf{1}_{\Xi}) M_i^{n+1,k} + \sigma^k |\xi| \mathbf{1}_{\Xi} M_{i\pm 1}^{n+1,k} \right) + Q_i^{n+1,k+1}, $$

with $i \pm 1 = i$ - sign $\xi$. Applying Lemma 2.4 for $a = M_i^{n+1,k+1}$ and $b = M_{i\pm 1}^{n+1,k}$, $M_i$, $M_i^{n+1,k}$, we respectively get:

$$ (3.16) \qquad H_0(M_{i\pm 1}^{n+1,k}) = H_0(M_i^{n+1,k+1}) + \partial_f H_0(M_i^{n+1,k+1})(M_{i\pm 1}^{n+1,k} - M_i^{n+1,k+1}) \\ \qquad \qquad + \Psi(M_i^{n+1,k+1}, M_{i\pm 1}^{n+1,k}), $$

$$ (3.17) \qquad H_0(M_i^{n+1,k}) = H_0(M_i^{n+1,k+1}) + \partial_f H_0(M_i^{n+1,k+1})(M_i^{n+1,k} - M_i^{n+1,k+1}) \\ \qquad \qquad + \Psi(M_i^{n+1,k+1}, M_i^{n+1,k}), $$

$$ (3.18) \qquad H_0(M_i) = H_0(M_i^{n+1,k+1}) + \partial_f H_0(M_i^{n+1,k+1})(M_i - M_i^{n+1,k+1}) \\ \qquad \qquad + \Psi(M_i^{n+1,k+1}, M_i). $$

Performing the linear combination

$$ \frac{1}{1+\alpha} \left( (3.18) + (\alpha - \sigma^k |\xi| \mathbf{1}_{\Xi})(3.17) + \sigma^k |\xi| \mathbf{1}_{\Xi}(3.16) \right), $$

and using (3.15) we obtain

$$ \frac{1}{1+\alpha} \left( H_0(M_i) + (\alpha - \sigma^k |\xi| \mathbf{1}_{\Xi}) H_0(M_i^{n+1,k}) + \sigma^k |\xi| \mathbf{1}_{\Xi} H_0(M_{i\pm 1}^{n+1,k}) \right) = \\ H_0(M_i^{n+1,k+1}) - \partial_f H_0(M_i^{n+1,k+1}) Q_i^{n+1,k+1} + \frac{1}{1+\alpha} \Psi(M_i^{n+1,k+1}, M_i) \\ + \frac{\alpha - \sigma^k |\xi| \mathbf{1}_{\Xi}}{1+\alpha} \Psi(M_i^{n+1,k+1}, M_i^{n+1,k+1}) + \frac{\sigma^k |\xi| \mathbf{1}_{\Xi}}{1+\alpha} \Psi(M_i^{n+1,k+1}, M_{i\pm 1}^{n+1,k}), $$

which corresponds to equality (3.14) after rearranging the terms.

Next we proceed by induction to show that the kinetic entropy is dissipated at every iteration assuming the CFL condition (3.10) holds for any integer $k$. The key argument is that under this CFL condition, the term $D_i^{n+1,k+1}$ defines a convex combination of negative quantities, and is thus negative. The initialization is obvious since we have $M_i^{n+1,0} = M_i$, so we focus on the recurrence. We want to show that (3.11) holds at some rank $k \ge 1$ assuming that it is satisfied at rank $k-1$. Under this assumption we can develop (3.9) as

$$ H_0(M_i^{n+1,k+1}) \le \frac{1}{1+\alpha} \left( H_0(M_i) + \alpha \left( H_0(M_i) - \sigma^k \xi \left( \tilde{H}_{0,i+1/2}^{n+1,k-1} - \tilde{H}_{0,i-1/2}^{n+1,k-1} \right) + \tilde{D}_i^{n+1,k} + \tilde{Z}_i^{n+1,k} \right) \right) \\ - \frac{\sigma^k \xi}{1+\alpha} \left( H_{0,i+1/2}^{n+1,k} - H_{0,i-1/2}^{n+1,k} \right) + \eta'(U_i^{n+1,k+1}) \cdot \left( \frac{1}{\xi} \right) Q_i^{n+1,k+1} + D_i^{n+1,k+1}, $$

with $\tilde{Z}_i^{n+1,k}$ and $\eta'(U_i^{n+1,k+1}) \cdot \left( \frac{1}{\xi} \right) Q_i^{n+1,k+1}$ macroscopically zero as per Remark 3.7. Therefore we have

$$ H_0(M_i^{n+1,k+1}) \le H_0(M_i) - \frac{\sigma^k \xi}{1+\alpha} \left( H_{0,i+1/2}^{n+1,k} + \alpha \tilde{H}_{0,i+1/2}^{n+1,k-1} - H_{0,i-1/2}^{n+1,k} - \alpha \tilde{H}_{0,i-1/2}^{n+1,k-1} \right) \\ + \frac{\alpha}{1+\alpha} \tilde{Z}_i^{n+1,k} + \eta'(U_i^{n+1,k+1}) \cdot \left( \frac{1}{\xi} \right) Q_i^{n+1,k+1} + \frac{\alpha}{1+\alpha} \tilde{D}_i^{n+1,k} + D_i^{n+1,k+1} $$
does not cause any issue as it vanishes upon integration over $\xi \in \mathbb{R}$. This is due to the collision term $Q_i^{n+1,k+1}$ satisfying the conservation constraints (1.7), meaning that its integral against $(1, \xi)^T$ vanishes. Therefore in (3.12) the only problematic terms are contained in the last line, as their sign can be positive since we do not assume $\sigma^k|\xi|1_\Xi \le \alpha$ anymore. Nevertheless, by regularity of $H_0(\cdot, \xi)$ and by definition (2.6) of $\Psi$, these terms write as a $O(M_i^{n+1,k+1} - M_i^{n+1,k})$ and vanish as $k \to \infty$ owing to the convergence of the method. As a consequence, from some rank $k$ these two terms become negligible compared to $-\Psi(M_i^{n+1,k+1}, M_i) < 0$ which remains bounded away from zero, and we recover a dissipation with negative sign. (Note that if $M^{n+1,k}$ was converging to $M$ as $k \to \infty$, it would imply that $M$ solves the fixed point problem and thus $M^{n+1,k} = M$ for all $k$; putting aside this trivial case, this is why $\Psi(M_i^{n+1,k+1}, M_i)$ remains bounded away from zero).

**Proof of Proposition 3.5.** First, we remark that for any $\xi$ there holds

$$ \partial_f H_0(M_i^{n+1,k+1}, \xi) Q_i^{n+1,k+1} \le \eta'(U_i^{n+1,k+1}) \cdot \frac{1}{\xi} Q_i^{n+1,k+1}, \quad (3.13) $$

which is implied by the relation (1.15) and by the fact that $Q_i^{n+1,k+1} = M_i^{n+1,k+1} - f_i^{n+1,k+1}$ is negative for any $\xi \notin \text{supp } M_i^{n+1,k+1}$. As a consequence of (3.13), inequality (3.9) (and equivalently (3.12)) is verified if there holds

$$ H_0(M_i^{n+1,k+1}) = \frac{H_0(M_i) + \alpha H_0(M_i^{n+1,k})}{1+\alpha} - \frac{\sigma^k \xi}{1+\alpha} (H_{0,i+1/2}^{n+1,k} - H_{0,i-1/2}^{n+1,k}) \\ + \partial_f H_0(M_i^{n+1,k+1}) Q_i^{n+1,k+1} + D_i^{n+1,k+1}. \quad (3.14) $$

To prove equality (3.14) we write the sub-iteration (3.2) as

$$ M_i^{n+1,k+1} = \frac{1}{1+\alpha} \left( M_i + \left( \alpha - \sigma^k |\xi| 1_\Xi \right) M_i^{n+1,k} + \sigma^k |\xi| 1_\Xi M_{i\pm 1}^{n+1,k} \right) + Q_i^{n+1,k+1}, \quad (3.15) $$

with $i \pm 1 = i$ — sign $\xi$. Applying Lemma 2.4 for $a = M_i^{n+1,k+1}$ and $b = M_{i\pm 1}^{n+1,k}$, $M_i, M_i^{n+1,k}$, we respectively get:

$$ H_0(M_{i\pm 1}^{n+1,k}) = H_0(M_i^{n+1,k+1}) + \partial_f H_0(M_i^{n+1,k+1})(M_{i\pm 1}^{n+1,k} - M_i^{n+1,k+1}) + \Psi(M_i^{n+1,k+1}, M_{i\pm 1}^{n+1,k}), \quad (3.16) $$

$$ H_0(M_i^{n+1,k}) = H_0(M_i^{n+1,k+1}) + \partial_f H_0(M_i^{n+1,k+1})(M_i^{n+1,k} - M_i^{n+1,k+1}) + \Psi(M_i^{n+1,k+1}, M_i^{n+1,k}), \quad (3.17) $$

$$ H_0(M_i) = H_0(M_i^{n+1,k+1}) + \partial_f H_0(M_i^{n+1,k+1})(M_i - M_i^{n+1,k+1}) + \Psi(M_i^{n+1,k+1}, M_i). \quad (3.18) $$

Performing the linear combination

$$ \frac{1}{1+\alpha} \left( (3.18) + \left( \alpha - \sigma^k |\xi| 1_\Xi \right) (3.17) + \sigma^k |\xi| 1_\Xi (3.16) \right), $$

and using (3.15) we obtain

$$ \begin{aligned} & \frac{1}{1+\alpha} \left( H_0(M_i) + \left( \alpha - \sigma^k |\xi| 1_\Xi \right) H_0(M_i^{n+1,k}) + \sigma^k |\xi| 1_\Xi H_0(M_{i\pm 1}^{n+1,k}) \right) \\ & = H_0(M_i^{n+1,k+1}) - \partial_f H_0(M_i^{n+1,k+1}) Q_i^{n+1,k+1} + \frac{1}{1+\alpha} \Psi(M_i^{n+1,k+1}, M_i) \\ & \quad + \frac{\alpha - \sigma^k |\xi| 1_\Xi}{1+\alpha} \Psi(M_i^{n+1,k+1}, M_i^{n+1,k}) + \frac{\sigma^k |\xi| 1_\Xi}{1+\alpha} \Psi(M_i^{n+1,k+1}, M_{i\pm 1}^{n+1,k}), \end{aligned} $$

which corresponds to equality (3.14) after rearranging the terms.
and the proof is complete by setting

$$ \tilde{H}_{0,i+1/2}^{n+1,k} = \frac{1}{1+\alpha} (H_{0,i+1/2}^{n+1,k} + \alpha \tilde{H}_{0,i+1/2}^{n+1,k-1}), \quad \tilde{D}_{i}^{n+1,k+1} = \frac{\alpha}{1+\alpha} \tilde{D}_{i}^{n+1,k} + D_{i}^{n+1,k+1}, \\ \tilde{Z}_{i}^{n+1,k+1} = \frac{\alpha}{1+\alpha} \tilde{Z}_{i}^{n+1,k} + \eta'(U_{i}^{n+1,k+1}) \cdot \left(\frac{1}{\xi}\right) Q_{i}^{n+1,k+1}. $$

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

## 3.2 Case of a non flat topography

In presence of a varying bathymetry, one difficulty is to design well-balanced schemes, that is to say schemes that preserve the lake at rest equilibrium given by $(h, q) = (-z, 0)$, where we recall $z(x)$ is a parametrization of the topography at the bottom. The hydrostatic reconstruction technique introduced in [3] is a well-known strategy to circumvent this difficulty. More recently [4], it has been established that when combined with the explicit kinetic scheme, this reconstruction step can prevent the entropy from being dissipated in some testcases featuring a non-flat bathymetry, no matter how refined the time and spatial steps are. This is due to a positive error term appearing in the discrete entropy inequality and which can dominate the dissipation. In this section our goal is to show that an iterative version of this explicit kinetic scheme does not suffer from this defect provided enough sub-iterations are performed. Therefore the iterative strategy, which approximates a fully implicit kinetic scheme, can be considered an improvement over the explicit method in terms of stability. We also emphasize that compared to the fully implicit kinetic scheme proposed in Section 2, the better stability properties of the iterative scheme does not translate in the possibility to use arbitrary large time steps. In fact a CFL condition is required to ensure the convergence of the sub-iterations, as already pointed by the Propositions 3.1 and 3.4 for iterative schemes applied to the flat bottom case.

Let us briefly recall the principle of the hydrostatic reconstruction, which is based on the reconstruction of the water height at every interface. Let $U_i = (h_i, h_i u_i)^T \in \mathbb{R}^2$ denote the vector of quantities of interest over cell $1 \le i \le P$, with $P$ the number of interior cells and with ghost cells corresponding to indices 0 and $P+1$. The reconstructed states are vectors from $\mathbb{R}^+ \times \mathbb{R}$ defined on the left and right neighborhoods of each cell interface as follows:

$$ (3.19) \qquad \forall 1 \le i \le P, \quad U_{i+1/2-} = \begin{pmatrix} h_{i+1/2-} \\ h_{i+1/2-} - u_i \end{pmatrix}, \quad U_{i-1/2+} = \begin{pmatrix} h_{i-1/2+} \\ h_{i-1/2+} - u_i \end{pmatrix}. $$

The reconstructed interfacial water heights are given by

$$ (3.20) \qquad h_{i-1/2+} = (h_i + z_i - z_{i-1/2})_+, \quad h_{i+1/2-} = (h_i + z_i - z_{i+1/2})_+ $$

with the interfacial bathymetry $z_{i+1/2} = \max(z_i, z_{i+1})$. The truly implicit kinetic scheme we are considering reads as below

$$ (3.21) \qquad U_i^{n+1} = U_i^n - \sigma(F_{i+1/2-}^{n+1} - F_{i-1/2+}^{n+1}), $$

with $\sigma = \Delta t/\Delta x$ and the numerical fluxes decomposed as:

$$ (3.22) \qquad \left\{ \begin{array}{l} F_{i+1/2-}^{n+1} = \mathcal{F}(U_{i+1/2-}^{n+1}, U_{i+1/2+}^{n+1}) + \frac{g}{2} \begin{pmatrix} 0 \\ (h_i^{n+1})^2 - (h_{i+1/2-}^{n+1})^2 \end{pmatrix} \\ F_{i-1/2+}^{n+1} = \mathcal{F}(U_{i-1/2-}^{n+1}, U_{i-1/2+}^{n+1}) + \frac{g}{2} \begin{pmatrix} 0 \\ (h_i^{n+1})^2 - (h_{i-1/2+}^{n+1})^2 \end{pmatrix} \end{array} \right. . $$

We recall that in our case the upwinding of the numerical flux $\mathcal{F}$ is induced at the kinetic level according to definition (3.4). In (3.22), the terms in factor of $g/2$ look like pressure variations and are consistent with the topography source term when we subtract them. In fact, owing to (3.20) we formally have

$$ \frac{g}{2} \left( (h_{i+1/2-}^{n+1})^2 - (h_{i-1/2+}^{n+1})^2 \right) = g \frac{h_{i+1/2-}^{n+1} + h_{i-1/2+}^{n+1}}{2} (h_{i+1/2-}^{n+1} - h_{i-1/2+}^{n+1}) \approx g h_i^{n+1} (z_{i+1/2} - z_{i-1/2}). $$
Next we proceed by induction to show that the kinetic entropy is dissipated at every iteration assuming the CFL condition (3.10) holds for any integer $k$. The key argument is that under this CFL condition, the term $D_i^{n+1,k+1}$ defines a convex combination of negative quantities, and is thus negative. The initialization is obvious since we have $M_i^{n+1,0} = M_i$, so we focus on the recurrence. We want to show that (3.11) holds at some rank $k \ge 1$ assuming that it is satisfied at rank $k-1$. Under this assumption we can develop (3.9) as

$$H_0(M_i^{n+1,k+1}) \le \frac{1}{1+\alpha} \left( H_0(M_i) + \alpha \left( H_0(M_i) - \sigma^k \xi \left( \tilde{H}_{0,i+1/2}^{n+1,k-1} - \tilde{H}_{0,i-1/2}^{n+1,k-1} \right) + \tilde{D}_i^{n+1,k} + \tilde{Z}_i^{n+1,k} \right) \right) \\ - \frac{\sigma^k \xi}{1+\alpha} \left( H_{0,i+1/2}^{n+1,k} - H_{0,i-1/2}^{n+1,k} \right) + \eta' \left( U_i^{n+1,k+1} \right) \cdot \left( \frac{1}{\xi} \right) Q_i^{n+1,k+1} + D_i^{n+1,k+1},$$

with $\tilde{Z}_i^{n+1,k}$ and $\eta'(U_i^{n+1,k+1}) \cdot \left(\frac{1}{\xi}\right) Q_i^{n+1,k+1}$ macroscopically zero as per Remark 3.6. Therefore we have

$$H_0(M_i^{n+1,k+1}) \le H_0(M_i) - \frac{\sigma^k \xi}{1+\alpha} \left( H_{0,i+1/2}^{n+1,k} + \alpha \tilde{H}_{0,i+1/2}^{n+1,k-1} - H_{0,i-1/2}^{n+1,k} - \alpha \tilde{H}_{0,i-1/2}^{n+1,k-1} \right) \\ + \frac{\alpha}{1+\alpha} \tilde{Z}_i^{n+1,k} + \eta'(U_i^{n+1,k+1}) \cdot \left(\frac{1}{\xi}\right) Q_i^{n+1,k+1} + \frac{\alpha}{1+\alpha} \tilde{D}_i^{n+1,k} + D_i^{n+1,k+1}$$

and the proof is complete by setting

$$\tilde{H}_{0,i+1/2}^{n+1,k} = \frac{1}{1+\alpha} \left( H_{0,i+1/2}^{n+1,k} + \alpha \tilde{H}_{0,i+1/2}^{n+1,k-1} \right), \quad \tilde{D}_i^{n+1,k+1} = \frac{\alpha}{1+\alpha} \tilde{D}_i^{n+1,k} + D_i^{n+1,k+1}, \\ \tilde{Z}_i^{n+1,k+1} = \frac{\alpha}{1+\alpha} \tilde{Z}_i^{n+1,k} + \eta'(U_i^{n+1,k+1}) \cdot \left(\frac{1}{\xi}\right) Q_i^{n+1,k+1}.$$

□

## 3.2. Case of a non flat topography

In presence of a varying bathymetry, one difficulty is to design well-balanced schemes, that is to say schemes that preserve the lake at rest equilibrium given by $(h, q) = (-z, 0)$, where we recall $z(x)$ is a parameterization of the topography at the bottom. One of the issue stems from the numerical diffusion which, when the bathymetry is non flat, can induce a transfer of water from one cell to its neighbors even in the case of the rest state. To illustrate this we can consider the toy update below

$$\frac{h_i^{n+1} - h_i^n}{\Delta t} = - \frac{q_{i+1/2}^n - q_{i-1/2}^n}{\Delta x} + \nu \frac{h_{i+1}^n - 2h_i^n + h_{i-1}^n}{\Delta x} = - \nu \frac{z_{i+1} - 2z_i + z_{i-1}}{\Delta x},$$

where $\nu > 0$ is a numerical viscosity coefficient, and where the second equality holds if we assume a discrete lake at rest such that $q_{i+1/2}^n = 0$ and $h_i^n + z_i = \text{Cst}$ for all $i$. For general bathymetries the right hand side doesn't vanish, and thus the water height fails to remain constant from one iteration to the other. A second problem lies in the momentum equation, which in the case of a lake at rest implies the balancing between the spatial derivative of the pressure and the source term:

$$\partial_x \left( \frac{g}{2} h^2 \right) = -gh \partial_x z.$$

One has to recover this equality at the discrete level by approximating the source term accordingly with the choice of numerical flux.

The hydrostatic reconstruction technique introduced in [5] is a well-known strategy to circumvent these difficulties. More recently [6], it has been established that when combined with the explicit kinetic scheme, this reconstruction step can prevent the entropy from being dissipated in some test cases featuring a non flat
This discretization of the topography source term can be interpreted at the kinetic level through the relations

$$ \int_{\mathbb{R}} \begin{pmatrix} 1 \\ \xi \end{pmatrix} (\xi - u_i)(M(U_i, \xi) - M(U_{i+1/2-}, \xi)) \mathrm{d}\xi = \begin{pmatrix} 0 \\ \frac{g}{2}(h_i^2 - h_{i+1/2-}^2) \end{pmatrix}, $$

$$ \int_{\mathbb{R}} \begin{pmatrix} 1 \\ \xi \end{pmatrix} (\xi - u_i)(M(U_i, \xi) - M(U_{i-1/2+}, \xi)) \mathrm{d}\xi = \begin{pmatrix} 0 \\ \frac{g}{2}(h_i^2 - h_{i-1/2+}^2) \end{pmatrix}, $$

and we introduce the following notations for convenience

$$ (3.23) \qquad \begin{cases} \delta M_{i+1/2-}(\xi) = (\xi - u_i)(M(U_i, \xi) - M(U_{i+1/2-}, \xi)), \\ \delta M_{i-1/2+}(\xi) = (\xi - u_i)(M(U_i, \xi) - M(U_{i-1/2+}, \xi)). \end{cases} $$

An issue with the nonlinear update (3.21) is that it cannot be solved analytically. Instead we will approximate it by an iterative process with a relaxation parameter $\alpha > 0$ similar to the one from Section 3.1, and which reads

$$ (3.24) \qquad \forall 1 \le i \le P, \quad (1+\alpha)U_i^{n+1,k+1} = U_i + \alpha U_i^{n+1,k} - \sigma^k (F_{i+1/2-}^{n+1,k} - F_{i-1/2+}^{n+1,k}), $$

where we can choose to take $U^{n+1,0} = U^n$ as the initialization. If the sequence defined through (3.24) converges, we recover the implicit scheme (3.21) by setting the macroscopic update as $U^{n+1} = \lim_{k \to \infty} U^{n+1,k}$ and $\sigma = \lim_{k \to \infty} \sigma^k$. In practice we will stop the sub-iterations for some $k$ large enough such that $U^{n+1,k} \approx U^{n+1}$. We also comment on the fact that when the bathymetry is flat ($z \equiv \text{Cst}$) the hydrostatic reconstruction becomes transparent, in the sense that the scheme (3.24) coincides with (3.3).

At the kinetic level, the recurrence relation (3.24) consists in introducing for any real $\xi$ the sequence $(f^{n+1,k}(\xi))_{k \in \mathbb{N}} \subset \mathbb{R}_+^P$ initialized with $f^{n+1,0}(\xi) = M(U^n, \xi)$ and defined recursively as:

$$ (3.25) \qquad \left\{ \begin{aligned} (1+\alpha)f_i^{n+1,k+1} &= M_i + \alpha M_i^{n+1,k} - \sigma^k \left( \xi(M_{i+1/2}^{n+1,k} - M_{i-1/2}^{n+1,k}) + \delta M_{i+1/2-}^{n+1,k} - \delta M_{i-1/2+}^{n+1,k} \right), \\ M_i^{n+1,k+1} &= f^{n+1,k+1} + \Delta t^k Q^{n+1,k+1} \end{aligned} \right. $$

where the interfacial Maxwellians $M_{i \pm 1/2}$ are defined through the following upwinding

$$ (3.26) \qquad M_{i \pm 1/2}(\xi) = \mathbb{1}_{\xi>0} M(U_{i \pm 1/2-}, \xi) + \mathbb{1}_{\xi<0} M(U_{i \pm 1/2+}, \xi) \quad \forall 1 \le i \le P, $$

and where we have used the notations (3.23).

We now consider the question of whether or not the iterative kinetic scheme with hydrostatic reconstruction (3.25)-(3.26) is structure preserving. We recall that in our context, structure preserving means the ability to keep the water height positive, to preserve the lake at rest steady state and to satisfy a discrete entropy inequality. When considering a flat bathymetry, this latter point was proved by making use of Lemma 2.4, providing in particular an estimate of the discrete spatial variation of the flux, namely

$$ \partial_f H(M_i)(M_i - M_{i \pm 1}) = H(M_i) - H(M_{i \pm 1}) + \Psi(M_i, M_{i \pm 1}), $$

which involves a conservative difference and the positive function $\Psi$ defined in (2.6) responsible for the dissipation. In presence of a varying bathymetry however, fluxes are computed after performing the hydrostatic reconstruction, and we now need an estimate for the following quantity

$$ \partial_f H(M_i)(\xi(M_{i+1/2} - M_{i-1/2}) + \delta M_{i+1/2-} - \delta M_{i-1/2+}), $$

which appears when multiplying the first line of (3.25) by $\partial_f H(M_i^{n+1,k}, \xi)$. Because the derivative of the kinetic entropy is multiplied by quantities all different from $M_i$, we can no longer use Lemma 2.4. Instead, we will need the following inequality that corresponds to the Proposition 3.1 from [4].
bathymetry, no matter how refined the time and spatial steps are. This is due to a positive error term appearing in the discrete entropy inequality and which can dominate the dissipation. In this section our goal is to show that an iterative version of this explicit kinetic scheme does not suffer from this defect provided enough sub-iterations are performed. Therefore the iterative strategy, which approximates a fully implicit kinetic scheme, can be considered an improvement over the explicit method in terms of stability. We also emphasize that compared to the fully implicit kinetic scheme proposed in Section 2, the stability properties of the iterative scheme does not translate in the possibility to use arbitrary large time steps. In fact a CFL condition is required to ensure the convergence of the sub-iterations, as already pointed by the Propositions 3.1 and 3.4 for iterative schemes applied to the flat bottom case.

Let us briefly recall the principle of the hydrostatic reconstruction, which is based on the reconstruction of the water height at every interface. Let $U_i = (h_i, h_i u_i)^T \in \mathbb{R}^2$ denote the vector of quantities of interest over cell $1 \le i \le P$, with $P$ the number of interior cells and with ghost cells corresponding to indices 0 and $P+1$. The reconstructed states are vectors from $\mathbb{R}^+ \times \mathbb{R}$ defined on the left and right neighborhoods of each cell interface as follows:

$$ \forall 1 \le i \le P, \quad U_{i+1/2-} = \begin{pmatrix} h_{i+1/2-} \\ h_{i+1/2-} - u_i \end{pmatrix}, \quad U_{i-1/2+} = \begin{pmatrix} h_{i-1/2+} \\ h_{i-1/2+} - u_i \end{pmatrix}. \qquad (3.19) $$

The reconstructed interfacial water heights are given by

$$ h_{i-1/2+} = (h_i + z_i - z_{i-1/2})_+, \quad h_{i+1/2-} = (h_i + z_i - z_{i+1/2})_+, \qquad (3.20) $$

with the interfacial bathymetry $z_{i+1/2} = \max(z_i, z_{i+1})$. The truly implicit kinetic scheme we are considering reads as below

$$ U_i^{n+1} = U_i^n - \sigma (F_{i+1/2-}^{n+1} - F_{i-1/2+}^{n+1}), \qquad (3.21) $$

with $\sigma = \Delta t/\Delta x$ and the numerical fluxes decomposed as:

$$ \left\{ \begin{aligned} F_{i+1/2-}^{n+1} &= \mathcal{F}(U_{i+1/2-}^{n+1}, U_{i+1/2+}^{n+1}) + \frac{g}{2} \begin{pmatrix} 0 \\ (h_i^{n+1})^2 - (h_{i+1/2-}^{n+1})^2 \end{pmatrix} \\ F_{i-1/2+}^{n+1} &= \mathcal{F}(U_{i-1/2-}^{n+1}, U_{i-1/2+}^{n+1}) + \frac{g}{2} \begin{pmatrix} 0 \\ (h_i^{n+1})^2 - (h_{i-1/2+}^{n+1})^2 \end{pmatrix}. \end{aligned} \right. \qquad (3.22) $$

We recall that in our case the upwinding of the numerical flux $\mathcal{F}$ is induced at the kinetic level according to definition (3.4). In (3.22), the terms in factor of $g/2$ look like pressure variations and are consistent with the topography source term when we subtract them. In fact, owing to (3.20) we formally have

$$ \frac{g}{2} \left( \left(h_{i+1/2-}^{n+1}\right)^2 - \left(h_{i-1/2+}^{n+1}\right)^2 \right) = g \frac{h_{i+1/2-}^{n+1} + h_{i-1/2+}^{n+1}}{2} \left(h_{i+1/2-}^{n+1} - h_{i-1/2+}^{n+1}\right) \approx g h_i^{n+1} (z_{i+1/2} - z_{i-1/2}). $$

This discretization of the topography source term can be interpreted at the kinetic level through the relations

$$ \int_{\mathbb{R}} \left( \frac{1}{\xi} \right) (\xi - u_i) (M(U_i, \xi) - M(U_{i+1/2-}, \xi)) \, d\xi = \begin{pmatrix} 0 \\ \frac{g}{2} (h_i^2 - h_{i+1/2-}^2) \end{pmatrix}, $$

$$ \int_{\mathbb{R}} \left( \frac{1}{\xi} \right) (\xi - u_i) (M(U_i, \xi) - M(U_{i-1/2+}, \xi)) \, d\xi = \begin{pmatrix} 0 \\ \frac{g}{2} (h_i^2 - h_{i-1/2+}^2) \end{pmatrix}, $$
**Proposition 3.8** (Audusse, Bouchut, Bristeau, and Sainte-Marie [4]). One has

$$ (3.27) \quad - \partial_f H(M_i^{n+1,k}, z_i) [\xi(M_{i+1/2}^{n+1,k} - M_{i-1/2}^{n+1,k}) + \delta M_{i+1/2-}^{n+1,k} - \delta M_{i-1/2+}^{n+1,k}] \le \tilde{G}_{i+1/2-}^{n+1,k} - \tilde{G}_{i-1/2+}^{n+1,k}, \\ \text{with } \tilde{G}_{i+1/2-}^{n+1,k} \text{ and } \tilde{G}_{i-1/2+}^{n+1,k} \text{ defined by} $$

$$ (3.28) \quad \tilde{G}_{i \pm 1/2 \mp} = \xi \mathbb{1}_{\xi < 0} H(M_{i \pm 1/2 \mp}, z_{i \pm 1/2}) + \xi \mathbb{1}_{\xi > 0} H(M_{i \pm 1/2 \mp}, z_{i \pm 1/2}) \\ \qquad + \xi H(M_i, z_i) - \xi H(M_{i \pm 1/2 \mp}, z_{i \pm 1/2}) + (\eta'(U_i) \cdot \begin{pmatrix} 1 \\ \xi \end{pmatrix} + gz_i)(\xi M_{i \pm 1/2 \mp} - \xi M_i + \delta M_{i \pm 1/2 \mp}). $$

Moreover, the difference $\tilde{G}_{i+1/2-}^{n+1,k} - \tilde{G}_{i-1/2+}^{n+1,k}$ is macroscopically conservative

$$ \int_{\mathbb{R}} (\tilde{G}_{i+1/2-}^{n+1,k} - \tilde{G}_{i-1/2+}^{n+1,k}) d\xi = \overline{G}_{i+1/2}^{n+1,k} - \overline{G}_{i-1/2}^{n+1,k}, $$

with the macroscopic entropy flux

$$ (3.29) \quad \overline{G}_{i+1/2} = \int_{\mathbb{R}} \xi (\mathbb{1}_{\xi>0} H(M_{i+1/2-}, z_{i+1/2}) + \mathbb{1}_{\xi<0} H(M_{i+1/2+}, z_{i+1/2})) d\xi. $$

We are now able to state our result.

**Proposition 3.9** The scheme (3.24) is structure preserving, in the sense that we have the three properties below:

(i) the sub-iterations are well balanced, that is to say $U^{n+1,k} = U^n$ for all $k \in \mathbb{N}$ if $U^n$ is a discrete lake at rest;

(ii) assuming that the water height vectors $h^n$ and $h^{n+1,k}$ are positive, the update $h^{n+1,k+1}$ defined in the iterative scheme (3.24) is also positive if for all $1 \le i \le P$ the CFL condition $\sigma^k|\xi| \le \alpha + M_i/M_i^{n+1,k}$ holds for any $\xi$ belonging to supp $M^{n+1,k}$;

(iii) the kinetic entropy of the iterative process (3.25) verifies the following kinetic entropy inequality

$$ (3.30) \quad H(M_i^{n+1,k+1}, z_i) \le H(M_i, z_i) - \sigma^k (\tilde{G}_{i+1/2-}^{n+1,k} - \tilde{G}_{i-1/2+}^{n+1,k}) + (1+\alpha)(\eta'(U_i^{n+1,k+1}) \cdot \begin{pmatrix} 1 \\ \xi \end{pmatrix} + gz_i)Q_i^{n+1,k+1} \\ \qquad + \alpha(H(M_i^{n+1,k}, z_i) - H(M_i^{n+1,k+1}, z_i)) \\ \qquad + (1+\alpha)(\Psi(M_i^{n+1,k}, f_i^{n+1,k+1}) - \Psi(M_i^{n+1,k+1}, f_i^{n+1,k+1})) - \Psi(M_i^{n+1,k}, M_i), $$

where $Q_i^{n+1,k+1} = M_i^{n+1,k+1} - f_i^{n+1,k+1}$ is a collision term verifying the conservation constraints (1.7), where $\Psi$ defined in (2.6) is positive and where $\tilde{G}_{i+1/2-}^{n+1,k}$ and $\tilde{G}_{i-1/2+}^{n+1,k}$ are defined by (3.28).

Before giving the proof we make the following remark.

**Remark 3.10** We reiterate the comments made in Remark 3.7 which are to say that in (3.30) the term

$$ (\eta'(U_i^{n+1,k+1}) \cdot \begin{pmatrix} 1 \\ \xi \end{pmatrix} + gz_i)Q_i^{n+1,k+1} $$

is macroscopically zero since $Q_i^{n+1,k+1}$ is a collision term satisfying the conservation constraints (1.7). Besides, assuming the method converges as $k \to \infty$, the quantity

$$ \alpha(H(M_i^{n+1,k}, z_i) - H(M_i^{n+1,k+1}, z_i)) + (1+\alpha)(\Psi(M_i^{n+1,k}, f_i^{n+1,k+1}) - \Psi(M_i^{n+1,k+1}, f_i^{n+1,k+1})) $$
and we introduce the following notations for convenience

$$ \left\{ \begin{aligned} \delta M_{i+1/2-}(\xi) &= (\xi - u_i)(M(U_i, \xi) - M(U_{i+1/2-}, \xi)), \\ \delta M_{i-1/2+}(\xi) &= (\xi - u_i)(M(U_i, \xi) - M(U_{i-1/2+}, \xi)). \end{aligned} \right. \qquad (3.23) $$

An issue with the nonlinear update (3.21) is that it cannot be solved analytically. Instead we will approximate it by an iterative process with a relaxation parameter $\alpha > 0$ similar to the one from Section 3.1, and which reads

$$ \forall 1 \le i \le P, \qquad (1 + \alpha) U_i^{n+1,k+1} = U_i + \alpha U_i^{n+1,k} - \sigma^k (F_{i+1/2-}^{n+1,k} - F_{i-1/2+}^{n+1,k}), \qquad (3.24) $$

where we can choose to take $U^{n+1,0} = U^n$ as the initialization. If the sequence defined through (3.24) converges, we recover the implicit scheme (3.21) by setting the macroscopic update as $U^{n+1} = \lim_{k \to \infty} U^{n+1,k}$ and $\sigma = \lim_{k \to \infty} \sigma^k$. In practice we will stop the sub-iterations for some $k$ large enough such that $U^{n+1,k} \approx U^{n+1}$. We also comment on the fact that when the bathymetry is flat ($z \equiv \text{Cst}$) the hydrostatic reconstruction becomes transparent, in the sense that the scheme (3.24) coincides with (3.3).

At the kinetic level, the recurrence relation (3.24) consists in introducing for any real $\xi$ the sequence $(f^{n+1,k}(\xi))_{k \in \mathbb{N}} \subset \mathbb{R}_+^P$ initialized with $f^{n+1,0}(\xi) = M(U^n, \xi)$ and defined recursively as:

$$ \left\{ \begin{aligned} (1 + \alpha) f_i^{n+1,k+1} &= M_i + \alpha M_i^{n+1,k} - \sigma^k \left( \xi (M_{i+1/2}^{n+1,k} - M_{i-1/2}^{n+1,k}) + \delta M_{i+1/2-}^{n+1,k} - \delta M_{i-1/2+}^{n+1,k} \right) \\ M_i^{n+1,k+1} &= f^{n+1,k+1} + \Delta t^k Q^{n+1,k+1}, \end{aligned} \right. \qquad (3.25) $$

where the interfacial Maxwellians $M_{i \pm 1/2}$ are defined through the following upwinding

$$ M_{i \pm 1/2}(\xi) = 1_{\xi>0} M(U_{i \pm 1/2-}, \xi) + 1_{\xi<0} M(U_{i \pm 1/2+}, \xi) \quad \forall 1 \le i \le P, \qquad (3.26) $$

and where we have used the notations (3.23).

We now consider the question of whether or not the iterative kinetic scheme with hydrostatic reconstruction (3.25) and (3.26) is structure preserving. We recall that in our context, structure preserving means the ability to keep the water height positive, to preserve the lake at rest steady state and to satisfy a discrete entropy inequality. When considering a flat bathymetry, this latter point was proved by making use of Lemma 2.4, providing in particular an estimate of the discrete spatial variation of the flux, namely

$$ \partial_f H(M_i)(M_i - M_{i \pm 1}) = H(M_i) - H(M_{i \pm 1}) + \Psi(M_i, M_{i \pm 1}), $$

which involves a conservative difference and the positive function $\Psi$ defined in (2.6) responsible for the dissipation. In presence of a varying bathymetry however, fluxes are computed after performing the hydrostatic reconstruction, and we now need an estimate for the following quantity

$$ \partial_f H(M_i)(\xi(M_{i+1/2} - M_{i-1/2}) + \delta M_{i+1/2-} - \delta M_{i-1/2+}), $$

which appears when multiplying the first line of (3.25) by $\partial_f H(M_i^{n+1,k}, \xi)$. Because the derivative of the kinetic entropy is multiplied by quantities all different from $M_i$, we can no longer use Lemma 2.4. Instead, we will need the following inequality that corresponds to the Proposition 3.1 from [6].

**Proposition 3.7 ([6]).** One has

$$ -\partial_f H(M_i^{n+1,k}, z_i) [\xi(M_{i+1/2}^{n+1,k} - M_{i-1/2}^{n+1,k}) + \delta M_{i+1/2-}^{n+1,k} - \delta M_{i-1/2+}^{n+1,k}] \le -\tilde{G}_{i+1/2-}^{n+1,k} + \tilde{G}_{i-1/2+}^{n+1,k}, \qquad (3.27) $$

with $\tilde{G}_{i+1/2-}^{n+1,k}$ and $\tilde{G}_{i-1/2+}^{n+1,k}$ defined by

$$ \tilde{G}_{i \pm 1/2 \mp} = \xi 1_{\xi<0} H(M_{i \pm 1/2+}, z_{i \pm 1/2}) + \xi 1_{\xi>0} H(M_{i \pm 1/2-}, z_{i \pm 1/2}) + \xi H(M_i, z_i) $$
will eventually become negligible compared to $-\Psi(M_i^{n+1,k}, M_i) < 0$ from some rank $k$. Integrating inequality (3.30) over $\xi \in \mathbb{R}$, this implies that there exists $K \in \mathbb{N}$ such that for any $k \ge K$ the fully discrete entropy inequality

$$ (3.31) \qquad \eta(U_i^{n+1,k+1}) \le \eta(U_i^n) - \sigma^k (\overline{G}_{i+1/2}^{n+1,k} - \overline{G}_{i-1/2}^{n+1,k}) $$

is satisfied at the macroscopic level, with $\overline{G}_{i+1/2}^{n+1,k}$ given in (3.29). Summing inequality (3.31) over every cell $1 \le i \le P$ we obtain the dissipation of the total energy up to boundary fluxes

$$ (3.32) \qquad \frac{1}{\Delta t^k} \sum_{i=1}^{P} \left( \eta(U_i^{n+1,k+1}) - \eta(U_i^n) \right) + \frac{1}{\Delta x} \left( \int_{\mathbb{R}} \xi H_{P+1/2}^{n+1,k}(\xi) \, d\xi - \int_{\mathbb{R}} \xi H_{1/2}^{n+1,k}(\xi) \, d\xi \right) \le 0. $$

In addition to the usual tolerance criterion where the iterations are stopped whenever two successive iterates are sufficiently close to each other, we can use (3.32) as a complementary condition to ensure the dissipation of total energy.

**Proof of prop. 3.9.** The proof makes use of the kinetic writing (3.25) of scheme (3.24).

(i) The well-balancedness is a consequence of the hydrostatic reconstruction, see [3].

(ii) Remarking that the quantity $\delta M_{i+1/2-}^{n+1,k} - \delta M_{i-1/2+}^{n+1,k}$ appearing in the last line of (3.25) defines an odd function of $\xi - u_i^{n+1,k}$, its integral over $\xi \in \mathbb{R}$ vanishes and we have at the macroscopic level

$$ (1 + \alpha)h_i^{n+1,k+1} = \int_{\mathbb{R}} \left( M_i + \alpha M_i^{n+1,k} - \sigma^k \xi (M_{i+1/2}^{n+1,k} - M_{i-1/2}^{n+1,k}) \right) d\xi. $$

Thus it is enough to prove the positivity of the integrand, whose developed form is

$$ M_i + \alpha M_i^{n+1,k} - \sigma^k \xi (1_{\xi>0} M_{i+1/2-}^{n+1,k} - 1_{\xi<0} M_{i-1/2+}^{n+1,k}) + \sigma^k \xi (1_{\xi>0} M_{i-1/2-}^{n+1,k} - 1_{\xi<0} M_{i+1/2+}^{n+1,k}). $$

By definition of the water height reconstruction (3.20), we have the inequalities $h_{i+1/2-}^{n+1,k} \le h_i^{n+1,k}$ and $h_{i-1/2+}^{n+1,k} \le h_i^{n+1,k}$. As a consequence $M_{i+1/2-}^{n+1,k} \le M_i^{n+1,k}$ and $M_{i-1/2+}^{n+1,k} \le M_i^{n+1,k}$, which allows us to bound the integrand from below by

$$ M_i + \alpha M_i^{n+1,k} - \sigma^k |\xi| M_i^{n+1,k}. $$

If $\xi$ does not belong to supp $M_i^{n+1,k}$ this quantity equals $M_i$ which is positive. Otherwise, it is made positive under the condition $\sigma^k |\xi| \le \alpha + M_i^0/M_i^{n+1,k}$ which gives the desired result.

(iii) We start to rewrite (3.25) as

$$ (3.33) \qquad (1 + \alpha)(f_i^{n+1,k+1} - M_i^{n+1,k}) = (M_i - M_i^{n+1,k}) - \sigma^k \xi (M_{i+1/2}^{n+1,k} - M_{i-1/2}^{n+1,k} + \delta M_{i+1/2-}^{n+1,k} - \delta M_{i-1/2+}^{n+1,k}). $$

The strategy is to multiply (3.33) by $\partial_f H(M_i^{n+1,k}, z_i)$ and to separate the discrete time variation from the flux and source contributions as

$$ (3.34) \qquad \partial_f H(M_i^{n+1,k}, z_i) [(1 + \alpha)(f_i^{n+1,k+1} - M_i^{n+1,k}) - (M_i - M_i^{n+1,k})] = - \sigma^k \partial_f H(M_i^{n+1,k}, z_i) [\xi(M_{i+1/2}^{n+1,k} - M_{i-1/2}^{n+1,k}) + \delta M_{i+1/2-}^{n+1,k} - \delta M_{i-1/2+}^{n+1,k}]. $$

We apply Lemma 2.4 to the left hand side of (3.34) to get

$$ (3.35) \qquad \partial_f H(M_i^{n+1,k}, z_i) [(1 + \alpha)(f_i^{n+1,k+1} - M_i^{n+1,k}) - (M_i - M_i^{n+1,k})] = $$
$$ - \xi H(M_{i \pm 1/2 \mp}, z_{i \pm 1/2}) + \left( \eta'(U_i) \cdot \left(\frac{1}{\xi}\right) + gz_i \right) \left( \xi M_{i \pm 1/2 \mp} - \xi M_i + \delta M_{i \pm 1/2 \mp} \right). \quad (3.28) $$

Moreover, the difference $\tilde{G}_{i+1/2-}^{n+1,k} - \tilde{G}_{i-1/2+}^{n+1,k}$ is macroscopically conservative

$$ \int_{\mathbb{R}} (\tilde{G}_{i+1/2-}^{n+1,k} - \tilde{G}_{i-1/2+}^{n+1,k}) d\xi = \overline{G}_{i+1/2}^{n+1,k} - \overline{G}_{i-1/2}^{n+1,k}, $$

with the macroscopic entropy flux

$$ \overline{G}_{i+1/2} = \int_{\mathbb{R}} \xi (1_{\xi>0} H(M_{i+1/2-}, z_{i+1/2}) + 1_{\xi<0} H(M_{i+1/2+}, z_{i+1/2})) d\xi. \quad (3.29) $$

We are now able to state our result.

**Proposition 3.8.** The scheme (3.24) is structure preserving, in the sense that we have the three properties below:

(i) the sub-iterations are well balanced, that is to say $U^{n+1,k} = U^n$ for all $k \in \mathbb{N}$ if $U^n$ is a discrete lake at rest;

(ii) assuming that the water height vectors $h^n$ and $h^{n+1,k}$ are positive, the update $h^{n+1,k+1}$ defined in the iterative scheme (3.24) is also positive if for all $1 \le i \le P$ the CFL condition $\sigma^k|\xi| \le \alpha + M_i/M_i^{n+1,k}$ holds for any $\xi$ belonging to supp $M^{n+1,k}$;

(iii) the kinetic entropy of the iterative process (3.25) verifies the following kinetic entropy inequality

$$ H(M_i^{n+1,k+1}, z_i) \le H(M_i, z_i) - \sigma^k (\tilde{G}_{i+1/2-}^{n+1,k} - \tilde{G}_{i-1/2+}^{n+1,k}) + (1+\alpha) \left( \eta'(U_i^{n+1,k+1}) \cdot \left(\frac{1}{\xi}\right) + gz_i \right) Q_i^{n+1,k+1} \\ + \alpha \left( H(M_i^{n+1,k}, z_i) - H(M_i^{n+1,k+1}, z_i) \right) + (1+\alpha) \left( \Psi(M_i^{n+1,k}, f_i^{n+1,k+1}) - \Psi(M_i^{n+1,k+1}, f_i^{n+1,k+1}) \right) - \Psi(M_i^{n+1,k}, f_i^{n+1,k+1}) - \Psi(M_i^{n+1,k}, M_i), \quad (3.30) $$

where $Q_i^{n+1,k+1} = M_i^{n+1,k+1} - f_i^{n+1,k+1}$ is a collision term verifying the conservation constraints (1.7), where $\Psi$ defined in (2.6) is positive and where $\tilde{G}_{i+1/2-}^{n+1,k}$ and $\tilde{G}_{i-1/2+}^{n+1,k}$ are defined by (3.28).

Before giving the proof we make the following remark.

**Remark 3.9.** We reiterate the comments made in Remark 3.6 which are to say that in (3.30) the term

$$ \left( \eta'(U_i^{n+1,k+1}) \cdot \left(\frac{1}{\xi}\right) + gz_i \right) Q_i^{n+1,k+1} $$

is macroscopically zero since $Q_i^{n+1,k+1}$ is a collision term satisfying the conservation constraints (1.7). Besides, assuming the method converges as $k \to \infty$, the quantity

$$ \alpha \left( H(M_i^{n+1,k}, z_i) - H(M_i^{n+1,k+1}, z_i) \right) + (1+\alpha) \left( \Psi(M_i^{n+1,k}, f_i^{n+1,k+1}) - \Psi(M_i^{n+1,k+1}, f_i^{n+1,k+1}) \right) $$

will eventually become negligible compared to $-\Psi(M_i^{n+1,k}, M_i) < 0$ from some rank $k$. Integrating inequality (3.30) over $\xi \in \mathbb{R}$, this implies that there exists $K \in \mathbb{N}$ such that for any $k \ge K$ the fully discrete entropy inequality

$$ \eta(U_i^{n+1,k+1}) \le \eta(U_i^n) - \sigma^k (\overline{G}_{i+1/2}^{n+1,k} - \overline{G}_{i-1/2}^{n+1,k}) \quad (3.31) $$
$$ (1 + \alpha)(H(f_i^{n+1,k+1}, z_i) - H(M_i^{n+1,k}, z_i) - \Psi(M_i^{n+1,k}, f_i^{n+1,k+1})) \\ - (H(M_i, z_i) - H(M_i^{n+1,k}, z_i) - \Psi(M_i^{n+1,k}, M_i)). $$

Furthermore, an upper bound on the right hand side of (3.34) is obtained by applying Proposition 3.8 which directly yields

$$ (3.36) \qquad - \partial_f H(M_i^{n+1,k}, z_i) [\xi(M_{i+1/2}^{n+1,k} - M_{i-1/2}^{n+1,k}) + \delta M_{i+1/2-}^{n+1,k} - \delta M_{i-1/2+}^{n+1,k}] \le \tilde{G}_{i+1/2-}^{n+1,k} - \tilde{G}_{i-1/2+}^{n+1,k}, $$

with $\tilde{G}_{i+1/2-}^{n+1,k}$ and $\tilde{G}_{i-1/2+}^{n+1,k}$ defined by (3.28). Injecting equality (3.35) and inequality (3.27) into (3.34) we obtain

$$ (3.37) \qquad (1 + \alpha)H(f_i^{n+1,k+1}, z_i) \le H(M_i, z_i) - \sigma^k(\tilde{G}_{i+1/2-}^{n+1,k} - \tilde{G}_{i-1/2+}^{n+1,k}) + \alpha H(M_i^{n+1,k}, z_i) \\ \qquad \qquad \qquad \qquad + (1 + \alpha)\Psi(M_i^{n+1,k}, f_i^{n+1,k+1}) - \Psi(M_i^{n+1,k}, M_i) $$

Using again Lemma 2.4 we can also write

$$ H(f_i^{n+1,k+1}, z_i) = H(M_i^{n+1,k+1}, z_i) + \partial_f H(M_i^{n+1,k+1}, z_i)(f_i^{n+1,k+1} - M_i^{n+1,k+1}) \\ + \Psi(M_i^{n+1,k+1}, f_i^{n+1,k+1}) \\ (3.38) \qquad \ge H(M_i^{n+1,k+1}, z_i) - \left(\eta'(U_i^{n+1,k+1}) \cdot \frac{1}{\xi} + gz_i\right)(M_i^{n+1,k+1} - f_i^{n+1,k+1}) \\ + \Psi(M_i^{n+1,k+1}, f_i^{n+1,k+1}), $$

where we used $\partial_f H = \partial_f H_0 + gz$ and (3.13) to get (3.38). Combining this inequality with (3.37) we finally get

$$ (1 + \alpha)H(M_i^{n+1,k+1}, z_i) \le H(M_i, z_i) - \sigma^k(\tilde{G}_{i+1/2-}^{n+1,k} - \tilde{G}_{i-1/2+}^{n+1,k}) + \alpha H(M_i^{n+1,k}, z_i) \\ + (1 + \alpha)\Psi(M_i^{n+1,k}, f_i^{n+1,k+1}) - \Psi(M_i^{n+1,k}, M_i) - (1 + \alpha)\Psi(M_i^{n+1,k+1}, f_i^{n+1,k+1}) \\ - (1 + \alpha)\left(\eta'(U_i^{n+1,k+1}) \cdot \frac{1}{\xi} + gz_i\right)(f_i^{n+1,k+1} - M_i^{n+1,k+1}). $$

After rearranging the terms and using $Q_i^{n+1,k+1} = -(f_i^{n+1,k+1} - M_i^{n+1,k+1})$ we obtain the desired kinetic entropy inequality (3.30).

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

The iterative kinetic scheme with hydrostatic reconstruction is validated through the numerical experiments performed in Section 5.2.

# 4 Towards the two dimensional Saint-Venant system

Being able to simulate the two dimensional Saint-Venant system is important if one wishes to deal with more relevant applications. In this Section, we investigate the possibility to extend the previous methods to the higher dimension, while retaining the good stability properties obtained in the one dimensional case. The following lines should be regarded as more exploratory, and in particular we leave for a future work the in-depth study of the two dimensional schemes that we briefly discuss below. In Section 4.1 we give the two dimensional Saint-Venant system and its kinetic interpretation. Then in Section 4.2 we derive a fully implicit scheme over flat bathymetry, and finally in Section 4.3 we propose an iterative kinetic scheme with hydrostatic reconstruction. In all cases a Cartesian mesh is considered.
is satisfied at the macroscopic level, with $\overline{G}_{i+1/2}^{n+1,k}$ given in (3.29). Summing inequality (3.31) over every cell $1 \le i \le P$ we obtain the dissipation of the total energy up to boundary fluxes

$$ \frac{1}{\Delta t^k} \sum_{i=1}^{P} \left( \eta(U_i^{n+1,k+1}) - \eta(U_i^n) \right) + \frac{1}{\Delta x} \left( \int_{\mathbb{R}} \xi H_{P+1/2}^{n+1,k}(\xi) \, d\xi - \int_{\mathbb{R}} \xi H_{1/2}^{n+1,k}(\xi) \, d\xi \right) \le 0. \quad (3.32) $$

In addition to the usual tolerance criterion where the iterations are stopped whenever two successive iterates are sufficiently close to each other, we can use (3.32) as a complementary condition to ensure the dissipation of total energy.

*Proof of Proposition 3.8.* The proof makes use of the kinetic writing (3.25) of scheme (3.24).

(i) The well-balancedness is a consequence of the hydrostatic reconstruction, see [5].

(ii) Remarking that the quantity $\delta M_{i+1/2-}^{n+1,k} - \delta M_{i-1/2+}^{n+1,k}$ appearing in the last line of (3.25) defines an odd function of $\xi - u_i^{n+1,k}$, its integral over $\xi \in \mathbb{R}$ vanishes and we have at the macroscopic level

$$ (1 + \alpha)h_i^{n+1,k+1} = \int_{\mathbb{R}} \left( M_i + \alpha M_i^{n+1,k} - \sigma^k \xi (M_{i+1/2}^{n+1,k} - M_{i-1/2}^{n+1,k}) \right) d\xi. $$

Thus it is enough to prove the positivity of the integrand, whose developed form is

$$ M_i + \alpha M_i^{n+1,k} - \sigma^k \xi (\mathbb{1}_{\xi>0} M_{i+1/2-}^{n+1,k} - \mathbb{1}_{\xi<0} M_{i-1/2+}^{n+1,k}) + \sigma^k \xi (\mathbb{1}_{\xi>0} M_{i-1/2-}^{n+1,k} - \mathbb{1}_{\xi<0} M_{i+1/2+}^{n+1,k}). $$

By definition of the water height reconstruction (3.20), we have the inequalities $h_{i+1/2-}^{n+1,k} \le h_i^{n+1,k}$ and $h_{i-1/2+}^{n+1,k} \le h_i^{n+1,k}$. As a consequence $M_{i+1/2-}^{n+1,k} \le M_i^{n+1,k}$ and $M_{i-1/2+}^{n+1,k} \le M_i^{n+1,k}$, which allows us to bound the integrand from below by

$$ M_i + \alpha M_i^{n+1,k} - \sigma^k |\xi| M_i^{n+1,k}. $$

If $\xi$ does not belong to supp $M_i^{n+1,k}$ this quantity equals $M_i$ which is positive. Otherwise, it is made positive under the condition $\sigma^k |\xi| \le \alpha + M_i^0/M_i^{n+1,k}$ which gives the desired result.

(iii) We start to rewrite (3.25) as

$$ (1 + \alpha) (f_i^{n+1,k+1} - M_i^{n+1,k}) = (M_i - M_i^{n+1,k}) - \sigma^k \xi (M_{i+1/2}^{n+1,k} - M_{i-1/2}^{n+1,k} + \delta M_{i+1/2-}^{n+1,k} - \delta M_{i-1/2+}^{n+1,k}). \quad (3.33) $$

The strategy is to multiply (3.33) by $\partial_f H(M_i^{n+1,k}, z_i)$ and to separate the discrete time variation from the flux and source contributions as

$$ \begin{aligned} & \partial_f H(M_i^{n+1,k}, z_i) \left[ (1 + \alpha) (f_i^{n+1,k+1} - M_i^{n+1,k}) - (M_i - M_i^{n+1,k}) \right] \\ & = -\sigma^k \partial_f H(M_i^{n+1,k}, z_i) \left[ \xi (M_{i+1/2}^{n+1,k} - M_{i-1/2}^{n+1,k}) + \delta M_{i+1/2-}^{n+1,k} - \delta M_{i-1/2+}^{n+1,k} \right]. \end{aligned} \quad (3.34) $$

We apply Lemma 2.4 to the left hand side of (3.34) to get

$$ \begin{aligned} & \partial_f H(M_i^{n+1,k}, z_i) \left[ (1 + \alpha) (f_i^{n+1,k+1} - M_i^{n+1,k}) - (M_i - M_i^{n+1,k}) \right] \\ & = (1 + \alpha) \left( H(f_i^{n+1,k+1}, z_i) - H(M_i^{n+1,k}, z_i) - \Psi(M_i^{n+1,k}, f_i^{n+1,k+1}) \right) \\ & \quad - \left( H(M_i, z_i) - H(M_i^{n+1,k}, z_i) - \Psi(M_i^{n+1,k}, M_i) \right). \end{aligned} \quad (3.35) $$
## 4.1 The two dimensional system and its kinetic representation

We consider the two dimensional Saint-Venant system written, with obvious notations, under the form

$$ (4.1) \qquad \left\{ \begin{aligned} \frac{\partial h}{\partial t} + \frac{\partial hu}{\partial x} + \frac{\partial hv}{\partial y} &= 0 \\ \frac{\partial hu}{\partial t} + \frac{\partial}{\partial x} \left(hu^2 + \frac{g}{2}h^2\right) + \frac{\partial}{\partial y} (huv) &= -gh \frac{\partial z}{\partial x} \\ \frac{\partial hv}{\partial t} + \frac{\partial}{\partial x} (huv) + \frac{\partial}{\partial y} \left(hv^2 + \frac{g}{2}h^2\right) &= -gh \frac{\partial z}{\partial y} \end{aligned} \right. . $$

with $(u, v)$ the horizontal velocity vector. A straightforward extension of Lemma 1.1 allows to get a kinetic interpretation of system (4.1) as studied in [6, 1]. To build the two dimensional Gibbs equilibrium, we define the function

$$ (4.2) \qquad \chi(z_1, z_2) = \frac{1}{4\pi} \mathbb{1}_{z_1^2+z_2^2 \le 4}, $$

and we have

$$ (4.3) \quad \forall U = (h, hu, hv)^T, \quad \forall \xi \in \mathbb{R}^2, \quad M(U, \xi) = \frac{h}{c^2} \chi\left(\frac{\xi_1 - u}{c}, \frac{\xi_2 - v}{c}\right) = \frac{1}{2g\pi} \mathbb{1}_{|\xi - (u, v)^T|_2 \le \sqrt{2gh}}, $$

with the velocity $c = \sqrt{\frac{g}{2}h}$. Let us remark that when averaging this Maxwellian function in the direction $\xi_1$ (resp. $\xi_2$), we get a function of $\xi_2$ (resp. $\xi_1$) that coincides with the half-disk Maxwellian (1.4) used in the one dimensional case. We then state the following lemma.

**Lemma 4.1** If the topography $z(x, y)$ is Lipschitz continuous, then $U = (h, hu, hv)^T$ is a weak solution to the Saint-Venant system (4.1) if and only if $M(U, \xi)$ satisfies the kinetic equation

$$ (4.4) \qquad \partial_t M + \xi \cdot \nabla_{(x,y)} M - g \nabla_{(x,y)} z \cdot \nabla_{(\xi_1, \xi_2)} M = Q, $$

for some "collision term" $Q(t, x, y, \xi)$ that satisfies, for a.e. $(t, x, y)$,

$$ (4.5) \qquad \int_{\mathbb{R}^2} Q \, d\xi_1 d\xi_2 = \int_{\mathbb{R}^2} \xi_1 Q \, d\xi_1 d\xi_2 = \int_{\mathbb{R}^2} \xi_2 Q \, d\xi_1 d\xi_2 = 0. $$

**Proof of Lemma 4.1.** The proof relies on simple computations. Classically, the integral of Equation (4.4) over $\xi \in \mathbb{R}^2$ gives the first line of the system (4.1) whereas the last two lines of this system are obtained by taking the scalar product of (4.4) against $\xi \in \mathbb{R}^2$ and by integrating over $\mathbb{R}^2$. $\square$

As in the one dimensional case, the Maxwellian (4.3) is associated with a kinetic entropy $H(f, \xi)$ satisfying the two dimensional analogous of Lemma 1.2, see [5]. This kinetic entropy is given by

$$ (4.6) \qquad H(f, \xi) = \frac{|\xi|^2}{2} f + gzf. $$

## 4.2 Fully implicit kinetic scheme over a flat bathymetry

Let us consider a bounded rectangular domain $\Omega = (0, L_x) \times (0, L_y)$ and its Cartesian discretization using $P$ cells in both $x$- and $y$-directions, such that the total number of cells is $P^2$. We denote $(P_{i,j})_{0 \le i,j \le P}$ the vertices with coordinates $(x_i, y_j)^T$ given by

$$ x_i = (i + 1/2)\Delta x, \quad y_j = (j + 1/2)\Delta y, $$

where $\Delta x = L_x/P$, $\Delta y = L_y/P$. We use the following notations:

* for all $(i, j) \in [1, P]^2$, $C_{i,j}$ is the rectangular cell centered on $P_{i,j}$ with area $|C_{i,j}| = \Delta x \Delta y$,