A mathematical analysis of the Kakinuma model for interfacial gravity waves. Part I: Structures and well-posedness

We consider a model, which we named the Kakinuma model, for interfacial gravity waves. As is well-known, the full model for interfacial gravity waves has a variational structure whose Lagrangian is an extension of Luke's Lagrangian for surface gravity waves, that is, water waves. The Kakinuma model is a system of Euler-Lagrange equations for approximate Lagrangians, which are obtained by approximating the velocity potentials in the Lagrangian for the full model. In this paper, we first analyze the linear dispersion relation for the Kakinuma model and show that the dispersion curves highly fit that of the full model in the shallow water regime. We then analyze the linearized equations around constant states and derive a stability condition, which is satisfied for small initial data when the denser water is below the lighter water. We show that the initial value problem is in fact well-posed locally in time in Sobolev spaces under the stability condition, the non-cavitation assumption and intrinsic compatibility conditions in spite of the fact that the initial value problem for the full model does not have any stability domain so that its initial value problem is ill-posed in Sobolev spaces. Moreover, it is shown that the Kakinuma model enjoys a Hamiltonian structure and has conservative quantities: mass, total energy, and in the case of the flat bottom, momentum.


Introduction
We are concerned with the motion of interfacial gravity waves at the interface between two layers of immiscible waters in a domain of the (n + 1)-dimensional Euclidean space in the rigid-lid case. Let t be the time, x = (x 1 , . . . , x n ) the horizontal spatial coordinates, and z the vertical spatial coordinate. We assume that the interface, the rigid-lid of the upper layer, and the bottom of the lower layer are represented as z = ζ(x, t), z = h 1 , and z = −h 2 + b(x), respectively, where ζ(x, t) is the elevation of the interface, h 1 and h 2 are mean thicknesses of the upper and lower layers, and b(x) represents the bottom topography. The only external force applied to the system is the constant and vertical gravity, and interfacial tension is neglected. Moreover, we assume that the waters in the upper and the lower layers are both incompressible and inviscid fluids with constant densities ρ 1 and ρ 2 , respectively, and that the flows are both irrotational. See Figure 1.1. Then, the motion of the waters is described by the velocity potentials Φ 1 and Φ 2 and the pressures P 1 and P 2 in the upper and the lower layers, respectively, satisfying the basic equations in the theory of fluid dynamics, which will be referred as the full model for interfacial gravity waves throughout in this paper. As shown by J. C. Luke [23], the basic equations for the surface gravity waves, that is, the water wave problem has a variational structure, whose Lagrangian is written in terms of the surface elevation of the water and the velocity potential, and the Lagrangian density is given by the vertical integral of the pressure in the water region. The full model for interfacial gravity waves has also a variational structure and the Lagrangian density L (Φ 1 , Φ 2 , ζ) is again given by the vertical integral of the pressure in both water regions. T. Kakinuma [17,18,19] proposed a model for interfacial gravity waves and applied his model for k = 1, 2, where {Z 1,i } and {Z 2,i } are appropriate function systems in the vertical coordinate z and may depend onh 1 (x) andh 2 (x), respectively, which are thickness of the upper and the lower layers in the rest state, whereas φ k = (φ k,0 , φ k,1 , . . . , φ k,N ) T , k = 1, 2, are unknown variables. Then, he derived an approximate Lagrangian density L app (φ 1 , φ 2 , ζ) = L (Φ app 1 , Φ app 2 , ζ) for unknowns (φ 1 , φ 2 , ζ). The Kakinuma model is a corresponding system of Euler-Lagrange equations for the approximated Lagrangian density L app (φ 1 , φ 2 , ζ). Different choices of the function systems {Z 1,i } and {Z 2,i } give different Kakinuma models and we have to carefully choose the function systems for the Kakinuma model to provide good approximations for interfacial gravity waves.
The Kakinuma model is an extension to interfacial gravity waves of the so-called Isobe-Kakinuma model for the surface gravity waves, that is, the water waves. In the case of the surface gravity waves, the basic equations are known to have a variational structure with Luke's Lagrangian density L Luke (Φ, ζ), where ζ is the surface elevation and Φ is the velocity potential of the water. The Isobe-Kakinuma model is a system of Euler-Lagrange equations for the approximated Lagrangian density L app (φ, ζ) = L Luke (Φ app , ζ), where Φ app is an approximate velocity potential and φ = (φ 0 , φ 1 , . . . , φ N ) T are unknown variables. The model was first proposed by M. Isobe [15,16] and then applied by T. Kakinuma to simulate numerically the water waves. We note that a similar model was derived by G. Klopman, B. van Groesen, and M. W. Dingemans [21], and used to simulate the water waves. See also Ch. E. Papoutsellis and G. A. Athanassoulis [29]. Recently, this model was analyzed from mathematical point of view. One possible choice of the function system {Z i } is a set of polynomials in z, for example, Z i (z; b(x)) = (z + h − b(x)) p i with integers p i satisfying 0 = p 0 < p 1 < · · · < p N . Under this choice of the function system {Z i }, the initial value problem to the Isobe-Kakinuma model was analyzed by Y. Murakami and T. Iguchi [27] in a special case and by R. Nemoto and T. Iguchi [28] in the general case. The hypersurface t = 0 in the space-time R n × R is characteristic for the Isobe-Kakinuma model, so that one needs to impose some compatibility conditions on the initial data for the existence of the solution. Under these compatibility conditions and a sign condition −∂ z P app ≥ c 0 > 0 on the water surface, they showed the well-posedness of the initial value problem locally in time, where P app is an approximate pressure in the Isobe-Kakinuma model calculated from Bernoulli's equation. Moreover, T. Iguchi [12,13] showed that under the choice of the function system in the case of the flat bottom, (z + h − b(x)) i in the case of a variable bottom, the Isobe-Kakinuma model is a higher order shallow water approximation for the water wave problem in the strongly nonlinear regime. Furthermore, V. Duchêne and T. Iguchi [8] showed that the Isobe-Kakinuma model also enjoys a Hamiltonian structure analogous to the one exhibited by V. E. Zakharov [32] on the full water wave problem. Our aim in the present paper is to extend these results on the surface gravity waves to interfacial gravity waves.
In view of these results on the Isobe-Kakinuma model, in the present paper we consider the Kakinuma model under the choice of the approximate velocity potentials in (1.1) as where N, N * and p 0 , p 1 , . . . , p N * are nonnegative integers satisfying 0 = p 0 < p 1 < · · · < p N * . In applications of the Kakinuma model, it would be better to choose N * = N and p i = 2i in the case of the flat bottom, and N * = 2N and p i = i in the case of a variable bottom. In the case N = N * = 0, that is, if we choose the approximation Φ app k (x, z, t) = φ k (x, t) for k = 1, 2 the functions independent of the vertical coordinate z, then the corresponding Kakinuma model is reduced to the shallow water equations. In the case N + N * > 0, the Kakinuma model is classified into a system of nonlinear dispersive equations.
As a shallow water limit h 1 |ξ|, h 2 |ξ| → 0, we have where c SW is the phase speed of infinitely long and small interfacial gravity waves. In Section 3, we will analyze the linear dispersion relation of the Kakinuma model and calculate the phase speed c K (ξ) of the plane wave solution related to the wave vector ξ. Under the choice N * = N and p i = 2i, or N * = 2N and p i = i in the approximation (1.4) of the velocity potentials, it turns out that which indicates that the Kakinuma model may be a good approximation of the full model for interfacial gravity waves in the shallow water regime h 1 |ξ|, h 2 |ξ| ≪ 1. We note that Miyata-Choi-Camassa model derived by M. Miyata [26] and W. Choi and R. Camassa [4] is a model for interfacial gravity waves in the strongly nonlinear regime and can be regarded as a generalization of the Green-Naghdi equations for water waves into a two-layer system. Let c MCC (ξ) be the phase speed of the plane wave solution related to the wave vector ξ for the linearized equations of the Miyata-Choi-Camassa model around the rest state. Then, we have so that the Kakinuma model gives a better approximation of the full model than the Miyata-Choi-Camassa model in the shallow water regime, at least, at the linear level. A rigorous analysis for the consistency of the Kakinuma model in the shallow water regime will be analyzed in the subsequent paper V. Duchêne and T. Iguchi [9]. On the other hand, in the deep water limit we have lim which is not consistent with the limit of the full model We notice that the Miyata-Choi-Camassa model is only apparently consistent with the full model in this deep water limit since lim h 1 |ξ|,h 2 |ξ|→∞ c MCC (ξ) 2 = 0, but we note also that We refer to V. Duchêne, S. Israwi, and R. Talhouk [10] for further discussion and the derivation of modified Miyata-Choi-Camassa models having either the same dispersion relation as the full model, or the same behavior as the Kakinuma model in the deep water limit. As we discuss below, thanks to the high-frequency behavior of the linearized equations, and contrarily to both the full model and the Miyata-Choi-Camassa model, the Kakinuma model has a non-trivial stability domain and, as a result, the initial value problem to the Kakinuma model is well-posed locally in time in Sobolev spaces under appropriate assumptions on the initial data.
As we have already seen, the roots ω ∈ C of the dispersion relation of the linearized equations of the full model around the rest state are always real, so that the corresponding initial value problem is well-posed. However, as for the nonlinear problem, even if the initial velocity is continuous on the interface, a discontinuity of the velocity in the tangential direction on the interface would be created instantaneously in general, so that the Kelvin-Helmholtz instability appears locally in space. As a result, the initial value problem for the full model turns out to be ill-posed. For more details, we refer to T. Iguchi, N. Tanaka, and A. Tani [14]. See also V. Kamotski and G. Lebeau [20] and D. Lannes [22]. In Section 4 we consider the linearized equations of the Kakinuma model around an arbitrary flow. After freezing the coefficients and neglecting lower order terms of the linearized equations, we calculate the linear dispersion relation and derive a stability condition, which is equivalent to on the interface, where P app 1 and P app 2 are approximate pressures of the waters in the upper and the lower layers in the Kakinuma model calculated from Bernoulli's equations, H 1 and H 2 are thickness of the upper and the lower layers, respectively, α 1 is a constant depending only on N , α 2 is a constant determined from {p 0 , p 1 , . . . , p N * }, and ∇ = (∂ x 1 , . . . , ∂ xn ) T is the nabla with respect to the horizontal spatial coordinates x = (x 1 , . . . , x n ). If ρ 1 = 0, then (1.8) coincides with the stability condition for the Isobe-Kakinuma model for water waves derived by R. Nemoto and T. Iguchi [28].
As in the case of the Isobe-Kakinuma model, the hypersurface t = 0 in the space-time R n ×R is characteristic for the Kakinuma model, so that one needs to impose some compatibility conditions on the initial data for the existence of the solution. Under these compatibility conditions, the non-cavitation assumption H 1 ≥ c 0 > 0 and H 2 ≥ c 0 > 0, and the stability condition (1.8), we will show in this paper that the initial value problem to the Kakinuma model is well-posed locally in time in Sobolev spaces. Here, we note that the coefficients α 1 and α 2 in the stability condition (1.8) converge to 0 as N, N * → ∞, so that the domain of stability diminishes as N and N * grow. This fact is consistent with the aforementioned properties of the full model.
Let us further comment on the significance of approximating an ill-posed system with well-posed systems. Firstly, while the initial value problem for the full model is ill-posed in Sobolev spaces, analytic solutions do exist, as shown by C. Sulem, P.-L. Sulem, C. Bardos, and U. Frisch [31] and C. Sulem and P.-L. Sulem [30] in the case where upper and lower boundaries are absent, and we expect that the corresponding solutions to the Kakinuma model provide valid approximations. Secondly, it should be recalled that the full model itself is a simplified model that discards effects that would stabilize the flow, especially vertical mixing across the pycnocline. In [22], D. Lannes considered another stabilizing effect, namely interfacial tension, and showed the existence and uniqueness of solutions with finite regularity to the corresponding initial value problem over a long time in the shallow water regime. The key physical mechanism at stake is that the Kelvin-Helmholtz instability which is responsible for ill-posedness issues occurs at sufficiently small spatial scale, so that it is possible to regularize the equations while being almost transparent to the behavior of the flow at large spatial scale, which is of practical interest for applications. Our results demonstrate that the Kakinuma model inherently incorporates such a stabilizing effect whose strength diminishes as N and N * grow, consistently with the expectation that the accuracy with respect to the full model increases.
As is well-known that the full model for interfacial gravity waves has a conserved energy where Ω 1 (t) and Ω 2 (t) are the upper and the lower layers, respectively. This is the total energy, that is, the sum of the kinetic energies of the waters in the upper and the lower layers and the potential energy due to the gravity. Moreover, T. B. Benjamin and T. J. Bridges [1] found that the full model can be written in Hamilton's canonical form where the canonical variable φ is defined by and the Hamiltonian H is the total energy E written in terms of the canonical variables (ζ, φ). Their result can be viewed as a generalization into interfacial gravity waves of Zakharov's Hamiltonian [32] for water waves. For mathematical treatments of the Hamiltonian for interfacial gravity waves, we refer to W. Craig and M. D. Groves [5] and W. Craig, P. Guyenne, and H. Kalisch [6]. The Kakinuma model has also a conserved energy E K , which is the total energy given by (1.9) with Φ 1 and Φ 2 replaced by Φ app 1 and Φ app 2 . Moreover, we will show that the Kakinuma model enjoys a Hamiltonian structure with a Hamiltonian H K the total energy in terms of canonical variables ζ and φ, where φ is defined by (1.10) with Φ 1 and Φ 2 replaced by Φ app 1 and Φ app 2 . This fact can be viewed as a generalization to the Kakinuma model for interfacial gravity waves of a Hamiltonian structure of the Isobe-Kakinuma model for water waves given by V. Duchêne and T. Iguchi [8].
The contents of this paper are as follows. In Section 2 we begin with reviewing the full model for interfacial gravity waves and derive the Kakinuma model. Then, we state one of the main results of this paper, that is, Theorem 2.1 about the well-posedness of the initial value problem to the Kakinuma model locally in time. In Section 3 we analyze the linear dispersion relation of the linearized equations of the Kakinuma model around the rest state in the case of the flat bottom and show (1.7). In Section 4 we derive the stability condition (1.8) by analyzing the linearized equations of the Kakinuma model around an arbitrary flow. In Section 5 we derive an energy estimate for the linearized equations with frozen coefficients and then transform the equations into a standard positive symmetric system by introducing an appropriate symmetrizer. In Section 6 we introduce several differential operators related to the Kakinuma model and derive elliptic estimates for these operators. In Section 7 we prove one of our main result, Theorem 2.1, by using the method of parabolic regularization of the equations. In Section 8 we prove another main result Theorem 8.4, which ensures that the Kakinuma model enjoys a Hamiltonian structure. Finally, in Section 9 we derive conservation laws of mass, momentum, and energy to the Kakinuma model together with the corresponding flux functions.
Notation. We denote by W m,p (R n ) the L p Sobolev space of order m on R n and H m = W m,2 (R n ). The norm of a Banach space B is denoted by · B . The L 2 -inner product is denoted by (·, ·) L 2 . We put ∂ t = ∂ ∂t , ∂ j = ∂ x j = ∂ ∂x j , and ∂ z = ∂ ∂z .
[P, Q] = P Q − QP denotes the commutator and [P ; u, v] = P (u·v)−(P u)·v −u·(P v) denotes the symmetric commutator. For a matrix A we denote by A T the transpose of A. For a vector φ = (φ 0 , φ 1 , . . . , φ N ) T we denote the last N components by φ ′ = (φ 1 , . . . , φ N ) T . We use the notational convention 0 0 = 0. We denote by C(a 1 , a 2 , . . .) a positive constant depending on a 1 , a 2 , . . .. f g means that there exists a non-essential positive constant C such that f ≤ Cg holds. f ≃ g means that f g and g f hold.

Kakinuma model and well-posedness
We begin with formulating mathematically the full model for interfacial gravity waves. In what follows, the upper layer, the lower layer, the interface, the rigid-lid of the upper layer, and the bottom of the lower layer, at time t, are denoted by Ω 1 (t), Ω 2 (t), Γ(t), Σ t , and Σ b , respectively. Then, the motion of the waters is described by the velocity potentials Φ 1 and Φ 2 and the pressures P 1 and P 2 in the upper and the lower layers satisfying the equations of continuity where ∆ = ∂ 2 1 + · · · + ∂ 2 n is the Laplacian with respect to the horizontal spatial coordinates x = (x 1 , . . . , x n ), and Bernoulli's equations The dynamical boundary condition on the interface is given by The kinematic boundary conditions on the interface, the rigid-lid, and the bottom are given by These are the basic equations for interfacial gravity waves. We can remove the pressures P 1 and P 2 from these basic equations. In fact, it follows from Bernoulli's equations (2.3)-(2.4) and the dynamical boundary condition (2.5) that Then, the basic equations consist of (2.1)-(2.2) and (2.6)-(2.10), and we can regard Bernoulli's equations (2.3)-(2.4) as the definition of the pressures P 1 and P 2 .
In the case of surface gravity waves, as shown by J. C. Luke [23], the basic equations have a variational structure and Luke's Lagrangian density is given by the vertical integral of the pressure P −P atm in the water region, where P atm is a constant atmospheric pressure. Therefore, it is natural to expect that even in the case of interfacial gravity waves the vertical integral of the pressure in the water regions would give a Lagrangian density L , so that we first define L pre by By Bernoulli's equations (2.3)-(2.4), this can be written in terms of the velocity potentials Φ 1 , Φ 2 , and the elevation of the interface ζ as The last two terms do not contribute to the calculus of variations of this Lagrangian, so that we define the Lagrangian density L (Φ 1 , Φ 2 , ζ) by and the action function J (Φ 1 , Φ 2 , ζ) by It is not difficult to check that the corresponding system of Euler-Lagrange equations is exactly the same as the basic equations (2.1)-(2.2) and (2.6)-(2.10) for interfacial gravity waves. We proceed to derive the Kakinuma model for interfacial gravity waves. Let Φ app 1 and Φ app 2 be approximate velocity potentials defined by (1.4) and define an approximate Lagrangian density L app (φ 1 , φ 2 , ζ) for φ 1 = (φ 1,0 , φ 1,1 , . . . , φ 1,N ) T , φ 2 = (φ 2,0 , φ 2,1 , . . . , φ 2,N * ) T , and ζ by which can be written explicitly as for i = 0, 1, . . . , N * , and Here and in what follows we use the notational convention 0 0 = 0. This system of equations is the Kakinuma model that we are going to consider in this paper. We consider the initial value problem to the Kakinuma model (2.14)-(2.16) under the initial condition For notational convenience, we decompose φ k as φ k = (φ k,0 , φ ′ k ) T for k = 1, 2 with φ ′ 1 = (φ 1,1 , . . . , φ 1,N ) and φ ′ 2 = (φ 2,1 , . . . , φ 2,N * ). Accordingly, we decompose the initial data φ k(0) as φ k(0) = (φ k,0(0) , φ ′ k(0) ) T for k = 1, 2. The hypersurface t = 0 in the space-time R n × R is characteristic for the Kakinuma model (2.14)-(2.16), so that the initial value problem (2.14)-(2.17) is not solvable in general. In fact, by eliminating the time derivative ∂ t ζ from the equations, we see that if the problem has a solution (ζ, φ 1 , φ 2 ), then the solution has to satisfy the N + N * + 1 relations (1 + |∇b| 2 )φ 2,j = 0 for i = 1, 2, . . . , N * , and Therefore, as a necessary condition the initial date (ζ (0) , φ 1(0) , φ 2(0) ) and the bottom topography b have to satisfy the relation (2.18)-(2.20) for the existence of the solution. These necessary conditions will be referred as the compatibility conditions.
The following theorem is one of our main results in this paper, which guarantees the wellposedness of the initial value problem to the Kakinuma model (2.14)-(2.17) locally in time.
Theorem 2.1. Let g, ρ 1 , ρ 2 , h 1 , h 2 , c 0 , M 0 be positive constants and m an integer such that m > n 2 + 1. There exists a time T > 0 such that for any initial data (ζ (0) , φ 1(0) , φ 2(0) ) and bottom topography b satisfying the compatibility conditions (2.18)-(2.20), the stability condition (1.8), and Remark 2.2. The term (∂ z (P app 2 − P app 1 )) | z=ζ in the stability condition (1.8) is explicitly given in (4.4). It includes the terms ∂ t φ k (x, 0) for k = 1, 2. Although the hypersurface t = 0 is characteristic for the Kakinuma model, we can uniquely determine them in terms of the initial data and b. For details, we refer to Remark 7.1. Under the condition (ρ 2 − ρ 1 )g > 0 and if the initial data and the bottom topography are suitably small, the stability condition (1.8) is automatically satisfied at t = 0. Remark 2.3. In the case N = N * = 0, that is, if we approximate the velocity potentials in the Lagrangian by functions independent of the vertical spatial variable z as Φ app k (x, z, t) = φ k (x, t) for k = 1, 2, then the Kakinuma model (2.14)-(2.16) is reduced to the nonlinear shallow water equations The compatibility conditions (2.18)-(2.20) are reduced to and the stability condition (1.8) is reduced to Therefore, we recover the conditions for the well-posedness in Sobolev spaces of the initial value problem to the nonlinear shallow water equations (2.22) proved by D. Bresch Given the initial data (ζ (0) , φ (0) ) for the canonical variables (ζ, φ) and the bottom topography b, the compatibility conditions (2.18)-(2.20) and the relation (2.23) determine the initial data (φ 1(0) , φ 2(0) ) for the Kakinuma model, which is unique up to an additive constant of the form (Cρ 2 , Cρ 1 ) to (φ 1,0(0) , φ 2,0(0) ). In fact, we have the following proposition, which is a simple corollary of Lemma 6.4 given in Section 6.

Linear dispersion relation
In this section we consider the linearized equations of the Kakinuma model (2.14)-(2.16) around the flow (ζ, φ 1 , φ 2 ) = (0, 0, 0) in the case of the flat bottom. The linearized equations have the form Putting we can rewrite the above equations as the following simple matrix form  where 1 = (1, . . . , 1) T and matrices A k,0 and A k,1 for k = 1, 2 are given by Therefore, the linear dispersion relation is given by where ξ ∈ R n is the wave vector, ω ∈ C is the angular frequency, and A k (ξ) = |ξ| 2 A k,0 + A k,1 for k = 1, 2. We can expand this dispersion relation as Here and in what follows, we use the notatioñ for a matrix A. Concerning the determinants appearing in the above dispersion relation, we have the following proposition, which was proved by R. Nemoto and T. Iguchi [28].
Thanks of this proposition and the dispersion relation (3.2), the linearized system (3.1) is classified into the dispersive system in the case N + N * > 0, so that the Kakinuma model (2.14)-(2.16) is a nonlinear dispersive system of equations.
Therefore, we can define the phase speed c K (ξ) of the plane wave solution to (3.1) related to the wave vector ξ ∈ R n by It follows from Proposition 3.1 that which is not consistent with the linear interfacial gravity waves However, as is shown by the following theorems the Kakinuma model gives a very precise approximation in the shallow water regime h 1 |ξ|, h 2 |ξ| ≪ 1 under an appropriate choice of the indices p i for i = 0, 1, . . . , N * .
Theorem 3.2. If we choose N * = N and p i = 2i for i = 0, 1, . . . , N * or N * = 2N and p i = i for i = 0, 1, . . . , N * , then for any ξ ∈ R n and any h 1 , h 2 , g > 0 we have where C is a positive constant depending only on N .

Stability condition
In this section, we will derive the stability condition (1.8) by analyzing a system of linearized equations to the Kakinuma model (2.14)-(2.16). We linearize the Kakinuma model around an arbitrary flow (ζ, φ 1 , φ 2 ) and denote the variation by (ζ,φ 1 ,φ 2 ). After neglecting lower order terms, the linearized equations have the form are approximate horizontal velocities in the upper and the lower layers on the interface, are approximate vertical velocities in the upper and the lower layers on the interface, and Here, P app 1 and P app 2 are approximate pressures in the upper and the lower layers calculated from Bernoulli's equations (2.3)-(2.4), that is, for k = 1, 2. Now, we freeze the coefficients in the linearized equations (4.1) and put Then, (4.1) can be written in the form  Therefore, the linear dispersion relation for (4.1) is given by where ξ ∈ R n is the wave vector and ω ∈ C the angular frequency. The left-hand side can be expanded as so that the linear dispersion relation is given simply as for k = 1, 2. The discriminant of this quadratic equation in ω is Therefore, the solutions ω to the dispersion relation (4.6) are real for any wave vector ξ ∈ R n if and only if Otherwise, the roots of the linear dispersion relation (4.6) have the form ω = ω r (ξ) ± iω i (ξ) satisfying ω i (ξ) → +∞ as ξ = (u 2 − u 1 )ξ and ξ → +∞, which leads to an instability of the problem. These consideration leads us to the following stability condition Here, we note that α 1 and α 2 are positive constants depending only on N and {p 0 , p 1 , . . . , p N * } and converge to 0 as N, N * → ∞. Therefore, as N and N * go to infinity the domain of stability diminishes.

Analysis of the linearized system
In this section, we still analyze the system of linearized equations (4.1) with frozen coefficients. We first derive an energy estimate for solutions to the linearized system by defining a suitable energy function, and then transform the linearized system into a standard symmetric form, for which the hypersurface t = 0 in the space-time R n × R is noncharacteristic. These results motivate the subsequent analysis on the nonlinear equations.

Energy estimate
With the notation (4.5), the linearized system (4.1) with frozen coefficients can be written in a symmetric form as We note that A 0 is symmetric in L 2 (R n ) whereas A 1 is skew-symmetric. Therefore, by taking L 2 -inner product of (5.1) with ∂ tU we have for any regular solutionU to (5.1), so that (U , A 0U ) L 2 would give a mathematical energy function to the linearized system (5.1) if we show the positivity of the symmetric operator A 0 in L 2 (R n ). We proceed to check the positivity. For simplicity, we consider first the case N = N * = 0 so that A 1,0 = A 2,0 = 1. Then, we see that Therefore, it is sufficient to analyze the positivity of this (2n + 1) × (2n + 1) matrix. The characteristic polynomial of this matrix is given by Therefore, the eigenvalues of the matrix are ρ 1 H 1 and ρ 2 H 2 of multiplicity n − 1 and λ 1 , λ 2 , λ 3 , which are the roots of the polynomial Here, we see that which is not necessarily positive even if u 1 = u 2 . Therefore, for the positivity of the symmetric operator A 0 we need a smallness of the horizontal velocities u 1 and u 2 . Such a condition is, of course, stronger restriction than the stability condition (4.8). This means that (U , A 0U ) L 2 is not an optimal energy function and we proceed to find out another one.
We are now considering the linearized system (5.1) with frozen coefficients, that is, Applying ∆ to the last equation in (5.2) we have Plugging the first and the second equations in (5.2) into (5.3) to removeψ 1 andψ 2 , we obtain In view of the relation following from Cramer's rule for k = 1, 2, the above equation forζ can be written as where u is an averaged horizontal velocity on the interface defined by Taking (5.4) into account, we consider the following constant coefficient second order partial differential equation where c 1 and c 2 are positive constants. By taking L 2 -inner product of (5.6) with (∂ t + u · ∇)ζ and using integration by parts, we see that for any regular solutionζ to (5.6). Here, we have Under this assumption, we obtain an energy estimate for the solutions to (5.6). Applying this consideration to (5.4), we see that the positivity condition is exactly the same as the stability condition (4.8), under which we can obtain an energy estimate for (5.4).
In [3] (see also [2]), D. Bresch and M. Renardy rewrote the the nonlinear shallow water equations (2.22), corresponding to the case N = N * = 0, as a scalar second order partial differential equation analogous to (5.4), and then used the abstract theory of T. J. R. Hughes, T. Kato, and J. E. Marsden [11] to obtain the local well-posedness of the initial-value problem under sharp hyperbolicity conditions, as mentioned in Remark 2.3. Our strategy is different as we rely on the symmetrization of the system and parabolic regularization to prove Theorem 2.1.
In view of (5.4) and the subsequent observation we rewrite the linearized system (5.1) with frozen coefficients in the form and v = u 2 − u 1 . By taking L 2 -inner product of this equation with (∂ t + u · ∇)U and using integration by parts, we see that for any regular solution to (5.1). We proceed to check the positivity of the symmetric operator A mod 0 in L 2 (R n ) under the stability condition (4.8). We see that On the other hand, the matrixÃ k,0 is nonsingular and its inverse matrix can be written as with a symmetric matrix Q k,0 for k = 1, 2. Moreover, q k,0 = det A k,0 detÃ k,0 = α k is positive and Q k,0 is nonnegative. In fact, for any ψ putting ζ and φ by We note that Q k,0 is not positive because it has a zero eigenvalue with an eigenvector 1. Now, for any φ, putting η = 1 · φ and ψ = A k,0 φ, we haveÃ k,0 from which we deduce the identity By using the decomposition (5.7) we see that Here, I 1 ≥ 0 since Q 1,0 and Q 2,0 are nonnegative, and so that it is sufficient to show the positivity of the matrix From Sylvester's criterion and since ρ k H k α k is positive for k = 1, 2, the positivity of the matrix A 0 is equivalent to Since v = u 2 − u 1 , under the stability condition (4.8) we have the positivity of A 0 , so that in view of (5.7) and the positivity of the matrix A k,0 for k = 1, 2 we finally obtain the equivalence Therefore, (A mod 0U ,U ) L 2 would provide a useful mathematical energy function.

Symmetrization of the linearized equations
We still consider the linearized equations (4.1) with frozen coefficients. However, for later use we defineφ 1 andφ 2 in place of (4.5) by Then, the linearized equations have the form In the following, we abbreviate simply l k (H k ) and A k (H k ) as l k and A k for k = 1, 2. We are going to show that the system can be transformed into a positive symmetric system of the form is the positive operator defined in the previous section with slight modification, and A is a skew-symmetric operator in L 2 (R n ). As before, we put v = u 2 − u 1 and define u by (5.5). Furthermore, we introduce the notation where α 1 and α 2 are positive constants defined by (4.7). Then, we have u = θ 2 u 1 + θ 1 u 2 and θ 1 + θ 2 = 1. We can also express u 1 and u 2 in terms of u and v as Applying ∆ to the third equation in (5.8) and differentiating the first and the second equations with respect to t, we obtain  In view of this, we introduce a symmetric matrix as and by Cramer's rule, Using these notations, we have Here, we see that On the other hand, taking the Euclidean inner product of the first and the second equations in (5.8) with −ρ 1 q 1 and ρ 2 q 2 , respectively, we obtain θ 2 (∂ tζ + u 1 · ∇ζ) + q 0 ρ 1 l 1 · ∆φ 1 = 0, θ 1 (∂ tζ + u 2 · ∇ζ) − q 0 ρ 2 l 2 · ∆φ 2 = 0, which are equivalent to (5.15) ∂ tζ + u · ∇ζ + q 0 ∆(ρ 1 l 1 ·φ 1 − ρ 2 l 2 ·φ 2 ) = 0, θ 1 θ 2 v · ∇ζ − q 0 ∆(θ 1 ρ 1 l 1 ·φ 1 + θ 2 ρ 2 l 2 ·φ 2 ) = 0.
It follows from the second equation in (5.15) that . We proceed to symmetrize the second term in the right-hand side of (5.14).
In the above calculation, we used the second equation in (5.15). Therefore, Summarizing the above calculations, if we define the symmetrizer A mod 0 by (5.17) then we obtain Therefore,U satisfies the symmetric system (5.11) with a skew-symmetric operator A defined by For the positive symmetric system (5.11), we can apply the standard theory for partial differential equations to show its well-posedness of the initial value problem. Moreover, these considerations help us to analyze the nonlinear problem (2.14)-(2.16).
Then, we have (L k,ij ) * = L k,ji for k = 1, 2, where (L k,ij ) * is the adjoint operator of L k,ij in L 2 (R n ). We also use u k and w k for k = 1, 2 defined by (4.2) and (4.3), which represent approximately the horizontal and the vertical components of the velocity field on the interface from the water region Ω k (t), respectively. Then, the Kakinuma model (2.14)-(2.16) can be written simply as Moreover, introducing φ 1 = (φ 1,0 , φ 1,1 , . . . , φ 1,N ) T , φ 2 = (φ 2,0 , φ 2,1 , . . . , φ 2,N * ) T , and we can write the Kakinuma model (2.14)-(2.16) more simply as By eliminating ∂ t ζ from the Kakinuma model, we obtain N + N * + 1 scalar relations the compatibility conditions can be written simply as We proceed to derive evolution equations for φ 1 and φ 2 . To this end, we differentiate the above compatibility conditions with respect to t and use equations of the Kakinuma model to eliminate ∂ t ζ. Then, we obtain (6.7) Here, we note that F 3 can be written in divergence form as On the other hand, the last equation in the Kakinuma model can be written as In view of these evolution equations (6.7)-(6.8) for φ 1 and φ 2 , we will consider the following equations for ϕ 1 and ϕ 2 .
Otherwise, for example, in the case p i = i (i = 0, . . . , N ) we obtain A similar estimate holds in other cases. These estimates give the desired one.
Although this lemma gives an a priori bound of the solution to (6.9), the equations in (6.9) do not have a good symmetry. In order to give an existence theorem to (6.9) with robust elliptic estimates, it is better to rewrite them in a symmetric form by introducing a good unknown variable. We introduce scalar functions ϕ 1 and ϕ 2 by (6.13) We also introduce second order differential operators P 1,i (H 1 ) (i = 1, . . . , N ) and Q 1 (H 1 ) acting on R N -valued functions ϕ ′ 1 = (ϕ 1,1 , . . . , ϕ 1,N ) T and P 2,i (H 2 , b) (i = 1, . . . , N * ) and Q 2 (H 2 ) acting on R N * -valued functions ϕ ′ 2 = (ϕ 2,1 , . . . , ϕ 2,N * ) T by (6.14) respectively, and put Then, we see easily that P 1 (H 1 ) and P 2 (H 2 , b) are symmetric in L 2 (R n ) and that where Q * denotes an adjoint operator of Q in L 2 (R n ). Therefore, we can rewrite (6.9) as These equations for (ϕ ′ 1 , ϕ 1 , ϕ ′ 2 , ϕ 2 ) do not have yet any good symmetry. But, it follows from the last equation that Using this we can remove ϕ 2 from the equations and obtain where (6.16) These equations for (ϕ ′ 1 , ϕ ′ 2 , ϕ 1 ) have a good symmetry and can be written in the matrix form which is symmetric in L 2 (R n ). Moreover, P(ζ, b) is positive in L 2 (R n ) as shown in the following lemma.
By this lemma, the explicit expression (6.18) of the operator P(ζ, b), and the standard theory of elliptic partial differential equations, we can obtain the following lemma. Lemma 6.3. Let ρ 1 , ρ 2 , h 1 , h 2 , c 0 , M be positive constants and m an integer such that m > n 2 +1. There exists a positive constant C = C(ρ 1 , ρ 2 , h 1 , h 2 , c 0 , m) such that if ζ and b satisfy Moreover, the solution is unique up to an additive constant to ϕ 1 .

Construction of the solution
In this section, we will prove Theorem 2.1 one of the main theorems in this paper. One possible strategy to construct the solution of the initial value problem to the Kakinuma model (2.14)-(2.16) would consist in firstly transforming the equations into a quasilinear positive symmetric system, that is, a quasilinear version of the positive symmetric system (5.11), secondly applying the method of parabolic regularization to construct the solution of the transformed system, and finally to show that the solution to the transformed system is in fact the solution of the Kakinuma model if we further impose the compatibility conditions (2.18)-(2.18) on the initial data. Here, in order to avoid the heavy computations that would be involved when following this strategy, we find it more convenient to instead apply the method of parabolic regularization to the Kakinuma model directly.

Parabolic regularization of the equations
We remind that the Kakinuma model (2.14)-(2.16) can be written compactly as (6.4), that is, 0 , φ 1,1 , . . . , φ 1,N ) T , φ 2 = (φ 2,0 , φ 2,1 , . . . , φ 2,N * ) T , l k and L k for k = 1, 2 are defined in (6.3), and Here u k and w k for k = 1, 2 are defined by (4.2) and (4.3) respectively. We regularize the Kakinuma model by adding artificial viscosity terms as We are going to show the existence of the solution to the initial value problem for this regularized Kakinuma model under the initial conditions For this regularized Kakinuma model, the compatibility conditions for the existence of the solution have the same form as the original Kakinuma model, that is, where L 1,i (H 1 ) for i = 0, 1, . . . , N and L 2,i (H 2 , b) for i = 0, 1, . . . , N * are defined in (6.5). Here, we note the identities We also note that f 3 (ζ, φ 1 , φ 2 , b) can be written in a divergence form as Therefore, applying the operator ∂ t − ε∆ to (7.5) we obtain for i = 1, 2, . . . , N * , On the other hand, we have N + N * + 2 evolution equations for one scalar function ζ. To select an appropriate evolution equation for ζ, we will use the notation defined by (5.13). We note that they depend on the unknown functions H 1 and H 2 . Taking Euclidean inner products of the first and the second equations in (7.3) with ρ 1 q 1 and ρ 2 q 2 , respectively, adding the resulting equations, and using the relation −ρ 1 l 1 · q 1 + ρ 2 l 2 · q 2 = 1, we obtain Plugging this into (7.6) and noting the last equation in (7.3), we have Therefore, thanks to Lemma 6.4 we obtain where G 1 = (G 1,0 , G 1,1 , . . . , G 1,N ) T and G 2 = (G 2,0 , G 2,1 , . . . , G 2,N * ) T are defined as a solution to the following equations (7.10) Precisely speaking, (G 1 , G 2 ) are defined uniquely up to an additive constant of the form (Cρ 2 , Cρ 1 ) to (G 1,0 , G 2,0 ). However, this indeterminacy does not cause any difficulties in the following arguments.
Proof. By the construction of the solution, we easily see that it satisfies (7.8) and in particular the last equation in (7.3). Therefore, it is sufficient to show that it satisfies also the first two equations in (7.3). By (7.7) and (7.8), we have so that by the uniqueness of the solution to the initial value problem of the heat equation, if the initial data satisfy the compatibility conditions (7.5), then the solution satisfies also (7.5) for all t ∈ [0, T ε ). Particularly, we obtain so that by the last equation in the compatibility conditions (7.5) we have Taking Euclidean inner products of the first and the second equations with ρ 1 q 1 and ρ 2 q 2 , respectively, adding the resulting equations, and using the relation −ρ 1 l 1 · q 1 + ρ 2 l 2 · q 2 = 1, we obtain which together with (7.7) implies Plugging this into (7.12), we see that the solution satisfies the first two equations in (7.3).
We have also Therefore, where we used the identity (5.7). Since Q 1,0 and Q 2,0 are nonnegative and in view of and the analysis in Section 5.1, we can show the desired equivalence.
the stability condition (7.19) with c 0 replaced by 2c 0 , and the compatibility conditions (7.5), then for any ε ∈ (0, 1] the solution (ζ ε , φ ε 1 , φ ε 2 ) constructed in Lemmas 7.2 and 7.3 satisfies Proof. Once again we simply write U = (ζ, φ 1 , φ 2 ) T in place of (ζ ε , φ ε 1 , φ ε 2 ) T . We define an energy function E m (t) by We assume that the solution (ζ(t), φ 1 (t), φ 2 (t)) satisfies (7.18) and the stability condition (7.19) for 0 ≤ t ≤ T . Then, the energy function E m (t) is equivalent to E m (t) = (ζ(t), ∇φ 1,0 (t), ∇φ 2,0 (t)) 2 Furthermore, we assume that for 0 ≤ t ≤ T , where the constant M 1 and the time T will be determined later. In the following we simply write the constants depending only on (g, ρ 1 , ρ 2 , h 1 , h 2 , c 0 , C 0 , M 0 ) by C 1 and the constants depending also on M 1 by C 2 . They may change from line to line. Then, it holds that for j = 0, 1, 2, . . .. We are going to evaluate the evolution of the energy function E m (t). To this end, we take the L 2 -inner product of (7.17) with ∂ t U β − ε∆U β + (u · ∇)U β and use integration by parts to get Here, we see that so that for 1 ≤ |β| ≤ m we have Similar estimate can be obtained in the case |β| = 0 more directly. On the other hand, it follows from (7.9) that Therefore, we obtain which yields Putting M 1 = 2C 1 and taking T > 0 so that C 2 T ≤ C 1 , we obtain by a continuity argument that (7.20) holds for 0 ≤ t ≤ T . It remains to show that (ζ(t), φ 1 (t), φ 2 (t)) satisfies (7.18) and the stability condition (7.19) for 0 ≤ t ≤ T . By the Sobolev embedding theorem, (7.7), and (7.9), we see that which yields (7.18), except for the estimate for a, by taking T > 0 sufficiently small. We now turn to the stability condition (7.19). In order to evaluate ∂ t a, we need to obtain estimates for ∂ t G ′ k for k = 1, 2. Differentiating (7.10) with respect to t, we have Therefore, by Lemma 6.4 with k = m − 2 we obtain On the other hand, it follows from (7.7) and (7.9) that Thus,

Hamilton's canonical form
We proceed to show that the Kakinuma model (8.1) is equivalent to Hamilton's canonical form with the Hamiltonian defined by (8.9). In the following, we fix b ∈ W m,∞ with m > n 2 + 1 and put U m b = {ζ ∈ H m ; inf x∈R n (h 1 − ζ(x)) > 0 and inf x∈R n (h 2 + ζ(x) − b(x)) > 0}, which is an open set in H m . We also use the function spaceH k = {φ ; ∇φ ∈ H m−1 }. For Banach spaces X and Y , we denote by B(X ; Y ) the set of all linear and bounded operators from X into Y . By Lemma 6.4, we see easily the following lemma.
Proof. Let us calculate Fréchet derivatives of the Hamiltonian H K (ζ, φ). Let us consider first U m b × H 2 ∋ (ζ, φ) → H K (ζ, φ). For anyφ ∈ H 2 , we see that where we used (8.5) and Lemma 8.1. The above calculations are also valid when (φ,φ) ∈ H 1 ×H 1 , provided we replace the L 2 inner products with the X ′ -X duality product where X =H 1 × (H 1 ) N or X =H 1 × (H 1 ) N * for the first lines, and X =H 1 for the last line. This gives the first equation of the lemma.
Similarly, for any (ζ, φ) ∈ U m b ×H 2 andζ ∈ H m we see that Here, we have where we used the identity stemming from (8.10). Again, the above identities are still valid for (ζ, φ) ∈ U m b ×H 1 provided we replace the L 2 inner products with suitable duality products. This concludes the proof of the Fréchet differentiability, and the second equation of the lemma. Now, we are ready to show another main result in this paper.
where F is given by (8.2) and