Hamiltonian aspects of 3-layer stratified fluids

The theory of 3-layer density stratified ideal fluids is examined with a view towards its generalization to the n-layer case. The focus is on structural properties, especially for the case of a rigid upper lid constraint. We show that the long-wave dispersionless limit is a system of quasi-linear equations that do not admit Riemann invariants. We equip the layer-averaged one-dimensional model with a natural Hamiltonian structure, obtained with a suitable reduction process from the continuous density stratification structure of the full two-dimensional equations proposed by Benjamin. For a a laterally unbounded fluid between horizontal rigid boundaries, the paradox about the non-conservation of horizontal total momentum is revisited, and it is shown that the pressure imbalances causing it can be intensified by three-layer setups with respect to their two-layer counterparts. The generator of the x-translational symmetry in the n-layer setup is also identified by the appropriate Hamiltonian formalism. The Boussinesq limit and a family of special solutions recently introduced by de Melo Virissimo and Milewski are also discussed.


Introduction
Density stratification is an important aspect of fluid dynamics, being inherent to a variety of phenomena concerning both the ocean and the atmosphere. In particular, displacement of fluid parcels from their neutral buoyancy position within a stratified flow can result in internal wave motion, whose evolution is of fundamental importance to energy propagation and distribution in both oceanic and atmospheric settings. Simplified one-dimensional models (in particular, their quasi-linear limit) have been introduced to isolate key elements in the dynamics of these phenomena, and have been the subject of a number of investigations, e.g., from [34] to the more recent [26,14,22], through [16,17,31,18,19].
Such simplified models are the focus of this paper. Two main classes of configurations have been examined in the literature: that of a free surface under constant pressure, and that of a fluid confined in a vertical channel by rigid horizontal plates. We first frame our discussion in more general terms by using a formalism, employed in [17] for the case n = 2, encompassing both classes. Within this setting, we first recall the derivation of the dispersionless n-layer equations by means of the so-called hydrostatic approximation, which allows to express the pressure in terms of layer thicknesses, plus a reference (interfacial) pressure. The free-surface case is then briefly addressed, in particular, for the sake of concreteness, when n = 2. For a comprehensive approach see, e.g., [15], and [29,1] for a discussion of a variational approach. Extension of these results, considering dispersive terms (as in (2.6) below), can be found in [3,35]. In this respect, our main goal is to highlight the following facts: (i) the non-existence of Riemann invariants, and (ii) the existence of a natural Hamiltonian structure.
To illustrate our approach, we focus on a stratified fluid composed of three homogeneous layers of constant density ρ 1 < ρ 2 < ρ 3 , confined in a vertical channel of fixed height h. We show that the apparently paradoxical feature of non-conservation of horizontal momentum, already hinted at in [4] and thoroughly discussed in [9,10] for the 2D Euler equations, can be detected, not unexpectedly, in the multiply stratified case, as also noted in [22]. We notice that the added degree of freedom afforded by the 3-layer setting makes the class of initial conditions leading to lack of horizontal momentum conservation much larger than its two-layer counterpart, since one can have a nonvanishing time-variation of the horizontal momentum at the first order in the density differences ρ j+1 − ρ j even with vanishing initial velocities. Hamiltonian aspects of such models have also been considered, notably in [5,21,20], with an emphasis on the 2-layer case. We approach this problem from the setup of Hamiltonian reduction methods, and show that the Hamiltonian structure of the resulting reduced equations naturally arises via a process of reduction from the Hamiltonian structure introduced by Benjamin [4] in the study of incompressible, density-stratified Euler system in two spatial dimensions. The Hamiltonian formalism allows us to derive the equations of motion as a system of conservation laws, and address the paradox mentioned above by identifying the generator of the x-translational symmetry of the problem. We present explicit initial conditions that illustrate the lack of horizontal momentum conservation, and note that, unlike the 2-layer case, even with zero initial velocity the horizontal momentum is not conserved. Next, we show explicitly that the reduced system does not admit Riemann invariants. We extend our study to the so-called Boussinesq approximation of retaining density differences for the buoyancy terms only (while disregarding density differences for the inertial terms). This limit is readily obtained in our canonical formalism, together with its effective Hamiltonian density and "symmetric" solutions.

Sharply stratified n-layered Euler fluids
The incompressible Euler equations for the velocity field u = (u, w) and non-constant density ρ(x, z, t), in the presence of gravity −gk and with horizontal fluid domain bottom at z = 0, are Dρ Dt = 0 , ∇ · u = 0 , D(ρu) Dt with boundary conditions u(x = ±∞, z, t) = 0 for all z ≥ 0 and w(x, 0, t) = 0 for all x ∈ R. As usual, D/Dt = ∂/∂t + u · ∇ is the material derivative. A classical way to reduce the dimensionality of the model is to study the spatial averages of certain fields, such as velocity or density, along suitable regions, as set forth by T. Wu [38]. In the case of fluids stratified by gravity, where the vertical direction plays a distinguished role with respect to the other two directions, one can study the evolution of vertical means of the fields. In particular, in [38] it was established that for a density stratified fluid it holds, for each layer of thickness where f (x, z, t) is any quantity, and the elevations ζ α are related with the thicknesses of the fluid's layers via Due to the intrinsic nonlinearity of the Euler equations, by applying the averaging procedure to (2.1) we obtain in general a non-closed system for the layer thicknesses η i and the layer mean velocities u i defined as Indeed, additional approximations must be introduced to close the system. In particular, the layer averaged equations of motion for a two-layer incompressible Euler fluid in an infinite channel were derived from the corresponding Euler equations in [17] and extended the case of n layers in [15]. Motions of typical wavelength L were considered, under the assumptions that the ratios can be considered small. Here h is the total height of the channel, while η i is the thickness of the i-th fluid homogeneous layer. By using the so-called columnar ansatz (see, e.g., [38]), and by noticing that equation (2.5), together with the incompressibility of each layer, implies that the ratio of vertical and horizontal velocities scales as , in [17] it was shown that the 2 + 1-dimensional Euler equations (2.1) reduce to the 1 + 1 dimensional equations x , the u i 's are the layer-mean velocities, ρ i is the density of the i-th layer, and P (x, t) is the interfacial pressure. These results were generalized in [15,2] to the case of n > 2 layers, together with the discussion of suitable conditions leading to systems of quasi-linear equations, which we summarize hereafter for the reader's convenience.
In fact, in this paper we shall be mainly interested in the leading order approximation in the long-wave small parameter series expansion, the so-called "hydrostatic limit," which discards the higher order dispersive terms D i and, consistently, gives the following simple expression for the pressure p i (x, z) with z in the i-th layer, i.e., for ζ i < z < ζ i−1 : (2.7) In these relations, P 0 denotes an x and t-dependent reference pressure which, without loss of generality, we can set from now on to be the pressure at the bottom of the channel.
In the hydrostatic approximation, the equations of motion for each layer contain the averaged x-derivative of the pressure p ix (which, owing to the noncommutative property of exchanging derivative with layer-averages, is in general different from the x-derivative of the averaged pressure Notice that in the hydrostatic approximation the x-derivative of the pressure does not depend on z and therefore its layer mean is readily computed as (2.9) We thus obtain the set of 2n equations for the 2n + 1 dependent variables η i , u i , P 0 .
3 The free surface equations: the case n = 2 A standard way to close system (2.10) consists of considering the free surface case (see, e.g. [34,26,14,20,7,15]). In this instance, the reference pressure P 0 can be eliminated by observing that Figure 1: Free-surface, two-layer fluid setup and relevant notation: P (x, t) = 0 is the free surface (air) pressure, ζ 0 and ζ 1 are the free surface and interface locations, and η 1 and η 2 are the layer thicknesses of the fluids with densities ρ 1 and ρ 2 , respectively. the pressure at the fluid's domain upper boundary z = ζ 0 can be consistently set to vanish, as depicted in Figure 1. Thus P 0 can be expressed, thanks to eq. (2.7) for i = 1, as So, for i = 1, . . . , n, the pressures p i become (linear) functions of the thicknesses η i and the system closes. An example of this procedure in the 2-layer case is as follows (the generalisation to arbitrary n being straightforward): the pressure in the lowest layer is p 2 (x, z) = P 0 − gρ 2 z. In the upper layer, we have since ζ 1 = η 2 and ζ 0 − ζ 1 = η 1 . Requiring the vanishing of the pressure at the free boundary z = ζ 0 = η 1 + η 2 yields, as expected, P 0 = g(ρ 1 η 1 + ρ 2 η 2 ), and so Thus, the equations of the motion deduced from (2.10) are These equations can be written in the standard matrix form for quasilinear systems as where the characteristic matrix reads The question of the existence of Riemann invariants for this quasi-linear system can be easily solved by computing the so-called Haantjes tensor H of the matrix A 2 [30,6,28], whose vanishing ensures, in the hyperbolic case, the existence of the Riemann invariants, and is not granted for systems with more than 2 quasi-linear equations. The computation can be performed by means of standard computer algebra programs. We recover the non-existence of Riemann invariants, conjectured in [34] and proved in [14], since it can be easily checked that the Haantjes tensor of A 2 has three non-vanishing components, Equations (3.4) do admit a natural Hamiltonian formulation. To show this, we pass from averaged velocities u i to the averaged momentum coordinates µ i defined by µ i = ρ i u i , with i = 1, 2. The evolution equations translate into A direct inspection shows that (3.8) admit the Hamiltonian formulation where the vector δ η 1 H, δ η 2 H, δ µ 1 H, δ µ 2 H is the differential of the effective energy functional The two-layer system (3.8) is a non-diagonalizable Hamiltonian system of conservation laws in Tsarev's sense [36,37]. The non-existence of Riemann invariants has deep consequences on the existence of conservation laws and on solutions of quasi-linear systems, and as such is a topic much studied in the literature. Notable for our purposes are the results of [36,37], showing that six quantities are guaranteed to be conserved by strong solutions: the Hamiltonian functional (3.10), +∞ −∞ η j dx, +∞ −∞ µ j dx, j = 1, 2 (these are Casimir functionals of the Poisson brackets), and the total horizontal momentum Π (x) = +∞ −∞ (η 1 µ 1 + η 2 µ 2 ) dx, this being the generator of the x-translational symmetry in the free surface case. Furthermore, in [33] and [1] the authors show that these conserved quantities are the only ones whose densities do not explicitly depend on x. Also notable in this regard are the results and conjectures in [27], about the complete integrability of quasilinear systems which are linearly degenerate (termed weakly nonlinear therein) but lack a complete set of Riemann invariants. Last but not least, as far as strong solutions of the quasilinear system are concerned, the non-existence of Riemann invariants removes the possibility of constructing, through their level sets, a coordinate system in the hodograph space. Existence of such a system provides a powerful tool for assuring the persistence of solutions in the hyperbolic region and for estimating times of possible shock formation. 4 The rigid lid constraint: the case n = 3 Another way of closing system (2.10) is by enforcing the so-called rigid lid constraint. Its basic physical motivation relies on the fact that surface waves effectively decouple from internal ones. This is a subtle process, since, from a physical point of view, the interfacial hydrostatic pressure cannot be computed a priori as in the free surface case, but rather it has to be eliminated by means of a subtraction process.
The rigid upper lid translates in the (geometrical ) constraint where h is the total height of the vertical fluid "channel". The sum of all layer-mass conservation that can be termed the dynamical constraint.
The system of equations (2.10) written for the momentum densities η k u k is The sum of all these equations gives the expression of the horizontal gradient of the pressure at the bottom of the channel, This explicit form for P 0 togetheer with the geometric and dynamical constraints (4.1),(4.4) guarantee the closure of (2.10) as a system of 2n − 2 equations for the 2n − 2 variables, say, (η α , u α ) with α = 2, . . . , n, that govern the evolution of a sharply n-layer stratified fluid with upper rigid lid, i.e,, confined in a vertical channel.
In analogy with the procedure we followed for the free surface case with 2 layers, rather than discussing the general case, we now concentrate our analysis to the case with four dependent variables, that is, to the rigid lid case with 3 layers.

Three-layer rigid lid case
Our geometrical setting is described by Figure 2. From (2.7) it is clear that the explicit expression for the hydrostatic pressures in the layers is and therefore The mass conservation equations read which imply, together with the vanishing asymptotic conditions on the velocities, the constraint as seen in (4.4), since The Euler equations lead to The mass conservation equations and the Euler equations imply the following set of evolution equations (see (4.5)): (4.12) By taking the sum of these three evolution equations and by using the constraint η 3 u 3 + η 2 u 2 + η 1 u 1 = 0, we can solve for the x-derivative of the pressure P 0 as (see (4.6)) where we used again the fact that ζ 2 = η 3 and (4.14) Hence we have the following equations of motion for (η 2 , η 3 , u 2 , u 3 ), where P 0 x is given by (4.14) and

Pressure imbalances and non-conservation of horizontal momentum
Horizontal momentum conservation can be violated in the dynamics of an ideal stratified fluid with a rigid lid constraint. This phenomenon, already suggested by Benjamin in [4], was highlighted and substantiated in [9,10]. The lack of horizontal momentum conservation is surprising, as the only acting body-force field is the vertical gravity, constraint forces generated through pressure against horizontal plates are necessarily vertical, and the fluid is free to move laterally. In the horizontal channel set-up, whenever hydrostatic conditions apply at infinity, the violation of momentum conservation is proportional (up to terms that arise from possibly different configurations at x = ±∞) to the difference P ∆ of the layer-averaged pressure at the far ends of the channel. In turn, in the one-dimensional models for sharply stratified configurations such as those considered here, the layer averaged pressure imbalances are eventually those given by the quantity P 0 (x, t).
Pressure imbalances at the far-ends of the channel (see also [22]), can be detected looking at equation (4.13). Indeed, its right-hand side in not a total x-derivatives, and so +∞ −∞ P 0 x dx need not vanish even for configurations with asymptotically flat interfaces ζ i at the far-ends of the channel. A very simple example of such a configuration with non-vanishing P ∆ is given in Figure 3 by splicing together arcs of parabolae for the interface locations. We choose, along with u i = 0 for i = 1, 2, 3, A specific feature of 3-(and, a fortiori, n-) layered fluid configurations is that this "paradox" about non-conservation of horizontal momentum is an effect of order one in the single density differences ρ 2 − ρ 1 and ρ 3 − ρ 2 , even with initial null velocities, as opposed to 2-layer, rigid lid configurations, Figure 3: A three-layer configuration with a non-vanishing initial pressure imbalance. Setting ρ 1 = 0.5, ρ 2 = 0.75, ρ 3 = 1, and g = h = 1, gives P ∆ −0.0037395.
where non-conservation of the horizontal total momentum is, with null-velocities, quadratic in the density difference ρ 2 − ρ 1 (see [10]). Indeed, equation (4.14) reduces to The numerator is equal, up to total derivatives, to g(ρ 2 − ρ 1 )(ρ 3 − ρ 2 )η 3x η 2 /(ρ 1 ρ 2 ), thus revealing the leading order dependence of P ∆ on the density differences for general initial layer thicknesses η i . Finally, a quick glance at the expression of P 0 x for the n-layer case given by equation (4.6) shows that the lack of conservation of the horizontal momentum persists for n > 3.
Remark 4.1 Non-dispersive, quasilinear equations governing stratified flows are generically systems of mixed type [19,8], i.e., exhibit hyperbolic to elliptic transitions corresponding to eigenvalues of the characteristic matrix becoming complex across some "sonic" surfaces. The presence of elliptic domains is associated with instabilities, and, as well known, the evolution from initial values in the elliptic domain is ill-posed. One can prove that the initial data with initial null velocities, and with initial profiles of the kind exemplified by equation (4.16), with density ratios ρ 2 /ρ 3 = 3/4, ρ 1 /ρ 3 = 1/2, lie inside the hyperbolic region, as the characteristic matrix has four real eigenvalues. By continuity, the system's evolution must remain hyperbolic for a possibly small but non-zero finite time. In fact, numerical integration of this initial value problem (not reported here) by using the reduced system derived below in §5.3, indicates that the evolution remains hyperbolic for a relatively long time, at least till the final time t = 5 we have tried, with a shock forming at around time t 2.8 (with the nondimensional parameters as in the example of Figure 3).

The 3-layer rigid lid case: the Hamiltonian structure
In this section, we extend our study of the 3-layer rigid-lid system considered in §4.1 with the aim of endowing its equations of motion with a Hamiltonian structure. In the free surface 2-layer case (as well as in the free surface 3-layer case, see the Appendix) the resulting equations, when written in the momentum coordinates µ k = ρ k u k , k = 1, 2, 3, are sufficiently simple to be cast in Hamiltonian form at a glance. This is not the case in the presence of a rigid top lid, basically owing to the nonlinear nature of the dynamical constraint (4.4). A natural avenue of attack to account for this is by means of a suitable Hamiltonian reduction of the Poisson structure for continuously stratified flows, originally introduced in [4], by extending to the multiple layer case a technique introduced in [12].

The 2D Benjamin model for heterogeneous fluids in a channel
Benjamin [4] proposed and discussed a set-up for the Hamiltonian formulation of an incompressible stratified Euler system in 2 spatial dimension, which we hereafter summarize for the reader's convenience.
Benjamin's idea was to consider, as basic variables for the evolution of a perfect inviscid and incompressible but heterogeneous fluid in 2D, subject to gravity, the density ρ together with the "weighted vorticity" Σ defined by The equations of motion for these two fields, ensuing from (2.1), are They can be written in the form is simply given by the sum of the kinetic and potential energy, D being the fluid domain. The most relevant feature of this coordinate choice is that (ρ, Σ) are physical variables. Their use, though confined to the 2D case with the above definitions, allows one to avoid the introduction of Clebsch variables (and the corresponding subtleties associated with gauge invariance of the Clebsch potentials) which are often used in the Hamiltonian formulation of both 2D and the general 3D case (see, e.g., [39]). As shown by Benjamin, equations (5.3) are a Hamiltonian system with respect to a Lie-theoretic Hamiltonian structure, that is, they can be written as for the Poisson bracket defined by the Hamiltonian operator (5.5)

The reduction process
By means of the Heaviside θ and Dirac δ generalized functions, a three-layer fluid configuration can be introduced with constant densities ρ i and velocity components u i (x, z), w i (x, z), i = 1, 2, 3 (for the upper i = 1, middle i = 2, and lower layer i = 3, respectively), with interfaces ζ 1 and ζ 2 . The global density and velocity variables can be written as Thus, the density-weighted vorticity Σ = (ρw) where Ω i = w ix − u iz for i = 1, 2, 3. Next, we assume the motion in each layer to be irrotational, so that Ω i = 0 for all i = 1, 2, 3. Therefore the density weighted vorticity acquires the form In the long-wave asymptotics described in Section 2, we have, at the leading order in the long wave expansion parameter, = h/L 1, i.e., we can neglect the vertical velocities w i and trade the horizontal velocities u i with their layer-averaged counterparts. Thus, from (5.7) and recalling the first of (5.6), we obtain and

(5.10)
To invert the map (5.8) we can integrate along the vertical direction z. Contrary to 2-layer case described in [12], now we need to integrate on two vertical slices of the channel in order to obtain four x-dependent fields. To this end, we extend the Benjamin's manifold M (parametrized by ρ(x, z) and Σ(x, z)) by the space F of isopycnals z = f (x), and trivially extend by 0 the Poisson structure J B given by (5.5). As shown in [11], the operator (5.5) provides a proper Poisson structure under the condition that the density ρ be constant at the physical boundaries z = 0 and z = h. Once equipped with the isopycnal z = f (x), we can define a projection π : M ≡ M × F → (C ∞ (R)) 4 by means of We choose to include the manifold of three-layer fluid configurations in the space M × F by adding to (5.8) the following choice for the isopycnal f : hereafter denoted by ζ. Equations (5.8) and (5.12) define a submanifold I of M, in one-to-one correspondence with the manifold S of 3-layer fluid configurations considered in Section 4.1, as it is parametrised by four functions of the horizontal coordinate x, and this set is equivalent to the 4-tuple (η 1 , η 2 , u 1 , u 2 ). For further reference, we explicitly remark that on I the projection π is defined by the 4-tuple (5.13) Further, as we shall see below, theζ component of the tangent vector (ρ,Σ,ζ) can be lineraly expressed in terms ofρ.
To obtain a Hamiltonian structure on the manifold S by reducing Benjamin's parent structure (5.5), we perform the following steps: 1. Starting from a 1-form on the manifold S, represented by the 4-tuple (α where π * is the tangent map to (5.13).
2. We apply Benjamin's operator (5.5) to the lifted one form α M to get the vector field This construction essentially works as in the two-layer case considered in [12], provided one point is taken into account, namely that the integrals in the second and fourth component of π have a variable upper bound. We have, for (ρ,Σ) ∈ T M I , We remark that, since the inequalities ζ 2 < ζ ≡ ζ 1 + ζ 2 2 < ζ 1 hold in the strict sense, the second term of this vector's fourth component vanishes. The same cannot be said for the analogous term in the vector's second component, since ρ(x, ζ) − ρ 3 = ρ 2 − ρ 3 = 0. As anticipated above, on I the mean of the tangent vector component coming from the ζ's,ζ = (ζ 1 +ζ 2 )/2, can be expressed in terms of the tangent vector componentρ. To this end, we can use the analogue of relations (5.9), which generically giveρ .
Thus, formula (5.16) for the components of the push forward π * does not explicitly depend on the last componentζ, and can be written as Applying the Poisson tensor (5.5) to this 1-form yields, after some manipulations that crucially use the fact that products of Dirac's δ supported at different locations can be consistently set to zero, the vector field S,x .

(5.25)
The push-forward under the map π * of this vector field gives the following four components: S,x , where we used definition (5.21) of ∆ ρ and the fact that the terms with the z-derivatives of the Dirac's δ give a vanishing contribution since they are integrated against functions of x only. From (5.26) we obtain the expression of the reduced Poisson tensor P on S, in the coordinates (ξ 1 , ξ 2 , τ 1 , τ 2 ), as (5.27) The variables (ξ 1 , ξ 2 , τ 1 , τ 2 ) are related to (ζ 1 , ζ 2 , σ 1 = ρ 2 u 2 − ρ 1 u 1 , σ 2 = ρ 3 u 3 − ρ 2 u 2 ) by (5.28) Solving these relations with respect to the ζ i 's and the σ i 's gives: (5.29) A straightforward computation shows that in these coordinates the Poisson operator (5.27) acquires the particularly simple form (5.30) Remark 5.1 According to the terminology favored by the Russian school, for Hamiltonian quasilinear systems of PDEs (see, e.g., [37]) the coordinates (ξ 1 , ξ 2 , τ 1 , τ 2 ) and, a fortiori, the coordinates (ζ 1 , ζ 2 , σ 1 , σ 2 ), are "flat" coordinates for the system. In view of the particularly simple form of the Poisson tensor (5.30), the latter set could be called a system of flat Darboux coordinates.

Remark 5.3
The steps performed here are borrowed from the Poisson reduction theory of Marsden, Ratiu and Weinstein (see, e.g., [32]); however, a more elaborate exposition of our procedure within such a general framework falls beyond the scope of the present work and we omit it here.

The reduced Hamiltonian
The full energy (per unit length) of the 2D fluid in the channel is just the sum of the kinetic and potential energy, The potential energy is readily reduced, using the first of (5.6), to When the layer thicknesses are not asymptotically zero, both energies can be appropriately renormalized subtracting the far field contributions of ζ i . To obtain the reduced kinetic energy density, we use the fact that at order O( 2 ) we can disregard the vertical velocity w, and trade the horizontal velocities with their layer-averaged means. Thus the x-density is computed as so that the reduced kinetic energy is We now use the dynamical constraint for localized solutions, whereby velocities vanish at infinty, to obtain Next, we express u 2 , u 3 in terms of σ 1 , σ 2 as The kinetic energy density turns out to be, in the new set of variables, so that the Hamiltonian functional is Explicitly, the equations of motion can be written as conservation laws, where the gradient of the Hamiltonian is, explicitly, As can be seen from the above formulae, even with the simple, constant Hamiltonian operator (5.30) in (5.43), this Hamiltonian gradient leads to rather lengthy (albeit explicit) expressions for the evolution equations, which are not particularly illuminating and hence are omitted here. Suffices to say that the conciseness of the Hamiltonian formalism allows to show quickly the existence of at least six conserved quantities, In particular, the quantity K indicates how the momentum paradox mentioned in Section 4.2 can be viewed from the canonical formulation of the evolution equations. In fact, we see that, rather than the horizontal momentum Π (x) , it is K that plays the role of generator of x-translations. By expressing its density K = ζ 1 σ 1 + ζ 2 σ 2 in terms of the horizontal mean velocities we obtain is the total horizontal momentum density. Thus, the total momentum Π (x) = +∞ −∞ π (x) dx is not conserved, while K = +∞ −∞ Kdx is.

The Boussinesq approximation
A dramatic simplification of the problem is provided by the so-called Boussinesq approximation, that is, the double scaling limit Since the Poisson tensor (5.30) is independent of the densities, this double scaling limit can be implemented most simply within the Hamiltonian formulation. While the potential energy is unchanged, from (5.41) the kinetic energy acquires the form so that the Hamiltonian energy functional in this Boussinesq limit is The ensuing equations of motion are with P the "canonical" Poisson tensor (5.30). As a system of quasi-linear equations, they can be cast in the form  where the charateristic matrix reads andg = g h ρ is the reduced gravity. As in the free surface case, the question of existence of Riemann invariants for this quasi-linear system can be easily answered by computing (by means of standard computer algebra programs) the Haantjes tensor H of the matrix A. A lengthy computation shows that H has 12 (out of 24) non-vanishing components, namely where r ij = ρ i − ρ j . As should be expected from the free surface case, even in the Boussinesq approximation the model for 3-layer fluid confined between rigid surfaces does not admit Riemann invariants. As can also be expected, this non-existence extends a fortiori to the general non-Boussinesq system, as well as to n-layered models with a rigid lid for n > 3. The implications on the structural properties of quasi-linear systems that do not admit Riemann invariants briefly discussed in Remark 3.1 naturally apply to all these cases as well.

Symmetric solutions
In the recent paper [22] the authors have focussed on the symmetric solutions defined by the equality of the upper (i = 1) and lower (i = 3) layer thicknesses, i.e., ζ 2 = h − ζ 1 , and the averaged horizontal velocities, u 1 = u 3 . In the Boussinesq approximation, our variables (σ 1 , σ 2 ) are actually proportional to the velocity shears, so that the symmetric solutions found in [22] are given by the relations A straightforward computation confirms that the submanifold defined by these relations is invariant under the flow (6.5) if and only if the relation is fulfilled among the density differences. In this case system (6.5) reduces to a system with 2 "degrees of freedom," parametrized, e.g., by the pair (ζ 2 ≡ ζ, σ 2 ≡ σ). The reduced "symmetric" equations of the motion are where we have defined ρ ∆ = ρ 3 − ρ 2 . These equations follow from the Hamiltonian functional with the "standard" Poisson tensor where the reference level for the potential energy in the Hamiltonian density is chosen at the midpoint of the channel. From equations (6.11) the characteristic velocities of the system can be read off immediately, since in this case ζ ∈ (0, h/2), the domain of hyperbolicity coincides with the rectangular region Alternatively, the Hamiltonian formulation shows that the simple map ζ → η, σ → hσ, ρ → ρ 2 , ρ ∆ → ρ 2 − ρ 1 turns the Hamiltonian in (6.12) and the operator in (6.13) into that of the two-layer Boussinesq system in a channel of half width h/2 reported in [13] (as also remarked in [22], though only reported explicitly from the motion equation viewpoint in the Ph.D. thesis [23]), thereby translating all the properties discussed therein to the present three-layer Boussinesq symmetric case. In fact, these symmetric solutions have a clear meaning in the Hamiltonian setting. Consider the involution J defined by Since its Jacobian is the Poisson tensor (5.30) is preserved by the involution. As far as the Hamiltonian (6.3) is concerned, the kinetic energy density is invariant under the involution (6.15), while the potential energy density, when ρ 3 − ρ 2 = ρ 2 − ρ 1 ≡ ρ ∆ , changes by a constant term added to the linear term H ∆ = −gρ ∆ ( ζ 1 + ζ 2 ). Since +∞ −∞ H ∆ dx is a Casimir function for the Poisson tensor (5.30), the equations of motion are left invariant by the involution (6.15). In other words, symmetric solutions correspond to the manifold of fixed points of the canonical involution (6.15), and H B,S is just the restriction of the Hamiltonian H B (up to a factor 1/2) to the space defined by relations (6.9) under conditions (6.10) for the density differences.

Summary and conclusions
In this paper we examined some structural properties of multi-layered incompressible Euler flows in the long-wave regime. In order to present explicit computations, we mainly focussed on the three-layer case. The aim was to point out similarities and differences with the two-layer setups which has received much attention in the literature (see, e.g., [19,9]).
At first, we briefly recalled how effective equations in one space dimension are obtained by layer-means and the hydrostatic approximation for the pressure in the dispersionless case. We then recalled some known facts about the free surface case, highlighting the natural Hamiltonian structure of the ensuing equations, and the non-existence of the Riemann invariants for the quasilinear resulting system, a fact proven in [14].
Next, we moved onto the main focus of the present paper, the case of stratified fluid in a vertical channel, a configuration that mathematically translates to the enforcement of the rigid lid upper constraint. The phenomenon of effective pressure differentials implying the "paradox" of nonconservation of the horizontal momentum, highlighted in [9] for the two-layer case, has been shown to persist for an n-layered situation, n ≥ 3, and in fact even be enhanced, for zero initial velocities, by scaling linearly with density differences, as opposed to quadratically as in the two layer case. A natural Hamiltonian structure on the configuration space of effective (in one space dimension) 3-layered fluid motions was derived by means of a geometrical reduction process from the full 2-dimensional Hamiltonian structure introduced in [4]. This allowed us to write the equations of motion in the form of a system of conservation laws, and to recover the correct conserved quantity (the impulse) associated with the translational Noether symmetry of the system. The Hamiltonian formulation led to a simple derivation of the Boussinesq limit of the motion equations. The (expected) non-existence of Riemann invariants in the rigid lid case was proved, and the geometrical meaning of the "symmetric" configurations recently studied in [22] is pointed out. The determination of invariant submanifolds of the full equations ensuing from the Hamiltonian (5.42), conditioned by suitable constraints on the density differences, and the study of the properties of such reduced systems is a non trivial question. Its analysis is left for future investigations.