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

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


Introduction.
We are interested in the numerical approximation of the wellknown Shallow Water Equations (SWE), given by ( 1) where z(x) denotes a given smooth topography and g > 0 is the gravity constant.The primitive variables are the water depth h ≥ 0 and its velocity u, which both depend on the space and time variables, respectively x ∈ R and t ∈ [0, ∞).At time t = 0, we assume that the initial water depth h(x, t = 0) = h 0 (x) and velocity u(x, t = 0) = u 0 (x) are given.In order to shorten the notations, we will use the following condensed form of (1), namely (2) where U = (h, hu) T , F(U) = (hu, hu 2 + gh 2 /2) T and S(U, z) = (0, −gh∂ x z) T .In case of necessity and since the topography is assumed to be constant in time in this work, we will also make use of the following equivalent form of (1), namely and V = (h, hu, z) T .At last, recall that (1) is strictly hyperbolic over the phase space Ω = {(h, hu) T ∈ R 2 | h > 0} and that the eigenstructure is composed by two genuinely non linear characteristic fields associated with the eigenvalues {u − c, u + c} where the sound speed is defined by c = √ gh.We recall also that the regions where u 2 < c 2 (resp. u 2 > c 2 ) are called subcritical (resp.supercritical).
As already stated, we are especially interested in the design of a numerical scheme satisfying the so-called fully well-balanced property, meaning that it should be able to preserve discrete approximations of the smooth stationary solutions of (1).These steady states are governed by the ordinary differential system ∂ x F (U) = S(U, z), namely Recall that the well-known "lake at rest" solution corresponds to (5) h + z = constant, u = 0, and that numerical schemes preserving these particular stationary solutions are said to be well-balanced (see for instance the recent book [26] for a review).Here we are also interested in preserving equilibriums (4) with a non-zero velocity, leading to the so-called fully well-balanced property.
In order to reach this objective, we propose to use a strategy based on a Lagrange-Projection (or equivalently Lagrange-Remap) decomposition, which allows to naturally decouple the acoustic and transport terms of the model.Such a decomposition proved to be useful and very efficient to deal with subsonic or near low-Froude number flows.
In this case, the usual CFL time step limitation of Godunov-type schemes is indeed driven by the acoustic waves and can thus be very restrictive.The Lagrange-Projection strategy allows to design a very natural implicit-explicit and large time-step schemes with a CFL restriction based on the (slow) transport waves and not on the (fast) acoustic waves.In this paper, we will stay focused on explicit schemes but the reader is referred to the pioneering work [20] and the more recent ones [14], [15], [16], [17] for more details.
Note that the large time step Lagrange-Projection scheme proposed in [17] for the SWE is well-balanced, but not fully well-balanced.
There is a huge amount of works about the design of numerical schemes for the SWE, but most of them intend to satisfy the well-balanced property, but not the full well-balanced property.To mention only a few of them, we refer for instance the reader to the following well-known contributions [3], [27], [30], [22], [23], [1], [19], [6], [7], [2].
The outline of the paper is as follows.In Section 2 we give a general introduction of the Lagrange-Projection approach in the case of flat topography.This technique will be generalized in Section 3 for the case of arbitrary topography in order to develop a fully well-balanced scheme.Finally, in Section 4, some numerical results will be shown in order to test this new scheme.

The Lagrange-Projection decomposition.
The aim of this section is to briefly present the Lagrange-Projection splitting strategy considered in this paper and to introduce the basics of the numerical approximation.
Let us begin with the Lagrange-Projection decomposition.By using the chain rule for the space derivatives, one first writes system (1) under the following equivalent form for smooth solutions, namely Then, the splitting simply consists in first solving the Lagrangian form of this system, namely which gives in Lagrangian coordinates, with τ = 1/h and τ ∂ and then the transport system (7) Remark 2.1.This splitting procedure can be interpreted as a two-step algorithm consisting in solving first the system in Lagrangian coordinates and then projecting the results in the Eulerian coordinates.See [24] for further details.

The Lagrange-Projection numerical algorithm.
Space and time will be discretized using a space step ∆x and a time step ∆t into a set of cells [x j−1/2 , x j+1/2 ) and instants t n+1 = n∆t, where x j+1/2 = j∆x and x j = (x j−1/2 + x j+1/2 )/2 are respectively the cell interfaces and cell centers, for j ∈ Z and n ∈ N.For a given initial condition x → U 0 (x), we will consider a discrete initial data U 0 j which approximates x j−1/2 U 0 (x) dx, for j ∈ Z.Therefore, the proposed algorithm aims at computing an approximation U n j of 1 ∆x x j+1/2 x j−1/2 U(x, t n ) dx where x → U(x, t n ) is the exact solution of (1) at all time t n , n ∈ N. Given the sequence {U n j } j , it is a matter of defining the sequence {U n+1 j } j , n ∈ N, since {U 0 j } j is assumed to be known.
Using these notations, the overall Lagrange-Projection algorithm can be described as follows: for a given discrete state U n j = (h, hu) n j , j ∈ Z that describes the system at instant t n , the computation of the approximation U n+1 j = (h, hu) n+1 j at the next time level is a two-step process defined by 1. Update U n j to U n+1− j by approximating the solution of (6), 2. Update U n+1− j to U n+1 j by approximating the solution of (7).

The case of a flat topography.
In order to prepare the forthcoming developments, we first consider the simple situation of a flat topography, which corresponds to the case ∂ x z = 0.In this case and using standard notations, the Lagrangian step writes (see [24]) and p n j = g(h n j ) 2 /2 for all j, or equivalently (8) Here, the constant a j+1/2 has to be chosen sufficiently large for the sake of stability, and more precisely larger than the Lagrangian sound speed hc according to the well-known subcharacteristic condition.In practice it is required that a i+1/2 is greater than the values h √ gh at the interface.The definition of a j+1/2 will be given later on and we also refer for instance the reader to [9,24,14,16] for more details.Moreover, the Lagrangian step is stable under the Courant-Friedrichs-Lewy (CFL) condition ∆t max As far as the the transport step is concerned, an upwind scheme is applied where the value u * j+1/2 is the velocity at the interface where u * ,+ j±1/2 = max(u * j±1/2 , 0) and u * ,− j±1/2 = min(u * j±1/2 , 0) for all j.
The CFL condition associated with the transport step reads ∆t max Easy calculations show that the system can be written in the equivalent form where which gives that the whole scheme is conservative in the usual sense of finite volume methods and writes These formulas will be useful in the next sections to deal with a generalization of this scheme to the case of a non flat topography and impose the fully well-balanced property.

An equivalent formulation.
In this section, we propose an equivalent formulation of the previous scheme based on the Reynolds' transport theorem which can be expressed for any scalar value X as where dV is a volume element at point x (the integration variable), and v(x, t) is the velocity at which the domain Ω(t) is moving.In the one dimension case considered in this paper, and considering that Ω(t) is an interval with limits a(t) and b(t), the Reynolds's theorem reduces to This theorem expresses the time variation of the amount of X contained inside the volume Ω(t) which moves at velocity v(x, t).
The Lagrangian coordinates.Considering that the velocity v(x, t) equals the fluid velocity u(x, t) in the SWE, the Reynolds' theorem allows to describe the fluid motion in the Lagrangian coordinates, assuming that the observer moves with the fluid flow.
Considering now successively that X = 1, h, hu and denoting by V (t) the infinitesimal control volume, we get where each point x(t; ξ) of the control volume is given by the solution of the following ordinary differential equation, with initial condition ξ, Note that with the above notations, and so that for an infinitesimal control volume one can write If we now denote α(ξ, t) = α(x(t; ξ), t) for any scalar quantity α, the previous relations also write The Lagrangian step.Let us now turn back to the numerical setting and consider for a given j that the control volume is defined at time t = 0 by the cell (x j−1/2 , x j+1/2 ).
The endpoints x ± (t; x j±1/2 ) of the control volume at a given time t are thus defined by the solutions of the following ordinary differential equations ( 12) Denoting h n+1− j and (hu) n+1− j constant approximate values of h and (hu) at time t n+1 on the moving domain whose volume is V j (∆t), a simple way to discretize (9) would be ( 13) Therefore, one clearly sees that the discretizations ( 13) and ( 14) of ( 9) (or equivalently (11)) allow to recover (8).At last, note now that from (12), , and by subtraction we also have and therefore The transport step.Let us now turn to the transport or remap step.In the framework of the present equivalent reformulation of our algorithm, we have now at hand some constant approximations denoted h n+1− j and (hu) n+1− j of h and (hu) at time t n+1 on the new cells (x − (∆t; x j−1/2 ), x + (∆t; x j+1/2 )) and whose volume is V j (∆t).In this context, the aim of the remap step is just to project these values on the original mesh, which simply writes with clear notations and X = h, hu, where we have set for all j It is easy to check that these final values coincide with the ones of the previous formulation of the scheme.
3. The Lagrange-Projection scheme with a non flat topography.Now we would like to extend this approach to the case of a non flat bottom paying a particular attention to the well-balanced property.Here the Lagrange-Projection approach aims at solving successively the Lagrangian form of our system, namely + gh∂ m z = 0, and the transport system At this stage, it is important to have in mind that in order to get the fully well-balanced property, we are going to consider in-cell reconstructions in both Lagrangian and Transport steps.Moreover, these reconstructions will make both steps of the strategy to be fully well-balanced for (hu), in the sense that if the solution at time t n is an equilibrium state, we will have (hu) n+1 j = (hu) n+1− j = (hu) n j for all j.This is clearly sufficient but not necessary as the fully well-balanced property simply asks for (hu) n+1 j = (hu) n j for all j.As far as h is concerned, we will have h n+1− j = h n j in general, but the overall algorithm will be fully well-balanced, namely h n+1 j = h n j for all j.Let us now describe the two steps in details.
and where At this stage, it remains to define u * j+1/2 , p * j+1/2 and {gh∂ x z} j for all j in such a way that the fully well-balanced property holds true for this step.
Reconstruction procedure.In order to get this property, we propose to perform incell reconstructions of the variables at time t n , that is to say just before starting the Lagrangian step.With this in mind, we suggest to understand the cell values h n j and (hu) n j as the averages of stationary solutions defined in the corresponding cell C j .Thus, we aim at defining the functions x → h(x; K 1,j , K 2,j ) and x → u(x; K 1,j , K 2,j ) such that with two constants K 1,j and K 2,j such that In practice, and since we are only interested in a first order scheme, it is more convenient to approximate the above integrals by the mid-point quadrature formula, which allows to identify the cell averages h n j (respectively (hu) n j ) and the mid-point values h(x j ; K 1,j , K 2,j ) (resp.x → u(x j ; K 1,j , K 2,j )).We are thus led to define the functions x → h(x; K 1,j , K 2,j ) and x → u(x; K 1,j , K 2,j ) by ( 16) for all j.Using the U notation, we have thus been able to reconstruct smooth stationary solutions x → U(x; K 1,j , K 2,j ) of the SWE within each cell and the average of which coincides with the solution at time t n for the conservative variables and up to second-order accuracy.These stationary solutions are then used to define left and right extrapolated values at the mesh interfaces by setting for all j Remark that, to close the definition of the reconstructed values, the bottom topography should be defined at the interface.Following the ideas in [1] and [11], let us consider a reconstructed bottom at z j+1/2 , consistent with z(x j+1/2 ), whose precise definition will be discussed later.Thanks to (16) and using the notation q = hu, we consider the reconstructed states h j+1/2± , q j+1/2± given by ( 17) It is important to notice that the values h i+1/2± may not be uniquely defined or not defined at all.Let us recall that this idea of using in-cell reconstruction using equilibria has been used, for instance, in [28,10,11].In particular, let us consider a similar notation and the following proposition in [10]: . There exists a solution (h * , u * ) to (19) if and only if (20) δ + δ 0 (h 0 , q 0 ) ≥ 0, where Moreover, i) If we have equality in (20), there is only one solution (h * , u * ) to (19), it is given by ii) If we have a strict inequality in (20), then there are exactly two different solutions (19), with h * sup ≤ h s < h * sub , and h * sup < h s for q 0 = 0.
iii) A solution (h * , u * ) to (19), with h * u * = 0 satisfies: iv) Assume that δ is such that (20) is satisfied strictly so that there are exactly two solutions h sub and h sup for (19).
This means that when U j and U j+1 are not critical points (δ 0 (h j , q j ) > 0 and δ 0 (h j+1 , q j+1 ) > 0) and the bottom is continuous, we can define the reconstructed values provided that the spatial resolution is high enough.Nevertheless, some problems may arise when we are near a critical point at the discrete level.
b) When the reconstructed states exist, we have in general two possible choices: a subcritical and supercritical value.We will select in general the solution that preserves the same character, that is, for U j subcritical (resp.supercritical) we choose U j+1/2− and U j−1/2+ subcritical (resp.supercritical).This is true far from sonic points and for ∆z sufficiently small.c) In the general case, the choice of z j+1/2 should be taken with care.We require that z j − z j+1/2 ≥ −δ 0 (h j , q j ) and z j+1 − z j+1/2 ≥ −δ 0 (h j+1 , q j+1 ) so that the reconstructed states exist.Moreover, we shall impose that the reconstructed states verify and which will play an essential role for the positivity of the scheme.
To conclude this paragraph, notice that the proposed reconstruction procedure can be understood as a generalization of the well-known hydrostatic reconstruction associated with the lake at rest equilibrium and defined by Let us now define {gh∂ x z} j such that {gh∂ x z} j ≈ 1 ∆x Since we have reconstructed stationary solutions in each cell, we clearly have for all j which, by integrating in the interior of the cell, gives that is to say where of course p n j+1/2− = g(h n j+1/2− ) 2 /2 and p n j−1/2+ = g(h n j−1/2+ ) 2 /2.Since, (hu) n j+1/2− = (hu) n j−1/2+ = (hu) n j , we finally get which is a simple way to account for the reconstruction procedure introduced in the Lagrangian step.It will also guarantee the fully well-balanced property.
Note that following the same idea introduced in Section 2. Regarding now the definition of (hu) n+1− j+1/2 , it is a matter of defining again a suitable approximation of (hu) at time t n+1− .Here, we propose to keep the same definition used in the case of a flat topography, namely which will be sufficient to get the fully well-balanced property, as we show now.

Positivity and fully well-balanced property.
In this section, we prove that the Lagrangian-Projection scheme defined by ( 15)-( 24) is fully well-balanced and that it preserves positivity of water height.
Well-balanced property.We shall check that the scheme defined previously preserves all smooth stationary states.Let us assume that the solution at time t n corresponds to a general steady state of the SWE, that is to say where K 1 and K 2 are two given constants.The reconstruction procedure gives U j+1/2− = U j+1/2+ and in particular As a consequence, the Lagrangian step (15) gives Under the usual CFL restriction, we thus have (hu) n+1− j = (hu) n j so that the Lagrangian step is fully well-balanced for (hu).Note that in general, h n+1− j = h n j , except if u j+1/2 = u j−1/2 gets true, which happens in the particular case of the lake at rest steady state (u n j = 0 for all j).We also have (25) gives Therefore, the scheme is thus fully well-balanced since h n+1 j = h n j and (hu) n+1 j = (hu) n j for all j.
Positivity property.It is clear that the Lagrangian step (15) preserves positivity of water thickness provided that L n j > 0, which will be true under a suitable CFL condition.The question is then whether the transport step preserves positivity or not.Following the ideas presented in [1], given a numerical scheme for water height in the form the scheme will preserve the positivity of h if whenever h n j = 0, one has The fully well-balanced Lagrange-Projection for variable h writes Now, definitions given in Remark 3.2.c) grants that when Conversely, when U n j is supercritical, we have It follows then that whenever h n j = 0 we get In that case, we get Now, it is easy to check that in this case and the result follows.
The simulation is carried out with 200 cells and open boundary conditions.As expected, this steady state is preserved by the scheme with an error 10 −12 .Figure 1 shows the surface elevation as well as the discharge at time t = 1 which corresponds just to the given initial condition.

Generation of subcritical steady state. The objective of this test is to
study the convergence of the scheme to steady states.To do so let us consider the bottom topography (27) and we consider as initial condition a lake at rest situation: h + z = 1, u = 0. We impose at the left boundary the condition hu(x = 0, t) = 0.5 and at the right boundary the condition h(x = 1, t) = 1.It is expected that, for long time, the solution will converge to the steady state presented previously.The simulation is carried out with 200 cells and the time evolution of the surface and discharge are shown in Figure 2.Here a comparison with the solutions computed with 800 cells are also shown.Figure 3 shows the steady state reached for t = 50 and comparison with the exact reference solution.Moreover, 10 −12 .Figure 4 shows the deviation from the expected constant values for q and u 2 2 + g(h + z) at the steady state.We see that the for the second.

Transcritical continuous steady state.
As it is known, the generalized hydrostatic reconstruction procedure may present some problems at critical points where the reconstructed states may not be well defined.We propose now to take as initial condition a transcritical steady state to check how this situation is handled by the scheme.We consider again the bottom topography (27)   .
We consider the subcritical solution for x < 0.5 and the supercritical solution for x > 0.5.
Previous values guarantee that these solutions exist and there is a critical point at x = 0.5.
Figure 5 shows the surface elevation for time t = 1 using 201 cells.We remark that the scheme behaves well at x = 0.5 where the critical point is located.Let us recall that critical points may be problematic for the reconstruction using equilibria.This is due to the fact that at discrete level we represent the exact smooth functions by a piece-wise constant function consisting on the cell averages which gives in practice a jump on the bottom at the interfaces.Nevertheless we check in Figure 5(b) that the scheme captures correctly the solution.Figure 6 shows the deviation from the expected constant values for q and u 2 2 + g(h + z) at the steady state.We see that the error is of order 10 −14 for the first one and 10 −13 for the second which shows that the steady state is preserved.Figures 7 and 8 show the evolution of the surface and the discharge at the beginning of the simulation, where the perturbation propagates along the domain.Figure 9 shows the surface elevation and discharge at time t = 10 when the steady state is reached.
Compared to the prescribed initial steady state, we remark that the perturbation has not affected the reached steady state.Nevertheless, a small artifact is seen at the critical point in Figure 9(b) although it is very small.This small artifact reduces as the number of cells is increased.
Consider now a smaller perturbation of 0.0001 for 0.15 < x < 0.2 and let us compare the scheme provided here with a similar based on just an hydrostatic reconstruction [1].
While the former is fully well-balanced, the latter is only well-balanced for lake at rest steady states.The simulations are done using 3201 volume cells.Figure 10 shows at time t = 0.05 the difference between the computed values h and q with respect to the unperturbed steady-state h 0 and q 0 described in Subsection 4.3.As we see, just from the beginning, the non fully well-balanced scheme is producing a perturbation in the   middle of the domain which is even bigger than the perturbation introduced on the initial condition.Nevertheless, the fully well-balanced scheme introduced here does not present this problem.

4.5.
Transcritical steady state with shock.We consider now a numerical test taken from [9] consisting of a transcritical flow with a shock over a bump.The interval is [0, 25], the initial condition and depth bottom is given by (28) h(x, t = 0) = 0.33, hu(x, t = 0) = 0.18, and the boundary conditions are given by hu(x = 0, t) = 0.18, h(x = 25, t) = 0.33.
We consider a uniform mesh composed by 201 cells.
The initial condition will evolve until a transcritical solution with a stationary shock is developed.Figure 11 shows the surface and discharge evolution obtained with the scheme proposed here.In Figure 12 we see the surface profile obtained and a comparison for a reference solution with increasing number of cells.We remark there is a small problem for the scheme at the critical point near x = 10 and the position of the shock   does not exactly coincide with reference solution.Nevertheless those problems disappear when the mesh is refined and the numerical solutions converge to the reference solution.
Figure 13 shows the steady state conditions and comparison with a reference solution for different number of cells.4.6.Formal convergence of the scheme.We intend now to study the formal convergence of the scheme proposed here.To do so, let us consider the initial condition given by z(x) = 0.1 cos(2πx), h(x, 0) = 1.1 + 0.1 sin(4πx) − z(x), hu(x, 0) = 0.0.
We compute a reference solution using 3200 points in the interval [0, 1] and using periodic boundary conditions up to time t = 0.2.Then we make a comparison with the solutions given by different number of cells.The L 1 errors are shown in Table 1.As expected, the scheme presented here is a first order scheme.

Traveling wave over an obstacle.
We consider now bottom topography and initial condition given by (29) z(x) = 0.5, if 0.6 < x < 0.7, 0, otherwise, h(x, t = 0) = 0.5, if 0.6 < x < 0.7, 1, otherwise, and hu(x, t = 0) = 0.This will produce a rarefaction followed by a shock that travels towards the obstacle at the bottom.We compute the solution using 201 cells in the   14 where dashed lines correspond to the results using the finer grid and continuous line to the coarser grid.The final time is shown in 15 once the waves have left the domain and the steady state is reached.We see that scheme is able to reproduce the evolution of the waves in the flat region as well as when the shock encounters the obstacle.

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

3. 1 .
The Lagrangian step.. Considering first the Lagrangian step, following the ideas used in the flat topography case, we propose the following generalization(15)

4. Numerical results. 4 . 1 .
Subcritical steady state.In this first test we will verify the fully well balanced property of the scheme.To do so, let us consider a bottom topography(27) z(x) = 0.5 exp − (x − 0.5) 2 0.005 , and the subcritical steady state solution verifying

Figure 2 .Figure 3 .
Figure 2. Generation of subcritical steady state from stationary initial condition: surface and discharge evolution.Continuous lines correspond to the solution computed with 200 cells and discontinuous to the one computed with 800 cells

Figure 9 .
Figure 9. Perturbation of a steady state at times t = 10

Figure 15 .
Figure 15.Travelling wave over an obstacle: surface elevation and discharge for time t = 1 3, we may interpret again the transport step as the projection of the values in Lagrangian coordinates into the Figure 10.Perturbation of a steady state.Comparison with the unperturbed steady state at time t = 0.05 for fully and non fully well-balanced schemes

Table 1 L
Figure 13.Transcritical steady state with shock: steady state condition for different number of cells 1 errors and numerical orders of accuracy.
(b) Zoom on the surface for different number of cells Figure 12.Transcritical steady state with shock: surface elevation