Christophe Berthon

Manuel J Castro Díaz

Tomás Arnaud Duran

Tomás Morales De Luna

Khaled Saleh

Tom Ás

Morales De

Artificial viscosity to get both robustness and discrete entropy inequalities

Keywords: AMS subject classifications. 65M08, 65M12 Evolution laws, finite volume schemes, CFL-like condition, robustness, discrete entropy inequalities

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1. Introduction. The present work concerns the stability of numerical schemes to approximate the weak solutions of a N × N evolution law systems in the following form:

∂ t w + ∂ x f (w) + A(w)∂ x w = 0, x ∈ R, t > 0, (1.1)
where (x, t) ∈ R + × R → w(x, t) denotes the unknown vector assumed to belong to an open convex set Ω ∈ R N . Here, f : Ω → R N is a smooth enough flux function and A : Ω → R N × R N denotes a given smooth matrix. From now on, it is worth noticing that we do not impose the matrix ∇ w f (w) + A(w) to be diagonalizable in R. As a consequence, in this work, the system (1.1) is not necessarily hyperbolic over Ω. Now, an important assumption must be put on the exact Riemann solutions of (1.1), denoted W ex x t ; w L , w R , for an initial data given by

w 0 (x) = w L , if x < 0, w R , if x > 0,
where w L and w R are two given constant states in Ω. Indeed, we assume the existence of λ > 0 such that the exact Riemann solution is in the form W ex x t ; w L , w R = w L for all x < -λ t, W ex x t ; w L , w R = w R for all x > λ t.

(1.2)

In fact, such an assumption is immediately satisfied by any hyperbolic system of conservation laws and it is, in general, assumed for non-conservative systems with real or complex eigenvalues.

In addition, we assume the system (1.1) to be endowed with an entropy inequality given by ∂ t η(w) + ∂ x G(w) ≤ 0, in the weak sense, (1.3) where the entropy function η : Ω → R is convex and the entropy flux function G : Ω → R is defined by ∇ w G(w) = ∇ w η(w) • (∇ w f (w) + A(w)). We recall that the existence of entropy inequalities is strongly related to the hyperbolicity of the system as long as (1.1) is in conservation form, namely A ≡ 0. However, by adopting a non-conservative system, with A ≡ 0, the system (1.1) may simultaneously get eigenvalues in C and admit entropy inequalities. These entropy inequalities are essential to correctly define the discontinuous shock solutions (see [START_REF] Dafermos | Hyperbolic conservation laws in continuum physics, volume 325 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences[END_REF][START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Lax | Shock waves and entropy[END_REF][START_REF] Lax | Hyperbolic systems of conservation laws and the mathematical theory of shock waves[END_REF][START_REF] Serre | Systems of conservation laws[END_REF]) as long as (1.1) is given in conservation form. Now, after [START_REF] Berthon | Why many theories of shock waves are necessary: kinetic relations for non-conservative systems[END_REF][START_REF] Castro | Why many theories of shock waves are necessary: convergence error in formally path-consistent schemes[END_REF][START_REF] Maso | Definition and weak stability of nonconservative products[END_REF][START_REF] Hou | Why nonconservative schemes converge to wrong solutions: error analysis[END_REF][START_REF] Lefloch | Shock waves for nonlinear hyperbolic systems in nonconservative form[END_REF], it is well-known that entropy inequalities, in general, are not sufficient to suitably define shock waves as long as the matrix ∇ w f (w)+A(w) never recasts as a Jacobian matrix. Indeed, because of the non-conservative product A(w)∂ x w, the system (1.1) presents a major difficulty because the discontinuous solutions have no sense in the distributional framework. However, to give a sense to the weak solutions, Dal Maso, LeFLoch and Murat [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF] have introduced a suitable theory, the so-called path-theory. Now, for both conservative and non-conservative systems, the entropy inequalities produce important estimations since we get, in a sense to be prescribed, η(w(•, t)) L 2 (R) ≤ η(w 0 ) L 2 (R) where w 0 stands for a given initial data.

During the last half century, numerous numerical schemes have been designed to approximate the weak solutions of (1.1). Clearly, it is not possible to give an exhaustive bibliography, but the reader is referred to the pioneer work by Parés [START_REF] Parés | Numerical methods for nonconservative hyperbolic systems: a theoretical framework[END_REF] where a path-consistent scheme is derived in order to produce relevant approximations of the weak solutions of (1.1) according to the path-theory [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF] (see also [14,[START_REF] Godlewski | Hyperbolic systems of conservation laws, volume 3/4 of Mathématiques & Applications[END_REF][START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws: Exponential-fit, Well-balanced and Asymptotic-preserving[END_REF][START_REF] Kröner | Numerical schemes for conservation laws[END_REF][START_REF] Lefloch | Shock waves for nonlinear hyperbolic systems in nonconservative form[END_REF][START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF] and references therein for a review about the numerical approximation of hyperbolic systems of conservation laws, balance laws and hyperbolic non-conservative systems). In the present work, we focus our attention on the well-known 3-point finite volume explicit schemes. These numerical techniques read

w n+1 i = w n i - ∆t ∆x f ∆ (w n i , w n i+1 ) -f ∆ (w n i-1 , w n i ) - ∆t 2∆x A L ∆ (w n i , w n i+1 ) • (w n i+1 -w n i ) + A R ∆ (w n i-1 , w n i ) • (w n i -w n i-1 ) , (1.4
) where f ∆ : Ω × Ω → R N denotes the Lipschitz-continuous numerical flux function such that f ∆ (w, w) = f (w) for all w in Ω. Moreover, A L,R

∆ : Ω × Ω → R N × R N is a matrix consistent with A(•), namely A L,R
∆ (w, w) = A(w). Here, w n i approximates w(x, t) for all x in a cell (x i- 1 2 , x i+ 1 2 ) of size ∆x at time t n . For the sake of simplicity in the forthcoming developments, both space and time increments ∆x and ∆t are constant, and we set x i = (x i- 1 2 + x i+ 1 2 )/2 the middle of the cell (x i-1 2 , x i+ 1 2 ). We underline that the definition of f ∆ and A L,R ∆ are free and the reader may refer to any scheme derivation proposed in the literature. Now, equipped with these numerical approximations, three natural questions arise:

1. How is the time increment ∆t restricted? 2. Do the updated states (w n+1 i ) i∈Z be in Ω as soon as w n i ∈ Ω for all i ∈ Z? 3. How can we restore, at the numerical level, an entropy inequality (1.3)?

Concerning the first question, the well-known CFL restriction produces an interesting response. Indeed, after [START_REF] Neumann | A method for the numerical calculation of hydrodynamic shocks[END_REF], stability criterion can be imposed as follows:

∆t ∆x max i∈Z, 1≤ ≤N |λ (w n i )| ≤ 1 2 , (1.5)
where (λ (w)) 1≤ ≤N denote the eigenvalues of the matrix ∇ w f (w) + A(w). According to the definition of f ∆ (•, •) and A L,R ∆ (•, •), such a time restriction may be insufficient to ensure a stability condition in a sense to be prescribed. Moreover, in (1.5), the eigenvalues are expected to be real. As a consequence, it is worth noticing that the CFL-like condition (1.5) is suitable for hyperbolic systems with explicit eigenvalues. Unfortunately, systems issuing from sophisticated physics do not necessarily admit real eigenvalues. For instance, let us consider the bi-layer shallow water model [10-13, 17,18] or sediment transport model [16,[START_REF] Cordier | Bedload transport in shallow water models: Why splitting (may) fail, how hyperbolicity (can) help[END_REF][START_REF] González-Aguirre | A robust model for rapidly varying flows over movable bottom with suspended and bedload transport: Modelling and numerical approach[END_REF][START_REF] Morales De Luna | A duality method for sediment transport based on a modified meyer-peter & müller model[END_REF] for which the associated eigenvalues do not have an explicit simple expression and belong to R or C.

Next, concerning the second question, this problem is frequently called realizability or robustness of the scheme. This issue is, in general, tackled in the derivation of a new scheme. Indeed, if the numerical robustness is lost, no longer the existence of the numerical approximation can be ensured. As a consequence, the robustness (or realizability) must be preserved to validate a scheme.

Finally, concerning the last question, the derivation of discrete entropy inequalities remains a very difficult task. Indeed, to enforce the required entropy stability, the updated states (w n+1 i ) i∈Z , given by (1.4), are expected to satisfy for all i ∈ Z

η(w n+1 i ) ≤ η(w n i ) - ∆t ∆x G ∆ (w n i , w n i+1 ) -G ∆ (w n i-1 , w n i ) , (1.6)
where G ∆ (•, •) is a numerical entropy flux function such that G ∆ (w, w) = G(w) for all w in Ω. As long as N ≥ 2, up to our knowledge, very few schemes are entropy preserving and satisfy a discrete entropy inequality in the form (1.6). The original Godunov scheme [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF] (see also [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF]) and the HLL scheme [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] are proved to be entropy preserving. The establishment of this stability property comes from the entropy inequality (1.3) satisfied by the exact solution. Some kinetic schemes [START_REF] Bouchut | Introduction to the mathematical theory of kinetic equations[END_REF][START_REF] Khobalatte | Maximum principle on the entropy and second-order kinetic schemes[END_REF] and relaxation schemes [7-9, 19, 20], including the well-known HLLC scheme [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Toro | Restoration of the contact surface in the hll-riemann solver[END_REF], are also proved to preserve the entropy stability, by adopting some convex minimization principles. Several works propose more general conditions to get the expected discrete entropy inequalities when considering conservative hyperbolic systems. For instance, in [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF], the authors introduce a local entropy condition per interface. Excepted for the Godunov scheme and the one-intermediate state HLL scheme (see also [START_REF] Berthon | A local entropy minimum principle for deriving entropy preserving schemes[END_REF]), this condition is hardly reachable. We also mention the important work by Tadmor [START_REF] Tadmor | Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems[END_REF][START_REF] Tadmor | Entropy stable schemes. Handbook of Numerical Methods for Hyperbolic Problems: Basic and Fundamental Issues[END_REF] where, arguing the entropy variables, discrete entropy inequalities (1.6) are derived according to a suitable control of the numerical viscosity. We also refer to the work by Bouchut [START_REF] Bouchut | Entropy satisfying flux vector splittings and kinetic bgk models[END_REF] where the discrete entropy inequalities (1.6) are obtained for hyperbolic system of conservation laws by arguing a suitable kinetic reformulation. In fact, most of the designed schemes perform good approximations, without spurious numerical perturbations, so that they are certainly entropy preserving but the establishment of the required entropy stability is not reachable. However, it is important to notice that several schemes are entropy violating (for instance, see [START_REF] Gallouët | Some recent finite volume schemes to compute Euler equations using real gas EOS[END_REF][START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF][START_REF] Masella | On a rough godunov scheme[END_REF][START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF]) and may capture wrong shock solutions. In the literature [START_REF] Dubois | A non-parameterized entropy correction for roe's approximate riemann solver[END_REF][START_REF] Helluy | A simple parameter-free entropy correction for approximate Riemann solvers[END_REF], entropy corrections have been introduced but fully discrete entropy inequalities (1.6) are, in general, not established.

In the present work, we address the three asked questions by adopting artificial viscosity. Following ideas introduced by Tadmor [START_REF] Tadmor | Numerical viscosity and the entropy condition for conservative difference schemes[END_REF][START_REF] Tadmor | The numerical viscosity of entropy stable schemes for systems of conservation laws[END_REF][START_REF] Tadmor | Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems[END_REF][START_REF] Tadmor | Entropy stable schemes. Handbook of Numerical Methods for Hyperbolic Problems: Basic and Fundamental Issues[END_REF], and its generalization to non-conservative systems in [15], the numerical viscosity is a essential tool to recover the required entropy estimation (1.6). As a consequence, the scheme improved by artificial viscosity under consideration is here given by

w n+1 i = w n i - ∆t ∆x f ∆ (w n i , w n i+1 ) -f ∆ (w n i-1 , w n i ) - ∆t 2∆x A L ∆ (w n i , w n i+1 ) • (w n i+1 -w n i ) + A R ∆ (w n i-1 , w n i ) • (w n i -w n i-1 ) + γ 2 ∆t ∆x w n i+1 -2w n i + w n i-1 , (1.7
) where γ ≥ 0 is a parameter to be fixed according to the expected stability properties. To judiciously control γ, we here suggest to derive a local Godunov-type reformulation [START_REF] Berthon | Improvement of the hydrostatic reconstruction scheme to get fully discrete entropy inequalities[END_REF][START_REF] Berthon | Continuum Mechanics, Applied Mathematics and Scientific Computing : Godunov's Legacy, chapter An easy control of the artificial numerical viscosity to get discrete entropy inequalities when approximating hyperbolic systems of conservation laws[END_REF][START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]. We adopt the local integral consistency stated by Harten, Lax and van Leer [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]Theorem 3.1] to get the required control of the artificial viscosity and to obtain the expected estimation (1.6) in the same spirit of Tadmor in [START_REF] Tadmor | Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems[END_REF] and its generalization to non-conservative systems [15] . Moreover, this control of the artificial viscosity comes with a relevant restriction of the time increment ∆t by exhibiting a CFL-like criterion. In fact, the local Godunov-type reformulation of the scheme (1.4) also leads to a condition to enforce the updated state w n+1 i to belong to Ω. This condition is, once again, controlled by the artificial viscosity.

The present paper is organized as follows. In the next section, we present local Godunov-type reformulations [START_REF] Berthon | Improvement of the hydrostatic reconstruction scheme to get fully discrete entropy inequalities[END_REF][START_REF] Berthon | Continuum Mechanics, Applied Mathematics and Scientific Computing : Godunov's Legacy, chapter An easy control of the artificial numerical viscosity to get discrete entropy inequalities when approximating hyperbolic systems of conservation laws[END_REF][START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] and the introduction of judicious approximate Riemann solvers. We also give assumptions to be satisfied by the exact Riemann solution for (1.1) in order to state our main result where we show a suitable CFL condition to get the expected robustness and the discrete entropy inequalities. Next, Section 3 and Section 4 are devoted to the establishment of the main result. Both robustness (Section 3) and discrete entropy inequalities (Section 4) are obtained by involving a relevant Godunov-type reformulation. The last section is devoted to numerical experiments in order to illustrate the relevance of the proposed stability enforcing technique. First, we perform numerical approximations of the solutions of the isentropic gas dynamic model and the shallow water model. We show the ability of the here designed technique to stabilize entropy violating schemes. Next, we consider the bi-layer shallow-water model. Because the eigenvalues of the model remain unknown, the interest of the stabilizing technique is twofold since, in addition to both robustness and entropy preserving properties, we obtain a natural CFL-like condition.

2. Godunov-type reformulation and main results. We here reformulate the viscous scheme (1.7) as a Godunov-type scheme. We address the reader to [14] and the references therein for further detail. To do so, we first introduce an approximate Riemann solver as follows:

WR x t ; w n i , w n i+1 =                          w n i , if x < -(λ i+ 1 2 + γ)t, wi+ 1 2 , if -(λ i+ 1 2 + γ)t < x < -λ i+ 1 2 t, w L i+ 1 2 , if -λ i+ 1 2 t < x < 0, w R i+ 1 2 , if 0 < x < λ i+ 1 2 t, wi+ 1 2 , if λ i+ 1 2 t < x < (λ i+ 1 2 + γ)t, w n i+1 , if x > (λ i+ 1 2 + γ)t, (2.1)
where w L i+ 1 2 and w R i+ 1 2 are defined by

w L i+ 1 2 = w n i - 1 λ i+ 1 2 f ∆ (w n i , w n i+1 ) -f (w n i ) - 1 2λ i+ 1 2 A L ∆ (w n i , w n i+1 ) • (w n i+1 -w n i ), (2.2) w R i+ 1 2 = w n i+1 + 1 λ i+ 1 2 f ∆ (w n i , w n i+1 ) -f (w n i+1 ) - 1 2λ i+ 1 2 A R ∆ (w n i , w n i+1 ) • (w n i+1 -w n i ), (2.3)
where we have set

wi+ 1 2 = 1 2 (w n i + w n i+1 ). (2.4)
At this level, the wave speed parameter λ i+ 1 2 is only imposed to be positive but it must be selected according to forthcoming stability conditions. Now, we are able to exhibit the expected Godunov-type reformulation. Lemma 2.1. Under the CFL-like condition ∆t ∆x max

i∈Z λ i+ 1 2 + γ ≤ 1 2 , (2.5)
the viscous scheme (1.7) equivalently reformulates as follows:

w n+1 i = 1 ∆x xi x i-1 2 WR x -x i-1 2 ∆t ; w n i-1 , w n i dx + 1 ∆x x i+ 1 2 xi WR x -x i+ 1 2 ∆t ; w n i , w n i+1 dx. (2.6)
From now on, let us remark that, as soon as γ = 0, the intermediate states wi+ 1 2 no longer appear within the approximate Riemann solver (2.1) so that the Godunov-type reformulation (2.6) coincides with the original scheme (1.4). Moreover, we emphasize that the artificial viscosity involved in (1.7) is entirely contained in the two intermediate states given by wi+ 1 2 . Such a splitting within the intermediate states will allow to control the entropy dissipation rate issuing from the artificial viscosity independently from the entropy dissipation rate of the initial scheme (1.4). Moreover, it is important to notice that this splitting results from the introduction of the free parameter λ i+ 1 2 > 0. In fact, to establish the entropy preserving property, a lower bound of λ i+ 1 2 will be exhibited. Now, we establish the scheme reformulation stated Lemma 2.1. Proof. Direct computations give the following sequence of equalities:

1 ∆x xi x i-1 2 WR x -x i-1 2 ∆t ; w n i-1 , w n i dx = 1 ∆x ∆x/2 0 WR x ∆t ; w n i-1 , w n i dx, = 1 2 w n i + γ ∆t ∆x ( wi+ 1 2 -w n i ) + λ i-1 2 ∆t ∆x (w R i+ 1 2 -w n i ), 1 ∆x x i+ 1 2 xi WR x -x i+ 1 2 ∆t ; w n i , w n i+1 dx. = 1 ∆x 0 -∆x/2 WR x ∆t ; w n i , w n i+1 dx, = 1 2 w n i+1 + γ ∆t ∆x ( wi+ 1 2 -w n i+1 ) + λ i+ 1 2 ∆t ∆x (w L i+ 1 2 -w n i+1 ).
Arguing the definition of the intermediate states, given by (2.2), (2.3) and (2.4), the updated state w n+1 i , defined by (1.7), rewrites as follows:

w n+1 i = 1 ∆x ∆x 2 0 WR x ∆t ; w n i-1 , w n i dx + 1 ∆x 0 -∆x 2 WR x ∆t ; w n i , w n i+1 dx.
Next, under the CFL restriction (2.5), two successive approximate Riemann solver,

namely WR x-x i-1 2 ∆t ; w n i-1 , w n i and WR x-x i+ 1 2 ∆t
; w n i , w n i+1 , never interact. As a consequence, the expected reformulation (2.6) is obtained by a change of variable in the integrals. The proof is thus completed.

In fact, the Godunov-type reformulation is not unique. It will be convenient to present a second reformulation based on the following approximate Riemann solver:

W R x t ; w n i , w n i+1 =              w n i , if x < -(λ i+ 1 2 + γ)t, w L i+ 1 2 , if -(λ i+ 1 2 + γ)t < x < 0, w R i+ 1 2 , if 0 < x < (λ i+ 1 2 + γ)t, w n i+1 , if x > (λ i+ 1 2 + γ)t, (2.7)
where

w L i+ 1 2 = 1 - λ i+ 1 2 λ i+ 1 2 + γ wi+ 1 2 + λ i+ 1 2 λ i+ 1 2 + γ w L i+ 1 2 , (2.8)
w R i+ 1 2 = 1 - λ i+ 1 2 λ i+ 1 2 + γ wi+ 1 2 + λ i+ 1 2 λ i+ 1 2 + γ w R i+ 1 2 . (2.9)
Equipped with this second approximate Riemann solver, we now give a new equivalent Godunov-type reformulation. Lemma 2.2. Under the CFL-like condition (2.5), the viscous scheme (1.7) equivalently reformulates as follows:

w n+1 i = 1 ∆x xi x i-1 2 W R x -x i-1 2 ∆t ; w n i-1 , w n i dx + 1 ∆x x i+ 1 2 xi W R x -x i+ 1 2 ∆t ; w n i , w n i+1 dx.
(2.10)

The establishment of the above result is similar to Lemma 2.1 and the proof is left to the reader.

We are now able to state our main result where we claim that the viscous scheme (1.7) is robust and entropy satisfying as soon as λ i+ 1 2 and γ are fixed large enough. Theorem 2.3. Let (w n i ) i∈Z be given in Ω. Let w n+1 i be given by the viscous scheme (1.7). Assume that

1 2 A L ∆ (w n i , w n i+1 ) + A R ∆ (w n i , w n i+1 ) • (w n i+1 -w n i ) = 1 ∆t ∆t 0 λ i+ 1 2 ∆t -λ i+ 1 2 ∆t A W ex x t ; w n i , w n i+1 ∂ x W ex x t ; w n i , w n i+1 dx dt, (2.11
)

Let λ i+ 1 2 > λ . Let us set E 0 i+ 1 2 = λ i+ 1 2 η(w L i+ 1 2 ) + η(w R i+ 1 2 ) -η(w n i ) -η(w n i+1 ) + G(w n i+1 ) -G(w n i ) , (2.12
)

D i+ 1 2 = 2η( wi+ 1 2 ) -η(w n i ) -η(w n i+1 ), (2.13)
and define γ as follows:

γ = max i∈Z 0, - E 0 i+ 1 2 D i+ 1 2 . ( 2

.14)

Let ∆t be restricted by the CFL-like condition (2.5). Assume that the Hessian matrix ∇ 2 w η(w) is positive definite. Then γ ≥ 0 is bounded and the viscous scheme (1.7) is 1. robust:

w n+1 i
∈ Ω, as long as γ satisfies in addition

γ ≥ max λ i+ 1 2 r i+ 1 2 w L i+ 1 2 -wi+ 1 2 -λ i+ 1 2 , λ i+ 1 2 r i+ 1 2 w R i+ 1 2 -wi+ 1 2 -λ i+ 1 2 ,
(2.15) where r i+ 1 2 > 0 is such that the ball B( wi+ 1 2 , r i+ 1 2 ) centered on wi+ 1 2 with radius r i+ 1 2 is entirely contained within Ω. 2. entropy preserving: w n+1 i verifies the discrete entropy inequality (1.6). The proof of this statement is the purpose of the two next sections. To conclude this section, let us comment the assumptions imposed in the above result. First, by enforcing λ i+ 1 2 to be large enough is a very natural assumption. Formally, the condition λ i+ 1 2 > λ can be found in [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]. Indeed, as long as the eigenvalues (λ (w)) 1≤ ≤N of the matrix ∇ w f (w) + A(w) are known, the assumption stated on λ i+ 1 2 reads

λ i+ 1 2 ≥ max 1≤ ≤N |λ (w n i )|, |λ (w n i+1 )|) .
As a consequence, this parameter cannot be as small as we want and that is coherent with the CFL-like restriction (2.5). Next, we emphasize that the assumption (2.11) only concerns systems (1.1) in non-conservative form [START_REF] Castro | Why many theories of shock waves are necessary: convergence error in formally path-consistent schemes[END_REF][START_REF] Castro | A two-layer finite volume model for flows through channels with irregular geometry: Computation of[END_REF]16,17,[START_REF] Parés | Numerical methods for nonconservative hyperbolic systems: a theoretical framework[END_REF]. In addition, let us mention that the entropies issuing from problems of physical interest admit, in general, a positive definite Hessian matrix.

3. Artificial viscosity and robustness. In this section, we show that the artificial viscosity may enforce the robustness of the numerical scheme.

Lemma 3.1. Let (w n i ) i∈Z be a given sequence in Ω. Let w n+1 i be given by the viscous scheme (2.6). Let ∆t be restricted according to the CFL-like condition (2.5). As soon as γ ≥ 0 is large enough according to (2.15), then w n+1 i ∈ Ω for all i ∈ Z. Proof. Since Ω is convex, in view of (2.10) and thanks to the CFL condition (2.5), it is sufficient to prove that w

L i+ 1 2 , w R i+ 1 2
∈ Ω whenever w n i , w n i+1 ∈ Ω. Now, clearly the intermediate state wi+ 1 2 , given by (2.4), stays in Ω. Moreover, since Ω is an open set, then there exists r i+ 1 2 > 0 such that w ∈ Ω for every w satisfying w -wi+ 1 2 < r i+ 1 2 . Therefore, it is a sufficient condition to have

w L i+ 1 2 -wi+ 1 2 = λ i+ 1 2 λ i+ 1 2 + γ w L i+ 1 2 -wi+ 1 2 < r i+ 1 2 , w R i+ 1 2 -wi+ 1 2 = λ i+ 1 2 λ i+ 1 2 + γ w R i+ 1 2 -wi+ 1 2 < r i+ 1 2 ,
where w L i+ to belong to Ω and the proof is achieved.

Let us emphasize that the admissible states, here, belong to an open set. According to some particular physics, Ω is not necessarily open (for instance, see the shallow-water model [START_REF] Berthon | Improvement of the hydrostatic reconstruction scheme to get fully discrete entropy inequalities[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Morales | Reliability of first order numerical schemes for solving shallow water system over abrupt topography[END_REF]). The transition from states in Ω to states in ∂Ω, in general, needs a specific attention and it is not considered in the present work.

In addition to this general result, a constructive method can be proposed in the particular case where Ω is of the form:

Ω = w ∈ R d , p j (w) > 0 , j = 1 , • • • , L , (3.1)
where (p j ) j=1 ,••• ,L are continuous functions. According to the convexity of Ω, we recall that a sufficient condition to have robustness is to enforce the quantities (2.8), (2.9) to belong to Ω, which means that p j ( w L, * i+1/2 ) > 0 and

p j ( w R, * i+1/2 ) > 0 for j = 1 , • • • , L. Note that w n i , w n i+1 ∈ Ω ⇒ p j ( w i+1/2 ) > 0 for j = 1 , • • • , L, and lim γ→+∞ w L, * i+1/2 = lim γ→+∞ w R, * i+1/2 = wi+1/2 .
According to the continuity of the functions p j , this implies the existence of a γ large enough such that p j ( w L, * i+1/2 ) > 0 and p j ( w R, * i+1/2 ) > 0 for all j = 1, • • • L, which ensures the robustness property and is consistent with the previous Lemma. Of course, this limit value of γ may depend on the functions (p j ) j=1,••• ,L involved in the definition of Ω. However, an explicit treshold can be exhibited under some additional regularity hypothesis. More precisely, let assume that the functions p j admit local Taylor expansions at order N j ∈ N * such that:

p j (x + h) -   p j (x) + Nj k=1 1 k! dp k j,x (h)   ≤ M j (N j + 1)! h Nj +1 , ∀x ∈ R d ,
with M j > 0 and where dp k j,x (h) is the k-th order differential of p j at point x. Note that, as will be illustrated later, all the models considered in this paper enter this frame.

In what follows we will drop the subscript "j" for simplification purposes. Also, when no confusion is possible, we will drop the subscript "i + 1/2" and rewrite (2.8), (2.9) as:

w * = w + ν δw , (3.2)
where ν = λ λ + γ and ( w * , δw) being either equal to ( w L, * i+1/2 , w L, * i+1/2 -wi+1/2 ) or ( w R, * i+1/2 , w R, * i+1/2 -wi+1/2 ). Then, according to (3.2), we have:

p( w * ) -p( w) + N k=1 1 k! ν k dp k w(δw) ≤ M (N + 1)! |ν| N +1 δw N +1 , (3.3)
from which we deduce a sufficient condition to have p( w * ) > 0:

p( w) + N k=1 1 k! ν k dp k w(δw) - M (N + 1)! |ν| N +1 δw N +1 > 0 . (3.4)
Recalling that λ ≥ 0 and γ ≥ 0, we set µ = λ + γ > 0. Then, multiplying (3.4) by µ N +1 , we can express the previous condition through the positivity of a polynomial of order N + 1 in µ:

Q(µ) := µQ 1 (µ) - M (N + 1)! |λ| N +1 δw N +1 > 0 , (3.5)
where

Q 1 (µ) = µ N p( w) + N k=1 1 k! λ k µ N -k dp k w(δw).
(3.6)

Note that the polynomial (3.5) only depends on w and δw, which are explicit interface quantities independent from γ, according to (3.2). From this, we remark that the dominant term of Q is µ N +1 p( w) ≥ 0, so that the study of the roots of Q will allow to identify a real value µ + such that Q(µ) > 0 for all µ > µ + , or in an equivalent way, for all γ ≥ 0 satisfying:

γ > µ + -λ . (3.7)
We have thus established the following result: Lemma 3.2. Assume that Ω is under the form (3.1), and consider (w n i ) i∈Z a given sequence in Ω. Let w n+1 i be given by the viscous scheme (2.6). Let ∆t be restricted according to the CFL-like condition (2.5). Then, if the viscous constant γ ≥ 0 satisfies the condition:

γ > max i∈Z max j=1,••• ,L µ + j -λ i+1/2 ,
where µ + j is the largest root of the polynomial Q defined in (3.5), we have w n+1 i

∈ Ω for all i ∈ Z.

Remark 1. We propose here two examples based on the models implied in Section 5 dedicated to numerical experiments. Example 1: Shallow Water equations. We have Ω = w ∈ R d , p(w) > 0 with p(w) = h and d = 2, 3. We have (3.3) with N = 1 and M = 0:

p( w * ) = p( w) + νδw.∇p = h + νδh ,
Since M = 0, the sufficient condition (3.5) reduces to Q 1 (µ) > 0, that is:

(λ + γ) h + λδh > 0 ,
from which we deduce the condition:

γ > -λ h * h ,
with h * = δh + h. One may remark that if the original scheme is robust, i.e. the intermediate states h L, * i+1/2 , h R, * i+1/2 given by (2.2), (2.3) are positive, then h * > 0 and there is no condition on γ. Example 2: Euler equations.

We have Ω = w = t (ρ, ρu, ρE) ∈ R 3 , p 1 (w) > 0 , p 2 (w) > 0 with p 1 (w) = ρ and

p 2 (w) = ρ 2 E - (ρu) 2 2
, the second constraint corresponding to the positivity of the internal energy. We refer to the previous example for the condition associated with p 1 . Regarding the second condition, noting that p 2 is a second order polynomial in the variables ρ, ρu, ρE, we have (3.3) with N = 2 and M = 0. Considering (3.5), we are left with the study of the second order polynomial:

Q 1 (µ) = µ 2 p( w) + µλ δw.∇p( w) + 1 2 λ 2 dp 2 w(δw) .
According to the concavity of p, the quantity dp 2 w(δw) is negative (equal to -(δρu) 2 ), so that the discriminant ∆ = λ 2 (δw.∇p( w))

2 + 2p( w)(δρu) 2 is non-negative. As a result Q 1 admits the following positive root:

µ + = -λ δw.∇p( w) + √ ∆ 2p( w)
,

and the resulting condition on γ follows from (3.7).

Artificial viscosity and discrete entropy inequalities.

In order to get the discrete entropy inequality (1.6), we apply for an entropy preserving sufficient condition established by Harten, Lax and van Leer in [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]. For the sake of completeness, we recall this result and the reader is referred to [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] for the proof.

Lemma 4.1 (Harten, Lax and van Leer [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]). Let w n i and w n i+1 be given in Ω. Assume that the approximate Riemann solver satisfies

1 ∆x ∆x/2 -∆x/2 η WR x ∆t ; w n i , w n i+1 dx ≤ 1 2 η(w n i ) + η(w n i+1 ) - ∆t ∆x G(w n i+1 ) -G(w n i ) , (4.1)
for a given entropy pair (η, G). Then, the Godunov-type scheme (2.6) satisfies a discrete entropy inequality (1.6) where the numerical entropy flux function reads as follows:

G ∆ (w n i , w n i+1 ) = G(w n i+1 ) - ∆x 2∆t η(w n i+1 ) + 1 ∆x ∆x/2 0 η WR x ∆t ; w n i , w n i+1 dx, = G(w n i+1 ) + λ i+ 1 2 η(w R i+ 1 2 ) + (λ i+ 1 2 + γ)(η( wi+ 1 2 ) -η(w n i+1 )).
Now, equipped with such a result, the expected discrete entropy inequality (1.6) is obtained as soon as the inequality (4.1) is proved to be satisfied. Since the approximate Riemann solver is given by (2.1), we easily get

1 ∆x ∆x/2 -∆x/2 η WR x ∆t ; w n i , w n i+1 dx = 1 2 -(λ i+ 1 2 + γ) ∆t ∆x η(w n i ) + η(w n i+1 ) + 2γ ∆t ∆x η( wi+ 1 2 ) + λ i+ 1 2 ∆t ∆x η(w L i+ 1 2 ) + η(w R i+ 1 2
) , so that the inequality (4.1) recasts as follows:

E 0 i+ 1 2 + γD i+ 1 2 ≤ 0, (4.2)
where the quantities E 0 i+ 1 2 and D i+ 1 2 are defined by (2.12) and (2.13). It is worth noticing that E 0 i+ 1 2 is nothing but an entropy dissipation rate associated with the original scheme (1.4) while E 0 i+ 1 2 + γD i+ 1 2 stands for the entropy dissipation rate of the viscous scheme (1.7). Since η is a convex function, we immediately get D i+ 1 2 ≤ 0, which coincides with an entropy dissipation rate associated with the artificial viscosity. Now, we remark that E 0 i+ 1 2 is non-positive for an entropy preserving scheme (1.4). Then, the expected inequality (4.1) is immediately satisfied. But, as soon as the original scheme (1.4) is entropy violating, we have E 0 i+ 1 2 > 0. However, we notice that neither E 0 i+ 1 2 nor D i+ 1 2 depends on γ. As a consequence, the objective now is to fix γ ≥ 0 large enough such that the required inequality (4.2) holds true. In fact, a particular attention must be paid on a possible blow-up of γ. Indeed, after (2.14), γ may eventually go to infinity as D i+ 1 2 goes to zero. In the next statement, we establish that γ remains bounded. Lemma 4.2. Let (w n i ) i∈Z be given in Ω. Let w n+1 i be given by the viscous scheme (2.6), or equivalently (1.7). Assume that the discretization of the non-conservation product A(w)∂ x w satisfies (2.11). Assume λ i+ 1 2 > λ . Let ∆t be restricted according to the CFL-like restriction (2.5). Moreover, let the Hessian matrix ∇ 2 w η(w) be positive definite. Then, γ ≥ 0, given by (2.14), is bounded and the discrete entropy inequality (1.6) is satisfied.

Proof. After Lemma 4.1, the expected discrete entropy inequality (1.6) holds as soon as (4.2) is verified. As a consequence, the proof is completed when the quantity

-E 0 i+ 1 2 /D i+ 1
2 is proved to be bounded for E 0

i+ 1 2 ≥ 0.
∂ t W ex + ∂ x f (W ex )+A(W ex )∂ x W ex = 0.
The integration of the above relation over (-λ i+ 1 2 ∆t, λ i+ 1 2 ∆t) × (0, ∆t), with λ i+ 1 2 large enough and according (1.2), gives

1 2λ i+ 1 2 ∆t 2 ∆t 0 λ i+ 1 2 ∆t -λ i+ 1 2 ∆t W ex x ∆t ; w n i , w n i+1 dx dt = 1 2 (w n i + w n i+1 ) - 1 2λ i+ 1 2 f (w n i+1 ) -f (w n i+1 ) - 1 2λ i+ 1 2 ∆t ∆t 0 λ i+ 1 2 ∆t -λ i+ 1 2 ∆t A W ex x t ; w n i , w n i+1 ∂ x W ex x t ; w n i , w n i+1 dx dt.
Arguing (2.11), we obtain

w HLL i+ 1 2 = 1 2 (w n i + w n i+1 ) - 1 2λ i+ 1 2 f (w n i+1 ) -f (w n i+1 ) - 1 4λ i+ 1 2 A L ∆ (w n i , w n i+1 ) + A R ∆ (w n i , w n i+1 ) • (w n i+1 -w n i ).
Now, let us notice that

w HLL i+ 1 2 = 1 2 w L i+ 1 2 + w R i+ 1 2 , (4.4)
where w L,R i+ 1 2 are given by (2.2) and (2.3).

Moreover, since W ex , is defined as an exact solution of (1.1), necessarily it satisfies the entropy inequality (1.3) so that we have

∂ t η(W ex ) + ∂ x G(W ex ) ≤ 0.
Once again, integrating the above inequality over (-λ i+ 1 2 ∆t, λ i+ 1 2 ∆t) × (0, ∆t), because of (1.2), we get

1 2λ i+ 1 2 ∆t λ i+ 1 2 ∆t -λ i+ 1 2 ∆t η W ex x t ; w n i , w n i+1 dx ≤ 1 2 η(w n i ) + η(w n i+1 ) - 1 2λ i+ 1 2 G(w n i+1 ) -G(w n i ) .
Arguing the well-known Jensen's inequality, from (4.3) the above inequality reads

η(w HLL i+ 1 2 ) ≤ 1 2 (w n i + w n i+1 ) - 1 2λ i+ 1 2 G(w n i+1 ) -G(w n i ) .
Next, with E 0 i+ 1 2 defined by (2.12), we obtain the following estimation:

E 0 i+ 1 2 ≤ λ i+ 1 2 η(w L i+ 1 2 ) + η(w R i+ 1 2 ) -2η(w HLL i+ 1 2
) .

Since η is a convex function with ∇ 2 w η(w) a positive definite matrix, because of (4.4), we have near

w L i+ 1 2 = w R i+ 1 2 η(w L i+ 1 2 ) + η(w R i+ 1 2 ) -2η(w HLL i+ 1 2 ) = O w R i+ 1 2 -w L i+ 1 2 2 .

Involving (2.2) and (2.3), we get

O w R i+ 1 2 -w L i+ 1 2 2 = O w n i+1 -w n i 2 so that - E 0 i+ 1 2 D i+ 1 2 = λ i+ 1 2 O(1),
and the proof is achieved.

Numerical experiments.

In this section we present several examples to illustrate the efficiency of the proposed stability enforcing technique. First, we perform numerical approximations of the solutions of the isentropic gas dynamic model. We show the ability of the here designed technique to stabilize entropy violating schemes. Next, we present an application of the proposed technique to the numerical approximation of the one layer shallow-water system with flat bathymetry to correct the VF-Roe scheme and, finally, an application to the bi-layer shallow-water model.

Compressible Euler equations.

We consider the 1-D evolution of a compressible inviscid fluid of density ρ, velocity u, pressure p and specific total energy E and we denote by e = E -1 2 u 2 the specific internal energy. This evolution is governed by the 1-D Euler equations which read:

   ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p(ρ, e)) = 0, ∂ t (ρE) + ∂ x (ρEu + p(ρ, e)u) = 0. (5.1)
The pressure function is given by an ideal gas equation of state: p(ρ, e) = (Γ -1)ρe where the heat capacity ratio is Γ = 1.4. Denoting w = (ρ, ρu, ρE) T and f (w) = (ρu, ρu 2 + p(ρ, e), ρEu + p(ρ, e)u) T system (5.1) can be written:

∂ t w + ∂ x f (w) = 0.
(5.2)

This system admits a family of entropy/entropy flux functions given by

η(w) = ρϕ(p/ρ Γ ), G(w) = η(w)u,
where ϕ is a smooth function which satisfies the restrictions given in [START_REF] Tadmor | A minimum entropy principle in the gas dynamics equations[END_REF]. In the following, we will consider the entropies η i (w) = ρϕ i (p/ρ Γ ) given by the following three functions:

ϕ 1 (θ) = ln(θ), ϕ 2 (θ) = -θ 1/(Γ+1) , ϕ 3 (θ) = θ -2 Γ .
The adopted initial scheme (1.6) is given by the VF-Roe method [START_REF] Gallouët | Some recent finite volume schemes to compute Euler equations using real gas EOS[END_REF] which is known to be entropy violating for some test-cases. The corrected scheme (1.7) is obtained by defining the numerical viscosity parameter γ ensuring a discrete entropy inequality for the entropy η 1 (w) = ρϕ 1 (p/ρ Γ ). Consequently, an interesting issue is to know whether the entropy inequality for the other two entropies is satisfied or not.

A smooth solution.

We consider the smooth initial data given by:

u 0 (x) = 1, p 0 (x) = 1, ρ 0 (x) = 1 + exp -1 1-x 2 , |x| ≤ 1, 1, |x| ≥ 1.
for which the exact solution is smooth and given by: u(x, t) = 1, p(x, t) = 1, ρ(x, t) = ρ 0 (x-t). The computations have been run on the time interval [0, 0.5] with grids given by the following space steps : ∆x = 2/N with N = 100 × 2 j with j ∈ {0, .., 5}. For this smooth solution, the VF-Roe scheme is entropy satisfying. Hence the corrected scheme rigorously coincide with the original scheme, the artificial viscosity γ is always zero and the discrete entropy (in)equalities are satisfied for all the entropies η i (w) for i = 1, 2, 3. Figure 5.1 displays the density at the final time and the numerical convergence of the L 1 -norm of the error between the numerical solution and the exact solution. We observe a first order (with respect to ∆x) convergence which is what is expected for this smooth solution.

A Riemann problem.

We consider the Riemann problem given by the following initial data:

ρ 0 (x) = 1.0, x < 0, 0.25, x > 0. u 0 (x) = 0.0, p 0 (x) = 1.0, x < 0, 0.1, x > 0.
The computations have been run on the time interval [0, 0.2] with grids given by the following space steps : ∆x = 1/N with N = 100 × 2 j with j ∈ {0, .., 5}. For this Riemann-type solution, the VF-Roe scheme is still expected to be entropy satisfying.

In practice, we observe that, for all mesh sizes, the artificial viscosity is activated only in the first two time steps near the initial discontinuity. On the first time step, the maximum value of γ is γ max = 1.53 and the ratio γ/λ i+ 1 2 satisfies (γ/λ i+ 1 2 ) max = 0.75. On the second time step, we have respectively γ max = 0.48 and (γ/λ i+ 1 2 ) max = 0.26, and these values do not depend on the mesh size. For all the following time steps, we observe that γ = 0 and that the VF-Roe scheme is entropy satisfying not only for η 1 but also for η 2 and η 3 . Figure 5.2 displays the density, velocity and pressure at the final time and the numerical convergence of the L 1 -norm of the error for these variables. We observe a convergence with an order slightly larger than 0.5 (with respect to ∆x) which is the expected order of convergence for this solution with a contact discontinuity.

5.1.3.

A Riemann problem with a sonic rarefaction wave. We consider the Riemann problem given by the following initial data:

ρ 0 (x) = 1.0, x < 0, 0.25, x > 0. u 0 (x) = 0.0, p 0 (x) = 1.0, x < 0, 0.01, x > 0.
The computations have been run on the time interval [0, 0.25] with grids given by the following space steps : ∆x = 1/N with N = 100 × 2 j with j ∈ {0, .., 10}. For this Riemann-type solution, the VF-Roe scheme is known to produce a steady non entropic discontinuity at the sonic point of the rarefaction wave. This discontinuity generates a production of (mathematical) entropy at the discrete level which can be balanced thanks to the artificial viscosity term in the corrected scheme. Figure 5.3 displays the density, velocity and pressure at the final time and the numerical convergence of the L 1 -norm of the error for these variables. We observe that without the artificial viscosity correction, the VF-Roe scheme does not converge because of the steady non-entropic discontinuity at x = 0. When the scheme is corrected with the artificial viscosity, the amplitude of the non entropic discontinuity gets smaller as the mesh is refined enabling the approximate solution to converge towards the exact solution with an order slightly larger than 0.5.

Contrary to the previous test-case (with no sonic point in the rarefaction wave), we observe here that the artificial viscosity correction is regularly activated during the whole computation. However, the maximum values of γ and γ/λ i+ 1 2 are obtained at the second time step with γ max = 0.6 and (γ/λ i+ 1 2 ) max = 0.28 and these maximum values do not depend on the mesh size.

By construction, the numerical method is entropy-satisfying for η 1 . However, one can check that the approximate solution also satisfies discrete entropy inequalities for η 2 and η 3 . In Table 5.1 we present the value of the L 1 -norm of the positive part of the entropy budget, defined as the left hand side of inequality (1.6), obtained with both the original VF-Roe scheme and with the corrected scheme. We can see that, even for η 2 and η 3 , while the VF-Roe scheme produces entropy, the corrected scheme preserves a local discrete entropy inequality. 1.05e-06 < e-14 4.40e-07 < e-14 1.49e-06 < e-14 400 3.29e-07 < e-14 1.37e-07 < e-14 4.69e-07 < e-14 1600 3.43e-08 < e-14 1.43e-08 < e-14 4.901e-08 < e-14 6400 4.33e-09 < e-14 1.80e-09 < e-14 6.20e-09 < e-14 25600 1.53e-09 < e-14 6.38e-10 < e-14 2.19e-09 < e-14 102400 1.24e-09 < e-14 5.14e-10 < e-14 1.77e-09 < e-14 Table 5.1: L 1 -norm of the positive part of the entropy budget for the three mathematical entropies η 1 , η 2 , η 3 . Black: VF-Roe scheme. Red: VF-Roe scheme with artificial viscosity 5.2. Shallow-water equations. We now consider the well-known shallow water equations on a flat bottom:

       ∂ t h + ∂ x q = 0, ∂ t q + ∂ x q 2 h + 1 2 gh 2 = 0 , (5.3)
where h = h(x, t) and q = q(x, t) refer to the water height and discharge respectively.

Gathering the flow variables in the vector w = (h, q) T and setting f (w) = (q, q 2 /h + gh 2 /2) T , the system can be written on the form (5.2). We recall that the mechanical energy η(w) = 1 2 hu 2 + 1 2 gh 2 plays the role of a mathematical entropy associated to the entropy flux G(w) = η(w) + 1 2 gh 2 u.

Dam-break problem.

As previously, we consider the entropy-violating VF-Roe scheme to exhibit the regularizing virtues of the proposed methodology. We consider the following dam-break problem: h(x, 0) = 1.5, x < 12.5, 0.02, x > 12.5, u(x, 0) = 0.0 , on the computational domain [0, 25] with ∆x = 0.0025 m, with a final time T = 0.1 s.

As in Test-case 5.1.3, this initial condition leads to an entropy-violating discontinuity at the level of the rarefaction wave. This phenomenon can be observed in Figure 5.4, which displays the water height and velocity profiles associated with the original VF-Roe scheme and the corrected version of the scheme, together with the energy budget, computed on the basis of (1.6). We observe that without correction the VF-Roe scheme produces an entropyviolating solution characterized by a sudden growth of entropy at the sonic point x = 12.5m. As expected, the artificial viscosity is able to correct this failure.

Concerning convergence analysis, depicted in Figure 5.4 (bottom right), computations have been run on a series of regular meshes with space step ∆x = 1/N with N = 100 × 2 j with j ∈ {0, .., 8}. We clearly see that the VF-Roe scheme does not converge and that the introduction of artificial viscosity allows to recover a proper convergence for both flow variables.

The time evolution of the ratio (γ/λ i+ 1 2 ) max is proposed in Figure 5.5 for ∆x = 0.0025 m. It can be observed that it is especially during the first iterations that artificial viscosity must be introduced to ensure entropy stability, while a small correction is still required throughout the simulation. Nevertheless, the maximum value of this ratio is around 6.935e -2 and induces a minor increase of the computational time, according to the time step restriction (2.5). Note that this time behavior do not depend on the space resolution, which is in accordance with the previous test-cases. 5.2.2. Double rarefaction wave in the presence of vacuum. We now show that the proposed method is compatible with occurrence of dry states. To that purpose we use the test case employed in [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF] consisting of a Riemann problem which solution is composed of two rarefaction waves, with a dry zone occurring between the two waves. This test case stands for a relevant benchmark to highlight the difficulty for solvers to deal with vacuum (see for instance [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]), and [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] [START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF] for an adaptation of We observe a proper propagation of the rarefaction waves, implying a water height close to zero at the center of the domain as expected. Based on the space distribution of the local ratio γ/λ, we note that a subsequent amount of viscosity is needed to ensure stability, notably at the reconnection points with the initial constant states. Note that without correction the scheme is unstable and did not provide any exploitable result. Similar observations were made with the VFRoe and HLL schemes.

5.3. Bi-layer shallow-water system. We consider the homogeneous bi-layer 1-D shallow water system (see [START_REF] Castro | A Q-scheme for a class of systems of coupled conservation laws with source term. Application to a two-layer 1-D shallow water system[END_REF]):

                     ∂ t h 1 + ∂ x q 1 = 0, ∂ t q 1 + ∂ x q 2 1 h 1 + 1 2 gh 2 1 = -gh 1 ∂ x h 2 , ∂ t h 2 + ∂ x q 2 = 0, ∂ t q 2 + ∂ x q 2 2 h 2 + 1 2 gh 2 2 = - ρ 1 ρ 2 gh 2 ∂ x h 1 .
(5.4) Index 1 refers to the upper layer while index 2 refers to the lower layer. This system uses the following notation:

• h j = h j (x, t) ≥ 0 is the thickness of the j-th layer at the section of coordinate

x at time t. • q j = q j (x, t) is the discharge of the j-th layer at the section of coordinate x at time t and is related with the averaged velocity at each layer by the following relation: q j = u j h j , j = 1, 2. • g is the gravitational constant.

• ρ j refers to the constant density of the j-th layer with ρ 1 < ρ 2 .

The bottom is assumed to be flat. This system can be written in the form

∂ t w + ∂ x f (w) + A(w)∂ x w = 0, ( 5
w =     h 1 q 1 h 2 q 2     , f (w) =            q 1 q 2 1 h 1 + 1 2 gh 2 1 q 2 q 2 2 h 2 + 1 2 gh 2 2            , A(w) =     0 0 0 0 0 0 gh 1 0 0 0 0 0 grh 2 0 0 0     , with r = ρ 1 /ρ 2 .
An entropy-entropy flux pair for the bi-layer shallow water system is given by

η(w) = 2 j=1 ρ j h j u 2 j 2 + g h 2 j 2 + gρ 1 h 1 h 2 (5.6a) G(w) = 2 j=1 ρ j h j u 2 j 2 + gh 2 j + gh j b u j + ρ 1 gh 1 h 2 (u 1 + u 2 ).
(5.6b)

The wave speeds are the eigenvalues of A(w) = J(w) + A(w) = ∇ w f (w) + A(w), that is, the roots of the characteristic polynomial:

p(λ) = (λ 2 -2u 1 λ + u 2 1 -gh 1 )(λ 2 -2u 2 λ + u 2 2 -gh 2 ) -rgh 1 gh 2 .
(5.7)

When r ∼ = 1, first order approximations of the eigenvalues were given in [START_REF] Schijf | Theoretical considerations on the motion of salt and fresh water[END_REF]:

λ ± ext = u 1 h 1 + u 2 h 2 h 1 + h 2 ± g(h 1 + h 2 ), (5.8)
λ ± int = u 1 h 2 + u 2 h 1 h 1 + h 2 ± g h 1 h 2 h 1 + h 2 1 - (u 1 -u 2 ) 2 g (h 1 + h 2 ) , (5.9)
where g = (1 -r)g. The exact expression of the eigenvalues can be obtained by using Ferrari's method to find an analytical solution for quartic equations. For each eigenvalue λ, an associated eigenvector is given by:

R i =     1 λ µ λµ     , (5.10)
where:

µ = (λ -u 1 ) 2 gh 1 -1.
An important difficulty of system (5.4) is related to the loss of hyperbolicity: for r ∼ = 1 this situation occurs approximately when the following inequality is satisfied: .11) This loss of hyperbolicity is related to the appearance of shear instabilities that may lead, in real flows, to intense mixing of the two layers. While, in practice, this mixture partially dissipates the energy, in numerical experiments these interface disturbances may grow and overwhelm the solution. Obviously, a simple model based on two layer of immiscible fluids is not able to simulate the mixing processes due to the development of shear instabilities: a more complex multilayer model or a continuously stratified model would be required. In [START_REF] Castro-Díaz | Numerical treatment of the loss of hyperbolicity of the two-layer shallow-water system[END_REF] authors propose a simple strategy consisting on adding an extra friction term in order to get rid of the related instabilities and go beyond of them by reaching again the hyperbolic character.

F 2 r = (u 1 -u 2 ) 2 g (h 1 + h 2 ) > 1. ( 5
Here, we consider a standard path-conservative Roe solver for system (5.5) based on the family of straight segments

Φ(s; w L , w R ) = w L + s(w R -w L )
described in [START_REF] Castro | A Q-scheme for a class of systems of coupled conservation laws with source term. Application to a two-layer 1-D shallow water system[END_REF]. The numerical scheme reads as follows:

w n+1 i = w n i - ∆t ∆x f ∆ (w n i , w n i+1 ) -f ∆ (w n i-1 , w n i ) - ∆t 2∆x A L ∆ (w n i , w n i+1 ) • (w n i+1 -w n i ) + A R ∆ (w n i-1 , w n i ) • (w n i -w n i-1 ) , (5.12
) where the numerical flux is given by

f ∆ (w n i , w n i+1 ) = 1 2 f (w n i ) + f (w n i+1 ) - 1 2 |A n ∆ | • (w n i+1 -w n i ), (5.13)
where A n ∆ is a Roe matrix according to the choice of segment-path and is defined as follows

A ∆ =     0 1 0 0 gh ∆,1 -u 2 ∆,1 2u ∆,1 gh ∆,2 0 0 0 0 1 rgh ∆,1 0 gh ∆,2 -u 2 ∆,2 2u ∆,2     , (5.14)
where we have dropped the time dependency for simplicity. In (5.14), h ∆,j and u ∆,j , j = 1, 2 are the Roe averages and are defined as

h ∆,j = h i,j + h i+1,j 2 , u ∆,j = h i,j u i,j + h i+1,j u i+1,j h i,j + h i+1,j , j = 1, 2.
As usual, |A n ∆ | is the matrix that has the same eigenvectors than matrix A n ∆ and whose eigenvalues are the absolute value of those of A n ∆ . Finally,

A L ∆ (w n i , w n i+1 ) = A R ∆ (w n i , w n i+1 ) =     0 0 0 0 0 0 gh ∆,1 0 0 0 0 0 grh ∆,2 0 0 0     .
(5.15)

When A ∆ has complex eigenvalues, this matrix is no more a Roe linearization in the usual sense. Nevertheless, the numerical scheme (5.12)-(5.15) can still be applied by redefining A ∆ . Following [START_REF] Castro-Díaz | Numerical treatment of the loss of hyperbolicity of the two-layer shallow-water system[END_REF], we consider the real Jordan decomposition of A ∆ , that is,

A ∆ = K ∆ • L ∆ • K -1 ∆ ,
where L ∆ is a block diagonal matrix whose diagonal blocks are either the real eigenvalues or 2 × 2 blocks of the form: α β -β α (5.16) associated to every pair of conjugate complex eigenvalues α ± iβ. K ∆ is the real matrix corresponding to the change of basis. Now, |A ∆ | can be formally defined by setting |L ∆ | as the diagonal matrix obtained from the Jordan matrix by taking the absolute values of the real eigenvalues and by replacing the diagonal blocks (5.16) corresponding to a pair of conjugate complex eigenvalues by the diagonal block:

α 2 + β 2 0 0 α 2 + β 2 .
5.3.1. Riemann problem. Let us consider the following Riemann problem:

h 1 (x, 0) = 0.95 if x < 0 0.05 otherwise, h 2 (x, 0) = 1 -h 1 (x, 0), u 1 (x, 0) = u 2 (x, 0) = 0,
in the interval [-1, 1] with ∆x = 0.005 and CF L = 0.4. Open boundary conditions are set and r = 0.2. We approximate the solution of this Riemann problem with the usual Roe scheme and with Roe scheme with the artificial viscosity method described in this paper. Figure 5.7 shows the free surface and interface computed with both schemes at times t = 1 and 2 s. At it can be observed, a non-entropic shock is created by Roe scheme at x = 0, while it disappears when the artificial viscosity scheme is used. Similar behaviour is observed when plotting the velocities (see Figure 5.8).

As we have done previously, we show the time evolution of the ratio (γ/λ i+ 1 2 ) max in Figure 5.9 for ∆x = 0.005 m. Again, one can observe that it is during the first iterations that artificial viscosity must be introduced to ensure entropic stability, while a small correction is still required throughout the simulation. The maximum value of this ratio is around 0.0315. 5.3.2. Non-hyperbolic regime. This test was inspired by Test 1 introduced in [START_REF] Castro-Díaz | Numerical treatment of the loss of hyperbolicity of the two-layer shallow-water system[END_REF]. It consists on the evolution of a perturbation of shear two-layer fluid that is close to the unstable region (appearance of complex eigenvalues). The simulation is carried out on a flat channel described by the interval [-5, 5]. The initial condition is given by h 1 (x, 0) = 0.4 -0.1e -16x 2 , h 2 (x, 0) = 1.0 -h 1 (x, 0), u 1 (x, 0) = 0.15, u 2 (x, 0) = -0.15.

Free boundary conditions are imposed and the system is simulated during T = 10 s. The CFL parameter is set to 0.4 and ∆x = 0.01 and r = 0.99. Note that

F 2 r = (u 1 -u 2 ) 2 g (h 1 + h 2 ) ≈ 0.917
at every point at time t = 0 s. Recall that values of F 2 r close to 1 result in the appearance of complex eigenvalues. Due to the perturbation on the interface, the bi-layer system becomes unstable near the region where the perturbation is located. On the one hand, these instabilities propagate along the channel if the previous Roe solver is used and the numerical scheme drives towards a completely useless result. On the other hand, adding properly some artificial viscosity to control the entropy of the system, allows to perform a stable simulation. Of course, it is impossible to properly simulate such complex fluid with a simple bi-layer shallow water system, but the artificial viscosity allows to provide some stable simulation that could be seen as the best approximation that could be obtained with the bi-layer shallow-water system. Figure 5.10 shows the free surface and the interface evolution at t = 5 and 10 s obtained with Roe scheme previously described (left) and with the Roe scheme with the artificial viscosity technique to control the entropy (right). As expected, the initial perturbation grows in time when Roe scheme is applied, producing non-physical waves, while the Roe scheme combined with the artificial viscosity technique is able to produce stable simulations. Similar results can be observed for the velocities (see Figure 5.11). Figure 5.12 shows the evolution of the stability number F 2 r (5.11) at the same time steps. It is clear that the strong oscillations obtained with the Roe scheme are nonphysical, while the solution obtained with the Roe scheme with artificial viscosity is the best approximation that one could obtain with the bi-layer shallow-water system for such complex shear flows. In this case, the time evolution of the ratio (γ/λ i+ 1 2 ) max in Figure 5. [START_REF] Castro | A two-layer finite volume model for flows through channels with irregular geometry: Computation of[END_REF] shows large values and it is present in the whole simulation, although the maximum values occur during the initial iterations.

Fig. 5 . 1 :

 51 Fig. 5.1: Left: space variations of the density at the final time t = 0.5. Mesh size: 100 cells. Right: L 1 -error with respect to ∆x.

Fig. 5 . 2 :

 52 Fig. 5.2: Space variations at the final time t = 0.2 of the density (upper-left), velocity (upper-right) and pressure (lower-left). Mesh size: 100 cells. Lower-right: L 1 -error with respect to ∆x.

Fig. 5 . 3 :

 53 Fig. 5.3: Space variations at the final time t = 0.2 of the density (upper-left), velocity (upper-right) and pressure (lower-left). Mesh size: 100 cells. Lower-right: L 1 -error with respect to ∆x.

Fig. 5 . 4 :

 54 Fig. 5.4: Space variations at the final time T = 0.1s of the water height (upperleft), energy budget (upper-right) and velocity (lower-left). Mesh size: 10 000 cells. Lower-right: L 1 -error with respect to ∆x.

Fig. 5 . 5 :

 55 Fig. 5.5: Time evolution of the maximum ratio γ/λ.

Fig. 5 . 6 :

 56 Fig. 5.6: Space variations at the final time T = 0.125s of the velocity height (upperleft), energy budget (upper-right) and water height (lower-left), together with the ratio γ/λ (lower-right) Mesh size: 10 000 cells.

Fig. 5 . 7 :

 57 Fig. 5.7: Test 1: Free surface and interface: comparison between Roe scheme (blue and green lines) and Roe with artificial viscosity (red and brown lines).

Fig. 5 . 8 :

 58 Fig. 5.8: Test 1: Velocities at each layers: comparison between Roe scheme (blue and green lines) and Roe with artificial viscosity (red and brown lines).

Fig. 5 . 9 :

 59 Fig. 5.9: Time evolution of the maximum ratio γ/λ.

Scientific Computing, 80:924-956, 2019.

Acknowledgements. Manuel J. Castro acknowledges financial support from the Spanish Government and FEDER through the coordinated Research project RTI2018-096064-B-C21 and the Andalusian Government Research projects UMA18-FEDERJA-161 and P18-RT-3163. Arnaud Duran acknowledges financial support from the French National Research Agency project NABUCO, grant ANR-17-CE40-0025 and from the French National program INSU-CNRS (Institut National des Sciences de l´Univers -Centre National de la Recherche Scientifique) program LEFE-MANU (Les Enveloppes Fluides et Environnement -Méthodes Mathématiques et Numériques), project DWAVE. Tomás Morales acknowledges financial support from the Spanish Government and FEDER through the coordinated Research project RTI2018-096064-B-C22.

Data availability statement. The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Declarations. The authors have no relevant financial or non-financial interests to disclose