2. The 2-D Boussinesq–Peregrine system and reduced models

We consider surface waves in a homogeneous fluid within the scope of the non-dimensional 2-D Boussinesq–Peregrine system (Peregrine Reference Peregrine1967, see also the review Lannes Reference Lannes2020 and references therein for the discussion of asymptotically equivalent systems and rigorous justification), assuming a flat bottom:

(2.1) \begin{align} & \eta _t + {\nabla}\boldsymbol{\cdot }\{(1 + \epsilon \eta ) \boldsymbol{u} \} = 0, \end{align}

(2.2) \begin{align} & \boldsymbol{u}_t + \epsilon (\boldsymbol{u} \boldsymbol{\cdot }{\nabla}) \boldsymbol{u} + {\nabla}\eta - \frac{\delta ^2}{3} {\nabla}({\nabla}\boldsymbol{\cdot }\boldsymbol{u}_t ) = \boldsymbol{0}. \end{align}

Here $\boldsymbol{u} = (u,v)^{\text{T}}$ are the depth-averaged horizontal velocities, $\eta$ is the free-surface elevation, $\epsilon$ is the small-amplitude parameter and $\delta$ is the long-wavelength parameter.

The 2-D Boussinesq–Peregrine system conserves mass exactly, while momentum and energy are conserved asymptotically (the details of our calculations can be found in Appendix A). The integral forms for these conservation laws can be written as

(2.3) \begin{align} & \frac{\textrm{d} \mathcal{M}}{\textrm{d} t} = \frac{\textrm{d}}{\textrm{d} t} \iint _D \eta \, \textrm{d}x \textrm{d}y = 0, \end{align}

(2.4) \begin{align} & \frac{\textrm{d} \mathcal{P}}{\textrm{d} t} = \frac{\textrm{d}}{\textrm{d} t} \iint _D \left ( 1 + \epsilon \eta \right ) \boldsymbol{u} \, \textrm{d}x \textrm{d}y = \textit{O}(\epsilon ^2,\epsilon \delta ^2,\delta ^4), \end{align}

(2.5) \begin{align} & \frac{\textrm{d} \mathcal{E}}{\textrm{d} t} = \frac{\textrm{d}}{\textrm{d} t} \iint _D \frac{1}{2} \left ( |\boldsymbol{u}|^2 + \eta ^2 \right ) \, \textrm{d}x \textrm{d}y = \textit{O}(\epsilon ,\delta ^2), \end{align}

respectively, where $D= \{(x, y)| x_1 \leqslant x \leqslant x_2, y_1 \leqslant y \leqslant y_2\}$ is the doubly periodic rectangular computational domain. These results generalise the 1-D results obtained by Katsaounis et al. (Reference Katsaounis, Mitsotakis and Sadaka2020), Khorbatly & Kalisch (Reference Khorbatly and Kalisch2024). In order to match the derivation of the reduced amplitude equations we choose the balance $\epsilon = \delta ^2$ . We therefore expect mass to be conserved in numerical simulations to machine precision, while the deviations in momentum and energy should scale as $\epsilon ^2$ and $\epsilon$ , respectively. This was indeed the case in our simulations (see the discussion in Appendix A). We numerically solve the system (2.1) and (2.2) using the pseudospectral method outlined in Appendix B.

Next, we derive the KdV $\theta$ and cKdV equations directly from the system (2.1) and (2.2), in the slow radius form. We apply the non-axisymmetric change of variables $(t,x,y) \rightarrow (t,r,\theta )$ to give

(2.6) \begin{align} & \eta _t + U_r + \frac{U}{r} + \frac{V_{\theta }}{r} + \epsilon \left [U\!\eta _r + \eta U_r + \frac{U\!\eta + V\!\eta _{\theta } + \eta V_{\theta }}{r} \right ] = 0, \end{align}

(2.7) \begin{align} & U_t + \eta _r + \frac{\epsilon }{3} \left [\frac{U_t + V_{t\theta }}{r^2} - \frac{U_{rt} - 3{\textit{VU}}_{\theta } + 3V^2 - V_{rt\theta }}{r} - U_{rrt} + 3{\textit{UU}}_r \right ] = 0, \end{align}

(2.8) \begin{align} & V_t + \frac{\eta _{\theta }}{r} + \frac{\epsilon }{3}\left [-\frac{U_{t\theta } + V_{t\theta \theta }}{r^2} - \frac{3UV + 3VV_{\theta } - U_{tr\theta }}{r} + 3{\textit{UV}}_r \right ] = 0, \end{align}

where the variables $U$ and $V$ are the cylindrically projected velocity components of $\boldsymbol{u}$ :

(2.9) \begin{eqnarray} u = U \cos \theta - V \sin \theta ,\quad v = U \sin \theta + V \cos \theta . \end{eqnarray}

To leading order, we obtain

(2.10) \begin{align} & \eta _{tt} - \eta _{rr} - \frac{\eta _r}{r} - \frac{\eta _{\theta \theta }}{r^2} = 0, \end{align}

(2.11) \begin{align} & U_t + \eta _r = 0, \end{align}

(2.12) \begin{align} & V_t + \frac{\eta _{\theta }}{r} = 0, \\[6pt] \nonumber \end{align}

where $\eta = \eta (\xi )$ , given $\xi = rk(\theta ) - t$ , is an exact solution of (2.10) if $k(\theta )$ has the form $k(\theta ) = a\cos \theta + b(a)\sin \theta$ and $a^2 + b(a)^2 = 1$ (coinciding with the general solution (1.14) for $u_0(z) = 0$ ), e.g. $a = \cos \varphi$ and $b = \sin \varphi$ , where $\varphi$ then has the meaning of the angle between the normal to the wavefront and the positive $x$ direction. In particular, if $\varphi = 0$ , i.e. waves propagate in the $x$ direction, then $\xi = r \cos \theta - t = x - t$ .

We then consider solutions in the far field such that $R = \epsilon rk(\theta ) \sim \textit{O}(1)$ and apply the change of variables $(t,r,\theta ) \rightarrow (R,\xi ,\theta )$ to obtain

(2.13) \begin{align} & \eta _{\xi } - kU_{\xi } - k'V_{\xi } - \epsilon \biggl [ kU_{\!R} + k\frac{U}{R} + k'V_R + k\frac{V_{\theta }}{R} + k(U\!\eta )_{\xi } + k'(V\!\eta )_{\xi } \biggr ]=O(\epsilon ^2), \end{align}

(2.14) \begin{align} & U_{\xi } - k\eta _{\xi } - \frac{\epsilon }{3}\bigl [3k\eta _R + 3k'{\textit{VU}}_{\xi } + kk' V_{\xi \xi \xi } + k^2 U_{\xi \xi \xi } + 3k{\textit{UU}}_{\xi } \bigr ]=O(\epsilon ^2), \end{align}

(2.15) \begin{align} & V_{\xi } - k'\eta _{\xi } - \frac{\epsilon }{3}\biggl [3k'\eta _R + 3k\frac{\eta _{\theta }}{R} + 3k'VV_{\xi } + kk' U_{\xi \xi \xi } + 3k{\textit{UV}}_{\xi } - k^{\prime 2} V_{\xi \xi \xi } \biggr ]=O(\epsilon ^2). \end{align}

We seek a solution of (2.13)–(2.15) in the form of an asymptotic multiple-scale expansion

(2.16) \begin{equation} \eta (R,\xi ,\theta ) = \eta ^{(0)}(R,\xi ,\theta ) + \epsilon \eta ^{(1)}(R,\xi ,\theta ) + \textit{O}(\epsilon ^2) \quad \mbox{as} \quad \epsilon \rightarrow 0, \end{equation}

and similarly for $U$ and $V$ . At leading order, $O(1)$ , we obtain

(2.17) \begin{align} & \eta ^{(0)} = \eta ^{(0)}(R,\xi ,\theta ), \end{align}

(2.18) \begin{align} & U^{(0)} = k\eta ^{(0)}, \end{align}

(2.19) \begin{align} & V^{(0)} = k'\eta ^{(0)}, \end{align}

assuming that any disturbances are generated only by the propagating waves. At next order, $O(\epsilon )$ , we obtain

(2.20) \begin{align} & \eta _{\xi }^{(1)} \!-\! kU_{\xi }^{(1)} \!-\! k'V_{\xi }^{(1)} = kU_{\!R}^{(0)} + k\frac{U^{(0)}}{R} + k'V_R^{(0)} + k\frac{V_{\theta }^{(0)}}{R} + k(U^{(0)}\eta ^{(0)})_{\xi } + k'(V^{(0)}\eta ^{(0)})_{\xi }, \end{align}

(2.21) \begin{align} & U_{\xi }^{(1)} = k\eta _{\xi }^{(1)} + k\eta _R^{(0)} + k'V^{(0)}U^{(0)}_{\xi } + \frac{1}{3}kk'V_{\xi \xi \xi }^{(0)} + \frac{1}{3}k^2U_{\xi \xi \xi }^{(0)} + kU^{(0)}U_{\xi }^{(0)}, \end{align}

(2.22) \begin{align} & V_{\xi }^{(1)} = k'\eta _{\xi }^{(1)} + k'\eta _R^{(0)} + k\frac{\eta _{\theta }^{(0)}}{R} + k'V^{(0)}V_{\xi }^{(0)} + \frac{1}{3}kk' U_{\xi \xi \xi }^{(0)} + k^{\prime 2}V_{\xi \xi \xi }^{(0)} + kU^{(0)}V_{\xi }^{(0)}, \end{align}

from which substituting (2.18) and (2.19), and (2.21) and (2.22) into (2.20), to eliminate $\eta ^{(1)}$ , $U^{(0)}$ , $U^{(1)}$ , $V^{(0)}$ and $V^{(1)}$ , gives the KdV $\theta$ equation for perturbed surface plane waves propagating at an angle $\varphi$ to the positive $x$ direction as

(2.23) \begin{align} 2\eta _R^{(0)} + 3 \eta ^{(0)}\eta _{\xi }^{(0)} + \frac{1}{3}\eta _{\xi \xi \xi }^{(0)} + 2kk' \frac{\eta _{\theta }^{(0)}}{R} = 0. \end{align}

The corresponding leading-order Cartesian velocities are calculated from (2.9) and (2.18), (2.19) to be

(2.24) \begin{equation} u^{(0)} = \cos \varphi \, \eta ^{(0)} \qquad \text{and} \qquad v^{(0)} = \sin \varphi \, \eta ^{(0)}. \end{equation}

If $x$ is chosen to be the direction of propagation then $\varphi = 0, \ k = \cos \theta$ , and (2.23) takes the form

(2.25) \begin{equation} 2\eta _R^{(0)} + 3 \eta ^{(0)}\eta _{\xi }^{(0)} + \frac{1}{3}\eta _{\xi \xi \xi }^{(0)} -2\sin \theta \cos \theta \frac{\eta _{\theta }^{(0)}}{R} = 0, \end{equation}

where $ R = \epsilon x$ , $\xi = x-t$ and the corresponding Cartesian velocities are $u^{(0)} = \eta ^{(0)}$ and $v^{(0)} = 0$ . This is a particular form of the KdV $\theta$ equation (3.3), derived directly from the 2-D Boussinesq–Peregrine system (2.1) and (2.2).

The cKdV equation can be derived using a simpler asymptotic expansion (see, for example, Miles Reference Miles1978; Johnson Reference Johnson1997; Sidorovas et al. Reference Sidorovas, Tseluiko, Choi and Khusnutdinova2024 and references therein). We apply the axisymmetric change of variables $(t,x,y) \rightarrow (t,r)$ to (2.1) and (2.2), which yields

(2.26) \begin{align} & \eta _t + U_r + \frac{U}{r} + \epsilon \left [U\!\eta _r + \eta U_r + \frac{U\!\eta }{r} \right ] = 0, \end{align}

(2.27) \begin{align} & U_t + \eta _r + \frac{\epsilon }{3} \left [\frac{U_t}{r^2} - \frac{U_{rt}}{r} - U_{rrt} + 3{\textit{UU}}_r \right ] = 0, \end{align}

where $U$ is the radial cylindrical velocity component of $\boldsymbol{u}$ . To leading order, we therefore obtain

(2.28) \begin{align} & \eta _{tt} - \eta _{rr} - \frac{\eta _r}{r} = 0, \end{align}

(2.29) \begin{align} & u_t + \eta _r = 0, \end{align}

where $\eta = \eta (\xi )$ , given $\xi = r - t$ , is an exact solution of (2.28). We then perform the change of variables $(t,r) \rightarrow (R,\xi )$ , where $R = \epsilon r$ is the slow radius variable and $\xi = r-t$ is the fast spatial variable in a moving coordinate frame, giving

(2.30) \begin{align} & \eta _{\xi } - U_{\xi } - \epsilon \left [U_{\!R} + \frac{U}{R} + U\!\eta _{\xi } + \eta U_{\xi } \right ] = O(\epsilon ^2), \end{align}

(2.31) \begin{align} & U_{\xi } - \eta _{\xi } - \frac{\epsilon }{3}\left [3\eta _R + 3{\textit{UU}}_{\xi } + U_{\xi \xi \xi } \right ] = O(\epsilon ^2). \end{align}

We now seek a solution of (2.30) and (2.31) in the form

(2.32) \begin{equation} \eta (R,\xi ) = \eta ^{(0)}(R,\xi ) + \epsilon \eta ^{(1)}(R,\xi ) + \textit{O}(\epsilon ^2) \quad \mbox{as} \quad \epsilon \rightarrow 0, \end{equation}

and similarly for the variable $U$ . We then consider solutions in the far field such that $R \sim \textit{O}(1)$ , which to leading order yields

(2.33) \begin{equation} \eta _{\xi }^{(0)} - U_{\xi }^{(0)} = 0, \end{equation}

giving that $U^{(0)} = \eta ^{(0)}$ , assuming that any disturbances are generated only by the propagating waves. At next order, we find that

(2.34) \begin{align} & \eta _{\xi }^{(1)} - U_{\xi }^{(1)} = U_{\!R}^{(0)} + \frac{U^{(0)}}{R} + 2\eta ^{(0)}\eta ^{(0)}_{\xi }, \end{align}

(2.35) \begin{align} & U_{\xi }^{(1)} = \eta _{\xi }^{(1)} + \eta _R^{(0)} + U^{(0)}U^{(0)}_{\xi } + \frac{1}{3}U_{\xi \xi \xi }^{(0)}, \end{align}

which after substituting (2.33) and (2.35) into (2.34) gives the axisymmetric cKdV equation for surface waves as

(2.36) \begin{equation} 2\eta _R^{(0)} + 3 \eta ^{(0)}\eta _{\xi }^{(0)} + \frac{1}{3}\eta _{\xi \xi \xi }^{(0)} + \frac{\eta ^{(0)}}{R} = 0. \end{equation}

Recalling that the depth-averaged Cartesian velocities are given by a linear combination of the polar velocity projections $U$ and $V$ by (2.9), this gives the Cartesian velocity components

(2.37) \begin{eqnarray} u^{(0)} = \frac{x}{\sqrt {x^2 + y^2}} \, \eta ^{(0)}, \quad v^{(0)} = \frac{y}{\sqrt {x^2 + y^2}} \, \eta ^{(0)}. \end{eqnarray}

Finally, the KPII equation is derived using the change of variables $(t,x,y) \rightarrow (T,\xi ,Y)$ , where $T = \epsilon t$ , $\xi = x - t$ and $Y = \sqrt {\epsilon } y$ (see, for example, Ablowitz & Segur Reference Ablowitz and Segur1979; Johnson Reference Johnson1997 and references therein). We also scale the Cartesian velocity in the $y$ direction as $v \rightarrow \sqrt {\epsilon } v$ to remain consistent with the coordinate scaling. Applying this change of variables to (2.1) and (2.2) gives

(2.38) \begin{align} & u_{\xi } - \eta _{\xi } + \epsilon \left (v_Y + \eta u_{\xi } + u\eta _{\xi } + \eta _T \right ) = O(\epsilon ^2), \end{align}

(2.39) \begin{align} & \eta _{\xi } - u_{\xi } + \epsilon \left (uu_{\xi } + \frac{1}{3}u_{\xi \xi \xi } + u_T \right ) = O(\epsilon ^2), \end{align}

(2.40) \begin{align} & \eta _Y - v_{\xi } + \epsilon \left (uv_{\xi } + \frac{1}{3} u_{\xi \xi Y} + v_T \right ) = O(\epsilon ^2). \end{align}

We now seek a solution of (2.38)–(2.40) in the form

(2.41) \begin{equation} \eta (T,\xi ,Y) = \eta ^{(0)}(T,\xi ,Y) + \epsilon \eta ^{(1)}(T,\xi ,Y) + \textit{O}(\epsilon ^2) \quad \mbox{as} \quad \epsilon \rightarrow 0, \end{equation}

and similarly for the variables $u$ and $v$ . We then consider solutions in the far field such that $T \sim \textit{O}(1)$ , and to leading order, $O(1)$ , we obtain

(2.42) \begin{align} u^{(0)}_{\xi} & = \eta ^{(0)}_{\xi }, \end{align}

(2.43) \begin{align} v^{(0)}_{\xi} &= \eta ^{(0)}_Y, \end{align}

such that $u^{(0)} = \eta ^{(0)}$ , assuming that any disturbances are generated only by the propagating waves. At the following order, $O(\epsilon )$ , we obtain

(2.44) \begin{align} \eta _{\xi }^{(1)} - u_{\xi }^{(1)} & = v_Y^{(0)} + \eta ^{(0)}u_{\xi }^{(0)} + u^{(0)}\eta _{\xi }^{(0)} + \eta _T^{(0)}, \end{align}

(2.45) \begin{align} u_{\xi }^{(1)} - \eta _{\xi }^{(1)} & = u^{(0)}u_{\xi }^{(0)} + \frac{1}{3} u_{\xi \xi \xi }^{(0)} + u_T^{(0)}, \end{align}

from which, taking the sum of (2.44) and (2.45) and substituting (2.42) gives

(2.46) \begin{equation} 2\eta ^{(0)}_T + 3\eta ^{(0)}\eta ^{(0)}_{\xi } + \frac{1}{3}\eta ^{(0)}_{\xi \xi \xi } + v^{(0)}_Y = 0. \end{equation}

Taking the $\xi$ derivative of (2.46) and substituting the $Y$ derivative of (2.43) gives the KPII equation for surface waves:

(2.47) \begin{equation} \left ( 2\eta ^{(0)}_T + 3\eta ^{(0)}\eta ^{(0)}_{\xi } + \frac{1}{3}\eta ^{(0)}_{\xi \xi \xi }\right )_{\xi } + \eta ^{(0)}_{\textit{YY}} = 0. \end{equation}

The slow time variable $T$ is more commonly used for the KPII equation and in the later sections makes comparisons to the solutions of the 2-D Boussinesq–Peregrine equations simpler. For the KdV $\theta$ and cKdV equations, the $R$ -variable version could be more natural (see, for example, the discussion concerning the cKdV equation in Sidorovas et al. Reference Sidorovas, Tseluiko, Choi and Khusnutdinova2024), but the slow time $T = \epsilon t$ version is easily recovered by the change of variables $(R, \xi , \theta ) \to (T, \xi , \theta )$ .

3. The KdV $\boldsymbol{\theta}$ equation and perturbed plane waves

In this section we first derive the KdV $\theta$ equation from the $(2 + 1)$ -dimensional amplitude equation (1.7) for a stratified fluid with a parallel shear flow, which was previously derived from the full set of Euler equations with the free-surface and flat-bottom boundary conditions by Khusnutdinova & Zhang (Reference Khusnutdinova and Zhang2016). Since the general solution of the equation for the function $k$ is associated with plane waves (Hooper et al. Reference Hooper, Khusnutdinova and Grimshaw2021), assuming that the waves exist, it will be in the form

(3.1) \begin{equation} k(\theta ; a) = a \cos \theta + b(a) \sin \theta , \end{equation}

with a problem-specific relation $b = b(a)$ . Indeed, the waves associated with the general solution (3.1) have the fast variable

(3.2) \begin{equation} \xi = r k(\theta ) - s t = ax + b(a) y - st, \end{equation}

with the constants $a$ and $b(a)$ having the meaning of the two wavenumbers of the plane wave propagating at an angle $\displaystyle \tan \varphi = b(a) / a$ to the flow (see figure 2). The wavefronts are tangent to the ring wave associated with the corresponding singular solution, obtained as a geometrical envelope of the general solution (see figure 3 for $u_0(z) = \gamma z, \ \gamma = \textrm {const}$ ).

Figure 2. Wavefront and wave vector of a plane wave propagating at an angle $\varphi$ to the direction of the shear flow $u_0(z)$ .

Figure 3. $(a)$ The general solution (1.14) for $a = 0.7$ (blue), $a=0.8$ (dashed blue) and $a=0.9$ (dashed dot blue) and its envelope (the singular solution, solid red) all for $u_0(z) = \gamma z$ , $\gamma = 0.2$ . $(b)$ The wavefronts corresponding to the general and singular solutions from $(a)$ described by $rk(\theta ) = 5$ .

Next, on the general solution (3.1), we obtain $k^2 + k^{\prime 2} = a^2 + b(a)^2$ and $F = -s + [u_0(z) - c] a$ . Hence, the modal equations (1.3)–(1.5) become $\theta$ independent, and therefore, the modal function is also $\theta$ independent, i.e. $\phi =\phi (z)$ . It then implies, using (1.9) and (1.17), that $\mu _4 = 0$ for any stratification and any parallel shear flow $u_0(z)$ . Hence, we obtain the amplitude equation, which has the form

(3.3) \begin{equation} \mu _1 A_R + \mu _2 {\textit{AA}}_{\xi } + \mu _3 A_{\xi \xi \xi } + \mu _5 \frac{A_\theta }{R} = 0, \end{equation}

and the coefficients of this equation are computed using the general solution (3.1) in (1.8) and (1.10). Generally, the coefficient $\mu _5$ has a non-trivial dependence on $\theta$ , while $\mu _1, \mu _2, \mu _3$ are some constants. For example, for surface waves in a homogeneous fluid over the linear shear flow $u_0(z) = \gamma z$ with $\gamma = \textrm {const}$ , we obtain

(3.4) \begin{align} \mu _1 & = \frac{2s - \gamma a}{s (s-\gamma a)^2}, \ \mu _2 = \frac{[a^2+b(a)^2]^2(3s-3s\gamma a + \gamma ^2 a^2)}{s^3(s-\gamma a)^3}, \ \mu _3 = \frac{[a^2+b(a)^2]^2}{3s^2}, \end{align}

(3.5) \begin{align} \mu _5 & = \frac{[a \cos \theta + b(a)\sin \theta ][2b(a)(s-\gamma a) \cos \theta - (2as + a^2 \gamma - b(a)^2 \gamma )\sin \theta ]}{s(s-\gamma a)^2}. \end{align}

When there is no shear flow, the coefficients reduce to $\mu _1 = 2$ , $\mu _2 = 3$ , $\mu _3 = 1/3$ and $\mu _5 = -2\sin \theta \cos \theta$ , resulting in the equation

(3.6) \begin{equation} A_R + \frac 32 {\textit{AA}}_{\xi } + \frac 1 6 A_{\xi \xi \xi } - \frac 12 \sin 2 \theta \frac{A_\theta }{R} = 0, \end{equation}

which coincides with (2.25) derived in § 2 directly from the 2-D Boussinesq–Peregrine system. If the wave has a plane wavefront then $A_{\theta } = 0$ (there is no change in the tangential direction along the wavefront), and the amplitude (3.3) (and hence (3.6)) has a reduction to the KdV equation (see Appendix C). However, the full equation (3.3), and its particular case (3.6), are more general. To the best of our knowledge, (3.3) and (3.6) have not been derived before and we will refer to both as the KdV $\theta$ equation.

In this section we consider plane waves subject to long transverse perturbations. As an example, we apply long transverse amplitude perturbations of finite strength to the line-soliton solution of the KdV $\theta$ equation (3.6), where we have changed the variables from $(R,\xi ,\theta )$ to $(T=\epsilon t,\xi ,\theta )$ for a more accurate comparison of numerical solutions of the reduced amplitude equation to those of the 2-D Boussinesq–Peregrine system. To leading order, the change of variable $R \rightarrow T$ does not change the form of the amplitude equation  (3.3), provided we map $\mu _1 \to \mu _1 /\! s$ and $\mu _5 \to \mu _5 /\! s$ :

(3.7) \begin{equation} \mu _1 A_T + s \mu _2 {\textit{AA}}_{\xi } + s \mu _3 A_{\xi \xi \xi } + \mu _5 \frac{A_{\theta }}{T} = 0. \end{equation}

The KdV $\theta$ equation (3.7) can be analysed analytically, using a mapping similar to that noted by Johnson (Reference Johnson1990) for ring waves. Indeed, for every value of $\theta$ , the change of variables

(3.8) \begin{equation} \tilde {\theta } = T \exp \left ( - \int \frac{\mu _1}{\mu _5} \, \textrm{d}\theta \right ) \end{equation}

maps the KdV $\theta$ equation (3.7) to the usual KdV equation

(3.9) \begin{equation} A_T + \frac{s \mu _2}{\mu _1} {\textit{AA}}_{\xi } + \frac{s \mu _3}{\mu _1} A_{\xi \xi \xi } = 0, \end{equation}

since (1.8) and (1.10) imply that the only $\theta$ -dependent coefficient in (3.7), and the slow $R$ version (3.3), is $\mu _5$ . However, the initial condition for the KdV equation (3.9) will have parametric dependence on $\tilde \theta$ , and hence, on the polar angle $\theta$ .

Here, we consider surface waves without a shear flow. Then $s = 1$ , and the KdV $\theta$ equation (3.7), for the waves propagating in the $x$ direction, takes the form

(3.10) \begin{eqnarray} &&A_T + \frac{3}{2} {\textit{AA}}_{\xi } + \frac{1}{6} A_{\xi \xi \xi } - \frac 12 \sin 2 \theta \frac{A_{\theta }}{T} = 0, \end{eqnarray}

where

(3.11) \begin{eqnarray} &&T = \epsilon t, \quad \xi = x-t, \quad \theta = \arctan \frac{y}{x}. \end{eqnarray}

A long transverse amplitude perturbation is added to the initial condition, at some $T = T_0$ , by multiplying the initial condition by a suitable function of $\theta$ . The unperturbed initial condition is the $\theta$ -independent exact line-soliton solution of the KdV $\theta$ equation (3.10), given by

(3.12) \begin{equation} A_0(T_0,\xi ,\theta ) = 2 \tilde v \operatorname {sech}^2 \left (\frac{1}{2} \sqrt {6 \tilde v} \left [\xi - \tilde v T_0 - \xi _0 \right ] \right )\!, \end{equation}

where $\tilde v \gt 0$ and $\xi _0$ are some constants. We take the perturbed initial condition to be in the form

(3.13) \begin{equation} A(T_0,\xi ,\theta ) = \left (1 + \alpha \operatorname {sech}^2 \beta \theta \right ) A_0(T_0,\xi ,\theta ), \end{equation}

although we could choose to multiply by any other function of $\theta$ . Here, the constants $\alpha$ and $\beta$ control the amplitude and width of the perturbation, $\tilde v$ is the speed of the soliton and $\xi _0$ is an initial phase shift. The multiplication of the original $\operatorname {sech}^2$ soliton by this $\theta$ -dependent function generates a smooth, localised perturbation in the central region of the wavefront.

The mapping (3.8) yields

(3.14) \begin{equation} \tilde \theta = T \tan \theta , \end{equation}

and the KdV equation (3.9) simplifies to become

(3.15) \begin{equation} A_T + \frac 32 {\textit{AA}}_{\xi } + \frac 16 A_{\xi \xi \xi } = 0, \end{equation}

which we now need to solve for the initial condition

(3.16) \begin{equation} A(T_0,\xi , \tilde \theta ) = \left [1 + \alpha \operatorname {sech}^2 \left (\beta \arctan \frac{\tilde \theta }{T_0}\right ) \right ]\! A_0 \!\left ( T_0,\xi , \arctan \frac{\tilde \theta }{T_0} \right )\!. \end{equation}

Next, to use the results of the IST (see Gardner et al. Reference Gardner, Green, Kruskal and Miura1967), we cast the KdV equation (3.15) into the canonical form

(3.17) \begin{equation} U_{\tau } - 6{\textit{UU}}_X + U_{\textit{XXX}} = 0, \end{equation}

via the scalings $\displaystyle U = - 3A/2, \ \tau = T/6, \ X = \xi .$ The wavefield is then defined by the spectrum of the Schrödinger equation

(3.18) \begin{equation} \varPsi _{\textit{XX}} + \left [\lambda - U(\tau _0,X;\tilde \theta ) \right ]\varPsi = 0, \end{equation}

where the potential is given by the perturbed initial condition (3.13), i.e.

(3.19) \begin{eqnarray} &&U(\tau _0,X;\tilde \theta ) = - \varLambda \operatorname {sech}^2 \left (\frac{X - 6 \tilde v \tau _0 - \xi _0}{L} \right )\!, \end{eqnarray}

where

(3.20) \begin{eqnarray} &&\Lambda = 3 \tilde v \Delta , \ \ L = \frac{2}{\sqrt {6 \tilde v}}, \ \ \varDelta (\tau _0;\tilde \theta ) = 1 + \alpha \operatorname {sech}^2\left ( \beta \arctan \frac{\tilde {\theta }}{6\tau _0} \right )\!, \end{eqnarray}

i.e. we have different initial conditions for different values of $\tilde \theta$ and, hence, for different values of the polar angle $\theta$ . The number of generated solitons depends on $\theta$ and is given by the greatest integer satisfying the inequality

(3.21) \begin{equation} N (\tau _0;\tilde \theta) \lt \frac{1}{2} \left [\big(1 + 4\varLambda L^2 \big)^{\frac{1}{2}} + 1 \right ] \end{equation}

(see, e.g., Landau & Lifshitz (Reference Landau and Lifshitz1959)). The discrete eigenvalues are given by

(3.22) \begin{eqnarray} \lambda = -k_n^2, \end{eqnarray}

where

(3.23) \begin{eqnarray} k_n (\tau _0;\tilde {\theta }) = \frac{1}{2L}\left [\big(1 + 4\varLambda L^2\big)^{\frac{1}{2}} - (2n-1) \right ] \gt 0, \quad n = 1,2,\ldots ,N. \end{eqnarray}

For positive $\alpha$ and $\tilde v$ values, the central perturbed region fissions into at least two solitons and a dispersive wave train. However, for small $\alpha$ values, the secondary soliton is small and for a very long time it blends with the radiation. The outer, unperturbed section of the line soliton propagates as a single curved soliton. We therefore expect a localised 2-D region of radiation behind a smaller soliton in the wake of the main perturbed soliton. The long-time asymptotics takes the form

(3.24) \begin{equation} U (\tau ,X;\tilde {\theta}) \simeq - \sum _{n=1}^N 2k_n^2 \operatorname {sech}^2 \big( k_n(X - 4k_n^2\tau - X_n)\big) + \text{radiation} \end{equation}

(e.g. Drazin Reference Drazin1983), where the phase shifts are given by

(3.25) \begin{eqnarray} X_1 = \frac{1}{2k_1} \ln \left (\frac{c_1}{2k_1} \right )\!, \quad X_n = \frac{1}{2k_n} \ln \left (\frac{c_n}{2k_n} \prod _{m=1}^{n-1} \left ( \frac{k_n - k_m}{k_n + k_m} \right )^2 \right )\!, \quad n \gt 1, \end{eqnarray}

and the constants $c_n$ are found as

(3.26) \begin{equation} c_n = \left (\int _{-\infty }^{\infty } \varPsi _n^2(x) \text{ d}x \right )^{-1}, \end{equation}

where $\varPsi _n(X)$ is the eigenfunction of (3.18) corresponding to the $n{\text{th}}$ eigenvalue, $-k_n^2$ :

(3.27) \begin{equation} \varPsi _n = \text{const}\left (1 - \tanh ^2 \frac{x}{L} \right )^{\dfrac{k_nL}{2}} {}_{2}F_{1}\left (1-n, 2k_nL + n,k_nL+1, \frac{1 - \tanh \dfrac{x}{L}}{2} \right )\!. \end{equation}

Here, ${}_{2}F_{1}(\boldsymbol{\cdot} )$ is the hypergeometric function and the constant should be chosen to normalise the eigenfunction at infinity:

(3.28) \begin{equation} \varPsi _n \sim \exp ^{-k_nx} \quad \text{as} \quad x \rightarrow + \infty . \end{equation}

We can therefore analytically describe the main soliton and the extra solitons produced for every direction defined by the polar angle $\theta$ .

When there is no perturbation, i.e. $\varDelta = 1$ , the solution is simply the exact line-soliton solution. When $\varDelta \gt 1$ , there is a localised perturbation to the wavefront, detectable in some region given by $-\tilde {\theta }_1 \lesssim \tilde {\theta } \lesssim \tilde {\theta }_1$ such that the amplitude $A(T,\xi ,\theta ) \gtrsim A(T_0,\xi ,\theta ) + \epsilon$ . The estimate for the region determining where the perturbation lies is therefore given by the formula $4 k_1(\tilde \theta )^2/3 \gtrsim - 2 U(\tau _0,X;\tilde \theta )/3 + \epsilon$ arising from (3.24). The small secondary solitons are then detectable when $A(T,\xi ,\theta ) \gtrsim \epsilon$ and, hence, when $4 k_n(\tilde \theta )^2/3 \gtrsim \epsilon$ . These regions are shown in figure 4, where we compare the performance of the KdV $\theta$ equation with the results of direct numerical simulations for the parent system.

Figure 4. Two-dimensional (a,b) and contour (c,d,e, f) plots of the numerical solution of both the 2-D Boussinesq–Peregrine system (a,b,d,f) and the KdV $\theta$ equation (c,e) for a perturbed line soliton, where $\alpha = 2$ , $\beta = 5$ , $\tilde v = 0.5$ , $\xi _0 = 180$ , $\epsilon = 0.05$ and $T\in [1,6]$ . The dark grey (dash-dotted) and the light grey (dashed) lines illustrate theoretical (IST) predictions for the perturbed regions of the main and secondary solitons, respectively.

The corresponding initial condition in the $xy$ plane for the 2-D Boussinesq–Peregrine system (2.1) and (2.2) is obtained using the initial condition to the KdV $\theta$ equation (3.13):

(3.29) \begin{align} & \eta (t_0,x,y) = 2 \tilde v \left (1 + \alpha \operatorname {sech}^2 \left ( \beta \arctan \frac{y}{x} \right ) \right ) \operatorname {sech}^2 \left (\frac{1}{2} \sqrt {6 \tilde v} \left (x - (1 + \epsilon \tilde v)t_0 - \xi _0 \right ) \right )\!, \end{align}

(3.30) \begin{align} & u(t_0,x,y) = \eta (t_0,x,y), \quad v(t_0,x,y) = 0. \end{align}

Here $\xi _0 = \textrm {const}$ . The computation parameters can be found in table 1 and the corresponding numerical results are shown in figure 4.

We solve the 2-D Boussinesq–Peregrine system for $t \in [20,120]$ and $\epsilon = 0.05$ , where comparison is made to the KdV $\theta$ equation that is solved for $T \in [1,6]$ . Comparisons can therefore be made at the same time. There is good quantitative agreement between the two solutions around the slow time value $T = 2$ . The second soliton is emerging and there is also some extra radiation in the 2-D Boussinesq–Peregrine system. In longer runs, shown for the slow time value $T = 6$ , the agreement between the KdV $\theta$ equation and the 2-D Boussinesq–Peregrine equations worsens. The main perturbed soliton starts to split sideways, spreading out across the wavefront and in doing so reduces in amplitude in the centre. This feature is not seen in the KdV $\theta$ equation. The second soliton has now fissioned from the original soliton, and the analytical estimates based on the IST give good predictions for the perturbed regions of the main and secondary solitons. In the 2-D Boussinesq–Peregrine equations the second soliton has radiated more than in the KdV $\theta$ regime, leading to the noticeable difference. However, there is still good qualitative agreement between the full parent system and the reduced model, despite the strong perturbation and long propagation time.

Eventually, all strong, long transverse amplitude perturbations split and travel sideways along the wavefront with two ring-wave arcs being radiated as the perturbation travels bidirectionally along the wavefront. For long enough perturbations, the $L^{\infty }$ norm $||d||_\infty$ equal to the maximum absolute value of the difference between solutions of the 2-D Boussinesq–Peregrine system and the KdV $\theta$ equation is consistent with the model and scales as $\textit{O}(\epsilon )$ . The computation parameters in table 2 correspond to figure 5. In figure 5, for the strong perturbation of $\alpha = 0.5$ considered here and the perturbation width parameter $\beta = 20$ , the fitted line has the gradient $0.76$ . When the perturbation is wider, $\beta = 10$ , the fitted line has a greater gradient of 1.1. Therefore, for sufficiently long perturbations, the KdV $\theta$ equation is accurate to $\textit{O}(\epsilon )$ . We note that the KPII equation is valid too in this regime, as well as for narrower perturbations. However, the KdV $\theta$ equation is a much simpler model, which can be used to describe the intermediate asymptotics of plane waves subject to sufficiently long transverse perturbations of finite amplitude.

Table 1. Computation parameters for the KdV $\theta$ equation and 2-D Boussinesq–Peregrine system to generate figure 4.

Table 2. Computation parameters for the KdV $\theta$ equation and 2-D Boussinesq–Peregrine system to generate figure 5.

Figure 5. Convergence of the solution of the KdV $\theta$ equation to the solution of the 2-D Boussinesq–Peregrine system for different values of parameter $\beta$ of the transverse perturbation and small parameter $\epsilon$ , where $||d||_\infty$ is the $L^\infty$ norm of the difference between solutions of the 2-D Boussinesq–Peregrine system and KdV $\theta$ equation, where $\alpha = 0.5$ , $\tilde v = 0.5$ , $\xi _0 = 200$ and $T \in [1,2]$ .

4. Perturbed ring waves

In this section we extend the studies of concentric ring waves by Chwang & Wu (Reference Chwang and Wu1977), Sidorovas et al. (Reference Sidorovas, Tseluiko, Choi and Khusnutdinova2024) to non-axisymmetric perturbed ring waves, within the framework of the 2-D Boussinesq–Peregrine system (2.1) and (2.2). We use the cKdV equation (2.36) and the related formulae (2.37) for velocity components to define the initial conditions for our numerical runs, and consider both localised and periodic perturbations of concentric waves.

We consider a ring wave that is placed sufficiently far away from the origin. Then, the cylindrical divergence/convergence term in the cKdV equation (2.36) can be treated as a perturbation of the KdV equation, and one can consider the axisymmetric solution initiated by the KdV soliton bent into a ring wave (see Ko & Kuel Reference Ko and Kuel1979; Dorfman, Pelinovsky & Stepanyants Reference Dorfman, Pelinovsky and Stepanyants1981; Johnson Reference Johnson1999; Sidorovas et al. Reference Sidorovas, Tseluiko, Choi and Khusnutdinova2024):

(4.1) \begin{equation} \eta (t_0,x,y) = 2 \tilde v \operatorname {sech}^2 \left ( \frac{1}{2} \sqrt {6 \tilde v} \left [\sqrt {x^2 + y^2} - (1 + \epsilon \tilde v)t_0 - r_0 \right ] \right )\!. \end{equation}

Here, $\tilde v \gt 0$ and $r_0$ are some constants. The Cartesian velocities $u$ and $v$ are defined using the related approximations given by (2.37).

Next, we consider the initial condition with a localised perturbation for the 2-D Boussinesq–Peregrine system (2.1) and (2.2) in the form

(4.2) \begin{align} \eta (t_0,x,y) & = 2 \tilde v \big(1 + \alpha H(x) \operatorname {sech}^2 \beta y \big) \operatorname {sech}^2\left ( \frac{1}{2} \sqrt {6 \tilde v} \left [\sqrt {x^2 + y^2} - (1 + \epsilon \tilde v)t_0 - r_0 \right ] \right )\!, \end{align}

(4.3) \begin{align} u(t_0,x,y) & = \frac{x}{\sqrt {x^2 + y^2}} \, \eta (t_0,x,y), \quad v(t_0,x,y) = \frac{y}{\sqrt {x^2 + y^2}} \, \eta (t_0.x,y), \end{align}

where $\alpha$ and $\beta$ control the amplitude and width of the perturbation, respectively, and $H = H(x)$ is the Heaviside function to only give the perturbation on one side of the ring. At the origin of the $xy$ plane, we set $\eta (t_0,0,0) = u(t_0,0,0) = v(t_0,0,0) = 0$ to avoid division by zero, and, for all ring-wave simulations, we use the computation parameters in table 3.

Table 3. Computation parameters for the 2-D Boussinesq–Peregrine system simulations of ring waves.

For the outward-propagating ring wave, we solve the 2-D Boussinesq–Peregrine system for $t \in [0,150]$ , $\epsilon = 0.01$ , $r_0 = 20$ , $\tilde v = 0.5$ and $\alpha = 0.5$ for the two values of $\beta$ shown, and the computation parameters in table 3 give figure 6.

Figure 6. Two-dimensional plots of an outward-propagating perturbed ring wave plotted at the times $t = 0,50,100,150$ for $\beta = 2$ (a) and $\beta = 20$ (b), where $\alpha = 0.5$ , $\tilde v = 0.5$ , $r_0 = 20$ , $\epsilon = 0.01$ and $t \in [0,150]$ .

We consider two perturbations. The wider perturbation, with $\beta = 2$ , has a width comparable to the radius of the initial ring. It expands and decreases in amplitude as the wave radially diverges. The narrower perturbation, with $\beta = 20$ , has a width comparable to the wavelength of the concentric wave. It instead rapidly splits and travels sideways along the wavefront, eventually becoming unnoticeable. Both outward-propagating perturbed ring waves are stable.

To generate inward-propagating ring waves in the 2-D Boussinesq–Peregrine system, we multiply the velocities (4.3) by $-1$ , to reverse the direction of propagation. We then solve this system for $t \in [0,250]$ , $\epsilon = 0.01$ , $r_0 = 180$ , $\tilde v = 0.5$ and $\alpha = 0.5$ for the two values of $\beta$ shown, and the computation parameters listed in table 3. The results are shown in figures 7 and 8.

Figure 7. Two-dimensional plots of an inward-propagating perturbed ring wave plotted at the times $t=0,50,100,150$ (a) and $t = 200$ (b), where $\alpha = 0.5$ , $\beta = 2$ , $\tilde v = 0.5$ , $r_0 = 180$ , $\epsilon = 0.01$ and $t \in [0,250]$ .

Figure 8. Two-dimensional plots of an inward-propagating perturbed ring wave plotted at the times $t=0,50,100,150$ (a) and $t = 200$ (b), where $\alpha = 0.5$ , $\beta = 20$ , $\tilde v = 0.5$ , $r_0 = 180$ , $\epsilon = 0.01$ and $t \in [0,250]$ .

The wider perturbation, with $\beta = 2$ , comparable to the radius of the ring, focuses to a point without splitting. The narrower perturbation, with $\beta = 20$ , comparable to the wavelength, splits, taking longer to do so than its outward-propagating counterpart. In both figures 7 and 8 the wave propagates through the origin and reforms the general shape as before, with the perturbation seemingly remaining on the same path, now as an outward-propagating perturbed ring wave. The case of the narrower ( $\beta = 20$ ) perturbation has more radiation on the trailing edge of the wave compared with the wider ( $\beta = 2$ ) perturbation case. When the perturbation is narrower, and once it has passed through the origin, the perturbation starts to deform the shape of the ring, as seen in figure 9. However, when the wave is at the origin, with the maximum amplitude $||\eta ||_{\infty } = 44.12$ occurring at $t = 177.5$ , the wave is almost concentric. Once the wave has passed through the origin, it retains the form of the perturbation, now as an outward-propagating perturbed ring wave. The wavefront is constantly radiating small waves, and the moustache-shaped radiation from the perturbation is clearly visible in both figures 8 and 9. Once the wave has propagated sufficiently far outward, it retains its concentric shape.

Figure 9. Contour plots of the inward-propagating perturbed ring wave from figure 8 before (a), during (b) and after (c,d) the propagation through the origin, where $\alpha = 0.5$ , $\beta = 20$ , $\tilde v = 0.5$ , $r_0 = 180$ , $\epsilon = 0.01$ .

We now apply a periodic perturbation to the inward-propagating ring wave. Theoretical considerations by Ostrovskii & Shrira (Reference Ostrovskii and Shrira1976), Shrira & Pesenson (Reference Shrira and Pesenson1984), Pesenson (Reference Pesenson1991) suggest that there exists a critical wavelength of the perturbation, such that above this wavelength, the perturbed ring wave is unstable, whereas below this critical value, it is stable. However, an analytical formula for this critical wavelength is not available at present, and deriving it is a separate problem beyond the scope of the present study. Here, in the first numerical study of the problem, we aim to provide an illustrative example of this qualitative difference. Taking the initial condition

(4.4) \begin{align} \eta (t_0,x,y) &= 2 \tilde v \!\left (1.25 + \alpha \cos \left ( \beta \arctan \frac{y}{x} \right )\right ) \operatorname {sech}^2\left ( \frac{1}{2} \sqrt {6 \tilde v} \left [\sqrt {x^2 + y^2} - x_0 \right ] \right )\!, \end{align}

(4.5) \begin{align} u(t_0,x,y) & = -\frac{x}{\sqrt {x^2 + y^2}} \, \eta (t_0,x,y), \quad v(t_0,x,y) = -\frac{y}{\sqrt {x^2 + y^2}} \, \eta (t_0,x,y), \end{align}

we solve the 2-D Boussinesq–Peregrine system (2.1) and (2.2) for $t \in [0,150]$ , $\epsilon = 0.01$ , $\tilde v = 0.5$ and $\alpha = 0.5$ for $\beta = 20$ and $\beta =4$ , and the computation parameters listed in table 3. The results are shown in figure 10.

Figure 10. Two-dimensional (a,b) and contour (c,d) plots of an inward-propagating periodically perturbed ring wave for $\beta = 20$ (a,c) and $\beta = 4$ (b,d), where $\alpha = 0.5$ , $\tilde v = 0.5$ , $r_0 = 180$ , $\epsilon = 0.01$ and $t \in [0,150]$ . The 2-D plots (a,b) are plotted at the times $t = 0, 50, 100, 150$ and the corresponding contour plots (c) and (d) are plotted at $t = 150$ .

It is clear that the higher-frequency periodic perturbation does not grow faster than the rate of the radial convergence, whereas the low-frequency counterpart does. Both cases produce periodic radiation patterns around the ring. However, the higher-frequency perturbation radiates considerably more, as seen in figure 10, but remains largely concentric, in contrast to the low-frequency perturbation, which notably deforms the ring. Hence, our simulations show significant qualitative differences in the behaviour of these two perturbations.

5. Pure and perturbed hybrid waves

In this section we construct and model hybrid waves, generated by initial conditions formed by an arc of a ring wave and two tangent plane waves (Ostrovskii & Shrira Reference Ostrovskii and Shrira1976; Hooper et al. Reference Hooper, Khusnutdinova and Grimshaw2021). The two regions are constructed using the same description for plane and ring waves as above, and the wavefront has a continuous derivative. Surface signatures of similarly looking internal waves generated in the Strait of Gibraltar can be seen in the SAR images shown in figure 1, as well as figures 13 and 14 of Apel (Reference Apel2003). As discussed in the introduction, all considerations of our present work can be extended to internal waves and to parallel shear flows.

We use the exact solution of the KdV $\theta$ equation given by

(5.1) \begin{eqnarray} &&\eta (T,\xi ,\theta ) = 2\tilde v \operatorname {sech}^2 \left (\frac{1}{2} \sqrt {6\tilde v} \left [ \xi - \tilde v T - \xi _0 \right ] \right )\!, \end{eqnarray}

where $\xi = \textit{rk}(\theta ) - \textit{st}$ , $T = \epsilon t$ , to define the plane-wave sections of a hybrid wave. We match the two plane waves to an arc of a ring wave centred at the origin at the polar angles $\theta = \pm \varphi$ .

The plane-wave sections are defined for $|y| \geqslant x \tan \varphi$ and we use the general solution of (1.12), $k(\theta ) = a \cos \theta + b(a) \sin \theta$ , and hence $\xi = ax + b(a)y - t$ , with $a = \cos \varphi$ and $b(a) = \sin \varphi$ . In the $xy$ plane, using the Cartesian velocities (2.24), the initial condition becomes

(5.2) \begin{align} \eta (t_0,x,y) & = 2\tilde v \operatorname {sech}^2\left (\frac{1}{2} \sqrt {6\tilde v} \left [ ax \pm b(a)y - (1 + \epsilon \tilde v)t_0 - \xi _0 \right ] \right )\!, \end{align}

(5.3) \begin{align} u(t_0,x,y) & = \cos \varphi \, \eta (t_0,x,y), \quad v(t_0,x,y) = \pm \sin \varphi \, \eta (t_0,x,y_0), \end{align}

where the upper and lower signs in $\eta$ and the velocity $v$ are for the upper and lower plane-wave sections, respectively.

For $|y| \lt x \tan \varphi$ , we take $k(\theta )$ to be the singular solution of (1.12), i.e. the geometrical envelope of the general solution, giving $k(\theta ) = 1$ , and hence $\xi = r - t$ . Therefore, in the $xy$ plane, the ring-wave part of the initial condition is given by

(5.4) \begin{align} \eta (t_0,x,y) & = 2\tilde v \operatorname {sech}^2\left (\frac{1}{2} \sqrt {6\tilde v} \left [ \sqrt {x^2 + y^2} - (1 + \epsilon \tilde v)t_0 - r_0 \right ] \right )\!, \end{align}

(5.5) \begin{align} u(t_0,x,y) & = \frac{x}{\sqrt {x^2+y^2}} \, \eta (t_0,x,y), \quad v(t_0,x,y) = \frac{y}{\sqrt {x^2+y^2}} \, \eta (t_0,x,y). \end{align}

The initial condition is therefore defined in the entire $xy$ plane, and we match the sections by choosing the constants $\xi _0 = r_0$ , so that at the initial moment of time there is a smooth transition between plane and ring parts at $\theta = \pm \varphi$ .

We modify the initial conditions to bring them to zero near the boundaries by multiplying them by the double $\tanh$ -function $M_x(x) M_y(y)$ , where

(5.6) \begin{align} M_x(x) & = \frac{1}{2} \left [\tanh \kappa (x - x_{min} - x_{\textit{span}}) - \tanh \kappa (x - x_{max} + x_{\textit{span}}) \right ], \end{align}

(5.7) \begin{align} M_y(y) & = \frac{1}{2} \left [\tanh \kappa (y - y_{min} - y_{\textit{span}}) - \tanh \kappa (y - y_{max} + y_{\textit{span}}) \right ], \end{align}

to maintain the periodic boundary conditions during our pseudospectral numerical runs. Typically, we take the parameter values $\kappa = 0.25$ , $x_{\textit{span}} = (x_{\textit{max}} - x_{\textit{min}})/20$ and $y_{\textit{span}} = (y_{\textit{max}} - y_{\textit{min}})/20$ . Such modification to the amplitude generates cylindrical radiation at the ends of the legs, similar to that described by Yuan & Wang (Reference Yuan and Wang2022), and we ensure that the computational domain is sufficiently large, so that these artefacts do not affect the central region. We solve the 2-D Boussinesq–Peregrine system (2.1) and (2.2) for $t \in [0,150]$ , $\epsilon = 0.01$ , $\tilde v = 0.5$ , $r_0 = 150$ and computation parameters listed in table 4. The results are shown in figure 11.

Table 4. Computation parameters for the 2-D Boussinesq–Peregrine system and KPII equation simulation of hybrid waves.

Figure 11. Two-dimensional (a,b,e,f) and contour (c,d) plots of an outward-propagating hybrid wave in the 2-D Boussinesq–Peregrine system, where $\tilde v = 0.5$ , $r_0 = 150$ , $\varphi = \pm \pi /8$ , $\epsilon = 0.01$ and $t \in [0,150]$ .

The initial condition evolves into a wavefront with three distinct sections and a small amount of radiation, not visible here in figure 11. The plane-wave legs propagate with constant amplitude, retaining their shape, and are well described by the exact soliton solution (5.2) and (5.3), obtained using the KdV $\theta$ equation (3.7). The amplitude of the ring-wave section gradually decreases due to the radial divergence and, provided the wave is sufficiently far from the origin, it is well described by the cKdV equation (2.36) (Sidorovas et al. Reference Sidorovas, Tseluiko, Choi and Khusnutdinova2024). The third section is the transition region between the legs and the central ring wave. It is not described by either the KdV $\theta$ or cKdV equations, but is described by the KPII equation (2.47) and, potentially, the modulation theory of Biondini, Hoefer & Moro (Reference Biondini, Hoefer and Moro2020), Ryskamp et al. (Reference Ryskamp, Maiden, Biondini and Hoefer2021) could also be extended to describe this region.

This outward-propagating hybrid wave can be described using the KPII equation (2.47). To do so, one must rewrite the initial conditions and relevant regions using the KPII variables $T = \epsilon t$ , $\xi = x - t$ and $Y = \sqrt {\epsilon } y$ . The plane-wave section ( $|y| \geqslant x\tan \varphi$ ) is mapped to $|Y| \geqslant \sqrt {\epsilon } (\xi + T_0/\epsilon ) \tan \varphi$ , where the surface elevation (5.2) is given by

(5.8) \begin{equation} \tilde \eta (T_0,\xi ,Y) = 2\tilde v \operatorname {sech}^2 \left (\frac{1}{2} \sqrt {6\tilde v} \left [ a \xi + \frac{b(a)}{\sqrt {\epsilon }}Y - \left (\tilde v + \frac{1 - a}{\epsilon } \right ) T_0 - r_0 \right ] \right )\!. \end{equation}

The ring-wave section ( $|y| \lt x\tan \varphi$ ) is mapped to $|Y| \lt \sqrt {\epsilon } (\xi + T_0/\epsilon ) \tan \varphi$ , where the surface elevation (5.4) is given by

(5.9) \begin{equation} \tilde \eta (T_0,\xi ,Y) = 2\tilde v \operatorname {sech}^2 \left (\frac{1}{2} \sqrt {6\tilde v} \left [ \sqrt { \left (\xi + \frac{T_0}{\epsilon } \right )^2 + \frac{Y^2}{\epsilon }} - \left (\tilde v + \frac{1}{\epsilon } \right ) T_0 - r_0 \right ] \right )\!. \end{equation}

Initial conditions for the KPII equation must have zero mass (Klein, Sparber & Markowich Reference Klein, Sparber and Markowich2007; Biondini, Bivolcic & Hoefer Reference Biondini, Bivolcic and Hoefer2025). To abide by this constraint, we modify the initial condition $\tilde \eta$ by adding, for each $Y_i$ in the $Y$ domain, a pedestal of the form

(5.10) \begin{equation} p_i = \frac{p_0}{2} \left (\tanh \kappa (\xi - \xi _{\textit{min}} - \xi _{\textit{span}}) - \tanh \kappa (\xi - \xi _{\textit{max}} - \xi _{\textit{span}}) \right )\!, \end{equation}

where $p_0$ and $\kappa$ are some suitable constants. One can then define the massless initial condition to be

(5.11) \begin{eqnarray} &&\eta (T_0,\xi ,Y_i) = \tilde \eta (T_0,\xi ,Y_i) - p_i(\xi ), \end{eqnarray}

where

(5.12) \begin{eqnarray} && \int _{\xi _{\textit{min}}}^{\xi _{\textit{max}}} \tilde \eta (T_0,\xi ,Y_i) \, \textrm{d}\xi = \int _{\xi _{\textit{min}}}^{\xi _{\textit{max}}} p_i(\xi ) \, \textrm{d}\xi . \end{eqnarray}

We usually take $\kappa = 0.5$ , $\xi _{\textit{span}} = (\xi _{\textit{max}} - \xi _{\textit{min}})/20$ . In practice, due to the size of the domain, this leads to a pedestal of amplitude less than $\epsilon$ , and therefore, the effect on the numerical solution is negligible. Further discussion of the numerical solution to the KPII equation is given in Appendix B.

As the outward-propagating hybrid wave moves away from the origin, the effects of cylindrical divergence are weak and the wavefront propagates in a quasi-stationary manner. The description of the central part of the wavefront can therefore be given analytically using the asymptotic solution constructed by Johnson (Reference Johnson1999). It is given for the initial-value problem

(5.13) \begin{eqnarray} && 2 \eta _T + 3\eta \eta _{\xi } + \frac{1}{3} \eta _{\xi \xi \xi } + \frac{\eta }{T} = 0, \quad \eta \!\left (T_0,\xi \right ) = A \operatorname {sech}^2 \left (\frac{\sqrt {3A}}{2} \xi \right )\!, \end{eqnarray}

which is posed sufficiently far away from the origin, and then the radial-divergence term can be viewed as a perturbation of the KdV equation. The three components of the wave are the primary wave, the shelf and the transition region leading back to the undisturbed medium. The variables

(5.14) \begin{align} \sigma & = T_0^{-1}, \quad X = \sigma T, \quad \mathcal{T} = \frac{1}{2}\sqrt {3{\textit{AX}}^{-\frac{2}{3}}}, \quad f(X) = \frac{3}{2}A\!\left (X^{\frac{1}{3}} - 1 \right )\!, \end{align}

(5.15) \begin{align} \varTheta & = \xi - \sigma ^{-1} f(X), \quad \bar {s} = \operatorname {sech}(\mathcal{T}\varTheta ), \quad \bar {t} = \tanh (\mathcal{T} \varTheta ) \end{align}

are required to construct the three components. The primary wave is given by

(5.16) \begin{equation} \eta _{\textit{primary}} = {\textit{AX}}^{-\frac{2}{3}}\bar {s}^2 + \frac{2\sigma }{3}\frac{X^{-\frac{2}{3}}}{\sqrt {3A}} \left [\!-1 + \bar {t} + (3 + 2\mathcal{T}\varTheta )\bar {s}^2 - \left (\frac{35}{12} + 3\mathcal{T}\varTheta + \mathcal{T}^2\varTheta ^2 \right )\bar {t}\bar {s}^2\! \right ], \end{equation}

the shelf is described by

(5.17) \begin{equation} \eta _{\textit{shelf}} = {\textit{AX}}^{-\frac{2}{3}}\bar {s}^2 - \frac{4\sigma }{3} \left (3A + 2 \sigma \xi \right )^{-\frac{1}{2}} \end{equation}

and, finally, the transition region is given in the form

(5.18) \begin{equation} \eta _{\textit{transition}} = {\textit{AX}}^{-\frac{2}{3}}\bar {s}^2 - \frac{4\sigma }{3} \left (3A \right )^{-\frac{1}{2}} \left (1 - \int _{\xi (2\sigma /\! X)^{1/3}}^{\infty } {\textit{Ai}}(\xi ') \, \textrm{d}\xi ' \right )\!, \end{equation}

where $ {\textit{Ai}}$ is the Airy function. The three components can be used to describe the central part of the hybrid wave in $(T,\xi )$ space for suitably small $\sigma$ and, hence, suitably large $T_0$ . This asymptotic solution was tested for a pure ring wave by Sidorovas et al. (Reference Sidorovas, Tseluiko, Choi and Khusnutdinova2024), and therefore, is applicable in the central part of the hybrid wave, which is well described by the cKdV equation. Hence, one can describe analytically the three main parts of the outward-propagating hybrid wave: the two legs and the central region.

We now model the inward-propagating hybrid wave. The wave is generated by the same set of initial conditions, (5.2)–(5.5), in the same regions. However, the velocities are multiplied by $-1$ , to reverse the direction of propagation, and similarly brought to zero near the boundaries. We solve the 2-D Boussinesq–Peregrine system (2.1) and (2.2) for $t \in [0,150]$ , $\epsilon = 0.01$ , $\tilde v = 0.5$ , $r_0 = 100$ and the computation parameters listed in table 4. The results are shown in figures 12–14.

Figure 12. Two-dimensional plots of an inward-propagating hybrid wave in the 2-D Boussinesq–Peregrine equations, where $\tilde v = 0.5$ , $r_0 = 100$ , $\varphi = \pm \pi /4$ , $\epsilon = 0.01$ and $t \in [0,150]$ .

Figure 13. The $L^\infty$ norm of $\eta$ for the inward-propagating hybrid wave, where $\tilde v = 0.5$ , $r_0 = 100$ , $\varphi = \pi /4$ , $\epsilon = 0.01$ and $t\in [0,150]$ . The results are shown for the varying computational grid sizes $\varDelta _x = \varDelta _y = 0.35$ (black), $\varDelta _x = \varDelta _y = 0.275$ (blue), $\varDelta _x = \varDelta _y = 0.2$ (green) and $\varDelta _x = \varDelta _y = 0.125$ (magenta).

Figure 14. Cross-section plot along $x = 0$ of the domain from figure 12 plotted five times from $t = 0$ to $t = 200$ .

The intermediate evolution of the wavefront is dominated by the ring-wave section focusing in on itself, increasing the amplitude dramatically, as shown by the $L^{\infty }$ norm in figure 13. However, this is not a finite-time blow-up of the numerical solution, since, as the resolution of the mesh is increased, the amplitude converges to a consistent value for the same $t$ value. The converging wavefront generates a large rogue wave (lump), which subsequently disintegrates into three pieces and eventually generates a stable ‘X-type’ formation, similar to that observed within the scope of the KPII equation in figure 1 of Ablowitz & Baldwin (Reference Ablowitz and Baldwin2012) and studied in Chakravarty & Kodama (Reference Chakravarty and Kodama2009, Reference Chakravarty and Kodama2014); Chakravarty, Lewkow & Maruno (Reference Chakravarty, Lewkow and Maruno2010). There also exists an arc of opposite curvature to the initial condition connecting two lump-like structures on the legs behind the X. Previous studies have focused on ‘V-type’ initial data for the KPII equation rather than our hybrid wave formation. These studies have shown similar X-type structures obtained numerically without a noticeable phase shift, see, for example, Chakravarty & Kodama (Reference Chakravarty and Kodama2009), Ryskamp et al. (Reference Ryskamp, Maiden, Biondini and Hoefer2021), Ryskamp, Hoefer & Biondini (Reference Ryskamp, Hoefer and Biondini2022). A weaker arc of opposite curvature to the initial condition connecting the two legs was observed by Yuan et al. (Reference Yuan, Grimshaw, Johnson and Wang2018a ), Ryskamp et al. (Reference Ryskamp, Maiden, Biondini and Hoefer2021) in the KPII equation and from the cross-section in figure 14 has a shape similar to the plane-wave derivative. The effect of curvature on the amplitude of the generated waves was noted in Porubov et al. (Reference Porubov, Truji, Lavrenov and Oikawa2005). To the best of our knowledge, the overall behaviour of the inward-propagating hybrid wave, including the generation of the large rogue wave and its disintegration into the three pieces, with the subsequent formation of the Xwave, has not been reported in any previous studies. Although the KPII equation does not support the exact lump solutions, unlike the KPI equation, our study shows that lumps, understood as localised 2-D waves, can exist as transient (emerging and then disappearing) features in the evolution of curved solitons, contributing to the mechanisms of generation of rogue waves (see, for example, Shrira & Pesenson Reference Shrira and Pesenson1984; Kharif, Pelinovsky & Slunyaev Reference Kharif, Pelinovsky and Slunyaev2008; Onorato et al. Reference Onorato, Residori, Bortolozzo, Montina and Arecchi2013; He & Chabchoub Reference He and Chabchoub2025; Nirunwiroj, Tseluiko & Khusnutdinova Reference Nirunwiroj, Tseluiko and Khusnutdinova2025 and references therein).

It turns out that the X-type wave formed at the end can be described by the solution of the KPII equation discussed, in a different context, by Ablowitz & Baldwin (Reference Ablowitz and Baldwin2012):

(5.19) \begin{equation} \eta = \frac{4}{3} \frac{\partial ^2}{\partial \xi ^2} \log \left ( F_2 \right )\!. \end{equation}

Here

(5.20) \begin{align}& F_2 = 1 + e^{\eta _1} + e^{\eta _2} + e^{\eta _1 + \eta _2 + A_{12}}, \end{align}

(5.21) \begin{align} & \eta _i = 6k_i \!\left ( \xi + 6P_iY - 6\!\left ( k_i^2 + 3 P_i^2 \right )\! T \!\right ) + c_i, \quad i = 1,2, \end{align}

(5.22) \begin{align} & e^{A_{12}} = \frac{(k_1 - k_2)^2 - (P_1 - P_2)^2}{(k_1 + k_2)^2 - (P_1 - P_2)^2}, \end{align}

such that $k_i$ and $P_i$ are constants determining the propagation angle and $c_i$ is an arbitrary constant shifting the wavefront. For our numerical run, the crossing angle is $\pi /4$ , and good agreement is obtained for the parameter values $6k_1 = 6k_2 = 1.1314$ , $6P_1 = -6P_2 = \epsilon ^{-1/2}$ and $c_1 = c_2 = 71.75k_1$ . The value of $e^{A_{12}}$ is indicative of the behaviour of the wavefront where the two solitons cross. For the case seen here, $e^{A_{12}} = 1.013$ , implying an X-type interaction with a very short stem between the legs, as seen from the evolved state of the inward-propagating hybrid wavefront. Indeed, there is no observable phase shift and the comparison between the KPII two-soliton solution and the inward-propagating hybrid wave is given in figure 15. The amplitude of the two-soliton cross depends on the propagation angle, and a discussion on how to obtain the maximum of the cross is discussed by Chakravarty et al. (Reference Chakravarty, Lewkow and Maruno2010).

Figure 15. Contour plots of (a) the inward-propagating hybrid wave at $t = 150$ from figure 12 and (b) the parameter-matched ‘X-type’ two-soliton solution of the KPII equation.

To model the entire process of the inward-propagating hybrid wave seen in figure 12 with the KPII equation, we must first reverse the direction of propagation. Therefore, for $T$ reducing, we map $t \rightarrow -t$ and, hence, $\partial\! /\!\partial T = - \partial\! /\!\partial T$ and $\xi \rightarrow \zeta = x + t$ , which yields the inward-propagating KPII equation

(5.23) \begin{equation} \left ( 2\eta _T - 3\eta \eta _{\zeta } - \frac{1}{3} \eta _{\zeta \zeta \zeta } \right )_{\zeta } - \eta _{\textit{YY}} = 0. \end{equation}

The initial condition for the plane-wave section ( $|y| \geqslant x \tan \varphi$ ) is defined when $|Y| \geqslant \sqrt {\epsilon } (\zeta - T_0/\epsilon ) \tan \varphi$ and is given by

(5.24) \begin{equation} \tilde \eta (T_0,\zeta ,Y) = 2\tilde v \operatorname {sech}^2 \left (\frac{1}{2} \sqrt {6\tilde v} \left [ a \zeta + \frac{b(a)}{\sqrt {\epsilon }}Y - \left (\tilde v + \frac{1 + a}{\epsilon } \right ) T_0 - \xi _0 \right ] \right )\!, \end{equation}

and the ring-wave section ( $|y| \lt x \tan \varphi$ ), defined when $|Y| \lt \sqrt {\epsilon } (\zeta - T_0/\epsilon ) \tan \varphi$ , is given by

(5.25) \begin{equation} \tilde \eta (T_0,\zeta ,Y) = 2\tilde v \operatorname {sech}^2 \left (\frac{1}{2} \sqrt {6\tilde v} \left [ \sqrt { \left (\zeta - \frac{T_0}{\epsilon } \right )^2 + \frac{Y^2}{\epsilon }} - \left ( \tilde v + \frac{1}{\epsilon } \right ) T_0 - \xi _0 \right ] \right )\!, \end{equation}

where, before solving, the mass of the solution should be removed to satisfy the zero-mass constraint of the KPII equation as before, such that $\eta = \tilde \eta - p(\xi )$ for the pedestal $p(\xi )$ . Solving (2.47) and (5.23) with the initial data of (5.8), (5.9) and (5.24), (5.25), respectively, for the parameter values given to match figures 11 and 12, we obtain the results that are shown in figure 16.

Figure 16. Two-dimensional plots of the numerical solution to the KPII equation for both the outward- (a) and inward- (b) propagating hybrid waves for $\varphi = \pi /8$ , $r_0 = 150$ (a) and $\varphi = \pi /4$ , $r_0 = 100$ (b) at $|T|= 1.5$ , where $\tilde v = 0.5$ , $\epsilon = 0.01$ and $|T| \in [0,1.5]$ .

We see that the computations give all of the same structures as those observed in the 2-D Boussinesq–Peregrine system, up to the accuracy of the reduced model. Notably, the transient lump-like structures in the middle and at the ends of the legs are present, as in figure 12, despite the fact that the KPII equation for weak surface tension used here does not possess lump soliton solutions. These transient lumps constantly decrease in amplitude, converting their energy to extend the legs behind the X. We also note that some radiation escapes and re-enters the periodic domain of the numerical solver, causing some small defects. The inclusion of a sponge layer can remove these artefacts.

Finally, in this section we model the evolution of an inward-propagating perturbed hybrid wave. We choose to apply a periodic perturbation to the ring-wave section, altering (5.4) and (5.5) but keeping (5.2) and (5.3) unchanged. We again take the velocities to be negative for inward propagation. We perturb the arc of a ring wave in the region $|y| \lt x \tan \varphi$ to be given by

(5.26) \begin{align} \eta (t_0,x,y) & = \left ( \frac 32 + \alpha \cos \left ( \beta \arctan \frac yx \right ) \right ) \eta _0(t_0,x,y), \end{align}

(5.27) \begin{align} u(t_0,x,y) & = - \frac{x}{\sqrt {x^2+y^2}} \, \eta (t_0,x,y), \quad v(t_0,x,y) = - \frac{y}{\sqrt {x^2+y^2}} \, \eta (t_0,x,y). \end{align}

Solving the 2-D Boussinesq–Peregrine system (2.1) and (2.2) for the perturbed initial condition given by (5.2) and (5.3), and (5.26) and (5.27) for $t \in [0,150]$ , $\epsilon = 0.01$ , $\tilde v = 0.5$ , $r_0 = 100$ , $\alpha = 0.5$ , $\beta = 12$ and the computation parameters listed in table 4, yields the results shown in figure 17.

Figure 17. Two-dimensional plots of a perturbed inward-propagating hybrid wave in the 2-D Boussinesq–Peregrine system, where $\alpha = 0.5$ , $\beta = 12$ , $\tilde v = 0.5$ , $r_0 = 100$ , $\varphi = \pm \pi /4$ , $\epsilon = 0.01$ and $t \in [0,200]$ .

The evolution of the inward-propagating perturbed hybrid wave is remarkably similar to that of the unperturbed case, despite the large amplitude and wavelength of the perturbation applied. The perturbed hybrid wave also generates a large rogue wave (larger than in the unperturbed case) and still eventually forms a stable X-type cross. The initial radiation behind the ring-wave section is qualitatively similar to that observed in the ring-wave study in § 4. The main difference is that the transient rogue wave disintegrates into more pieces than in the pure hybrid wave evolution.

Finally, we note that the $(2+1)$ -dimensional cKdV-type equation (1.7) can be used to initiate hybrid waves in the general case of a stratified fluid with a shear flow (see Hooper et al. Reference Hooper, Khusnutdinova and Grimshaw2021).