3.1. Hamiltonian form of the equations: equations for the moving surfaces

This subsection is devoted to a reformulation of the governing equations and of the boundary conditions (2.3)–(2.8) by collapsing the dynamical variables into a one-dimensional representation. More precisely, we will rewrite (2.3)–(2.8) as a Hamiltonian system of partial differential equations of the form (cf. Olver Reference Olver1993)

(3.1) \begin{align} \partial _t\mathfrak{m}=J\frac {\delta H}{\delta \mathfrak{m}}, \end{align}

where $t\mapsto \mathfrak{m}(t)$ denotes a path in a Hilbert space $\mathcal{H}$ equipped with an inner product $({\cdot },{\cdot })$ , and the associated Hamiltonian functional $H:\mathcal{D}\subset \mathcal{H}\rightarrow \mathbb{R}$ is defined on a dense subspace $\mathcal{D}$ of $\mathcal{H}$ . In our considerations $\mathfrak{m}=(\eta _1,\eta _2,\eta _3,\xi _1,\xi _2,\xi _3)^{\rm T}$ , the Hilbert space is an appropriate Cartesian product of spaces of square-integrable functions and $\mathcal{D}=\mathcal{S}(\mathbb{R})\times \mathcal{S}(\mathbb{R})\times \mathcal{S}(\mathbb{R})\times \mathcal{S}(\mathbb{R})\times \mathcal{S}(\mathbb{R})\times \mathcal{S}(\mathbb{R})$ , where $\mathcal{S}(\mathbb{R})$ is the subspace of Schwartz functions in the variable $x$ . Moreover, $J=\begin{pmatrix}0 & I_3\\-I_3 &0 \end{pmatrix}$ (where $I_3$ is the $3\times 3$ identity matrix) defines a Poisson bracket for functionals by the formula

(3.2) \begin{align} \{F_1,F_2\}=\left (\frac {\delta F_1}{\delta \mathfrak{m}},J\frac {\delta F_2}{\delta \mathfrak{m}}\right )\!. \end{align}

We also denote by ${\delta H}/{\delta \mathfrak{m}}$ the variational derivative of $H$ at $\mathfrak{m}\in \mathcal{D}$ through the formula

(3.3) \begin{align} \lim _{\varepsilon \to 0}\frac {H(\mathfrak{m}+\varepsilon \mathfrak{m}_0)-H(\mathfrak{m})}{\varepsilon }=\left (\mathfrak{m}_0,\frac {\delta H}{\delta \mathfrak{m}}\right )\quad \textrm {for} \quad \mathfrak{m}_0\in \mathcal{D}, \end{align}

where $({\cdot },{\cdot })$ is the inner product in the Hilbert space $L^2(\mathbb{R})\times L^2(\mathbb{R})\times L^2(\mathbb{R})\times L^2(\mathbb{R})\times L^2(\mathbb{R})\times L^2(\mathbb{R})$ .

3.2. The Hamiltonian functional

The natural candidate for the Hamiltonian functional is the total energy of the flow given as $E_K+E_P$ , where $E_K$ stands for the kinetic energy, and $E_P$ denotes the potential energy. In terms of the quantities of the flow, the kinetic energy is written as

(3.4) \begin{align} E_K=\iint_{\varOmega _1\cup \varOmega _2\cup \varOmega _3}\hat {\rho }\!\left (\frac {\hat {u}^2+\hat {v}^2}{2}\right ) {\rm d}y\,{\rm d}x, \end{align}

which, due to stratification and (2.10), can be detailed as

(3.5) \begin{align} E_K & = \int _{\mathbb{R}}\int _{-h}^{-h_3+\eta _3(x,t)}\rho _3\! \left (\frac {u_3^2+v_3^2}{2}\right ){\rm d}y\, {\rm d}x +\int _{\mathbb{R}}\int _{-h_3+\eta _3(x,t)}^{\eta _2(x,t)}\rho _2\! \left (\frac {u_2^2+v_2^2}{2}\right ){\rm d}y\, {\rm d}x \nonumber\\& \quad +\int _{\mathbb{R}}\int _{\eta _2(x,t)}^{h_1+\eta _1(x,t)}\rho _1\! \left (\frac {u_1^2+v_1^2}{2}\right ){\rm d}y\, {\rm d}x \nonumber\\\quad & = \frac {\rho _3}{2}\int _{\mathbb{R}}\int _{-h}^{-h_3+\eta _3(x,t)}|\boldsymbol{\nabla }\varphi _3|^2\,{\rm d}y\, {\rm d}x+ \frac {\rho _2}{2}\int _{\mathbb{R}}\int _{-h_3+\eta _3(x,t)}^{\eta _2(x,t)}|\boldsymbol{\nabla }\varphi _2|^2\,{\rm d}y\, {\rm d}x \nonumber\\& \quad +\frac {\rho _1}{2}\int _{\mathbb{R}}\int _{\eta _2(x,t)}^{h_1+\eta _1(x,t)}|\boldsymbol{\nabla }\varphi _1|^2\,{\rm d}y\, {\rm d}x. \end{align}

In defining the potential energy, we act on the realisation that not all the latter energy can be converted into kinetic energy. However, the difference in the potential energy between the perturbed wave state and the pure background state is available for this conversion. Hence we have

(3.6) \begin{align} E_P & = g\rho _3\int _{\mathbb{R}}\left (\int _{-h}^{-h_3+\eta _3}y\,{\rm d}y-\int _{-h}^{-h_3}y\,{\rm d}y\right )+g\rho _2\int _{\mathbb{R}}\left (\int _{-h_3+\eta _3}^{\eta _2}y\,{\rm d}y-\int _{-h_3}^0y\,{\rm d}y\right )\nonumber\\& \quad + g\rho _1\int _{\mathbb{R}}\left (\int _{\eta _2}^{h_1+\eta _1}y\,{\rm d}y-\int _0^{h_1}y\,{\rm d}y\right )\nonumber\\& = \frac {g\rho _1}{2}\int _{\mathbb{R}}\eta _1^2\,{\rm d}x+\frac {g(\rho _2-\rho _1)}{2}\int _{\mathbb{R}}\eta _2^2\,{\rm d}x + \frac {g(\rho _3-\rho _2)}{2}\int _{\mathbb{R}}\eta _3^2\,{\rm d}x. \end{align}

We note that the above choice for the potential energy allows us to avoid any non-physical constant terms in the energy density and hence infinite energy. More specifically, the potential energy functional (3.6) is not convergent if the potential energy of the pure current background state is not subtracted.

3.3. Computations of variations of $H$

We first list a few useful formulas for computing the variations of $H$ .

We proceed now to compute the variations of the terms that make up the kinetic energy $E_K$ from (3.5). Utilising the divergence theorem (Evans Reference Evans2010), and (3.12) as well as (3.7), (3.8) and (3.9), we obtain

(3.13) \begin{align} \delta \left (\int _{\mathbb{R}}\int _{-h}^{-h_3+\eta _3}|\boldsymbol{\nabla }\varphi _3|^2\,{\rm d}y\,{\rm d}x\right )={}&2 \int _{\mathbb{R}}(\varphi _{3,y}-\eta _{3,x}\varphi _{3,x})_{s_3} \big (\delta \varPhi ^3-(\delta \eta _3)(\varphi _{3,y})_{s_3}\big ){\rm d}x \nonumber\\ & \quad +\int _{\mathbb{R}}|\boldsymbol{\nabla }\varphi _3|^2_{s_3}\delta \eta _3\, {\rm d}x. \end{align}

Analogously, by the same token we have

(3.14) \begin{align} \delta \left (\int _{\mathbb{R}}\int _{-h_3+\eta _3}^{\eta _2}|\boldsymbol{\nabla }\varphi _2|^2\,{\rm d}y\,{\rm d}x\right ) & = 2 \int _{\mathbb{R}} (\varphi _{2,y}-\eta _{2,x}\varphi _{2,x})_{s_2} \big (\delta \varPhi ^2_2-(\delta \eta _2)(\varphi _{2,y})_{s_2}\big ) {\rm d}x \nonumber\\& \quad -2 \int _{\mathbb{R}} (\varphi _{2,y}-\eta _{3,x}\varphi _{2,x})_{s_3} \big (\delta \varPhi ^2_3-(\delta \eta _3)(\varphi _{2,y})_{s_3}\big ) {\rm d}x \nonumber\\& \quad +\int _{\mathbb{R}}|\boldsymbol{\nabla }\varphi _2|^2_{s_2}\delta \eta _2\,{\rm d}x -\int _{\mathbb{R}}|\boldsymbol{\nabla }\varphi _2|^2_{s_3}\delta \eta _3\,{\rm d}x \end{align}

and

(3.15) \begin{align} \delta \left (\int _{\mathbb{R}}\int _{\eta _2}^{h_1+\eta _1}|\boldsymbol{\nabla }\varphi _1|^2\,{\rm d}y\,{\rm d}x\right ) & = 2 \int _{\mathbb{R}} (\varphi _{1,y}-\eta _{1,x}\varphi _{1,x})_{s_1} \big (\delta \varPhi ^1_1-(\delta \eta _1)(\varphi _{1,y})_{s_1}\big ) {\rm d}x \nonumber\\& \quad-2 \int _{\mathbb{R}} (\varphi _{1,y}-\eta _{2,x}\varphi _{1,x})_{s_2} \big (\delta \varPhi ^1_2-(\delta \eta _2)(\varphi _{1,y})_{s_2}\big ) {\rm d}x \nonumber\\& \quad +\int _{\mathbb{R}}|\boldsymbol{\nabla }\varphi _1|^2_{s_1}\delta \eta _1\,{\rm d}x -\int _{\mathbb{R}}|\boldsymbol{\nabla }\varphi _1|^2_{s_2}\delta \eta _2\,{\rm d}x. \end{align}

Collecting together the variations of $H$ (from (3.13)–(3.15) and (3.6)) with respect to $\eta _1$ , $\eta _2$ and $\eta _3$ , respectively, we find that

(3.16) \begin{align} \frac {\delta H}{\delta \eta _1}&=-\rho _1\,\eta _{1,t}(\varphi _{1,y})_{s_1}+\frac {\rho _1}{2}|\boldsymbol{\nabla }\varphi _1|^2_{s_1}+g\rho _1\eta _1 \nonumber\\&=-\rho _1\,\eta _{1,t}(\varphi _{1,y})_{s_1}-\rho _1(\varphi _{1,t})_{s_1} \nonumber\\&=-\frac {{\rm d}}{{\rm d}t}\left (\rho _1\,\varphi _1(x,h_1+\eta _1(x,t),t)\right ) \nonumber\\&=-\xi _{1,t}, \\[-28pt] \nonumber \end{align}

(3.17) \begin{align} \frac {\delta H}{\delta \eta _2}&=\rho _1\,\eta _{2,t}(\varphi _{1,y})_{s_2}-\frac {\rho _1}{2}|\boldsymbol{\nabla }\varphi _1|^2_{s_2}-\rho _2\,\eta _{2,t}(\varphi _{2,y})_{s_2}+\frac {\rho _2}{2}|\boldsymbol{\nabla }\varphi _2|^2_{s_2}+(\rho _2-\rho _1)g\eta _2 \nonumber \\&=-\rho _2\left ((\varphi _{2,t})_{s_2}+\eta _{2,t}(\varphi _{2,y})_{s_2}\right )+\rho _1\left ((\varphi _{1,t})_{s_2}+\eta _{2,t}(\varphi _{1,y})_{s_2}\right ) \nonumber\\&=-\frac {{\rm d}}{{\rm d}t}\left (\rho _2\,\varphi _2(x,\eta _2(x,t),t)\right )+\frac {{\rm d}}{{\rm d}t}\left (\rho _1\,\varphi _1(x,\eta _2(x,t),t)\right ) \nonumber\\&=-\xi _{2,t} \\[2pt] \nonumber \end{align}

and

(3.18) \begin{align} \frac {\delta H}{\delta \eta _3}&=\rho _2\,\eta _{3,t}(\varphi _{2,y})_{s_3}-\frac {\rho _2}{2}|\boldsymbol{\nabla }\varphi _2|^2_{s_3}-\rho _3\,\eta _{3,t}(\varphi _{3,y})_{s_3}+\frac {\rho _3}{2}|\boldsymbol{\nabla }\varphi _3|^2_{s_3}+(\rho _3-\rho _2)g\eta _3 \nonumber \\&=-\rho _3\left ((\varphi _{3,t})_{s_3}+\eta _{3,t}(\varphi _{3,y})_{s_3}\right )+\rho _2\left ((\varphi _{2,t})_{s_3}+\eta _{3,t}(\varphi _{2,y})_{s_3}\right ) \nonumber \\&=-\frac {{\rm d}}{{\rm d}t}\left (\rho _3\,\varphi _3(x,-h_3+\eta _3(x,t),t)\right )+\frac {{\rm d}}{{\rm d}t}\left (\rho _2\,\varphi _2(x,-h_3+\eta _3(x,t),t)\right )\nonumber \\&=-\xi _{3,t}. \end{align}

Moreover, glancing at (2.17) and utilising (3.13)–(3.15), we infer that

(3.19) \begin{align} \frac {\delta H}{\delta \xi _1}=\eta _{1,t},\quad \frac {\delta H}{\delta \xi _2}=\eta _{2,t},\quad \frac {\delta H}{\delta \xi _3}=\eta _{3,t}. \end{align}

We summarise the conclusions regarding the variations of $H$ computed in (3.16), (3.17), (3.18) and (3.19) in the following theorem.

Theorem 3.3. The governing equations and boundary conditions (2.3)–(2.8) possess the Hamiltonian formulation

(3.20) \begin{align} \begin{cases} \dfrac {\delta H}{\delta \eta _1}=-\xi _{1,t},\quad \dfrac {\delta H}{\delta \xi _1}=\eta _{1,t},\\[9pt] \dfrac {\delta H}{\delta \eta _2}=-\xi _{2,t},\quad \dfrac {\delta H}{\delta \xi _2}=\eta _{2,t},\\[9pt] \dfrac {\delta H}{\delta \eta _3}=-\xi _{3,t},\quad \dfrac {\delta H}{\delta \xi _3}=\eta _{3,t}. \end{cases} \end{align}

3.4. The DN operators

The DN operator corresponding to the uppermost layer is defined as

(3.21) \begin{align} \begin{pmatrix}G_{11} & G_{12}\\ G_{21} & G_{22}\end{pmatrix}\begin{pmatrix}(\varphi _{1})_{s_2}\\ (\varphi _{1})_{s_1}\end{pmatrix}= \begin{pmatrix} -\sqrt {1+\eta _{2,x}^2}\left (\dfrac{\partial \varphi _1}{\partial \boldsymbol{n}_2} \right )_{s_2} \\[6pt] \sqrt {1+\eta _{1,x}^2}\left (\dfrac{\partial \varphi _1}{\partial \boldsymbol{n}_1} \right )_{s_1} \end{pmatrix} \!, \end{align}

where $\boldsymbol{n}_2={(-\eta _{2,x},1)}/{\sqrt {1+\eta _{2,x}^2}}$ and $\boldsymbol{n}_1={(-\eta _{1,x},1)}/{\sqrt {1+\eta _{1,x}^2}}$ are the unit normal vectors to the upper interface $y=\eta _2(x,t)$ and to the top surface $y=-h_1+\eta _1(x,t)$ , respectively. The DN operator corresponding to the middle layer is

(3.22) \begin{align} \begin{pmatrix}\varGamma _{11} & \varGamma _{12}\\ \varGamma _{21} &\varGamma _{22}\end{pmatrix}\begin{pmatrix}(\varphi _2)_{s_2}\\ (\varphi _2)_{s_3}\end{pmatrix}= \begin{pmatrix} \sqrt {1+\eta _{2,x}^2}\left (\dfrac{\partial \varphi _2}{\partial \boldsymbol{n}_2} \right )_{s_2} \\[6pt] -\sqrt {1+\eta _{3,x}^2}\left (\dfrac{\partial \varphi _2}{\partial \boldsymbol{n}_3} \right )_{s_3} \end{pmatrix} \!, \end{align}

where $\boldsymbol{n}_3={(-\eta _{3,x},1)}/{\sqrt {1+\eta _{3,x}^2}}$ denotes the outward pointing unit normal vector corresponding to the lower interface $y=-h_3+\eta _3(x,t)$ . Moreover, the DN operator corresponding to the bottom layer is

(3.23) \begin{align} G(\eta _3)(\varphi _3)_{s_3}=\sqrt {1+\eta _{3,x}^2}\left ({\partial \varphi _3}/{\partial \boldsymbol{n}_3}\right )_{s_3}=(\boldsymbol{\nabla }\varphi _3)_{s_3}\boldsymbol{\cdot }(-\eta _{3,x},1). \end{align}

These DN operators for each domain $\varOmega_1,\varOmega_2,\varOmega_3$ could be obtained recursively in the form of asymptotic series in $\eta _1, \eta _2, \eta_3$ following the technique described in the Appendix of Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005); see also Appendix A of Ivanov, Martin & Todorov (Reference Ivanov, Martin and Todorov2022).

Appealing to the divergence theorem (Evans Reference Evans2010), the kinetic energy of the system (3.5) can be rewritten as

(3.24) \begin{align} {\rm E}_K & = \frac {\rho _1}{2}\int (\varphi _1)_{s_1}\left (\boldsymbol{\nabla }\varphi _1\right )_{s_1}\boldsymbol{\cdot }(-\eta _{1,x},1)\,{\rm d}x -\frac {\rho _1}{2}\int (\varphi _1)_{s_2}\left (\boldsymbol{\nabla }\varphi _1\right )_{s_2}\boldsymbol{\cdot }(-\eta _{2,x},1)\,{\rm d}x \nonumber\\& \quad + \frac {\rho _2}{2}\int (\varphi _2)_{s_2}\left (\boldsymbol{\nabla }\varphi _2\right )_{s_2}\boldsymbol{\cdot }(-\eta _{2,x},1)\,{\rm d}x -\frac {\rho _2}{2}\int (\varphi _2)_{s_3}\left (\boldsymbol{\nabla }\varphi _2\right )_{s_3}\boldsymbol{\cdot }(-\eta _{3,x},1)\,{\rm d}x \nonumber\\& \quad +\frac {\rho _3}{2}\int (\varphi _3)_{s_3}\left (\boldsymbol{\nabla }\varphi _3\right )_{s_3}\boldsymbol{\cdot }(-\eta _{3,x},1)\,{\rm d}x. \end{align}

With the help of the DN operators, the kinetic energy can be written as

(3.25) \begin{align} {\rm E}_K & = \frac {\rho _1}{2}\int (\varphi _1)_{s_1}\left [G_{21}(\varphi _1)_{s_2}+G_{22}(\varphi _1)_{s_1}\right ]\,{\rm d}x \nonumber\\& \quad +\frac {\rho _1}{2}\int (\varphi _1)_{s_2}\left [G_{11}(\varphi _1)_{s_2}+G_{12}(\varphi _1)_{s_1}\right ]\,{\rm d}x \nonumber\\& \quad +\frac {\rho _2}{2}\int (\varphi _2)_{s_2}\left [\varGamma _{11}(\varphi _2)_{s_2} +\varGamma _{12}(\varphi _2)_{s_3}\right ]\,{\rm d}x \nonumber\\& \quad +\frac {\rho _2}{2}\int (\varphi _2)_{s_3}\left [\varGamma _{21}(\varphi _2)_{s_2} +\varGamma _{22}(\varphi _2)_{s_3}\right ]\,{\rm d}x \nonumber\\& \quad +\frac {\rho _3}{2}\int (\varphi _3)_{s_3} G(\eta _3)(\varphi _3)_{s_3}\,{\rm d}x. \end{align}

Attempting to simplify (3.25) we note first that from the kinematic condition at the upper interface $y=\eta _2(x,t)$ , we have

(3.26) \begin{align} \eta _{2,t}=(\boldsymbol{\nabla }\varphi _2)_{s_2}\boldsymbol{\cdot }(\eta _{2,x},-1)=(\boldsymbol{\nabla }\varphi _1)_{s_2}\boldsymbol{\cdot }(\eta _{2,x},-1), \end{align}

which, in terms of the DN operators, becomes

(3.27) \begin{align} \eta _{2,t}=-(\varGamma _{11}(\varphi _2)_{s_2}+\varGamma _{12}(\varphi _2)_{s_3})=G_{11}(\varphi _1)_{s_2}+G_{12}(\varphi _1)_{s_1}. \end{align}

Moreover,

(3.28) \begin{align} \eta _{3,t}=(\varphi _{2,y}-\varphi _{2,x}\eta _{3,x})_{s_3}=(\varphi _{3,y}-\varphi _{3,x}\eta _{3,x})_{s_3}, \end{align}

which implies

(3.29) \begin{align} \eta _{3,t}=G(\eta _3)(\varphi _3)_{s_3}=- \big[\varGamma _{21}(\varphi _2)_{s_2} +\varGamma _{22}(\varphi _2)_{s_3}\big]. \end{align}

Recalling that

(3.30) \begin{align} \xi _1=\rho _1(\varphi _1)_{s_1},\quad\xi _2=\rho _2 (\varphi _2)_{s_2}-\rho _1 (\varphi _1)_{s_2},\quad\xi _3=\rho _3 (\varphi _3)_{s_3}-\rho _2(\varphi _2)_{s_3}, \end{align}

and employing (3.27) and (3.29), we see that the kinetic energy can be simplified as

(3.31) \begin{align} \begin{split} {E}_{K}={}&\frac {1}{2}\int \xi _1\left [G_{21}(\varphi _1)_{s_2}+G_{22}(\varphi _1)_{s_1}\right ]\,{\rm d}x \\ &{}+\frac {1}{2}\int \xi _2\left [\varGamma _{11}(\varphi _2)_{s_2} +\varGamma _{12}(\varphi _2)_{s_3}\right ]\,{\rm d}x\\ &{}+\frac {1}{2}\int \xi _3\, G(\eta _3)(\varphi _3)_{s_3}\,{\rm d}x. \end{split} \end{align}

Moreover, from (3.27) and (3.29), we also obtain the system

(3.32) \begin{align} \xi _2 &=\rho _2 (\varphi _2)_{s_2}+\rho _1 G_{11}^{-1}[\varGamma _{11}(\varphi _2)_{s_2} +\varGamma _{12}(\varphi _2)_{s_3}+\rho _1^{-1}G_{12}\xi _1], \\[-9pt] \nonumber \end{align}

(3.33) \begin{align} \xi _3&=-\rho _2 (\varphi _2)_{s_3}-\rho _3\, G^{-1}(\eta _3)\,[\varGamma _{21}(\varphi _2)_{s_2} +\varGamma _{22}(\varphi _2)_{s_3}], \\[9pt] \nonumber \end{align}

which can be rearranged as

(3.34) \begin{align} \begin{split} (\rho _2+\rho _1 G_{11}^{-1}\varGamma _{11})(\varphi _2)_{s_2}+(\rho _1 G_{11}^{-1}\varGamma _{12})(\varphi _2)_{s_3} &=\xi _2-G_{11}^{-1}G_{12}\xi _1,\\ (-\rho _3\, G^{-1}(\eta _3)\,\varGamma _{21})(\varphi _2)_{s_2} -(\rho _2 +\rho _3\, G^{-1}(\eta _3)\,\varGamma _{22})(\varphi _2)_{s_3} &=\xi _3, \end{split} \end{align}

where the superscript $-1$ in the previous system (and in the subsequent equations) is used to refer to the inverse of the operator in question.

Taking into account the fact that (3.34) is a linear system of equations with operator-valued coefficients, we obtain

(3.35) \begin{align} \begin{split} &(\varphi _2)_{s_2}=\mathcal{D}^{-1}\big ( \rho _1 G \xi _3 + (\rho _2 G+\rho _3\varGamma _{22})\varGamma _{12}^{-1}(G_{11}\xi _2-G_{12}\xi _1) \big )\equiv \varPhi _2^2(\eta _i, \xi _i), \\ &(\varphi _2)_{s_3}=-(\mathcal{D}^*)^{-1}\big ( \rho _3(G_{11}\xi _2-G_{12}\xi _1) + (\rho _2 G_{11}+\rho _1\varGamma _{11})\varGamma _{21}^{-1} G \xi _3 \big ), \end{split} \end{align}

where $\mathcal{D}:= ( (\rho _2 G+\rho _3\varGamma _{22})\varGamma _{12}^{-1}(\rho _2 G_{11}+\rho _1\varGamma _{11})-\rho _1 \rho _3 \varGamma _{21} )$ , and the star $(*)$ notation denotes the operation of conjugation. (The conjugation is with respect to the standard real bilinear form $(f,g)=\int _{\mathbb{R}} f(x)\, g(x)\, {\rm d}x$ .) We used also the fact that $G_{\textit{ij}}=G^*_{\textit{ji}}$ and $\varGamma _{\textit{ij}}=\varGamma ^*_{\textit{ji}}$ ; see e.g. Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005). Moreover, from (3.30), we also have

(3.36) \begin{align} (\varphi _1)_{s_2} & =\frac {1}{\rho _1}\mathcal{D}^{-1}\big ( \rho _1 (\rho _2 G \xi _3 +\rho _3\varGamma _{21} \xi _2)- (\rho _2 G+\rho _3\varGamma _{22})\varGamma _{12}^{-1}(\rho _2 G_{12}\xi _1+\rho _1 \varGamma _{11}\xi _2) \big ) \nonumber \\ & \equiv \varPhi ^1_2(\eta _i, \xi _i), \nonumber \\(\varphi _3)_{s_3}& =-(\mathcal{D}^*)^{-1}\big ( \rho _2 G_{11}\xi _2-\rho _2 G_{12}\xi _1+\rho _1 \varGamma _{12}\xi _3 {-} (\rho _2 G_{11}+\rho _1\varGamma _{11})\varGamma _{21}^{-1} \varGamma _{22} \xi _3 \big ). \end{align}

The previous expressions of the surface values of the potentials together with (3.6) allow us to express the Hamiltonian as a functional of the quantities $\eta _i,\xi _i\ (i=1,2,3)$ , i.e. we have

(3.37) \begin{align} \begin{split} H={}&{{E}_{{K}}(\xi _i)}+V(\eta _i)\\ ={}&\frac {1}{2}\int \xi _1 A_{11} \xi _1 \,{\rm d}x + \frac {1}{2}\int \xi _2 A_{22} \xi _2 \,{\rm d}x+\frac {1}{2}\int \xi _3 A_{33} \xi _3 \,{\rm d}x \\ &{}+ \int \xi _1 A_{12} \xi _2 \, {\rm d}x + \int \xi _1 A_{13} \xi _3 \, {\rm d}x + \int \xi _2 A_{23} \xi _3 \, {\rm d}x \\ &{}+\frac {1}{2} \int \big ( g\rho _1 \eta _1^2 +g(\rho _2-\rho _1) \eta _2^2 + g(\rho _3-\rho _2) \eta _3^2\big ) \, {\rm d}x, \end{split} \end{align}

where the last integral is the potential energy $V$ , and the operators $A_{\textit{ij}}\ (i,j=1,2,3)$ , written in terms of the DN operators, are

(3.38) \begin{align} &A_{11}=\frac {1}{\rho _1}G_{22}- \frac {\rho _2}{\rho _1} G_{21}\mathcal{D}^{-1}(\rho _2 G +\rho _3 \varGamma _{22})\varGamma _{12}^{-1}G_{12} , \nonumber\\&A_{22}= \varGamma _{11} \mathcal{D}^{-1} (\rho _2 G + \rho _3 \varGamma _{22})\varGamma _{12}^{-1}G_{11}-\rho _3 G_{11} \mathcal{D}^{-1} \varGamma _{21}, \nonumber\\& A_{33}= \big ( \varGamma _{22} \varGamma _{12}^{-1}(\rho _2 G_{11}+\rho _1 \varGamma _{11}) -\rho _1 \varGamma _{21} \big )\mathcal{D}^{-1} G,\nonumber\\&A_{12}=\rho _3 G_{21} \mathcal{D}^{-1} \varGamma _{21} - \frac {1}{2} G_{21}\big ( \mathcal{D}^{-1}(\rho _2 G + \rho _3 \varGamma _{22})\varGamma _{12}^{-1}+ \varGamma _{21}^{-1} (\rho _2 G + \rho _3 \varGamma _{22}) (\mathcal{D}^*)^{-1} \big )\varGamma _{11},\nonumber\\&A_{13}=\frac {1}{2}\rho _2 G_{21} \big (\mathcal{D}^{-1}+(\mathcal{D}^*)^{-1}\big )G,\nonumber\\&A_{23}=\frac {1}{2}\big ( (\rho _1 \varGamma _{11}-\rho _2 G_{11}) \mathcal{D}^{-1} - \varGamma _{12} (\mathcal{D}^*)^{-1} (\rho _2 G_{11}+\rho _1 \varGamma _{11})\varGamma _{21}^{-1}\big ) G. \end{align}

The Hamiltonian equations can be represented, via the DN operators, with the help of the quantities $(\varphi _{1,x}(x,y,t))_{s_1}$ and $(\varphi _{1,y}(x,y,t))_{s_1}$ . From (3.37) and because the operators $A_{\textit{ij}}$ do not depend on $\xi _i$ , we obtain (in terms of the Hamiltonian variables) the evolution equation

(3.39) \begin{align} \eta _{1,t}=\frac {\delta H}{\delta \xi _1}=A_{11}\xi _1+A_{12}\xi _2+A_{13}\xi _3:=F_1(\eta _i,\xi _i). \end{align}

In addition, by (3.37) and (3.20), $(\varphi _{1,x}(x,y,t))_{s_1}$ and $(\varphi _{1,y}(x,y,t))_{s_1}$ satisfy the system

(3.40) \begin{align} \begin{split} &(\varphi _{1,x})_{s_1}+(\varphi _{1,y})_{s_1} \eta _{1,x}= \frac {1}{\rho _1}\xi _{1,x},\\ &(\varphi _{1,y})_{s_1}-(\varphi _{1,x})_{s_1}\eta _{1,x}=\eta _{1,t}=F_1(\eta _i,\xi _i). \end{split} \end{align}

Solving the system, we find

(3.41) \begin{align} \begin{split} & (\varphi _{1,x})_{s_1}= \frac {\xi _{1,x}-\rho _1\eta _{1,x} F_1}{\rho _1 \big(1+\eta _{1,x}^2 \big)} , \\&(\varphi _{1,y})_{s_1}= \frac {\eta _{1,x}\xi _{1,x}+\rho _1 F_1}{\rho _1 \big(1+\eta _{1,x}^2 \big)}, \end{split} \end{align}

which allow us to write (2.18) through the Hamiltonian variables. More precisely, we have

(3.42) \begin{align} \xi _{1,t} + \frac {\xi _{1,x}^2 -2\rho _1\eta _{1,x}\xi _{1,x}F_1- \rho _1^2F_1^2}{2\rho _1 \big(1+\eta _{1,x}^2 \big)} +g\rho _1\eta _1 =0. \end{align}

The counterpart equation for $\eta _1$ is of course (3.39),

(3.43) \begin{align} \eta _{1,t} = F_1(\xi _i,\eta _i). \end{align}

From (3.37), where the operators $A_{\textit{ij}}$ do not depend on $\xi _i$ , we have

(3.44) \begin{align} \eta _{2,t}=\frac {\delta H}{\delta \xi _2}= A_{12}^* \xi _1+A_{22}\xi _2+A_{23} \xi _3 =: F_2(\eta _i,\xi _i). \end{align}

Thus $F_2(\eta _i,\xi _i)$ is defined through the canonical variables. Next, from (3.35) we have $(\varphi _2)_{s_2}\equiv \varPhi _2^2(\eta _i, \xi _i)$ as a known quantity, defined through the canonical variables. Since

(3.45) \begin{align} \frac {{\rm d}}{{\rm d}x} \varPhi _2^2= (\varphi _{2,x} )_{s_2}+(\varphi _{2,y})_{s_2}\eta _{2,x}, \end{align}

and from (2.12)

(3.46) \begin{align} F_2=\eta _{2,t}=(\varphi _{2,y})_{s_2}-(\varphi _{2,x} )_{s_2}\eta _{2,x}, \end{align}

we have a system

(3.47) \begin{align} & (\varphi _{2,x} )_{s_2}+ \eta _{2,x}(\varphi _{2,y})_{s_2}= (\varPhi ^2_2 )_x, \\[-9pt] \nonumber \end{align}

(3.48) \begin{align} -&\eta _{2,x} (\varphi _{2,x} )_{s_2}+(\varphi _{2,y})_{s_2}=F_2, \\[9pt] \nonumber \end{align}

which gives

(3.49) \begin{align} (\varphi _{2,x} )_{s_2}&= \frac {(\varPhi ^2_2)_x-\eta _{2,x}F_2}{1+\eta _{2,x}^2 }, \\[-12pt] \nonumber \end{align}

(3.50) \begin{align} (\varphi _{2,y})_{s_2}&= \frac {\eta _{2,x}(\varPhi ^2_2)_x +F_2}{1+\eta _{2,x}^2 }, \\[9pt] \nonumber \end{align}

and similarly,

(3.51) \begin{align} (\varphi _{1,x})_{s_2}&= \frac {(\varPhi ^1_2)_x-\eta _{2,x}F_2}{1+\eta _{2,x}^2 }, \\[-12pt] \nonumber \end{align}

(3.52) \begin{align} (\varphi _{1,y})_{s_2}&= \frac {\eta _{2,x}(\varPhi ^1_2)_x +F_2}{1+\eta _{2,x}^2 }, \\[9pt] \nonumber \end{align}

where $\varPhi ^1_2$ is defined in (3.36). Therefore, $\xi _2$ satisfies the equation

(3.53) \begin{align} \xi _{2,t}+\frac {\rho _2 (\varPhi ^2_{2,x} )^2-\rho _1 (\varPhi _{2,x}^1 )^2-2F_2\xi _{2,x}\eta _{2,x}+(\rho _1-\rho _2)F_2^2 }{2 \big(1+\eta _{2,x}^2 \big)}+(\rho _2-\rho _1)g \eta _2 =0. \end{align}

Similarly, for the last pair of equations, we obtain

(3.54) \begin{align} \eta _{3,t}=\frac {\delta H}{\delta \xi _3}=A_{13}^{*}\xi _1 +A_{23}^{*}\xi _2+A_{33}\xi _3:=F_3(\eta _i,\xi _i) \end{align}

and

(3.55) \begin{align} \xi _{3,t}+\frac {\rho _3 (\varPhi ^3_{x} )^2-\rho _2 (\varPhi _{3,x}^2 )^2-2F_3\xi _{3,x}\eta _{3,x}+(\rho _2-\rho _3)F_3^2 }{2\big(1+\eta _{3,x}^2 \big)}+(\rho _3-\rho _2)g \eta _3 =0. \end{align}

We note that (3.42), (3.53), (3.55), (3.43), (3.44), (3.54) render a closed system of six equations for the six variables $(\eta _i,\xi _i)$ , $i=1,2,3$ . All quantities in these equations are given in terms of the DN operators. Most importantly, (3.42), (3.53), (3.55), (3.43), (3.44), (3.54) also suggest the pattern for more superimposed layers.

4.1. The derivation of the dispersion relation

A significant role in the derivation of the dispersion relation will be played by the DN operators introduced in (3.21), (3.22) and (3.23). The DN operators are analytic in their dependence on $(\eta _1,\eta _2)$ , $(\eta _2,\eta _3)$ and $\eta _3$ , respectively, and have convergent Taylor series expansions

(4.1) \begin{align} \begin{pmatrix} G_{11}(\eta _1,\eta _2) & G_{12}(\eta _1,\eta _2)\\[5pt] G_{21}(\eta _1,\eta _2) & G_{22}(\eta _1,\eta _2) \end{pmatrix}=\sum _{k_1,k_2=0}^{\infty }\begin{pmatrix} G_{11}^{(k_1 k_2)}(\eta _1,\eta _2) & G_{12}^{(k_1 k_2)}(\eta _1,\eta _2)\\[7pt] G_{21}^{(k_1 k_2)}(\eta _1,\eta _2) & G_{22}^{(k_1 k_2)}(\eta _1,\eta _2) \end{pmatrix} \!, \\[-28pt] \nonumber \end{align}

(4.2) \begin{align} \begin{pmatrix} \varGamma _{11}(\eta _2,\eta _3) & \varGamma _{12}(\eta _2,\eta _3)\\[5pt] \varGamma _{21}(\eta _2,\eta _3) & \varGamma _{22}(\eta _2,\eta _3) \end{pmatrix}=\sum _{p_1,p_2=0}^{\infty }\begin{pmatrix} \varGamma _{11}^{(p_2 p_3)}(\eta _2,\eta _3) & \varGamma _{12}^{(p_2 p_3)}(\eta _2,\eta _3)\\[7pt] \varGamma _{21}^{(p_2 p_3)}(\eta _2,\eta _3) & \varGamma _{22}^{(p_2 p_3)}(\eta _2,\eta _3) \end{pmatrix}\!, \\[-28pt] \nonumber \end{align}

(4.3) \begin{align} G(\eta _3)=\sum _{l=0}^{\infty }G^{(l)}(\eta _3), \end{align}

with each linear operator $G_{\textit{ij}}^{(k_1 k_2)}(\eta _1,\eta _2)$ homogeneous of degree $k_1$ in $\eta _1$ , and of degree $k_2$ in $\eta _2$ , for $i,j=1,2$ , each linear operator $\varGamma _{\textit{ij}}^{(p_2 p_3)}(\eta _2,\eta _3)$ homogeneous of degree $p_2$ in $\eta _2$ , and of degree $p_3$ in $\eta _3$ , for $i,j=1,2$ , and each linear operator $G^{(l)}(\eta _3)$ homogeneous of degree $l$ in $\eta _3$ ; cf. Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005, Reference Craig, Guyenne and Sulem2012). Each of the operators $G_{\textit{ii}}^{(k_1 k_2)}(\eta _1,\eta _2)$ with $i=1,2$ and $k_1,k_2\geq 0$ , $\varGamma _{\textit{ii}}^{(p_2 p_3)}(\eta _2,\eta _3)$ with $i=1,2$ and $p_1,p_2\geq 0$ , and $G^{(l)}(\eta _3)$ with $l\geq 0$ , is self-adjoint.

Remark 4.1. The DN operators can be understood as certain pseudo-differential operators. More precisely, letting $m$ be a complex-valued function of a real variable whose derivatives of any order have polynomial growth, and setting

(4.4) \begin{align} D:=-i\partial _x, \end{align}

we define

(4.5) \begin{align} (m(D)f)(x):=\frac {1}{2\pi }\iint \mathrm{e}^{\mathrm{i}k(x-y)}\,m(k)\,f(y)\,{\rm d}y\,{\rm d}k. \end{align}

The operator $m(D)$ is called a Fourier multiplier operator and maps $\mathcal{S}(\mathbb{R})$ into $\mathcal{S}(\mathbb{R})$ . Furthermore,

1. (i) $m(D)$ extends to a self-adjoint operator in $L^2(\mathbb{R})$ if and only if $m$ is real valued (cf. Reed & Simon Reference Reed and Simon1975)

2. (ii) $m(D)$ is bounded if and only if $m\in L^{\infty }(\mathbb{R})$ .

Then the leading-order terms of the DN operators are given (cf. Craig et al. Reference Craig, Guyenne, Nicholls and Sulem2005) by means of Fourier multipliers as

(4.6) \begin{align} \begin{pmatrix} G_{11}^0 & G_{12}^0\\[5pt] G_{21}^0 & G_{22}^0 \end{pmatrix} = \begin{pmatrix} D\coth (h_1 D) & -D\,{\textrm {csch}}(h_1 D)\\[5pt] -D\,{\textrm {csch}}(h_1 D) & D\coth (h_1 D) \end{pmatrix}\!, \\[-28pt] \nonumber \end{align}

(4.7) \begin{align} \begin{pmatrix} \varGamma _{11}^0 & \varGamma _{12}^0\\[5pt] \varGamma _{21}^0 & \varGamma _{22}^0 \end{pmatrix} =\begin{pmatrix}D\coth (h_3 D) & -D\,{\textrm {csch}}(h_3 D)\\[5pt] -D\,{\textrm {csch}}(h_3 D) & D\coth (h_3 D) \end{pmatrix} \\[2pt] \nonumber \end{align}

and

(4.8) \begin{align} G^0=D\tanh (h_b D), \end{align}

where $h_b:=h-h_3$ is the mean depth of the bottom layer, for $h_3$ the mean depth of the middle layer.

Under consideration will be monochromatic solutions of the leading-order linear equations, meaning that all components of these solutions are proportional to $\exp (ikx)$ , where $k:={2\pi }/{L}$ represents the wavenumber, $L$ being the corresponding wavelength. The latter means that each component is an eigenfunction of the operator $D$ with eigenvalue $k$ . These facts allow us to replace the operator $D$ by its eigenvalue $k$ when acting on the previously mentioned monochromatic solutions. Consequently, leading-order approximations of the operators $A_{\textit{ij}}(k)$ $(i,j=1,2,3)$ , which are denoted by $A_{\textit{ij}}^0(k)$ $(i,j=1,2,3)$ , acquire the form

(4.9) \begin{align} &A_{11}^0(k)= \frac {k}{\rho _1}\coth (kh_1)-\frac {k\rho _2\big ( \rho _3\cosh (kh_3)\cosh (kh_b)+\rho _2 \sinh (kh_3)\sinh (k h_b) \big )}{\rho _1 \sinh (kh_1)\, \Delta (k)} , \nonumber\\ &A_{22}^0(k)= \frac {k}{\Delta (k)} \cosh (k h_1) \big (\rho _2 \cosh (k h_3) \sinh (k h_b) + \rho _3\sinh (k h_3) \cosh (k h_b) \big ), \nonumber\\ & A_{33}^0(k)= \frac {k}{\Delta (k)}\sinh (k h_b) \big (\rho _1 \sinh (k h_1) \sinh (k h_3) + \rho _2 \cosh (k h_3) \cosh (k h_1)\big ),\nonumber\\ &A_{12}^0(k)= \frac {k}{\Delta (k)}\big ( \rho _2\cosh (kh_3)\sinh (kh_b)+\rho _3 \sinh (kh_3)\cosh (k h_b) \big ) ,\nonumber\\ &A_{13}^0(k)= \frac { k \rho _2 }{ \Delta (k)} \sinh (k h_b),\nonumber\\ &A_{23}^0(k)= \frac {k \rho _2 }{\Delta (k)} \sinh (k h_b) \cosh (k h_1), \end{align}

where

(4.10) \begin{align} \begin{split} \Delta (k)={}& \rho _2\big (\rho _3 \cosh (k h_3) \cosh (k h_b)+\rho _2 \sinh (kh_3)\sinh (kh_b)\big )\cosh (kh_1) \\ & {}+ \rho _1\big (\rho _2 \cosh (k h_3) \sinh (k h_b)+\rho _3 \sinh (kh_3)\cosh (kh_b)\big )\sinh (kh_1). \\ \end{split} \end{align}

Note that

(4.11) \begin{align} \begin{split} \Delta (k)\gt{}& \rho _1^2\big ( \cosh (k h_3) \cosh (k h_b)+ \sinh (kh_3)\sinh (kh_b)\big )\cosh (kh_1) \\ &{}+ \rho _1^2\big (\cosh (k h_3) \sinh (k h_b)+ \sinh (kh_3)\cosh (kh_b)\big )\sinh (kh_1) \\ ={}&\rho _1^2\big (\cosh (kh_3+kh_b)\cosh (kh_1)+\sinh (kh_3+kh_b)\sinh (kh_1)\big )\\ ={}&\rho _1^2\cosh (k(h_3+h_b+h_1))=\rho _1^2\cosh (k(h+h_1))\gt 0. \end{split} \end{align}

We observe that in the leading order, all $A^0_{\textit{ij}}(k)$ are even functions of $k$ , therefore the operators $A^0_{\textit{ij}}(D)$ are self-conjugate operators. Thus from (3.39), (3.44) and (3.54), we have

(4.12) \begin{align} \begin{split} \eta _{1,t}&=A^0_{11}(D)\,\xi _1+A^0_{12}(D)\,\xi _2+A^0_{13}(D)\,\xi _3,\\ \eta _{2,t}&=A_{12}^0 (D)\,\xi _1+A^0_{22}(D)\,\xi _2+A^0_{23}(D)\, \xi _3,\\ \eta _{3,t}&=A_{13}^{0}(D)\,\xi _1 +A_{23}^{0}(D)\,\xi _2+A^0_{33}(D)\,\xi _3. \end{split} \end{align}

On the other hand, retaining only the linear terms in (3.42), (3.53) and (3.55), we obtain the equations

(4.13) \begin{align} \begin{split} \xi _{1,t}&=-g\rho _1\eta _1,\\ \xi _{2,t}&=-g(\rho _2-\rho _1)\eta _2,\\ \xi _{3,t}&=-g(\rho _3-\rho _2)\eta _3. \end{split} \end{align}

Denoting now $\mathfrak{m}=(\eta _1,\eta _2,\eta _3,\xi _1,\xi _2,\xi _3)^{\rm T}$ , we will solve the linear system (4.12)–(4.13) by utilising the Fourier transform: setting $\hat {f}(k)=({1}/{\sqrt {2\pi }})\int _{\mathbb{R}}\mathrm{e}^{-\mathrm{i}kx}f(x)\,{\rm d}x$ for $f\in \mathcal{S}(\mathbb{R})$ , in each component of $\mathfrak{m}$ , we see (after taking into account (3.20), (3.37) and (4.13)) that the system (3.20) is transformed for any fixed $k\in \mathbb{R}$ into the linear autonomous system of ordinary differential equations

(4.14) \begin{align} \partial _t\hat {\mathfrak{m}}(k,t)=\mathcal{M}(k)\,\hat {\mathfrak{m}}(k,t), \end{align}

where

(4.15) \begin{align} \mathcal{M}(k)=\begin{pmatrix} 0 & 0 & 0 & A^0_{11}(k) & A^0_{12}(k) & A^0_{13}(k)\\[5pt] 0 & 0 & 0 & A^0_{12}(k) & A^0_{22}(k) & A^0_{23}(k)\\[5pt] 0 & 0 & 0 & A^0_{13}(k) & A^0_{23}(k) & A^0_{33}(k)\\[5pt] -g\rho _1 & 0 & 0 & 0 & 0 & 0\\[5pt] 0 & g(\rho _1-\rho _2) & 0 & 0 & 0 & 0\\[5pt] 0 & 0 & g(\rho _2-\rho _3) & 0 & 0 &0 \end{pmatrix}\!. \end{align}

Remark 4.2. According to Constantin & Ivanov (Reference Constantin and Ivanov2019), the unique solution to (4.14) is given by $\hat {\mathfrak{m}}(k,t)=e^{\mathcal{M}(k)\,t}\,\hat {\mathfrak{m}}_0(k)$ (where $\mathfrak{m}_0(k)=({1}/{\sqrt {2\pi }})\int _{\mathbb{R}}\mathfrak{m}_0(x)\,\mathrm{e}^{\mathrm{i}kx}\,{\rm d}x$ are the initial data) and corresponds (via the inverse Fourier transform) to the solution

(4.16) \begin{align} \mathfrak{m}(x,t)=\frac {1}{2\pi }\int _{\mathbb{R}}e^{\mathcal{M}(k)\,t}\,\hat {\mathfrak{m}}_0(k)\,\mathrm{e}^{\mathrm{i}kx}\,{\rm d}k,\quad t\geq 0, \end{align}

of (3.20) with initial data $\mathfrak{m}_0\in \mathcal{S}(\mathbb{R})$ . We note now that if $\lambda (k)=-ikc\ (c\neq 0)$ is a purely imaginary eigenvalue of $\mathcal{M}(k)$ with eigenvector $\mathfrak{v}(k)\neq 0$ , then $e^{-ikc t}$ is an eigenvalue of $e^{\mathcal{M}(k)\,t}$ with eigenvector $\mathfrak{v}(k)$ . Therefore, we see from (4.16) that the solution $\mathfrak{m}(x,t)$ of (3.20) acquires the form

(4.17) \begin{align} \frac {1}{2\pi }\int _{\mathbb{R}} \mathrm{e}^{\mathrm{i}k(x-ct)}\,\hat {\mathfrak{m}}_0(k)\,{\rm d}k. \end{align}

This shows that a purely imaginary eigenvalue $\lambda (k)=-ikc$ of $\mathcal{M}(k)$ corresponds to the fundamental oscillation mode $\mathrm{e}^{\mathrm{i}k(x-ct)}\,\mathfrak{m}_0(k)$ , with frequency $ {|k|}/{2\pi }$ , and constant speed $c$ .

Searching now for solutions $\hat {\mathfrak{m}}=\hat {\mathfrak{m}}_0\, \mathrm{e}^{\mathrm{i}k(x-ct)}$ in (4.14) yields

(4.18) \begin{align} (\mathcal{M}(k)+ikcI_6)\hat {\mathfrak{m}}=0, \end{align}

which implies $\textrm {det} (\mathcal{M}(k)+ikcI_6)=0$ , since it is meaningful to require $\hat {\mathfrak{m}}\neq 0$ . Since for $\lambda \in \mathbb{R}$ it holds that (and accounting for $A^0_{\textit{ij}}:=A^0_{\textit{ji}}$ when $i\gt j$ )

(4.19) \begin{align} \begin{split} &\det(\mathcal{M}(k)-\lambda I_6)\\ &={}\lambda ^6-(b_1 A^0_{11}+b_2 A^0_{22}+b_3 A^0_{33})\lambda ^4\\ &\quad{}+[b_1 b_2(A^0_{11} A^0_{22}-(A^0_{12})^2)+b_2 b_3 (A^0_{22} A^0_{33}-(A^0_{23})^2)+b_1 b_3(A^0_{11} A^0_{33}-(A^0_{13})^2)]\lambda ^2\\ &\quad{}-b_1 b_2 b_3\det(A^0_{i,j})_{i,j=\overline {1,3}}, \end{split} \end{align}

we have that the wave speed $c$ satisfies the equation

(4.20) \begin{align} \begin{split} c^6&+\frac {b_1 A^0_{11}+b_2 A^0_{22}+b_3 A^0_{33}}{k^2}c^4\\ &+\frac {b_1 b_2(A^0_{11} A^0_{22}-(A^0_{12})^2)+b_2 b_3 (A^0_{22} A^0_{33}-(A^0_{23})^2)+b_1 b_3(A^0_{11} A^0_{33}-(A^0_{13})^2)}{k^4}c^2\\ &+\frac {b_1 b_2 b_3\textrm {det}(A^0_{\textit{ij}})_{i,j=\overline {1,3}} }{k^6}=0, \end{split} \end{align}

where $b_1=-g\rho _1,\ b_2=g(\rho _1-\rho _2),\ b_3=g(\rho _2-\rho _3)$ . To ease further calculations, we find it useful to introduce the notation

(4.21) \begin{align} \begin{split} &\varDelta _1(k):= \rho _2 \cosh (k h_3) \sinh (k h_b)+\rho _3 \sinh (kh_3)\cosh (kh_b) ,\\&\varDelta _2(k):= \rho _2 \sinh (kh_3)\sinh (kh_b)+ \rho _3 \cosh (k h_3) \cosh (k h_b),\\& \varDelta _3(k):=\rho _2 \cosh (k h_3)\cosh (k h_1) + \rho _1 \sinh (k h_1)\sinh (k h_3),\\&\varDelta _4(k):=\rho _1 \cosh (k h_3)\sinh (k h_1) + \rho _2 \sinh (k h_3)\cosh (k h_1),\\&\varDelta _5(k):=\rho _1 \cosh (k h_1)\sinh (k h_3) + \rho _2 \sinh (k h_1)\cosh (k h_3), \end{split} \end{align}

which entails

(4.22) \begin{align} \Delta (k)= \rho _1 \sinh (k h_1)\, \varDelta _1(k) + \rho _2\cosh (kh_1)\, \varDelta _2(k). \end{align}

Setting also

(4.23) \begin{align} \tilde {\Delta }(k)= \rho _3 \cosh (k h_b)\,\varDelta _3(k)+\rho _2 \sinh (k h_b)\, \varDelta _4(k), \end{align}

we have

(4.24) \begin{align} \begin{split} M_{12}:= &A^0_{11} A^0_{22}-(A^0_{12})^2=\frac {k^2\, \tilde {\Delta }(k)\, \varDelta _1(k)\sinh (k h_1)}{\rho _1\,\varDelta ^2(k)},\\ M_{13}:= &A^0_{11} A^0_{33}-(A^0_{13})^2= \frac {k^2\, \tilde {\Delta }(k)\, \varDelta _5(k) \sinh (k h_b)}{\rho _1\, \varDelta ^2(k)} ,\\ M_{23}:= & A^0_{22} A^0_{33}-(A^0_{23})^2 =\frac {k^2\, \tilde {\Delta }(k) \cosh (k h_1) \sinh (k h_b)\sinh (k h_3) }{\varDelta ^2(k)} . \end{split} \end{align}

The coefficients of (4.20) can now be written explicitly as

(4.25) \begin{align} \begin{split} a_4:={}& \frac {b_1 A^0_{11}+b_2 A^0_{22}+b_3 A^0_{33}}{k^2}\\ ={}& \frac {g\big ( (\rho _2-\rho _1) (\rho _3-\rho _2)\sinh (k h_1) \sinh (k h_b) \sinh (k h_3) - \rho _2 \rho _3 \sinh (k(h_1+h_b+h_3))\big )}{k\, \Delta (k)} ,\\ a_2:={} &\frac {b_1 b_2 M_{12}+b_2 b_3 M_{23}+b_1 b_3 M_{13}}{k^4} \\ ={} &\frac {g^2\, \tilde {\Delta }(k)}{k^2\, \varDelta ^2(k)}\big ((\rho _2-\rho _1)(\rho _3-\rho _2)\cosh (k h_1) \sinh (k h_b) \sinh (k h_3) \big . \\ & \big . \phantom {*******} +(\rho _2-\rho _1)\, \varDelta _1 \sinh (k h_1) + (\rho _3-\rho _2)\, \varDelta _5 \sinh (k h_b) \big ),\\ a_0: =&\frac {b_1 b_2 b_3\det(A^0_{\textit{ij}})_{i,j=\overline {1,3}} }{k^6}= -\frac {g^3(\rho _2-\rho _1)(\rho _3-\rho _2)\sinh (k h_1)\sinh (k h_b) \sinh (k h_3)}{k^3\, \Delta (k)}. \end{split} \end{align}

4.2. Short-wave limit

The short-wave limit for the solutions of (4.20) can be obtained by using the asymptotics $\sinh (k h_i) \sim \pm (1/2)\exp (kh_i)$ and $ \cosh (k h_i) \sim (1/2)\exp (k h_i)$ for $k\to \pm \infty $ . Thus

(4.26) \begin{align} a_4 & \sim - \frac {g}{|k|}\left (1+ \frac {\rho _2 - \rho _1}{\rho _2 + \rho _1}+ \frac {\rho _3-\rho _2}{\rho _3+\rho _2} \right )\! , \\[-12pt] \nonumber \end{align}

(4.27) \begin{align} a_2 & \sim \frac {g^2}{k^2} \left (\frac {\rho _2 - \rho _1}{\rho _2 + \rho _1}+ \frac {\rho _3-\rho _2}{\rho _3+\rho _2} + \frac {\rho _2 - \rho _1}{\rho _2 + \rho _1} \boldsymbol{\cdot }\frac {\rho _3-\rho _2}{\rho _3+\rho _2} \right )\!, \\[-12pt] \nonumber \end{align}

(4.28) \begin{align} a_0 & \sim - \frac {g^3}{|k|^3} \left (\frac {\rho _2 - \rho _1}{\rho _2 + \rho _1}\right ) \left (\frac {\rho _3-\rho _2}{\rho _3+\rho _2}\right )\!.\end{align}

Thus the roots of (4.20) are real and given by the expressions

(4.29) \begin{align} c_1^2 \sim \frac {g}{|k|} , \quad c_{2}^2 \sim \frac {g}{|k|} \left (\frac {\rho _2 - \rho _1}{\rho _2 + \rho _1}\right )\! , \quad c_{3}^2 \sim \frac {g}{|k|} \left (\frac {\rho _3 - \rho _2}{\rho _3 + \rho _2}\right )\!. \end{align}

5.1. Bounds for the roots of the dispersion relation

To study (5.1) in greater detail, we denote $c^2:=X$ , so we are led to consider the polynomial equation

(5.2) \begin{align} \begin{split} &P(X):=X^3-g H X^2\\[1em] &+g^2\left [\rho _1(\rho _2-\rho _1)\frac {h_1 h_3}{\rho _1\rho _2}+\rho _2(\rho _3-\rho _2)\frac {h_3 h_b }{\rho _2\rho _3} +\rho _1(\rho _3-\rho _1)\frac {h_1 h_b }{\rho _1\rho _3}\right ]X\\[1em] &-g^3\rho _1(\rho _2-\rho _1)(\rho _3-\rho _2)\frac {h_1 h_3 h_b }{\rho _1\rho _2\rho _3}=0. \end{split} \end{align}

The following result gives the interval where the roots of $P$ lie.

Proof. We readily note from above that the hypotheses $\rho _3\gt \rho _2\gt \rho _1$ and $h\gt h_3$ yield

(5.3) \begin{align} P(0)\lt 0. \end{align}

Furthermore, it is also easy to see that if $H_m\geq H=h+h_1$ , then we have

(5.4) \begin{align} \begin{split} P(gH_m)={}&g^3 H_m^2(H_m-H)+g^3 H_m\left [\rho _1(\rho _2-\rho _1)\frac {h_1 h_3}{\rho _1\rho _2}+\rho _1(\rho _3-\rho _1)\frac {h_1 h_b }{\rho _1\rho _3}\right ]\\ &{}+\frac {g^3(\rho _3-\rho _2)h_3 h_b }{\rho _2\rho _3}\underbrace {(\rho _2 H_m- (\rho _2-\rho _1)h_1)}_{\geq \rho _2 h +\rho _1 h_1} \gt 0. \end{split} \end{align}

Inequalities (5.3) and (5.4) yield that there is $X_1\in (0, gH )$ such that $P(X_1)=0$ . Close inspection of the coefficients of $P$ shows that $P(X)\lt 0$ for all $X\leq 0$ . This implies that all the real roots of $P$ lie in $(0, gH )$ .

Although explicit formulas are available for the roots of $P$ , the intricacy of these formulas makes them difficult to use. Instead, using a general result of Prodanov (Reference Prodanov2022a,Reference Prodanovb), we provide bounds for the roots, defined in terms of the coefficients of $P$ . First, we establish some notation. Let

(5.5) \begin{align} \mathcal{A} & :=-gH, \\[-12pt] \nonumber \end{align}

(5.6) \begin{align} \mathcal{B} &:=g^2\left [\frac {\rho _2-\rho _1}{\rho _2}h_1 h_3+\frac {\rho _3-\rho _2}{\rho _3}h_3 h_b +\frac {\rho _3-\rho _1}{\rho _3}h_1 h_b \right ]\!, \\[-12pt] \nonumber \end{align}

(5.7) \begin{align} \mathcal{C}& :=-g^3(\rho _2-\rho _1)(\rho _3-\rho _2)\frac {h_1 h_3 h_b }{\rho _2\rho _3} \\[9pt] \nonumber \end{align}

be the coefficients of the polynomial from (5.2).

The locations of the three roots of the polynomial $P(X)$ depend (cf. Prodanov Reference Prodanov2022a ,Reference Prodanovb) on the ratio $\mathcal{B}/\mathcal{A}^2$ , which will be analysed in the sequel. For a fixed value of $\mathcal{A}$ (i.e. for a fixed total depth $H$ ), the maximum of $\mathcal{B}$ can be found, for example, from the Lagrangian function

(5.8) \begin{align} L(h_1,h_3,h_b,\lambda )=\mathcal{B}(h_1,h_3,h_b) - \lambda (h_1+h_3+h_b - H). \end{align}

The maximal value of $\mathcal{B}$ is

(5.9) \begin{align} \mathcal{B}_{\textit{max}}= \frac {\rho _2 (\rho _3 -\rho _1) g^2 H^2}{\rho _2^2 + (3\rho _3-\rho _1)\rho _2 +\rho _3\rho _1}= \frac {\nu _2 (1-\nu _1 ) }{\nu _2^2 + (3-\nu _1)\nu _2 +\nu _1} \mathcal{A}^2 \end{align}

(where $\nu _1={\rho _1}/{\rho _3}\lt 1$ and $\nu _2= {\rho _2}/{\rho _3}\lt 1$ ), attained for

(5.10) \begin{align} h_1^0&= \frac {(\rho _2 + \rho _3)\rho _2}{\rho _2^2 + (3\rho _3-\rho _1)\rho _2 +\rho _3\rho _1 }H, \\[-12pt] \nonumber \end{align}

(5.11) \begin{align} h_3^0&= \frac {(\rho _3 - \rho _1)\rho _2} {\rho _2^2 + (3\rho _3-\rho _1)\rho _2 +\rho _3\rho _1 }H, \\[-12pt] \nonumber \end{align}

(5.12) \begin{align} h_b^0&= \frac {(\rho _1 + \rho _2) \rho _3} {\rho _2^2 + (3\rho _3-\rho _1)\rho _2 +\rho _3\rho _1 }H. \\[9pt] \nonumber \end{align}

It is easy to check that

(5.13) \begin{align} \frac {\nu _2 (1-\nu _1 ) }{\nu _2^2 + (3-\nu _1)\nu _2 +\nu _1}\lt \frac {1}{3}, \end{align}

thus we always have $ 0\lt \mathcal{B}\lt ({1}/{3}) \mathcal{A} ^2$ . In Prodanov (Reference Prodanov2022a,Reference Prodanovb), three possible scenarios for the ratio $ {\mathcal{B}}/{\mathcal{A} ^2}$ can occur:

(5.14) \begin{align} 0\lt \mathcal{B}\lt \frac {2}{9} \mathcal{A} ^2, \\[-12pt] \nonumber \end{align}

(5.15) \begin{align} \frac {2}{9} \mathcal{A} ^2 \lt \mathcal{B}\lt \frac {1}{4} \mathcal{A} ^2 , \\[-12pt] \nonumber \end{align}

(5.16) \begin{align} \frac {1}{4} \mathcal{A} ^2 \lt \mathcal{B}\lt \frac {1}{3} \mathcal{A} ^2 . \\[9pt] \nonumber \end{align}

In each of the above scenarios, bounds detailing the locations of the three positive roots of $P(X)$ are available; cf. (5.22) below. Before passing to a discussion on these bounds, we examine the term on the left-hand side of (5.13), so we identify some conditions for the fulfilment of (5.14), (5.15) and (5.16). These conditions are given in terms of the ratio between the density of the top layer and suitable linear combinations of the densities of the other layers as follows:

(5.17) \begin{align} \text{(i)}\ 0\lt \mathcal{B}\lt \frac {2}{9} \mathcal{A} ^2 \quad \text{when} \quad \nu _1& \gt \frac {\nu _2(3-2\nu _2)}{7\nu _2+2}, \quad \text{i.e.} \quad \rho _1 \gt \rho _U=\frac {3\rho _3-2\rho _2}{7\rho _2+2\rho _3} \rho _2 ; \\[-12pt] \nonumber \end{align}

(5.18) \begin{align} \text{(ii)}\ 0 \lt \mathcal{B}\lt \frac {1}{4} \mathcal{A} ^2 \quad \text{when} \quad \nu _1 &\gt \frac {\nu _2(1-\nu _2)}{3\nu _2+1}, \quad \text{i.e.} \quad \rho _1 \gt \rho _L= \frac {\rho _3-\rho _2}{3\rho _2+\rho _3} \rho _2 . \\[9pt] \nonumber \end{align}

Of course, we always have $\rho _U\gt \rho _L$ . The above conditions in terms of the densities are only sufficient conditions. Unfortunately, there is no simple criterion (in terms of densities) that discriminates between (5.14), (5.15) and (5.16). For example, a necessary (but not sufficient) condition to have the situation in (5.16) is

(5.19) \begin{align} \rho _1 \lt \rho _L= \frac {\rho _3-\rho _2}{3\rho _2+\rho _3} \rho _2 . \end{align}

In most practical-relevant situations, the relative differences in densities are relatively small. It can be shown that this assumption leads to the situation in (5.14). More precisely, assuming that

(5.20) \begin{align} \frac {\rho _2-\rho _1}{\rho _2}\lt \frac {4}{9},\quad \frac {\rho _3-\rho _2}{\rho _3}\lt \frac {4}{9},\quad \frac {\rho _3-\rho _1}{\rho _3}\lt \frac {4}{9}, \end{align}

and using the Cauchy–Schwartz inequality, we infer from (5.6) that

(5.21) \begin{align} \mathcal{B}\lt \frac {4}{9} g^2(h_1 h+h_3(h-h_3))\lt \frac {4}{9} g^2\left(h_1 h+\frac {h^2}{4}\right)\lt \frac {4}{9} g^2\frac {(h+h_1)^2}{2}=\frac {2}{9}\mathcal{A}^2. \end{align}

Therefore, (5.20) is another sufficient condition that leads to the case (5.14). In practical applications, the internal waves occur between layers of one type of liquid (or very similar liquids) with slightly different properties (and densities), thus the condition (5.20) covers a wide range of practically possible situations.

According to Prodanov (Reference Prodanov2022a,Reference Prodanovb), when $ \mathcal{A}\lt 0$ for the three real roots of $P(X)$ , denoted by $0\lt X_3\leq X_2\leq X_1$ , in the case (5.14) the following bounds are valid:

(5.22) \begin{align} 0\lt -\mathcal{C} &\leq -\mathcal{C}_2, \nonumber\\ -\frac {\mathcal{C}}{\mathcal{B}} \leq X_3 &\leq \mu _2, \nonumber\\ \mu _2\leq X_2 &\leq \lambda _2, \nonumber\\ \lambda _1\leq X_1 &\leq \sigma _2, \end{align}

where the following notations are used:

(5.23) \begin{align} \mathcal{C}_1& :=-\frac {2\mathcal{A}^3}{27}+\frac {\mathcal{A}\mathcal{B}}{3}+\frac {2}{27}\sqrt {(\mathcal{A}^2-3\mathcal{B})^3}, \quad \mathcal{C}_2 :=-\frac {2\mathcal{A}^3}{27}+\frac {\mathcal{A}\mathcal{B}}{3}-\frac {2}{27}\sqrt {(\mathcal{A}^2-3\mathcal{B})^3}, \nonumber \\ \lambda _{1}&:=-\frac {\mathcal{A}}{2}+\sqrt {\frac {\mathcal{A}^2}{4}-\mathcal{B}},\quad \lambda _{2}:=-\frac {\mathcal{A}}{2} -\sqrt {\frac {\mathcal{A}^2}{4}-\mathcal{B}}, \nonumber \\ \mu _1&:=-\frac {\mathcal{A}}{3}+\frac {\sqrt {3}}{3}\sqrt {\frac {\mathcal{A}^2}{3}-\mathcal{B}}, \quad \mu _2:=-\frac {\mathcal{A}}{3}-\frac {\sqrt {3}}{3}\sqrt {\frac {\mathcal{A}^2}{3}-\mathcal{B}}, \end{align}

and

(5.24) \begin{align} \sigma _1:=-\frac {\mathcal{A}}{3}-2\frac {\sqrt {3}}{3}\sqrt {\frac {\mathcal{A}^2}{3}-\mathcal{B}}, \quad \sigma _2:=-\frac {\mathcal{A}}{3}+2\frac {\sqrt {3}}{3}\sqrt {\frac {\mathcal{A}^2}{3}-\mathcal{B}}. \end{align}

5.2. Approximate solutions for the wave speeds

We provide here some further insights resulting from (5.1). More precisely, we take into account several propagation regimes defined by the mutual relations existing between monomials in (5.1). We distinguish the following cases.

1. (i) If $a_2 c^2$ and $c_0$ are much smaller than $c^6$ and $a_4 c^4$ , then the approximate equation that we obtain has the roots $\pm \sqrt {g H}$ . However, from estimates (5.22) for the largest root, we observe that if $\mathcal{B}/\mathcal{A}^2\ll 1$ , then

   (5.25) \begin{align} \sigma _2 &= |\mathcal{A}|\left [1- \frac {\mathcal{B}}{\mathcal{A}^2}+ \mathcal{O}\left ( \left ( \frac {\mathcal{B}}{\mathcal{A}^2}\right )^2\right ) \right ]\!, \\[-12pt] \nonumber \end{align}
   (5.26) \begin{align} \lambda _1 &= |\mathcal{A}|\left [1- \frac {\mathcal{B}}{\mathcal{A}^2}+ \mathcal{O}\left ( \left ( \frac {\mathcal{B}}{\mathcal{A}^2}\right )^2\right ) \right ]\,\end{align}
   and the next correction to the approximate expression $X_1\approx gH$ is
   (5.27) \begin{align} X_1\approx |\mathcal{A}|- \frac {\mathcal{B}}{|\mathcal{A}|}= gH-\frac {g}{H} \left [\frac {\rho _2-\rho _1}{\rho _2}h_1 h_3+\frac {\rho _3-\rho _2}{\rho _3}h_3 h_b +\frac {\rho _3-\rho _1}{\rho _3}h_1 h_b \right ]\!. \end{align}
   Thus the propagation speed on the surface is approximately
   (5.28) \begin{align} c_{1,a}=\pm \sqrt {gH-\frac {g}{H} \left [\frac {\rho _2-\rho _1}{\rho _2}h_1 h_3+\frac {\rho _3-\rho _2}{\rho _3}h_3 h_b +\frac {\rho _3-\rho _1}{\rho _3}h_1 h_b \right ]}. \end{align}

2. (ii) The dominant terms are $a_4 c^4$ and $a_2 c^2$ . Indeed, from (5.22), for $X_2$ if $\mathcal{B}/\mathcal{A}^2\ll 1$ , then we have

   (5.29) \begin{align} \mu _2 &= \frac {\mathcal{B}}{2\,|\mathcal{A}|}+ |\mathcal{A}|\,\mathcal{O}\left ( \left ( \frac {\mathcal{B}}{\mathcal{A}^2}\right )^2\right )\!, \\[-12pt] \nonumber \end{align}
   (5.30) \begin{align} \lambda _2 &= \frac {\mathcal{B}}{|\mathcal{A}|}+ |\mathcal{A}|\,\mathcal{O}\left ( \left ( \frac {\mathcal{B}}{\mathcal{A}^2}\right )^2\right )\!, \\[9pt] \nonumber \end{align}
   therefore $X_2$ is of order $\mathcal{B}/|\mathcal{A}|\ll X_1$ . Moreover, again from (5.22), $(|\mathcal{C}|/\mathcal{B})\lt \mu _2\approx {\mathcal{B}}/{2|\mathcal{A}|}$ , thus $\mathcal{B}^2\gt 2\,|\mathcal{A}\mathcal{C}|$ . This justifies the further assumption $\mathcal{B}^2 \gg |\mathcal{A}\mathcal{C}|$ , under which it can be seen that when $c^2=X_2$ , the dominant terms in the equation are $a_4 c^4$ and $a_2 c^2$ . Then the solution of the approximate equation gives $ c_{2,a}^2=\mathcal{B}/|\mathcal{A}|$ or
   (5.31) \begin{align} c_{2,a}=\pm \sqrt {\frac {g}{H} \left [\frac {\rho _2-\rho _1}{\rho _2}h_1 h_3+\frac {\rho _3-\rho _2}{\rho _3}h_3 h_b +\frac {\rho _3-\rho _1}{\rho _3}h_1 h_b \right ] }. \end{align}

3. (iii) The dominant terms are $a_2 c^2$ and $a_0$ . From Viète formulas, we have

   (5.32) \begin{align} X_1X_2X_3=-\mathcal{C}=|\mathcal{C}|.\end{align}
   In the leading order, $X_1=|A|$ and $X_2=\mathcal{B}/|\mathcal{A}|$ , therefore $X_3= |\mathcal{C}/\mathcal{B}|$ , thus
   (5.33) \begin{align} c_{3,a}=\pm \sqrt {\frac {g(\rho _3 - \rho _2)(\rho _2-\rho _1)h_1 h_b h_3} {h_1h_b\rho _2(\rho _3-\rho _1) + h_b h_3\rho _2(\rho _3-\rho _2)+h_1h_3\rho _3(\rho _2-\rho _1)}}. \end{align}
   In these derivations, we made the assumptions $\mathcal{B}/\mathcal{A}^2\ll 1$ and $\mathcal{B}^2 \gg |\mathcal{A}\mathcal{C}|$ , which are compatible with the inequalities of (5.14), and of course with all kind of realistic physical data. Therefore, it is expected that the approximate formulas perform best in this case. Nevertheless, the results in table 1 show that (5.28), (5.31) and (5.33) perform quite well also in (5.15) and (5.16).

Table 1.

The approximate and exact wave speeds for some values of the parameters. Rows 1 and 2 correspond to internal waves, row 3 refers to equatorial internal waves, rows 4 and 5 pertain to three-layer lakes, and rows 6 and 7 refer to laboratory experiments. The dimensionalities of all quantities are given in the International System of Units (SI), $g=9.81$  m s⁻². All velocities are given in non-dimensional units by division by $\sqrt {gH}$ . The quantities $c_i$ are the positive roots of (5.1).

Remark 5.2. The two-layer limit case $\rho _3 =\rho _2$ , with depths $h_1$ and $h$ in (5.28) leads to

(5.34) \begin{align} c_{1,a}^*=\pm \sqrt {g(h_1+h)\left (1-\frac { h_1 h}{(h_1+h)^2} \frac {\rho _2-\rho _1}{\rho _2} \right )}. \end{align}

This appears also in e.g.Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005)and Kodaira et al. (Reference Kodaira, Waseda, Miyata and Choi2016).

Remark 5.3. The two-layer limit case $\rho _3 =\rho _2$ , with depths $h_1$ and $h$ in (5.31), leads to the well-known formula for the phase speed of an internal wave, which is related to the so-called Brunt–Väisälä frequency,

(5.35) \begin{align} c^*_{2,a}=\pm \sqrt {\frac {g (\rho _2 - \rho _1) h h_1 }{\rho _2(h+h_1)} }. \end{align}

As a matter of fact, the precise formula for two layers in the rigid lid approximation is formula (4.5) of Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005), but because in the ocean the two densities are very close ( $\rho _1 \approx \rho _2$ ), in many situations (5.35) is used instead (Massel Reference Massel2015 ; Kodaira et al. Reference Kodaira, Waseda, Miyata and Choi2016 ).

The propagation speed (5.31) is associated with the disturbance described by $\eta _2$ because the internal wave formed at this interface usually has a significant amplitude moving with the mentioned speed. Clearly, $\sqrt {gH}$ is the speed of the surface waves (over the total depth $H$ ) which due to the influence of the interfaces is modified to (5.28), while (5.33) is the speed of internal waves, observable at the interface $y=\eta _3$ . In general, waves with all possible speeds have their traces on the surface and on both interfaces. However, their amplitudes can be of quite different magnitudes.

Table 1 displays a comparison between the solutions of the approximate formulas above and the (positive) numerical solutions of (5.1). It is remarkable that the approximate formulas (5.28)–(5.33) give a very good approximation (for a range of physically feasible values of the parameters) to the actual solutions of (5.1). Rows 1–6 of table 1 correspond to (5.14), and the approximation formulas provide very good approximations. Row 7 illustrates (5.15). It does not seem that there are any real-life examples satisfying case (5.16). The cases (5.15) and (5.16) exhibit very large differences in the densities, and are not realistic in an oceanographic context. The approximation formulas work less well in the examples that exhibit large density contrasts, but nevertheless the approximate values are very good as a guide to the locations of the exact ones.

The first example in table 1 reflects typical parameter values for interfacial ocean internal waves. The next example (shown as no. 2 in table 1) concerns oceanographic and limnological systems, where stratified water columns are common due to variations in temperature, salinity and pressure. A three-layer structure with densities approximately 1000, 1020 and 1100 kg m⁻³, and respective thicknesses 10, 100 and 3000 m, can be justified by combining observations from coastal, open ocean and hypersaline environments; cf. https://oceanexplorer.noaa.gov/technology/ctd/ctd.html, https://www.bodc.ac.uk/data/bodc_database/ctd/, https://podaac.jpl.nasa.gov/dataset/OMG_L2_CTD.

The third example (no. 3 in table 1) is with parameter values, typical for the internal equatorial waves in the Pacific Ocean; see e.g. Constantin & Ivanov (Reference Constantin and Ivanov2019). The total depth of the Pacific Ocean is approximately 4000 m on average, and the upper layer is several hundred metres deep. We note that standard seawater density rarely exceeds 1050  $\rm kg\,m^{- 3}$ .

The next two examples (nos 4 and 5 in table 1) of a three-layer system concern temperate lakes, in which – particularly during the summer stratification period – the water column consistently organises into three distinct layers: the epilimnion, metalimnion and hypolimnion. The epilimnion, which typically extends from the surface to depths of 5–20 m, is a warm, well-mixed layer influenced by wind and solar heating. Beneath it lies the metalimnion, or thermocline, a transition zone several metres thick, where temperature and density change rapidly with depth, effectively inhibiting vertical mixing. The hypolimnion, the coldest and densest layer, lies below the thermocline and remains largely isolated from surface influences, often becoming oxygen-depleted over time because of biological activity and lack of mixing. This three-layer structure is not just a convenient simplification, but a physically necessary framework that accurately captures the thermal, chemical and biological dynamics of stratified lakes. Unlike oceanographic systems, where stratification is governed by both temperature and salinity over much greater depths, lake stratification is shallow, thermally driven and seasonally reversible, reflecting the unique density behaviour of freshwater and the limited depth of most inland water bodies (cf. Huttula Reference Huttula, Bengtsson, Herschy and Fairbridge2012).

To investigate flow behaviour in sharply stratified systems, we consider (in rows 6 and 7 of table 1) a three-layer fluid column with engineered density contrasts that far exceed those found in natural oceanographic environments (see National Institute of Standards and Technology, US Department of Commerce, NIST Chemistry WebBook, SRD 69). Unlike ocean stratification, which typically involves density differences of only 1–3 $\%$ due to temperature and salinity gradients, this system spans nearly an order of magnitude in density, using immiscible or semi-miscible fluids to maintain stable, well-defined interfaces.

The top layer typically consists of light hydrocarbons such as pentane or hexane. These fluids have densities in the range 600–700 $\rm kg\,m^{- 3}$ and form a buoyant, sharply defined upper boundary. In experimental set-ups, this layer is often kept relatively thin, with height approximately 1–2 cm, to minimise mixing and maintain interface clarity.

The middle layer is composed of dense aqueous solutions, such as $\rm ZnCl_2$ brines or glycerol–water mixtures, with densities in the range 1000–2200 $\rm kg\,m^{- 3}$ . This layer serves as a transitional buffer and is typically the thickest, with heights 3–5 cm, allowing for stable stratification and tunable flow dynamics. Its density can be adjusted to fine-tune the separation between the top and bottom layers.

The bottom layer features ultra-dense fluids such as perfluorodecalin, or liquid metals such as gallium alloys, which can reach densities 6400 $\rm kg\,m^{- 3}$ . This layer is usually 1–3 cm thick and provides a heavy, immiscible foundation that anchors the system and prevents upward mixing.

7.1. The DN operators and evolution equations

We use the techniques developed in the previous sections and repeat the same steps. The Hamiltonian formulation can be established in this case as well; however, we are going to employ only techniques related to the DN operators. For the middle and lower layers, we have the DN operators (3.22) and (3.23). The DN operator for the upper layer is as in (2.12) of Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005) and is defined as

(7.1) \begin{align} \left ( \frac {\partial \varphi _1}{\partial (-\boldsymbol {n}_2)} \right )_{s_2}\sqrt {1+\eta _{2,x}^2}=-{\boldsymbol {n}_2}\boldsymbol{\cdot }(\boldsymbol{\nabla }\varphi _1)_{s_2} \equiv - G_1(\eta _2) (\varphi _1)_{s_2}. \end{align}

More details are provided in (A.10) of Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005), e.g.

(7.2) \begin{align} G_1(\eta _2)=\sum _{l=0}^{\infty }G_1^{(l)}(\eta _2), \\[-28pt] \nonumber \end{align}

(7.3) \begin{align} G_1^{(0)}= D \tanh (h_1 D), \quad G_1^{(1)}=- D \eta _2(x)\, D + G_1^{(0)}\, \eta _2(x)\, G_1^{(0)} +\cdots . \\[2pt] \nonumber \end{align}

Then obviously

(7.4) \begin{align} \eta _{2,t}=- G_1(\eta _2) (\varphi _1)_{s_2}=\varGamma _{11}(\varphi _2)_{s_2}+\varGamma _{12}(\varphi _2)_{s_3}, \end{align}

and as in (3.27), we obtain

(7.5) \begin{align} \eta _{3,t}=G(\eta _3)(\varphi _3)_{s_3}=-[\varGamma _{21}(\varphi _2)_{s_2}+\varGamma _{22}(\varphi _2)_{s_3}]. \end{align}

Equations (7.4) and (7.5) together with

(7.6) \begin{align} \xi _2:=\rho _2(\varphi _2)_{s_2}-\rho _1(\varphi _1)_{s_2}, \quad \xi _3:=\rho _3(\varphi _3)_{s_3}-\rho _2(\varphi _2)_{s_3} \end{align}

can be solved for $(\varphi _1)_{s_2}$ , $(\varphi _2)_{s_2}$ , $(\varphi _2)_{s_3}$ and $(\varphi _3)_{s_3}$ in terms of $\xi _2,\xi_3$ and the corresponding DN operators (that depend on $\eta _2$ and $\eta _3$ ). We obtain the following expressions (different from those in (3.35) and (3.36)):

(7.7) \begin{align} (\varphi _2)_{s_2}&=B^{-1}\big [ \rho _1 G \xi _3 + (\rho _3 \varGamma _{22}+\rho _2 G) \varGamma _{12}^{-1} G_1 \xi _2 \big ]=:\varPhi ^2_2(\eta _i, \xi _i) , \\[-12pt] \nonumber \end{align}

(7.8) \begin{align} (\varphi _2)_{s_3}&=-(B^*)^{-1}\big [ \rho _3 G_1 \xi _2 + (\rho _1 \varGamma _{11}+\rho _2 G_1) \varGamma _{21}^{-1} G \xi _3 \big ]=:\varPhi ^2_3(\eta _i, \xi _i), \\[9pt] \nonumber \end{align}

where

(7.9) \begin{align} B= (\rho _3 \varGamma _{22}+\rho _2 G) \varGamma _{12}^{-1}(\rho _2 G_1+\rho _1 \varGamma _{11})-\rho _1\rho _3 \varGamma _{21}. \end{align}

This allows to express (7.4) and (7.5) in the form

(7.10) \begin{align} \eta _{2,t}&= A_{22}\xi _2+A_{23} \xi _3, \\[-12pt] \nonumber \end{align}

(7.11) \begin{align} \eta _{3,t}&= A_{32}\xi _2+A_{33} \xi _3 , \\[9pt] \nonumber \end{align}

where

(7.12) \begin{align} A_{22}&=-\rho _3 \varGamma _{12} (B^*)^{-1} G_1+ \varGamma _{11} B^{-1}(\rho _3 \varGamma _{22}+\rho _2 G) \varGamma _{12}^{-1} G_1, \\[-12pt] \nonumber \end{align}

(7.13) \begin{align} A_{23}&=\rho _1 \varGamma _{11} B^{-1} G -\varGamma _{12} (B^*)^{-1}(\rho _1 \varGamma _{11}+\rho _2 G_1) \varGamma _{21}^{-1} G, \\[-12pt] \nonumber \end{align}

(7.14) \begin{align} A_{32}&= \rho _3 \varGamma _{22} (B^*)^{-1} G_1 -\varGamma _{21} B^{-1}(\rho _3 \varGamma _{22}+\rho _2 G) \varGamma _{12}^{-1} G_1, \\[-12pt] \nonumber \end{align}

(7.15) \begin{align} A_{33}&=-\rho _1 \varGamma _{21} B^{-1} G +\varGamma _{22} (B^*)^{-1}(\rho _1 \varGamma _{11}+\rho _2 G_1) \varGamma _{21}^{-1} G \\[9pt] \nonumber \end{align}

are operators that depend only on $\eta _2$ and $\eta _3$ . These operators, of course, are not the same as those in (3.38). Thus defining

(7.16) \begin{align} F_2(\eta _i, \xi _i):= A_{22}\xi _2+A_{23} \xi _3, \quad F_3(\eta _i, \xi _i):=A_{32}\xi _2+A_{33} \xi _3 , \quad i=2,3, \end{align}

we have the nonlinear dynamics for $\eta _2$ and $\eta _3$ from (7.4) in the form

(7.17) \begin{align} \eta _{2,t}= F_2(\eta _i, \xi _i), \quad \eta _{3,t}=F_3(\eta _i, \xi _i), \end{align}

and also

(7.18) \begin{align} (\varphi _1)_{s_2}&=-G_1^{-1}\eta _{2,t}=-G_1(\eta _2)^{-1}F_2(\eta _i, \xi _i)=:\varPhi ^1_2(\eta _i, \xi _i), \\[-12pt] \nonumber \end{align}

(7.19) \begin{align} (\varphi _3)_{s_3}&=G^{-1}\eta _{3,t}=G(\eta _3)^{-1}F_3(\eta _i, \xi _i)=:\varPhi ^3(\eta _i, \xi _i). \\[9pt] \nonumber \end{align}

The time evolution of $\xi _2$ and $\xi _3$ is given by (3.53) and (3.55) as

(7.20) \begin{align} \xi _{2,t}&+\frac {\rho _2 (\varPhi ^2_{2,x})^2-\rho _1 (\varPhi _{2,x}^1)^2-2F_2\xi _{2,x}\eta _{2,x}+(\rho _1-\rho _2)F_2^2 }{2(1+\eta _{2,x}^2)}+(\rho _2-\rho _1)g \eta _2 =0, \\[-12pt] \nonumber \end{align}

(7.21) \begin{align} \xi _{3,t}&+\frac {\rho _3 (\varPhi ^3_{x})^2-\rho _2 (\varPhi _{3,x}^2)^2-2F_3\xi _{3,x}\eta _{3,x}+(\rho _2-\rho _3)F_3^2 }{2(1+\eta _{3,x}^2)}+(\rho _3-\rho _2)g \eta _3 =0, \\[9pt] \nonumber \end{align}

with the new definitions of $F_2,F_3$ in (7.16), and the quantities in (7.7), (7.8), (7.18) and (7.19).

Thus we have the full nonlinear system of the time evolution in terms of DN operators given by (7.17), (7.20) and (7.21). We will illustrate how this system produces model equations in different propagation regimes.

7.2. Leading-order linear approximation and dispersion relation

In the leading-order approximation, we are using the DN operators from (7.3), (4.8) and (4.7) to obtain

(7.22) \begin{align} A_{22}^0(k)&= \frac {k[\rho _2 \cosh (kh_3)\sinh (kh_b) + \rho _3 \sinh (kh_3)\cosh (kh_b)]\sinh (k h_1)}{ \boldsymbol {D}(k) }, \\[-12pt] \nonumber \end{align}

(7.23) \begin{align} A_{23}^0(k)&=A_{32}^0(k)= \frac {k \rho _2 \sinh (k h_b)\sinh (k h_1) }{ \boldsymbol {D}(k) }, \\[-12pt] \nonumber \end{align}

(7.24) \begin{align} A_{33}^0(k)&= \frac {k[\rho _1 \cosh (k h_1) \sinh (kh_3) + \rho _2 \cosh (kh_3)\sinh (kh_1)]\sinh (k h_b)}{ \boldsymbol {D}(k) }, \\[9pt] \nonumber \end{align}

where

(7.25) \begin{align} \begin{split} \boldsymbol {D}(k)={}& \rho _2\Big (\rho _3 \cosh (k h_3) \cosh (k h_b)+\rho _2 \sinh (kh_3)\sinh (kh_b)\Big )\sinh (kh_1) \\ & {}+ \rho _1\Big (\rho _2 \cosh (k h_3) \sinh (k h_b)+\rho _3 \sinh (kh_3)\cosh (kh_b)\Big )\cosh (kh_1). \end{split} \end{align}

Note that for $k\ne 0$ , the sign of $\boldsymbol {D}(k)$ coincides with the sign of $k$ , thus $\boldsymbol {D}(k)\ne 0$ . The equations for $\xi _i$ are as in (4.13), and the linearised equations in the leading-order approximation therefore are

(7.26) \begin{align} \eta _{2,t}&=A^0_{22}(D)\,\xi _2+A^0_{23}(D)\, \xi _3, \\[-12pt] \nonumber \end{align}

(7.27) \begin{align} \eta _{3,t}&=A_{32}^{0}(D)\,\xi _2+A^0_{33}(D)\,\xi _3, \\[-12pt] \nonumber \end{align}

(7.28) \begin{align} \xi _{2,t}&=-g(\rho _2-\rho _1)\eta _2, \\[-12pt] \nonumber \end{align}

(7.29) \begin{align} \xi _{3,t}&=-g(\rho _3-\rho _2)\eta _3. \\[9pt] \nonumber \end{align}

The analogue of the matrix $\mathcal{M}$ in (4.14) is now

(7.30) \begin{align} \mathcal{M}(k)=\begin{pmatrix} 0 & 0 & A^0_{22}(k) & A^0_{23}(k)\\[5pt] 0 & 0 & A^0_{23}(k) & A^0_{33}(k)\\[5pt] g(\rho _1-\rho _2) & 0 & 0 &0 \\[5pt] 0 & g(\rho _2-\rho _3) & 0 & 0 \end{pmatrix}\!. \end{align}

The characteristic equation $\textrm {det} (\mathcal{M}(k)+ikcI_4)=0$ leads to the following bi-quadratic equation for the propagation speeds $c$ :

(7.31) \begin{align} &c^4 + a_2 c^2 +a_0 =0, \\[-12pt] \nonumber \end{align}

(7.32) \begin{align} a_2&=-\frac {g(A^0_{22}(\rho _2 - \rho _1) + A^0_{33}(\rho _3-\rho _2)) }{k^2}=-\frac {g\, \boldsymbol {D}_1(k)}{k\,\boldsymbol {D}(k)}, \\[-12pt] \nonumber \end{align}

(7.33) \begin{align} a_0&= \frac {g^2(\rho _3 - \rho _2)(\rho _2 - \rho _1)(A^0_{22}A^0_{33} - (A^0_{23})^2)}{k^4} \nonumber \\ & = \frac {g^2(\rho _3 - \rho _2)(\rho _2 - \rho _1)\sinh (kh_1)\sinh (kh_3)\sinh (k h_b)}{k^2\, \boldsymbol {D}(k)}, \\[9pt] \nonumber \end{align}

where

(7.34) \begin{align} \boldsymbol {D}_1(k)={}& [\rho _2(\rho _3-\rho _1)\cosh (k h_3) \sinh (k h_b) \nonumber \\&{}+ \rho _3(\rho _2-\rho _1) \sinh (kh_3)\cosh (kh_b)]\sinh (k h_1) \nonumber \\ &{}+\rho _1(\rho _3-\rho _2)\cosh (kh_1) \sinh (kh_3)\sinh (k h_b). \end{align}

7.3. Short-wave limit

The short-wave limit for the solutions of (7.31) can be obtained in a similar way by the asymptotics $\sinh (k h_i) \sim \pm (1/2)\exp (kh_i)$ and $ \cosh (k h_i) \sim (1/2)\exp (k h_i)$ for $k\to \pm \infty $ . Thus

(7.35) \begin{align} a_2 & \sim -\frac {g}{|k|} \left (\frac {\rho _2 - \rho _1}{\rho _2 + \rho _1}+ \frac {\rho _3-\rho _2}{\rho _3+\rho _2} \right )\!, \\[-12pt] \nonumber \end{align}

(7.36) \begin{align} a_0 & \sim \frac {g^2}{k^2} \left (\frac {\rho _2 - \rho _1}{\rho _2 + \rho _1}\right ) \left (\frac {\rho _3-\rho _2}{\rho _3+\rho _2}\right ) \!. \\[9pt] \nonumber \end{align}

The roots of (7.31) are real as well, and coincide with the two already known values in (4.29):

(7.37) \begin{align} c_{1}^2 \sim \frac {g}{|k|} \left (\frac {\rho _2 - \rho _1}{\rho _2 + \rho _1}\right )\! , \quad c_{2}^2 \sim \frac {g}{|k|} \left (\frac {\rho _3 - \rho _2}{\rho _3 + \rho _2}\right )\!. \end{align}

7.4. Long-wave limit

Passing to the long-wave limit $k\to 0$ gives

(7.38) \begin{align} a_2 & \to -g\frac {\rho _2(\rho _3-\rho _1) h_1 h_b + \rho _3(\rho _2-\rho _1) h_1 h_3 + \rho _1 (\rho _3-\rho _2) h_3 h_b}{\rho _2 \rho _3 h_1 + \rho _1 \rho _2 h_b + \rho _1 \rho _3 h_3}, \\[-12pt] \nonumber \end{align}

(7.39) \begin{align} a_0 & \to \frac {g^2(\rho _3 - \rho _2)(\rho _2 - \rho _1) h_1 h_3 h_b}{\rho _2 \rho _3 h_1 + \rho _1 \rho _2 h_b + \rho _1 \rho _3 h_3}, \\[9pt] \nonumber \end{align}

therefore we obtain the equation

(7.40) \begin{align} &(\rho _2 \rho _3 h_1 + \rho _1 \rho _2 h_b + \rho _1 \rho _3 h_3) c^4 \nonumber \\ &{}-g[\rho _2(\rho _3-\rho _1) h_1 h_b + \rho _3(\rho _2-\rho _1) h_1 h_3 + \rho _1 (\rho _3-\rho _2) h_3 h_b] c^2 \nonumber \\ & {}+g^2(\rho _3 - \rho _2)(\rho _2 - \rho _1) h_1 h_3 h_b=0, \end{align}

which appeared also in de Melo Viríssimo & Milewski (Reference de Melo Viríssimo and Milewski2019) and Ratliff (Reference Ratliff2022), where it is obtained by other methods. We will prove in the following that (7.40) has four real roots for all values of the densities and all values of the (mean) depths $h_1,h_3,h_b$ . To begin, we denote $c^2=Z$ in (7.40) and obtain the second-degree equation in $Z$

(7.41) \begin{align} (\rho _2 \rho _3 h_1 & + \rho _1 \rho _2 h_b + \rho _1 \rho _3 h_3) Z^2 \nonumber \\ & -g[\rho _2(\rho _3-\rho _1) h_1 h_b + \rho _3(\rho _2-\rho _1) h_1 h_3 + \rho _1 (\rho _3-\rho _2) h_3 h_b] Z \nonumber \\ & +g^2(\rho _3 - \rho _2)(\rho _2 - \rho _1) h_1 h_3 h_b=0. \end{align}

Taking into account that $\rho _3\gt \rho _2\gt \rho _1$ , and applying the Cauchy-Schwarz inequality, we have first

(7.42) \begin{align} &\rho _2(\rho _3-\rho _1) h_1 h_b + \rho _3(\rho _2-\rho _1) h_1 h_3 + \rho _1 (\rho _3-\rho _2) h_3 h_b \nonumber \\ &\geq 2h_1\sqrt {h_3 h_b}\sqrt {\rho _2\rho _3(\rho _3-\rho _1)(\rho _2-\rho _1)}\nonumber \\ & +2h_b\sqrt {h_1 h_3}\sqrt {\rho _1\rho _2(\rho _3-\rho _1)(\rho _3-\rho _2)}\nonumber \\ & +2h_3\sqrt {h_1 h_b}\sqrt {\rho _1\rho _3 (\rho _2-\rho _1)(\rho _3-\rho _2)}\nonumber \\ &\gt 2\sqrt {(\rho _3-\rho _2)(\rho _2-\rho _1)}\Big [\sqrt {\rho _2\rho _3}\,h_1\sqrt {h_3 h_b}+\sqrt {\rho _1\rho _2}\,h_b\sqrt {h_1 h_3}+\sqrt {\rho _1\rho _3}\,h_3\sqrt {h_1 h_b} \Big ]\nonumber \\ &=2\sqrt {(\rho _3-\rho _2)(\rho _2-\rho _1)}\sqrt {h_1 h_3 h_b}\Big [\sqrt {\rho _2\rho _3 h_1}+\sqrt {\rho _1\rho _2 h_b}+\sqrt {\rho _1\rho _3 h_3}\Big ]\nonumber \\ &\geq 2\sqrt {(\rho _3-\rho _2)(\rho _2-\rho _1)}\sqrt {h_1 h_3 h_b} \sqrt {\rho _2\rho _3 h_1+\rho _1\rho _2 h_b+\rho _1\rho _3 h_3}, \end{align}

where for the last inequality we have used that $\rho _i\ (i=1,2,3),\ h_1,h_3,h_b$ are positive. By squaring in the obtained inequality, we see that the discriminant of (7.41) is strictly positive, i.e. (7.41) has two distinct real roots, which, by the Viète formulas, are positive. From the latter finding we infer that the fourth-order equation (7.40) has four real roots: two positive and two negative.

We obtain an approximation for the smaller root by neglecting $c^4$ in (7.40),

(7.43) \begin{align} c_{-}:=\pm \sqrt {\frac {g(\rho _3 - \rho _2)(\rho _2 - \rho _1) h_1 h_3 h_b}{\rho _2(\rho _3-\rho _1) h_1 h_b + \rho _3(\rho _2-\rho _1) h_1 h_3 + \rho _1 (\rho _3-\rho _2) h_3 h_b }}, \end{align}

which in turn is approximately equal to (5.33) when $\rho _1 \approx \rho _2$ in the expression $ \rho _1 (\rho _3-\rho _2) h_3 h_b \approx \rho _2 (\rho _3-\rho _2) h_3 h_b$ .

By neglecting the last term in (7.40), we obtain an approximation for the larger root,

(7.44) \begin{align} c_{+}:=\pm \sqrt {\frac {g[\rho _2(\rho _3-\rho _1) h_1 h_b + \rho _3(\rho _2-\rho _1) h_1 h_3 + \rho _1 (\rho _3-\rho _2) h_3 h_b]}{\rho _2 \rho _3 h_1 + \rho _1 \rho _2 h_b + \rho _1 \rho _3 h_3}}. \end{align}

Under the further assumption that $\rho _1\approx \rho _2 \approx \rho _3 = \rho$ in the denominator and in the ratios $\rho _i/\rho _{\!j}$ , we have

(7.45) \begin{align} c_{+} &\approx \pm \sqrt {\frac {g[\rho _2(\rho _3-\rho _1) h_1 h_b + \rho _3(\rho _2-\rho _1) h_1 h_3 + \rho _1 (\rho _3-\rho _2) h_3 h_b]}{\rho ^2( h_1 + h_b + h_3)}} \nonumber \\ & \approx \pm \sqrt {\frac {g}{h_1+h} \left (\frac {(\rho _2-\rho _1) h_1 h_3}{\rho _2}+ \frac {(\rho _3-\rho _2) h_3 h_b}{\rho _3}+\frac {(\rho _3-\rho _1) h_1 h_b}{\rho _3} \right )} , \end{align}

which is (5.33).

Therefore the approximation formulas for both models agree under the further assumption of small relative differences between densities, $|\rho _i - \rho _{\!j}|/\rho _n \ll 1$ .

Finally, we point out that each propagation speed has its trace at every surface. For the largest propagation speed $c_+$ (mode 1), the ratio of the wave amplitudes $\eta _2,\eta_3$ of the upper and lower interfaces in the linear regime is $ ({\eta _2}/{\eta _3} )_{c=c_+} \gt 0$ . In this case, both waves have the same polarity. For the slow-propagating mode $c_-$ (mode 2), $ ({\eta _2}/{\eta _3} )_{c=c_-} \lt 0,$ and in this case the waves have the opposite polarity. See e.g. the details in de Melo Viríssimo & Milewski (Reference de Melo Viríssimo and Milewski2019).

7.5. Boussinesq regime for the rigid lid system

The next example of evolution equations arising from (7.17), (7.20) and (7.21) is the Boussinesq regime. As a result, we obtain a coupled system of nonlinear evolution equations. We recall the expansions of the DN operators in the case under consideration (Craig et al. Reference Craig, Guyenne, Nicholls and Sulem2005; Ivanov et al. Reference Ivanov, Martin and Todorov2022):

(7.46) \begin{align} G_1(\eta _2)={}& D\tanh (h_1D)-D\eta _2 D+D\tanh (h_1D)\, \eta _2 D\tanh (h_1D)+\cdots , \\[-12pt] \nonumber \end{align}

(7.47) \begin{align} G(\eta _3)={}& D\tanh (h_bD)+D\eta _3 D-D\tanh (h_b D)\, \eta _3 D\tanh (h_bD)+\cdots ,\\[-12pt] \nonumber \end{align}

(7.48) \begin{align} \varGamma _{11}(\eta _2, \eta _3)={}& D \coth (h_3D) + D\,{\textrm {csch}}(h_3D)\,\eta _3 D\,{\textrm {csch}} (h_3D) \nonumber \\ & {}+ D\eta _2D-D\coth (h_3D)\,\eta _2D\coth (h_3D)+\cdots ,\\[-12pt] \nonumber \end{align}

(7.49) \begin{align} \varGamma _{12}(\eta _2,\eta _3)={}& -D\,{\textrm {csch}}(h_3D)-D\,{\textrm {csch}}(h_3D)\,\eta _3 D\coth (h_3D) \nonumber \\ & {}+D \coth (h_3D)\,\eta _2 D\, {\textrm {csch}}(h_3D)+\cdots, \\[-12pt] \nonumber \end{align}

(7.50) \begin{align} \varGamma _{21}(\eta _2,\eta _3)={}& -D\,{\textrm {csch}}(h_3D) -D\coth (h_3D)\,\eta _3 D\, {\textrm {csch}}(h_3D) \nonumber \\ & {}+D\,{\textrm {csch}}(h_3D)\, \eta _2 D \coth (h_3D)+\cdots ,\\[-12pt] \nonumber \end{align}

(7.51) \begin{align} \varGamma _{22}(\eta _2, \eta _3)={}& D \coth (h_3D) + D\coth (h_3D)\,\eta _3 D\coth (h_3D) - D \eta _3 D \nonumber \\ & {}-D\,{\textrm {csch}}(h_3D)\, \eta _2 D\, {\textrm {csch}}(h_3D)+\cdots . \\[9pt] \nonumber\end{align}

We use these expansions in the Boussinesq regime, i.e. assuming that the small-amplitude parameters for each combination $\varepsilon _{\textit{ij}} = {|\eta _{\!j}|_{{max}}}/{h_i}$ are of the same order of magnitude, say $\varepsilon $ , for all $i,j$ , and also that the long-wave parameters $\delta _{\!j}={2 \pi h_{\!j}}/{\lambda }=k h_{\!j}$ are of the same order $\delta $ , and moreover $\varepsilon \ll 1$ and $ \delta ^2 \ll 1$ are of the same order as well. Furthermore, we recall the fact that the eigenvalue of $D$ is the wavenumber $k$ , so the action of $D$ will generate a quantity of order $\delta $ . On the other hand, the multiplication by $\eta _{\!j}$ will generate a quantity of order $\varepsilon$ . Keeping this in mind, the relevant expansion of (7.12)–(7.15) in terms of the DN operators (7.46)–(7.51) is

(7.52) \begin{align} A_{22}=A_{220}D^2 + A_{221}D^4+ B_{222} D\eta _2D+B_{223} D\eta _3 D +\cdots , \\[-12pt] \nonumber \end{align}

(7.53) \begin{align} A_{23}=A_{32}=A_{230}D^2 + A_{231}D^4+ B_{232} D\eta _2D+B_{233} D\eta _3 D +\cdots , \\[-12pt] \nonumber \end{align}

(7.54) \begin{align} A_{33}=A_{330}D^2 + A_{331}D^4+ B_{332} D\eta _2D+B_{333} D\eta _3 D +\cdots , \\[9pt] \nonumber \end{align}

where the expansion coefficients are

(7.55) \begin{align} & A_{220} = \frac {h_1(h_3\rho _3 + h_b\rho _2) }{ h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2}, \\[-9pt] \nonumber \end{align}

(7.56) \begin{align} & A_{221}=- \dfrac {\begin{array}{c} h_1^2 \big(h_1h_3^2\rho _1\rho _3^2 + 2h_1h_3h_b\rho _1\rho _2\rho _3 + h_1h_b^2\rho _1\rho _2^2 + h_3^3\rho _2\rho _3^2 + 3h_3^2h_b\rho _2^2\rho _3 \big.\\\big. + 3h_3h_b^2\rho _2^3 + h_b^3\rho _2^2\rho _3\big)\end{array}}{3(h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2)^2}, \\[-9pt] \nonumber \end{align}

(7.57) \begin{align} & B_{222}= \frac {-(h_3\rho _3 + h_b\rho _2)^2\rho _1 + h_1^2\rho _2\rho _3^2}{(h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2)^2}, \\[-9pt] \nonumber \end{align}

(7.58) \begin{align} & B_{223}= \frac {(\rho _2-\rho _3)h_1^2\rho _2\rho _3 }{(h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2)^2}, \\[-9pt] \nonumber \end{align}

(7.59) \begin{align} & A_{230}=\frac { h_1 h_b \rho _2}{h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2} , \\[-9pt] \nonumber \end{align}

(7.60) \begin{align} & A_{231} = - \dfrac { \begin{array}{c} h_1 h_b \rho _2\big[6h_1h_bh_3 \rho _2^2+ h_b(2h_1^2 + 3 h_3^2)\rho _1\rho _2 + (3 h_3^2 + 2h_b^2)h_1\rho _2\rho _3 \big.\\ + \big. h_3(2h_1^2 + h_3^2 +2 h_b^2)\rho _1\rho _3 \big]\end{array}}{6(h_1\rho _2\rho _3 + h_3\rho _1 \rho _3 + h_b\rho _1\rho _2)^2} , \\[-9pt] \nonumber \end{align}

(7.61) \begin{align} & B_{232}= -\frac {h_b \rho _1 \rho _2 [h_b \rho _2 + (h_1 + h_3)\rho _3] }{( h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2 )^2 } , \\[-9pt] \nonumber \end{align}

(7.62) \begin{align} & B_{233}=\frac {h_1 \rho _3\rho _2 [(h_3 + h_b)\rho _1 + h_1 \rho _2 ]}{( h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2)^2}, \\[9pt] \nonumber\end{align}

(7.63) \begin{align} & A_{330}=\frac {h_b(h_1\rho _2 + h_3\rho _1)}{h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2} \\[-9pt] \nonumber \end{align}

(7.64) \begin{align} & A_{331}\nonumber \\ &=-\frac {\begin{array}{c} h_b^2 \big(h_1^3\rho _1\rho _2^2 + 3h_1^2h_3\rho _2^3 + h_1^2h_b\rho _2^2\rho _3 + 3h_1h_3^2\rho _1\rho _2^2 + 2h_1h_3h_b\rho _1\rho _2\rho _3 + h_3^3\rho _1^2\rho _2 \big.\\ \big. + h_3^2h_b\rho _1^2\rho _3 \big)\end{array}}{3(h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2)^2}, \\[-9pt] \nonumber \end{align}

(7.65) \begin{align} & B_{332}=\frac {h_b^2(\rho _1 -\rho _2)\rho _1\rho _2}{(h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2)^2}, \\[-9pt] \nonumber \end{align}

(7.66) \begin{align} & B_{333}=\frac {(h_1\rho _2 + h_3\rho _1)^2\rho _3 - h_b^2\rho _1^2\rho _2}{(h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2)^2}. \\[9pt] \nonumber \end{align}

Since in the Boussinesq regime the equations are maximum quadratic in the field variables, (7.20) and (7.21) suggest that we need only the leading order of the quantities $\varPhi ^i_{\!j}$ and $F_{\!j}$ :

(7.67) \begin{align} \varPhi ^2_2=&\frac {h_1\rho _3\xi _2 - h_b\rho _1\xi _3}{h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2}, \\[-12pt] \nonumber \end{align}

(7.68) \begin{align} \varPhi ^1_2=&-\frac {h_b \rho _2 (\xi _2 + \xi _3) + h_3\rho _3\xi _2}{h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2}, \\[-12pt] \nonumber \end{align}

(7.69) \begin{align} \varPhi ^3=&\frac {h_1 (\xi _2 +\xi _3)\rho _2 + h_3 \xi _3\rho _1}{h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2} , \\[-12pt] \nonumber \end{align}

(7.70) \begin{align} \varPhi ^2_3=&\frac {h_1\rho _3\xi _2 - h_b\rho _1\xi _3}{h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2}. \\[9pt] \nonumber \end{align}

Introducing the variables $\mathfrak{u}_{\!j}=(\xi _{\!j})_x$ , we first notice that $F_{\!j}$ in leading order are linear combinations of $\mathfrak{u}_{j,x}$ , while $\varPhi ^i_{j,x}$ are linear combinations of $\mathfrak{u}_{j}$ only. However, we notice that $\mathcal{O}(h_{\!j} \mathfrak{u}_{j,x}) = \mathcal{O}(kh_{\!j} \mathfrak{u}_{j}) = \delta\, \mathcal{O}(\mathfrak{u}_{j})\ll \mathcal{O}(\mathfrak{u}_{j})$ . Thus the terms with $F_{\!j}$ in (7.20) and (7.21) can be neglected in comparison to the other terms with $\varPhi ^i_{j,x}$ (as in the water wave theory of a single layer). Then we obtain the system

(7.71) \begin{align} \mathfrak{u}_{2,t} & + \frac {[\rho _2(h_1\rho _3\mathfrak{u}_2 - h_b\rho _1\mathfrak{u}_3)^2-\rho _1(h_b \rho _2 (\mathfrak{u}_2 + \mathfrak{u}_3) + h_3\rho _3\mathfrak{u}_2)^2]_x}{2(h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2)^2} +(\rho _2-\rho _1)g\eta _{2,x}=0 , \nonumber \\\mathfrak{u}_{3,t} &+ \frac {[\rho _3(h_1 \rho _2 (\mathfrak{u}_2 +\mathfrak{u}_3) + h_3 \rho _1 \mathfrak{u}_3 )^2-\rho _2(h_1\rho _3\mathfrak{u}_2 - h_b\rho _1\mathfrak{u}_3)^2]_x}{2(h_1\rho _2\rho _3+ h_3\rho _1 \rho _3 + h_b\rho _1\rho _2)^2} +(\rho _3-\rho _2)g\eta _{3,x}=0, \nonumber\\\eta _{2,t} & = -A_{220}\mathfrak{u}_{2,x} + A_{221}\mathfrak{u}_{2,\textit{xxx}}- B_{222} [\eta _2 \mathfrak{u}_{2}]_x -B_{223} [\eta _3 \mathfrak{u}_{2} ]_x \nonumber\\&{}-A_{230} \mathfrak{u}_{3,x} + A_{231}\mathfrak{u}_{3,\textit{xxx}}- B_{232} [\eta _2 \mathfrak{u}_{3}]_x -B_{233} [\eta _3 \mathfrak{u}_{3} ]_x ,\nonumber\\\eta _{3,t} & = -A_{230}\mathfrak{u}_{2,x} + A_{231} \mathfrak{u}_{2,\textit{xxx}} - B_{232} [\eta _2 \mathfrak{u}_{2}]_x -B_{233} [\eta _3 \mathfrak{u}_{2} ]_x \nonumber\\& {}-A_{330}\mathfrak{u}_{3,x} + A_{331} \mathfrak{u}_{3,\textit{xxx}} - B_{332} [\eta _2 \mathfrak{u}_{3}]_x -B_{333} [\eta _3 \mathfrak{u}_{3} ]_x . \end{align}

Letting $\rho _1 \to 0$ , we obtain the two-layer domain with surface $y=\eta _2(x,t)$ and interface $y=\eta _3(x,t)$ , whereby the coefficients $A_{ijk}, B_{ijk}\ (i, j\in \{1,2\},\ k\in \{0,1,2,3\})$ are given by

(7.72) \begin{align} A_{220}& =\frac {h_3\rho _3 + h_b\rho _2}{ \rho _2\rho _3}, \quad A_{221} = - \frac { h_3^3 \rho _3^2 + 3h_3^2h_b\rho _2 \rho _3 + 3h_3h_b^2\rho _2^2 + h_b^3\rho _2 \rho _3}{3\rho _2 \rho _3 ^2}, \nonumber \\ B_{222}& = \frac {1}{\rho _2 },\quad B_{223} = \frac {\rho _2-\rho _3 }{\rho _2\rho _3 }, \nonumber \\ A_{230}& = \frac { h_b }{\rho _3 } ,\quad A_{231} = - \frac { h_b \big[6h_bh_3 \rho _2+ \big(3 h_3^2 + 2h_b^2 \big) \rho _3 \big]}{6 \rho _3^2} , \quad B_{232} = 0 ,\quad B_{233} = \frac { 1 }{\rho _3 }, \nonumber \\ A_{330}& = \frac {h_b }{\rho _3}, \quad A_{331} = -\frac {h_b^2 ( 3 h_3\rho _2 + h_b \rho _3 )}{3 \rho _3 ^2},\quad B_{332} = 0,\quad B_{333} =\frac {1}{\rho _3 }. \end{align}

Consequently, system (7.71) becomes

(7.73) \begin{align} \begin{split} \mathfrak{u}_{2,t}+{} &\frac {1}{ \rho _2 } \mathfrak{u}_2 \mathfrak{u}_{2,x} +\rho _2 g\eta _{2,x}=0 , \\ \mathfrak{u}_{3,t}+{} &\frac {\rho _2-\rho _3}{\rho _2 \rho _3 }\mathfrak{u}_2 \mathfrak{u}_{2 ,x} + \frac {1}{\rho _3}(\mathfrak{u}_2\mathfrak{u}_3)_x+ \frac {1}{\rho _3}\mathfrak{u}_3 \mathfrak{u}_{3 ,x} +(\rho _3-\rho _2)g\eta _{3,x}=0 ,\\ \eta _{2,t}={}&- \frac {h_3\rho _3 + h_b\rho _2}{ \rho _2\rho _3} \mathfrak{u}_{2,x} - \frac { h_3^3 \rho _3^2 + 3h_3^2h_b\rho _2 \rho _3 + 3h_3h_b^2\rho _2^2 + h_b^3\rho _2 \rho _3}{3\rho _2 \rho _3 ^2} \mathfrak{u}_{2,\textit{xxx}}\\ &{}-\frac {1}{\rho _2} [\eta _2 \mathfrak{u}_{2}]_x - \frac {\rho _2-\rho _3 }{\rho _2\rho _3 } [\eta _3 \mathfrak{u}_{2} ]_x \\ &{}- \frac { h_b }{\rho _3 } \mathfrak{u}_{3,x} - \frac { h_b \big[6h_bh_3 \rho _2+ \big(3 h_3^2 + 2h_b^2\big) \rho _3 \big]}{6 \rho _3^2} \mathfrak{u}_{3,\textit{xxx}} - \frac { 1 }{\rho _3 } [\eta _3 \mathfrak{u}_{3} ]_x ,\\ \eta _{3,t}={}&- \frac { h_b }{\rho _3 } \mathfrak{u}_{2,x} - \frac { h_b \big[6h_bh_3 \rho _2+ \big(3 h_3^2 + 2h_b^2\big) \rho _3 \big]}{6 \rho _3^2} \mathfrak{u}_{2,\textit{xxx}} - \frac { 1 }{\rho _3 } [\eta _3 \mathfrak{u}_{2} ]_x \\ & {}-\frac {h_b }{\rho _3} \mathfrak{u}_{3,x} -\frac {h_b^2 ( 3 h_3\rho _2 + h_b \rho _3 )}{3 \rho _3 ^2} \mathfrak{u}_{3,\textit{xxx}} - \frac {1}{\rho _3 }[\eta _3 \mathfrak{u}_{3} ]_x , \end{split} \end{align}

which recovers the system obtained by Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005) on p. 1626, derived by employing similar methods involving expansions of the DN operators for two layers with a free surface.

The procedure utilised to obtain the Boussinesq system (7.73) highlights a general method for the derivation of nonlinear models for water waves propagation starting from the general system (7.17), (7.20) and (7.21).

7.6. Limit to the full equations for two layers with a free surface

The limit to the full equations for two layers with a free surface as they appear (however, in Hamiltonian form) in Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005) could be done by simply setting up $\rho _1=0$ in all formulas in this section. Furthermore, given a depth $h^*$ , there is a relation between the DN operators $G_{\textit{ij}}(h^*, \eta _l, \eta _u)$ and $\varGamma _{\textit{ij}}(h^*, \eta _l, \eta _u)$ , which reads

(7.74) \begin{align} G_{11}(h^*, \eta _l, \eta _u)& =\varGamma _{22}(h^*, \eta _l, \eta _u) , \quad G_{12}(h^*, \eta _l, \eta _u)=\varGamma _{21}(h^*, \eta _l, \eta _u) , \\[-12pt] \nonumber \end{align}

(7.75) \begin{align} G_{21}(h^*, \eta _l, \eta _u)& =\varGamma _{12}(h^*, \eta _l, \eta _u) , \quad G_{22}(h^*, \eta _l, \eta _u)=\varGamma _{11}(h^*, \eta _l, \eta _u), \\[9pt] \nonumber \end{align}

where $y=\eta _l(x,t)$ is the lower boundary of the domain, and $y=\eta _u(x,t)$ is the upper boundary of the domain.

From (7.6), we obtain

(7.76) \begin{align} \xi _2:=\rho _2(\varphi _2)_{s_2} , \quad \xi _3:=\rho _3(\varphi _3)_{s_3}-\rho _2(\varphi _2)_{s_3}. \end{align}

We define the operator (essentially the operator $B$ in Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005) for the free surface case.)

(7.77) \begin{align} \mathcal{B}:=\rho _3 \varGamma _{22}+\rho _2 G, \end{align}

which is self-adjoint. Then we have

(7.78) \begin{align} B^{-1}= \frac {1}{\rho _2}G_1^{-1}\varGamma _{12} \mathcal{B} ^{-1 }, \quad ( B^{-1})^*= \frac {1}{\rho _2} \mathcal{B} ^{-1 } \varGamma _{21}G_1^{-1 }. \end{align}

From (7.7) and (7.8), we obtain

(7.79) \begin{align} (\varphi _2)_{s_2}&=\frac {1}{\rho _2} \xi _2 =:\varPhi ^2_2(\xi _2) , \\[-12pt] \nonumber \end{align}

(7.80) \begin{align} (\varphi _2)_{s_3}&=-\mathcal{B} ^{-1 } \left [ G\xi _3 + \frac {\rho _3}{\rho _2} \varGamma _{21} \xi _2 \right ]=:\varPhi ^2_3(\eta _i, \xi _i). \\[9pt] \nonumber \end{align}

Then from (7.5), (7.79) and (7.80) after some calculations, we have

(7.81) \begin{align} (\varphi _3)_{s_3}=-G^{-1}[\varGamma _{21}(\varphi _2)_{s_2}+\varGamma _{22}(\varphi _2)_{s_3}]=\mathcal{B} ^{-1 } \left [-\varGamma _{21}\xi _2+\varGamma _{22}\xi _3 \right ]=:\varPhi ^3(\eta _i, \xi _i). \end{align}

In view of (7.74)–(7.75), these formulas reproduce (2.26)–(2.28) of Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005). Next, from (7.12)–(7.15), we obtain the expressions for the operators:

(7.82) \begin{align} A_{22}&= \frac {1}{\rho _2}\varGamma _{11}-\frac {\rho _3}{\rho _2}\varGamma _{12}\mathcal{B} ^{-1 }\varGamma _{21}, \\[-12pt] \nonumber \end{align}

(7.83) \begin{align} A_{23}&= -\varGamma _{12}\mathcal{B} ^{-1 }G, \\[-12pt] \nonumber \end{align}

(7.84) \begin{align} A_{32}&= \frac {1}{\rho _2}\left (-\mathcal{B} +{\rho _3}\varGamma _{22}\right )\mathcal{B} ^{-1 }\varGamma _{21}=-G\mathcal{B} ^{-1 } \varGamma _{21}=A_{23}^* , \\[-12pt] \nonumber \end{align}

(7.85) \begin{align} A_{33}&= \varGamma _{22}\mathcal{B} ^{-1 }G. \\[9pt] \nonumber \end{align}

The operator $A_{33}$ is self-adjoint since $G$ and $\varGamma _{22}$ are self-adjoint operators (Craig et al. Reference Craig, Guyenne, Nicholls and Sulem2005). Indeed, we can write it formally as

(7.86) \begin{align} A_{33}=\frac {1}{\rho _3} \varGamma _{22}\left ( 1+\frac {\rho _2}{\rho _3} \varGamma _{22}^{-1} G \right )^{-1}\varGamma _{22}^{-1}G. \end{align}

In the leading order, $ ( {\rho _2}/{\rho _3})\,| \varGamma _{22}^{-1} G|= ({\rho _2}/{\rho _3})\, |{\tanh (kh_b)} \tanh (k h_3)|\lt 1$ , so, at least in some domain for the values of $(k,\eta _2,\eta _3)$ in the neighbourhood of $(0,0,0)$ , we have an expansion

(7.87) \begin{align} A_{33}=\frac {1}{\rho _3} \left ( G - \frac {\rho _2}{\rho _3} G \varGamma _{22}^{-1 }G + \frac {\rho _2^2}{\rho _3^2} G \varGamma _{22}^{-1 }G \varGamma _{22}^{-1 }G+\cdots \right )\!, \end{align}

from which it is evident that $A_{33}=A_{33}^*$ , thus $A_{33}= G \mathcal{B} ^{-1 } \varGamma _{22}$ . The operator $A_{22}$ is also self-adjoint since $\varGamma _{11}$ is self adjoint and $\varGamma _{12}^*=\varGamma _{21}$ ; cf. Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005).

The functionals $F_2$ and $F_3$ are as in (7.16); then the full equations of motion are the same as in (7.17), (7.20) and (7.21) (with $\rho _1=0$ ), i.e. explicitly,

(7.88) \begin{align} \eta _{2,t}&=\left ( \frac {1}{\rho _2}\varGamma _{11}-\frac {\rho _3}{\rho _2}\varGamma _{12}\mathcal{B} ^{-1 }\varGamma _{21} \right ) \xi _2 -\varGamma _{12}\mathcal{B} ^{-1 }G \xi _3 =: F_2(\eta _i,\xi _i) , \end{align}

(7.89) \begin{align} \eta _{3,t}&= G\mathcal{B} ^{-1 }[ -\varGamma _{21} \xi _2+ \varGamma _{22} \xi _3]=:F_3(\eta _i, \xi _i), \end{align}

(7.90) \begin{align} \xi _{2,t}&=-\frac {\rho _2 (\varPhi ^2_{2,x})^2-2F_2\xi _{2,x}\eta _{2,x}-\rho _2F_2^2 }{2(1+\eta _{2,x}^2)}-\rho _2g \eta _2 , \end{align}

(7.91) \begin{align} \xi _{3,t}&=-\frac {\rho _3 (\varPhi ^3_{x})^2-\rho _2 (\varPhi _{3,x}^2)^2-2F_3\xi _{3,x}\eta _{3,x}+(\rho _2-\rho _3)F_3^2 }{2(1+\eta _{3,x}^2)}-(\rho _3-\rho _2)g \eta _3 . \end{align}

Equation (7.88) follows also from (7.4), (7.79) and (7.80); and (7.89) follows from (7.5) and (7.81). The DN operators depend on $\eta _2$ and $\eta _3$ in a complicated nonlinear way; however, as we have mentioned, there is a systematic approach for their derivation. In the Boussinesq limit, the above system produces the system (7.73) (see the derivation in Craig et al. Reference Craig, Guyenne, Nicholls and Sulem2005), although we have obtained it as a limit from (7.71). Figure 2 shows the relations between the derived models in this section.

Figure 2.

Relations between the derived models: limits and approximations.