Integrable dispersive generalisation of Ovsyannikov two-layer shallow water model

Abstract

Ovsyannikov equations describe long-wavelength dynamics of nonlinear waves in two-layer shallow water systems. We suggest their simple generalisations which model qualitatively dispersive properties of such systems. These generalised equations are completely integrable by the known methods of nonlinear physics and we present their periodic and soliton solutions. The theory is illustrated by concrete examples.

1. Introduction

Internal gravity waves at the interface of two fluids with different densities have been studied for a long time in a number of papers (see, e.g., Ref. [42] and references therein). As is well known, a surface separating two layers is stable when the density $\rho_2$ of the upper layer is smaller than the density $\rho_1$ of the lower layer, $\rho_2 \lt \rho_1$, so the propagation of the linear interfacial waves is also stable. It is clear that in such a two-layer system there exist two modes of wave propagation – one with both fluids moving mainly in phase and another one with counter-phase motion of two fluids. The first mode reduces to the usual surface gravity waves in the limit $\rho_1-\rho_2\ll\rho_1$ and the second mainly interfacial mode is expected to demonstrate relatively slow motions since the difference between gravity forces acting on each layer is very small compared with the total gravity force. This means that we can separate approximately these two motions and study the interfacial waves independently of the surface waves in the main approximation with respect to the small parameter $\varepsilon=(\rho_1-\rho_2)/\rho_1\ll1$. Actually, existence of this small parameter allows one to introduce the slow time variable $T=\epsilon^{1/2} t$ and derive equations for the slow dynamics which govern the propagation of the interfacial waves. This idea was realized by Ovsyannikov in Ref. [39] for his “model-III” of shallow water dynamics of two-layer flows in the long wavelength (dispersionless) limit (see also Refs. [4, 6]). On the basis of this model, the theory of turbulent bores propagating in two-layer flows was developed in Refs. [34, 35].

As is known, besides turbulent bores, there exist solitons as well as undular bores (or dispersive shock waves (DSWs)) which connect two states of the flow with different parameters, where a transition region of strong dissipation is replaced by an expanding with time region of strong nonlinear oscillations (see, e.g., review articles [14, 22] and references therein). Gurevich and Pitaevskii suggested the general approach [19], based on the Whitham theory [44, 45] of modulations of nonlinear waves, to the analytical description of DSWs. This theory was successfully applied to the theory of undular bores described by different models of surface gravity waves (see, e.g., [7, 9, 11–13, 17, 31, 32, 43, 46]) and to some models of the two-layer flows (see, e.g., [15, 28, 36–38]). However, it was not discussed so far for any dispersive generalisation of the Ovsyannikov model. The goal of this paper is to introduce a completely integrable generalisation of the Ovsyannikov ‘model III’ of the two-layer flows and to study its soliton solutions. This is the first step to developing the theory of DSWs in the two-layer shallow water model.

To obtain the desirable generalisation, we should notice the following. In a simpler case of shallow water equations, there are several physically reasonable dispersive generalisations. First, the shallow water equations represent the dispersionless limit of the Gross-Pitaevskii (or defocusing nonlinear Schrödinger (NLS)) equation [3, 16, 40] for Bose-Einstein condensates (BECs) flows, which can be considered as one of their dispersive generalizations. Its one-dimensional version is completely integrable [47] and the corresponding analytical theory of DSWs was developed in Refs. [10, 18]. Second, the shallow water equations are the well-known dispersionless limit of the Boussinesq family of equations for surface waves dynamics (see, e.g., [45]), and a particular case of the Kaup-Boussinesq equations [30] has an advantage of being the completely integrable system which is the most convenient for development of the analytical theory of DSWs [7, 11, 12]. It is important that both these dispersive generalisations correspond to the same dispersion relation for linear harmonic waves.

If we turn to the Ovsyannikov model, then we notice that it is the dispersionless limit of the Landau-Lifshitz equation for a magnetic with an easy-plane anisotropy [20] or an equivalent model for a two-component BEC with very small difference of inter-species interaction constants [8, 41], and the corresponding theory of DSWs was developed in [21]. This is the completely integrable generalisation of the Ovsyannikov equations analogous to the transition from the shallow water equations to the NLS equation. Another generalisation of the Ovsyannikov equations similar to the transition from the shallow water equations to the Kaup-Boussinesq system has been found recently in Ref. [26]. The dispersion relation for linear harmonic waves does not coincide exactly with the real dispersion relation for interfacial waves in two-layer flows, so this model only gives a qualitative description of solitons or undular bores in such systems. It is also worth mentioning that the two main Ovsyannikov’s model approximations – the shallow water approximation and separation of two modes of wave propagation – exclude some important physical effects which exist in the two-fluid dynamics. In particular, this model does not describe large amplitude solitons for which ‘fully-nonlinear’ theories have been developed (see, e.g., [2, 5, 15]). Since we neglect interaction between the two modes of propagation, the resonance effects appearing in case of equality of the corresponding phase velocities are also excluded from our consideration (see, e.g., [1]). In spite of these reservations, the presented here model has an advantage of being completely integrable. So it is amenable to a quite detailed analytical treatment and can describe some specific properties of two-fluid systems qualitatively, as it will be demonstrated in this paper.

2. Ovsyannikov equations

Let $h_1=h_1(x,t)$ and $h_2=h_2(x,t)$ denote local depths of the lower and upper layers of our fluid system, correspondingly. We suppose that these fluids are incompressible and the density $\rho_2$ of the upper layer is slightly smaller than the density $\rho_1$ of the lower layer, so we have a small parameter

(2.1)\begin{equation} \varepsilon=1-\rho_2/\rho_1\ll1. \end{equation}

The wavelength of waves is supposed to be much longer than both depths, so in the standard shallow water approximation the flows are characterised by the averaged across each layer horizontal velocities $u_1=u_1(x,t)$ and $u_2=u_2(x,t)$, hence, we have four dynamical variables $h_1,h_2,u_1,u_2$ in total. The conservation of mass in each layer gives evident equations

(2.2)\begin{equation} h_{1,t}+(h_1u_1)_x=0,\qquad h_{2,t}+(h_2u_2)_x=0, \end{equation}

and it is not difficult to find that conservation of momenta in the layers yields the equations

(2.3)\begin{equation} \begin{aligned} & u_{1,t}+u_1u_{1,x}+gh_{1,x}+(1-\varepsilon)gh_{2,x}=0,\\ & u_{2,t}+u_2u_{2,x}+gh_{1,x}+gh_{2,x}=0, \end{aligned} \end{equation}

where $g$ is the gravity acceleration constant. To simplify the notation, we will assume in what follows that $g=1$.

It is easy to find from the system of Eqs. (2.2) and (2.3), that the characteristic velocity $v_c$ must satisfy the algebraic equation

(2.4)\begin{equation} (u_1-v_c)^2(u_2-v_c)^2-(u_1-v_c)^2h_2-(u_2-v_c)^2h_1+\varepsilon h_1h_2=0. \end{equation}

To understand its physical consequences, let us consider first a particular case of waves propagating along the background with $u_1=u_2=0$. Then the equation $v_c^4-(h_1+h_2)v_c^2+\varepsilon h_1h_2=0$ gives in the main approximation with respect to the small parameter $\varepsilon$ two roots

(2.5)\begin{equation} |v_{c,\text{fast}}|\approx\sqrt{h_1+h_2}\quad\text{and}\quad |v_{c,\text{slow}}|\approx\sqrt{\frac{\varepsilon h_1h_2}{h_1+h_2}}. \end{equation}

It is clear that the fast velocity $v_{c,\text{fast}}$ corresponds to the in-phase motions in the layers equivalent to shallow waves in a single-layer liquid with the total depth $h_1+h_2$. In this case there is no exchange of momentum between the layers in the main approximation. The slow velocity $v_{c,\text{slow}}$ corresponds to the counter-phase motions in the layers with exchange of the momentum across the interface between the layers. This is the characteristic velocity of the interfacial waves in this particular case with $u_1=u_2=0$. Since we are interested in these interfacial waves, we assume that in this case all velocities have the same order of magnitude, that is $(u_1-v_c)^2\sim(u_2-v_c)^2\sim\varepsilon$. Then the first term in Eq. (2.4) can be neglected in the main approximation and we obtain

(2.6)\begin{equation} v_c=\frac{1}{h_1+h_2}\left\{u_1h_2+u_2h_1\pm\sqrt{h_1h_2[\varepsilon(h_1+h_2)-(u_1-u_2)^2]}\right\}. \end{equation}

These characteristic velocities are real, that is the motion is stable, if

(2.7)\begin{equation} |u_1-u_2| \lt \sqrt{\varepsilon(h_1+h_2)}. \end{equation}

The obtained above formulas define the characteristic scales of the interfacial wave phenomena. Therefore, we introduce some characteristic length $H$ (its exact value will be determined later), so

(2.8)\begin{equation} h_1=\overline{h}_1 H,\qquad h_2=\overline{h}_2H, \end{equation}

define the scaled velocities

(2.9)\begin{equation} u_1=\sqrt{\varepsilon H}\,\overline{u}_1,\qquad u_2=\sqrt{\varepsilon H}\,\overline{u}_2, \end{equation}

and the corresponding space and time variables

(2.10)\begin{equation} \overline{x}=\frac{x}{H},\qquad \overline{t}=\sqrt{\frac{\varepsilon}{H}}\,t. \end{equation}

In new variables, Eqs. (2.3) take the form

(2.11)\begin{equation} \begin{aligned} & \varepsilon\left(\overline{u}_{1,\overline{t}}+\overline{u}_1\overline{u}_{1,\overline{x}}-\overline{h}_{2,\overline{x}}\right)+(\overline{h}_1+\overline{h}_2)_{\overline{x}}=0,\\ & \varepsilon\left(\overline{u}_{2,\overline{t}}+\overline{u}_2\overline{u}_{2,\overline{x}}\right)+(\overline{h}_1+\overline{h}_2)_{\overline{x}}=0. \end{aligned} \end{equation}

Their difference gives for small but finite $\varepsilon$ the equation

(2.12)\begin{equation} (\overline{u}_1-\overline{u}_2)_{\overline{t}}+\frac12(\overline{u}_1^2-\overline{u}_2^2-2\overline{h}_2)_{\overline{x}}=0. \end{equation}

In the main approximation in the limit of small $\varepsilon$ each equation (2.11) reduces to $(\overline{h}_1+\overline{h}_2)_{\overline{x}}=0$, that is $\overline{h}_1+\overline{h}_2=a(\overline{t})$. The sum of Eqs. (2.2) gives

(2.13)\begin{equation} a_{\overline{t}}+(\overline{h}_1\overline{u}_1+\overline{h}_2\overline{u}_2)_{\overline{x}}=0. \end{equation}

We assume that there is no sources or sinks of the liquids, so the total flux is constant, $\overline{h}_1\overline{u}_1+\overline{h}_2\overline{u}_2=b=\mathrm{const}$, and then $a=\mathrm{const}$. We choose its value equal to $a=1$, that is

(2.14)\begin{equation} \overline{h}_1+\overline{h}_2=1\quad\text{or}\quad h_1+h_2=H. \end{equation}

This means that in this approximation the interfacial waves do not disturb the total depth of the two-layer system and it remains equal to $H=\mathrm{const}$. This explains the meaning of the parameter $H$. Now, we can use the reference frame where the constant total flux $b$ is equal to zero,

(2.15)\begin{equation} \overline{h}_1\overline{u}_1+\overline{h}_2\overline{u}_2=0 \quad \text{or} \quad h_1u_1+h_2u_2=0. \end{equation}

Thus, we have two constrains (2.14) and (2.15), so the interfacial waves are described by only two variables. We define them by the formulas

(2.16)\begin{equation} v=\overline{u}_1-\overline{u}_2,\qquad w=\overline{h}_2-\overline{h}_1, \end{equation}

so that

(2.17)\begin{equation} \begin{aligned} & \overline{h}_1=\frac12(1-w),\qquad \overline{h}_2=\frac12(1+w),\\ & \overline{u}_1=\frac12v(1+w),\qquad \overline{u}_2=-\frac12v(1-w), \end{aligned} \end{equation}

and the conditions (2.14), (2.15) are satisfied automatically. To find the dynamical equations for the variables $v$ and $w$, we substitute Eqs. (2.17) into the first Eq. (2.3), written in terms of $\overline{h}_1,\overline{u}_1$, and into Eq. (2.12). After evident simplifications, we arrive at the equations

(2.18)\begin{equation} w_t-\frac12[v(1-w^2)]_x=0,\qquad v_t-\frac12[w(1-v^2)]_x=0, \end{equation}

where we use for $\overline{t}$ and $\overline{x}$ the previous notation $t$ and $x$. Eqs. (2.18) are the Ovsyannikov equations that govern the wave dynamics of interfacial waves in the long wavelength (dispersionless) limit.

In new variables, the characteristic velocities (2.6) can be written as

(2.19)\begin{equation} v_c=\sqrt{\varepsilon H}\,v_{\pm},\quad v_{\pm}=vw\pm\frac12\sqrt{(1-v^2)(1-w^2)}. \end{equation}

Then the system (2.18) can be diagonalised and in terms of the Riemann invariants

(2.20)\begin{equation} r_{\pm}=\frac12vw\pm\frac12\sqrt{(1-v^2)(1-w^2)}, \end{equation}

it takes the form

(2.21)\begin{equation} \frac{\partial r_+}{\partial t}+v_+\frac{\partial r_+}{\partial x}=0,\qquad \frac{\partial r_-}{\partial t}+v_-\frac{\partial r_-}{\partial x}=0, \end{equation}

where

(2.22)\begin{equation} v_+=\frac12(3r_++r_-),\qquad v_-=\frac12(r_++3r_-). \end{equation}

Eq. (2.18) are hyperbolic in the region

(2.23)\begin{equation} -1\leq v\leq1,\qquad -1\leq w\leq1, \end{equation}

so the system (2.20) can be solved by the standard hodograph method with reservation that the mapping $(v,w)\to(r_+,r_-)$ is degenerate at the diagonals of this square, and therefore the hodograph transform solution yields different physical solutions in four triangles obtained after cutting the square by its diagonals. This leads to a quite complicated classification of different DSWs appearing in such not genuinely nonlinear systems (see, e.g., Refs. [15, 21, 24]).

Now we will discuss the completely integrable dispersive generalisation of the Ovsyannikov equations.

3. Dispersive generalisation

We want to generalise Eqs. (2.18) in such a way, that the resulting equations are completely integrable and reproduce in the long wavelength limit the dispersion relation for linear interfacial waves at least qualitatively. In particular, the main term in the series expansion of the dispersion relation $\omega=\omega(k)$ in the limit $k\to0$, $k$ being the wave number, must be

(3.1)\begin{equation} \omega\approx k\left(vw\pm\frac12\sqrt{(1-v^2)(1-w^2)}\right), \end{equation}

so that the phase velocity $\omega/k$ becomes equal to the dispersionless characteristic velocities (2.19). The method of finding such generalisations was suggested in Refs. [26, 27]. This method is based on the condition of asymptotic integrability of wave equations which can be formulated in the following way. The wave equation is called asymptotically integrable, if the Hamilton equations governing the propagation of short-wavelength wave packets with a carrier wave number $k$ along a large-scale background wave $r_+=r_+(x,t),r_-=r_-(x,t)$ have an integral $Q(x,k,r_+,r_-)=\mathrm{const}$ or $k=k(x,r_+,r_-,q)$, where $q$ is some constant. If such an integral exists, then the function $k$ satisfies the equations

(3.2)\begin{equation} \frac{\partial k}{\partial r_+}=\frac{\partial\omega/\partial r_+}{v_+-\partial\omega/\partial k},\qquad \frac{\partial k}{\partial r_-}=\frac{\partial\omega/\partial r_-}{v_--\partial\omega/\partial k}. \end{equation}

This condition of asymptotic integrability imposes strong restrictions upon the form of the dispersion relation $\omega=\omega(k,r_+,r_-)$. In fact, only a few dispersion relation which satisfy Eqs. (3.2) are known (see Refs. [26, 27]), and we choose for our generalisation of the Ovsyannikov equations the dispersion relation

(3.3)\begin{equation} \omega=k\left(r_++r_-\pm\frac12\sqrt{(r_+-r_-)^2+\sigma k^2}\right), \end{equation}

where $\sigma=\pm1$. The corresponding integral, that is the solution of Eqs. (3.2), is given by the formula

(3.4)\begin{equation} k^2=4\sigma(q-r_+)(q-r_-). \end{equation}

In case of the Ovsyannikov equations with the Riemann invariants given by Eqs. (2.20) the dispersion relation reads

(3.5)\begin{equation} \omega=k\left(vw\pm\frac12\sqrt{(1-v^2)(1-w^2)+\sigma k^2}\right), \end{equation}

and the integral (3.4) takes the form

(3.6)\begin{equation} k^2=\sigma[(\lambda-vw)^2-(1-v^2)(1-w^2)], \end{equation}

where we replaced $q$ by a more common notation $\lambda=2q$.

The dispersion relation (3.5) coincides exactly (for $\sigma=1$) with the dispersion relation for linear waves in the two-component Bose-Einstein condensate [8, 21], so the knowledge of the integral (3.5) leads to some important consequences, for example, to the theory of motion of so-called ‘magnetic’ solitons in such condensates [25, 41]. In our case of the two-layer flows the dispersion relation obviously satisfies the condition (3.1), and this guarantees the correct description of the evolution in the dispersionless approximation. The dependence of the expression under the square root sign on $k^2$ models the dispersive properties. The actual dispersion relation for linear waves in our two-layer model is more complicated (see Appendix A). Nevertheless, it has the same general structure, so we hope that our model provides a qualitatively correct description of properties of nonlinear waves in two-layer shallow water systems. At the same time it should be noticed that Ovsyannikov’s scaling (2.9), (2.10) implies that $\varepsilon$ is the smallest parameter of the theory, so we cannot find the dependence of nonlinear wave structure on $\varepsilon$, as it was done in the theories without taking this limit $\varepsilon\to0$; see, e.g., Ref. [1] and references therein. Actually, we consider the limiting case of waves at the critical value of the total depth $H=1$ (in non-dimensional variables) corresponding to the negligibly small difference of fluids’ densities $\varepsilon\to0$.

To get the nonlinear wave equations for wave dynamics in our model, we assume that the integral (3.6) is just a result of the quasiclassical approximation applied to the linear spectral problem [26, 27]

(3.7)\begin{equation} \phi_{xx}=\mathcal{A}\phi, \end{equation}

where

(3.8)\begin{equation} \mathcal{A}=-\frac14k^2=\frac{\sigma}{4}[(1-v^2)(1-w^2)-(\lambda-vw)^2]. \end{equation}

Since $v$ and $w$ depend also on time $t$, the function $\phi$ depends on $t$, too, and, as was shown in Refs. [23, 26, 27], this dependence is determined by the equation

(3.9)\begin{equation} \phi_t=-\frac12\mathcal{B}_x\phi+\mathcal{B}\phi_x, \end{equation}

where

(3.10)\begin{equation} \mathcal{B}=-\frac{\omega}{k}=-\frac12(\lambda+vw). \end{equation}

The compatibility condition $(\phi_{xx})_t=(\phi_t)_{xx}$ of Eqs. (3.7) and (3.9),

(3.11)\begin{equation} \mathcal{A}-2\mathcal{B}_x\mathcal{A}-\mathcal{B}\mathcal{A}_x+\frac12\mathcal{B}_{xxx}=0, \end{equation}

in our particular case (3.8) and (3.10) yields the system

(3.12)\begin{equation} \begin{aligned} & v_t-\frac12[w(1-v^2)]_x-\frac{\sigma v(vw)_{xxx}}{2(v^2-w^2)}=0,\\ & w_t-\frac12[v(1-w^2)]_x+\frac{\sigma w(vw)_{xxx}}{2(v^2-w^2)}=0. \end{aligned} \end{equation}

This is the desirable dispersive generalisation of the Ovsyannikov system (2.18). Thus, we have obtained the system of Eqs. (3.12) and the corresponding Lax pair (3.8), (3.10), so we can apply for its investigation the well-developed methods.

As was shown in Ref. [39], the Ovsyannikov system (2.18) can be transformed to the shallow water equations

(3.13)\begin{equation} \rho_t+(\rho u)_x=0,\qquad u_t+uu_x+\rho_x=0. \end{equation}

In fact, the same transformation casts the system (3.12) to the Kaup-Boussinesq system [30]. Indeed, we identify the Riemann invariants $r_{\pm}=u/2\pm\sqrt{\rho}$ with the Riemann invariants (2.20) of the Ovsyannikov equations and obtain the relationships

(3.14)\begin{equation} u=vw,\qquad \rho=\frac14(1-v^2)(1-w^2). \end{equation}

Then it is easy to show that $u$ and $\rho$ must satisfy the Kaup-Boussinesq equations

(3.15)\begin{equation} \rho_t+(\rho u)_x-\frac{\sigma}4u_{xxx}=0,\qquad u_t+uu_x+\rho_x=0, \end{equation}

which corresponds to the Lax pair (3.7), (3.9) with

(3.16)\begin{equation} \mathcal{A}=\frac{\sigma}{4}[\rho-(\lambda-u)^2],\qquad \mathcal{B}=-\frac12(\lambda+u). \end{equation}

Periodic and soliton solutions of the Kaup-Boussinesq system (3.15) for both signs $\sigma=\pm1$ of the dispersion term were discussed in Refs. [7, 11, 12], and they can be transformed to the Ovsyannikov wave variables. To present the results in convenient for us form, we will do it independently in the next section.

4. Periodic and soliton solutions

To be definite, we confine ourselves to the upper branch of the dispersion relation (3.5) which corresponds to the negative dispersion (see Appendix), so we take the model (3.12) with $\sigma=-1$. For finding its periodic and soliton solutions, we will use the ‘reduced’ finite-gap integration method (see, e.g., [24]).

The second-order linear spectral problem (3.7) has two-independent basis solutions $\phi_+,\phi_-$, so we define the ‘squared basis function’ $g=\phi_+\phi_-$. It is well known that it satisfies the differential equations

(4.1)\begin{equation} g_{xxx}-2\mathcal{A}_xg-4\mathcal{A}g_x=0,\qquad g_t=\mathcal{B}g_x-\mathcal{B}_xg, \end{equation}

which follow from Eqs. (3.7) and (3.9), correspondingly. The first equation is easily integrated to give

(4.2)\begin{equation} \frac12gg_{xx}-\frac14g_x^2-\mathcal{A}g^2=-P(\lambda), \end{equation}

and (quasi)periodic solutions are distinguished in the finite-gap integration method by the condition that the integration constant $P(\lambda)$ is a polynomial function of the spectral parameter $\lambda$. In our case, $\mathcal{A}$ is a second-degree polynomial in $\lambda$, so a single-phase solution is obtained, when we choose a linear in $\lambda$ function $g$,

(4.3)\begin{equation} g=2(\lambda-\mu), \end{equation}

and then $P(\lambda)$ must be a 4th degree polynomial in $\lambda$. As a result, Eq. (4.2) becomes the identity

(4.4)\begin{equation} (\lambda-\mu)\mu_{xx}+\mu_x^2+\left[(1-v^2)(1-w^2)-(\lambda-vw)^2\right](\lambda-\mu)^2=P(\lambda), \end{equation}

where

(4.5)\begin{equation} P(\lambda)=\prod_{i=1}^4(\lambda-\lambda_i)=\lambda^4-s_1\lambda^3+s_2\lambda^2-s_3\lambda+s_4, \end{equation}

and $s_k$ are symmetrical functions of the zeros $\lambda_i$,

(4.6)\begin{equation} s_1=\sum\lambda_i,\quad s_2=\sum_{i \lt j}\lambda_i\lambda_j,\quad s_3=\sum_{i \lt j \lt k}\lambda_i\lambda_j\lambda_k,\quad s_4=\lambda_1\lambda_2\lambda_3\lambda_4. \end{equation}

Comparing coefficients before different degrees of the spectral parameter $\lambda$, we obtain useful formulas relating our wave variables $v,w$ with the function $\mu$. In particular, we get

(4.7)\begin{equation} vw=\frac12s_1-\mu\equiv a(\mu),\qquad v^2+w^2=1+s_2-2s_1\mu+3\mu^2\equiv b(\mu), \end{equation}

so $v^2$ and $w^2$ are expressed in terms of $\mu$ by the formulas

(4.8)\begin{equation} v^2,w^2=\frac12\left[b(\mu)\pm\sqrt{b^2(\mu)-4a^2(\mu)}\right]. \end{equation}

Consequently, if $\mu=\mu(x,t)$ is found, then $v$ and $w$ are also known.

The second Eq. (4.1) gives with account of Eqs. (3.16) and (4.3) $\mu_t=-\frac14s_1\mu_x$, that is $\mu$ is only a function of $\xi=x-Vt$, $V=s_1/4$. Eq. (4.4) after substitution $\lambda=\mu$ yields $\mu_x^2=P(\mu)$, so we arrive at

(4.9)\begin{equation} \mu_{\xi}=\pm\sqrt{P(\mu)}. \end{equation}

This equation can be solved by standard methods. If the zeroes $\lambda_i$ are ordered according to inequalities

(4.10)\begin{equation} -1\leq\lambda_1\leq\lambda_2\leq\lambda_3\leq\lambda_4\leq1, \end{equation}

where the conditions $-1\leq\lambda_1$ and $\lambda_4\leq1$ are imposed in order to ensure that $0\leq v^2,w^2\leq1$, then a periodic solution corresponds to oscillations of $\mu$ in the interval

(4.11)\begin{equation} \lambda_2\leq\mu\leq\lambda_3, \end{equation}

where $P(\mu)\geq0$. The integral $\xi=\pm\int d\mu/\sqrt{P(\mu)}$ is equal to the elliptic integral of the first kind and its inversion yields the periodic solution $\mu=\mu(\xi)$ expressed in terms of elliptic functions. It is convenient to represent the result in two different forms corresponding to different initial conditions for the function $\mu(\xi)$.

The first form corresponds to the condition $\mu(0)=\lambda_2$. Then standard calculation yields

(4.12)\begin{equation} \mu=\frac{\lambda_2(\lambda_3-\lambda_1)-\lambda_1(\lambda_3-\lambda_2)\mathrm{sn}^2(W,m)}{\lambda_3-\lambda_1-(\lambda_3-\lambda_2)\mathrm{sn}^2(W,m)}, \end{equation}

where

(4.13)\begin{equation} W=\frac12\sqrt{(\lambda_3-\lambda_1)(\lambda_4-\lambda_2)}\,\xi,\qquad m=\frac{(\lambda_4-\lambda_1)(\lambda_3-\lambda_2)}{(\lambda_4-\lambda_2)(\lambda_3-\lambda_1)}. \end{equation}

This form of the solution is convenient for transition to the soliton limit $\lambda_3\to\lambda_4$, when $m\to1$, so we get

(4.14)\begin{equation} \mu=\frac{\lambda_2(\lambda_4-\lambda_1)-\lambda_1(\lambda_4-\lambda_2)\tanh^2W}{\lambda_4-\lambda_1-(\lambda_4-\lambda_2)\tanh^2W}, \qquad W=\frac12\sqrt{(\lambda_4-\lambda_1)(\lambda_4-\lambda_2)}\,\xi. \end{equation}

This is a dark soliton solution for the variable $\mu$ which reaches its minimal value at $\xi=0$.

Another form of the solution corresponds to the initial condition $\mu(0)=\lambda_3$. Then we get

(4.15)\begin{equation} \mu=\frac{\lambda_3(\lambda_4-\lambda_2)-\lambda_4(\lambda_3-\lambda_2)\mathrm{sn}^2(W,m)}{\lambda_4-\lambda_2-(\lambda_3-\lambda_2)\mathrm{sn}^2(W,m)}. \end{equation}

Actually, this is the same solution (4.12), but shifted by a half-wavelength, so the origin $\xi=0$ is shifted from the wave trough to the wave crest. Consequently, $W$ and $m$ are given by the same formulas (4.13). Now the soliton limit $\lambda_2\to\lambda_1$, $m\to1$, gives a bright soliton solution for the variable $\mu$:

(4.16)\begin{equation} \mu=\frac{\lambda_3(\lambda_4-\lambda_2)-\lambda_4(\lambda_3-\lambda_1)\tanh^2W}{\lambda_4-\lambda_1-(\lambda_3-\lambda_2)\tanh^2W}, \qquad W=\frac12\sqrt{(\lambda_4-\lambda_1)(\lambda_3-\lambda_1)}\,\xi. \end{equation}

Solitons with different polarities in a two-layer systems were considered in Ref. [1] without restriction to the limit of very small difference of fluids’ densities.

Besides soliton solutions, there is also a kink solution for $\lambda_2=\lambda_1,\lambda_3=\lambda_4$:

(4.17)\begin{equation} \mu=\lambda_4-\frac{\lambda_4-\lambda_1}{\exp[\pm(\lambda_4-\lambda_1)\xi]+1}, \end{equation}

which connects two different plateau states with $\mu=\lambda_1$ and $\mu=\lambda_4$. In physical variables, this kink solution corresponds to shock waves similar to those considered in Ref. [15] for the Choi-Camassa model [5].

The solutions given by Eqs. (4.12)–(4.17) look very similar to solutions found previously for other completely integrable equations, as, for example, the Kaup-Boussinesq equation [7, 11, 12] whose Lax pair (3.15), (3.16) is related to the Lax pair (3.7)–(3.10). However, in our case the relationships (4.8) between the physical variables $v,w$ and the function $\mu=\mu(\xi)$ are much more complicated than similar relationships for the Kaup-Boussinesq system. In particular, to get the real positive values of $v^2,w^2$, we must satisfy the condition that

(4.18)\begin{equation} b^2(\mu)-4a^2(\mu)\geq0\qquad\text{for}\qquad \lambda_2\leq\mu\leq\lambda_3. \end{equation}

This condition defines some region in four-dimensional space (4.10) of the parameters $\lambda_1,\lambda_2,\lambda_3,\lambda_4$. The form of this region is quite complicated, so we will confine ourselves to studying this condition for the soliton solutions (4.14). The case of the solution (4.16) can be considered in a similar way.

So, we consider the solution (4.14) with $\lambda_3=\lambda_4$ and demand that $b^2(\mu)-4a^2(\mu)$ does not vanish inside the interval

(4.19)\begin{equation} \lambda_2 \lt \mu \lt \lambda_4. \end{equation}

Let us find the locations of the zeroes of the equation $b^2(\mu)-4a^2(\mu)=0$. Obviously, it splits into two equations $b(\mu)\pm2a(\mu)=0$, so we define the quadratic functions of $\mu$:

(4.20)\begin{equation} \begin{aligned} & f_1(\mu)=b(\mu)+2a(\mu)=3\mu^2-2(s_1+1)\mu+1+s_2+s_1,\\ & f_2(\mu)=b(\mu)-2a(\mu)=3\mu^2-2(s_1-1)\mu+1+s_2-s_1, \end{aligned} \end{equation}

where

(4.21)\begin{equation} s_1=\lambda_1+\lambda_2+2\lambda_4,\qquad s_2=\lambda_1\lambda_2+2(\lambda_1+\lambda_2)\lambda_4+\lambda_4^2. \end{equation}

At the end points of the interval (4.19) the function $f_1(\mu)$ takes the values

(4.22)\begin{equation} f_1(\lambda_2)=(1+\lambda_4)^2-(\lambda_2-\lambda_1)(1-\lambda_2+2\lambda_4),\qquad f_1(\lambda_4)=(1+\lambda_1)(1+\lambda_2)\geq0. \end{equation}

The expression for $f_1(\lambda_2)$, considered as a quadratic function of $1+\lambda_4$, has a negative discriminant $-(\lambda_2-\lambda_1)(1+\lambda_1) \lt 0$, and, since

(4.23)\begin{equation} \left.f_1(\lambda_2)\right|_{\lambda_4=\lambda_2}=(1+\lambda_1)(1+\lambda_2)\geq0, \end{equation}

the function $f_1(\mu)$ has the same signs at both sides of the interval (4.19). Therefore, it can have zeroes inside this interval only if its negative minimum is located within (4.19). The function $f_1(\mu)$ has a minimum at

(4.24)\begin{equation} \mu_m^{(1)}=\frac13(\lambda_1+\lambda_2+2\lambda_4+1), \end{equation}

equal to

(4.25)\begin{equation} f_1(\mu_m^{(1)})=(1+\lambda_1)(1+\lambda_2)+2(1+\lambda_1+\lambda_2)\lambda_4+\lambda_4^2-\frac13(1+\lambda_1+\lambda_2+2\lambda_4)^2. \end{equation}

We have $\mu_m^{(1)}-\lambda_2=[2(\lambda_4-\lambda_2)+(1+\lambda_1)]/3 \gt 0$ and $\lambda_4-\mu_m^{(1)}=[\lambda_4-(1+\lambda_1+\lambda_2)]/3$, is negative for

(4.26)\begin{equation} \lambda_4 \lt 1+\lambda_1+\lambda_2, \end{equation}

that is $\mu_m^{(1)}$ is located outside the interval (4.19), so $f_1(\mu)$ does not vanish inside it in this case. But if $1+\lambda_1+\lambda_2 \lt \lambda_4$, then $\mu_m^{(1)}$ belongs to this interval, and it is easy to show that $f_1(\mu_m^{(1)}) \gt 0$ for

(4.27)\begin{equation} 1+\lambda_1+\lambda_2 \lt \lambda_4 \lt 1+\lambda_1+\lambda_2+\sqrt{3(1+\lambda_1)(1+\lambda_2)}. \end{equation}

Combining the inequalities (4.25) and (4.27), we arrive at the condition

(4.28)\begin{equation} \lambda_4 \lt 1+\lambda_1+\lambda_2+\sqrt{3(1+\lambda_1)(1+\lambda_2)}, \end{equation}

when the function $f_1(\mu)$ does not vanish inside the interval (4.19). Naturally, this condition is satisfied for all values of $\lambda_4$, $\lambda_2 \lt \lambda_4 \lt 1$, if the right-hand side of this inequality is greater than unity.

In a similar way, we find the values of $f_2(\mu)$ at the end points of the interval (4.19):

(4.29)\begin{equation} f_2(\lambda_2)=(1-\lambda_4)^2+(\lambda_2-\lambda_1)(1+\lambda_2-2\lambda_4),\qquad f_2(\lambda_4)=(1-\lambda_1)(1-\lambda_2)\geq0. \end{equation}

In this case, $f_2(\lambda_2)$ can be either positive, or negative, depending on the values of $\lambda_4$. We have

(4.30)\begin{equation} \left.f_2(\lambda_2)\right|_{\lambda_4=1}=-(\lambda_2-\lambda_1)(1-\lambda_2) \lt 0, \end{equation}

and $f_2(\lambda_2)$, considered as a function of $\lambda_4$, becomes positive for

(4.31)\begin{equation} \lambda_4 \lt 1+\lambda_2-\lambda_1-\sqrt{(\lambda_2-\lambda_1)(1-\lambda_1)}, \end{equation}

that is for these values of $\lambda_4$ the function $f_2(\mu)$ has the same signs at both sides of the interval (4.19). It is easy to find that $f_2(\mu)$ takes its minimal value at

(4.32)\begin{equation} \mu_m^{(2)}=\frac13(\lambda_1+\lambda_2+2\lambda_4-1), \end{equation}

and this minimal value is equal to

(4.33)\begin{equation} f_2(\mu_m^{(2)})=(1-\lambda_1)(1-\lambda_2)-2(1-\lambda_1-\lambda_2)\lambda_4+\lambda_4^2-\frac13(1-\lambda_1-\lambda_2-2\lambda_4)^2. \end{equation}

Now we have two possibilities: (i) the point $\mu_m^{(2)}$ belongs to the interval (4.19) and then we have to find the condition that $f_2(\mu_m^{(2)})\geq0$; (ii) the point $\mu_m^{(2)}$ is located outside the interval (4.19), and then the value of $f_2(\mu_m^{(2)})$ is irrelevant, provided the condition (4.31) is fulfilled.

We have $\mu_m^{(2)}-\lambda_2=[2\lambda_4-(1-\lambda_1+2\lambda_2)]/3$ and $\lambda_4-\mu_m^{(2)}=[(\lambda_4-\lambda_2)+(1-\lambda_1)]/3 \gt 0$. Consequently, $\mu_m^{(2)}$ belongs to the interval for $\lambda_4 \gt \lambda_2+(1-\lambda_1)/2$. Combining this inequality with another condition $f_2(\mu_m^{(2)})\geq0$, we arrive for the case (i) at the condition

(4.34)\begin{equation} \lambda_2+(1-\lambda_1)/2 \lt \lambda_4 \lt -1+\lambda_1+\lambda_2+\sqrt{3(1-\lambda_1)(1-\lambda_2)}. \end{equation}

Such a region of $\lambda_4$ exists for

(4.35)\begin{equation} \lambda_1 \lt \lambda_2 \lt \frac14(1+3\lambda_1). \end{equation}

If $\lambda_4 \lt \lambda_2+(1-\lambda_1)/2$, then the point $\mu_m^{(2)}$ is located outside the interval (4.19). In this case (ii) we combine this condition with inequality (4.31) and arrive at inequalities

(4.36)\begin{equation} \lambda_4 \lt 1+\lambda_2-\lambda_1-\sqrt{(\lambda_2-\lambda_1)(1-\lambda_1)} \lt \lambda_2+(1-\lambda_1)/2, \end{equation}

which define the region with positive values of $f_2(\mu)$ in the interval (4.19). Such a region of $\lambda_4$ exists for

(4.37)\begin{equation} \frac14(1+3\lambda_1) \lt \lambda_2. \end{equation}

This completes the analysis of soliton solutions (4.14) with single-valued dependence of $v^2$ and $w^2$ on $\xi$.

A typical soliton solution is shown in Fig. 1, where we have chosen $\lambda_1=-0.2,\lambda_2=-0.1$ and $\lambda_4=0.689$. The condition (4.28) is obviously fulfilled and the value of $\lambda_4$ is taken close to the upper admissible limit $-1+\lambda_1+\lambda_2+\sqrt{3(1-\lambda_1)(1-\lambda_2)}\approx0.68998$ of the condition (4.34). Above this bifurcation point the topology of curves changes and they reconnect into three separate two-valued plots defined on three disconnected segments of the $\xi$-axis. Similar results can be obtained for the soliton solution (4.16).

Figure 1. Plots of the curves (4.8) for the soliton solution (4.14) with $\lambda_1=-0.2,\lambda_2=-0.1,\lambda_4=0.689$. Either of the curves can be taken as a plot of $v^2$, then the other curve represents the plot of $w^2$.

Now let us turn to the “kink” solution (4.17). In this case with $\lambda_2=\lambda_1,\lambda_3=\lambda_4$ we have

(4.38)\begin{equation} s_1=2(\lambda_1+\lambda_4),\qquad s_2=\lambda_1^2+\lambda_4^2+4\lambda_1\lambda_4, \end{equation}

and we get

(4.39)\begin{equation} \begin{aligned} & f_1(\mu)=3\mu^2-2[2(\lambda_1+\lambda_4)+1]\mu+1+\lambda_1^2+\lambda_4^2+4\lambda_1\lambda_4+2(\lambda_1+\lambda_4),\\ & f_2(\mu)=3\mu^2-2[2(\lambda_1+\lambda_4)-1]\mu+1+\lambda_1^2+\lambda_4^2+4\lambda_1\lambda_4-2(\lambda_1+\lambda_4). \end{aligned} \end{equation}

These two functions are simultaneously positive for $\lambda_4 \gt \lambda_1$ in the shadowed region shown in Fig. 2, where the separating straight lines are given by the expressions

(4.40)\begin{equation} \lambda_4=1\pm\sqrt{3}+(2\pm\sqrt{3})\lambda_1\quad\text{and}\quad \lambda_4=-1\pm\sqrt{3}+(2\mp\sqrt{3})\lambda_1. \end{equation}

Figure 2. The shadowed region corresponds to the admissible values of the parameters $\lambda_4 \gt \lambda_1$ for “kink” solutions (4.17). The point $A$ has coordinates $A=(-1/\sqrt{3},1/\sqrt{3})$.

They cross each other at the point with coordinates $A=(-1/\sqrt{3},1/\sqrt{3})$. Typical plots of $v^2,w^2$ are shown in Fig. 3 for $\lambda_1=-0.4,\lambda_4=0.6$. If we assume, for example, that the upper curve corresponds to $v^2$ and the lower curve to $w^2$, then this solution describes a specific dispersive shock wave joining the flows:

(4.41)\begin{equation} \begin{aligned} & v^2\to1,\quad w^2\to\lambda_4^2\quad\text{for}\quad \xi\to-\infty,\\ & v^2\to1,\quad w^2\to\lambda_1^2\quad\text{for}\quad \xi\to+\infty. \end{aligned} \end{equation}

Figure 3. Plots of the curves (4.8) for the “kink” solution (4.17) with $\lambda_1=-0.4,\lambda_4=0.6$. Either of the curves can be taken as a plot of $v^2$, then the other curve represents the plot of $w^2$.

This shock corresponds to a jump of the relative thickness of the layers propagating with velocity $V=(\lambda_1+\lambda_4)/2$. This formula agrees with the Rankine-Hugoniot condition which follows from the first Eq. (2.18). Denoting jumps by square brackets, in case $v=1$ at both sides of the shock, we obtain $-V[w]+\frac12[w^2]=0$ and, consequently,

(4.42)\begin{equation} V=\frac{[w^2]}{2[w]}=\frac12(\lambda_1+\lambda_4). \end{equation}

In the opposite case, when we interpret the upper curve in Fig. 3 as a plot of $w^2$ and the lower curve as a plot of $v^2$, we obtain a shock which joints two uniform flows of a lighter fluid with different velocities through a “puddle” of a heavier fluid in between. It is interesting that in case of $\lambda_4=-\lambda_1$ we have $V=0$, that is this solution describes a stationary “puddle” of a heavier fluid whose profile depends on the value of the relative velocity $v=|\lambda_1|$ at infinity.

Investigation of stability of these solutions is outside the scope of this paper.

5. Conclusion

In this paper, we have studied the completely integrable generalisation of Ovsyannikov model of two-layer shallow water hydrodynamics. This generalisation is based on the fact that the large-scale Ovsyannikov hydrodynamics is compatible with propagation of short-wavelength wave packets with a certain dispersion relation in the sense of asymptotic integrability [26], that is the Hamilton equations which describe the packet’s propagation have an integral of motion. It turns out that this integral can be considered as a quasiclassical approximation to the equation of some Lax pair [27, 29]. Then the other equation of this Lax pair is related to the dispersion relation, so the transition from the quasiclassical approximation to the appropriate wave theory (‘quantisation’) yields the Lax pair for the desired generalisation [26, 27]. Equations obtained in this way can be solved by standard methods of the theory of completely integrable nonlinear wave equations. However, in case of our model the relationships between the physical variables and the most natural mathematical variables of the finite-gap integration method are quite complicated and they only provide physical real values for some sets of parameters which characterize the solutions. We have described these sets for the soliton and kink solutions and illustrated the theory by concrete examples. Although Ovsyannikov’s scaling approach excludes consideration of some features of two-fluid systems, which depend on the ratio of fluids’ densities, nevertheless our model describes qualitatively some essential properties of real two-layer shallow water systems.