Global stability of weak solutions for a multilayer Saint–Venant model with interactions between layers

Nonlinear Analysis 163 (2017) 177–200

Contents lists available at ScienceDirect

Nonlinear Analysis www.elsevier.com/locate/na

Global stability of weak solutions for a multilayer Saint–Venant model with interactions between layers Bernard Di Martinoa,d , Boris Haspotb,d , Yohan Penelc,d, * a

Université de Corse, SPE, UMR CNRS 6134, Campus Grimaldi, BP 52, 20250 Corte, France Université Paris Dauphine, PSL Research University, Ceremade, UMR CNRS 7534, Place du Maréchal De Lattre De Tassigny, 75775 Paris cedex 16, France c Sorbonne Universités, UMPC, Univ. Paris 06, LJLL, UMR CNRS 7958, 75005 Paris, France d ANGE project-team (Inria, CEREMA, UPMC, CNRS), 2 rue Simone Iff, CS 42112, 75589 Paris cedex 12, France b

article

info

Article history: Received 5 December 2016 Accepted 26 July 2017 Communicated by Enzo Mitidieri MSC 2010: 35B65 35Q35 76D03

abstract This paper concerns the global stability of weak solutions for the multilayer system introduced by Audusse et al., which models incompressible free surface flows. To do this, it is proven that this model admits the so called BD-entropy and a gain of integrability on the velocity. It allows to obtain enough compactness estimates in order to show the stability of global weak solutions. © 2017 Elsevier Ltd. All rights reserved.

Keywords: Navier–Stokes equations Free surface flow Degenerate viscosities Global weak solutions Saint–Venant model

1. Introduction The issue of modelling and simulating free-surface flows is extensively addressed in the literature. It is of major interest for a large amount of engineering applications such as the design of harbours, the protection of coasts, the production of energy or the prevention of natural hazards. Depending on the wavelengths of hydrodynamic processes at stake, several models of reduced complexity have been designed. A renowned simplified model implemented in many industrial codes is the system of viscous Shallow Water (SW) equations [16,24] which consists of a first-order hyperbolic partial differential equation (PDE) modelling the conservation of volume and of a second-order parabolic PDE for the momentum. The SW

*

Corresponding author at: Sorbonne Universit´ es, UMPC, Univ. Paris 06, LJLL, UMR CNRS 7958, 75005 Paris, France E-mail addresses: [email protected] (B. Di Martino), [email protected] (B. Haspot), [email protected] (Y. Penel). http://dx.doi.org/10.1016/j.na.2017.07.010 0362-546X/© 2017 Elsevier Ltd. All rights reserved.

178

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

Fig. 1. Multilayer approach.

equations are dedicated to a specific regime of water flows, namely when dispersion effects can be neglected and for water heights small compared to the characteristic longitudinal length of the domain. For such flows, the SW equations turn out to provide reliable numerical results. From SW equations to the Navier–Stokes (NS) equations with a free surface, there exists in the literature a hierarchy of models of increasing complexity including Boussinesq type models [17,23,29–31] with higher order derivatives to account for dispersion effects (necessary for modelling shoaling) or non-hydrostatic models [9,11] with a larger amount of unknowns (like the hydrodynamic pressure). The derivation of intermediate models is aimed at widening the range of applications of hydrodynamic models. For the specific regime addressed by the SW equations, another technique consists in splitting the flow into horizontal layers similar to a discretisation procedure along the vertical axis in order to improve the accuracy of the results. In this framework, the SW equations correspond to a coarse vertical mesh with a single layer. As a consequence, this multilayer approach is still relevant for non shallow flows. First, such models have been introduced with 2 or 3 layers for immiscible multifluid flows [10,27,28]. They were then extended to an arbitrary number of layers without [2,4] or with [3,14,15,22] mass transfer between layers. A major consequence is the noticeable increase of the number of unknowns related to the number of layers. In the inviscid case, open questions like the hyperbolicity of the model still hold (let us mention that recently Aguillon et al. in [1] proved the well-posedness of the Riemann problem for a two layer model). In the present work, we focus on the viscous case. We are interested in proving the stability of global weak solutions (in a future work, we shall consider the construction of global approximate solution which will imply the existence of global weak solutions). In the sequel of this section, the equations under study are detailed (Section 1.1) and a review of classical techniques to obtain global existence of weak solutions is presented (Section 1.2). The main result is stated in Section 1.3 (Theorem 1.1). Elements necessary for its proof are given in subsequent sections. 1.1. The multilayer Saint–Venant model We consider in this paper a multilayer description of a geophysical flow with a free surface and a varying topography issue from [3] that we are going to recall. N is the number of layers which might correspond to physical discontinuities but in the present approach layers are predetermined elements of the discretisation. Horizontal layers ℓα are separated by given surfaces z = zα+1/2 (t, x) where x ∈ Td , with d ∈ {1, 2, 3}.1 See Fig. 1 for notations. Without loss of generality, we assume that all layers have the same thickness 1 The model considered here results from an averaging process over the vertical axis applied to the Navier–Stokes equations. Hence it has a smaller dimension: d = 1 corresponds to the 2D Navier–Stokes equations and d = 2 to the 3D equations.

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

179

hα = zα+ 1 − zα− 1 = h/N , so that we have 2

2

α h(t, x). N The multilayer approach amounts to approximating the velocity field by a layer-wise constant function through a Galerkin discretisation procedure. More precisely, uα denotes an approximation of the average velocity over the layer ℓα ∫ z 1 (t,x) α+ 2 N uα (t, x) ≈ u(t, x, z) dz h(t, x) z 1 (t,x) zα+ 1 (t, x) = zb (x) + 2

α− 2

where u satisfies the Navier–Stokes equations. Let us introduce the notations u = (u1 , . . . , uN ) ∈ RN ,

u ¯=

and

N 1 ∑ uα . N α=1

The multilayer Saint–Venant model proposed by Audusse et al. [3] is obtained by integrating the hydrostatic Navier–Stokes equations over each layer, which reads with α ∈ {1, . . . , N }: ⎧ ∂t h + div(h¯ u) = 0, ⎪ ⎪ ) ( ⎪ g ⎨ ∂t (huα ) + div(huα ⊗ uα ) + ∇h2 = −gh∇zb + N uα+ 1 Gα+ 1 − uα− 1 Gα− 1 2 2 2 2 (1.1) 2 ( ) ⎪ +div 4νhD(uα ) + κ(uα+1 − uα ) − κ(uα − uα−1 ), ⎪ ⎪ ( ) ⎩ h(0, ·), u1 (0, ·), . . . , uN (0, ·) = (h0 , u1,0 , . . . , uN,0 ), with u0 = u1 and uN +1 = uN and g the acceleration of gravity. The model has for initial data (h0 , uα,0 ) with α ∈ {1, . . . , N }. The term zb (x) denotes the bottom topography assumed sufficiently smooth and stationary. In the sequel we assume that zb ∈ W 1,∞ (Td ). The term D(u) = 12 (∇u + ∇T u) corresponds to the linearised stress tensor. Gα+ 1 represents the mass 2 flux through the interface z = zα+ 1 from ℓα+1 to ℓα . It is defined by 2

Gα+ 1 2

α N ( ) 1 ∑ ∑ = 2 div h(uj − ui ) . N j=1 i=α+1

(1.2)

Let us consider standard kinematic boundary conditions at the surface and at the bottom which read G 1 = GN + 1 = 0. 2

2

Finally, we define the value of the interface velocity uα+ 1 by means of an upwind formula 2 { uα , if Gα+ 1 ≤ 0, 2 uα+ 1 = 2 uα+1 , otherwise. Let us observe that using the definition (1.3), we have: 1 1 Gα+ 1 uα+ 1 = Gα+ 1 (uα + uα+1 ) − |Gα+ 1 |(uα − uα+1 ). 2 2 2 2 2 2

(1.3)

(1.4)

Remark 1. In [14], the authors consider the following choice: 1 uα+ 1 = (uα + uα+1 ). 2 2 From (1.4), we notice that our choice for uα+ 1 gives the same value for Gα+ 1 uα+ 1 as in [14] up to an 2 2 2 additional term − 21 |Gα+ 1 |(uα − uα+1 ). 2

Remark 2. Let us notice that the definition of Gα+1/2 (1.2) implies the following partial mass law ∂t h + div (huα ) = N (Gα+ 1 − Gα− 1 ). 2

2

(1.5)

180

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

1.2. Main results of existence of global solutions for the compressible Navier–Stokes equations A large amount of papers in the literature has been devoted to the study of existence of global weak solutions for the compressible Navier–Stokes system. The first result of existence is due to P.-L. Lions in [21] where the mono-layer counterpart of (1.1) with constant viscosity coefficients is considered for a gamma law P (ρ) = ργ (with P the pressure and ρ the density which replaces the height h) with γ ≥ 59 in dimension 3 and γ ≥ 32 in dimension 2 (we refer also to [13] and [8]). The case of the compressible Navier–Stokes equations with degenerate viscosity coefficients is completely different in terms of analysis. Indeed one of the main issues is due to the fact that there is no equivalent to the so-called “effective pressure” (or in other words we cannot invert the viscous stress tensor). Recently several authors obtained significant progress on the existence of global weak solutions with degenerate viscosity coefficients. Bresch and Desjardins [6] introduced a new entropy (the so-called BD entropy) which gives new estimates on the gradient of the density provided that the viscosity coefficients µ and λ satisfy the following algebraic relation λ(ρ) = 2µ′ (ρ) − 2µ(ρ).

(1.6)

Let us mention that a particular choice of viscosity coefficients λ(ρ) = 0 and µ(ρ) = µρ satisfying (1.6) leads to the so-called viscous shallow water system which corresponds to our problem in the mono layer ′ framework. At least heuristically Bresch and Desjardins observed that the quantity µ√(ρ) ρ ∇ρ is conserved ∞ 2 d in L (0, T ; L (R )) norm for any T > 0. This allows to prove the existence of global weak solutions [6] with either a drag friction or a cold pressure term (a pressure that is singular at the vacuum). The addition of a friction term allows to get a gain of integrability on the velocity which provides enough compactness estimates in order to deal with the stability of the term ρu ⊗ u. Indeed compared with the constant viscosity case there is no control on the gradient of the velocity ∇u in L2 (R+ , L2 (Rd )) and it is not possible to apply classical Sobolev embeddings to deal with the term ρu ⊗ u. This is related to the fact that the viscosity coefficients are degenerate (see the relation (1.6)). The same remark holds when a cold pressure is added. We refer also to [7,35] for more developments on the existence of global weak solutions with a cold pressure or with a drag friction. The problem of stability of global weak solutions for the classical γ law (when 1 < γ < +∞ for d = 2 and 1 < γ < 3 for d = 3) has been solved by Mellet and Vasseur [26]. To do this they introduced a new energy estimate allowing a gain of integrability on the velocity. However the problem of existence of global weak solutions remains open. Indeed it remains to prove the existence of global approximate solutions to the system satisfying uniformly energy estimates, BD entropy and the gain of integrability ` a la Mellet–Vasseur which is tricky. However recently for the particular case of the shallow water system, the proof has been completed simultaneously and independently by Vasseur and Yu [32,33] and Li and Xin [20] using different methods. Concerning the existence of global strong solutions with large initial data for degenerate viscosity coefficients, the problem remains completely open in dimensions greater than 1. We can however mention some results in the case d = 1. For viscosity coefficients of the form µ(ρ) = ρα with 0 < α < 21 , the BD entropy allows to bound the density from below. It allowed Mellet and Vasseur [25] to prove the existence of global strong solutions for initial density far away from vacuum. Indeed the BD entropy gives a bound 1 on ∂x (ρα− 2 ) in L∞ (0, T ; L2 (R)) for all T > 0 and a control on ρ−1 in L∞ (0, T ; L∞ (R)) from Sobolev embeddings. Next it is classical to propagate any regularity on the density and the velocity in order to prove the uniqueness. This result has been recently extended by the second author in [18] to the case of general degenerate viscosity coefficients α ≥ 12 and in particular the shallow water system (α = 1) which corresponds to System (1.1) for d = 1. The main idea was to rewrite the system by introducing a suitable effective velocity v and apply a maximum principle.

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

181

Let us also recall some results on multilayer systems. To our knowledge most existence results concern immiscible fluid flows. In other words it means that there is no mass flux between each layer at the interface, in particular Gα+ 1 = 0 for any α. In [12,27] the authors obtained existence results of weak solutions for the 2 bilayer case with a viscous term of the form ν∆uα . In [34], it is proven stability of global weak solutions for viscous terms like in (1.1) with surface tensions and with test functions depending on the density itself. When mass transfer is involved, let us mention the work from Fern´andez Nieto et al. [14] where the authors construct numerical solutions of finite element type satisfying the classical energy inequality. In our paper, we prove the stability of global weak solutions of System (1.1). The main difficulty comes from the terms describing the transfer of flux between the layers which are not taken into account in the immiscible case. In particular it makes the analysis more difficult when we wish to prove the BD entropy and the gain of integrability ` a la Mellet–Vasseur which gives enough compactness informations to deal with the term huα ⊗ uα . These two estimates are the cornerstone of the proof of stability of global weak solutions following the arguments developed in [26]. However the lack of compactness for the mass flux terms prevents from recovering the expected limit. This is due to the fact that we cannot prove the convergence almost everywhere of the terms Gnα+ 1 . 2 In a future work, we shall prove the existence of global weak solutions. It remains essentially to construct global approximate solutions satisfying uniformly all the entropy inequalities (in order to do this, we shall follow the method developed in [32,33]. 1.3. Main results Before stating the main result (Theorem 1.1), we define the notion of weak solutions in the following way. Definition 1. Let d ∈ {1, 2, 3} be the space dimension. The solution (h, u1 , . . . , uN ) is said to be a global weak solution of (1.1) supplemented with initial conditions h(0, ·) = h0 , (huα )(0, ·) = mα,0 ,

(1.7)

√ √ h0 ∈ L1 (Td ), h0 ∇ h0 ∈ L2 (Td ), h0 ≥ 0, √ h0 |uα,0 | ∈ L2 (Td ), √ √ 2 h0 |uα,0 | log(1 + |uα,0 | ) ∈ L2 (Td ),

(1.8)

such that for any α ∈ {1, . . . , N }:

if the following smoothness assumptions are satisfied for any α ∈ {1, . . . , N }: ( ) √ ( ) √ ( ) • h ∈ L∞ 0, T ; L1 (Td ) , ∇ h ∈ L∞ 0, T ; L2 (Td ) , huα ∈ L∞ 0, T ; L2 (Td ) , √ √ ( ) √ ( ) 2 • h∇uα ∈ L2 (0, T ) × Td , h|uα | log(1 + |uα | ) ∈ L∞ 0, T ; L2 (Td ) , with h ≥ 0 satisfying in the sense of distributions over [0, T ] × Td for any α ∈ {1, . . . , N }: {

∂t h + div(huα ) = N (Gα+ 1 − Gα− 1 ), 2 2 h(0, ·) = h0 ,

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

182

and if the following equality holds for all smooth test functions φ(t, x) with compact support such that φ(T, ·) = 0, we have: ∫ ∫ T ∫ [ √ (√ ] ) √ √ g mα,0 · φ(0, ·) dx + h huα ∂t φ + huα ⊗ huα : ∇φ + h2 divφ dxdt 2 Td 0 Td ∫ T ∫ [ (G 1 M 1 √ α+ 2 α+ 2 + N (uα + uα+1 ) − h(uα − uα+1 ) 2 2 d 0 T ) Gα− 1 Mα− 1 √ 2 2 − (uα−1 + uα ) + h(uα−1 − uα ) 2 2 ] − gh∇zb + κ(uα+1 − 2uα + uα−1 ) · φ dxdt − ⟨4νh D(uα ), ∇φ⟩ = 0, Gn d

2

α+ 1

where Mα+ 1 is the weak limit in L ((0, T ) × T ) of hn 2 . Moreover, we give sense to the diffusion term 2 and the flux term: ∑∫ T ∫ √ √ ) (√ huα,i h∂jj φi + 2∂j φi ∂j h dxdt ⟨4νh D(uα ), ∇φ⟩ = − i,j

0

Td

and Gα+ 1 defined by (1.2). 2

) ( Remark 3. In the previous definitions, the sequences (hn )n∈N and Gnα+ 1 n

(h

, unα )n∈N

2

of global weak solutions defined in Theorem 1.1.

are related to the sequence n∈N

Let us now state our main result about global weak solutions for the multilayer system (1.1). Theorem 1.1. Let 1 ≤ d ≤ 3 and (h0 , m1,0 , . . . , mN,0 ) satisfying the assumption (1.8). Let us assume that there exists a sequence of global weak solutions (hn , un1 , . . . , unN )n∈N for System (1.1) such that the energy inequalities (2.1), (2.3) and (2.4) are uniformly verified. In particular the corresponding initial data are chosen such that: hn0 > 0, hn0 −−−−−→ h0 in L1 (Td ), hn0 un0,α −−−−−→ h0 u0,α n→+∞

n→+∞

and satisfy the following bounds (where C > 0 is independent from n): ) ∫ (∑ ∫ ⏐ √ ⏐2 N n 2 ⏐ ⏐ n |u0 | n 2 + (h0 ) dx < C, h0 ⏐∇ hn0 ⏐ dx < C, 2 Td Td α=1 and ∫

N ∑

Td α=1

2

hn0

( ) 1 + |un0 | 2 log 1 + |un0,α | dx < C. 2

In addition we assume that hn is a continuous function on R+ × Td such that for any (t, x) ∈ R+ × Td , we have: hn (t, x) > 0. √ √ √ Then√up to a subsequence, (hn , hn un1 , . . . , hn unN ) converges strongly to a global weak solution (h, hu1 , . . . , huN ) of System (1.1) in the sense of Definition 1. More precisely, hn converges strongly in √ 3 d n n C((0, T ); L 2 (T )), h uα converges strongly in L2 ((0, T ); L2 (Td )) and the momentum mnα = hn unα converges strongly in L1 ((0, T ); L1 (Td )) for any T > 0.

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

183

Remark 4. The previous theorem can be generalised to the case Rd with 1 ≤ d ≤ 3 provided that the previous convergence are local in space. Furthermore in a future work inspired of [32], we shall prove the existence of such a sequence (hn , un )n∈N of global regular approximate solutions satisfying uniformly (2.1), (2.3) and (2.4) only in the case Td for technical reasons. Remark 5. In the paper of Audusse et al. [3] the authors claim that the modelling is physically relevant since they exhibit a classical energy. In the present work, in order to prove the stability of global weak solutions, we need in addition to obtain the so-called BD entropy which is another hint of the physical interest of the model. Remark 6. Let us emphasise that in Theorem 1.1, it seems difficult to deal with the mass transfer flux, essentially because we are not able to prove that |Gnα+ 1 |(unα+1 − unα ) converges in the sense of distributions 2 to |Gα+ 1 |(uα+1 − uα ). Indeed it is not clear to prove the convergence almost everywhere of Gnα+ 1 . 2

2

In [14], with the choice uα+ 1 = 21 (uα + uα+1 ), the additional term |Gnα+ 1 |(unα+1 − unα ) does not appear. 2 2 Then it is easy to deal with the mass transfer term. For this specific choice for uα+ 1 , we are also able to 2 prove the BD entropy but it seems tricky to obtain a gain of integrability ` a la Mellet–Vasseur. For this reason we do not have enough compactness information to treat the convection term. We could obtain global weak solutions for the system proposed in [14] if we consider friction terms of the 1+ϵ form h|uα | uα with ϵ > 0 in each layer. Indeed in this case the friction terms ensure directly a gain of integrability on the velocity. The paper unfolds as follows. In Section 2, we give new estimates for System (1.1) involving the BD entropy and some gain of integrability on the velocity uα . In Section 3, we show the stability of global weak solutions following the arguments developed in [26]. Proofs of the BD entropy and of the gain of integrability on the velocities uα are postponed to Appendix. 2. A priori energy estimates In this section, we are interested in proving at least heuristically different energy estimates: the classical energy of the system, the BD entropy (see [5]) which is less obvious and an equivalent of the Mellet–Vasseur estimate from [26]. All the proofs are transferred in Appendix. 2.1. Classical energy Proposition 2.1. Let (h, u1 , . . . , uN ) be a classical solution of System (1.1). Then, the following equality holds: ∫ N ∫ N ∫ ∑ ∑ d 2 2 E dx + 4νh|D(uα )| dx + N κ|uα+1 − uα | dx dt Td d d α=1 T α=1 T N ∫ ∑ N 2 + (2.1) |uα+1 − uα | |Gα+ 1 | dx = 0 2 2 α=1 Td with 1 E= 2

( 2

N gh +

N ∑ α=1

) 2

h|uα |

+ N ghzb .

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

184

2.2. BD entropy Unfortunately the energy estimate (2.1) is not sufficient in order to prove the existence of global weak solutions. Indeed we need additional compactness information to deal with the pressure term and the convection terms. As in [5], we would like to prove that a BD entropy estimate is satisfied. Let us introduce as in [18,19] the effective velocity vα = uα + 4ν∇ log h. Then System (1.1) can be written ⎧ ∂ h + div(h¯ u) = 0, ⎪ ⎪ t g ⎪ 2 ⎪ ⎪ ⎨h [∂t vα + (uα · ∇)vα ] − 2ν div (h curl vα ) + 2 ∇h = −gh∇zb ( ) ( ) (2.2) Gα+ 1 − Gα− 1 ⎪ 2 2 ⎪ + N G (u − u ) − G (u − u ) + 4νN h∇ 1 1 1 1 α α ⎪ α+ 2 α− 2 α− 2 α+ 2 ⎪ h ⎪ ⎩ + κ(uα+1 − uα ) − κ(uα − uα−1 ), with curl vα = (∇vα − ∇T vα ) the vorticity. Proposition 2.2 (BD Entropy). If we assume that (h, u1 , . . . , uN ) is a smooth solution of System (1.1), then ∫ ∫ ∫ ∫ N ∑ 1 d Ng d d 2 2 h |vα | dx + h2 dx + N g zb h dx + 2ν h|curl vα | dx 2 dt Td α=1 2 dt Td dt Td d T ∫ ∫ N ∫ ∑ 2 2 + 4N νg |∇h| dx + 4N νg ∇zb · ∇h dx + κ |vα+1 − vα | dx Td

Td

α=1

Td

N ∫ N −1 ∫ N ∑ ∑ ( )2 1 N ∑ 2 |Gα+ 1 | |vα+1 − vα | dx + div(h(uα − uj )) dx = 0. + 2 2 2 α=1 Td hN d α=1 T j=α+1

(2.3)

Remark 7. From this entropy √ we deduce two new pieces of information which are√essential to obtain √ convergence results. Firstly, hvα is bounded in L∞ (0, T ; L2 (Td )). We deduce that h∇ log h = 2∇ h is bounded in L∞ (0, T ; L2 (Td )). This is the crucial point ensured by the BD entropy. On the other hand, √ thanks to (1.2), the last term of this estimate also gives a bound for Gα+ 1 / h in L2 (0, T ; L2 (Td )) that 2 enables to give sense to the term uα+ 1 Gα+ 1 . 2

2

Remark 8. We mention that we can also obtain energy and BD entropy with the choice for uα+ 1 = 2 1 2 (uα + uα+1 ) used in [14] (see Appendix A.4). 2.3. Mellet–Vasseur logarithmic estimate In order to deal with the convection term huα ⊗ uα which is only bounded in L∞ (0, T ; L1 (Td )), it is important to get a gain of integrability on the velocity as in [26]. Let us mention that in [5], in order to overcome this difficulty, the authors need to work with a friction term. We have the following result. Proposition 2.3. If we assume that (h, u1 , . . . , uN ) is a smooth solution of System (1.1), then ] ) ∫ [ ∫ N ( 2 ( ) [ ( )] ∑ d 1 + |uα | 2 2 2 h log 1 + |uα | dx + 3ν h 1 + log 1 + |uα | |D(uα )| dx dt Td 2 Td α=1 2−δ ( (∫ ) (∫ ) (∫ ) 2δ N [ ( )] 2 ∑ 2 6−δ δ 2 2 2−δ ≤ C h|∇uα | dx + C h dx × h 2 + log 1 + |uα | dx Td Td Td α=1 ) ∫ 2 ] 1 + |uα | [ 2 +g h 1 + log(1 + |uα | ) |∇zb | dx (2.4) 2 Td for any δ ∈ (0, 2) and for some constant C ≥ 0.

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

185

Remark 9. Let us mention that it seems difficult to obtain a similar result with the choice for uα+ 1 used 2 in [14].

3. Stability of global weak solutions Let us assume that there exists a sequence of global approximate solutions (hn , unα )n∈N satisfying uniformly in n all the following estimates for every α and all T > 0 with C > 0 depending on the initial data (h0 , uα,0 ): hn (t, x) > 0 almost everywhere on (0, +∞) × Td ,

(3.1)

and  √  n n  h uα 

≤ C,

(3.2a)

∥hn ∥L∞ (0,T ;L2 (Td )) ≤ C, √   n n ≤ C,  h ∇uα  2 2 d

(3.2b)

∥∇hn ∥L2 (0,T ;L2 (Td )) ≤ C,  √    ≤ C. ∇ hn  ∞ 2 d

(3.2d)

L∞ (0,T ;L2 (Td ))

(3.2c)

L (0,T ;L (T ))

L

(3.2e)

(0,T ;L (T ))

The initial data satisfy the following conditions: hn0 is bounded in L2 (Td ), hn0 ≥ 0 a.e. in Td , 2

2

hn0 |unα,0 | = |mnα,0 | /hn0 is bounded in L1 (Td ), √ ∇ hn0 is bounded in L2 (Td ), ∫ 2 1 + |unα,0 | 2 log(1 + |unα,0 | )dx ≤ C. hn0 2 Td

(3.3)

The proof of the stability of the sequence (hn , un )n∈N follows the same arguments as in [26]. We adapt them to our case. Step 1: Convergence of

√

hn .

Lemma 3.1. We have that for any T > 0: √

( ) hn is bounded uniformly in L∞ 0, T ; H 1 (Td ) , √ ( ) ∂t hn is bounded uniformly in L2 0, T ; H −1 (Td ) . As a consequence, up to a subsequence, √

(√ √

) ( ) hn converges a.e. and strongly in L2 0, T ; L2 (Td ) . We write

h a.e. and L2 ((0, T ) × Td ) strong. ( ) 3 Moreover, (hn ) converges strongly to h in C [0, T ]; L 2 (Td ) . hn −→

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

186

√ √ Proof . hn is uniformly bounded in L∞ ((0, T ), H 1 (Td )) due to (3.2b) and (3.2e). The estimate on ∂t hn can be deduced from the mass equation since ) (√ √ 1√ n ∂t hn = h div¯ un − div hn u ¯n . 2 The first term on the right hand side √ is bounded in L2 (0, T ; L2 (Td )) and the second term is bounded in ∞ −1 d L (0, T ; H (T )), so it implies that ∂t hn is uniformly bounded in L2 ((0, T ), H −1 (Td )). The Aubin–Lions’ lemma gives directly the strong convergence in L2 ((0, T ) × Td ). √ 3 To prove the convergence in C([0, T ]; L 2 (Td )) we first deduce by Sobolev embeddings that hn is bounded in the space L∞ (0, T ; L6 (Td )). We deduce that √ N 3 hn ∑ (√ n n ) h uα is bounded in L∞ (0, T ; L 2 (Td )). hn u ¯n = N α=1 √ √ 3 The continuity equation thus yields ∂t hn bounded in L∞ (0, T ; W −1, 2 (Td )) and since ∇hn = 2 hn ∇ hn 3 we have hn bounded in L∞ ((0, T ), W 1, 2 (Td )). From the Aubin–Lions’ lemma we deduce that hn converges 3 to h up a subsequence in C([0, T ]; L 2 (Td )). □

Step 2: Convergence of the pressure. Lemma 3.2. The pressure (hn )2 is bounded in Lr ((0, T )×Td ) for all r ∈ [1, 2). In particular, (hn )2 converge to h2 strongly in Lr ((0, T ) × Td ) for every T > 0 and 1 ≤ r < 53 . Proof . From inequalities (3.2b) and (3.2d) we deduce that hn ∈ L2 (0, T ; H 1 (Td )). Hence hn ∈ L2 (0, T ; L6 (Td )). In addition hn is bounded in L∞ ((0, T ), L2 (Td )). By interpolation hn is bounded in 10 L 3 ((0, T ) × Td ). We conclude recalling that (hn )2 converges almost everywhere to h2 and is uniformly 5 bounded in L 3 ((0, T ) × Td ), it implies then the strong convergence of (hn )2 to h2 in Lr ((0, T ) × Td ) for 1 ≤ r < 53 . □

Step 3: Bound for

√

hn unα .

2

2

Lemma 3.3. hn |unα | log(1 + |unα | ) is bounded in L∞ (0, T ; L1 (Td )). Proof . We use the inequality given by the Mellet–Vasseur approach and it implies that : ∫ ∫ N 2 ∑ d 1 + |unα | 2 2 2 hn log(1 + |unα | )dx + 3νhn (1 + log(1 + |unα | ))|D(unα )| dx dt 2 d d T T α=1 2−δ (∫ ) (∫ ) 2δ N ∑ 2 6−δ n 2−δ n 2 2 n ≤C +C (h ) dx × [2 + log(1 + |uα | )] δ h dx +g

α=1 N ∫ ∑ Td

α=1

Td

hn

Td

1 + (unα )2 2 (1 + log(1 + |unα | ))∇zb dx 2

for any δ ∈ (0, 2). Since hn is uniformly bounded in L (∫ C Td

6−δ

(hn ) 2−δ dx

) 2−δ 2

(∫ × Td

10 3

((0, T ) × Td ) we deduce that for δ small enough: 2

2

[2 + log(1 + |unα | )] δ hn dx

) 2δ ≤ C.

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

187

Since zb is bounded in W 1,∞ (Td ), applying Gr¨ onwall’s lemma we deduce that N ∫ ∑ α=1

Td

2

hn

1 + |unα | 2 log(1 + |unα | )(t)dx ≤ CT , 2

∀t ∈ (0, T ).

□

Step 4: Convergence of the momentum. Lemma 3.4. Up to a subsequence, the momentum mnα = hn unα converges strongly in L2 (0, T ; Lp (Td )) to some mα (t, x) for all p ∈ [1, 2). In particular hn unα −−−−−→ mα almost everywhere in Td × (0, T ). n→+∞

√ √ √ ∞ 6 d n n T, L we Proof . We have hn unα = hn hn unα . Since hn is bounded √ √in L (0, √ (T )) √ deduce that h uα 3 ∞ d n n n n n n n n is bounded in L (0, T ; L 2 (T )). Next, since ∇(h uα ) = h h ∇uα + 2 h uα ∇ h , we obtain that ∇(hn unα ) is bounded in L2 ((0, T ), L1 (Td )). In particular, we have hn unα is bounded in L2 (0, T ; W 1,1 (Td )). Let us bound now ∂t (hn unα ) in order to apply the Aubin–Lions lemma. Let us consider the momentum equation (1.1)-b,c,d). We have then the uniform bounded estimates using in particular the fact that hn unα is bounded in L2 (0, T ; W 1,1 (Td )): √ √ div(hn unα ⊗ unα ) = div( hn unα ⊗ hn unα ) ∈ L∞ (0, T ; W −1,1 (Td )) ( ) ) 1( n div(hn Dunα ) = div D(hn unα ) − uα ⊗ ∇hn + (unα ⊗ ∇hn )T ∈ L2 (0, T ; W −1,1 (Td )) 2 ∇(hn )2 ∈ L∞ (0, T ; W −1,1 (Td )). Now using Remark 7, we have: unα+ 1 Gnα+ 1 2 2

Gnα+ 1 √ = √ 2 hn unα+ 1 ∈ L2 (0, T ; L1 (Td )). 2 hn

In addition we know that unα − unα+1 is uniformly bounded in L2 ((0, T ) × Td ) and hn ∇zb is bounded in L∞ ((0, T ), L2 (Td )). In conclusion ∂t (hn unα ) is uniformly bounded in L2 (0, T ; W −1,1 (Td )) and using the Aubin–Lions’ lemma we deduce that hn unα converges strongly in L2 ((0, T ), W s,1 (Td )) for −1 ≤ s < 1. By Sobolev embeddings we obtain what we wish. □ Note that we can define uα (t, x) = mα (t, x)/h(t, x) over E = {(t, x); h(t, x) > 0} but uα (t, x) is not uniquely defined in the vacuum set E c . In order to define properly uα over {h = 0} we have to study the weak limit of the term unα+1 − unα . Step 5: Convergence of unα+1 − unα . We know via the energy estimate (2.1) that (unα+1 − unα ) is uniformly bounded in L2 ((0, T ) × Td ). Then − unα ) converges weakly in L2T (L2 ) to µα up to a subsequence. For any d ≥ 1, we have due to ∞ 1{h=0} ∈ L∞ T (L ) (unα+1

(unα+1 − unα )1{h=0} −−−−−→ µα 1{h=0} in D′ ((0, T ) × Td ). n→+∞

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

188

Since unα converges almost everywhere to uα on {h > 0}. We have then: (unα+1 − unα )1{h>0} −−−−−→ (uα+1 − uα )1{h>0} = µα 1{h>0} in D′ ((0, T ) × Td ). n→+∞

We assume now that u1 = 0 on {h = 0}. This choice defines all the values of uα over {h = 0}. Indeed we have u1 = 0 over {h = 0} and by iteration u2 = µ1 over {h = 0} and so on. Step 6: Convergence of

√

hn unα .

√ √ Lemma 3.5. The quantity hn unα converges strongly in L2 ((0, T ) × Td ) to mα / h (defined to be zero when h = 0). In particular, we have mα (t, x) = 0 a.e. over E c and there exists a function uα (t, x) such that mα (t, x) = h(t, x)uα (t, x) and √ √ hn unα −→ huα strongly in L2 ((0, T ) × Td ). √ Proof . Since mnα / hn is bounded in L∞ (0, T ; L2 (Td )), Fatou’s Lemma yields for almost every t ∈ (0, T ) ∫ ∫ (mnα )2 (mnα )2 (t)dx ≤ lim inf (t)dx < ∞. lim inf n n→+∞ h hn Td n→+∞ Td In particular, we have mα (t, x) = 0 a.e. over {h(t, x) = 0}. So, if we define the limit velocity uα (t, x) by setting uα (t, x) = mα (t, x)/h(t, x) when h(t, x) ̸= 0 and uα (t, x) = 0 when h(t, x) = 0, we have mα (t, x) = h(t, x)uα (t, x) and ∫ Td

m2α dx = h

∫

2

h|uα | dx < ∞. Td

Moreover, Fatou’s lemma yields that for almost every t ∈ (0, T ) ∫ ∫ 2 2 2 2 h|uα | log(1 + |uα | )dx = h|uα | log(1 + |uα | )dx d T {h>0} ∫ ∫ 2 2 2 n n n 2 = lim inf h |uα | log(1 + |uα | )dx ≤ lim inf hn |unα | log(1 + |unα | )dx. {h>0} n→+∞

n→+∞

mn

Td

mn

Let us point out that since unα = hnα has a limit over {h > 0} which is uα and in addition hnα is well defined 2 2 because hn > 0 almost everywhere. We deduce that h|uα | log(1 + |uα | ) is in L∞ (0, T ; L1 (Td )). √ n n n n Next, everywhere, it√is readily seen that over {h(t, √ since mα and h converge almost √ √ x) ̸= 0}, h uα = n n n n mα / h converges almost everywhere to huα = mα / h (we observe that mα / h has a sense since hn > 0 almost everywhere). Moreover, we have √ √ hn unα 1{|unα |≤M } −−−−−→ huα 1{|uα |≤M } almost everywhere. n→+∞

As a matter of fact, the convergence holds almost everywhere over {h(t, x) ̸= 0}, and over {h(t, x) = 0}, we √ √ have hn unα 1{|unα |≤M } ≤ M hn → 0. To conclude the proof of the lemma, for M > 0, there exists C > 0 such that: ∫ T∫ √ ∫ T∫ √ √ √ 2 2 | hn unα − huα | dx ≤ C | hn unα 1{|unα |≤M } − huα 1{|uα |≤M } | dx d d 0 T 0 T ∫ T∫ √ ∫ T∫ √ 2 2 n +C | hn uα 1{|unα |≥M } | dx + C | huα 1{|uα |≥M } | dx. 0

Td

0

Td

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

We observe that: ∫ T∫ Td

0

√ 2 | hn unα 1|unα |≥M | dx ≤

1 log(1 + M 2 )

T

∫

∫

2

Td

0

189

2

hn |unα | log(1 + |unα | )dx

and T

∫

√ 2 | huα 1|uα |≥M | dx ≤

∫ Td

0

∫

1 log(1 + M 2 )

2

2

h|uα | log(1 + |uα | )dx. Td

We have now: ∫ 0

T

√ √ 2 | hn unα 1{|unα |≤M } − huα 1{|uα |≤M } | dx d T (∫ ∫ T √ √ 2 | hn unα 1{|un |≤M,√hn ≤M } − huα 1{|uα |≤M,√h≤M } | dx ≤C

∫

∫

T

∫ |

+

√

Td

0

(∫

T

∫

≤C ∫

T

∫ | T

√

∫

≤C

2

√ √ 2 | hn unα 1{|un |≤M,√hn ≤M } − huα 1{|uα |≤M,√h≤M } | dx α

Td

T

√ √ |( hn − h)unα 1{|un |≤M,√hn >M } |dx

T

∫

Td

α

|

+

√

Td

0 T

√

∫

+ Td

)

− uα 1{|uα |≤M,√h>M } )| dx

∫

0

0

h(unα 1{|un |≤M,√hn >M } α

0

+

dx

α

Td

(∫

)

α

√ √ |( hn − h)unα 1{|un |≤M,√hn >M } |dx

Td

2

√ √ 2 | hn unα 1{|un |≤M,√hn ≤M } − huα 1{|uα |≤M,√h≤M } | dx

Td

0

∫

huα 1{|uα |≤M,√h>M } |

∫

+

∫

−

√

0

0

∫

hn unα 1{|un |≤M,√hn >M } α

T

+ ∫

α

Td

0

h1{ h>M } (unα 1{|unα |≤M } √

2

)

− uα 1{|uα |≤M } )| dx

h|unα 1{|unα |≤M } | |1{|√h⟩M } − 1{|√hn ⟩M } |dx.

The first term on the right hand side converges to 0 when M goes to +∞ by dominated convergence. The √ √ second term converges to 0 when n goes to +∞ since hn converges strongly to h in L2 ((0, T ) × Td ). The third and fourth term converge to 0 when M goes to +∞ when we apply the Chebyshev lemma. We deduce that: ∫ √ √ 2 lim sup | hn unα − huα | dx = 0. □ n→+∞

Td

Step 7: Convergence of the diffusion terms. Lemma 3.6. We have div(hn ∇unα ) −→ div(h∇u) in D′ ((0, T ) × Td ) div(hn ∇T unα ) −→ div(h∇T u) in D′ ((0, T ) × Td ).

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

190

Proof . Let ϕ(t, x) a test function, then ∫ T∫ ∫ T∫ ∫ div(hn ∇unα )ϕ dx = − hn ∇unα : ∇ϕ dx = 0

Td

Td

0

T

0

hn unα

∫ Td p

2

(∇hn · ∇ϕ) · unα dx +

∫ 0

T

∫ Td

hn unα · ∆ϕdx.

d

Thanks to Lemma 3.4, converges strongly in L (0, T ; L (T )) for 1 ≤ p <√2. This √ is nenough to prove n n n· the convergence of the second term. For the first term, we have ∇h · u = 2∇ h hn uα , we know that α √ √ n 2 d 2 n n h uα converges strongly in L ((0, T ) × T ) and ∇ h converges weakly in L ((0, T ) × Td ) then ∇hn · unα converges in the sense of distributions to ∇h · uα . Step 8: Convergence of Gnα+ 1 unα+ 1 2 2 Let us recall that we have: Gnα+ 1 = 2

α N ( ) 1 ∑ ∑ div hn (unj − uni ) . 2 N j=1 i=α+1

√ √ √ √ Since we know that hn unj converges strongly to huj in L2t,x and hn converges strongly to h we deduce that Gnα+ 1 converges in the sense of distributions to Gα+ 1 with: 2

2

Gα+ 1 = 2

α N ( ) 1 ∑ ∑ div h(uj − ui ) . 2 N j=1 i=α+1

We recall that we have: Gnα+ 1 unα+ 1 = 2

2

1 n 1 G 1 (un + unα+1 ) − |Gnα+ 1 |(unα − unα+1 ). 2 2 α+ 2 α 2

(3.4)

Let us consider the first term on the right hand side. We are going to show that Gnα+ 1 (unα + unα+1 ) converges 2 in the sense of distributions to Gα+ 1 (uα + uα+1 ). Let us take φ a C ∞ function with compact support in 2 (0, T ) × Td , we have then: ∫ 0

T

∫ Td

=−

Gnα+ 1 unα · φ dxdt = 2

∫ T∫ N α ( ) 1 ∑ ∑ div hn (unj − uni ) unα · φ dxdt 2 N j=1 i=α+1 0 Td

∫ T∫ √ α N √ √ √ ( ) 1 ∑ ∑ hn (unj − uni ) · hn ∇unα φ − hn (unj − uni ) · ∇φ hn unα dxdt. 2 N j=1 i=α+1 0 Td

√ √ √ ∞ 2 Using the fact that hn unα converges strongly to hu in L (L ), that hn ∇un converges weakly up to α T √ a subsequence in L2T (L2 ) to h∇u (indeed we have this convergence also in the sense of distributions), we deduce that: ∫ T∫ Gnα+ 1 unα · φ dxdt 0

2

Td

∫ T∫ √ α N √ √ √ ) ( 1 ∑ ∑ h(uj − ui ) · h∇uα φ − h(uj − ui ) · ∇φ huα dxdt 2 n→+∞ N j=1 i=α+1 0 Td ∫ T∫ −−−−−→ Gα+ 1 uα · φ dxdt. −−−−−→ −

n→+∞

0

Td

2

We proceed similarly for the term Gnα+ 1 unα+1 . Let us write now the second term on the right hand side of 2 (3.4) as follows: |Gnα+ 1 |(unα − unα+1 ) = 1{h=0} |Gnα+ 1 |(unα − unα+1 ) + 1{h>0} |Gnα+ 1 |(unα − unα+1 ). 2

2

2

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

191

⏐ n ⏐ ⏐G 1 √ ⏐ 3 ⏐ √α+ 2 ⏐ = ⏐ hn hn ⏐ is uniformly bounded in L2T (L 2 ). We are going to prove now that We know that ⏐ ⏐ 1{h=0} |Gnα+ 1 | converges strongly to 0 in L1T (L1 ). We have then: |Gnα+ 1 | 2

2

T

∫

∫

∫

1{h=0} |Gnα+ 1 | 2 Td

0

≤

α ∑

N ∑

∫

j=1 i=α+1

0

T

∫

T

α N ∑ ∑

∫

dxdt ≤ 0

(

Td

Td

( ) 1{h=0} |div hn (unj − uni ) | dxdt

j=1 i=α+1

√ √ √ √ ) 1{h=0} hn | hn div(unj − uni )| + 2 1{h=0} |∇ hn | hn |unj − uni | dxdt. (3.5)

√ √ Since | hn div(unj −uni )| is uniformly bounded in L2T (L2 ) and hn 1{h=0} converges strongly to 0 in LpT (L6−ϵ ) for any p ≥ 2 and ϵ > 0, we obtain using H¨ older’s inequality that: ∫ T∫ α N ∑ ∑ √ √ 1{h=0} hn | hn div(unj − uni )| dxdt −−−−−→ 0. j=1 i=α+1

0

n→+∞

Td

Let us estimate now the second term on the right hand side of (3.5), we have then: ∫ T∫ ∫ T∫ √ √ √ √ 1{h=0} |∇ hn | hn |unj − uni |dxdt ≤ 1{h=0} |∇ hn | hn (|unj | + |uni |)dxdt. 0

Td

0

Td

unj ,

Let us consider simply the term in we have then: ∫ T∫ ∫ T∫ √ √ √ √ n n n 1{h=0} |∇ hn | hn |uni |1{|uni |≤M } dxdt 1{h=0} |∇ h | h |ui |dxdt ≤ 0 Td 0 Td ∫ T∫ √ √ + 1{h=0} |∇ hn | hn |uni |1{|uni |>M } dxdt 0 Td ∫ T∫ √ √ ≤M 1{h=0} |∇ hn | hn dxdt 0 Td ∫ T∫ √ √ 1 2 1 + 1{h=0} |∇ hn | hn |uni | log(1 + |uni | ) 2 dxdt. 1 2 (log(1 + M )) 2 0 Td √ The first term on the √ right hand side goes to 0 when n goes to +∞ since hn 1{h=0} converges strongly to 2 0 in LpT (L6−ϵ ) and ∇ hn is uniformly bounded in L∞ T (L ). The second term goes also to 0 when M goes to +∞ and because the integral is uniformly bounded using the Mellet–Vasseur inequality. It proves that: |Gnα+ 1 |1{h=0} −−−−−→ 0 in L1 ((0, T ) × Td ). 2

(3.6)

n→+∞

We have seen that: 1 n 1 G 1 (un + unα+1 ) −−−−−→ Gα+ 1 (uα + uα+1 ) in D′ ((0, T ) × Td ). 2 n→+∞ 2 2 α+ 2 α |Gn

α+ 1

|

In addition we know that √hn2 is uniformly bounded in L2T (L2 ), it implies that up to a subsequence it √ converges to Mα+ 1 in L2T (L2 ). In addition we know that hn (unα − unα+1 ) converges strongly in L2T (L2 ) to 2 Gα+ 1 (uα − uα+1 ). We have then: 2

√ 1 1 n |Gα+ 1 |(unα + unα+1 ) −−−−−→ Mα+ 1 h(uα + uα+1 ) in D′ ((0, T ) × Td ). 2 n→+∞ 2 2 2 Finally we proved that: √ 1 1 Gnα+ 1 unα+ 1 −−−−−→ Gα+ 1 (uα + uα+1 ) − Mα+ 1 h(uα − uα+1 ) in D′ ((0, T ) × Td ). 2 2 2 2 n→+∞ 2 2

□

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

192

Acknowledgements The work presented in this paper was supported in part by CNRS-INSU, TelluS-INSMI-MI program, project CORSURF. It was realised during the secondment of the first and second authors in the ANGE Inria team. The second author has been partially funded by the ANR project INFAMIE ANR-15-CE40-0011. Appendix A Hereafter, A : B =

2

∑

denotes the scalar product upon matrices and |A| = A : A.

i,j Aij Bij

A.1. Proof of Proposition 2.1 We follow here the arguments of [3]. The main difficulty concerns the coupling between the different equalities through the flux terms. Simplifications arise only after summing the equations. Multiplying the momentum equations (1.1) by uα and summing over α, we obtain: ) ∑N ∫ ∑N ∫ ( G − u G · uα dx = 1. [∂ (hu ) + div(hu ⊗ u )] · u dx − N u 1 1 1 1 t α α α α d d α=1 T α=1 T α− α− α+ 2 α+ 2 ∑N ∫ ∑N ∫ ∑N2 ∫ 2( 2 2 N 1 d − Gα− 1 ) dx − N α=1 Td uα+ 1 Gα+ 1 − 1 α=1 Td h|uα | dx + α=1 Td |uα | (Gα+ 2 2 dt 2 2 2 2 ) uα− 1 Gα− 1 · uα dx; ∫ ∑N2 ∫ 2 d 2. 12 α=1 Td huα · ∇h2 dx = N g2 dt h2 dx; Td ∫ ∑N ∫ ∑N 2 3. )| dx; d uα · div (4νhD(uα )) dx = − Td 4νh α=1 α=1 |D(u T ∫ ∑N ∑αN ∫ 2 4. α − uα−1 ) · uα ] dx = − α=1 ∫Td [(uα+1 − uα ) · uα − (u α=1 Td |uα+1 − uα | dx; ∫ ∑N d 5. α=1 Td huα · ∇zb dx = N dt Td zb h dx. Let us observe now that: N ∫ N ∫ ( ) ( ) ∑ N ∑ 2 2 Gα+ 1 |uα | − |uα+1 | dx − N uα+ 1 Gα+ 1 − uα− 1 Gα− 1 · uα dx 2 2 2 2 2 2 α=1 Td d α=1 T ∫ N N ∑ 2 = |G 1 | |uα+1 − uα | dx. 2 α=1 Td α+ 2

(A.1)

Combining the different previous estimates and (A.1), we obtain the energy estimate (2.1). A.2. Proof of Proposition 2.2 A.2.1. A useful identity We aim at proving the identity n−1 i ∑∑

∀ n ≥ 2,

n ∑

(uj − uk )(ui − ui+1 ) =

i=1 j=1 k=i+1

n−1 ∑

n ∑

(ui − uk )2 .

i=1 k=i+1

This is equivalent to showing En :=

n−1 ∑

n ∑

i=1 k=i+1

⎡ ⎤ i ∑ ⎣ (uj − uk )(ui − ui+1 ) − (ui − uk )2 ⎦ = 0.

(A.2)

j=1

Let us first notice that (k ≥ i) due to a telescoping procedure (ui − uk )2 = (ui − uk )

k−1 ∑ j=i

(uj − uj+1 ).

(A.3)

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

193

Inserting (A.3) into (A.2) and switching series twice, we get ⎡ ⎤ n−1 n i k−1 ∑ ∑ ∑ ∑ ⎣ (uj − uk )(ui − ui+1 ) − (ui − uk ) En = (uj − uj+1 )⎦ i=1 k=i+1

j=1

j=i

⎡ ⎤ n k−1 i k−1 ∑ ∑ ∑ ∑ ⎣ (uj − uk )(ui − ui+1 ) − (ui − uk ) (uj − uj+1 )⎦ = k=2 i=1 n ∑

=

j=1

j=i

⎡

k−1 i ∑∑

⎣

(uj − uk )(ui − ui+1 ) −

k−1 ∑

i=1 j=1

k=2

(ui − uk )

i=1

k−1 ∑

⎤ (uj − uj+1 )⎦

j=i

⎡ ⎤ n k−1 k−1 ∑ ∑ k−1 ∑ ∑ k−1 ∑ ⎣ = (uj − uk )(ui − ui+1 ) − (ui − uk )(uj − uj+1 )⎦ j=1 i=j

k=2

i=1 j=i







=0

which ends the proof. A.2.2. Proof of Proposition 2.2 Multiplying the momentum equations of (2.2) by vα , integrated over Td and summing over α we get: ] ∑N ∫ [ d ∑N ∫ 2 2 2 1 h|v | + N G (|v | − |v | ) dx; 1. h (∂ v + (u · ∇)v ) · v dx = 1 α α+1 α t α α α α d d α=1 T α=1 T α+ 2 2 dt ∫ ∫ ∑N g ∫ 2 Ng d 2 2 2. α=1 2 Td vα · ∇h dx = 2 dt Td h dx + 4N νg Td |∇h| dx; ∫ ∫ ∫ ∑N 3. α=1 g Td h∇zb · vα dx = N g Td zb h dx + 4gN ν Td ∇zb · ∇hdx; ∑N ∫ ∑N ∫ 2 4. α=1 Td vα · 2νdiv(h curl vα ) dx = − α=1 Td 2νh|curl vα | dx; 5. Since for all α ∈ {1, . . . , N } we have vα − vα−1 = uα − uα−1 we deduce that: N ∫ N ∫ ∑ ∑ 2 κ (uα+1 − uα ) · vα − (uα − uα−1 ) · vα dx = −κ |vα+1 − vα | dx. Td

α=1

α=1

Td

Combining the previous equality we have: 1 d 2 dt

∫ h Td

N ∑

2

|vα | dx +

α=1

∫

∫

Ng d 2 dt

h2 dx + N g

Td

∫

Td

Td

∇zb · ∇h dx + κ Td

N ∑

∫

∫

Td

N ∫ ∑ α=1

N ∑

2

2νh|curl vα | dx

zb h dx +

∫

2

|∇h| dx + 4N νg

+ 4N νg

d dt

2

|vα+1 − vα | dx

Td

) 1 2 2 N Gα+ 1 (|vα+1 | − |vα | ) dx − N Gα+ 1 (uα+ 1 − uα ) − Gα− 1 (uα− 1 − uα ) · vα dx 2 2 2 2 2 2 α=1 d α=1 T ( ) N ∫ Gα+ 1 − Gα− 1 ∑ 2 2 · vα dx = 0. (A.4) − 4νN h∇ h d α=1 T +

(

Next we have due to G 1 = GN + 1 = 0: 2

N ∑

2

N ∑

( ) 1 2 2 Gα+ 1 (uα+ 1 − uα ) − Gα− 1 (uα− 1 − uα ) · vα Gα+ 1 (|vα+1 | − |vα | ) − 2 2 2 2 2 2 α=1 α=1 =

N N N −1 ∑ ∑ 1∑ 2 2 Gα+ 1 (|vα+1 | − |vα | ) − Gα+ 1 (uα+ 1 − uα ) · vα − Gα+ 1 (uα+1 − uα+ 1 ) · vα+1 2 2 2 2 2 2 α=1 α=1 α=0

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

194

N −1 ∑

[ ] 1 2 2 −Gα+ 1 uα+ 1 · (vα − vα+1 ) + (uα+1 · vα+1 − uα · vα ) − (|vα+1 | − |vα | ) 2 2 2 α=1 [ ] N −1 ∑ 1 2 2 2 2 = −Gα+ 1 uα+ 1 · (uα − uα+1 ) + |uα+1 | − |uα | + 4ν(uα+1 − uα ) · ∇ log h − (|vα+1 | − |vα | ) ; 2 2 2 α=1 [ ] N −1 ∑ 1 2 2 = −Gα+ 1 uα+ 1 · (uα − uα+1 ) + (|uα+1 | − |uα | ) , 2 2 2 α=1 ( ) N −1 ∑ 1 = −Gα+ 1 (uα − uα+1 ) · uα+ 1 − (uα+1 + uα ) , 2 2 2 α=1 =

since 2

2

|vα+1 | − |vα | = (uα+1 − uα ) · (uα+1 + uα + 8ν∇ log h) 2

2

= |uα+1 | − |uα | + 8ν(uα+1 − uα ) · ∇ log h. Then N −1 ∑

( −Gα+ 1 (uα − uα+1 ) · uα+ 1 2

α=1

2

1 − (uα+1 + uα ) 2

) =

N −1 ⏐ 1 ∑ ⏐⏐ ⏐ 2 ⏐Gα+ 1 ⏐ |uα+1 − uα | . 2 2 α=1

Finally we have proven that: N N ( ) ∑ 1∑ 2 2 Gα+ 1 (|vα+1 | − |vα | ) − Gα+ 1 (uα+ 1 − uα ) − Gα− 1 (uα− 1 − uα ) · vα 2 2 2 2 2 2 α=1 α=1

=

N |G ∑ α+ 1 | 2

α=1

2

2

|uα+1 − uα | .

From the relation (1.2) and by integration by parts, we have ( ) ( ) N ∫ N −1 ∫ Gα+ 1 − Gα− 1 Gα+ 1 ∑ ∑ 2 2 2 h vα · ∇ dx = h(vα − vα+1 ) · ∇ dx h h d d T T α=1 α=1 N −1 ∫ Gα+ 1 ∑ 2 div (h(uα − uα+1 )) dx =− h d α=1 T N −1 ∫ α N ∑ 1 ∑ ∑ =− div(h(uj − ui ))div(h(uα − uα+1 )) dx 2 d hN α=1 T j=1 i=α+1 N −1 ∫ N ∑ ∑ [ ]2 1 =− div(h(uα − uj )) dx. 2 d hN α=1 T j=α+1

(A.5)

(A.6)

A detailed proof of the latter relation is given in Appendix A.2.1. Combining (A.4), (A.5) and (A.6) implies the estimate (2.3). A.3. Proof of Proposition 2.3 Let us first rewrite System (1.1) under the non-conservative form ⎧ u) = 0, (a) ⎪ ⎨∂t h + div(h¯ g h [∂t uα + (uα · ∇)uα ] + ∇h2 = −gh∇zb + N Gα+ 1 (uα+ 1 − uα ) 2 2 2 ⎪ ⎩ + N Gα− 1 (uα − uα− 1 ) + div (4νhD(uα )) + κ(uα+1 − uα ) − κ(uα − uα−1 ). (b) 2

2

(A.7)

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

195

Let us set Φ(x) =

1 + x2 log(1 + x2 ) 2

[ ] ϕ(x) = Φ ′ (x) = x 1 + log(1 + x2 )

and

(A.8)

and notice that ) 2 ( ) [ ( )] 1 + |u| 2 2 log 1 + |u| = (u · ∂u) 1 + log 1 + |u| . ∂ 2 )] [ ( 2 We multiply Eq. (A.7)(b) by uα 1 + log 1 + |uα | and we integrate over Td . Hence each term becomes [ ] ) ] ∫ ∫ [ ( 2 1. Td huα · ∂t uα 1 + log(1 + |uα | ) dx = Td ∂t hΦ(|uα |) − Φ(|uα |)∂t h dx. 2. Using (1.5), we obtain by integration by parts: ∫ ∫ [ ] 2 huα · (uα · ∇)uα 1 + log(1 + |uα | ) dx = huα · ∇Φ(|uα |) dx Td ∫ Td [ ] = Φ(|uα |) ∂t h − N (Gα+1/2 − Gα−1/2 ) dx. (

Td

3. For the pressure term we apply the same approach as in [26] ⏐ ⏐∫ [ ] ⏐ ⏐ 2 2 ⏐ ⏐ ⏐ d 1 + log(1 + |uα | ) uα · ∇h dx ⏐ T ⏐ ⏐ ⏐∑ ∫ ⏐ ⏐ ⏐∫ [ ( )] ⏐ ⏐ ⏐ ⏐ 2u u αi αk 2 2 2 ⏐ ⏐ ⏐ ≤⏐ h h 1 + log 1 + |uα | (divuα )⏐⏐ dx 2 ∂i uαk dx ⏐ + ⏐ Td ⏐ i,j Td 1 + |uα | ⏐ ⏐ (∫ ) 21 (∫ ) 12 ⏐∫ [ ( )] ⏐ ⏐ 2 2 ≤2 h|∇uα | dx h2 1 + log 1 + |uα | h3 dx + ⏐⏐ (divuα )dx ⏐⏐ . Td

Td

Td

Since (divu)2 =

∑∑ i

∂i ui ∂j uj ≤

∑∑ 1 (

j

i

j

2

) 2 (∂i ui )2 + (∂j uj )2 ≤ d |D(u)|

where d is the dimension of the space, we have ⏐ ⏐∫ )] [ ( ⏐ ⏐ 2 2 ⏐ (divuα )dx ⏐⏐ ⏐ d h 1 + log 1 + |uα | T ) 12 (∫ (∫ [ ( )] [ ( )] ) 21 2 2 2 3 (divuα ) dx ≤ h 1 + log 1 + |uα | h 1 + log 1 + |uα | dx d d T ∫T [ ∫ )] )] ( [ ( 2 2 2 ≤ν h 1 + log 1 + |uα | |D(uα )| dx + Cν h3 1 + log 1 + |uα | dx Td

Td

according to the Young’s inequality with Cν = It follows that for Cν′ large enough: ⏐∫ [ ⏐ ] ⏐ ⏐ 2 2 ⏐ ⏐ ⏐ d 1 + log(1 + |uα | ) uα · ∇h dx ⏐ T (∫ ) ∫ [ ( )] 2 2 2 ≤ν h|∇uα | dx + ν h 1 + log 1 + |uα | |D(uα )| dx d d T ∫ T [ ( )] 2 ′ 3 + Cν h 1 + log 1 + |uα | dx. d 4ν .

Td

Finally, for any δ ∈ (0, 2), the last term is bounded by means of the H¨older’s inequality ∫ Td

(∫ [ ( )] 2 h 1 + log 1 + |uα | dx ≤ 3

Td

h

6−δ 2−δ

) 2−δ 2 dx

(∫ × Td

[ ( )] 2 ) 2δ δ 2 h 1 + log 1 + |uα | dx .

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

196

4. Likewise we have ∫ ∫ [ ] 2 huα · ∇zb 1 + log(1 + |uα | ) dx ≤ Td

Td

h

2 ] 1 + |uα | [ 2 1 + log(1 + |uα | ) |∇zb |dx. 2

1,∞

Since zb is assumed bounded in W this term can be treated by the Gr¨onwall’s lemma. 5. For the viscous terms we have ∫ [ ] 2 uα · div (4νhD(uα )) 1 + log(1 + |uα | ) dx Td ∫ [ ] ∑ ∫ uα uα · ∂j uα 2 2 = −4ν h 1 + log(1 + |uα | ) |D(uα )| dx − 8ν h i 2 Dij (uα )dx 1 + |uα | Td Td i,j and we have for Cα > 0 large enough: ∫ ∫ [ ] [ ] 2 2 2 uα · div (4νhD(uα )) 1 + log(1 + |uα | ) dx + 4ν h 1 + log(1 + |uα | ) |D(uα )| dx d d T T ∫ 2 ≤ Cα h|∇uα | dx. Td

6. For the friction terms (since by definition u0 = u1 and uN = uN +1 ) we have: N ∫ ∑ α=1

[ ] 2 [(uα+1 − uα ) · uα − (uα − uα−1 ) · uα ] 1 + log(1 + |uα | ) dx

Td

=−

N −1 ∫ ∑ α=1

[ [ ] [ ]] 2 2 (uα − uα+1 ) · uα 1 + log(1 + |uα | ) − uα+1 1 + log(1 + |uα+1 | ) dx ≤ 0

Td

since the function ϕ defined by (A.8) is increasing. Combining all the previous estimates, we have: ] ) ∫ [ ∫ N ( 2 ( ) [ ( )] ∑ 1 + |uα | d 2 2 2 h log 1 + |uα | dx + 3νh 1 + log 1 + |uα | |D(uα )| dx dt Td 2 Td α=1 ( ∫ N [ ] ∑ 2 + [Gα+ 1 (uα+ 1 − uα ) + Gα− 1 (uα − uα− 1 )] · uα 1 + log(1 + |uα | ) 2

Td

α=1

2

2

2

) 2 ( ) 1 + |uα | 2 log 1 + |uα | (Gα+1/2 − Gα−1/2 ) dx + 2 ) (∫ ) 2−δ (∫ N ( (∫ ∑ 2 6−δ 2 2−δ ≤ C h|∇uα | dx + C h dx × Td

α=1

Td

) 2δ )] 2 [ ( δ 2 h 2 + log 1 + |uα | dx

Td

2

) ] 1 + |uα | [ 2 +g h 1 + log(1 + |uα | ) |∇zb | dx . 2 Td ∫

We have since G 1 = GN + 1 = 0: 2

N ∑ α=1

∫ Td

2

[ ] 2 [Gα+ 1 (uα+ 1 − uα ) + Gα− 1 (uα − uα− 1 )] · uα 1 + log(1 + |uα | ) 2

2

2

2

2

( ) 1 + |uα | 2 log 1 + |uα | (Gα+ 1 − Gα− 1 )dx 2 2 2 N ∫ [ ] ∑ 2 = Gα+ 1 (uα+ 1 − uα ) · uα 1 + log(1 + |uα | ) dx +

α=1

Td

2

2

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

+

+

N ∫ ∑ d α=1 T N ∫ ∑

α=1

=

N ∑ α=1

+

∫ Td

Td

−

=

α=0 T N ∫ ∑

d α=1 T ∫ N −1 ∑

α=0 N −1 ∫ ∑ α=1

2

2

2

( ) 1 + |uα | 2 log 1 + |uα | (Gα+ 1 − Gα− 1 )dx 2 2 2

2

N −1 ∫ ∑

Td

[ ] 2 Gα− 1 (uα − uα− 1 ) · uα 1 + log(1 + |uα | ) dx

[ ] 2 Gα+ 1 (uα+ 1 − uα ) · uα 1 + log(1 + |uα | ) dx

d

+

197

Td

2

[ ] 2 Gα+ 1 (uα+1 − uα+ 1 ) · uα+1 1 + log(1 + |uα+1 | ) dx 2

2

2

( ) 1 + |uα | 2 log 1 + |uα | Gα+ 1 dx 2 2 2

( ) 1 + |uα+1 | 2 log 1 + |uα+1 | Gα+ 1 dx 2 2

Gα+ 1 1{G 2

α+ 1 2

≤0} (uα+1

2

− uα ) · uα+1 (1 + log(1 + |uα+1 | ))

[ ] 2 + Gα+ 1 1{G 1 ≥0} (uα+1 − uα ) · uα 1 + log(1 + |uα | ) α+ 2 2 ) ( 2 2 ( ) 1 + |u 1 + |uα | α+1 | 2 2 log 1 + |uα | − log(1 + |uα+1 | ) dx + Gα+ 1 2 2 2 N −1 ∫ ∑ 2 =− |Gα+ 1 |1{G 1 ≤0} (uα+1 − uα ) · uα+1 (1 + log(1 + |uα+1 | )) α=1

2

Td

α+ 2

[ ] 2 (uα − uα+1 ) · uα 1 + log(1 + |uα | ) 2 ( ) 2 2 ( ) 1 + |u 1 + |uα | α+1 | 2 2 + |Gα+ 1 |1{G 1 ≤0} log 1 + |uα | − log(1 + |uα+1 | ) dx α+ 2 2 2 2 ) ( 2 2 ( ) 1 + |uα | 1 + |uα+1 | 2 2 log 1 + |uα+1 | − log(1 + |uα | ) dx + |Gα+ 1 |1{G 1 ≥0} α+ 2 2 2 2 [ ] N −1 ∫ ∑ =− |Gα+ 1 | 1{G 1 ≤0} Ψ (uα , uα+1 ) + 1{G 1 ≥0} Ψ (uα+1 , uα ) dx ≤ 0. + |Gα+ 1 |1{G

≥0} α+ 1 2

α=1

Td

2

α+ 2

α+ 2

[ ] Indeed, the function Ψ : (x, y) ↦→ y(y − x) 1 + log(1 + y 2 ) + Φ(x) − Φ(y) satisfies Ψ (y, y) = 0 and ∂x Ψ (x, y) = ϕ(x) − ϕ(y). As the function ϕ defined by (A.8) is increasing, it shows that Ψ (x, y) ≥ 0. Finally, we obtain the estimate ∫ ∫ N [ ( )] ∑ d 2 2 hΦ(|uα |)dx + 3ν h 1 + log 1 + |uα | |D(uα )| dx dt d d T T α=1 (∫ (∫ ) (∫ ) 2−δ N [ ( )] 2 ) 2δ ∑ 2 6−δ δ 2 2 ′′ ′ 2−δ × h 1 + log 1 + |uα | ≤ Cν h|∇uα | dx + Cν h dx dx α=1

Td

Td

Td

] 1 + u2α [ 2 +g h 1 + log(1 + |uα | ) |∇zb | dx. 2 Td ∫

Remark 10. Let us notice that Ψ satisfies the following estimate ( )] ( )] (y − x)2 [ (y − x)2 [ 1 + log 1 + min(x, y)2 ≤ Ψ (x, y) ≤ 3 + log 1 + max(x, y)2 . 2 2

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

198

Indeed given (A.8) we check that ∫ ∫ 1 [ ′ ( )] Φ (y) − Φ ′ y + s(x − y) ds = (y − x)2 Ψ (x, y) = (y − x)

1

0

0

∫

1

( ) sΦ ′′ y + s(1 − t)(x − y) dsdt

0

2z 2 1+z 2

with Φ ′′ (z) = 1 + + log(1 + z 2 ). We obtain the bound above inserting min(x, y) ≤ y + s(1 − t)(x − y) ≤ max(x, y) for all (t, s) ∈ [0, 1]2 into the integral. We observe that the function Ψ (uα , uα+1 ) measures the distance between uα and uα+1 in the Mellet– Vasseur inequality (2.4). In particular, this provides a better estimate for uα+1 − uα in L2T (L2 ) according to ( ) the classical energy inequality. Indeed we get here a factor log 1 + min(uα , uα+1 )2 .

A.4. Energy and BD entropy when uα+ 1 = 21 (uα + uα+1 ) 2

We recall that G 1 = GN + 1 = 0. With this definition of uα+ 1 , the energy due to the term in Gα+ 1 2 2 2 2 computed in Section 2 gives N ∫ N ∫ ( ) ( ) ∑ N ∑ 2 2 Gα+ 1 |uα | − |uα+1 | dx − N uα+ 1 Gα+ 1 − uα− 1 Gα− 1 · uα dx 2 2 2 2 2 2 α=1 Td d α=1 T N ∫ ) ( N ∑ 2 2 = Gα+ 1 |uα | + |uα+1 | dx 2 2 α=1 Td N ∫ ) ( N ∑ − (uα+1 + uα )Gα+ 1 − (uα−1 − uα )Gα− 1 · uα dx 2 2 2 α=1 Td N ∫ N ∫ ( ) N ∑ N ∑ 2 2 Gα+ 1 |uα | − |uα+1 | dx − uα+1 · uα Gα+ 1 = 2 2 2 α=1 Td 2 α=1 Td 2

2

+ |uα | Gα+ 1 − uα−1 · uα Gα− 1 − |uα | Gα− 1 dx 2 2 2 N ∫ N ∫ ( ) ∑ N N ∑ 2 2 2 = Gα+ 1 |uα | − |uα+1 | dx − uα+1 · uα Gα+ 1 + |uα | Gα+ 1 dx 2 2 2 2 α=1 Td 2 α=1 Td N −1 ∫ N ∑ 2 + uα · uα+1 Gα+ 1 + |uα+1 | Gα+ 1 dx = 0 2 2 2 α=0 Td and the energy is then d dt

∫ E dx + Td

N ∫ ∑ α=1

2

4νh|D(uα )| dx + N

Td

N ∫ ∑ α=1

κ(uα+1 − uα )2 dx = 0

Td

with 1 E= 2

( 2

N gh +

N ∑

) hu2α

+ N gzb h.

α=1

For the BD-entropy, we have N N ( ) ∑ 1∑ 2 2 Gα+ 1 (uα+ 1 − uα ) − Gα− 1 (uα− 1 − uα ) · vα Gα+ 1 (|vα+1 | − |vα | ) − 2 2 2 2 2 2 α=1 α=1

=

N N N −1 ∑ ∑ 1∑ 2 2 Gα+ 1 (|vα+1 | − |vα | ) − Gα+ 1 (uα+ 1 − uα ) · vα + Gα+ 1 (uα+1 − uα+ 1 ) · vα+1 2 2 2 2 2 2 α=1 α=1 α=0

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

199

N −1 ∑

[ ] 1 2 2 −Gα+ 1 uα+ 1 · (vα − vα+1 ) + (uα+1 · vα+1 − uα · vα ) − (|vα+1 | − |vα | ) 2 2 2 α=1 [ ] N −1 ∑ 1 2 2 2 2 = −Gα+ 1 uα+ 1 · (uα − uα+1 ) + |uα+1 | − |uα | + 4ν(uα+1 − uα ) · ∇ log h − (|vα+1 | − |vα | ) 2 2 2 α=1 [ ] N −1 ∑ 1 2 2 = −Gα+ 1 uα+ 1 · (uα − uα+1 ) + (|uα+1 | − |uα | ) 2 2 2 α=1 [ ( )] N −1 ∑ 1 = −Gα+ 1 (uα − uα+1 ) uα+ 1 − (uα+1 + uα ) = 0. 2 2 2 α=1 =

Then we obtain ∫

N ∑

∫ ∫ ∫ Ng d d 2 2 |vα | dx + h h dx + N g zb h dx + 2νh|curl vα | dx 2 dt dt d d d d T T T T α=1 ∫ ∫ ∫ N ∑ 2 2 + 4N νg |∇h| dx + 4N νg ∇zb · ∇h dx + κ |vα+1 − vα | dx

1 d 2 dt

2

Td

+

N −1 ∫ ∑ α=1

Td

1 hN 2

Td N ∑

(

div(h(uα − uj ))

α=1

)2

Td

dx = 0.

j=α+1

References

[1] N. Aguillon, E. Audusse, E. Godlewski, M. Parisot, Hyperbolicity of the layer wise discretized hydrostatic Euler equations: The bilayer case, preprint, 2016. [2] E. Audusse, A multilayer Saint–Venant model: derivation and numerical validation, Discrete Contin. Dyn. Syst. Ser. B 5 (2) (2005) 189–214. [3] E. Audusse, M.-O. Bristeau, B. Perthame, J. Sainte-Marie, A multilayer Saint–Venant system with mass exchanges for shallow water flows. Derivation and numerical validation, Math. Model. Numer. Anal. 45 (1) (2011) 169–200. [4] F. Bouchut, V. Zeitlin, A robust well-balanced scheme for multi-layer shallow water equations, Discrete Contin. Dyn. Syst. Ser. B 13 (4) (2010) 739–758. [5] D. Bresch, B. Desjardins, Existence of weak solutions for 2d viscous shallow water equations and convergence to the quasi-geostrophic model, Comm. Math. Phys. 238 (1–2) (2003) 211–223. [6] D. Bresch, B. Desjardins, On the construction of approximate solutions for the 2D viscous shallow water model and for compressible Navier-Stokes models, J. Math. Pures Appl. 86 (2006) 362–368. [7] D. Bresch, B. Desjardins, E. Zatorska, Two-velocity hydrodynamics in fluid mechanics: Part II existence of global κ-entropy solutions to compressible Navier-Stokes systems with degenerate viscosities, J. Math. Pures Appl. 104 (4) (2015) 801–836. [8] D. Bresch, E. Jabin, Global existence of weak solutions for compresssible Navier–Stokes equations: Thermodynamically unstable pressure and anisotropic viscous stress tensor, arXiv:1507.04629. [9] M.-O. Bristeau, J. Sainte-Marie, Derivation of a non-hydrostatic shallow water model; Comparison with Saint–Venant and Boussinesq systems, Discrete Contin. Dyn. Syst. Ser. B 10 (4) (2008). [10] M. Castro, J. Mac´ıas, C. Par´ es, A Q-scheme for a class of systems of coupled conservation laws with source term. Application to a two-layer 1-D shallow water system, ESAIM Math. Model. Numer. Anal. 35 (1) (2001) 107–127. [11] V. Casulli, A semi-implicit finite difference method for non-hydrostatic, free-surface flows, Int. J. Numer. Methods Fluids 30 (4) (1999) 425–440. [12] B. Di Martino, M. Peybernes, P. Orenga, On a bi-layer shallow water problem with the rigid-lid hypothesis, Math. Models Methods Appl. Sci. 15 (2005) 843–869. [13] E. Feireisl, A. Novotny, Singular Limits in Thermodynamics of Viscous Fluids, in: Adv. Math. Fluid Mech., Birkha¨ user, 2009. [14] E.D. Fern´ andez-Nieto, E.H. Kon´ e, T. Chac´ on-Rebollo, A multilayer method for the hydrostatic Navier-Stokes equations: a particular weak solution, J. Sci. Comput. 60 (2) (2014) 408–437. [15] E.D. Fern´ andez-Nieto, E.H. Kon´ e, T. Morales De Luna, R. B¨ urger, A multilayer shallow water system for polydisperse sedimentation, J. Comput. Phys. 238 (1) (2013) 281–314. [16] J.-F. Gerbeau, B. Perthame, Derivation of viscous Saint–Venant system for laminar shallow water; numerical validation, Discrete Contin. Dyn. Syst. Ser. B 1 (1) (2001) 89–102.

200

B. Di Martino et al. / Nonlinear Analysis 163 (2017) 177–200

[17] A.E. Green, P.M. Naghdi, A derivation of equations for wave propagation in water of variable depth, J. Fluid Mech. 78 (02) (1976) 237–246. [18] B. Haspot, Existence of global strong solution for the compressible Navier-Stokes equations with degenerate viscosity coefficients in 1D, hal-01082319, arXiv:1411.5503. [19] B. Haspot, New formulation of the compressible Navier-Stokes equations and parabolicity of the density, hal-01081580. [20] J. Li, Z. Xin, Global existence of weak solutions to the barotropic compressible Navier-Stokes flows with degenerate viscosities, 2015, arXiv:1504.06826. [21] P.L. Lions, Mathematical Topics in Fluid Mechanics, in: Compressible Models, vol. 2, Oxford Science Publications, 1998. [22] P.J. Lynett, P.L.-F. Liu, Linear analysis of the multi-layer model, Coastal Eng. 51 (5) (2004) 439–454. [23] P.A. Madsen, R. Murray, O.R. Sørensen, A new form of the Boussinesq equations with improved linear dispersion characteristics, Coastal Eng. 15 (4) (1991) 371–388. [24] F. Marche, Derivation of a new two-dimensional viscous shallow water model with varying topography, bottom friction and capillary effects, Eur. J. Mech. B Fluids 26 (1) (2007) 49–63. [25] A. Mellet, A. Vasseur, Existence and uniqueness of global strong solutions for one-dimensional compressible Navier-Stokes equations, SIAM J. Math. Anal. 39 (4) (2007) 1344–1365. [26] A. Mellet, A. Vasseur, On the barotropic compressible Navier-Stokes equations, Comm. Partial Differential Equations 32 (1–3) (2007) 431–452. [27] M.L. Mu˜ noz-Ruiz, F.J. Chatelon, P. Orenga, On a bi-layer shallow-water problem, Nonlinear Anal. RWA 4 (1) (2003) 139–171. [28] G. Narbona-Reina, J.D.D. Zabsonr´ e, E.D. Fern´ andez-Nieto, D. Bresch, Derivation of a bilayer model for shallow water equations with viscosity. numerical validation, Comput. Model. Eng. Sci. 43 (1) (2009) 27–72. [29] O. Nwogu, Alternative form of Boussinesq equations for nearshore wave propagation, J. Waterway, Port, Coastal, Ocean Eng. 119 (6) (1993) 618–638. [30] D.H. Peregrine, Long waves on a beach, J. Fluid Mech. 27 (04) (1967) 815–827. [31] F. Serre, Contribution ` a l’´ etude des ´ ecoulements permanents et variables dans les canaux, La Houille Blanche (6) (1953) 830–872. [32] A. Vasseur, C. Yu, Existence of global weak solution for 3D degenerate compressible Navier-Stokes equations, Invent. Math. 206 (3) (2016) 935–974. [33] A. Vasseur, C. Yu, Global weak solutions to compressible quantum Navier-Stokes equations with damping, SIAM J. Math. Anal. 48 (2) (2016) 1489–1511. [34] J.D.D. Zabsonr´ e, G. Narbona-Reina, Existence of a global weak solution for a 2D viscous bi-layer Shallow Water model, Nonlinear Anal. RWA 10 (2009) 2971–2984. [35] E. Zatorska, On the flow of chemically reacting gaseous mixture, J. Differential Equations 253 (2012) 3471–3500.