Derivation of a bilayer model for shallow water equations with viscosity. Numerical validation

Derivation of a bilayer model for Shallow Water equations with viscosity. Numerical validation G. Narbona-Reina∗ , J.D.D. Zabsonr´e † , E.D. Fern´andez-Nieto1 and D. Bresch ‡ Abstract In this work we present a new two-dimensional bilayer Shallow-Water model including viscosity and friction effects on the bottom and interface level. It is obtained following [[6]] from an asymptotic analysis of non-dimensional and incompressible Navier-Stokes equations with hydrostatic approximation. In order to obtain the viscosity effects into the model we must have into account a second order approximation. To evaluate this model we perform two numerical tests consisting of an internal dam-break problem for both, one and two dimensional cases. In the first one we make a comparison between the model obtained and the Navier- Stokes simulation. Shallow Water equations, bilayer models, viscosity, friction, capillarity, Finite Volume methods. 1 Introduction The goal of this paper is the derivation of a new viscous bilayer Shallow Water model. We also present some numerical test with the aim of checking its validity. So we shall compare the solution obtained by this new model with the Navier-Stokes solution. The Shallow Water (SW) equations are usually used to simulate a large number of geophysical and engineering applications as ocean circulation, coastal areas, rivers, etc. But sometimes these equations are not sufficient to model specific situations as, for instance, the flow involved in the Strait of Gibraltar. In this physical domain two layers of water with different properties are founded, the denser Mediterranean and the Atlantic water. So in this case we must consider a model of at least two layers. For this purpose we can find several derivations of one and two layers SW models to tackle these kind of situations in one and two dimensional case. Usually there are two effects that are neglected in Shallow Water models: the viscosity and the Coriolis force. On the one hand it is important to remark that the viscous effects are essential to obtain an accurate approximation in dam-breaks problems or hydraulic jumps situations, as we can see in [[6, 9]]. On the other hand the Coriolis force plays an important role in geophysical fluid dynamics applications. To include this force in equations does not add any meaningful difficulty for the deduction of the model. Nevertheless its effects are of major importance in these kind of applications, see for example [[8, 14]]. The objective is to find a bilayer model as complete as possible that takes the Coriolis force into consideration, includes the viscosity effect, friction -at the bottom and at the interface- and tension effect -on the surface and the interface-. ∗Dpto. Matem´atica Aplicada I, E.T.S. Arquitectura. Universidad de Sevilla. Avda. Reina Mercedes 2. 41012 Sevilla, Spain, ([email protected], [email protected]) †Universit´e Polytechnique de Bobo-Dioulasso, 01BP 1091 Bobo 01, [email protected] ‡LAMA, UMR5127 CNRS, Universit´e de Savoie, 73376 Le Bourget du Lac, France, [email protected] 1 The pioneer work performed in [[6]] has been considered as a basis to develop the deduction of our model. In this work a viscous Shallow Water model is obtained for the one-dimensional one-layer case by performing the asymptotic analysis of the Navier-Stokes equations where friction effect at the bottom has been taken into account. When only first order approximation is considered the viscous terms doesn’t appear in the equations, so moving on to the second order is needed to get the viscosity effects. The authors also put in evidence the difference between two approximations through an application to a dam-break problem. In [[9]], a viscous one layer 2D Shallow-Water system is derived. The originality of that work is the introduction of a surface-tension term through the capillary effects at the free surface and quadratic friction term at the bottom. These terms have been useful to establish the existence of global weak solutions in [[2]]. With regard to bilayer models we must mention the work performed in [[12]] and those developed in [[1]]. In [[12]], it is deduced a bilayer viscous Shallow-Water model which takes into account the friction at the interface. But instead of asymptotic analysis development, several simplifications are used in the boundary conditions to deduce the final system. The energy of the system is also obtained under restrictive hypothesis on the data. In [[1]], a derivation of a multi-layer Shallow-Water model is performed to extend the case of one layer established in [[6]]. In this work, using the hydrostatic pressure and the kinematic boundary conditions, the author derives momentum equations of the form: ∂tZHα Hα−1 udz +∂xZHα Hα−1 u2dz +ghα∂xh=ν0 ǫ∂zu(Hα(t, x)) −ν0 ǫ∂zu(Hα−1(t, x)) and use at the leading order a finite difference method with respect to the vertical variable when the equation is an interface equation to deduce the friction term: µ∂zu(Hα) = µUα+1 −Uα hα+1 +hα .(1.1) Another works related to the derivation of 2D Shallow-Water model can be find in [[4]] and in [[5]]. In [[4]] the authors include the atmospheric pressure in the derivation. In [[5]] a non-viscous two-layer Shallow Water system is deduced following [[6]]. Linear friction conditions have been taken into account on the interface and on the bottom. The deduction of the bilayer model developed in the present work has been obtained by integrating the three dimensional Navier-Stokes equations with Coriolis force and by using the asymptotic analysis to get the viscosity effects. We have also considered the friction effects on the surface and at the bottom and the tension effect on the surface and at the interface. We would like to remark that the friction term usually carries some difficulties in its treatment (both in the model deduction and for the proof of the existence of solution). So we can find several definitions to avoid these troubles. Often one takes non linear friction terms, for example in [[13]] the tangential displacement has been considered in the definition of the friction traction or a more complex expression can be found in [[10]] where a microscopic study has been tackled to set it. In [[9]] a non linear friction under the form −γ|v|αvis considered at the bottom. So, in a equivalent way we could define the interface friction as −γ|v1−v2|(v1−v2). In [[16]] we have performed the theoretical ana-lysis of the model presented in this work, obtaining the existence of weak solutions. In order to get it we have not been able to control the interface friction of quadratic form, so we have defined fric(v1, v2) = −ξB(h1, h2)(v1−v2), B(h1, h2) = h1h2 ρ1 ρ2h1+ρ2 ρ1h2 ,(1.2) with ρithe density of each layer and ξa positive constant. So to be in agreement with these results we have considered in this paper the linear definition of the friction term (1.2). The drag coefficient Bis also used in [[7]] and allows us to control the friction term. Note that in [[3]] the authors study a system of 3D Navier-Stokes equations in a two-layer thin 2 domain with an interface condition (νi∂3ui j−k(u1 j−u2 j))|x3=0 = 0 i, j = 1,2.(1.3) This condition is the same kind condition appearing in the Primitive Equations of the Coupled Atmosphere and Ocean which describes the atmosphere/ocean interaction. The paper is organized as follows: In Section 2 we develop the derivation of the model. First we write the equations in non-dimensional variables, later we perform the hydrostatic approximation to obtain the Shallow Water equations and finally we study the asymptotic analysis of two layers. In the last part of the section we state the final systems found and we include some remarks about them. Section 3 is devoted to show some numerical test in which we notice the improvement obtained in the solution when considering viscous Shallow Water model in front of considering the first order one. In Test 1 we compare them with a solution of the Navier-Stokes equations for an internal dam-break problem in one-dimensional case. In the second test we present a 2D circular dam-break problem comparing the solutions of the models deduced in this work. 2 Derivation of the bilayer viscous Shallow Water model In this section we perform the derivation of the model proposed in this paper. We start from the Navier-Stokes equations in a periodic domain D(t)∈R3. We consider a two layer environment of inmiscible fluids including three boundary regions. We assume that the bottom is defined by the function b(x), and we denote by η1,2(t, x) the interface and the free surface is given by η(t, x). The vertical direction is denoted by zand by xwe denote a point in a domain χ⊂R2. So the global domain is D(t) = D1(t)∪D2(t)∪Γb∪ Γ1,2(t)∪Γs(t), being: D1(t) = {(x, z)∈R3/x ∈χ, b(x)< z < η1,2(t, x)}; D2(t) = {(x, z)∈R3/x ∈χ, η1,2(t, x)< z < η(t, x)}; Γb={(x, z)∈R3/x ∈χ, z =b(x)}; Γ1,2={(x, z)∈R3/x ∈χ, z =η1,2(t, x)}; Γs={(x, z)∈R3/x ∈χ, z =η(t, x)}; (2.1) From now on, subscript 1 will correspond to the layer located below and subscript 2 to those located on the top. We denote by h1(t, x) = η1,2(t, x)−b(x) the thickness of the layer 1 and by h2(t, x) = η(t, x)−η1,2(t, x) the thickness of the second one. See Fig. 1. We consider ui= (vi, wi) the velocity of each layer, ρithe density, µidenote the dynamic viscosity and piis the pressure, for i= 1,2. With this notation, the Navier-Stokes equations for each layer i= 1,2 state as: ρi∂tui+ (ρiui∇)ui−div (σi) + 2ρi−→ Ω×ui=−ρigez; div (ui) = 0.(2.2) We have included the Coriolis force given by the term 2ρi−→ Ω×uiwhere −→ Ω = Ω(0,cos θ, sin θ), being θthe latitude. The stress tensor is defined as σi= 2µiD(ui)−piId where D(u) = ∇u+∇⊥u 2is the strain tensor. Finally gis the gravity constant. We denote by subscript nthe normal component and by τthe tangent component, that is, d=dnn+dττfor all d∈R2. In order to obtain a well-posed system we impose the following conditions on the boundaries: •On the free surface z=η(t, x): We assume the pressure to be constant. One usually neglect the atmospheric pressure effect but here we have considered it. If we denote by nsthe unit normal vector to η(t, x) towards the increasing zand by κthe mean curvature of the surface (κ= div (ns)), the tension effect on the surface is given by: σ2·ns=α2κ·ns,(2.3) 3 h h 2 1 b(x) η1,2 η= b+h = b+h1 +h 12 z x Figure 1: Domain. being α2constant. Finally we impose the kinematic condition for the surface: ∂tη+v2·∇xη=w2.(2.4) •At the interface z=η1,2(t, x): First we consider the conditions related to the movement of the interface that is advected by the two flows, i.e., the kinematic condition that we write as: ∂tη1,2+vi·∇xη1,2=wi;i= 1,2.(2.5) At the interface we consider the friction effects between the two layers and the continuity of the tension force. These conditions concern the tangent and normal components of the stress term respectively in the following sense. We consider the friction term between the two layers with coefficient γdone by −γ(u1−u2) , so we impose: (σi·n1,2)τ=−γ(u1−u2)τ;i= 1,2.(2.6) We now consider κ1,2= div (n1,2) the mean curvature at the interface being n1,2the unit normal vector to Γ1,2pointing from layer 1 to layer 2, so we have: (σ1·n1,2)n= (σ2·n1,2)n+ ((α1−α2)κ1,2·n1,2)n,(2.7) being α1constant. •At the bottom z=b(x): We consider a Navier condition with a friction coefficient α: (σ1·nb)τ=α(u1)τ,(2.8) and a no-penetration condition: u1·nb= 0; (2.9) being nbthe unit normal vector to Γbpointing to the increasing z. To obtain the model, first we shall write these equations under a non-dimensional form, secondly we shall develop the vertical integration in each layer to obtain the Shallow Water system. Finally we shall perform the asymptotic analysis studding both first and second order approximations. Therefore, two models are proposed depending if the viscous and friction terms are included or not. 4 2.1 Dimensionless equations Before changing equations to non-dimensional variables, we shall write the Navier-Stokes equations and the boundary conditions in a explicit form. First we look at the Coriolis term that reads: ~ Ω×ui= Ω(wicos θ−(visin θ)e2, visin θe1,−vicos θe1).(2.10) If we write the equations for each component of the velocity, we have: ρi∂tvi+ρivi∇xvi+ρiwi∂zvi+ 2ρiΩ(wicos θe1+ (vi)⊥sin θ)− −2µidiv x(Dx(vi)) −µi∂2 zvi−µi∇x(∂zwi) + ∇xpi= 0, i = 1,2; (2.11) ρi∂twi+ρivi∇xwi+ρiwi∂zwi−2ρiΩvicos θe1−µi∆xwi− −µi∂z(div xvi)−2µi∂2 zwi+∂zpi=ρig, i = 1,2; (2.12) and div xvi+∂zwi= 0 i= 1,2.(2.13) Now, we explicitly write the boundary conditions concerning tension and friction terms using the definition of the stress tensor when needed. •Free surface. We take the normal vector ns=1 p1 + (∇xη)2„−∇xη 1«, so the tension condition state at follows: 8 < : (−2µ2Dx(v2) + p2+α2κ)∇xη+µ2(∇xw2+∂zv2) = 0; −µ2(∇xw2+∂zv2)∇xη+ 2µ2∂zw2−p2−α2κ= 0. (2.14) •Interface. In this case the normal vector is n1,2=1 p1 + (∇xη1,2)2„−∇xη1,2 1«, obtaining for the tension condition: (2µ1Dx(v1)−p1)|∇xη1,2|2−2µ1(∇xw1+∂zv1)∇xη1,2+ +2µ1∂zw1−p1= (2µ2Dx(v2)−p2+ (α1−α2)κ1,2)|∇xη1,2|2− −µ2(∇xw2+∂zv2)∇xη1,2+ 2µ2∂zw2−p2+ (α1−α2)κ1,2. (2.15) For the friction condition we must consider the tangent vector ~τ = (τ1, τ2), that we define as follows: τ1=1 |∇xη1,2|„∇⊥ xη1,2 0«, τ2=1 p|∇xη1,2|2+|∇xη1,2|4„−∇xη1,2 −|∇xη1,2|2«.(2.16) So, the conditions for i= 1,2 are: 1 p1 + |∇xη1,2|2µi(∇xwi+∂zvi) = −γ(v1−v2); 1 p1 + |∇xη1,2|2`−2µiDx(vi)∇xη1,2+µi(∇xwi+∂zvi)(1 −|∇xη1,2|2)+ +2µi∂zwi∇xη1,2´=γ((v1−v2) + (w1−w2)∇xη1,2). (2.17) •Bottom. In the same way, we consider the normal vector nb=1 p1 + |∇xb|2„−∇xb 1«. So for the no-penetration condition we obtain: −v1∇xb+w1= 0.(2.18) 5 For the other one we take ~τ = (τ1, τ2) being τ1=1 |∇xb|„∇⊥ xb 0«, τ2=1 p|∇xb|2+|∇xb|4„−∇xb −|∇xb|2«,(2.19) and the condition state as: 1 p1 + |∇xb|2µ1(∇xw1+∂zv1) = αv1; (2.20) 1 p1 + |∇xb|2`2µ1Dx(v1)∇xb−µ1(∇xw1+∂zv1)(1 −|∇xb|2)− −2µ1∂zw1∇xb´=−α(v1+w1∇xb). (2.21) We introduce now a small parameter ǫ=H Lwhere Hand Lare two characteristics dimensions along the edges OZ and OX respectively. We also introduce some others characteristic dimensions: Vfor the horizontal velocity, W=ǫV for the vertical component of the velocity and P=V2for the pressure. Next, we consider the following dimensionless variables, we use the overline notation to denote them: x=L¯x z =H¯z vi=Vviwi=ǫV wi t=L V¯ t pi=V2pi Rei=V L µi Ro =V 2LΩFr =V √gH (2.22) γ=Vγ α =V α κ=ǫ Lκ κ1,2=ǫ Lκ1,2 αi=V2L ǫαi;i= 1,2 b=Hb. (2.23) where we have denoted by Re the Reynolds number, Ro the Rossby number and Fr the Froude number. Thus, the equations get as follows (we drop the “overline” notation for the sake of clarity): ρi∂tvi+ρivi∇xvi+ρiwi∂zvi+ρi 1 Roǫwicos θe1+ρi 1 Ro(vi)⊥sin θ− −2 Rei div x(Dx(vi)) −1 Rei 1 ǫ2∂2 zvi−1 Rei∇x(∂zwi) + ∇xpi= 0; (2.24) ρi∂twi+ρivi∇xwi+ρiwi∂zwi−ρi 1 ǫ 1 Rovicos θe1−1 Rei ∆xwi− −1 ǫ2 1 Rei ∂z(div xvi)−21 ǫ2 1 Rei ∂2 zwi+1 ǫ2∂zpi=−ρi 1 ǫ2 1 Fr2; (2.25) div xvi+∂zwi= 0.(2.26) In the same way the boundary conditions must be modified, we specify them next. 6 •Conditions at the free surface: ∂tη+v2·∇xη=w2; (2.27) „−21 Re2 Dx(v2) + p2+α2κ«∇xη+1 Re2∇xw2+1 ǫ2 1 Re2 ∂zv2= 0; (2.28) „−ǫ21 Re2∇xw2−1 Re2 ∂zv2«∇xη+ 2 1 Re2 ∂zw2−p2−α2κ= 0.(2.29) •Conditions at the interface: ∂tη1,2+v2·∇xη1,2=w2; (2.30) ∂tη1,2+v1·∇xη1,2=w1; (2.31) −„2ǫ21 Re1∇xw1+ 2 1 Re1 ∂zv1«∇xη1,2+ 2 1 Re1 ∂zw1−p1= =−„2ǫ21 Re2∇xw2+ 2 1 Re2 ∂zv2«∇xη1,2+ 2 1 Re2 ∂zw2−p2+ +(α1−α2)κ1,2; (2.32) 1 Rei (∇xwi+1 ǫ2∂zvi) = −1 ǫγ(v1−v2)p1 + ǫ2|∇xη1,2|2; (2.33) −21 Rei Dx(vi)∇xη1,2+1 Rei (∇xwi+1 ǫ2∂zvi)(1 −ǫ2|∇xη1,2|2)+ +2 1 Rei ∂zwi∇xη1,2= =1 ǫγ((v1−v2) + ǫ2(w1−w2)∇xη1,2)p1 + ǫ2|∇xη1,2|2. (2.34) •Conditions at the bottom: −v1∇xb+w1= 0; (2.35) 1 Re1 (∇xw1+1 ǫ2∂zv1) = 1 ǫα v1p1 + ǫ2|∇xb|2; (2.36) 21 Re1 Dx(v1)∇xb−1 Re1 (∇xw1+1 ǫ2∂zv1)(1 −ǫ2|∇xb|2)− −21 Re1 ∂zw1∇xb=−α1 ǫ(v1+ǫ2w1∇xb)p1 + ǫ2|∇xb|2. (2.37) 2.2 Hydrostatic approximation We assume ǫto be small and we take the formal expression of the system at O(ǫ2), keeping the terms of order zero and one. Thus, the hydrostatic system state as follows: ρi∂tvi+ρidiv x(vi⊗vi) + ρi∂z(viwi) + ǫρi 1 Rowicos θe1+ρi 1 Ro(vi)⊥sin θ− −2 Rei div x(Dx(vi)) −1 ǫ2 1 Rei ∂2 zvi−1 Rei∇x(∂zwi) + ∇xpi= 0; (2.38) −ǫρi 1 Rovicos θe1−1 Rei ∂z(div xvi)−21 Rei ∂2 zwi+∂zpi=−ρi 1 Fr2; (2.39) div xvi+∂zwi= 0.(2.40) And boundary conditions: 7 •Conditions at the free surface: ∂tη+v2∇xη=w2; (2.41) „−21 Re2 Dx(v2) + p2+α2κ«∇xη+1 Re2∇xw2= 0; (2.42) ∂zv2=O(ǫ); (2.43) −1 Re2 ∂zv2∇xη+ 2 1 Re2 ∂zw2−p2−α2κ= 0.(2.44) •Conditions at the interface: ∂tη1,2+v2·∇xη1,2=w2; (2.45) ∂tη1,2+v1·∇xη1,2=w1; (2.46) −21 Re1 ∂zv1∇xη1,2+ 2 1 Re1 ∂zw1−p1= =−21 Re2 ∂zv2∇xη1,2+ 2 1 Re2 ∂zw2−p2+ (α1−α2)κ1,2; (2.47) 1 Rei ∂zvi=ǫfric(v1, v2); (2.48) −21 Rei Dx(vi)∇xη1,2+1 Rei (∇xwi+1 ǫ2∂zvi) + 2 1 Rei ∂zwi∇xη1,2− −1 Rei ∂zvi|∇xη1,2|2=−1 ǫfric(v1, v2). (2.49) Where we have denoted by fric(v1, v2) = −γ(v1−v2) the friction term between the two layers. •Conditions at the bottom: −v1∇xb+w1= 0; (2.50) 1 Re1 ∂zv1=ǫα v1; (2.51) 21 Re1 Dx(v1)∇xb−1 Re1 (∇xw1+1 ǫ2∂zv1) + 1 Re1 ∂zv1|∇xb|2− −21 Re1 ∂zw1∇xb=−1 ǫα v1. (2.52) 2.3 Shallow Water system and asymptotic analysis To obtain the Shallow Water equations, we assume that the height is small with respect to the length of the domain, that is, ǫ≪1. We first integrate the equations for each layer. Then we shall perform the asymptotic analysis of the system by introducing an asymptotic regime hypotheses over the physical data. We shall obtain the system at first order, but we will must analyze the second order to obtain a system with viscosity. We first perform the integration of the layer 1. We want to obtain the pressure value, so from equation (2.39) for i= 1: ∂zp1=−ρ11 Fr2+ǫρ11 Ro cos θv1e1+1 Re1 ∂z(div xv1) + 2 1 Re1 ∂2 zw1.(2.53) 8 To get p1, we integre this equations from zto η1,2, with z∈(b, η1,2), p1(z)−p1(η1,2) = −ρ11 Fr2(z−η1,2) + ǫρ11 Ro cos θe1Zz η1,2 v1dz+ +1 Re1 (div xv1−div xv1|z=η1,2) + 2 1 Re1 (∂zw1−∂zw1|z=η1,2). (2.54) By using the divergence free condition we obtain the following expression for p1: p1(z) = p1(η1,2)−ρ11 Fr2(z−η1,2) + ǫρ11 Ro cos θe1Zz η1,2 v1dz− −1 Re1 (div xv1−div xv1|z=η1,2). (2.55) Now, we integrate the system equations from bto η1,2. For the equation (2.40) and using conditions (2.46) and (2.50), we have: ∂tη1,2+ div xZη1,2 b v1dz = 0.(2.56) If we integrate the equation for the horizontal velocity (2.38), and we use the condition (2.46) it gives: ρ1∂tZη1,2 b v1dz +ρ1div xZη1,2 b (v1⊗v1)dz +ǫρ1 1 Ro cos θe1Zη1,2 b w1dz+ +ρ11 Ro sin θZη1,2 b (v1)⊥dz −21 Re1 div xZη1,2 b Dx(v1)dz +∇xZη1,2 b p1dz+ +1 ǫ2 1 Re1 ∂zv1|z=b+1 Re1∇xw1|z=b−p1(η1,2)∇xη1,2+ +2 1 Re1 Dx(v1)|z=η1,2∇xη1,2−1 ǫ2 1 Re1 ∂zv1|z=η1,2−1 Re1∇xw1|z=η1,2− −21 Re1 Dx(v1)|z=b∇xb+p1(b)∇xb= 0. (2.57) Due to condition (2.49), we can write 21 Re1 Dx(v1)|z=η1,2∇xη1,2−1 ǫ2 1 Re1 ∂zv1|z=η1,2−1 Re1∇xw1|z=η1,2= =1 ǫfric(v1, v2) + 2 1 Re1 ∂zw1|z=η1,2∇xη1,2−1 Re1 ∂zv1|z=η1,2|∇xη1,2|2. (2.58) And thanks to condition (2.52), we have −21 Re1 Dx(v1)|z=b∇xb+1 ǫ2 1 Re1 ∂zv1|z=b+1 Re1∇xw1|z=b= =1 ǫα v1|z=b+1 Re1 ∂zv1|z=b|∇xb|2−21 Re1 ∂zw1|z=b∇xb. (2.59) So, we get for the first layer the following equation: ρ1∂tZη1,2 b v1dz +ρ1div xZη1,2 b (v1⊗v1)dz +ǫρ1 1 Ro cos θe1Zη1,2 b w1dz+ +ρ11 Ro sin θZη1,2 b (v1)⊥dz −21 Re1 div xZη1,2 b Dx(v1)dz +∇xZη1,2 b p1dz− −p1(η1,2)∇xη1,2+p1(b)∇xb+1 Re1 (2∂zw1−∂zv1∇xη1,2)|z=η1,2∇xη1,2− −1 Re1 (2∂zw1−∂zv1∇xb)|z=b∇xb+1 ǫfric(v1, v2) + 1 ǫα v1|z=b= 0. (2.60) 9 Subtracting the last two expressions we get: v1−v2=−ǫ 6µ01 h0 1α0ev1|z=b+„1 + ǫ 3γ0„h0 1 µ01 +h0 2 µ02 ««(ev1−ev2)|z=η1,2+O(ǫ2).(2.114) We solve the following system for obtaining ev1|z=band ( ev1−ev2)|z=η1,2: 8 > > < > > : v1=„1 + ǫ 3µ01 α0h0 1«ev1|z=b−ǫ 6µ01 h0 1γ0(ev1−ev2)|z=η1,2; v1−v2=−ǫ 6µ01 h0 1α0ev1|z=b+„1 + ǫ 3γ0„h0 1 µ01 +h0 2 µ02 ««(ev1−ev2)|z=η1,2. (2.115) The results are: ev1|z=b=Dδv1+Dǫγ0h0 1 6µ01 (v1−v2),(2.116) and (ev1−ev2)|z=η1,2=Dǫα0h0 1 6µ01 v1+Dβ−1(v1−v2),(2.117) being β=„1 + ǫ 3µ01 α0h0 1«−1 , δ = 1 + ǫ 3γ0„h0 1 µ01 +h0 2 µ02 «(2.118) and D=„β−1δ−ǫ2 36µ2 01 (h0 1)2α0γ0«−1 .(2.119) As we can check that β−1δ > 1 and having into account that we have a second order approximation, we set: D=βδ−1.(2.120) So we can rewrite the solutions as: ev1|z=b=βv1+βδ−1ǫγ0h0 1 6µ01 (v1−v2),(2.121) and (ev1−ev2)|z=η1,2=βδ−1ǫα0h0 1 6µ01 v1+δ−1(v1−v2).(2.122) We take these values to equations (2.60) and (2.67): ρ1∂t(f h1v1) + ρ1div x(f h1v1⊗v1)−ǫ 2ρ11 Ro cos θe1(f h1)2div xv1+ +ρ11 Ro sin θf h1(v1)⊥−2ǫµ01div x(f h1Dx(v1)) + 1 2ρ11 Fr2∇x(f h1)2− −ǫ 2ρ11 Ro cos θe1∇x((f h1)2v1) + ρ11 Fr2f h1∇xeb+ +ǫρ11 Ro cos θf h1ˆ(v1·∇xeb)e1−(v1·e1)∇xeb˜+ +δ−1γ0„βǫα0h0 1 6µ01 v1+ (v1−v2)«+βα0„v1+δ−1ǫγ0h0 1 6µ01 (v1−v2)«− −2ǫµ01∇xf h1div xv1+f h1∇x(ep1(gη1,2)) = 0. (2.123) 16 ρ2∂t(f h2v2) + ρ2div x(f h2v2⊗v2)−ǫ 2ρ2 1 Ro cos θe1(f h2)2div xv2+ +ρ21 Ro sin θf h2(v2)⊥−2ǫµ02div x(f h2Dx(v2)) + 1 2ρ21 Fr2∇x(f h2)2− −ǫ1 2ρ2 1 Ro cos θe1∇x((f h2)2v2)−2ǫµ02∇x(f h2div xv2)+ +ǫρ2 1 Ro cos θe1h0 2w0 2|z=η1,2−δ−1γ0„βǫα0h0 1 6µ01 v1+ (v1−v2)«+ +∇xgη1,2(ep2(gη1,2) + 2ǫµ02div xv2)− −ǫα02∇x(h0 2∆xη0) + ǫα02∇xη0∆xη0= 0, (2.124) where we can check now how the coefficients for both friction terms have been corrected for two layers. In order to simplify the last two terms in (2.123), we use the interfaz condition (2.32) to write: ep1(gη1,2) + 2ǫµ01div xv1=ep2(gη1,2) + 2ǫµ02div xv2−ǫ(α01 −α02)κ1,2.(2.125) So if we also include the expression of κ1,2we get: −2ǫµ01∇xf h1div xv1+f h1∇x(ep1(gη1,2)) = −2ǫµ01∇xf h1div xv1+ +f h1∇x(−2ǫµ01div xv1+ep2(gη1,2) + 2ǫµ02div xv2−ǫ(α01 −α02)∆xgη1,2) = =−2ǫµ01∇x(f h1div xv1) + f h1∇x(ep2(gη1,2) + 2ǫµ02div xv2)− −ǫ(α01 −α02)f h1∇x(∆xgη1,2). (2.126) Finally the equation for layer 1 reads as: ρ1∂t(f h1v1) + ρ1div x(f h1v1⊗v1)−ǫ 2ρ11 Ro cos θe1(f h1)2div xv1+ +ρ1 1 Ro sin θf h1(v1)⊥−2ǫµ01div x(f h1Dx(v1)) + 1 2ρ1 1 Fr2∇x(f h1)2− −ǫ 2ρ11 Ro cos θe1∇x((f h1)2v1) + ρ11 Fr2f h1∇xeb−2ǫµ01∇x(f h1div xv1)+ +ǫρ1 1 Ro cos θf h1ˆ(v1·∇xeb)e1−(v1·e1)∇xeb˜+ +δ−1γ0„βǫα0h0 1 6µ01 v1+ (v1−v2)«+βα0„v1+δ−1ǫγ0h0 1 6µ01 (v1−v2)«+ +f h1∇x(ep2(gη1,2) + 2ǫµ02div xv2)−ǫ(α01 −α02)f h1∇x(∆xgη1,2) = 0. (2.127) In the same way we work on the last terms in (2.124) by using (2.62) to rewrite them as: ep2(gη1,2) + 2ǫµ02div xv2=ρ21 Fr2f h2−ǫρ21 Ro cos θe1f h2v2−ǫα02κ. (2.128) So: ∇gη1,2(ep2(gη1,2) + 2ǫµ02div xv2) = =∇gη1,2(ρ21 Fr2f h2)−ǫρ21 Ro cos θe1∇x(gη1,2)f h2v2−ǫα02∇xη0 1,2∆xη0.(2.129) 17 We can finally write the equation for layer 2 as ρ2∂t(f h2v2) + ρ2div x(f h2v2⊗v2)−ǫ 2ρ2 1 Ro cos θe1(f h2)2div xv2+ +ρ21 Ro sin θf h2(v2)⊥−2ǫµ02div x(f h2Dx(v2)) + ρ21 Fr2(f h2∇xf h2+f h2∇xgη1,2)− −ǫ1 2ρ2 1 Ro cos θe1∇x((f h2)2v2)−2ǫµ02∇x(f h2div xv2)+ +ǫρ2 1 Ro cos θe1h0 2w0 2|z=η1,2−δ−1γ0„βǫα0h0 1 6µ01 v1+ (v1−v2)«− −ǫρ2 1 Ro cos θe1f h2v2∇xgη1,2−ǫα02h0 2∇x(∆xη0) = 0. (2.130) Remark 2.2. We can use the equation (2.128) to write (2.127) as: ρ1∂t(f h1v1) + ρ1div x(f h1v1⊗v1)−ǫ 2ρ11 Ro cos θe1(f h1)2div xv1+ +ρ1 1 Ro sin θf h1(v1)⊥−2ǫµ01div x(f h1Dx(v1)) + 1 2ρ1 1 Fr2∇x(f h1)2− −ǫ 2ρ11 Ro cos θe1∇x((f h1)2v1) + ρ11 Fr2f h1∇xeb−2ǫµ01∇x(f h1div xv1)+ +ǫρ1 1 Ro cos θf h1ˆ(v1·∇xeb)e1−(v1·e1)∇xeb˜+ +δ−1γ0„βǫα0h0 1 6µ01 v1+ (v1−v2)«+βα0„v1+δ−1ǫγ0h0 1 6µ01 (v1−v2)«− −ǫα01f h1∇x∆xf h1−ǫα01f h1∇x∆xb+ρ2 1 Fr2f h1∇xf h2− −2ǫρ21 Ro cos θe1f h1∇x(f h2v2)−ǫα02f h1∇x∆xf h2= 0. (2.131) 2.4 Final models In this section we write the final equations for the two models obtained with dimension and dropping the cosines terms, having into account that 1 Rei =ǫµ0i, α =ǫα0, αi=ǫα0i, γ =ǫγ0.(2.132) We also divide the second and fourth equations in the systems by ρ1and ρ2respectively. For a good writing of equations we introduce some notation about the coefficients involved in the system. First we define the density relation by r=ρ2 ρ1, we explicit the dynamic viscosity as µi=ρiνifor i= 1,2, being νithe kinematic viscosity and finally we take the following definition for the friction and tension coefficients: γ= ˜γρ2,α= ˜αρ1,αi= ˜αiρiwith ˜γ, ˜αand ˜αibeing positive constants. Next we introduce some remarks about the approximations obtained mainly related to the friction terms. First, we state the system without viscosity, from equations (2.74), (2.86), (2.82) and 18 (2.85), we found the following system that we denoted by (BL1): (BL1) 8 > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > : ∂th2+ div x(h2v2) = 0; ∂t(h2v2) + div x(h2v2⊗v2) + 2Ω sin θh2(v2)⊥+ +1 2g∇xh2 2+gh2∇xη1,2= ˜γ(v1−v2); ∂th1+ div x(h1v1) = 0; ∂t(h1v1) + div x(h1v1⊗v1) + 2Ω sin θh1(v1)⊥+1 2g∇xh2 1+ +gh1∇xb+rgh1∇xh2=−r˜γ(v1−v2)−˜αv1. (2.133) In the same way, we consider equations (2.88), (2.131), (2.100) and (2.130) to write the viscous model with correction on the bottom and interface friction, (BL2): (BL2) 8 > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > : ∂th2+ div x(h2v2) = 0; ∂t(h2v2) + div x(h2v2⊗v2) + 2Ω sin θh2(v2)⊥+1 2g∇xh2 2+ +gh2∇xη1,2=δ−1˜γ„β˜αh1 6ν1 v1+ (v1−v2)«+ ˜α2h2∇x∆xh2+ +˜α2h2∇x∆xη1,2+ 2ν2div x(h2Dx(v2)) + 2ν2∇x(h2div xv2); ∂th1+ div x(h1v1) = 0; ∂t(h1v1) + div x(h1v1⊗v1) + 2Ω sin θh1(v1)⊥+1 2g∇xh2 1+ +gh1∇xb+rgh1∇xh2=−δ−1˜γ r „β˜αh1 6ν1 v1+ (v1−v2)«− −β˜α„v1+δ−1r˜γh1 6ν1 (v1−v2)«+ ˜α1h1∇x∆xh1+ +˜α1h1∇x∆xb+ 2ν1div x(h1Dx(v1)) + 2ν1∇x(h1div xv1), (2.134) being β=„1 + ˜α 3ν1 h1«−1 , δ = 1 + ˜γ 3„rh1 ν1 +h2 ν2«.(2.135) Remark 2.3. We want to remark that friction terms that we have obtained in (2.134) have the same order of the friction and viscosity coefficients. That is to say, if we suppose that the coefficients νi,˜αi,˜αand ˜γare of order ε, with ε∼10−3, then we next prove that the terms δ−1˜γr „β˜αh1 6ν1 v1+ (v1−v2)«, β˜α„v1+δ−1r˜γh1 6ν1 (v1−v2)«(2.136) have order εtoo. We write them as follows: δ−1˜γrβ ˜α 6ν1 h1v1+δ−1˜γr (v1−v2), β˜α v1+δ−1˜γrβ ˜α 6ν1 h1(v1−v2).(2.137) So, it is enough to prove that the coefficients given by δ−1˜γβ ˜α ν1 |{z } [1] , δ−1˜γ |{z} [2] ,and β˜α |{z} [3] (2.138) 19 have order ε. First we develop βand δ−1: β=3ν1 3ν1+ ˜αh1 , δ−1=3ν1ν2 3ν1ν2+ ˜γrν2h1+ ˜γν1h2 (2.139) and we observe that they have order ε0because νi,˜αand ˜γhave the same order. Now we study each term in (2.138). Term [1]: Since βand δ−1have order ε0and ν1,˜αand ˜γhave order εwe deduce that this first term have order ε. Term [2]: this term has order εbecause δ−1has order ε0and ˜γhas order ε. Term [3]: for the same reason that for the second one, we find the εorder for this term, because βhas order ε0and ˜αhas order ε. Remark 2.4. We have performed the deduction of a bilayer Shallow Water equations following the work developed in [[6]] for the one-layer case. Then, if we throw out the layer on top, we should get the model for one-layer obtained by them in that work. So, taking h2=ρ2= 0 in the final system (2.134): 8 > > > > > < > > > > > : ∂th1+div x(h1v1) = 0; ∂t(h1v1) + div x(h1v1⊗v1) + 2Ω sin θh1(v1)⊥+1 2g∇xh2 1+gh1∇xb+ +rgh1∇xh2=−β˜αv1+ 2ν1div x(h1Dx(v1)) + 2ν1∇x(h1div xv1); (2.140) that is just the same correction for the friction that Gerbeau and Perthame have found but in 2d case. Remark 2.5. In [[16]] a theoretical study of a simplified (BL2) model is performed proving the existence of global weak solution for the system above but in the particular case when b= 0. To obtain this result, the following form for the friction coefficient is taken: γ=h1h2 ν1 ν2h1+ν2 ν1h2 (2.141) Remark 2.6. The proof of the global weak solution for the (BL2) system is in course, it shall appear in a forthcoming paper. 3 Numerical assessment This section is devoted to check the validity of the new viscous bilayer model that we have derived in the previous section. In the first test we solve a 1D internal dam-break problem following the work [[6]] and we compare the numerical solution obtained by solving (BL1) and (BL2) with the Navier-Stokes equations. In Test 2 we make a comparison of the solution of the models for a 2D dam-break problem. The results obtained show us that the new viscous model improves the no viscous one for both unknowns, height and discharge. But as it is already confirmed in precedent works (see [[6, 9, 1, 12]]) we notice that the more significant difference relies on the discharge. Test 1: An internal dam-break problem. We present a test for which we compare the solutions obtained for Navier-Stokes equations with variable density with those given by systems (BL1) and (BL2). In the numerical discretization we have defined the friction term at the interface as follows: fric(v1, v2) = −κB(h1, h2)(v1−v2),with B(h1, h2) = h1h2 ρ1 ρ2h1+ρ2 ρ1h2 and κ > 0.(3.1) In [[6]] a dam-break test is calculated to validate the viscous model obtained. In this work we present a similar case but solving an internal dam-break problem where two flows with 20 different densities are involved, with the aim of emphasizing the importance of the friction term between layers. To make the comparison between the approximated Shallow-water systems and Navier- Stokes equations, we have computed the non-dimensional problem in each case without tension terms (i.e. ˜αi= 0, for i= 1,2). For the sake of clarity we specify these problems below. So for the first order approximation we have: (BL1adim) 8 > > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > > : ∂th2+ div x(h2v2) = 0; ∂t(h2v2) + div x(h2v2⊗v2) + 1 2 1 Fr2∇xh2 2+1 Fr2h2∇xη1,2= = ˜γ0(v1−v2); ∂th1+ div x(h1v1) = 0; ∂t(h1v1) + div x(h1v1⊗v1) + 1 2 1 Fr2∇xh2 1+ +1 Fr2h1∇xb+r1 Fr2h1∇xh2=−r˜γ0(v1−v2)−˜α0v1. (3.2) For the second order approximation we solve: (BL2adim) 8 > > > > > > > > > > > > < > > > > > > > > > > > > : ∂th2+ div x(h2v2) = 0; ∂t(h2v2) + div x(h2v2⊗v2) + 1 2 1 Fr2∇xh2 2+ +1 Fr2h2∇xη1,2=δ−1˜γ0„β˜α0h1 6ǫ2Re1v1+ (v1−v2)«+ +2 1 Re2 div x(h2Dx(v2)) + 2 1 Re2∇x(h2div xv2); (3.3) for the upper layer and for the lower one: (BL2adim) 8 > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > : ∂th1+ div x(h1v1) = 0; ∂t(h1v1) + div x(h1v1⊗v1) + 1 2 1 Fr2∇xh2 1+1 Fr2h1∇xb+ +r1 Fr2h1∇xh2=−δ−1˜γ0r„β˜α0h1 6ǫ2Re1v1+ (v1−v2)«− −β˜α0„v1+δ−1r˜γ0h1 6ǫ2Re1(v1−v2)«+ +2 1 Re1 div x(h1Dx(v1)) + 2 1 Re1∇x(h1div xv1), (3.4) being β=„1 + ˜α0 3ǫ2Re1h1«−1 , δ = 1 + ˜γ0 3ǫ2(rh1Re1+h2Re2).(3.5) 21 And finally for the Navier-Stokes equations we take the following problem: (NSadim) 8 > > > > > > > > > > > > > > < > > > > > > > > > > > > > > : ρ∂tv+ρv∇xv+ρw∂zv− −2 Rediv x(Dx(v)) −1 Re 1 ǫ2∂2 zv−1 Re∇x(∂zw) + ∇xp= =−γ[v]δε(ρ); ρ∂tw+ρv∇xw+ρw∂zw−1 Re∆xw− −1 ǫ2 1 Re∂z(div xv)−21 ǫ2 1 Re∂2 zw+1 ǫ2∂zp=−ρ1 ǫ2 1 Fr2−γ[w]δε(ρ); div xv+∂zw= 0. (3.6) with boundary conditions at the bottom: (1 Re(∇xw+1 ǫ2∂zv) = 1 ǫα v; w= 0.(3.7) And δεbeing an approximation of the Dirac mass on the interface, defined as follows: δε=1 εξ“ρ ε”|∇ρ|,with ξ(ϑ) = 1 2(1 + cos(πϑ)) if |ϑ|<1 0 otherwise (3.8) So the term involving δεin (3.6) is an approximation of the friction condition at the interface that we have in the two layers case. We take a domain Dwith length L= 1 meters and we set the friction coefficient as ν=ν0ζ, γ =γ0ζ, α =α0ζwith ζbeing a quantity depending on the jump height, given by ζ=h1L−h1R L. The height before the jump is taken as h1L= 0.9 and h1R= 0.1 after the jump. If we denote by ρi,i= 1,2 the densities associated with each layer, we have defined r=ρ2 ρ1= 0.98 and the upper layer density, ρ2= 1. We have also set the following constant data: ǫ= 0.04,1 Rei =1 Re =ǫ 10 (for i= 1,2), 1 Fr2= 1,and ζ= 10 to put in evidence the influence of the friction term. The resolution of problems (BL1adim) and (BL2adim) has been performed by using a Finite Volume method of Roe’s type; see [[11]]. We take the following initial conditions: h(t= 0) = h1Lx < 0; h1Rx > 0,q(t= 0) = 0.(3.9) In order to solve the Navier-Stokes equations for this problem, we have solved it in a two-dimensional domain D×Υ, Υ with length 1, in the terms that we specify next. To calculate the density ρ, we have solved the transport problem for the salinity Sgiven by: ∂tS+∇(uS)−ς∆S= 0,(3.10) ςbeing the molecular diffusion. Then, the density is updated by the state equation ρ=ρ0[1 + F(S)], where ρ0is a reference density and F(S) is a function of the salinity S. For the time being we are not able to solve the bilayer problem using Navier-Stokes equations so in order to simulate this situation, we have taken the following constant piecewise function for the initial density value, (see Fig. 3) related to initial condition (3.9): ρ(t= 0) = ρ2if {x < 0, y > h1L}and {x > 0, y > h1R} ρ1otherwise, (3.11) and the corresponding initial data for the salinity. At the initial time, we have taken u= 0. Regarding the constants involving the Navier-Stokes problem, we have fixed the reference 22 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1(BL1) (BL2) 1 1.002 1.004 1.006 1.008 1.01 1.012 1.014 1.016 1.018 1.02 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 1.002 1.004 1.006 1.008 1.01 1.012 1.014 1.016 1.018 1.02 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 1.002 1.004 1.006 1.008 1.01 1.012 1.014 1.016 1.018 1.02 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 1.002 1.004 1.006 1.008 1.01 1.012 1.014 1.016 1.018 1.02 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 1.002 1.004 1.006 1.008 1.01 1.012 1.014 1.016 1.018 1.02 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 1.002 1.004 1.006 1.008 1.01 1.012 1.014 1.016 1.018 1.02 Figure 2: Evolution in time of the interface from t=1 to t=6. 23 −5 0 5 −2 −1.8 −1.6 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 ρ2 ρ1 Figure 3: Initial density for Navier-Stokes problem. density ρ0= 1 and in our case, we have considered linear functions F,F(S) = bS for a positive constant b, concretely we set b= 1. Finally, the molecular diffusion for the salinity problem is taken as ς= 10−5. We numerically solve this problem using a Finite Element discretization in the stable pair of spaces (P2,P1). The computational work has been performed by using the software Freefem++ (http://www.freefem.org). Computing the test at time T= 6 seconds for the three problems, we compare the solutions obtained for the interface and the velocity. In Fig. 3 we show the interface level obtained in each case for times from t= 1 to t= 6. For problems (BL1adim) and (BL2adim) we plot the height of the lower layer, h1, but for Navier-Stokes equations, this value must be obtained as a function of the density profile, so we show the isolines of the density (colored lines). As we can see in the figure this approximation for Navier-Stokes is not too meaningful to compare it with the models studied here. Anyway we can see that the solution given by the two systems keep on the profile of the Navier-Stokes solution. Finally, in Fig. 3 we show the velocities of the lower layer obtained by solving problems (BL1adim), (BL2adim) and (NSadim) at times t= 1 to 6. We notice that the solutions of the Shallow-Water systems are getting further when the time increases and that the second order approximation gives us a closer solution to the Navier-Stokes velocity. Test 2: Circular dam-break problem in a 2D domain. We consider a circular dam-break problem in both, surface and interface with a no constant bottom. The domain is the square D= [0,2] ×[0,2], the bottom is given by the following function: b(x, y) = 1 8(1 + cos(2πx))(1 + cos(2πy)) (x−1)2+ (y−1)2≤0.12; 0 otherwise. (3.12) The initial condition is given by: h1(t= 0) + b(x, y) = 1.1 (x−0.9)2+ (y−1)2≤0.22; 0.6 otherwise. (3.13) h2(t= 0) + h1(t= 0) + b(x, y) = 1.7 (x−1.6)2+ (y−1)2≤0.12; 1.2 otherwise; (3.14) and q1(t= 0) = q2(t= 0) = 0. A longitudinal section in y= 1 of the height initial condition is shown in Fig. 5. The CFL is set to 0.7 and we consider a partition with ∆x= ∆y= 0.02, the final time is T= 2. The friction coefficients and the kinematic viscosity has been taken as ˜γ= ˜α=ν1= ν2= 10−3and the density ratio is set to 0.8. 24 −0.5 0 0.5 −2 0 2 4 6 8 10 12 14 16 x 10−3 Navier−Stokes (BL1) (BL2) −0.5 0 0.5 −2 0 2 4 6 8 10 12 x 10−3 −0.5 0 0.5 0 1 2 3 4 5 6 7 8 9 10 x 10−3 −0.5 0 0.5 −1 0 1 2 3 4 5 6 7 8x 10−3 −0.5 0 0.5 −1 0 1 2 3 4 5 6 7x 10−3 −0.5 0 0.5 −1 0 1 2 3 4 5 6x 10−3 Figure 4: Evolution in time of the velocity of the layer below for times t=1 to t=6. 25 4 Conclusions In this work we propose a bilayer Shallow Water 2D model, taking into account viscosity and tension effects on the surface and the interface. The model is obtained from the Navier-Stokes equations through a second order development in the asymptotic analysis and the integration process, following [[6]]. The main difficulty is related to the correction of the friction term at the interface. Usually to make this correction we must write the velocities at the interface in function of the average velocities. Due to that the friction term depends on the velocity difference of the two layers, we have a coupled problem and a second order friction correction cannot be performed layer by layer. To solve this problem we set out a linear system of equations where the unknowns are the velocities at the interface and at the bottom, that we write in function of the average velocities. So we obtain the correction for both friction terms: at the interface and at the bottom. And, in particular we can observe in the model the influence of the friction at the bottom in the upper layer. Finally we present two numerical tests to check the influence of the viscosity terms and friction corrections in the model. In the first test a one-dimensional internal dam-break problem is presented. We make a comparison between the solution of two models (first and second order) and the solution of Navier-Stokes equations with variable density. In this case we observe that the interface position and velocities computed by the models are comparable with the solution of the Navier-Stokes equations with variable density. Moreover, we can see that the velocity obtained from the second order model is closer to the velocity computed for Navier-Stokes problem. In the second test, a problem with higher velocities is considered. We set a double circular dam-break problem with a bump in the bottom. The test is designed in order to obtain a great difference between the velocities of the two layers, consequently we find an important influence of friction terms. As motivated in the first test, we show that the effects added in the second order model are significants. The research of E.D. Fern´andez-Nieto and G. Narbona-Reina to carry on this work was partially supported by the Spanish Government Research project MTM 2006-01275. The authors wish to thank Manuel J. Castro D´ıaz for interesting discussions about the numerical tests. References [1] E. Audusse. A multilayer Saint-Venant model: derivation and numerical validation. Discrete Contin. Dyn. Syst. Ser. B, 5(2):189–214, 2005. [2] D. Bresch and B. Desjardins. Existence of global weak solutions for a 2D viscous shallow water equations and convergence to the quasi-geostrophic model. Comm. Math. Phys., 238(1-2):211–223, 2003. [3] I. D. Chueshov, G. Raugel, and A. M. Rekalo. Interface boundary value problem for the Navier-Stokes equations in thin two-layer domains. J. Differential Equations, 208(2):449– 493, 2005. [4] S. Ferrari and F. Saleri. A new two-dimensional shallow water model including pressure effects and slow varying bottom topography. M2AN Math. Model. Numer. Anal., 38(2):211–234, 2004. [5] J. GarcRodrez. Paralelizaci esquemas de volmenes finitos: aplicacia resoluci sistemas de tipo aguas someras. Thesis of the University of Mga (Spain), 2005. [6] J. Gerbeau and B. Perthame. Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation. Discrete Contin. Dyn. Syst. Ser. B, 1(1):89–102, 2001. [7] M. Ishii and T. Hibiki. Thermo-fluid dynamics of two-phase flow. Springer, New York, 2006. With a foreword by Lefteri H. Tsoukalas. 32 [8] C. Lucas. Effet cosinus sur un mod`ele visqueux de type Saint-Venant et ses ´equations limites de type quasi-g´eostrophique et lacs. C. R. Math. Acad. Sci. Paris, 345(6):313–318, 2007. [9] F. Marche. Theoretical and numerical study of Shallow Water models. Application to Nearshore hydrodinamics. Thesis of the University of Bordeaux (France), 2005. [10] S. Ozaki, K. Hashiguchi, T. Okayasu, and D. H. Chen. Finite element analysis of particle assembly-water coupled frictional contact problem. CMES Comput. Model. Eng. Sci., 18(2):101–119, 2007. [11] C. Par´es and M. Castro. On the well-balance property of Roe’s method for nonconservative hyperbolic systems. Applications to shallow-water systems. M2AN Math. Model. Numer. Anal., 38(5):821–852, 2004. [12] M. Peybernes. Analyse de probl mathtiques de la mnique des fluides de type bi-couche et onti libre. Thesis of the University of Pascal Paoli (France), 2006. [13] Y. Song, D. McFarland, L. Bergman, and A. Vakakis. Stick-Slip-Slap Interface Response Simulation: Formulation and Application of a General Joint/Interface Element. CMES Comput. Model. Eng. Sci., 10(2):153–170, 2005. [14] A. Wacher and D. Givoli. Remeshing and refining with moving finite elements. Application to nonlinear wave problems. CMES Comput. Model. Eng. Sci., 15(3):147–164, 2006. [15] J. ZabsonrMods visqueux en smentation et stratification. Obtention formelle, Stabilitique et schs volumes finis bien ilibr . Thesis of the University of Savoie (France), 2008. [16] J. Zabsonrd G. Narbona-Reina. Existence of global weak solution for a 2D viscous bilayer Shallow-Water model. Nonlin. Anal. Real World Appl., 10(5):2971–2984, 2009. 33