Christophe Berthon

Victor Michel-Dansac

A fully well-balanced hydrodynamic reconstruction

The present work focuses on the numerical approximation of the weak solutions of the shallow water model over a non-flat topography. In particular, we pay close attention to steady solutions with nonzero velocity. The goal of this work is to derive a scheme that exactly preserves these stationary solutions, as well as the commonly preserved lake at rest steady solution. To address this issue, we propose an extension of the well-known hydrostatic reconstruction. By appropriately defining the reconstructed states at the interfaces, any numerical flux function, combined with a relevant source term discretization, produces a well-balanced scheme that preserves both moving and non-moving steady solutions. This eliminates the need to construct specific numerical fluxes. Additionally, we prove that the resulting scheme is consistent with the homogeneous system on flat topographies, and that it reduces to the hydrostatic reconstruction when the velocity vanishes. To increase the accuracy of the simulations, we propose a linear well-balanced high-order procedure. Several numerical experiments demonstrate the effectiveness of the numerical scheme.

Introduction

In this work, we consider the numerical approximation of the weak solutions of the shallow water equations, given as follows:

     ∂ t h + ∂ x q = 0, ∂ t q + ∂ x q 2 h + 1 2 gh 2 = -gh∂ x Z, (1.1)
where h ≥ 0 is the water height and q ∈ R stands for the discharge. The given topography function Z is assumed to be smooth enough, and g > 0 is the gravity constant. For the sake of convenience in the notation, as long as h > 0, we define the velocity u as follows:

u = q h . (1.2)
In fact, defining the velocity in dry regions, when the water height h vanishes, requires special attention. In this paper, we impose the following conventions: lim h→0 q h = 0 and q 2 h 3 is bounded when h → 0.

(1.

3)

The first convention is commonly used when simulating dry areas. The second convention, however, is less usual. It is associated with the behavior of the Froude number in dry areas. * Université de Nantes, CNRS UMR 6629, Laboratoire de Mathématiques Jean Leray, 2 rue de la Houssinière, BP 92208, 44322 Nantes, France; email address: christophe.berthon@univ-nantes.fr; homepage: https://www.math.sciences.univ-nantes. fr/~berthon † Université de Strasbourg, CNRS, Inria, IRMA, F-67000 Strasbourg, France; email address: victor.michel-dansac@inria.fr; homepage: https://irma.math.unistra.fr/~micheldansac ‡ corresponding author 1

For the sake of clarity, we introduce the following condensed notation

W = h q , F (W ) =    q q 2 h + 1 2 gh 2    and S(W ) = 0 -gh∂ x Z , (1.4)
so that the shallow water model (1.1) rewrites

∂ t W + ∂ x F (W ) = S(W ). (1.5)
This system is endowed with an initial data

W (x, t = 0) = W 0 (x), (1.6)
where W 0 is defined according to prescribed physics. When deriving numerical schemes to approximate the weak solutions of (1.1), special attention must be given to the steady solutions, see [START_REF] Bermúdez | Upwind methods for hyperbolic conservation laws with source terms[END_REF][START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF]. Such solutions of major interest are governed by the system ∂ x F (W ) = S(W ). Solving this time-independent system yields smooth steady solutions that only depend on x, given by    q =: q 0 , q 2 2h 2 + g(h + Z) =: B 0 , (1.7) where q 0 and B 0 are given constant values and with a Froude number Fr 2 = q 2 gh 3 different from 1. Over the last three decades, a significant amount of research has been devoted to the derivation of numerical schemes able to exactly preserve the steady solutions. After pioneering work by Bermúdez and Vázquez [START_REF] Bermúdez | Upwind methods for hyperbolic conservation laws with source terms[END_REF], and next by Greenberg and Leroux [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF], it is now well-known that steady solutions play a crucial role when designing numerical discretizations. Indeed, small errors in the approximation of a steady solution may accumulate over time and render the simulation irrelevant. To avoid such an issue, well-balanced schemes have been introduced (see [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF], as well as [START_REF] Bermúdez | Upwind methods for hyperbolic conservation laws with source terms[END_REF] for the C-property). A scheme is said to be well-balanced if it exactly preserves at least some steady solutions. In [START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF], Gosse obtains a well-balanced scheme by solving the nonlinear Bernoulli equation (1.7) in each cell and at each time step. Next, in [START_REF] Jin | A steady-state capturing method for hyperbolic systems with geometrical source terms[END_REF], Jin derives an approximately well-balanced scheme, where accurately approximating the steady solutions, rather than exactly preserving them, results in good scheme behavior.

However, these approaches have several drawbacks. Indeed, exactly solving Bernoulli's equation is computationally costly, and the approximately well-balanced property may be insufficient in some cases. An interesting alternative was introduced by Audusse et al. in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF], where the hydrostatic reconstruction was designed to exactly preserve the lake at rest steady state. The lake at rest is governed by particular case of (1.7), where q 0 is fixed equal to 0, thus leading to h + Z =: H 0 , with H 0 a given constant free surface. The simplicity of this approach makes it very attractive.

Indeed, let us briefly recall the hydrostatic reconstruction. We introduce a space discretization made of cells (x i-1/2 , x i+1/2 ) of constant size ∆x > 0, such that x i+1/2 = x i-1/2 +∆x. We also introduce the time discretization t n+1 = t n + ∆t, where the time increment ∆t > 0 is restricted by a CFL-type condition. To define the hydrostatic reconstruction, the authors of [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] merely have to adopt a numerical flux function F W n i , W n i+1 that is consistent with the exact flux function F (W ) given by (1.4), i.e. F(W, W ) = F (W ). Next, at each interface x i+ 1 2 , they introduce the following reconstructed water heights:

h n i+ 1 2 ,-= max 0, h n i + Z i -Z i+ 1 2 and h n i+ 1 2 ,+ = max 0, h n i+1 + Z i+1 -Z i+ 1 2 ,
with Z i+ 1 2 = max(Z i , Z i+1 ), where

Z i = 1 ∆x x i+ 1 2 x i-1 2 Z(x)dx
is a finite volume discretization of the topography function Z. Equipped with such a reconstruction, with W n i a given approximation of W (x, t n ) over the cell (x i-1 2 , x i+ 1 2 ), the approximate solution of (1.1) at time t n+1 is given by

W n+1 i = W n i - ∆t ∆x F W n i+ 1 2 ,-, W n i+ 1 2 ,+ -F W n i-1 2 ,-, W n i-1 2 ,+ + ∆t S n i , (1.8)
where

W n i+ 1 2 ,-= h n i+ 1 2 ,- h n i+ 1 2 ,-u n i , W n i+ 1 2 ,+ = h n i+ 1 2 ,+ h n i+ 1 2 ,+ u n i+1 and u n i = (hu) n i h n i .
Concerning the source term S n i , it now defined as follows:

S n i =   0 g 2 (h n i+ 1 2 ,-) 2 -(h n i-1 2 ,+ ) 2 ∆x   .
Thanks to its simplicity and versatility, the hydrostatic reconstruction is frequently used (for instance, see [START_REF] Morales De Luna | Reliability of first order numerical schemes for solving shallow water system over abrupt topography[END_REF][START_REF] Natalini | A well-balanced numerical scheme for a one dimensional quasilinear hyperbolic model of chemotaxis[END_REF][START_REF] Chen | A New Hydrostatic Reconstruction Scheme Based on Subcell Reconstructions[END_REF][START_REF] Berthon | Improvement of the Hydrostatic Reconstruction Scheme to Get Fully Discrete Entropy Inequalities[END_REF], for applications and extensions) to get a scheme that exactly preserves, specifically, the lake at rest steady solution. Several other techniques have also been developed to yield lake at rest-preserving schemes (for instance, see the non-exhaustive list [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations[END_REF][START_REF] Xing | High-order finite volume WENO schemes for the shallow water equations with dry states[END_REF][START_REF] Berthon | Efficient well-balanced hydrostatic upwind schemes for shallow-water equations[END_REF][START_REF] Delestre | A 'well-balanced' finite volume scheme for blood flow simulation[END_REF][START_REF] Chen | A New Hydrostatic Reconstruction Scheme Based on Subcell Reconstructions[END_REF][START_REF] Li | Well-balanced discontinuous Galerkin methods with hydrostatic reconstruction for the Euler equations with gravitation[END_REF][START_REF] Duran | Energy-stable staggered schemes for the Shallow Water equations[END_REF]).

Yet, merely considering the lake at rest preservation may be insufficient in some simulations, see [START_REF] Xing | On the advantage of well-balanced schemes for moving-water equilibria of the shallow water equations[END_REF]. As a consequence, more recently, new approaches have been proposed to deal with still or moving steady states, which are described by (1.7) with q = 0 or q ̸ = 0. For a few examples of schemes that exactly preserve moving steady states, the reader is referred to [START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equation with moving water[END_REF][START_REF] Xing | Exactly well-balanced discontinuous Galerkin methods for the shallow water equations with moving water equilibrium[END_REF] for high-order schemes, to [START_REF] Bouchut | A subsonic-well-balanced reconstruction scheme for shallow water flows[END_REF] for a Suliciu relaxation scheme dealing with the subcritical case, to [START_REF] Berthon | A fully well-balanced, positive and entropy-satisfying Godunov-type method for the shallow-water equations[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography or Manning friction[END_REF][START_REF] Michel-Dansac | A two-dimensional high-order well-balanced scheme for the shallow water equations with topography and Manning friction[END_REF] for Godunov-type schemes that exactly capture moving steady states, or to [START_REF] Gómez-Bueno | High-order well-balanced methods for systems of balance laws: a control-based approach[END_REF] for a control-based approach.

In this work, we are interested in a linear extension of the hydrostatic reconstruction [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] in order to exactly capture moving steady states. This issue is addressed by designing new reconstructed states W n i+1/2,± at the interface x i+1/2 . Unlike the generalized hydrostatic reconstruction from [START_REF] Castro | Well-balanced numerical schemes based on a generalized hydrostatic reconstruction technique[END_REF], where the nonlinear Bernoulli equation in solved in each cell and where the positivity of the water height may be lost, we here introduce a suitable linearization. The relevant properties of the hydrostatic reconstruction, namely the ability to deal with transitions between wet and dry areas, as well as its versatility, are also satisfied by the technique derived in this paper.

To address such issues, the paper is structured as follows. Firstly, Section 2 contains a few necessary comments related to the steady states (1.7). Indeed, the smoothness of the steady solutions (1.7) is a priori governed by the smoothness of the topography function Z. However, because of the topography discretization, the required smoothness is lost. As a consequence, a discontinuous extension of (1.7) must be considered. In fact, the discontinuous extension of the steady solution turns out to be the main ingredient of the forthcoming hydrodynamic extension. Then, in Section 3, we introduce an easy interface reconstruction technique that preserves the moving and non-moving steady states, building on the versatility of the hydrostatic reconstruction [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. Moreover, we prove that the derived hydrodynamic reconstruction can handle transitions between wet and dry areas. Since the hydrodynamic reconstruction is defined according to a sequence of properties, in Section 4, we specify the reconstruction and we present an explicit definition. Next, in Section 5, we suggest a linear high-order accurate technique which preserves the moving and non-moving steady states. Finally, in Section 6, we present numerical experiments that illustrate the relevance of the hydrodynamic reconstruction.

Discretization of steady solutions

This section is devoted to some comments about the discretization of the smooth steady solutions given by (1.7). First, let us emphasize that arguing the smoothness of the topography function Z(x), we have

Z i+1 -Z i = O(∆x).
Now, we consider a non-dry region; namely with h n i > 0. For a given pair (q 0 , B 0 ) ∈ R 2 and an arbitrary time t n , according to (1.7) a discrete steady solution is naturally given by, for all i ∈ Z,

     q n i =: q 0 , (q n i ) 2 2(h n i ) 2 + g(h n i + Z i ) =: B 0 .
(2.1)

In fact, we now show that this natural choice of the approximate steady solutions may introduce inconsistencies coming from a loss of smoothness in (1.7). Indeed, from (2.1), we easily obtain a local per interface definition of a steady solution, given as follows for all i ∈ Z:

     q n i = q n i+1 , (q n i ) 2 2(h n i ) 2 + g(h n i + Z i ) = (q n i+1 ) 2 2(h n i+1 ) 2 + g(h n i+1 + Z i+1 ).
(2.2)

Since, at the interface located at x i+1/2 , the pair (Z i , Z i+1 ) defines a discontinuous topography function with a small jump in O(∆x), the local definition (2.2) is not sufficient to ensure the smoothness required to derive (1.7). Actually, even if the error to the smoothness is small and controlled by O(∆x), this failure may create non-physical steady solutions. Indeed, let us momentarily consider the equivalent augmented system [START_REF] Lefloch | The Riemann problem for the shallow water equations with discontinuous topography[END_REF] as follows:

         ∂ t h + ∂ x q = 0, ∂ t q + ∂ x q 2 h + 1 2 gh 2 + gh∂ x Z = 0, ∂ t Z = 0. (2.3)
We easily see that this system contains a stationary wave endowed with the following Riemann invariants:

q and q 2 2h 2 + g(h + Z).

(2.4)

Naturally, from these Riemann invariants, we easily recover the steady solutions (1.7). As a consequence, (2.2) exactly coincides with the preservation of the Riemann invariants for the stationary wave. In other words, (2.2) is relevant in defining discrete steady solutions independently of the smoothness of the pair (Z i , Z i+1 ). However, it is essential to note that (2.2) can produce physically inconsistent solutions. Indeed, as demonstrated in [START_REF] Lefloch | The Riemann problem for the shallow water equations with discontinuous topography[END_REF], the Riemann problem of the shallow water equations (1.1) -or equivalently (2.3) -with discontinuous topography may admit non-unique solutions. While the steady solution preserving the Riemann invariants (2.4) is a possible solution, it is not the only one and does not seem to be the expected physical solution [START_REF] Aleksyuk | The exact Riemann solver for the shallow water equations with a discontinuous bottom[END_REF][START_REF] Lefloch | The Riemann problem for the shallow water equations with discontinuous topography[END_REF].

At this level, we conjecture that the steady solution made of the Riemann invariants (2.4) is unstable, whereas the physical unsteady solution, with the same initial data made of constant Riemann invariants (2.4), is stable in a sense to be prescribed. Such an assertion will be illustrated by numerical experiments displayed in Section 6.

Another important comment about the steady solutions concerns the wet/dry transitions. Indeed, after [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF], non-moving steady solutions and dry regions cannot coexist. As a consequence, as soon as a dry area is present, smooth steady solutions must be at rest.

A family of moving steady states-preserving schemes

In this section, we present a simple state interface reconstruction method that enables the numerical method to preserve both moving and non-moving steady solutions as given by (2.2), following the approach by Audusse et al. [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. With a scheme given by (1.8), in the spirit of the usual hydrostatic reconstruction, we here design the reconstructed states W n i+1/2,± at the interface x i+1/2 and the source term discretization S n i in order to preserve the expected moving steady solutions (2.1).

The interface reconstructed states are now given by W n i+ 1 2 ,± = (h n i+ 1 2 ,± , q n i+ 1 2 ,± ) ⊺ , where we have set

               h n i+ 1 2 ,-= max 0, h n i + Z i -Z i+ 1 2 + 2 Fr 2 (h n i , h n i+ 1 2 , q n i ) H h n i , h n i+ 1 2 , q n i , Z i+ 1 2 -Z i , h n i+ 1 2 ,+ = max 0, h i+1 + Z i+1 -Z i+ 1 2 + 2 Fr 2 (h i+1 , h i+ 1 2 , q i+1 ) H h i+1 , h i+ 1 2 , q i+1 , Z i+ 1 2 -Z i+1 , q n i+ 1 2 ,-= q n i , q i+ 1 2 ,+ = q n i+1 , (3.1)
with an intermediate reconstruction (W i+ 1 2 , Z i+ 1 2 ) at the interface x i+ 1 2 defined by

(W n i+ 1 2 , Z i+ 1 2 ) = (W n i , Z i ) if Z i > Z i+1 , (W n i+1 , Z i+1 ) otherwise. (3.2)
Moreover, we have introduced the approximate Froude number as follows:

Fr 2 (h n i , h n i+ 1 2 , q n i ) = (q n i ) 2 (h n i + h n i+ 1 2 ) 2g(h n i ) 2 (h n i+ 1 2 ) 2 . (3.3)
Next, concerning the source term discretization, S n i = (0, (S q ) n i ) ⊺ , we have adopted the following definition:

∆x(S q ) n i = -g 2h n i-1 2 ,+ h n i+ 1 2 ,- h n i-1 2 ,+ + h n i+ 1 2 ,- Z i+ 1 2 -Z i-1 2 + 4g h n i-1 2 ,+ + h n i+ 1 2 ,- H h n i-1 2 ,+ , h n i+ 1 2 ,-, q n i , Z i+ 1 2 -Z i-1 2 3 . (3.4)
At this level, it is worth noticing that the hydrodynamic reconstruction cannot be fully characterized until the function H is defined. A possible choice of the particular function is detailed in Section 4.

From now on, let us emphasize that the interface reconstruction (3.1) is nothing but the standard hydrostatic reconstruction from [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF], augmented with an additional term governed by the function H. Of course, this new term serves to perturb the hydrostatic reconstruction in order to provide a hydrodynamic reconstruction that is capable of preserving the moving steady state solutions. Now, we impose suitable hypotheses to be satisfied by the perturbation H so that the hydrodynamic reconstruction scheme (1.8) is consistent, preserves the steady solutions (2.1), and efficiently handles wet/dry transitions.

To address this issue, the function H must be endowed with several properties.

(H-1) In order to recover the required consistency, the perturbation H : R ⋆ + ×R ⋆ + ×R×R → R must be continuous. Moreover, it must satisfy the following asymptotic behavior for all h L > 0, h R > 0 and q ∈ R:

lim ∆Z→0 H(h L , h R , q, ∆Z) = 0.
(H-2) Next, the well-balanced property will be recovered by imposing

H (h L , h R , q, ∆Z) = 1 2 (h R -h L ), for all (h L , h R , q, ∆Z) ∈ R ⋆ + × R ⋆ + × R × R such that ∆Z = -(h R -h L ) 1 -Fr 2 (h L , h R , q) ,

and with Fr

2 (h L , h R , q 0 ) ̸ = 1 so that h R -h L = O(∆Z).
(H-3) The last property we enforce concerns the wet/dry transition. To that end, we have to make sure that all the involved quantities are well-defined in dry regions. To address such an issue, we require that the following limit holds in order for h n i-1 2 ,+ and h n i+ 1 2 ,-to be bounded:

lim h n i →0 Fr 2 (h n i , h n i± 1 2 , q n i ) H h n i , h n i± 1 2 , q n i , Z i± 1 2 -Z i = 0.
Next, to make sure that the source term vanishes in dry areas, we also have to impose

lim h n i →0 1 h n i-1 2 ,+ + h n i+ 1 2 ,- H h n i-1 2 ,+ , h n i+ 1 2 ,-, q n i , Z i+ 1 2 -Z i-1 2 3 = 0, lim h n i →0 h n i-1 2 ,+ h n i+ 1 2 ,- h n i-1 2 ,+ + h n i+ 1 2 ,- = 0.
Before we state the main properties satisfied by the scheme (1.8)-(3.1)-(3.4), let us underline once again that H is nothing but a small perturbation. This is due to the fact that the topography function is assumed to be smooth, resulting in Z i+1 -Z i being of order ∆x. As a consequence of (H-1), we have

H(h n i , h n i+1/2 , q n i , Z i+1/2 -Z i ) = O(∆x) and H(h n i+1 , h n i+1/2 , q n i+1 , Z i+1/2 -Z i+1 ) = O(∆x) for all i ∈ Z.
In this sense, the application H is clearly a small perturbation of the original hydrostatic reconstruction designed in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. Moreover, it is also worth mentioning that obtaining an explicit definition of H, such that assumptions (H-1), (H-2) and (H-3) hold, is not a trivial task. Section 4 is devoted to exhibiting an admissible perturbation. More specifically, we will show that the set of admissible perturbation according to assumptions (H-1), (H-2) and (H-3) is not empty, and we will design suitable approximations of the admissible functions H.

In addition, let us emphasize that assumptions (H-3) are necessary to define a wet/dry transition. In dry regions, where h n i-1 = h n i = h n i+1 = 0, the hydrodynamic reconstruction (3.1) and the source term (3.4) are well-defined, and they vanish due to conventions (1.3).

We are now able to state our main result, which outlines the properties of the scheme (1.8) endowed with the hydrodynamic reconstruction (3.1) and the source term discretization (3.4).

Theorem 1. Let H(h L , h R , q 0 , ∆Z) be a function which satisfies the assumptions (H-1), (H-2) and (H-3). For non-negative water heights h n i ≥ 0 for all i ∈ Z, the scheme (1.8)-(3.1)-(3.4) satisfies the following properties:

(1-a) it is consistent with the shallow water equations (1.1);

(1-b) it preserves the steady states with nonzero velocity, in the sense that if

(W n i ) i∈Z satisfy (2.2) for all i ∈ Z and Fr 2 (h n i , h n i± 1 2 , q n i ) ̸ = 1 then W n+1 i = W n i ; (1-c) it is non-negativity-preserving, i.e., if h n i ≥ 0, then h n+1 i ≥ 0.
Now, in order to establish this main result, we first need to prove the following two lemmas. They state some needed properties to be satisfied by the hydrodynamic reconstruction. The first lemma deals with wet areas, while the second one addresses the wet/dry and dry cases.

Lemma 2. For positive water heights h n i > 0 for all i ∈ Z far away from dry areas, the hydrodynamic reconstruction (3.1), with assumptions (H-1) and (H-2), satisfies the following properties:

(2-a) if Z i = Z i+1 , then h n i+ 1 2 ,-= h n i and h n i+ 1 2 ,+ = h n i+1 ;
(2-b) the hydrodynamic reconstruction (3.1) degenerates to the standard hydrostatic reconstruction from [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] as soon as q n i = q n i+1 = 0;

(2-c) if two consecutive states (W n i , Z i ) and (W n i+1 , Z i+1 ) satisfy the local per interface steady state definition (2.2), then h n i+ 1 2 ,-= h n i+ 1 2 and h n i+ 1 2 ,+ = h n i+ 1 2 .
Proof. The proof of (2-a) is immediate, since it relies on property (H-1).

Similarly, by inspection, we note that if

q n i = q n i+1 = 0 in (3.1), then Fr 2 (h i , h i± 1 2 , 0) = 0, to get      h n i+ 1 2 ,-= max 0, h n i + Z i -Z i+ 1 2 , h n i+ 1 2 ,+ = max 0, h n i+1 + Z i+1 -Z i+ 1 2 ,
which is nothing but the standard hydrostatic reconstruction, and thus (2-b) holds. Regarding (2-c), since (W i , Z i ) and (W i+1 , Z i+1 ) define a local per interface steady state according to (2.2), then we have

     q n i = q n i+1 =: q 0 , q 2 0 2(h n i ) 2 + g(h n i + Z i ) = q 2 0 2(h n i+1 ) 2 + g(h n i+1 + Z i+1 ) =: B 0 .
Next, by definition of (h i+ 1 2 , Z i+ 1 2 ) given by (3.2), we get

q 2 0 2(h n i+ 1 2 ) 2 + g(h n i+ 1 2 + Z i+ 1 2 ) = B 0 .
As a consequence, we have

q 2 0 2(h n i ) 2 +g(h n i +Z i ) = q 2 0 2(h n i+ 1 2 ) 2 +g(h n i+ 1 2 +Z i+ 1 2 ) and q 2 0 2(h n i+1 ) 2 +g(h n i+1 +Z i+1 ) = q 2 0 2(h n i+ 1 2 ) 2 +g(h n i+ 1 2 +Z i+ 1 2 ), to get Z i+ 1 2 -Z i = -(h n i+ 1 2 -h n i ) 1-Fr 2 (h n i , h n i+ 1 2 , q n i ) and Z i+ 1 2 -Z i+1 = -(h n i+ 1 2 -h n i+1 ) 1-Fr 2 (h n i+1 , h n i+ 1 2 , q n i+1 ) .
Next, arguing property (H-2), we obtain

H h n i , h n i+ 1 2 , q n i , Z i+ 1 2 -Z i = 1 2 h n i+ 1 2 -h n i and H h n i+1 , h n i+ 1 2 , q n i+1 , Z i+ 1 2 -Z i+1 = 1 2 h n i+ 1 2 -h n i+1 .
Plugging these values into the definition (3.1) yields the following chain of equalities:

h n i+ 1 2 ,-= max 0, h n i + Z i -Z i+ 1 2 + (q n i ) 2 (h n i + h n i+ 1 2 ) g(h n i ) 2 (h n i+ 1 2 ) 2 1 2 h n i+ 1 2 -h n i , = max 0, h n i + Z i -Z i+ 1 2 + q 2 0 2g 1 (h n i ) 2 - 1 (h n i+ 1 2 ) 2 , = max 0, B 0 g -Z i+ 1 2 - q 2 0 2g 1 (h n i+ 1 2 ) 2 , = max 0, h n i+ 1 2 , = h n i+ 1 2 .
Similar relations lead to h i+ 1 2 ,+ = h i+ 1 2 , which concludes the proof of Lemma 2. Next, in addition to the above result, it is also necessary to establish the behavior of the hydrodynamic reconstruction (3.1) in dry areas. This is the object of the following lemma. Lemma 3. Assume H satisfies the assumptions (H-1), (H-2) and (H-3). Then, the hydrodynamic reconstruction (3.1) verifies:

(3-a) if h n i = 0, then h n i-1 2 ,+ = h n i+ 1 2 ,-= 0; (3-b) if h n i = 0, then S n i = 0.
Proof. Let us first focus on property (3-a). Because of assumptions (H-3), arguing the water height reconstruction (3.1) in the cell (x i-1/2 , x i+1/2 ) with h n i = 0, we immediately get

h n i-1 2 ,+ = max(0, Z i-1 -Z i-1 2 ) and h n i+ 1 2 ,-= max(0, Z i -Z i+ 1 2 ). Since we have Z i+ 1 2 = max(Z i , Z i+1 ) for all i ∈ Z, we then obtain h n i-1 2 ,+ = 0 and h n i+ 1
2 ,-= 0. Next, concerning (3-b), once again assumptions (H-3) immediately enforce S n i = 0 as soon as h n i = 0, and the proof is achieved.

From now on, it should be noted that properties (3-a) and (3-b), which comprise Lemma 3, are satisfied by the standard hydrostatic reconstruction from [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF].

Using the intermediate results established in Lemma 2 and Lemma 3, we can now proceed to prove Theorem 1.

Proof of Theorem 1. Regarding the consistency in (1-a), let us recall that the numerical flux F is assumed to be consistent with the exact flux function F (W ) given by (1.4). Since this consistent flux function F is applied without modification to the reconstructed states given by (3.1), the consistency of the flux is maintained, based on the proof given in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. We still need to prove the consistency of the source term (3.4). More specifically, we have to show that (S q ) n i is consistent with -gh∂ x Z. The approximate source term reads

(S q ) n i = -g 2h n i-1 2 ,+ h n i+ 1 2 ,- h n i-1 2 ,+ + h n i+ 1 2 ,- Z i+ 1 2 -Z i-1 2 ∆x + 4g h n i-1 2 ,+ + h n i+ 1 2 ,- H h n i-1 2 ,+ , h n i+ 1 2 ,-, q n i , Z i+ 1 2 -Z i-1 2 3
∆x .

(3.5)

Arguing (H-1), we obtain that the right part of (S q ) n i vanishes when ∆x approaches 0 due to the smoothness assumption on the topography. Then, it is obvious that the remainder is consistent with -gh∂ x Z, which concludes the proof of (1-a).

Next, let us establish the well-balanced property (1-b). Assume that (W n i-1 , Z i-1 ), (W n i , Z i ) and (W n i+1 , Z i+1 ) define the same steady state (2.1), with constant discharge q 0 and Bernoulli's constant B 0 . We now have to prove that

(W n+1 i , Z i ) = (W n i , Z i ).
According to Lemma 2 (property (2-c)), the hydrodynamic reconstruction (3.1) becomes

h n i± 1 2 ,-= h n i± 1 2 ,+ = h n i± 1 2 , q n i± 1 2 ,-= q n i± 1 2 ,+ = q n i± 1 2 = q 0 . Let us set W n i± 1 2 = (h n i± 1 2
, q 0 ) ⊺ so that the scheme (1.8) now reads

W n+1 i = W n i - ∆t ∆x F W n i+ 1 2 , W n i+ 1 2 -F W n i-1 2 , W n i-1 2 + ∆t S n i . (3.6)
Since F is a consistent numerical flux, we know that F(W, W ) = F (W ), where F is the flux of the shallow water equations defined by (1.4). Thus, (3.6) writes

h n+1 i q n+1 i = h n i q n i - ∆t ∆x    q n i+ 1 2 -q n i-1 2 (q n i+ 1 2 ) 2 h n i+ 1 2 + 1 2 g(h n i+ 1 2 ) 2 - (q n i-1 2 ) 2 h n i-1 2 - 1 2 g(h n i-1 2 ) 2    + ∆t 0 (S q ) n i .
(3.7)

We immediately note that W n+1 i = W n i = (h n i , q 0 ) ⊺ for all i ∈ Z as soon as the source term (S q ) n i satisfies

∆x(S q ) n i = q 2 0 h n i+ 1 2 + 1 2 g(h n i+ 1 2 ) 2 - q 2 0 h n i-1 2 - 1 2 g(h n i-1 2 ) 2 . (3.8)
Now, to establish the above relation, since we have

h n i+ 1 2 ,± = h n i+ 1 2 and h n i-1 2 ,± = h n i-1 2
, the source term (S q ) n i defined by (3.5) reads

(S q ) n i = -g 2h n i-1 2 h n i+ 1 2 h n i-1 2 + h n i+ 1 2 Z i+ 1 2 -Z i-1 2 ∆x + 4g h n i-1 2 + h n i+ 1 2 H h n i-1 2 , h n i+ 1 2 , q 0 , Z i+ 1 2 -Z i-1 2 3
∆x .

(3.9)

Next, we remark that the states

(W n i-1 2 , Z i-1 2 ) and (W n i+ 1 2 , Z i+ 1 2
) satisfy the local per interface steady state condition (2.2). We then have

Z i+ 1 2 -Z i-1 2 = -(h n i+ 1 2 -h n i-1 2 ) 1 -Fr 2 (h n i-1 2 , h n i+ 1 2 , q 0 ) , (3.10)
and, as a consequence of hypothesis (H-2), we get

H(h n i-1 2 , h n i+ 1 2 , q 0 , Z i+ 1 2 -Z i-1 2 ) = 1 2 (h n i+ 1 2 -h n i-1 2
).

Therefore, (3.9) now reads ∆x(S q ) n i = -g

2h n i-1 2 h n i+ 1 2 h n i-1 2 + h n i+ 1 2 Z i+ 1 2 -Z i-1 2 + g h n i+ 1 2 -h n i-1 2 3 2 h n i-1 2 + h n i+ 1 2 . (3.11)
Plugging (3.10) into (3.11), the above relation reformulates

∆x(S q ) n i = q 2 0 h n i+ 1 2 + 1 2 g(h n i+ 1 2 ) 2 - q 2 0 h n i-1 2 - 1 2 g(h n i-1 2 ) 2 , (3.12)
which concludes the proof of the well-balanced property (1-b).

The final property (1-c) is a direct consequence of the definition of h n i+1/2,± , which involves taking a maximum with 0. Then, to prove the non-negativity preservation satisfied by the scheme, we exactly follow the proof given in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. The establishment of Theorem 1 is thus completed.

One possible choice for the function H

The goal of this section is to propose an expression for the function H that satisfies the required properties (H-1) through (H-3). Recall that H is a function from

R ⋆ + × R ⋆ + × R × R to R, applied for instance on (h n i , h n i+1/2 , q n i , Z i+1/2 -Z i ) to compute h i±1/2,-.
For the sake of clarity, throughout this section, H will be written as a function H(h L , h R , q, ∆Z). Its arguments may be omitted for more concise notation.

As a first step, consider H given as a solution to the following polynomial equation of degree five, in the spirit of [START_REF] Berthon | A fully well-balanced, positive and entropy-satisfying Godunov-type method for the shallow-water equations[END_REF]:

2 H g h2 -H 2 2 -q2 h = -g∆Z h2 -H 2 2 , (4.1)
where h = (h L + h R )/2. This expression allows us to state the following result. Proof. For this proof, we consider data far from a dry area. Let us start with property (H-2). To that end, we assume a steady solution with Fr 2 (h L , h R , q) ̸ = 1, and we take ∆Z = -(h R -h L ) 1 -Fr 2 (h L , h R , q) . Given the expression (3.3) of the squared Froude number, plugging the above value of ∆Z in (4.1) leads to:

2 H g h2 -H 2 2 -q2 h -(h R -h L ) g h2 -H 2 2 -q2 h h2 -H 2 2 h 2 L h 2 R = 0. (4.2)
Now, we prove that H = (h R -h L )/2 is a solution to the fifth-degree equation (4.2). First, let us note that

1 h L h R h2 - (h R -h L ) 2 4 = 1 h L h R (h R + h L ) 2 -(h R -h L ) 2 4 = 1.
Then, plugging H = (h R -h L )/2 in the left-hand side of (4.2) yields:

2 (h R -h L ) 2 g h2 - (h R -h L ) 2 4 2 -q2 h -(h R -h L ) g h2 - (h R -h L ) 2 4 2 -q2 h = 0,
which proves that H = (h R -h L )/2 is one of the multiple solutions to (4.2). This proves property (H-2). We turn to property (H-1). The solution to (4.2) we have just exhibited satisfies property (H-2), i.e. it satisfies

H = (h R -h L )/2 = O(h R -h L ) as soon as ∆Z = -(h R -h L ) 1 -Fr 2 (h L , h R , q) for Fr 2 (h L , h R , q) ̸ = 1.
However, note that, when ∆Z = 0, this condition is also satisfied, and we get, in this case, h R -h L = O(∆Z). This immediately proves that the previously exhibited solution H satisfies H = ∆Z→0 O(∆Z), thereby proving property (H-1).

The proof is thus concluded, since we have exhibited a solution to (4.2) that satisfies both properties (H-1) and (H-2). Therefore, according to Lemma 4, the nonlinear equation (4.1) has a solution that satisfies both properties (H-1) and (H-2). However, finding the correct solution is a complex process that would negate all the benefits of our linear approach. Indeed, we would need to use Newton's method at each interface and each time step to compute the solutions to (4.1). Moreover, we would then have to choose the correct solution from among (at most) five possible ones. Thus, we elect not to pursue this nonlinear direction. Instead, we provide a relevant linearization of part of (4.1).

Assuming that H ̸ = ± h, (4.1) rewrites

2 H 1 - q2h h2 -H 2 2 = -∆Z. (4.3)
We temporarily assume that h L ̸ = h R , and we set ∆h = h R -h L . We suggest the following linearization around H = ∆h/2 of the expression in brackets in (4.3), thus modifying the equation satisfied by H, to get a quadratic equation in H:

2 H 1 - q2 (h L + h R ) 2gh 2 L h 2 R + 4 sgn(∆Z) |∆Z| |∆h| 3 (∆h -2 H) = -∆Z,
This expression can be simplified by remarking that

1 -Fr 2 (h L , h R , q) = 1 - q2 (h L + h R ) 2gh 2 L h 2 R , to get 2 H 1 -Fr 2 +4 sgn(∆Z) |∆Z| |∆h| 3 (∆h -2 H) = -∆Z.
Solving this quadratic equation for H leads to

H = 1 4     ∆h + 1 -Fr 2 4 sgn(∆Z) |∆h| 3 |∆Z| ± ∆h + 1 -Fr 2 4 sgn(∆Z) |∆h| 3 |∆Z| 2 + |∆Z||∆h| 3     . (4.4)
Choosing the correct sign for the ± in (4.4) makes it possible to state, and prove, the following result.

Lemma 5. Let

H = 1 4 E -sgn(1 -Fr 2 ) sgn(∆Z) E 2 + |∆Z||∆h| 3 , with E = ∆h + 1 -Fr 2 4 sgn(∆Z) |∆h| 3 |∆Z| . (4.5)
Then H satisfies the required properties (H-1) (if ∆Z and 1-Fr 2 do not simultaneously vanish), (H-2) and (H-3).

Proof. We prove the three properties (H-1), (H-2) and (H-3) in order.

(H-1) Let us first note that (4.5) contains a division by ∆Z. Nevertheless, this expression turns out to be infinitely continuously differentiable around ∆Z = 0.

To prove (H-1), let us compute the limit of H when ∆Z goes to 0. To that end, for the sake of simplicity, we consider the case where ∆Z > 0, ∆h > 0 and 1 -Fr 2 > 0. In this case, the Taylor expansion provided in Appendix A proves (H-1). An immediate consequence is that H tends to 0 as ∆Z tends to 0, despite the a priori indeterminate division by |∆Z|. Investigating the other cases (1 -Fr 2 < 0 or ∆h ≤ 0) yields the same limit.

(H-2) We now prove that the expression of H satisfies property (H-2), i.e. that H = ∆h/2 when a steady solution has been reached. We therefore assume that the solution is steady, i.e. that ∆Z = -(1 -Fr 2 )∆h, so that the expression of E is simplified as follows:

E = ∆h + 1 -Fr 2 4 sgn(-(1 -Fr 2 )∆h) |∆h| 3 |(1 -Fr 2 )||∆h| = ∆h 1 - 1 4 |1 -Fr 2 | ,
and E 2 + |∆Z||∆h| 3 becomes:

E 2 + |∆Z||∆h| 3 = |∆h| 2 1 + 1 4 |1 -Fr 2 | 2 .
Plugging the two expressions above into H, we get the following chain of equalities:

H = 1 4   ∆h 1 - 1 4 |1 -Fr 2 | + sgn(∆h) |∆h| 2 1 + 1 4 |1 -Fr 2 | 2   = ∆h 4 1 - 1 4 |1 -Fr 2 | + 1 + 1 4 |1 -Fr 2 | = ∆h 2 ,
which completes the proof of (H-2).

(H-3) Finally, we turn to (H-3). We first make the general remark that H, given by (4.5), rewrites as

H(h L , h R , q, ∆Z) = (h R -h L )B(h L , h R , q, ∆Z), (4.6)
with B a bounded function. The boundedness of B comes from the boundedness assumption on the Froude number, given in (1.3).

To prove the first equality of (H-3), we consider the function

C(h L , h R , q L ) = 2 Fr 2 (h L , h R , q L ) H(h L , h R , q L , ∆Z). (4.7)
We have to prove that C goes to zero when h L tends to 0, in the two cases where h R = h L and h R = 0. Technical Taylor expansions of C are provided in Appendix B. They show that this property is satisfied, thereby proving the first equality of (H-3).

To prove the second equality of (H-3), we have to show that lim

h n i →0 + H h i-1 2 ,+ , h i+ 1 2 ,-, q i , Z i+ 1 2 -Z i-1 2 3 h i-1 2 ,+ + h i+ 1 2 ,- = 0. (4.8)
To that end, recall from property (3-a) that, if h i = 0, then h i-1 2 ,+ = h i+ 1 2 ,-= 0. This follows from the fact that H satisfies (H-3). Therefore, proving (4.8) requires proving that

lim h L →0 + h R →0 + H(h L , h R , q, ∆Z) 3 h L + h R = 0. (4.9)
Arguing (4.6), we get

H 3 h L + h R ≤ |h R -h L | 3 |h L | + |h R | |B| 3 ≤ (|h R | + |h L |) 3 |h L | + |h R | max |B| 3 = (|h R | + |h L |) 2 max |B| 3 .
Taking the limit of the above expression when h L and h R go to zero proves that (4.9) holds. Therefore, the third equality of (H-3) is satisfied.

Finally, the proof of the third equality of (H-3) is obtained by arguing the same arguments as above, in addition to the well-known limit lim

x→0 + y→0 + xy x + y = 0.
The proof of (H-3) is therefore concluded.

We have therefore proven the three properties, which concludes the proof.

As a conclusion, the expression (4.5) of H satisfies the required properties, as proven in Lemma 5. The only limit where this expression is ill-defined is when both Fr 2 tends to 1 and ∆Z tends to 0. This corresponds to a well-known issue, since the steady equations are themselves ill-defined when Fr 2 = 1 and ∆Z = 0. This resonant case is well-documented in the literature, as seen for instance in [START_REF] Castro-Orgaz | Minimum Specific Energy and Transcritical Flow in Unsteady Open-Channel Flow[END_REF][START_REF] Gómez-Bueno | Collocation Methods for High-Order Well-Balanced Methods for Systems of Balance Laws[END_REF] and references therein. In the numerical experiments, we set H = 0 when Fr 2 = 1 and ∆Z = 0.

Linear extension to high-order accuracy

In Section 3, we proposed a linear reconstruction technique able to turn any first-order, non-well-balanced scheme into a fully well-balanced one. The present section is dedicated to its high-order extension. Since the first-order scheme is linear, we aim to develop a linear high-order and well-balanced extension, i.e. one that does not require solving Bernoulli's relations. Indeed, typically, fully well-balanced high-order extensions require solving some nonlinear equations, which can be computationally expensive and may fail in the case of multiple or non-existent solutions. For instance, we mention [START_REF] Castro | Well-balanced high order extensions of Godunov's method for semilinear balance laws[END_REF][START_REF] Xing | Exactly well-balanced discontinuous Galerkin methods for the shallow water equations with moving water equilibrium[END_REF][START_REF] Britton | High Order Still-Water and Moving-Water Equilibria Preserving Discontinuous Galerkin Methods for the Ripa Model[END_REF][START_REF] Gómez-Bueno | Collocation Methods for High-Order Well-Balanced Methods for Systems of Balance Laws[END_REF], where nonlinear equations have to be numerically solved.

To avoid the computational expense associated with a nonlinear solver, we adopt the method described in [START_REF] Berthon | A Very Easy High-Order Well-Balanced Reconstruction for Hyperbolic Systems with Source Terms[END_REF], which was initially introduced on a generic hyperbolic system, and later applied to the 2D shallow water case in [START_REF] Michel-Dansac | A two-dimensional high-order well-balanced scheme for the shallow water equations with topography and Manning friction[END_REF]. For the sake of completeness, the application of this technique to the current scheme is summarized below. For the remainder of this section, we build a scheme of order d + 1.

To that end, we start by considering the following generic polynomial reconstruction of degree d in space (see for instance [START_REF] Diot | Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF][START_REF] Diot | The multidimensional optimal order detection method in the threedimensional case: very high-order finite volume method for hyperbolic systems[END_REF]):

W n i (x) = W n i + Π n i (x -x i ), (5.1)
where the degree d polynomial Π n i is defined such that, for all x ∈ (x i-1 2 , x i+ 1 2 ), the following relations hold:

W n i (x) = W (x, t n ) + O(∆x d+1 ) and 1 ∆x x i+ 1 2 x i-1 2 W n i (x) dx = W n i . (5.2)
This reconstruction of degree d naturally defines a scheme of space order d + 1, see for instance [START_REF] Diot | Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF][START_REF] Diot | The multidimensional optimal order detection method in the threedimensional case: very high-order finite volume method for hyperbolic systems[END_REF].

Given the polynomial reconstruction (5.1), we can define a high-order, non-well-balanced scheme. To that end, consider the following modification of (1.8):

W n+1 i = W n i - ∆t ∆x F W n i+ 1 2 ,-, W n i+ 1 2 ,+ -F W n i-1 2 ,-, W n i-1 2 ,+ + ∆t S n i , (5.3)
where the polynomial reconstruction is applied at each inner interface, to get

W n i+ 1 2 ,-= W n i + Π n i ∆x 2 and W n i+ 1 2 ,+ = W n i+1 + Π n i+1 - ∆x 2 ,
and where the source term is nothing but a high-order approximation of the averaged source term:

S n i = 1 ∆x x i+ 1 2 x i-1 2 S(W (x, t n ), x) dx + O(∆x d+1 ).
To compute S n i in practice, one can use a quadrature formula of order d + 1, see for instance [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]. Equipped with the high-order scheme (5.3), we make it well-balanced, without having to solve nonlinear equations. We accomplish this by leveraging the procedure from [START_REF] Berthon | A Very Easy High-Order Well-Balanced Reconstruction for Hyperbolic Systems with Source Terms[END_REF]. To address this issue, we first introduce a steady solution detector.

A steady state detector

Let us define the following indicator, which will help us write the well-balanced extension of the high-order scheme:

θ n i+ 1 2 = ε n i+ 1 2 ε n i+ 1 2 + ∆x C n i+ 1 2 d+1 , with (5.4)
ε n i+ 1 2 =   q n i+1 (q n i+1 ) 2 2(h n i+1 ) 2 + g(h n i+1 + Z i+1 )   -   q n i (q n i ) 2 2(h n i ) 2 + g(h n i + Z i )   , (5.5)
where

C n i+ 1 2 ̸ = 0 is independent of ∆x. An expression of C n i+ 1 2
will be proposed before performing the numerical experiments, in (6.1). The properties enjoyed by this indicator are summarized in the following result. Proposition 6. The expressions (5.4) -(5.5) ensure the following properties:

(6-a) if (W n i , Z n i ) and (W n i+1 , Z n i+1
) define a steady state, i.e. satisfy the local per interface steady state relation (2.2), then θ n

i+ 1 2 = 0; (6-b) otherwise, θ n i+ 1 2 = 1 + O(∆x d+1 ).
Proof. The proof is straightforward and present in [START_REF] Berthon | A Very Easy High-Order Well-Balanced Reconstruction for Hyperbolic Systems with Source Terms[END_REF]. The first property is verified by inspection of (5.5), and the second one through a Taylor expansion of θ n i+ 1 2 .

The high-order well-balanced reconstruction

Equipped with the steady solution detector (5.4) -(5.5) and the high-order, non-well balanced scheme (5.3), we are able to fully state the linear, high-order, well-balanced scheme.

To that end, we modify the high-order scheme (5.3) as follows:

W n+1 i = W n i - ∆t ∆x F W n i+ 1 2 ,-, W n i+ 1 2 ,+ -F W n i-1 2 ,-, W n i-1 2 ,+ + ∆t S n i , (5.6)
where the polynomial reconstruction W has been replaced with the modified polynomial reconstruction W , defined by

W n i+ 1 2 ,-= W n i + θ n i+ 1 2 Π n i ∆x 2 and W n i+ 1 2 ,+ = W n i+1 + θ n i+ 1 2 Π n i+1 - ∆x 2 ,
and where the high-order source term S n i has been replaced with a convex combination S n i between the high-order source term S n i and the order one source term S n i , defined in (3.4):

S n i = 1 - θ n i-1 2 + θ n i+ 1 2 2 S n i + θ n i-1 2 + θ n i+ 1 2 2 S n i .
Note that Proposition 6 implies that, if there is a steady state at the interface

x i+ 1 2 , then W n i+ 1 2 ,-= W n i and W n i+ 1 2 ,+ = W n i+1 . Otherwise, W n i+ 1 2 ,-and W n i+ 1
2 ,+ are high-order approximations of the solution at the interface x i+ 1 2 .

Equipped with this modified high-order reconstruction, the high-order hydrodynamic reconstruction is simply computed by applying (3.1) to W n i+ 1 2 ,-and W n i+ 1 2 ,+ at the interface x i+ 1 2 , instead of applying it to W i and W i+1 . Thanks to these definitions, we are able to state the following result, which represents a high-order extension of Theorem 1.

Theorem 7. The scheme (5.3) enjoys the following properties:

(7-a) it is high-order accurate, i.e. consistent with the shallow water equations (1.1) up to O(∆x d+1 ), (7-b) it is non-negativity-preserving, (7-c) it is fully well-balanced, in the sense that it exactly preserves the steady states (2.1).

Proof. The proof of Theorem 7 is present in [START_REF] Berthon | A Very Easy High-Order Well-Balanced Reconstruction for Hyperbolic Systems with Source Terms[END_REF]. For the sake of conciseness, it is omitted here.

Numerical experiments

This last section comprises several numerical experiments, designed to validate the properties of the scheme. Firstly, in Section 6.1, the consistency and order of accuracy are assessed. Secondly, in Section 6.2, we perform several experiments to test the well-balanced property of the scheme (for a lake at rest in Section 6.2.1 and for moving steady states in Section 6.2.2). Thirdly, two dam-break problems are computed in Section 6.3. Lastly, we simulate an unstable steady contact wave in Section 6.4.

Recall that any consistent and non-negativity-preserving numerical flux F can be used in the scheme (1.8). For the numerical experiments, we use the simple HLL flux from [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]. This flux imposes a standard CFL restriction on the time step, as discussed for instance in [START_REF] Vila | Simplified Godunov Schemes for 2 × 2 Systems of Conservation Laws[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics. A practical introduction[END_REF]. When it is required, in Section 6.1 and in Section 6.3, the reference solution is given by the uncorrected HLL scheme with a naive centred discretization of the source term.

For clarity, we label the schemes under consideration as follows:

• the HSR scheme is the first-order accurate hydrostatic reconstruction from [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF],

• the HDRp scheme is the hydrodynamic reconstruction (3.1) for the well-balanced scheme of order p, constructed with the method described in Section 5 if p ≥ 2. For the second-order scheme, we use the minmod limiter (see for instance [START_REF] Van Leer | Towards the Ultimate Conservative Difference Scheme, V. A Second Order Sequel to Godunov's Method[END_REF]) in conjunction with the SSPRK2 time discretization. For the third-order scheme, we use the reconstruction from [START_REF] Schmidtmann | Relations Between WENO3 and Third-Order Limiting in Finite Volume Methods[END_REF] with the SSPRK3 time discretization. Both time discretizations are presented in, for instance, [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF][START_REF] Gottlieb | Strong stability-preserving high-order time discretization methods[END_REF]. Unless otherwise mentioned, the space domain is (0, 1) and the gravity constant is g = 9.81. As prescribed in [START_REF] Berthon | A Very Easy High-Order Well-Balanced Reconstruction for Hyperbolic Systems with Source Terms[END_REF], the coefficient

C n i+ 1 2 in (5.4) is given by C 0 i+ 1 2 = 1 and, for n ≥ 1, C n i+ 1 2 = C θ 1 2 w n i+1 -w n-1 i+1 ∆t + w n i -w n-1 i ∆t , (6.1)
where C θ is a constant depending on the numerical experiment. Unless mentioned in a specific experiment, we take C θ = 1. Before proceeding with the numerical experiments themselves, we briefly mention how dry zones are numerically treated. Let ε m = 2 -52 be the machine epsilon. According to (1.3), the velocity u is computed as follows:

u =    q h if h > ε m , 0 otherwise.
In addition, for the high-order scheme, the non-negativity-preserving limiting procedure from [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF] is used on the reconstructed water height. The remaining potential divisions by h or ∆Z are handled by leveraging properties (H-1), (H-2) and (H-3). Further, to tackle steady solutions with an emerged bottom, it should be noted that they are not solution to Bernoulli's relation (1.7). Indeed, after [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF], such steady states must be at rest, which means that they must satisfy q 0 = 0 in (1.7). In addition, the height and topography satisfy either h i + Z i < Z i+1 and h i+1 = 0 (Figure 1, left panel), or Z i > h i+1 + Z i+1 and h i = 0 (Figure 1, right panel). Since such steady solutions are not solution to Bernoulli's relation, H cannot capture them without a modification. To address this issue, we impose the following additional properties on H:

     h L < ∆Z =⇒ H(h L , 0, 0, ∆Z) = ∆h 2 , h R < -∆Z =⇒ H(0, h R , 0, ∆Z) = ∆h 2 .
These properties ensure that steady states at rest with an emerged bottom are correctly captured by the scheme. Note that a similar technique was already used in [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography or Manning friction[END_REF][START_REF] Michel-Dansac | A two-dimensional high-order well-balanced scheme for the shallow water equations with topography and Manning friction[END_REF] to make smooth-steady-state-capturing schemes able to handle such non-smooth steady solutions at rest.

h i Z i Z i+1 h i+1 Z i+1 Z i
Figure 1: Non-smooth, emerged steady state at rest not governed by (1.7). Left panel: lake at rest with h i+1 = 0 and h i + Z i < Z i+1 . Right panel: lake at rest with h i = 0 and Z i > h i+1 + Z i+1 .

In Figure 2, we display the reference solution and the approximations given by the HSR1, HDR1 and HDR3 schemes with 40 cells. We observe good agreement with the exact solution. To obtain a more precise assessment of the error, error lines are shown in Figure 3. We observe that the schemes exhibit the expected orders of accuracy. As expected, neither the high-order well-balanced procedure nor the hydrodynamic reconstruction impedes the order of accuracy. For the sake of completeness, the values of the errors are reported in Table 1 (we only present errors on h, but the results for q are similar). Figure 3: Experiment from Section 6.1: error in L 2 norm on h (left panel) and on q (right panel) with respect to the number of cells.

Well-balanced property

We now turn to experiments that assess the well-balanced property. We first examine submerged and emerged lake at rest steady solutions in Section 6.2.1. Next, in Section 6.2.2, we tackle the experiments from [START_REF] Goutal | Proceedings of the 2 nd Workshop on Dam-Break Wave Simulation[END_REF], which involve moving steady solutions that are reached after a transient, unsteady state.

Lake at rest

We begin by studying steady states at rest, taking Z(x) = ω(x) once again. The initial discharge is zero everywhere (q(x, 0) = 0), and the initial height is given in Table 2. Note that the resulting initial condition is nothing but a steady state at rest of the shallow water system (1.1). Therefore, since the HSR and HDR schemes are wellbalanced, we expect them to exactly preserve this initial condition. We fix the final time at t end = 1, take 50 cells, and prescribe the exact steady solution as inhomogeneous Dirichlet boundary conditions. We conduct two Table 1: Experiment from Section 6.1: errors (in L 2 norm) and orders of accuracy on h. The errors on q are displayed in Figure 3.

experiments: the first one has a submerged bottom (no dry zones), and the second one has an emerged bottom (with a dry area).

experiment figure h(x, 0)
submerged bottom Figure 4 2 -Z(x) emerged bottom Figure 5 max(0, 0.5 -Z(x))

Table 2: Setup of the experiments from Section 6.2.1.

Submerged bottom. First, in Figure 4, we display the results of the lake at rest with a submerged bottom. As expected, the initial condition is exactly preserved (up to machine precision) by all the schemes under consideration (HSR scheme, HDR scheme, and its high-order extensions). These conclusions are confirmed by the values of the errors reported in Table 3. x q HSR1 HDR1 HDR2 HDR3

Figure 4: Lake at rest with submerged bottom from Section 6.2.1: free surface h + Z (left panel) and discharge q (right panel), at time t end = 1 with 50 cells.

Emerged bottom. Next, the results of the steady state at rest with emerged bottom are depicted in Figure 5, and the errors are collected in Table 4. Similarly to the submerged bottom case, the dry zones did not negatively impact the well-balanced property. error on q 1.24 • 10 -16 5.17 Table 4: Lake at rest with emerged bottom from Section 6.2.1: errors (in L 2 norm) between the initial condition and the approximate solutions at time t end = 1, using 50 cells.

Moving steady solutions

To assess the ability of the HDR scheme to capture moving steady solutions, we now examine the well-known test cases from [START_REF] Goutal | Proceedings of the 2 nd Workshop on Dam-Break Wave Simulation[END_REF]. Namely, we run three test cases: a subcritical flow, a transcritical flow without shock and a transcritical flow with a shock. Each of these test cases follows the same principle: the initial condition consists in a steady state at rest, which is then perturbed by an inflow boundary condition at the left of the domain. After a transient state, the resulting flow becomes a moving steady state (with nonzero velocity). For the subcritical flow and the transcritical flow without shock, this moving steady state satisfies

q = cst = q 0 and B = q 2 2h 2 + g(h + Z) = cst = q 2 0 2H 2 0 + gH 0 .
This final steady state therefore depends on two parameters: the inflow discharge, denoted by q 0 , and the initial free surface, denoted by H 0 . We expect the HDR scheme and its high-order extensions to exactly capture the final steady state, and the HSR scheme to have a nonzero approximation error. However, the transcritical flow with shock, being a non-smooth steady state, is not expected to be exactly captured by the HDR scheme. The space domain in (0, [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF], where the function Z(x) = max(0, 0.05(x -8)(12 -x)) is considered. The initial conditions are defined as h(x, 0)+Z(x) = H 0 and q(x, 0) = q 0 , with H 0 and q 0 provided in Table 5. The final times are also provided in Table 5, and we take 75 discretization cells. At the left boundary, we prescribe homogeneous Neumann boundary conditions on h, and we impose q(0, t) = q 0 . At the right boundary, we prescribe homogeneous Neumann boundary conditions on q, and we impose h (25, Subcritical flow. The first experiment, which converges towards a subcritical flow, is illustrated in Figure 6. Table 6 contains the values of the errors to the underlying steady state. As expected, we observe that the HDR1 scheme, contrary to the HSR1 scheme, exactly captures the resulting moving steady state. In addition, the high-order extensions HDR2 and HDR3 also exactly capture the moving steady state. Table 6: Subcritical flow from Section 6.2.2: errors (in L 2 norm) between the exact steady state and the approximate solutions at time t end = 500, using 75 cells.

Transcritical flow. The results of the second experiment, which involves a transcritical steady flow, are displayed in Figure 7 and Table 7. Similar to the previous case, we observe that the steady state is exactly captured by the HDR scheme, unlike the HDR scheme. However, we observe a small kink near x = 10. This defect arises because, at the critical point x = 10, the Froude number is 1 and the topography derivative vanishes. Note that similar defects were already observed in earlier work, see for instance [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF][START_REF] Berthon | A well-defined moving steady states capturing Godunov-type scheme for Shallow-water model[END_REF]. This kink disappears when the mesh is sufficiently refined. Transcritical flow with shock. The results of the final experiment are displayed in Figure 8. This experiment involves a transcritical flow with a shock; as expected, since it is not smooth, it is not exactly captured by the HDR scheme, let alone by the HSR scheme. However, note that the loss of precision of the HDR scheme only occurs in the vicinity of the shock (around x = 12), while the continuous steady states before and after the shock are exactly preserved.

Dam-break problems

The purpose of these experiments is to evaluate the performance of the scheme on two standard dam-break problems: the first one without dry areas, and the second one with a dry area. For both problems, we take Z(x) = x/2 and set the initial discharge to q(x, 0) = 0. The initial water height is determined according to Table 8, which also contains the values of t end and θ C . Homogeneous Neumann boundary conditions are prescribed, and we take 50 discretization cells.

Wet dam-break. Figure 9 depicts the solutions of the wet dam-break experiment. We observe no difference between the HSR1 and HDR1 schemes, and it is worth noting that the HDR scheme's high-order extensions provide a more accurate approximation of the exact solution, despite minor oscillations on the discharge for the HDR3 scheme. These oscillations are solely due to the high-order polynomial reconstruction. Dry dam-break. The results of the dry dam-break problem are presented in Figure 10. Two significant differences between the HSR1 and HDR1 schemes are noteworthy. First, the HDR1 scheme produces a small kink near the critical point x = 0.5. However, that this kink disappears when the mesh is refined, or when increasing the order of accuracy by using the HDR2 or HDR3 schemes. Second, the HDR scheme produces a more accurate approximation of the wet/dry transition, than the HSR scheme, despite having the same number of cells.

Stationary contact wave

This last experiment corresponds to the situation described in Section 2. The discontinuous topography is given by

Z(x) = Z L = 0 if x < 0.5, Z R = 0.01 otherwise.
We consider a Riemann problem with initial data

W (x, 0) = (h L , q L ) ⊺ if x < 0.5, (h R , q R ) ⊺ otherwise,
with q L = q R = 1, h L = 1 and where h R ≃ 0.2545853624828563 is defined such that, up to machine precision, the Riemann invariants are constant, i.e.

q 2 L 2h 2 L + g(h L + Z L ) = q 2 R 2h 2 R + g(h R + Z R ).
In addition, Neumann boundary conditions are prescribed, and the final time is t end = 0.075. As discussed in Section 2, we conjectured that such Riemann data, although it is solution to the discrete form of Bernoulli's equation, is not a stable steady solution. Therefore, it should not be exactly preserved by the numerical scheme.

In Figure 11, we compare the numerical solution of the four schemes to a reference solution computed with 5000 cells. The numerical solutions with 100 cells are displayed in the top panels, and we observe good agreement between the numerical and reference solutions, especially for the higher order schemes.

However, there is still a kink present around x = 0.5, which corresponds to the of the topography discontinuity. To further analyze this issue, in the bottom panels of Figure 11, we provide a zoom on the interval (0.46, 0.54), computed with 2000 cells. We observe in the bottom left panel that the HSR1 and HDR1 schemes display a sharp water height discontinuity around x = 0.5. Using the higher order HDR2 and HDR3 schemes, the amplitude of this discontinuity decreases. In the bottom right panel, these findings are confirmed, although the discharge remains continuous with the HDR1 scheme, contrary to the HSR1 scheme where a sharp oscillation is present. Like in the case of the water height, the higher order HDR2 and HDR3 schemes allow a better approximation of the reference solution.

Conclusion and outlook

In this paper, we have presented an extension (3.1)-(3.4) to the hydrostatic reconstruction from [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. Applied to a numerical scheme with any consistent numerical flux function, this hydrodynamic reconstruction possesses the following properties:

(i) consistency with the shallow water system (1.1), (ii) preservation of moving steady solutions (1.7) as well as of the lake at rest, (iii) handling of transitions between wet and dry areas. These properties are summarized in Theorem 1. The hydrodynamic reconstruction depends on the choice of a function H, which has to satisfy properties (H-1), (H-2) and (H-3). We have exhibited such a function, and proven that is satisfies the required properties, in Lemma 5. Numerical experiments have confirmed that the numerical scheme endowed with the hydrodynamic reconstruction is indeed consistent, well-balanced, and able to treat dry/wet transitions.

Nevertheless, there are some potential improvements to the method. First, one could design a function H with a more compact expression, without losing the properties outlined in Lemma 5. Second, one could modify the function H to try and remove any kinks appearing when the Froude number approaches unity. To accomplish this, L'Hôpital's rule could be used instead of the standard Bernoulli relations to compute the correct slope of the water height close to the critical point.

A Taylor expansions of H

The goal here is to provide a Taylor expansion of the function H given by (4.5), in the case where ∆Z > 0, ∆h > 0 and 1 -Fr 2 > 0, when ∆Z approaches zero. The computations are performed below, where we have temporarily set F = 1 -Fr 2 in order to save some space.

B Taylor expansions of C

In this section, we give the Taylor expansions of the function C(h L , h R , q L , ∆Z), given by (4.7), when either h L or h R go to 0. The goal is to prove (H-3), i.e., prove that C is continuous when either h L tends to 0 and h R ̸ = 0, and when h L tends to 0 and h R = h L . Recall that C(h L , h R , q L , ∆Z) = 2 Fr 2 (h L , h R , q L , ∆Z) H(h L , h R , q L , ∆Z).

Since q L = h L u L , with u the velocity, we also note that

Fr 2 (h L , h R , h L u L , ∆Z) = h 2 L u 2 L (h L + h R ) 2gh 2 L h 2 R = u 2 L 2gh R 1 + h L h R .
First, we consider the case where h L goes to 0 and h R ̸ = 0. According to assumptions (1.3), in this case, u L also goes to zero. To model this phenomenon, we assume that u L = u(h L ), where the function u is such that u(0) = 0. In this case, we get, again using symbolic computation software,

C = h L =0 + ± 4g 2 h 2 R + 16g h R |∆Z| (±u(0) 2 -2g(h R ± 2∆Z)) + 4gh R (16g∆Z -u(0) 2 ) + u(0) 4 32g 2 h R |∆Z| u(0) 2 ± √ h R u(0) 2 -2gh R √ h R ± 4 |∆Z| 32g 2 h R |∆Z| u(0) 2 + O(h L ).
In the above expression, for the sake of clarity, the ± symbols correspond to sgn(∆Z). In any case, since u(0) = 0, this Taylor expansion shows that lim

h L →0 + h R ̸ =0 C(h L , h R , q L , ∆Z) = 0,
which is what we had set out to prove. Second, we have to prove that lim h L →0 + C(h L , h L , q L , ∆Z) = 0.

However, recall from (4.6) that H = (h R -h L )B, with B a bounded function. The above result is established by arguing the boundedness (1.3) of the Froude number. Therefore, property (H-3) is satisfied by H, when H is given by (4.5).

Lemma 4 .

 4 If the correct solution H of the polynomial equation (4.1) is chosen, then H satisfies the consistency property (H-1) and the well-balanced property (H-2).

Figure 2 :

 2 Figure 2: Experiment from Section 6.1: values of h (left panel) and q (right panel) at time t end with 40 cells.

- 4 • 10 - 15 - 2 • 10 - 15 0 2 •

 41015210152 10 -15

• 10 -

 10 17 1.24 • 10 -16 3.59• 10 -17

Figure 6 :

 6 Figure 6: Subcritical flow from Section 6.2.2: free surface h+Z (left panel) and discharge q (right panel), displayed at time t end = 500 with 75 cells.

q 3 .

 3 74 • 10 -2 4.73 • 10 -14 5.15 • 10 -14 5.21 • 10 -14 error on B 1.45 • 10 -1 4.50 • 10 -14 5.12 • 10 -14 5.92 • 10 -14

Figure 7 :

 7 Figure 7: Transcritical flow without shock from Section 6.2.2: free surface h + Z (left panel) and discharge q (right panel), displayed at time t end = 125 with 75 cells.

Figure 8 :

 8 Figure 8: Transcritical flow with shock from Section 6.2.2: free surface h + Z (left panel) and discharge q (right panel), displayed at time t end = 1000 with 75 cells.

Figure 9 :

 9 Figure 9: Wet dam-break from Section 6.3: free surface h + Z (left panel) and discharge q (right panel), displayed at time t end = 0.05 with 50 cells.

Figure 10 :

 10 Figure 10: Dry dam-break from Section 6.3: free surface (left panel) and discharge q (right panel), displayed at time t end = 0.075 with 50 cells.

Figure 11 :

 11 Figure11: Stationary contact from Section 6.4: free surface h + Z (left panels) and discharge q (right panels), displayed at time t end = 0.075 with 50 cells (top panels) and 2000 cells (bottom panels). The bottom panels are zoomed in on the interval (0.46, 0.54).

  P2error on h 8.88 • 10 -17 2.01 • 10 -16 1.09 • 10 -16 4.44 • 10 -17 error on q 5.25 • 10 -16 1.42 • 10 -15 2.32 • 10 -15 1.61 • 10 -15

	HSR, P0	HDR, P0	HDR, P1	HDR,

Table 3 :

 3 Lake at rest with submerged bottom from Section 6.2.1: errors (in L 2 norm) between the initial condition and the approximate solutions at time t end = 1, using 50 cells.

	h + Z
	1
	0.5								2 • 10 -16
									0
	0	0	0.2	0.4	0.6	0.8	1	x	-2 • 10 -16	0	0.2	0.4	0.6	0.8	1

  Lake at rest with emerged bottom from Section 6.2.1: free surface h + Z (left panel) and discharge q (right panel), displayed at time t end = 1 with 50 cells. P2 error on h 1.85 • 10 -17 2.75 • 10 -17 3.07 • 10 -17 1.32• 10 -17

			q
			4 • 10 -16
					x
	HSR1	HDR1	HDR2	HDR3
	Figure 5: HSR, P0	HDR, P0	HDR, P1	HDR,

Table 5 :

 5 t) = H 0 if the flow is subcritical; otherwise, homogeneous Neumann boundary conditions are prescribed on h. Setup of the experiments from Section 6.2.2.

	experiment	figure	q 0	H 0	t end
	subcritical Figure 6 4.42	2	500
	transcritical Figure 7 1.53 0.66	125
	transcritical with shock Figure 8 0.18 0.33 1000

Table 7 :

 7 Transcritical flow without shock from Section 6.2.2: errors (in L 2 norm) between the exact steady state and the approximate solutions at time t end = 125, using 75 cells.

Table 8 :

 8 Setup of the experiments from Section 6.3.

				experiment	figure			h(x, 0) + Z(x)	t end	θ C
				wet dam-break Figure 9			1.5 if x < 0.5 1 otherwise	0.05	0.15
				dry dam-break Figure 10		1 Z(x) otherwise if x < 0.5	0.075	0.1
	h + Z								q
	1.5								0.8
									0.6
	1
									0.4
	0.5								0.2
	0	0	0.2	0.4	0.6	0.8	1	x	0	0	0.2	0.4	0.6	0.8	1	x
					reference	HSR1			HDR1	HDR2	HDR3

Order of accuracy

The first round of numerical experiments consists in measuring the order of accuracy. To that end, we introduce this useful compactly supported C ∞ bump function:

Note that ω(x) vanishes when |x -1/2| ≥ 1/4. We take Z(x) = ω(x), and the initial condition is given by h(x, 0) = 2 -Z(x) + cos 2 (2πx) and q(x, 0) = sin(2πx).

These initial data are evolved until the final time t end = 5 • 10 -3 , chosen in order to avoid the formation of discontinuities, so as to be able to compute the order of accuracy. A reference solution, to which the results of the schemes are compared, is computed with 20 • 2 12 = 81920 cells. Periodic boundary conditions are prescribed.