Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface

In this paper we propose a stable and robust strategy to approximate the 3d incompressible hydrostatic Euler and Navier-Stokes systems with free surface. Compared to shallow water approximation of the Navier-Stokes system, the idea is to use a Galerkin type approximation of the velocity ﬁeld with piecewise constant basis functions in order to obtain an accurate description of the vertical proﬁle of the horizontal velocity. Such a strategy has several advantages. It allows


Introduction
In this paper we present layer-averaged Euler and Navier-Stokes models for the numerical simulation of incompressible free surface flows over variable topographies. We are mainly interested in applications to geophysical water flows such as tsunamis, lakes, rivers, estuarine waters, hazardous flows in the context either of advection dominant flows or of wave propagation.
The simulation of these flows requires stable, accurate, conservative schemes able to sharply resolve stratified flows, to handle efficiently complex topographies and free surface deformations, and to capture robustly wet/dry fronts. In addition, the application to realistic three-dimensional problems demands efficient methods with respect to computational cost. The present work is aimed at building a simulation tool endowed with these properties.
Due to computational issues associated with the free surface Navier-Stokes or Euler equations, the simulations of geophysical flows are often carried out with shallow water type models of reduced complexity. Indeed, for vertically averaged models such as the Saint-Venant system [12], efficient and robust numerical techniques (relaxation schemes [21], kinetic schemes [7,47], . . . ) are available and avoid to deal with moving meshes. In order to describe and simulate complex flows where the velocity field cannot be approximated by its vertical mean, multilayer models have been developed [3,5,8,20,22,30,31,46,53]. Unfortunately these models are physically relevant for non miscible fluids. In [9,10,27,35,36,49], some authors have proposed a simpler and more general formulation for multilayer model with mass exchanges between the layers. The obtained model has the form of a conservation law with source terms and presents remarkable differences with respect to classical models for non miscible fluids. In the multilayer approach with mass exchanges, the layer partition is merely a discretization artefact, and it is not physical. Therefore, the internal layer boundaries do not necessarily correspond to isopycnic surfaces. A critical distinguishing feature of our model is that it allows fluid circulation between layers. This changes dramatically the properties of the model and its ability to describe flow configurations that are crucial for the foreseen applications, such as recirculation zones.
Compared to previous works of some of the authors [9,10,27], that handled only the 2D configurations, this paper deals with the 3D case on unstructured meshes reinforcing the need of efficient numerical schemes. The key points of this paper are the following -A formulation of the 3D Navier-Stokes system under the form of a set of conservation laws with source terms on a fixed 2D domain. -Compared to previous works of some of the authors [9,10], we propose a complete kinetic interpretation of the model allowing to derive a robust and accurate numerical scheme. Notice that the kinetic interpretation is valid for the implicit treatment of the vertical exchange terms arising in the multilayer description. -Choosing a Newtonian rheology for the fluid, we propose energy-consistent -at the continuous and discrete levels -models extending previous results [27] in the 3D context. -Compared to 2D (x, z) situations, the 3D approximation with unstructured meshes raises new difficulties (boundary conditions, computational costs, upwinding, . . . ). The fact that the proposed strategy is also efficient in 3D proves -to some extent -its applicability. We give a complete description of all the ingredients of the numerical scheme. Even if some parts have been already published in 2D, the objective is to have a self-contained paper for 3D applications. -The numerical approximation of the 3D Navier-Stokes system is endowed with strong stability properties (consistency, well-balancing, positivity of the water depth, wet/dry interfaces treatment, . . . ). -Using academic examples, we prove the accuracy of the proposed numerical procedure especially convergence curves towards a 3D non-stationary analytical solution with wet-dry interfaces have been obtained (see Sect. 6.2.1).
Most of the numerical models in the literature for environmental stratified flows use finite difference or finite element schemes solving the free surface Navier-Stokes equations. We refer in particular to [32,39] and references therein for a partial review of these methods. Since the layer-averaged model has the form of a conservation law with source terms, we single out a finite volume scheme. Moreover, the kinetic interpretation of the continuous model leads to a kinetic solver endowed with strong stability properties (well-balancing, domain invariant, discrete entropy [11]). The viscous terms are discretized using a finite element approach. Considering various analytical solutions we emphasize the accuracy of the discrete model and we also show the applicability of the model to real geophysical situations. The numerical method is implemented in Freshkiss3D [2] and other various academic tests are documented on the web site.
The outline of the paper is as follows. In Section 2, we recall the incompressible and hydrostatic Navier-Stokes equations and the associated boundary conditions. The layer-averaged system obtained by a vertical discretization of the hydrostatic model is described in Section 3. The kinetic interpretation of the model is given in Section 4 allowing to derive a numerical scheme presented in Section 5. Numerical validations and application to real test cases are shown in Section 6.

The hydrostatic Navier-Stokes system
We consider the three-dimensional hydrostatic Navier-Stokes system [40] describing a free surface gravitational flow moving over a bottom topography z b (x, y). For free surface flows, the hydrostatic assumption consists in neglecting the vertical acceleration, see [24,25,37,43] for justifications of such hydrostatic models.
The incompressible and hydrostatic Navier-Stokes system consists in the model where U(t, x, y, z) = (u, v, w) T is the velocity, u(t, x, y, z) = (u, v) T is the horizontal velocity, σ = −pI d + µ∇ x,y u = −pI d + Σ is the total stress tensor, p is the fluid pressure, g represents the gravity acceleration and ρ 0 is the fluid density. The quantity ∇ denotes ∇ = ∂ ∂x , ∂ ∂y , ∂ ∂z T , ∇ x,y corresponds to the projection of ∇ on the horizontal plane i.e. ∇ x,y = ∂ ∂x , ∂ ∂y T . We assume a Newtonian fluid, µ is the viscosity coefficient and we will make use of ν = µ/ρ 0 . We consider a free surface flow (see Fig. 1a), 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. Due to the hydrostatic assumption in equation (2.3), the pressure gradient in equation (2.2) reduces to ρ 0 g∇ x,y η.

Bottom and free surface
Let n b and n s be the unit outward normals at the bottom and at the free surface respectively defined by (see Fig. 1a) , and n s = 1 1 + |∇ x,y η| 2 −∇ x,y η 1 .
The system (2.1)-(2.3) is completed with boundary conditions. On the bottom we prescribe an impermeability condition U.n b = 0, (2.4) whereas on the free surface, we impose the kinematic boundary condition ∂η ∂t + u(t, x, y, η).∇ x,y η − w(t, x, y, η) = 0. Concerning the dynamical boundary conditions, at the bottom we impose a friction condition given e.g. by a Navier law with κ a Navier coefficient. For some applications, one can choose κ = κ(h, u| b ). At the free surface, we impose the no stress condition whereũ = (u, 0) T , p a (t, x, y) and W (t, x, y) are two given quantities, p a (resp. W ) mimics the effects of the atmospheric pressure (resp. the wind blowing at the free surface) and t s is a given unit horizontal vector. Throughout the paper p a = 0 except in Section 6.2.2 where the effects of the atmospheric pressure is considered.

Fluid boundaries and solid walls
On solid walls, we prescribe an impermeability condition U.n = 0, (2.8) coupled with an homogeneous Neumann condition n being the outward normal to the considered wall.
In this paper we consider fluid boundaries on which we prescribe zero, one or two of the following conditions depending on the type of the flow (fluvial or torrential): Water level h + z b (x, y) given, flux hU given.

Energy balance
The smooth solutions of the system (2.1)-(2.7) satisfy the energy balance

The hydrostatic Euler system
In the case of an inviscid fluid, the system (2.1)-(2.7) consists in the incompressible and hydrostatic Euler equations with free surface and reads coupled with the two kinematic boundary conditions (2.4), (2.5) and p(t, x, y, η(t, x, y)) = 0. We recall the fundamental stability property related to the fact that the hydrostatic Euler system admits, for smooth solutions, an energy conservation that can be written under the form with E defined by (2.11).

The layer-averaged Euler system
The layer-averaging process for the 2D hydrostatic Euler and Navier-Stokes systems is precisely described in the paper [27] with a general rheology, the reader can refer to it. In the following, we present a Galerkin type approximation of the 3D Euler system also leading to a layer-averaged version of the Euler system, the obtained model reduces to [27] in the 2D context.
Using the notations (3.1), let us consider the space P N,t 0,h of piecewise constant functions defined by 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 for X ∈ (u, v, w). The three following propositions hold.
instead of (3.7) leads to a vanishing right hand side in equation (3.11). But the centered choice (3.12) does not allow to obtain an energy balance in the variable density case and does not give a maximum principle, at the discrete level, see [9]. Simple calculations show that any other choice than (3.7) or (3.12) leads to a non negative r.h.s. in (3.11). Indeed, let us consider the energy balance given in Proposition 3.2. In order to ensure either conservation or dissipation of the energy, it is necessary to ensure for any values of G α+1/2 , u α and u α+1 . In other words, one should have for any value of G α+1/2 , u α and u α+1 . And this gives the two possible choices (3.7) or (3.12) for u α+1/2 .
It is noticeable that, thanks to the kinematic boundary condition at each interface, the vertical velocity is no more a variable of the system (3.5). This is an advantage of this formulation over the hydrostatic model where the vertical velocity is needed in the momentum equation (2.2) and is deduced from the incompressibility condition (2.1). Even if the vertical velocity w no more appears in the model (3.4) and (3.5), it can be obtained as follows.
The quantities {w α } N α=1 are obtained only using a post-processing of the variables governing the system (3.4) and (3.5).
Proof of Proposition 3.1. Considering the divergence free condition (2.12), using the decomposition (3.3) and the space of test functions (3.2), we consider the quantity with G α±1/2 defined by The sum for α = 1, . . ., N of the above relations gives equation (3.4) where the kinematic boundary conditions (2.4) and (2.5) corresponding to leading, after simple computations, to equation (3.5).
Proof of Proposition 3.2. In order to obtain (3.8) we multiply (3.14) by g(h + z b ) − |u α | 2 /2 and (3.5) by u α then we perform simple manipulations. More precisely, the momentun equation along the x axis multiplied by u α gives Considering first the left hand side of the preceding equation excluding the pressure terms, we denote and using (3.14) we have Now we consider the contribution of the pressure terms over the energy balance i.e.
and it comes Performing the same manipulations over the momentum equation along y and adding the terms I u,α ,I v,α ,I p,u,α ,I p,v,α and (3.14) multiplied by g(h + z b ) − |u α | 2 /2 gives the result.
Proof of Proposition 3.4. Using the boundary condition (2.4), an integration from z b to z of the divergence free condition (2.1) easily gives Replacing formally in the above equation u (resp. w) by u N (resp. w N ) defined by (3.3) and performing an integration over the layer L 1 of the obtained relation yields i.e. w 1 = ∇ x,y .(z b u 1 ) − z 1 ∇ x,y .u 1 , corresponding to (3.13) for α = 1. A similar computation for the layers L 2 , . . ., L N proves the result (3.13) for α = 2, . . ., N . A more detailed version of this proof is given in [27].

The layer-averaged Navier-Stokes system
In Section 3.1, we have applied the layer-averaging to the Euler system, we now use the same process for the hydrostatic Navier-Stokes system. First, we consider the Navier-Stokes system (2.1)-(2.7) for a Newtonian fluid and then with a simplified rheology.

Complete model
The layer-averaging process applied to the Navier-Stokes system (2.1)-(2.7) leads to the following proposition.
Notice that in (3.23), we use the notation Proof of Proposition 3.5. The proof is given in Appendix A.
Remark 3.6. Notice that in the definition (3.21), since we consider viscous terms we use a centered approximation.

Simplified rheology
The viscous terms in the layer-averaged Navier-Stokes system are difficult to discretize especially when a discrete version of the energy balance has to be preserved. Hence, we propose a modified version of the model given in Proposition 3.5.
Proposition 3.7. The layer-averaged Navier-Stokes can be written under the form and
The proof of the energy balance (3.34) is similar to the one given in Proposition 3.5.
Remark 3.8. The layer-averaged Navier-Stokes system defined by (3.24) and (3.25) has the form We denote with F (U ) the fluxes of the conservative part, and with S e (U, ∂ t U, ∂ x U ) and S v,f (U ) the source terms, representing respectively the momentum exchanges and the viscous, wind and friction effects. The numerical scheme for the system (3.35) will be given in Section 5.

Kinetic description for the Euler system
In this section we give a kinetic interpretation of the system (3.4)-(3.8). The numerical scheme for the system (3.4), (3.5) and (3.13) will be deduced from the kinetic description.

Preliminaries
We begin this section by recalling the classical kinetic approach -used in [48] for example -for the 1D Saint-Venant system with the water depth h(t, x) ≥ 0, the water velocity u(t, x) ∈ R and a slowly varying topography z b (x). The kinetic Maxwellian is given by where U = (h, hu) T , ξ ∈ R and x + ≡ max(0, x) for any x ∈ R. It satisfies the following moment relations, These definitions allow us to obtain a kinetic representation of the Saint-Venant system.
The standard way to use Lemma 4.1 is to write a kinetic relaxation equation [7,16,17], like x, ξ)dξ, and > 0 is a relaxation time. In the limit → 0 we recover formally the formulation (4.4), (4.5). We refer to [16] for general considerations on such kinetic relaxation models without topography, the case with topography being introduced in [48]. Note that the notion of kinetic representation as (4.4), (4.5) differs from the so called kinetic formulations where a large set of entropies is involved, see [47]. For systems of conservation laws, these kinetic formulations include non-advective terms that prevent from writing down simple approximations. In general, kinetic relaxation approximations can be compatible with just a single entropy. Nevertheless this is enough for proving the convergence as ε → 0, see [15].
To build the Gibbs equilibria, we choose the function This choice corresponds to the 2D version of the kinetic maxwellian used in 1D (see Rem. 4.2) and we have and where (ξ, γ) ∈ R 2 . In other words, we have Starting from the 2D maxwellian in the single layer case i.e.
and computing its integral w.r.t. the variable γ yields that is exactly the expression (4.2).
The quantity M α satisfies the following moment relations The interest of the function χ 0 and hence the particular form (4.8) lies in its link with a kinetic entropy. Consider the kinetic entropy Then one can check the relations Let us introduce the Gibbs equilibria N α+1/2 defined by for α = 0, . . ., N by where G α+1/2 is defined by (3.6) and u α+1/2 ,v α+1/2 are given by (3.7). The quantity N α+1/2 satisfies the following moment relations Notice that from (3.6), we can give a kinetic interpretation on the exchange terms under the form for α = 1, . . ., N . Then we have the two following results.
Proof of Proposition 4.4. The proof is obtained multiplying (4.17) by So finally, equation (4.17) multiplied by H α (M α , ξ, γ, z b ) gives It remains to calculate the sum of the preceding relations from α = 1, . . ., N and to integrate the obtained relation in ξ, γ over R 2 that completes the proof.
for i, j = 1, . . ., N with δ i,j the Kronecker symbol. Then, using Proposition 4.3, we can write Hence, using the above notations, the layer-averaged Euler system (3.4) and (3.5) can be written under the form

Numerical scheme
The numerical scheme for the model (3.35) proposed in this section extends the results presented by some of the authors in [4,7,9,10]. Compared to these previous results, it has the following advantages -it gives a 3D approximation of the Navier-Stokes system whereas 2D situations (x, y) and (x, z) where considered in [4,9,10], -the implicit treatment of the vertical exchanges terms gives a bounded CFL condition even when the water depth vanishes, -the kinetic interpretation, on which is based the numerical scheme, is also valid for the vertical exchange terms -that was not the case in [9,10] -and allows to derive a robust and accurate numerical scheme, -the numerical approximation of the system given in (3.35) is endowed with strong stability properties (wellbalancing, positivity of the water depth, . . . ), -convergence curves towards a 3D non-stationary analytical solution with wet-dry interfaces have been obtained (see Sect. 6.2.1).
First, we focus on the Euler part of the system (3.35) then in Section 5.6, a numerical scheme for the viscous terms is proposed.
Notice that, as a consequence of the layer-averaged discretization, the system (3.35) and the Boltzmann type equation (4.17) are 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.

Semi-discrete (in time) scheme
We consider discrete times t n with t n+1 = t n + ∆t n . For the time discretisation of the layer-averaged Navier-Stokes system (3.35) we adopt the following scheme where the superscript n has been omitted and the integer p = 0, 1/2, 1 will be precised below. Using the expressions (3.24) and (3.25) for the layer averaged model, the semi-discrete in time scheme (5.1) writes for α = 1, . . ., N . The vertical velocities {w α } N α=1 are defined by (3.13). The first two equations (5.2) and (5.3) consist in an explicit time scheme where the horizontal fluxes and the topography source term are taken into account whereas in equation (5.4) an implicit treatment of the exchange terms between layers is proposed. The implicit part of the scheme requires to solve a linear problem (see Lem. 5.1) but, on the contrary of previous work of some of the authors [10], it implies that the CFL condition (5.31) no more depends on the exchange terms.
When ν = κ = 0 in equation (5.4), equations (5.2)-(5.4) correspond to the layer-averaged of the Euler system. The choice p = 1 (resp. p = 1/2) in equation (5.4) corresponds to an implicit (resp. semi-implicit) treatment of the viscous and friction terms whereas the choice p = 0 implies an explicit treatment and requires a CFL condition. Notice that -the implicit or semi-implicit treatment of the viscous terms requires to solve a system that is linear since h n+1 is known from the mass conservation equation, -the explicit discretisation is simpler to implement but implies a restrictive CFL when the spatial discretisation (vertical and horizontal) becomes fine. In the context of geophysical flows, this limitation is rarely severe. In particular, the CFL condition induced by the explicit treatment of the viscous terms is less restrictive than the one coming from the advection part (see Eq. (5.31)) in all the numerical tests of Section 6.

Space discretization
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.
For the space discretization of the system (5.2)-(5.4), we use a finite volume technique for the Euler partthat is described below -and a finite element approach -P 1 on T h -for the viscous part that is described in Section 5.6.

Finite volume formalism for the Euler part
In this section and in Section 5.4, we propose a space discretization for the model (5.2)-(5.4) without the viscous and friction terms i.e. the system completed with (5.5).
We recall now the general formalism of finite volumes on unstructured meshes. 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. 2): 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. 2) 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 and z b (x, y) as We will also use the notation with U α defined by (4.9). A finite volume scheme for solving the system (5.6) and (5.7) is a formula of the form where using the notations of (5.1) Here we consider first-order explicit schemes where and and for the boundary nodes Relation (5.10) tells how to compute the values U n+1/2 i knowing U i and discretized values z b,i of the topography. Following (5.11), the term F i,j in (5.10) 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 , z b,i −z b,j , n i,j ) ∈ R 2N +1 are the numerical fluxes, see [19].
In the next section we define F(U i , U j , z b,i − z b,j , n i,j ) using the kinetic interpretation of the system. The computation of the value U i,e , which denotes a value outside C i (see Fig. 2b), 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 Section 5.7. Notice that we assume a flat topography on the boundaries i.e. z b,i = z b,i,e .

Discrete kinetic equation
The choice of a kinetic scheme is motivated by several arguments. First, the kinetic interpretation is a suitable starting point for building a stable numerical scheme. We will prove in Section 5.4 that the proposed kinetic scheme preserves positivity of the water depth and ensures a discrete local maximum principle for a tracer concentration (temperature, salinity, . . .). Second, the construction of the kinetic scheme does not need the computation of the system eigenvalues. This point is very important here since these eigenvalues are not available in explicit analytical form, and they are hardly accessible even numerically. Furthermore, as previously mentioned, hyperbolicity of the multilayer model may not hold, and the kinetic scheme allows overcoming this difficulty.

Without topography
In a first step we consider a situation with flat bottom. Following Proposition 4.1, the model (3.4) and (3.5) reduces, for each layer, to a classical Saint-Venant system with exchange terms and its kinetic interpretation (see Eq. (4.17)) is given by with the notations defined in Section 4.2.
Let C i be a cell, see Fig. 2. The integral over C i of the convective part of the kinetic equation (5.15) gives with M α,i = M (U α,i , ξ, γ), n i,j being the outward normal to the cell C i . The quantity M α,i,j is defined by the classical kinetic upwinding Therefore, the kinetic scheme applied for equation (5.15) is given by Notice that the previous definition is consistent with (4.14). From (4.16), we get By analogy with the computations in (4.11), we can recover the macroscopic quantities U n+1 α,i at time t n+1 by integration of the relation (5.18) The scheme (5.17) and the definition (5.20) allow to complete the definition of the macroscopic scheme (5.10), (5.13) and (5.14) with the numerical flux given by the flux vector splitting formula [17] where I N is the identity matrix of size N and G N,i is defined by Hence, the resolution of the discrete kinetic equation ( Remark 5.2. Compared to an explicit treatment of the vertical exchanges terms as presented in [9,10], the implicit scheme (5.18) requires to invert for each cell a small matrix whose size corresponds to the number of layers. Depending on the type of simulation carried out, it increases the computational costs e.g. for the tsunami simulation the difference is around 20%.
But it is worth noticing that the explicit treatment of the vertical exchanges terms can lead to severe constraints on the CFL condition since the quantity is not bounded, see Proposition 5.2 from [10]. (ii) Denoting G d N,i (resp. G nd N,i ) the diagonal (resp. non diagonal) part of G N,i we can write where all the entries of the matrix J N,i = (I N + ∆tG d N,i ) −1 (−∆tG nd N,i ), are non negative and less than 1. And hence, we can write proving all the entries of (I N + ∆tG N,i ) −1 are non negative. (iii) Let us consider the vector 1 whose entries are all equal to 1. Since we have we also have 1 = (I N + ∆tG N,i ) −t 1. Now let T be a vector whose entries {T α } 1≤α≤N are non negative, then that completes the proof.

With topography
The hydrostatic reconstruction scheme (HR scheme for short) for the Saint-Venant system has been introduced in [6] in the 1D case and described in 2D for unstructured meshes in [4]. The HR in the context of the kinetic description for the Saint-Venant system has been studied in [7].
In order to take into account the topography source and to preserve relevant equilibria, the HR leads to a modified version of (5.10) under the form where We would like here to propose a kinetic interpretation of the HR scheme, which means to interpret the above numerical fluxes as averages with respect to the kinetic variables of a scheme written on a kinetic function f . More precisely, we would like to approximate the solution to (4.17) by a kinetic scheme such that the associated macroscopic scheme is exactly (5.22) and (5.23) with homogeneous numerical flux F given by (5.21). We denote M * α,i,j = M (U * α,i,j , ξ, γ) for any α = 1, . . ., N and we consider the scheme For the exchange terms, by analogy with (5.19) we define and using (4.16) we get It is easy to see that in the previous formula, we have the moment relations This means equivalently that |u α,i | + |v α,i | + √ 2gh i ≤ v m . We consider a CFL condition strictly less than one, where σ i = ∆t n j∈Ki L i,j /|C i |, and β is a given constant. Then the following proposition holds.

Macroscopic scheme
The numerical scheme for the system (5.6)-(5.8) is given by (5.20), (5.25), (5.26) and requires to calculate fluxes having the form with M given by (4.10), n x and n y being the components of a normal unit vector n. Defining the change of variables we can write where χ 0 is defined by (4.7). A second change of variables y 1 = n x z 1 + n y z 2 , y 2 = n x z 2 − n y z 1 ,ũ = n x u + n y v gives

The discrete layer-averaged Navier-Stokes system
In this section, we detail the space discretization of the viscous terms. Several expressions have been obtained for the viscous terms, see Section 3.2. In this section, we give a numerical scheme for the model (5.4), rewriting it under the form with the definitions (3.26), (3.27), (3.32), (3.33) and (3.30) for T 0 α , Γ α±1/2 . It remains to give a fully discrete scheme for the viscous and friction terms {S v,f,α }.
The discretization of (5.34) is done using a finite element/finite difference approximation obtained as follows. We depart form the triangulation defined in Section (5.2) and we use the cells values of the variables -inherited from the finite volume framework -to define a P 1 approximation of the variables.
Notice that, compared to the advection and pressure terms, the discretization of the viscous terms raises less difficulties and we propose a stable scheme that will be extended to more general rheology terms [27] and more completely analyzed in a forthcoming paper.
Using a classical P 1 finite element type approximation with mass lumping of equation (5.34), we get with the matrices where ϕ i , ϕ j are the basis functions. We have presented an explicit in time version of (5.35) that is stable under a classical CFL condition. An implicit or semi-implicit version of (5.35) can also be used.
The main purpose of this paper is to propose a stable and robust numerical approximation of the incompressible Euler system with free surface. Voluntarily, we give few details concerning the numerical approximation of the dissipative terms: -the viscous and friction terms are dissipative and hence a reasonable approximation leads to a stable numerical scheme. -In this paper, we consider a simplified Newtonian rheology for the fluid, the numerical approximation of the general (layer-averaged) rheology [27] will be studied in a forthcoming paper.

Boundary conditions
The contents of this section slightly differ from previous works of one of the authors [26] and valid for the classical Saint-Venant system. First, we focus on the boundary conditions for the layer-averaged Euler system i.e. the system (5.2) and (5.3) for ν = 0, κ = 0 and then for the viscous part.

Layer-averaged Euler system
In this section we detail the computation of the boundary flux F(U i , U e,i , n i ) appearing in (5.10), (5.13) and (5.14). The variable U n i,e can be interpreted as an approximation of the solution in a ghost cell adjacent to the boundary. As before we introduce the vector Solid wall. If we consider a node i 0 belonging to a solid wall, we prescribe a slip condition written u α · n i0 = 0, (5.37) for α = 1, . . ., N . We assume the continuity of the water depth h i0,e = h i0 and of the tangential component of velocity.
From (5.33) with (5.37) we obtain .t i0 = 0, for α = 1, . . ., N for a vector t i0 orthogonal to n i0 . The condition (5.37) is therefore prescribed weakly but a posteriori, in order to be sure that (5.37) is satisfied, we can apply Fluid boundary. Even if the considered model is more complex than the Shallow water system, we can consider that the type of the flow depends, for each layer, on the value of the Froud number F r α = |u α |/ √ gh, a flow is said torrential, for |u α | > √ gh and fluvial, for |u α | < √ gh. Generally, for the fluid boundaries, the conditions prescribed by the user depend on the type of the flow defined by this criterion.
We have also to notice that with n the outward unit normal to the boundary edge, an inflow boundary corresponds to u α · n < 0 and an outflow one to u α · n > 0.
We will treat the following cases: for a fluvial flow boundary, we distinguish the cases where the flux or the water depth are given, while for a torrential flow we distinguish the inflow or outflow boundaries.
Fluvial boundary. Flux given. We consider first a fluvial boundary, so we assume that |u α | < gh. (5.38) If for each layer, the flux q g,α is given, then we wish to impose with n i .t i = 0. The value q g,α depends on the value of the prescribed flux along the vertical axis. If one directly imposes (5.39), it leads to instabilities (especially because the numerical values are not necessarily in the regime of validity of this condition). We propose to discretize it in a weak form. We denote If a 1 ≥ 0, we prescribe F h (U e,α,i ) = 0, F hu (U e,α,i ) = 0, and F hv (U e,α,i ) = 0.
If a 1 < 0, we have to write a third equation to be able to compute the three components of U e,α and by analogy with what is done for the Saint-Venant system -where the Riemann invariant related to the outgoing characteristic is preserved -we assume the quantity u α .n is constant though the interface, i.e.
In practice, we use a Newton-Raphson algorithm to solve an equivalent form of equation (5.46), namely Once the above equation has been solved, from (5.42) to (5.43) we deduce Remark 5.5. Notice that in the procedure proposed to calculate u e,α,i , h e,i , even if h e,i represents a total water depth, a different value of h e,i is calculated for each layer α. h e,i is only used to ensure (5.39).
Fluvial boundary. Water depth given. We verify that the flow is actually fluvial, i.e.
(u α,i .n i − gh i )(u α,i .n i + gh i ) ≤ 0. (5.47) Since the water depth is given, we write h e,i = h g,i . (5.48) We assume the continuity of the tangential component with t i .n i = 0. To define completely u e,α,i , we assume, as in the previous case, that the Riemann invariant is constant along the outgoing characteristic (5.41), so we obtain Sometimes it appears that the numerical values do not satisfy the condition (5.47), then the flow is in fact torrential and -if u α,i .n i > 0 , the condition (5.48) cannot be satisfied (see Sect. 6.2.4), -if u α,i .n i < 0 , one condition is missing and we prescribe u e,α,i .n i = u α,i .n i . Torrential inflow boundary. For a torrential inflow boundary we assume that the water depth and the flux are given, then we prescribe In this case we have to compute (hu) e,α,i .n i or u e,α,i .n i . We consider an inflow boundary, so (hu) g,α,i .n i < 0 therefore using the notation (5.40) we have a 1 < 0. By analogy with the previous section we denote As in the section entitled Flux given, the above equation has a unique solution m < 2 for a 1 < 0.
Torrential outflow boundary. In the case of a torrential outflow boundary, we do not prescribe any condition. We assume that the two Riemann invariants are constant along the outgoing characteristics leading to u e,α,i .n i − 2 gh e,i = u α,i .n i − 2 gh i , u e,α,i .n i + 2 gh e,i = u α,i .n i + 2 gh i , and we deduce h e,i = h i , u e,α,i .n i = u α,i .n i . We assume that we also have (hu) e,α,i .t i = (hu) α,i .t i .

Layer-averaged Navier-Stokes system
Because of the fractional step we use, the boundary conditions for the layer-averaged Euler system are, to some extent, independent from the one used for the rheology terms.
For the resolution of equation (5.34), boundary conditions associated with the operator have to be specified and usually we prescribe homogeneous Neumann boundary conditions (corresponding to an imposed stress). Of course, in particular cases, Dirichlet or Robin type boundary conditions can also be considered.

Toward second order schemes
In order to improve the accuracy of the results the first-order scheme defined in Sections 5.3-5.5 can be extended to a formally second-order one using a MUSCL like extension (see [52]).

Second order reconstruction for the layer-averaged Euler system
In the definition of the flux (5.23), we replace the piecewise constant values U i,j , U j,i by more accurate reconstructions deduced from piecewise linear approximations, namely the valuesŨ i,j ,Ũ j,i reconstructed on both sides of the interface. The reconstruction procedure is similar to the one used and described in Section 5.1 from [4].
The second order reconstruction is only applied for the horizontal fluxes. For the exchange terms along the vertical axis involving the quantities G α±1/2 , we keep the first order approximation. Despite this, we recover over the simulations (see Sects. 6.1 and 6.2) a second order type convergence curve. For this reason, we call this reconstruction "second order".

Modified Heun scheme
The explicit time scheme (5.2) and (5.3) used in the previous paragraphs corresponds to a first order explicit Euler scheme. The second-order accuracy in time is usually recovered by the Heun method [18] that is a slight modification of the second order Runge-Kutta method. More precisely, for a dynamical system written under the form ∂y ∂t = f (y), (5.51) the Heun scheme consists in defining y n+1 by y n+1 = y(t n + ∆t n ) = y n +ỹ n+2 2 , (5.52) withỹ n+1 = y n + ∆t n f (y n , t n ),ỹ n+2 =ỹ n+1 + ∆t n f (ỹ n+1 , t n+1 ). To overcome this difficulty, we propose an improvement of the Heun scheme Proposition 5.6. The scheme defined by y n+1 = (1 − γ)y n + γỹ n+2 with y n+1 = y n + ∆t n 1 f (y n ),ỹ n+2 =ỹ n+1 + ∆t n 2 f (ỹ n+1 ), and ∆t n = 2∆t n 1 ∆t n is second order and compatible with a CFL constraint. Since γ ≥ 0, y n+1 is a convex combination of y n and y n+2 so the scheme preserves the positivity. For the previous relations ∆t n 1 and ∆t n 2 respectively satisfy the CFL conditions associated with y n andỹ n+1 .

Numerical applications
In this section, we use the numerical scheme to simulate several test cases: analytical solutions or in situ measurements, stationary or non-stationary solutions, for the Euler and Navier-Stokes systems. The obtained results emphasize the accuracy of the numerical procedure in a wide range of typical applications and its applicability to a real tsunami case. We also propose simulations of the hydrodynamic regime in a raceway agitated by a paddlewheel and confront the results with experimental measurements.
The numerical simulations presented in this section have been obtained with the code Freshkiss3D [2] where the numerical scheme presented in this paper is implemented.

Stationary analytical solution
First, we compare our numerical model with stationary analytical solutions for the free surface Euler system proposed by some of the authors in [23].
We consider as geometrical domain a channel (x, y) ∈ [0, x max ] × [0, 2]. The analytical solution given in Proposition 3.1 from [23] and defined by with α = 1 m 2 s −1 , β = 1 m −1 , z b = cst, x max = 20 m and h 0 (x, y) = 1 2 is a stationary regular analytical solution of the incompressible and hydrostatic Euler system with free surface (2.12), (2.13), (2.4) and (2.5) with p a = 0. In order to obtain the simulated solution, we consider the topography defined by (6.1), (6.3) and we impose the following boundary conditions -solid wall for the two boundaries y = 0 m and y = 2 m, -given water depth h 0 (x max , y) at x = x max = 20 m, -given flux defined by (6.2) at x = 0 m.
We have performed the simulations for several unstructured meshes having 290 nodes and 2 layers, 597 nodes and 4 layers, 1010 nodes and 8 layers, 2112 nodes and 17 layers, see Remark 6.1.
Remark 6.1. In each case where a convergence curve towards an analytical solution is presented, we have proceeded as follows. First, we choose a sequence of unstructured meshes for the considered horizontal geometrical domain. Then the number of layers is adapted so that each 3D element of the mesh can be approximatively considered as a regular polyhedron.  On Figure 3a, we have depicted the features of the analytical solution we use for the convergence test, it clearly appears on Figure 3a that the velocity profile of the chosen analytical solution varies along the z axis. The L 2 errors for the convergence test and the corresponding convergence rate are given in Table 1. Figure 3b gives the convergence curve towards the analytical solution i.e. the log(L 2 -error) of the water depth -at time T = 300 s when the stationary regime is reached -versus log(h a0 /h a ) for the first and second-order schemes and they are compared to the theoretical order (we denote by h a the average edge length and h a0 the average edge length of the coarser mesh).
Remark 6.2. Following the results given in Section 3.4 from [23], it is possible to obtain stationary analytical solutions with discontinuities for the Euler system. In this case, the unknowns are not given by algebraic expressions but are obtained through the resolution of an ODE involving only the water depth h.
The numerical scheme has been used in the context of such a discontinuous analytical solution. As planned, for the first and second order schemes, we recover a first order convergence of the simulated solution towards the analytical one because of the discontinuity of the reference solution.

Non-stationary analytical solutions
In a recent paper [28], some of the authors have proposed time-dependent 3D analytical solutions for the Euler and Navier-Stokes equations, some of them concern hydrostatic models. We confront our numerical scheme to these situations where analytical solutions are available.

Radially-symmetrical parabolic bowl
The Thacker' analytical solution [50], corresponds to a periodic oscillation in a parabolic bowl. In [28] an extension of the Thacker' radially-symmetrical solution to the situation where the velocity field depends on the vertical coordinate is proposed. This means the proposed solution, described hereafter in Proposition 6.3, is analytical for the 3D incompressible hydrostatic Euler system but does not correspond to a shallow water regime. Proposition 6.3. For some t 0 ∈ R, (α, β, γ) ∈ R 3 + * such that γ < 1 let us consider the functions h, u, v, w, p defined for t ≥ t 0 by h(t, x, y) = max 0, 1 , with ω = √ 4αg, r = x 2 + y 2 and with a bottom topography defined by 8) and the function f given by c being a negative constant such that c ≤ 4g 2 /(γ − 1). From equation (2.1), the vertical velocity w can be expressed under the form Then h, u, v, w, p as defined previously satisfy the 3D hydrostatic Euler system (2.12) and (2.13) completed with (2.4) and (2.5).
The geometrical domain is defined by (x, y) ∈ [−L/2, L/2] 2 and the chosen parameters are α = 2, β = 1, γ = 0.3, c = −1, L = 1, the considered analytical solution is depicted on Figure 4. The initial conditions correspond to (6.4)-(6.7) at time t = t 0 = 0 s. In order to evaluate the convergence rate of the simulated solution h sim towards the analytical one h anal , we have performed a convergence test, the errors and convergence rates appear over Table 2. The five unstructured meshes we have considered have respectively 1273 nodes and a single layer, 11 104 nodes and 6 layers, 30 441 nodes and 15 layers, 59 473 nodes and 30 layers and 98 137 nodes and 50 layers, see Remark 6.1. We have plotted (see Fig. 5) the log(L 2 -error) over the water depth at time T = 2π/ω s versus log(h a0 /h a ) for the first and second-order schemes and they are compared to the theoretical order.
The analytical solution of Proposition 6.3 is non stationary and hence, the errors due to the time scheme are combined with the one induced by the space discretization. Moreover this test case has a lot of wet/dry interfaces where the second order reconstruction in space cannot be applied. This can explained the differences between the theoretical and observed slopes for the convergence tests of the second order schemes, see Table 2 and Figure 5.

Draining of a tank
Considering the Navier-Stokes system (2.1)-(2.3) completed with the boundary conditions (2.5)-(2.8), the following proposition holds, see [28] for more details about the proposed analytical solution.
Then h, u, v, w, p as defined previously satisfy the 3D hydrostatic Navier-Stokes system ( Figure 6b gives the convergence curve towards the analytical solution i.e. the log(L 2 -error) of the water depth -at time T = 1 s -versus log(h a0 /h a ) for the first and second-order schemes and they are compared to the theoretical order (we denote by h a the average edge length and h a0 the average edge length of the coarser mesh). Notice that in this test case, the errors due to the space and time discretization are combined, this explains why the theoretical orders of convergence are not exactly obtained. Moreover, the boundary conditions (inflow prescribed) play an important role and since their numerical treatment is only at the first order in space, this also explains the difference between the theoretical and observed orders of convergence. With the mesh having 2781 nodes, we have tested the influence of the number of layers, see Figure 6c. When the numbers of layers increase, we recover the analytical velocity profile.

Simulation of a tsunami
In this section, we test our discrete model in the case of a real tsunami propagation for which field measurements are available (free surface variations recorded by buoys). Even if in such cases, involving long wave propagation, 2D shallow water models can be used instead of 3D description, and we test here the capacity of our model and of the numerical procedure to handle this complex situation.  Simulation of tsunami waves generated by earthquakes is very important in Earth science for hazard assessment and for recovering earthquakes characteristics. Indeed, tsunami waves can be analyzed to recover the earthquake source that generated the tsunami and are now classically used in joint inversion methods. It has been shown that tsunami waves provide strong constraints on the spatial distribution of the source, especially in the case of shallow slip [38]. In some cases, far-field tsunami gauges may help constrain the earthquake source process even though they are affected by the compressibility of the water column and of the Earth [38,54].
The 2014/04/01 Iquique earthquake struck off the coast of Chile at 20:46 local time (23:46 UTC), with a moment magnitude of 8.1. The epicenter of the earthquake was approximately 95 kilometers (59 mi) northwest of Iquique, as shown in Figure 7.
We have carried out simulations of the tsunami induced by the earthquake using -a topography obtained from the National Oceanic and Atmospheric Administration (NOAA, [13]) using the ETOPO1 data (1 arcmin global relief model), -an unstructured mesh whose dimensions -a square of 2224.2 km 2 -correspond to the domain covered by Figure 7, -a source corresponding to the seafloor displacement induced by the earthquake (Fig. 7). This source is obtained by computing the 3D final displacements of the seafloor generated by the earthquake coseismic slip. This coseismic slip has been itself retrieved by inversion of numerous geodetic and seismic data, according to the model determined by Vallée et al. [51]. The source is activated at time t 0 , just after the earthquake occurrence (t 0 is here 2014/04/01, 23h47mn25s) We did not consider here the Coriolis force, the tides and the ocean currents. The results shown in Figure 8 have been obtained with a mesh containing 545 821 nodes and 5 layers (computation time was 35 min with a Mac book air 1.7 GHz Intel core i7). We compare the numerical solutions -provided by the first order scheme (space and time) and the second order scheme (space and time) -with the DART measurements (obtained from the NOAA website http://www.ndbc.noaa.gov/dart.shtml). A series of simulations have been performed using several meshes and we present "converged" results in the sense that a finer mesh would give the same results. This is illustrated in Figure 8d, where we plot the simulation results obtained with three meshes having respectively 311 687 nodes (coarse mesh), 545 821 nodes (fine mesh), and 985 327 nodes (very fine mesh): the curves corresponding to the fine (cyan) and very fine (blue) curves are very similar. Figures 8a-8c shows that the second order scheme significantly improves the results both for the amplitude and phase of the water waves. The second order scheme is able to very accurately reproduce the shape of the first wave at the closest DART buoy 32 401, located at 287 km from the epicenter. The two following peaks in the waveform are quite well reproduced up to about 1.755 × 10 5 s. This is also the case at the DART buoy 32 402, located 853 km from the epicenter. The arrival time of the first wave is very well reproduced at the three DART buoys, slightly better with the second order scheme. At the most distant buoy 32 412 (1650 km from the source), the first order scheme is not able to reproduce the recorded wave. The second order scheme reproduces the first wave quite well but not the rest of the waveform, possibly due to Earth curvature effects that are not taken into account here. Globally, the low frequency content of the signals is better explained by the model than the high frequency fluctuations. These high frequency fluctuations may be related to effects not accounted for here, such as spatio-temporal heterogeneity of the real source, small wavelength fluctuation of the topography, and possibly non-hydrostatic effects [1,29].

Monai valley benchmark
In 2004, as part of a workshop organized by the US National Science Foundation, an experiment has been set up, reproducing the impact of a tsunami wave on the shore of the Okushiri island, in the Monai village area. The objective of the experiment was to provide a set of well organized data, reproducing the 1993 tsunami event, to validate numerical codes for Tsunami simulation [41,44].
This test case has been simulated by various numerical tools, it is well suited to test the numerical treatment of wet/dry interfaces. We reproduce hereafter the results obtained with Freshkiss3D [2].
The geometrical domain of 5.448 m × 3.402 m is depicted over Figure 9 where the bathymetry and the free surface elevation at time t = 16 s are presented. The topography, the input wave and the time series of surface elevation at different gauges are available here [45]. For the three gauges located respectively at (x = 4.521, y = 1.196), (x = 4.521, y = 1.696), and (x = 4.521, y = 2.196), we compare the simulation results and the experimental data, see Figure 10. The unstructured mesh used has 8 layers and 35 000 nodes for the horizontal mesh (corresponding to a mean edge length of 0.028 m), notice that it is significantly coarser than the recommended grid sizes of ∆x = ∆y = 0.014 m. Since the test case corresponds to a shallow water flow, the results we obtain are similar to those provided with the resolution of the classical Saint-Venant system, see (Clawpack [33], Hysea [42], . . . ).

Hydrodynamics in a raceway
In this section, we test our model in the case of raceways used as High Rate Algal Ponds developed to produce microalgae biomass [34] or treat wastewater as enhanced stabilization ponds. Figure 10. Comparison between the free surface variations recorded by the 3 gauges and the corresponding simulations (2nd order schemes in space and time).

Experimental measurements
Raceways are annular shaped ponds where the water is mixed with a paddlewheel, see Figure 11a. In this experiment, the raceway has a length of 4.2 m and a width of 1.2 m with perfect circular shape at each extremity. The height of water is 0.5 m for a total volume of water circulated of 2.37 m 3 , see Figure 11b. The paddlewheel has a diameter of 1.36 m and a width of 0.5 m with 6 blades equally distributed. Based on plastic material, blades are reinforced by two circular plates at each side with holes in between each blade in order to avoid air trapping. The paddlewheel was placed at the beginning of the straight part of the raceway (rotating axe at 0.66 m after the end of the curve) and a depth of 0.2 m leaving 0.3 m between the downiest part of the paddlewheel and the bottom of the raceway. The rotation frequency of the paddlewheel is maintained in order to obtain a speed at the circumference of 0.6 m s −1 .
In order to measure the water velocities at different localized points of the pond, a correlation wedge flow sensor from NIVUS (model POA-V2XXK) was used. It measures continuously the speed in a window of 5 • (inclination of the measure is 40 • ) of 16 layers from the bottom until a maximum height of 1 m. The accuracy alleged by the supplier is 0.5% for speed between 0.05 to 0.5 m s −1 and 1% over (max 6 m s −1 ). Hence, for the points A in , A out , B in , B out , C in , C mid , C out , P in , P out , Q in , Q out , R in and R out depicted over Figure 11b and having, in the plane (O, x, y), the coordinates (given in meter) Such experimental measurements allow to describe finely the hydrodynamic regime within the raceway. The measured average horizontal velocity is 0.21 m s −1 . A gradient of velocity can be observed from the bottom to the top of water as well as from the inner to the outer border. A dead zone is observed at 2.5 m close to the outer border (B out ). Profile is much more linear in the second straight part (return) with an acceleration of water close to the outer border (P out ) and a deceleration in the inner border with a dead zone after the curve close to the inner border (P in ). Finally, an ascendant flow is observed at the beginning of the second straight part (P out ) with high speed coming from the bottom to the top of the profile along the outer border (Q out ).

Simulation results
The simulation has been carried out with a mesh having 3330 nodes and 25 layers. Starting from a raceway at rest, after 50 s a stationary regime is reached -especially far from the paddlewheel. The viscosity of the fluid is ν = 0.005 m 2 s −1 and the bottom friction κ = 0.002 m 2 s −1 . The paddlewheel and its modeling in the multilayer framework is described in [14] and not presented here.
A global view of the hydrodynamics in the raceway at time T = 50 s is given over Figure 12 and the comparisons between the simulation results and the experimental measurements are given over Figure 13. Figure 13 enables to formulate three main comments -the velocity field in the raceway is really 3D in the sense that two closed points can have very different velocity fields, see e.g. points A out and A in or P out and P in . -The results are in good agreement with the measurements.
-At some points of the raceway (A in near the paddlewheel or C out near the turn), the non-hydrostatic effects can be significant and this can explain the discrepancy between the simulations and the measurements.

Conclusion
In this paper, we have presented a layer-averaged version of the 3D incompressible, hydrostatic Euler and Navier-Stokes systems with free surface. Compared to previous works of some of the authors, the numerical scheme is improved (implicit treatment of the vertical exchanges, description of the topography source term, . . . ) and hence, based on a kinetic interpretation of the system for the Euler part, we have derived a stable, robust and efficient numerical scheme in a finite volume/finite element framework on fixed unstructured meshes. The numerical scheme is endowed with strong stability properties (domain invariant, well-balancing, wet/dry interfaces treatment, . . . ).
The numerical scheme is successfully validated with analytical solutions and is shown to be applicable to simulate complex test cases, like a tsunami propagation over a real bathymetry, proving its accuracy and efficiency. and using the closure relation (3.21), it comes for α = 2, . . ., N − 1 R x,α u α = ∂ (u α h α Σ xx,α ) ∂x + ∂ (u α h α Σ xy,α ) ∂y An analogous relation can be easily obtained for R y,α v α and computing the sum over the layers of the obtained quantity, we get proving the result.