1. Introduction
The Korteweg–de Vries- and Ostrovsky-type family of models plays an important role in understanding the behaviour of long nonlinear internal waves commonly observed in coastal seas, narrow straits and river–sea interaction areas (see Grimshaw et al. Reference Grimshaw, Ostrovsky, Shrira and Stepanyants1998; Lamb Reference Lamb2005; Helfrich & Melville Reference Helfrich and Melville2006; Bona, Lannes & Saut Reference Bona, Lannes and Saut2008; Ostrovsky et al. Reference Ostrovsky, Pelinovsky, Shrira and Stepanyants2015; Khusnutdinova & Zhang Reference Khusnutdinova and Zhang2016; Stastna Reference Stastna2022; Ostrovsky et al. Reference Ostrovsky, Pelinovsky, Shrira and Stepanyants2024, and references therein). The Ostrovsky equation (Ostrovsky Reference Ostrovsky1978; see also Grimshaw Reference Grimshaw1985),
\begin{equation} (A_T + \nu A A_\xi + \lambda A_{\xi \xi \xi })_\xi = \gamma A ,\end{equation}
with constant coefficients
$\nu$
,
$\lambda$
and
$\gamma$
, is a rotationally modified extension of the integrable Korteweg–de Vries (KdV) equation (Boussinesq Reference Boussinesq1871; Korteweg & de Vries Reference Korteweg and de Vries1895; Gardner et al. Reference Gardner, Greene, Kruskal and Miura1967), accounting for the leading-order balance of weak nonlinear, dispersive and rotational effects. In the general setting of a density stratified fluid described by the Euler equations on an
$f$
-plane, where
$f$
stands for the Coriolis parameter, with boundary conditions appropriate for oceanic applications, (1.1) is written for the amplitude
$A(\xi , T)$
of a single plane internal mode
$\phi (z)$
in a reference frame moving with the linear long-wave speed
$c_0$
. Here,
$\xi$
and
$T$
denote the characteristic and slow-time variables, respectively. In physical variables, where the bottom is at
$z = - h$
and unperturbed surface is at
$z=0$
, the modal equations have the form
\begin{align} &\displaystyle (\rho _0 W^2 \phi _{z} )_z + \rho _0 N^2 \phi = 0 , \end{align}
\begin{align} \displaystyle \phi = 0 \quad \text{at} \quad & z =-h, \quad \hbox{and} \quad W^2 \phi _{z }= g\phi \quad \text{at} \quad z =0 . \end{align}
Here,
$\rho _{0} (z)$
is the stable background density profile,
$g$
is gravity,
$N^2 = -g\rho _{0z} / \rho _0$
,
$W=c_0-u_{0}$
, where
$u_{0}(z)$
is the background shear flow supported by a body force, and it is assumed that there are no critical levels, that is
$W\ne 0$
for any
$z$
in the flow domain. The nonlinear, dispersive and rotational coefficients
$\nu , \lambda$
and
$\gamma$
, respectively, are given by
\begin{eqnarray} I\nu = 3\, \int _{-h}^0\rho _0 W^2 \phi _z^3\, {\rm d}z , \quad I\lambda = \int _{-h}^0\rho _0 W^2 \phi ^2\, {\rm d}z , \quad I\gamma = f^2 \int ^{0}_{-h} \, \rho _{0} \varPhi \phi _{z} \, {\rm d}z , \end{eqnarray}
where
\begin{eqnarray} I = 2\, \int _{-h}^0 \rho _0 W \phi _z^2 \, {\rm d}z , \quad \rho _0 W\varPhi = \rho _0 W\phi _{z} -(\rho _{0}u_{0})_z \phi , \end{eqnarray}
and
$f$
is the Coriolis parameter (a single mode reduction of the bi-modal system derived in Alias, Grimshaw & Khusnutdinova Reference Alias, Grimshaw and Khusnutdinova2013). Note that when there is no shear flow, that is
$u_0 (z) \equiv 0$
, then
$\varPhi \equiv \phi _z$
and
$\gamma = f^2/2c$
; in this case
$\lambda \gamma \gt 0$
(normal dispersion), but sufficiently strong shear near a pycnocline may lead to a situation where
$\lambda \gamma \lt 0$
(anomalous dispersion, see Alias, Grimshaw & Khusnutdinova Reference Alias, Grimshaw and Khusnutdinova2014).
The Ostrovsky equation became a paradigm forming model for studying the effects of rotation on the evolution of internal waves with the natural initial conditions in the form of KdV solitons (see Grimshaw & Helfrich Reference Grimshaw and Helfrich2008; Grimshaw et al. Reference Grimshaw R.H., Helfrich and Johnson2013; Ostrovsky et al. Reference Ostrovsky, Pelinovsky, Shrira and Stepanyants2015; Stepanyants Reference Stepanyants2020, and references therein) and cnoidal waves, with an emphasis on modulational instability (see Whitfield & Johnson Reference Whitfield and Johnson2014, Reference Whitfield and Johnson2017; Johnson Reference Johnson, Oregano and Perkins2025, and references therein), as well as the related qualitative analysis of the long-time asymptotics of strongly interacting internal modes described by solutions of coupled Ostrovsky equations (Alias et al. Reference Alias, Grimshaw and Khusnutdinova2014). One of the aims of our current study is to extend the modelling to situations when the initial conditions are close to a cnoidal wave (approximately a chain of KdV solitons), but they are not perfectly periodic, and rather have some local defects within the computational domain (hence, they can be viewed as being quasi-periodic), and to show, for the case of normal dispersion
$\lambda \gamma \gt 0$
, that this class of initial conditions also leads to the generation of slower-moving bursts, qualitatively similar to the previously modelled unsteady wavepackets (see Grimshaw & Helfrich Reference Grimshaw and Helfrich2008; Whitfield & Johnson Reference Whitfield and Johnson2017; Stepanyants Reference Stepanyants2020, and references therein). The motivation for that study stems from the observational data of the type shown in figure 1, where we can see formation of a wavetrain of internal solitary waves close to an imperfect cnoidal wave. Indeed, given that the waves generally propagate in a variable environment, quasi-periodic initial conditions seem to be a more natural choice than a single (pure or only slowly modulated) cnoidal wave. Recently, there appeared renewed interest in the possibility of generating rogue waves and breathers by various localised perturbations of cnoidal waves both in integrable and non-integrable settings (see Kuznetsov & Mikhailov Reference Kuznetsov and Mikhailov1975; Onorato et al. Reference Onorato, Residori, Bortolozzo, Montina and Arecchi2013; Chabchoub et al. Reference Chabchoub, Hoffman, Onorato, Genty, Dudley and Akhmediev2013; Kedziora, Ankiewicz & Akhmediev Reference Kedziora, Ankiewicz and Akhmediev2014; Maiden & Hoefer Reference Maiden and Hoefer2016; Chabchoub et al. Reference Chabchoub, Slunyaev, Hoffman, Genty, Dias, Kibler, Dudley and Akhmediev2021; He et al. Reference He, Witt, Trillo, Chabchoub and Hoffmann2022; Bertola, Jenkins & Tovbis Reference Bertola, Jenkins and Tovbis2023; Grava et al. Reference Grava, McLaughlin, Girotti, Jenkins and Minakov2023; Hoefer, Mucalica & Pelinovsky Reference Hoefer, Mucalica and Pelinovsky2023; Mao et al. Reference Mao, Chandramouli, Xu and Hoefer2023; Chandramuli et al. Reference Chandramouli, Mao and Hoefer2024, and references therein), among other possible mechanisms (e.g. Kharif, Pelinovsky & Slunyaev Reference Kharif, Pelinovsky and Slunyaev2009; Zakharov Reference Zakharov2009; Pelinovsky & Slunyaev Reference Pelinovsky and Slunyaev2016; Bokaeeyan, Ankiewicz & Akhmediev Reference Bokaeeyan, Ankiewicz and Akhmediev2019; Choi et al. Reference Choi, Bokhove, Kalogirou and Kelmanson2022; Slunyaev & Shrira Reference Slunyaev and Shrira2023; Agafontsev et al. Reference Agafontsev, Gelash, Randoux and Suret2024; Congy et al. Reference Congy, El, Roberti, Tovbis, Randoux and Suret2024; Flamarion, Pelinovsky & Didenkulova Reference Flamarion, Pelinovsky and Didenkulova2024; Slunyaev Reference Slunyaev2024, and references therein). Moreover, internal rogue waves registered in the oceans have been linked to the KdV solitons (Osborne Reference Osborne2010). Our study extends the line of research related to localised perturbations by adding the effect of rotation, as well as discussing the generalised (shock-like) travelling waves in the absence of rotation (for the KdV equation).
Figure 1. The first
$1.7$
hours of the colour contour time series of temperature profiles off Northern Oregon from the surface to 35 m depth. The figure is adapted from Stanton & Ostrovsky (Reference Stanton and Ostrovsky1998).
Unlike the KdV equation, the Ostrovsky equation has a constraint on the mean value of its regular solutions. For example, for periodic solutions on the interval
$[-L, L]$
, any regular solution should have zero mean
\begin{equation} \int _{-L}^L A\ {\rm d}\xi = 0. \end{equation}
Additional constraints appear when solutions are considered in the class of functions vanishing at infinity (Benilov Reference Benilov1992). The existing derivations of the Ostrovsky equation from the Euler equations do not allow one to consider the arbitrary initial conditions for the field variables of the parent system, but only those which agree with this constraint. This restriction on the choice of initial conditions of the Cauchy problem for the parent system became known as the ‘zero-mean contradiction’. Hence, another aim of our study is to generalise the construction of weakly nonlinear solutions leading to the Ostrovsky equation in such a way that the zero-mean contradiction is avoided. Indeed, this has been previously done within the scope of the derivation of Ostrovsky-type models in the simpler settings of the Boussinesq–Klein–Gordon and coupled Boussinesq equations (Khusnutdinova, Moore & Pelinovsky Reference Khusnutdinova, Moore and Pelinovsky2014; Khusnutdinova & Tranter Reference Khusnutdinova and Tranter2019, Reference Khusnutdinova and Tranter2022). In contrast to the previous work, in the fluids context the mean-field equations appear to be coupled to equations for deviations from the mean values, presenting a new challenge, which is addressed in our present paper.
Hence, the aim of our study is twofold, and the rest of the paper is organised as follows. In § 2 we introduce Helfrich’s rotation-modified two-layer Miyata–Maltseva–Choi–Camassa (MMCC-f) model (Helfrich Reference Helfrich2007) and obtain its simpler weakly nonlinear reduction which is used to develop the subsequent (weakly nonlinear) derivations of the reduced model (for the original MMCC model see Miyata Reference Miyata1988; Maltseva Reference Maltseva1989; Choi & Camassa 1996, Reference Choi and Camassa1999). While this setting is simpler than the full Euler equations, our derivations reveal that it retains the key complexity: the mean-field equations are coupled with the equations for the deviations. Therefore, this model provides an appropriate framework for our developments. Next, in § 3 we refine the derivation of the uni-directional model by considering the simultaneous evolution of mean fields and deviations from them. This approach introduces the fast characteristic variable
$\xi = x - c_0 t$
and two slow-time variables
$\tau = \sqrt {\alpha }t$
and
$T = \alpha t$
, rather than just one, where
$\alpha$
is the small-amplitude parameter (we assume the maximum balance conditions for the weak nonlinearity, dispersion and rotation). We find a way to bypass the zero-mean contradiction by developing simultaneous asymptotic expansions of both the mean-field variables and deviations from them. This allows us to construct a more general class of solutions allowing the initial conditions for all fluid variables to have arbitrary (and generally time-dependent) mean values, while the emerging Ostrovsky equations have zero mean by construction. In § 4 we use the constructed weakly nonlinear solution combined with numerical modelling using the KdV and Ostrovsky equations in order to investigate the effect of rotation on the evolution of cnoidal waves of the KdV equation close to their solitonic limit and having local amplitude and/or periodicity defects. We begin the section by considering the effect of rotation on simple soliton and cnoidal wave solutions, as well as dark and bright breathers of the KdV equation and expansion and contraction periodicity defects in order to set up the phenomenological framework for the discussion of our main numerical results concerning generic localised perturbations. The expansion/contraction defects are introduced by cutting the cnoidal wave at the trough and symmetrically inserting a piece of a straight line/extracting a small symmetric piece around the trough and gluing the remaining parts together, respectively, see the first two rows of figure 17 in Appendix A. In all cases we take a sufficiently large computational domain and impose periodic boundary conditions to model the resulting quasi-periodic solutions. We show that, combined with the effect of rotation, initial conditions in the form of cnoidal waves with local defects can lead to bursts of large-amplitude internal waves and shear currents at relatively small slow-time values and tens of times faster than the generation of wavepackets from a soliton and a slowly modulated cnoidal wave initial conditions considered in earlier studies (Grimshaw & Helfrich Reference Grimshaw and Helfrich2008; Whitfield & Johnson Reference Whitfield and Johnson2017). This significantly increases the likelihood of emergence of such rotation-induced bursts under real oceanic conditions (see the relevant discussions and oceanic synthetic-aperture radar images in Stastna et al. Reference Stastna, Poulin, Rowe and Subich2009; Whitfield & Johnson Reference Whitfield and Johnson2017; Grimshaw, da Silva & J. Reference Grimshaw, da Silva and Magalhaes2017). For the generic local defects introduced by adding a localised perturbation, this can be attributed to the formation of a pair of bright and dark KdV breathers and contraction and expansion periodicity defects, with the subsequent effect of rotation. For the pure expansion and contraction defects the bursts do not take place without rotation. In fact, in pseudospectral numerical runs the smoothed counterparts of such initial conditions evolve almost like travelling wave solutions of the KdV equation, for a very long time. We argue that this happens because these functions (and some other functions constructed from the known KdV solutions by similar procedures) satisfy all (infinitely many) conservation laws of the KdV equation (see Appendix A for the details of this and related discussions). Moreover, a cnoidal wave with an expansion defect has a continuous first derivative, satisfying the Weirstrass–Erdmann corner condition for broken (non-smooth) extremals (e.g. Fox Reference Fox1954) and it is a generalised (‘shock-like’) travelling wave, using the terminology introduced by Gavrilyuk & Shyue (Reference Gavrilyuk and Shyue2022). We also make remarks about a possible weak formulation for the generalised travelling waves of the KdV equation. The Cauchy problem for the KdV equation with periodic boundary conditions in
$L^2$
(the square-integrable functions on a torus) is globally well posed, including uniqueness and continuous dependence with respect to the initial data (Bourgain Reference Bourgain1993). The bursts observed in our modelling with pure periodicity defects are then attributed solely to the effect of rotation. Moreover, we show that the effects discussed in the paper are structurally stable with respect to the natural (compatible with the period of the background cnoidal wave) variations in the size of the computational domain. We finish with a discussion in § 5, where we also show an example where a local perturbation of a cnoidal wave of the type considered in our paper has led to the generation of a rogue wave. A pseudospectral scheme used to solve the Ostrovsky equation is discussed in Appendix B.
2. The MMCC-
$\boldsymbol{f}$
model
We consider the
$f$
-plane extension of the MMCC model for an inviscid, incompressible two-layer fluid with the rigid lid introduced by Helfrich (Reference Helfrich2007)
\begin{align} h_{it} + (h_i \overline {u}_i)_x =0, \end{align}
where
\begin{align} \overline {u}_{it} + \overline {u}_i \overline {u}_{\textit{ix}} -f\overline {v}_i=-g\eta _x+\dfrac {1}{\rho _i}P_x + D_i, \end{align}
\begin{align} D_i = h_i^{-1}\left \{ \dfrac {h_i^3}{3}\big [ \overline {u}_{\textit{ixt}} + \overline {u}_i \overline {u}_{\textit{ixx}} -(\overline {u}_{\textit{ix}})^2 \big ]\right \}_x \!, \end{align}
\begin{align} \overline {v}_{it} + \overline {u}_i \overline {v}_{\textit{ix}} + f\overline {u}_{i} = 0. \end{align}
Here,
$h_i$
,
$\rho _i$
are the layer depths and densities,
$\overline {u}_i, \overline {v}_i$
denote the depth-averaged, over each layer, horizontal velocities in the
$x$
and
$y$
directions with
$i=1$
and
$2$
referring to the upper and lower layers, respectively;
$f$
is the Coriolis parameter,
$g$
is gravity and
$P$
is the pressure at the interface. The subscripts
$t$
and
$x$
denote partial derivatives. In the absence of motion
$h_1=h_0$
, and
$h_1+h_2=H$
(total depth) and
$\eta$
is the interfacial displacement (
$\eta = h_0 - h_1$
, see figure 2).
Figure 2. Schematic of a two-layer fluid with rotation in the rigid-lid approximation.
In the Boussinesq approximation
$((\rho_2 - \rho_1)/\rho_1 \ll 1)$
, the equations can be simplified by eliminating
$P_x$
. Using
$\sqrt {g'H},\, H,\, l$
(a typical wavelength) and
${l}/{\sqrt {g'H}}$
to non-dimensionalise
$(u_i,\,v_i),\, h_i,\, x$
and
$t$
, respectively, and changing the variables to
$s=\overline {u}_2-\overline {u}_1$
,
$v=\overline {v}_2-\overline {v}_1$
,
$h_2=1-h_1$
,
$U = \overline{u}_1 h_1 + \overline{u}_2 h_2$
and
$V = \overline {v}_1h_1 + \overline {v}_2h_2$
with the barotropic transport in
$x$
direction
$U = F(t) \neq 0$
in general, we obtain
\begin{align} &\qquad\,\, \eta _{t} + (c_0^2 s +\sigma s\eta - s\eta ^2 + F\eta )_x = 0, \end{align}
\begin{align}& s_t + \left(\dfrac {\sigma }{2}s^2 -s^2\eta +\eta +Fs\right)_x - \tilde \gamma v = \tilde \beta (D_2-D_1), \end{align}
\begin{align}&\,\, v_t { -} sv\eta _{x} + \sigma sv_x -2 s\eta v_x + Fv_x + sV_x + \tilde \gamma s = 0, \end{align}
\begin{align}&\quad V_t + ( c_0^2 sv + \sigma sv\eta - sv\eta ^2 + FV )_x +\tilde \gamma F = 0, \\[9pt] \nonumber \end{align}
where
${ c_0^2 = h_0 (1-h_0)}, \sigma = 2 h_0 - 1, g' = g \Delta \rho /\rho _1, L_R = \sqrt {g'H}/f, \tilde \gamma = {l}/{L_R}, \tilde \beta = ({H}/{l} )^2$
.
Next, we consider the Ostrovsky equation regime and assume that
$\tilde \beta =\mathcal{O}(\alpha )$
and
$\tilde \gamma =\mathcal{O}(\alpha ^{1/2})$
, where
$\alpha$
is the small-amplitude parameter, and scale
$\tilde \beta = \alpha \beta ,\, \tilde \gamma = \sqrt {\alpha } \gamma$
, where
$\beta , \gamma = \mathcal{O}(1)$
. Then, considering the asymptotic expansions
\begin{align} (\eta , s, v, V) &= \alpha (\eta , s, v, V)^{(1)} + \alpha ^2 (\eta , s, v, V)^{(2)} + O(\alpha ^3), \end{align}
\begin{align} F &= \alpha F^{(1)} + \alpha ^{2} F^{(2)} + O(\alpha ^3), \end{align}
and dropping the terms of
$O(\alpha ^3)$
(not used in our subsequent derivations), gives the simpler weakly nonlinear equations (all variables have been scaled by
$\alpha$
, but for simplicity we keep the same notations for the new
$O(1)$
variables)
\begin{align} \eta _{t} + c_0^2 s_x = - \alpha (\sigma s\eta + F\eta )_x, \end{align}
\begin{align} s_t + \eta _x = \sqrt {\alpha } \gamma v - \alpha \left(\sigma {\textit{ss}}_x + {\textit{Fs}}_x - \dfrac { \beta c_0^2}{3} s_{\textit{xxt}}\right), \end{align}
\begin{align} v_t = - \sqrt {\alpha } \gamma s - \alpha (\sigma {sv}_x + {Fv}_x + {sV}_x), \end{align}
\begin{align} V_t = - \sqrt {\alpha } \gamma F - \alpha \big(c_0^2 {sv} + {FV} \big)_x. \end{align}
For
$F = 0$
, the equations reduce to the model derived by Gerkema (Reference Gerkema1996), Helfrich (Reference Helfrich2007).
We consider the periodic solutions on the interval
$[-L, L]$
(a typical setting for numerical runs using pseudospectral methods). Alongside (2.11)–(2.14) we will consider the equations for the evolving mean fields by averaging the system with respect to
$x$
over this interval and denoting the mean values by hats
\begin{align} \hat {\eta }_t = 0, \end{align}
\begin{align} \hat {s}_t = \sqrt {\alpha } \gamma \hat {v}, \end{align}
\begin{align} \hat {v}_t = - \sqrt {\alpha } \gamma \hat {s} - \alpha \frac {1}{2L} \int _{-L}^{L} (\sigma s v_x +sV_x)\, {\rm d}x, \end{align}
\begin{align} \hat {V}_t = -\sqrt {\alpha } \gamma F. \end{align}
Here, the mean fields are generally time-dependent. Hence, the equations for the evolving mean values and deviations from the mean values are generally coupled because of the integral term in (2.17). This differs from our previous derivations of Ostrovsky-type equations free from zero-mean contradiction within the scope of the Boussinesq–Klein–Gordon equation (Khusnutdinova et al. Reference Khusnutdinova, Moore and Pelinovsky2014; Khusnutdinova & Tranter Reference Khusnutdinova and Tranter2019) and coupled Boussinesq equations (Khusnutdinova & Tranter Reference Khusnutdinova and Tranter2022), where the mean fields were described by independent equations.
3. Weakly nonlinear uni-directional solution free from zero-mean contradiction
In this section, we construct a large class of solutions of the system (2.11)–(2.18) describing uni-directionally propagating waves on top of evolving mean fields by presenting each variable in the form of the sum of its mean value (generally, time-dependent) and deviation from this mean value
\begin{equation} \eta = \hat \eta + \tilde \eta , \ s = \hat s + \tilde s, \ v = \hat v + \tilde v, \ V = \hat V + \tilde V, \end{equation}
where the unknowns
$\hat \eta , \hat s, \hat v, \hat V$
and the given function
$F$
are assumed to be functions of
$t, \tau = \sqrt \alpha t, T = \alpha t$
, while the deviations
$\tilde \eta , \tilde s, \tilde v, \tilde V$
are assumed to depend on the fast characteristic variable
$\xi = x - c_0 t$
, where
$c_0$
is the linear long-wave speed, and two slow-time variables
$\tau$
and
$T$
. Both the averages and deviations from the averages are sought in the form of asymptotic multiple-scale expansions in powers of
$\sqrt {\alpha }$
\begin{align} (\hat \eta , \hat s, \hat v, \hat V) &= (\hat \eta , \hat s, \hat v, \hat V)^{(0)} + \sqrt {\alpha } (\hat \eta , \hat s, \hat v, \hat V)^{(1)}+ \alpha (\hat \eta , \hat s, \hat v, \hat V)^{(2)} + O\big(\alpha ^{3/2}\big), \end{align}
\begin{align} (\tilde \eta , \tilde s, \tilde v, \tilde V) &= (\tilde \eta , \tilde s, \tilde v, \tilde V)^{(0)} + \sqrt {\alpha } (\tilde \eta , \tilde s, \tilde v, \tilde V)^{(1)} + \alpha (\tilde \eta , \tilde s, \tilde v, \tilde V)^{(2)} + O\big(\alpha ^{3/2}\big). \end{align}
The prescribed function
$F(t)$
defining the barotropic transport in the
$x$
direction is also assumed to be given in the form
\begin{equation} F = F^{(0)} + \sqrt {\alpha } F^{(1)} + \alpha F^{(2)} + O\big(\alpha ^{3/2}\big). \end{equation}
The average interfacial displacement
$\hat \eta$
is a constant, and if
$\hat \eta ^{(0)} \ne 0$
, then without loss of generality we assume that
$\hat \eta = \hat \eta ^{(0)}$
. If
$\hat \eta ^{(0)} = 0$
, but
$\hat \eta ^{(1)} \ne 0$
, then without loss of generality
$\hat \eta = \sqrt {\alpha } \hat \eta ^{(1)}$
, etc. Similarly, if
$F^{(0)} \ne 0$
, then we can assume
$F = F^{(0)}$
, etc.
At
$O(1)$
we obtain the equations
\begin{align} &\qquad\qquad\quad -c_0 \tilde \eta ^{(0)}_{\xi } + c_0^2 \tilde s^{(0)}_{\xi } =0, \quad \hat s^{(0)}_t - c_0 \tilde s^{(0)}_{\xi } + \tilde \eta ^{(0)}_{\xi } =0, \end{align}
\begin{align} &\hat v^{(0)}_t - c_0 \tilde v^{(0)}_\xi = 0, \quad \hat V^{(0)}_t - c_0 \tilde V^{(0)}_\xi = 0, \quad \hat s^{(0)}_t =0,\quad \hat v^{(0)}_t=0, \quad \hat V^{(0)}_t = 0, \end{align}
implying
\begin{align} & \hat s^{(0)} = A^{(0)}(\tau , T), \quad \hat v^{(0)} = B^{(0)}(\tau , T), \quad \hat V^{(0)} = C^{(0)}(\tau , T), \end{align}
\begin{align} &\qquad\qquad \displaystyle \tilde s^{(0)} = \frac {1}{c_0} \tilde \eta ^{(0)}, \quad \tilde v^{(0)} = 0, \quad \tilde V^{(0)} = 0. \end{align}
Here, the functions
$(A, B, C)^{(0)}$
are arbitrary functions of their variables, and we have used that all deviations should have zero mean, by construction of the solution.
Next, collecting the terms of
$O(\sqrt {\alpha })$
, using relations (3.7), (3.8) and simplifying the resulting system by virtue of its averaged members, we obtain the system
\begin{align} &\qquad\quad \displaystyle -c_0 \tilde \eta ^{(1)}_{\xi } + c_0^2 \tilde s^{(1)}_{\xi } = - \tilde \eta ^{(0)}_\tau , \quad - c_0 \tilde s^{(1)}_{\xi } + \tilde \eta ^{(1)}_{\xi } = - \frac {1}{c_0} \tilde \eta ^{(0)}_{\tau }, \end{align}
\begin{align} &\qquad\qquad\qquad\qquad\qquad \displaystyle \tilde v^{(1)}_\xi = \frac {\gamma }{c_0^2} \tilde \eta ^{(0)}, \quad \tilde V^{(1)}_\xi = 0, \end{align}
\begin{align} &\displaystyle \hat s^{(1)}_t = \gamma B^{(0)} - A^{(0)}_{\tau }, \quad \hat v^{(1)}_t= - \gamma A^{(0)} - B^{(0)}_{\tau }, \quad \hat V^{(1)}_t = - \gamma F^{(0)} - C^{(0)}_\tau . \end{align}
Equations (3.9) imply
\begin{equation} \tilde \eta ^{(0)}_{\tau } = 0, \quad \tilde s^{(1)} = \frac {1}{c_0} \tilde \eta ^{(1)}, \end{equation}
again using that all deviations have zero mean.
To avoid secular growth in
$\hat s^{(1)}$
and
$\hat v^{(1)}$
, we require that
\begin{equation} A^{(0)}_{\tau } = \gamma B^{(0)}, \quad B^{(0)}_{\tau } = - \gamma A^{(0)}, \end{equation}
which implies
$A^{(0)}_{\tau \tau } = - \gamma ^2 A^{(0)},$
yielding
\begin{equation} A^{(0)} = A(T) \cos \gamma [\tau + \phi (T)], \quad B^{(0)} = -A(T) \sin \gamma [\tau + \phi (T)], \end{equation}
where
$A(T)$
and
$\phi (T)$
are arbitrary functions of
$T$
.
To avoid secular growth in
$\hat V^{(1)}$
one can consider two options:
-
(i)
$F^{(0)} = F^{(0)}(\tau , T)$
, then
$C^{(0)}_{\tau } = -\gamma F^{(0)}$
, implying(3.15)assuming that
\begin{equation} \hat V^{(0)} = C^{(0)} = -\gamma \int _0^{\tau } F^{(0)}(\tilde \tau , T) {\rm d} \tilde \tau , \end{equation}
$F^{(0)}$
is such that the
$C^{(0)}$
is bounded;
-
(ii)
$F^{(0)} = F^{(0)}(t, \tau , T)$
such that
$\hat V^{(1)} = - \int _0^t [\gamma F^{(0)} (t, \tau , T) + C^{(0)}_{\tau }(\tau , T)] {\rm d} \tilde t$
is bounded, e.g.
${ \gamma } F^{(0)} = - C^{(0)}_{\tau } + \sin \omega t, \ \mbox{where} \ \omega = \omega (\tau , T)$
.
If
$F^{(0)} \ne 0$
, we can assume that
$F^{(1)} = 0$
, etc.
In what follows, we consider the first case. Then, summarising, we have
\begin{align} &\qquad\displaystyle \hat s^{(0)} = A(T) \cos \gamma [\tau + \phi (T)], \quad \hat v^{(0)} = - A(T) \sin \gamma [\tau + \phi (T)], \end{align}
\begin{align} &\,\,\,\displaystyle \hat V^{(0)} = -\gamma \int _0^{\tau } F^{(0)}(\tilde \tau , T) {\rm d} \tilde \tau , \ \ \tilde s^{(0)} = \frac {1}{c_0} \tilde \eta ^{(0)}(\xi , T), \ \ \tilde v^{(0)} = 0, \ \ \tilde V^{(0)} = 0; \end{align}
\begin{align} &\displaystyle \hat s^{(1)} = A^{(1)}(\tau , T), \quad \hat v^{(1)} = B^{(1)}(\tau , T), \quad \hat V^{(1)} = C^{(1)}(\tau , T), \quad \tilde s^{(1)} = \frac {1}{c_0} \tilde \eta ^{(1)}, \end{align}
\begin{align} &\qquad \displaystyle \tilde v^{(1)} = \frac {\gamma }{c_0^2} \left (\int _{-L}^\xi \tilde \eta ^{(0)}(\tilde \xi , T) {\rm d} \tilde \xi - \left \langle \int _{-L}^\xi \tilde \eta ^{(0)}(\tilde \xi , T) {\rm d} \tilde \xi \right \rangle \right ), \quad \tilde V^{(1)} = 0. \end{align}
Here,
$\left \langle \int _{-L}^\xi \tilde \eta ^{(0)}(\tilde \xi , T) {\rm d} \tilde \xi \right \rangle = ({1}/{2L}) \int _{-L}^{L} ( \int _{-L}^\xi \tilde \eta ^{(0)}(\tilde \xi , T) {\rm d} \tilde \xi ) {\rm d}\xi$
is the mean value.
Finally, collecting the terms at
$O(\alpha )$
, using (3.16)–(3.19) and simplifying the resulting system by virtue of its averaged members, we obtain
\begin{align} & -c_0 \tilde \eta ^{(2)}_{\xi } + c_0^2 \tilde s^{(2)}_{\xi } = - \frac {\sigma }{c_0} (\hat \eta ^{(0)} +\tilde \eta ^{(0)}) \tilde \eta ^{(0)}_\xi \nonumber \\ &\quad -\sigma \left ( A(T) \cos \gamma [\tau + \phi (T)] + \frac {1}{c_0} \tilde \eta ^{(0)} \right ) \tilde \eta ^{(0)}_\xi - F^{(0)} \tilde \eta ^{(0)}_\xi - \tilde \eta ^{(1)}_\tau - \tilde \eta ^{(0)}_T, \end{align}
\begin{align} &- c_0 \tilde s^{(2)}_{\xi } + \tilde \eta ^{(2)}_{\xi } = \frac {{ \gamma ^2}}{c_0^2} \left (\int _{-L}^\xi \tilde \eta ^{(0)}(\tilde \xi , T) {\rm d} \tilde \xi - \left \langle \int _{-L}^\xi \tilde \eta ^{(0)}(\tilde \xi , T) {\rm d} \tilde \xi \right \rangle \right ) - \frac {\beta c_0^2}{3} \tilde \eta ^{(0)}_{\xi \xi \xi } \nonumber \\ &\quad - { \frac {\sigma }{c_0}} \left (A(T) \cos \gamma [\tau + \phi (T)] + \frac {1}{c_0} \tilde \eta ^{(0)}\right ) \tilde \eta ^{(0)}_\xi - \frac {1}{c_0} \left(F^{(0)} \tilde \eta ^{(0)}_\xi { +} \tilde \eta ^{(1)}_\tau { +} \tilde \eta ^{(0)}_T\right), \end{align}
\begin{align} &\qquad\qquad\qquad\qquad\qquad\qquad\qquad { \tilde v^{(2)}_\xi = \frac {\gamma }{c_0^2} \tilde \eta ^{(1)},} \end{align}
\begin{align} &\qquad\qquad\qquad\qquad { \tilde V^{(2)}_\xi = - A(T) \sin \gamma [\tau + \phi (T)] \tilde \eta ^{(0)}_\xi ,} \end{align}
\begin{align} &\,\, \hat s^{(2)}_t = \gamma B^{(1)} - A^{(1)}_{\tau } - A_T \cos \gamma [\tau + \phi (T)] + \gamma A \phi _T \sin \gamma [\tau + \phi (T)], \end{align}
\begin{align} & \hat v^{(2)}_t= - \gamma A^{(1)} - B^{(1)}_{\tau } + A_T \sin \gamma [\tau + \phi (T)] + \gamma A \phi _T \cos \gamma [\tau + \phi (T)], \end{align}
\begin{align} &\qquad\qquad\qquad \hat V^{(2)}_t = - \gamma F^{(1)} - C^{(1)}_\tau + \gamma \int _0^{\tau } F^{(0)}_T(\tilde \tau , T) {\rm d} \tilde \tau . \end{align}
Equations (3.20) and (3.21) yield the equation
\begin{align} - 2 \tilde \eta ^{(1)}_\tau &= 2 \tilde \eta ^{(0)}_T + \left [\frac {\sigma \hat \eta ^{(0)}}{c_0} + 2 \sigma A \cos \gamma (\tau + \phi ) + 2 F^{(0)}(\tau , T) \right ] \tilde \eta ^{(0)}_\xi + \frac {3 \sigma }{c_0} \tilde \eta ^{(0)} \tilde \eta ^{(0)}_\xi \nonumber \\ &\quad +\frac {\beta c_0^3}{3} \tilde \eta ^{(0)}_{\xi \xi \xi } - \frac {\gamma ^2}{c_0} \left ( \int _{-L}^\xi \tilde \eta ^{(0)}(\tilde \xi , T) {\rm d} \tilde \xi - \left\langle \int _{-L}^\xi \tilde \eta ^{(0)}(\tilde \xi , T){\rm d} \tilde \xi \right\rangle \right ). \end{align}
Integrating (3.27) with respect to
$\tau$
and avoiding secular growth with
$\tau$
, we obtain
\begin{eqnarray} &&2 \tilde \eta ^{(0)}_T + \frac {\sigma \hat \eta ^{(0)}}{c_0} \tilde \eta ^{(0)}_\xi + \frac {3 \sigma }{c_0} \tilde \eta ^{(0)} \tilde \eta ^{(0)}_\xi + \frac { \beta c_0^3}{3} \tilde \eta ^{(0)}_{\xi \xi \xi } \nonumber \\ &&\quad -\frac {\gamma ^2}{c_0} \left (\int _{-L}^\xi \tilde \eta ^{(0)}(\tilde \xi , T) {\rm d} \tilde \xi - \left \langle \int _{-L}^\xi \tilde \eta ^{(0)}(\tilde \xi , T){\rm d} \tilde \xi \right \rangle \right ) = 0, \end{eqnarray}
and
\begin{equation} \tilde \eta ^{(1)} = - \left [\frac {\sigma A}{\gamma } \sin \gamma (\tau + \phi ) + \int _0^{\tau } F^{(0)}(\tilde \tau , T) {\rm d}\tilde \tau \right ] \tilde \eta _\xi ^{(0)}, \end{equation}
where we did not add an arbitrary zero-mean-field function since this just redefines
$\tilde \eta ^{(0)} (\xi , T)$
.
Differentiating (3.28) with respect to
$\xi$
and changing the variable
$\xi$
to
$\zeta = \xi - ({\sigma \hat \eta ^{(0)}}/{2c_0}) T$
leads to the traditional form of the Ostrovsky equation
\begin{eqnarray} \left (\tilde \eta ^{(0)}_T + \alpha _1 \tilde \eta ^{(0)} \tilde \eta ^{(0)}_\zeta + \beta _1 \tilde \eta ^{(0)}_{\zeta \zeta \zeta } \right )_\zeta = \gamma _1 \tilde \eta ^{(0)}, \end{eqnarray}
where
\begin{eqnarray} \alpha _1 = \frac {3 \sigma }{2 c_0}, \quad \beta _1 = \frac { \beta c_0^3}{6}, \quad \gamma _1 = \frac {\gamma ^2}{2 c_0}. \end{eqnarray}
Next, substituting (3.29) for
$\tilde \eta ^{(1)}$
into (3.22) for
$\tilde v^{(2)}_{\xi }$
we conclude that
\begin{equation} \tilde v^{(2)} = - \frac {1}{c_0^2} \left [\sigma A \sin \gamma (\tau + \phi ) + \gamma \int _0^\tau F^{(0)}(\tilde \tau , T) {\rm d} \tilde \tau \right ] \tilde \eta ^{(0)}, \end{equation}
using the zero-mean condition once again. Similarly, (3.23) yields
\begin{equation} \tilde V^{(2)} = - A \sin \gamma (\tau + \phi ) \tilde \eta ^{(0)}. \end{equation}
Finally, the equations for the mean values yield
$A_T = \phi _T = 0$
in order to avoid secular growth, and then simplify to take the form
\begin{align} \hat s^{(2)}_t &= \gamma B^{(1)} - A^{(1)}_{\tau }, \quad \hat v^{(2)}_t = -\gamma A^{(1)} - B^{(1)}_{\tau }, \nonumber\\ \hat V^{(2)}_t &= -\gamma F^{(1)} - C^{(1)}_{\tau } + \gamma \int _0^{\tau } F_T^{(0)} (\tilde \tau , T) {\rm d} \tilde \tau . \end{align}
Assuming that
$F^{(1)} = F^{(1)}(\tau , T)$
, we then have
\begin{equation} A^{(1)}_{\tau } = \gamma B^{(1)}, \quad B^{(1)}_{\tau } = -\gamma A^{(1)}, \quad C^{(1)}_{\tau } = - \gamma F^{(1)}+ \gamma \int _0^{\tau } F_T^{(0)} (\tilde \tau , T) {\rm d} \tilde \tau . \end{equation}
In what follows, we consider
$(F, A, B)^{(1)} = 0$
, while
\begin{equation} C^{(1)} = \gamma \int _0^{\tau } \left ( \int _0^{\hat \tau } F_T^{(0)} (\tilde \tau , T) {\rm d}\tilde \tau \right ) {\rm d} \hat \tau , \end{equation}
provided this function is bounded. A sufficient condition for the latter is given by
$F_T^{(0)} = 0$
, implying that
$F^{(0)} = F^{(0)} (\tau )$
, which we assume here.
Hence, considering the equations up to
$O(\alpha )$
, i.e. up to the accuracy of the governing equations (2.11)–(2.18), allows us to fully define all terms at
$O(1)$
and
$O(\sqrt {\alpha })$
. The procedure can be continued to any order, but instead of using the truncated weakly nonlinear formulation, we would need to use the original strongly nonlinear equations.
To summarise, up to the accuracy of the governing equations, we obtained the following large class of uni-directional waves described by the Ostrovsky equation and propagating over the non-zero, generally evolving, mean fields:
\begin{align} \eta &= \hat \eta ^{(0)} + \tilde \eta ^{(0)} - \sqrt {\alpha } \left [\frac {\sigma A}{\gamma } \sin \gamma (\tau + \phi ) + \int _0^\tau F^{(0)}(\tilde \tau ) {\rm d} \tilde \tau \right ] \tilde \eta ^{(0)}_{\zeta } + O(\alpha ), \end{align}
\begin{align} s &= A \cos \gamma (\tau + \phi ) + \frac {1}{c_0} \tilde \eta ^{(0)} -\sqrt {\alpha } \frac {1}{c_0} \left [\frac {\sigma A}{\gamma } \sin \gamma (\tau + \phi ) + \!\int _0^\tau \!F^{(0)}(\tilde \tau ) {\rm d} \tilde \tau \right ] \tilde \eta ^{(0)}_{\zeta } + O(\alpha ), \end{align}
\begin{align} v &= - A \sin \gamma (\tau + \phi ) +\sqrt {\alpha } \frac {\gamma }{c_0^2} \left [\int _{-L}^\zeta \tilde \eta ^{(0)}(\tilde \zeta , T) {\rm d} \tilde \zeta - \left \langle \int _{-L}^\zeta \tilde \eta ^{(0)}(\tilde \zeta , T) d \tilde \zeta \right \rangle \right ] + O(\alpha ), \end{align}
\begin{align} V &= -\gamma \int _0^{\tau } F^{(0)}(\tilde \tau ) {\rm d} \tilde \tau + O(\alpha ), \end{align}
where
$\hat \eta ^{(0)},\, A,\, \phi$
are arbitrary constants,
$F^{(0)}(\tilde \tau )$
is a function such that
$\int _0^\tau F^{(0)}(\tilde \tau ) {\rm d} \tilde \tau$
is bounded (e.g.
$\sin \omega \tilde \tau$
) and
$\tilde \eta ^{(0)} = \tilde \eta ^{(0)}(\zeta , T)$
satisfies the Ostrovsky equation (3.30), where
\begin{equation} \zeta = \xi - \frac {\sigma \hat \eta ^{(0)}}{2 c_0} T = x - \left (c_0 + \frac {\alpha \sigma \hat \eta ^{(0)}}{2 c_0}\right ) t, \quad T = \alpha t. \end{equation}
We note that this weakly nonlinear solution can be rewritten in a more convenient and asymptotically equivalent form
\begin{align}&\qquad\qquad\qquad\qquad\qquad \eta = \hat \eta ^{(0)} + \tilde \eta ^{(0)}\left (T, \theta \right ) + O(\alpha ), \end{align}
\begin{align}&\qquad\qquad\qquad\quad s = A \cos \gamma (\tau + \phi ) + \frac {1}{c_0} \tilde \eta ^{(0)}\left (T, \theta \right ) + O(\alpha ), \end{align}
\begin{align}& v = - A \sin \gamma (\tau + \phi ) +\sqrt {\alpha } \frac {\gamma }{c_0^2} \left [\int _{-L}^\theta \tilde \eta ^{(0)}(\tilde \theta , T) {\rm d} \tilde \theta - \left \langle \int _{-L}^\theta \tilde \eta ^{(0)}(\tilde \theta , T) {\rm d} \tilde \theta \right \rangle \right ] + O(\alpha ), \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad V = -\gamma \int _0^{\tau } F^{(0)}(\tilde \tau ) {\rm d} \tilde \tau + O(\alpha ), \end{align}
where
\begin{equation} \theta = \zeta - \frac {\sqrt {\alpha } \sigma A}{\gamma } \sin \gamma (\tau + \phi ) - \sqrt {\alpha } \int _0^\tau F^{(0)}(\tilde \tau ) {\rm d} \tilde \tau , \end{equation}
and
\begin{eqnarray} \left (\tilde \eta ^{(0)}_T + \alpha _1 \tilde \eta ^{(0)} \tilde \eta ^{(0)}_\theta + \beta _1 \tilde \eta ^{(0)}_{\theta \theta \theta } \right )_\theta = \gamma _1 \tilde \eta ^{(0)}, \end{eqnarray}
with
$ \alpha _1 = ({3 \sigma }/{2 c_0}), \; \beta _1 = ({ \beta c_0^3}/{6}), \; \gamma _1 = ({\gamma ^2}/{2 c_0}).$
Note that we should not add any additional transport terms to (3.47) since we simply replace
$\zeta$
with
$\theta$
in the solution of (3.30) in order to combine the
$O(1)$
and
$O(\sqrt {\alpha })$
terms in the asymptotic expansions. This representation, where the effective phase
$\theta$
now depends on
$\zeta$
and
$\tau$
, gives a clear description of the main effects related to non-zero mass: generally, there appear
$O(1)$
inertial oscillations in both shear variables, as well as
$O(\sqrt {\alpha })$
oscillations with the same frequency in the phase of the interfacial displacement and shear. Non-zero barotropic transport
$F^{(0)}$
in the propagation direction also gives an
$O(\sqrt {\alpha })$
shift in the phase of the interfacial displacement and shear in the propagation direction, as well as giving an
$O(1)$
contribution to the transverse barotropic transport.
Thus, the solution
$\tilde \eta ^{(0)}$
of the Ostrovsky equation (3.47), has zero mean by construction,
$\eta$
has an arbitrary constant mean value, while the shear variables
$s$
and
$v$
generally have time-dependent mean values. The barotropic transport variables
$F$
and
$V$
are also generally time-dependent.
We finish the section by considering the limit of the constructed weakly nonlinear solution (3.42)–(3.47) to the case when there is no rotation. Assuming
$\phi = 0$
, consider a finite value of
$\tau$
and let
$\gamma \to 0$
. Then
\begin{align} \eta &= \hat \eta ^{(0)} + \tilde \eta ^{(0)}\left (T, \theta \right ) + O(\alpha ), \end{align}
\begin{align} s &= A + \frac {1}{c_0} \tilde \eta ^{(0)}\left (T, \theta \right ) + O(\alpha ), \end{align}
\begin{align} v &= O(\alpha ), \quad V = O(\alpha ), \end{align}
where
\begin{align} \theta = \zeta - \sigma A T - \sqrt {\alpha } \int _0^\tau F^{(0)}(\tilde \tau ) {\rm d} \tilde \tau , \end{align}
and
\begin{eqnarray} \tilde \eta ^{(0)}_T + \alpha _1 \tilde \eta ^{(0)} \tilde \eta ^{(0)}_\theta + \beta _1 \tilde \eta ^{(0)}_{\theta \theta \theta } = 0 ,\end{eqnarray}
with
$ \alpha _1 = {3 \sigma }/{2 c_0}, \; \beta _1 = { \beta c_0^3}/{6}$
, is the KdV equation.
In the rest of this paper we consider the case
$F^{(0)} = 0$
(i.e. the barotropic transport is absent). A particular choice
$A = {\hat \eta ^{(0)}}/{c_0}$
corresponds to the case when there is no background shear flow, and then, in the KdV equation (3.52),
\begin{equation} \theta = \xi - \alpha _1 \hat \eta ^{(0)} T = x - \left ( c_0 + \frac {3 \alpha \sigma \hat \eta ^{(0)}}{2 c_0} \right ) t. \end{equation}
In general,
$A = ({\hat \eta ^{(0)}}/{c_0}) + s_0$
, where
$s_0 = \textrm {const}$
represents a constant background shear, and then we have, again in the absence of rotation,
\begin{equation} \theta = \xi - (\alpha _1 \hat \eta ^{(0)} + \sigma s_0) T = x - \left ( c_0 + \frac {3 \alpha \sigma \hat \eta ^{(0)}}{2 c_0} + \alpha \sigma s_0 \right ) t. \end{equation}
Suppose one allows for the presence of a weak piecewise-constant shear flow in the two layers and assume that
$s_0 \gt 0$
. In the absence of rotation, sufficient conditions allowing one to avoid the appearance of the long-wave instability and critical levels are given, in the non-dimensional variables used in this paper, by the conditions
\begin{equation} { s_0 \lt \left (\frac {\rho _1 (1-h_0) + \rho _2 h_0}{ \rho _2} \right ) ^{1/2} \quad \mbox{and} \quad s_0 \lt \left (\frac {\rho _1 h_0 (1-h_0)}{\rho _1 (1-h_0) + \rho _2 h_0} \right )^{1/2}} \end{equation}
(see Ovsyannikov Reference Ovsyannikov1979, Reference Ovsyannikov1985; Bontozoglou Reference Bontozonglou1991; Boonkasame & Milewski Reference Boonkasame and Milewski2011; Lannes & Ming Reference Lannes and Ming2015; Khusnutdinova & Zhang Reference Khusnutdinova and Zhang2016 for the details and the necessary conditions), where the second condition is stronger than the first. In the Boussinesq approximation used in this paper it implies
\begin{equation} { s_0 \lt c_0 = \left [h_0 (1-h_0)\right ]^{1/2}.} \end{equation}
In the following sections we consider the basic case
$s_0 = 0$
.
The constructed weakly nonlinear solution can be used to describe waves of small amplitude. The conservation of energy to leading order was verified using the asymptotic approximation of Helfrich (Reference Helfrich2007). Fully nonlinear waves in the absence of rotation are significantly different from the weakly nonlinear description provided by the KdV equation (see Jo & Choi Reference Jo and Choi2002, Reference Jo and Choi2008; Camassa et al. Reference Camassa, Choi, Michallet, Rusas and Sveen2006; Choi Reference Choi2006; Camassa et al. Reference Camassa, Rusas, Saxena and Tiron2010; Barros, Choi & Milewski Reference Barros, Choi and Milewski2020; Doak, Barros & Milewski Reference Doak, Barros and Milewski2022 and references therein). The extended weakly nonlinear models can be used to describe MMCC waves of greater amplitude than that described by the KdV equation (Sidorovas et al. Reference Sidorovas, Tseluiko, Choi and Khusnutdinova2025), paving the way for similar extensions in the presence of rotation.
4. The effect of rotation on the evolution of cnoidal waves with defects
The MMCC system admits periodic nonlinear travelling wave solutions (Camassa et al. Reference Camassa, Rusas, Saxena and Tiron2010) (see also Sidorovas et al. Reference Sidorovas, Tseluiko, Choi and Khusnutdinova2025 for the weakly nonlinear modelling of the cnoidal waves of small and moderate amplitude). The aim of this section is to model the effect of rotation on the evolution of the waves generated by initial conditions close to the cnoidal waves of the KdV equation but having local amplitude and/or periodicity defects. We use the weakly nonlinear solution constructed in § 3 and first model the effect of rotation on the exact classical solutions of the KdV equation: solitons and cnoidal waves (Boussinesq Reference Boussinesq1871; Korteweg & de Vries Reference Korteweg and de Vries1895), as well as bright and dark breathers (Kuznetsov & Mikhailov Reference Kuznetsov and Mikhailov1975; Hoefer et al. Reference Hoefer, Mucalica and Pelinovsky2023). Then, we consider cnoidal waves with the local expansion and contraction periodicity defects introduced by cutting at the trough and symmetrically inserting a piece of a straight line, and by symmetrically cutting away a part close to the trough between the two neighbouring peaks and gluing together the remaining parts of the solution, respectively (see the first two rows in figure 17 of Appendix A). This allows us to set up a framework for the discussion of the effects observed in our subsequent modelling of cnoidal waves with generic localised perturbations. As a by-product of our study involving initial conditions with local periodicity defects we came across an observation that cnoidal waves with expansion defects can be viewed as generalised (‘shock-like’) travelling waves of the KdV equation. In Appendix A, we prove that all conservation laws of the KdV equation (infinitely many) are identically satisfied for them, provided that we understand the conserved quantities as the natural sum of integrals, treating the points of discontinuity similarly to shocks. Moreover, a cnoidal wave with an expansion defect has a continuous first derivative, satisfying the Weirstrass–Erdmann corner condition for the associated variational problem. We also make comments about a possible weak formulation.
In our numerical experiments, we use two sets of parameters. In the first set, which is used for the majority of our simulations (except the runs shown in figures 13, 14(b), 15(b) and 16), the pycnocline is closer to the surface, and internal waves are the waves of depression. Hence, we show the bottom view in the majority of our figures in this section. The depths of the upper and lower layers are
$37.5$
and
$112.5$
m, respectively, and the total depth of the system, denoted by
$H$
, is
$150.0$
m. The densities of the upper and lower layers are chosen to be
$\rho _{1} = 1000.0 ,\, \rho _{2} = 1003.1$
kg m
$^{-3}$
, respectively. The small parameters are defined as
$\alpha = a/H$
,
$\tilde \beta = (H/l)^2$
and
$\tilde \gamma =(lf)/\sqrt {(g'H)}$
, where
$a$
is a wave amplitude,
$l$
is a wavelength,
$f$
is the Coriolis parameter,
$g' = g\Delta \rho /\rho _1$
is the reduced gravity. We let
$\alpha = 0.005$
and
$\tilde \beta = 0.030$
. We assume midlatitude oceanic values for the Coriolis parameter
$f = 5\times 10^{-5} \text{ s}^{-1}$
, and reduced gravity
$g'= 0.030 \text{ m s}^{-2}$
. Hence,
$\tilde \gamma \approx 0.020$
. The non-dimensional unperturbed upper layer depth
$h_0=0.250$
gives values of the linear long-wave speed
$c_0 = (h_0-h_0^{2})^{1/2} \approx 0.433$
and coefficient
$\sigma = 2h_0-1 = -0.500$
. The scaled
$\mathcal{O}(1)$
parameters are
$\beta = \tilde \beta /\alpha$
and
$\gamma = \tilde \gamma /\sqrt {\alpha }$
. Then, the coefficients of the Ostrovsky equation (3.47), are given by
$\alpha _1 \approx -1.732, \beta _1 \approx 0.081, \gamma _1 \approx 0.096$
. This regime is close to one of the regimes considered by Helfrich (Reference Helfrich2007), and the internal waves are waves of depression. The second set of parameters is used to model the internal waves of elevation shown in figures 13, 14(b), 15(b) and 16, and here the pycnocline is closer to the bottom: the depths of the upper and lower layers are
$120.0$
and
$30.0$
m, respectively. The total depth of the system is again
$H=150.0$
m. The other parameters are unchanged except the non-dimensional unperturbed upper layer depth
$h_0=0.800$
, giving values of the linear long-wave speed
$c_0 = (h_0-h_0^{2})^{1/2} = 0.400$
and coefficient
$\sigma = 2h_0-1 = 0.600$
. Then, the coefficients of the Ostrovsky equation (3.30) are found to be
$\alpha _1 = 2.250, \beta _1 = 0.064, \gamma _1 \approx 0.104$
. The Ostrovsky equation is solved using a pseudospectral method described in Appendix B. The constant shear is absent,
$s_0=0$
, and the initial phase of the inertial oscillations is
$\phi =0$
.
4.1. The effect of rotation on solitons and cnoidal waves
We begin by modelling the evolution of initial conditions in the form of solitary and cnoidal wave solutions of the KdV equation (3.52) associated with the Ostrovsky equation (3.47), which we write in the form
\begin{eqnarray} \eta ^{(0)}_T + \alpha _1 (\eta ^{(0)} - \hat \eta ^{(0)}) \eta ^{(0)}_\theta + \beta _1 \eta ^{(0)}_{\theta \theta \theta } = 0, \end{eqnarray}
where
$\eta ^{(0)} = \hat \eta ^{(0)} + \tilde \eta ^{(0)}$
. We consider the initial condition for
$\eta$
in the form of a soliton solution
\begin{align} \eta ^{(0)} \Big |_{T=0} = \dfrac {3 v_s}{\alpha _1} \text{sech}^2 \left [ \dfrac {1}{2} \sqrt {\dfrac {v_s}{\beta _1}}(\theta + \alpha _1 \hat \eta ^{(0)} T -v_s T-\zeta _0) \right ]_{T=0}, \end{align}
where
$\zeta _0$
and
$v_s\gt 0$
are arbitrary constants, and the cnoidal wave solution
\begin{align} \eta ^{(0)} = \dfrac {6\beta _1}{\alpha _1}\left \{ u_2 + (u_3 - u_2) \text{ cn}^2\left [ (\theta + \alpha _1 \hat \eta ^{(0)} T - v_c T - \zeta _0) \sqrt {\dfrac {u_3 -u_1}{2}} ; m \right ] \right \}_{T=0}, \end{align}
given in terms of the Jacobi elliptic functions, where
$u_1 \lt u_2 \lt u_3$
are real,
$v_c = 2 \beta _1 (u_1 + u_2 + u_3)$
and
$m ={(u_3-u_2)}/{(u_3-u_1)}$
is the elliptic parameter. We recall that, in the case under study
\begin{equation} \theta = \zeta - \frac {2 \alpha _1 \hat \eta ^{(0)}}{3} T ,\end{equation}
and, accounting for the transport term in (4.1), we get the phase
\begin{equation} \theta + \alpha _1 \hat \eta ^{(0)} T - v_{s/c} T - \zeta _0 = \zeta + \frac {\alpha _1 \hat \eta ^{(0)}}{3} T - v_{s/c} T - \zeta _0, \end{equation}
choosing the variable
$\zeta$
as a convenient variable defining the speed of the moving reference frame in both cases, with and without rotation. Here,
$ \zeta + ({\alpha _1 \hat \eta ^{(0)}}/{3}) T = \xi = x - c_0 t$
. In the limit
$m \to 1$
, the cnoidal wave approaches a soliton (generally, on a non-zero pedestal). The period of the cnoidal wave is given by
\begin{align} \lambda = \dfrac {2K(m)}{\sqrt {(u_3 - u_1)/2}}, \end{align}
where
$K(m)$
is the complete elliptic integral of the first kind (see, for example, Kamchatnov Reference Kamchatnov2000). We use the same initial conditions for
$\eta ^{(0)}$
when rotation is present, and obtain the corresponding weakly nonlinear solutions using (3.42)–(3.46) and the Ostrovsky equation (3.47).
Figure 3. The effect of rotation on the interfacial displacement for the KdV soliton initial condition in simulations with periodic boundary conditions. (a) view from above, (b) view from below, (c) in simulations without the sponge layers, (d) in simulations with the sponge layers, (e,f,g) comparison of the lead wavepacket in simulations without the sponge layers (black) and with the sponge layers (red dashed) at different moments of time.
In figure 3, we compare the evolution of a single soliton in a two-layer fluid system with rotation in numerical experiments with periodic boundary conditions, either with or without the sponge layers near the boundaries (see, for example, Alias et al. Reference Alias, Grimshaw and Khusnutdinova2013 for the discussion of the sponge layers). The sponge layers near the boundaries act as a filter absorbing the radiation. The computational domain is
$2L = 100$
, with the number of modes
$M = 998$
, the spatial step
$\Delta x \approx 10^{-1}$
, the total simulation time
$T_{max} = 200$
and the temporal step
$\Delta T = 10^{-2}$
. Other constants are
$v_s \approx 0.487$
and
$\zeta _0 = 0$
.
Figure 4. The effect of rotation on the interfacial displacement for the KdV cnoidal wave initial condition: (a) view from above and (b) view from below.
Figure 3(a,b) shows three-dimensional plots of the time evolution of the interfacial displacement for initial conditions in the form of a single soliton, the views from above (a) and below (b). Figure 3(c,d) compares the evolution of the same initial conditions without (c) and with (d) the sponge layers near the boundaries of the periodic domain. The soliton evolves into an unsteady wavepacket that has been intensively studied previously (see Grimshaw & Helfrich Reference Grimshaw and Helfrich2008; Grimshaw et al. Reference Grimshaw R.H., Helfrich and Johnson2013; Ostrovsky et al. Reference Ostrovsky, Pelinovsky, Shrira and Stepanyants2015; Stepanyants Reference Stepanyants2020, and references therein), and for a long time the radiation penetrating through the boundaries in the absence of sponge layers does not have a large effect on the shape or speed of the main wavepacket, which is shown in more detail in figure 3(e,f,g). The overall wave pattern remains stable and consistent throughout the simulation.
Figure 5. The effect of rotation on the evolution of the maximum, minimum and amplitude of the interfacial displacement for a soliton (a) and cnoidal wave (b) initial condition.
Similarly, in figure 4 we show the effect of rotation for the cnoidal wave initial condition. The computational domain is
$2L = 26.40$
, with the number of modes
$M = 262$
, the spatial step
$\Delta x \approx 10^{-1}$
, the total simulation time
$T_{max} = 100$
and the temporal step
$\Delta T = 10^{-2}$
. The parameters characterising the initial condition are
$u_1 = -10^{-3}, u_2 = 0$
and
$u_3 = 3$
. It is evident that rotation leads to formation of a rather regular quasi-periodic wave pattern shown in figure 4, without the formation of bursts of large waves. Overall, this wave looks similar to a cnoidal wave, but it has a slowly oscillating amplitude, resembling the modulationally stable periodic waves in Whitfield & Johnson (Reference Whitfield and Johnson2014).
We also note that, both for the soliton and cnoidal wave initial conditions, the time evolution of the maximum and minimum of the free-surface elevation is quasi-periodic, with no significant bursts, as shown in figure 5. All computational parameters are the same as in the previous figures. Hence, we conclude that rotation alone does not trigger the formation of large waves in these simulations, but it does lead to the emergence of a left-propagating (i.e. slower moving) wavepacket for initial condition in the form of a soliton.
4.2. The effect of rotation on breathers on a cnoidal wave background
The exact breather on a cnoidal wave solution of the KdV equation was studied by Kuznetsov & Mikhailov (Reference Kuznetsov and Mikhailov1975) and Hoefer et al. (Reference Hoefer, Mucalica and Pelinovsky2023). There are two types of breathers: bright and dark (this terminology refers to the canonical form of the KdV equation, where solitons are waves of elevation). For a bright breather propagating as a dislocation of a cnoidal wave, the exact solution of the KdV equation (4.1) takes the form
\begin{align} \eta ^{(0)} = \frac {12 \beta _1}{\alpha _1} \left \{ k^2 -1 +\frac {E(k)}{K(k)} +\frac {\partial ^2}{\partial \zeta ^2} \left [ \log \left ( \tau ({ \zeta + \frac {\alpha _1 \hat \eta ^{(0)}}{3} T}, \beta _1 T) \right ) \right ] \right \}\!, \end{align}
where
$k\in (0,1)$
is the elliptic modulus (the elliptic parameter
$m = k^2$
),
$K(k)$
is the complete elliptic integral of the first kind,
$E(k)$
is the complete elliptic integral of the second kind and the
$\tau\hbox{-}$
function is given by
\begin{align} \tau (x,t) &:= \varTheta (x-c_0t +\alpha _b) \exp \{\kappa _b(x-c_bt+x_0)\}\nonumber\\&\quad + \varTheta (x-c_0t -\alpha _b) \exp \{-\kappa _b(x-c_bt+x_0)\} ,\end{align}
with
$\kappa _b\gt 0,\, c_b \gt c_0$
and
$\alpha _b\in (0,K(k))$
. Here,
$\varTheta (x) = \theta _4 ({\pi x}/{2K(k)} )$
, where
$\theta _4$
is the Jacobi theta function of type four, given by
$\theta _4(u) = 1+2 \sum \nolimits _{n=1}^{\infty } (-1)^n q^{n^2}\cos (2nu)$
. The solution is parameterised by
$\lambda \in (-\infty ,-k^2)$
. The parameters are defined as follows (Hoefer et al. Reference Hoefer, Mucalica and Pelinovsky2023):
\begin{align} \varphi _\gamma &= \arcsin \left ( \frac {\sqrt {-\lambda -k^2}}{\sqrt {1-2k^2-\lambda }} \right ), \end{align}
\begin{align} \alpha _b &= F(\varphi _\gamma ,k),\quad \kappa _b = \frac {\sqrt {1-\lambda -k^2}\sqrt {-\lambda -k^2}}{\sqrt {1-2k^2-\lambda }} - Z(\varphi _\gamma ,k), \end{align}
\begin{align} c_0 &= 4(2k^2-1),\quad c_b = c_0 + \frac {4\sqrt {1-\lambda -2k^2} \sqrt {1-\lambda -k^2} \sqrt {-\lambda -k^2}}{\kappa _b}, \end{align}
where
$F(\varphi ,k):= \int _0^\varphi {{\rm d}\theta }/{\sqrt {1-k^2 \sin ^2\theta }}$
is the elliptic integral of the first kind, and
$Z(\varphi ,k):= \int _0^\varphi \sqrt {1-k^2 \sin ^2\theta }{\rm d}\theta$
is the elliptic integral of the second kind.
For a dark breather propagating as a dislocation of a cnoidal wave the exact solution is again given by (4.7), where the solution is now parameterised by
$\lambda \in (1-2k^2, 1-k^2)$
, and the
$\tau $
-function takes the form
\begin{align} \tau (x,t) &:= \varTheta (x-c_0t +\alpha _d) \exp \{-\kappa _d(x-c_dt+x_0)\} \nonumber\\&\quad + \varTheta (x-c_0t -\alpha _d) \exp \{\kappa _d(x-c_dt+x_0)\} ,\end{align}
with
$\kappa _d\gt 0,\, c_d \lt c_0$
and
$\alpha _d\in (0,K(k))$
. The parameters are given by
\begin{align}&\qquad\qquad\qquad\qquad\quad \varphi _\alpha = \arcsin \Big ( \frac {\sqrt {1-k^2-\lambda }}{k} \Big ),\quad \end{align}
\begin{align}&\qquad\qquad\qquad\qquad \alpha _d = F(\varphi _\alpha ,k),\quad \kappa _d = Z(\varphi _\alpha ,k),\quad \end{align}
\begin{align}& c_0 = 4(2k^2-1),\quad c_d = c_0 - \frac {4\sqrt {(k^2+\lambda ) (\lambda -1+2k^2) (1-k^2 -\lambda )}}{\kappa _d}. \end{align}
Figure 6. Numerical solution for a bright (a,c,e,g) and dark (b,d,f,h) breather on a cnoidal wave initial condition (view from below): (a,b) interfacial displacement in the absence of rotation, (c,d) interfacial displacement under the effect of rotation, (e,f) / (g,h) shear in the direction of wave propagation / orthogonal direction, under the effect of rotation.
Figure 7. The effect of rotation on the evolution of the maximum, minimum and amplitude of the interfacial displacement for a bright (a) and a dark (b) breather on a cnoidal wave initial condition. Black solid, red dashed and blue dot-dashed lines correspond to 5, 7 and 9 peaks in the domain, respectively.
Numerical solutions initiated with the bright and dark breather on a cnoidal wave initial condition at
$T=0$
in the absence of rotation (i.e. in the KdV regime) are shown in the figure 6(a,b). For the bright-breather simulations, shown in the figure 6(a,c,e,g), the computational domain length is
$2L = 72.86$
, with the number of modes
$M=728$
, the spatial step
$\Delta x \approx 10^{-1}$
, the total simulation time
$T_{max} = 100$
and the temporal step
$\Delta T = 10^{-2}$
. The parameters used in the initial condition are
$[k,\lambda ] \approx [0.9998, -1.30]$
. For the dark-breather simulations, shown in the figure 6(b,d,f,h), the computational domain length is
$2L = 80.64$
, with the number of modes
$M=804$
, the spatial step
$\Delta x \approx 10^{-1}$
, the total simulation time
$T_{max} = 100$
and the temporal step
$\Delta T = 10^{-2}$
. The other parameters are
$[k,\lambda ] \approx [0.9998, -0.50]$
. In the absence of rotation, the figure 6(a,b) shows stable propagation of the bright- and dark-breather solutions on top of the cnoidal wave, in good agreement with the available analytical solution (4.7).
Next, the effect of rotation on the evolution of the bright and dark breather on a cnoidal wave initial condition is shown in the subsequent rows of figure 6, for the interfacial displacement (figure 6 c,d), shear in the direction of wave propagation (figure 6 e,f) and shear in the orthogonal direction (figures 6 g,h). We notice a striking difference from the results shown in the figure 6(a,b) (no rotation): there emerges a rather strong left-propagating localised burst clearly visible both in the free-surface elevation (figure 6 c,d) and the shear in the direction of wave propagation (figure 6 e,f). The counterpart of the wave is also present in the shear in the orthogonal direction (figure 6 g,h), but this signal is weak. It must be noted that the wave forms soon after the initiation of the simulation, and it continues to grow for a long time.
In these runs, the effect of rotation on a dark breather on a cnoidal wave is qualitatively similar to that on the bright breather, although less pronounced. In both cases, the background cnoidal wave continues propagating to the right (in the moving reference frame), and there appears a significant burst both in the interfacial displacement and the shear in the direction of wave propagation, moving to the left. For the dark breather, the signal in the shear in the orthogonal direction is also present, but it is barely noticeable to the naked eye.
Figure 7 shows the time evolution of the maximum, minimum and the amplitude of the interfacial displacement for both cases of a bright (figure 7
a) and a dark (figure 7
b) breather on a cnoidal wave initial condition. We see that the observed effects of rotation are structurally stable with regard to the size of the computational domain. We experimented with 5, 7 and 9 major peaks in the domain. In the bright-breather case, the computational domains are
$2L= 37.42$
(5 peaks),
$2L = 55.02$
(7 peaks),
$2L = 72.86$
(9 peaks). For the dark-breather case, the computational domains are
$2L= 45.44$
(5 peaks),
$2L = 63.04$
(7 peaks),
$2L = 80.64$
(9 peaks). The remaining numerical parameters are the same as before. It is again evident that the wave amplitude grows. The burst forms soon after the initiation of the numerical runs. The close-up views of the large waves visible in the free-surface elevation are shown in figure 8 for the time around
$T=40$
(bright-breather case, figure 8
a) with approximately
$29\,\%$
increase in the amplitude compared with the initial condition, and around
$T=9$
(dark-breather case, figure 8
b) with approximately
$15\,\%$
increase in the amplitude compared with the initial condition. Hence, we conclude that, under the effect of rotation, the moving dislocation on top of the otherwise regular cnoidal wave can lead to the emergence of strong bursts of interfacial waves and shear currents in the direction of propagation of the cnoidal wave.
Figure 8. Close-up view from below of the large burst in the interfacial displacement for a bright (a) and a dark (b) breather on a cnoidal wave initial condition around
$T = 40$
(a) and
$T = 9$
(b), respectively.
We also note that we performed preliminary simulations with higher-amplitude breathers, and these suggest a more complex dynamics: the burst initially propagates to the right, but part of the mass is subsequently emitted into the background cnoidal wave, and the burst subsequently decreases in amplitude and changes direction. A detailed investigation of this behaviour is beyond the scope of the present work and is left for future work.
4.3. The effect of rotation on cnoidal waves with periodicity defects
Recent research related to the wave packets described by the Schrödinger equation has shown that localised phase defects can lead to the emergence of rogue waves (He et al. Reference He, Witt, Trillo, Chabchoub and Hoffmann2022). Here, we investigate whether periodicity defects introduced into the long cnoidal waves in the KdV–Ostrtovsky regime can also lead to the emergence of large localised bursts of energy, under the effect of rotation. Also, can it happen already in the absence of rotation, i.e. in the KdV regime?
Motivated by figure 1, we consider cnoidal waves close to their solitonic limit and introduce two types of periodicity defects (see figure 17 in Appendix A): contraction and expansion defects, depending on whether the distance between the two neighbouring peaks is shorter or longer than the period of the cnoidal wave. The first defect is introduced by symmetrically cutting away a small part close to the trough between the two neighbouring peaks, and gluing together the remaining parts of the solution. Naturally, the resulting function has discontinuous first derivative at one point within the computational domain, but the jump in the derivative is small because we have cut close to extremum (smoothed in numerical simulations). The second defect is introduced by cutting the graph at the trough and symmetrically inserting a piece of a straight line. The resulting function has continuous first derivative, and discontinuous second derivative at two points within the computational domain. To our surprise, in the absence of rotation (KdV regime), these functions evolved almost like travelling waves of the KdV equation. Both types of waves were long lived, and a cnoidal wave with an expansion defect, smoothed in a pseudospectral simulation, was stable, with no visible changes at the end of the long run (the difference between the numerical solution and the initial condition travelling at the speed of the cnoidal wave was of order
$10^{-4}$
). The numerical results initiated with such initial conditions in the absence of rotation are shown in the first row of figure 9 for the contraction defect (figure 9
a) and the expansion defect (figure 9
b). Intrigued by this observation, we managed to prove that, for the constructed functions (and their natural generalisations), all (infinitely many) conservation laws of the KdV equation, understood in the sense of a sum of integrals, are satisfied exactly. This approach was inspired by the recent work by Gavrilyuk et al. (Reference Gavrilyuk, Nkonga, Shyue and Truskinovsky2020), Gavrilyuk & Shyue (Reference Gavrilyuk and Shyue2022) and Gavrilyuk, Nkonga & Shyue (Reference Gavrilyuk, Nkonga and Shyue2024), where interesting generalised solutions with singularities satisfying the Weierstrass–Erdmann corner conditions (e.g. Fox Reference Fox1954) and requiring the continuity of the derivative at the junction, were constructed in the context of non-integrable Benjamin–Bona–Mahoney and conduit equations. To the best of our knowledge, the discontinuous generalised (shock-like) travelling wave solutions of the KdV equation in the form of cnoidal waves with periodicity defects constructed in our paper have not been discussed before. A cnoidal wave with an expansion defect satisfies the corner condition exactly, it’s smoothed counterpart was extremely stable in our numerical runs, but a smoothed counterpart of a cnoidal wave with a contraction defect also turned out to be long lived. Further discussion of solutions of the KdV equation with defects (both stable and unstable) can be found in Appendix A. Importantly, we conclude that, for the duration of our long simulations, such natural periodicity defects do not lead to the focusing of energy in the absence of rotation.
Figure 9. Numerical solution for a cnoidal wave with a contraction (a,c,e,g) and expansion (b,d,f,h) defect initials condition (view from below): (a,b) interfacial displacement in the absence of rotation, (c,d) interfacial displacement under the effect of rotation, (e,f) / (g,h) shear in the direction of wave propagation / orthogonal direction, under the effect of rotation.
Figure 10. The effect of rotation on the evolution of the maximum, minimum and amplitude of the interfacial displacement for a cnoidal wave with a contraction (a) and expansion (b) defect initial condition. Black solid, red dashed and blue dot-dashed lines correspond to 5, 7 and 9 peaks in the domain, respectively.
Next, we use the cnoidal waves with expansion and contraction defects as the initial conditions in the presence of rotation. The effect of rotation on the cnoidal wave with a contraction and expansion defects is shown in the subsequent rows of figure 9, in the left and right columns, respectively, for the interfacial displacement (figure 9 c,d), shear in the direction of wave propagation (figure 9 e,f) and shear in the orthogonal direction (figure 9 g,h). Here, we see the formation of a strong burst of the interfacial waves and shear in the direction of wave propagation. In the moving reference frame the burst propagates to the left. Qualitatively, this is similar to the behaviour observed in the previous section, where numerical runs were initiated using the breather on a cnoidal wave initial condition, but the important difference is that this effect is solely due to rotation and no bursts of any kind are present in the absence of rotation.
Figure 10 shows the time evolution of the maximum, minimum and amplitude of the interfacial displacement for both cases of a contraction (figure 10
a) and expansion (figure 10
b) periodicity defect. We show that the observed effects of rotation are structurally stable with regard to the size of the computational domain. We experimented with 5, 7 and 9 major peaks in the domain. In the contraction defect case, the computational domains are
$2L= 37.42$
(5 peaks),
$2L = 55.02$
(7 peaks),
$2L = 72.86$
(9 peaks). In the expansion defect case, the computational domains are
$2L= 45.44$
(5 peaks),
$2L = 63.04$
(7 peaks),
$2L = 80.64$
(9 peaks). The remaining numerical parameters are the same as before. The burst forms soon after the initiation of the numerical runs, and the wave continues to grow for a long time after that. The close-up view of the large wave visible in the free-surface elevation is shown in figure 11 for the time around
$T=25$
(contraction defect, figure 11
a), with approximately 49 % increase in the amplitude compared with the initial condition and around
$T=32$
(expansion defect, figure 11
b), with an approximately
$31\,\%$
increase in the amplitude compared with the initial condition. Hence, we conclude that, under the effect of rotation, both contraction and expansion defects present in the otherwise regular cnoidal wave can lead to the emergence of strong bursts of interfacial waves and shear currents in the direction of propagation of the cnoidal wave.
Figure 11. Close-up view from below of the large burst in the interfacial displacement for a cnoidal wave with a contraction (a) and expansion (b) defect around
$T = 25$
(a) and
$T = 32$
(b), respectively.
4.4. The effect of rotation on cnoidal waves with generic localised defects
Finally in this section, motivated by the recent numerical and experimental generation of breathers in a fluid conduit during the interaction of a cnoidal wave with a soliton by Maiden & Hoefer (Reference Maiden and Hoefer2016) and Mao et al. (Reference Mao, Chandramouli, Xu and Hoefer2023) and generalising this scenario, we initiate numerical runs for our problem with the initial condition which can be described as a cnoidal wave with a generic localised perturbation. Namely, the initial condition is given by the function
\begin{align} \eta ^{(0)}|_{T=0} &= \dfrac {6\beta _1}{\alpha _1}\left \{ u_2 + (u_3 - u_2) \text{ cn}^2\left [ \left(\zeta { + \frac {\alpha _1 \hat \eta ^{(0)}}{3} T} - v_c T\right) \sqrt {\dfrac {u_3 -u_1}{2}} ; m \right ] \right \}_{T=0} \nonumber \\&\quad + A_1 \text{sech}^2 \left [ A_2 (\zeta +\zeta _0) \right ]\!, \end{align}
and
$u_1 \lt u_2 \lt u_3$
are real,
$v_c = 2 \beta _1 (u_1 + u_2 + u_3)$
and
$m ={(u_3-u_2)}/{(u_3-u_1)}$
, where
$A_1$
,
$A_2$
and
$\zeta _0$
are arbitrary constants.
Figure 12. Numerical solution for a cnoidal wave with a generic localised defect initial condition (waves of depression, view from above (a,c,e,g) and below (b,d,f,h)): (a,b) interfacial displacement in the absence of rotation, (c,d) interfacial displacement under the effect of rotation, (e,f) / (g,h) shear in the direction of wave propagation / orthogonal direction, under the effect of rotation.
We experiment with both sets of physical parameters discussed at the beginning of § 4. For the first set, numerical solutions initiated with the initial condition with a generic localised defect are shown in figure 12. The computational domain is
$2L = 79.20$
, with the number of modes
$M=790$
, the spatial step
$\Delta x \approx 10^{-1}$
, the total simulation time
$T_{max} = 100$
and the temporal step
$\Delta T = 10^{-2}$
. The parameters used in the initial condition are
$u_1 = -10^{-3},\, u_2=0,\, u_3=3,\, A_1 = -0.80,\, A_2 = 1.00$
and
$\zeta _0 = -1.30$
.
The initial condition introduces a localised perturbation to the otherwise regular cnoidal wave. In our runs, in the absence of rotation (KdV regime), the introduction of such a localised defect leads to the formation of a pair of bright and dark breathers, as well as the clearly noticeable expansion defect, which is shown in the first row of figure 12, showing both the view from above (figure 12 a) and the view from below (figure 12 b). Qualitatively, this is similar to the waves generated in the interaction of a cnoidal wave with a solitary wave in a conduit in Mao et al. (Reference Mao, Chandramouli, Xu and Hoefer2023). In the moving reference frame, the bright breather propagates to the right, while the dark breather propagates to the left, and the expansion defect moves with the speed of the cnoidal wave. The effect of rotation on the evolution of the same initial condition is shown in the subsequent rows of the same figure, for the top (figure 12 c,e,g) and bottom (figure 12 d,f,h) views of the interfacial displacement (figure 12 c,d), shear in the direction of wave propagation (figure 12 e,f) and shear in the orthogonal direction (figure 12 g,h). There are strong bursts both in the interfacial displacement and shear in the direction of wave propagation, which can be associated with the previously considered types of ‘defects’ in the otherwise regular cnoidal wave. The signal in the shear in the orthogonal direction is again weak. Hence, we conclude that the mechanism of formation of the bursts of strong interfacial displacements and shear in the direction of wave propagation in these simulations can be interpreted as formation of KdV-type breathers and the expansion defect on a cnoidal wave background, with the subsequent effect of rotation.
In order to test our theoretical framework further, we also modelled the effect of rotation on internal waves of elevation, using our second set of parameters discussed at the beginning of § 4. In figure 13, the computational domain is
$2L = 79.20$
, with the number of modes
$M=790$
, the spatial step
$\Delta x \approx 10^{-1}$
, the total simulation time
$T_{max} = 100$
and the temporal step
$\Delta T = 10^{-2}$
. The parameters used in the initial condition are
$u_1 = -10^{-3},\, u_2=0,\, u_3=3,\, A_1 = 0.80,\, A_2 = 1.00$
and
$\zeta _0 = -1.30$
. Here, the pycnocline is closer to the bottom, the upper layer depth takes
$80\,\%$
of the total depth, whereas in previous cases it was
$25\,\%$
. The results are similar to our previous simulations for internal waves of depression, but here we see, in the absence of rotation (i.e. in the KdV regime, see figure 13
a,b), the formation of a pair of bright and dark breathers, as before, but also the noticeable formation of both expansion and contraction defects. The comparison of the evolution of maxima, minima and amplitude of the waves shown in figures 12 and 13, under the effect of rotation, is given in figure 14, illustrating the growth of the amplitude of the waves and further supporting our interpretation of the observed (in all cases) formation of bursts of strong interfacial displacements and shear in the direction of wave propagation. For both sets, the computational domains are
$2L= 44.00$
(5 peaks),
$2L = 61.60$
(7 peaks) and
$2L = 79.20$
(9 peaks). The remaining numerical parameters are the same as before. The results are structurally stable. The large waves formed at the early stage of our computations are shown in figure 15. At around
$T=19$
, the first and second sets have approximately
$33\,\%$
and
$36\,\%$
increases, respectively, in the amplitudes relative to the initial conditions.
Figure 13. Numerical solution for a cnoidal wave with a generic localised defect initial condition (waves of elevation, view from above (a,c,e,g) and below (b,d,f,h)): (a,b) interfacial displacement in the absence of rotation, (c,d) interfacial displacement under the effect of rotation, (e,f) / (g,h) shear in the direction of wave propagation / orthogonal direction, under the effect of rotation.
Figure 14. The effect of rotation on the evolution of the maximum, minimum and amplitude of the interfacial displacement for a cnoidal wave with a generic localised defect initial condition: waves of depression (a) and waves of elevation (b). Black solid, red dashed and blue dot-dashed lines correspond to 5, 7 and 9 peaks in the domain, respectively.
Figure 15. Close-up view from below (a) and above (b) the large bursts in the interfacial displacement for a cnoidal wave with a generic localised defect initial condition around
$T = 19$
for the waves of depression (a) and waves of elevation (b).
5. Discussion
In this paper, we addressed several issues related to the modelling of internal waves in the Ostrovsky regime, i.e. in the KdV regime taking account of rotation. We chose Helfrich’s
$f$
-plane extension (MMCC-f) of the two-layer MMCC model as our parent system. The MMCC-f model has similar properties to the full Euler equations with rotation with respect to construction of the weakly nonlinear solution, but the technical details related to the derivation are more manageable. In our derivation, we represented all field variables as the sums of time-dependent mean values and both spatially and time-dependent deviations from these evolving means. We considered periodic solutions, and mean-field equations were obtained by averaging the equations over the period of the problem. In contrast to our previous research within the scope of the Boussinesq-type equations, the resulting equations for the mean fields turned out to be coupled to the equations for the deviations, which presented a significant new challenge. We managed to construct a large class of solutions for uni-directional waves by introducing two slow-time variables, and simultaneously constructing asymptotic expansions in powers of the square root of the amplitude parameter (as opposed to the traditional derivation of the Ostrovsky equation with just one slow time and in powers of the amplitude parameter). Since the resulting reduced model, the Ostrovsky equation, has been derived for the zero-mean deviations by construction, the so-called ‘zero-mean contradiction’ meaning the existence of the zero-mean constraints on all field variables of the parent system has been by-passed. Importantly, the Ostrovsky equation is decoupled from the variation of mean fields, with the latter acting as additive inertial oscillations in the shear variables as well as additive oscillations of the same frequency in the phase of both the interfacial displacement and shear variables. Non-zero barotropic transport in the propagation direction contributed to the shift in the phase of the interfacial displacement and shear in the propagation direction, and to the transverse barotropic transport.
Next, we used the constructed weakly nonlinear solution combined with extensive numerical modelling using the Ostrovsky equation in order to study the effect of rotation on the evolution of internal waves with initial conditions in the form of KdV cnoidal waves, but with various local defects. We built the phenomenological framework by considering pure bright and dark breathers on cnoidal wave initial conditions and cnoidal waves with contraction and expansion periodicity defects. The latter defects were introduced ‘by hand’, by modifying the cnoidal wave solution near the trough in between the two peaks, and we showed that such functions satisfy all (infinitely many) conservation laws of the KdV equation, where the integration is understood as the sum over the natural subintervals (see Appendix A). Moreover, the cnoidal waves with an expansion defect also have a continuous first derivative, satisfying the so-called ‘corner condition’ necessary for it to be a non-smooth extremal of the relevant variational problem. It must be noted that in our numerical runs, in the absence of rotation (i.e. in the KdV regime), initial conditions with both the contraction and expansion defects led to long-lived states, and a cnoidal wave with an expansion defect behaved very closely to a travelling wave of the KdV equation. The important difference between the defects represented by dislocation (bright and dark breather) and periodicity (contraction and expansion) perturbations of the cnoidal wave is that the second type does not lead to the formation of bursts in the absence of rotation. Qualitatively, both types of defects lead to the formation of bursts of large interfacial displacements and shear in the direction of wave propagation, under the effect of rotation.
Finally, we considered initial conditions in the form of cnoidal waves with generic localised perturbations. Our modelling has shown that, in the absence of rotation (in the KdV regime), such initial conditions typically produce a pair of a fast-moving bright and slow-moving dark breathers and expansion and contraction periodicity defects, moving with the speed of the cnoidal wave. Under the effect of rotation, the splitting of a generic localised perturbation into these ‘basic’ defects is followed by formation of several bursts of interfacial displacements and shear in the direction of wave propagation, which can be associated with the effect of rotation on breathers and the periodicity defects. In all simulations the large bursts propagated to the left in our moving reference frame, i.e. slower than the speed of our frame. Although the detailed description of the observed interfacial bursts goes beyond the scope of our present study, we note that they are qualitatively similar to the unsteady wavepackets observed in the previous studies of the effect of rotation on the KdV solitary waves (Grimshaw & Helfrich Reference Grimshaw and Helfrich2008) and slowly modulated waves (Whitfield & Johnson Reference Whitfield and Johnson2017), but they appear tens of times faster (comparing the evolution for the same physical parameters and different initial conditions). We also note that we compared the evolution of a cnoidal wave with a localised perturbation in the simulations with periodic boundary conditions with the simulations in a bigger domain and with zero boundary conditions and sponge layers near the boundaries. In the latter case the initial condition coincided with the former one in the central part of the domain, and then was gradually reduced to zero towards the boundaries. Our results have shown good agreement between these two simulations in the central part of the latter simulation. Hence, for the problems of that type, simulations with periodic boundary conditions in a smaller domain can be used instead of simulations in the large domain.
Figure 16. Rogue wave generation: three- and two-dimensional views of the interfacial displacement from above (a) and initial condition at
$T=0$
vs interfacial displacement at
$T = 8.8$
(b).
In this study we did not aim to systematically investigate the probability of generation of rogue waves due to this scenario (in the sense of the classical definition of rogue waves as waves with an amplitude more than twice the significant wave height), but it is worth mentioning the generation of the rogue wave shown in figure 16 using the same physical parameters as in figure 13 and initial condition (4.16) with
$u_1 = -10^{-3},\, u_2=0, u_3=3,\, A_1 = 0.60,\, A_2 = 0.60$
and
$\zeta _0 = -2.60$
. The computational domain is
$2L=272.81$
with the number of modes
$M=2726$
, the spatial step
$\Delta x \approx 10^{-1}$
, the total simulation time
$T_{max} = 100$
, with the temporal step
$\Delta T = 10^{-2}$
. This could be an interesting direction of further research.
We hope that our study sheds light both on the typical structure of nearly cnoidal waves (with defects) observed in oceanic observations (see, for example, figure 1 in the Introduction), and on the role of rotation in the formation of bursts of large interfacial displacements and shear in the direction of wave propagation. Further developments could concern bi-directional propagation and effects related to barotropic transport, which were discussed in the construction of the weakly nonlinear solution, but were left out of scope of the subsequent modelling. Another issue which was left out of scope of the present study was the effect of the depth-dependent shear currents (see Hooper, Khusnutdinova & Grimshaw Reference Hooper, Khusnutdinova and Grimshaw2021, Tseluiko et al. Reference Tseluiko, Alharthi, Barros and Khusnutdinova2023 for the recent developments related to the presence of a long-wave instability, in the absence of rotation). As a by-product of our study we constructed some curious generalised travelling waves of the KdV equation, whose smoothed counterparts naturally emerged in the evolution of initial conditions in the form of cnoidal waves with local perturbations in the absence of rotation (KdV regime). It would be interesting to investigate whether similar solutions can be found in other integrable models, as well as the non-integrable strongly nonlinear parent models, and what is their meaning and role in the relevant physical contexts.
Acknowledgements
K. Khusnutdinova is grateful to R. Barros, R. Camassa, E. Ferapontov, S. Gavrilyuk, B. Kruglikov, D. Pelinovsky and V. Shrira for useful discussions, and to Z. Makridin and N. Makarenko for sending the paper Maltseva (Reference Maltseva1989). We also thank the referees for their helpful comments, and J.-C. Saut for the references to Kato (Reference Kato1983), Kruzhkov & Faminskii (Reference Kruzhkov and Faminskii1984) and Bourgain (Reference Bourgain1993). K. Nirunwiroj would like to acknowledge the Royal Thai Government Scholarship supporting his postgraduate studies at Loughborough University. K. Khusnutdinova would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme ‘Emergent phenomena in nonlinear dispersive waves’, where partial work on this paper was undertaken.
Funding
This work was supported by EPSRC grant EP/R014604/1.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Generalised (shock-like) travelling waves of the KdV equation
The cnoidal waves with expansion/contraction defects have a jump in the second/first derivative at some points within the domain. Our analysis below follows the recent line of research developed in the papers by Gavrilyuk et al. (Reference Gavrilyuk, Nkonga, Shyue and Truskinovsky2020), Gavrilyuk & Shyue (Reference Gavrilyuk and Shyue2022) and Gavrilyuk et al. (Reference Gavrilyuk, Nkonga and Shyue2024). The key idea is to treat the points of discontinuity similarly to shocks in the theory of hyperbolic conservation laws, i.e. to split the solution into several parts where the solution and all derivatives are smooth, and to match them at the points of discontinuity in such a way that the ‘generalised Rankine–Hugoniot (gRH)’ conditions following from the dispersive conservation laws are satisfied. In the previous work developed within the scope of nonintegrable equations, e.g. the Benjamin–Bona–Mahoney equation considered by Gavrilyuk & Shyue (Reference Gavrilyuk and Shyue2022), the conservation laws have led to some non-trivial gRH conditions, and numerical experiments have supported the importance of these conditions. Our observation reported in this appendix is that the cnoidal waves of the KdV equation with the expansion or contraction defects introduced in our paper do not lead to any non-trivial jump conditions – all conservation laws are satisfied, in the same sense as in Gavrilyuk & Shyue (Reference Gavrilyuk and Shyue2022). Moreover, the cnoidal wave with an expansion defect satisfies the Weirstrass–Erdmann corner condition and is a broken extremal of the associated variational problem.
The KdV equation
\begin{equation} u_t - 6 u u_x + u_{xxx} = 0 ,\end{equation}
can be represented in the form of a conservation law
\begin{equation} \frac {\partial u}{\partial t}+\frac {\partial }{\partial x}\big ({-}3u^2+u_{xx}\big )=0. \end{equation}
Equation (A2) implies conservation of ‘mass’
\begin{equation} \frac {\rm d}{{\rm d} t}\int ^{L}_{-L} u {\rm d}x=0, \end{equation}
provided the function
$u(x,t)$
either vanishes together with its spatial derivatives as
$L \to \infty$
, or it is spatially periodic on the interval
$[-L, L]$
. The latter case is considered in our paper.
One also has conservation laws for ‘momentum’
\begin{equation} \frac {\partial }{\partial t} \frac {u^2}{2}+\frac {\partial }{\partial x}\left (uu_{xx}-\frac {1}{2}u_x^2 -2u^3 \right )=0 ,\end{equation}
and ‘energy’
\begin{equation} \frac {\partial }{\partial t} \left ( u^3+\frac {1}{2}u^2_x \right )+\frac {\partial }{\partial x}\left ({-}\frac {9}{2}u^4+3u^2u_{xx}-6uu_x^2 + u_x u_{xxx}-\frac {1}{2}u_{xx}^2\right )=0. \end{equation}
It is well known that the KdV equation has infinitely many polynomial conservation laws, which can be proven using a generating function (Miura, Gardner & Kruskal Reference Miura, Gardner and Kruskal1968). We reproduce the main part of the proof below, since we need some formulae from it for our subsequent considerations. The Gardner transformation
\begin{equation} u=w+\epsilon w_x+ \epsilon ^2 w^2, \end{equation}
where
$\epsilon$
is the grading parameter, yields
\begin{equation} u_t-6uu_x+u_{xxx}=\big(1+2\epsilon ^2 w+\epsilon \partial _x\big)\big[w_t+\big(w_{xx} - 3w^2-2\epsilon ^2 w^3\big)_x\big]=0. \end{equation}
Therefore, if
$w(x,t)$
satisfies the conservation law
\begin{equation} w_t+(w_{xx} - 3w^2-2\epsilon ^2 w^3)_x=0, \end{equation}
then
$u(x,t)$
satisfies the KdV equation. Next, one represents
$w$
in the form of an asymptotic expansion in powers of
$\epsilon$
\begin{equation} w = \sum _{n=0}^{\infty }\epsilon ^n w_n. \end{equation}
Collecting the coefficients of equal powers of
$\epsilon$
in (A6) gives us
\begin{align} &\displaystyle w_0=u \, , \qquad w_1=-w_{0x} =-u_x \, , \qquad w_2=-w_{1x} -w_0^2 =u_{xx}-u^2, \end{align}
\begin{align} &\qquad\qquad\quad \displaystyle w_n = - w_{n-1, x} - \sum _{k=0}^{n-2} w_k w_{n- 2 - k} \quad \mbox{for} \quad n \geqslant 2. \end{align}
Now, substituting expansion (A9) into (A8) one obtains an infinite series of the KdV conservation laws as the coefficients at even powers of
$\epsilon$
. The coefficients of odd powers are exact differentials (referred to as ‘trivial conservation laws’, see Miura et al. Reference Miura, Gardner and Kruskal1968, for details). The integrals
\begin{align} I_n=\int ^{L}_{-L}w_{2n} {\rm d}x, \quad n=0,\,1,\,2,\, \ldots ,\end{align}
are called ‘Kruskal integrals’. In particular, the first three integrals take the form
\begin{align} I_1=\int _{-L}^L u {\rm d}x, \quad I_2 = \int _{-L}^L \frac 12 u^2 {\rm d}x, \quad I_3= \int _{-L}^L \left (u^3 + \frac 12 u_x^2 \right ) {\rm d}x. \end{align}
Next, to prove that cnoidal waves with our contraction and expansion defects satisfy all conservation laws of the KdV equation, we also need a recurrence formula for the fluxes. Substituting expansion (A9) into (A8) we obtain
\begin{eqnarray} w_{xx} - 3 w^2 + 2 \epsilon ^2 w^3 = \sum _{n=0}^\infty \varepsilon ^n f_n, \end{eqnarray}
where
\begin{align} &\displaystyle f_0 = w_{0, xx} - 2 w_0^2, \quad f_1 = w_{1, xx} - 6 w_0 w_1, \quad f_2 = w_{2, xx} - 3 w_1^2 - 6 { w_0} w_2 + 2 w_0^3, \end{align}
\begin{align} &\qquad\quad \displaystyle f_n = w_{n, xx} - 3 \sum _{k=0}^n w_k w_{n-k} + 2 \sum _{k+l+m = n-2} w_k w_l w_m \quad \mbox{for} \quad n \geqslant 2 ,\\[9pt] \nonumber \end{align}
with the last summation being over all triples of integers
$(k, l, m)$
such that
$k+l+m = n-2$
. The Kruskal integrals and the recurrence relations for the fluxes are invariant with respect to the Galilean transformations.
Let us now consider a cnoidal wave with an expansion defect, as shown in the figure 17(a). This wave is constructed by cutting at a trough of the cnoidal wave between two peaks and symmetrically adding a segment of a straight line at the corresponding level. In the original
$(x,t)$
frame, the entire wave moves with the constant speed of the cnoidal wave. We now consider it in the reference frame moving with this constant speed, then the wave is stationary. We place the origin in the middle of the defect. Then, the function describing the wave profile is even, see the figure 17(a). Let us denote the period of the cnoidal wave by
$\ell$
. If we consider
$m$
peaks in the domain, then without the defect
$L=m \ell /2$
. With the defect, given that the length of the inserted symmetric interval is
$d$
, then
$L = (m \ell +d)/2$
. The resulting function has a corner at two points in the periodic domain
$[-L, L]$
, and the graph is symmetric about the vertical axis passing through the middle of the interval which has been assumed to be
$x=0$
, with the function
$w_0 = u$
being even in that frame. From the recurrence relation (A11), we observe that functions
$w_{2n}$
are even and
$w_{2n+1}$
are odd. Let the
$x$
-coordinates of the left and right endpoints of the line segment be denoted by
$x_1$
and
$x_2$
, respectively. With our choice of the reference frame we have
$x_1 = - x_2 = - x_d$
. Consequently, treating the points of discontinuity similarly to shocks in the theory of hyperbolic equations (see Gavrilyuk et al. Reference Gavrilyuk, Nkonga, Shyue and Truskinovsky2020, Gavrilyuk & Shyue Reference Gavrilyuk and Shyue2022, Gavrilyuk et al. Reference Gavrilyuk, Nkonga and Shyue2024), we represent the integral as the sum over three subintervals, and obtain
\begin{align} \frac {\rm d}{{\rm d}t} I_n &= - \int _{-L}^{-x_d} {f_{2n}}_ x {\rm d}x - \int _{-x_d}^{x_d} {f_{2n}}_x{\rm d}x - \int _{x_d}^{L} {f_{2n}}_x {\rm d}x \end{align}
\begin{align} &= - \left [f_{2n}\right ]_{-L}^ {-x_d} - \left [f_{2n} \right ]_{-x_d}^ {x_d} - \left [f_{2n}\right ]_{x_d}^ {L} = 0, \end{align}
by periodicity, the constant value of
$f_{2n}$
on the interval
$[-x_d, x_d]$
, and since the flux
$f_{2n}$
is an even function. We note that this can be generalised by first extracting a small part close to the trough, symmetrically inserting a segment of a straight line and gluing the remaining parts to it. The resulting functions will also satisfy all conservation laws of the KdV equation. Note that unlike the examples constructed in Gavrilyuk & Shyue (Reference Gavrilyuk and Shyue2022) in the context of the Benjamin–Bona–Mahoney equation, there appear no non-trivial gRH conditions.
Figure 17. Schematic of construction of the generalised travelling waves of the KdV equation in the form of a long-lived cnoidal wave with an expansion and contraction periodicity defects (a and b, respectively), and a short-lived soliton with similar amplitude defects (c and d, respectively). The waves have negative polarity: the peaks have negative amplitude.
Similarly we can consider a cnoidal wave with a contraction defect, as shown in the figure 17(b), again in the reference frame where the wave is stationary. This wave is constructed by extracting a symmetric interval around a trough between two neighbouring peaks and gluing the remaining parts together. Let us again denote the period of the cnoidal wave by
$\ell$
and consider
$m$
peaks in the domain, then without the defect
$L=m \ell /2$
. With the defect, given that the length of the extracted symmetric interval is
$d$
, then
$L = (m \ell -d)/2$
. The resulting function has a corner at one point in the periodic domain
$[-L, L]$
, and the graph is symmetric about the vertical axis passing through that point. Let us denote this point by
$x_0$
. We can then associate the origin of the moving reference frame with
$x_0$
, i.e. assume
$x_0 = 0$
, and the function
$w_0 = u$
will be even in that frame. From the recurrence relation (A11), we observe that functions
$w_{2n}$
are even and
$w_{2n+1}$
are odd. Recalling that non-trivial conservation laws correspond to even values
$2n$
, we deduce from the recurrence relation (A16) that the corresponding fluxes
$f_{2n}$
can only contain terms representing even functions. Therefore,
\begin{eqnarray} \frac {\rm d}{{\rm d}t} I_n &=& - \int _{-L}^{0}{ f_{2n}}_x {\rm d}x - \int _{0}^{L}{ f_{2n}}_x {\rm d}x \end{eqnarray}
\begin{eqnarray} &=& - \left [f_{2n} \right ]_{-L}^ {0} - \left [f_{2n} \right ]_{0}^ {L} = 0 ,\end{eqnarray}
by periodicity and since the flux
$f_{2n}$
is an even function.
Another way to prove the above relies on using that all non-trivial fluxes
$f_{2n}$
are polynomial functions of even weight with respect to the scaling symmetries of the KdV equation. This also implies that the fluxes
$f_{2n}$
are even functions, and the rest follows from that.
Next, the KdV equation (A1) can be written in Lagrangian form
\begin{equation} \delta \int {\mathcal L} {\rm d}x {\rm d}t = 0, \quad \mbox{where} \quad {\mathcal L} = \frac 12 \phi _t \phi _x - \phi _x^3 - \frac 12 \phi _{xx}^2, \quad u = \phi _x. \end{equation}
The cnoidal waves with expansion defects can be viewed as generalised (shock-like) travelling waves (see Gavrilyuk & Shyue Reference Gavrilyuk and Shyue2022 for a relevant discussion in the context of the Benjamin–Bona–Mahoney equation). Then Lagrangian density (A21) becomes a function of
$u$
and
$u_x$
:
${\mathcal L} = {\mathcal L}(u, u_x)$
. The Weirstrass–Erdmann corner conditions for broken extremals require the continuity of
\begin{equation} \frac {\partial \mathcal L}{\partial u_{x}} \quad \mbox{and} \quad {\mathcal L} - u_{x} \frac {\partial {\mathcal L}}{\partial u_{x}} ,\end{equation}
at each junction (e.g. Fox Reference Fox1954) , which leads to the requirement of continuity of
$u_x$
, since
$u$
is continuous. This condition is satisfied exactly for a cnoidal wave with an expansion defect, and approximately for a cnoidal wave with a contraction defect, provided we cut close to extremum. The smoothed counterparts of cnoidal waves with both expansion and contraction defects were long lived in our numerical simulations, with the former function generating states behaving very closely to a travelling wave.
Having observed that the constructed cnoidal waves with expansion defects generate long-lived states close to travelling waves, it is also temping to try to construct similar solutions by cutting around the peak rather than the trough. In contrast to the previous cases, this can be done for a single peak as well, starting from the exact soliton solution, as shown in figure 17(c,d). However, in numerical runs, the smoothed counterparts of these solutions turned out to be unstable and short lived, fissioning and giving rise to the usual solitons, in full agreement with the inverse scattering transform predictions (Gardner et al. Reference Gardner, Greene, Kruskal and Miura1967, see also Drazin & Johnson Reference Drazin and Johnson1989) for the respective smoothed initial conditions. Similar behaviour is observed in the cnoidal wave with an insertion at a peak rather than though. Such smoothed initial conditions do not produce long-lived states close to travelling waves in pseudospectral simulations, they are unstable.
Finally, let us also make some remarks concerning a possible weak formulation. Consider a
$2L$
-periodic regular travelling wave described by a function
$u = u(\xi ), \xi = x - c t$
of the KdV equation (A1)
\begin{equation} \big({-}c u - 3 u^2 + u_{\xi \xi }\big)_{\xi } = 0. \end{equation}
Let us multiply (A23) by a test function
$\phi$
, which is assumed to be smooth and
$2L$
-periodic, and integrate by parts. Then,
\begin{equation} \int _{-L}^{L} \big({-}c u - 3 u^2 + u_{\xi \xi }\big)_{\xi } \phi {\rm d} \xi = - \int _{-L}^L \big({-}c u + u^2 + u_{\xi \xi }\big) \phi _{\xi } {\rm d} \xi = 0, \end{equation}
where we used periodicity of the functions. Consider a cnoidal wave with an expansion defect. It has a jump in the second derivative at two points in the periodic domain and does not satisfy the strong formulation (A23). However, it does satisfy (A24). Indeed, let us consider the wave in the reference frame moving with the speed
$c$
. The wave is stationary and even in this reference frame. Then, it is natural to require the test function
$\phi$
to be even too. Equation (A23) is invariant with respect to the Galilean transformation. The integral in the right-hand side of (A24) is well defined. Represent it as the sum of three integrals, over the intervals
$[-L, -x_d], [-x_d, x_d], [x_d, L]$
. In the first and third integrals,
$ (- c u + u^2 + u_{\xi \xi }) = \textrm { const}$
(cnoidal wave), and then the sum of these integrals is equal to zero using that
$\phi$
is even and periodic. In the middle integral,
$u = \textrm {const}$
(expansion defect), and then again the integral vanishes using that
$\phi$
is an even function. Hence, a cnoidal wave with an expansion defect satisfies a weak formulation, in the above sense.
If we relax the conditions on the test function, and require only the periodicity but do not require
$\phi$
to be an even function, then we can still make the above argument work by requiring that the length of the defect is commensurate with the period
$\ell$
of the cnoidal wave, i.e.
$ x_d = n \ell$
, where
$n$
is a natural number. We note that we experimented with both commensurate and non-commensurate lengths of the defect, and there was no significant difference in the states produced at the end of long pseudospectral runs.
The global well posedness of the Cauchy problem with periodic boundary conditions in
$L^2$
, including uniqueness and continuous dependence with respect to the initial data, was proven by Bourgain (Reference Bourgain1993) (the global existence of the weak
$L^2$
solutions on the infinite line was proven by Kato Reference Kato1983, Kruzhkov & Faminskii Reference Kruzhkov and Faminskii1984). We believe that our considerations provide strong evidence that cnoidal waves with expansion defects considered in this paper are weak solutions of the KdV equation. The states close to them naturally emerged in numerical modelling of the evolution of initial conditions given by cnoidal waves with localised perturbations (see, for example, the figures 12(a,b) and 13(a,b), showing the evolution in the absence of rotation, i.e. in the KdV regime).
Appendix B. Numerical scheme for the Ostrovsky equation
The Ostrovsky equation
\begin{align} \left ( u_{T} + \alpha _1 u u_{\zeta } + \beta _1 u_{\zeta \zeta \zeta } \right )_\zeta = \gamma _1 u ,\end{align}
is numerically solved by pseudospectral methods (see Fornberg Reference Fornberg1996; Trefethen Reference Trefethen2000), which usually offer the highest accuracy and computational efficiency for smooth data on periodic domains. By implementing the fast Fourier transform algorithm for spatial derivatives and the fourth-order Runge–Kutta scheme for time stepping, the method yields accurate approximation of the solution to the differential equation.
The spatial domain is discretised by
$M$
equidistant points with spacing
$\Delta x = 2\pi /M$
. Then, the discrete Fourier transform of the equation with respect to
$\zeta$
gives
\begin{align} \hat {u}_T - i\left (k^3 \beta _1 -\dfrac {\gamma _1}{k}\right )\hat {u} = - \dfrac {ik\alpha _1}{2} \hat {(u^2)}, \end{align}
where
$k$
is the scaled wavenumber. We use the fourth-order Runge–Kutta scheme for temporal integration. By the integrating factor method of Kassam and Trefethen (e.g. Trefethen Reference Trefethen2000), we multiply the equation by
$K=\exp [-i(k^3\beta _1-({\gamma _1}/{k}))T]$
to obtain
\begin{align} \hat {U}_T = -\dfrac {ik\alpha _1}{2} K \mathcal{F}\left\{ \left( \mathcal{F}^{-1}\left [ \dfrac {\hat {U}}{K} \right ]\right)^2 \right\}, \end{align}
where
$\hat {U} = \exp [-i(k^3\beta _1-({\gamma _1}/{k}))T] \, \hat {u} = K \hat {u}\,$
and
$\mathcal{F}$
is the Fourier transform. Discretising the time domain as
$T_n=n \Delta T$
, and introducing the function
\begin{align} E = \exp \left[\frac {i}{2}\left(k^3\beta _1-\frac {\gamma _1}{k}\right)\Delta T\right], \end{align}
the optimised Runge–Kutta time stepping has the form
\begin{align} \hat {u}^{(n+1)} = E^2 \hat {u}^{(n)} + \dfrac {1}{6}\left[E^2 k_1 + 2E (k_2 + k_3) + k_4\right], \end{align}
where
\begin{align} k_1 &= -\dfrac {ik\alpha _1}{2} \Delta T \; \mathcal{F}\big\{ \big(\mathcal{F}^{-1} [ \hat {u}^{(n)} ] \big)^2\big\}, \nonumber\\k_2 &= -\dfrac {ik\alpha _1}{2} \Delta T \; \mathcal{F}\big\{ \big(\mathcal{F}^{-1} [ E (\hat {u}^{(n)} +k_1/2) ] \big)^2\big\},\nonumber\\k_3 &= -\dfrac {ik\alpha _1}{2} \Delta T \; \mathcal{F}\big\{ \big(\mathcal{F}^{-1} [ E \hat {u}^{(n)} +k_2/2]\big)^2\big\},\nonumber\\k_4 &= -\dfrac {ik\alpha _1}{2} \Delta T \; \mathcal{F}\big\{ \big(\mathcal{F}^{-1} [ E^2 \hat {u}^{(n)} + E k_3 ]\big)^2\big\}. \end{align}
Open access
