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. [Reference Esler and Pearce15, Reference Ivanov, Kamchatnov, Congy and Pavloff21, Reference Kamchatnov24]).

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. [Reference Kamchatnov26, Reference Kamchatnov27]. 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. [Reference Kamchatnov26, Reference Kamchatnov27]), 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 [Reference Congy, Kamchatnov and Pavloff8, Reference Ivanov, Kamchatnov, Congy and Pavloff21], 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 [Reference Kamchatnov25, Reference Qu, Pitaevskii and Stringari41]. 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. [Reference Barros and Gavrilyuk1] 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 [Reference Kamchatnov26, Reference Kamchatnov27]

(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. [Reference Kamchatnov23, Reference Kamchatnov26, Reference Kamchatnov27], 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. [Reference Ovsyannikov39], 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 [Reference Kaup30]. 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. [Reference Congy, Ivanov, Kamchatnov and Pavloff7, Reference El, Grimshaw and Kamchatnov11, Reference El, Grimshaw and Pavlov12], 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., [Reference Kamchatnov24]).

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. [Reference Barros and Gavrilyuk1] 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. [Reference Esler and Pearce15] for the Choi-Camassa model [Reference Choi and Camassa5].

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 [Reference Congy, Ivanov, Kamchatnov and Pavloff7, Reference El, Grimshaw and Kamchatnov11, Reference El, Grimshaw and Pavlov12] 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 [Reference Kamchatnov26], 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 [Reference Kamchatnov27, Reference Kamchatnov and Shaykin29]. 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 [Reference Kamchatnov26, Reference Kamchatnov27]. 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.