# A semi-implicit exactly fully well-balanced relaxation scheme for the Shallow Water Linearized Moment Equations

Celia Caballero, Irene Gómez-Bueno, Alessia Del Grosso, Julian Koellermeier, Tomás Morales de Luna

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

To cite this version:

Celia Caballero, Irene Gómez-Bueno, Alessia Del Grosso, Julian Koellermeier, Tomás Morales de Luna. A semi-implicit exactly fully well-balanced relaxation scheme for the Shallow Water Linearized Moment Equations. 2025. hal-04885181

HAL Id: hal-04885181

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

Preprint submitted on 14 Jan 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A semi-implicit exactly fully well-balanced
relaxation scheme for the Shallow Water
Linearized Moment Equations

C. Caballero-Cárdenas, I. Gómez-Bueno, A. Del Grosso,
J. Koellermeier, T. Morales de Luna

## Abstract

When dealing with shallow water simulations, the velocity profile is often assumed to be constant along the vertical axis. However, since in many applications this is not the case, modeling errors can be significant. Hence, in this work, we deal with the Shallow Water Linearized Moment Equations (SWLME), in which the velocity is no longer constant in the vertical direction, where a polynomial expansion around the mean value is considered. The linearized version implies neglecting the non-linear terms of the basis coefficients in the higher order equations. As a result, the model is always hyperbolic and the stationary solutions can be more easily computed. Then, our objective is to propose an efficient, accurate and robust numerical scheme for the SWLME model, specially adapted for low Froude number situations. Hence, we describe a semi-implicit second order exactly fully well-balanced method. More specifically, a splitting is performed to separate acoustic and material phenomena. The acoustic waves are treated in an implicit manner to gain in efficiency when dealing with subsonic flow regimes, whereas the second order of accuracy is achieved thanks to a polynomial reconstruction and Strang-splitting method. We also exploit a reconstruction operator to achieve the fully well-balanced character of the method. Extensive numerical tests demonstrate the well-balanced properties and large speed-up compared to traditional methods.

## 1 Introduction

Shallow water models have many applications in atmospheric and oceanic dynamics, which have motivated extensive research over the past decades: [3, 6, 7, 14, 19, 22, 26, 27, 36, 50, 51, 56, 57], among many others. These models are commonly represented by the shallow water equations, or the Saint-Venant system in one-dimensional form. The equations can be derived by performing a depth-averaging procedure on the Navier-Stokes equations, which assumes a constant horizontal velocity across the vertical axis. While this assumption can simplify calculations in many cases, there are numerous situations where a more
complex approach is required. The vertical averaging involved leads to a loss of
vertical detail, introducing modelling errors that may be significant and should
not be neglected [31, 34, 44].

Looking at the literature, a first possibility in this regard is to decompose
the fluid into two or more layers, each with its own constant horizontal velocity
[1, 4, 11, 19, 20, 29, 31, 32, 39].

In another direction, a second possibility is provided by shallow water mo-
ment models, in which the velocity is no longer constant in the vertical direction
as a polynomial expansion around the mean value is considered [45]. The expan-
sion coefficients, also called moments, allow for the representation of complex
vertical velocity profiles. The greater the number of moments, the better the ver-
tical approximation of the velocity. However, the model may not be hyperbolic
if two (or more) moments are considered. We note that a number of alternative
models has been derived in the literature for the 1D case [17, 35, 41, 44, 53]
and more recently also for the 2D extension [5, 55]. The technique can also
be combined with a multilayer approach by proposing a vertical decomposition
into layers that use withing a moment approach [34].

Let us remark that the derivation of steady states for those models is intri-
cate due to the combination of non-conservative products and non-linear terms.
A solution to this problem is provided in [43], where the authors propose a par-
tially linearized version of the model which allows for small deviations from the
constant velocity profile. The new Shallow Water Linearized Moment Equations
(SWLME) are obtained neglecting non-linear terms of the basis coefficients in
the higher order equations. Not only this new model is always hyperbolic, re-
gardless of the number of moments, but it also allows for a simpler description
of the model’s stationary solutions, whereas before they were complicated both
analytically and numerically. This fact is of no small importance as it is es-
sential to be able to preserve steady states to deal with small deviations from
them, otherwise numerical artifacts and spurious oscillations may appear in the
numerical results.

Numerical preservation of steady states has recently gained a lot of interest. A numerical scheme with this property is defined either *well-balanced* [2, 3, 21, 30, 50] if only zero-velocity steady states are preserved or *fully well-balanced* [7, 8, 9, 18, 21] if all steady states can be preserved. In fact, one may distinguish between well-balanced or exactly well-balanced schemes depending on whether an approximation of the exact steady states are preserved [37].

For all the previous reasons, in this paper, we aim to numerically simulate
shallow water flows with non-constant vertical velocity by means of the SWLME
model. Then, our objective is to propose an efficient, accurate and robust
numerical scheme for the selected model. Concerning the last property, we aim
for a scheme which is not only consistent, but also exactly fully well-balanced.
Then, the numerical method should be efficient, namely numerical simulations
should be fast by neglecting restrictive CFL conditions for the time step. A large
part of applications and simulations for shallow water flows occur in subsonic
regimes, i.e. when material waves are much slower than acoustic waves. In
this type of situation, the CFL condition may be restrictive due to the fast
acoustic waves, which are not of main interest as they quickly leave the domain. Hence, we aim to deal with this problem by exploiting a numerical scheme whose CFL condition is based on the slow material waves. This also implies that the transport waves will be better approximated as the CFL condition is adapted to them. For this purpose, we look for an approach that treats implicitly the acoustic waves and explicitly the material waves. Different possibilities arise for this purpose. One option would be to use an IMEX approach as in [10, 23, 46]. Another alternative is the use of Lagrange-projection methods, which provide a natural decomposition of acoustic and transport waves [13, 15, 24, 25, 28]. However, even if this strategy was extended to two-layer shallow water flows [29], this was far from trivial and the moment models seem to be even more difficult to deal with in that framework. A different method for splitting of acoustic and transport phenomena can be found in [42]. It was originally developed for the compressible Euler equations and then extended to the Ripa model [33] and the shallow water system [16]. It combines a relaxation approach with an splitting of acoustic and transport waves that avoids several difficulties when compared to Lagrange-projection approach.

It is important to take into account that implicit schemes can be in general highly diffusive, hence at least second order of accuracy is needed. However, combining implicit methods with a higher order of accuracy could become a complex task, since a naive extension could lead to non-linear systems to be solved iteratively and thus a higher computational cost. Of course, this would be in contradiction with the goal of having an efficient method. One of the advantages of this technique, as shown in [16], is precisely the fact that, thanks to the relaxation, only linear systems need to be solved.

In conclusion, we aim to develop a numerical scheme, which is second order accurate and exactly fully well-balanced for the Shallow Water Linearized Moment Equations (SWLME), that splits the acoustic and transport phenomena. The former are then solved implicitly in a efficient way using a relaxation approach.

*Outline of the paper.* In the next section 2, we present the mathematical model and splitting we exploit for the numerical simulations and schemes. Then, both first- and second order methods are described in section 3 along with details for the exactly fully well-balanced character. Numerical simulations are presented in section 4 to show the good behaviour of the schemes. Finally, we draw conclusions in section 5.

## 2 Model description

The classical shallow water equations (SWE) model the evolution of the water height $h$ and the mean velocity $u_0$ under shallow conditions. They can be derived by depth-averaging of the incompressible Navier-Stokes equations. For the frictionless case with gravitational constant $g$ and bottom topography $z(x)$,
the evolution equations are given by

$$ \partial_t \left( \frac{h}{hu_0} \right) + \partial_x \left( hu_0 + \frac{1}{2} gh^2 \right) = \begin{pmatrix} 0 \\ -gh\partial_x z \end{pmatrix}, \quad (1) $$

The simplicity of the model allows to easily compute analytical properties of the model such as its eigenvalues, given by $\lambda_{1,2} = u_0 \pm \sqrt{gh}$, and its steady states given by

$$ hu_0 = \text{const}, \quad (2) $$

$$ \frac{1}{2} u_0^2 + g(h+z) = \text{const.} \quad (3) $$

However, the SWE suffer from the inability to model vertical variations of the velocity, since they only include a depth-averaged velocity $u_0$. It has been shown in many applications, that vertical velocity variations are an important feature in shallow flows and need to be taken into account for an accurate representation, see also [45]. To overcome this deficiency, different extended shallow water models have been derived. The model used in this work is based on the Shallow Water Moment Equations (SWME), which assume an expansion of the velocity profile $u(t, x, \xi)$ around the depth-averaged value:

$$ u(t, x, \xi) = u_0(t, x) + \sum_{i=1}^{N} u_i(t, x) \phi_i \left( \frac{\xi - z(x)}{h(t, x)} \right), \quad \xi \in [z(x), z(x) + h(t, x)] \quad (4) $$

where $N$ is referred to as the order of the model, i.e., how many expansion terms are used, $u_i$ are the corresponding expansion coefficients, also called *moments*, and $\phi_i: [0, 1] \to \mathbb{R}$ are *scaled Legendre polynomials* of degree $i$, also called basis functions, defined as

$$ \phi_i(\zeta) = \frac{1}{i!} \frac{d^i}{d\zeta^i} (\zeta - \zeta^2)^i. \quad (5) $$

The expansion (4) allows to include polynomial velocity profiles at the expense of including the additional moments $u_i$ as variables. Evolution equations for these moments are derived by computing higher order projections of the Navier Stokes equations [45]. This results in the SWME, a closed system of $N+2$ evolution equations for the variables $h, u_0, u_1, \dots, u_N$. Nevertheless, SWME present two main drawbacks: (1) the model is not always hyperbolic, leading to potential instabilities as analysed in detail in [44], and (2) it is not possible to derive an easy and general expression, using algebraic relations, for steady states as in (2), due to the multitude of non-linear terms originating from the expansion and projection.

In this paper, we consider the so-called Shallow Water Linearized Moment Equations (SWLME), derived first in [43] and further analyzed in [41]. The model overcomes both deficiencies of the original SWME model from [45] using the following assumption: the velocity profile $u(t, x, \xi)$ only has small deviations from a constant profile, i.e., $u_i = O(\epsilon)$ for $i \ge 1$. In the higher-order equations,
all term of order $O(\epsilon^2)$ are then neglected, while the conservation of mass and the momentum equation are left unchanged.

The SWLME are given by

$$ \partial_t \begin{pmatrix} h \\ hu_0 \\ hu_1 \\ \vdots \\ hu_N \end{pmatrix} + \partial_x \begin{pmatrix} hu_0 \\ hu_0 + g \frac{h^2}{2} + \frac{1}{3}hu_1^2 + \cdots + \frac{1}{2N+1}hu_N^2 \\ 2hu_0u_1 \\ \vdots \\ 2hu_0u_N \end{pmatrix} = Q\partial_x \begin{pmatrix} h \\ hu_0 \\ hu_1 \\ \vdots \\ hu_N \end{pmatrix} + \begin{pmatrix} 0 \\ -gh\partial_x z \\ 0 \\ \vdots \\ 0 \end{pmatrix}, \quad (6) $$

where the non-conservative right-hand side term includes $Q = \text{diag}(0, 0, u_0, \dots, u_0)$.

As proven in [43], the eigenvalues are given by

$$ \lambda_{1,2} = u_0 \pm \sqrt{gh + \sum_{i=1}^{N} \frac{3u_i^2}{2i+1}} \quad \text{and} \quad \lambda_{i+2} = u_0, \quad \text{for } i=1, \dots, N, \quad (7) $$

meaning that the system is hyperbolic for positive water heights.

The Froude number, defined by

$$ Fr(U) = \frac{|u_0|}{\sqrt{gh + \sum_{i=1}^{N} \frac{3u_i^2}{2i+1}}}, \quad (8) $$

determines the fluid regime: if its value is less than 1, it is subcritical; if it is 1 it represents a critical regime, whereas a Froude number bigger than 1 characterises a supercritical flow. In this work, we are especially interested in the subcritical regime, as this is where the improvements in computational efficiency of implicit versus explicit schemes are observed. However, our schemes could be used for other regimes too.

**Remark 1** The standard definition of the Froude number in equation (8) does not take into account the higher moments $u_i$, for $i > 0$. In [43], this has been first discussed by inclusion of additional dimensionless numbers that indicate the deviation from the constant profile. However, since we assume that the expansion coefficients are small, we shall consider here the standard definition of the Froude number from equation (8). Let us remark that one could define an extended Froude number by

$$ \widetilde{Fr}(U) = \frac{|u_0|}{\sqrt{gh + \sum_{i=1}^{N} \frac{3u_i^2}{2i+1}}} $$

Steady states can be computed by solving

$$ hu_0 = C_1, \quad (9) $$

$$ \frac{1}{2}u_0^2 + g(h+z) + \frac{3}{2} \sum_{i=1}^{N} \frac{1}{2i+1} u_i^2 = C_2, \quad (10) $$

$$ \frac{u_i}{h} = C_k, \quad i \ge 1, k \ge 3, \quad (11) $$
and we refer to [43] for analysis and discussion on the solvability for existing steady states.

The hyperbolicity and existence of easily computable steady states allows the model to be used in standard numerical simulations. As already said, here we aim at defining a numerical scheme that is well-balanced for the steady states and efficient, specially in the case of low Froude numbers. To this end, we shall consider a relaxation approach for system (6), given by

$$ 
\begin{aligned}
\partial_t \begin{pmatrix} h \\ hu_0 \\ hu_1 \\ \vdots \\ hu_N \\ h\pi \end{pmatrix} + \partial_x \begin{pmatrix} hu_0 \\ hu_0^2 + \pi + \frac{1}{3}hu_1^2 + \cdots + \frac{1}{2N+1}hu_N^2 \\ 2hu_0u_1 \\ \vdots \\ 2hu_0u_N \\ a^2u_0 \end{pmatrix} = \\
= \tilde{Q}\partial_x \begin{pmatrix} h \\ hu_0 \\ hu_1 \\ \vdots \\ hu_N \\ h\pi \end{pmatrix} + \begin{pmatrix} 0 \\ -gh\partial_x z \\ 0 \\ \vdots \\ 0 \\ 0 \end{pmatrix} + \begin{pmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ 0 \\ \frac{1}{\varepsilon} \left(g\frac{h^2}{2} - \pi\right) \end{pmatrix},
\end{aligned} 
\quad (12) $$

where $\tilde{Q} = \text{diag}(0, 0, u_0, \dots, u_0, 0)$, the constant $a$ has to fulfill the subcharacteristic condition (see [12]) given by:

$$ h\sqrt{gh} \le a. \quad (13) $$

Moreover, when $\varepsilon \to 0$, one has $\pi = g\frac{h^2}{2}$. Here, an instantaneous relaxation of the pressure will be considered as usual. Consequently, the $\varepsilon$-term will be neglected and $\pi$ is then set to an initial pressure value $g\frac{h^2}{2}$ at each time step. Therefore, the resulting model can be written in a compressed way as follows

$$ \partial_t U + \partial_x (F(U) + F_P(U)) + B(U)\partial_x U + S(U)\partial_x z = 0, \quad (14) $$

where

$$ F(U) = (hu_0, hu_0^2 + \frac{1}{3}hu_1^2 + \cdots + \frac{1}{2N+1}hu_N^2, 2hu_0u_1, \dots, 2hu_0u_N, 0)^T, \quad (15) $$

$$ F_P(U) = (0, \pi, 0, \dots, 0, a^2 u_0)^T, \quad (16) $$

$$ B = -\tilde{Q}, \quad (17) $$

$$ S(U) = (0, gh, 0, \dots, 0)^T. \quad (18) $$

Now, following the ideas in [16, 33], we perform a splitting, obtaining a pressure system and a transport one. The pressure system to be solved is written as

$$ \partial_t U + \partial_x F_P(U) + S(U)\partial_x z = 0, \quad (19) $$
and the one for the transport is given by

$$ \partial_t U + \partial_x F(U) + B(U) \partial_x U = 0. \qquad (20) $$

Note that the eigenvalues of system (19) are

$$ \lambda_{1,2}^{P} = \pm \frac{a}{h} \quad \text{and} \quad \lambda_{i+2}^{P} = 0, \quad \text{for } i = 1, \dots, N, \qquad (21) $$

while the ones of system (20) are

$$ \lambda_{1,2}^{T} = u_0 \pm \sqrt{\sum_{i=1}^{N} \frac{3u_i^2}{2i+1}} \quad \text{and} \quad \lambda_{i+2}^{T} = u_0, \quad \text{for } i = 1, \dots, N. \qquad (22) $$

This means that stability restriction for the pressure system comes from the pressure term, while in the transport system the stability is governed by velocity. Therefore, we shall consider an implicit approach for the pressure system, which will allow to overcome the restrictive CFL condition that comes from the pressure in the case of small Froude number.

## 3 Numerical scheme

We will consider a finite volume approach, that is, the spatial domain is divided into a series of cells, denoted by intervals $[x_{i-1/2}, x_{i+1/2})$ for each integer $i \in \mathbb{Z}$. For simplicity, we assume that each cell has a uniform length, represented by $\Delta x = x_{i+1/2} - x_{i-1/2}$. Here, the points $x_{i+1/2}$ mark the interfaces between cells, while the cell centers are designated as $x_i = \frac{x_{i-1/2} + x_{i+1/2}}{2}$. For the moment, the time variable will remain continuous. Moreover, the approximation of the cell averages will be denoted by

$$ U_i(t) \approx \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} U(x, t) dx. $$

Since we will consider here schemes up to second order, the cell averages will be identified with the cell center values, which corresponds to a mid point quadrature formula for computing the cell averages. The same quadrature rule will also be used to approximate the integrals of the source term and the non-conservative ones as it will explained later.

In order to obtain fully well-balanced schemes, we will follow the strategy described in [21], in which the well-balanced property is transferred to the reconstruction operator. A fully exactly well-balanced reconstruction operator $P_i(x)$ is a reconstruction operator that satisfies

$$ P_i(x, \{U_j^e\}_{j \in S_i}) = U^e(x), \quad \forall x \in [x_{i-1/2}, x_{i+1/2}], \quad \forall i, $$

for a given continuous stationary solution $U^e$, and where $S_i$ denotes the stencil corresponding to cell $i$.
Notice that a stationary solution $U^e$ must satisfy:

$$ \partial_x(F(U^e) + F_P(U^e)) + B(U^e)\partial_x U^e + S(U^e)\partial_x z = 0. \quad (23) $$

The stationary solutions can be computed by solving (9)-(11).

Specifically, being $t_0$ the initial time at each time step, we compute at every time step $[t_0, t_0 + \Delta t]$ and at every cell $i$, the steady state $U_i^{e,t_0}$ that satisfies (23) such that $U_i^{e,t_0}(x_i) = U_i^{t_0}$.

For each cell $i$, if the ODE (23) is subtracted from the system (14), and a similar splitting to the one described above is applied, we can write the semi-discrete version of the pressure and transport systems as follows:

$$ \begin{aligned} U'_i(t) = & -\frac{1}{\Delta x} (F_{P,i+1/2}(t) - F_P(U_i^{e,t_0}(x_{i+1/2}))) - F_{P,i-1/2}(t) + F_P(U_i^{e,t_0}(x_{i-1/2})) \\ & -\frac{1}{\Delta x} \left( \int_{x_{i-1/2}}^{x_{i+1/2}} (S(P_{i,U}(x)) - S(U_i^{e,t_0}(x))) z'(x) dx \right), \end{aligned} \quad (24) $$

and

$$ \begin{aligned} U'_i(t) = & -\frac{1}{\Delta x} (F_{i+1/2}(t) - F(U_i^{e,t_0}(x_{i+1/2}))) - F_{i-1/2}(t) + F(U_i^{e,t_0}(x_{i-1/2})) \\ & -\frac{1}{\Delta x} (\mathcal{B}_{i+1/2-}(t) + \mathcal{B}_{i-1/2+}(t)) \\ & -\frac{1}{\Delta x} \left( \int_{x_{i-1/2}}^{x_{i+1/2}} B(P_{i,U}(x, t_0)) \frac{dP_{i,U}}{dx} dx - \int_{x_{i-1/2}}^{x_{i+1/2}} B(U_i^{e,t_0}(x)) \frac{dU_i^{e,t_0}}{dx} dx \right), \end{aligned} \quad (25) $$

where

* $F_{P,i+1/2} = (0, \pi_{i+1/2}^*, 0, \dots, 0, a^2 u_{0,i+1/2}^*)^T$, being $\pi_{i\pm\frac{1}{2}}^*(t)$ and $u_{0,i\pm\frac{1}{2}}^*(t)$ approximations of the pressure and velocity at the interface, respectively, where $\pi_{i\pm\frac{1}{2}}^*(t) \approx \pi(x_{i\pm\frac{1}{2}}, t)$ and $u_{0,i\pm\frac{1}{2}}^*(t) \approx u_0(x_{i\pm\frac{1}{2}}, t)$.

* $P_{i,X}(x, t)$ is the fully well-balanced reconstruction operator of any variable $X$, that we will describe in subsection 3.1.

* $F_{i\pm 1/2} = \mathcal{F}(U_{i\pm 1/2-}, U_{i\pm 1/2+})$, where $\mathcal{F}$ is a consistent numerical flux, i.e. $\mathcal{F}(U, U) = F(U)$. In this work, the HLL numerical flux is considered (see [40]), that is,

$$ F_{i+1/2} = \frac{F(U_{i+1/2-}) + F(U_{i+1/2+})}{2} - \frac{1}{2} \left( \alpha_{i+1/2}^0 (U_{i+1/2+} - U_{i+1/2-}) + \alpha_{i+1/2}^1 (F(U_{i+1/2+}) - F(U_{i+1/2-})) \right) \quad (26) $$
where the coefficients $\alpha_{i+1/2}^0$ and $\alpha_{i+1/2}^1$ are defined as

$$ \alpha_{i+1/2}^0 = \frac{S_{i+1/2}^R |S_{i+1/2}^L| - S_{i+1/2}^L |S_{i+1/2}^R|}{S_{i+1/2}^R - S_{i+1/2}^L}, \quad \alpha_{i+1/2}^1 = \frac{|S_{i+1/2}^R| - |S_{i+1/2}^L|}{S_{i+1/2}^R - S_{i+1/2}^L}. \qquad (27) $$

Here, $S_{i+1/2}^L$ and $S_{i+1/2}^R$ are, respectively, approximations of the minimum and maximum wave propagation speeds, computed from equation (7).

• $B_{i+1/2\pm} = \frac{1}{2}B_{i+1/2}(U_{i+1/2+} - U_{i+1/2-}) \pm \frac{1}{2}\alpha_{i+1/2}^1 B_{i+1/2}(U_{i+1/2+} - U_{i+1/2-})$, where $B_{i+1/2}$ is a consistent matrix with the original matrix $B$. Note that the term $B(U)\partial_x U$ corresponds to a non-conservative product. Therefore, one must give a precise definition of weak solutions. One way to do so is by means of the path-conservative theory (see [52]). Nevertheless, we will not focus on this problem here, which is out of the scope of the paper. In practice, we define

$$ B_{i+1/2} = \int_{0}^{1} B(sU_{i+1/2+} + (1-s)U_{i+1/2-})ds \approx B((U_{i+1/2+} + U_{i+1/2-})/2) $$

In the previous definitions, the time dependency has been omitted for the sake of simplicity. Moreover, notice that, in system (24), only the variables $hu_0$ and $h\pi$ change with time. Then, using that $h_i$ is constant in this system, we freeze it at each time step when we solve the pressure system. Therefore, for first and second order schemes, taking into account that the mid point rule is applied to approximate the numerical source term, the equations corresponding to the variables $u_0$ and $\pi$ read as follow

$$ \left\{ \begin{aligned} (u_0)_i'(t) &= -\frac{1}{h_i(t_0)\Delta x} \left( \pi_{i+1/2}^*(t) - \pi_{i-1/2}^*(t) - \pi_i^{e,t_0}(x_{i+1/2}) + \pi_i^{e,t_0}(x_{i-1/2}) \right), \\ \pi_i'(t) &= -\frac{a^2}{h_i(t_0)\Delta x} \left( u_{0,i+1/2}^*(t) - u_{0,i-1/2}^*(t) - u_{0,i}^{e,t_0}(x_{i+1/2}) + u_{0,i}^{e,t_0}(x_{i-1/2}) \right). \end{aligned} \right. \qquad (28) $$

Denoting as $\vec{w} = \pi + au_0$ and $\overleftarrow{w} = \pi - au_0$ the Riemann invariants, system (28) can be rewritten as two transport equations:

$$ \left\{ \begin{aligned} \dot{\vec{w}}_i'(t) &= -\frac{a}{h_i(t_0)\Delta x} \left( \vec{w}_{i+1/2}(t) - \vec{w}_{i-1/2}(t) - \overleftarrow{w}_i^{e,t_0}(x_{i+1/2}) + \overleftarrow{w}_i^{e,t_0}(x_{i-1/2}) \right) \\ \dot{\overleftarrow{w}}_i'(t) &= \frac{a}{h_i(t_0)\Delta x} \left( \overleftarrow{w}_{i+1/2}(t) - \overleftarrow{w}_{i-1/2}(t) - \overleftarrow{w}_i^{e,t_0}(x_{i+1/2}) + \overleftarrow{w}_i^{e,t_0}(x_{i-1/2}) \right), \end{aligned} \right. \qquad (29) $$

where $\vec{w}_{i+1/2}(t) \approx \vec{w}(x_{i+1/2}, t)$ and $\overleftarrow{w}_{i+1/2}(t) \approx \overleftarrow{w}(x_{i+1/2}, t)$ are the numerical fluxes at the intercells, and will be computed by means of a reconstruction operator. In addition, the variables corresponding to the pressure and the velocity can be easily obtained from the following relations:

$$ \pi = \frac{\vec{w} + \overleftarrow{w}}{2}, \quad u_0 = \frac{\vec{w} - \overleftarrow{w}}{2a}. \qquad (30) $$
Then, the idea is to first solve (29) and then we can use (30) to define $\pi_{i+1/2}^*(t)$ and $u_{0,i+1/2}^*(t)$ as

$$ \pi_{i+1/2}^*(t) = \frac{P_{i, \vec{w}}(x_{i+1/2}, t) + P_{i+1, \vec{w}}(x_{i+1/2}, t)}{2}, \quad (31) $$

$$ u_{0,i+1/2}^*(t) = \frac{P_{i, \vec{w}}(x_{i+1/2}, t) - P_{i+1, \vec{w}}(x_{i+1/2}, t)}{2a}. \quad (32) $$

**Remark 2** As far as the stability of the scheme is concerned, for the pressure step the explicit schemes must satisfy the following stability restriction:

$$ \Delta t_P \le \frac{\Delta x \max_i \{h_i\}}{a}, \quad (33) $$

while the one for the transport step the restriction can be written as

$$ \Delta t_T \le \frac{\Delta x}{2 \max\{\max_i\{|\lambda_{1,i}^T|\}, \max_j\{|\lambda_{2,j}^T|\}\}}. \quad (34) $$

Of course, we will choose the maximum $\Delta t$ satisfying both restrictions. The main advantage of using semi-implicit schemes is that they do not require restriction (33), which is specially interesting in low Froude number situations.

**Remark 3** In practice, the value of the constant $a$ will be computed as

$$ a = \max_i \{h_i \sqrt{gh_i}\}. $$

## 3.1 Fully well-balanced reconstruction operator

We now illustrate the process for reconstructing the variables, presented in a general form for a variable $X$, which could represent $\vec{w}$, $\vec{w}$, $h$, or $hu_j$, for $j = 0, \dots, N$.

In order to reconstruct a variable $X$ at a given time $t_0$, we apply the strategy described in [21], consisting on, at each cell, computing the stationary solution $X_i^{e,t_0}(x)$ and applying a standard reconstruction operator $Q$ to the fluctuations with respect to the stationary solution, giving

$$ P_{i,X}^{t_0}(x) = X_i^{e,t_0}(x) + Q_{i,X}^{t_0}(x). \quad (35) $$

Following [38], the reconstruction operator at time $t \in [t_0, t_0 + \Delta t]$ is expressed as the sum of the well-balanced reconstruction operator at time $t_0$ and a standard reconstruction operator applied to the time fluctuations. Thus, for any variable $X$, the reconstruction operator is defined as:

$$ \begin{aligned} P_{i,X}(x,t) &= P_{i,X}^{t_0}(x) + \tilde{Q}_{i,X}(x,t) \\ &= X_i^{t_0,e}(x) + Q_{i,X}^{t_0}(x) + \tilde{Q}_{i,X}(x,t), \quad t \in [t_0, t_0 + \Delta t], x \in [x_{i-1/2}, x_{i+1/2}). \end{aligned} $$
Here, $P_{i,X}^{t_0}$ and $Q_{i,X}^{t_0}$ denote the reconstruction operators defined in (35) for variable $X$ at time $t_0$ and $\tilde{Q}_{i,X}(x,t)$ is a reconstruction operator based on time fluctuations, given by

$$ \tilde{Q}_{i,X}(x,t) = \tilde{Q}_i(x; \{X_j^{t,f}\}_{j \in S_i}), \text{ where } X_j^{t,f} = X_j(t) - X_j^{t_0}, \quad j \in S_i. $$

Note that stationary solutions are recalculated at each initial time $t_0$.

Thus, the first-order well-balanced reconstruction operator is

$$ P_{i,X}^{o1}(x,t) = X_i^{e,t_0}(x) + X_i(t) - X_i^{e,t_0}(x_i), \quad (36) $$

while for second order schemes, it is given by:

$$ P_{i,X}^{o2}(x,t) = X_i^{e,t_0}(x) - X_i^{e,t_0}(x_i) + \Delta X_i^{t,f}(x-x_i) + X_i(t) + \Delta X_i^{t,f}(x-x_i). \quad (37) $$

Here, employing the avg or harmod function limiter from [49], we set

$$ \Delta X_i^{t_0,f} = \frac{1}{\Delta x} \left( \phi_{i+}^{t_0}(X_i^{t_0,f} - X_{i-1}^{t_0,f}) + \phi_{i-}^{t_0}(X_{i+1}^{t_0,f} - X_i^{t_0,f}) \right), $$

where

$$ \phi_{i-}^{t_0} = \begin{cases} \frac{|d_{i-}|}{|d_{i-}|+|d_{i+}|} & \text{if } |d_{i-}| + |d_{i+}| > 0, \\ 0 & \text{otherwise,} \end{cases} $$

and

$$ \phi_{i+}^{t_0} = \begin{cases} \frac{|d_{i+}|}{|d_{i-}|+|d_{i+}|} & \text{if } |d_{i-}| + |d_{i+}| > 0, \\ 0 & \text{otherwise,} \end{cases} $$

where $d_{i-} = X_i^{t_0,f} - X_{i-1}^{t_0,f}$ and $d_{i+} = X_{i+1}^{t_0,f} - X_i^{t_0,f}$, with

$$ X_j^{t_0,f} = X_j^{t_0} - X_i^{e,t_0}(x_j) $$

for a given cell $i$. Moreover, we define

$$ \Delta X_i^{t,f} = \frac{1}{\Delta x} \left( \tilde{\phi}_{i+}^{t_0}(X_i^{t,f} - X_{i-1}^{t,f}) + \tilde{\phi}_{i-}^{t_0}(X_{i+1}^{t,f} - X_i^{t,f}) \right), $$

where $\tilde{\phi}_{i\pm}^{t_0} = \phi_{i\pm}^{t_0}$ and $X_i^{t,f} = X_i(t) - X_i^{t_0}$. Note that both $\Delta X_i^{t_0,f}$ and $\Delta X_i^{t,f}$ utilize the same limiters calculated at time $t_0$ to avoid introducing nonlinearities.

## 3.2 First order schemes

We will now describe two different versions of our fully well-balanced first order schemes: one in which both the pressure and the transport systems are solved explicitly, and another one in which the pressure system is solved implicitly, while the transport one is solved explicitly.

Let us assume that the cell averages at time $t_n$, denoted by $U_i^n$, are known. Our goal is to compute the cell averages at the next time step, $U_i^{n+1}$. To accomplish this, we begin by solving the pressure system over a time step $\Delta t$, yielding
an intermediate solution, which we denote as $U_i^{n+1-}$. Then, this intermediate solution $U_i^{n+1-}$ serves as the initial condition for solving the transport system, again over a time step $\Delta t$, which provides the final approximation for $U_i^{n+1}$.

**Remark 4** Note that we could also solve the transport system first, and then the pressure one. However, we have only considered the previously described version of the scheme since in previous works it has been observed that this one performs best in terms of stability (see [16]).

### 3.2.1 Explicit scheme

For the explicit scheme, following the semi-discrete scheme (24), the pressure system updates the variable $hu_0$ as follows:

$$ (hu_0)_i^{n+1-} = (hu_0)_i^n - \frac{\Delta t}{\Delta x} \left( \pi_{i+1/2}^{*,n} - \pi_{i-1/2}^{*,n} - \pi_i^{e,n}(x_{i+1/2}) + \pi_i^{e,n}(x_{i-1/2}) \right), \quad (38) $$

where $\pi_{i\pm 1/2}^{*,n}$ is obtained from the expression (31) using the first order reconstruction operator (36) to get the values of the Riemann invariants at the intercells. Note that, within this step $h$ and $(hu_j)$ for $j = 1, \dots, N$ do not change, that is, $h_i^{n+1-} = h_i^n$ and $(hu_j)_i^{n+1-} = (hu_j)_i^n$ for $j > 0$.

Once we have $U_i^{n+1-}$, we can perform the transport step considering (25)

$$ \begin{aligned} U_i^{n+1} = & U_i^{n+1-} - \frac{\Delta t}{\Delta x} \left( F_{i+1/2}^{n+1-} - F(U_i^{e,t_{n+1-}}(x_{i+1/2})) - F_{i-1/2}^{n+1-} + F(U_i^{e,t_{n+1-}}(x_{i-1/2})) \right) \\ & - \frac{\Delta t}{\Delta x} \left( B_{i+1/2-}^{n+1-} + B_{i-1/2+}^{n+1-} \right) \\ & - \frac{\Delta t}{\Delta x} B(U_i^{n+1-}) \left( U_{i+1/2-}^{n+1-} - U_i^{e,n+1-}(x_{i+1/2}) - U_{i-1/2+}^{n+1-} + U_i^{e,n+1-}(x_{i-1/2}) \right), \end{aligned} \quad (39) $$

where the reconstructed states are always computed using the well-balanced first order reconstruction operator (36), the terms $F_{i+1/2}^{n+1-}$ and $B_{i+1/2-}^{n+1-}$ are evaluated at the reconstructed states at time $t_{n+1-}$ and the derivative in the integral has been approximated using

$$ \frac{dP_{i,U}(x_i)}{dx} \approx \frac{P_{i,U}^{01}(x_{i+1/2}) - P_{i,U}^{01}(x_{i-1/2})}{\Delta x} = \frac{U_{i+1/2-} - U_{i-1/2+}}{\Delta x}. $$

### 3.2.2 Semi-implicit scheme

In the semi-implicit case, the only difference with the explicit one happens in the pressure step. Here, the variable $hu_0$ is updated using

$$ (hu_0)_i^{n+1-} = (hu_0)_i^n - \frac{\Delta t}{\Delta x} \left( \pi_{i+1/2}^{*,n+1-} - \pi_{i-1/2}^{*,n+1-} - \pi_i^{e,n}(x_{i+1/2}) + \pi_i^{e,n}(x_{i-1/2}) \right). \quad (40) $$
In practice, this value is obtained as

$$ (hu_0)_i^{n+1-} = h_i^n \cdot \frac{\vec{w}_i^{n+1-} - \overleftarrow{w}_i^{n+1-}}{2a}, $$

where the variables $\vec{w}_i^{n+1-}$ and $\overleftarrow{w}_i^{n+1-}$ are given by system (29) which is discretised as follows

$$ \vec{w}_i^{n+1-} = \vec{w}_i^n - \frac{a\Delta t}{h_i^n \Delta x} \left( \vec{w}_{i+1/2}^{n+1} - \vec{w}_{i-1/2}^{n+1} - \vec{w}_i^{e,n}(x_{i+1/2}) + \vec{w}_i^{e,n}(x_{i-1/2}) \right), $$

and

$$ \overleftarrow{w}_i^{n+1-} = \overleftarrow{w}_i^n + \frac{a\Delta t}{h_i^n \Delta x} \left( \overleftarrow{w}_{i+1/2}^{n+1} - \overleftarrow{w}_{i-1/2}^{n+1} - \overleftarrow{w}_i^{e,n}(x_{i+1/2}) + \overleftarrow{w}_i^{e,n}(x_{i-1/2}) \right), $$

where (36) is used to compute $\vec{w}_{i\pm 1/2}^{n+1}$ and

$$ \vec{w}_i^{e,n}(x_{i\pm 1/2}) = \frac{1}{2}g(h_i^e)^2(x_{i\pm 1/2}) + au_{0,i}^e(x_{i\pm 1/2}). $$

### 3.3 Second order schemes

To achieve second order accuracy, we employ the Strang splitting method (see [47, 48, 54]) for time-stepping. Specifically, we proceed as follows:

1. Take a half-step for the transport system using a time step of $\Delta t/2$.

2. Advance the pressure system over a full time step $\Delta t$.

3. Conclude with another half-step for the transport system, again using a time step of $\Delta t/2$, to obtain the final values $U_i^{n+1}$.

This could also be written in terms of the approximate solution operators of the systems as

$$ U(x, t + \Delta t) = S_T^{\Delta t/2} \circ S_P^{\Delta t} \circ S_T^{\Delta t/2} (U(x, t)), \quad (41) $$

where the terms $S_T^{\Delta t}$ and $S_T^{\Delta t}$ denote the approximate solution operators over the interval $[t, t + \tau]$ for the pressure system and the transport one, respectively.

In each step, second order spatial approximations are required. While each sub-step is first order in time, second order accuracy overall is achieved through the Strang splitting method.

**Remark 5** As mentioned for the first order schemes, in this second order case, we could also consider another variant of the scheme:

$$ U(x, t + \Delta t) = S_P^{\Delta t/2} \circ S_T^{\Delta t} \circ S_P^{\Delta t/2} (U(x, t)). $$

However, following again the work [16], we have only considered the approach (41) due to stability reasons.

For the sake of readability, the explicit description of the discrete version of the second order schemes are not detailed here. The procedure in this case is analogous to that already described in the first order case using (41) as well as second order reconstructions given in (37).
# 4 Numerical experiments

In this section, we present several numerical experiments designed to evaluate the performance of our proposed schemes. In each experiment we consider open boundary conditions and gravitational acceleration fixed at $g = 9.812$. Unless stated otherwise, the number of moments is set to $N = 8$.

In the following tests, when we mention a CFL number greater than 1, we are specifically referring to the stability constraint related to the sound speed (33). Let us remark that our semi-implicit methods also present a stability constraint linked to the transport step, which depends on the flow velocity (34). This means that we cannot increase the sound-speed-related CFL indefinitely, as it remains bounded by the transport CFL limit.

## 4.1 Test 1. Well-balanced property check: lake-at-rest

In order to check the well-balanced property, let us consider a lake-at-rest initial condition given in [43]. The bottom topography is given by

$$z(x) = \begin{cases} 2 - x^2 & \text{if } -0.5 \le x \le 0.5, \\ 1.75 & \text{otherwise.} \end{cases} \qquad (42)$$

The water depth is $h(x) = 3 - z(x)$ and a zero velocity profile is considered $u_j(x) = 0$, for $j = 0, 1, \dots, 8$.

The spatial domain is $[-1, 1]$ and the experiment is run until the final time $t = 0.5$, considering a 400-cell mesh. We take the CFL value to be 0.9 for the explicit schemes whereas the implicit methods use the CFL parameter set to 10. Table 1 shows the $L^1$-errors at $t = 0.5$ between the initial water-at-rest state and the final numerical solutions for the different schemes. For simplicity, in the previous table, only the variables $h$ and $hu_0$ are shown: the errors for the other variables are exactly zero for this test. As expected, all methods give machine accuracy.

<table><caption>Table 1: Test 1. Well-balanced property check: lake-at-rest. $L^1$-errors at time $t = 0.5$. The errors of all schemes are machine precision, indicating the well-balanced property.</caption><thead><tr><th>Variable</th><th>EXP O1</th><th>EXP O2</th><th>IMP O1</th><th>IMP O2</th></tr></thead><tbody><tr><td>$h$</td><td>$1.33 \cdot 10^{-15}$</td><td>$4.44 \cdot 10^{-16}$</td><td>$3.10 \cdot 10^{-15}$</td><td>$2.22 \cdot 10^{-15}$</td></tr><tr><td>$hu_0$</td><td>$7.48 \cdot 10^{-15}$</td><td>$2.06 \cdot 10^{-14}$</td><td>$7.60 \cdot 10^{-15}$</td><td>$7.60 \cdot 10^{-15}$</td></tr></tbody></table>

## 4.2 Test 2. Well-balanced property check: subcritical stationary flow with zero moments

Let us check the fully well-balanced property, considering now an initial condition which corresponds to a subcritical stationary solution given in [30]. The
<figure><img src="image_3.png" /><figcaption>Figure 1: Test 2. Well-balanced property check: subcritical stationary flow with zero moments. Initial condition: free surface and bottom (left) and velocity $u_0$ (right).</figcaption></figure>

bottom topography is considered to be

$$z(x) = \begin{cases} \frac{1}{4}(\cos((x + \frac{1}{2})5\pi) + 1) & \text{if } 1.3 \le x \le 1.7, \\ 0 & \text{otherwise.} \end{cases} \qquad (43)$$

We consider the subcritical steady state characterised by the constants

$$C_1 = 3.5, \quad C_2 = 21.15525, \quad C_k = 0, \ k \ge 3, \qquad (44)$$

in (9)-(10)-(11), which is shown in Figure 1. The spatial domain is $[0, 3]$ and the experiment is run until $t = 0.5$, considering a 400-cell mesh. For the explicit and the implicit schemes, we take the CFL value to be 0.9 and 1.26, respectively. Recall that, in practice, the CFL number is restricted by the transport step of the scheme in the implicit case. Since the Froude number varies between 0.4 and 0.78, no larger CFL condition is possible in this test.

Table 2 shows the $L^1$-errors at time $t = 0.5$ between the initial subcritical initial condition and the final numerical solutions for the different schemes. As expected, all methods give machine accuracy.

<table><caption>Table 2: Test 2. Well-balanced property check: subcritical stationary flow with zero moments. $L^1$-errors at time $t = 0.5$ for the subcritical initial condition given by the constants (44). The errors of all schemes are machine precision, indicating the well-balanced property.</caption><thead><tr><th>Variable</th><th>EXP O1</th><th>EXP O2</th><th>IMP O1</th><th>IMP O2</th></tr></thead><tbody><tr><td>$h$</td><td>$2.90 \cdot 10^{-14}$</td><td>$3.41 \cdot 10^{-14}$</td><td>$3.66 \cdot 10^{-14}$</td><td>$3.24 \cdot 10^{-14}$</td></tr><tr><td>$hu_0$</td><td>$6.26 \cdot 10^{-14}$</td><td>$1.82 \cdot 10^{-14}$</td><td>$2.35 \cdot 10^{-14}$</td><td>$2.08 \cdot 10^{-14}$</td></tr></tbody></table>

As mentioned before, the numerical scheme described here is specially designed for situations where the Froude number is low, as it then makes sense
<table><caption>Table 3: Test 2. Subcritical stationary flow with zero moments. CPU runtime in s for the different schemes using a 400-cell mesh. The CFL is equal to 0.9 for the explicit schemes and 10 for the implicit ones. As expected by the relative CFL numbers, we observe a speedup of about 10 for the respective implicit schemes.</caption><thead><tr><th>EXP O1</th><th>EXP O2</th><th>IMP O1</th><th>IMP O2</th></tr></thead><tbody><tr><td>23.58</td><td>90.51</td><td>2.79</td><td>8.28</td></tr></tbody></table>

to make implicit the pressure part of the system in order to avoid an overly restrictive stability condition. Let us now consider a subcritical steady state scenario with a lower Froude number, characterized by the following constants:

$$C_1 = 0.5, \quad C_2 = 21.15525, \quad C_k = 0, \quad k \ge 3, \qquad (45)$$

Again, the final time is $t = 0.5$. In this case the Froude number is between 0.05 and 0.07, which allows us to consider CFL up to 10 for the implicit scheme. This initial condition is preserved up to machine precision as well, but we are now interested in the efficiency of the scheme. Table 3 shows the different CPU times needed for the explicit and implicit schemes. Notice that the first order implicit scheme shows a speedup of 8.45 compared to the explicit one, while for the second order case the speedup of the implicit scheme is 10.93 compared to the explicit one. These numbers are close to the relative increase in the CFL number of 11.1.

## 4.3 Test 3. Well-balanced property check: subcritical stationary flow with non-zero moments

Let us consider an initial condition which corresponds to a stationary solution in the subcritical regime with non-zero moments inspired by an experiment given in [43]. The bottom topography is given again by (43), the spatial domain is $[0, 3]$ and the final time is $t = 0.5$, considering a 400-cell mesh. We take the CFL value to be 0.9 for the explicit schemes and 9.15 for the implicit ones (which is the maximum allowed for this particular solution). We consider the subcritical steady state with non-zero moments characterised by the constants:

$$C_1 = 0.5, \quad C_2 = 21.15525, \quad C_k = 0.005, \quad k \ge 3, \qquad (46)$$

which is shown in Figure 2. The Froude number is small, between 0.05 and 0.07. Figure 3 shows the velocity profile at $x = 1.5$. Table 4 shows the $L^1$-errors at $t = 0.5$ between the initial subcritical state and the final numerical solutions for the different schemes. As expected, all methods give machine accuracy. The CPU times for the different schemes are presented in Table 5 with the purpose of analysing the efficiency of the methods. Notice that the first and second order implicit schemes show a speedup of 9.4 and 9.5, respectively, over the explicit ones. This corresponds directly to the increase of CFL number by a factor of about 10.
<figure><img src="image_4.png" /><figcaption>Figure 2: Test 3. Well-balanced property check: subcritical steady state with non-null moments. Initial condition: free surface $\eta$ (up-left), $u_0$ (up-right) and $u_j$, $j \ge 1$ (down).</figcaption></figure>
<table><caption>Table 5: Test 3. Subcritical stationary flow with non-zero moments. CPU times for the different schemes using a 400-cell mesh. The CFL is equal to 0.9 for the explicit schemes and 9.15 for the implicit ones. Again, the speedup is about 10, in line with the increase of the CFL number.</caption><thead><tr><th>EXP O1</th><th>EXP O2</th><th>IMP O1</th><th>IMP O2</th></tr></thead><tbody><tr><td>23.07</td><td>74.85</td><td>2.46</td><td>7.86</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 3: Test 3. Subcritical stationary flow with non-zero moments: velocity profile at $x = 1.5$.</figcaption></figure>

<table><caption>Table 4: Test 3. Well-balanced property check: subcritical stationary flow with non-zero moments. $L^1$-errors at time $t = 0.5$ for the subcritical initial condition given by the constants (46). The errors of all schemes are machine precision, indicating the well-balanced property.</caption><thead><tr><th>Variable</th><th>EXP O1</th><th>EXP O2</th><th>IMP O1</th><th>IMP O2</th></tr></thead><tbody><tr><td>$h$</td><td>2.66 · 10<sup>-14</sup></td><td>1.68 · 10<sup>-14</sup></td><td>3.46 · 10<sup>-14</sup></td><td>2.02 · 10<sup>-14</sup></td></tr><tr><td>$q_0$</td><td>8.03 · 10<sup>-14</sup></td><td>3.73 · 10<sup>-14</sup></td><td>1.46 · 10<sup>-13</sup></td><td>2.99 · 10<sup>-14</sup></td></tr><tr><td>$hu_j, j \ge 1$</td><td>4.78 · 10<sup>-16</sup></td><td>3.36 · 10<sup>-16</sup></td><td>7.49 · 10<sup>-16</sup></td><td>3.22 · 10<sup>-16</sup></td></tr></tbody></table>

## 4.4 Test 4. Second order accuracy check

This experiment is devoted to check the accuracy of the explicit and implicit second order schemes. Let us consider a perturbation of the subcritical stationary solution $U^{0,*}$ considered in Test 3, characterised by the constants (46), with
the bottom (43). The initial perturbation is given by

$$U^0(x) = U^{0,*}(x) + (10^{-4}e^{-200(x-2)^2}, 0, \dots, 0). \quad (47)$$

The spatial domain is $[0, 3]$ and the final time is $t = 0.1$. We consider four grids consisting of 25, 50, 100 and 200 uniform points in order to compute the errors, comparing with a reference solution, which is computed with a mesh of 800 cells using the same numerical scheme. The CFL number is set to 0.9 for explicit schemes, whereas a CFL value equal to 2 is considered for implicit ones. The results are shown in Tables 6 and 7: the expected second order of accuracy is achieved in all variables.

<table><caption>Table 6: Test 4. Second order accuracy check. Explicit scheme. CFL= 0.9.</caption><thead><tr><th rowspan="2">No. of cells</th><th colspan="2">h</th><th colspan="2">u<sub>0</sub></th><th colspan="2">u<sub>j</sub>, j ≥ 1</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>25</td><td>7.08e-03</td><td>-</td><td>1.36e-03</td><td>-</td><td>3.54e-05</td><td>-</td></tr><tr><td>50</td><td>2.09e-03</td><td>1.76</td><td>3.07e-04</td><td>2.14</td><td>1.05e-05</td><td>1.76</td></tr><tr><td>100</td><td>6.08e-04</td><td>1.78</td><td>9.14e-05</td><td>1.75</td><td>3.05e-06</td><td>1.78</td></tr><tr><td>200</td><td>1.37e-04</td><td>2.15</td><td>2.17e-05</td><td>2.07</td><td>6.87e-07</td><td>2.15</td></tr></tbody></table>

<table><caption>Table 7: Test 4. Second order accuracy check. Implicit scheme. CFL=2.</caption><thead><tr><th rowspan="2">No. of cells</th><th colspan="2">h</th><th colspan="2">u<sub>0</sub></th><th colspan="2">u<sub>j</sub>, j ≥ 1</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>25</td><td>7.08e-03</td><td>-</td><td>1.36e-03</td><td>-</td><td>3.54e-05</td><td>-</td></tr><tr><td>50</td><td>2.10e-03</td><td>1.75</td><td>3.16e-04</td><td>2.10</td><td>1.05e-05</td><td>1.75</td></tr><tr><td>100</td><td>6.14e-04</td><td>1.77</td><td>1.04e-04</td><td>1.61</td><td>3.08e-06</td><td>1.77</td></tr><tr><td>200</td><td>1.41e-04</td><td>2.13</td><td>2.93e-05</td><td>1.82</td><td>7.08e-07</td><td>2.12</td></tr></tbody></table>

## 4.5 Test 5. Perturbation of a subcritical steady state

We want to perform a test, similar to the one proposed in Test 4, but now considering a more physical initial condition. With that aim, we consider now $N = 2$ moments, so that we have a parabolic velocity profile. The spatial domain is $[0, 3]$, the bottom topography is again (43) and a final time $t = 0.1$ is considered in order to observe the propagation of the perturbation before it has left the domain. As initial condition, we consider a perturbation of a subcritical stationary solution. That subcritical steady state $U^{0,*}$ is characterised by the following constants:

$$C_1 = 0.5, \quad C_2 = 21.15525, \quad C_3 = -0.005, \quad C_4 = -0.001, \quad (48)$$

and the initial perturbation is then given by

$$U^0(x) = U^{0,*}(x) + (10^{-4}e^{-200(x-2)^2}, 0, \dots, 0), \quad (49)$$
and shown in Figure 4.

In Figure 5, we have plotted the solution obtained with the implicit schemes considering different CFL values. We have zoomed in the free surface to see the differences between the schemes, observing, as expected that the first order schemes are, in general, more diffusive, and that the diffusion is increased when the CFL is bigger. This behaviour is also observed in the variable $hu_0$, while no major differences can be perceived in $hu_1$ and $hu_2$. As previously mentioned, the velocity profile is now parabolic. As an example, Figure 6 shows the vertical profile at the final time at point $x = 2.25$ using the second order implicit scheme with CFL=5.

<figure><img src="image_6.png" /><figcaption>Figure 4: Test 5. Perturbation of a subcritical steady state. Initial condition: free surface $\eta$ (up-left), $hu_0$ (up-right), $hu_1$ (down-left) and $hu_2$ (down-right).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: Test 5. Perturbation of a subcritical steady state. Numerical solution for first and second order implicit schemes using different CFL values: free surface $\eta$ (up-left), zoom in the free surface (up-right), $hu_0$ (middle-left), $hu_1$ (middle-right) and $hu_2$ (down).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Test 5. Perturbation of a subcritical stationary flow: velocity profile at $x = 2.25$ using the second order implicit scheme with CFL=5.</figcaption></figure>

## 4.6 Test 6. Dam-break problem

Even though the aim of our schemes is to use them in low Froude number situations, since it is in that case when they outperform explicit schemes, we are now going to check that in fact they can also be used in other cases, such as in dam-break situations. In order to do so, we consider a problem similar to one shown in [44]. The bottom topography is considered to be flat, i.e., $z(x) = 0$ and an initial condition given by:

$$h^0(x) = \begin{cases} 2 & \text{if } x \le 0, \\ 1 & \text{otherwise,} \end{cases} \qquad (50)$$

$$u_0^0(x) = 0.25; u_1^0(x) = -0.005; u_8^0(x) = 0.005; u_j^0(x) = 0, j = 2, \dots, 7.$$

The spatial domain is [-0.4, 0.4] and the final time is $t = 0.1$.

Figures 7 and 8 display the solutions provided by the methods of order 1 and 2 at $t = 0.01$ and $t = 0.1$ for $\eta$, the velocity $u_0$ and first and last moments $u_1$ and $u_8$. The CFL value is set to 0.9 for the explicit schemes, whereas the values 1 and 2 (the second being the maximum allowed) are taken for the implicit schemes.

The images that depict moments can be reinterpreted as a vertical profile for velocity. As an example, in Figure 9 we present the profiles obtained at the points $x = 0$ and $x = 0.15$ for the second order implicit scheme with a CFL number equal to 2 at the final time.
<figure><img src="image_9.png" /><figcaption>(a) Free surface $\eta$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Velocity $u_0$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) First moment $u_1$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(d) Last moment $u_8$</figcaption></figure>

Figure 7: Test 6. Standard dam-break test. Zoom of the variables $\eta$, $u_0$, $u_1$ and $u_8$ at time $t = 0.01$.
<figure><img src="image_13.png" /><figcaption>(a) Free surface $\eta$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Velocity $u_0$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) First moment $u_1$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) Last moment $u_8$</figcaption></figure>

Figure 8: Test 6. Standard dam-break test. Variables $\eta$, $u_0$, $u_1$ and $u_8$ at time $t = 0.1$.
<figure><img src="image_17.png" /><figcaption>Figure 9: Test 6. Standard dam-break test: velocity profile at $x = 0$ and $x = 0.15$ at the final time $t = 0.1$ using the second order implicit scheme with CFL= 2.</figcaption></figure>

# 5 Concluding remarks

This work considers the development of numerical schemes for the Shallow Water Linearized Moment Equations (SWLME), which allow to obtain a vertical distribution of the horizontal velocity. This is in contrast to the standard Shallow Water Equations (SWE) where the horizontal velocity is constant in the vertical direction. The improved description of the velocity profile is useful in many practical applications, where a constant vertical profile could be very restrictive. Compared to other Shallow Water Moment Equations (SWME), the SWLME assume that higher order moments are small, which results in an explicit algebraic expression for steady states and explicit computation of the eigenvalues.

SWLME have been solved by means of a semi-implicit first and second order exactly fully well-balanced method. A splitting approach allowed to separate acoustic and material phenomena. This is especially interesting for the case of low Froude regimes, where the CFL restriction due to the pressure term is too restrictive when compared to the velocity stability restriction. Indeed, using an implicit approach for the pressure system made it possible to use a larger CFL condition and a significant speed-up comparable to this increase is obtained. The resulting scheme is therefore efficient in such situations and exactly preserve the stationary solutions, as shown in the several numerical experiments that have been performed.

# Acknowledgments

This work is supported by proyect PID2022-137637NB-C21 funded by MCIN/AEI/10.13039/501100011033/ and ERDF A way of making Europe, and proyect PDC2022-133663-C21 funded by MCIN/AEI/10.13039/501100011033 and European Union NextGenerationEU/PRTR.
References

[1] E. Audusse et al. “Approximation of the hydrostatic Navier–Stokes system for density stratified flows by a multilayer model: Kinetic interpretation and numerical solution”. In: *Journal of Computational Physics* 230.9 (2011), pages 3453–3478. ISSN: 0021-9991. DOI: 10.1016/j.jcp.2011.01.042.

[2] E. Audusse, C. Chalons, and P. Ung. “A simple well-balanced and positive numerical scheme for the shallow-water system”. In: *Communications in Mathematical Sciences* 13 (2015), pages 1317–1332. DOI: 10.4310/CMS.2015.v13.n5.a11.

[3] E. Audusse et al. “A Fast and Stable Well-Balanced Scheme with Hydrostatic Reconstruction for Shallow Water Flows”. In: *SIAM Journal on Scientific Computing* 25.6 (2004), pages 2050–2065. DOI: 10.1137/S1064827503431090.

[4] E. Audusse et al. “A multilayer Saint-Venant system with mass exchanges for Shallow Water flows. Derivation and numerical validation”. In: *ESAIM Mathematical Modelling and Numerical Analysis* 45 (2009). DOI: 10.1051/m2an/2010036.

[5] M. Bauerle et al. *On the rotational invariance and hyperbolicity of shallow water moment equations in two dimensions*. 2023. arXiv: 2306.07202 [math.NA].

[6] C. Berthon and F. Foucher. “Efficient well-balanced hydrostatic upwind schemes for shallow-water equations”. In: *Journal of Computational Physics* 231.15 (2012), pages 4993–5015. ISSN: 0021-9991. DOI: https://doi.org/10.1016/j.jcp.2012.02.031.

[7] C. Berthon and V. Michel-Dansac. “A simple fully well-balanced and entropy preserving scheme for the shallow-water equations”. In: *Applied Mathematics Letters* 86 (2018), pages 284–290. ISSN: 0893-9659. DOI: 10.1016/j.aml.2018.07.013.

[8] C. Berthon and C. Chalons. “A fully well-balanced, positive and entropy-satisfying Godunov-type method for the shallow-water equations”. In: *Mathematics of Computation* 85.299 (2015), pages 1281–1307. DOI: 10.1090/mcom3045.

[9] C. Berthon et al. “Fully well-balanced, positive and simple approximate Riemann solver for shallow water equations”. In: *Bulletin of the Brazilian Mathematical Society, New Series* 47.1 (2016), pages 117–130. DOI: 10.1007/s00574-016-0126-1.

[10] G. Bispen et al. “IMEX Large Time Step Finite Volume Methods for Low Froude Number Shallow Water Flows”. In: *Communications in Computational Physics* 16.2 (2014), pages 307–347. DOI: 10.4208/cicp.040413.160114a.
[11] L. Bonaventura et al. "Multilayer shallow water models with locally variable number of layers and semi-implicit time discretization". In: *Journal of Computational Physics* 364 (2018), pages 209-234. DOI: 10.1016/j.jcp.2018.03.017.

[12] F. Bouchut. *Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources*. Frontiers in Mathematics. Basel: Birkhäuser Verlag, 2004, pages viii+135. ISBN: 3-7643-6665-6.

[13] R. Bourgeois et al. "Recasting an operator splitting solver into a standard finite volume flux-based algorithm. The case of a Lagrange-projection-type method for gas dynamics". In: *Journal of Computational Physics* 496 (2024), page 112594. ISSN: 0021-9991. DOI: https://doi.org/10.1016/j.jcp.2023.112594.

[14] S. Busto and M. Dumbser. "A staggered semi-implicit hybrid finite volume / finite element scheme for the shallow water equations at all Froude numbers". In: *Applied Numerical Mathematics* 175 (2022), pages 108-132. DOI: 10.1016/j.apnum.2022.02.005.

[15] C. Caballero-Cárdenas et al. "Implicit and implicit-explicit Lagrange-projection finite volume schemes exactly well-balanced for 1D shallow water system". In: *Applied Mathematics and Computation* 443 (2023), page 127784. DOI: 10.1016/j.amc.2022.127784.

[16] C. Caballero-Cárdenas et al. "A Semi-Implicit Fully Exactly Well-Balanced Relaxation Scheme for the Shallow Water System". In: *SIAM Journal on Scientific Computing* 46 (2024), A2503-A2527. DOI: 10.1137/23M1621289.

[17] F. N. Cantero-Chinchilla, O. Castro-Orgaz, and A. A. Khan. "Depth-integrated nonhydrostatic free-surface flow modeling using weighted-averaged equations". In: *International Journal for Numerical Methods in Fluids* 87.1 (2018), pages 27-50. ISSN: 1097-0363. DOI: 10.1002/fld.4481.

[18] M. Castro, C. Chalons, and T. Morales de Luna. "A Fully Well-Balanced Lagrange-Projection-Type Scheme for the Shallow-Water Equations". In: *SIAM Journal on Numerical Analysis* 56 (2018), pages 3071-3098. DOI: 10.1137/17M1156101.

[19] M. Castro, J. Macías, and C. Parés. "A Q-scheme for a class of systems of coupled conservation laws with source term. Application to a two-layer 1-D shallow water system". In: *M2AN Math. Model. Numer. Anal.* 35.1 (2001), pages 107-127. ISSN: 0764-583X.

[20] M. Castro et al. "Numerical Treatment of the Loss of Hyperbolicity of the Two-Layer Shallow-Water System". In: *Journal of Scientific Computing* 48 (2011), pages 16-40. DOI: 10.1007/s10915-010-9427-5.

[21] M. J. Castro and C. Parés. "Well-Balanced High-Order Finite Volume Methods for Systems of Balance Laws". In: *Journal of Scientific Computing* 82.2 (2020). DOI: 10.1007/s10915-020-01149-5.
[22] V. Casulli and R. A. Walters. "An unstructured grid, three-dimensional model based on the shallow water equations". In: *International Journal for Numerical Methods in Fluids* 32.3 (2000), pages 331-348. DOI: 10.1002/(SICI)1097-0363(20000215)32:3<331::AID-FLD941>3.0.CO;2-C.

[23] V. Casulli and E. Cattani. "Stability, accuracy and efficiency of a semi-implicit method for three-dimensional shallow water flow". In: *Computers & Mathematics with Applications* 27.4 (1994), pages 99-112. DOI: 10.1016/0898-1221(94)90059-0.

[24] C. Chalons, M. Girardin, and S. Kokh. "An All-Regime Lagrange-Projection Like Scheme for the Gas Dynamics Equations on Unstructured Meshes". In: *Communications in Computational Physics* 20.1 (2016), pages 188-233. ISSN: 1815-2406, 1991-7120. DOI: 10.4208/cicp.260614.061115a.

[25] C. Chalons et al. "A large time-step and well-balanced Lagrange-Projection type scheme for the shallow-water equations". In: *Communications in Mathematical Sciences* 15 (2016). DOI: 10.4310/CMS.2017.v15.n3.a9.

[26] A. Chertock et al. "Locally divergence-free well-balanced path-conservative central-upwind schemes for rotating shallow water MHD". In: *Journal of Computational Physics* 518 (2024), page 113300. ISSN: 0021-9991. DOI: https://doi.org/10.1016/j.jcp.2024.113300.

[27] A. Chinnayya, A.-Y. Leroux, and N. Seguin. "A Well-balanced Numerical Scheme for the Approximation of the Shallow-water Equations with Topography: The Resonance Phenomenon". In: *International Journal on Finite Volumes* 1 (2004).

[28] A. Del Grosso et al. "On Lagrange-Projection Schemes for Shallow Water Flows Over Movable Bottom with Suspended and Bedload Transport". In: *Numerical Mathematics: Theory, Methods and Applications* 41 (2023). DOI: 10.4208/nmtma.OA-2023-0082.

[29] A. Del Grosso et al. "On well-balanced implicit-explicit Lagrange-projection schemes for two-layer shallow water equations". In: *Applied Mathematics and Computation* 442 (2023), page 127702. DOI: 10.1016/j.amc.2022.127702.

[30] M. J. C. Díaz, J. A. López-García, and C. Parés. "High order exactly well-balanced numerical methods for shallow water systems". In: *Journal of Computational Physics* 246 (2013), pages 242-264.

[31] E. D. Fernández-Nieto, E. H. Koné, and T. C. Rebollo. "A Multilayer Method for the Hydrostatic Navier-Stokes Equations: A Particular Weak Solution". en. In: *Journal of Scientific Computing* 60.2 (2014), pages 408-437. ISSN: 0885-7474, 1573-7691. DOI: 10.1007/s10915-013-9802-0.

[32] E. Fernández-Nieto, M. Castro, and C. Madroñal. "On an Intermediate Field Capturing Riemann Solver Based on a Parabolic Viscosity Matrix for the Two-Layer Shallow Water System". In: *Journal of Scientific Computing* 48 (2011), pages 117-140. DOI: 10.1007/s10915-011-9465-7.
[33] E. Franck and L. Navoret. "Semi-implicit Two-Speed Well-Balanced Relaxation Scheme for Ripa Model". In: *Finite Volumes for Complex Applications IX - Methods, Theoretical Aspects, Examples*. Edited by R. Klöfkorn et al. Cham: Springer International Publishing, 2020, pages 735–743. ISBN: 978-3-030-43651-3.

[34] J. Garres-Díaz et al. "A general vertical decomposition of Euler equations: Multilayer-moment models". In: *Applied Numerical Mathematics* 183 (2023), pages 236–262. DOI: 10.1016/j.apnum.2022.09.004.

[35] J. Garres-Díaz et al. "Shallow Water Moment Models for Bedload Transport Problems". In: *Communications of Computational Physics* 30.3 (2021), pages 903–941.

[36] J.-F. Gerbeau and B. Perthame. "Derivation of viscous Saint-Venant system for laminar shallow water; Numerical validation". In: *Discrete and Continuous Dynamical Systems - B* 1.1 (2001), pages 89–102. ISSN: 1531-3492. DOI: 10.3934/dcdsb.2001.1.89.

[37] I. Gómez-Bueno et al. "Collocation Methods for High-Order Well-Balanced Methods for Systems of Balance Laws". In: *Mathematics* 9.15 (2021), page 1799. DOI: 10.3390/math9151799.

[38] I. Gómez-Bueno et al. "Implicit and semi-implicit well-balanced finite-volume methods for systems of balance laws". In: *Applied Numerical Mathematics* 184 (2023), pages 18–48. DOI: 10.1016/j.apnum.2022.09.016.

[39] E. Guerrero Fernández, M. J. Castro-Díaz, and T. Morales de Luna. "A Second-Order Well-Balanced Finite Volume Scheme for the Multilayer Shallow Water Model with Variable Density". en. In: *Mathematics* 8.5 (2020), page 848. DOI: 10.3390/math8050848.

[40] A. Harten, P. D. Lax, and B. van Leer. "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws". In: *Upwind and High-Resolution Schemes*. Berlin, Heidelberg: Springer Berlin Heidelberg, 1997, pages 53–79. ISBN: 978-3-642-60543-7. DOI: 10.1007/978-3-642-60543-7_4.

[41] Q. Huang, J. Koellermeier, and W.-A. Yong. "Equilibrium stability analysis of hyperbolic shallow water moment equations". In: *Mathematical Methods in the Applied Sciences* 45.10 (2022), pages 6459–6480. DOI: 10.1002/mma.8180.

[42] D. Iampietro et al. "A Mach-sensitive implicit-explicit scheme adapted to compressible multi-scale flows". In: *Journal of Computational and Applied Mathematics* 340 (2018). DOI: 10.1016/j.cam.2018.02.019.

[43] J. Koellermeier and E. Pimentel-García. "Steady states and well-balanced schemes for shallow water moment equations with topography". In: *Applied Mathematics and Computation* 427 (2022), page 127166. DOI: 10.1016/j.amc.2022.127166.
[44] J. Koellermeier and M. Rominger. "Analysis and numerical simulation of hyperbolic shallow water moment equations". In: *Communications in Computational Physics* 28.3 (2020), pages 1038–1084.

[45] J. Kowalski and M. Torrilhon. "Moment Approximations and Model Cascades for Shallow Flow". In: *Communications in Computational Physics* 25 (2019). DOI: 10.4208/cicp.OA-2017-0263.

[46] E. Macca. "Shock-Capturing methods: Well-Balanced Approximate Taylor and Semi-Implicit schemes". PhD thesis. 2022. DOI: 10.13140/RG.2.2.11244.85128.

[47] G. I. Marchuk. *Metody rasshchepleniya*. Russian. Moskva: Nauka, 1988. ISBN: 5-02-013783-9.

[48] G. I. Marchuk. "Splitting and alternating direction methods". In: (1990), pages 197–462. DOI: 10.1016/S1570-8659(05)80035-3.

[49] A. Marquina and S. Serna. "Shock-capturing schemes: high accuracy versus total-variation boundedness". In: *PAMM* 7.1 (2007), pages 1024101–1024102. ISSN: 1617-7061. DOI: 10.1002/pamm.200700672.

[50] V. Michel-Dansac et al. "A well-balanced scheme for the shallow-water equations with topography". In: *Computers & Mathematics with Applications* 72.3 (2016), pages 568–593.

[51] S. Noelle et al. "Well-balanced finite volume schemes of arbitrary order of accuracy for shallow water flows". In: *Journal of Computational Physics* 213.2 (2006), pages 474–499. DOI: 10.1016/j.jcp.2005.08.019.

[52] C. Parés. "Numerical methods for nonconservative hyperbolic systems: a theoretical framework." In: *SIAM Journal on Numerical Analysis* 44.1 (2006), pages 300–321. DOI: 10.1137/050628052.

[53] P. M. Steffler and J. Yee-Chung. "Depth averaged and moment equations for moderately shallow free surface flow". In: *Journal of Hydraulic Research* 31.1 (1993), pages 5–17. DOI: 10.1080/00221689309498856.

[54] G. Strang. "On the Construction and Comparison of Difference Schemes". In: *SIAM Journal on Numerical Analysis* 5.3 (1968), pages 506–517. DOI: 10.1137/0705041.

[55] R. Verbiest and J. Koellermeier. *Hyperbolic Axisymmetric Shallow Water Moment Equations*. 2023. arXiv: 2302.07952 [math.NA].

[56] C. Vreugdenhil. "Numerical Methods for Shallow Water Flow". In: *In VKI, Computational Fluid Dynamics, Volume 1 48 p (SEE N90-27989 22-34)* 1 (1990).

[57] X. C. Xi Chen and A. K. Alexander Kurganov. "A Well-Balanced Partial Relaxation Scheme for the Two-Dimensional Saint-Venant System". In: *Communications in Computational Physics* 33.5 (2023), pages 1432–1465. ISSN: 1815-2406. DOI: https://doi.org/10.4208/cicp.OA-2022-0319.