Hydroelastic lumps in shallow water

We consider a three-dimensional incompressible and inviscid fluid of density $\rho_{w}$ covered by a deformable ice sheet of density $\rho_{i}$ and the thickness $d$. The displacement of the ice cover is denoted by $z=\eta(x,y,t)$, where $x$ and $y$ are horizontal coordinates, and the $z$-axis points upwards with $z=0$ the undisturbed ice sheet. The fluid is bounded below by a rigid bed located at $z=-h$, where $h$ is constant. The motion of the fluid is assumed to be irrotational; therefore, we can introduce velocity potential, $\phi$, which satisfies the Laplace equation

$$\phi_{xx}+\phi_{yy}+\phi_{zz}=0\,,\quad\textrm{for}\;-h<z<\eta\left(x,y,t\right)\,.$$

(1)

At the free surface $z=\eta\left(x,y,t\right)$, the kinematic and dynamic boundary conditions read

$$\eta_{t}=\phi_{z}-\left(\eta_{x}\phi_{x}+\eta_{y}\phi_{y}\right)$$

(2)

and

$$\phi_{t}+\frac{1}{2}\left(\phi_{x}^{2}+\phi_{y}^{2}+\phi_{z}^{2}\right)+g\eta+\frac{P}{\rho_{w}}=p\left(x,y,t\right)\,,$$

(3)

where $g$ is the acceleration due to gravity, $P$ the restoring force due to elastic bending, and $p(x,y,t)$ the external pressure exerted on the ice sheet (if $p<0$, the pressure acts downwards). The pressure resulting from the linear elastic bending model takes the form

$$P=D\Delta^{2}\eta+\rho_{i}d\eta_{tt}\,,$$

(4)

where $D=\frac{Ed^{3}}{12\left(1-\nu^{2}\right)}$ is the flexural rigidity ($E$ is Young’s modulus and $\nu$ is the Poisson ratio), and $\Delta^{2}=\partial_{xxxx}+\partial_{yyyy}+2\partial_{xxyy}$ is the bi-Laplacian operator acting on horizontal variables. Hence the dynamic condition can be rewritten as

$$\phi_{t}+\frac{1}{2}\left(\phi_{x}^{2}+\phi_{y}^{2}+\phi_{z}^{2}\right)+g\eta+\frac{D}{\rho_{w}}\Delta^{2}\eta+\frac{\rho_{i}d}{\rho_{w}}\eta_{tt}=p\,.$$

(5)

Finally, the impermeability boundary condition at the bottom,

$$\phi_{z}=0\,,\quad\text{at $z=-h$}\,,$$

(6)

completes the whole system.

Ablowitz et al. [23] introduced an explicit non-local formulation of the classical equations of water waves in both two and three dimensions. Following their work, we will derive a Benney-Luke-type equation for hydroelastic waves in the shallow-water regime in the subsequent analyses.

We reformulate these equations in terms of the surface displacement and surface velocity potential. If we denote by $q(x,y,t)=\phi(x,y,\eta(x,y,t),t)$ the value of $\phi$ at the free surface, then the dynamic boundary condition can be recast to

$$q_{t}+\frac{1}{2}|\nabla{q}|^{2}+g\eta-\frac{\left(\eta_{t}+\nabla q\cdot\nabla\eta\right)^{2}}{2\left(1+|\nabla\eta|^{2}\right)}+\frac{D}{\rho_{w}}\Delta^{2}\eta+\frac{\rho_{i}d}{\rho_{w}}\eta_{tt}=p\,,$$

(7)

where $\nabla$ is the horizontal gradient operator. For arbitrary harmonic functions $\phi$ and $\psi$, it is straightforward to verify that

$$\left(\phi_{z}\psi_{x}+\psi_{z}\phi_{x}\right)_{x}+\left(\phi_{z}\psi_{y}+\psi_{z}\phi_{y}\right)_{y}+\left(\phi_{z}\psi_{z}-\psi_{x}\phi_{x}-\psi_{y}\phi_{y}\right)_{z}=0\,.$$

(8)

We choose $\psi=e^{\text{i}\left(k_{1}x+k_{2}y\right)+kz}$, where $k=\sqrt{k_{1}^{2}+k_{2}^{2}}$ is the magnitude of the wavenumber vector. Substituting this solution into Eq. (8), integrating over the space occupied by the fluid, and applying the divergence theorem, one obtains the global relation

	$\displaystyle 0=$	$\displaystyle\,\int_{\partial D}e^{\text{i}(k_{1}x+k_{2}y)+kz}\left[(\text{i}k_{1}\phi_{z}+k\phi_{x})n_{1}+(\text{i}k_{2}\phi_{z}+k\phi_{y})n_{2}\right.$		(9)
		$\displaystyle\,\left.+(k\phi_{z}-\text{i}k_{1}\phi_{x}-\text{i}k_{2}\phi_{y})n_{3}\right]\mathrm{d}S\,,$

where $(n_{1},n_{2},n_{3})^{\top}$ is the outward unit normal vector of $\partial D$, the boundary of the fluid body, and $\mathrm{d}S$ is the surface element. Using the kinematic boundary condition at the free surface and the impermeability condition at the bottom, Eq. (9) can be recast to

$\displaystyle\iint_{\mathbb{R}^{2}}e^{\text{i}(k_{1}x+k_{2}y)}\Big{[}e^{k(\eta+h)}(k\eta_{t}-\text{i}k_{1}q_{x}-\text{i}k_{2}q_{y})+(\text{i}k_{1}\phi_{x}+\text{i}k_{2}\phi_{y})_{z=-h}\Big{]}\mathrm{d}x\mathrm{d}y=0\,.$

(10)

In the same vein, substituting $\psi=e^{\text{i}(k_{1}x+k_{2}y)-kz}$ into Eq. (8) yields

$\displaystyle\iint_{\mathbb{R}^{2}}e^{\text{i}(k_{1}x+k_{2}y)}\Big{[}-e^{-k(\eta+h)}(k\eta_{t}+\text{i}k_{1}q_{x}+\text{i}k_{2}q_{y})+(\text{i}k_{1}\phi_{x}+\text{i}k_{2}\phi_{y})_{z=-h}\Big{]}\mathrm{d}x\mathrm{d}y=0\,.$

(11)

Subtracting Eq. (11) from Eq. (10) gives

$\displaystyle\iint_{\mathbb{R}^{2}}\left\{\text{i}\eta_{t}\cosh[k(\eta+h)]+(\textbf{k}\cdot\nabla{q})\frac{\sinh\left[k(\eta+h)\right]}{k}\right\}e^{\text{i}(k_{1}x+k_{2}y)}\mathrm{d}x\mathrm{d}y=0\,,$

(12)

where we denote by $\mathbb{k}=(k_{1},k_{2})^{\top}$ the wavenumber vector.

It is convenient to non-dimensionalize the problem using the Boussinesq scaling to proceed with the derivation. In this respect, we introduce the typical wavelength $l$, the small amplitude $a$, and the shallow-water speed $c_{0}=\sqrt{gh}$, and rescale the variables as

		$\displaystyle x=lx^{*}\,,\quad y=ly^{*}\,,\quad k_{1}=\frac{k_{1}^{*}}{l}\,,\quad k_{2}=\frac{k_{2}^{*}}{l}\,,$
		$\displaystyle t=\frac{lt^{*}}{c_{0}}\,,\quad q=\frac{gla}{c_{0}}q^{*}\,,\quad\eta=a\eta^{*}\,,\quad p=\frac{p^{*}ga^{2}}{h}\,.$

Then the dimensionless dynamic boundary condition and global relation read

$$q_{t}+\frac{\epsilon}{2}|\nabla q|^{2}+\eta-\frac{\epsilon\mu^{2}}{2}\frac{(\eta_{t}+\epsilon\nabla q\cdot\nabla\eta)^{2}}{1+\epsilon^{2}\mu^{2}|\nabla\eta|^{2}}+\delta\Delta^{2}\eta+\gamma\eta_{tt}=\epsilon{p}$$

(13)

and

$\displaystyle 0=\iint_{\mathbb{R}^{2}}e^{\text{i}(k_{1}x+k_{2}y)}\left\{\text{i}\eta_{t}\cosh\left[\mu k(1+\epsilon\eta)\right]+(\textbf{k}\cdot\nabla q)\frac{\sinh\left[\mu k(1+\epsilon\eta)\right]}{\mu k}\right\}\mathrm{d}x\mathrm{d}y\,,$

(14)

respectively, where asterisks have been dropped for the ease of notations, and the dimensionless parameters read

$$\epsilon=\frac{a}{h}\,,\quad\mu=\frac{h}{l}\,,\quad\delta=\frac{Ed^{3}}{12\rho_{w}g(1-\nu^{2})l^{4}}\,,\quad\gamma=\frac{{\rho_{i}}dh}{{\rho_{w}}l^{2}}\,.$$

Let $\epsilon,\delta,\gamma=O(\mu^{2})$ and this assumption corresponds to $a=0.1\,\mathrm{m}$, $d=1\,\mathrm{m}$, $h=10\,\mathrm{m}$, $l=100\,\mathrm{m}$, and $E=10^{9}\,\mathrm{N/m}^{2}$. It is noted that in the real situation, $\gamma=\frac{\rho_{i}}{\rho_{w}}\frac{d}{h}\mu^{2}$ is usually much smaller than $\mu^{2}$ due to the smallness of $d/h$; however, it does not qualitatively change the model equation if we involve the inertia effect in the current problem.

We expand the dynamic condition (13) and the global relation (14) with respect to small parameters and retain terms valid up to $O(\epsilon,\delta,\gamma,\mu^{2})$ to achieve a balance between dispersion and nonlinearity. After some algebra, one obtains the following Boussinesq-type equations

$$\eta=-q_{t}-\frac{\epsilon}{2}|\nabla q|^{2}-\delta\Delta^{2}\eta-\gamma\eta_{tt}+\epsilon{p}\,,$$

(15)

$$\left(1-\frac{\mu^{2}}{2}\Delta\right)\eta_{t}+\left(\Delta-\frac{\mu^{2}}{6}\Delta^{2}\right)q+\epsilon\left(\nabla q\cdot\nabla\eta\right)+\epsilon\eta\Delta q=0\,,$$

(16)

where $\Delta=\partial_{xx}+\partial_{yy}$ is the horizontal Laplace operator, and Eq. (16) is obtained by taking the inverse Fourier transform. By combining Eqs. (15) and (16) to eliminate $\eta$ and $\eta_{t}$, one arrives at

$$\left[1-\left(\frac{\mu^{2}}{2}\Delta+\delta\Delta^{2}\right)\right]q_{tt}-\left(\Delta-\frac{\mu^{2}}{6}\Delta^{2}\right)q-\gamma q_{tttt}+\epsilon\left(\partial_{t}|\nabla q|^{2}+q_{t}\Delta q\right)-\epsilon{p_{t}}=0\,.$$

(17)

Upon noticing $q_{tt}\sim\Delta{q}$ to leading order, Eq. (17) can be reformulated to

$$q_{tt}-\Delta{q}-\left(\frac{\mu^{2}}{3}+\gamma\right)\Delta^{2}q-\delta\Delta^{3}q+\epsilon\left(\partial_{t}|\nabla{q}|^{2}+q_{t}\Delta{q}\right)=\epsilon{p_{t}}\,,$$

(18)

a Benney-Luke-type equation with a time-dependent external force.

By letting $p=0$ in Eq. (18), we can derive the fifth-order Kadomtsev-Petviashvili (KP) equation, a single evolution equation describing uni-directional waves with weak variations in the transverse direction. For this purpose, we seek a solution of the form $q(x,y,t)=f(\bar{x},\bar{y},\bar{\tau})$, where $\bar{\tau}=\frac{1}{2}\mu^{2}t$, $\bar{x}=x-t$, and $\bar{y}=\mu y$. In this asymptotic limit, we can rewrite Eq. (18) as

$\displaystyle f_{\bar{x}\bar{\tau}}+f_{\bar{y}\bar{y}}+\left(\frac{1}{3}+\frac{\gamma}{\mu^{2}}\right)\partial_{\bar{x}}^{4}f+\frac{\delta}{\mu^{2}}\partial_{\bar{x}}^{6}f+\frac{3\epsilon}{\mu^{2}}f_{\bar{x}}f_{{\bar{x}{\bar{x}}}}=0$

upon neglecting higher-order terms. Differentiating the above equation with respect to $\bar{x}$ and substituting $\eta=f_{\bar{x}}+O\left(\mu^{2}\right)$, we obtain the fifth-order KP equation

$$\left[\eta_{\bar{\tau}}+\left(\frac{1}{3}+\frac{\gamma}{\mu^{2}}\right)\partial_{\bar{x}}^{3}\eta+\frac{\delta}{\mu^{2}}\partial_{\bar{x}}^{5}\eta+\frac{3\epsilon}{\mu^{2}}\eta\eta_{\bar{x}}\right]_{\bar{x}}+\eta_{\bar{y}\bar{y}}=0\,.$$

(19)

In the context of hydroelastic waves, the fifth-order KP equation (19) was first obtained by Hărăgus-Courcelle and Andrej [24] via a regular asymptotic expansion and re-derived by Guyenne and Părău [19] based on the Taylor expansion of the Dirichlet-Neumann operator in the Hamiltonian framework. If the $y$-dependence is negligible, Eq. (19) reduces to the fifth-order Korteweg-de Vries equation, which is known as a reduced model for hydroelastic waves in channels (see [25] for more details).

Finally, we show that without the external forcing, Eq. (18) has a Hamiltonian structure with the Hamiltonian functional

$\displaystyle\mathcal{H}[q,\xi]=\frac{1}{2}\iint_{\mathbb{R}^{2}}\left[\left(\xi-\frac{\epsilon}{2}|\nabla{q}|^{2}\right)^{2}+|\nabla{q}|^{2}+\delta|\nabla\Delta{q}|^{2}-\left(\frac{\mu^{2}}{3}+\gamma\right)(\Delta{q})^{2}\right]\mathrm{d}x\mathrm{d}y\,.$

(20)

The canonical variables are $q$ and $\xi$, where $\xi$ is defined as

$$\xi=q_{t}+\frac{\epsilon}{2}|\nabla{q}|^{2}\,.$$

(21)

Calculating the variational derivatives with respect to the canonical variables gives

$\displaystyle\frac{\delta\mathcal{H}}{\delta\xi}=\xi-\frac{\epsilon}{2}|\nabla{q}|^{2}=q_{t}$

(22)

and

	$\displaystyle\frac{\delta\mathcal{H}}{\delta{q}}=$	$\displaystyle\,\epsilon\nabla\cdot\left[\left(\xi-\frac{\epsilon}{2}|\nabla{q}|^{2}\right)\nabla{q}\right]-\Delta{q}-\left(\frac{\mu^{2}}{3}+\gamma\right)\Delta^{2}{q}-\delta\Delta^{3}{q}$		(23)
	$\displaystyle=$	$\displaystyle\,\epsilon\left[q_{t}\Delta{q}+\frac{1}{2}\left(|\nabla{q}|^{2}\right)_{t}\right]-\Delta{q}-\left(\frac{\mu^{2}}{3}+\gamma\right)\Delta^{2}{q}-\delta\Delta^{3}{q}$
	$\displaystyle=$	$\displaystyle\,-q_{tt}-\frac{\epsilon}{2}\left(|\nabla{q}|^{2}\right)_{t}=-\xi_{t}\,.$

Traditionally, weakly nonlinear wavepacket dynamics are investigated using the cubic nonlinear Schrödinger (NLS) equation, which describes the envelope evolution of a monochromatic carrier wave in deep water. In the case of finite depth, the wave envelope is strongly coupled with the induced mean flow, and their interactions are governed by the Benney-Roskes-Davey-Stewartson (BRDS) system. It is well known that the envelope equations play a key role in understanding the existence of solitary wave solutions and their stability properties for both gravity-capillary and flexural-gravity problems (see [16, 26] for example). In the context of 3D hydroelastic waves in the shallow-water regime, we expect to gain some theoretical predictions from the BRDS system derived from Eq. (18) in the absence of external forcing. For this purpose, we expand the solution of Eq. (18) as

$$q=\alpha{q_{0}}+\alpha^{2}{q_{1}}+\alpha^{3}{q_{2}}+\cdots\,,$$

(24)

where the assumption of small amplitude is taken through the expansion in powers of $\alpha$. Then the equation to leading order reads

$$q_{0tt}-\Delta{q_{0}}-\left(\frac{\mu^{2}}{3}+\gamma\right)\Delta^{2}{q_{0}}-\delta\Delta^{3}{q_{0}}=0\,,$$

whose real solution is given by

$$q_{0}=A(X,Y,T)e^{\text{i}\theta}+\text{c.c.}+M(X,Y,T)\,,$$

(25)

where c.c. stands for complex conjugation, $T=\alpha{t}$, $X=\alpha{x}$ and $Y=\alpha{y}$ are slow variables, $\theta=kx-\omega{t}$ is the fast variable with the dispersion relation $\omega^{2}=k^{2}-\left(\frac{\mu^{2}}{3}+\gamma\right)k^{4}+\delta{k^{6}}$, and $M(X,Y,T)$ is a real slow-varying, mean term.

We substitute $\partial_{t}=-\omega\partial_{\theta}+\alpha\partial_{T}$, $\partial_{x}=k\partial_{\theta}+\alpha\partial_{X}$, and $\partial_{y}=\alpha\partial_{Y}$ into the Benney-Luke-type equation. At second order, one obtains

$\displaystyle\left[2\text{i}\omega{A_{T}}+2\text{i}kA_{X}-4\text{i}k^{3}\left(\frac{\mu^{2}}{3}+\gamma\right)A_{X}+6\text{i}\delta{k^{5}}A_{X}\right]e^{\text{i}\theta}=0\,,$

(26)

and solve $q_{1}$ for what remains, i.e., solve

$$\mathcal{L}q_{1}=-3\text{i}\epsilon{k^{2}}\omega{A^{2}}e^{2i\theta}+\text{c.c.}$$

with $\mathcal{L}=\left(\omega^{2}-k^{2}\right)\partial^{2}_{\theta}-\left(\frac{\mu^{2}}{3}+\gamma\right)k^{4}\partial^{4}_{\theta}-\delta{k^{6}\partial^{6}_{\theta}}$. The former equation implies that $A(X,Y,T)=A(X-c_{\text{g}}T,Y)$, where $c_{\text{g}}$ is called the group velocity defined as

$$c_{\text{g}}=\frac{\rm{d}\omega}{\rm{d}k}=\frac{1}{\omega}\left[k-2\left(\frac{\mu^{2}}{3}+\gamma\right)k^{3}+3\delta k^{5}\right]\,.$$

The latter equation can be solved by letting $q_{1}=Be^{2\text{i}\theta}+\text{c.c.}$, which gives

$$B=\frac{-3\text{i}\epsilon\omega{k^{2}A^{2}}}{-4\omega^{2}+4k^{2}-16\left(\frac{\mu^{2}}{3}+\gamma\right)k^{4}+64\delta{k^{6}}}\,.$$

To obtain the BRDS system, we transform to a moving reference frame $X-c_{\text{g}}T$ and denote by $\tau=\alpha{T}$ a much slower time scale. Then to the third order, we can get

$$\left(1-c_{\text{g}}^{2}\right)M_{XX}+M_{YY}=\chi(|A|^{2})_{X}\,,$$

(27)

$$\text{i}A_{\tau}+\frac{\omega^{{}^{\prime\prime}}}{2}A_{XX}+\frac{c_{\text{g}}}{2k}A_{YY}-\epsilon{kAM_{X}}+\Xi|A|^{2}A=0\,,$$

(28)

where

$$\chi=-\epsilon(k^{2}c_{\text{g}}+2k\omega)\,,\qquad\Xi=\frac{3\epsilon^{2}\omega}{-4\left(\frac{\mu^{2}}{3}+\gamma\right)+20\delta{k^{2}}}\,.$$

The BRDS system (27)–(28) is valid for arbitrary $k>0$; however, in the current context, we focus on the case $k=k_{\text{min}}$, where $k_{\min}$ is the wavenumber corresponding to the phase speed minimum. Since $c_{\text{g}}^{2}<1$, $\omega^{{}^{\prime\prime}}>0$, and $c_{\text{g}}/2k>0$, the system is of ‘elliptic-elliptic’ type as the coefficients of $M_{XX}$, $M_{YY}$, $A_{XX}$, and $A_{YY}$ are of the same sign. Following [27], the system’s behavior then depends on the sign of the parameter $\Theta$ defined as

$$\Theta=\Xi-\frac{\epsilon{k}\chi}{1-c_{\text{g}}^{2}}\,.$$

It was proved by Cipolatti [27] that $\Theta>0$ is a sufficient condition for the BRDS system to have a localized ground state solution, and the system is therefore called ‘focussing’. For the case that transverse modulations are absent (i.e., $\partial_{YY}=0$), the BRDS system (27)–(28) degenerates to the one-dimensional NLS equation $\text{i}A_{\tau}+\frac{\omega^{\prime\prime}}{2}A_{XX}+\Theta|A|^{2}A=0$. It is well known that the simplest solution to the one-dimensional focussing NLS equation takes an exact formula: $A=\widetilde{a}\,\text{sech}(\widetilde{b}X)e^{\text{i}\Omega\tau}$ with $\widetilde{a}=\pm\sqrt{2\Omega/\Theta}$ and $\widetilde{b}=\sqrt{2\Omega/\omega^{\prime\prime}}$.

In infinite depth, the effect of the mean flow vanishes, and the evolution is governed by the two-dimensional NLS equation

$$\text{i}A_{\tau}+A_{XX}+A_{YY}+|A|^{2}A=0\,,$$

(29)

where coefficients have all been normalized to unity. It was shown by Alfimov et al. [28] that for the nonlinear eigenvalue problem of Eq. (29),

$$\left\{\begin{aligned} &\rho_{XX}+\rho_{YY}-\rho+\rho^{3}=0\,,\\ &\rho(X,Y+l)=\rho(X,Y)\,,\\ &\lim_{X\to\pm\infty}\rho(X,Y)=0\,,\end{aligned}\right.$$

there exists a one-parameter family of solutions $\rho_{l}(X,Y)$ for $l>2\pi/\sqrt{3}$, where $l$ represents the basic period in the $Y$-direction. In contrast to plane solitons, this new type of solution features a nontrivial variation in the transverse direction. As the basic period $l\rightarrow{2\pi/\sqrt{3}}+0$, the solution $\rho_{l}$ degenerates to the plane soliton solution, and hence the emergence of transversally periodic solitons is called a dimension-breaking bifurcation. Following the new branch by increasing the bifurcation parameter $l$, a fully localized solution with circular symmetry can be obtained as $l\to\infty$; that is to say, these new solutions are a bridge between plane solitons and lumps. Inspired by the results of the asymptotic model, Milewski and Wang [29] found transversally periodic solitary waves in the full gravity-capillary wave problem corresponding to its NLS counterpart. Additionally, these authors predicted that the threshold wavelength of transverse instability of a plane solitary wave coincides with the dimension-breaking bifurcation point.

Motivated by their work, we numerically explore various types of localized steady solutions to the BRDS system through the dimension-breaking bifurcation, which provide a solid underpinning for the existence of solitary waves in the Benney-Luke-type equation. Without loss of generality, we set $\mu=0.3$, $\epsilon=\mu^{2}$, $\gamma=\delta=\mu^{3}$, and as a result, the BRDS system is of focussing type since the coefficient $\Theta>0$. To find localized time harmonic solutions to the BRDS system, we assume $A=e^{\text{i}\Omega\tau}\rho(X,Y)$, and then the steady BRDS system in the non-local form can be written as

$\displaystyle-\Omega\rho+\frac{\omega^{{}^{\prime\prime}}}{2}\rho_{XX}+\frac{c_{\text{g}}}{2k}\rho_{YY}-\vartheta\epsilon{k}\chi\widetilde{\Delta}^{-1}\left[\rho^{2}\right]_{XX}\rho+\Xi\rho^{3}=0\,,$

(30)

where $\widetilde{\Delta}=\left(1-c_{\text{g}}^{2}\right)\partial_{XX}+\partial_{YY}$ and $\vartheta=1$. The parameter $\vartheta$ varying from 0 to 1 is introduced to assist with computing steady solutions to Eq. (30). The numerical process can be divided into three steps in general. First of all, by letting $\vartheta=0$, the problem can be reduced to the steady focussing cubic NLS equation. This equation, after normalization, has countably many localized radial solutions (see [30] for example), but we focus only on computing the ground state solution. Secondly, the obtained ground state solution is taken as the initial guess, and the non-local term in Eq. (30) is included by using a numerical continuation in $\vartheta$ up to $\vartheta=1$, resulting in a lump solution to the BRDS system. The lump is computed using the standard pseudo-spectral method in a horizontally double periodic domain with small wavenumbers $k_{1}$ and $k_{2}$ in the $X$- and $Y$-direction, respectively. Finally, the dimension-breaking bifurcation curve can be obtained by using $k_{2}$ as the continuation parameter, and the interested readers are referred to [16] for more details on the numerical scheme.

By gradually increasing $k_{2}$ (equivalently, decreasing the basic period $\bar{l}$ in the $Y$-direction), a one-parameter family of solutions $\rho_{\bar{l}}(X,Y)$ for $\bar{l}>2\pi/1.205$ is obtained. As the basic period $\bar{l}\to{2\pi/1.205+0}$, $\rho_{\bar{l}}$ degenerates to a plane soliton, while $\bar{l}\to\infty$, a lump solution emerges. The intermediate states which bridge the two limiting cases correspond to transversally periodic solitons. Fig. 1 shows the typical wave profiles on this new branch: 1 and 1 are the two limiting cases, while 1 is the intermediate case. It is noted that the bifurcation wavenumber $k_{2}\approx 1.205$ is closely related to the transverse instability of the plane soliton.

As shown in Fig. 1, three types of solitons (plane soliton, transversally periodic soliton, and lump) have been found in the BRDS system. Therefore we can expect their counterparts in the Benney-Luke-type equation (Eq. (18) without external forcing). To find solitary waves and the dimension-breaking bifurcation diagram, we start with computing lumps in the Benney-Luke-type equation. The numerical scheme is an extension of Petviashvili’s method [31], whose essence is to perform iterations in the Fourier space supplemented by a normalization factor. For this purpose, we first assume a lump is moving in an arbitrary horizontal direction, namely

$$q(x,y,t)=q(x-c_{1}t,y-c_{2}t)\,,$$

where $c_{1}$ and $c_{2}$ are constants. Substituting this ansatz into Eq. (18) and performing the Fourier transform, one obtains

$$\hat{q}=\epsilon\frac{\text{i}(c_{1}k_{1}+c_{2}k_{2})\widehat{|\nabla{q}|^{2}}+c_{1}\widehat{q_{x}\Delta{q}}+c_{2}\widehat{q_{y}\Delta{q}}}{D}=\mathcal{P}\left[\hat{q}\right]\,,$$

(31)

where

$\displaystyle D=k_{1}^{2}+k_{2}^{2}-\left(c_{1}k_{1}+c_{2}k_{2}\right)^{2}-\left(\frac{\mu^{2}}{3}+\gamma\right)\left(k_{1}^{2}+k_{2}^{2}\right)^{2}+\delta\left(k_{1}^{2}+k_{2}^{2}\right)^{3}\,,$

(32)

and hats denote the Fourier transform in $x$ and $y$. Following Ablowitz et al. [23], we introduce a multiplier in every iteration to prevent unlimited growth or reduction in amplitude, and hence the numerical scheme can be proposed as

$$\hat{q}_{n+1}=\alpha_{n}\mathcal{P}[\hat{q}_{n}]\,,\quad\text{with }\alpha_{n}=\dfrac{\displaystyle{\iint}{|\hat{q}_{n}|^{2}\mathrm{d}k_{1}\mathrm{d}k_{2}}}{\displaystyle{\iint}{{\hat{q}_{n}}^{*}}\mathcal{P}[\hat{q}_{n}]\mathrm{d}k_{1}\mathrm{d}k_{2}}\,.$$

(33)

In most numerical experiments, we compute solitary waves propagating in the positive $x$-direction, that is, $c_{2}=0$ and $c_{1}=c$. In Fig. 2, we present bifurcation curves and typical wave profiles of solitary waves in the Benney-Luke-type equation for $\mu=0.3$, $\epsilon=\mu^{2}$, and $\gamma=\delta=\mu^{3}$. The bifurcation diagrams 2 and 2 show the relation between the translating speed $c$ and the center of the free-surface displacement for plane solitary waves (locally confined in the direction of wave propagation without variations in the transverse direction) and lumps (localized in all horizontal directions), respectively. It is shown that both plane solitary waves and lumps bifurcate from $c\approx 0.9848$, and only depression waves, which feature a negative free-surface elevation at their center, can be found for positive $c$. The dashed lines in 2 and 2 represent the bifurcation curves that take the leading-order approximation of amplitude $\eta\sim-{q_{t}}$ in the BRDS system, indicating that the theoretical predictions are valid only for a very narrow range of speed in the vicinity of the bifurcation point. The dimension-breaking bifurcation curves are plotted in Fig. 2, demonstrating the relationship between the basic wavenumber in the transverse direction and $\eta(0,0)$ for fixed $c=0.9798$, $0.9811$, $0.9823$ (from right to left). Typical examples of three types of solitary waves are plotted on the right column: plane solitary wave (2), lump (2), and transversally periodic solitary wave, which are localized in the propagation direction and periodic in the transverse direction (2).

In the gravity-capillary (GC) wave problem, the KP-I equation proposed for unidirectional waves with weak transverse variation has rational solutions, which can help to understand the algebraic decay of lumps in the Benney-Luke equation with strong surface tension. However, regarding hydroelastic waves, no analytic travelling-wave solutions are known to the fifth-order KP equation. Alternatively, we manage to understand the far-field behavior of lumps of the Benney-Luke-type equation with the aid of the associated BRDS system. Using the envelope equations, Kim and Akylas [26] showed that GC lumps feature algebraically decaying tails at infinity regardless of water depth owing to the algebraic decay of the induced mean flow. Following their argument, we first take the Fourier transform of the mean flow equation (27), which yields

$\displaystyle\frac{\partial{M}}{\partial{X}}=\chi\iint_{\mathbb{R}^{2}}\widehat{|A|^{2}}\,\frac{k_{1}^{2}\,e^{\text{i}(k_{1}X+k_{2}Y)}}{(1-c_{\text{g}}^{2})k_{1}^{2}+k_{2}^{2}}\,\mathrm{d}k_{1}\mathrm{d}k_{2}\,.$

where $\widehat{|A|^{2}}$ is the Fourier transform of $|A|^{2}$ in $X$ and $Y$. Therefore, as $X^{2}+Y^{2}\to\infty$,

$\displaystyle\frac{\partial{M}}{\partial{X}}\sim-\frac{I_{0}\chi}{2\pi\sqrt{1-c_{\text{g}}^{2}}}\frac{\partial}{\partial{Y}}\left[\frac{Y}{X^{2}+(1-c_{\text{g}}^{2})Y^{2}}\right]\,,$

where

$$I_{0}=\iint_{\mathbb{R}^{2}}|A|^{2}\mathrm{d}X\mathrm{d}Y\,.$$

This indicates that the induced mean flow decays algebraically at infinity, and the envelope $A$ of primary harmonic features exponential decaying tails according to Eq. (28). As a result, the induced mean flow controls the behavior at the tails of a lump. The free surface elevation at the tails is obtained $\eta\sim\alpha^{2}cM_{X}$ and decays algebraically as well,

$\displaystyle\eta\sim-\frac{\alpha^{2}cI_{0}\chi}{2\pi\sqrt{1-c_{\text{g}}^{2}}}\frac{\partial}{\partial{Y}}\left[\frac{Y}{X^{2}+(1-c_{\text{g}}^{2})Y^{2}}\right]\,.$

The numerical solutions to the Benney-Luke-type equation are presented in Fig. 3 to confirm the asymptotic results from the BRDS system. The plots of $\ln|\eta|$ against $\ln|X|$ and $\ln|Y|$ for the lump with the speed $c=0.983$ are shown. Note that, in both plots, the free-surface elevation approaches a straight line with slope $-2$ to a good approximation, which is consistent with the prediction that the tails ultimately decay algebraically like $X^{-2}$ and $Y^{-2}$ in the $x$- and $y$-direction, respectively.

The stability of line solitary waves subject to transverse long-wavelength perturbations was explored extensively in gravity-capillary water-wave problems by different groups. Kim and Akylas [32] provided a long wave stability analysis and concluded that the leading-order instability growth rate was associated with the change rate of the mechanical energy relative to the wave speed. However, their argument could not give the threshold wavelength of the transverse instability. On the other hand, the linear stability analysis of plane solitary waves of the 2D NLS equation indicated the critical wavenumber of transverse instability [33]. Motivated by this result, Akers and Milewski [34] and Wang and Milewski [30] obtained the threshold wavelength of perturbations and performed direct numerical simulations to verify the theoretical prediction. The studies mentioned above are related to depression GC solitary waves, while the elevation ones were discussed by Wang and Vanden-Broeck [35], and similar results were obtained. In the subsequent analyses, we first perform an asymptotic analysis of the instability of plane solitary waves subject to transverse long-wavelength perturbations and then specify the critical wavenumber of triggering instability.

We assume that $U({\sigma};c)$ is a plane solitary wave solution characterized by the translating speed $c$, where $\sigma=x-ct$. Substituting this solution into the unforced Benney-Luke-type equation yields

$$(c^{2}-1)U_{\sigma\sigma}-\left(\frac{\mu^{2}}{3}+\gamma\right)\partial_{\sigma}^{4}U-\delta\partial_{\sigma}^{6}U-3c\epsilon U_{\sigma}U_{\sigma\sigma}=0\,.$$

(34)

Taking the derivative of Eq. (34) with respect to $c$ yields

$$-2cU_{\sigma\sigma}+3\epsilon{U_{\sigma}U_{\sigma\sigma}}=\left(c^{2}-1\right)U_{c\sigma\sigma}-\left(\frac{\mu^{2}}{3}+\gamma\right)\partial_{\sigma}^{4}U_{c}-\delta\partial_{\sigma}^{6}U_{c}-3c\epsilon\left({U_{\sigma}U_{c\sigma}}\right)_{\sigma}\,.$$

(35)

We perturb the plane solitary wave solution in the transverse direction using a cosine longwave, namely $q=U(\sigma;c)+\tilde{q}(\sigma;c)e^{\text{i}\beta{y}+\lambda{t}}+\text{c.c.}$ with a small wavenumber $\beta$ in the direction transverse to wave propagation. Substituting the ansatz into the Benney-Luke-type equation and collecting the coefficients of $e^{\text{i}\beta{y}+\lambda{t}}$, one can obtain

	$\displaystyle 0=$	$\displaystyle\,\left[\lambda^{2}+\beta^{2}-\left(\frac{\mu^{2}}{3}+\gamma\right)\beta^{4}+\delta\beta^{6}+\epsilon\beta^{2}cU_{\sigma}+\epsilon\lambda{U_{\sigma\sigma}}\right]\tilde{q}+\left(2\epsilon\lambda{U_{\sigma}}-3c\epsilon{U_{\sigma\sigma}}-2c\lambda\right){\tilde{q}}_{\sigma}$		(36)
		$\displaystyle\,+\left[c^{2}-1+2\beta^{2}\left(\frac{\mu^{2}}{3}+\gamma\right)-3\delta\beta^{4}-3c\epsilon{U_{\sigma}}\right]\tilde{q}_{\sigma\sigma}+\left[3\delta\beta^{2}-\left(\frac{\mu^{2}}{3}+\gamma\right)\right]\partial_{\sigma}^{4}\tilde{q}-\delta\partial_{\sigma}^{6}\tilde{q}\,.$

Furthermore, we assume

	$\displaystyle\lambda$	$\displaystyle=$	$\displaystyle\beta\lambda^{(1)}+\beta^{2}\lambda^{(2)}+\beta^{3}\lambda^{(3)}+\cdots\,,$		(37a)
	$\displaystyle\tilde{q}$	$\displaystyle=$	$\displaystyle q^{(0)}+\beta{q^{(1)}}+\beta^{2}{q^{(2)}}+\cdots\,.$		(37b)

By substituting expansions (37a) and (37b) into Eq. (36) and equating like powers of $\beta$, we obtain, at O(1),

$$0=\left(c^{2}-1\right)q^{(0)}_{\sigma\sigma}-\left(\frac{\mu^{2}}{3}+\gamma\right)\partial_{\sigma}^{4}q^{(0)}-\delta\partial_{\sigma}^{6}q^{(0)}-3c\epsilon\left(U_{\sigma}q^{(0)}_{\sigma}\right)_{\sigma}\,.$$

(38)

It is evident that $q^{(0)}=U_{\sigma}$ is a solution to Eq. (38). To $O(\beta)$, $q^{(1)}$ satisfies the forced problem

$$\lambda^{(1)}\left(2cU_{\sigma\sigma}-3\epsilon{U_{\sigma}U_{\sigma\sigma}}\right)=\left(c^{2}-1\right)q^{(1)}_{\sigma\sigma}-\left(\frac{\mu^{2}}{3}+\gamma\right)\partial_{\sigma}^{4}{q^{(1)}}-{\delta}\partial_{\sigma}^{6}{q^{(1)}}-3c\epsilon\left(U_{\sigma}q^{(1)}_{\sigma}\right)_{\sigma}\,.$$

(39)

Recalling (35), the solution to Eq. (39) is $q^{(1)}=-\lambda^{(1)}U_{c}$. At the next order, the equation for $q^{(2)}$ reads

		$\displaystyle\,\left[\left(c^{2}-1\right)\partial_{\sigma}^{2}-\left(\frac{\mu^{2}}{3}+\gamma\right)\partial_{\sigma}^{4}-{\delta}\partial_{\sigma}^{6}\right]{q^{(2)}}-3c\epsilon\left(U_{\sigma}q^{(2)}_{\sigma}\right)_{\sigma}$		(40)
	$\displaystyle=$	$\displaystyle\,-\left({\lambda^{(1)}}^{2}+c\epsilon{U_{\sigma}}+1+\epsilon\lambda^{(2)}U_{\sigma\sigma}\right)q^{(0)}-\epsilon\lambda^{(1)}U_{\sigma\sigma}q^{(1)}-2\epsilon\lambda^{(2)}U_{\sigma}q^{(0)}_{\sigma}$
		$\displaystyle\,-2\epsilon\lambda^{(1)}U_{\sigma}q^{(1)}_{\sigma}+2c\lambda^{(1)}q^{(1)}_{\sigma}+2c\lambda^{(2)}q^{(0)}_{\sigma}-2\left(\frac{\mu^{2}}{3}+\gamma\right)q^{(0)}_{\sigma\sigma}-3{\delta}\partial_{\sigma}^{4}q^{(0)}\,.$

The adjoint operator of the left-hand side of Eq. (40) is

$$\mathcal{L^{+}}=(c^{2}-1)\partial^{2}_{\sigma}-\left(\frac{\mu^{2}}{3}+\gamma\right)\partial^{4}_{\sigma}-\delta\partial^{6}_{\sigma}-3c\epsilon(U_{\sigma}\partial_{\sigma})_{\sigma}\,,$$

(41)

and it is obvious that $\mathcal{L^{+}}U_{\sigma}=0$. If we denote by $R$ the right-hand side of (40), then the solvability condition for the inhomogeneous equation (40) is

$$\int_{-\infty}^{\infty}{R}U_{\sigma}\mathrm{d}\sigma=0\,,$$

which gives

$$\int_{-\infty}^{\infty}\left[U_{\sigma}^{2}+c\epsilon{U_{\sigma}^{3}}-2\left(\frac{\mu^{2}}{3}+\gamma\right)U_{\sigma\sigma}^{2}+3\delta{U_{\sigma\sigma\sigma}^{2}}\right]\mathrm{d}\sigma=\left({\lambda^{(1)}}\right)^{2}\int_{-\infty}^{\infty}\left(\frac{1}{2}\epsilon{U_{\sigma}^{3}-cU_{\sigma}^{2}}\right)_{c}\mathrm{d}\sigma\,.$$

(42)

The integrals on the right and left sides are denoted by $M_{1}$ and $M_{2}$, respectively. We infer from Eq. (42) that the plane solitary wave $U$ is transversely unstable if

$$M_{1}M_{2}>0\,.$$

By letting

$$E=\int_{-\infty}^{\infty}\left(cU_{\sigma}^{2}-\frac{1}{2}\epsilon{U_{\sigma}^{3}}\right)\mathrm{d}\sigma\,,$$

the values of $E$ and $M_{2}$ can be calculated numerically. As shown in Fig. 4, $E$ and $M_{2}$ are both positive and decrease monotonically along with the translating speed $c$. Consequently, the instability condition is satisfied and plane solitary waves are proved to be transversally unstable under long-wavelength disturbances.

It is mentioned in §2 that the critical wavenumber for the onset of transverse instability of a plane solitary wave is coincident with the dimension-breaking bifurcation point. It is natural to establish the relation between the amplitude of plane solitary waves and the critical perturbation wavenumber $k_{\mathrm{c}}$ for the transverse instability in the Benney-Luke-type equation. It was shown in the last paragraph of §2 that for the BRDS plane solitons, when the wavenumber in the $Y$-direction increases to 1.205, denoted as $K_{\mathrm{c}}$, the dimension-breaking phenomenon occurs. Using the relations $k_{\mathrm{c}}\approx\alpha{K_{\mathrm{c}}}$ and $||\eta||_{\infty}\approx 2\alpha\omega|A|$ and eliminating $\alpha$ by using wave amplitude, we then obtain

$\displaystyle k_{\mathrm{c}}\approx 0.4167||\eta||_{\infty}\,.$

(43)

On the other hand, the critical wavenumber for the onset of transverse instability in the Benney-Luke-type equation can be obtained by direct numerical computations. We start with a lump solution and then gradually shrink the computational domain in the $y$-direction (by fixing the wave speed) until the variations in the transverse direction become trivial. Fig. 5 shows the relation between the amplitude of plane solitary waves and the dimension-breaking bifurcation point (solid line), indicating that the asymptotic prediction (dashed line) given by Eq. (43) is valid only for small-amplitude waves. It then follows that the parameter space ($||\eta||_{\infty}$, $k_{2}$) shown in Fig. 5 can be divided into two regions: below and above the solid curve. For a given plane solitary wave (which can be determined by its amplitude), the transverse instability occurs when it is superimposed with a harmonic perturbation, orthogonal to the direction of wave propagation, with the wavenumber in the lower region. In contrast, it is stable when the perturbation wavenumber is in the upper region.