Alexandre Chiapolino
email: alexandre.chiapolino@rs2n.eu

Richard Saurel
email: richard.saurel@univ-amu.fr

Models and methods for two-layer shallow water flows

Keywords: multiphase flows, shallow water, hyperbolic systems, Riemann problem, non-conservative, relaxation

Models and methods for two-layer shallow water flows

Introduction

Two-layer (and multi-layer) shallow water models are particularly useful in some limit cases of multi-fluid and variable density flows separated by nearly horizontal interfaces. These models govern the dynamics of incompressible fluids spreading under gravity effects. It can be for example:

-Flows of the same liquid but at different temperatures, resulting in density differences, such situation being typical of oceanic flows; -Flows of two liquids of different densities; -Flows of two gases evolving at low Mach number.

The two-layer approach is particularly interesting compared to multidimensional approaches, that consider vertical motion, as it enables much faster computations. It is also helpful when the height of one of the phases is arbitrarily small, as there is no need to spatially resolve it.

Thereby, no numerical diffusion of the nearly horizontal interface is present and no interface tracking is needed. However, there are obviously some limitations with this approach:

-The vertical velocity component is neglected; -The velocity is assumed uniform in cross sections of each layer. Such type of modeling also involves serious difficulties. Indeed, most models are not hyperbolic, this issue having serious consequences both for propagation phenomenon, which becomes ill-posed, and for the design of numerical methods. A second serious difficulty appears as nonconservative terms are present in the momentum equations. The present paper addresses these two difficulties and provides solutions.

In the frame of averaged (or homogenized) equations in fluid mechanics, the issue related to the lack of hyperbolicity appears in different type of models, such as those of non-equilibrium two-phase flows. Only a few models seem well-posed with this respect, [START_REF] Marble | Dynamics of a containing small solid particles[END_REF] [START_REF] Marble | Dynamics of a containing small solid particles[END_REF], [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF], [START_REF] Saurel | Modelling compressible dense and dilute two-phase flows[END_REF] [START_REF] Saurel | Modelling compressible dense and dilute two-phase flows[END_REF]). There are mainly two types of remedy to cure this issue: 2 -Consider compressibility of the phases and deal with pressure relaxation [START_REF] Lallemand | Pressure relaxation procedures for multiphase compressible flows[END_REF]. This approach involves sound propagation in the phases and is particularly efficient in many situations.

It has been adopted in the last two above-mentioned references.

-Consider turbulent effects in the phases, as they result in the appearance of a "turbulent sound speed" [START_REF] Forestier | Solveur de type Godunov pour simuler les écoulements turbulents compressibles[END_REF] [START_REF] Forestier | Solveur de type Godunov pour simuler les écoulements turbulents compressibles[END_REF], [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shock-bubble interaction[END_REF] [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shock-bubble interaction[END_REF], Lhuillier et al. (2013) [START_REF] Lhuillier | On the quest for a hyperbolic effective-field model of disperse flows[END_REF]). In the frame of shallow water flows, these effects have been studied in [START_REF] Richard | A new model of roll waves: comparison with Brock's experiments[END_REF] [START_REF] Richard | A new model of roll waves: comparison with Brock's experiments[END_REF] and [START_REF] Gavrilyuk | Spilling breakers in shallow water: Applications to Favre waves and to the shoaling and breaking of solitary waves[END_REF] [START_REF] Gavrilyuk | Spilling breakers in shallow water: Applications to Favre waves and to the shoaling and breaking of solitary waves[END_REF].

In the present work, the first method is adopted and the fluids are considered weakly compressible. The resulting model is strictly hyperbolic and in the limit of stiff pressure relaxation, the conventional (non-hyperbolic) two-layer model is recovered. This approach is reminiscent of the model of [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF], except that extra pressure terms are present in the momentum equations of the new formulation. It also gives another interpretation of the relaxation approach, now based on compressibility and pressure effects.

The second issue is addressed as well and is related to the presence of non-conservative terms in the momentum equations. By examining the Riemann problem structure, it appears that local constants are present, at locations where the derivative of the Heaviside function emerges.

Consequently the non-conservative products become well-defined. Also, local conservation laws are obtained and used in the frame of HLL-type Riemann solver.

The accuracy of the new solver is checked against results of [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] as well as results obtained with a flow solver based on the VFRoe method of [START_REF] Gallouet | Un schéma Godunov approché[END_REF] [START_REF] Gallouet | Un schéma Godunov approché[END_REF] as it is able to deal, to some extent, with both conservative and non-conservative systems. The new method, based on HLL-type solver shows results of high accuracy and is oscillation free.

The paper is organized as follows. The two-layer hyperbolic model is presented in Section 2 and its stiff mechanical relaxation limit is examined. Both approximate VFRoe solver and non-conservative HLL solver are considered in Section 3. A Godunov-type scheme is derived in the same section. Results and validations are addressed in Sections 4 and 5. Conclusions are given in Section 6.

Hyperbolic two-layer shallow water model

The conventional two-layer shallow water model [START_REF] Ovsyannikov | Two-layer shallow water model[END_REF] [START_REF] Ovsyannikov | Two-layer shallow water model[END_REF]) reads,

                       ∂(h 1 ρ 1 ) ∂t + ∂ (h 1 ρ 1 u 1 ) ∂x = 0, ∂ (h 1 ρ 1 u 1 ) ∂t + ∂ h 1 ρ 1 u 2 1 + 1 2 ρ 1 gh 2 1 + gρ 2 h 1 h 2 ∂x = ρ 2 gh 2 ∂h 1 ∂x , ∂(h 2 ρ 2 ) ∂t + ∂ (h 2 ρ 2 u 2 ) ∂x = 0, ∂ (h 2 ρ 2 u 2 ) ∂t + ∂ h 2 ρ 2 u 2 2 + 1 2 ρ 2 gh 2 2 ∂x = -ρ 2 gh 2 ∂h 1 ∂x .
(2.1) h 1 and h 2 denote the heights of the two layers, ρ 1 and ρ 2 represent the densities of the fluids, considered constant at this level, u 1 and u 2 denote the fluid velocities, averaged in each layer and g represents the gravity constant. Topography effects have been omitted for the sake of simplicity as well as friction with the bottom and between layers.

System (2.1) has been examined in [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF], [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF] [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF] and [START_REF] Monjarret | Local well-posedness of the two-layer shallow water model with free surface[END_REF] [START_REF] Monjarret | Local well-posedness of the two-layer shallow water model with free surface[END_REF] and appeared hyperbolic for small velocity drift only,

(u 1 -u 2 ) 2 < (h 1 + h 2 )g 1 - ρ 2 ρ 1 . (2.2)
Moreover the wave speeds can hardly be computed, rendering the system intricate to solve numerically. A method is given in [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF] [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF] to overcome this difficulty.

In the present approach, pressure non-equilibrium effects result in an unconditionally hyperbolic formulation:

                                       ∂h 1 ∂t + u 1 ∂h 1 ∂x = µ(p 1 -p 0 ) ρ 1 c 2 1 , ∂(h 1 ρ 1 ) ∂t + ∂ (h 1 ρ 1 u 1 ) ∂x = 0, ∂ (h 1 ρ 1 u 1 ) ∂t + ∂ h 1 ρ 1 u 2 1 + h 1 p 1 (ρ 1 ) + 1 2 ρ 1 gh 2 1 + gρ 2 h 1 h 2 ∂x = ρ 2 gh 2 ∂h 1 ∂x + p 0 ∂h 1 ∂x , ∂h 2 ∂t + u 2 ∂h 2 ∂x = µ(p 2 -p 0 ) ρ 2 c 2 2 , ∂(h 2 ρ 2 ) ∂t + ∂ (h 2 ρ 2 u 2 ) ∂x = 0, ∂ (h 2 ρ 2 u 2 ) ∂t + ∂ h 2 ρ 2 u 2 2 + h 2 p 2 (ρ 2 ) + 1 2 ρ 2 gh 2 2 ∂x = -ρ 2 gh 2 ∂h 1 ∂x + p 0 ∂h 2 ∂x .
(

p k (ρ k ) = p (0) k + c 2 k ρ k -ρ (0) k , (2.4)
with k = 1, 2. Other options, such as Tait EOS for instance are possible. We will see that the choice of the EOS is not important, only the related sound speed c k has importance. p 0 denotes the (constant) atmospheric pressure and p

(0) k = p 0 .
The assumption of constant atmospheric pressure is quite realistic when dealing with flows of gases having different densities. When dealing with liquids and large hydrostatic effects, the model can be reformulated with extra pressure terms as a function of heights of the fluid layers, and does not pose extra fundamental issues than those already addressed in the present contribution. Such extension is given in Appendix, but for the sake of simplicity of the presentation, this extended model is not considered in the rest of the paper.

The pressure relaxation parameter µ is related to the fluid sound speeds and heights of layers.

It controls the rate at which pressure equilibrium is reached. Following [START_REF] Saurel | Modelling compressible dense and dilute two-phase flows[END_REF] [START_REF] Saurel | Modelling compressible dense and dilute two-phase flows[END_REF], the first equation of System (2.3) can be written as,

d 1 h 1 dt = h 1 τ p 1 -p 0 ρ 1 c 2 1 ,
where

d 1 dt = ∂ ∂t + u 1 ∂
∂x and τ is the pressure relaxation time,

τ = h 1 c 1 , (2.5)
corresponding to the following pressure relaxation parameter estimate:

µ ≃ h 1 τ ≃ c 1 . (2.6)
In most situations, this relaxation time is of the order of 1 100 second, meaning that the relaxation parameter µ is large:

µ ≃ Max τ -1 1 , τ -1 2
or alternatively µ ≃ Min (c 1 , c 2 ). In practical computations, the relaxation time τ will be assumed of the same order as the computational time step and stiff pressure relaxation will be done at the end of each time step. Therefore, there is no need of precise knowledge of the pressure relaxation parameter µ.

This system is reminiscent of Baer and Nunziato's (1986) [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] model widely used in two-phase flow modeling. It is also reminiscent of Abgrall and Karni's (2009) [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] relaxation model, except that pressure terms have been added to the momentum equations (h 1 p 1 (ρ 1 ) and h 2 p 2 (ρ 2 )). To maintain mechanical equilibrium, extra non-conservative terms have been added in the righthand side (p 0 ∂h k ∂x ). These terms are not in contradiction with the total momentum conservation that reads:

∂ (h 1 ρ 1 u 1 + h 2 ρ 2 u 2 ) ∂t + ∂ h 1 ρ 1 u 2 1 + h 1 p 1 (ρ 1 ) + 1 2 ρ 1 gh 2 1 + gρ 2 h 1 h 2 + h 2 ρ 2 u 2 2 + h 2 p 2 (ρ 2 ) + 1 2 ρ 2 gh 2 2 -p 0 (h 1 + h 2 ) ∂x = 0
Let us now examine some relevant properties to check validity of this formulation.

Hyperbolicity

System (2.3) is expressed in primitive-variable formulation (in the absence of source terms) as,

∂W ∂t + A(W ) ∂W ∂x = 0, (2.7)      u 1 0 0 0 0 0 0 u 1 ρ 1 0 0 0 p 1 -p 0 h 1 ρ 1 + g c 2 1 + 1 2 gh 1 ρ 1 u 1 g ρ 2 ρ 1 gh 2 ρ 1 0 0 0 0 u 2 0 0 0 0 0 0 u 2 ρ 2 g 0 0 p 2 -p 0 h 2 ρ 2 + g c 2 2 + 1 2 gh 2 ρ 2 u 2               . (2.8)
The wave speeds are solutions of det A(W ) -λI = 0 resulting in,

(u 1 -λ) (u 2 -λ) (u 2 -λ) 2 -c 2 2 + 1 2 gh 2 (u 1 -λ) 2 -c 2 1 + 1 2 gh 1 = 0. (2.9)
Six real and distinct eigenvalues appear as:

       λ 1 = u 1 , λ 2 = u 1 + c 2 1 + 1 2 gh 1 , λ 3 = u 1 -c 2 1 + 1 2 gh 1 , λ 4 = u 2 , λ 5 = u 2 + c 2 2 + 1 2 gh 2 , λ 6 = u 2 -c 2 2 + 1 2 gh 2 .
(2.10)

Those eigenvalues correspond to the wave speeds emerging at a given initial discontinuity, as schematized in Fig. 1.

System (2.
3) is consequently strictly hyperbolic. This model is however relevant with respect to the physics expressed in (2.1) if it tends to the same equations when pressure relaxation is stiff. This limit is examined hereafter.

Stiff pressure relaxation limit

The mass and height equations of a given phase are analyzed:

       ∂h k ∂t + u k ∂h k ∂x = µ(p k -p 0 ) ρ k c 2 k , ∂(h k ρ k ) ∂t + ∂ (h k ρ k u k ) ∂x = 0, with k = 1, 2.
Their combination results in,

d k ρ k dt + ρ k ∂u k ∂x = - ρ k h k µ (p k -p 0 ) ρ k c 2 k .
Inserting the phase k equation of state p k (ρ k ), the following pressure evolution equation is obtained,

d k p k dt + ρ k c 2 k ∂u k ∂x = - µ (p k -p 0 ) h k .
As the atmospheric pressure p 0 is constant, the last equation can be expressed as,

d k (p k -p 0 ) dt + ρ k c 2 k ∂u k ∂x = - µ (p k -p 0 ) h k .
The phase pressure is expressed around the equilibrium state with the following expansion,

p k = p (0) k + ǫp (1) k + . . . ,
where ǫ is of the order of the inverse of pressure relaxation parameter (ǫ ≃ µ -1 ≃ τ ), tending to zero in most situations (ǫ → 0 + ) as discussed earlier (see also [START_REF] Kapila | Two-phase modeling of deflagrationto-detonation transition in granular materials: Reduced equations[END_REF] [START_REF] Kapila | Two-phase modeling of deflagrationto-detonation transition in granular materials: Reduced equations[END_REF] for estimates in the context of granular flows), p

k and p

k represent respectively the leading and first-order pressure terms of the Taylor expansion.

Inserting these definitions in the pressure evolution equation,

d k p (0) k + ǫp (1) k + . . . -p 0 dt + ρ k c 2 k ∂u k ∂x = - p (0) k + ǫp (1) k + . . . -p 0 h k ǫ ,
the following results are obtained:

-At leading order (ǫ -1 ):

p (0) k = p 0 ; -At first order p (1)
k = -ρ k c 2 k h k ∂u k ∂x .
Inserting this last result in the height equations,

∂h k ∂t + u k ∂h k ∂x = µ (p k -p 0 ) ρ k c 2 k ≃ p (1) k ρ k c 2 k ≃ -h k ∂u k ∂x ,
they become,

∂h k ∂t + ∂ (h k u k ) ∂x ≃ 0.
The mass equations are unchanged while modifications in the momentum equations appear as a consequence of the equilibrium condition (p k = p 0 ). They finally result at leading order in,

       ∂ (h 1 ρ 1 u 1 ) ∂t + ∂ h 1 ρ 1 u 2 1 + 1 2 ρ 1 gh 2 1 + gρ 2 h 1 h 2 ∂x = ρ 2 gh 2 ∂h 1 ∂x , ∂ (h 2 ρ 2 u 2 ) ∂t + ∂ h 2 ρ 2 u 2 2 + 1 2 ρ 2 gh 2 2 ∂x = -ρ 2 gh 2 ∂h 1 ∂x .
System (2.1) is recovered, complemented by two conservation equations for the heights, that are in agreement with the two mass equations as soon as the densities are constants.

It thus appears that System (2.3) tends to System (2.1) when pressure relaxation is stiff.

As System (2.3) is hyperbolic, it is a good candidate to approximate (2.1) numerically with a two-step procedure:

-Solve the hyperbolic system (2.3) without source terms;

-Relax the pressures onto the atmospheric one and reset the heights. This is similar to the method of Saurel and Abgrall (1999) [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF] to compute flows with interfaces separating fluids. Before entering in the details of the hyperbolic solver, let us present the pressure relaxation process, that is particularly simple in the present context.

Stiff pressure relaxation solver

Let us consider for example EOS (2.4). Consequently, the densities as functions of pressures are given by:

ρ k = ρ (0) k + p k -p (0) k c 2 k . (2.11)
As the pressures relax to the atmospheric one (p k = p (0) k = p 0 ), the densities at relaxed pressure are just,

ρ * k = ρ (0) k , (2.12)
where the superscript * denotes the relaxed pressure state. As the masses of each layer are computed by associated mass balance equations and are constant during the relaxation process,

m k = h k ρ k = h * k ρ * k , (2.13)
the heights at relaxed states are reset as,

h * k = h k ρ k ρ (0) k . (2.14)
The stiff pressure relaxation solver just consists in the reset of the heights of the fluids h k → h * k and is independent of the equations of state. At this level, the relaxation method of Abgrall and Karni (2009) [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] is recovered. We now address the design of hyperbolic solvers.

Approximate Riemann solvers

Two different approximate methods are considered to solve the Riemann problem of System (2.3), the VFRoe method [START_REF] Gallouet | Un schéma Godunov approché[END_REF] [START_REF] Gallouet | Un schéma Godunov approché[END_REF]) and a new HLL-type Riemann solver [START_REF] Harten | On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws[END_REF] [START_REF] Harten | On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws[END_REF]). As System (2.3) contains non-conservative terms, an approach dealing with both conservative and non-conservative systems is examined first.

VFRoe solver

The VFRoe method considers the equations in non-conservative formulation,

∂W ∂t + A(W ) ∂W ∂x = 0, (3.1)
with

W = (h 1 , ρ 1 , u 1 , h 2 , ρ 2 , u 2 ) T and W = W L + W R 2 ,
where W L and W R are respectively the left and right-state vectors at a given cell boundary.

The VFRoe method considers the exact Riemann problem solution of (3.1). Note that (3.1) is a local linearization of the non-linear flow model (2.3) around state W . The VFRoe solution is thus the exact solution of an approximate problem.

The exact solution of (3.1) may be found in many textbooks related to hyperbolic systems (LeVeque (2002) [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF], Toro (2013) [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: A practical introduction[END_REF]) and can be summarized as follows,

W * = W L + λ i <0 a i R i = W R - λ i >0 a i R i , (3.2)
where the wave strengths a i are the coefficients resulting from the decomposition of the eigenvectors.

W R -W L = λ i a i R i . (3.3)
For the sake of space restriction, the right eigenvectors R i and the wave strengths a i are not detailed, associated formulas being considerably large. The main weakness of this method is related to the average W which can be far from the solution of the non-linear problem, resulting in positivity issues especially when large amplitude waves are present.

With the help of the Riemann problem solution (3.2), the various equations of System (2.3) are updated with a Godunov-type method (stable under the conventional CF L condition) as,

h n+1 k = h n k - ∆t ∆x (h k u k ) * i+ 1 2 -(h k u k ) * i-1 2 + ∆t ∆x h n k,i u * k,i+ 1 2 -u * k,i-1 2 , (3.4)
(h k ρ k ) n+1 i = (h k ρ k ) n i - ∆t ∆x (h k ρ k u k ) * i+ 1 2 -(h k ρ k u k ) * i-1 2 , (3.5)
(h 1 ρ 1 u 1 ) n+1 i = (h 1 ρ 1 u 1 ) n i - ∆t ∆x h 1 ρ 1 u 2 1 + h 1 (p 1 -p 0 ) + 1 2 ρ 1 h 2 1 * i+ 1 2 -h 1 ρ 1 u 2 1 + h 1 (p 1 -p 0 ) + 1 2 ρ 1 h 2 1 * i-1 2 + ∆t ∆x h n 1,i (-ρ 2 gh 2 ) * i+ 1 2 -(-ρ 2 gh 2 ) * i-1 2 , (3.6)
(h 2 ρ 2 u 2 ) n+1 i = (h 2 ρ 2 u 2 ) n i - ∆t ∆x h 2 ρ 2 u 2 2 + h 2 (p 2 -p 0 ) + 1 2 ρ 2 h 2 2 + ρ 2 gh 1 h 2 * i+ 1 2 -h 2 ρ 2 u 2 2 + h 2 (p 2 -p 0 ) + 1 2 ρ 2 h 2 2 + ρ 2 gh 1 h 2 * i-1 2 + ∆t ∆x h n 1,i (ρ 2 gh 2 ) * i+ 1 2 -(ρ 2 gh 2 ) * i-1 2 , (3.7)
where n + 1 and n denote two consecutive time steps and superscript * denotes the VFRoe Riemann problem solution given by Eq. (3.2). Indexes i and i ± 1 2 denote respectively the center of the current numerical cell and its corresponding boundaries.

HLL-type Riemann solver

Let us consider a simplified solver, based on Rankine-Hugoniot conditions, such as the HLL solver. In this frame, the two extreme waves S L and S R are approximated following [START_REF] Davis | Simplified second-order Godunov-type methods[END_REF] [20] as,

           S L,k = min u L,k -c 2 L,k + 1 2 gh L,k , u R,k -c 2 R,k + 1 2 gh R,k , S R,k = max u L,k + c 2 L,k + 1 2 gh L,k , u R,k + c 2 R,k + 1 2 gh R,k , (3.8)
with k = 1, 2. The indexes L and R denote respectively the left and right states at a given cell boundary. The two extreme waves are considered as,

S L = min (S L,1 , S L,2 ) , S R = max (S R,1 , S R,2 ) . (3.9)
The two contact waves u 1 and u 2 are considered as well for the transport of the heights h 1 and h 2 , as depicted in Fig. 2. Regarding the transport equations, the exact Riemann problem solution is straightforward:

x t S R S L u 1 u 2 W R W L W * W * L W * R
     h * 1 x t < u * 1 = h 1,L , h * 1 x t > u * 1 = h 1,R , h * 2 x t < u * 2 = h 2,L , h * 2 x t > u * 2 = h 2,R . (3.10)
These solutions indicate that the non-conservative terms have contributions between the two extreme waves S R and S L , at points where h 1 and h 2 are discontinuous. More precisely, only the discontinuity in h 1 needs attention, as the non-conservative terms involving the atmospheric pressure (considered constant) transform to fluxes,

p 0 ∂h k ∂x = ∂(p 0 h k ) ∂x .
It thus remains to analyze only the non-conservative term,

ρ 2 gh 2 ∂h 1 ∂x .
The solution states for (ρ 2 h 2 ) are given by,

(ρ 2 h 2 ) * L = (ρ 2 h 2 ) L u 2,L -S L u * 2 -S L and (ρ 2 h 2 ) * R = (ρ 2 h 2 ) R u 2,R -S R u * 2 -S R .
A schematic representation is given in Fig. 3. These formulas need u * 2 for their practical use that is unknown at this level. However, according to the sign of the velocity difference u * 1u * 2 , only two instances may occur:

x t S R S L u 1 u 2 W R W L (h 2 ρ 2 ) * L (h 2 ρ 2 ) * R
-If u * 1 > u * 2 , the ρ 2 gh 2 ∂h 1 ∂x term becomes locally g (ρ 2 h 2 ) * R ∂h 1 ∂x . As (ρ 2 h 2 ) * R is constant at the point where ∂h 1
∂x is discontinuous, the non-conservative term becomes locally;

ρ 2 gh 2 ∂h 1 ∂x = ∂ g (ρ 2 h 2 ) * R h 1 ∂x .
-If u * 1 < u * 2 , the same reasoning yields,

ρ 2 gh 2 ∂h 1 ∂x = ∂ g (ρ 2 h 2 ) * L h 1 ∂x .
It thus appears that the momentum equations are locally conservative. However, their explicit determination and use require knowledge of both u * 1 and u * 2 that are themselves solutions of integration of the momentum equations.

To simplify the algorithm, a single solution state is considered for the apparent densities

(ρ k h k ) * instead of the two (ρ k h k ) *
L and (ρ k h k ) * R in the same mind as in the HLL solver for the Euler equations:

(h k ρ k ) * = (h k ρ k ) R (u k,R -S R ) -(h k ρ k ) L (u k,L -S L ) S L -S R . (3.11)
Thanks to this approximation, the momentum equations become locally,

       ∂ (h 1 ρ 1 u 1 ) ∂t + ∂ h 1 ρ 1 u 2 1 + h 1 p 1 (ρ 1 ) -p 0 + 1 2 ρ 1 gh 2 1 + gρ 2 h 1 h 2 -g(h 2 ρ 2 ) * h 1 ∂x = 0, ∂(h 2 ρ 2 u 2 ) ∂t + ∂ h 2 ρ 2 u 2 2 + h 2 p 2 (ρ 2 ) -p 0 + 1 2 ρ 2 gh 2 2 + g(h 2 ρ 2 ) * h 1 ∂x = 0. (3.12)
Denoting the momentum fluxes by,

     F 1,mom =h 1 ρ 1 u 2 1 + h 1 p 1 (ρ 1 ) -p 0 + 1 2 ρ 1 gh 2 1 + gρ 2 h 1 h 2 -g(h 2 ρ 2 ) * h 1 , F 2,mom =h 2 ρ 2 u 2 2 + h 2 p 2 (ρ 2 ) -p 0 + 1 2 ρ 2 gh 2 2 + g(h 2 ρ 2 ) * h 1 , (3.13)
the momentum numerical fluxes are then given by,

F * k,mom = F k,mom,R S L -F k,mom,L S R + S L S R (U k,mom,L -U k,mom,R ) S L -S R , (3.14)
with U k,mom = h k ρ k u k . The mass numerical fluxes are computed by the HLL approximation as well and read,

F * k,mass = (h k ρ k ) R S L (u k,R -S R ) -(h k ρ k ) L S R (u k,L -S L ) S L -S R . (3.15)
System (2.1) being non-conservative, the conservative variable-state vector and in particular the fluid velocities are also needed for the computations. With the help of Eqs. (3.12) and (3.13), the momentum variables are computed with the HLL approximation as,

(h k ρ k u k ) * = F k,mom,R -F k,mom,L -S R U k,mom,R + S L U k,mom,L S L -S R . (3.16)
Using Eqs. (3.11) and (3.16), the speeds of the fluids are given by,

u * k = (h k ρ k u k ) * (h k ρ k ) * . (3.17)
The associated Godunov-type method now reads,

                     h n+1 k = h n k - ∆t ∆x (hu) * k,i+ 1 2 -(hu) * k,i-1 2 + ∆t ∆x h n k,i u * k,i+ 1 2 -u * k,i-1 2 , (h k ρ k ) n+1 i = (h k ρ k ) n i - ∆t ∆x F * k,mass,i+ 1 2 -F * k,mass,i-1 2 , (h 1 ρ 1 u 1 ) n+1 i = (h 1 ρ 1 u 1 ) n i - ∆t ∆x F * 1,mom,i+ 1 2 -F * 1,mom,i-1 2 + ∆t ∆x h n 1,i -g (h 2 ρ 2 ) * i+ 1 2 -(h 2 ρ 2 ) * i-1 2 , (h 2 ρ 2 u 2 ) n+1 i = (h 2 ρ 2 u 2 ) n i - ∆t ∆x F * 2,mom,i+ 1 2 -F * 2,mom,i-1 2 + ∆t ∆x h n 1,i g (h 2 ρ 2 ) * i+ 1 2 -(h 2 ρ 2 ) * i-1 2 . (3.18)
Efficiency of both VFRoe and HLL solvers are now investigated on various test problems of the literature.

Results and validations

It is important to address in priority the effects of the fluid EOS with the present relaxation approach. Indeed the model is hyperbolic as a consequence of compressibility terms in the momentum equations. Extra tests, where VFRoe and HLL solvers are compared, are addressed subsequently.

Effects of the artificial sound speed

The EOS (2.4) involves sound speed c k that has influence on computed results, as shown hereafter in Figs. [START_REF] Lallemand | Pressure relaxation procedures for multiphase compressible flows[END_REF], 5 and 6. The examined configuration consists in a limit case where the initial height of the first fluid (lower layer) is as low as numerically admissible, so that only the upper layer (second fluid) evolves significantly. With this specific configuration, the solution of the two-layer shallow water system (2.3) is meant to be compared to the exact solution of the one-layer Saint-Venant equations.

A dam-break problem is used to this end. The following test is proposed in LeVeque's textbook (2002) [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] to illustrate behavior of the solution of the Saint-Venant equations. It consists in a dam, separating two levels of fluids, that bursts at time t = 0. All variables of the current test problem are in dimensionless units as done in [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. This test is the shallow water equivalent of the shock-tube problem of gas dynamics and appears to be an excellent benchmark as the flow deals with shock and expansion waves that create arduous conditions. The constant gravity is normalized and reads g = 1. The numerical domain has a length set to 10 with a height discontinuity initially located in the middle. On the left of this discontinuity, the fluid is initially at h = 3 and h = 1 on the right. The fluid is initially at rest on either side of the discontinuity.

To mimic the Saint-Venant system with the two-layer model (2.3), the height of the first fluid is initially set to h 1 = ǫ = 10 -6 throughout the whole numerical domain. Its density is set to ρ 1 = 1.2 and its velocity is set to u 1 = 0. The second fluid, placed above the first one, has initial heights h 2 = 3 at left and h 2 = 1 at right. Its density is set to ρ 2 = 1 and its velocity is set to u 2 = 0.

Figure 4 shows the results with constant sound speed set to c k = 100. The same sound speed has been taken for both fluids for the sake of simplicity in this illustration. Also, two mesh resolutions are used in Fig. 4, a coarse one made of 100 cells (A) and a fine one made of 10, 000 cells (B).

It appears that the speed of sound influences computed results. The wave speeds of the two-layer model (2.3) involve the effective sound speeds given by c 2 k + 1 2 gh k while the singlelayer wave speed is √ gh. When c k is significantly greater than 1 2 gh k (of the order of unity in the present example), excessive numerical diffusion is present, as shown in Plot A of Fig. 4 where c k = 100 for both fluids. Indeed, at the current time, the left and right-facing waves are considerably dissipated and even exit the domain. However, this feature is purely numerical and the system does converge to the exact solution, as shown in Plot B of Fig. 4 that uses a fine grid made of 10, 000 cells.

c k > 1 2 gh k . (4.1)
In the upcoming computations, the following sound speed is used:

c k = θ k 1 2 gh k , with θ k > 1, k = 1, 2. (4.2)
θ k is a numerical parameter that controls the numerical diffusion as illustrated in Fig. 5.

h lef t 1 = h right 1 = 10 -6 , u lef t 1 = u right 1 = 0, ρ 1 = 1.2, h lef t 2 = 3, h right 2 = 1, u lef t 2 = u right 2 = 0, ρ 2 = 1. Final time: t ≈ 2.
All results use a 100-cell mesh. First-order Godunov-type numerical scheme is used with CF L = 0.9. For the sake of clarity, 50 symbols are plotted for the exact solution. An optimum appears for θ k = 2. All variables are in dimensionless units.

In order to unambiguously fulfill the above-mentioned subcharacteristic condition (4.1), θ k must be greater than unit. θ k ∈ [2 , 5] seems to be a fair choice as it is low enough to control numerical diffusion and high enough to ensure stability. Indeed, as seen in Plot D of Fig. 5, θ k = 2 provides accurate results whereas θ k = 10 (Plot E) and θ k = 50 (Plot F) show excessive numerical dissipation.

Figure 6 repeats the same test with θ k = 2 and a 1000-cell grid. The Godunov method (3.18) including non-conservation terms is extended to second order with the MUSCL-type method detailed for example in Toro (2013) [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: A practical introduction[END_REF] (see also [START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF] [START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF] when non-conservative terms are present).

The results show excellent agreement with the exact solution. Besides, they also reveal that:

-Incompressible behavior is recovered as the densities ρ k are constant; -The two-layer shallow water model (2.3) tends to the single-layer Saint-Venant equations in the limit h 1 → ǫ; -Second-order extension of the Godunov-type scheme (3.18) and associated non-conservative terms does not pose specific difficulties. These various computations have been done with the HLL solver while the VFRoe one failed immediately, as it was unable to preserve positivity of the height h 1 . It is also important to note that the method does not require any fluid EOS, nor relaxation parameter, as Eq. (2.14) determines efficiently the heights at relaxed states. The only "thermodynamic" information is the sound speed, and more precisely θ k in Eq. (4.2). The method is robust and accurate with 2 ≤ θ k ≤ 5.

Effects of the fluid densities

The preceding dam-break problem showed that the two-layer shallow water model (2.3) is able to recover the single-layer Saint-Venant system in the limit h 1 → 0. When the density ratio r = ρ 2 ρ 1 is small, the effects of the surrounding fluid (upper layer) are expected to become insignificant and the one-layer Saint-Venant solution is meant to be recovered as well.

h lef t 1 = h right 1 = 10 -6 , u lef t 1 = u right 1 = 0, ρ 1 = 1.2, h lef t 2 = 3, h right 2 = 1, u lef t 2 = u right 2 = 0, ρ 2 = 1. Final time: t ≈ 2.
Second-order MUSCL-type numerical scheme using van Leer's limiter (see [START_REF] Van Leer | Towards the ultimate difference scheme. II. Monotonicity and conservation combined in a second-order scheme[END_REF], [START_REF] Sweby | High resolution schemes using flux limiters for hyperbolic conservation laws[END_REF]) is considered with CF L = 0.5 and 1000 cells. For the sake of clarity, only 50 symbols out of 1000 are plotted for the HLL-type computations. All variables are in dimensionless units.

The forthcoming tests analyze the effects the fluid densities on a configuration presenting initially a Heaviside function regarding the height of the first fluid (lower layer) located in the middle of the numerical domain. The first layer is initially at height h 1 above the flat ground and the top of the plateau is located at height h ′ 1 . The second fluid surrounds the lower layer and is set initially at constant height h 2 . The initial configuration is schematically depicted in Fig. 7 with data summarized in Table 1. Test In all following tests, the gravity constant is set to g = 10 m.s -2 and the top of the Heaviside plateau is at h ′ 1 = 10 m. The results are shown at time t = 1 s and computed on a 1000-cell mesh with second-order MUSCL-type scheme and CF L = 0.5. Figure 8 Small density ratio: r = ρ 2 ρ 1 ≪ 1 When the density ratio is small, such as the situation of Fig. 8 G and H considering a waterair-type configuration, the solutions of the two-layer system (2.3) are in excellent agreement with the exact solutions of the Saint-Venant model. For these two tests, h ′ 1 = 10 m, h 2 = 20 m and h 1 = 4 m for Plot G and h 1 = 10 -6 m for Plot H.

h 1 h ′ 1 h 2 x 1 = 20 m x 2 = 10 m x 3 = 20 m ρ 1 ρ 2 g x y h 2
h 1 (m) h ′ 1 (m) h 2 (m) ρ 1 (kg.m -3 ) ρ 2 (kg.m -3 ) G 4 10
In Plot G of Fig. 8, the Heaviside profile of the lower layer gives rise to four waves moving in each direction. The extreme waves steepen through compression waves into two shocks, while the back waves spread out as rarefaction waves.

However, the solution in Plot H of Fig. 8 is significantly different as the initial height of the lower layer is as low as numerically acceptable (outside the Heaviside profile). Thereby the solution evolves continuously throughout the entire numerical domain as there is not enough fluid, regarding the lower layer (h 1 ), to observe a compression process: expansion waves only are present.

Comparable densities: r = ρ 2 ρ 1 → 1 Plots I and J of Fig. 8 x (m) x (m) gh k with θ k = 2. The dashed lines represent the initial conditions. Final time: t ≈ 1 s. All results use a 1000-cell mesh. Second-order MUSCL-type numerical scheme is used with Sweby's limiter (φ = 1.1, see [START_REF] Sweby | High resolution schemes using flux limiters for hyperbolic conservation laws[END_REF]) and CF L = 0.5. For the sake of clarity, 50 symbols are plotted for the exact solution of the single-layer model. On the column at left, the initial height is h 1 = 4 m and on the column at right h 1 = 10 -6 m. When the density ratio is small: r = ρ2 ρ1 ≪ 1, as in configurations G and H, the two-layer and single-layer models are in perfect agreement. Large differences appear when the density ratio increases as shown in configurations I, J, K and L.

h 1 (m) h 1 (m) h 1 (m) h 1 (m) h 1 (m) h 1 (m) G H I J K L
interaction of the two fluids influences significantly the computed results as an interesting wave structure appears in Plot I of Fig. 8. Right and left-facing shock waves propagate faster than those of the single-layer system. These shocks induce height increase of the fluid layer.

They are followed by expansion waves that decrease these heights. Contact waves follow these expansion waves, followed by extra expansion waves that decrease the initial height h ′ 1 . The solution is quite different when h 1 = 10 -6 m initially (Plot J) where only two expansion waves are observable. The interaction of the two fluids influences the flow, as the first fluid moves more difficultly into the second one as a result of comparable densities.

Large upper layer

Plots K and L of Fig. 8 keep on analyzing the present density ratio (water-oil) but with different height for the second fluid that is now set to h 2 = 100 m. Doing so, the domain is mainly filled with the upper layer and Plots K and L show that the first fluid moves into the second one difficultly. The difference of initial height for the first fluid h 1 = 4 m (Plot K) and h 1 = 10 -6 m (Plot L) is minor compared to the effect of the large layer of second fluid (h 2 = 100 m). As the entire domain is mainly filled with heavy fluids, the flow is slowly set to motion under gravity effects.

Concluding remarks

Those last results reveal that the two-fluid model (2.3) is able to recover the one-layer Saint-Venant system when the effects of the surrounding fluid are negligible, as expected. This behavior appears when the density ratio between the lightest fluid and the heaviest one is small: r =

ρ 2 ρ 1 ≪ 1.
They also reveal the importance of the two-layer model when the density ratio is arbitrary.

Indeed, the two-layer system is able to deal with interactions between fluids unlike the conventional one-fluid Saint-Venant model. Note that the previous tests have been computed with the HLL-type Riemann solver presented in Section 3.2, the VFRoe method being unable to keep positivity of the heights. Note also that drag effects are absent in these computations.

Comparison of VFRoe and HLL

two solvers considered in the present paper are tested on a flow configuration examined in Abgrall and Karni (2009) [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] and [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF] [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF]. The following tests set gravity constant to g = 10 m.s -2 and density ratio to r = ρ 2 ρ 1 = 0.98. Hereby, ρ 1 = 1200 kg.m -3 and ρ 2 = 1176 kg.m -3 are used. Initially, different heights are present from either side of the initial discontinuity and results in the creation of a flow under gravity effects. The numerical domain is 1 m long and the initial discontinuity is located at x = 0.5 m. On the left of this discontinuity h 1 = 0.5 m and h 2 = 0.5 m. On the right, h 1 = 0.45 m and h 2 = 0.55 m. The initial conditions are schematically depicted in Fig. 9. The first-order Godunov-type scheme is used with CF L = 0.7 in the following tests. Doing so, the comparison between solvers is free of extra ingredients, such as gradient limiters. However, the pressure relaxation step restores incompressibility as shown in Fig. 10 where the densities ρ 1 and ρ 2 remain constant. The numerical solutions consist in 4 waves emerging from the initial discontinuity. Both methods predict the same solution qualitatively, but the VFRoe solver produces spurious oscillations while the HLL one is monotonic.

h lef t 1 h lef t
x (m) x (m) h 2 (m) u 2 (m.s -1 ) h 1 (m) u 1 (m.s -1 ) ρ 1 (kg.m -3 ) ρ 2 (kg.m -3 )
Figure 10: Comparison of the computed solutions with the present HLL-type Riemann solver (thick lines) versus the computed solutions with the VFRoe method (thin lines and symbols). The dashed lines represent the initial conditions:

h lef t 1 = 0.5 m, h right 1 = 0.45 m, h lef t 2 = 0.5 m, h right 2 = 0.55 m, u lef t 1 = u right 1 = u lef t 2 = u right 2 = 0 m.s -1 , ρ 1 = 1200 kg.m -3 , ρ 2 = 1176 kg.m -3 .
For the first fluid, the numerical parameter reads θ 1 = 3.5, for the second fluid, θ 2 = 3. These parameters are minimum values for successful computations with VFRoe. Final time: t ≈ 0.12 s. First-order Godunov-type numerical scheme is used with CF L = 0.7 and 100 cells. For the sake of clarity, only 50 symbols out of 100 are plotted for the VFRoe method.

As initially u lef t = 0 m.s -1 , ρ 1 = 1200 kg.m -3 , ρ 2 = 1176 kg.m -3 . For the first fluid, the numerical parameter reads θ 1 = 3.5, for the second fluid, θ 2 = 3. Final time: t ≈ 0.12 s. First-order Godunov-type numerical scheme is used with CF L = 0.7 and 10, 000 cells. For the sake of clarity, only 50 symbols out of 10, 000 are plotted for the VFRoe method. Both methods converge to the same solution.

x (m) x (m) h 2 (m) u 2 (m.s -1 ) h 1 (m) u 1 (m.s -1 )
next test repeats the previous one with non-zero initial velocities. Those are set to u 1 = u 2 = 2.5 m.s -1 throughout the entire domain. This test was examined in [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] and computed with both 400 and 10, 000-cell meshes. The same grids are used here to compare the present model and HLL solver with the results given in [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF]. Figure 12 shows the results at time t ≈ 0.07 s with 400-cell mesh.

x (m) x (m) h 2 (m) u 2 (m.s -1 ) h 1 (m) u 1 (m.s -1 ) ρ 1 (kg.m -3 ) ρ 2 (kg.m -3 )
Figure 12: Comparison of the computed solutions with the present HLL-type Riemann solver (thick lines) versus the computed solutions with the VFRoe algorithm (thin lines and symbols). The dashed lines represent the initial conditions:

h lef t 1 = 0.5 m, h right 1 = 0.45 m, h lef t 2 = 0.5 m, h right 2 = 0.55 m, u lef t 1 = u right 1 = u lef t 2 = u right 2 = 2.5 m.s -1 , ρ 1 = 1200 kg.m -3 , ρ 2 = 1176 kg.m -3
. For the first fluid, the numerical parameter reads θ 1 = 3.5, for the second fluid, θ 2 = 3. Final time: t ≈ 0.07 s. First-order Godunov-type numerical scheme is used with CF L = 0.7 and 400 cells. For the sake of clarity, only 50 symbols out of 400 are plotted for the VFRoe method.

ditions are easier for the VFRoe-type that does not oscillate. The expected behavior is recovered. The heights and velocity profiles are transported to the right and the effects of gravity seen in Figs. 10 and 11 are still present. As expected the densities remain constant as a consequence of pressure relaxation. Figure 13 shows the results of the same test with 10, 000-cell mesh.

x (m) x (m) h 2 (m) u 2 (m.s -1 ) h 1 (m) u 1 (m.s -1 )
Figure 13: Comparison of the computed solutions with the present HLL-type Riemann solver (thick lines) versus the computed solutions with the VFRoe algorithm (thin lines and symbols). The dashed lines represent the initial conditions:

h lef t 1 = 0.5 m, h right 1 = 0.45 m, h lef t 2 = 0.5 m, h right 2 = 0.55 m, u lef t 1 = u right 1 = u lef t 2 = u right 2 = 2.5 m.s -1 , ρ 1 = 1200 kg.m -3 , ρ 2 = 1176 kg.m -3 .
For the first fluid, the numerical parameter reads θ 1 = 3.5, for the second fluid, θ 2 = 3. Final time: t ≈ 0.07 s. First-order Godunov-type numerical scheme is used with CF L = 0.7 and 10, 000 cells. For the sake of clarity, only 100 symbols out of 10, 000 are plotted for the VFRoe method.

The HLL-type solver and the VFRoe-type method are in excellent agreement and are in excellent agreement with the results given in [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] as well.

It is interesting to examine the rate of convergence of the present method (HLL-type solver) and give comparison to existing methods. In [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF], only first-order were carried out while in [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF] [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF], higher-order computational results are provided. Comparisons are consequently done with the results of figures 2.2, 2.3 and 2.4 of this last reference. The results are given in Fig. 14.

Figure 14 shows results comparable to those of [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF] [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF] with all mesh sizes. The second-order MUSCL-type scheme provides velocity profiles accurate enough with both 400-cell and 800-cell meshes as observed in [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF]. However, the velocity spike with the 800-cell mesh seems slightly greater than the reference solution and the close-up view on the height profile reveals that only the 800-cell mesh cancels out sufficiently numerical dissipation and allows to observe properly the constant plateau. The same conclusions are obtained in [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF] [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF]. Also, a slight oscillation is observed with the coarse mesh made of 100 cells. The same observation holds in [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF]. However, the left and right-facing shocks, seen on the velocity profile, differ from those observed in [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF] where initial data seem to be not exactly the same as those of [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF]. The present results are nonetheless in agreement with the ones of [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF] [START_REF] Abgrall | Two-layer shallow water system: A relaxation approach[END_REF]. The overall height is consequently different from the results of [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF] [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF]. Nevertheless, velocity profiles indicate that the present method seems to have a similar rate of convergence as the one of [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF].

Note that for this test problem, θ 1 = θ 2 = 2 induces spurious oscillations, similar to those observed in Fig. 5, Plot C. Consequently, θ 1 = θ 2 = 3 is used for both fluids.

Comparison of two-layer shallow water solutions versus two-dimensional two-fluid computations

The averaged (or homogenized) solution computed by the present one-dimensional two-layer shallow water system (2.3) is now compared to the solution of a multidimensional model involving material interfaces. As mentioned in the introduction, the two-layer approach is expected to provide comparable results with considerable computational savings.

In the following, solutions computed with the compressible two-phase flow model of Saurel et al. ( 2009) [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF] are used as reference solutions. This model is a pressure disequilibrium system which tends, in its asymptotic limit of stiff pressure relaxation, to the model of Kapila et al. x (m) x (m) . For the two fluids, the numerical parameters read θ 1 = θ 2 = 3. Final time: t ≈ 0.07 s. The CFL number is 0.7. A close-up view of the spanning plateau is provided at top with the h 1 profile. The overall height h 1 +h 2 +B(x) is given in the middle. B(x) represents topography of the ground and is considered constant in the present work, B(x) = 0 m. For proper comparison with the results of [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF] [START_REF] Kurganov | Central-upwind schemes for two-layer shallow water equations[END_REF] who considered B(x) = -1 m, unit is subtracted from h 1 + h 2 . Finally the velocity profile u 1 is displayed at bottom. [START_REF] Kapila | Two-phase modeling of deflagrationto-detonation transition in granular materials: Reduced equations[END_REF], able to compute fluid interfaces as diffuse numerical zones. Interface sharpening can be achieved with the method of [START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF] [START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF].

h 1 (m) h 1 (m) h 1 + h 2 -1 (m) h 1 + h 2 -1 (m) u 1 (m.s -1 ) u 1 (m.s -1 )
To compare the solutions computed by both approaches (multidimensional interface model and one-dimensional two-layer shallow water one), the test configuration schematically depicted in Fig. 15 is used with data summarized in Table 2. The accurate resolution of interfaces in computations is improved with a second-order numerical method. The MUSCL-type method with "Overbee" limiter at interfaces is used as detailed in [START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF]. The very same second-order method is used with the two-layer shallow water system, computed on a 1000-cell mesh with van Leer's limiter [START_REF] Van Leer | Towards the ultimate difference scheme. II. Monotonicity and conservation combined in a second-order scheme[END_REF]. Non-reflecting boundary conditions are used for the shallow water computations.

h 1 h ′ 1 h 2 x 1 = 40 m x 2 =
Results are shown in Fig. 16 at times t = 5 s and t = 8 s. The computed averaged height and averaged velocity of the multidimensional computation are determined by the integration on the two-dimensional numerical domain as,

h 1 = h 2 0 α 1 dy, (5.1)
and

u x = h 2 0 (α 1 ρ 1 u x ) dy h 2 0 (α 1 ρ 1 )dy , (5.2)
where α 1 denotes the volume fraction of krypton and u x denotes the averaged velocity of the two-phase mixture in the x-axis direction. The one-fluid shallow water solution is shown as well in Fig. 16 for comparison.

As expected, the two-layer model provides better results than the one-fluid model. The density ratio r = ρ 2 ρ 1 = 1.29 3.506 ≃ 0.37 being moderate, the interaction between the two fluids is meant to be significant. The results of the two-layer shallow water system present a large zone where the agreement with the two-dimensional simulation is very good. Beyond this zone, the results computed by the multidimensional interface model present oscillations. These oscillations are due to the presence of Kelvin-Helmholtz instabilities as seen for example in Fig. 17 showing the 2D results at time t = 2 s. Indeed, the multidimensional solution involves hydrodynamic instabilities that cannot be accounted for with the present two-layer shallow water model. The overall qualitative behavior of the one-dimensional approach is correct but the krypton is spread too far ahead with the two-layer model.

To improve agreement between 1D and 2D computations, drag effects are added in the two- layer formulation. Pressure (or "acoustic") drag considered only and is modeled through the velocity relaxation terms that appear in the right-hand side of the momentum equations,

       ∂ (h 1 ρ 1 u 1 ) ∂t + ∂ h 1 ρ 1 u 2 1 + h 1 p 1 (ρ 1 ) + 1 2 ρ 1 gh 2 1 + gρ 2 h 1 h 2 ∂x = ρ 2 gh 2 ∂h 1 ∂x + p 0 ∂h 1 ∂x + Z 1 Z 2 Z 1 + Z 2 A I (u 2 -u 1 ), ∂ (h 2 ρ 2 u 2 ) ∂t + ∂ h 2 ρ 2 u 2 2 + h 2 p 2 (ρ 2 ) + 1 2 ρ 2 gh 2 2 ∂x = -ρ 2 gh 2 ∂h 1 ∂x + p 0 ∂h 2 ∂x - Z 1 Z 2 Z 1 + Z 2 A I (u 2 -u 1 ).
(5.3) The numerical advantages of this system are significant. In addition to its simple resolution through the HLL-type Riemann solver (Section 3.2), computational time saving is tremendous.

Z k = ρ k c k
For instance, the computation of the preceding test problem required about 30 hours for the two-dimensional simulation (computing only half of the domain) with a parallel code (MPI) running on 64 cores, while the one-dimensional computation needed only a couple of seconds in sequential implementation.

Conclusion

A pressure relaxation model 6 equations has been built, especially devoted to twolayer shallow water flows. The mathematical structure of the new formulation is well-posed and results in a strictly hyperbolic model. The system considers weak compressibility of the fluids, which is responsible for its hyperbolic behavior, and is shown to tend to the conventional, but conditionally hyperbolic, two-layer shallow water model in the stiff pressure relaxation limit.

A simple, efficient and robust HLL-type Riemann solver has been derived to solve the corresponding non-conservative system. Computational examples have shown capabilities of the present formulation.

Compared to multi-D computations of gravity-driven interfacial flows, the new model offers tremendous numerical advantages and computational savings. This is done at the price of a single parameter in the drag force model.

This research work can be continued in many directions. Among them, the consideration of variable topography, friction with the ground and interfacial area creation through turbulence modeling seem important.

AppendixA. Extended formulation with variable external pressure

To account for variable external pressure, System (2.3) is generalized as,

                                     

. 3 )

 3 Two equations have been added and express the transport of the heights of the fluid layers that are assumed to vary as a function of pressure differentials (p kp 0 ). p k denotes the thermodynamic pressure of fluid k, given by barotropic (and convex) equations of state p k (ρ k ). Example of such equation of state (EOS) is,

4 Figure 1 :

 41 Figure 1: Schematic representation in the (x,t) diagram of the six waves (2.10) present in the flow model (2.3) and emerging at a given initial discontinuity.

Figure 2 :

 2 Figure 2: Schematic representation in the (x,t) diagram of the two extreme waves and the two contact waves considered for the transport of the two heights.

Figure 3 :

 3 Figure 3: Schematic representation in the (x,t) diagram of the two levels (ρ 2 h 2 ) * L,R in the Riemann problem solution.

Figure 4 : 1 2 1 2

 411 Figure4: Comparison of the computed solutions with the present HLL-type Riemann solver (thick lines) versus the exact solution of the one-layer Saint-Venant equations (thin lines and symbols). Results on the left (A) use a 100-cell mesh while results on the right (B) use a 10, 000-cell mesh. The two plots use constant sound speed set to c k = 100 for both fluids. The dashed lines represent the initial conditions. For the sake of clarity, 50 symbols are plotted for the exact solution, shown at time t ≈ 2. First-order Godunov-type numerical scheme is used with CF L = 0.9. Computed results are shown at the same time in full lines. Results at left (A) show a curved line where the shock and expansion waves have already exited the domain. Same computation is rerun with 10, 000 cells (B) restoring the two waves inside the domain. All variables are in dimensionless units.

Figure 5 : 1 2

 51 Figure 5: Influence of the artificial speed of sound. The computed solutions of the present HLL-type Riemann solver for the two-layer system are displayed in thick lines. The exact solution of the one-layer Saint-Venant equations is shown in thin lines and symbols. System (2.3) is solved with various sound speeds c k = θ k 1 2 gh k . Figures C, D, E and F use respectively θ k = 0.1, θ k = 2, θ k = 10 and θ k = 50. The dashed lines represent the initial conditions: h lef t

2 Figure 6 : 1 2

 261 Figure 6: Comparison of the computed solutions with the present HLL-type Riemann solver (symbols) versus the exact solution of the one-layer Saint-Venant equations (thick lines). The numerical system uses an artificial sound speed reading c k = θ k 1 2 gh k with θ k = 2. The dashed lines represent the initial conditions: h lef t

Figure 7 :

 7 Figure 7: Schematic representation of the test problem analyzing the effects of the fluid densities. A fluid layer with a Heaviside profile is set to motion under gravity effects and interacts with the lighter fluid initially above with a Heaviside profile as well.

  examines two different ratios. Plots G and H consider ρ 1 = 1000 kg.m -3 and ρ 2 = 1 kg.m -3 , this situation being typical of water-air configurations while Plots I, J, K and L consider ρ 1 = 1000 kg.m -3 and ρ 2 = 990 kg.m -3 , a situation reminiscent of water-oil flows. Besides, two different values of h 1 are used. The left column of Fig. 8 considers initially h 1 = 4 m while the right column uses h 1 = 10 -6 m. We will see that this initial data influences significantly computed results. All plots of Fig. 8 use h 2 = 20 m with the exception of Plots K and L that use a much larger height for the upper layer, h 2 = 100 m. This data also influences the results.

  repeat the same test (h ′ 1 = 10 m, h 2 = 20 m, h 1 = 4 m (Plot I) and h 1 = 10 -6 m (Plot J)) but with a water-oil-type configuration. As expected, the solutions are different from the one-fluid solutions, as the two fluids are dense and interact each other, this interaction being taken into account by the two-fluid model only. The two-layer solutions are shown with full lines and the single-layer Saint-Venant solution is shown with symbols, just to compare the limit solutions.

Figure 8 :

 8 Figure 8: Influence of the density ratio between the two fluids. Two different density ratios are used: r = ρ2 ρ1 ≪ 1 and r = ρ2 ρ1 → 1. The computed solutions of the present HLL-type Riemann solver for the two-layer system are displayed in thick lines. The exact solutions of the one-layer Saint-Venant equations are shown in thin lines and symbols for comparison. System (2.3) uses an artificial sound speed reading ck = θ k 1 2 gh k with θ k = 2.The dashed lines represent the initial conditions. Final time: t ≈ 1 s. All results use a 1000-cell mesh. Second-order MUSCL-type numerical scheme is used with Sweby's limiter (φ = 1.1, see[START_REF] Sweby | High resolution schemes using flux limiters for hyperbolic conservation laws[END_REF]) and CF L = 0.5. For the sake of clarity, 50 symbols are plotted for the exact solution of the single-layer model. On the column at left, the initial height is h 1 = 4 m and on the column at right h 1 = 10 -6 m. When the density ratio is small: r = ρ2 ρ1 ≪ 1, as in configurations G and H, the two-layer and single-layer models are in perfect agreement. Large differences appear when the density ratio increases as shown in configurations I, J, K and L.

Figure 9 : 1 = 0.5 m, h right 1 = 1 = u right 1 = u lef t 2 = u right 2 = 1 = u right 1 = u lef t 2 = u right 2 = 2 .

 911112211222 Figure 9: Schematic representation of the test problem comparing the computed solutions with the present HLLtype Riemann solver and those computed with the VFRoe method. The initial conditions are h lef t 1 = 0.5 m, h right 1

Figure 10

 10 Figure 10 displays the results obtained with the HLL-type solver and the VFRoe method on a 100-cell mesh. The initial velocity is set to u 1 = u 2 = 0 m.s -1 throughout the entire domain. The results are shown at time t ≈ 0.12 s. The hyperbolic model (2.3) considers compressible fluids during the wave propagation stage.

1 u right 1 = u lef t 2 = u right 2 =

 1122 0 m.s -1 , the velocity u = 0 m.s -1 and the full state vector W in the VFRoe solver induces singularities in the decomposition of the eigenvectors. To lower this effect, different sound speeds (through θ 1 and θ 2 ) have been considered in the various fluids. In the absence of such correction, the VFRoe solver fails immediately. Consequently, θ 1 = θ 2 is used as well with the HLL-type solver for proper comparison.

Figure 11

 11 Figure 11 addresses the same test problem on a 10, 000-cell mesh, showing convergence to the same solution of both VFRoe and HLL. With refined mesh, the interfacial waves are clearly connected to a constant plateau spanning as time evolves. Besides the four-wave structure is clearly observable in the velocity plots.

Figure 11 :h lef t 1 = 0 . 5 m, h right 1 = 0 .45 m, h lef t 2 = 0 . 5 m, h right 2 = 0 .55 m, u lef t 1 = u right 1 = u lef t 2 = u right 2

 1110510205201122 Figure 11: Comparison of the computed solutions with the present HLL-type Riemann solver (thick lines) versus the computed solutions with the VFRoe method (thin lines and symbols). The dashed lines represent the initial conditions: h lef t 1 = 0.5 m, h right 1

Figure 14 : 1 =u lef t 1 = u right 1 = u lef t 2 = u right 2 = 2 .

 14111222 Figure 14: Rate of convergence of the present HLL-type solver. The second-order Godunov-type method is used with Minmod limiter and with four mesh resolutions. The thin solid lines represent results obtained with 100 cells (left column) and 400 cells (right column). The dashed-dotted lines represent results obtained with 200 cells (left column) and 800 cells (right column). The thick solid lines represent the results obtained with 10, 000 cells and considered as "reference" solution. The thick dashed lines represent the initial conditions: h lef t 1 = 0.5 m, h right 1

20 m x 3 Figure 15 :

 315 Figure 15: Schematic representation of the test problem comparing the multidimensional interface approach and the one-dimensional two-layer shallow water model. A fluid layer with a Heaviside profile is set to motion under gravity effects and interacts with the lighter fluid initially above.

  Krypton and air are considered at rest and at atmospheric conditions initially. The flow is set to motion by the gravity acceleration |g| = 10 m.s -2 and consists of two ideal gases evolving at low Mach number. Due to symmetric boundary conditions, only half of the domain is computed with the multidimensional approach. This latter uses an unstructured mesh made of 510, 000 triangular elements with spatial discretization varying from about 1.5 cm in the zone of interest ( ∼ = 480, 000 elements located between y = 0 m and y = 1 m) to about 1 m in the upper far field. The numerical boundaries are considered as atmospheric outflows except for the left side (symmetric condition) and bottom one (flat ground) considered as walls.

Figure 16 :

 16 Figure 16: Comparison of the one-dimensional solution (thin lines and symbols) of the present two-layer shallow water model (HLL-type solver, θ k = 2) versus the two-dimensional computation of the diffuse interface model of Saurel et al. (2009) [24] (thick lines). The exact solution of the one-layer Saint-Venant equations is plotted (dotted lines) as well for comparison. The dashed lines represent the initial conditions. For the sake of clarity, 50 symbols are plotted for the two-layer solution. The column at left shows the results at time t = 5 s and the column at right shows the same results at time t = 8 s. For symmetry reasons, only half of the numerical domain is computed with the two-dimensional simulation. Both computations (two-layer and diffuse interface systems) use CF L = 0.8.

Figure 17 :

 17 Figure 17: Kelvin-Helmholtz instabilities observed during the descent of krypton due to gravity effects, g = -10 m.s -2 . The figure presents krypton volume fraction contours. The black rectangle represents the initial position of the gas. The results are shown at time t = 2 s and are computed with the diffuse interface model of Saurel et al. (2009)[START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF] on an unstructured mesh made of 510, 000 triangular elements. MUSCL-type method is used with the "sharpening-interface" method of[START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF] [START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF] and CF L = 0.8. For symmetry reasons, only half of the numerical domain is computed.

  denotes the acoustic impedance of fluid k and A I denotes the specific interfacial area. ρ 1 = 3.506 kg.m -3 , c 1 = 218 m.s -1 for krypton and ρ 2 = 1.29 kg.m -3 , c 2 = 340 m.s -1 for air are used in this work. This acoustic drag effect modeling was developed in Saurel et al. (2003) [6] and Chinnayya et al. (2004) [25], pages 504 and 510. It is obtained by local interfacial pressure integration over the surface of a piece of interface. The interfacial pressure is estimated through an approximate Riemann solver for the Euler equations of gas dynamics. Let us mention that upon integration over height, as done in Eqs. (5.1) and (5.2), the specific interfacial area A I becomes dimensionless. shown in Figs. 18 and 19, computed results are significantly improved when drag effects are considered. The two figures show the solutions with respectively A I = 0.0003 and A I = 0.00015. A comparison of the different results at time t = 8 s is presented in Fig. 20, showing the influence of the A I parameter. Indeed, the overall solution presents very good agreement with the 2D simulation using the diffuse interface two-phase flow model. The results show the ability of the new two-layer shallow water model to predict the spreading and dispersal of two gases evolving at low Mach number.

Figure 18 :Figure 19 :- 1 )A

 18191 Figure 18: Comparison of the one-dimensional solution (thin lines and symbols) of the present two-layer shallow water model (HLL-type solver, θ k = 2) versus the two-dimensional computation of the diffuse interface model of Saurel et al. (2009)[START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF] (thick lines). Drag effects are included in the shallow water system with A I = 0.0003. For the sake of clarity, 50 symbols are plotted for the two-layer solution. The column at left shows the results at time t = 5 s and the column at right shows the same results at time t = 8 s. For symmetry reasons, only half of the numerical domain is computed with the two-dimensional simulation. Both computations (two-layer and diffuse interface systems) use CF L = 0.8.

Figure 20 :

 20 Figure 20: Influence of the interfacial area in drag force between the two layers. The solutions of the present onedimensional two-layer shallow water system (HLL-type solver, θ k = 2) are displayed in thin lines and symbols. The solution of the two-dimensional computation of the diffuse interface model of Saurel et al. (2009)[START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF] is shown in thick lines. Drag effects are included in the shallow water system with various values of the specific interfacial area A I . For the sake of clarity, 50 symbols are plotted for the two-layer solution. The results are shown at time t = 8 s. For symmetry reasons, only half of the numerical domain is computed with the twodimensional simulation. Both computations (two-layer and diffuse interface systems) use CF L = 0.8. 1D and 2D computations are in good agreement with A I = 0.0002.

Table 1 :

 1 Initial conditions of the test problem analyzing the effects of the fluid densities.

	20	1000	1

Table 2 :

 2 Initial conditions of the test problem comparing the multidimensional interface approach and the onedimensional two-layer shallow water model. γ k represents the thermodynamic polytropic coefficient of fluid k used in the compressible two-phase flow model.

  ∂ (h 1 ρ 1 u 1 ) ∂t + ∂ h 1 ρ 1 u 2 1 + h 1 p 1 (ρ 1 , ρ 2 , h 2 ) + 1 2 ρ 1 gh 2 ∂ (h 2 ρ 2 u 2 ) ∂t + ∂ h 2 ρ 2 u 2 2 + h 2 p 2 (ρ 2 ) + 1 2 ρ 2 gh 2With this formulation, the equation of state of the first (heaviest) fluid transforms to,p 1 = p 0 + ρ 2 gh 2 + c 2 1 ρ 1ρ The gρ 2 h 1 h 2 term present in System (2.3) is now considered in the equation of state of the first fluid. The hydrostatic pressure of the first fluid has been consequently reformulated in the relaxation term of the height equation as well. It is important to check that the hyperbolicity of the reformulated system is still valid. Also, the stiff pressure relaxation limit is meant to recover Ovsyannikov's[START_REF] Ovsyannikov | Two-layer shallow water model[END_REF] system (Eq. (2.1)). These points are examined hereafter.System (A.1) is expressed in primitive-variable formulation (without source terms) as,

									∂W ∂t	+ A(W )	∂W ∂x	= 0,	(A.4)
	∂h 1 ∂t	+ u 1	∂h 1 ∂x	=	µ(p 1 -p 0 -ρ 2 gh 2 ) 1 ρ 1 c 2	,
	∂(h 1 ρ 1 ) ∂t	+	∂ (h 1 ρ 1 u 1 ) ∂x	= 0,
	∂h 2 ∂t	+ u 2	∂h 2 ∂x	=	2 ρ 2 c 2 µ(p 2 -p 0 )	,	∂x			1	= ρ 2 gh 2	∂h 1 ∂x	+ p 0	∂h 1 ∂x	,	(A.1)
	∂(h 2 ρ 2 ) ∂t	+	∂ (h 2 ρ 2 u 2 ) ∂x	= 0,
									∂x				2	= -ρ 2 gh 2	∂h 1 ∂x	+ p 0	∂h 2 ∂x	.
															(0) 1	,	(A.2)
	while it is unchanged for the upper layer,
									p 2 = p 0 + c 2 2 ρ 2 -ρ (0) 2	.	(A.3)

Acknowledgements

The authors are very grateful to Jeaniffer Vides for participating in multiple valuable discussions that definitely helped to improve the quality of this work.

Part of this work has been carried out in the framework of the Labex MEC (ANR-10-LABX-0092) and of the A*MIDEX project (ANR-11-IDEX-0001-02) funded by ANR. Support from CEA Gramat and especially Emmanuel Labedie are particularly acknowledged.

The authors are grateful to the referees for their interesting and constructive comments.

The wave speeds are solutions of det A(W ) -λI = 0 resulting in,

The six real and distinct eigenvalues of System (2.3) are recovered,

(A.7) System (A.1), like System (2.3), is therefore strictly hyperbolic.

Asymptotic limit

Let us now consider the expansion,

The first equation of System (A.1) becomes,

Hence, at leading order,

This last result is introduced in the momentum equation of the first fluid of System (A.1) and yields,

The momentum equation of System (2.1) is thus recovered. Similar manipulations on the second layer momentum equation of the present model lead to,

The second momentum equation of System (2.1) is recovered as well. The present formulation thus tends in the limit of stiff pressure relaxation to the conventional two-layer Saint-Venant model. Thereby, System (A.1) allows to compute solutions of the non-hyperbolic model [START_REF] Ovsyannikov | Two-layer shallow water model[END_REF] with a hyperbolic step using the preceding wave speeds (Eq. (A.7)) followed by a correction step regarding the heights.

This latter is supposed to be isentropic. The equation of state for the second fluid (lightest) reads:

. At the end of the relaxation step, p * 2 = p 0 and implies ρ * 2 = ρ (0) 2 . As the mass quantity is known from the hyperbolic evolution and is constant during the relaxation process, m 2 = h 2 2 = h * 2 ρ * 2 , the relaxation step just consists in reseting the height:

For the first fluid, the EOS now reads:

. At the end of the relaxation process, the pressure reads

1 and the relaxation step reduces to h

as before with System (2.3). The only difference is that ρ

1 no longer represents the density of fluid 1 at atmospheric pressure p 0 but is corrected by hydrostatic effects as ρ

, where ρ


Fenemedre Qaeze
email: fenemedre.qaeze@irit.fr

Romain Guillaume
email: romain.guillaume@irit.fr

Caroline Thierry
email: caroline.thierry@univ-tlse2.fr

ROBUST OPTIMIZATION FOR MINE PLANNING

Keywords: uncertainty, scenario, mine planning, MaxMin criterion 1

In this paper, we are interested in mine planning under uncertainty on the capacity of extraction in a mine complex. Indeed, in real life there is always a gap between planned activities and activities carried out during the period. In this paper, we determine the quantity of loaded mining materials, transported mining materials and processed mining materials in order to minimize the impact of uncertainty into the mining production chain. In order to achieve this challenge, we proposed a model to represent the mining complex. This model takes into account the uncertainty on the capacity by the use of scenarios. In order to face the uncertainty on the capacity of extraction we propose a robust approach with the MaxMin decision criterion. To the best of our knowledge, while the bulk of the literature treats the optimization in mining complex with stochastic approaches, this paper uses a robust approach under uncertainty. A case study using data from a nickel laterite company is used in order to implement the proposed model under uncertainty on capacity.

INTRODUCTION

A mining complex could be considered as a supply chain system where material is transformed from one processing plant to another. Indeed, a typical mining complex "will include a number of mines, a number of processing plants and a number of products" [START_REF] Whittle | Not for the faint-hearted[END_REF]. In this context different kinds of uncertainty in the mining complex have been identified in the literature: the geological uncertainty [START_REF] Montiel | Optimizing mining complexes with multiple processing and transportation alternatives: An uncertainty-based approach[END_REF], which includes the uncertainty on the ore grade [START_REF] Lagos | Robust Planning for an Open-Pit Mining Problem under Ore-Grade Uncertainty[END_REF], the uncertainty on the cut-off grade [START_REF] Azimi | An uncertainty based multi-criteria ranking system for open pit mining cut-off grade strategy selection[END_REF] and the uncertainty on metal [START_REF] Lamghari | A diversified Tabu search approach for the open-pit mine production scheduling problem with metal uncertainty[END_REF]. Besides the uncertainty around the element grade, the uncertainty on the metal price, operating cost [START_REF] Dehghani | Evaluation of the mining projects under economic uncertainties using multidimensional binomial tree[END_REF] and the uncertainty on supply and demand [START_REF] Asad | Implementing a parametric maximum flow algorithm for optimal open pit mine design under uncertain supply and demand[END_REF] have been considered. Different kinds of planning decisions in the optimization of the smelting furnace impact the mine complex in various levels horizons. There are short-term, middle-term and long-term decisions. In this paper, we focus on middle-term decisions. The determination of the open pit limit involves the decision of considering a bloc as a waste or not [START_REF] Dagdelen | Open Pit Optimization -Strategies Of Improving Economics Of Mining Projects Through Mine Planning[END_REF]. The determination of the sequence of blocs involves the decision of considering one or another period for the extraction of blocs [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF]. Finally, the common objective of these decisions is to maximize the profit or to maximize the production.

As we can see, in the bulk of the literature, ore grade uncertainty is well studied. Indeed, one of the principal problem in the optimization of a mining complex is to determine the design of the open pit limit [START_REF] Dagdelen | Open Pit Optimization -Strategies Of Improving Economics Of Mining Projects Through Mine Planning[END_REF], [START_REF] Leite | Stochastic optimization model for open pit mine planning: application and risk analysis at copper deposit[END_REF], [START_REF] Amaya | A scalable approach to optimal block scheduling[END_REF]. The determination of the open pit limit is dependent on the valuation of a mining bloc. The valuation of a bloc is computed from the planned profit on a bloc. This planned profit value takes into account the quantity and quality of products that a bloc could produce and the cost of the operations (extraction, transportation and the transformation). After the determination of the open pit limit, the principal problem is the determination of the sequence of the extracted blocs [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF], [START_REF] Askari-Nasab | Open pit optimisation using discounted economic block values[END_REF], [START_REF] Dimtrakopoulos | Stochastic optimization of mineral value chainsdevelopments and applications for the global optimisation of mining complexes with uncertainty[END_REF]. In this context, the problem is to determine the best period for a set of blocs to be extracted and the best extraction order of these blocs. The period and the order of extraction are factors that could increase generated profit. In this paper, the objective is to maximize the production of the processing plant more precisely of the smelting furnace, which is the critical resource of the processing plant, with minimal cost. This maximization is done knowing the sequence of extracted blocs. In this paper the uncertainty is integrated by taking into account different scenarios. Stochastic approaches have been proposed to address the problem of uncertainty on ore grade [START_REF] Montiel | Optimizing mining complexes with multiple processing and transportation alternatives: An uncertainty-based approach[END_REF], [START_REF] Lagos | Robust Planning for an Open-Pit Mining Problem under Ore-Grade Uncertainty[END_REF] and [START_REF] Azimi | An uncertainty based multi-criteria ranking system for open pit mining cut-off grade strategy selection[END_REF]. In this context of uncertainty, a probability distribution on the state of the world is assessable. Hence this approach consists in considering that the different scenarios have equivalent probabilities. In such, so called, stochastic approaches, the optimization consists in choosing the maximum value in average. To the best of our knowledge, the bulk of the literature goes ahead optimization with stochastic approaches. In this paper, a robust approach is proposed. In decision theory, the following classical decision criterions are: the criterion of Wald (1950) or Maximin criterion, the Maximax criterion, the Hurwicz criterion, the regret minimization criterion or MinMax criterion, the Leximin and Leximax criterions, and the Laplace criterion [START_REF] Kast | La théorie de la décision[END_REF]. Our study focuses on the MaxMin criterion, which minimizes the maximal cost between planned solutions and real solutions. Indeed, in this paper we focus on the uncertainty of capacity of extraction and we try to reduce the impact of the uncertainty on the mining complex.

The remainder of this paper is organized as follows. In the following section, we introduce the context of the study. In the next section, we introduce the proposed deterministic model in deterministic context then the model under uncertainty on capacity. In this section, the notations and assumptions for the models are provided. In the following section, a numerical test is conducted to show the efficiency of the robust approach with the use of the MaxMin criterion. Finally, the conclusions and future works are outlined.

CONTEXT OF THE STUDY

In this study we are interested on the exploitation of the mine complex. The mine complex is composed of a set of mines, a set of processing plants with stock areas, a calcination plant and a smelting furnace. In the figure 1, we illustrate the general context of this study. The extracted blocs from each mine are stored at the seaside. At this place, the mining material could be considered as blended and it is also stated that it is not possible to extract more than the capacity of extraction, which is composed of specifics trucks and specific workforce. From this place, transportation is scheduled in order to supply the processing plant. Such as for the extraction process, it is not possible to transport more mining material than the capacity of transportation which is composed of a set of boats and trucks. The transportation is the link between mines and the processing plant which is composed of stock areas, a calcination plant and a smelting furnace. In this study each mine has an order of bloc extraction, so called: "sequence of extracted blocs" and the smelting furnace is considered as the customer of the system under study.

As we said before, this study focuses on the uncertainty of capacity at the mine. Indeed, in real life, there is a gap between the quantities planned by the mines decision makers and the real extracted quantities. More precisely, we are interested on the uncertainty due to a hurricane which could happen in a short and well known period: for instance in New Caledonia from January to March. Hurricanes stop the extraction (of one mine) and the transportation from this mine to the furnace. We cannot predict the hurricanes so the uncertainty on the occurrence of hurricanes induces an uncertainty on the capacity of the mining process (only the periods are known). Thus, the problem is: how to plan the extraction and the transportation to guarantee a realistic production of the furnace? In other words, the problem is to plan the production of the furnace such that the production of furnace by period will be realizable for all scenarios (due to the uncertainty). There is the possibility to increase the extraction capacity after a hurricane but this induces extra cost. Hence, we look for an extracting, transportation and production plan, which can be adaptable to uncertainty for a minimal cost. The solution will be a balance between safety stocks and corrective decisions (that could increase the extraction capacity).

We have presented the context of the paper: planning optimization under uncertainty on capacity. In the next sections, we propose a deterministic model to represent the proposed optimization problem without uncertainty then a model under uncertainty on the capacity. In this context, we use the MaxMin criterion to optimize the mining complex.

Equations

represents the link between the inventory at the end of a period t, the inventory from precedent period ( the transport to the processing plant inventory (

) and the quantity of calcined mining material (

).

Equations compute the quantity of lost mining material during the calcination process.

Equation [START_REF] Dagdelen | Open Pit Optimization -Strategies Of Improving Economics Of Mining Projects Through Mine Planning[END_REF] represents the capacity of extraction.

Equation ( 9) represents the capacity constraint of the transportation.

Equation

represents the capacity constraint of the smelting furnace.

Approach under uncertainty

In this paper we deal with uncertainty on the extraction capacity. This capacity may decreased due to the meteorological conditions. For instance in New Caledonia from January to march, hurricanes often stop the extraction (of one mine) and the transportation from this mine to the furnace. Thus, the problem is "how to plan the extraction and the transportation to guarantee the well production of the furnace?". In this study we assume that the occurrence of the hurricane cannot be predicted (only the period where it is possible to have one are known). Moreover the extraction and transportation plan can only be adjusted after the hurricane is passed: for instance it is possible to increase the capacity after a hurricane. Indeed, during the hurricane there is a work stoppage, and after the hurricane overtime or temporary staff can be used. Due to the uncertain nature of a hurricane this increase of working hours cannot be planned before the occurrence of the hurricane. Thus, to take the uncertainty into account, different scenarios are considered according to the possible periods of occurrence of the hurricanes and their impacts on the capacity. So in the robust model, we want to have a stable flow through the furnace with maximal value and minimal cost. Hence the decision variable is not scenario dependent but all other variables are scenario dependent. For instance, , becomes as the stock value depends on the scenario. A scenario is a variation of extraction capacity due to a hurricane. Hence we need to distinguish the real quantity of extracted materials and the maximal capacity of extraction. Indeed, in robust model we need to:

-Distinguish the bloc and the period so become where b is the block number as well as the period of extraction of this bloc for the nominal plan ( without damage), -Introduce a new decision variable which is the quantity of bloc b of mines m extracted at period t in scenario s -Add parameter which is the real capacity of extraction at mine m at period t in scenario s.

In the figure 3, we represent the flow graph of the robust model.

In this study, the robust approach consists in choosing the best decision on the extraction, loading, transportation, storing, in order to catch a stable quantity in input for the smelting furnace. Concerning the scenarios in which we have damage and cause an increase of capacity or a decrease of capacity we have levers: to make more stock or to increase workforce with temporary work contract and overtime. The last solution is not use as a rule because it is too expensive for a society. Thus, the levers, which increase capacity, are considered. In the proposed model, the extra capacity is modeled at constraint with the variable: .

In this optimization context the optimization consists in computing solution that takes into account the scenarios. An apparition of damage characterizes a scenario during the horizon. Before the observation of the damage, the mine production system follows the nominal plan without damage (the different scenario dependent variables are equal for all scenarios). After undergo of a damage, a new plan is computed. Hence, the scenario dependent variables for a scenario in which the damage appears are now different than the ones of the nominal plan. The objective is to propose a realizable plan for the smelting furnace for all the scenarios so that the variable linked to the smelting furnace is not scenario dependent.

The figure 4 bellow illustrates the process of computation of new values of the scenario dependent variables at the moment of observation of a damage that is taking into account in the single proposed robust MIP model.

As shown in figure 3, we find two main parts in the model. The first part is dependent of the scenario; this is the part, which deals with the extraction, transportation and blending process. The part, which is not scenario dependent, is the one, which deals with the smelting furnace.

In deterministic context, we only find a planned mining material. Indeed, in deterministic context, planned activity is equal to realized activity.

Model under uncertainty

We have presented the general context of the optimization under uncertainty, in the previous

The application that is discussed is for a lateritemining complex. The laterite-mining complex is localized in a zone where the risk of occurrence of a hurricane is high. This risk is the consequence of the meteorological uncertainty. The purpose of this example is to implement the Maximax criterion.

Description of the input data

In order to implement the model, we took into account a set of twelve periods and two mines.

The proposed robust approach is tested through a series of seven scenarios. The first three scenarios concern the case where the first mine is impacted. The first scenario takes into account a hurricane emerging in March, the second scenario with a hurricane emerging in April and the third scenario takes into account a hurricane emerging in May. As we consider a set of two mines for the tests, there are also three scenarios for the second mine according to the period of occurrence of the hurricane (March, April or May). Moreover, we consider a scenario in which there is no hurricane observed.

Description of the experimental process

The objective of the experimental process is to compare the robust approach and the classical approach (re-planning).

In real life, we usually find the classical approach. This approach consists in a replanning when damage occurs.

The robust approach consists in the computation of a production plan that takes into account a set of scenarios. The tests are done as follow. First, we compute a production plan under uncertain context by taking into account a set of scenarios. Secondly, we generate a set of production plans, on the one hand with re-planning and on the other hand by taking into account each scenario. The re-planning has been done as follow. We consider a set of damages with a period of occurrence of these damages. As shown in figure 4, first, we generate a production plan in a context without damages, then we consider the occurrence of a first damage. Before this first damage, the variables have the same values as in a production plan without damages. After the occurrence of damage we re-planned the production and we have new values for the variables. The variables linked directly to the smelting furnace are not scenario dependent, they are stable. The first period of occurrence of a damage is t=5. We save the generated solution for this first replanning process. Another re-planning process is done for a damage at period t=6 ant at period t=7.

Finally, we choose the production plan that generates the maximum of profit and we compare this solution to the robust solution.

In the next section, we analyze the obtained results.

Analysis of the result

The different values of the input data are presented is figure 5. We have compared the profits generated by a robust approach and by the different re-planning (that takes into account each different damages).

The results show that for damage at period t=5 (hurricane in December), robust approach is better at 77%, for a damage at t=6 (hurricane in January) robust approach is better at 52% and for a damage at t=7 (hurricane in February) robust approach is better at 37%.

As we said before the best result is obtained with damage in February (t=7). We compare the classical approach (with re-planning, figure 6) and robust approach (MaxiMin criterion, figure 7). The figures bellow, show the evolution of the stock at the processing plant and the evolution of the extra-capacity for a damage in February. We can see that for the robust approach, we have a better stability in the consumption of the stock. This paper is the first paper under our knowledge that applies such a robust approach in order to optimize a mining complex. In this purpose, the MaxMin criterion has been used. The proposed approach is a new way of thinking the optimization of a mining complex. Indeed, with this decision criterion, decision maker belief is taken into account. Indeed, the different scenarios that have been tested result from the industry expertise. In the future, we propose to develop a model that deals with the grade of element. Indeed, in real life, the concentration of an element influences the fusion speed in the smelting furnace. As other perspective, we will apply other decision criterions as min max Regret or leximin in order to optimize a mining complex and will make comparison between them.

Figure 1 ::

 1 Figure 1: Context of the study

Figure 2 :

 2 Figure 2: Deterministic approach

Figure 4 :

 4 Figure 4: Illustration of the re-planning

Figure 5 :

 5 Figure 5: Input data

Figure 6 :

 6 Figure 6: Stock and extra capacity in February; classical approach

Figure 7 :

 7 Figure 7: Stock and extra capacity in February; robust approach

section; we introduce now the proposed robust model. The presentation of the model proceeds as follow. First, we give the context of the optimization by introducing the notations of new variables and parameters and finally we introduce the objective function and the constraints of the different constraints.

Constraint represents the link between the quantity (

) of removable mining material and the effectively extracted in several periods for a scenario s.

Constraints ( 15) and (17) express that before the observation of damage there is the same extracted and loaded materials.

represents the constraint of capacity for the real extraction of mining material for a scenario s.

Hurricane. Constraint (18) shows the link between loaded mining materials (

) and transported mining materials ( ) from seaside to the processing plant. Indeed, all the loaded mining materials are transported.

Constraints show the link between the quantity (

) of extracted mining material , inventory of the end of the precedent period ( and the outflow ( at the loading zone for a scenario s.

Constraints

represent the link between the inventory ( at the end of a period t and for a scenario s, the inventory ( from precedent period the transport ( ) to the processing plant inventory and the quantity ( ) of calcined mining material.

Constraints

compute the quantity of lost mining material during the calcination process at period t.

Equation

represents the capacity of the smelting furnace at period t.

We have presented in the previous section the deterministic model and the model under uncertainty. In the next section, we illustrate the implementation of these models on data from a nickel industry. We describe the input data and finally we analyze the output data and the results.