Shallow water asymptotic models for the propagation of internal waves

We are interested in asymptotic models for the propagation of internal waves at the interface between two shallow layers of immiscible fluid, under the rigid-lid assumption. We review and complete existing works in the literature, in order to offer a unified and comprehensive exposition. Anterior models such as the shallow water and Boussinesq systems, as well as unidirectional models of Camassa-Holm type, are shown to descend from a broad Green-Naghdi model, that we introduce and justify in the sense of consistency. Contrarily to earlier works, our Green-Naghdi model allows a non-flat topography, and horizontal dimension d = 2. Its derivation follows directly from classical results concerning the one-layer case, and we believe such strategy may be used to construct interesting models in different regimes than the shallow-water/shallow-water studied in the present work.

oceanography, as it takes into account the dispersive effects neglected by the shallow-water (Saint-Venant) model and allows waves of greater amplitude than the Boussinesq model.
However, the previously mentioned works are restricted to the formal level, and the rigorous, mathematical justification of asymptotic models received a satisfactory answer only recently. We say that a model is fully justified (using the terminology of [37]) if the Cauchy problem for both the full Euler system and the asymptotic model is well-posed for a given class of initial data, and over the relevant time scale; and if the solutions with corresponding initial data remain close. The full justification of a system (S) follows from: • (Consistency) One proves that families of solutions to the asymptotic model, existing and controlled over the relevant time scale satisfies the full Euler system up to a small residual.
• (Existence) One proves that solutions of the full Euler system and solutions of the the model (S) with corresponding initial data do exist.
• (Convergence) One proves that the solutions of the full Euler system, and the ones of the asymptotic model, with corresponding initial data, remain close over the relevant time scale.
A result of Alvarez-Samaniego and Lannes [3] provides the existence and uniqueness of a solution to the full Euler system over the relevant time scale, uniformly with respect to the dimensionless parameters at stake, as well as a stability result with respect to perturbation of the equations. As a consequence, any consistent and well-posed asymptotic model is automatically justified in the sense described above. This strategy may be applied to the (quasilinear, hyperbolic) shallow-water model, to various Boussinesq models [33,5,6,11,43] and to Green-Naghdi models [3,29,31]. We let the reader refer to [37, Appendix C] for a reader's digest of the numerous known results on this aspect. All the aforementioned models coincide at the lower order of precision as a simple wave equation, which in dimension d = 1 predicts that any initial perturbation of the free surface will split up into two counter-propagating waves. An important family of models is dedicated to the precise study of the evolution of one of the two waves when higher order terms are included. The most famous example of such model is the Korteweg-de Vries [34] equation, but various extensions and generalizations have been proposed; see [32,15,30] and references therein. Again, the full justification of such models is recent: see [33,44,7,12,15,30].
All the aforementioned works are concerned with the case of a single layer of homogeneous fluid. Such assumption of may seem too crude for applications to oceanographic problems, as variations of salinity induce variation of density. In the present work, we are interested in the simplest setting that models such a variation of density: we consider a system of two layers of homogeneous, immiscible fluid, and we are interested in the evolution of the interface between the two layers. A considerable amount of interests has been given to such bi-fluidic systems; see [27] for a comprehensive review of the ins and outs on this topic.
The governing equations of the bi-fluidic systems share many aspects and properties of the aforementioned water-wave system, and its study has often been carried out in parallel. In particular, one can derive asymptotic models in analogy with the ones presented above. It is out of the scope of this introduction to present an exhaustive review of all the different models, as many settings are of interests. Here, and in the present work, we restrict ourselves to the case of a surface delimited by a flat, rigid lid (as deformation of the surface is in practice small compared to the deformation of the interface), and to the so-called shallow water/shallow water regime. In this regime, the two layers are assumed to be of comparable depth, and both small when compared to the typical horizontal wavelength of the flow. In that case, the models corresponding to the shallow water and Boussinesq systems have been derived in [13,18], and justified in the sense of consistency in [8] (where, incidentally, a much larger range of scaling regimes are studied). Green-Naghdi type models where obtained in the one-dimensional case in [39], and in the two-dimensional in [14]. An extensive study of scalar models has been provided by one of the authors in [22]. Let us note that all the aforementioned works are restricted to the case of a flat bottom, contrarily to the present work.
As attested earlier, bi-fluidic models have a similar structure as the ones in the one-layer case. As a matter of fact, one recovers the latter from the former when we assume that the mass density 1. The role of surface tension. Contrarily to the water wave case, the Cauchy problem for the full Euler system is ill-posed in Sobolev spaces in the absence of surface tension. However, surface tension is very small in practical cases, so that its effect is systematically negligible in all present asymptotic models. In [36], Lannes shows that a small amount of surface tension is sufficient to guarantee the well-posedness over times consistent with observations, provided that a stability estimate holds; see the somewhat more precise description in Section 2.2.
2. Absence of stability result. An equivalent result as the stability result (with respect to perturbations of the equation) for the full Euler system as described above is not known in the bi-fluidic case, partly due to the difficulties described in the previous item. In order to deal with this, a strategy consists in proving such stability result on the model itself. Thus the full justification of the model is a consequence of its well-posedness, and the full Euler system's consistency with the asymptotic model (and not the other way around). This strategy has been applied by the authors in [23], and we recall these results in Section 4.2.

3.
A non-local operator. As noticed in [8], shallow water models for internal waves with a rigid lid contain a non-local operator, which involves in particular the projection into the space of gradient functions; see Definition 3.11 in Section 3.3. This operator appears only in the case d = 2, and under the rigid-lid configuration (see [20] for shallow water/shallow water models with a free surface). The precise effect and meaning of this non-local term is yet to be fully understood.

4.
A critical-ratio. There exists a critical ratio for the depth of the two layers for which the first order (quadratic) nonlinearities vanish; δ 2 = γ in, e.g., (4.11). This phenomenon does not occur in the one-layer case, and motivates a precise study of unidirectional asymptotic models with stronger nonlinearities than in the classical long wave regime, and especially in the Camassa-Holm regime, for which first order dispersion and nonlinearities are formally of same magnitude in the critical case. This study has been carried out in [22], completed in [23], and the results are presented in Section 4.4.
In the present work, we report and complete recent mathematical results concerning the bifluidic system under the rigid lid assumption; from the well-posedness of the full Euler system to the justification of various models in the shallow water regime. Our aim is to provide a unified and comprehensive exposition of the existing theory. The above concerns and remarks appear spontaneously in the course of the study. We would like to mention the work of Saut, which pursues similar objectives than ourselves in [42].
Organization of the paper. The present paper is organized as follows. In Section 2, we introduce the non-dimensionalized full Euler equations describing the evolution of the two-fluid system with a rigid-lid we consider. We roughly describe in Theorem 2.2 its well-posedness result, obtained by Lannes in [36] (for a flat topography). Section 3 is dedicated to the construction and justification (in the sense of consistency) of the Green-Naghdi models. This result, in dimension d = 2, and allowing non-flat topography, is new to our knowledge. Several lower order models, for which stronger results have been recently obtained, are shown to descend directly from our Green-Naghdi system, and are described in Section 4. More precisely: • Section 4.1: The shallow water (Saint Venant) model, introduced in [8] and studied in details in [26]; • Section 4.2: A very recent Green-Naghdi model in the Camassa-Holm regime, introduced and rigorously justified in [23]; • Section 4.3: Boussinesq models, whose study follows from results in [7,21], and that we adapt to our case; • Section 4.4: Unidirectional (scalar) models generalizing the classical Korteweg-de Vries equation, whose rigorous justification has been investigated in [22].
Notations. In the following, C 0 denotes any nonnegative constant whose exact expression is of no importance. The notation a b means that a ≤ C 0 b.
We denote by C(λ 1 , λ 2 , . . . ) a nonnegative constant depending on the parameters λ 1 , λ 2 ,. . . and whose dependence on the λ j is always assumed to be nondecreasing. Let p be any constant with 1 ≤ p < ∞. We denote L p = L p (R d ) the space of all Lebesguemeasurable functions f with the standard norm The real inner product of any functions f 1 and f 2 in the Hilbert space L 2 (R d ) is denoted by For any real constant s ≥ 0, For convenience, we will make use of the following notation for given h 0 , s, t 0 ≥ 0:

For a vector-valued function
. The function spaces are endowed with canonical norms: For any functions u = u(t, X) and v(t, X) defined on [0, T ) × R d with T > 0, we denote the inner product, the L p -norm and especially the L 2 -norm, as well as the Sobolev norm, with respect to the spatial variable X, by u, v = u(t, ·), v(t, ·) , u L p = u(t, ·) L p , and |u| H s = |u(t, ·)| H s , respectively. We denote L ∞ ([0, T ); H s ) the space of functions such that u(t, ·) is controlled in H s , uniformly for t ∈ [0, T ): Finally, C k (R d ) denote the space of k-times continuously differentiable functions.
We conclude this section by the nomenclature that we use to describe the different regimes that appear in the present work. A regime is defined through restrictions on the set of admissible dimensionless parameters of the system, which are precisely defined in (2.2), below. Definition 1.1 (Regimes). We designate by shallow water regime the set of parameters with fixed 0 < µ max , δ min , δ max < ∞. We designate by Camassa-Holm regime (see [15]) the set and by long wave regime the set with some fixed M ∈ (0, ∞). Additionally, unless otherwise indicated, it is assumed that Bo −1 ∈ [0, Bo −1 min ] and β ∈ [0, 1]. The dependency of the constants on Bo −1 min is not displayed (Bo −1 min ≤ 1 in any oceanographic application).
2 The full Euler system 2.1 Construction of the full Euler system The system we study consists in two layers of immiscible, homogeneous, ideal, incompressible fluid under the only influence of gravity. The two layers are infinite in the horizontal dimension, and delimited above by a flat rigid lid and below by a non-necessarily flat bottom. The derivation of the governing equations of such a system is not new. We briefly recall it below, and refer to [8,4,22] for more details. We assume that the interface and bottom are given as the graph of a function (resp. ζ(t, X) and b(X)) which expresses the deviation from their rest position (resp. (X, 0),(X, −d 2 )) at the spatial coordinate X ∈ R d (d = 1 or d = 2) and at time t. Therefore, at each time t ≥ 0, the domains of the upper and lower fluid (denoted, respectively, Ω t 1 and Ω t 2 ), are given by We assume that the two domains are strictly connected, that is there exists h > 0 such that We denote by (ρ 1 , v 1 ) and (ρ 2 , v 2 ) the mass density and velocity fields of, respectively, the upper and the lower fluid. The two fluids are assumed to be homogeneous and incompressible, so that the mass densities ρ 1 , ρ 2 are constant, and the velocity fields v 1 , v 2 are divergence free.
Here and thereafter, we use the notation ∇ X,z to designate the gradient operator with respect to the variable (X, z), while the ∇ and ∆ simply denote the gradient and the Laplacian operators with respect to the variable X.
As we assume the flows to be irrotational, one can express the velocity field as gradients of a potential: v i = ∇ X,z φ i (i = 1, 2) , and the velocity potentials satisfy Laplace's equation The fluids being ideal, they satisfy the Euler equations. Integrating the momentum equations yields Bernoulli equations, written in terms of the velocity potentials: where P denotes the pressure inside the fluid. From the assumption that no fluid particle crosses the surface, the bottom or the interface, one deduces kinematic boundary conditions, and the set of equations is closed by the continuity of the stress tensor at the interface, which reads where k(ζ) = −∇ · 1 √ 1+|∇ζ| 2 ∇ζ denotes the mean curvature of the interface, and σ is the surface tension coefficient.
Altogether, the governing equations of our problem are given by the following where n denotes the unit upward normal vector at the surface at stake.
The next step consists in nondimensionalizing the system. Thanks to an appropriate scaling, the two-layer full Euler system (2.1) can be written in dimensionless form. The study of the linearized system (see [36] for example), which can be solved explicitly, leads to a well-adapted rescaling.
Let a (resp. a b ) be the maximum amplitude of the deformation of the interface. We denote by λ a characteristic horizontal length (that we assume to be identical in any of the directions if d = 2; see [37] for a treatment of the anisotropic case when d = 2), say the wavelength of the interface. Then the typical velocity of small propagating internal waves (or wave celerity) is given by Consequently, we introduce the dimensionless variables 1 the dimensionless unknowns as well as the following dimensionless parameters We conclude by remarking that the system can be reduced into two evolution equations coupling Zakharov's canonical variables [51,19], namely (withdrawing the tildes for the sake of readability) the deformation of the free interface from its rest position, ζ, and the trace of the dimensionless upper potential at the interface, ψ, defined as follows: ψ ≡ φ 1 (t, X, ζ(t, X)).
Indeed, φ 1 and φ 2 are uniquely deduced from (ζ, ψ) as solutions of the following Laplace's problems: More precisely, we define the so-called Dirichlet-Neumann operators.
where φ 1 and φ 2 are uniquely defined (up to a constant for φ 2 ) as the solutions in H 2 (R d ) of the Laplace's problems (2.3)-(2.4).
The existence and uniqueness of a solution to (2.3)-(2.4), and therefore the well-posedness of the Dirichlet-Neumann operators follow from classical arguments detailed, for example, in [37].
Using the above definition, and after straightforward computations, one can rewrite the nondimensionalized version of (2.1) as a simple system of two coupled evolution equations, namely where we denote We will refer to (2.5) as the full Euler system, and solutions of this system will be exact solutions of our problem.

A well-posedness theorem on the full Euler system
We mention here that Lannes [36] recently ensured that the Cauchy problem for (2.5) (with a flat bottom: β = 0) is locally well-posed in Sobolev spaces, with an existence time consistent with observations. Earlier results showed that the problem was ill-posed in the absence of surface tension, outside of the analytic framework. It was subsequently proved that taking into account the surface tension restores the local well-posedness of the equations, but with a very small existence time of the solution when the surface tension is small, which is the case in the oceanographic setting. It has to be noted that none of the asymptotic models presented in the following sections (and as a matter of fact, no asymptotic model known to us) share the same property, and that the surface tension term could be withdrawn from the equations (by setting Bo −1 = 0) without hurting their well-posedness. The reason for this apparent paradox is that the positive role of surface tension is to regularize Kelvin-Helmholtz instabilities that appear at high frequencies, while the main part of the wave, which is captured by the asymptotic models, is located at low frequencies and is thus unaffected by surface tension.
In [36], Lannes introduces a stability criterion, whose role is to ensure that the aforementioned frequency threshold is high enough, and shows that, under this condition, the combined effect of surface tension and gravity is sufficient to control the regularity of the flow.
Somewhat more precisely, one has the following result; see [36,Theorems 5 and 6] for the precise statements.
Theorem 2.2. Let p ∈ P SW and the initial data U 0 ≡ (ζ 0 , ψ 0 ) satisfy the following: 1. U 0 belongs to an energy space of sufficiently smooth, bounded functions (in particular, the following is required: 2. ζ 0 satisfies the non-vanishing depth condition: 3. A stability criterion is satisfied, which can be roughly expressed by Then there exists a unique solution to (2.5) (with flat bottom: β = 0) with initial data U | t=0 ≡ U 0 , bounded in the same energy space (no loss of derivatives). The flow is continuous with respect to time, and defined for t ∈ [0, − T ], where T > 0 depends only on the quantities defined through the three above conditions, and in particular can be chosen independent of the parameters p ∈ P SW .

The Green-Naghdi/Green-Naghdi model
In the following, we construct Green-Naghdi type models for the system (2.5), that is asymptotic models with precision O(µ 2 ), in the sense of consistency. As we shall see, and contrarily to earlier works, our construction relies only on asymptotic expansions which can be straightforwardly deduced from known results on the one-layer case. Thus we start by recalling below these results, which can be found in particular in [37]. We then deduce equivalent asymptotic expansions in the bi-fluidic setting in Section 3.2, and finally use these expansions to construct our asymptotic models in Section 3.3.

Asymptotic expansions in the water-wave case
The proof of the following statements may be found in [37] (with depth D = 1, but the general case is obtained by straightforward change of variables). For simplicity, and without lack of generality, we set = β = 1 in this section.
Let us now recall that the Dirichlet-Neumann operator may be equivalently defined through the vertically averaged mean velocity, thanks to the following Proposition.
Then one has the identity Proof. This striking result is a consequence of a simple calculation, that we recall. Let ϕ ∈ C ∞ c (R d ) be a test function. Then one has where we used Green's identity, and the Laplace's equation satisfied by φ. Since this result is valid for any test function 2) holds in the strong sense.
Let us conclude with the asymptotic expansion of the quantities defined above. Here and in the following, we denote, for convenience, 2

Asymptotic expansions in the bi-fluidic case
Our specific operators may be deduced from the classical Dirichlet-Neumann operator used in the water-wave problem, and Defined in Definition 3.1. Thus the following results are easily deduced from the ones of the previous section. Let us first define u 1 (resp. u 2 ) the vertically averaged mean velocity of the upper layer (resp. lower layer): where φ 1 and φ 2 are uniquely defined (up to a constant for φ 2 ) as the solutions in H 2 (R d ) of the Laplace's problems (2.3)-(2.4).
Proposition 3.5. Let ζ, b, ψ satisfy the hypothesis of Definition 3.4. Then one has the identity The proof of these identities is identical as the one of Proposition 3.2 (when considering the upper and lower potential respectively, and using that Thus we omit the proof, and continue with the asymptotic expansions of the above quantities. Proof. Expansions (3.9),(3.10), simply follow from Proposition 3.
It follows that one has the identity , and the expansions (3.9),(3.10) follow. Expansion (3.11) is a straightforward consequence of (3.9),(3.10), and Proof. As above, the expansions can be deduced from Proposition 3.3, once we remark that by definition, φ 2 (X, z) satisfies In other words, one has the identity Thus one deduces from Proposition 3.3 the following estimates: 14) Furthermore, one has from [8, Proposition 3] that so that estimate (3.12) is now straightforward. Finally, estimate (3.13) is easily deduced from the previous estimates. 3 Note that by definition of the Dirichlet-Neumann operators G µ and G µ,δ −1 , this identity yields In other words, and as remarked in [36], one has the identity In particular, the bound (3.16) is not optimal; see [36, Proposition 1 and Remark 6].

Construction of the Green-Naghdi/Green-Naghdi model
Let us recall the full Euler system (2.5): By Proposition 3.5, the first equation of (3.17) yields When we plug the expansions of Propositions 3.6 and 3.7 into the second equation of (3.17), and withdrawing O(µ 2 ) terms, one obtains Remark 3.8. Equations (3.18) and (3.19) are very similar to the system obtained in [14]. It may also be recovered from system (60) in [20] when setting α = 0 (notation therein), and after straightforward adjustments (in particular, we use a different scaling in the non-dimensionalizing step). (3.19), the latter up to a remainder term, R, bounded by Proof. The fact that (ζ p , u p 1 , u p 2 ) satisfy (3.18) has been expressed earlier in Proposition 3.5. That (3.19) approximately holds is a consequence of the asymptotic expansions of Propositions 3.6 and 3.7. Let us detail briefly the argument.
Subtracting (3.19) to the last equation in (2.5) yields (withdrawing the explicit reference to p ∈ P SW for the sake of readability) We now show how to estimate each of these terms.
Recall (3.13) in Proposition 3.7. It follows from tame product estimates in H s (see [23, Appendix A] for example) Here and below, we denote t 0 > d/2. Identically, using (3.11) in Proposition 3.6, one obtains It is now clear that one can choose N sufficiently large so that the following holds: with C as in the Proposition (note that one can deduce a control in H s of u 1 from a control in H s+2 of ∇ψ, thanks to (3.9) -being non optimal).
The estimate on R III is obtained similarly. Using identity (3.8) as well as first order expansions (3.9),(3.14), one obtains so that one deduces from the above estimates that with C as in the Proposition, and for N sufficiently large.
The estimate on R I requires a control of the time derivatives. One can obtain equivalent results as in Propositions 3.6 an 3.7, and in particular for s ≥ t 0 + 1; following the same method, but after differentiating the equations (with respect to the time variable, t). We do not detail the proof, and refer to [20,37] for examples of applications of this strategy.
Finally, note that one can control ∂ t ∇ψ and ∂ t ζ using only a control on ∇ψ and ζ, using that (ζ, ψ) satisfies the full Euler system (2.5); allowing for a loss of derivatives. At the end of the day, one sees that if N is sufficiently large, then one has with C as in the Proposition.
Our aim is now to approximate (3.18),(3.19) as a system of coupled evolution equations (thus directly comparable with (3.17)). In order to do so, we introduce a new velocity variable, v, which shall satisfy In dimension d = 1, identity (3.24) (assuming that v → 0 as |x| → ∞) uniquely defines v as the shear mean velocity: However, such an explicit expression is not available in dimension d = 2. In that case, we make use of the following operator, defined in [8].
and ΠV = V , where Π = ∇∇ |D| 2 is the orthogonal projector onto the gradient vector fields of This allows to define v as the unique gradient solution to (3.24).
Definition 3.11. Let ζ ∈ L ∞ (R d ) be such that |ζ| L ∞ < 1 and | ζ − βb| L ∞ < δ −1 , so that Then we define v as the unique gradient solution to (3.24); or, in other words, Note that the condition |ξ| L ∞ < 1 is ensured by the following: Remark 3.12. In the one dimensional space (d = 1), one has Π = Id, and one can check that the operator Q simply reduces to Q[ξ]W = 1 1+ξ W , so that one recovers v = u 2 − γu 1 and Note also that in that case, conditions |ζ| L ∞ < 1 and | ζ − βb| L ∞ < δ −1 may be replaced by the slightly less stringent non-vanishing depth condition: Let us emphasize again that in the case d = 2, there is no reason to think that v = u 2 −γu 1 holds, even with precision O(µ); and in particular that u 2 − γu 1 is a gradient vector fields. In the same way, one would like to write, seeing (3.24), u 2 = δQ[ δζ] h1h2 h1+γh2 v and u 1 = −Q[− ζ] h1h2 h1+γh2 v ; but unfortunately, it is not clear that u 1 , u 2 are gradient vector fields (as a matter of fact, their second order expansion tends to show that it is not true). However, one has the following expansion: Then one has h1+γh2 v . Proof. The first estimate follows from where we used identity (3.24). Consequently, Lemma 3.10 yields The control of ∇ψ − u 1 H s as in estimate (3.27) is now a consequence of Proposition 3.6, and the continuity of the operator Q expressed in Lemma 3.10. The control of ∇ψ − u 1 H s is then deduced, using once again Proposition 3.6 and triangular inequality. Estimate (3.27) is proved. Estimate (3.28) is obtained in the same way, but using the control of u 2 − ∇H µ,δ ψ H s displayed in Proposition 3.7.
Following the same strategy, one order further, yields The last term in the identity above is estimated in part thanks to Proposition 3.6, and in part thanks to the first order estimate (3.27). Estimate (3.29) then follows as above from the definition and continuity of the the operator Q; see Lemma 3.10. Estimate (3.30) is obtained in the same way, and we omit the detailed calculations.
One deduces from Proposition 3.13 the following approximate equation equivalent to (3.19) (withdrawing again O(µ 2 ) terms): Proposition 3.14 (Consistency). Let U p ≡ (ζ p , ψ p ) be a family of solutions to the full Euler system (2.5) such that there exists T > 0, s ≥ 0 for which (ζ p , ∇ψ p ) is bounded in L ∞ ([0, T ); H s+N ) d+1 (N sufficiently large), uniformly with respect to (µ, , δ, γ) ≡ p ∈ P SW ; see Definition 1.1. Moreover, assume that b ∈ H s+N and Define v p through Definitions 3.4 and 3.11. Then (ζ p , v p ) satisfy (3.25) and (3.31), the latter up to a remainder term, R, bounded by The Proposition is obtained as in the proof of Proposition 3.9, but using the asymptotic expansions of Proposition 3.13; we omit its proof.
Unidimensional case (d = 1). Recall that in the one dimensional space, one has simply for any V ∈ L 2 (R), and denoting h 1 ≡ 1 − ζ and h 2 ≡ δ −1 + ζ. In particular, one can check that u 1 = u 1 = −h2 v h1+γh2 and u 2 = u 2 = h1 v h1+γh2 . The system (3.25),(3.31) thus becomes where we denote If, additionally, one assumes the bottom is flat (by setting β = 0), then one recovers the following system, as already introduced in [22]: where we define: Proposition 3.14 thus generalizes the consistency result obtained in [22,23] to the case d = 2, and to non-flat topography.

Lower order models
The system of equations (3.25),(3.31) is very broad, in the sense that it has been obtained with minimal assumptions: allowing d = 1 and d = 2, non-flat topography, and in the shallow water regime of Definition 1.1. It is justified by a consistency result (Proposition 3.14). As argued in the introduction, the consistency result alone is not sufficient to fully justify a model. In particular, its well-posedness should be confirmed. As a matter of fact, contrarily to the water-wave case [3], the well-posedness of the Green-Naghi model in the bi-fluidic case is not clear, and similar systems have been proved to be ill-posed; see [38] and discussion in [16]. In the following subsections, we show that existing models in the literature directly descend from our Green-Naghdi model (3.25),(3.31), after additional assumptions (typically, restricting to d = 1, flat bottom, and/or more stringent regimes), or with a lower precision. Their justification in the sense of consistency is therefore a direct application of Proposition 3.14, and stronger results (well-posedness, convergence) are stated when available.

The Shallow water (Saint Venant) model
In this section, we consider only the first order terms in equation (3.31) (equivalently, we set µ = 0; this corresponds to the assumption that the horizontal velocity is constant across the vertical layers). The system (3.25),(3.31), withdrawing O(µ) terms, is now simply where we recall that h1+γh2 v , where the operator Q is defined in Definition 3.10.
A similar system is obtained when using a different velocity variable, such as the shear velocity at the interface V ≡ ∇H µ,δ ψ − γ∇ψ. In that case, one obtains the following [8]: where R is defined similarly as Q: R[ ζ]W is the only gradient solution to System (4.1) and system (4.2) are equivalent up to order O(µ). In fact from proposition 3.13, one has V = u 2 − γ u 1 + O(µ), and In one dimension (i.e. d = 1), both (4.1) and (4.2) read System (4.2) has been derived and justified in the sense of consistency in [8], in the case of a flat bottom, and without surface tension. An equivalent consistency result clearly holds for (4.1) (as a consequence of Proposition 3.14 in particular). More precisely, one has Define v as in definition 3.11, and V ≡ ∇H µ,δ ψ − γ∇ψ. Then (ζ, v) satisfies (4.1), up to a remainder (0, R 1 ) ; and (ζ, V ) satisfies (4.2), up to a remainder (0, R 2 ) ; with , uniform with respect to the parameters p ∈ P SW .

A Green-Naghdi type model in the Camassa-Holm regime
The present section is limited to the so-called Camassa-Holm regime (that is using additional assumptions = O( √ µ); see Definition 1.1), flat bottom case (β = 0) and one dimensional space d = 1. Moreover, we neglect the surface tension component in our model. All the estimates are now understood uniformly with respect to (µ, , δ, γ) ∈ P CH .
In that case, one can easily check that the following approximations hold for Q, R defined in section 3.3 Plugging these expansions into system (3.33) yields a simplified model, precise with the same order of magnitude as the original model (that is O(µ 2 )) in the Camassa-Holm regime. Furthermore, after several additional transformations, ones may produce an equivalent model (again, in the sense of consistency) which possesses a structure similar to symmetrizable quasilinear systems, thus allowing its full justification. The following system of equations has been introduced and justified by the authors in [23].
with q i (X) ≡ 1 + κ i X (i = 1, 2) and κ 1 , κ 2 , ς are defined by : We recall that the shear mean velocity is uniquely defined in the case of one dimension, d = 1 (see (3.26) and (3.8)) by System (4.6) is fully justified as an asymptotic model by the following results (see [23] for a more precise statement): Proposition 4.4 (Consistency). Let U p ≡ (ζ p , ψ p ) be a family of solutions of the full Euler system (2.5) for which (ζ p , ∂ x ψ p ) is bounded in L ∞ ([0, T ); H s+N ) 2 with s ≥ 0 and sufficiently large N , and uniformly with respect to (µ, , δ, γ) ≡ p ∈ P CH ; see Definition 1.1. Moreover, assume Define v p as in (4.10). Then (ζ p , v p ) satisfies (4.6), up to a remainder (0, R) , bounded by , uniform with respect to the parameters p ∈ P CH . System (4.6) is well-posed (in the sense of Hadamard) in the energy space X s = H s (R)×H s+1 (R), endowed with the norm provided that the following ellipticity condition (for the operator T) holds: Theorem 4.5 (Existence and uniqueness). Let s 0 > 1/2, s ≥ s 0 + 1, and let U 0 = (ζ 0 , v 0 ) ∈ X s satisfy (H1),(H2). Then there exists a maximal time T max > 0, uniformly bounded from below with respect to p ∈ P CH , such that the system of equations (4.6) admits a unique solution and preserving the conditions (H1),(H2) (with different lower bounds) for any t ∈ [0, T max / ). Moreover, for any 0 ≤ T < T max , there exists C 0 , λ T = C(h −1 01 , h −1 02 , δ max , δ −1 min , M, T, U 0 X s ), independent of p ∈ P CH , such that one has the energy estimate If T max < ∞, one has or one of the two conditions (H1),(H2) ceases to be true as t −→ T max / .
Proposition 4.6 (Stability). Let s ≥ s 0 + 1 with s 0 > 1/2, and let U 0,1 = (ζ 0,1 , v 0,1 ) ∈ X s and U 0,2 = (ζ 0,2 , v 0,2 ) ∈ X s+1 . Under the assumptions of Theorem 4.5, let U j be the solution of system (4.6) with U j | t=0 = U 0,j .Then there exists T, λ, Finally, the following "convergence result" states that the solutions of our system approach the solutions of the full Euler system, with as good a precision as µ (and Bo −1 ) is small. The above results hold on time interval t ∈ [0, T / ] with T bounded by below, independently of p ∈ P CH , provided that a stronger criterion is satisfied by the initial data. This corresponds to setting = 1 in Theorem 2.2; see criterion (5.5) and Theorem 6 in [36] for the precise statement.
Remark 4.8. The new model allows to fully justify any well-posed system, consistent with our model (4.6), thanks to an a priori estimate between two approximate solutions of our system, established in [23] (see Proposition 7.1 and Theorem 7.5 therein). This is used in particular to obtain the convergence results of the unidirectional and decoupled approximations stated in Section 4.4.
Remark 4.9 (The model with surface tension). The previous results concerning (4.6), where the surface tension effects are neglected, still hold when surface tension is taken into account; see (2.5). These results will be precisely stated in a forthcoming paper.

Boussinesq models
In this section, we restrict ourselves to the case of flat bottom (β = 0), and unidimensional case (d = 1). Moreover, we restrict the regime under study to the so-called long wave regime, where = O(µ); see Definition 1.1. We also assume that the surface tension term is at most of size Bo −1 = O(µ), with the following: In that case, when withdrawing O(µ 2 ) terms, one can easily check (see the proof of Proposition 4.10, below for detailed calculations) that the system (3.33) becomes a simple quasilinear system, with additional linear dispersive terms: The full Euler system is consistent with this model, with the same precision as the Green-Naghdi models, provided that the assumptions of the long wave regime, and in particular ≤ M µ, hold (in addition to the aforementioned β = 0, d = 1); see Proposition 4.10, below.
Boussinesq systems with improved frequency dispersion. Let us emphasize that system (4.11) is only one of a large family of Boussinesq-type models, that are consistent with precision O(µ 2 ). Briefly, one can make use of • Near-identity change of variable. Define, for θ 1 , θ 2 ≥ 0, the following When rewriting (4.11) with respect to this new variable, v θ1,θ2 , and withdrawing O(µ 2 ) terms, one obtains with U θ1,θ2 ≡ (ζ, v θ1,θ2 ) and • The Benjamin-Bona-Mahony trick. Using that ∂ t U θ1,θ2 + A 0 ∂ x U θ1,θ2 = O(µ), one has for any λ 1 , λ 2 ∈ R: Plugging this approximation into (4.12) yields, withdrawing again O(µ 2 ) terms, The above transformations are useful, for example, with the aim of improving the frequency dispersion, that is choosing the coefficients so that the dispersion relation fits the one of the full Euler system with high precision, even for relatively large µ (using truncated Taylor series or Padé approximant). It may also be useful for mathematical purposes to generate such a large family of models, which are all equivalent in the sense of consistency, but may have very different properties (well-posedness, integrability, etc.). We let the reader refer to [37] for a more detailed account.
In [43], Saut and Xu offer an in-depth study of the well-posedness of such Boussinesq systems in the water-wave setting (thus with different coefficients). It would be interesting, but out of the scope of the present work, to adapt the techniques developed therein to our bi-fluidic systems (4.13).
A fully justified symmetric Boussinesq model. Another strategy for constructing a model with improved properties, that has been used in [7,21] and that we develop in the following, consists in symmetrizing the original model (4.11), up to precision O(µ 2 ). Define Multiplying (4.11) with (S 0 + S[U ] − µT ∂ 2 x ) and withdrawing O(µ 2 ) terms yields with the following symmetric matrices and operators: This new model is fully justified, in the sense described in the introduction. Let us detail the consistency, well-posedness and convergence results below. Proposition 4.10 (Consistency). Let U p ≡ (ζ p , ψ p ) be a family of solutions to the full Euler system (2.5) with β = 0 and d = 1, such that there exists T > 0, s ≥ 0 for which (ζ p , ∂ x ψ p ) is bounded in L ∞ ([0, T ); H s+N ) 2 (N sufficiently large), uniformly with respect to (µ, , δ, γ) ≡ p ∈ P LW ; see Definition 1.1, and bo −1 ≤ bo −1 min , bo min > 0. Moreover, assume that Define v p ≡ u 2 − γu 1 ; see Definition 3.4. Then (ζ p , v p ) satisfy (4.11) (resp. (4.14)) up to a remainder term, R B (resp. R S ), bounded by Proof. Let us first recall that since we assume d = 1, then defining v p ≡ u 2 − γu 1 is equivalent as defining v p through Definition 3.11, but requires only the non-vanishing depth condition h 1 , h 2 ≥ h 0 > 0, instead of the more stringent condition of Definition 3.11. Let us also emphasize that system (3.25), (3.31) is exactly (that is, without any approximation) system (3.33), in the case β = 0 and d = 1. Thus, by Proposition 3.14, we know that (ζ p , v p ) satisfies (3.33) up to a remainder term, R, satisfying with C as in the Proposition 4.10. Thus one only needs to check that the neglected terms from (3.33) in (4.11) (resp. (4.14)), when using that ≤ M µ, are estimated in the same way. Let us detail briefly.
We first claim that the following holds, for any s ≥ 0 and with t 0 > 1/2: The formal expansion (as powers of ζ) is easily checked, so that the estimate in L ∞ norm is straightforward: The control in H s is slightly more elaborate, as h 1 , h 2 are not bounded in H s , since they do not decrease at infinity. We let the reader refer to [23,Lemma 4.5] to see how this technical difficulty may be faced. It follows from (4.15) that ). One obtains in the same way the following estimates: ). Let us recall that estimates on ∂ t ζ and ∂ t v may be deduced from estimates on ζ, v as they satisfy the full Euler system (2.5), allowing for a loss of derivatives.
Altogether, choosing N sufficiently large and since ≤ M µ when p ∈ P LW , it follows that (ζ, v) satisfy (4.11) up to a remainder term, denoted R B , and bounded by One obtains similarly the estimate concerning the symmetric system (4.11), after controlling the extra error terms: This is easily checked if N is sufficiently large, and using again that ≤ M µ. Proposition 4.10 is proved.
In addition to its justification in the sense of consistency, one is able to fully justify the symmetric Boussinesq model (4.14), thanks to its following properties: 1. The matrices S 0 , Σ 0 , S 1 , Σ 1 are 2-by-2 symmetric matrices; 2. S[·] and Σ[·] are linear mappings with values into 2-by-2 symmetric matrices; 3. S 0 and S 2 are definite positive.
These properties allow to control a natural energy of the system: which is equivalent to the scaled Sobolev norm X s+1 µ : provided that U is sufficiently small, (for the nonlinear terms to be controlled).
The following Propositions are a direct consequence from [21, Proposition 2.6 and 2.8] (where the author deal with the case of internal waves with a free surface, so that the system has four equations, but possess the exact same structure), and we do not detail the proof further on. Moreover, one has the following estimate for t ∈ [0, T / ]: A stability result similar to Proposition 4.7 holds, so that the symmetric Boussinesq system is well-posed in the sense of Hadamard.
Theorem 4.12 (Convergence). Let p ∈ P LW and U 0 ≡ (ζ 0 , ψ 0 ) ∈ H s+N (R) 2 , with N sufficiently large, satisfy the hypotheses of Theorem 5 in [36]; see Theorem 2.2. Then there exists C, T, 0 > 0, independent of p ∈ P LW and Bo −1 ≤ µ bo −1 min , such that for any 0 ≤ ≤ 0 , • There exists a unique solution U ≡ (ζ, ψ) to the full Euler system (2.5), with β = 0 and d = 1, defined on [0, T ] and with initial data (ζ 0 , ψ 0 ) (provided by Theorem 5 in [36]); • The non-vanishing depth condition is satisfied for ≤ 0 , so that one can define v ≡ u 2 − γu 1 , with u 1 , u 2 as in Definition 3.4; • There exists a unique solution U B ≡ (ζ B , v B ) to the symmetric Boussinesq model (4.14), defined on [0, T ] and with initial data (ζ 0 , v | t=0 ) (provided by Theorem 4.11); • The difference between the two solutions is controlled as The above results hold on time interval t ∈ [0, T / ] with T bounded by below, independently of p ∈ P LW , provided that a stronger criterion is satisfied by the initial data. This corresponds to setting = 1 in Theorem 2.2; see criterion (5.5) and Theorem 6 in [36] for the precise statement.

Scalar models
In this section, we are interested in the justification of scalar asymptotic models for the propagation of internal waves (as opposed to all aforementioned models, which consist in a system of evolution equations). The derivation and study of such models have a very rich and ancient history, starting with the work of Boussinesq [9] and Korteweg-de Vries [34] which introduced the famous Korteweg-de Vries equation ∂ t u + c∂ x u + αu∂ x u + ν∂ 3 x u = 0, as a models for the propagation of surface gravity waves, in the long wave regime. However, the complete rigorous justification of such model is much more recent [33,44,7] (see [21] for the bi-fluidic case).
One obvious discrepancy between scalar models such as the KdV equation and coupled models, is that the former selects a direction of propagation (right or left, depending on the sign of c). On the contrary, the first order approximation of coupled models is a linear wave equation, which predicts that any initial perturbation of the flow will split into two counter-propagating waves. This yields two very different possible justifications of scalar models: • Unidirectional approximation. One proves that if the initial perturbation (that is the deformation of the interface as well as shear layer-mean velocity) is carefully chosen, then the flow can be approximated as a solution of a scalar equation. Physically speaking, this means that we focus our attention on only one of the two counter-propagating waves after they have split.
• Decouled approximation. For a generic initial perturbation of the flow, we approximate the flow as the superposition of two (uncoupled) waves, each one driven by a scalar equation. This justification is of course more general as for the admissible initial data, but its precision is often much worse, as controlling the coupling effects between the two counter-propagating waves (and in particular their secular growth), which are neglected in the decoupled approximation, is arduous.
A major difference between the water-wave case and the bi-fluidic case is that in the latter, there exists a critical ratio (δ 2 = γ) for which the first order (quadratic) nonlinearity of our models vanishes; see the Boussinesq system (4.11) for example. This allows to consider a regime with greater nonlinearities than the long wave regime (and in particular the Camassa-Holm regime; see Definition 1.1), and thus motivates higher order models than the Korteweg-de Vries equations, in order to recover a O(µ 2 ) precision. Here, we focus on the so-called Constantin-Lannes equation: (4.17) where β, α i (i = 1, 2, 3), ν, κ 1 , κ 2 , are fixed, given parameters . This equation has been studied and justified as a model for the propagation of unidirectional surface gravity waves in [15]. In particular, the well-posedness of the Cauchy problem for (4.17) in Sobolev spaces is proved, provided β > 0. The justification of both the unidirectional and decoupled approximation in the bi-fluidic case, in the sense of consistency, have been worked out in [22]. In what follows, we restrict to the Camassa-Holm regime and to the Constantin-Lannes equations for the sake of simplicity, but more general results may be found therein. The full justification, and in particular the convergence results stated below is then a consequence of the properties of the Green-Naghdi type model introduced by the authors in [23], and recalled in section 4.2. We let the reader refer to the aforementioned references for more details, and disclose the statements below without further comments.
If the initial data satisfies a stronger stability hypothesis (see [36,Theorem 6]), then the result holds for large time t ∈ [0, T / ], with T independent of p.
• Moreover, assume that the initial data is sufficiently localized in space, that is more precisely (1 + | · | 2 )∂ k ζ 0 and (1 + | · | 2 )∂ k v | t=0 ∈ H s (R), k ∈ {0, . . . , 7}, then one has the improved estimate The first three items hold for large time t ∈ [0, T / ], with T independent of p, provided that a stronger criterion is satisfied by the initial data (see hypotheses in [36,Theorem 6]). In that case, the last item is valid for time t ∈ [0, T / max{ , µ}].