Formal derivation of a bilayer model coupling shallow water and Reynolds lubrication equations: evolution of a thin pollutant layer over water

Formal derivation of a bilayer model coupling shallow water and Reynolds lubrication equations: evolution of a thin pollutant layer over water E.D. Fern´andez-Nieto∗, G. Narbona-Reina∗, J.D. Zabsonr´e† Abstract In this paper a bilayer model is derived to simulate the evolution of a thin film flow over water. This model is derived from the incompressible Navier-Stokes equations together with suitable boundary conditions including friction and capillary effects. The derivation is based on the different properties of the fluids, thus, we perform a multiscale analysis in space and time, and a different asymptotic analysis to derive a system coupling two different models: the Reynolds lubrication equation for the upper layer and the shallow water model for the lower one. We prove that the model is provided of a dissipative entropy inequality, up to a second order term. Moreover, we propose a correction of the model −by taking into account the second order extention for the pressure−that admits an exact dissipative entropy inequality. Two numerical tests are presented. In the first one we compare the numerical results with the viscous bilayer shallow water model proposed in [G. Narbona-Reina, J.D.D. Zabsonr´ e, E.D. Fern´ andez-Nieto, D. Bresch, CMES Comput. Model. Eng. Sci., 2009]. In the second test the objective is to show some of the characteristic situations that can be studied with the proposed model. We simulate a problem of pollutant dispersion near the coast. For this test the influence of the friction coefficient on the coastal area affected by the pollutant is studied. 1 Introduction In this paper, we derive a bilayer model of two immiscible fluids where the upper layer can be represented by a Reynolds lubrication model and the lower one by a shallow water model. This model can be used to simulate the evolution of a pollutant viscous fluid over water. Our purpose is to study the evolution of a system that consists of two layers of Newtonian viscous fluids with different properties. For the upper layer we consider a thin liquid ∗Dpto. Matem´atica Aplicada I, Universidad de Sevilla, Spain. [email protected], gnarb[email protected] †Institut des sciences exactes et appliqu´ees, Universit´e polytechnique de Bobo-Dioulasso, Burkina Faso. [email protected] 1 film of slow flow −Reynolds number of unity order−and for the lower one we assume a fluid like water −high Reynolds number−. Under these assumptions it is not possible to choose the same type of equations to model each layer. Thus, for the upper layer we follow the Reynolds’ theory of lubrication to deduce the equation that defines its behaviour. In this equation the pressure is the one which plays the main role. The Reynolds lubrication equation is classically used to model a fluid between two very close surfaces or a very thin film. This equation was first derived in [31] where several hypotheses are used in order to describe the pressure of the film. In fact, in this derivation, inertial, gravitational, viscosity variation, slip condition at the bottom, surface roughness and thermal effects are neglected. In addition, the author restricts his analysis to an isoviscous, incompressible fluid. Later, many works related to the derivation of generalized Reynolds lubrication equation have been made in order to get a better modelling of a thin film (see for example [11, 15, 35, 37]). In [28], a review of long-scale evolution of thin liquid films is presented. A general mathematical theory of Reynolds lubrication equations is introduced. In their analysis, the authors use the slip condition at the bottom and take into account surface tension effects at the free surface. A general nonlinear evolution equation or equations are then derived and several particular cases are considered. The condition on the free surface with surface tension effects revealed the presence of a term of the form ∂x(σ∂2 xh), being σthe surface tension coefficient and hthe height of the thin layer. This term appears at the leading order due to the scaling used for this purpose. In another way, the shallow water system, which is considered to model the lower layer, is applied to study a large number of geophysical and engineering applications as ocean circulation, coastal areas, rivers, etc. Many derivations of this model had been made in order to model shallow flows, (see for example [27], [42], [30]). The pioneering work [19] has been considered as a basis to develop the deduction of the shallow water layer of the model proposed in this paper. This derivation takes into account laminar friction at the bottom and viscous effects. Viscous and capillary effects are useful to obtain an existence result of global weak solutions in [41]. In [24], a viscous two dimensional one-layer shallow water system, taking into account surface tension, capillary effects and quadratic friction terms, has been derived. In [27], the authors derived a bilayer shallow water model where friction and surface tension at the interface and free surface are introduced with a second order approximation. In the model considered in this work, coupling Reynolds lubrication and shallow water equations, the shallow-flow assumption is taken, that is, the height (H) of the layer is much smaller than its length (L). Thus, for the derivation of the model all the variables are written in terms of this aspect ratio ε=H L, assumed small. As we have mentioned before, the obtention of the model is inspired on the simulation of the transport of a pollutant over water, so following this idea we consider that this smallness ratio εis not the same for the two layers. This means that the order of all the characteristic variables, namely, velocities, pressures, viscosities, space and time, are different for each layer. Due to these hypotheses, the idea of the present work is to make a multiscale analysis in space 2 and time for the incompressible bilayer Navier-Stokes equations and get a simplified system of three equations. In the literature one can find several works related to the multiscale analysis in time in order to model the evolution of the topography in oceans submitted to tidal effects, see for example [12, 20, 32, 34]. The novelty of this work rests on two features, that, to the best of our knowledge, have not been tackled before in the modelling of multilayer systems. The first one is the coupling of two different equations on a bilayer model and the second one is the multiscale analysis –in space and time– developed in the two layers. Both of them provide also the main difficulties to derive the proposed model. The paper is organized as follows: Section 2 is devoted to the derivation of the model. First we write the equations in non-dimensional variables, taking into account a different scale for each layer. Next, to deduce the shallow water equation we first perform the hydrostatic approximation and we use an asymptotic analysis to deduce the shallow water system. We also use an asymptotic analysis to deduce a Reynolds lubrication equation for the upper layer. In Section 3 we present the obtained model and in Subsection 3.1 we study the energy of the model. Finally, in Section 4 we describe a numerical scheme to discretize the proposed model and two numerical tests. 2 The thin film - shallow water model In this section we give a formal derivation of the model. To deduce the model we start from the incompressible Navier-Stokes equations for both layers and the steps to follow are: first, we nondimensionalize equations, then we take into account the asymptotic regime for the physical parameters and finally we integrate the equations to get the averaged model. Following the goal of this work we search for a bilayer model where each one of the fluids has different properties. So we consider that they have not only different physical properties, but also a different flow behaviour. As we have mentioned before, we would like to use this model to simulate the evolution of a pollutant over water, so we also want to consider different thickness for each layer. All these considerations imply some difficulties to develop the derivation of the model. First, due to the different flow behaviour, we must follow different ways to obtain the equations that model each layer. For the lower layer we derive a shallow water model and a Reynolds lubrication equation is used for the upper layer. The fact of considering different thicknesses leads us to take different characteristic space variables. Moreover, in order to take into account the “slow flow” property for the fluid on the upper layer we also consider different characteristic velocities. Thus, we must develop a two scale analysis of the problem in space and time. First, we introduce the domain of study and we set out the incompressible Navier-Stokes equations together with appropriate boundary conditions. In particular, the interaction 3 z h2 h1 b Γs Γ1,2 x Figure 1: Schematic representation of the two layers over the solid boundary b, with free surface denoted by Γs, interface between the two fluids, Γ1,2and heights labeled h1and h2. between the two fluids is held in the interface boundary conditions through the friction and capillary effects. We also take into account capillary effects at the free surface and friction force at the bottom. We consider the two-dimensional domain Ω(t) = Ω1(t)∪Ω2(t)∪Γb∪Γ1,2(t)∪Γs(t), where: Ω1(t) = {(x, z)∈R2:x∈ω, b(x)< z < I(x, t)}; Ω2(t) = {(x, z)∈R2:x∈ω, I(x, t)< z < η(x, t)}; Γ1,2(t) = {(x, z)∈R2:x∈ω, z =I(x, t)}; Γs(t) = {(x, z)∈R2:x∈ω, z =η(x, t)}; Γb={(x, z)∈R2:x∈ω, z =b(x)}; and ωis a bounded domain in R. We denote by Ithe interface level, I(x, t) = b(x)+h1(x, t), and by ηthe free surface, η(x, t) = b(x) + h1(x, t) + h2(x, t) (see Figure 1). We consider the incompressible Navier-Stokes equations for each layer:    ρi(∂tui+ui∂xui+wi∂zui) = −∂xpi+ρiνi(∂2 xui+∂2 zui); ρi(∂twi+ui∂xwi+wi∂zwi) = −∂zpi+ρiνi(∂2 xwi+∂2 zwi)−ρig; ∂xui+∂zwi= 0; i= 1,2.(1) where we denote by vi= (ui, wi) the velocity field for each layer, ρithe density, νiis the kinematic viscosity and pithe pressure, for i= 1,2; gis the constant gravity. As a general notation rule, the subscript 1 corresponds to the lower layer and the subscript 2 to the upper layer. We also introduce the ratio of densities r=ρ2 ρ1 4 and the stress tensors σi= 2ρiνiD(vi)−piId where D(vi) = ∇vi+∇vT i 2and Id denotes the identity matrix. To complete the problem, we impose the following boundary conditions: •At the free surface, z=η(x, t): –The normal stress balance: (σ2·ns) = −δ κns(2) where δis the surface tension coefficient, κ=−divxnsis the mean curvature and nsthe unitary outward normal vector to the free surface, ns=1 √1+|∂xη|2(−∂xη, 1). –The kinematic condition: ∂tη+u2∂xη=w2.(3) •At the interface, z=I(x, t): –The kinematic conditions for each velocity: ∂th1+u2∂xI=w2; ∂th1+u1∂xI=w1.(4) –The normal stress balance: (σ1·nI)n−(σ2·nI)n= (δIκInI)n,(5) with δIthe interfacial tension coefficient, κI=−divxnIis the mean curvature of the interface and nI=1 p1 + |∂xI|2(−∂xI,1), the unitary normal vector to the interface pointing from layer 1 to layer 2. The subscript ndenotes the normal component of the vector. –The friction condition (Navier-slip boundary condition): (σi·nI)τ=−cρ2(v1−v2)τfor i= 1,2; (6) The positive friction coefficient is denoted by c. The subscript τdenotes the tangential component of the vector. •At the bottom, z=b(x): We consider the Navier-slip boundary conditions: –The no penetration condition: v1·nb= 0,(7) where nb=1 p1 + |∂xb|2(−∂xb, 1). 5 –The friction condition: (σ1·nb)τ=α(u1)τ,(8) where αis the positive friction coefficient. Remark 2.1 A linear friction term between the two layers is considered here. It coincides with the friction law used in [27] that is of the form: c(v1−v2). In [9], this type of friction term is used to study a system of 3D Navier-Stokes equations in a two-layer thin domain with an interface condition of the form (νi∂zui j−k(u1 j−u2 j)) = 0, i, j = 1,2. 2.1 Dimensionless equations To derive the model it is suitable to write the system (1)-(8) in dimensionless form. First, let us divide the conservation equations by densities ρito get          ∂tui+ui∂xui+wi∂zui=−1 ρi ∂xpi+νi(∂2 xui+∂2 zui); ∂twi+ui∂xwi+wi∂zwi=−1 ρi ∂zpi+νi(∂2 xwi+∂2 zwi)−g; ∂xui+∂zwi= 0. (9) Now we set the dimensionless variables, where we must take into account the different properties of the fluids in two layers, so we make it separately. Nevertheless the characteristic variables of both layers must be related because we study the coupled system. We must indicate that the Reynolds layer is not only thin and slow, but thinner and slower than the shallow water layer. Thus, we will make an assumption relating the aspect ratios of the two layers. Nondimensionalization for layer 1: First, we establish the dimensionless variables for the lower layer following a shallow water nondimensionalization (see [19]). We denote by H,Land Uthe characteristic height, length and velocity respectively. Thus, the characteristic time is T=L U. To impose the shallow flow condition we assume that the aspect ratio between the characteristic height and length is small, as commonly we denote it by ε=H L. We denote with the superscript asterisk (∗) the dimensionless variables: x=Lx∗, z1=Hz∗ 1, h1=Hh∗ 1 u1=Uu∗ 1, w1=εUw∗ 1, t =Tt∗ 1 p1=ρ1U2p∗ 1, Fr1=U √gH ;Re1=UL ν1 ; where we denoted by Fr1the Froude number and by Re1the Reynolds number. Nondimensionalization for layer 2: For the upper layer we take a nondimensionalization suitable for a Reynolds lubrication 6 equation following [28]. It mainly affects to the characteristic pressure that are larger due to the lubrication effect. Now we denote by H2,U2and T2=L U2the characteristic height, velocity and time for this layer respectively. Thus, we also have a different ratio height-length, namely, ε2=H2 L. The dimensionless variables in this case read: x=Lx∗, z2=H2z∗ 2, h2=H2h∗ 2 u2=U2u∗ 2, w2=ε2U2w∗ 2, t =T2t∗ 2 p2=ρ2ν2U2 ε2H2 p∗ 2, Fr2=U2 √gH2 ;Re2=U2L ν2 . Since we study the coupled system we must take into account that the layer 2 is thinner and slower than the layer 1. These related aspects lead us to search for a relationship between the characteristic height and velocity of the two layers that express these properties. Thus, we assume that H2=εH;U2=ε2U, so consequently ε2=ε2and T2=1 ε2T. We write the dimensionless variables of the layer 2 in terms of H,Uand ε: x=Lx∗, z2=εHz∗ 2, h2=εHh∗ 2 u2=ε2Uu∗ 2, w2=ε4Uw∗ 2, t =1 ε2Tt∗ 2 p2=ρ2ν2U εH p∗ 2, Fr2=ε2U √gεH ;Re2=εUH ν2 . Remark 2.2 We would like to give some remarks about the dimensionless procedure. First, in order to justify our hypothesis on the aspect ratio ε2=ε2we focus on the oceanic circulation case. In accordance with [24], the average thickness Hof the oceans in coastal domain is nearly 100 m whereas their horizontal characteristic value is about 100 km. So the aspect ratio ε=H/L is about 10−3. Thus, we are in the case that the pollutant layer has an aspect ratio of 2∼10−6and a thickness 10−3times smaller than the layer 1, which is equivalent to a thickness of H2=εH ∼10cm, what seems reasonable. Secondly, we would like to clarify the existence of two dimensionless times t∗ 1and t∗ 2. It just comes from considering two characteristic velocities which automatically gives two different characteristic times Tand T2=1 ε2T. So the times do not have the same order, that’s why we also consider two dimensionless times. Regarding the friction coefficients we take α=Uα∗, c =Uc∗, 7 and the Capillary numbers −at the free surface (C) and the interface (CI)−as: C=ε2Uρ2ν2 δ, CI=Uρ1ν1 δI . We assume Has the characteristic height for the bottom, so b=Hb∗. Due to the different characteristic heights we must pay attention to the ranges for the dimensionless vertical component. Thus, for the layer 1, since z1∈[b, b +h1] we obtain directly that z∗ 1∈[b∗, b∗+h∗ 1]. For the layer 2, the rank of the dimensional variable is z2∈[b+h1, b +h1+h2] so we deduce that z∗ 2∈[1 ε(b∗+h∗ 1),1 ε(b∗+h∗ 1) + h∗ 2]. In what follows –for the second layer– we denote ηε=1 ε(b∗+h∗ 1) + h∗ 2the free surface level and Iε=1 ε(b∗+h∗ 1) the interface level. Next we write the equations and the boundary conditions in dimensionless form, we omit the superscript (∗) in the notation for the sake of clarity. •Layer 1: ∂t1u1+u1∂xu1+w1∂zu1−1 Re1∂2 xu1+1 ε2∂2 zu1+∂xp1= 0; (10) ∂t1w1+u1∂xw1+w1∂zw1−1 Re1∂2 xw1+1 ε2∂2 zw1+1 ε2∂zp1=−1 ε2 1 Fr2 1 ; (11) ∂xu1+∂zw1= 0.(12) •Layer 2: ε4Re2(∂t2u2+u2∂xu2+w2∂z2u2) = −∂xp2+∂2 z2u2+ε4∂2 xu2; ε8Re2(∂t2w2+u2∂xw2+w2∂z2w2) = −∂z2p2−ε4Re2 Fr2 2 +ε4(∂2 z2w2+ε4∂2 xw2); ∂xu2+∂z2w2= 0. (13) •Conditions at the free surface: ∂t1(b+h1) + ε3∂t2h2+ε2u2∂x(b+h1) + ε3u2∂xh2=ε3w2; (14) (∂z2u2+ε4∂xw2)1−ε2(∂x(b+h1+εh2))2−4ε3∂xu2∂x(b+h1+εh2) = 0; (15a) −p2+ 2ε4∂zw2−ε3∂x(b+h1+εh2)(∂z2u2+ε4∂xw2) =−ε5C−1∂2 x(b+h1+εh2) 1 + ε2(∂x(b+h1+εh2))23/2.(15b) 8 •Conditions at the interface: ∂t1(b+h1) + ε2u2∂x(b+h1) = ε3w2; (16) ∂t1(b+h1) + u1∂x(b+h1) = w1; (17) −2ρ1 1 Re1 (ε2∂xw1+∂zu1)∂x(b+h1) + ε2ρ121 Re1 ∂xu1−p1∂x(b+h1)2+ρ121 Re1 ∂zw1−p1 =−2ρ2 1 Re2 ε3ε4∂xw2+∂z2u2∂x(b+h1) + ε21 Re2 ρ2(2ε4∂xu2−p2)∂x(b+h1)2 +ρ2 1 Re2 (2ε4∂zw2−p2) + ερ1 C−1 I Re1 ∂2 x(b+h1)(1 + ε2∂x(b+h1)2); (18) −41 Re1 ∂xu1∂x(b+h1) + 1 Re1 (∂xw1+1 ε2∂zu1)(1 −ε2∂x(b+h1)2) =−r εcu1−ε2u2+ε2(w1−ε3w2)∂x(b+h1); (19) −4ε41 Re2 ∂xu2∂x(b+h1) + ε1 Re2 (ε4∂xw2+∂z2u2)(1 −ε2∂x(b+h1)2) =−1 εcu1−ε2u2+ε2(w1−ε3w2)∂x(b+h1). (20) •Conditions at the bottom: −u1∂xb+w1= 0; (21) −ρ1 4 Re1 ∂xu1∂xb+ρ1 1 Re1 (∂xw1+1 ε2∂zu1) −ε2ρ1 1 Re1 (∂xw1+1 ε2∂zu1) (∂xb)2=α εu1(1 + ε2(∂xb)2). (22) In the next subsections we develop the derivation of the equations for both layers. We begin with the shallow water layer following [19, 27] and secondly we deduce the Reynolds lubrication layer following [28]. In order to get the viscous effect and the influence of the pressure of the thin film flow on the system we derive a second order approximation, that is, we keep the terms of order ε0and εin the system. 2.2 Layer 1: shallow water flow To obtain the shallow water model we first take the hydrostatic approximation and then develop the asymptotic analysis of (10). 9 and the final equation for the layer 1 reads: ∂t1(h1¯u1) + ∂x(h1¯u2 1) + 1 2 1 Fr2 1 ∂xh2 1+1 Fr2 1 h1∂xb−4εµ01∂x(h1∂x¯u1) +r Re2 h1C−1 0∂3 x(b+h1+εh2) + β0∂x(εh2)−h1µ01C−1 I0∂3 x(b+h1) +r Re2 εh2C−1 0∂3 x(b+h1) + β0∂x(b+h1)+1 ρ1 α0γ(h1)¯u1= 0. (52) 3 Final model In this section we expose and discuss the final model obtained in the previous section as a formal second order approximation of the initial problem defined by equations (1)-(8). First we will write this system in dimensional variables, and secondly we study the energy balance of the model. The final deduced model is given in non-dimensional variables by equations (32), (52), (48) and (43). Notice that the model contains three equations, mass and momentum for the shallow water flow and the mass equation for the thin film flow. Next we present the system written in dimensional variables: (M1)≡                                      ∂th1+∂x(h1u1)=0, ∂t(h1u1) + ∂x(h1u12) + 1 2g∂xh2 1+gh1∂xb−4ν1∂x(h1∂xu1) +h1δ ρ1 ∂3 x(b+h1+h2) + rg∂xh2−h1 δI ρ1 ∂3 x(b+h1) +h2δ ρ1 ∂3 x(b+h1) + rg∂x(b+h1)+α ρ1 γ(h1)u1= 0. ∂th2+∂x(h2u1) + ∂x−h221 ρ21 c+1 3ν2 h2∂xp2= 0. (53) with ∂xp2=δ ∂3 x(b+h1+h2) + ρ2g∂x(b+h1+h2) and γ(h1) = 1 + α 3ν1 h1−1 .(54) Remark 3.1 The continuity of the tangential stress is usually considered for the interface separating two immiscible viscous fluids and it may correspond to the most realistic case for fluid-fluid interfaces. In our case −due to the different properties of the fluids involved in our system−we found it more appropriate to impose a Navier-slip condition. This 16 condition incorporates the possibility of fluid slip at the interface with the thin lubrication layer [4, 33, 40]. Nevertheless, we would like to know what is the influence of this change into the model. To introduce this condition into the model, we must replace the friction condition (6) by the following one (σ1·nI)τ=−(σ2·nI)τ.(55) Where by nIwe denote the unitary normal vector to the interface pointing from layer 1 to layer 2. But we also need an additional condition. When the continuity of the tangential stress is considered, the continuity of the velocities at the interface is also usually imposed [38, 36, 22], u1=u2at z=I(x, t).(56) By imposing this condition, we obtain the model (M1)without the term ∂x(−h2 21 ρ2 1 c∂xp2)in the third equation of (53). Note that formally it corresponds to impose an infinity friction coefficient, that implies the continuity of the velocity at the interface. 3.1 Energy of the model In this subsection we prove that the model (M1) admits a dissipative energy inequality up to a second order term. Then, we propose a variation of this model, named (M2). The objective is to obtain a model provided of an exact energy balance, as we expose later in Proposition 3.1. The idea of the modification is to change the value of the velocity of layer 1 at the interface written in terms of the pressure p2. In order to complete the equation for layer 1, the last step on the derivation has been to replace in the integrated momentum equation the values of u1|Iand u1|b, according to equations (50) and (51). As a variation of the model we propose to change the value of u1|Ito also consider the terms of order ε2coming from ∂xp2. Thus, from (43) we write: ∂xp2=C−1 0∂3 x(b+h1+εh2) + β0∂x(b+h1+εh2). Then, from (49) we deduce the value for the velocity u1|I: u1|I=εh2 c0Re2C−1 0∂3 x(b+h1+εh2) + β0∂x(b+h1+εh2). 17 The final system obtained in this case reads as follows: (M2)≡                                      ∂th1+∂x(h1u1)=0, ∂t(h1u1) + ∂x(h1u12) + 1 2g∂xh2 1+gh1∂xb−4ν1∂x(h1∂xu1) +h1δ ρ1 ∂3 x(b+h1+h2) + rg∂xh2−h1 δI ρ1 ∂3 x(b+h1) +h2δ ρ1 ∂3 x(b+h1+h2) + rg∂x(b+h1+h2)+α ρ1 γ(h1)u1= 0, ∂th2+∂x(h2u1) + ∂x−h221 ρ21 c+1 3ν2 h2∂xp2= 0. (57) With ∂xp2=δ ∂3 x(b+h1+h2) + ρ2g∂x(b+h1+h2) and γ(h1) = 1 + α 3ν1 h1−1 .(58) We have the following result: Proposition 3.1 The models (M1)and (M2)defined by (53) and (57) respectively, admit an entropy inequality: ∂tu2 1 2+gh1(b+h1 2) + rgh2b+h1+h2 2 −δ ρ1 ∂t(h1+h2)∂2 x(b+h1+h2) + 1 2(∂x(h1+h2))2 +δI ρ1 ∂th1∂2 x(b+h1) + 1 2(∂xh1)2 +∂xh1u1u2 1 2+g(h1+b)+rgh2u1(b+ 2h1+h2) −∂x4ν1h1∂xu2 1 2+h2 2 1 ρ21 c+1 3ν2 h2∂x1 2(δ∂2 x(b+h1+h2) + ρ2g(b+h1+h2))2 −δ ρ1 ∂x(h1+h2)u1∂2 x(b+h1+h2)−(h1+h2)∂t(∂x(h1+h2)) +δI ρ1 ∂xh1u1∂2 x(b+h1)−h1∂t(∂xh1) ≤R. (59) For model (M2)we have an exact dissipative entropy inequality, with R=R2=−4ν1h1(∂xu1)2−rg2h2 21 c+1 3ν2 h2∂x(b+h1+h2)2−α ρ1 γ(h1)u2 1. 18 For model (M1)we have an approximated dissipative entropy inequality, with R=R2+rgu1∂xh2 2+δ ρ1 u1h2∂3 x(h2),  Proof. First we multiply the momentum equation by u1and use the mass conservation equation of the first layer for simplifications. Then, we obtain ∂tu2 1 2+gh1b+h1 2+∂xh1u1u2 1 2+g(h1+b)−∂x4ν1h1∂xu2 1 2 +δ ρ1 h1u1∂3 x(b+h1+h2)−δI ρ1 h1u1∂3 x(b+h1) + δ ρ1 h2u1∂3 x(b+h1+ξh2) +δ ρ1 ∂th2∂2 x(b+h1+h2) + δ ρ1 ∂x(h2u1)∂2 x(b+h1+h2) +rgh1u1∂x(h2) + rgh2u1∂x(b+h1+ξ h2) = −4ν1h1(∂xu1)2−α ρ1 γ(h1)u2 1.(60) where the coefficient ξ= 0 for model (M1) and ξ= 1 for model (M2). Secondly, we multiply the equation for the thin film flow by δ∂2 x(b+h1+h2) + ρ2g(b+ h1+h2) to obtain ∂trgb+h2 2+rgh1∂th2+rg(b+h1+h2)∂x(h2u1) +δ ∂th2∂2 x(b+h1+h2) + δ ∂x(h2u1)∂2 x(b+h1+h2) −∂xh2 2 1 ρ21 c+1 3ν2 h2∂x1 2(δ∂2 x(b+h1+h2) + ρ2g(b+h1+h2))2 =−h2 2 1 ρ21 c+1 3ν2 h2δ∂3 x(b+h1+h2) + ρ2g∂x(b+h1+h2)2.(61) We use the mass conservation equation to write u1h1∂xh2+h1∂th2=∂x(rgu1h1h2)−∂x(h1u1)h2+h1∂th2=∂x(rgu1h1h2)+ ∂t(h1h2),(62) and to develop the following product affecting the terms with δI: ∂x(h1u1)∂2 xh1=−∂t(h1∂2 xh1) + h1∂t(∂2 xh1) = −∂t(h1∂2 xh1) + ∂x(h1∂t∂xh1)−1 2∂t((∂xh1)2). (63) We make similar calculations for the term with δ. Finally, by adding equations (60) and (61), and taking into account (62) and (63), we obtain the entropy inequality (59).  19 Remark 3.2 For model (M1)we obtain a dissipative entropy inequality, up to the term rg u1∂xh2 2+δ ρ1 u1h2∂3 x(h2). It is of order of ε2because we assumed h2to be of order of ε.  4 Numerical tests We present in this section two academic tests (subsections 4.2 and 4.3). In the first one, a comparison with the numerical results obtained by the viscous bilayer shallow water model proposed in [27] is presented. In the second test the objective is to show some of the characteristic situations that can be studied with the proposed model. Concretely, we simulate the problem of a pollutant dispersion near the coast. We study the influence of the friction coefficient in order to determine the coastal area affected by the pollutant. Previously, in Subsection 4.1 we present a numerical scheme to discretize the proposed model (53)-(54). 4.1 Numerical scheme Before describing the numerical scheme we observe that the proposed model (53)-(54) can be written under the following structure: ∂tW+∂xF(W) = S1(W)∂xb+∂x(D(W)∂xS2(W)) + Sδ,1(W) + Sδ,2(W) + SF(W) (64) where, if we denote q1=u1h1, then W=  h1 q1 h2 , F =          q1 q2 1 h1 +1 2gh2 1+rgh1h2+ξrg 2h2 2 h2 q1 h1          , S1(W) =   0 −g(h1+r h2) 0 , D(W) =     0 0 0 0 4ν1h10 0 0 gh2 2(1 c+h2 3ν2 )     , S2(W) =    0 q1 h1 b+h1+h2   , Sδ,1(W) =     0 −¯ δ r(h1+h2)∂3 x(b+h1)+(h1+ξ h2)∂3 x(h2)+¯ δIh1∂3 x(b+h1) 0     , 20 Sδ,2(W) =     0 0 ¯ δ∂xh2 2(1 c+h2 3ν2 )∂3 x(b+h1+h2)    , SF(W) =    0 −γ1(h1)q1 h1 0   , and γ1(h1) = 3 ¯α ν1 3ν1+ ¯αh1 ,¯α=α ρ1 ,¯ δ=δ ρ2 ,¯ δI=δI ρ1 . Where the coefficient ξ= 0 corresponds to model (M1) and ξ= 1 to model (M2). The Jacobian matrix of F(W) is A(W) =     0 1 0 −u2 1+g(h1+rh2) 2u1gr(h1+ξh2) −u1 h2 h1 h2 h1 u1     .(65) Note that when h1= 0 we set that u1=q1/h1= 0. Then, the third component of the flux, (h2q1/h1), is zero and consequently their derivatives are also zero. That is, in the case that h1= 0, the third row of A(W) is set to zero. The eigenvalues of A(W) are λ1=u1−C, λ2=u1, λ3=u1+C, with C=sgh1(1 −r) + rgh11 + h2 h12+rgh2 2 h1 (ξ−1). Remark 4.1 Let us remark that for the characteristic variables we have that h2=εh1. If we denote ε=h2/h1, then, C=pgh1(1 −r) + rgh1(1 + ε)2+rgh1ε2(ξ−1). Thus, when εtends to zero, we have that Ctends to √gh1. That is, the eigenvalues coincide with the ones of the transport matrix for the shallow water equations with a passive scalar transport equation.  For the discretization of the system we use a finite volume method. Computing cells Ii= [xi−1/2, xi+1/2] are considered. For simplicity, we suppose that these cells have constant size ∆x. Let us define xi+1 2=i∆xand xi= (i−1/2)∆x, the center of the cell Ii. Let ∆t be the time step and define tn+1 =tn+ ∆t, being t0= 0. Wn idenotes the approximation of the cell averages of the exact solution provided by the numerical scheme: Wn i∼ =1 ∆xZxi+1/2 xi−1/2 W(x, tn)dx. (66) 21 The source term S(W)∂xb(x) is discretized following the ideas introduced in [8] and [29]. The discretization of B(W)∂xWfirst requires us to interpret this term as a Borel measure (see [13]), depending on the choice of a family of paths linking two given states. Here the family of segments are considered as in [29]. We also consider a semi-discretization in time of the friction term between the fluid and the bottom. Let us suppose that the values Wn iare known. In order to advance in time we proceed as follows: •First Step. We define W∗ i= [h∗ 1,i q∗ 1,i h∗ 2,i]Tas W∗ i=Wn i−∆t ∆xDFn,+ i−1/2+DFn,− i+1/2,(67) where DFn,± i+1/2=DF± i+1/2(Wn i, Wn i+1) are the generalized Roe flux difference computed using the family of segments: DF± i+1/2(Wi, Wi+1) = 1 2F(Wi+1)−F(Wi)−S1(Wi+1/2)(bi+1 −bi)−Sδ,1,i+1/2∆x(68) ±Sδ,2,i+1/2+1 2|Ai+1/2|∆IW|i+1/2+D(Wi+1/2)S2(Wi+1)−S2(Wi) ∆x being |Ai+1/2|the absolute value of the Roe matrix Ai+1/2, Wi+1/2=Wi+Wi+1 2,∆IW|i+1/2=  bi+1 +h1i+1 −bi−h1i q1i+1 −q1i h2i+1 −h2i , Sδ,1,i+1/2=    0 −¯ δ r(h2,i+1/2+h1,i+1/2)φi+1/2(b+h1)+(h1,i+1/2+ξ h2,i+1/2)φi+1/2(h2) 0     +  0 ¯ δIh1,i+1/2φi+1/2(b+h1) 0 , Sδ,2,i+1/2=    0 0 ¯ δh2 2,i+1/2(1 c+h2,i+1/2 3ν2 )φi+1/2(b+h1+h2)     , and φi+1/2(ω) = 1 ∆x3(ωi+2 −3ωi+1 + 3ωi−ωi−1). •Second step. Semidiscretization in time of the friction term. We define 22 Wn+1 i= [h∗ 1,i qn+1 1,i h∗ 2,i]T, where qn+1 1,i =q∗ 1,i 1+∆tγ1(h∗ 1)/hn 1 . The scheme is L∞-stable, for small values of δ, under the CFL condition ∆t ∆xmax i|q1i| h1i +pgh1i+kD(Wi)kL∞/∆x≤1. 4.2 Test 1: comparison with a viscous bilayer shallow water model In this first test we compare the numerical results of the model proposed in this paper with the two-layer viscous shallow water model proposed in [27]. The objective of this test is to show the difference between the proposed model and the two-layer shallow water equations. To deduce this bilayer shallow water model two immiscible flows with different physical properties have been considered. It includes viscosity and friction effects on the bottom and at the interface level. It is obtained from an asymptotic analysis of non-dimensional and incompressible Navier-Stokes equations with hydrostatic approximation. In order to obtain the viscous effect into the model, a second order approximation is also considered. The model proposed in [27] reads as:                                                  ∂th1+∂x(h1u1) = 0, ∂t(h1u1) + ∂x(h1u12) + 1 2g∂xh2 1+gh1∂xb+rgh1∂xh2 =−γ−1c rβ¯αh1 6ν1 u1+ (u1−u2)+¯ δIh1∂3 x(b+h1) −β¯αu1+γ−1r¯αh1 6ν1 (u1−u2)+ 4ν1∂x(h1∂xu1), ∂th2+∂x(h2u2) = 0, ∂t(h2u2) + ∂x(h2u22) + 1 2g∂xh2 2+gh2∂x(b+h1) =γ−1cβ¯αh1 6ν1 u1+ (u1−u2)+ 4ν2∂x(h2∂xu2) + ¯ δh2∂3 x(b+h1+h2). (69) Where ¯α=α ρ1 , β =1 + ¯α 3ν1 h1−1 , γ = 1 + c 3rh1 ν1 +h2 ν2. Remember that cis the coefficient linked to the friction term between both layers (see equation (6)), αdefines the friction coefficient between the first layer and the bottom (see 23 equation (8)), ¯ δ=δ/ρ2and ¯ δI=δI/ρ1, with δand δIthe surface and interfacial tension coefficients respectively. In order to make a comparison of the numerical results of both models we consider a very simple test. The domain is set to [0,25] and it is discretized with 200 points. The final time is t= 2.5 s. A flat bottom is considered, b(x) = 0, and the initial conditions are (see Figure 2): q1(x, 0) = q2(x, 0) = 0, h1(x, 0) = 1, h2(x, 0) = 0.04 if x∈[12,13], 0 otherwise. 0 5 1 0 1 5 2 0 2 5 −0.2 0 0.2 0.4 0.6 0.8 1 b b+h1 b+h1+h2 q1 Figure 2: Test 1.- Initial conditions for Test 1 with flat solid surface, a local perturbation in h2and null fluxes for both fluids. No flux boundary conditions are considered, i.e. ∂xhi= 0, ∂x(hiui) = 0 for i= 1,2. From a numerical point of view they are imposed by considering a ghost cell where the value of the unknowns coincide with the value at the boundary cell. Let us remark that at the contact points between the pollutant drop and the liquid is not neccesary to use any special numerical treatment. We have try with several numerical tests and several numerical treatments and we conclude that the numerical results are not sensitive. Note that this is a consequence of the different nature of the two fluids, with different characteristic length and velocity. Moreover, the pressure terms are more significative that the surface tension terms in the framework of the model. Nevertheless, if we consider the interaction between thin film flows at the contact points (the three-phase gas/liquid/liquid points) a special treatment must be considered. For example, a balance of local forces at contact points is used to study the hydrodynamics of the spreading of one thin liquid film on the surface of another in [26] (see also [14], [39]). For this test the friction coefficients are set to c= 1 and ¯α= 10−3. We set ρ1= 1027, corresponding to the sea water, and ρ2= 920, δ= 0.033, corresponding to the density and 24 tension surface coefficient of a marine residual fuel. The interfacial tension coefficient is set to δI= 0.027, corresponding to an oil/sea-water interface (see [21]). In the areas where h2= 0 we set δI= 0.072, corresponding to the sea-water surface tension. Then, we have r= 0.8958,¯ δ=δ ρ2 = 0.3587 ×10−4,¯ δI=δI ρ1 =0.2629 ×10−4if h2>0, 0.7010 ×10−4if h2= 0.(70) The water viscosity and the pollutant viscosity are set to ν1= 10−6, ν2= 5.9783 ×10−4.(71) In Figure 3 the evolution of the layers for the proposed model, coupling Saint Venant and Reynolds lubrication equations (SVR in what follows) is compared with the viscous bilayer shallow water system (2SW in what follows). We can observe that the evolution of the upper layer is completely different for these two models. The layer on the top is supposed to be a thin film flow, so the velocity of this layer must be lower than the velocity of the water fluid considered on the shallow water layer. This hypothesis on the velocity has been taken into account to deduce the SVR model. This is the behavior that we observe in this test (see Figure 3, left column). In Figure 4, a zoom of the evolution of the free surface near the pollutant is presented. If we look at the 2SW model solutions (Figure 3, right column) we can see that the velocity of the pollutant layer is very similar to the velocity of the water layer. Let us remark that to derive the 2SW model, the same shallow water assumptions for both layers, with different density and viscosity, have been considered. With this numerical test we see that the 2SW is not well adapted to study the evolution of a thin pollutant layer over water. It is necessary to consider a Reynolds lubrication theory to model the evolution of the thin film flow. 4.3 Test 2: pollutant dispersion near the coast The aim of this test is to show some of the characteristic situations in which the proposed model can be used. Concretely, we study the dispersion of a pollutant near the coast, for an academic case. The domain is [0,12] and it is discretized with 200 points. The bottom is defined by the following function: b(x) =      e−(x−8)2 10 if x≤8, 1+e−(x−20)2 50 −e−122 50 if x > 8; As initial condition we set a horizontal free surface and two pollutant slicks (see Figure 5). Concretely, q1(x, 0) = q2(x, 0) = 0, h1(x, 0) = max(0.78 −b(x),0) −h2(x, 0), 25 Conclusions A new bilayer model is presented in this paper to simulate the transport of a viscous thin layer of fluid over water. For this aim, two kinds of equations has been considered for each layer: Reynolds lubrication equation to model the upper layer and a shallow water model to describe the evolution of the water layer. An analysis of two different scales in space and time is done in the derivation of the coupling model. The model can be applied to simulate a transport of a viscous pollutant over water. We have proved that the model verifies, up to a second order term, a dissipative entropy inequality. Moreover, we have proposed a correction of the model that takes into account the second order extention for the pressure law. This version of the model verifies an exact dissipative entropy inequality. Finally, some academic numerical tests are presented. The objective of the first test is to show the difference between the proposed model and the bilayer shallow water equations. We have observed that the bilayer shallow water model is not well adapted to study the evolution of a thin pollutant layer over water. It is necessary to consider a Reynolds lubrication theory to model the evolution of the thin film flow, as we consider in the proposed model. In the second test we simulate the problem of a pollutant dispersion near the coast. In this test we have studied the influence of the friction coefficient on the amplitude of the coastal area affected by the pollutant. Acknowledgments This research has been partially supported by the Spanish Government Research project MTM2009-07719. The authors would like to thank to the Instituto Universitario de Investigaci´on de Matem´aticas de la Universidad de Sevilla (IMUS) for the grant awarded to J.D.Z. References [1] C. Ancey. Plasticity and geophysical flows: A review. J. Non-Newtonian Fluid Mech. 142, 4-35, 2007. [2] M. Arghir, A. Alsayed, D. Nicolas. The finite volume solution of the Reynolds equation of lubrication with film discontinuities. Int. J. Mech. Sci. 44, 2119–2132, 2002. [3] E. Audusse. A multilayer Saint-Venant model. Disc. Cont. Dyn. System. Series B, 5(2), 189-214, 2005. [4] M.G. Blyth, C. Pozrikidis. Stagnation-point flow against a liquid film on a plane wall. Acta Mech. 180, 203-219, 2005. [5] G. Capriz. On the Vibrations of Shaft Rotating on Lubricated Bearings. An. Math. Pure. Appl. 50, 223–248, 1960. 32 [6] M.J. Castro, A.M. Ferreiro, J.A. Garc´ıa, J.M. Gonz´alez, J. Mac´ıas, C. Par´es, M.E. V´azquez. On the numerical treatment of wet/dry fronts in shallow flows: applications to one-layer and two-layer systems. Math. Comp. Model. 42 (3-4): 419–439, 2005. [7] M. J. Castro D´ıaz, E. D. Fern´andez-Nieto, A. M. Ferreiro. Sediment transport models in Shallow Water equations and numerical approach by high order finite volume methods. Computers & Fluids 37, 299–316, 2008. [8] T. Chac´on, A. Dom´ınguez, E.D. Fern´andez. Asymptotically balanced schemes for non-homoegeneous hyperbolic systems – application to the Shallow Water equations. C.R. Acad. Sci. Paris, Ser. I 338: 85–90, 2004. [9] I.D. Chueskov , G. Raugel and A.M. Rekalo. Interface boundary value problem for the Navier-Stokes equations in thin domains. J. Differential Equations 208, 449-493, 2005. [10] G. Cimatti. A Rigorous Justification of the Reynolds Equation. Quart. Appl. Math. XLV (4), 627–644, 1987. [11] W. F. Cope. The hydrodynamic theory of film lubrication. Proc. R. Soc. London, Ser. A 197, 201–217, 1949. [12] S. Cordier, C. Lucas, J. D. D. Zabsonr´e. A two time-scale model for tidal bed-load transport. Commun. Math. Sci. (in press) [13] G. Dal Maso, P.G. LeFloch, F. Murat. Definition and weak stability of nonconservative products. J. Math. Pures Appl. 74: 483–548, 1995. [14] N.D. DiPietro, C. Huh, R.G. Cox. The hydrodynamics of the spreading of one liquid on the surface of another. J. Fluid Mech. 84, 529–549, 1978. [15] D. Dowson. A generalized Reynolds equation for fluid-film lubrication. Int. J. Mech. Sci. 4, 159–170, 1962. [16] H. G. Elrod. A Derivation of the Basic Equations for Hydrodynamics Lubrication with a Fluid Having Constant Properties. Quart. Appl. Math. 27, 349–385, 1960. [17] S. Ferrari, F. Saleri. A new to dimensional Shallow-Water model including pressure effects and slow varying bottom topography. M2AN 38(2), 211–234, 2004. [18] J.A. Garc´ıa-Rodr´ıguez. Paralelizaci´on de esquemas de volumenes finitos: aplicaci´on a la resoluci´on de sistemas de tipo aguas someras.Ph.D. thesis, University of M´alaga, 2005. [19] F. Gerbeau, B. Perthame. Derivation of viscous Saint-Venant system for laminar Shallow-Water; Numerical validation. Disc. Cont. Dynam. Syst. Series B 1(1), 89– 102, 2001. 33 [20] S. Hulscher. Formation and migration of large-scale, rhythmic sea-bed patterns: a stability approach, Ph.D. thesis, Utrecht University, 1996. [21] A. Konno, K. Izumiyama. On the relationship of the oil/water interfacial tension and the spread of oil slick under ice cover. Proceedings of the 17th International Symposium on Okhotsk Sea & Sea Ice, 275–282, 2002. [22] J.N. Koster. Multilayer fluid dynamics of inmiscible fluids. NASA. Lewis Research Center, Second Microgravity Fluid Physics Conference, 65–71, 1994. [23] T. Koyama, I. Neretnieks, L. Jing. A numerical study on differences in using Navier- Stokes and Reynolds equations for modeling the fluid flow and particle transport in single rock fractures with shear. Int. J. Rock Mech. & Min. Sci. 45, 1082–1101, 2008. [24] F. Marche. Theoretical and numerical study of Shallow Water models. Applications to Nearshore hydrodynamics.Ph. D. Thesis, University of Bordeaux, 2005. [25] E. Marusic-Paloka, M. Starcevic. Derivation of reynolds equation for gas lubrication via asymptotic analysis of the compressible navier-stokes system. Nonlinear Analysis: Real World Applications, 4565–4571, 2009. [26] M.J. Miksis, J.M. Vanden-Broeck. Motion of a triple junction. J. Fluid Mech. 437, 385–394, 2001. [27] G. Narbona-Reina, J. D. D. Zabsonr´e, E.D. Fern´andez-Nieto, D. Bresch. Derivation of a bi-layer Shallow-Water model with viscosity. Numerical validation. CMES 43(1), 27-71, 2009. [28] A. Oron, S. H. Davis, S. G. Bankoff. Long scale evolution of thin films. Rev. Mod. Phys. 69(3), 931–980, 1997. [29] C. Par´es, M.J. Castro. On the well-balance property of Roe’s method for nonconservative hyperbolic systems. Applications to shallow-water systems. ESAIM-Math. Model. Num. 38 (5), 821–852, 2004. [30] M. Peybernes. Analyse de probl`eme math´ematiques de la m´ecanique des fluides de type bi-couche et `a fronti`ere libre.Ph.D. thesis, University of Pascal Paoli, 2006. [31] O. Reynolds. On the theory of lubrication and its application to Mr. Beauchamp Tower’s experiment. Phil. Trans. R. Soc. London Part I, 228–310, 1886. [32] P. C. Roos. Seabed Pattern Dynamics and Offshore Sand Extraction, Ph.D. thesis, University of Twente, 2004. [33] B. Santra, B. S. Dandapat and H. I. Andersson. Axisymmetric stagnation-point flow over a lubricated surface. Acta Mechanica, 194(1-4), 1–10, 2007. 34 [34] H. Schuttelaars. Evolution and Stability Analysis of Bottom Patterns in Tidal Embayements, Ph.D. thesis, Utrecht University, 1997. [35] J. B. Shukla, S. Kumar, P. Chandra. Generalized Reynolds equation with slip at bearing surfaces: Multiple-layer lubrication theory . Wear 60, 253–268, 1980. [36] H.A. Stone, J.W.M. Bush. Time-dependent drop deformation in a rotation high viscosity fluid. Quart. App. Math. 54(3), 551–556, 1996. [37] G.H. Wannier. A Contribution to the Hydrodynamics of Lubrication. Quart. Appl. Math. 8, 1–32, 1950. [38] J. Wilkes. Fluid Mechanics for Chemical Engineers with Microfluidics and CFD. Prentice Hall, 2006. [39] S.D.R. Wilson, J. Williams. The flow of a liquid film on the inside of a rotating cylinder, and some related problems. Phys. Fluids 9, 2184–2190, 1997. [40] A. Yeckel, L. Strong, S. Middleman. Viscous film flow in the stagnation region of the jet impinging on planar surface. AIChE J. 40, 1611–1617, 1994. [41] J.D.D. Zabsonr´e, G. Narbona-Reina. Existence of a global weak solution for a 2D viscous bi-layer Shallow Water model. Nonlinear Analysis: Real World Applications 10, 2971–2984, 2009. [42] J.D.D. Zabsonr´e. Mod`eles visqueux en s´edimentation et stratification: obtention formelle, stabilit´e th´eorique et sch´emas volumes finis bien ´equilibr´es. Ph.D. thesis, University of Savoie, 2008. 35