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

$$ - \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] \leq -\tilde{G}_{i+1/2-}^{n+1,k} + \tilde{G}_{i-1/2+}^{n+1,k}, \quad (3.36) $$

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

$$ (1 + \alpha)H(f_i^{n+1,k+1}, z_i) \leq 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). \quad (3.37) $$

Using again Lemma 2.4 we can also write

$$ \begin{aligned} 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}) \\ &\quad + \Psi(M_i^{n+1,k+1}, f_i^{n+1,k+1}) \\ &\geq H(M_i^{n+1,k+1}, z_i) - \left(\eta'(U_i^{n+1,k+1}) \cdot \left(\frac{1}{\xi}\right) + gz_i\right) (M_i^{n+1,k+1} - f_i^{n+1,k+1}) \\ &\quad + \Psi(M_i^{n+1,k+1}, f_i^{n+1,k+1}), \end{aligned} \quad (3.38) $$

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) \leq 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 \left(\frac{1}{\xi}\right) + gz_i\right)\left(f_i^{n+1,k+1} - M_i^{n+1,k+1}\right). $$

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_3.png" /></figure>

Compared to the flat topography case, proving the convergence of the iterative kinetic scheme with hydrostatic reconstruction (3.24) as $k \to \infty$ is significantly more difficult owing to the reconstruction step. Instead we investigate this question through numerical experiments in Section 5, where in all considered test cases the sub-iterations seem to converge with no issue.

## 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.
• $\partial C_{i,j}$ is the boundary of $C_{i,j}$,

• indices $(i, j)$ such that $i \in \{0, P+1\}$ or $j \in \{0, P+1\}$ indicate a ghost cell outside of $\Omega$.

We define the piecewise constant functions $U^n(x, y)$ and $z(x, y)$ on cells $C_{i,j}$ as

$$ (4.7) \qquad U^n(x, y) = U_{i,j}^n, \quad z(x, y) = z_{i,j}, \quad \text{for } (x, y) \in C_{i,j}, $$

with $U_{i,j}^n = (h_{i,j}^n, (hu)_{i,j}^n, (hv)_{i,j}^n)^T$ such that they approximate the cell averages as below

$$ U_{i,j}^n \approx \frac{1}{|C_{i,j}|} \int_{C_{i,j}} U(t^n, x, y) \, dxdy, \quad z_{i,j} \approx \frac{1}{|C_{i,j}|} \int_{C_{i,j}} z(x, y) \, dxdy. $$

In the case of a flat topography, the integral over $C_{i,j}$ of the convective part of the kinetic equation (4.4) gives

$$ (4.8) \qquad \int_{C_{i,j}} \left( \frac{\partial M}{\partial t} + \xi \cdot \nabla_{(x,y)} M \right) dxdy \approx |C_{i,j}| \frac{\partial M_{i,j}}{\partial t} + \int_{\partial C_{i,j}} (\xi \cdot n \partial_{C_{i,j}}) M_{i,j} \, d\ell \\ \approx |C_{i,j}| \frac{\partial M_{i,j}}{\partial t} + \xi_1 \Delta y (1_{\xi_1>0} (M_{i,j} - M_{i-1,j}) + 1_{\xi_1<0} (M_{i+1,j} - M_{i,j})) \\ + \xi_2 \Delta x (1_{\xi_2>0} (M_{i,j} - M_{i,j-1}) + 1_{\xi_2<0} (M_{i,j+1} - M_{i,j})), $$

with $M_{i,j}(t) = M(U_{i,j}(t), \xi)$ and $n_{\partial C_{i,j}}(\ell)$ the normal to $\partial C_{i,j}$ at the contour abscissa $\ell \in \partial C_{i,j}$ taken outward to $C_{i,j}$. From the spatial semi-discretization (4.8) of the convection operator, we deduce an implicit Euler scheme for the kinetic interpretation (4.4)

$$ (4.9) \qquad f_{i,j}^{n+1} = M_{i,j}^n - \sigma_x \xi_1 (f_{i+1/2,j}^{n+1} - f_{i-1/2,j}^{n+1}) - \sigma_y \xi_2 (f_{i,j+1/2}^{n+1} - f_{i,j-1/2}^{n+1}), $$

where we used $\sigma_x = \Delta t^n / \Delta x$, $\sigma_y = \Delta t^n / \Delta y$ and the interfacial kinetic densities

$$ f_{i+1/2,j}^{n+1} = 1_{\xi_1>0} f_{i,j}^{n+1} + 1_{\xi_1<0} f_{i+1,j}^{n+1}, \qquad f_{i,j+1/2}^{n+1} = 1_{\xi_2>0} f_{i,j}^{n+1} + 1_{\xi_2<0} f_{i,j+1}^{n+1}. $$

Denoting $f$ the vector of interior values from $\mathbb{R}^{P^2}$

$$ f = (f_{1,1}, f_{2,1}, \dots, f_{P,1}, f_{1,2}, \dots, f_{P,P})^T, $$

the kinetic scheme (4.9) also writes

$$ (4.10) \qquad (\mathbf{I}_{P2} + \mathbf{L}_{P2}) f^{n+1} = M + B^{n+1}, $$

where we have used the particular geometry of the mesh and with $\mathbf{I}_{P2}$ is the identity matrix of length $P^2$, and the block matrix $\mathbf{L}_{P2}$ is defined by

$$ \mathbf{L}_{P2} = \begin{pmatrix} \mathbf{D} & \mathbf{N}^{+} & 0 & \dots & 0 \\ \mathbf{N}^{-} & \mathbf{D} & \mathbf{N}^{+} & \dots & \vdots \\ 0 & \dots & \dots & \dots & 0 \\ \vdots & \dots & \mathbf{N}^{-} & \mathbf{D} & \mathbf{N}^{+} \\ 0 & \dots & 0 & \mathbf{N}^{-} & \mathbf{D} \end{pmatrix}, $$

where $\mathbf{D}, \mathbf{N}^{\pm}$ are $P \times P$ matrices defined by $\mathbf{N}^{+} = -\sigma_y \xi_2 1_{\xi_2<0} \mathbf{I}_P$, $\mathbf{N}^{-} = \sigma_y \xi_2 1_{\xi_2>0} \mathbf{I}_P$ and

$$ \mathbf{D} = \sigma_x |\xi_1| \begin{pmatrix} 1 & -1_{\xi_1<0} & 0 & \dots & 0 \\ -1_{\xi_1>0} & 1 & -1_{\xi_1<0} & \dots & \vdots \\ 0 & \dots & \dots & \dots & 0 \\ \vdots & \dots & -1_{\xi_1>0} & 1 & -1_{\xi_1<0} \\ 0 & \dots & 0 & -1_{\xi_1>0} & 1 \end{pmatrix} + \sigma_y |\xi_2| \mathbf{I}_P. $$
## 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

$$ 
\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( h u^2 + \frac{g}{2} h^2 \right) + \frac{\partial}{\partial y} (h u v) &= -g h \frac{\partial z}{\partial x} \\
\frac{\partial hv}{\partial t} + \frac{\partial}{\partial x} (h u v) + \frac{\partial}{\partial y} \left( h v^2 + \frac{g}{2} h^2 \right) &= -g h \frac{\partial z}{\partial y}
\end{aligned}
\right. 
\qquad (4.1) $$

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 [1, 4]. To build the two dimensional Gibbs equilibrium, we define the function

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

and we have

$$ \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 \le \sqrt{2gh}}, \qquad (4.3) $$

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

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

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

$$ \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. \qquad (4.5) $$

*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 [3]. This kinetic entropy is given by

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

## 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:
Likewise, the vector $B^{n+1}(\xi) \in \mathbb{R}^P$ accounting for the boundary conditions is given by

$$ B^{n+1} = \begin{pmatrix} \mathbf{N}^{-} & \mathbf{D}_{1}^{-} & \mathbf{D}_{1}^{+} & 0 \\ 0 & \mathbf{D}_{2}^{-} & \mathbf{D}_{2}^{+} & \vdots \\ \vdots & \vdots & \vdots & 0 \\ 0 & \mathbf{D}_{P}^{-} & \mathbf{D}_{P}^{+} & \mathbf{N}^{+} \end{pmatrix} \begin{pmatrix} M_{\text{bottom}} \\ M_{\text{left}} \\ M_{\text{right}} \\ M_{\text{top}} \end{pmatrix}, \quad \left\{ \begin{array}{ll} M_{\text{bottom}} & = (M_{1,0}^{n+1}, \dots, M_{P,0}^{n+1})^{T} \\ M_{\text{left}} & = (M_{0,1}^{n+1}, \dots, M_{0,P}^{n+1})^{T} \\ M_{\text{right}} & = (M_{P+1,1}^{n+1}, \dots, M_{P+1,P}^{n+1})^{T}, \\ M_{\text{top}} & = (M_{1,P+1}^{n+1}, \dots, M_{P,P+1}^{n+1})^{T} \end{array} \right. $$

with $(\mathbf{D}_k^-)_{i,j} = \delta_{1,i} \delta_{k,j} \sigma_x |\xi_1| \mathbf{1}_{\xi_1>0}$ and $(\mathbf{D}_k^+)_{i,j} = \delta_{P,i} \delta_{k,j} \sigma_x |\xi_1| \mathbf{1}_{\xi_1<0}$ for all $(i, j, k) \in [1, P]^3$. Since the matrix $\mathbf{I}_{P2} + \mathbf{L}_{P2}$ has the same structure as the matrix $\mathbf{I} + \sigma \mathbf{L}$ studied in Lemma 2.1, and since the kinetic representation (4.4) admits the kinetic entropy (4.6), the results of Propositions 2.2 and 2.3 are still valid. The practical computation of the inverse of the matrix $\mathbf{I}_{P2} + \mathbf{L}_{P2}$, of the numerical fluxes at the macroscopic level, and the implementation of the scheme is left for a later work.

## 4.3 Iterative kinetic scheme with hydrostatic reconstruction

The iterative kinetic scheme (3.2) and its version with hydrostatic reconstruction (3.25) can be extended to the two dimensional case. First we need to detail the reconstruction step over a Cartesian mesh. Note that it is also possible to apply it on unstructured meshes as detailed in [6]. The idea is to apply formulas (3.19)-(3.20) along the normal of a given edge separating two neighboring cells. More precisely we have

$$ (4.11) \qquad \forall 1 \le i, j \le P, \quad U_{i \pm 1/2 \mp, j} = \begin{pmatrix} 1 \\ u_{i,j} \\ v_{i,j} \end{pmatrix} h_{i \pm 1/2 \mp, j}, \quad U_{i, j \pm 1/2 \mp} = \begin{pmatrix} 1 \\ u_{i,j} \\ v_{i,j} \end{pmatrix} h_{i, j \pm 1/2 \mp}, $$

with the reconstructed interfacial water heights given by

$$ (4.12) \qquad h_{i \pm 1/2 \mp, j} = (h_{i,j} + z_{i,j} - z_{i \pm 1/2, j})_{+}, \quad h_{i, j \pm 1/2 \mp} = (h_{i,j} + z_{i,j} - z_{i, j \pm 1/2})_{+}, $$

and with the interfacial bathymetry $z_{i+1/2,j} = \max(z_{i,j}, z_{i+1,j})$ and $z_{i,j+1/2} = \max(z_{i,j}, z_{i,j+1})$. It is then possible to consider the Maxwellian densities associated with the reconstructed interfacial states (4.11) through an upwinding as follows

$$ \begin{aligned} M_{i+1/2,j}(\xi) &= \mathbf{1}_{\xi_1>0} M(U_{i+1/2-,j}, \xi) + \mathbf{1}_{\xi_1<0} M(U_{i+1/2+,j}, \xi), \\ M_{i,j+1/2}(\xi) &= \mathbf{1}_{\xi_2>0} M(U_{i,j+1/2-,j}, \xi) + \mathbf{1}_{\xi_2<0} M(U_{i,j+1/2+,j}, \xi). \end{aligned} $$

An implicit hydrostatic reconstruction kinetic scheme stemming from (4.9) consists to find macroscopic states $(U_{i,j}^{n+1})_{1 \le i,j \le P}$ such that

$$ \begin{aligned} M_{i,j}^{n+1} = M_{i,j}^{n} &- \sigma_x \xi_1 (M_{i+1/2,j}^{n+1} - M_{i-1/2,j}^{n+1}) - \sigma_y \xi_2 (M_{i,j+1/2}^{n+1} - M_{i,j-1/2}^{n+1}) + Q_{i,j}^{n+1} \\ &+ \sigma_x (\xi_1 - u_{i,j}^{n+1,k}) (M_{i+1/2-,j}^{n+1,k} - M_{i-1/2+,j}^{n+1,k}) + \sigma_y (\xi_2 - v_{i,j}^{n+1,k}) (M_{i,j+1/2-}^{n+1,k} - M_{i,j-1/2+}^{n+1,k}), \end{aligned} $$

where $Q_{i,j}^{n+1}$ is a collision term verifying (4.5), and where the last line of the above expression accounts for the topography source term. As in the one dimensional case, it is not possible to compute explicitly this update, instead we approximate it by a fixed point method

$$ (4.13) \qquad (1 + \alpha) f_{i,j}^{n+1,k+1} = M_{i,j}^{n} + \alpha M_{i,j}^{n+1,k} - \sigma_x \xi_1 (M_{i+1/2,j}^{n+1,k} - M_{i-1/2,j}^{n+1,k}) - \sigma_y \xi_2 (M_{i,j+1/2}^{n+1,k} - M_{i,j-1/2}^{n+1,k}) \\ \qquad \qquad + \sigma_x (\xi_1 - u_{i,j}^{n+1,k}) (M_{i+1/2-,j}^{n+1,k} - M_{i-1/2+,j}^{n+1,k}) + \sigma_y (\xi_2 - v_{i,j}^{n+1,k}) (M_{i,j+1/2-}^{n+1,k} - M_{i,j-1/2+}^{n+1,k}), $$

where $\alpha > 0$ is a relaxation parameter, combined with the projection step

$$ (4.14) \qquad M_{i,j}^{n+1,k+1} = f_{i,j}^{n+1,k+1} + Q_{i,j}^{n+1,k}. $$
- 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$,

- $\partial C_{i,j}$ is the boundary of $C_{i,j}$,

- indices $(i, j)$ such that $i \in \{0, P+1\}$ or $j \in \{0, P+1\}$ indicate a ghost cell outside of $\Omega$.

We define the piecewise constant functions $U^n(x, y)$ and $z(x, y)$ on cells $C_{i,j}$ as

$$ U^n(x, y) = U_{i,j}^n, \quad z(x, y) = z_{i,j}, \quad \text{for } (x, y) \in C_{i,j}, \qquad (4.7) $$

with $U_{i,j}^n = (h_{i,j}^n, (hu)_{i,j}^n, (hv)_{i,j}^n)^T$ such that they approximate the cell averages as below

$$ U_{i,j}^n \approx \frac{1}{|C_{i,j}|} \int_{C_{i,j}} U(t^n, x, y) \, dx \, dy, \quad z_{i,j} \approx \frac{1}{|C_{i,j}|} \int_{C_{i,j}} z(x, y) \, dx \, dy. $$

In the case of a flat topography, the integral over $C_{i,j}$ of the convective part of the kinetic equation (4.4) gives

$$ \begin{aligned} \int_{C_{i,j}} \left( \frac{\partial M}{\partial t} + \xi \cdot \nabla_{(x,y)} M \right) dx dy &\approx |C_{i,j}| \frac{\partial M_{i,j}}{\partial t} + \int_{\partial C_{i,j}} (\xi \cdot n \partial C_{i,j}) M_{i,j} d\ell \\ &\approx |C_{i,j}| \frac{\partial M_{i,j}}{\partial t} + \xi_1 \Delta y (1_{\xi_1>0} (M_{i,j} - M_{i-1,j}) + 1_{\xi_1<0} (M_{i+1,j} - M_{i,j})) \\ &\quad + \xi_2 \Delta x (1_{\xi_2>0} (M_{i,j} - M_{i,j-1}) + 1_{\xi_2<0} (M_{i,j+1} - M_{i,j})), \end{aligned} \qquad (4.8) $$

with $M_{i,j}(t) = M(U_{i,j}(t), \xi)$ and $n_{\partial C_{i,j}}(\ell)$ the normal to $\partial C_{i,j}$ at the contour abscissa $\ell \in \partial C_{i,j}$ taken outward to $C_{i,j}$. From the spatial semi-discretization (4.8) of the convection operator, we deduce an implicit Euler scheme for the kinetic interpretation (4.4)

$$ f_{i,j}^{n+1} = M_{i,j}^{n} - \sigma_{x} \xi_{1} (f_{i+1/2,j}^{n+1} - f_{i-1/2,j}^{n+1}) - \sigma_{y} \xi_{2} (f_{i,j+1/2}^{n+1} - f_{i,j-1/2}^{n+1}), \qquad (4.9) $$

where we used $\sigma_x = \Delta t^n / \Delta x$, $\sigma_y = \Delta t^n / \Delta y$ and the interfacial kinetic densities

$$ f_{i+1/2,j}^{n+1} = 1_{\xi_1>0} f_{i,j}^{n+1} + 1_{\xi_1<0} f_{i+1,j}^{n+1}, \quad f_{i,j+1/2}^{n+1} = 1_{\xi_2>0} f_{i,j}^{n+1} + 1_{\xi_2<0} f_{i,j+1}^{n+1}. $$

Denoting $f$ the vector of interior values from $\mathbb{R}^{P^2}$

$$ f = (f_{1,1}, f_{2,1}, \dots, f_{P,1}, f_{1,2}, \dots, f_{P,P})^T, $$

the kinetic scheme (4.9) also writes

$$ (\mathbf{I}_{P2} + \mathbf{L}_{P2}) f^{n+1} = M + B^{n+1}, \qquad (4.10) $$

where we have used the particular geometry of the mesh and with $\mathbf{I}_{P2}$ is the identity matrix of length $P^2$, and the block matrix $\mathbf{L}_{P2}$ is defined by

$$ \mathbf{L}_{P2} = \begin{pmatrix} \mathbf{D} & \mathbf{N}^{+} & 0 & \dots & 0 \\ \mathbf{N}^{-} & \mathbf{D} & \mathbf{N}^{+} & \dots & \vdots \\ 0 & \dots & \dots & \dots & 0 \\ \vdots & \ddots & \mathbf{N}^{-} & \mathbf{D} & \mathbf{N}^{+} \\ 0 & \dots & 0 & \mathbf{N}^{-} & \mathbf{D} \end{pmatrix}, $$
At the macroscopic level, an iterative scheme is obtained by integrating (4.13) against $(1, \xi)$ which is implemented in practice. It reads

$$ (4.15) \qquad (1 + \alpha) U_{i,j}^{n+1,k+1} = U_{i,j}^{n} + \alpha U_{i,j}^{n+1,k} - \sigma_x (F_{x,i+1/2,-,j}^{n+1,k} - F_{x,i-1/2+,j}^{n+1,k}) - \sigma_y (F_{y,i,j+1/2-}^{n+1,k} - F_{y,i,j-1/2+}^{n+1,k}), $$

where we have the following definition of the numerical fluxes

$$ F_{x,i\pm1/2\mp,j}^{n+1,k} = \int_{\mathbb{R}^2} \xi_1 \begin{pmatrix} 1 \\ \xi_1 \\ \xi_2 \end{pmatrix} M_{i\pm1/2,j}^{n+1,k} d\xi_1 d\xi_2 + \frac{g}{2} \begin{pmatrix} 0 \\ (h_{i,j}^{n+1,k})^2 - (h_{i\pm1/2\mp,j}^{n+1,k})^2 \\ 0 \end{pmatrix}, $$

$$ F_{y,i,j\pm1/2\mp}^{n+1,k} = \int_{\mathbb{R}^2} \xi_2 \begin{pmatrix} 1 \\ \xi_1 \\ \xi_2 \end{pmatrix} M_{i,j\pm1/2}^{n+1,k} d\xi_1 d\xi_2 + \frac{g}{2} \begin{pmatrix} 0 \\ 0 \\ (h_{i,j}^{n+1,k})^2 - (h_{i,j\pm1/2\mp}^{n+1,k})^2 \end{pmatrix}. $$

The results obtained in Section 3 for the 1d Saint-Venant remain valid in the two dimensional setting since, the kinetic representation (4.4) of the Saint-Venant equations with Maxwellian (4.3) possesses a convex kinetic entropy $H$ given by (4.6).

The complete analysis of the 2d scheme is left for a forthcoming paper but a numerical experiment consisting of the two dimensional parabolic bowl [17] is presented in Section 5.3.

# 5 Numerical simulations

In this section we evaluate the behavior and the efficiency of the proposed kinetic schemes through numerical experiments. In Section 5.1, we focus on the fully implicit kinetic scheme over a flat bathymetry that corresponds to (2.1). Then, in Section 5.2, the iterative kinetic scheme with hydrostatic reconstruction is investigated numerically. In particular, we exhibit a testcase where our iterative method succeeds to dissipate the total energy while its explicit version increases it. Finally in Section 5.3 a two dimensional test case is proposed.

## 5.1 Fully implicit kinetic scheme in the flat bottom case

<figure><img src="image_5.png" /><figcaption>Figure 2: Slow moving shock approximated by various kinetic schemes, including explicit, implicit and iterative strategies. The initial condition is given by a Riemann data with discontinuity at position $x = 0$.</figcaption></figure>
where $\mathbf{D}$, $\mathbf{N}^{\pm}$ are $P \times P$ matrices defined by $\mathbf{N}^{+} = -\sigma_y \xi_2 \mathbf{1}_{\xi_2<0} \mathbf{I}_P$, $\mathbf{N}^{-} = \sigma_y \xi_2 \mathbf{1}_{\xi_2>0} \mathbf{I}_P$ and

$$ \mathbf{D} = \sigma_x |\xi_1| \begin{pmatrix} 1 & -\mathbf{1}_{\xi_1<0} & 0 & \dots & 0 \\ -\mathbf{1}_{\xi_1>0} & 1 & -\mathbf{1}_{\xi_1<0} & \dots & \vdots \\ 0 & \dots & \dots & \dots & 0 \\ \vdots & \dots & -\mathbf{1}_{\xi_1>0} & 1 & -\mathbf{1}_{\xi_1<0} \\ 0 & \dots & 0 & -\mathbf{1}_{\xi_1>0} & 1 \end{pmatrix} + \sigma_y |\xi_2| \mathbf{I}_P. $$

Likewise, the vector $B^{n+1}(\xi) \in \mathbb{R}^P$ accounting for the boundary conditions is given by

$$ B^{n+1} = \begin{pmatrix} \mathbf{N}^{-} & \mathbf{D}_{1}^{-} & \mathbf{D}_{1}^{+} & 0 \\ 0 & \mathbf{D}_{2}^{-} & \mathbf{D}_{2}^{+} & \vdots \\ \vdots & \vdots & \vdots & 0 \\ 0 & \mathbf{D}_{P}^{-} & \mathbf{D}_{P}^{+} & \mathbf{N}^{+} \end{pmatrix} \begin{pmatrix} M_{\text{bottom}} \\ M_{\text{left}} \\ M_{\text{right}} \\ M_{\text{top}} \end{pmatrix}, \quad \left\{ \begin{array}{l} M_{\text{bottom}} = (M_{1,0}^{n+1}, \dots, M_{P,0}^{n+1})^{\mathrm{T}} \\ M_{\text{left}} = (M_{0,1}^{n+1}, \dots, M_{0,P}^{n+1})^{\mathrm{T}} \\ M_{\text{right}} = (M_{P+1,1}^{n+1}, \dots, M_{P+1,P}^{n+1})^{\mathrm{T}} \\ M_{\text{top}} = (M_{1,P+1}^{n+1}, \dots, M_{P,P+1}^{n+1})^{\mathrm{T}} \end{array} \right. $$

with $(\mathbf{D}_k^-)_{i,j} = \delta_{1,i} \delta_{k,j} \sigma_x |\xi_1| \mathbf{1}_{\xi_1>0}$ and $(\mathbf{D}_k^+)_{i,j} = \delta_{P,i} \delta_{k,j} \sigma_x |\xi_1| \mathbf{1}_{\xi_1<0}$ for all $(i,j,k) \in \llbracket 1, P \rrbracket^3$. Since the matrix $\mathbf{I}_{P2} + \mathbf{L}_{P2}$ has the same structure as the matrix $\mathbf{I} + \sigma \mathbf{L}$ studied in Lemma 2.1, and since the kinetic representation (4.4) admits the kinetic entropy (4.6), the results of Propositions 2.2 and 2.3 are still valid. The practical computation of the inverse of the matrix $\mathbf{I}_{P2} + \mathbf{L}_{P2}$, of the numerical fluxes at the macroscopic level, and the implementation of the scheme is left for a later work.

### 4.3. Iterative kinetic scheme with hydrostatic reconstruction

The iterative kinetic scheme (3.2) and its version with hydrostatic reconstruction (3.25) can be extended to the two dimensional case. First we need to detail the reconstruction step over a Cartesian mesh. Note that it is also possible to apply it on unstructured meshes as detailed in [4]. The idea is to apply formulas (3.19) and (3.20) along the normal of a given edge separating two neighboring cells. More precisely we have

$$ \forall 1 \le i, j \le P, \quad U_{i \pm 1/2 \mp, j} = \begin{pmatrix} 1 \\ u_{i,j} \\ v_{i,j} \end{pmatrix} h_{i \pm 1/2 \mp, j}, \quad U_{i,j \pm 1/2 \mp} = \begin{pmatrix} 1 \\ u_{i,j} \\ v_{i,j} \end{pmatrix} h_{i,j \pm 1/2 \mp}, \qquad (4.11) $$

with the reconstructed interfacial water heights given by

$$ h_{i \pm 1/2 \mp, j} = (h_{i,j} + z_{i,j} - z_{i \pm 1/2, j})_{+}, \quad h_{i,j \pm 1/2 \mp} = (h_{i,j} + z_{i,j} - z_{i,j \pm 1/2})_{+}, \qquad (4.12) $$

and with the interfacial bathymetry $z_{i+1/2,j} = \max(z_{i,j}, z_{i+1,j})$ and $z_{i,j+1/2} = \max(z_{i,j}, z_{i,j+1})$. It is then possible to consider the Maxwellian densities associated with the reconstructed interfacial states (4.11) through an upwinding as follows

$$ M_{i+1/2,j}(\xi) = \mathbf{1}_{\xi_1>0} M(U_{i+1/2-,j}, \xi) + \mathbf{1}_{\xi_1<0} M(U_{i+1/2+,j}, \xi), $$

$$ M_{i,j+1/2}(\xi) = \mathbf{1}_{\xi_2>0} M(U_{i,j+1/2-}, \xi) + \mathbf{1}_{\xi_2<0} M(U_{i,j+1/2+}, \xi). $$

An implicit hydrostatic reconstruction kinetic scheme stemming from (4.9) consists to find macroscopic states $(U_{i,j}^{n+1})_{1 \le i,j \le P}$ such that

$$ M_{i,j}^{n+1} = M_{i,j}^{n} - \sigma_x \xi_1 (M_{i+1/2,j}^{n+1} - M_{i-1/2,j}^{n+1}) - \sigma_y \xi_2 (M_{i,j+1/2}^{n+1} - M_{i,j-1/2}^{n+1}) + Q_{i,j}^{n+1} $$
**Slow moving shock.** To assess the efficiency and interest of the implicit scheme (2.2), we perform a numerical test involving a Riemann problem with a slowly moving shock over a flat bottom. This configuration is achieved for a nearly transcritical flow where the material velocity $u$ is positive and satisfies $u - \sqrt{gh} \approx 0$ and $u + \sqrt{gh} \gg 1$. Hence the maximum eigenvalue severely constrains the time step, however a small time step might not be necessary to accurately resolve the slow shock. We set the gravitational acceleration $g$ to 10 and the Riemann problem corresponds to the initial state $U_L = (1, 4.75)^T$ for $x < 0$ and $U_R = (2, 4.75)^T$ for $x > 0$. In Figure 2 we compare several schemes with an explicit time step $\Delta t_{\text{exp}}$ given by the CFL condition $\Delta t_{\text{exp}} \le 0.45 \frac{\Delta x}{\lambda_{\text{max}}}$, as well as the implicit kinetic scheme using a time step $\Delta t_{\text{imp}} = 10\Delta t_{\text{exp}}$. We also use the iterative scheme (3.3) with parameter $\alpha = 1$, and plot the results at time $t = 0.5$. We notice that in the discharge profile, an oscillation appears downwind of the left-going shock, which is quite pronounced for the explicit and iterative kinetic schemes, and less so for the fully implicit ones. As expected the implicit scheme using $\Delta t_{\text{imp}}$ strongly diffuses the fast traveling rarefaction. On the other hand the slow shock seems to be slightly less impacted by the large time steps, however it is still less diffused when using $\Delta t_{\text{exp}}$. Despite requiring ten times less iterations to reach the final time, the use of large time increments for the implicit kinetic scheme only results in around two percents faster computations compared to the explicit strategy which is due to the high quadratic cost of the implicit method. We believe that it is not possible to lower this cost when it comes to unconditionally stable methods, because the associated stencil has to cover the entire computational domain.

<figure><img src="image_6.png" /><figcaption>Figure 3: Comparing implicit, iterative and explicit kinetic solvers on a Riemann problem.</figcaption></figure>

**Riemann problem.** We compare the fully implicit kinetic scheme and iterative kinetic scheme to explicit methods. The testcase is given by the Riemann problem with initial data $U^0(x) = \mathbb{1}_{x<0}U_L +$
$$ + \sigma_x (\xi_1 - u_{i,j}^{n+1,k}) (M_{i+1/2-,j}^{n+1,k} - M_{i-1/2+,j}^{n+1,k}) + \sigma_y (\xi_2 - v_{i,j}^{n+1,k}) (M_{i,j+1/2-}^{n+1,k} - M_{i,j-1/2+}^{n+1,k}), $$

where $Q_{i,j}^{n+1}$ is a collision term verifying (4.5), and where the last line of the above expression accounts for the topography source term. As in the one dimensional case, it is not possible to compute explicitly this update, instead we approximate it by a fixed point method

$$ (1 + \alpha) f_{i,j}^{n+1,k+1} = M_{i,j}^{n} + \alpha M_{i,j}^{n+1,k} - \sigma_{x} \xi_{1} (M_{i+1/2,j}^{n+1,k} - M_{i-1/2,j}^{n+1,k}) - \sigma_{y} \xi_{2} (M_{i,j+1/2}^{n+1,k} - M_{i,j-1/2}^{n+1,k}) \\ + \sigma_{x} (\xi_{1} - u_{i,j}^{n+1,k}) (M_{i+1/2-,j}^{n+1,k} - M_{i-1/2+,j}^{n+1,k}) + \sigma_{y} (\xi_{2} - v_{i,j}^{n+1,k}) (M_{i,j+1/2-}^{n+1,k} - M_{i,j-1/2+}^{n+1,k}), \quad (4.13) $$

where $\alpha > 0$ is a relaxation parameter, combined with the projection step

$$ M_{i,j}^{n+1,k+1} = f_{i,j}^{n+1,k+1} + Q_{i,j}^{n+1,k}. \quad (4.14) $$

At the macroscopic level, an iterative scheme is obtained by integrating (4.13) against $(1, \xi)$ which is implemented in practice. It reads

$$ (1 + \alpha) U_{i,j}^{n+1,k+1} = U_{i,j}^{n} + \alpha U_{i,j}^{n+1,k} - \sigma_{x} (F_{x,i+1/2-,j}^{n+1,k} - F_{x,i-1/2+,j}^{n+1,k}) - \sigma_{y} (F_{y,i,j+1/2-}^{n+1,k} - F_{y,i,j-1/2+}^{n+1,k}), \quad (4.15) $$

where we have the following definition of the numerical fluxes

$$ F_{x,i \pm 1/2, \mp, j}^{n+1,k} = \int_{\mathbb{R}^2} \xi_1 \begin{pmatrix} 1 \\ \xi_1 \\ \xi_2 \end{pmatrix} M_{i \pm 1/2, j}^{n+1,k} \mathrm{d}\xi_1 \mathrm{d}\xi_2 + \frac{g}{2} \begin{pmatrix} 0 \\ (h_{i,j}^{n+1,k})^2 - (h_{i \pm 1/2, \mp, j}^{n+1,k})^2 \\ 0 \end{pmatrix}, $$

$$ F_{y,i,j \pm 1/2 \mp}^{n+1,k} = \int_{\mathbb{R}^2} \xi_2 \begin{pmatrix} 1 \\ \xi_1 \\ \xi_2 \end{pmatrix} M_{i,j \pm 1/2}^{n+1,k} \mathrm{d}\xi_1 \mathrm{d}\xi_2 + \frac{g}{2} \begin{pmatrix} 0 \\ 0 \\ (h_{i,j}^{n+1,k})^2 - (h_{i,j \pm 1/2 \mp}^{n+1,k})^2 \end{pmatrix}. $$

The results obtained in Section 3 for the 1d Saint-Venant remain valid in the two dimensional setting since, the kinetic representation (4.4) of the Saint-Venant equations with Maxwellian (4.3) possesses a convex kinetic entropy $H$ given by (4.6).

The complete analysis of the 2d scheme is left for a forthcoming paper but a numerical experiment consisting of the two dimensional parabolic bowl [17] is presented in Section 5.3.

## 5. NUMERICAL SIMULATIONS

In this section we evaluate the behavior and the efficiency of the proposed kinetic schemes through numerical experiments. In Section 5.1, we focus on the fully implicit kinetic scheme over a flat bathymetry that corresponds to (2.1). Then, in Section 5.2, the iterative kinetic scheme with hydrostatic reconstruction is investigated numerically. In particular, we exhibit a test case where our iterative method succeeds to dissipate the total energy while its explicit version increases it. Finally in Section 5.3 a two dimensional test case is proposed. Unless specified otherwise, we use a CFL constant of 0.45 and the relaxation parameter $\alpha$ is set to one for the iterative kinetic schemes.

### 5.1. Fully implicit kinetic scheme in the flat bottom case

**Slow moving shock.** To assess the efficiency and interest of the implicit scheme (2.2), we perform a numerical test involving a Riemann problem with a slowly moving shock over a flat bottom. This configuration is achieved for a nearly transcritical flow where the material velocity $u$ is positive and satisfies $u - \sqrt{gh} \approx 0$ and $u + \sqrt{gh} \gg 1$.
$l_{x>0}U_R$ where we define

$$U_L = \begin{pmatrix} 2 \\ 1/2 \end{pmatrix}, \quad U_R = \begin{pmatrix} 1 \\ 1/2 \end{pmatrix}.$$

The gravity constant $g$ is set to 100 and the solution consists in a 1-rarefaction and a 2-shock. The iterative kinetic scheme uses the half-disk Maxwellian, and we choose the parameters $\alpha = 1$ and $\epsilon_{tol} = 10^{-9}$ for the stopping criterion. All the schemes use an explicit time step, and the results are given in Figure 3 at time $t = 0.025$. Three aspects have to be considered, namely the accuracy, the computational cost and the stability. Over Fig. 3, we see that in terms of efficiency both iterative and implicit kinetic schemes are at their disadvantage. Especially, the quadratic complexity of the fully implicit version results in a steeper slope of the efficiency curve. However this is only one part of the picture, and we know from Proposition 3.6 and Remark 3.7 that the iterative kinetic scheme (3.3) satisfies a discrete entropy inequality without restriction on the time step, assuming enough iterations are performed. Concretely the greater stability comes with a higher level of diffusion which is noticeable in the first two plots of Figure 3. This increased diffusion remains within acceptable margin, and is the price to pay to have better stability properties.

<figure><img src="image_7.png" /><figcaption>Figure 4: Parabolic bowl approximated by explicit and iterative kinetic schemes. First row: elevation and discharge at time 0.75, second row: convergence and efficiency curves. The stopping criteria used in the two kinetic iterative schemes combines the standard tolerance condition with tolerance $\epsilon = 10^{-9}$ and the entropy condition (3.32).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Slow moving shock approximated by various kinetic schemes, including explicit, implicit and iterative strategies. The initial condition is given by a Riemann data with discontinuity at position $x = 0$.</figcaption></figure>

Hence the maximum eigenvalue severely constrains the time step, however a small time step might not be necessary to accurately resolve the slow shock. We set the gravitational acceleration $g$ to 10 and the Riemann problem corresponds to the initial state $U_L = (1, 4.75)^T$ for $x < 0$ and $U_R = (2, 4.75)^T$ for $x > 0$. In Figure 2 we compare several schemes with an explicit time step $\Delta t_{\text{exp}}$ given by the CFL condition $\Delta t_{\text{exp}} \le 0.45 \frac{\Delta x}{\lambda_{\text{max}}}$, as well as the implicit kinetic scheme using a time step $\Delta t_{\text{imp}} = 10\Delta t_{\text{exp}}$. We also use the iterative scheme (3.3) with parameter $\alpha = 1$, and plot the results at time $t = 0.5$. We notice that in the discharge profile, an oscillation appears downwind of the left-going shock, which is quite pronounced for the explicit and iterative kinetic schemes, and less so for the fully implicit ones. As expected the implicit scheme using $\Delta t_{\text{imp}}$ strongly diffuses the fast traveling rarefaction. On the other hand the slow shock seems to be slightly less impacted by the large time steps, however it is still less diffused when using $\Delta t_{\text{exp}}$. Despite requiring ten times less iterations to reach the final time, the use of large time increments for the implicit kinetic scheme only results in around two percents faster computations compared to the explicit strategy which is due to the high quadratic cost of the implicit method. We believe that it is not possible to lower this cost when it comes to unconditionally stable methods, because the associated stencil has to cover the entire computational domain.

**Riemann problem.** We compare the fully implicit kinetic scheme and iterative kinetic scheme to explicit methods. The test case is given by the Riemann problem with initial data $U^0(x) = \mathbb{1}_{x<0}U_L + \mathbb{1}_{x>0}U_R$ where we define

$$U_L = \begin{pmatrix} 2 \\ 1/2 \end{pmatrix}, \quad U_R = \begin{pmatrix} 1 \\ 1/2 \end{pmatrix}.$$

The gravity constant $g$ is set to 100 and the solution consists in a 1-rarefaction and a 2-shock. The iterative kinetic scheme uses the half-disk Maxwellian, and we choose the parameters $\alpha = 1$ and $\varepsilon_{\text{tol}} = 10^{-9}$ for the stopping criterion. All the schemes use an explicit time step, and the results are given in Figure 3 at time $t = 0.025$. Three aspects have to be considered, namely the accuracy, the computational cost and the stability. In Figure 3, we see that in terms of efficiency both iterative and implicit kinetic schemes are at their disadvantage. Especially, the quadratic complexity of the fully implicit version results in a steeper slope of the efficiency curve. However this is only one part of the picture, and we know from Proposition 3.5 and Remark 3.6 that the iterative kinetic scheme (3.3) satisfies a discrete entropy inequality without restriction on the time step, assuming enough iterations are performed. Concretely the greater stability comes with a higher level of diffusion
## 5.2 Iterative kinetic scheme with hydrostatic reconstruction

**Parabolic bowl.** Now we consider the Thacker's testcase, also known as the parabolic bowl testcase, taken from [17]. The analytic solution in the domain $\Omega = (0, L)$ corresponds to

$$ egin{cases} h(t, x) = -\frac{h_0}{a} \left( \left[ \left(x - \frac{L}{2}\right) + \frac{1}{2} \cos(2Bt) \right]^2 - 1 \right) \mathbb{1}_{x \in W(t)} \\ u(t, x) = B \sin(2Bt) \mathbb{1}_{x \in W(t)} \\ z(x) = h_0 \left( \frac{1}{a^2} \left(x - \frac{L}{2}\right)^2 - 1 \right) \end{cases}, \quad W(t) = \frac{L}{2} - \frac{1}{2} \cos(2Bt) + (-a, a), $$

where we set $h_0 = 0.5$, $a = 1$, $L = 4$, $B = \frac{1}{2a}\sqrt{2gh_0}$ and $g = 10$. We plot the numerical solution at time $t = 0.75s$ in Figure 4. This testcase is relevant as it provides us with a non trivial analytical solution enabling to plot convergence curves, and it is known to be challenging numerically, as it presents a varying bottom together with an evolving wet/dry front and a discontinuous velocity profile. It is interesting to note that the different choice of Maxwellian used in the two iterative kinetic schemes has very little impact on the approximation. In both cases we obtain a convergence with first order accuracy, and unsurprisingly the numerical cost is higher than for fully explicit methods due to the number of sub-iterations required to update the solution. One should also note that the use of the half-disk Maxwellian is slightly more expensive than the simpler index Maxwellian. Besides, in this testcase the iterative kinetic scheme with index Maxwellian was always able to fulfill the entropy condition (3.32) after some iterations, which we only proved rigorously for the half-disk Maxwellian. Hence despite using the wrong Maxwellian, it seems that the iterative kinetic scheme in question still has better stability properties than fully explicit methods. This will be further corroborated with the next testcase.

<figure><img src="image_8.png" /><figcaption>Figure 5: Evolution of the relative total energy obtained for various explicit and iterative kinetic schemes.</figcaption></figure>

**Total energy dissipation.** To assess the good entropy stability property of our iterative kinetic scheme with hydrostatic reconstruction (3.22)-(3.24), we wish to compare it to its explicit counterpart in a testcase where this latter method fails to dissipate the entropy. To construct such a testcase,
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Comparing implicit, iterative and explicit kinetic solvers on a Riemann problem.</figcaption></figure>

which is noticeable in the first two plots of Figure 3. This increased diffusion remains within acceptable margin, and is the price to pay to have better stability properties.

**Total water height conservation.** We consider a double dam break problem over the domain $(-1, 1)$ with flat bottom, $g = 10$ and the following initial condition:

$$ \forall x \in (-1, 1), \quad h^0(x) = 2 \cdot 1_{|x|<0.2}, \quad q^0(x) = 0. $$

For sufficiently small times $0 \le t \le T_f$ the solution remains supported in $(-1, 1)$, hence we can enforce homogeneous Dirichlet boundary conditions on $h$ and $q$. In practice the simulation uses 512 cells and stops at $T_f = 0.05$. The relative total water height and total discharge are plotted in Figure 4 for the fully implicit kinetic scheme (2.1) and the iterative kinetic scheme (3.24). The latter achieves a deviation of less than $10^{-15}$ for both the total water height and total discharge. On the other hand, the maximum deviations for the fully implicit scheme are respectively of order $10^{-11}$ and $10^{-8}$. Hence when the conservation of mass and momentum is of critical importance, the iterative kinetic scheme should be preferred.
we can use the discrete entropy inequality verified by the explicit kinetic scheme established in [4], Theorem 3.6. It reads as

$$
\begin{aligned}
H(f_i^{n+1}, z_i) \le H(M_i^n, z_i) & - \sigma(\tilde{G}_{i+1/2-}^n - \tilde{G}_{i-1/2+}^n) \\
& - C_1 \sigma |\xi| (1_{\xi<0}(M_{i+1/2+} + M_{i+1/2-})(M_{i+1/2+} - M_{i+1/2-})^2 \\
& + 1_{\xi>0}(M_{i-1/2+} + M_{i-1/2-})(M_{i-1/2+} - M_{i-1/2-})^2) \\
& + C_2(\sigma \xi_{\max})^2 M_i ((M_i - M_{i+1/2-})^2 + (M_i - M_{i-1/2+})^2),
\end{aligned}
$$

with the kinetic entropy fluxes $\tilde{G}_{i+1/2-}^n$, $\tilde{G}_{i-1/2+}^n$ defined in (3.28) and with $C_1, C_2 > 0$ two constants. Since terms in factor of $C_1$ have negative sign, they play the role of a dissipation; on the other hand the ones in factor of $C_2$ are a positive error contribution. A possibility is thus to choose an initial condition such that the dissipative term cancels, but not the error one. In the dissipation, quantities of the form $M_{i\pm 1/2+} - M_{i\pm 1/2-}$ vanish when the free surface elevation $h+z$ and the velocity $u$ are constant across the mesh. At the same time, terms of the form $(M_i - M_{i\pm 1/2\mp})$ appearing in the positive error term vanish if the bathymetry is taken flat. Hence we consider an initial condition given by a flat free surface, a non-flat bathymetry and constant velocity which is nonzero (otherwise one would get a lake at rest steady state). More precisely we set the spatial domain to $\Omega = (0,1)$ with periodic boundary conditions, with $g=10$ and we have initially

$$ \forall x \in (0,1), \quad h^{\text{in}}(x) = -z(x), \quad u^{\text{in}}(x) = 1, \quad z(x) = -5 + \frac{1}{2} \left(1 + \cos\left(5\pi\left(x - \frac{1}{2}\right)\right)\right) 1_{|x - \frac{1}{2}| \le \frac{1}{5}} $$

The results can be seen in Figure 5, where we plot the time evolution of the total energy for various schemes over a mesh of $P=100$ cells.

Interestingly all the iterative methods manage to dissipate the total energy at each time step, even the scheme using the index Maxwellian (2.18), for which there is no proof of discrete entropy inequality. On the contrary, as expected the explicit kinetic scheme with half-disk Maxwellian increases the energy in the first few time steps, after what it decreases. The same goes for the explicit HLL scheme, and as a result these two explicit methods are not entropy stable for this testcase. For comparison we also added in dark blue the iterative kinetic scheme with $\alpha=0$ and whose sub-iterations stop as soon as the entropy condition (3.32) is verified. We can see that after some time this scheme becomes less dissipative than iterative kinetic methods using the standard tolerance condition; moreover this time roughly corresponds to the time at which the fully explicit schemes stop increasing the entropy.

## 5.3 Two dimensional simulation

We investigate the two dimensional iterative kinetic scheme with hydrostatic reconstruction given in Section 4.3. The exact solution in the domain $\Omega = (0, L)^2$ can be found in [17], and reads

$$
\left\{
\begin{array}{l}
h(t, x, y) = \frac{2\eta h_0}{a^2} \left( \left[ \left(x - \frac{L}{2}\right) \cos(\omega t) + \left(y - \frac{L}{2}\right) \sin(\omega t) \right] - z(x, y) \right)_+ \\
u(t, x, y) = -\eta \omega \sin(\omega t) \mathbf{1}_{h(t,x,y)>0} \\
v(t, x, y) = \eta \omega \cos(\omega t) \mathbf{1}_{h(t,x,y)>0}
\end{array}
\right.
, \quad z(x, y) = -h_0 \left(1 - \frac{r^2}{a^2}\right),
$$

with $r = \sqrt{(x-L/2)^2 + (y-L/2)^2}$, $L=1$, $h_0 = 0.1$, $a = 0.25$, $\eta = 0.2$, $\omega = \sqrt{2gh_0}/a$ and $g=4$. The initial condition is obtained by setting $t=0$ in the above expressions. The results are obtained with the 2D version of the iterative kinetic scheme (4.15) with the Maxwellian defined by (4.2)-(4.3) and are displayed in Figures 6 and 7. We see that when increasing the tolerance value to $\epsilon_{\text{tol}} = 10^{-5}$, the experimental order of convergence of the iterative scheme decreases, which illustrates that a large tolerance error prevents the sub-iterations to converge to the implicit update. On the other hand, the smaller $\epsilon_{\text{tol}}$ is, the more iterations are needed to reach the stopping criteria which translates to an increase in computational time.
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Total water height and discharge evolution. The light blue line stands for the fully implicit kinetic scheme (2.1) and the dark dotted line for the iterative kinetic scheme (3.24).</figcaption></figure>

## 5.2. Iterative kinetic scheme with hydrostatic reconstruction

**Parabolic bowl.** Now we consider the Thacker's test case, also known as the parabolic bowl test case, taken from [17]. The analytic solution in the domain $\Omega = (0, L)$ corresponds to

$$ \begin{cases} h(t, x) = -\frac{h_0}{a} \left( \left[ \left(x - \frac{L}{2}\right) + \frac{1}{2} \cos(2Bt) \right]^2 - 1 \right) \mathbf{1}_{x \in W(t)} \\ u(t, x) = B \sin(2Bt) \mathbf{1}_{x \in W(t)} \\ z(x) = h_0 \left( \frac{1}{a^2} \left(x - \frac{L}{2}\right)^2 - 1 \right), \end{cases} \quad W(t) = \frac{L}{2} - \frac{1}{2} \cos(2Bt) + (-a, a), $$

where we set $h_0 = 0.5$, $a = 1$, $L = 4$, $B = \frac{1}{2a}\sqrt{2gh_0}$ and $g = 10$. We plot the numerical solution at time $t = 0.75$ in Figure 5. This test case is relevant as it provides us with a non trivial analytical solution enabling to plot convergence curves, and it is known to be challenging numerically, as it presents a varying bottom together with an evolving wet/dry front and a discontinuous velocity profile. It is interesting to note that the different choice of Maxwellian used in the two iterative kinetic schemes has very little impact on the approximation. In both cases we obtain a convergence with first order accuracy, and unsurprisingly the numerical cost is higher than for fully explicit methods due to the number of sub-iterations required to update the solution. One should also note that the use of the half-disk Maxwellian is slightly more expensive than the simpler index Maxwellian. Besides, in this test case the iterative kinetic scheme with index Maxwellian was always able to fulfill the entropy condition (3.32) after some iterations, which we only proved rigorously for the half-disk Maxwellian. Hence despite using the wrong Maxwellian, it seems that the iterative kinetic scheme in question still has better stability properties than fully explicit methods. This will be further corroborated with the next test case.

**Total energy dissipation.** To assess the good entropy stability property of our iterative kinetic scheme with hydrostatic reconstruction (3.22)–(3.24), we wish to compare it to its explicit counterpart in a test case where this latter method fails to dissipate the entropy. To construct such a test case, we can use the discrete entropy
<figure><img src="image_9.png" /><figcaption>Figure 6: Convergence and efficiency curves obtained with the 2D parabolic bowl testcase with final time $t = 4\pi/\omega$. Different tolerances $\varepsilon$ are compared for the iterative kinetic scheme. A CFL constant of 1/2 was used, except for the case $\varepsilon = 10^{-9}$ where we set it to 1/10.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: Numerical approximation of the 2D parabolic bowl using the iterative kinetic scheme with $\varepsilon = 10^{-7}$ over a 100 × 100 mesh. From left to right and top to bottom: initial condition, approximation at time $t = 1.1708$, $t = 2.3416$ and $t = 3.5124$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Parabolic bowl approximated by explicit and iterative kinetic schemes. First row: elevation and discharge at time 0.75, second row: convergence and efficiency curves, third row: conservation of the total water height up to machine precision. The stopping criteria used in the two kinetic iterative schemes combines the standard tolerance condition with tolerance $\epsilon = 10^{-9}$ and the entropy condition (3.32).</figcaption></figure>
Acknowledgments

The authors wish to express their warmest thanks to François Bouchut for many fruitful discussions.
Antonin Leprevost and Bilal Al Taki have contributed preliminary versions of the work.

References

[1] Sebastien Allgeyer, Marie-Odile Bristeau, David Froger, Raouf Hamouda, V. Jauzein, Anne Mangeney, Jacques Sainte-Marie, Fabien Souillé, and Martin Valée, *Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface*, ESAIM, Math. Model. Numer. Anal. **53** (2019), no. 6, 1981–2024 (English).

[2] K.R. Arun, A.J. Das Gupta, and S. Samantaray, *Analysis of an asymptotic preserving low mach number accurate imex-rk scheme for the wave equation system*, Applied Mathematics and Computation **411** (2021), 126469.

[3] Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and Benoît Perthame, *A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows*, SIAM J. Sci. Comput. **25** (2004), no. 6, 2050–2065 (English).

[4] Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, and Jacques Sainte-Marie, *Kinetic entropy inequality and hydrostatic reconstruction scheme for the Saint-Venant system*, Math. Comput. **85** (2016), no. 302, 2815–2837 (English).

[5] Emmanuel Audusse and Marie-Odile Bristeau, *A 2d Well-balanced Positivity Preserving Second Order Scheme for Shallow Water Flows on Unstructured Meshes*, Research Report RR-5260, INRIA, 2004.

[6] Emmanuel Audusse and Marie-Odile Bristeau, *A well-balanced positivity preserving “second-order” scheme for shallow water flows on unstructured meshes*, J. Comput. Phys. **206** (2005), no. 1, 311–333 (English).

[7] F. Berthelin and F. Bouchut, *Relaxation to isentropic gas dynamics for a BGK system with single kinetic entropy*, Methods Appl. Anal. **9** (2002), no. 2, 313–327 (English).

[8] Georgij Bispen, Koottungal Revi Arun, Maria Lukàcovà-Medvid'ovà, and Sebastian Noelle, *IMEX Large Time Step Finite Volume Methods for Low Froude Number Shallow Water Flows*, Communications in Computational Physics **16** (2014), 307–347.

[9] F. Bouchut, *Construction of BGK models with a family of kinetic entropies for a given system of conservation laws*, J. Stat. Phys. **95** (1999), no. 1-2, 113–170 (English).

[10] François Bouchut, *Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources*, Front. Math., Basel: Birkhäuser, 2004 (English).

[11] François Bouchut and Xavier Lhébrard, *Convergence of the kinetic hydrostatic reconstruction scheme for the Saint Venant system with topography*, Math. Comput. **90** (2021), no. 329, 1119–1153 (English).

[12] Bouchut, F., *Entropy satisfying flux vector splittings and kinetic BGK models*, Numer. Math. **94** (2003), no. 4, 623–672 (English).

[13] Marie-Odile Bristeau and Benoît Coussin, *Boundary Conditions for the Shallow Water Equations solved by Kinetic Schemes*, Research Report RR-4282, INRIA, 2001, Projet M3N.

[14] Marie-Odile Bristeau, Nicole Goutal, and Jacques Sainte-Marie, *Numerical simulations of a non-hydrostatic shallow water model*, Comput. Fluids **47** (2011), no. 1, 51–64 (English).
<figure><img src="image_8.png" /><figcaption>FIGURE 6. Evolution of the relative total energy obtained for various explicit and iterative kinetic schemes.</figcaption></figure>

inequality verified by the explicit kinetic scheme established in [6], Theorem 3.6. It reads as

$$ H(f_i^{n+1}, z_i) \le H(M_i^n, z_i) - \sigma(\tilde{G}_{i+1/2-}^n - \tilde{G}_{i-1/2+}^n) \\ - C_1\sigma|\xi|(1_{\xi<0}(M_{i+1/2+} + M_{i+1/2-})(M_{i+1/2+} - M_{i+1/2-})^2 \\ + 1_{\xi>0}(M_{i-1/2+} + M_{i-1/2-})(M_{i-1/2+} - M_{i-1/2-})^2) \\ + C_2(\sigma\xi_{\max})^2 M_i((M_i - M_{i+1/2-})^2 + (M_i - M_{i-1/2+})^2), $$

with the kinetic entropy fluxes $\tilde{G}_{i+1/2-}^n, \tilde{G}_{i-1/2+}^n$ defined in (3.28) and with $C_1, C_2 > 0$ two constants. Since terms in factor of $C_1$ have negative sign, they play the role of a dissipation; on the other hand the ones in factor of $C_2$ are a positive error contribution. A possibility is thus to choose an initial condition such that the dissipative term cancels, but not the error one. In the dissipation, quantities of the form $M_{i\pm 1/2+} - M_{i\pm 1/2-}$ vanish when the free surface elevation $h+z$ and the velocity $u$ are constant across the mesh. At the same time, terms of the form $(M_i - M_{i\pm 1/2\mp})$ appearing in the positive error term vanish if the bathymetry is taken flat. Hence we consider an initial condition given by a flat free surface, a non flat bathymetry and constant velocity which is nonzero (otherwise one would get a lake at rest steady state). More precisely we set the spatial domain to $\Omega = (0, 1)$ with periodic boundary conditions, with $g = 10$ and we have initially

$$ \forall x \in (0, 1), \quad h^{\text{in}}(x) = -z(x), \quad u^{\text{in}}(x) = 1, \quad z(x) = -5 + \frac{1}{2}\left(1 + \cos\left(5\pi\left(x - \frac{1}{2}\right)\right)\right)1_{|x - \frac{1}{2}| \le \frac{1}{5}}. $$

The results can be seen in Figure 6, where we plot the time evolution of the total energy for various schemes over a mesh of $P = 100$ cells. Interestingly all the iterative methods manage to dissipate the total energy at each time step, even the scheme using the index Maxwellian (2.18), for which there is no proof of discrete
[15] F. Coron and B. Perthame, *Numerical passage from kinetic to fluid equations*, SIAM J. Numer. Anal. **28** (1991), no. 1, 26–42 (English).

[16] David Coulette, Emmanuel Franck, Philippe Helluy, Michel Mehrenberger, and Laurent Navoret, *High-order implicit palindromic discontinuous galerkin method for kinetic-relaxation approximation*, Computers & Fluids **190** (2019), 485–502.

[17] Olivier Delestre, Carine Lucas, Pierre-Antoine Ksinant, Frédéric Darboux, Christian Laguerre, T.-N.-Tuoi Vo, François James, and Stéphane Cordier, *SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies*, Int. J. Numer. Methods Fluids **72** (2013), no. 3, 269–300 (English).

[18] J.-F. Gerbeau and B. Perthame, *Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation*, Discrete Contin. Dyn. Syst., Ser. B **1** (2001), no. 1, 89–102 (English).

[19] Laurent Gosse, *Computing qualitatively correct approximations of balance laws. Exponential-fit, well-balanced and asymptotic-preserving*, SIMAI Springer Ser., vol. 2, Milano: Springer, 2013 (English).

[20] N. Goutal and J. Sainte-Marie, *A kinetic interpretation of the section-averaged Saint-Venant system for natural river hydraulics*, Int. J. Numer. Methods Fluids **67** (2011), no. 7, 914–938 (English).

[21] J. Haack, S. Jin, and J.-G. Liu, *An all-speed asymptotic-preserving method for the isentropic euler and navier-stokes equations*, Communications in Computational Physics **12**(4) (2012), 955–980.

[22] Shi Jin and Lorenzo Pareschi, *Asymptotic-preserving (ap) schemes for multiscale kinetic equations: a unified approach*, Hyperbolic Problems: Theory, Numerics, Applications (Basel) (Heinrich Freistühler and Gerald Warnecke, eds.), Birkhäuser Basel, 2001, pp. 573–582.

[23] B. Perthame, *Boltzmann type schemes for gas dynamics and the entropy property*, SIAM J. Numer. Anal. **27** (1990), no. 6, 1405–1421 (English).

[24] B. Perthame and C. Simeoni, *A kinetic scheme for the Saint-Venant system with a source term*, Calcolo **38** (2001), no. 4, 201–231 (English).

[25] Benoît Perthame, *Kinetic formulation of conservation laws*, Oxf. Lect. Ser. Math. Appl., vol. 21, Oxford: Oxford University Press, 2002 (English).

[26] Perthame, B., *Second-order Boltzmann schemes for compressible Euler equations in one and two space dimensions*, SIAM J. Numer. Anal. **29** (1992), no. 1, 1–19 (English).

[27] Saint-Venant, *Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit.*, C. R. Acad. Sci., Paris **73** (1871), 147–154 (French).

[28] Andrea Thomann, Gabriella Puppo, and Christian Klingenberg, *An all speed second order well-balanced imex relaxation scheme for the euler equations with gravity*, Journal of Computational Physics **420** (2020), 109723.

[29] Yulong Xing and Chi-Wang Shu, *A survey of high order schemes for the shallow water equations*, J. Math. Study **47** (2014), no. 3, 221–249 (English).

## A Expression of the numerical fluxes

The optimal choice for the Maxwellian is given by (1.4). Unfortunately the explicit expression for the numerical fluxes appearing in (2.20) is hardly possible with the choice (1.4) and the use of approximate quadrature formula for the integrals in (1.4) will degrade the accuracy of the scheme and increase the
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Water height and discharge $L^2$ errors obtained for the 2d parabolic bowl test case at final time $t = 2\pi/\omega$ with the HLL scheme, explicit kinetic scheme and iterative kinetic scheme.</figcaption></figure>

entropy inequality. On the contrary, as expected the explicit kinetic scheme with half-disk Maxwellian increases the energy in the first few time steps, after what it decreases. The same goes for the explicit HLL scheme, and as a result these two explicit methods are not entropy stable for this test case. For comparison we also added in dark blue the iterative kinetic scheme with $\alpha = 0$ and whose sub-iterations stop as soon as the entropy condition (3.32) is verified. We can see that after some time this scheme becomes less dissipative than iterative kinetic methods using the standard tolerance condition; moreover this time roughly corresponds to the time at which the fully explicit schemes stop increasing the entropy.

## 5.3. Two dimensional simulation

We investigate the two dimensional iterative kinetic scheme with hydrostatic reconstruction given in Section 4.3. The exact solution in the domain $\Omega = (0, L)^2$ can be found in [17], and reads

$$ egin{cases} h(t, x, y) = \frac{2\eta h_0}{a^2} \left( \left[ \left(x - \frac{L}{2}\right) \cos(\omega t) + \left(y - \frac{L}{2}\right) \sin(\omega t) \right] - z(x, y) \right)_+ \\ u(t, x, y) = -\eta\omega \sin(\omega t) \mathbf{1}_{h(t,x,y)>0} \\ v(t, x, y) = \eta\omega \cos(\omega t) \mathbf{1}_{h(t,x,y)>0}, \end{cases} \qquad z(x, y) = -h_0 \left( 1 - \frac{r^2}{a^2} \right), $$

with $r = \sqrt{(x - L/2)^2 + (y - L/2)^2}$, $L = 1$, $h_0 = 0.1$, $a = 0.25$, $\eta = 0.2$, $\omega = \sqrt{2gh_0}/a$ and $g = 10$. Setting $t = 0$ in the above expressions yields the initial condition, and the numerical approximation is computed with the 2d version of the iterative kinetic scheme (4.15) using the Maxwellian defined by (4.2) and (4.3), with $\alpha = 3$, $\varepsilon = 10^{-9}$ and a CFL constant of 0.2. The $L^2$ errors for the water height and the discharge are compared to that of the explicit HLL and explicit kinetic schemes in Figure 7, where we see similar results. A surface plot of the elevation and bathymetry is also provided in Figure 8 at different times.

## 6. CONCLUSION

In summary, using a kinetic solver, we built implicit schemes for the one- and two-dimensional Saint-Venant systems that are structure preserving and in particular entropy dissipating, therefore improving upon the previously investigated explicit scheme [6]. Indeed we showed that an implicit version of the scheme offers stronger stability results: in the flat bottom case we got rid of the CFL condition while retaining the positivity and the entropy dissipation properties; whereas in the varying bottom case we showed that it is possible to fix the creation of entropy induced by the hydrostatic reconstruction technique when using a fixed point sub-iteration
computational costs. Hence, we choose M defined by the first expression in (2.17) and relation (2.20) becomes

$$ U_i^{\text{int}} = \frac{1}{2\sqrt{3}} \left( \int_{\min(0, a_j^n)}^{\min(0, b_j^n)} \sum_{j=i}^{P} \sqrt{\frac{2h_j^n}{g}} \binom{1}{\xi} \frac{(-\sigma\xi)^{j-i}}{(1-\sigma\xi)^{j-i+1}} d\xi \right. \\ \left. \qquad + \int_{\max(0, a_j^n)}^{\max(0, b_j^n)} \sum_{j=1}^{i} \sqrt{\frac{2h_j^n}{g}} \binom{1}{\xi} \frac{(\sigma\xi)^{i-j}}{(1+\sigma\xi)^{i-j+1}} d\xi \right), $$

with $a_j^n = u_j^n - \sqrt{3}c_j^n$ and $b_j^n = u_j^n + \sqrt{3}c_j^n$. The expressions of $h_i^{\text{int}}$ and $(hu)_i^{\text{int}}$ are given by

$$ \begin{array}{l} h_i^{\text{int}} = \frac{1}{2\sqrt{3}} \left( \sum_{j=i}^{P} \sqrt{\frac{2h_j^n}{g}} \underbrace{\int_{\min(0, a_j^n)}^{\min(0, b_j^n)} \frac{(-\sigma\xi)^{j-i}}{(1-\sigma\xi)^{j-i+1}} d\xi}_{(Ah)_{i,j}} \right. \\ \left. \qquad + \sum_{j=1}^{i} \sqrt{\frac{2h_j^n}{g}} \underbrace{\int_{\max(0, a_j^n)}^{\max(0, b_j^n)} \frac{(\sigma\xi)^{i-j}}{(1+\sigma\xi)^{i-j+1}} d\xi}_{(Bh)_{i,j}} \right) \\[1.5em] (A.1) \qquad \qquad \qquad \qquad \qquad \qquad \qquad (hu)_i^{\text{int}} = \frac{1}{2\sqrt{3}} \left( \sum_{j=i}^{P} -\frac{1}{\sigma} \sqrt{\frac{2h_j^n}{g}} \underbrace{\int_{\min(0, a_j^n)}^{\min(0, b_j^n)} \frac{(-\sigma\xi)^{j-i+1}}{(1-\sigma\xi)^{j-i+1}} d\xi}_{(Ahu)_{i,j}} \right) \\[1.5em] (A.2) \qquad \qquad \qquad \qquad \qquad \qquad \qquad + \sum_{j=1}^{i} \frac{1}{\sigma} \sqrt{\frac{2h_j^n}{g}} \underbrace{\int_{\max(0, a_j^n)}^{\max(0, b_j^n)} \frac{(\sigma\xi)^{i-j+1}}{(1+\sigma\xi)^{i-j+1}} d\xi}_{(Bhu)_{i,j}} \end{array} $$

Now we need to compute analytically the integrals of both expressions using the following lemmas.

**Lemma A.1** If we denote $y = 1 - \frac{1}{1+x}$ for all $x \in \mathbb{R}\setminus\{-1\}$ and $C \in \mathbb{R}$ we have the following primitive:

$$ \int \frac{x^k}{(1+x)^{k+1}} dx = \ln(|1+x|) - \sum_{l=1}^{k} \frac{y^l}{l} + C. $$

**Lemma A.2** Using the same notation as in the previous lemma, we have

$$ \int \frac{x^k}{(1+x)^k} dx = -k \ln(|1+x|) + x + \sum_{l=1}^{k-1} l \frac{y^{k-l}}{k-l} + C'. $$

**Proof of Lemma A.1.** We have

$$ I = \int \frac{x^k}{(1+x)^{k+1}} dx = \int \frac{x^k}{(1+x)^k} \frac{1}{1+x} dx = \int (1 - \frac{x}{1+x})^k \frac{1}{1+x} dx $$

We pose $y = 1 - \frac{1}{1+x}$

$$ I = \int y^k (1-y) \frac{dy}{(1-y)^2} = \int \frac{y^k - 1}{1-y} + \frac{1}{1-y} dy $$

Now we use the formula $y^k - 1 = (y-1)(y^{k-1} + y^{k-2} + \dots + y + 1)$. And we obtain

$$ I = - \int \sum_{l=0}^{k-1} y^l dy - \ln(|1-y|) + C \qquad C \in \mathbb{R} $$
<figure><img src="image_10.png" /><figcaption>FIGURE 8. Numerical approximation of the 2d parabolic bowl using the iterative kinetic scheme with $\epsilon = 10^{-7}$ over a 100 × 100 mesh. From left to right and top to bottom: initial condition, approximation at time $t = 2\pi/3\omega$, $t = 4\pi/3\omega$ and $t = 2\pi/\omega$.</figcaption></figure>

process that approximates the implicit update. The benefits of our approach were then illustrated through numerical experiments.

In general the choice of implicit solvers comes with an increase in computational cost compared to explicit schemes. Yet in some regimes with multiscale behavior, such as the one corresponding to low Froude numbers, it can make sense efficiency-wise to use an implicit scheme with large time steps. In this regard, the ability to explicitly write the update of our fully implicit scheme when the bathymetry is flat makes it competitive compared to other implicit strategies resolved by a direct method. However when the bathymetry is not flat, this implicit kinetic scheme cannot be used, instead we can consider the iterative kinetic scheme with hydrostatic reconstruction, which should be more relevant for real-life applications.

Future works include an improvement of the hydrostatic reconstruction that would allow to preserve a broader family of steady states. This would be achieved by modifying the velocity in addition to the water height at each interface. Another point to investigate is whether a semi-implicit time discretization combined with a kinetic solver can be used in order to reduce the amount of diffusion and improve the efficiency while keeping a fully discrete entropy inequality.

## ACKNOWLEDGMENTS

The authors wish to express their warmest thanks to François Bouchut for many fruitful discussions. Antonin Leprevost and Bilal Al Taki have contributed preliminary versions of the work.

## FUNDING

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

## DATA AVAILABILITY STATEMENT

The research data associated with this article are included in the article.
$$ = \ln(|1 + x|) - \sum_{l=1}^{k} \frac{y^l}{l} + C' \qquad C' \in \mathbb{R} $$

□

**Proof of Lemma A.2.** We already have denoted $y = \frac{x}{1+x} = 1 - \frac{1}{1+x}$

$$ I = \int \left(\frac{x}{1+x}\right)^k dx = \int \frac{y^k dy}{(1-y)^2} = \int \left(\frac{y^k - 1}{(1-y)^2} + \frac{1}{(1-y)^2}\right) dy $$

where the formula $y^k - 1 = (y-1)(y^{k-1} + y^{k-2} + \dots + y + 1)$ has been used. Hence

$$ \begin{aligned} I &= - \int \sum_{l=0}^{k-1} \frac{y^l}{1-y} dy + x + C = - \int \sum_{l=0}^{k-1} \frac{y^l - 1}{1-y} dy - \int \frac{1}{1-y} \sum_{l=0}^{k-1} dy + x + C \\ &= \int \sum_{l=1}^{k-1} \frac{y^l - 1}{y-1} dy + k \ln(|1-y|) + x + C' = \int \sum_{l=1}^{k-1} \sum_{p=0}^{l-1} y^p dy - k \ln(|1+x|) + x + C' \\ &= \sum_{l=1}^{k-1} l \int y^{k-1-l} dy - k \ln(|1+x|) + x + C' = \sum_{l=1}^{k-1} l \frac{y^{k-l}}{k-l} - k \ln(|1+x|) + x + C'', \end{aligned} $$

with $(C, C', C'') \in \mathbb{R}^3$. □

We are now able to compute the quantities $Ah_{i,j}, Bh_{i,j}, Ahu_{i,j}, Bhu_{i,j}$

$$ \begin{aligned} Ah_{i,j} &= \int_{\min(0, a_j^n)}^{\min(0, b_j^n)} \frac{(-\sigma\xi)^{j-i}}{(1-\sigma\xi)^{j-i+1}} d\xi = -\frac{1}{\sigma} \int_{-\min(0, a_j^n)\sigma}^{-\min(0, b_j^n)\sigma} \frac{(x)^{j-i}}{(1+x)^{j-i+1}} dx \\ (A.3) \qquad &= \frac{1}{\sigma} \left[ \ln(|1+x|) - \sum_{l=1}^{j-i} \frac{y^l}{l} \right]_{-\min(0, b_j^n)\sigma}^{-\min(0, a_j^n)\sigma}. \end{aligned} $$

$$ \begin{aligned} Bh_{i,j} &= \int_{\max(0, a_j^n)}^{\max(0, b_j^n)} \frac{(\sigma\xi)^{i-j}}{(1+\sigma\xi)^{i-j+1}} d\xi = \frac{1}{\sigma} \int_{\max(0, a_j^n)\sigma}^{\max(0, b_j^n)\sigma} \frac{(x)^{i-j}}{(1+x)^{i-j+1}} dx \\ (A.4) \qquad &= \frac{1}{\sigma} \left[ \ln(|1+x|) - \sum_{l=1}^{i-j} \frac{y^l}{l} \right]_{\max(0, a_j^n)\sigma}^{\max(0, b_j^n)\sigma} \end{aligned} $$

And similarly we obtain the formulas for *Ahu* and *Bhu* under the form

$$ (A.5) \qquad Ahu_{i,j} = \frac{1}{\sigma} \left[ -(j-i+1) \ln(|1+x|) + x + \sum_{l=1}^{j-i} \frac{y^{j-i+1-l}}{j-i+1-l} \right]_{-\min(0, b_j^n)\sigma}^{-\min(0, a_j^n)\sigma} $$

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

To conclude this paragraph, we give the final expression of $U_i^{\text{int}}$

$$ (A.7) \qquad h_i^{\text{int}} = \frac{1}{2\sigma\sqrt{3}} \left( \sum_{j=i}^{P} \sqrt{\frac{2h_j^n}{g}} \left[ \ln(|1+x|) - \sum_{l=1}^{j-i} \frac{y^l}{l} \right]_{-\min(0, b_j^n)\sigma}^{-\min(0, a_j^n)\sigma} \right) $$

$$ (A.8) \qquad + \sum_{j=1}^{i} \sqrt{\frac{2h_j^n}{g}} \left[ \ln(|1+x|) - \sum_{l=1}^{i-j} \frac{y^l}{l} \right]_{\max(0, a_j^n)\sigma}^{\max(0, b_j^n)\sigma} $$

$$ (A.9) \qquad (hu)_i^{\text{int}} = \frac{1}{2\sigma^2\sqrt{3}} \left( -\sum_{j=i}^{P} \sqrt{\frac{2h_j^n}{g}} \left[ -(j-i+1)\ln(|1+x|) + x + \sum_{k=1}^{j-i} (j-i+1-k)\frac{y^k}{k} \right]_{-\min(0, b_j^n)\sigma}^{-\min(0, a_j^n)\sigma} \right) $$
REFERENCES

[1] S. Allgeyer, M.O. Bristeau, D. Froger, R. Hamouda, V. Jauzein, A. Mangeney, J. Sainte-Marie, F. Souillé and M. Valée, Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface. *ESAIM: Math. Model. Numer. Anal.* **53** (2019) 1981-2024.

[2] K. Arun, A. Das Gupta and S. Samantaray, Analysis of an asymptotic preserving low mach number accurate imex-rk scheme for the wave equation system. *Appl. Math. Comput.* **411** (2021) 126469. https://www.sciencedirect.com/science/article/pii/S0096300321005580.

[3] E. Audusse and M.O. Bristeau, A 2d well-balanced positivity preserving second order scheme for shallow water flows on unstructured meshes. Research Report RR-5260, INRIA (2004). https://inria.hal.science/inria-00070738.

[4] E. Audusse and M.O. Bristeau, A well-balanced positivity preserving “second-order” scheme for shallow water flows on unstructured meshes. *J. Comput. Phys.* **206** (2005) 311–333.

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

[6] E. Audusse, F. Bouchut, M.O. Bristeau and J. Sainte-Marie, Kinetic entropy inequality and hydrostatic reconstruction scheme for the Saint-Venant system. *Math. Comput.* **85** (2016) 2815–2837.

[7] F. Berthelin and F. Bouchut, Relaxation to isentropic gas dynamics for a BGK system with single kinetic entropy. *Methods Appl. Anal.* **9** (2002) 313–327.

[8] G. Bispen, K.R. Arun, M. Lukàcovà-Medvid'ovà and S. Noelle, IMEX large time step finite volume methods for low froude number shallow water flows. *Commun. Comput. Phys.* **16** (2014) 307–347.

[9] F. Bouchut, Construction of BGK models with a family of kinetic entropies for a given system of conservation laws. *J. Stat. Phys.* **95** (1999) 113–170.

[10] F. Bouchut, Entropy satisfying flux vector splittings and kinetic BGK models. *Numer. Math.* **94** (2003) 623–672.

[11] F. Bouchut, Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources. Springer Science & Business Media (2004).

[12] F. Bouchut and X. Lhébrard, Convergence of the kinetic hydrostatic reconstruction scheme for the Saint-Venant system with topography. *Math. Comput.* **90** (2021) 1119–1153. https://hal-upec-upem.archives-ouvertes.fr/hal-01515256v3/file/kin-hydrostoj conv.pdf.

[13] M.O. Bristeau and B. Coussin, Boundary conditions for the shallow water equations solved by kinetic schemes. Research Report RR-4282, INRIA (2001). https://hal.inria.fr/inria-00072305.

[14] M.O. Bristeau, N. Goutal and J. Sainte-Marie, Numerical simulations of a non-hydrostatic shallow water model. *Comput. Fluids* **47** (2011) 51–64. www.sciencedirect.com/science/article/pii/S0045793011000727.

[15] F. Coron and B. Perthame, Numerical passage from kinetic to fluid equations. *SIAM J. Numer. Anal.* **28** (1991) 26–42.

[16] D. Coulette, E. Franck, P. Helluy, M. Mehrenberger and L. Navoret, High-order implicit palindromic discontinuous galerkin method for kinetic-relaxation approximation. *Comput. Fluids* **190** (2019) 485–502. https://www.sciencedirect.com/science/article/pii/S0045793019301835.

[17] O. Delestre, C. Lucas, P.A. Ksinant, F. Darboux, C. Laguerre, T.N.T. Vo, F. James and S. Cordier, SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies. *Int. J. Numer. Methods Fluids* **72** (2013) 269–300.

[18] J.F. Gerbeau and B. Perthame, Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation. *Discrete Contin. Dyn. Syst. Ser. B* **1** (2001) 89–102.

[19] I. Gómez-Bueno, S. Boscarino, M. Castro, C. Parés and G. Russo, Implicit and semi-implicit well-balanced finite-volume methods for systems of balance laws. *Appl. Numer. Math.* **184** (2023) 18–48.

[20] L. Gosse, Computing Qualitatively Correct Approximations of Balance Laws. Exponential-fit, Well-balanced and Asymptotic-preserving. *SIMAI Springer Series* Vol. 2. Springer, Milano (2013).

[21] N. Goutal and J. Sainte-Marie, A kinetic interpretation of the section-averaged Saint-Venant system for natural river hydraulics. *Int. J. Numer. Methods Fluids* **67** (2011) 914–938.

[22] J. Haack, S. Jin and J.G. Liu, An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations. *Commun. Comput. Phys.* **12** (2012) 955–980.

[23] S. Jin and L. Pareschi, Asymptotic-preserving (AP) schemes for multiscale kinetic equations: a unified approach, in *Hyperbolic Problems: Theory, Numerics, Applications*, edited by H. Freistühler and G. Warnecke. Birkhäuser Basel, Basel (2001) 573–582.

[24] B. Perthame, Boltzmann type schemes for gas dynamics and the entropy property. *SIAM J. Numer. Anal.* **27** (1990) 1405–1421.
$$ (A.10) \qquad + \sum_{j=1}^{i} \sqrt{\frac{2h_j^n}{g}} \left[ -(i-j+1) \ln(|1+x|) + x + \sum_{k=1}^{i-j} (i-j+1-k) \frac{y^k}{k} \right]_{\max(0, a_j^n)\sigma}^{\max(0, b_j^n)\sigma} $$

## B Computations of the fluxes involving the boundary conditions

We assume the ghost quantities $U_0^{n+1}$ and $U_{P+1}^{n+1}$ at time $t^{n+1}$ to be known. The exterior contribution given in (2.19) also writes

$$ U_i^{\text{ext}} = \int_{\mathbb{R}^{-}} \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_{\mathbb{R}^{+}} \left(\frac{1}{\xi}\right) \frac{(\sigma\xi)^i}{(1+\sigma\xi)^i} M_0^{n+1} d\xi. $$

Using computations similar to what has been proposed in Appendix A, we get

$$ U_i^{\text{ext}} = \frac{1}{2\sqrt{3}} \left[ \sqrt{\frac{2h_{P+1}^{n+1}}{g}} \int_{\min(0, a_{P+1}^{n+1})}^{\min(0, b_{P+1}^{n+1})} \left(\frac{1}{\xi}\right) \frac{(-\sigma\xi)^{P-i+1}}{(1-\sigma\xi)^{P-i+1}} d\xi \right. \\ \left. \qquad + \sqrt{\frac{2h_0^{n+1}}{g}} \int_{\max(0, a_0^{n+1})}^{\max(0, b_0^{n+1})} \left(\frac{1}{\xi}\right) \frac{(\sigma\xi)^i}{(1+\sigma\xi)^i} d\xi \right], $$

or equivalently

$$ h_i^{\text{ext}} = \frac{1}{2\sqrt{3}} \left[ \sqrt{\frac{2h_{P+1}^{n+1}}{g}} \int_{\min(0, a_{P+1}^{n+1})}^{\min(0, b_{P+1}^{n+1})} \frac{(-\sigma\xi)^{P-i+1}}{(1-\sigma\xi)^{P-i+1}} d\xi \right. \\ \left. \qquad + \sqrt{\frac{2h_0^{n+1}}{g}} \int_{\max(0, a_0^{n+1})}^{\max(0, b_0^{n+1})} \frac{(\sigma\xi)^i}{(1+\sigma\xi)^i} d\xi \right], $$

$$ (hu)_i^{\text{ext}} = \frac{1}{2\sqrt{3}} \left[ \sqrt{\frac{2h_{P+1}^{n+1}}{g}} \int_{\min(0, a_{P+1}^{n+1})}^{\min(0, b_{P+1}^{n+1})} \xi \frac{(-\sigma\xi)^{P-i+1}}{(1-\sigma\xi)^{P-i+1}} d\xi \right. \\ \left. \qquad + \sqrt{\frac{2h_0^{n+1}}{g}} \int_{\max(0, a_0^{n+1})}^{\max(0, b_0^{n+1})} \xi \frac{(\sigma\xi)^i}{(1+\sigma\xi)^i} d\xi \right]. $$

As explained in Section 2.4, in practice we will replace the unknown values of $U_0^{n+1}$, $U_{P+1}^{n+1}$ with that of $U_0^n$, $U_{P+1}^n$. The expression of $h_i^{\text{ext}}$ can then be established by the mean of Lemma A.2. We have now to find an analytic expression for the quantity $\int \frac{x^{k+1}}{(1+x)^k} dx$ in order to obtain the final expression $(hu)_i^{\text{ext}}$. The following lemma holds.

**Lemma B.1** Let $k \in \mathbb{N}^*$. If we denote $y = 1 - \frac{1}{1+x}$ for all $x \in \mathbb{R} \setminus \{-1\}$ and $C \in \mathbb{R}$ we have the following expression

$$ \int \frac{x^{k+1}}{(1+x)^k} dx = \left( -\sum_{r=1}^{k-2} (k-r-1) \frac{k-r}{2} \frac{y^r}{r} \right) \mathbb{1}_{k \ge 3} - \left( \frac{k(k-1)}{2} \ln|1-y| \right) \mathbb{1}_{k \ge 2} \\ - \frac{k+1}{(1-y)} + \frac{1}{2(1-y)^2} - \left( \sum_{q=1}^{k-1} (k-q) \frac{y^q}{q} \right) \mathbb{1}_{k \ge 2} - k \ln|1-y| + C $$

**Proof.** We begin by performing the change of variable $y = 1 - \frac{1}{1+x}$

$$ \int \frac{x^{k+1}}{(1+x)^k} dx = \int y^k \left( \frac{1}{1-y} - 1 \right) \frac{dy}{(1-y)^2} = \int \frac{y^k}{(1-y)^3} dy - \int \frac{y^k}{(1-y)^2} dy. $$
[25] B. Perthame, Second-order Boltzmann schemes for compressible Euler equations in one and two space dimensions. *SIAM J. Numer. Anal.* **29** (1992) 1–19.

[26] B. Perthame, Kinetic Formulation of Conservation Laws. *Oxford Lecture Series in Mathematics and its Applications*. Vol. 21. Oxford University Press, Oxford (2002).

[27] B. Perthame and C. Simeoni, A kinetic scheme for the Saint-Venant system with a source term. *Calcolo* **38** (2001) 201–231.

[28] Saint-Venant, Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit. *C. R. Acad. Sci. Paris* **73** (1871) 147–154.

[29] A. Thomann, G. Puppo and C. Klingenberg, An all speed second order well-balanced imex relaxation scheme for the euler equations with gravity. *J. Comput. Phys.* **420** (2020) 109723. https://www.sciencedirect.com/science/article/pii/S0021999120304976.

[30] Y. Xing and C.W. Shu, A survey of high order schemes for the shallow water equations. *J. Math. Study* **47** (2014) 221–249. semanticscholar.org/paper/0b3a29cd823dceb85e3a0a1f9622a306bc62591c.

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

Please help to maintain this journal in open access!

This journal is currently published in open access under the *Subscribe to Open model* (S2O). We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.

## APPENDIX A. EXPRESSION OF THE NUMERICAL FLUXES

The optimal choice for the Maxwellian is given by (1.4). Unfortunately the explicit expression for the numerical fluxes appearing in (2.20) is hardly possible with the choice (1.4) and the use of approximate quadrature formula for the integrals in (1.4) will degrade the accuracy of the scheme and increase the computational costs. Hence, we choose *M* defined by the first expression in (2.17) and relation (2.20) becomes

$$U_i^{\text{int}} = \frac{1}{2\sqrt{3}} \left( \int_{\min(0, a_j^n)}^{\min(0, b_j^n)} \sum_{j=i}^{P} \sqrt{\frac{2h_j^n}{g}} \left(\xi\right) \frac{(-\sigma\xi)^{j-i}}{(1-\sigma\xi)^{j-i+1}} d\xi + \int_{\max(0, a_j^n)}^{\max(0, b_j^n)} \sum_{j=1}^{i} \sqrt{\frac{2h_j^n}{g}} \left(\xi\right) \frac{(\sigma\xi)^{i-j}}{(1+\sigma\xi)^{i-j+1}} d\xi \right),$$

with $a_j^n = u_j^n - \sqrt{3} c_j^n$ and $b_j^n = u_j^n + \sqrt{3} c_j^n$. The expressions of $h_i^{\text{int}}$ and $(hu)_i^{\text{int}}$ are given by

$$h_i^{\text{int}} = \frac{1}{2\sqrt{3}} \left( \sum_{j=i}^{P} \sqrt{\frac{2h_j^n}{g}} \underbrace{\int_{\min(0, a_j^n)}^{\min(0, b_j^n)} \frac{(-\sigma\xi)^{j-i}}{(1-\sigma\xi)^{j-i+1}} d\xi}_{(Ah)_{i,j}} + \sum_{j=1}^{i} \sqrt{\frac{2h_j^n}{g}} \underbrace{\int_{\max(0, a_j^n)}^{\max(0, b_j^n)} \frac{(\sigma\xi)^{i-j}}{(1+\sigma\xi)^{i-j+1}} d\xi}_{(Bh)_{i,j}} \right) \quad (A.1)$$

$$(hu)_i^{\text{int}} = \frac{1}{2\sqrt{3}} \left( \sum_{j=i}^{P} -\frac{1}{\sigma} \sqrt{\frac{2h_j^n}{g}} \underbrace{\int_{\min(0, a_j^n)}^{\min(0, b_j^n)} \frac{(-\sigma\xi)^{j-i+1}}{(1-\sigma\xi)^{j-i+1}} d\xi}_{(Ahu)_{i,j}} + \sum_{j=1}^{i} \frac{1}{\sigma} \sqrt{\frac{2h_j^n}{g}} \underbrace{\int_{\max(0, a_j^n)}^{\max(0, b_j^n)} \frac{(\sigma\xi)^{i-j+1}}{(1+\sigma\xi)^{i-j+1}} d\xi}_{(Bhu)_{i,j}} \right). \quad (A.2)$$

Now we need to compute analytically the integrals of both expressions using the following lemmas.
Making use of $y^k - 1 = (y-1)(y^{k-1} + y^{k-2} + \cdots + 1)$ as before, we remark the following relation for $k \ge 1$

$$ \frac{y^k}{1-y} = \frac{y^k - 1}{1-y} + \frac{1}{1-y} = - \sum_{p=0}^{k-1} y^p + \frac{1}{1-y} $$

Dividing this by $1-y$ leads to

$$ \begin{aligned} \frac{y^k}{(1-y)^2} &= - \sum_{p=0}^{k-1} \frac{y^p}{1-y} + \frac{1}{(1-y)^2} = - \sum_{p=0}^{k-1} \left( \frac{y^p - 1}{1-y} + \frac{1}{1-y} \right) + \frac{1}{(1-y)^2} \\ &= \left( \sum_{p=1}^{k-1} \sum_{q=0}^{p-1} y^q \right) \mathbb{1}_{k \ge 2} - \frac{k}{1-y} + \frac{1}{(1-y)^2} \end{aligned} $$

Iterating this one more time we find

$$ \begin{aligned} \frac{y^k}{(1-y)^3} &= \left( \sum_{p=1}^{k-1} \sum_{q=0}^{p-1} \frac{y^q}{1-y} \right) \mathbb{1}_{k \ge 2} - \frac{k}{(1-y)^2} + \frac{1}{(1-y)^3} \\ &= \left( \sum_{p=1}^{k-1} \sum_{q=0}^{p-1} \frac{y^q - 1}{1-y} + \frac{1}{1-y} \right) \mathbb{1}_{k \ge 2} - \frac{k}{(1-y)^2} + \frac{1}{(1-y)^3} \\ &= \left( - \sum_{p=2}^{k-1} \sum_{q=1}^{p-1} \sum_{r=0}^{q-1} y^r \right) \mathbb{1}_{k \ge 3} + \frac{k(k-1)}{2(1-y)} \mathbb{1}_{k \ge 2} - \frac{k}{(1-y)^2} + \frac{1}{(1-y)^3} \end{aligned} $$

As a consequence we get the following primitives up to a constant

$$ \int \frac{y^k}{(1-y)^2} dy = \left( \sum_{p=1}^{k-1} \sum_{q=1}^{p} \frac{y^q}{q} \right) \mathbb{1}_{k \ge 2} + k \ln|1-y| + \frac{1}{(1-y)} $$

$$ \int \frac{y^k}{(1-y)^3} dy = \left( - \sum_{p=2}^{k-1} \sum_{q=1}^{p-1} \sum_{r=1}^{q} \frac{y^r}{r} \right) \mathbb{1}_{k \ge 3} - \frac{k(k-1)}{2} \ln|1-y| \mathbb{1}_{k \ge 2} \\ - \frac{k}{(1-y)} + \frac{1}{2(1-y)^2} $$

Finally, we simplify the double and triple sums

$$ \sum_{p=1}^{k-1} \sum_{q=1}^{p} \frac{y^q}{q} = \sum_{q=1}^{k-1} \sum_{p=q}^{k-1} \frac{y^q}{q} = \sum_{q=1}^{k-1} (k-q) \frac{y^q}{q} $$

From this we deduce that

$$ \begin{aligned} \sum_{p=2}^{k-1} \sum_{q=1}^{p-1} \sum_{r=1}^{q} \frac{y^r}{r} &= \sum_{p=2}^{k-1} \sum_{r=1}^{p-1} (p-r) \frac{y^r}{r} \\ &= \sum_{p=1}^{k-2} \sum_{r=1}^{p} (p-r+1) \frac{y^r}{r} = \sum_{r=1}^{k-2} \sum_{p=r}^{k-2} (p-r+1) \frac{y^r}{r} \\ &= \sum_{r=1}^{k-2} \left( \frac{(k-r-1)(k+r-2)}{2} + (k-r-1)(1-r) \right) \frac{y^r}{r} \\ &= \sum_{r=1}^{k-2} (k-r-1) \frac{k-r}{2} \frac{y^r}{r} \end{aligned} $$
**Lemma A.1.** If we denote $y = 1 - \frac{1}{1+x}$ for all $x \in \mathbb{R}\setminus\{-1\}$ and $C \in \mathbb{R}$ we have the following primitive:

$$ \int \frac{x^k}{(1+x)^{k+1}} dx = \ln(|1+x|) - \sum_{l=1}^{k} \frac{y^l}{l} + C. $$

**Lemma A.2.** Using the same notation as in the previous lemma, we have

$$ \int \frac{x^k}{(1+x)^k} dx = -k \ln(|1+x|) + x + \sum_{l=1}^{k-1} l \frac{y^{k-l}}{k-l} + C'. $$

*Proof of Lemma A.1.* We have

$$ I = \int \frac{x^k}{(1+x)^{k+1}} dx = \int \frac{x^k}{(1+x)^k} \frac{1}{1+x} dx = \int \left(1 - \frac{x}{1+x}\right)^k \frac{1}{1+x} dx. $$

We pose $y = 1 - \frac{1}{1+x}$

$$ I = \int y^k (1-y) \frac{dy}{(1-y)^2} = \int \frac{y^k - 1}{1-y} + \frac{1}{1-y} dy. $$

Now we use the formula $y^k - 1 = (y-1)(y^{k-1} + y^{k-2} + \dots + y + 1)$. And we obtain

$$ \begin{aligned} I &= - \int \sum_{l=0}^{k-1} y^l dy - \ln(|1-y|) + C && C \in \mathbb{R} \\ &= \ln(|1+x|) - \sum_{l=1}^{k} \frac{y^l}{l} + C' && C' \in \mathbb{R}. \end{aligned} $$

*Proof of Lemma A.2.* We already have denoted $y = \frac{x}{1+x} = 1 - \frac{1}{1+x}$

$$ I = \int \left( \frac{x}{1+x} \right)^k dx = \int \frac{y^k dy}{(1-y)^2} = \int \left( \frac{y^k - 1}{(1-y)^2} + \frac{1}{(1-y)^2} \right) dy $$

where the formula $y^k - 1 = (y-1)(y^{k-1} + y^{k-2} + \dots + y + 1)$ has been used. Hence

$$ \begin{aligned} I &= - \int \sum_{l=0}^{k-1} \frac{y^l}{1-y} dy + x + C = - \int \sum_{l=0}^{k-1} \frac{y^l - 1}{1-y} dy - \int \frac{1}{1-y} \sum_{l=0}^{k-1} dy + x + C \\ &= \int \sum_{l=1}^{k-1} \frac{y^l - 1}{y-1} dy + k \ln(|1-y|) + x + C' = \int \sum_{l=1}^{k-1} \sum_{p=0}^{l-1} y^p dy - k \ln(|1+x|) + x + C' \\ &= \sum_{l=1}^{k-1} l \int y^{k-1-l} dy - k \ln(|1+x|) + x + C' = \sum_{l=1}^{k-1} l \frac{y^{k-l}}{k-l} - k \ln(|1+x|) + x + C'', \end{aligned} $$

with $(C, C', C'') \in \mathbb{R}^3$.

We are now able to compute the quantities $Ah_{i,j}, Bh_{i,j}, Ah_{u_{i,j}}, Bhu_{i,j}$

$$ Ah_{i,j} = \int_{\min(0, a_j^n)}^{\min(0, b_j^n)} \frac{(-\sigma \xi)^{j-i}}{(1 - \sigma \xi)^{j-i+1}} d\xi = -\frac{1}{\sigma} \int_{-\min(0, a_j^n)}^{-\min(0, b_j^n)} \frac{(x)^{j-i}}{(1+x)^{j-i+1}} dx $$
As a conclusion we have the expression

$$ \int \frac{x^{k+1}}{(1+x)^k} dx = \left( -\sum_{r=1}^{k-2} (k-r-1) \frac{k-r}{2} \frac{y^r}{r} \right) \mathbb{1}_{k \ge 3} - \left( \frac{k(k-1)}{2} \ln|1-y| \right) \mathbb{1}_{k \ge 2} \\ - \frac{k+1}{(1-y)} + \frac{1}{2(1-y)^2} - \left( \sum_{q=1}^{k-1} (k-q) \frac{y^q}{q} \right) \mathbb{1}_{k \ge 2} - k \ln|1-y| + C $$

where $C \in \mathbb{R}$ and with $y = x/(x+1)$. $\square$
$$ = \frac{1}{\sigma} \left[ \ln(|1+x|) - \sum_{l=1}^{j-i} \frac{y^l}{l} \right]_{-\min(0, b_j^n)\sigma}^{-\min(0, a_j^n)\sigma} . \qquad (A.3) $$

$$ \begin{aligned} B h_{i,j} &= \int_{\max(0, a_j^n)}^{\max(0, b_j^n)} \frac{(\sigma\xi)^{i-j}}{(1 + \sigma\xi)^{i-j+1}} d\xi = \frac{1}{\sigma} \int_{\max(0, a_j^n)\sigma}^{\max(0, b_j^n)\sigma} \frac{(x)^{i-j}}{(1+x)^{i-j+1}} dx \\ &= \frac{1}{\sigma} \left[ \ln(|1+x|) - \sum_{l=1}^{i-j} \frac{y^l}{l} \right]_{\max(0, a_j^n)\sigma}^{\max(0, b_j^n)\sigma} . \end{aligned} \qquad (A.4) $$

And similarly we obtain the formulas for Ahu and Bhu under the form

$$ Ahu_{i,j} = \frac{1}{\sigma} \left[ -(j-i+1) \ln(|1+x|) + x + \sum_{l=1}^{j-i} l \frac{y^{j-i+1-l}}{j-i+1-l} \right]_{-\min(0, b_j^n)\sigma}^{-\min(0, a_j^n)\sigma} \qquad (A.5) $$

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

To conclude this paragraph, we give the final expression of $U_i^{\text{int}}$

$$ h_i^{\text{int}} = \frac{1}{2\sigma\sqrt{3}} \left( \sum_{j=i}^{P} \sqrt{\frac{2h_j^n}{g}} \left[ \ln(|1+x|) - \sum_{l=1}^{j-i} \frac{y^l}{l} \right]_{-\min(0, b_j^n)\sigma}^{-\min(0, a_j^n)\sigma} + \sum_{j=1}^{i} \sqrt{\frac{2h_j^n}{g}} \left[ \ln(|1+x|) - \sum_{l=1}^{i-j} \frac{y^l}{l} \right]_{\max(0, a_j^n)\sigma}^{\max(0, b_j^n)\sigma} \right) \qquad (A.7) $$

$$ \begin{aligned} (hu)_i^{\text{int}} = \frac{1}{2\sigma^2\sqrt{3}} & \left( -\sum_{j=i}^{P} \sqrt{\frac{2h_j^n}{g}} \left[ -(j-i+1)\ln(|1+x|) + x + \sum_{k=1}^{j-i} (j-i+1-k)\frac{y^k}{k} \right]_{-\min(0, b_j^n)\sigma}^{-\min(0, a_j^n)\sigma} \right. \\ & \left. + \sum_{j=1}^{i} \sqrt{\frac{2h_j^n}{g}} \left[ -(i-j+1)\ln(|1+x|) + x + \sum_{k=1}^{i-j} (i-j+1-k)\frac{y^k}{k} \right]_{\max(0, a_j^n)\sigma}^{\max(0, b_j^n)\sigma} \right). \end{aligned} \qquad (A.8) $$

## APPENDIX B. COMPUTATIONS OF THE FLUXES INVOLVING THE BOUNDARY CONDITIONS

We assume the ghost quantities $U_0^{n+1}$ and $U_{P+1}^{n+1}$ at time $t^{n+1}$ to be known. The exterior contribution given in (2.19) also writes

$$ U_i^{\text{ext}} = \int_{\mathbb{R}^{-}} \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_{\mathbb{R}^{+}} \left(\frac{1}{\xi}\right) \frac{(\sigma\xi)^i}{(1+\sigma\xi)^i} M_0^{n+1} d\xi. $$

Using computations similar to what has been proposed in Appendix A, we get

$$ U_i^{\text{ext}} = \frac{1}{2\sqrt{3}} \left[ \sqrt{\frac{2h_{P+1}^{n+1}}{g}} \int_{\min(0, a_{P+1}^{n+1})}^{\min(0, b_{P+1}^{n+1})} \left(\frac{1}{\xi}\right) \frac{(-\sigma\xi)^{P-i+1}}{(1-\sigma\xi)^{P-i+1}} d\xi + \sqrt{\frac{2h_0^{n+1}}{g}} \int_{\max(0, a_0^{n+1})}^{\max(0, b_0^{n+1})} \left(\frac{1}{\xi}\right) \frac{(\sigma\xi)^i}{(1+\sigma\xi)^i} d\xi \right], $$

or equivalently

$$ h_i^{\text{ext}} = \frac{1}{2\sqrt{3}} \left[ \sqrt{\frac{2h_{P+1}^{n+1}}{g}} \int_{\min(0, a_{P+1}^{n+1})}^{\min(0, b_{P+1}^{n+1})} \frac{(-\sigma\xi)^{P-i+1}}{(1-\sigma\xi)^{P-i+1}} d\xi + \sqrt{\frac{2h_0^{n+1}}{g}} \int_{\max(0, a_0^{n+1})}^{\max(0, b_0^{n+1})} \frac{(\sigma\xi)^i}{(1+\sigma\xi)^i} d\xi \right], $$

$$ (hu)_i^{\text{ext}} = \frac{1}{2\sqrt{3}} \left[ \sqrt{\frac{2h_{P+1}^{n+1}}{g}} \int_{\min(0, a_{P+1}^{n+1})}^{\min(0, b_{P+1}^{n+1})} \xi \frac{(-\sigma\xi)^{P-i+1}}{(1-\sigma\xi)^{P-i+1}} d\xi + \sqrt{\frac{2h_0^{n+1}}{g}} \int_{\max(0, a_0^{n+1})}^{\max(0, b_0^{n+1})} \xi \frac{(\sigma\xi)^i}{(1+\sigma\xi)^i} d\xi \right]. $$

As explained in Section 2.4, in practice we will replace the unknown values of $U_0^{n+1}$, $U_{P+1}^{n+1}$ with that of $U_0^n$, $U_{P+1}^n$. The expression of $h_i^{\text{ext}}$ can then be established by the mean of Lemma A.2. We have now to find an analytic expression for the quantity $\int \frac{x^{k+1}}{(1+x)^k} dx$ in order to obtain the final expression $(hu)_i^{\text{ext}}$. The following lemma holds.
**Lemma B.1.** Let $k \in \mathbb{N}^*$. If we denote $y = 1 - \frac{1}{1+x}$ for all $x \in \mathbb{R} \setminus \{-1\}$ and $C \in \mathbb{R}$ we have the following expression

$$ \int \frac{x^{k+1}}{(1+x)^k} dx = \left( -\sum_{r=1}^{k-2} (k-r-1) \frac{k-r}{2} \frac{y^r}{r} \right) \mathbb{1}_{k \ge 3} - \left( \frac{k(k-1)}{2} \ln|1-y| \right) \mathbb{1}_{k \ge 2} \\ - \frac{k+1}{(1-y)} + \frac{1}{2(1-y)^2} - \left( \sum_{q=1}^{k-1} (k-q) \frac{y^q}{q} \right) \mathbb{1}_{k \ge 2} - k \ln|1-y| + C. $$

*Proof.* We begin by performing the change of variable $y = 1 - \frac{1}{1+x}$

$$ \int \frac{x^{k+1}}{(1+x)^k} dx = \int y^k \left( \frac{1}{1-y} - 1 \right) \frac{dy}{(1-y)^2} = \int \frac{y^k}{(1-y)^3} dy - \int \frac{y^k}{(1-y)^2} dy. $$

Making use of $y^k - 1 = (y-1)(y^{k-1} + y^{k-2} + \dots + 1)$ as before, we remark the following relation for $k \ge 1$

$$ \frac{y^k}{1-y} = \frac{y^k - 1}{1-y} + \frac{1}{1-y} = - \sum_{p=0}^{k-1} y^p + \frac{1}{1-y}. $$

Dividing this by $1-y$ leads to

$$ \frac{y^k}{(1-y)^2} = - \sum_{p=0}^{k-1} \frac{y^p}{1-y} + \frac{1}{(1-y)^2} = - \sum_{p=0}^{k-1} \left( \frac{y^p - 1}{1-y} + \frac{1}{1-y} \right) + \frac{1}{(1-y)^2} \\ = \left( \sum_{p=1}^{k-1} \sum_{q=0}^{p-1} y^q \right) \mathbb{1}_{k \ge 2} - \frac{k}{1-y} + \frac{1}{(1-y)^2}. $$

Iterating this one more time we find

$$ \begin{aligned} \frac{y^k}{(1-y)^3} &= \left( \sum_{p=1}^{k-1} \sum_{q=0}^{p-1} \frac{y^q}{1-y} \right) \mathbb{1}_{k \ge 2} - \frac{k}{(1-y)^2} + \frac{1}{(1-y)^3} \\ &= \left( \sum_{p=1}^{k-1} \sum_{q=0}^{p-1} \frac{y^q - 1}{1-y} + \frac{1}{1-y} \right) \mathbb{1}_{k \ge 2} - \frac{k}{(1-y)^2} + \frac{1}{(1-y)^3} \\ &= \left( - \sum_{p=2}^{k-1} \sum_{q=1}^{p-1} \sum_{r=0}^{q-1} y^r \right) \mathbb{1}_{k \ge 3} + \frac{k(k-1)}{2(1-y)} \mathbb{1}_{k \ge 2} - \frac{k}{(1-y)^2} + \frac{1}{(1-y)^3}. \end{aligned} $$

As a consequence we get the following primitives up to a constant

$$ \int \frac{y^k}{(1-y)^2} dy = \left( \sum_{p=1}^{k-1} \sum_{q=1}^{p} \frac{y^q}{q} \right) \mathbb{1}_{k \ge 2} + k \ln|1-y| + \frac{1}{(1-y)} \\ \int \frac{y^k}{(1-y)^3} dy = \left( - \sum_{p=2}^{k-1} \sum_{q=1}^{p-1} \sum_{r=1}^{q} \frac{y^r}{r} \right) \mathbb{1}_{k \ge 3} - \frac{k(k-1)}{2} \ln|1-y| \mathbb{1}_{k \ge 2} - \frac{k}{(1-y)} + \frac{1}{2(1-y)^2}. $$

Finally, we simplify the double and triple sums

$$ \sum_{p=1}^{k-1} \sum_{q=1}^{p} \frac{y^q}{q} = \sum_{q=1}^{k-1} \sum_{p=q}^{k-1} \frac{y^q}{q} = \sum_{q=1}^{k-1} (k-q) \frac{y^q}{q}. $$

From this we deduce that

$$ \sum_{p=2}^{k-1} \sum_{q=1}^{p-1} \sum_{r=1}^{q} \frac{y^r}{r} = \sum_{p=2}^{k-1} \sum_{r=1}^{p-1} (p-r) \frac{y^r}{r} $$
$$ 
\begin{aligned}
&= \sum_{p=1}^{k-2} \sum_{r=1}^{p} (p-r+1) \frac{y^r}{r} = \sum_{r=1}^{k-2} \sum_{p=r}^{k-2} (p-r+1) \frac{y^r}{r} \\
&= \sum_{r=1}^{k-2} \left( \frac{(k-r-1)(k+r-2)}{2} + (k-r-1)(1-r) \right) \frac{y^r}{r} \\
&= \sum_{r=1}^{k-2} (k-r-1) \frac{k-r}{2} \frac{y^r}{r}.
\end{aligned}
$$

As a conclusion we have the expression

$$ 
\begin{aligned}
\int \frac{x^{k+1}}{(1+x)^k} dx = & \left( -\sum_{r=1}^{k-2} (k-r-1) \frac{k-r}{2} \frac{y^r}{r} \right) \mathbb{1}_{k \ge 3} - \left( \frac{k(k-1)}{2} \ln|1-y| \right) \mathbb{1}_{k \ge 2} \\
& - \frac{k+1}{(1-y)} + \frac{1}{2(1-y)^2} - \left( \sum_{q=1}^{k-1} (k-q) \frac{y^q}{q} \right) \mathbb{1}_{k \ge 2} - k \ln|1-y| + C
\end{aligned}
$$

where $C \in \mathbb{R}$ and with $y = x/(x+1)$. $\square$