1 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 Ovsyannikov (1979) to the more recent (Duchêne et al. 2016; Chesnokov et al. 2017; de Melo Viríssimo and Milewski 2019), through (Choi and Camassa 1996, 1999; Lvov and Tabak 2001; Chumakova et al. 2009a, b).

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 Choi and Camassa (1999) 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., Choi (2000), and Gavrilyuk et al. (1998); Barros (2006) for a discussion of a variational approach. Extension of these results, considering dispersive terms (as in (2.6)), can be found in Barros et al. (2007), Percival et al. (2008). In this respect, our main goal is to highlight the following facts: (1) the non-existence of Riemann invariants and (2) 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 \(\rho _1<\rho _2<\rho _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 Benjamin (1986) and thoroughly discussed in Camassa et al. (2012), Camassa et al. (2013) for the 2D Euler equations, can be detected, not unexpectedly, in the multiply stratified case, as also noted in de Melo Viríssimo and Milewski (2019). We notice that the added degree of freedom afforded by the three-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 non-vanishing time-variation of the horizontal momentum at the first order in the density differences \(\rho _{j+1}-\rho _{j}\) even with vanishing initial velocities. Hamiltonian aspects of such models have also been considered, notably in Benjamin and Bridges (1997), Craig and Sulem (1993), Craig et al. (2005), with an emphasis on the two-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 (1986) 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 two-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.

2 Sharply Stratified \({\varvec{n}}\)-Layered Euler Fluids

The incompressible Euler equations for the velocity field \({\mathbf {u}}=(u,w)\) and non-constant density \(\rho (x,z,t)\), in the presence of gravity \(-g{\mathbf {k}}\) and with horizontal fluid domain bottom at \(z=0\), are

$$\begin{aligned} \frac{D \rho }{D t}=0\,, \qquad \nabla \cdot {\mathbf {u}} =0\,, \qquad \frac{D (\rho {\mathbf {u}}) }{D t} + \nabla p + \rho g {\mathbf {k}}=0\,, \end{aligned}$$
(2.1)

with boundary conditions \({\mathbf {u}}(x=\pm \infty ,z,t)={\mathbf {0}}\) for all \(z\ge 0\) and \(w(x,0,t)=0\) for all \(x\in {{\mathbb {R}}}\). As usual, \(D/Dt=\partial /\partial t+{\mathbf {u}}\cdot \nabla \) 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 Wu (1981). 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 Wu (1981), it was established that for a density-stratified fluid, for each layer of thickness \(\eta _i=\zeta _{i-1}-\zeta _i\),

$$\begin{aligned} \int \limits _{\zeta _{i}}^{\zeta _{i-1}} \frac{D}{Dt} f\, \mathrm {d}z= \frac{\partial }{\partial t} \int \limits _{\zeta _{i}}^{\zeta _{i-1}} f \, \mathrm {d}z+ \frac{\partial }{\partial x} \int \limits _{\zeta _{i}}^{\zeta _{i-1}} u f \, \mathrm {d}z\, , \qquad i=1,\dots , n, \end{aligned}$$
(2.2)

where f(xzt) is any quantity and the elevations \(\zeta _\alpha \) are related with the thicknesses of the fluid’s layers via

$$\begin{aligned} \left\{ \begin{aligned}&\zeta _n\equiv 0\quad (\zeta _n\, \text { is the flat bottom})\, ,\\&\zeta _{i}=\sum _{k=0}^{n-i-1} \eta _{n-k},\quad i=1,\ldots , n-1\, . \end{aligned} \right. \end{aligned}$$
(2.3)

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 \(\eta _i\) and the layer mean velocities \({{\overline{u}}}_i\) defined as

$$\begin{aligned} {{\overline{u}}}_i ={\displaystyle {\frac{1}{\eta _i}{\int \limits _{\zeta _{i}}^{\zeta _{i-1}} u(x,z)\, \mathrm {d}z}}}\, . \end{aligned}$$
(2.4)

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 Choi and Camassa (1999), and extended the case of n layers in Choi (2000). Motions of typical wavelength L were considered under the assumptions that the ratios

$$\begin{aligned} \epsilon =\frac{h}{L}\simeq \frac{\eta _i}{L}\,,\quad i=1,\ldots , n\,, \end{aligned}$$
(2.5)

can be considered small. Here, h is the total height of the channel, while \(\eta _i\) is the thickness of the i-th fluid homogeneous layer.

By using the so-called columnar ansatz (see, e.g., Wu 1981), 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 \(\epsilon \), in Choi and Camassa (1999) it was shown that the \(2+1\)-dimensional Euler equations (2.1) reduce to the \(1+1\)-dimensional equations

$$\begin{aligned} {\eta _i}_t+({{\overline{u}}}_i\eta _i)_x=0\,,\qquad {{{\overline{u}}}_i}_t+{{{\overline{u}}}_i}{{{\overline{u}}}_i}_x +(-1)^ig {\eta _i}_x + \frac{P_x}{\rho _i} + D_i =0\, ,\qquad i=1, 2\,,\nonumber \\ \end{aligned}$$
(2.6)

where the \(D_i\) are dispersive terms, which at \(O(\epsilon ^2)\) read \(D_i=\frac{1}{3 \eta _i} [\eta _i^3 ({{{\overline{u}}}_i}_{xt} +{{{\overline{u}}}_i}{{{\overline{u}}}_i}_{xx}-({{{\overline{u}}}_i}_x)^2)]_x \), the \({{\overline{u}}}_{i}\)’s are the layer-mean velocities, \(\rho _i\) is the density of the i-th layer, and P(xt) is the interfacial pressure. These results were generalized in Choi (2000), Barros et al. (2020) 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 \(\zeta _{i}<z<\zeta _{i-1}\):

$$\begin{aligned} \begin{aligned} p_i(x,z)&=p(x,\zeta _{i}(x))- \rho _i g (z-\zeta _{i}) \\&= P^0- g \sum _{k=0}^{n-i-1}\rho _{n-k} \eta _{n-k} - \rho _{i} g (z-\zeta _{i}) \, , \qquad i=1,2,\dots ,n-1\, ,\\ p_n(x,z)&=P^0-\rho _n g z\, . \end{aligned}\nonumber \\ \end{aligned}$$
(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 \(\overline{{p_i}_x}\) (which, owing to the non-commutative property of exchanging derivative with layer averages, is in general different from the x-derivative of the averaged pressure \((\overline{p_i})_x\)). Since by (2.3) \(\zeta _{i}=\sum _{k=0}^{n-i-1} \eta _{n-k}\) for \(i=1,\ldots , n-1\) while \(\zeta _n\equiv 0\), equation (2.7) implies that

$$\begin{aligned} \begin{aligned} {{p_i}_x}=&P^0_x- g \sum _{k=0}^{n-i-1}\rho _{n-k} \eta _{{n-k}\, , x} + g \rho _{i} (\zeta _{i})_x \\ =&P^0_x- g \sum _{k=0}^{n-i-1}(\rho _{n-k}-\rho _{i}) \eta _{{n-k}\, , x} \, , \quad i=1,\dots ,n-1\, ,\\ { p_n}_x=&P^0_x\, . \end{aligned} \end{aligned}$$
(2.8)

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

$$\begin{aligned} \begin{aligned} \overline{{p_i}_x}\equiv&\frac{1}{\eta _i}\int \limits _{\zeta _{i}}^{\zeta _{i-1}} {p_i}_x \, \mathrm {d}z \, = {p_i}_x \\&= P^0_x- g \sum _{k=0}^{n-i-1}(\rho _{n-k}-\rho _{i}) \eta _{{n-k}\, , x} \, , \qquad i=1,\dots ,n-1\, ,\\ \overline{ p_n}_x=&P^0_x\, . \end{aligned} \end{aligned}$$
(2.9)

We thus obtain the set of 2n equations

$$\begin{aligned} \begin{aligned}&{\eta _i}_t+ (\eta _i {{\overline{u}}}_i)_x=0 \, , \qquad i=1,\dots ,n-1\,,\\&{{\overline{u}}_i}_t + {{\overline{u}}}_i\ {{{\overline{u}}}_i}_x+ \frac{P^0_x}{\rho _i} - g \sum _{k=0}^{n-i-1}\frac{\rho _{n-k}-\rho _i}{\rho _i} \eta _{{n-k}\, , x}=0\, , \qquad i=1,\dots ,n-1\,, \\&{\eta _n}_t+ (\eta _n{{\overline{u}}}_n)_x=0\,, \\&{{\overline{u}}_n}_t + {{\overline{u}}}_n {{{\overline{u}}}_i}_x+ \frac{P^0_x}{\rho _n}=0\, , \end{aligned} \end{aligned}$$
(2.10)

for the \(2n+1\) dependent variables \(\big (\eta _i, {{\overline{u}}}_i, P^0\big )\).

3 The Free-Surface Equations: The Case \(\varvec{n=2}\)

A standard way to close system (2.10) consists of considering the free-surface case (see, e.g., Ovsyannikov 1979; Duchêne et al. 2016; Chesnokov et al. 2017; Craig et al. 2005; Bona et al. 2008; Choi 2000). In this instance, the reference pressure \(P^0\) can be eliminated by observing that the pressure at the fluid’s domain upper boundary \(z=\zeta _0\) can be consistently set to vanish, as shown in Fig. 1. Thus, \(P^0\) can be expressed, thanks to eq. (2.7) for \(i=1\), as

$$\begin{aligned} P^0=g\, \sum _{k=1}^n \rho _k\eta _k\,. \end{aligned}$$
(3.1)

So, for \(i=1,\ldots ,n\), the pressures \(p_{i}\) become (linear) functions of the thicknesses \(\eta _{i}\) and the system closes. An example of this procedure in the two-layer case is as follows (the generalization to arbitrary n being straightforward). The pressure in the lowest layer is \(p_2(x,z)=P^0-g \rho _2 z\). In the upper layer, we have

$$\begin{aligned} p_1(x,z)=p_2(x,\zeta _1)-g\rho _1(z-\zeta _1)=P^0-g\rho _2\eta _2-g\rho _1\eta _1, \end{aligned}$$
(3.2)

since \(\zeta _1=\eta _2\) and \(\zeta _0-\zeta _1=\eta _1\). Requiring the vanishing of the pressure at the free boundary \(z=\zeta _0=\eta _1+\eta _2\) yields, as expected, \(P^0=g(\rho _1\eta _1+\rho _2\eta _2)\), and so

$$\begin{aligned} \overline{{p_2}_x}=P^0_x=g(\rho _1{\eta _1}_x+\rho _2{\eta _2}_x),\quad \overline{{p_1}_x}=g(\rho _1{\eta _1}_x+\rho _1{\eta _2}_x)\, . \end{aligned}$$
(3.3)

Thus, the equations of the motion deduced from (2.10) are

$$\begin{aligned} \begin{array}{lcl} {\eta _1}_t+ (\eta _ 1{\overline{u}}_1)_x=0\,, &{}\quad &{} {\eta _2}_t+ (\eta _ 2{\overline{u}}_2)_x=0\,,\\ &{}&{}\\ {{\overline{u}}_1}_t + {\overline{u}}_1\ {{\overline{u}}_1}_x +g\left( {\eta _1}_x+{\eta _2}_x\right) =0\,, &{}\quad &{} {{\overline{u}}_2}_t + {\overline{u}}_2\ {{\overline{u}}_2}_x +g \left( {\displaystyle {\frac{\rho _1}{\rho _2}{\eta _1}_x+{\eta _2}_x}}\right) =0\,. \end{array} \end{aligned}$$
(3.4)

These equations can be written in the standard matrix form for quasi-linear systems as

$$\begin{aligned} \left( \begin{array}{c} {\eta _1}_t\\ {\eta _2}_t\\ {{{\overline{u}}}_1}_t\\ {{{\overline{u}}}_2}_t \end{array} \right) +{\mathbf {A_2}} \left( \begin{array}{c} {\eta _1}_x\\ {\eta _2}_x\\ {{{\overline{u}}}_1}_x\\ {{{\overline{u}}}_2}_x \end{array} \right) =0\,, \end{aligned}$$
(3.5)

where the characteristic matrix reads

$$\begin{aligned} {\mathbf {A_2}}=\left( \begin{array}{cccc} {{\overline{u}}}_{{1}}&{}0&{}\eta _{{1}}&{}0\\ 0&{}{{\overline{u}}}_{{2}}&{}0&{}\eta _{{2}}\\ g&{}g&{}{{\overline{u}}}_{{1}}&{}0\\ {\displaystyle {{\frac{g\rho _{{1}}}{\rho _{{2}}}}}}&{}g&{}0&{}{{\overline{u}}}_{{2}}\end{array} \right) \,. \end{aligned}$$
(3.6)

The question of the existence of Riemann invariants for this quasi-linear system can be easily solved by computing the so-called Haantjes tensor \(\mathcal {H}\) of the matrix \({\mathbf {A_2}}\) (Haantjes 1955; Bogoyavlenskij 1993; Ferapontov and Marshall 2007), 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 Ovsyannikov (1979) and proved in Chesnokov et al. (2017), since it can be easily checked that the Haantjes tensor of \({\mathbf {A_2}}\) has three non-vanishing components,

$$\begin{aligned} {\mathcal {H}^1}_{1\, 2}={\mathcal {H}^3}_{2\, 3} ={\frac{\eta _{{1}}{g}^{2} \left( \rho _{{1}}-\rho _{{2}} \right) }{\rho _ {{2}}}} \quad \text {and } {\mathcal {H}^4}_{1\, 3}= -{\frac{\eta _{{1}}{g}^{2}\rho _{{1}} \left( \rho _{{1}}-\rho _{{2}} \right) }{{\rho _{{2}}}^{2}}} \, . \end{aligned}$$
(3.7)
Fig. 1
figure 1

Free-surface, two-layer fluid setup and relevant notation: \(P(x,t)=0\) is the free-surface (air) pressure, \(\zeta _0\) and \(\zeta _1\) are the free surface and interface locations, and \(\eta _1\) and \(\eta _2\) are the layer thicknesses of the fluids with densities \(\rho _1\) and \(\rho _2\), respectively

System (3.4) does admit a natural Hamiltonian formulation. To show this, we pass from averaged velocities \({{\overline{u}}}_i\) to the averaged momentum coordinates \({{\overline{\mu }}}_i\) defined by \({{\overline{\mu }}}_i=\rho _i{{\overline{u}}}_i\), with \(i=1,2\). The evolution equations translate into

$$\begin{aligned} \begin{array}{lcl} {\eta _1}_t+{\displaystyle {\frac{1}{\rho _1}}} (\eta _ 1{{\overline{\mu }}}_1)_x=0\,, &{}\quad &{} {\eta _2}_t+{\displaystyle {\frac{1}{\rho _2}}} (\eta _ 2{{\overline{\mu }}}_2)_x=0\,, \\ {{{\overline{\mu }}}_1}_t + {\displaystyle {\frac{1}{\rho _1}}} {{\overline{\mu }}}_1{{{\overline{\mu }}}_1}_x +g\left( \rho _1{\eta _1}_x+\rho _1{\eta _2}_x\right) =0\,, &{}\quad &{} {{{\overline{\mu }}}_2}_t + {\displaystyle {\frac{1}{\rho _2}}} {{\overline{\mu }}}_2{{{\overline{\mu }}}_2}_x\\ &{}&{}+g\left( \rho _1{\eta _1}_x+\rho _2{\eta _2}_x\right) =0\,. \end{array} \end{aligned}$$
(3.8)

A direct inspection shows that (3.8) admits the Hamiltonian formulation

$$\begin{aligned} \left( \begin{array}{c} {\eta _1}_t\\ {\eta _2}_t\\ {{{\overline{\mu }}}_1}_t\\ {{{\overline{\mu }}}_2}_t \end{array}\right) =- \left( \begin{array}{cccc} 0&{}0&{}\partial _x&{}0\\ 0&{}0&{}0&{}\partial _x\\ \partial _x&{}0&{}0&{}0\\ 0&{}\partial _x&{}0&{}0 \end{array}\right) \left( \begin{array}{c} \delta _{\eta _{1}} H\\ \delta _{\eta _{2}}H\\ \delta _{{{\overline{\mu }}}_{1}}H\\ \delta _{{{\overline{\mu }}}_{2} }H \end{array}\right) \,, \end{aligned}$$
(3.9)

where the vector \(\left( \delta _{\eta _{1}} H, \delta _{\eta _{2}}H, \delta _{{{\overline{\mu }}}_{1}}H, \delta _{{{\overline{\mu }}}_{2}} H\right) \) is the differential of the effective energy functional

$$\begin{aligned} H= \int \limits _{-\infty }^{+\infty } \, \frac{1}{2}\left( \frac{\eta _1}{\rho _1}{{\overline{\mu }}}_1^2+\frac{\eta _2}{\rho _2}{{\overline{\mu }}}_2^2 +g (\rho _1 \eta _1^2+2\rho _1\eta _1\eta _2+\rho _2\eta _2^2)\right) \mathrm {d}x\, . \end{aligned}$$
(3.10)

Remark 3.1

The two-layer system (3.8) is a non-diagonalizable Hamiltonian system of conservation laws in Tsarev’s sense (Tsarev 1985, 1991). 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 Tsarev (1985, 1991), showing that six quantities are guaranteed to be conserved by strong solutions: the Hamiltonian functional (3.10), \(\int _{-\infty }^{+\infty } \eta _j\, \mathrm {d}x, \, \int _{-\infty }^{+\infty } {{\overline{\mu }}}_j\, \mathrm {d}x\), \(j=1,2\) (these are Casimir functionals of the Poisson brackets), and the total horizontal momentum \(\Pi ^{(x)}=\int _{-\infty }^{+\infty } (\eta _1{{\overline{\mu }}}_1+\eta _2{{\overline{\mu }}}_2)\, \mathrm {d}x\), this being the generator of the x-translational symmetry in the free-surface case. Furthermore, in Montgomery and Moodie (2001) and Barros (2006), 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 Ferapontov (1994), about the complete integrability of quasi-linear 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 quasi-linear 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 \(\varvec{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

$$\begin{aligned} \eta _1+\eta _2+\dots +\eta _n = h\, , \end{aligned}$$
(4.1)

where h is the total height of the vertical fluid “channel”. The sum of all layer-mass conservation equations \({\eta _k}_t+(\eta _k{{\overline{u}}}_k)_x=0\) yields

$$\begin{aligned} \left( \eta _1 {\overline{u}}_1 + \eta _2 {\overline{u}}_2 + \dots +\eta _n {\overline{u}}_n\right) _x=0 . \end{aligned}$$
(4.2)

Taking into account the boundary conditions

$$\begin{aligned} \lim _{x \rightarrow \infty } {\overline{u}}_i =0\,, \qquad i=1,2,\dots , n\,, \end{aligned}$$
(4.3)

equation (4.2) translates into the total flux independence of x

$$\begin{aligned} \eta _1 {\overline{u}}_1 + \eta _2 {\overline{u}}_2 + \dots +\eta _n {\overline{u}}_n=0\, , \end{aligned}$$
(4.4)

that can be termed the dynamical constraint.

The system of equations (2.10) written for the momentum densities \(\eta _k{{\overline{u}}}_k\) is

$$\begin{aligned} ({\eta _i {\overline{u}}_i})_t + (\eta _i {\overline{u}}_i^2)_x+ \frac{P^0_x \eta _i }{\rho _i} - g \sum _{k=0}^{n-i-1}\frac{\rho _{n-k}-\rho _i}{\rho _i} ({\eta _{n-k}})_x \eta _i=0\,, \qquad i=1,2,\dots ,n\,.\nonumber \\ \end{aligned}$$
(4.5)

The sum of all these equations gives the expression of the horizontal gradient of the pressure at the bottom of the channel,

$$\begin{aligned} P^0_x=-\frac{\sum _{i=1}^n(\eta _i {\overline{u}}_i^2)_x+g \sum _{i=1}^{n} \sum _{k=0}^{n-i-1}\frac{\rho _{n-k}-\rho _i}{\rho _i} ({\eta _{n-k}})_x \eta _i}{\sum _{i=1}^{n}\eta _i /\rho _i}\,. \end{aligned}$$
(4.6)

This explicit form for \(P^0\) together with the geometric and dynamical constraints (4.1),(4.4) guarantees the closure of (2.10) as a system of \(2n-2\) equations for the \(2n-2\) variables, say, \((\eta _\alpha , {{\overline{u}}}_\alpha )\) with \(\alpha =2,\ldots , 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 two 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 three layers.

4.1 Three-Layer Rigid Lid Case

Our geometrical setting is shown in Fig. 2. From (2.7), it is clear that the explicit expression for the hydrostatic pressures in the layers is

$$\begin{aligned} \begin{aligned} p_3=&P^0-g \rho _3 z \\ p_2=&P^0-g \rho _3 \zeta _2 - g \rho _2 (z-\zeta _2)=P^0\\&-g (\rho _3-\rho _2) \zeta _2 - g \rho _2 z \\ p_1=&P^0-g \rho _3 \zeta _2 - g \rho _2 (\zeta _1-\zeta _2)- g \rho _1 (z-\zeta _1)\\&=P^0-g (\rho _3-\rho _2) \zeta _2-g (\rho _2-\rho _1) \zeta _1 - g \rho _1 z \end{aligned} \end{aligned}$$
(4.7)

and therefore

$$\begin{aligned}&{p_3}_x=P^0_x\,, \quad {p_2}_x=P^0_x-g (\rho _3-\rho _2) {\zeta _2}_x\, , \nonumber \\&\quad \quad {p_1}_x=P^0_x-g (\rho _3-\rho _2) {\zeta _2}_x-g (\rho _2-\rho _1) {\zeta _1}_x \, . \end{aligned}$$
(4.8)

The mass conservation equations read

$$\begin{aligned}&{\zeta _2}_t+({{\overline{u}}}_3 \zeta _2)_x=0\, , \quad {(\zeta _1-\zeta _2)}_t+({{\overline{u}}}_2 (\zeta _1-\zeta _2))_x=0\, ,\nonumber \\&\quad -{\zeta _1}_t+({{\overline{u}}}_1 (h- \zeta _1))_x=0\, , \quad \end{aligned}$$
(4.9)

which imply, together with the vanishing asymptotic conditions on the velocities, the constraint

$$\begin{aligned} {{\overline{u}}}_3 \zeta _2 + {{\overline{u}}}_2 (\zeta _1-\zeta _2)+{{\overline{u}}}_1 (h- \zeta _1)=0\, , \end{aligned}$$
(4.10)

as seen in (4.4), since \(\eta _3=\zeta _2\), \(\eta _2=\zeta _1-\zeta _2\), \(\eta _1=h-\zeta _1\). The Euler equations lead to

$$\begin{aligned} \begin{aligned}&\rho _3 {{{\overline{u}}}_3}_t+ \rho _3 {{\overline{u}}}_3 {{{\overline{u}}}_3}_x + P^0_x=0 \, ,\\&\rho _2 {{{\overline{u}}}_2}_t+ \rho _2 {{\overline{u}}}_2 {{{\overline{u}}}_2}_x + P^0_x-g (\rho _3-\rho _2) {\zeta _2}_x=0 \, , \\&\rho _1 {{{\overline{u}}}_1}_t+ \rho _1 {{\overline{u}}}_1 {{{\overline{u}}}_1}_x + P^0_x-g (\rho _3-\rho _2) {\zeta _2}_x-g (\rho _2-\rho _1) {\zeta _1}_x=0 \, .\\ \end{aligned} \end{aligned}$$
(4.11)

The mass conservation equations and the Euler equations imply the following set of evolution equations (see (4.5)):

$$\begin{aligned} \begin{aligned}&(\eta _3 {{{\overline{u}}}_3})_t+ (\eta _3 {{{\overline{u}}}_3}^2 )_x + P^0_x \, \frac{\eta _3}{\rho _3}=0 \, ,\\&(\eta _2 {{{\overline{u}}}_2})_t+ (\eta _2 {{{\overline{u}}}_2}^2)_x + \left( P^0_x-g (\rho _3-\rho _2) {\zeta _2}_x\right) \frac{\eta _2}{\rho _2}=0 \, , \\&(\eta _1 {{{\overline{u}}}_1})_t+ (\eta _1 {{{\overline{u}}}_1}^2)_x + \left( P^0_x-g (\rho _3-\rho _2) {\zeta _2}_x-g (\rho _2-\rho _1) {\zeta _1}_x\right) \frac{\eta _1}{\rho _1}=0 \, .\\ \end{aligned} \end{aligned}$$
(4.12)

By taking the sum of these three evolution equations and by using the constraint \(\eta _3 {{\overline{u}}}_3+\eta _2 {{\overline{u}}}_2+\eta _1 {{\overline{u}}}_1=0\), we can solve for the x-derivative of the pressure \(P^0\) as (see (4.6))

$$\begin{aligned} P^0_x{=}\frac{-(\eta _3 {{{\overline{u}}}_3}^2+\eta _2 {{{\overline{u}}}_2}^2+\eta _1 {{{\overline{u}}}_1}^2 )_x+ g\left( \frac{\rho _3-\rho _2}{\rho _2} {\eta _3}_x \eta _2 {+} \frac{\rho _3-\rho _2}{\rho _1} {\eta _3}_x \eta _1 -\frac{\rho _2-\rho _1}{\rho _1} {\eta _1}_x \eta _1 \right) }{ \frac{\eta _3}{\rho _3}+\frac{\eta _2}{\rho _2}+\frac{\eta _1}{\rho _1}}\,,\nonumber \\ \end{aligned}$$
(4.13)

where we used again the fact that \(\zeta _2=\eta _3\) and \(\zeta _1=h-\eta _1\). Substituting \(\eta _1=h-\eta _2-\eta _3\) yields

$$\begin{aligned} P^0_x=\frac{-(\eta _3 {{{\overline{u}}}_3}^2+\eta _2 {{{\overline{u}}}_2}^2+\eta _1 {{{\overline{u}}}_1}^2 )_x+ g\left( \frac{\rho _3(\rho _1-\rho _2)}{\rho _1\rho _2} {\eta _3}_x \eta _2 + \frac{\rho _3-\rho _1}{\rho _1} {\eta _3}_x (h-\eta _3) +\frac{\rho _2-\rho _1}{\rho _1} {\eta _2}_x (h-\eta _2-\eta _3) \right) }{ \frac{1}{\rho _1}\left( h + \frac{\rho _1-\rho _2}{\rho _2}\eta _2+\frac{\rho _1-\rho _3}{\rho _3}\eta _3\right) }\,.\nonumber \\ \end{aligned}$$
(4.14)

Hence, we have the following equations of motion for \((\eta _2,\eta _3,{{\overline{u}}}_2,{{\overline{u}}}_3)\),

$$\begin{aligned} \begin{array}{lcl} {\eta _2}_t+ (\eta _ 2\overline{u_2})_x=0\,, &{}\quad &{} {\eta _3}_t+ (\eta _ 3\overline{u_3})_x=0\,, \\ \rho _2 {{{\overline{u}}}_2}_t+ \rho _2 {{\overline{u}}}_2 {{{\overline{u}}}_2}_x + P^0_x-g (\rho _3-\rho _2) {\eta _3}_x=0\,, &{}\quad &{} \rho _3 {{{\overline{u}}}_3}_t+ \rho _3 {{\overline{u}}}_3 {{{\overline{u}}}_3}_x + P^0_x=0\,, \end{array}\nonumber \\ \end{aligned}$$
(4.15)

where \(P^0_x\) is given by (4.14) and \({{\overline{u}}}_1=-(\eta _3 {{\overline{u}}}_3+\eta _2 {{\overline{u}}}_2)/(h-\eta _2-\eta _3)\).

Fig. 2
figure 2

Three-layer rigid lid setup and relevant notation: \(P^0\) is the bottom pressure, \(\zeta _i\) and \(\eta _i\), \(i=1,2\), are the interface locations and layer thicknesses, respectively.

4.2 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 in Benjamin (1986), was highlighted and substantiated in Camassa et al. (2012), Camassa et al. (2013). 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 setup, whenever hydrostatic conditions apply at infinity, the violation of momentum conservation is proportional (up to terms that arise from possibly different configurations at \(x=\pm \infty \)) to the difference \(P_{\Delta }\) 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 de Melo Viríssimo and Milewski 2019) can be detected looking at equation (4.13). Indeed, its right-hand side in not a total x-derivatives, and so \(\int \limits _{-\infty }^{+\infty } P^0_x\, \mathrm {d}x\) need not vanish even for configurations with asymptotically flat interfaces \(\zeta _i\) at the far-ends of the channel. A very simple example of such a configuration with non-vanishing \(P_{\Delta }\) is given in Fig. 3 by splicing together arcs of parabolae for the interface locations. We choose, along with \({{\overline{u}}}_i=0\) for \(i=1,2,3\),

$$\begin{aligned} \eta _3= & {} {\left\{ \begin{array}{ll} \frac{h}{4} (1-x^2)+\frac{h}{5}&{}\text{ if } -1<{x}<1\\ \frac{h}{5}&{}\text{ otherwise }\end{array}\right. }\,,\nonumber \\ \eta _2= & {} {\left\{ \begin{array}{ll} \frac{h}{3} (x^2-2\, x)+\frac{2h}{5}&{}\text{ if } 0<{x}<2\\ \frac{2h}{5}&{}\text{ otherwise }\end{array}\right. }\, . \end{aligned}$$
(4.16)

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 \(\rho _2-\rho _1\) and \(\rho _3-\rho _2\), even with initial null velocities, as opposed to two-layer, rigid lid configurations, where non-conservation of the horizontal total momentum is, with null-velocities, quadratic in the density difference \(\rho _2-\rho _1\) (see Camassa et al. 2013). Indeed, equation (4.14) reduces to

$$\begin{aligned} P^0_x=\frac{ g\left( \frac{\rho _3(\rho _1-\rho _2)}{\rho _1\rho _2} {\eta _3}_x \eta _2 + \frac{\rho _3-\rho _1}{\rho _1} {\eta _3}_x (h-\eta _3) +\frac{\rho _2-\rho _1}{\rho _1} {\eta _2}_x (h-\eta _2-\eta _3) \right) }{ \frac{1}{\rho _1}\left( h + \frac{\rho _1-\rho _2}{\rho _2}\eta _2+\frac{\rho _1-\rho _3}{\rho _3}\eta _3\right) }\,.\nonumber \\ \end{aligned}$$
(4.17)

The numerator is equal, up to total derivatives, to \( g (\rho _2-\rho _1) (\rho _3-\rho _2){\eta _3}_x\eta _2/(\rho _1\rho _2)\), thus revealing the leading order dependence of \(P_\Delta \) on the density differences for general initial layer thicknesses \(\eta _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, quasi-linear equations governing stratified flows are generically systems of mixed type (Chumakova et al. 2009b; Boonkasame and Milewski 2011), 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 \(\rho _2/\rho _3=3/4, \rho _1/\rho _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 nonzero 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\simeq 2.8\) (with the non-dimensional parameters as in the example of Fig. 3).

Fig. 3
figure 3

A three-layer configuration with a non-vanishing initial pressure imbalance. Setting \(\rho _1=0.5\), \(\rho _2=0.75\), \(\rho _3=1\), and \(g=h=1\), gives \(P_{\Delta }\simeq -0.0037395\)

5 The Three-Layer Rigid Lid Case: The Hamiltonian Structure

In this section, we extend our study of the three-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 two-layer case (as well as in the free-surface three-layer case, see Appendix), the resulting equations, when written in the momentum coordinates \({{\overline{\mu }}}_k=\rho _k{{\overline{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 Benjamin (1986), by extending to the multiple-layer case a technique introduced in Camassa et al. (2017).

5.1 The 2D Benjamin Model for Heterogeneous Fluids in a Channel

Benjamin (1986) proposed and discussed a setup for the Hamiltonian formulation of an incompressible stratified Euler system in two 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 \(\rho \) together with the “weighted vorticity" \(\Sigma \) defined by

$$\begin{aligned} \Sigma =\nabla \times (\rho {\varvec{u}})=(\rho w)_x-(\rho u)_z. \end{aligned}$$
(5.1)

The equations of motion for these two fields, ensuing from (2.1), are

$$\begin{aligned} \begin{array}{l} \rho _t+u\rho _x+w\rho _z =0\\ \Sigma _t+u\Sigma _x +w\Sigma _z +\rho _x\big (gz-\frac{1}{2}(u^2+w^2)\big )_z+\frac{1}{2}\rho _z\big (u^2+w^2\big )_x=0\ . \end{array} \end{aligned}$$
(5.2)

They can be written in the form

$$\begin{aligned} {\rho _t}=-\left[ \rho , {\displaystyle {\displaystyle {\frac{\delta H}{\delta \Sigma }}}}\right] \, , \qquad \Sigma _t= -\left[ \rho , {\displaystyle {\displaystyle {\frac{\delta H}{\delta \rho }}}}\right] -\left[ \Sigma , {\displaystyle {\displaystyle {\frac{\delta H}{\delta \Sigma }}}}\right] \, , \end{aligned}$$
(5.3)

where by definition, \([A, B] \equiv A_xB_z-A_zB_x\), and the functional

$$\begin{aligned} H= {\displaystyle {\int \limits _\mathcal {D} \rho \left( \frac{1}{2}|{\varvec{u}} |^2+g z\right) \,\mathrm{d}x\,\mathrm{d}z }} \end{aligned}$$
(5.4)

is simply given by the sum of the kinetic and potential energy, \(\mathcal {D}\) being the fluid domain. The most relevant feature of this coordinate choice is that \((\rho ,\Sigma )\) 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., Zakharov et al. 1985).

As shown by Benjamin, Eqs. (5.3) are a Hamiltonian system with respect to a Lie-theoretic Hamiltonian structure; that is, they can be written as

$$\begin{aligned} \rho _t=\{\rho , H\} ,\qquad \Sigma _t=\{\Sigma , H\} \end{aligned}$$

for the Poisson bracket defined by the Hamiltonian operator

$$\begin{aligned} J_B=- \left( \begin{array}{cc} 0 &{} \rho _x \partial _z -\rho _z \partial _x \\ \rho _x \partial _z -\rho _z \partial _x &{} \Sigma _x \partial _z -\Sigma _z \partial _x \end{array} \right) . \end{aligned}$$
(5.5)

5.2 The Reduction Process

By means of the Heaviside \(\theta \) and Dirac \(\delta \) generalized functions, a three-layer fluid configuration can be introduced with constant densities \(\rho _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 \(\zeta _1\) and \(\zeta _2\). The global density and velocity variables can be written as

$$\begin{aligned} \begin{aligned}&\rho (x,z)=\rho _{3}+(\rho _{2}-\rho _{3})\theta (z-\zeta _{2})+(\rho _{1}-\rho _{2})\theta (z-\zeta _{1})\\&u (x,z)=u_{3}+(u_{2}-u_{3})\theta (z-\zeta _{2})+(u_{1}-u_{2})\theta (z-\zeta _{1}) \\&w (x,z)=w_{3}+(w_{2}-w_{3})\theta (z-\zeta _{2})+(w_{1}-w_{2})\theta (z-\zeta _{1})\, . \end{aligned} \end{aligned}$$
(5.6)

Thus, the density-weighted vorticity \(\Sigma =(\rho w)_{x}-(\rho u)_{z}\) is

$$\begin{aligned} \begin{aligned} \Sigma =&\rho _3({w_3}_x-{u_3}_z)+\theta (z-\zeta _{2})\left( \rho _{2}{w_2}_x-\rho _{2}{u_2}_z+\rho _{3}{u_3}_z-\rho _{3}{w_3}_x\right) \\&+\theta (z-\zeta _{1})\left( \rho _{1}{w_1}_x-\rho _{1}{u_1}_z+\rho _{2}{u_2}_z-\rho _{2}{w_2}_x \right) \\&+\delta (z-\zeta _{2})\left( (\rho _{3}w_{3}-\rho _{2}w_{2}){\zeta _2}_x+(\rho _{3}u_{3}-\rho _{2}u_{2}) \right) \\&+\delta (z-\zeta _{1}) \left( (\rho _{2}w_{2}-\rho _{1}w_{1}){\zeta _1}_x+(\rho _{2}u_{2}-\rho _{1}u_{1}) \right) \\ =&\rho _{3}\Omega _{3}+\theta (z-\zeta _{2})(\rho _{2}\Omega _{2}-\rho _{3}\Omega _{3})+\theta (z-\zeta _{1})(\rho _{1}\Omega _{1}-\rho _{2}\Omega _{2})\\&+\left( (\rho _{3}u_{3}-\rho _{2}u_{2})+(\rho _{3}w_{3}-\rho _{2}w_{2}){\zeta _2}_x\right) \delta (z-\zeta _{2}) \\&+\left( (\rho _{2}u_{2}-\rho _{1}u_{1})+(\rho _{2}w_{2}-\rho _{1}w_{1}){\zeta _1}_x\right) \delta (z-\zeta _{1})\,, \end{aligned} \end{aligned}$$

where \(\Omega _{i}={w_i}_x-{u_i}_z\) for \(i= 1,2,3\). Next, we assume the motion in each layer to be irrotational, so that \(\Omega _{i}=0 \) for all \( i=1,2,3\). Therefore, the density-weighted vorticity acquires the form

$$\begin{aligned} \Sigma =&\left( (\rho _{3}u_{3}-\rho _{2}u_{2})+(\rho _{3}w_{3}-\rho _{2}w_{2}){\zeta _2}_x \right) \delta (z-\zeta _{2})\nonumber \\&+\left( (\rho _{2}u_{2}-\rho _{1}u_{1})+(\rho _{2}w_{2}-\rho _{1}w_{1})\zeta _{1x} \right) \delta (z-\zeta _{1})\,. \end{aligned}$$
(5.7)

In the long-wave asymptotics described in Sect. 2, we have, at the leading order in the long-wave expansion parameter \(\epsilon = h/L\ll 1\),

$$\begin{aligned} u_{i} \sim {{{\overline{u}}}_{i}}\,,\qquad w_{i} \sim 0\,, \end{aligned}$$

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

$$\begin{aligned} \begin{aligned} \rho (x,z)&=\rho _{3}+(\rho _{2}-\rho _{3})\theta (z-\zeta _{2})+(\rho _{1}-\rho _{2})\theta (z-\zeta _{1})\\ \Sigma (x,z)&=\left( \rho _{3}{{{\overline{u}}}_{3}}-\rho _{2}{{{\overline{u}}}_{2}} \right) \delta (z-\zeta _{2}) +\left( \rho _{2}{{{\overline{u}}}_{2}}-\rho _{1}{{{\overline{u}}}_{1}} \right) \delta (z-\zeta _{1}) \, . \end{aligned} \end{aligned}$$
(5.8)

The x and z-derivative of the Benjamin’s variables given by equations (5.8) are generalized functions supported at the surfaces \(\{z=\zeta _1\}\cup \{z=\zeta _2\}\) and are computed as

$$\begin{aligned} \begin{aligned} \rho _{x}=&-(\rho _{2}-\rho _{3}){\zeta _2}_x \delta (z-\zeta _{2})-(\rho _{1}-\rho _{2}){\zeta _1}_x \delta (z-\zeta _{1}) \\ \rho _{z}=&(\rho _{2}-\rho _{3})\delta (z-\zeta _{2})+(\rho _{1}-\rho _{2})\delta (z-\zeta _{1})\, , \end{aligned} \end{aligned}$$
(5.9)

and

$$\begin{aligned} \begin{aligned} \Sigma _{z}=&(\rho _{3}{{{\overline{u}}}_{3}}-\rho _{2}{{{\overline{u}}}_{2}}) \delta '(z-\zeta _{2}) +(\rho _{2}{{{\overline{u}}}_{2}}-\rho _{1}{{{\overline{u}}}_{1}})\delta '(z-\zeta _{1}) \\ \Sigma _{x}=&-(\rho _{3}{{{\overline{u}}}_{3}}-\rho _{2}{{{\overline{u}}}_{2}})\delta '(z-\zeta _{2}){\zeta _2}_x -(\rho _{2}{{{\overline{u}}}_{2}}-\rho _{1}{{{\overline{u}}}_{1}})\delta '(z-\zeta _{1}){\zeta _1}_x\\&+\left( \rho _{3}{{{{\overline{u}}}_3}_x}-\rho _{2}{{{{\overline{u}}}_2}_x} \right) \delta (z-\zeta _{2}) +\left( \rho _{2}{{{{\overline{u}}}_2}_x}-\rho _{1}{{{{\overline{u}}}_1}_x} \right) \delta (z-\zeta _{1})\, . \end{aligned} \end{aligned}$$
(5.10)

To invert the map (5.8), we can integrate along the vertical direction z. Contrary to the two-layer case described in Camassa et al. (2017), 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 \(\mathcal {M}\) (parameterized by \(\rho (x,z)\) and \(\Sigma (x,z)\)) by the space \(\mathcal {F}\) of isopycnals \(z=f(x)\) and trivially extend by 0 the Poisson structure \(J_B\) given by (5.5).

As shown in Camassa et al. (2014), the operator (5.5) provides a proper Poisson structure under the condition that the density \(\rho \) be constant at the physical boundaries \(z=0\) and \(z=h\). Once equipped with the isopycnal \(z=f(x)\), we can define a projection \(\pi : \widetilde{\mathcal {M}}\equiv \mathcal {M}\times \mathcal {F} \rightarrow (C^\infty ({{\mathbb {R}}}))^4\) by means of

$$\begin{aligned} \begin{aligned}&\pi \left( \rho (x,z), \Sigma (x,z), f \right) =\left( \xi _{1},\xi _2,\tau _1,\tau _2\right) \\&\quad = \left( \int \limits _{0}^{h} (\rho (x,z)-\rho _{max} )\, \mathrm {d}z, \int \limits _{0}^{f} (\rho (x,z) -\rho _{max} )\, \mathrm {d}z, \int \limits _{0}^{h} \Sigma (x,z) \mathrm {d}z,\int \limits _{0}^{f} \Sigma (x,z) \mathrm {d}z\right) \, . \end{aligned}\nonumber \\ \end{aligned}$$
(5.11)

We choose to include the manifold of three-layer fluid configurations in the space \(\mathcal {M}\times \mathcal {F}\) by adding to (5.8) the following choice for the isopycnal f:

$$\begin{aligned} f={\displaystyle {{\frac{\zeta _{1}+\zeta _{2}}{2}}}}\, , \end{aligned}$$
(5.12)

hereafter denoted by \( {{\overline{\zeta }}}\). Equations (5.8) and (5.12) define a submanifold \(\mathcal {I}\) of \(\widetilde{\mathcal {M}}\), in one-to-one correspondence with the manifold \(\mathcal {S}\) of three-layer fluid configurations considered in Sect. 4.1, as it is parameterized by four functions of the horizontal coordinate x, and this set is equivalent to the 4-tuple \((\eta _1,\eta _2,{{\overline{u}}}_1,{{\overline{u}}}_2)\). For further reference, we explicitly remark that on \( \mathcal {I}\) the projection \(\pi \) is defined by the four-tuple

$$\begin{aligned}&\pi \left( \rho (x,z), \Sigma (x,z), {{\overline{\zeta }}}\right) =\left( \int \limits _{0}^{h} (\rho (x,z)-\rho _3 )\, \mathrm {d}z, \int \limits _{0}^{{{\overline{\zeta }}}} (\rho (x,z) \right. \nonumber \\&\left. \quad -\rho _3 )\, \mathrm {d}z, \int \limits _{0}^{h} \Sigma (x,z) \mathrm {d}z,\int \limits _{0}^{{{\overline{\zeta }}}} \Sigma (x,z) \mathrm {d}z\right) \, \, . \end{aligned}$$
(5.13)

Further, as we shall see below, the \(\dot{{\bar{\zeta }}}\) component of the tangent vector \(({\dot{\rho }},{\dot{\Sigma }},\dot{{\bar{\zeta }}})\) can be linearly expressed in terms of \({\dot{\rho }}\).

To obtain a Hamiltonian structure on the manifold \(\mathcal {S}\) by reducing Benjamin’s parent structure (5.5), we perform the following steps:

  1. 1.

    Starting from a 1-form on the manifold \(\mathcal {S}\), represented by the four-tuple \((\alpha _S^{(1)}, \alpha _S^{(2)},\alpha _S^{(3)},\alpha _S^{(4)})\), we construct its pull-back to \(\mathcal {I}\), that is, the 1-form \(\varvec{\alpha }_M=(\alpha _M^{(1)},\alpha _M^{(2)}, 0)\) at \(T^*{\widetilde{M}}\vert _\mathcal {I}\) satisfying the relation

    $$\begin{aligned} \int \limits _{-\infty }^{+\infty }\int \limits _0^h( \alpha _{M}^{(1)} {\dot{\rho }}+\alpha _{M}^{(2)}{\dot{\Sigma }})\, \mathrm {d}x\, \mathrm {d}z = \int \limits _{-\infty }^{+\infty } \sum _{k=1}^4\alpha _{S}^{(k)} \cdot \left( \pi _{*} ({\dot{\rho }},{\dot{\Sigma }})\right) ^k \, \mathrm {d}x\>\>, \end{aligned}$$
    (5.14)

    where \(\pi _*\) is the tangent map to (5.13).

  2. 2.

    We apply Benjamin’s operator (5.5) to the lifted one form \(\varvec{\alpha }_M\) to get the vector field

    $$\begin{aligned} {\mathbf {Y}}\equiv \left( \begin{array}{c} Y_{M}^{(1)}\\ Y_{M}^{(2)}\\ 0 \end{array}\right) = J_B \cdot \left( \begin{array}{c} \alpha _M^{(1)}\\ \alpha _M^{(2)}\\ 0 \end{array}\right) \, . \end{aligned}$$
    (5.15)
  3. 3.

    Thanks to the form of Benjamin’s Poisson structure, it turns out that \({\mathbf {Y}}\) is still supported on the locus \(\{z=\zeta _1\}\cup \{z=\zeta _2\}\), and can be easily projected with \(\pi _*\) to obtain the vector field \((X_1,X_2,X_3,X_4)\) on \(\mathcal {S}\). The latter depends linearly on \(\{\alpha _S^{(i)}\}_{i=1,\ldots ,4}\) and defines the reduced Poisson operator \(\mathcal {P}\) on \(\mathcal {S}\).

This construction essentially works as in the two-layer case considered in Camassa et al. (2017), provided one point is taken into account, namely that the integrals in the second and fourth components of \(\pi \) have a variable upper bound. We have, for \(({\dot{\rho }},{\dot{\Sigma }})\in TM\big \vert _\mathcal {I}\),

$$\begin{aligned} \pi _{*}\left( \begin{array}{c} {\dot{\rho }}\\ {\dot{\Sigma }}\\ \dot{{\overline{\zeta }}} \end{array} \right) = \left( \begin{array}{l} \int _0^h {\dot{\rho }}\, \mathrm {d}z\\ \int _0^{{\overline{\zeta }}} {\dot{\rho }} \, \mathrm {d}z+\dot{{\overline{\zeta }}}\,(\rho (x, {\overline{\zeta }})-\rho _3)\\ \int _0^h {\dot{\Sigma }}\, \mathrm {d}z\\ \int _0^{{\overline{\zeta }}} {\dot{\Sigma }} \, \mathrm {d}z+\dot{{\overline{\zeta }}}\, \Sigma (x, {\overline{\zeta }}) \,. \end{array} \right) \,. \end{aligned}$$
(5.16)

We remark that since the inequalities \(\zeta _2<{\overline{\zeta }}\equiv {\displaystyle {\frac{\zeta _1+\zeta _2}{2}}}<\zeta _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 \(\rho (x, {\overline{\zeta }})-\rho _3=\rho _2-\rho _3\ne 0\). As anticipated above, on \(\mathcal {I}\) the mean of the tangent vector component coming from the \(\zeta \)’s, \(\dot{{\bar{\zeta }}}= ({\dot{\zeta }}_1+{\dot{\zeta }}_2)/2\) can be expressed in terms of the tangent vector component \({\dot{\rho }}\). To this end, we can use the analogue of relations (5.9), which generically give

$$\begin{aligned} {\dot{\rho }}=(\rho _3-\rho _2){\dot{\zeta }}_2\delta (z-\zeta _2)+(\rho _2-\rho _1){\dot{\zeta }}_1\delta (z-\zeta _1)\, . \end{aligned}$$
(5.17)

Integrating this with respect to z in [0, h] yields

$$\begin{aligned} \int \limits _0^h {\dot{\rho }}\, \mathrm {d}z=(\rho _3-\rho _2){\dot{\zeta }}_2+(\rho _2-\rho _1){\dot{\zeta }}_1\, , \end{aligned}$$
(5.18)

while by integrating over \([0, {\overline{\zeta }}]\), we obtain

$$\begin{aligned} \int \limits _0^{{\overline{\zeta }}} {\dot{\rho }}\, \mathrm {d}z=(\rho _3-\rho _2){\dot{\zeta }}_2\, . \end{aligned}$$
(5.19)

Solving the linear system given by (5.18) and (5.19) gives

$$\begin{aligned} (\rho _2-\rho _3)\dot{{\overline{\zeta }}}=\Delta _\rho \int \limits _0^h{\dot{\rho }}\, \mathrm {d}z-(\frac{1}{2}+\Delta _\rho )\int \limits _0^{{\overline{\zeta }}} {\dot{\rho }}\, \mathrm {d}z\, , \end{aligned}$$
(5.20)

where

$$\begin{aligned} \Delta _\rho :={\displaystyle { \frac{\rho _2-\rho _3}{2(\rho _2-\rho _1)}}}\, .\end{aligned}$$
(5.21)

Thus, formula (5.16) for the components of the push forward \(\pi _*\) does not explicitly depend on the last component \(\dot{{\overline{\zeta }}}\) and can be written as

$$\begin{aligned} \pi _{*}\left( \begin{array}{c} {\dot{\rho }}\\ {\dot{\Sigma }}\\ \dot{{\overline{\zeta }}}\end{array} \right) = \left( \begin{array}{l} \int _0^h {\dot{\rho }}\, \mathrm {d}z\\[2mm] \Delta _\rho \int _0^h {\dot{\rho }}\, \mathrm {d}z +(\frac{1}{2}-\Delta _\rho ) \int _0^{{\overline{\zeta }}} {\dot{\rho }} \, \mathrm {d}z\\[2mm] \int _0^h {\dot{\Sigma }}\, \mathrm {d}z\\[2mm] \int _0^{{\overline{\zeta }}} {\dot{\Sigma }} \, \mathrm {d}z \end{array} \right) \,. \end{aligned}$$
(5.22)

By substituting this result in relation (5.14), we find

$$\begin{aligned} \begin{aligned}&\int \limits _{-\infty }^{+\infty } \int \limits _0^h ({\dot{\rho }}\, \alpha ^{(1)}_M+{\dot{\Sigma }}\, \alpha ^{(2)}_M)\, \mathrm {d}x\, \mathrm {d}z\\&\quad = \int \limits _{-\infty }^{+\infty } \left( \int \limits _0^h {\dot{\rho }} (\alpha ^{(1)}_S+\Delta _\rho \alpha ^{(2)}_S)\, \mathrm {d}z+\right. \left. \int \limits _0^{{\overline{\zeta }}}\, {\dot{\rho }}\, (\frac{1}{2}-\Delta _\rho ) \alpha ^{(2)}_S) \,\mathrm {d}z\right. \\&\qquad \left. +\int \limits _0^h {\dot{\Sigma }}\alpha ^{(3)}_S\, \mathrm {d}z+\int \limits _0^{{\overline{\zeta }}} {\dot{\Sigma }}\alpha ^{(4)}_S \, \mathrm {d}z\, \right) \mathrm {d}x\, . \end{aligned} \end{aligned}$$
(5.23)

Hence, the non-vanishing components of the 1-form \({\varvec{\alpha }}_M\) pulled back to \(\mathcal {I}\) are

$$\begin{aligned} \alpha ^{(1)}_M= \alpha ^{(1)}_S+\Delta _\rho \alpha ^{(2)}_S\theta ({{\overline{\zeta }}}-z) +\left( \frac{1}{2}-\Delta _\rho \right) \alpha ^{(2)}_S \,,\qquad \alpha ^{(2)}_M= \alpha ^{(3)}_S+\alpha ^{(4)}_S\theta ({{\overline{\zeta }}}-z)\,.\nonumber \\ \end{aligned}$$
(5.24)

Applying the Poisson tensor (5.5) to this 1-form yields, after some manipulations that crucially use the fact that products of Dirac’s \(\delta \) supported at different locations can be consistently set to zero, the vector field

$$\begin{aligned} \begin{aligned} Y_{M}^{(1)}=&\big ((\rho _{2}-\rho _{3})\delta (z-\zeta _{2})+(\rho _{1}-\rho _{2})\delta (z-\zeta _{1})\big )\, \alpha _{S,x}^{(3)} +(\rho _{2}-\rho _{3})\delta (z-\zeta _{2})\, \alpha _{S,x}^{(4)}\\ Y_{M}^{(2)}=&\left( (\rho _{2}-\rho _{3})\delta (z-\zeta _{2})+(\rho _{1}-\rho _{2})\delta (z-\zeta _{1})\right) \alpha _{S,x}^{(1)}\\&+\big (\frac{1}{2}(\rho _2-\rho _3)\delta (z-\zeta _2)+(\rho _1-\rho _2){\Delta _\rho }\delta (z-\zeta _1)\big )\alpha _{S,x}^{(2)} \\&+\big ((\tau _1-\tau _2)\delta '(z-\zeta _{1}){\zeta _1}_x + \tau _2\delta '(z-\zeta _{2}){\zeta _2}_x \big )\alpha _{S,x}^{(3)} \,. \end{aligned} \end{aligned}$$
(5.25)

The push-forward under the map \(\pi _*\) of this vector field gives the following four components:

$$\begin{aligned} \begin{aligned} \int \limits _0^h Y_{M}^{(1)}\, \mathrm {d}z&=(\rho _1-\rho _3)\alpha _{S,x}^{(3)}+(\rho _{2}-\rho _{3})\alpha _{S,x}^{(4)}\\ \Delta _\rho \int \limits _0^h Y_{M}^{(1)}\,&\mathrm {d}z +\left( \frac{1}{2}-\Delta _\rho \right) \int \limits _0^{{\overline{\zeta }}} Y_1 \, \mathrm {d}z\\&={\Delta _\rho }\left( (\rho _1-\rho _3)\alpha _{S,x}^{(3)}+(\rho _{2}-\rho _{3})\alpha _{S,x}^{(4)}\right) \\&\quad +\left( \frac{1}{2}-{\Delta _\rho }\right) (\rho _2-\rho _3)(\alpha _{S,x}^{(3)}+\alpha _{S,x}^{(4)})\\&=\big ({\Delta _\rho }(\rho _1-\rho _3)+\left( \frac{1}{2}-{\Delta _\rho }\right) (\rho _2-\rho _3)\big )\alpha _{S,x}^{(3)}+\frac{1}{2}(\rho _2-\rho _3)\alpha _{S,x}^{(4)}\\&=\frac{1}{2}(\rho _2-\rho _3)\alpha _{S,x}^{(4)} \\ \int \limits _0^h Y_{M}^{(2)}\, \mathrm {d}z&=(\rho _1-\rho _2)\alpha _{S,x}^{(1)}+\big (\frac{\rho _2-\rho _3}{2}+{\Delta _\rho }\, (\rho _1-\rho _2)\big )\alpha _{S,x}^{(2)}=(\rho _1-\rho _3)\alpha _{S,x}^{(1)}\\ \int \limits _0^{{\overline{\zeta }}}Y_{M}^{(2)}\, \mathrm {d}z&=(\rho _2-\rho _3)\alpha _{S,x}^{(1)}+\frac{\rho _2-\rho _3}{2} \alpha _{S,x}^{(2)}\, , \end{aligned} \end{aligned}$$
(5.26)

where we used definition (5.21) of \({\Delta _\rho }\) and the fact that the terms with the z-derivatives of the Dirac’s \(\delta \) 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 \({{\mathcal {P}}}\) on \({{\mathcal {S}}}\), in the coordinates \((\xi _1,\xi _2, \tau _1,\tau _2)\), as

$$\begin{aligned} {{\mathcal {P}}}=\left( \begin{array}{cccc} 0&{}0&{}\rho _1-\rho _3&{}\rho _2-\rho _3\\ 0&{}0&{}0&{}\frac{1}{2}(\rho _2-\rho _3)\\ \rho _1-\rho _3&{}0&{}0&{}0\\ \rho _2-\rho _3&{}\frac{1}{2}(\rho _2-\rho _3)&{}0&{}0\end{array} \right) \partial _x\, . \end{aligned}$$
(5.27)

The variables \((\xi _1,\xi _2, \tau _1,\tau _2)\) are related to \((\zeta _1,\zeta _2,\sigma _1=\rho _2{{\overline{u}}}_2-\rho _1{{\overline{u}}}_1,\sigma _2=\rho _3{{\overline{u}}}_3-\rho _2{{\overline{u}}}_2)\) by

$$\begin{aligned} \begin{aligned} \xi _1&=\left( h-\zeta _{{2}} \right) \left( \rho _{{2}}-\rho _{{3}} \right) + \left( h-\zeta _{{1}} \right) \left( \rho _{{1}}-\rho _{{2}} \right) \,,\qquad \xi _2=\frac{1}{2}\, \left( \rho _{{2}}-\rho _{{3}} \right) \left( \zeta _{{1}}-\zeta _{{2 }} \right) \,, \\ \tau _1&=\sigma _1+\sigma _2\,,\qquad \tau _2=\sigma _2\,. \end{aligned}\nonumber \\ \end{aligned}$$
(5.28)

Solving these relations with respect to the \(\zeta _i\)’s and the \(\sigma _i\)’s gives:

$$\begin{aligned} \begin{aligned} \zeta _1&=-{\frac{\xi _{{1}}}{\rho _{{1}}-\rho _{{3}}}}+{\frac{2\xi _{{2}}}{\rho _{{1}}-\rho _{{3}}}}+h\,,\qquad \zeta _2=-{\frac{\xi _{{1}}}{\rho _{{1}}-\rho _{{3}}}}-{\frac{2 \left( \rho _{{ 1}}-\rho _{{2}} \right) \xi _{{2}}}{ \left( \rho _{{2}}-\rho _{{3}} \right) \left( \rho _{{1}}-\rho _{{3}} \right) }}+h\,,\\ \sigma _1&=\tau _1-\tau _2\,,\qquad \sigma _2=\tau _2 . \end{aligned}\nonumber \\ \end{aligned}$$
(5.29)

A straightforward computation shows that in these coordinates the Poisson operator (5.27) acquires the particularly simple form

$$\begin{aligned} {{\mathcal {P}}}=\left( \begin{array}{cccc} 0&{}0&{}-\partial _x&{}0\\ 0&{}0&{}0&{}-\partial _x\\ -\partial _x&{}0&{}0&{}0\\ 0&{}-\partial _x&{}0&{}0 \end{array} \right) \,. \end{aligned}$$
(5.30)

Remark 5.1

According to the terminology favored by the Russian school, for Hamiltonian quasi-linear systems of PDEs (see, e.g., Tsarev 1991), the coordinates \((\xi _1, \xi _2, \tau _1, \tau _2)\) and, a fortiori, the coordinates \((\zeta _1,\zeta _2, \sigma _1,\sigma _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.2

It should be clear how to proceed in the n-layered case, with a stratification given by densities \(\rho _1<\rho _2<\cdots <\rho _n\) and interfaces \(\zeta _1>\zeta _2>\cdots >\zeta _{n-1}\). We consider intervals

$$\begin{aligned} I_1=[0,h],\, I_2=\left[ 0, \frac{\zeta _1+\zeta _2}{2}\right] , I_3=\left[ 0, \frac{\zeta _2+\zeta _3}{2}\right] ,\,\dots \,, I_n=\left[ 0, \frac{\zeta _{n-2}+\zeta _{n-1}}{2}\right] \,,\nonumber \\ \end{aligned}$$
(5.31)

and define the analogue of the map (5.13) by integrating both \(\rho (z)-\rho _n\) and \(\sigma (x,z)\) over each of the intervals \(I_k\), with \(k=1,\ldots ,n\). The procedure can be repeated verbatim. We conjecture that in the long-wave limit, the quantities

$$\begin{aligned} (\zeta _1,\zeta _2,\ldots , \zeta _{n-1}, \sigma _1,\sigma _2,\ldots , \sigma _{n-1}), \end{aligned}$$
(5.32)

where \(\sigma _k=\rho _{k+1}{{\overline{u}}}_{k+1}-\rho _{k}{{\overline{u}}}_{k}\), for \(k=1,\ldots , n-1\), are flat Darboux coordinates for the reduced Poisson structure. We verified that this conjecture holds true in the four-layer case.

Remark 5.3

The steps performed here are borrowed from the Poisson reduction theory of Marsden, Ratiu and Weinstein (see, e.g., Marsden and Ratiu 1986); 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.

5.3 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,

$$\begin{aligned} H=\int \limits _{-\infty }^{+\infty } \int \limits _0^h \frac{1}{2}\,\rho \left( u^2+w^2\right) \, \mathrm {d}x\, \mathrm {d}z+\int \limits _{-\infty }^{+\infty } \int \limits _0^h g\rho z\, \mathrm {d}x\, \mathrm {d}z\,. \end{aligned}$$
(5.33)

The potential energy is readily reduced, using the first of (5.6), to

$$\begin{aligned} U= \int \limits _{-\infty }^{+\infty } \frac{1}{2}\left( \,g \left( \rho _{{2}}-\rho _{{1}} \right) {\zeta _{{1}}}^{2}+ g \left( \rho _{{3}}-\rho _{{2}} \right) {\zeta _{{2}}}^{2}\right) \, \mathrm {d}x\,. \end{aligned}$$
(5.34)

When the layer thicknesses are not asymptotically zero, both energies can be appropriately renormalized subtracting the far-field contributions of \(\zeta _i\). To obtain the reduced kinetic energy density, we use the fact that at order \(O(\epsilon ^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

$$\begin{aligned}&\mathcal {T}=\frac{1}{2}\left( \int \limits _0^{\zeta _2} \rho _3 {{\overline{u}}}_3^2 \, \mathrm {d}z+\int \limits _{\zeta _2}^{\zeta _1}\rho _2 {{\overline{u}}}_2^2\, \mathrm {d}z+\int \limits _{\zeta _1}^h \rho _1 {{\overline{u}}}_1^2\, \mathrm {d}z\right) \nonumber \\&\quad =\frac{1}{2}\left( \rho _3{\zeta _2}{{\overline{u}}}_3^2 + \rho _2({\zeta _1-\zeta _2}){{\overline{u}}}_2^2+ \rho _1{(h-\zeta _1)}{{\overline{u}}}_1^2\right) \, , \end{aligned}$$
(5.35)

so that the reduced kinetic energy is

$$\begin{aligned} T=\int \limits _{-\infty }^{+\infty } \frac{1}{2}\left( \rho _3{\zeta _2}{{\overline{u}}}_3^2 + \rho _2({\zeta _1-\zeta _2}){{\overline{u}}}_2^2+ \rho _1{(h-\zeta _1)}{{\overline{u}}}_1^2\right) \mathrm {d}x\, . \end{aligned}$$
(5.36)

We now use the dynamical constraint for localized solutions, whereby velocities vanish at infinity,

$$\begin{aligned} (h-\zeta _1){{\overline{u}}}_1+(\zeta _1-\zeta _2){{\overline{u}}}_2+\zeta _2{{\overline{u}}}_3=0\, , \end{aligned}$$
(5.37)

to obtain

$$\begin{aligned} {{\overline{u}}}_{{1}}={\frac{(\zeta _{{1}}-\zeta _{{2}}){{\overline{u}}}_{{2}}+\zeta _{{2 }}{{\overline{u}}}_{{3}}}{\zeta _{{1}}-h}}\, . \end{aligned}$$
(5.38)

Next, we express \({{\overline{u}}}_2,{{\overline{u}}}_3\) in terms of \(\sigma _1,\sigma _2\) as

$$\begin{aligned} \begin{aligned} {{\overline{u}}}_2&=\frac{\rho _{{3}} \left( h-\zeta _{{1}} \right) \sigma _{{1}}}{\Psi }-\frac{\zeta _{{2}}\rho _{{1 }}\sigma _{{2}}}{\Psi }\\ {{\overline{u}}}_3&=\frac{\rho _2(h-\zeta _1)\sigma _1}{\Psi }+\frac{\left( h\rho _{{2}}+( \rho _{{1}}-\rho _{{2}}) \zeta _{{1}}-\zeta _{{2} }\rho _{{1}}\right) \sigma _2}{\Psi }\, , \end{aligned} \end{aligned}$$
(5.39)

where

$$\begin{aligned} \Psi =h\rho _{{2}}\rho _{{3}}-\rho _{{3}} \left( \rho _{{2}}-\rho _{{1}} \right) \zeta _{{1}}-\rho _{{1}} \left( \rho _{{3}}-\rho _{{2}} \right) \zeta _{{2}}\,. \end{aligned}$$
(5.40)

The kinetic energy density turns out to be, in the new set of variables,

$$\begin{aligned} \begin{aligned} \mathcal {T}=&{1\over 2\Psi } \Big ( (h-\zeta _1)\left( \rho _3\zeta _1+ (\rho _2-\rho _3) \zeta _2\right) \sigma _1^2 \\&+2(h-\zeta _1)\rho _2\,\zeta _2\,\sigma _1\sigma _2+ \left( (\rho _1-\rho _2) \zeta _2\,\zeta _1+\rho _2\, \zeta _2 \,h-\rho _1\,\zeta _2^2 \right) \sigma _2^2 \Big )\,, \end{aligned} \end{aligned}$$
(5.41)

so that the Hamiltonian functional is

$$\begin{aligned} H=\int \limits _{-\infty }^{+\infty } \left( \mathcal {T}+\frac{1}{2}\,g \left( ( \rho _{{2}}-\rho _{{1}}) \zeta _1^{2}+ ( \rho _{{3}}-\rho _{{2}}) \zeta _2^{2}\right) \right) \mathrm {d}x\,. \end{aligned}$$
(5.42)

Explicitly, the equations of motion can be written as conservation laws,

$$\begin{aligned} \left( \begin{array}{c} {\zeta _1}_t\\ {\zeta _2}_t\\ {\sigma _1}_t\\ {\sigma _2}_t \end{array}\right) =- \left( \begin{array}{c} \left( \delta _{\sigma _{1}} H\right) _x\\ \left( \delta _{\sigma _{2}} H\right) _x\\ \left( \delta _{\zeta _{1}} H\right) _x\\ \left( \delta _{\zeta _{2} }H\right) _x \end{array} \right) \, , \end{aligned}$$
(5.43)

where the gradient of the Hamiltonian is, explicitly,

$$\begin{aligned} \begin{aligned} \delta _{\zeta _{1}} H&= \frac{1}{2\Psi }\left( \left( h\rho _{{3}}-2\,\zeta _{{1}}\rho _{{3}}+\zeta _{{2}}(\rho _{{3}}-\rho _{{2}})\right) {\sigma _{{1}}}^{2}-2\,\zeta _{{2}}\rho _{{2}} \sigma _{{2}}\sigma _{{1}}-\zeta _{{2}} \left( \rho _{{2}}-\rho _{{1}} \right) {\sigma _{{2}}}^{2}\right) \\&\quad -\frac{\rho _3(\rho _1-\rho _2)\, \mathcal {T}}{\Psi }+g \left( \rho _{{2}}-\rho _{{1}} \right) \zeta _{{1}}\,,\\ \delta _{\zeta _{2} }H&= \frac{1}{2\Psi } \left( \left( h-\zeta _{{1}} \right) \left( \rho _{{2}}-\rho _{{3}} \right) {\sigma _{{1}}}^{2}+2\, \left( h-\zeta _{{1}} \right) \rho _{{2}}\sigma _{{2 }}\sigma _{{1}}\right. \\&\quad \left. + \left( \rho _{{2}}h+\zeta _{{1}}\rho _{{1}}-\zeta _{{1}}\rho _{{2}}-2\,\zeta _{{2}}\rho _{{1}} \right) {\sigma _{{2}}}^{2} \right) \\&\quad -\frac{\rho _{{1}} \left( \rho _{{2}}-\rho _{{3}} \mathcal {T}\right) }{\Psi }+g \left( \rho _{{3}}-\rho _{{2}} \right) \zeta _{{2}}\,,\\ \delta _{\sigma _{1}} H&=\frac{1}{\Psi } \left( \left( \zeta _{{1}}\rho _{{3}}+\zeta _{{2}}\rho _{{2}}-\zeta _{{2}}\rho _{{3}} \right) \left( h-\zeta _{{1}} \right) \sigma _{{1}}+\rho _{{2}}\zeta _{{2} } \left( h-\zeta _{{1}} \right) \sigma _{{2}} \right) \,,\\ \delta _{\sigma _{2}} H&= \frac{1}{\Psi } \left( \rho _{{2}}\zeta _{{2}} \left( h-\zeta _{{1}} \right) \sigma _{{1}}+\zeta _{{2 }} \left( h\rho _{{2}}+\zeta _{{1}}(\rho _{{1}}-\rho _{{2}})-\zeta _{ {2}}\rho _{{1}} \right) \sigma _{{2}}\right) \,. \end{aligned} \end{aligned}$$
(5.44)

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,

$$\begin{aligned} \begin{aligned}&Z_j=\int \limits _{-\infty }^{+\infty } \zeta _j\, \mathrm {d}x,\qquad S_j=\int \limits _{-\infty }^{+\infty } \sigma _j\, \mathrm {d}x, \quad j=1,2\\&H\> \text { given by equation (5.42)},\quad K=\int \limits _{-\infty }^{+\infty } (\zeta _1\sigma _1+\zeta _2\sigma _2) \, \mathrm {d}x\, . \end{aligned} \end{aligned}$$
(5.45)

In particular, the quantity K indicates how the momentum paradox mentioned in Sect. 4.2 can be viewed from the canonical formulation of the evolution equations. In fact, we see that, rather than the horizontal momentum \(\Pi ^{(x)}\), it is K that plays the role of generator of x-translations ( cf. Benjamin 1986 for the case of continuous stratification). By expressing its density \(\mathcal {K}=\zeta _1\sigma _1+\zeta _2\sigma _2 \) in terms of the horizontal mean velocities, we obtain

$$\begin{aligned} \mathcal {K}=\rho _{{2}}{{\overline{u}}}_{{2}} \left( \zeta _{{1}}-\zeta _ {{2}} \right) +\rho _{{3}}{{\overline{u}}}_{{3}}\zeta _{{2}}-\rho _{{1}}{{\overline{u}}}_{{1}}\zeta _{{1}}=\pi ^{(x)}-h\rho _1{{\overline{u}}}_1\,, \end{aligned}$$
(5.46)

where

$$\begin{aligned} \pi ^{(x)}=\left( h-\zeta _{{1}} \right) \rho _{{1}}{{\overline{u}}}_{{1}}+\rho _{{2}}{{\overline{u}}}_{{2}} \left( \zeta _{{1}}-\zeta _{{2}} \right) +\rho _{{3}}\zeta _{{2}}{{\overline{u}}}_{{3}} \, \end{aligned}$$

is the total horizontal momentum density. Thus, the total momentum \(\Pi ^{(x)}=\int _{-\infty }^{+\infty } \pi ^{(x)}\, \mathrm {d}x\) is not conserved, while \(K=\int _{-\infty }^{+\infty } \mathcal {K} \mathrm {d}x\) is.

6 The Boussinesq Approximation

A dramatic simplification of the problem is provided by the so-called Boussinesq approximation, that is, the double scaling limit

$$\begin{aligned} \rho _i\rightarrow {{\overline{\rho }}}=\frac{\rho _1+\rho _2+\rho _3}{3},\quad i=1,2,3, \quad \text {with}\ g(\rho _1-\rho _2)\ \text {and}\ g(\rho _2-\rho _3)\ \text {both finite.}\nonumber \\ \end{aligned}$$
(6.1)

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

$$\begin{aligned} {\frac{ \left( h-{\zeta _{{1}}}\right) \zeta _{{1}} {\sigma _{{1}}}^{2}}{2h{{{\overline{\rho }}}}}}+{\frac{\zeta _{{2}} \left( h-\zeta _{{1}} \right) \sigma _{{2}}\sigma _{{1} }}{h{{\overline{\rho }}}}}+{\frac{\zeta _{{2}} \left( h-\zeta _{{2}} \right) {\sigma _{{2}}}^{2}}{2h{{{\overline{\rho }}}}}}\,, \end{aligned}$$
(6.2)

so that the Hamiltonian energy functional in this Boussinesq limit is

$$\begin{aligned} \begin{aligned} H_B=&\int \limits _{-\infty }^{+\infty } \frac{1}{2\, h{{\overline{\rho }}}} \left( { { \zeta _1(h-\zeta _{{1}}){\sigma _{{1}}}^{2}}}+ 2{ {\zeta _{{2}} \left( h-\zeta _{{1}} \right) \sigma _{{1}}\sigma _{{2}}}} +{ {\zeta _{{2}} \left( h-\zeta _{{2}} \right) {{\sigma _2}^2}}} \right) \, \\&+g(\left( \rho _{{2}}-\rho _{{1}} \right) {\zeta _{{1}}}^{2}+ \left( \rho _{{3}}-\rho _{{2}} \right) {\zeta _{{2}}}^{2}) \,\mathrm {d}x\,. \end{aligned} \end{aligned}$$
(6.3)

The ensuing equations of motion are

$$\begin{aligned} \left( \begin{array}{c} {\zeta _1}_t\\ {\zeta _2}_t\\ {\sigma _1}_t\\ {\sigma _2}_t \end{array}\right) =\mathcal {P} \left( \begin{array}{c} \delta _{\zeta _{1}} H_B\\ \delta _{\zeta _{2}} H_B\\ \delta _{\sigma _{1}}H_B\\ \delta _{\sigma _{2} }H_B \end{array} \right) \,, \end{aligned}$$
(6.4)

with \(\mathcal {P}\) the “canonical" Poisson tensor (5.30). As a system of quasi-linear equations, they can be cast in the form

$$\begin{aligned} \left( \begin{array}{c} {\zeta _1}_t\\ {\zeta _2}_t\\ {\sigma _1}_t\\ {\sigma _2}_t \end{array}\right) + {\mathbf {A}} \left( \begin{array}{c} {\zeta _1}_x\\ {\zeta _2}_x\\ {\sigma _1}_x\\ {\sigma _2}_x \end{array}\right) =0\, , \end{aligned}$$
(6.5)

where the characteristic matrix reads

$$\begin{aligned} {\mathbf {A}}= \frac{1}{h\, {{\overline{\rho }}}}\left( \begin{array}{cccc} \begin{array}{l} { \left( 2\,\zeta _1 -h\right) \sigma _1}\\ +{\sigma _2\zeta _2}\end{array} &{} { \left( \zeta _{{1}} -h \right) \sigma _{{2}} }&{}\zeta _{{1}} \left( \zeta _{{1}} -h \right) &{} \zeta _{{2}} \left( \zeta _{{1}} -h \right) \\ \zeta _{{2}}\sigma _1&{} \begin{array}{l}\left( \zeta _{{1}} -h\right) \sigma _{{1}}\\ +\left( 2\,\zeta _{{2}}-h \right) \sigma _{{2}}\end{array}&{} \zeta _{{2}} \left( \zeta _{{1}}-h \right) &{}\zeta _{{2}} \left( \zeta _{{2}} -h\right) \\ \sigma _{{1}}^{2}+{\tilde{g}} \left( \rho _{{1}}-\rho _{{2}} \right) &{} \sigma _{{1}}\sigma _{{2}}&{}\begin{array}{l}\left( 2\,\zeta _{{1}}-h \right) \sigma _{{1}}\\ +\sigma _{{2}}\zeta _{{2}}\end{array}&{}\zeta _{{2}}\sigma _{{1}}\\ \sigma _{{1}}\sigma _{{2}}&{}\sigma _{{2}}^{2}+{\tilde{g}} \left( \rho _{{2}}-\rho _{{3}} \right) &{}\left( \zeta _{{1}}-h \right) \sigma _{{2}} &{}\begin{array}{l} \left( \zeta _{{1}}-h \right) \sigma _{{1}}\\ +\left( 2\,\zeta _{{2}} -h\right) \sigma _{{2}}\end{array} \end{array}\right) \nonumber \\ \end{aligned}$$
(6.6)

and \({\tilde{g}} =g\, h\, {{\overline{\rho }}}\) 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 \(\mathcal {H}\) of the matrix \({\mathbf {A}}\). A lengthy computation shows that \(\mathcal {H}\) has 12 (out of 24) non-vanishing components, namely

$$\begin{aligned} \begin{aligned} {{{\mathcal {H}}}^1}_{1\, 2}&={{{\mathcal {H}}}^3}_{2\, 3}=-{{{\mathcal {H}}}^4}_{1\, 3}=-{\frac{g\,r_{{32}}\,\zeta _{{2}} \left( h-\zeta _{{1}} \right) {\sigma _{{1}} }^{2}}{{{{\bar{\rho }}}}^{3}{h}^{2}}}\\&\quad +{\frac{g\,r_{{21}}\, \left( \zeta _{{1}}-2\,\zeta _{{2}} \right) \left( h-\zeta _{{1}} \right) {\sigma _{{2}}}^{2}}{{{{\bar{\rho }}}} ^{3}{h}^{2}}}+{\frac{{g}^{2}r_{{21}}r_{{32}} \zeta _{{2}} \left( h-\zeta _ {{1}} \right) }{{{{\bar{\rho }}}}^{2}h}}\\ {{{\mathcal {H}}}^2}_{1\, 2}&={{{\mathcal {H}}}^3}_{2\, 4}=-{{{\mathcal {H}}}^4}_{1\, 4}=-{\frac{g\, r_{{32}}\,\zeta _{{2}} \left( \zeta _{{2}}+h-2\,\zeta _{{1}} \right) {\sigma _{{1}}}^{2}}{{{{\bar{\rho }}}}^{3}{h}^{2}}}\\&\quad -{\frac{g\, r_{{21}}\, \zeta _{{2}} \left( h-\zeta _{{1}} \right) {\sigma _{{2}}}^{2}}{{{{\bar{\rho }}}}^{3}{h}^{ 2}}}+{\frac{{g}^{2}\, r_{{21}}r_{{32}}\, \zeta _{{2}} \left( h-\zeta _{{1}} \right) }{{{{\bar{\rho }}}}^{2}h}}\\ {{{\mathcal {H}}}^1}_{1\, 4}&=-{{{\mathcal {H}}}^2}_{1\, 3}=-{{{\mathcal {H}}}^3}_{3\, 4}=2\,{\frac{g\, r_{{21}}\zeta _{{2}}\sigma _{{2}} \left( \zeta _{{1}}-\zeta _{{2} } \right) \left( h-\zeta _{{1}} \right) }{{{{\bar{\rho }}}}^{3}{h}^{2}}} \\ {{{\mathcal {H}}}^1}_{2\, 4}&=-{{{\mathcal {H}}}^2}_{2\, 3}=-{{{\mathcal {H}}}^4}_{3\, 4}= -2\,{\frac{g\, r_{{32}}\,\zeta _{{2}}\sigma _{{1}} \left( \zeta _{{1}}-\zeta _{{2 }} \right) \left( h-\zeta _{{1}} \right) }{{{{\bar{\rho }}}}^{3}{h}^{2}}}\, \end{aligned} \end{aligned}$$
(6.7)

where \(r_{ij}=\rho _i-\rho _j\). As should be expected from the free surface case, even in the Boussinesq approximation the model for three-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.

6.1 Symmetric Solutions

In the recent paper (de Melo Viríssimo and Milewski 2019), the authors have focused on the symmetric solutions defined by the equality of the upper (\(i=1\)) and lower (\(i=3\)) layer thicknesses, i.e., \(\zeta _2=h-\zeta _1\), and the averaged horizontal velocities, \({{\overline{u}}}_1={{\overline{u}}}_3\). In the Boussinesq approximation, our variables \((\sigma _1,\sigma _2)\) are actually proportional to the velocity shears,

$$\begin{aligned} \sigma _1={{\overline{\rho }}}({{\overline{u}}}_2-{{\overline{u}}}_1),\quad \sigma _2={{\overline{\rho }}}({{\overline{u}}}_3-{{\overline{u}}}_2)\, , \end{aligned}$$
(6.8)

so that the symmetric solutions found in de Melo Viríssimo and Milewski (2019) are given by the relations

$$\begin{aligned} \zeta _2=h-\zeta _1,\quad \sigma _2=-\sigma _1\,. \end{aligned}$$
(6.9)

A straightforward computation confirms that the submanifold defined by these relations is invariant under the flow (6.5) if and only if the relation

$$\begin{aligned} \rho _3-\rho _2=\rho _2-\rho _1 \end{aligned}$$
(6.10)

is fulfilled among the density differences. In this case system (6.5) reduces to a system with 2 “degrees of freedom," parameterized, e.g., by the pair \((\zeta _2\equiv \zeta , \sigma _2\equiv \sigma )\). The reduced “symmetric" equations of the motion are

$$\begin{aligned} \left( \begin{array}{c} \zeta _{t}\\ \sigma _t\end{array} \right) =\frac{1}{{{\overline{\rho }}}h} \left( \begin{array}{cc} \left( 4\,\zeta - \,h\right) \sigma &{} { { \zeta \left( 2\,\zeta -h\right) }}\\ 2\,{\sigma }^{2}-g{{\overline{\rho }}}\rho _{{}_\Delta }h &{} \left( 4\,\zeta - \,h\right) \sigma \end{array}\right) \left( \begin{array}{c} \zeta _{x}\\ \sigma _x\end{array} \right) \, , \end{aligned}$$
(6.11)

where we have defined \(\rho _{{}_\Delta }=\rho _3-\rho _2\). These equations follow from the Hamiltonian functional

$$\begin{aligned} H_{B,S}=\int \limits _{-\infty }^{+\infty } \left( {\frac{ \zeta \left( h- 2\,\zeta \right) {\sigma }^{2}}{2{{\overline{\rho }}}h}}+\frac{1}{2}\,g\rho _{{}_\Delta }\left( \zeta -{h\over 2} \right) ^2 \right) \mathrm {d}x\, , \end{aligned}$$
(6.12)

with the “standard" Poisson tensor

$$\begin{aligned} \mathcal {P}_{(2)}=\left( \begin{array}{cc} 0&{}-\partial _x\\ -\partial _x&{}0\end{array}\right) \, , \end{aligned}$$
(6.13)

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,

$$\begin{aligned} \lambda _\pm =\frac{1}{{{\overline{\rho }}}h} \left( \sigma (h-4\zeta )\pm \sqrt{ \zeta \left( h-2\,\zeta \right) \left( g h{{\overline{\rho }}}\rho _{{}_\Delta }-2\,{\sigma }^{2} \right) }\right) \, ; \end{aligned}$$
(6.14)

since in this case \(\zeta \in (0,{h}/2)\), and the domain of hyperbolicity coincides with the rectangular region

$$\begin{aligned} \left( 0,\frac{h}{2}\right) \times \left( {\displaystyle {-\frac{1}{\sqrt{2}}\sqrt{gh{{\overline{\rho }}}\rho _{{}_\Delta }}}}, {\displaystyle {\frac{1}{\sqrt{2}}\sqrt{gh{{\overline{\rho }}}\rho _{{}_\Delta }}}}\right) \, . \end{aligned}$$

Alternatively, the Hamiltonian formulation shows that the simple map \(\zeta \rightarrow \eta \), \(\sigma \rightarrow h \sigma \), \({{\overline{\rho }}}\rightarrow \rho _2\), \(\rho _{{}_\Delta }\rightarrow \rho _2-\rho _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 Camassa et al. (2019) (as also remarked in de Melo Viríssimo and Milewski 2019, though only reported explicitly from the motion equation viewpoint in the PhD thesis de Melo Viríssimo 2018), 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 \(\mathcal {J}\) defined by

$$\begin{aligned} {\widetilde{\zeta }}_1=h-\zeta _2\,,\quad {\widetilde{\zeta }}_2=h-\zeta _1\,,\quad {\widetilde{\sigma }}_1=-\sigma _2\,,\quad {\widetilde{\sigma }}_2=-\sigma _1\,. \end{aligned}$$
(6.15)

Since its Jacobian is

$$\begin{aligned} \mathcal {J}_*=\left( \begin{array}{cccc} 0&{}-1&{}0&{}0\\ -1&{}0&{}0&{}0\\ 0&{}0&{}0&{}-1\\ 0&{}0&{}-1&{}0 \end{array} \right) \, , \end{aligned}$$
(6.16)

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 \(\rho _3-\rho _2=\rho _2-\rho _1\equiv \rho _\Delta \), changes by a constant term added to the linear term \(H_\Delta =-g\rho _\Delta \, ({\widetilde{\zeta }}_1+{\widetilde{\zeta }}_2)\). Since \( \int \limits _{-\infty }^{+\infty } H_\Delta \, \mathrm {d}x \) 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.

7 Summary and Conclusions

In this paper, we examined some structural properties of multilayered incompressible Euler flows in the long-wave regime. In order to present explicit computations, we mainly focused 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., Chumakova et al. 2009b; Camassa et al. 2012).

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 quasi-linear resulting system, a fact proven in Chesnokov et al. (2017).

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 non-conservation of the horizontal momentum, highlighted in Camassa et al. (2012) for the two-layer case, has been shown to persist for an n-layered situation, \(n\ge 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) three-layered fluid motions was derived by means of a geometrical reduction process from the full two-dimensional Hamiltonian structure introduced in Benjamin (1986). 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 de Melo Viríssimo and Milewski (2019) 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.

Future developments will include the extension of our formalism to the next order in the long-wave parameter asymptotics, with the inclusion in this framework of the dispersive terms \(D_i\) of the n-layer equations (2.6). Framing this problem within Dubrovin’s approach of the expansion of Hamiltonian PDEs in the dispersion parameter (see, e.g., Dubrovin 2008; Dubrovin et al. 2015) is a task worth pursuing, as well as the comparison of our setting with that of Percival et al. (2008).