L Boittin

F Bouchut

M.-O Bristeau

A Mangeney

J Sainte-Marie

F Souillé

The Navier-Stokes system with temperature and salinity for free surface flows -Numerical scheme and validation

Keywords: Navier-Stokes equations, free surface flows, variable density flows, layeraveraged formulation, finite volume scheme

. These systems model hydrostatic free surface flows with density variations. We show that the finite volume scheme presented is well balanced with regards to the steady state of the lake at rest and preserves the positivity of the water height. A maximum principle on the density is also proved as well as a discrete entropy inequality in the case of the Euler system with variable density. Some numerical validations are finally shown with comparisons to 3D analytical solutions and experiments.

Introduction

In this paper we present a numerical scheme for the 3D incompressible Navier-Stokes-Fourier system with variable density and free surface, as well as numerical test cases. This model describes variable density flows with free surface, the density variations coming from differences in temperature and/or salinity. The model is presented in the companion paper [START_REF] Boittin | Low-Mach type approximation of the Navier-Stokes system with temperature and salinity for free surface flows[END_REF], in which a layer-averaged formulation is also given. The layer-averaged formulation suppresses the need for moving meshes [START_REF] Decoene | Sigma transformation and ALE formulation for threedimensional free surface flows[END_REF], [START_REF] Donea | Encyclopedia of Computational Mechanics[END_REF]. It allows to perform 3D simulations with a 2D fixed mesh.

Variable density flows are frequently studied by oceanographers. Different systems of coordinates exist, among which terrain-following coordinates and isopycnal coordinates. For a discussion of the advantages and disadvantages of the various coordinates frequently used in ocean models, the reader is referred to [START_REF] Griffies | Developments in ocean climate modelling[END_REF] and [START_REF] Song | Parametric vertical coordinate formulation for multiscale, Boussinesq, and non-Boussinesq ocean modeling[END_REF]. The layer-averaged model presented here is not a terrain-following coordinate model. Though the layer thicknesses are defined as fractions of the total water height, it is not a σ-coordinate system. The model also differs from isopycnal coordinate models because in the layer-averaged formulation, the layers exchange mass between themselves, which means that the internal layer boundaries are actually not physical.

For the Euler part of the Navier-Stokes-Fourier system with variable density, a finitevolume formalism is adopted. The hydrostatic reconstruction technique is used [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for Shallow Water flows[END_REF]. Therefore, the topography is accurately represented and the scheme is wellbalanced. Yet the discretization of the nonconservative pressure terms demands special care. For the viscosity terms, we use finite elements as in [START_REF] Allgeyer | Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface[END_REF].

In [START_REF] Audusse | Approximation of the hydrostatic Navier-Stokes system for density stratified flows by a multilayer model. Kinetic interpretation and numerical validation[END_REF]), a similar model was studied and simulated. The scheme presented in [START_REF] Audusse | Approximation of the hydrostatic Navier-Stokes system for density stratified flows by a multilayer model. Kinetic interpretation and numerical validation[END_REF]) was a 2D (x -z) scheme and relied on a kinetic interpretation. In the present work we present a fully 3D scheme which is more flexible in a certain sense, because the kinetic flux is only one of the possible choices for the numerical flux. With any flux consistent with the semi-discrete in time Euler system, the resulting scheme is well-balanced and preserves the nonnegativity of the water depth. A maximum principle on the density is satisfied. In order to prove an in-cell entropy inequality, we adopt a kinetic flux, already used in the context of the Shallow Water equations in [START_REF] Perthame | A kinetic scheme for the Saint-Venant system with a source term[END_REF], [START_REF] Audusse | Approximation of the hydrostatic Navier-Stokes system for density stratified flows by a multilayer model. Kinetic interpretation and numerical validation[END_REF]. The entropy inequality is satisfied for a constant topography and includes third-order rest terms. Moreover, the unknowns in [START_REF] Audusse | Approximation of the hydrostatic Navier-Stokes system for density stratified flows by a multilayer model. Kinetic interpretation and numerical validation[END_REF] were not the same, which resulted in a complicated numerical scheme -nonlinear systems were solved at each step and a Newton fixed-point method was used. The present scheme is simpler and does not involve nonlinear systems. Finally, the present scheme is more stable than the scheme in [START_REF] Audusse | Approximation of the hydrostatic Navier-Stokes system for density stratified flows by a multilayer model. Kinetic interpretation and numerical validation[END_REF], the CFL condition of which could actually degenerate and give a time step equal to zero. With the proposed scheme, the computational cost of the simulation of a non-Boussinesq flow is not greater than that of the simulation of a Boussinesq flow.

The proposed numerical scheme is validated on three test cases simulated with the Freshkiss3d code Freshkiss3d (2020). Each of these test cases allows to validate an aspect of the numerical scheme: wet/dry interface treatment, buoyancy terms, diffusion effects, second order extensions (space and time). . . A real application of hydrodynamics in a river with chlorides entries will be presented in a forthcoming paper (the research report is available [START_REF] Souillé | Remontées de chlorures dans la Vilaine[END_REF], in French).

The first test is a convergence test towards an analytical solution [START_REF] Bristeau | Various analytical solutions for the incompressible Euler and Navier-Stokes systems with free surface[END_REF] for the Euler system with variable density. In the second test, a lock exchange simulation is performed and the results are compared with experimental data available from the literature [START_REF] Adduce | Gravity Currents Produced by Lock Exchanges: Experiments and Simulations with a Two-Layer Shallow-Water Model with Entrainment[END_REF]. Finally, in two diffusion cases, the differences between the Navier-Stokes-Fourier with variable density and Boussinesq models are evidenced.

The paper is organized as follows. In section 2, the layer-averaged Navier-Stokes-Fourier system and the Euler with variable density introduced in (Boittin et al., 2020) are recalled. A numerical scheme for the layer-averaged Euler model with variable density is presented in Section 3, its properties are studied. An extension of this scheme for the layer-averaged Navier-Stokes-Fourier model with variable density is presented in section 4. The numerical test cases are presented in section 5.

2 The layer-averaged models minder_models From now on, we just call the Euler system, the Euler system with variable density. We briefly recall the features of the multilayer models presented in [START_REF] Boittin | Low-Mach type approximation of the Navier-Stokes system with temperature and salinity for free surface flows[END_REF]) and studied here. The multilayer Navier-Stokes-Fourier model is a layer-averaged version of the incompressible, hydrostatic Navier-Stokes-Fourier system

∇ • U = - ρ (T ) ρ 2 c p ∇ • (λ∇T ) + µ|∇ x,y u| 2 + µ ∂u ∂z 2 , (1) eq
:NSF_1_p2 The temperature T is linked to the density ρ via the state equation T = T (ρ). The viscosity coefficient is denoted by µ, the heat conductivity by λ and the specific heat capacity at constant pressure by c p .

∂ρ ∂t + ∇ • (ρU) = 0, (2)
The energy balance for model ( 1)-( 3) is We consider a free surface flow, therefore we assume z b (x, y) ≤ z ≤ η(t, x, y) := h(t, x, y) + z b (x, y), with z b (x, y) the bottom elevation and h(t, x, y) the water depth, see Fig. 1.

∂ ∂t ρ |u| 2 2 + ρe + ∇ x,y • u ρ |u| 2 2 + η z ρgdz 1 + ρe -µ∇ |u| 2 2 = ∇.(λ∇T ),
Let n b and n s be the unit outward normals at the bottom and at the free surface respectively defined by

n b = 1 1 + |∇ x,y z b | 2 ∇ x,y z b -1
, and n s = 1 1 + |∇ x,y η| 2 -∇ x,y η 1 .

On the bottom we prescribe an impermeability condition U.n b = 0, (4) eq:bottom_p2

and a friction condition given e.g. by a Navier law

µ 1 + |∇ x,y z b | 2 ∂u ∂n b = -κu, (5)
?eq:fric_p2?

with κ a Navier coefficient. On the free surface, the kinematic boundary condition

∂η ∂t + u(t, x, y, η) • ∇ x,y η -w(t, x, y, η) = 0, (6)
eq:free_surf_p is satisfied, along with the no stress condition

µ ∂ ũ ∂n s = 0, (7)
?eq:bound3ns_p

with ũ = (u, 0) T . On solid walls, we prescribe a slip condition

U • n = 0, (8)
?eq:slip?

coupled with an homogeneous Neumann boundary condition

µ ∂u ∂n = 0,
n being the outward normal to the considered wall. Boundary conditions for the temperature also have to be considered, we can choose either Neumann or Dirichlet conditions namely at the bottom where φ T | z b , φ T | η are two given temperature fluxes and T 0 b , T 0 s are two given temperatures. Since T = T (ρ), the boundary conditions for ρ naturally ensue from the boundary conditions for T .

λ ∂T ∂n b = φ T | z b , (9)
The system is completed with some initial conditions h(0, x, y) = h 0 (x, y), ρ(0, x, y) = ρ 0 (x, y), U(0, x, y, z) = U 0 (x, y, z).

The system (1)-(3) was derived from the compressible Navier-Stokes system in [START_REF] Boittin | Low-Mach type approximation of the Navier-Stokes system with temperature and salinity for free surface flows[END_REF]. More specifically, the derivation consisted in performing the incompressible limit. This model respects the second principle of thermodynamics (non-decreasing entropy).

2.1

The layer-averaged Navier-Stokes-Fourier model f_multilayer ? We consider a discretization of the fluid domain by layers (see Fig. 2) where the layer α contains the points of coordinates (x, y, z) with z ∈ L α (t, x, y) = (z α-1/2 , z α+1/2 ) and {z α+1/2 } α=1,...,N is defined by

z α+1/2 (t, x, y) = z b (x, y) + α j=1 h j (t, x, y), α ∈ [0, . . . , N ], h α (t, x, y) = z α+1/2 (t, x, y) -z α-1/2 (t, x, y) = l α h(t, x, y),
(13) eq:layer and N α=1 l α = 1. Using the notations (13), let us consider the space P N,t 0,h of piecewise constant functions defined by P N,t

0,h = 1 z∈Lα(t,x,y) (z), α ∈ {1, . . . , N } , (14)
?eq:P0_space?

Figure 2: Notations for the layerwise discretization.

fig: free_p2 where 1 z∈Lα(t,x,y) (z) is the characteristic function of the layer L α (t, x, y). Using this formalism, the projection of ρ, u, v and w on P N,t 0,h is a piecewise constant function defined by

X N (t, x, y, z, {z α }) = N α=1 1 z∈Lα(t,x,y) (z)X α (t, x, y), (15)
?eq:ulayer?

for X ∈ (ρ, u, v, w). When the quantities {ρ α (t, x, y)} α=1,...,N are known, if the function T = T (ρ) is known, it is possible to recover the temperature using the formula

T N (t, x, z, {z α }) = N α=1 1 z∈Lα(t,x,y) (z)T (ρ α (t, x, y)).
The layer-averaged Navier-Stokes-Fourier system introduced in (Boittin et al., 2020) reads

∂h ∂t + N α=1 ∇ x,y • (h α u α ) = - N α=1 ρ (T α ) ρ 2 α c p (S T,α -S µ,α ), ( 16
) eq:massesvml ∂ρ α h α ∂t + ∇ x,y • (ρ α h α u α ) = ρ α+1/2 G α+1/2 -ρ α-1/2 G α-1/2 , α = 1, ..., N, (17)
?eq:massesvm

∂ρ α h α u α ∂t + ∇ x,y • (ρ α h α u α ⊗ u α ) + ∇ x,y h α p α = p α+1/2 ∇ x,y z α+1/2 -p α-1/2 ∇ x,y z α-1/2 +u α+1/2 ρ α+1/2 G α+1/2 -u α-1/2 ρ α-1/2 G α-1/2 + ∇ x,y • (µh α ∇ x,y u α ) +Γ α+1/2 (u α+1 -u α ) -Γ α-1/2 (u α -u α-1 ) -κ α u α , α = 1, ..., N, ( 18
) eq:mvtsvml_n with G α+1/2 = - N j=1 α p=1 l p -1 j≤α ∇ x,y • (h j u j ) + α j=1 ρ (T j ) ρ 2 j c p (S T,j -S µ,j ), ( 19
) ?eq:Qalphabis_ κ α = κ if α = 1 0 if α = 1 , S T,α = λ∇ x,y • (h α ∇ x,y T α ) + 2λ α+1/2 T α+1 -T α h α+1 + h α -2λ α-1/2 T α -T α-1 h α + h α-1 , (20) S_T_alpha λ α+1/2 = λ for α = 1, . . . , N -1, T α = T (ρ α ).
For α = 0, 2λ α+1/2 T α+1 -Tα

h α+1 +hα = φ T | z b if the Neumann boundary condition (9) is chosen, or h 0 = h 1 , T 0 = T 0 b if the Dirichlet boundary condition (10) is chosen. Likewise, for α = N , 2λ α+1/2 T α+1 -Tα h α+1 +hα = φ T | η with the boundary condition (11), or h N +1 = h N , T N +1 = T 0 s
with the boundary condition [START_REF] Donea | Encyclopedia of Computational Mechanics[END_REF]. The dissipation term due to the viscous effects is

S µ,α = -h α µ|∇ x,y u α | 2 -Γ α+1/2 |u α+1 -u α | 2 2 -Γ α-1/2 |u α -u α-1 | 2 2 -κ α |u α | 2 , ( 21
) S_mu_alpha with Γ α+1/2 = 2µ α+1/2 h α+1 + h α , ( 22
) Gamma_alpha_pl µ α+1/2 =      0 if α = 0 µ if α = 1, . . . , N -1 0 if α = N. ( 23
) ?mu_alpha_plus
The term |∇ x,y u α | 2 actually denotes

|∇ x,y u α | 2 = (∇ x,y u α ) : (∇ x,y u α ) T = ∂u α ∂x 2 + ∂u α ∂y 2 + ∂v α ∂x 2 + ∂v α ∂y 2 .
The velocities u α+1/2 and the densities ρ α+1/2 at the interfaces are defined by

v α+1/2 = v α if G α+1/2 ≤ 0 v α+1 if G α+1/2 > 0 (24) ?eq:upwind_uT_ for v = u, ρ.
The pressure terms p α , p α+1/2 are given by

p α = g ρ α h α 2 + N j=α+1 ρ j h j and p α+1/2 = g N j=α+1 ρ j h j . ( 25
) ?eq:palpha1_p2
The pressure is hydrostatic. The terms G α+1/2 represent the mass exchanges between the layers. Notice that some of the viscous and diffusion terms have been simplified, see [START_REF] Boittin | Low-Mach type approximation of the Navier-Stokes system with temperature and salinity for free surface flows[END_REF]. Since the right-hand side of the total height conservation equation ( 16) is nonzero, we expect to observe dilatation and contraction due to the temperature diffusion and to the viscosity. We recall here the following result, obtained in [START_REF] Boittin | Low-Mach type approximation of the Navier-Stokes system with temperature and salinity for free surface flows[END_REF].

y_balance_ns ? Proposition 2.1 The system (16)-( 18) completed with the equation

∂ ∂t (ρ α h α e α ) + ∇ x,y • (ρ α h α u α e α ) = ρ α+1/2 e α+1/2 G α+1/2 -ρ α-1/2 e α-1/2 G α-1/2 + p α ρ (T α ) ρ 2 α c p (S T,α -S µ,α ) + S T,α -S µ,α
admits, for smooth solutions, the energy balance

∂ ∂t E α + ∇ x,y • (u α (E α + h α p α -µh α ∇ x,y u α )) +Γ α+1/2 |u α+1 | 2 -|u α | 2 2 -Γ α-1/2 |u α | 2 -|u α-1 | 2 2 = ρ α+1/2 u 2 α+1/2 2 + gρ α+1/2 z α+1/2 G α+1/2 + p α+1/2 G α+1/2 - ∂z α+1/2 ∂t -ρ α-1/2 u 2 α-1/2 2 + gρ α-1/2 z α-1/2 G α-1/2 -p α-1/2 G α-1/2 - ∂z α-1/2 ∂t - 1 2 ρ α+1/2 (u α+1/2 -u α ) 2 + gh α (ρ α+1/2 -ρ α ) G α+1/2 + 1 2 ρ α-1/2 (u α-1/2 -u α ) 2 -gh α (ρ α-1/2 -ρ α ) G α-1/2 + S T,α , (26) eq:energy_mcl_
with

E α = ρ α h α |u α | 2 2 + ρ α gh α z α 2 + e α . ( 27
) ?eq:energ_al_N
Note that in (26), we use the notation

u α ∇ x,y u α = u α ∂uα ∂x + v α ∂vα ∂x u α ∂uα ∂y + v α ∂vα ∂y = ∇ x,y |u 2 α | 2 .
The sum of Eqs. (26) for α = 1, . . . , N gives

∂ ∂t N α=1 E α + N α=1 ∇ x,y • u α (E α + h α p α ) = - N α=1 ρ α+1/2 |u α+1 -u α | 2 2 |G α+1/2 | + N α=1 S T,α - g 2 N α=1 h α (ρ α+1/2 -ρ α ) + h α+1 (ρ α+1/2 -ρ α+1 ) G α+1/2 .
The sum of S T,α over the layers gives

N α=1 S T,α = λ N α=1 ∇ x,y • (h α ∇ x,y T α ) -∇T | s • n s + ∇T | b • n b .
As explained in [START_REF] Boittin | Low-Mach type approximation of the Navier-Stokes system with temperature and salinity for free surface flows[END_REF], the terms on the last line of the right-hand side are third-order terms.

The layer-averaged Euler system

What we refer to hereafter as the layer-averaged Euler system is the system ( 16)-( 18) without viscosity and without diffusion terms, i.e.

∂h ∂t + N α=1 ∇ x,y • (h α u α ) = 0, ( 28
) eq:massesvml ∂ρ α h α ∂t + ∇ x,y • (ρ α h α u α ) = ρ α+1/2 G α+1/2 -ρ α-1/2 G α-1/2 , α = 1, ..., N, (29)
?eq:massesvm

∂ρ α h α u α ∂t + ∇ x,y • (ρ α h α u α ⊗ u α ) + ∇ x,y h α p α = p α+1/2 ∇ x,y z α+1/2 -p α-1/2 ∇ x,y z α-1/2 +u α+1/2 ρ α+1/2 G α+1/2 -u α-1/2 ρ α-1/2 G α-1/2 , α = 1, ..., N. ( 30
) eq:mvtsvml_p
The quantity G α+1/2 (resp. G α-1/2 ) corresponds to mass exchange accross the interface z α+1/2 (resp. z α-1/2 ) and G α+1/2 is defined here by

G α+1/2 = α j=1 ∂h j ∂t + ∇ x,y • (h j u j ) = - N j=1 α p=1 l p -1 j≤α ∇ x,y • (h j u j ), (
31) eq:Qalphabis_p for α = 1, . . . , N . Notice that the mass conservation for the layer α writes The energy balance verified by the system (28)-(30) is given in [START_REF] Boittin | Low-Mach type approximation of the Navier-Stokes system with temperature and salinity for free surface flows[END_REF]. It is very similar to the balance (26), obviously without the viscosity and diffusion terms. In the balance for the Euler system, the internal energy e α does not intervene. It is actually equal to 0 because there is no volume variation.

∂h α ∂t + ∇ x,y • (h α u α ) = G α+1/2 -G α-1/2 , (32) eq
3 Numerical scheme for the layer-averaged Euler system sec:euler_num In this section, a numerical scheme for the layer-averaged Euler system is designed and analyzed. It extends the work done by some of the authors in [START_REF] Allgeyer | Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface[END_REF][START_REF] Audusse | Approximation of the hydrostatic Navier-Stokes system for density stratified flows by a multilayer model. Kinetic interpretation and numerical validation[END_REF]. Before specifying the scheme for the Euler system, a common strategy for the time discretization of the Navier-Stokes-Fourier and Euler systems is presented. The discretization of the diffusion terms does not present any additional difficulty, however including these terms considerably lengthens the equations. This is why it seems preferable to explain the numerical scheme for the Euler system first. The advantages of the numerical scheme are the following

• it gives a 3D approximation of the Navier-Stokes-Fourier system, while only 2D situations were considered in [START_REF] Audusse | Approximation of the hydrostatic Navier-Stokes system for density stratified flows by a multilayer model. Kinetic interpretation and numerical validation[END_REF] • it can be implemented with any flux that is consistent with the homogeneous Saint-Venant system; the kinetic flux is used only for the discrete entropy property stated for a constant topography

• the scheme is endowed with strong stability properties (well-balanced, positivity of the water depth)

• convergence curves towards a 3D non-stationary analytical solution with wet-dry interfaces are obtained, see paragraph 5.1.

Strategy for the time discretization

The system ( 16)-( 18) has the form

∂U ∂t + ∇ x,y • F (U) = S p (U, z b ) + S e (U, ∂ t U, ∂ x U) + S v,f (U), ( 33
) eq:glo_p2
where the vector of unknowns is

U = (h, ρ 1 h 1 , . . . , ρ N h N , q
x,1 , . . . , q x,N , q y,1 , . . . , q y,N ) T , with q x,α = ρ α h α u α , q y,α = ρ α h α v α . We denote by and S v,f (U), representing respectively the mass and momentum exchanges and the viscous and friction effects. Notice that, as a consequence of the layer-averaged discretization, the system (33) is made of only 2d (x, y) partial differential equations with source terms. Hence, the spacial approximation of the considered PDEs is performed on a 2d planar mesh. We consider discrete times t n with t n+1 = t n + ∆t n . For the time discretisation of the layer-averaged Navier-Stokes-Fourier system (33) we adopt the following scheme

F (U) = (F x (U), F y (U)) T the
U n+1 = U -∆t n (∇ x,y • F (U) -S p (U, z b )) + ∆t n S n+1 e + ∆t n S n+l v,f , (34
) eq:glo_dis

where the integer l = 0, 1/2, 1 will be precised below. In (34) and wherever there is no ambiguity the superscript n has been omitted.

Semi-discrete (in time) scheme

From now on and until the end of this section, the system considered is the Euler system. Similarly to [START_REF] Allgeyer | Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface[END_REF], the semi-discrete in time scheme (34) with S n+l v,f = 0 writes

h n+1/2 α = h α -∆t n ∇ x,y (h α u α ) (35) eq:nscont_laye (ρ α h α ) n+1/2 = ρ α h α -∆t n ∇ x,y • (ρ α h α u α ) (36) eq:nsml22_dbis (ρ α h α u α ) n+1/2 = ρ α h α u α -∆t n ∇ x,y • (ρ α h α u α ⊗ u α ) + ∇ x,y h α p α -p α ∇ x,y h α + ρ α gh α ∇ x,y z α , (37)
eq:nsml22_d

h n+1 = h n+1/2 = N α=1 h n+1/2 α = h -∆t n N α=1 ∇ x,y • (h α u α ), ( 38
) eq:nsml11_d (ρ α h α ) n+1 = (ρ α h α ) n+1/2 + ∆t n ρ n+1 α+1/2 G α+1/2 -ρ n+1 α-1/2 G α-1/2 , ( 39
) eq:nsml33bis_d (ρ α h α u α ) n+1 = (ρ α h α u α ) n+1/2 + ∆t n u n+1 α+1/2 ρ n+1 α+1/2 G α+1/2 -u n+1 α-1/2 ρ n+1 α-1/2 G α-1/2 , ( 40
) eq:nsmlbis_d with G α+1/2 = - N j=1 α p=1 l p -1 j≤α ∇ x,y • (h j u j ).
(41) eq:nsf_6d

Since the relations

p α±1/2 = p α ∓ ρ α gh α 2 , z α = z α+1/2 + z α-1/2 2
hold, in (37) we often use the identity

p α+1/2 ∇ x,y z α+1/2 -p α-1/2 ∇ x,y z α-1/2 = p α ∇ x,y h α -ρ α gh α ∇ x,y z α . (
U n+1/2 = U -∆t n (∇ x,y • F (U) -S p (U, z b )) , (
U n+1 = U n+1/2 + ∆t n S n+1 e , (44)
?eq:glo_dis2?

this implicit step being precised in paragraph 3.4.3.

Eigenvalues for the advection and pressure part

In order to propose a finite volume discretisation for Eqs. ( 35)-(38), an estimation of the eigenvalues is necessary and given below. Without the exchange terms, the system (28)-(30) writes

∂h ∂t + N α=1 ∇ x,y • (h α u α ) = 0, ( 45
) eq:massesvml_s ∂(ρ α h α ) ∂t + ∇ x,y • (ρ α h α u α ) = 0, ( 46
) ?eq:massesvml1 ∂(ρ α h α u α ) ∂t + ∇ x,y • (ρ α h α u α ⊗ u α ) + p α ∇ x,y h α -ρ α gh α ∇ x,y z α = 0, ( 47
) eq:mvtsvml_se with p α ∇ x,y h α -ρ α gh α ∇ x,y z α = gh α N j=α+1 ∇ x,y (ρ j h j ) + 1 2 ∇ x,y (ρ α h α ) + ρ α gh α 2 ∇ x,y h α + ρ α gh α α-1 j=1 ∇ x,y h j + ∇ x,y z b .
This system is the continuous version of system (35)-(38). We rewrite the system (45)-(47) under the form

∂h ∂t + N j=1 (l j h∇ x,y • u j ) + N j=1
(l j u j ) • ∇ x,y h = 0, (48) eq:massesvml_s

∂ρ α ∂t + u α .∇ x,y ρ α = 0, (49)
?eq:massesvml1

∂u α ∂t + (u α • ∇ x,y )u α + 1 ρ α h α (p α ∇ x,y h α -ρ α gh α ∇ x,y z α ) = 0, ( 50
) eq:mvtsvml_se1
and the quasilinear form of the system (48)-(50) writes

∂ Ũ ∂t + A( Ũ)∇ x,y Ũ = s b ( Ũ), ( 51
) eq:quasi_lin with Ũ = (h, u 1 , . . . , u N , ρ 1 , . . . , ρ N ) T , and
A( Ũ) = A 1 ( Ũ) A 2 ( Ũ) A 3 ( Ũ) A 4 ( Ũ) , A 1 ( Ũ) =           N j=1 l j u j l 1 h . . . . . . . . . l N h p1 u 1 0 . . . . . . 0 p2 0 u 2 0 . . . 0 . . . 0 . . . u j . . . 0 . . . 0 . . . 0 . . . 0 pN 0 0 . . . 0 u N           , A 2 ( Ũ) =         0 . . . . . . 0 gh 2 1 2 0 . . . 0 gh 2 h 1 gh 2 2 2 0 0 . . . 0 gh 2 j 2 0 gh N h 1 0 0 gh 2 N 2         , A 3 ( Ũ) = 0 N ×N +1 , A 4 ( Ũ) = diag(u j ), pj = g ρ j ρ j l j + ρ j j-1 i=1 l i + j-1 i=1 ρ i l i .
For the sake of simplicity, the expression of the matrix A( Ũ) given below corresponds to the 1D case i.e. for v i = 0, i = 1, . . . , N . Notice that A 2 ( Ũ) and A 3 ( Ũ) T are rectangular matrices with N + 1 rows and N columns.

The following proposition holds, consisting in a version of the Cauchy's interlace theorem [START_REF] Hwang | Cauchy's interlace theorem for eigenvalues of hermitian matrices[END_REF] in the case of non symmetric matrix. Proposition 3.1 The system (51) is strictly hyperbolic for h > 0 and the eigenvalues of A( Ũ) belong to the interval (λ min , λ max ) with

λ min = min j {u j , v j } - max{ρ j } min{ρ j } gh, λ max = max{u j , v j } + max{ρ j } min{ρ j } gh. prop:hyper Proof of prop. 3.1 Since A( Ũ) is a block-matrix, its eigenvalues consist in the eigen- values of A 1 ( Ũ) completed with the set {u i } N i=1 .
Writing the characteristic polynomial of A 1 ( Ũ) under the form (development e.g. along the first row)

P A 1 = Π N i=1 (λ -u i ) λ - N j=1 l j u j - N i=1 l i hp i λ -u i , the eigenvalues of A 1 ( Ũ) satisfy λ - N j=1 l j u j = Q A 1 (λ),
:pol_carac with Q A 1 (λ) = N i=1 l i hp i λ-u i . For N = 3, the functions λ → Q A 1 (λ) and λ → λ -N
j=1 l j u j are depicted over Fig. 3 and it is easy to see that the four eigenvalues λ i exists with the interlacing

λ 1 < u 1 ≤ λ 2 ≤ . . . ≤ u 3 < λ 4 .

Moreover we have

Q A 1 (λ max ) = N i=1 l i hp i λ max -u i ≤ min{ρ j } max{ρ j } N i=1 l i hp i √ gh ≤ max{ρ j } min{ρ j } gh ≤ λ max - N j=1 l j u j ,
and likewise

Q A 1 (λ min ) = N i=1 l i hp i λ min -u i ≥ - min{ρ j } max{ρ j } N i=1 l i hp i √ gh ≥ - max{ρ j } min{ρ j } gh ≥ λ min - N j=1 l j u j , therefore the eigenvalues {λ i } N i=1 of A 1 ( Ũ) satisfy λ min ≤ λ i ≤ λ max , i = 1, . . . N,
proving the result.

Finite volume formalism for the Euler part

subsec:fv In this paragraph, we propose a space discretization for the model ( 35)-(40) completed with (41). We first recall the general formalism of finite volumes on unstructured meshes.

Let Ω denote the computational domain with boundary Γ, which we assume is polygonal.

Let T h be a triangulation of Ω for which the vertices are denoted by P i with S i the set of interior nodes and G i the set of boundary nodes. The dual cells C i are obtained by joining the centers of mass of the triangles surrounding each vertex P i . We use the following notations (see Fig. 4):

• K i , set of subscripts of nodes P j surrounding P i ,

• |C i |, area of C i , • Γ ij , boundary edge between the cells C i and C j , • L ij , length of Γ ij , • n ij , unit normal to Γ ij , outward to C i (n ji = -n ij ).
If P i is a node belonging to the boundary Γ, we join the centers of mass of the triangles adjacent to the boundary to the middle of the edge belonging to Γ (see Fig. 4) and we denote

• Γ i , the two edges of C i belonging to Γ,

• L i , length of Γ i (for sake of simplicity we assume in the following that L i = 0 if P i does not belong to Γ),

• n i , the unit outward normal defined by averaging the two adjacent normals. We define the piecewise constant functions U n (x, y) on cells C i corresponding to time

t n as U n (x, y) = U n i , for (x, y) ∈ C i , with U n i = (h n i , ρ n 1,i h n 1,i , . . . , ρ n N,i h n N,i , q n
x,1,i , . . . , q n x,N,i , q n y,1,i , . . . , q n y,N,i ) T i.e.

U n i ≈ 1 |C i | C i U(t n , x, y)dxdy.
For the topography, we choose a piecewise constant approximation under the form

z i ≈ 1 |C i | C i z b (x, y)dxdy.
We will also use the notation

U n α,i ≈ 1 |C i | C i U α (t n , x, y)dxdy, with U α defined by U α = (h α , ρ α h α , ρ α h α u α , ρ α h α v α ) T .
(52) eq:u_alpha

The horizontal fluxes and the pressure terms

ntal_discrete A finite volume scheme for solving the system (35)-( 38) is a formula of the form

U n+1/2 i = U i - j∈K i σ i,j F i,j -σ i F e,i + j∈K i σ i,j S p (U i , U j , z b,i , z b,j ), ( 53
) eq:upU0
where using the notations of ( 34)

j∈K i L i,j F i,j ≈ C i ∇ x,y • F (U)dxdy, (54) eq:flux_dis1
with

σ i,j = ∆t n L i,j |C i | , σ i = ∆t n L i |C i | .
Here we consider first-order explicit schemes where

F i,j =             F h i,j F ρ 1 h 1 i,j
. . .

F ρ N h N i,j F u 1 i,j
. . .

F u N i,j             . ( 55
) eq:flux_def and F h i,j = N α=1 F hα i,j , (56)
eq:flux and for the boundary nodes

F i,e =             F h i,e F ρ 1 h 1 i,e
. . .

F ρ N h N i,e F u 1 i,e
. . .

F u N i,e             . ( 57
) eq:fluxbis
The fluxes F hα i,j , F ραhα i,j

, F uα i,j appearing in expressions (55),( 56), (57) are numerical fluxes such that

F mα i,j = F mα (U α,i , U α,j , n i,j ), with m = h, ρh, u, α = 1, . . . , N .
Relation (53) tells how to compute the values U n+1/2 i knowing U i and discretized values z b,i of the topography. Following (54), the term F i,j in (53) denotes an interpolation of the normal component of the flux F (U).n i,j along the edge C i,j . The functions F (U i , U j , n i,j ) ∈ R 2N +1 are the numerical fluxes, see [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF].

Until now, the expression for the numerical fluxes is not detailed since any numerical fluxes (Rusanov, HLL,. . . ) can be used [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]. In paragraph 3.5 we define F(U i , n i,j ) using kinetic fluxes and we prove a discrete entropy inequality for the system. The computation of the value U i,e , which denotes a value outside C i (see Fig. 4-(b)), defined such that the boundary conditions are satisfied, and the definition of the boundary flux F (U i , U e,i , n i ) are described in [START_REF] Allgeyer | Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface[END_REF]. Notice that we assume a flat topography on the boundaries i.e. z b,i = z b,i,e .

For the discretization of the pressure source term S p (U, z b ), we adopt a strategy defined below.

The hydrostatic reconstruction technique

subsubsec:HR The hydrostatic reconstruction scheme (HR scheme for short) for the Saint-Venant system has been introduced in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for Shallow Water flows[END_REF] in the 1d case and described in 2d for unstructured meshes in [START_REF] Audusse | A well-balanced positivity preserving second-order scheme for Shallow Water flows on unstructured meshes[END_REF]. The HR in the context of the kinetic description for the Saint-Venant system has been studied in [START_REF] Audusse | Kinetic entropy inequality and hydrostatic reconstruction scheme for the Saint-Venant system[END_REF].

In order to take into account the topography variations and to preserve relevant equilibria, the HR leads to a modified version of (53) under the form

U n+1/2 i = U i - j∈K i σ i,j F * i,j -σ i F i,e + j∈K i σ i,j S * p,i,j , (58) eq:upU0_HR
where

F * i,j = F (U * i,j , U * j,i , n i,j ), ( 59
) ?eq:flux_HR? S * p,i,j = S p (U i , U * i,j , z b,i , z b,j , n i,j ) =         0 p * 1,i,j (h 1,i,j -h 1,i )n i,j -gρ 1,i h * 1,i,j (z 1,i,j -z 1,i )n i,j . . . p * α,i,j (h α,i,j -h α,i )n i,j -gρ α,i h * α,i,j (z α,i,j -z α,i )n i,j . . .         (60) eq:Sp_HR with z * b,i,j = max(z b,i , z b,j ), h * i,j = max(h i + z b,i -z * b,i,j , 0), U * i,j = (h * i,j , ρ 1,i l 1 h * i,j , . . . , ρ N,i l N h * i,j , ρ 1,i l 1 h * i,j u 1,i , . . . , ρ N,i l N h * i,j u N,i , . . .) T , z * α,i,j = z * b,i,j + lα 2 + α-1 j=1 l j h * i,j , z α,i,j = z α,j,i = z * α,i,j +z * α,j,i 2 , h α,i,j = h α,j,i = h * α,i,j +h * α,j,i 2 , h * α,i,j = h α,i +h * α,i,j 2 , p * α,i,j = p α,i +p * α,i,j 2 ,
(61) eq:state_HR and

z α = z b + l α 2 + α-1 j=1 l j h, p α = ρ α gh α 2 + N j=α+1 ρ j gh j .
Throughout this work, the * refers to the HR technique.

Remark 3.2 Since the quantity S(U, z b ) appearing in (33) contains non conservative terms, its integration over the cell C i is not straightforward and we have used the result proposed by Bouchut [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF], Proposition 5.3) to obtain the expression (60).

The vertical exchange terms

sec:exchanges We give the fully discrete expression of the step for the vertical exchanges, described by equations (39)-(40). The step for the vertical exchanges consists in

U n+1 i = U n+1/2 i + ∆t n G n+1 i , ( 62
) eq:upU1_HR with G n+1 i =                      0 ρ n+1 3/2,i G 3/2,i ρ n+1 5/2,i G 5/2,i -ρ n+1 3/2,i G 3/2,i . . . ρ n+1 N -1/2,i G N -1/2,i -ρ n+1 N -3/2,i G N -3/2,i -ρ n+1 N -1/2,i G N -1/2,i u n+1 3/2,i ρ n+1 3/2,i G 3/2,i u n+1 5/2,i ρ n+1 5/2,i G 5/2,i -u n+1 3/2,i ρ n+1 3/2,i G 3/2,i . . . u n+1 N -1/2,i ρ n+1 N -1/2,i G N -1/2,i -u n+1 N -3/2,i ρ n+1 N -3/2,i G N -3/2,i -u n+1 N -1/2,i ρ n+1 N -1/2,i G N -1/2,i                      .
The mass conservation for the layer α is governed by Eq. (31) and its discretization gives

h n+1 α,i = h n+1/2 α,i + ∆t(G α+1/2,i -G α-1/2,i ).
The previous equation multiplied by ρ n+1 α,i from equation (39) gives the relation

ρ n+1 α,i h n+1/2 α,i = ρ n+1/2 α,i h n+1/2 α,i + ∆t((ρ n+1 α+1/2,i -ρ n+1 α,i )G α+1/2,i -(ρ n+1 α-1/2,i -ρ n+1 α,i )G α-1/2,i
). Hence, relation (62) can be rewritten under the form

(H n+1/2 N,i + ∆tG N,i )ρ n+1 i = (ρ i h i ) n+1/2 , ( 63
) eq:matG where H n+1/2 N,i is a diagonal matrix of size N with coefficients (H N,i j ) 1≤j≤N = h n+1/2 j
and the matrix G N,i is given by

G N,i =                 |G 3/2,i | + h n+1/2 1,i - |G 3/2,i | + h n+1/2 1,i 0 0 . . . 0 |G 3/2,i | - h n+1/2 1,i . . . . . . 0 . . . 0 0 . . . . . . . . . 0 0 . . . 0 |G α-1/2,i | - h n+1/2 α,i |G α+1/2,i | + h n+1/2 α,i - |G α-1/2,i | - h n+1/2 α,i - |G α+1/2,i | + h n+1/2 α,i 0 . . . . . . 0 . . . . . . - |G N -1/2,i | + h n+1/2 N,i 0 . . . 0 0 |G N -1/2,i | - h n+1/2 N,i - |G N -1/2,i | - h n+1/2 N,i                 If h n+1/2 α,i
= 0 for all α, then we trivially have ρ n+1 α,i = 0 for all α. Let us now assume that there exists a layer α such that h n+1 α,i > 0. Then the matrix H n+1/2 N,i + ∆tG N,i is a strictly diagonally dominant matrix. Therefore, it is invertible and the entries of its inverse are all nonnegative -see the proof made in [START_REF] Audusse | Kinetic entropy for layer-averaged hydrostatic Navier-Stokes equations[END_REF]. This system and its numerical resolution are studied afterwards, in section 3.4.5.

The variable update

Hence the sum of relations (58) and (62) gives

U n+1 i = U i -. . . , (64) eq
:upU_HR

The space discretization of the system (38)-(40) and its numerical resolution allow to determine the quantities

h n+1 i , ρ n+1 α,i , (ρ α h α u α ) n+1 i ,
for any i ∈ I knowing the quantities {h j , (ρ α h α ) j , (ρ α h α u α ) j } j∈I . Thus, it is possible to recover h n+1 α,i , ρ n+1 α,i and u n+1 α,i from

h n+1 α,i = l α h n+1 i , ρ n+1 α,i = (ρ α h α ) n+1 i l α h n+1 i , u n+1 α,i = (ρ α h α u α ) n+1 i (ρ α h α ) n+1 i .
Concerning the temperature T n+1 α,i , we use the formula (ii ) preserves the steady state of a lake at rest, (iii ) is consistent with the system (28)-( 30),(31).

T n+1 α,i = ρ -1 (ρ n+1 α,i ).
n_properties ?

Proof (i) For α = 1, . . . , N , h n+1/2 α
is obtained using the fully discrete version of (35), that is to say,

h n+1/2 α,i = h α,i - j∈K i σ i,j F hα i,j
Since the flux F preserves the positivity for the shallow water equations, then h n+1/2 α,i is positive. As a sum of positive terms, h n+1 is positive. This proves (i).

(ii) In the constant density case and with a non-flat topography, (ii) is proved in [START_REF] Audusse | Approximation of the hydrostatic Navier-Stokes system for density stratified flows by a multilayer model. Kinetic interpretation and numerical validation[END_REF]. In the variable density case with flat topography, the proof is trivial because at equilibrium, the density is constant in each of the layers. In the general case, the continuous equations describing the static equilibrium are

u α = 0, α = 1, . . . , N, ∇ x,y η = 0, ∇ x,y (h α p α ) = p α+1/2 ∇ x,y z α+1/2 -p α-1/2 ∇ x,y z α-1/2 , α = 1, . . . , N.
(65) eq:static_eq_p

While the first two equations are easy to write at the discrete level, the third one is not.

There is no simple discrete formulation of (65). Yet, for u α,i = 0, ∀α, ∀i, and η i = η eq , ∀i (66) eq:static_equi with η eq a constant, due to the use of the HR technique, the proposed numerical scheme verifies

h n+1/2 α,i = h α,i , (ρh) n+1/2 α,i = (ρh) α,i , (ρh) n+1 α,i = (ρh) n+1/2 α,i , (ρhu) n+1 α,i = (ρhu) n+1/2 α,i
, for all α and all i. Therefore, starting from a situation described by (66) and a discrete version of (65), the discrete equilibrium is preserved. Starting from a situation near the equilibrium, the system will evolve towards equilibrium. We have empirically checked that a system initially far from equilibrium will reach an equilibrium, see [START_REF] Audusse | Approximation of the hydrostatic Navier-Stokes system for density stratified flows by a multilayer model. Kinetic interpretation and numerical validation[END_REF].

(iii) Since the flux F is consistent with the homogeneous (S p (U, z b ) = 0)), Eq. (60) is a consistent discretization of the non-conservative pressure terms, and (62) is a consistent discretization of the vertical exchange terms, then (iii) is true.

F ραhα i,j = ρ α,i,j F hα i,j , ( 67
) eq:flux_less_d with ρ α,i,j = ρ α,i if F hα i,j ≥ 0 ρ α,j if F hα i,j ≤ 0
(68) eq:upwind_rho then the numerical scheme (58)-( 62) satisfies a maximum principle on the density i.e. prop:max_prin for any α, i one has

ρ n+1 α,i ≤ max α {ρ α,i , ρ α,j }, ∀j ∈ K i .
Remark 3.5 The formula (67) was initially proposed in [START_REF] Larrouturou | How to preserve the mass fractions positivity when computing compressible multi-component flows[END_REF]) see also [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF].

Proof of prop 3.4 Let us first deal with the horizontal exchanges. Due to the choice of F ραhα i,j

, the numerical discretization of (36) can be decomposed into

(ρ α h α ) n+1/2 i = ρ α,i h n α,i - j∈K i σ i,j |F hα i,j | + - j∈K i σ i,j ρ α,j |F hα i,j | -
The right hand side of this expression is positive. Indeed, h n α,i -j∈K i σ i,j |F hα i,j | + is positive due to the CFL condition and

-j∈K i σ i,j ρ α,j |F hα i,j | -is positive because |F hα i,j | - is negative. The right hand side is factorized by max{ρ α,i , ρ α,j } (ρ α h α ) n+1/2 i ≤ max{ρ α,i , ρ α,j } h n α,i - j∈K i σ i,j |F hα i,j | + - j∈K i σ i,j |F hα i,j | -.
Therefore, dividing by h n+1/2 α,i positive, we obtain that + ∆tG N,i defined by (63) is an M -matrix. Let 1 be the vector the entries of which are all equal to 1. We notice that

ρ n+1/2 α,i ≤ max{ρ α,i , ρ α,j }, ∀j ∈ K i . (
(H n+1/2 N,i + ∆tG N,i )1 = h n+1/2 i , so we also have 1 = (H n+1/2 N,i + ∆tG N,i ) -1 h n+1/2 i
. Let (ρ i h i ) n+1/2 be the vector the entries of which are

(ρ i h i ) n+1/2 α = ρ n+1/2 α,i h n+1/2 α,i , α = 1, . . . , N . Then ||(H n+1/2 N,i + ∆tG N,i ) -1 (ρ i h i ) n+1/2 || ∞ ≤ ||ρ n+1/2 i || ∞ ||(H n+1/2 N,i + ∆tG N,i ) -1 h n+1/2 i || ∞ , which is exactly ||ρ n+1 i || ∞ ≤ ||ρ n+1/2 i || ∞ . ( 70
) maxpple_vertic To conclude, relationship (69) is applied to ρ n+1/2 α,i
for all α. Combining with (70) gives the maximum principle on the density. mark:max_prin Remark 3.6 Let us present a more accurate result for the maximum principle on the vertical exchanges. Let α, α 0 and α 1 such that 1 ≤ α 0 < α < α 1 ≤ N and

G α j +1/2,i < 0, G α j -1/2,i > 0 for j ∈ {0, 1}.
The coefficients of the lines α 0 and α 1 of matrix G N,i are respectively (G N,i ) α 0 ,j = δ α 0 ,j and (G N,i ) α 1 ,j = δ α 1 ,j with δ k,l the Kronecker symbol, which means that

ρ n+1 α 0 ,i = ρ n+1/2 α 0 ,i and ρ n+1 α 1 ,i = ρ n+1/2
α 1 ,i . The system can be solved using forward elimination and backward substitution. Denoting by ρ n+1/2 k-l,i the vector (ρ

n+1/2 k,i , ρ n+1/2 k+1,i , . . . , ρ n+1/2 l,i
) T , we have the following results

||ρ n+1 1-α 0 ,i || ∞ ≤ ||ρ n+1/2 1-α 0 ,i || ∞ , ||ρ n+1 α 0 -α 1 ,i || ∞ ≤ ||ρ n+1/2 α 0 -α 1 ,i || ∞ , ||ρ n+1 α 0 -N,i || ∞ ≤ ||ρ n+1/2 α 0 -N,i || ∞ .
The layers receiving no mass from the layers above and below them separate groups of layers which exchange mass between themselves.

Remark 3.7 For the following semi-implicit scheme for the vertical exchanges, the maximum principle on the density of proposition 3.4 is also verified:

(ρ α,i h α,i ) n+1 = (ρ α,i h α,i ) n+1/2 + ∆t 2 (ρ n+1 α+1/2,i G α+1/2,i -ρ n+1 α-1/2,i G α-1/2,i ) + ∆t 2 (ρ n+1/2 α+1/2,i G α+1/2,i -ρ n+1/2 α-1/2,i G α-1/2,i ).
The more accurate maximum principle presented in remark 3.6 is verified as well. However, in the rest of the paper, we work only with the fully implicit scheme for the vertical exchanges. More specifically, proposition 3.10 is stated only for the fully implicit scheme.

Remark 3.8

We have presented a first order in space discretization of the system. In practice, we apply a formally second order extension in space and time presented in [START_REF] Audusse | A well-balanced positivity preserving second-order scheme for Shallow Water flows on unstructured meshes[END_REF] and [START_REF] Allgeyer | Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface[END_REF]. More specifically, the modified Heun scheme which is used is presented in [START_REF] Allgeyer | Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface[END_REF]). For the obtained numerical scheme we are able to prove the consistence, the well-balancing and the nonnegativity of the water depth. But a discrete entropy inequality such as the one in proposition 3. [START_REF] Bristeau | Various analytical solutions for the incompressible Euler and Navier-Stokes systems with free surface[END_REF] has not yet been obtained. The proof of proposition 3.10 cannot be adapted for the second order scheme; a different strategy would be needed.

Remark 3.9 At each time step, to advance

• from h n i to h n+1/2 i • from (ρ α,i h α,i ) n to (ρ α,i h α,i ) n+1/2
convex combinations are used, which gives the scheme a certain stability. Then, to advance from (ρ α,i h α,i ) n+1/2 to (ρ α,i h α,i ) n+1 , the fact that the matrix of the system (matrix

H n+1/2 N,i
+ ∆tG N,i , defined in the proof of proposition 3.4) is an M -matrix gives stability to the computation.

Kinetic fluxes

ec:kin_fluxes Whereas the proposed numerical scheme can be adapted to any finite volume solver for the classical Saint-Venant system, in Section 5, the numerical simulations are performed using a kinetic solver and hence the numerical fluxes in (55) in the kinetic context now are specified.

To define the numerical fluxes, we introduce the functions χ 0 , M α

χ 0 (z 1 , z 2 ) = 1 4π 1 z 2 1 +z 2 2 ≤4 , M α = M (U α , ξ, γ) = h α c 2 α χ 0 ξ -u α c α , γ -v α c α ,
with c α = p α /ρ α , U α defined by (52) and where (ξ, γ) ∈ R 2 . We also define the quantity M ρ α by M ρ α = ρ α M α . The quantity M ρ α satisfies the following moment relations

R 2 1 ξ ρ α M (U α , ξ, γ) dξdγ =   ρ α h α ρ α h α u α ρ α h α v α   , ( 71
) eq:kinmom1 R 2   ξ 2 ξγ γ 2   ρ α M (U α , ξ, γ) dξdγ =   ρ α h α u 2 α + h α p α ρ α h α u α v α ρ α h α v 2 α + h α p α   . ( 72
) ?eq:kinmom2?
Hence in the context of the kinetic description, the fluxes appearing in (55) have the expressions

F hα i,j = R 2 M α,i,j ζ i,j dξdγ, F ραhα i,j = ρ α,i,j F hα i,j , F uα i,j = ρ α,i,j R 2 ξ γ M α,i,j ζ i,j dξdγ (73) eq:kin_fluxes with M α,i,j = M * α,i,j 1 ζ i,j ≥0 + M * α,j,i 1 ζ i,j ≤0 , ( 74
) eq:kin_dis22 where M * α,i,j = M (U * α,i,j , ξ, γ), U * α,i,j = (l α h * i,j , l α h * i,j u α,i , l α h * i,j v α,i ) T .
The * refers to the HR technique, see (61). The density ρ α,i,j is defined by (68) and

ζ i,j = ξ γ .n i,j .
We give here some details about Proposition 3.3 in the case of the kinetic flux. We consider the equation

f n+1/2- α,i = M α,i - j∈K i σ i,j ζ i,j 1 ζ i,j ≥0 M * α,i,j - j∈K i σ i,j M * α,j,i ζ i,j 1 ζ i,j ≤0
(75) eq:kin_horizon Using (71), we see that integrating equation (75) for α = 1, . . . , N with respect to ξ, γ gives the HR scheme for (35).

Let us now give the CFL condition for the kinetic flux. There exists a velocity v m ≥ 0 such that for all α, i

ξ -u α,i c α,i ≥ v m or γ -v α,i c α,i ≥ v m ⇒ M (U α,i , ξ, γ) = 0.
A CFL condition strictly less than one is considered

σi |u α,i | + |v α,i | + v m c α,i ≤ β < 1 2 for all i, α,
where σi = ∆t n Σ j∈K i L i,j /|C i |, and β is a given constant. More precisely, the CFL condition writes

∆t n ≤ 1 2 min α,i |C i | Σ j∈K i L i,j (|u α,i | + |v α,i | + v m c α,i )
.

Under this CFL condition, the kinetic function f n+1/2α,i remains non-negative, i.e.

f n+1/2- α,i ≥ 0, ∀(ξ, γ) ∈ R 2 , ∀i, ∀α.
The proof can be found in [START_REF] Allgeyer | Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface[END_REF]. Therefore, the water depth h n+1/2 is non-negative. Note that the CFL condition does not depend on the vertical exchange terms because they are treated implicitly.

Discrete entropy inequality

In this paragraph, a discrete entropy inequality is proved in the case of a flat topography. The crucial point of the numerical scheme is the treatment of the pressure source term S p,α written under the form (42). Indeed the other terms appearing in the numerical scheme are either conservative -and hence easily incorporated in the numerical fluxesor similar to terms appearing in the constant density case, see [START_REF] Allgeyer | Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface[END_REF]. The term S p,α is an extension of the topography term for the Saint-Venant system but in a far more complex setting. Proposition 3.10 is interesting since until now, the properties satisfied by the numerical scheme detailed in paragraph 3.4.5 concern Eqs. ( 35)-(41) except Eq. (37). Hence, for the momentum equation (37), only the equilibrium at rest is proved.

In the context of the kinetic description, the relation between the mass and momentum fluxes is simple (see (73)) and it is possible to slightly modify the definitions (60) in order to obtain an in cell discrete entropy inequality. The authors do not claim that only the kinetic fluxes allow to obtain such a result but, as in [START_REF] Audusse | Kinetic entropy inequality and hydrostatic reconstruction scheme for the Saint-Venant system[END_REF], it is not clear whether the result holds with other numerical fluxes in particular the definition (77) is related to the kinetic description and not easily available for others finite volume solvers (Rusanov, HLL,. . . ).

We consider in this paragraph a discrete form of (42) that is slightly different from (60) and defined by

S p,α,i,j = p α,i ( ĥα,i,j -h α,i )n i,j -gρ α,i h α,i (z α,i,j -z α,i )n i,j , ( 76
) eq:Sp_HR_mod with z α+1/2,i = z b,i + α l=1 h l,i , with z b,i = z b,j = cst, ∀j, z α,i = z α+1/2,i + z α-1/2,i 2 , z α,i,j = z α,i + z α,j 2 , ĥα,i,j = R 2 (M α,i 1 ζ i,j ≤0 + M α,j 1 ζ i,j ≥0 )dξdγ, = R 2
M α,i,j dξdγ (77) eq:h_alpha_i_j

p α,i = ρ α,i 2 R 2 ξ γ -u α,i 2 M α,i dξdγ R 2 M α,i dξdγ . ( 78
) eq:state_HR_mo
Since we consider a flat topography, the notations associated with the HR do not appear in the previous definitions. The discretization (76) is very similar to the discretization (60) written for a flat topography, but a bit of upwinding with respect to the advection is included in (76). Note that this kind of upwinding could help to reduce the error term due to the topography in [START_REF] Audusse | Kinetic entropy inequality and hydrostatic reconstruction scheme for the Saint-Venant system[END_REF]. For a non-flat topography, the HR induces upwinding with respect to the topography, not with respect to the advection.

The main interest of the following proposition is to justify the discretization (76) but for three reasons it is a partial result:

• it only concerns flat topography situations whereas it is well known that the numerical treatment of topography source terms is a very difficult issue,

• the extension of the expression (76) to the situation of a non flat topography is not natural since in the simple case of a single layer with a constant density i.e. the classical Saint-Venant system, the definition (76) does not exactly match with previous works of some of the authors (Audusse et al., 2016) (whereas definition (60) exactly reduces to the scheme studied in [START_REF] Audusse | Kinetic entropy inequality and hydrostatic reconstruction scheme for the Saint-Venant system[END_REF] in the Saint-Venant case),

• the numerical tests carried out with the two possible discretizations of S p,α,i,j , namely (60) and (76) give very similar results especially similar convergence curves, see paragraph 5.1.

For these reasons and even if the result proposed in the following proposition is an interesting stability property, the numerical simulations presented in Section 4 have been obtained using the discretization (60).

crete_entropy Proposition 3.10 When considering a flat bottom, the scheme (58),(62) with the fluxes defined by (73) satisfies an in cell fully discrete entropy inequality having the form

E n+1 α,i -E α,i + j∈K i σ i,j R 2 gz α,i,j + ξ 2 + γ 2 2 - 1 2 ξ γ -ũα,i,j 2 ρ α,i,j M α,i,j ζ i,j dξdγ -∆t n ρ n+1 α+1/2,i |u n+1 α+1/2,i | 2 2 + gz α+1/2,i G α+1/2,i +∆t n ρ n+1 α-1/2,i |u n+1 α-1/2,i | 2 2 + gz α-1/2,i G α-1/2,i -p n * α+1/2,i ∆t n G α+1/2,i -z n+1 α+1/2,i + z α+1/2,i + p n * α-1/2,i ∆t n G α-1/2,i -z n+1 α-1/2,i + z α-1/2,i = d α,i + e α,i + f α,i
where E α,i is the discrete energy

E α,i = ρ α,i h α,i |u α,i | 2 2 + ρ α,i gh α,i z α,i ,
and where d α,i is a sum of non-positive and hence dissipative terms whereas e α,i (resp. f α,i ) contains errors terms of magnitude O(diam(C i ) 3 ) (resp. O(∆t n ) 2 ). The expressions of d α,i and e α,i are given by

d α,i = - j∈K i σ i,j |ũ α,i,j -u α,i | 2 2 R 2 ρ α,i M α,i |ζ i,j |dξdγ + j∈K i σ i,j |u α,i -u α,j | 2 4 {ζ i,j ≤0} ρ α,j M α,j ζ i,j dξdγ (79) def e α,i = j∈K i σ i,j g R 2 ρ α,i h α,i u α,i -ρ α,i,j M α,i,j ξ γ (z α,i,j -z α,i ) • n i,j dξdγ - j∈K i σ i,j p α,i R 2 (M α,i,j -M α,i )(ũ α,i,j -u α,i ) • n i,j dξdγ - j∈K i σ i,j (ũ α,i,j -u α,i ) • {ζ i,j ≤0} ξ γ -u α,i ρ α,i M α,i - ξ γ -u α,j ρ α,j M α,j ζ i,j dξdγ - j∈K i σ i,j |2ũ α,i,j -u α,i -u α,j | 2 4 {ζ i,j ≤0} ρ α,j M α,j ζ i,j dξdγ. (80) def
The quantity f α,i is defined by

f α,i = ∆t n gh α,i 2 (ρ n+1 α,i -ρ n+1 α+1/2,i )G α+1/2,i + (ρ n+1 α,i -ρ n+1 α-1/2,i )G α-1/2,i +g (ρ α,i h α,i ) n+1 -ρ n+1 α,i h α,i (z n+1 α,i -z α,i ) + (ρ α,i h α,i ) n+1 2 |u n+1 α,i -u α,i | 2 -∆t n ρ n+1 α+1/2,i u α,i -u n+1 α,i 2 .(u n+1 α,i -2u n+1 α+1/2,i + u α,i )G α+1/2,i +∆t n ρ n+1 α-1/2,i u α,i -u n+1 α,i 2 .(u n+1 α,i -2u n+1 α-1/2,i + u α,i )G α-1/2,i ,
where the first line is the discrete version of the term appearing in the last line of the continuous energy balance given in Eq. (26) (see also [START_REF] Boittin | Low-Mach type approximation of the Navier-Stokes system with temperature and salinity for free surface flows[END_REF]) and the other lines come from the explicit time scheme and vanish in the semi-discrete (in space) case. The velocity ũα,i,j is defined by

ũα,i,j = R 2 M α,i,j ξ γ dξdγ R 2 M α,i,j dξdγ , ( 81
) eq:utilde
where M α,i,j is defined by (74) and the pressure p n * α±1/2,i is defined by

p n * α±1/2,i = p α,i ∓ ρ n+1 α,i g h α,i 2 .
:entropy_ineq Remark 3.11 Since we use an explicit time scheme it is natural to have error terms f α,i of order O(∆t n ) 2 . Concerning the error terms e α,i due to the space discretization, we point out that they are of order O(diam(C i ) 3 ) i.e. smaller than residuals with second terms.

Proof of prop. 3.10

The proof of this proposition is long but only contains simple computations. The authors have not found a simpler presentation.

Starting from the set of discrete equations (64), the energy balance for the cell i at the layer α is obtained by performing the sum of the two following quantities:

• the mass conservation equation over the layer α in (64) multiplied by gz α,i -|u α,i | 2 /2

• the momentum equation in (64) over the layer α multiplied by u α,i .

In other words, the energy balance comes from a rewriting of the quantity

E α,i := gz α,i - |u α,i | 2 2 (ρ α,i h α,i ) n+1 -ρ α,i h α,i + j∈K i σ i,j ρ α,i,j F hα i,j +∆t n ρ n+1 α+1/2,i G α+1/2,i -ρ n+1 α-1/2,i G α-1/2,i +u α,i . (ρ α,i h α,i u α,i ) n+1 -ρ α,i h α,i u α,i + j∈K i σ i,j (F hαu α,i i,j -S p,α,i,j ) -∆t n ρ n+1 α+1/2,i u n+1 α+1/2,i G α+1/2,i -ρ n+1 α-1/2,i u n+1 α-1/2,i G α-1/2,i .
(82) eq:energ_init

Since the manipulations necessary to obtain the result are, to some extent tedious, we proceed as follows: first we consider the terms involving time derivatives then those involving the horizontal fluxes and finally, we consider the vertical exchange terms.

The discrete time derivatives The terms appearing in (82) writes

E t α,i := gz α,i - |u α,i | 2 2 (ρ α,i h α,i ) n+1 -ρ α,i h α,i + u α,i . (ρ α,i h α,i u α,i ) n+1 -ρ α,i h α,i u α,i ,
and can be rewritten under the form

E t α,i = (ρ α,i h α,i ) n+1 |u n+1 α,i | 2 2 -(ρ α,i h α,i ) |u α,i | 2 2 + (ρ α,i gh α,i z α,i ) n+1 -(ρ α,i gh α,i z α,i ) -(ρ α,i gh α,i ) n+1 (z n+1 α,i -z α,i ) - (ρ α,i h α,i ) n+1 2 |u n+1 α,i -u α,i | 2 . ( 83
) eq:time_d
The last term in (83) is classical in the context of explicit in time schemes and give rise to a non-negative term in the energy balance.

The horizontal fluxes The quantities we are now considering are

E xy α,i,j := gz α,i - |u α,i | 2 2 ρ α,i,j R 2 M α,i,j ζ i,j dξdγ + u α,i . ρ α,i,j R 2 ξ γ M α,i,j ζ i,j dξdγ -S p,α,i,j .
And using (76) we rewrite E xy α,i,j under the form

E xy α,i,j := R 2 gz α,i,j + ξ 2 + γ 2 2 - 1 2 ξ γ -ũα,i,j 2 ρ α,i,j M α,i,j ζ i,j dξdγ + R 2 1 2 ξ γ -ũα,i,j 2 - ξ γ -u α,i 2 ρ α,i,j M α,i,j ζ i,j dξdγ -p α,i ( ĥα,i,j -h α,i )u α,i .n i,j +g R 2 ρ α,i,j M α,i,j z α,i,j -z α,i u α,i .n i,j -ζ i,j dξdγ (84) {?} +g R 2 (ρ α,i h α,i -ρ α,i,j (z α,i,j -z α,i )u α,i • n i,j dξdγ.
(85) eq:E_ij_xy

The last two lines of the previous equation reduce to

- j∈K i σ i,j p α,i R 2 (M α,i,j -M α,i )(ũ α,i,j -u α,i ) • n i,j dξdγ.
Now, using the definition (81), we rewrite the second line of (85), denoted E 2,xy α,i,j , under the form

E 2,xy α,i,j = - 1 2 (ũ α,i,j -u α,i ).n i,j R 2 ξ γ -u α,i 2 ρ α,i M α,i dξdγ - R 2 1 2 ξ γ -u α,i 2 ρ α,i,j M α,i,j ζ i,j dξdγ + R 2 1 2 ξ γ -ũα,i,j 2 ρ α,i,j M α,i,j ζ i,j dξdγ + 1 2 (ũ α,i,j -u α,i ).n i,j R 2 ξ γ -u α,i 2 ρ α,i M α,i dξdγ = -p α,i R 2 M α,i (ũ α,i,j -u α,i ).n i,j dξdγ -(ũ α,i,j -u α,i ). R 2 ξ γ - ũα,i,j + u α,i 2 ρ α,i,j M α,i,j ζ i,j dξdγ +(ũ α,i,j -u α,i ).n i,j R 2 1 2 ξ γ -u α,i 2 ρ α,i M α,i dξdγ, ( 86
) eq:e2_xy_fin
where the definition (78) has been used. Let n i,j = (n 1,i,j , n 2,i,j ) T , because the Gibbs equilibrium M α,i is an even function of the variables ξ -u α,i and γ -v α,i , we get

1 2 R 2 ξ γ -u α,i 2 ρ α,i M α,i dξdγn 1,i,j = R 2 (ξ -u α,i ) 2 + (γ -v α,i ) 2 2 ρ α,i M α,i dξdγn 1,i,j = R 2 (ξ -u α,i ) 2 ρ α,i M α,i dξdγn 1,i,j = R 2 (ξ -u α,i )ξn 1,i,j ρ α,i M α,i dξdγ = R 2 (ξ -u α,i )ζ i,j ρ α,i M α,i dξdγ, ( 87
) eq:even_odd_1
where we have used that since the function (ξ -u α,i )M α,i is even then

R 2 (ξ -u α,i )γn 2,i,j ρ α,i M α,i dξdγ = 0.
Likewise, we obtain

1 2 R 2 ξ γ -u α,i 2 ρ α,i M α,i dξdγn 2,i,j = R 2 (γ -v α,i )ζ i,j ρ α,i M α,i dξdγ.
(88) eq:even_odd_2

And therefore, using Eqs. (87),(88), we can rewrite the last line of relation (86) under the form

(ũ α,i,j -u α,i ). R 2 ξ γ -u α,i ρ α,i M α,i ζ i,j dξdγ,
leading to the following expression for E 2,xy α,i,j

E 2,xy α,i,j = -p α,i R 2 M α,i (ũ α,i,j -u α,i ).n i,j dξdγ -(ũ α,i,j -u α,i ). R 2 ξ γ - ũα,i,j + u α,i 2 ρ α,i,j M α,i,j ζ i,j dξdγ +(ũ α,i,j -u α,i ). R 2 ξ γ -u α,i ρ α,i M α,i ζ i,j dξdγ.
We rewrite the second line of E 2,xy α,i,j under the form

|ũ α,i,j -u α,i | 2 2 R 2 ρ α,i,j M α,i,j ζ i,j dξdγ-(ũ α,i,j -u α,i ). R 2 ξ γ -u α,i ρ α,i,j M α,i,j ζ i,j dξdγ, or equivalently |ũ α,i,j -u α,i | 2 2 R 2 ρ α,i,j M α,i,j ζ i,j dξdγ-(ũ α,i,j -u α,i ). R 2 ξ γ -u α,i,j ρ α,i,j M α,i,j ζ i,j dξdγ + (ũ α,i,j -u α,i ). R 2
(u α,i -u α,i,j )ρ α,i,j M α,i,j ζ i,j dξdγ.

Using the previous identities and considering the cases ζ i,j ≥ 0 and ζ i,j ≤ 0, simple computations give the following expression for E 2,xy α,i,j

E 2,xy α,i,j = -p α,i R 2 M α,i (ũ α,i,j -u α,i ).n i,j dξdγ + |ũ α,i,j -u α,i | 2 2 {ζ i,j ≥0} ρ α,i M α,i ζ i,j dξdγ +(ũ α,i,j -u α,i ). {ζ i,j ≤0} ξ γ -u α,i ρ α,i M α,i - ξ γ -u α,j ρ α,j M α,j ζ i,j dξdγ + 1 2 (ũ α,i,j -u α,i ). {ζ i,j ≤0} ũα,i,j -u α,i -2(u α,j -u α,i ) ρ α,j M α,j ζ i,j dξdγ,
and the last line of the previous expression can be written under the form

- 1 2 |ũ α,i,j -u α,i | 2 {ζ i,j ≤0} ρ α,j M α,j ζ i,j dξdγ + (ũ α,i,j -u α,i ).(ũ α,i,j -u α,j ) {ζ i,j ≤0} ρ α,j M α,j ζ i,j dξdγ, or equivalently - 1 2 |ũ α,i,j -u α,i | 2 {ζ i,j ≤0} ρ α,j M α,j ζ i,j dξdγ - |u α,i -u α,j | 2 4 {ζ i,j ≤0} ρ α,j M α,j ζ i,j dξdγ+ |2ũ α,i,j -u α,i -u α,j | 2 4 {ζ i,j ≤0} ρ α,j M α,j ζ i,j dξdγ,
where the identity ab = (a + b) 2 /4 -(a -b) 2 /4 has been used. Hence, the final expression for E 2,xy α,i,j is given by

E 2,xy α,i,j = -p α,i R 2 M α,i (ũ α,i,j -u α,i ).n i,j dξdγ + |ũ α,i,j -u α,i | 2 2 R 2 ρ α,i M α,i |ζ i,j |dξdγ - |u α,i -u α,j | 2 4 {ζ i,j ≤0} ρ α,j M α,j ζ i,j dξdγ +(ũ α,i,j -u α,i ). {ζ i,j ≤0} ξ γ -u α,i ρ α,i M α,i - ξ γ -u α,j ρ α,j M α,j ζ i,j dξdγ + |2ũ α,i,j -u α,i -u α,j | 2 4 {ζ i,j ≤0} ρ α,j M α,j ζ i,j dξdγ,
where in the previous expression, the second line is nonnegative and the last two lines are third order terms (i.e. of order O(diam(C i ) 3 )) when considering Lipschitz continuous solutions.

Then for the third line of (85), performing simple manipulations we have

P xy α,i,j = -p α,i ( ĥα,i,j -h α,i )u α,i .n i,j = -p α,i R 2 (M α,i,j -M α,i )dξdγu α,i .n i,j = -p α,i R 2 (M α,i,j ũα,i,j .n i,j -M α,i u α,i .n i,j )dξdγ + p α,i R 2 M α,i,j (ũ α,i,j -u α,i ).n i,j dξdγ = -p α,i R 2 (M α,i,j ζ i,j -M α,i u α,i .n i,j )dξdγ + p α,i R 2 M α,i,j (ũ α,i,j -u α,i ).n i,j dξdγ = -p α,i R 2 (M α,i,j ζ i,j -M α,i u α,i .n i,j )dξdγ + p α,i R 2 M α,i (ũ α,i,j -u α,i ).n i,j dξdγ + p α,i R 2 (M α,i,j -M α,i )(ũ α,i,j -u α,i ).n i,j dξdγ (89
) eq:p1_cons

where the definitions (81),(77) have been used. Notice that for the second term in the first line of (89), we have

j∈K i p α,i R 2 M α,i u α,i .n i,j dξdγ = p α,i R 2 M α,i u α,i dξdγ. j∈K i n i,j = 0,
and using the discrete form of the continuity equation, for the first term in the first line of (89) we get

- j∈K i σ i,j p α,i R 2 M α,i,j ζ i,j dξdγ = p α,i h n+1 α,i -h α,i -∆t n (G α+1/2,i -G α-1/2,i ) .
The vertical exchange terms It remains to examine the contribution of the vertical exchange terms over the energy balance, namely in Eq. (82) the quantity

V α,i := ∆t n gz α,i - |u α,i | 2 2 ρ n+1 α+1/2,i G α+1/2,i -ρ n+1 α-1/2,i G α-1/2,i +∆t n u α,i . ρ n+1 α+1/2,i u n+1 α+1/2,i G α+1/2,i -ρ n+1 α-1/2,i u n+1 α-1/2,i G α-1/2,i .
And we write

V α,i = ∆t n ρ n+1 α+1/2,i |u n+1 α+1/2,i | 2 2 G α+1/2,i -∆t n ρ n+1 α-1/2,i |u n+1 α-1/2,i | 2 2 G α-1/2,i +∆t n g ρ n+1 α+1/2,i z α+1/2,i G α+1/2,i -ρ n+1 α-1/2,i z α-1/2,i G α-1/2,i -∆t n ρ n+1 α+1/2,i |u n+1 α+1/2,i -u n+1 α,i | 2 2 + u α,i -u n+1 α,i 2 .(u n+1 α,i -2u n+1 α+1/2,i + u α,i ) G α+1/2,i +∆t n ρ n+1 α-1/2,i |u n+1 α-1/2,i -u n+1 α,i | 2 2 + u α,i -u n+1 α,i 2 .(u n+1 α,i -2u n+1 α-1/2,i + u α,i ) G α-1/2,i -∆t n g h α,i 2 ρ n+1 α+1/2,i G α+1/2,i + ρ n+1 α-1/2,i G α-1/2,i .
(90) eq:exchang Remark 3.12 When G α+1/2,i ≤ 0, then the third line of (90) is nonnegative namely,

-∆t n ρ n+1 α+1/2,i |u α,i -u n+1 α,i | 2 2 G α+1/2,i = O (∆t n ) 3 , whereas for G α+1/2,i > 0 |u n+1 α+1/2,i -u n+1 α,i | 2 2 - u α,i -u n+1 α,i 2 .(u n+1 α,i -2u n+1 α+1/2,i + u α,i ) = u n+1 α+1,i -u n+1 α,i 2 (u n+1 α+1,i -2u α,i + u n+1 α,i ) = O (∆t n ) 2 .
Hence, when it is not a dissipative term, the third line of (90) is a O (∆t n ) 3 term. The same result holds for the fourth line of (90).

All the contributions Now, summarizing the computations carried out in the previous paragraphs, we sum all the contributions, namely E t α,i , E xy α,i,j and V α,i leading to a new expression for E α,i under the form

E n+1 α,i -E α,i + j∈K i σ i,j R 2 gz α,i,j + ξ 2 + γ 2 2 - 1 2 ξ γ -ũα,i,j 2 ρ α,i,j M α,i,j ζ i,j dξdγ -∆t n ρ n+1 α+1/2,i |u n+1 α+1/2,i | 2 2 + gz α+1/2,i G α+1/2,i +∆t n ρ n+1 α-1/2,i |u n+1 α-1/2,i | 2 2 + gz α-1/2,i G α-1/2,i = d α,i + e α,i + F α,i
where d α,i are non-positive and hence dissipative terms whereas e α,i are errors terms. d α,i and e α,i are given by (79) and (80) respectively. The quantity F α,i is defined by

F α,i = (ρ α,i gh α,i ) n+1 (z n+1 α,i -z α,i ) -p α,i h n+1 α,i -h α,i -∆t n (G α+1/2,i -G α-1/2,i ) -∆t n ρ n+1 α+1/2,i g h α,i 2 G α+1/2,i -∆t n ρ n+1 α-1/2,i g h α,i 2 G α-1/2,i + (ρ α,i h α,i ) n+1 2 |u n+1 α,i -u α,i | 2 -∆t n ρ n+1 α+1/2,i |u n+1 α+1/2,i -u n+1 α,i | 2 2 + u α,i -u n+1 α,i 2 .(u n+1 α,i -2u n+1 α+1/2,i + u α,i ) G α+1/2,i +∆t n ρ n+1 α-1/2,i |u n+1 α-1/2,i -u n+1 α,i | 2 2 + u α,i -u n+1 α,i 2 .(u n+1 α,i -2u n+1 α-1/2,i + u α,i ) G α-1/2,i ,
and the first three lines of F α,i can be rewritten under the form

F 1 α,i = (ρ α,i gh α,i ) n+1 (z n+1 α,i -z α,i ) -p α,i h n+1 α,i -h α,i + (ρ α,i h α,i ) n+1 2 |u n+1 α,i -u α,i | 2 +∆t n p n * α+1/2,i G α+1/2,i -p n * α-1/2,i G α-1/2,i +∆t n gh α,i 2 (ρ n+1 α,i -ρ n+1 α+1/2,i )G α+1/2,i + (ρ n+1 α,i -ρ n+1 α-1/2,i )G α-1/2,i . (91
) eq:f_alpha

Now, we rewrite the first two terms of the first line of (91) under the form

F 2 α,i = (ρ α,i gh α,i ) n+1 2 (z n+1 α+1/2,i + z n+1 α-1/2,i -z α+1/2,i -z α-1/2,i ) -p α,i z n+1 α+1/2,i -z n+1 α-1/2,i -z α+1/2,i + z α-1/2,i = -p n * α+1/2,i (z n+1 α+1/2,i -z α+1/2,i ) + p n * α-1/2,i (z n+1 α-1/2,i -z α-1/2,i ) +g (ρ α,i h α,i ) n+1 -ρ n+1 α,i h α,i (z n+1 α,i -z α,i ).
Hence, we have for F α,i

F α,i = p n * α+1/2,i ∆t n G α+1/2,i -z n+1 α+1/2,i + z α+1/2,i -p n * α-1/2,i ∆t n G α-1/2,i -z n+1 α-1/2,i + z α-1/2,i -∆t n ρ n+1 α+1/2,i |u n+1 α+1/2,i -u n+1 α,i | 2 2 + u α,i -u n+1 α,i 2 .(u n+1 α,i -2u n+1 α+1/2,i + u α,i ) G α+1/2,i +∆t n ρ n+1 α-1/2,i |u n+1 α-1/2,i -u n+1 α,i | 2 2 + u α,i -u n+1 α,i 2 .(u n+1 α,i -2u n+1 α-1/2,i + u α,i ) G α-1/2,i +∆t n gh α,i 2 (ρ n+1 α,i -ρ n+1 α+1/2,i )G α+1/2,i + (ρ n+1 α,i -ρ n+1 α-1/2,i )G α-1/2,i +g (ρ α,i h α,i ) n+1 -ρ n+1 α,i h α,i (z n+1 α,i -z α,i ) + (ρ α,i h α,i ) n+1 2 |u n+1 α,i -u α,i | 2 ,
the first two lines being conservatives terms and the four following ones being error terms corresponding to f α,i . The fifth line is the discrete version of the term appearing in the last line of the continuous energy balance given in Eq. (26) (see also [START_REF] Boittin | Low-Mach type approximation of the Navier-Stokes system with temperature and salinity for free surface flows[END_REF]) and the last line of the previous relation comes from the time scheme and is of order

O (∆t n ) 2 .
In order to conclude the proof, it remains to prove that all the quantities appearing in e α,i are third order terms i.e. of magnitude O(diam(C i ) 3 ). The terms in each of the sums in e α,i are obviously second-order terms. Since the sum is made on the faces, gradients of second-order terms appear. These gradients are indeed third-order terms.

If one wishes to introduce even more upwinding in the discretization of S p,α by using the discretization S p,α,i,j = p α,i ( ĥα,i,j -h α,i )n i,j -gρ α,i,j ĥα,i,j (z α,i,j -z α,i )n i,j , Proposition 3.10 still holds but the expression of the rest term e α,i becomes

e α,i = j∈K i σ i,j g R 2 ρ α,i,j M α,i,j z α,i,j -z α,i (u α,i -u α,i,j ) • n i,j dξdγ + j∈K i σ i,j g R 2 ρ α,i,j M α,i,j z α,i,j -z α,i u α,i,j • n i,j -ζ i,j dξdγ - j∈K i σ i,j p α,i R 2 (M α,i,j -M α,i )(ũ α,i,j -u α,i ) • n i,j dξdγ - j∈K i σ i,j (ũ α,i,j -u α,i ) • {ζ i,j ≤0} ξ γ -u α,i ρ α,i M α,i - ξ γ -u α,j ρ α,j M α,j ζ i,j dξdγ - j∈K i σ i,j |2ũ α,i,j -u α,i -u α,j | 2 4 {ζ i,j ≤0} ρ α,j M α,j ζ i,j dξdγ.
This new expression of e α,i also contains only third-order error terms, thus the result is not deteriorated.

4 Numerical scheme for the layer-averaged Navier-Stokes-Fourier system sec:NS_num The discretization of the full layer-averaged Navier-Stokes-Fourier is presented. The main difficulties have already been tackled in section 3.

Semi-discrete (in time) scheme

The semi-discrete in time scheme (34) yields the following system

h n+1 = h n+1/2 = h -∆t n N α=1 ∇ x,y • (h α u α ) - N α=1 ρ (T α ) ρ 2 α c p (S T,α -S µ,α ), ( 92
) eq:NSF_sd_1 (ρ α h α ) n+1/2 = ρ α h α -∆t n ∇ x,y • (ρ α h α u α ) , (93)
?eq:NSF_sd_2?

(ρ α h α u α ) n+1/2 = ρ α h α u α -∆t n ∇ x,y • (ρ α h α u α ⊗ u α ) + ∇ x,y h α p α -p α ∇ x,y h α + ρ α gh α ∇ x,y z α , (94)
?eq:NSF_sd_3?

(ρ α h α ) n+1 = (ρ α h α ) n+1/2 -∆t n ρ n+1 α+1/2 G α+1/2 -ρ n+1 α-1/2 G α-1/2 , (95)
?eq:NSF_sd_4?

(ρ α h α u α ) n+1 = (ρ α h α u α ) n+1/2 -∆t n u n+1 α+1/2 ρ n+1 α+1/2 G α+1/2 -u n+1 α-1/2 ρ n+1 α-1/2 G α-1/2 + ∇ x,y • (µh n+1 α ∇ x,y u n+l α ) + Γ n+1 α+1/2 (u n+1 α+l -u n+l α ) -Γ n+1 α-1/2 (u n+l α -u n+l α-1 ) -κ α u n+l α , ( 96
) eq:NSF_sd_5
where l = 0, 1 and with

G α+1/2 = - N j=1 α p=1 l p -1 j≤α ∇ x,y • (h j u j ) + α j=1 ρ (T j ) ρ 2 j c p (S T,j -S µ,j ).
and Γ α+1/2 , S T,α , and S µ,α respectively defined by ( 22),( 20) and [START_REF] Thacker | Some exact solutions to the non-linear shallow-water wave equations[END_REF]. Note that the definition of the mass exchange terms G α+1/2 is different from the definition of the mass exchange terms for the Euler system given in (41).

Spatial discretization of the diffusion terms

The Euler part of the system is discretized in space as in section 3.4. We present here only the discretization of the diffusion terms. For the discretization of the viscosity terms in the momentum equation, we refer to [START_REF] Allgeyer | Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface[END_REF]. A classical P 1 finite element type approximation with mass-lumping is used. Let us define the number of cells N x as well as the vector of unknowns in layer α

U α = (h α,1 , . . . , h α,Nx , (ρ α h α ) 1 , . . . , (ρ α h α ) N x , (ρ α h α u α ) 1 , . . . , (ρ α h α u α ) N x , (ρ α h α v α ) 1 , . . . , (ρ α h α v α ) N x ) T
5 Numerical validation ec:validation In this section, we confront the proposed numerical scheme to three test cases. In the first one we present convergence results in the case of a 3d analytical solution for the Euler system. This test case is 3d, non-stationary with wet/dry interfaces and density variations and hence it is a challenging problem to obtain convergence curves towards the analytical solution. The second test deals with the simulation of the lock exchange phenomenon and the comparison with experimental results available in the literature. The third test consists in two simple diffusion cases for which we present a validation with an analytical solution and comparisons between the Navier-Stokes-Fourier and Boussinesq models. All the presented simulations have been obtained with the numerical code Freshkiss3d (2020).

Parabolic bowl

ubsec:thacker In [START_REF] Bristeau | Various analytical solutions for the incompressible Euler and Navier-Stokes systems with free surface[END_REF], the authors have proposed an analytic solution for the Euler system i.e. the system (1)-( 3) with λ = 0, µ = 0. It is an extension of the Thacker' analytical solution [START_REF] Thacker | Some exact solutions to the non-linear shallow-water wave equations[END_REF], corresponding to a periodic oscillation in a parabolic bowl.

The following proposition holds, its proof is detailed in [START_REF] Bristeau | Various analytical solutions for the incompressible Euler and Navier-Stokes systems with free surface[END_REF].

Proposition 5.1 For any nonnegative function s → ρ(s) and for some (a, α, η, h 0 ) ∈ R 3 × R + , let us consider the functions h, u, v, w, p, φ defined for (x, y) ∈ [-L/2, L/2] 2 , t ≥ t 0 by h(t, x, y) = max 0, h 0 -α (x -η cos(ωt)) 2 + (y -η sin(ωt)) 2 2 , u(t, x, y, z) = -ηω sin(ωt), v(t, x, y, z) = ηω cos(ωt), w(t, x, y, z) = -αηω (x sin(ωt) -y cos(ωt)) , p(t, x, y, z) = p a (t) +

h+z b z ρ(T (t, x, y, z 1 ))dz 1 , T (t, x, y, z) = a (h + z b -z) ,
with ω = √ αg and with a bottom topography defined by z b (x, y) = α 2 (x 2 + y 2 ), then h, u, v, w, p, φ as defined previously satisfy the 3D hydrostatic Euler system with variable density (system (1)-( 3) with λ = 0, µ = 0) completed with the kinematic boundary conditions (4), (6).

NS2d_var_dens

For the numerical validation the parameters are set to η = 0.1, h 0 = 0.1, a = 10, α = 1 and L = 4 and we consider a simplified state law given by ρ(T ) = ρ 0 + βT with ρ 0 = 1000 and β = 10. The free surface is plotted at different times in Figure 5 highlighting the planar motion of the fluid in the bowl. On Figure 6, the density in the density_slice slice plane (x, y=0, z) after a period for a mesh with 31316 triangles and 30 layers is plotted.

The convergence towards the analytic solution is assessed by plotting the logarithm of the cumulative error (in L 2 -norm) at time τ = 2π/ω versus the space discretization (i.e. log(l 0 /l i ) where l 0 is the average edge length of the mesh i and l 0 the average edge length of the coarsest mesh) for several unstructured meshes with 934, 2194, 4020, 6408, 9066, 12674 triangles. More precisely, the error in L 2 -norm is computed by summing on all the nodes of the mesh at each time step, for each layer. Then, the cumulative error is obtained by summing the errors at each time step, normalized by the time step.

Figures 7 and8 show the cumulative errors obtained with a constant number of layers equal to 10 (the mesh is refined in the horizontal direction but not in the vertical direc-tion). The analytic solution being non-stationary, errors in time and space accumulate over time and the theoretical rate of convergence is thus hard to obtain. show the cumulative error in L 2 norm for meshes with 934, 2194, 4020, 6408, 9066, 12674 triangles and 10, 12, 14, 16, 18, 20 layers respectively. Increasing the number of layers while refining the horizontal mesh is a reasonable idea because by doing so, the proportions of the 3D wedge cells are preserved. A super-convergence phenomenon can be observed for ρh when the number of layers is increased as the mesh is refined. A rate of convergence higher than the theoretical rate is also observed for ρhu and ρhv. The faster one refines the mesh in the vertical direction, the higher the convergence rates obtained. The results shown on figures 9 and 10 prove the stability of the numerical scheme for the Euler system. v_convergence

To summarize the results on the convergence of the schemes and to overcome the difficulty of interpreting the super-convergence cases, we propose another way of plotting the results in the case of multi-layer models. The error is plotted as a function of the horizontal space step l i and of the vertical proportion l p = 1/N at the same time. Figure 11 shows the horizontal convergence rates of the first-and second-order schemes, as well as the first order vertical convergence. Note that the vertical discretization proportion does not correspond to the vertical space step because the water depth varies (h α = l p h). Initially, fluids of different densities, ρ 1 = 1090kg.m -3 and ρ 0 = 1000kg.m -3 respectively, are at rest and separated by a wall located at x 0 = 0.3m. When the vertical barrier is removed, the denser fluid flows under the lighter one due to the difference in the hydrostatic pressure. The initial water height is h 0 = 0.3m and the length of the domain is L = 3m. The reduced gravity is defined by g * = g (1 -γ) where γ = ρ 0 /ρ 1 is the density ratio. We also define the buoyancy velocity as u b = √ g * h 0 . The Navier-Stokes equations are usually made dimensionless using the Grashof number defined by

Gr = u b h 0 ν 2 , (97)

Comparison with the Boussinesq assumption

In this section we mainly test the effects of the temperature diffusion and we compare the solutions obtained with the model ( 1)-(3) to the solutions with the Boussinesq assumption defined below. So, we introduce the system

∇ • U = 0, ( 98
) eq:div_boussi_ ρ 0 c p ∂T ∂t + ∇ • (T U) = ∇ • (λ∇T ), (99)
?eq:temp_bouss

ρ 0 ∂u ∂t + ∇ x,y • (u ⊗ u) + ∂(uw) ∂z + ∇ x,y η z ρgdz = 0, (
100) eq:mom_boussi_ i.e. where the Boussinesq assumption is madeρ 0 is a constant. We consider simple diffusion cases in which a basin is initially at rest i.e. U = (0, 0, 0) T and the free surface and bottom are flat and equal to h(t 0 , x, y) = h 0 and z b (t 0 , x, y) = 0 respectively. In the basin, the temperature is initially distributed such that T (t 0 , x, y, z) = T 0 (z), where T 0 = T 0 (z) is a given function.

Starting from the initial conditions described above, it is easy to see that the velocity and density for t ≥ t 0 are given by U = (0, 0, 0) T and ρ = ρ(T ), where T = T (t, z) is governed by the heat equation:

   ∂T ∂t = D 0 ∂ 2 T ∂z 2 , T (t 0 , z) = T 0 (z).
(101) eq:heat_boussi

The coefficient D 0 is the diffusivity defined by D 0 = λ ρ 0 cp . Now considering the system (1)-(3) with µ = 0 i.e. without the Boussinesq assumption and starting also from the initial conditions described above, it is easy to see that the velocity and temperature for t ≥ t 0 is given by U = (0, 0, w) T and T = T (ρ), where w is defined by:

w = - λ c p z z b ρ ρ 2 ∂ 2 T ∂z 2 dz,
where ρ = ∂ρ ∂T is deduced from the state law and ρ = ρ(t, z) is governed by the equation:

   ∂ρ ∂t + w ∂ρ ∂z = λ c p ρ ρ ∂ 2 T ∂z 2 , ρ(t 0 , z) = ρ 0 (z).
(102) eq:heat_non_bo

The velocity w is the result of the fluid dilatation. When w = 0 the two models are almost identical, the only difference being the diffusivity D. It is either a constant in the Boussinesq model or defined by D = λ ρcp in the case of the Navier-Stokes-Fourier model. In the following examples we show that for common water state laws, w is sufficiently small to have no noticeable effect on T and ρ. However, this term is essential in order to obtain rigorous mass conservation.

Temperature diffusion

In this test, the initial temperature T 0 (and density ρ 0 ) is constant in the domain and a Dirichlet boundary condition is applied at the bottom with a temperature of T b = 0 (cf. figure 15). At the free surface we impose a homogeneous Neumann boundary condition (φ T | η = 0, where φ T is the thermal flux). We compare the temperature obtained by a numerical resolution of the system (92)-(96) (with a simplified state law given by T (ρ) = ρ 0 -ρ β with ρ 0 = 1000 and β = 10) to the analytical solution (104). Notice that from Eq. (102) we obtained a slightly different dimensionless heat equation

∂ T ∂ t - ρ(T 0 -T b ) h 0 ∂ T ∂ z z z b ρ ρ 2 ∂ T ∂ z dz = ∂ 2 T ∂ z2 .
We obtain a good matching between the numerical results and the analytical solution both with the Navier-Stokes-Fourier system (with µ = 0) and the Boussinesq system (98)-(100) which validates the numerical treatment of the diffusion. The results are plotted on Figure 17 in the case of the Navier-Stokes-Fourier model. However, there is no noticeable difference between the analytical solution and the Navier-Stokes-Fourier numerical solution even though ρ = -β has been chosen higher than its real physical value (in the case of water 0.03 < ρ < 0.13). tion_analytic Note that the mass is strictly conserved in the Navier-Stokes-Fourier model, even though the density varies. Dilatation effects induce a variation of the water height such that the integral of the density over the total volume of fluid stays constant over time. This is not the case with the Boussinesq model, where the volume is not affected by temperature variation and mass conservation is violated (cf. Figure 18).

Thermal equilibrium

In this second test a well stratified fluid is considered. No exterior forcing is applied and zero thermal flux boundary conditions are considered so that only the internal diffusion due to gradients of temperature affects the evolution of the density in the fluid (cf. Figure 19). Given a non null thermal conductivity, the density converges towards a stationary and uniform solution. The evolution of density and temperature obtained with both Navier-Stokes-Fourier and Boussinesq is plotted in Fig. 20. The temperature converges rigorously towards an equilibrium temperature of T eq = T m = (T 0 + T 1 )/2 = 20 • C in the case of the Boussinesq model. For the Navier-Stokes-Fourier model, the equilibrium temperature is shifted below T m due to dilatation effects and is equal to T eq = 19.977 • C. Because of the diffusion effects, the equilibrium temperature is lower than T m . Note that for both models the density does not converge towards ρ m as the state law is not linear.

As in Fig. 18, Fig. 21 emphasizes the effects induced by the Boussinesq assumption over mass and volume variations.

Conclusion

In this work, we have proposed and analyzed a finite-volume scheme to solve the Navier-Stokes-Fourier equations, which describe free-surface variable density flows. With any flux consistent with the semi-discrete in time Euler system, the proposed scheme is wellbalanced and preserves the non-negativity of the water depth. In the case of a kinetic flux, a discrete entropy balance is proved for a flat topography. The numerical scheme is validated. The confrontation is made with results obtained with the simulation of the Boussinesq system. Notably, in a simple thermal diffusion case, the equilibrium temperature is not the same for the two systems.

A discrete entropy balance for a non-flat topography has yet to be obtained. Another challenge is the design and analysis of a numerical scheme for a system with a non-Newtonian rheology and in such a situation, we expect to be able to simulate complex interactions between the rheology terms and the temperature fluxes. Finally, simulations of the Navier-Stokes-Fourier system could be performed to investigate the propagation of internal waves in a stratified ocean.

ρgdz 1 =

 1 ?eq:NSF_2_p2? ∂(ρu) ∂t + ∇ x,y • (ρu ⊗ u) + ∂(ρuw) ∂z + ∇ x,y η z µ∆ x,y u + µ ∂ ∂z ∂u ∂z , (3) eq:NSF_3_p2 where U(t, x, y, z) = (u, v, w) T is the velocity, u = (u, v) T is the horizontal velocity vector and ρ is the density. The notation ∇ denotes ∇ = ∂ ∂x , ( ∂ ∂y , ( ∂ ∂z T , ∇ x,y corresponds to the projection of ∇ on the horizontal plane i.e. ∇ x,y = ∂ ∂x , ∂ ∂y T and the quantity |∇ x,y u| 2 means |∇ x,y u| 2 = (∇ x,y u) : (∇ x,y u) T .

2 +

 2 with e is the internal energy of the fluid governed by∂(ρe) ∂t + ∇ x,y • (ρeu) = η z ρgdz 1 ∇ x,y • u + µ|∇ x,y u| 2 + µ ∂u ∂z ∇.(λ∇T ).

Figure 1 :

 1 Figure 1: Flow domain with water height h(t, x, y), free surface η(t, x, y) and bottom z b (x, y).

  tations3d_rho

BC:neumann_bot or T b = T 0 b( 10 )

 010 BC:dirichlet_b and at the free surface λ ∂T ∂n s = φ T | η (11) BC:neumann_sur or T s = T 0 s (12) BC:dirichlet_s

  :exchange_te and the sum (total or partial) of the previous equations with G 1/2 = G N +1/2 = 0 (see (Boittin et al., 2020, Prop. 3.1)) gives Eqs. (28) and (31).

Figure 3 :

 3 Figure 3: The two functions λ → Q A 1 (λ) and λ → λ -N j=1 l j u j , each intersection of the two curves is an eigenvalue of A 1 ( Ũ).

  fig:pol_carac

Figure 4 :

 4 Figure 4: (a) Dual cell C i and (b) Boundary cell C i .

3. 4 . 5

 45 Properties of the numerical scheme op_num_scheme prop:num_flux Proposition 3.3 Consider a consistent numerical flux F for the homogeneous problem ( (35)-(38) i.e. Eq. (43) with S p (U, z b ) = 0) that preserves the non-negativity of the water depth h i under the corresponding CFL condition, then the finite volume scheme (58)-(62) (i ) preserves the non-negativity of the water depth,

Proposition 3 . 4

 34 Consider a consistent numerical flux F for the homogeneous problem (35)-(38) i.e. Eq. (43) with S p (U, z b ) = 0 that preserves nonnegativity of h i (t) and such that

Figure 5 :

 5 Figure 5: Analytical solution of prop. 5.1, 3D planar surface in a parabolic bowl: free surface at t = 0 (red), t = τ /4 (dark grey), t = τ /2 (blue), with the period τ defined by τ = 2π/ω.

Figure 6 :

 6 Figure 6: Numerical result of the parabolic bowl with variable density. Free surface and density contour in the slice plane (x, y=0, z) at initial time (left) and at time τ = 2π/ω with first order scheme (right).

Figure 8 :Figures 9 and 10

 810 Figure 8: Parabolic bowl: convergence of ρhu and ρhv in L 2 -norm towards the analytical solution, constant number of layers.

Figure 10 :

 10 Figure 10: Parabolic bowl: convergence of ρhu and ρhv in L 2 -norm towards the analytical solution, increasing number of layers.

Figure 11 :

 11 Figure 11: Parabolic bowl: error of ρh in L 2 -norm as a function of vertical and horizontal discretization for first order (a) and second order (b) numerical schemes.

ergence_plane 5 . 2

 52 Lock exchangeGravity currents triggered by lock-exchanges are encountered in many applications and their numerical simulation is a challenge. In this section we show the ability of our numerical scheme to properly simulate the propagation of lock-exchange induced density currents. The computed front position is compared to the experiments carried out by[START_REF] Adduce | Gravity Currents Produced by Lock Exchanges: Experiments and Simulations with a Two-Layer Shallow-Water Model with Entrainment[END_REF]. The results presented were obtained with the Navier-Stokes-Fourier model where the dissipation terms S µ,α due to the viscosity were neglected.

Figure 12 :

 12 Figure 12: Fluid domain of the lock-exchange test case

  {?} with ν the kinematic viscosity. Simulations have been carried out with a Grashof number of Gr = 2.53 × 10 8 on different meshes. The evolution of the density and the position of the front are presented in figures 13 and 14 respectively. The initial opening of the gate has been delayed (1s) to take into account the time of opening in the experiment. The figure14shows the convergence of the numerical scheme and we observe a good matching between the numerical simulation and the experimental data of[START_REF] Adduce | Gravity Currents Produced by Lock Exchanges: Experiments and Simulations with a Two-Layer Shallow-Water Model with Entrainment[END_REF].

Figure 13 :

 13 Figure 13: Lock exchange: computed density with the most refined mesh in the slice plane (x, y = 0, z) with Gr = 2.53 × 10 8 at times t = 3, 7, 9 and 11s (from top to bottom).

Figure 14 :

 14 Figure 14: Lock exchange: front position as a function of time for different meshes with comparison to Adduce et al. (2012) experimental results (where N t is the number of triangles and N is the number of layers).

  results_front

Figure 15 :

 15 Figure 15: Fluid domain of the diffusion test case with Dirichlet boundary condition at the bottom.n_diffusion_1

Figure 16 : 2 0 2 T

 1622 Figure 16: Temperature diffusion: evolution of the density in the slice plane (x, y=0, z) with the Navier-Stokes-Fourier model at time t = 0.03 (left) and t = 0.06 (right).

Figure 17 :

 17 Figure 17: Temperature diffusion: dimensionless temperature T as a function of z/h 0 at different times and comparison between analytical solution and numerical simulation with a number of layers equal to N = 20.

Figure 18 :

 18 Figure 18: Evolution of the mass ratio (m/m 0 ) and volume ratio (V /V 0 ) for the Navier-Stokes-Fourier and Boussinesq models.

  ass_vs_volume

Figure 19 : 2 0,Figure 20 :

 19220 Figure 19: Thermal equilibrium: fluid domain for the diffusion test case

Figure 21 :

 21 Figure 21: Thermal equilibrium: evolution of the mass ratio (m/m 0 ) and volume ratio (V /V 0 ) for the Navier-Stokes-Fourier and the Boussinesq models.

  ass_vs_volume

  69) eq:maxpple_hor

	We deal next with the vertical exchanges, see paragraph 3.4.3.
	n+1/2 N,i Note that the matrix H

  Parabolic bowl: convergence of h and ρh in L 2 -norm towards the analytical solution, constant number of layers.

					h							ρh
		ΫͷǤͷ						ΫͳǤͲ
		ΫǤͲ						ΫͳǤͷ
	L 2						L 2
	h	ΫǤͷ					ρh	ΫʹǤͲ
		ΫǤͲ	1 st					ΫʹǤͷ	1 st
			2 nd						2 nd
		ΫǤͷ	1 st 2 nd					Ϋ͵ǤͲ	1 st 2 nd
		ͲǤͲ	ͲǤʹ	ͲǤͶ	ͲǤ	ͲǤͺ	ͳǤͲ	ͲǤͲ	ͲǤʹ	ͲǤͶ	ͲǤ	ͲǤͺ	ͳǤͲ
				0.5 * log(l 0 /l i )					0.5 * log(l 0 /l i )
	Figure 7: ergence_cstNL
					ρhu					ρhv
								ΫͳǤͲ
		ΫͳǤͷ
								ΫͳǤͷ
		ΫʹǤͲ
	L 2						L 2	ΫʹǤͲ
	ρhu	ΫʹǤͷ					ρhv	ΫʹǤͷ
		Ϋ͵ǤͲ	1 st 2 nd					Ϋ͵ǤͲ	1 st 2 nd
		Ϋ͵Ǥͷ	1 st 2 nd					Ϋ͵Ǥͷ	1 st 2 nd
		ͲǤͲ	ͲǤʹ	ͲǤͶ	ͲǤ	ͲǤͺ	ͳǤͲ	ͲǤͲ	ͲǤʹ	ͲǤͶ	ͲǤ	ͲǤͺ	ͳǤͲ
				0.5 * log(l 0 /l i )					0.5 * log(l 0 /l i )

  Parabolic bowl: convergence of h and ρh in L 2 -norm towards the analytical solution, increasing number of layers.

							h						ρh
		ΫͷǤͷ						ΫͳǤͲ
		ΫǤͲ						ΫͳǤͷ
	h L 2	ΫǤͷ ΫǤͲ					ρh L 2	ΫʹǤͷ ΫʹǤͲ
		ΫǤͷ						Ϋ͵ǤͲ
		ΫͺǤͲ	1 st 2 nd					Ϋ͵Ǥͷ	1 st 2 nd
		ΫͺǤͷ	1 st 2 nd					ΫͶǤͲ	1 st 2 nd
			ͲǤͲͲ	ͲǤʹͷ	ͲǤͷͲ	ͲǤͷ	ͳǤͲͲ	ͳǤʹͷ	ͲǤͲͲ	ͲǤʹͷ	ͲǤͷͲ	ͲǤͷ	ͳǤͲͲ	ͳǤʹͷ
					log(l 0 /l i )					log(l 0 /l i )
	Figure 9: h_convergence
							ρhu						ρhv
			ΫͳǤͷ
									ΫͳǤͷ
			ΫʹǤͲ
									ΫʹǤͲ
			ΫʹǤͷ
	L 2						L 2	ΫʹǤͷ
	ρhu	Ϋ͵ǤͲ Ϋ͵Ǥͷ ΫͶǤͲ ΫͶǤͷ	1 st 2 nd 1 st 2 nd				ρhv	Ϋ͵ǤͲ Ϋ͵Ǥͷ ΫͶǤͲ ΫͶǤͷ	1 st 2 nd 1 st 2 nd
			ͲǤͲͲ	ͲǤʹͷ	ͲǤͷͲ	ͲǤͷ	ͳǤͲͲ	ͳǤʹͷ	ͲǤͲͲ	ͲǤʹͷ	ͲǤͷͲ	ͲǤͷ	ͳǤͲͲ	ͳǤʹͷ
					log(l 0 /l i )					log(l 0 /l i )

Acknowledgments

The authors acknowledge the Inria Project Lab "Algae in Silico" for its financial support. This research is also supported by the ERC SLIDEQUAKES ERC-CG-2013-PE10-617472.

and the vector containing the temperatures in layer α T α = (T α,1 , . . . , T α,Nx ) T .

The discretization of ∇ x,y •(µh

with the 4N x × 4N x block matrices

where the non-zero coefficients are given by

.

The ϕ i are the basis functions and δ i,j is the Kronecker symbol. The area A j is the area of the support of the test function ϕ j . The discretization of the terms S T,α due to temperature diffusion is similar. To discretize ρ (Tα) ρ 2 α cp S T,α , we propose the simplification

The term S µ,α ensures that the energy of the layer-averaged system is consistent, see [START_REF] Boittin | Low-Mach type approximation of the Navier-Stokes system with temperature and salinity for free surface flows[END_REF]. Indeed, the following identity holds

It is important to ensure the same consistency at the discrete level, i.e. an analogous identity should be verified at the discrete level. This gives guidelines for the discretization of S µ,α . A similar problem is studied in [START_REF] Grapsas | An unconditionally stable staggered pressure correction scheme for the compressible Navier-Stokes equations[END_REF] in the case of the compressible Navier-Stokes equations. However the numerical analysis of the scheme with the dissipation and diffusion terms has not been done yet. It will be performed in a future work.