2.1. Problem definition

We consider ocean surface waves propagating on deep water in the framework of potential-flow theory, thereby assuming incompressible inviscid flows, irrotational fluid motions and negligible effects of surface tension. A Cartesian coordinate system is chosen with the undisturbed water surface located at $z = 0$. The system can be described as a boundary value problem governed by the Laplace equation,

(2.1)\begin{equation} \nabla^2_3 \varPhi =0\quad \textrm{for}\ -\infty< z< \zeta(\boldsymbol{x},t), \end{equation}

where $\varPhi (\boldsymbol {x},z,t)$ denotes the velocity potential, $\zeta (\boldsymbol {x},t)$ is the free surface elevation, $\boldsymbol {x}$ is the position vector in the horizontal plane, $t$ is the time and $\nabla _3 = (\boldsymbol {\nabla },\partial _z)$, with $\boldsymbol {\nabla }=(\partial _x,\partial _y)$ denoting the gradient in the horizontal plane. Equation (2.1) should be solved subject to the nonlinear kinematic and dynamic boundary conditions (cf. Davey & Stewartson Reference Davey and Stewartson1974) at the free water surface $z= \zeta (\boldsymbol {x},t)$, as follows:

(2.2a)\begin{equation} \partial_t\zeta +\boldsymbol{\nabla} \varPhi\boldsymbol{\cdot} {\boldsymbol{\nabla}}\zeta - \partial_z\varPhi = 0 \end{equation}

and

(2.2b)\begin{equation} \varGamma \varPhi + \partial_t (\nabla_3\varPhi)^2+ \tfrac{1}{2} \nabla_3\varPhi\boldsymbol{\cdot} \nabla_3 (\nabla_3\varPhi)^2 = 0, \end{equation}

where the operator $\varGamma$ is defined as $\varGamma = \partial _{tt} +g\partial _z$ with $g$ the gravitational acceleration; a deep-water boundary condition is

(2.3)\begin{equation} \partial_z \varPhi =0\quad \textrm{for}\ z\to -\infty . \end{equation}

2.2. Stokes expansion and separation of harmonics

In order to solve the boundary value problem (2.1)–(2.3), we seek the solutions of unknown $\varPhi$ and $\zeta$ in a form of power series in wave steepness defined as $\epsilon = k_0A_0$ (a so-called Stokes expansion), with $k_0$ and $A_0$ denoting the characteristic wavenumber and wave amplitude, respectively,

(2.4a)\begin{equation} \varPhi = \epsilon \varPhi^{(1)} + \epsilon^2 \varPhi^{(2)}+ \epsilon^3 \varPhi^{(3)} + O(\epsilon^4) \end{equation}

and

(2.4b)\begin{equation} \zeta = \epsilon \zeta^{(1)} + \epsilon^2 \zeta^{(2)}+ \epsilon^3\zeta^{(3)} + O(\epsilon^4), \end{equation}

where we consider the first three orders and the superscripts denote the order in $\epsilon$. Substituting (2.4) into the boundary value problem (2.1)–(2.3) leads to the decomposition of the fully nonlinear system into different problems through a collection of the terms at the same order in $\epsilon$. The decomposed problems can be solved successively from the lowest to higher orders, as presented in the following from the first (§ 2.3) to the third order (§ 5).

Let linear surface elevation $\zeta ^{(1)}$ be expressed in two equivalent forms, as follows:

(2.5a)\begin{equation} \zeta^{(1)}(\boldsymbol{x},t) = \textstyle\frac{1}{2} A(\boldsymbol{x},t)\,\mathrm{e}^{\mathrm{i} (\boldsymbol{k}_0\cdot\boldsymbol{x}-\omega_0t)} +\text{c.c.} \end{equation}

or

(2.5b)\begin{equation} \zeta^{(1)}(\boldsymbol{x},t) = \frac{1}{2}\int_{-\infty}^{\infty} \hat{\zeta}(\boldsymbol{k}) \,\mathrm{e}^{\mathrm{i} (\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}-\omega(\boldsymbol{k})t)}\,\mathrm{d} \boldsymbol{k}^2 +\text{c.c.}, \end{equation}

in which $\text {c.c.}$ denotes the complex conjugates; $A$ denotes the complex wave envelope of the surface displacement of the carrier wave, with $\boldsymbol {k}_0=(k_0,0)$ the wavenumber vector chosen in the $x$ direction; $\omega _0$ denotes the angular frequency of the carrier wave that satisfies the dispersion relation $\omega _0=\omega (\boldsymbol {k}_0)$ given by $\omega (\boldsymbol {k}) = \sqrt {g k}$ (where $k=|\boldsymbol {k}|$); $\boldsymbol {k}=(k_x,k_y)$ denotes a wavevector in the horizontal plane. Equation (2.5a) denotes the first-order elevation being expressed in an envelope-type form and (2.5b) a form through linear superposition of monochromatic waves in the Fourier $\boldsymbol {k}$ plane. This paper focuses on the former which leads to a nonlinear evolution equation for $A$ that allows for the relaxation of a narrow banded assumption (cf. Chu & Mei Reference Chu and Mei1971; Davey & Stewartson Reference Davey and Stewartson1974).

Equating (2.5a) and (2.5b) we obtain a relation between $A$ and $\hat {\zeta }$, as follows:

(2.6)\begin{equation} A(\boldsymbol{x},t) = \int_{-\infty}^{\infty} \hat{\zeta}(\pmb{\kappa}+\boldsymbol{k}_0)\, \mathrm{e}^{\mathrm{i}[\pmb{\kappa}\boldsymbol{\cdot}\boldsymbol{x} -(\omega(\pmb{\kappa}+\boldsymbol{k}_0) -\omega_0)t]}\,\mathrm{d} \pmb{\kappa}^2, \end{equation}

where the integral variable was replaced with $\pmb {\kappa }$ through $\boldsymbol {k} = \pmb {\kappa }+\boldsymbol {k}_0$. Introducing a Fourier transform for $A$ defined as follows:

(2.7a)\begin{equation} A(\boldsymbol{x},t) = \int_{-\infty}^{\infty} \hat{A}(\boldsymbol{k},t)\,\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d} \boldsymbol{k}, \end{equation}

we obtain

(2.7b)\begin{equation} \hat{A}(\boldsymbol{k},t) = \hat{\zeta}(\boldsymbol{k}+\boldsymbol{k}_0)\,\mathrm{e}^{-\mathrm{i}(\omega(\boldsymbol{k}+\boldsymbol{k}_0)-\omega_0)t}. \end{equation}

A linear evolution equation for $A$ for waves of a broad bandwidth is derived by Trulsen et al. (Reference Trulsen, Kliakhandler, Dysthe and Velarde2000) in two equivalent forms, as follows:

(2.8a)\begin{equation} \partial_tA + \mathrm{i} \omega_0 \mathcal{L} A = 0, \end{equation}

with

(2.8b)\begin{equation} \mathcal{L}(\partial_x,\partial_y) = \left[ \left(1-\frac{\mathrm{i}\partial_x}{k_0}\right)^2 +\frac{(\mathrm{i}\partial_y)^2}{k^2_0} \right]^{1/4}-1 \end{equation}

and

(2.8c)\begin{equation} \partial_t \hat{A} + \mathrm{i}(\omega(\boldsymbol{k}+\boldsymbol{k}_0)-\omega_0)\hat{A}= 0, \end{equation}

where (2.8a) is used for constructing the solutions in §§ 3, 4 and 6.2; (2.8c) is in a form convenient for numerical implementations, as explained in § 4.4. It is clear that a narrow-banded assumption implies $k_0\mathcal {L} A\ll O(\epsilon )$, which, again, this paper aims to drop. For convenience and later reference, we introduce a dimensionless bandwidth parameter defined, as follows:

(2.9a,b)\begin{equation} \delta_{x,y}\sim \frac{\Delta k_{x,y} }{k_0}\quad \text{or}\quad \delta_{x,y}\sim O\left(\frac{\partial_{x,y}A}{\epsilon}\right), \end{equation}

where $\delta _x$ ($\Delta k_x$) and $\delta _y$ ($\Delta k_y$) denote the dimensionless (dimensional) bandwidth in the $x$ and $y$ direction, respectively. The bandwidth is not assumed small, in contrast to conventional, reduced-form evolution equations for $A$ by earlier works, e.g. Dysthe (Reference Dysthe1979) and Stiassnie (Reference Stiassnie1984). Specifically, $\delta _{x,y}>1$ are permitted in this paper. Similar to (2.5a), $\zeta$ and $\varPhi$ are expressed in an envelope-type form through the separation of wave harmonics, as follows:

(2.10a)\begin{align} \varPhi(\boldsymbol{x},z,t) &=\textstyle\frac{1}{2}\epsilon (B + \epsilon^2 B^{(31)})\,\mathrm{e}^{k_0z}\,\mathrm{e}^{\mathrm{i} (\boldsymbol{k}_0\cdot\boldsymbol{x}-\omega_0t)} +\text{c.c.} \nonumber\\ &\quad + \epsilon^2 [\varPhi^{(20)}(\boldsymbol{x},z,t)+ (\varPhi^{(22)}(\boldsymbol{x},z,t)+\text{c.c.})] \nonumber\\ &\quad + \epsilon^3(\varPhi^{(33)}(\boldsymbol{x},z,t)+\text{c.c.}), \end{align}

(2.10b)\begin{align} \zeta(\boldsymbol{x},t) &= \textstyle\frac{1}{2} \epsilon(A + \epsilon^2 A^{(31)})\,\mathrm{e}^{\mathrm{i} (\boldsymbol{k}_0\cdot\boldsymbol{x}-\omega_0t)} +\text{c.c.}\nonumber\\ &\quad + \epsilon^2 [\zeta^{(20)} + (\textstyle\frac{1}{2} A^{(22)}\,\mathrm{e}^{2\mathrm{i} (\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t)}+\text{c.c.})] \nonumber\\ &\quad + \textstyle\frac{1}{2} (\epsilon^3 A^{(33)}\,\mathrm{e}^{3\mathrm{i}(\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t)} + \text{c.c.}), \end{align}

in which the superscripts $(ij)$ denote $O(\epsilon ^i)$ and $j$th harmonic, $A^{(ij)}= A^{(ij)}(\boldsymbol {x},t)$ and $B^{(ij)}=B^{(ij)} (\boldsymbol {x},z,t)$ are the modulated amplitude and potential, respectively, and the potentials at second and third harmonics are given by, respectively,

(2.11a)\begin{equation} \varPhi^{(22)}(\boldsymbol{x},z,t) = \textstyle\frac{1}{2} B^{(22)}\,\mathrm{e}^{2k_0z}\,\mathrm{e}^{2\mathrm{i} (\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t)} \end{equation}

and

(2.11b)\begin{equation} \varPhi^{(33)}(\boldsymbol{x},z,t)= \textstyle\frac{1}{2} B^{(33)}\,\mathrm{e}^{3k_0z}\,\mathrm{e}^{3\mathrm{i}(\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t)}. \end{equation}

It is worth noting, that although (2.10a,b) are in a form similar to the so-called harmonic expansion, the wave fields at second and third order in such a form are introduced for convenience. They are due to the separation of wave harmonics arising mainly from the forcing equations (e.g. (3.1b) and (4.1)) at a still water surface, presented in § 3 at second order and in § 4 at third order.

2.3. Linear wave fields

Solutions to a linearized boundary value problem from (2.1)–(2.3) are obtained for linear wave fields. They are given here without detailed derivations (see § 13 in Mei, Stiassnie & Yue Reference Mei, Stiassnie and Yue2005), expressed in terms of $\hat {A}(\boldsymbol {k})$, as follows:

(2.12a)\begin{equation} \boldsymbol{V}^{(1)} \equiv [\boldsymbol{u}^{(1)},w^{(1)}] = \textstyle\frac{1}{2}\bar{\boldsymbol{V}}(\boldsymbol{x},z,t) \,\mathrm{e}^{k_0z}\,\mathrm{e}^{\mathrm{i} (\boldsymbol{k}_0\cdot\boldsymbol{x}-\omega_0t)} +\text{c.c.}, \end{equation}

with

(2.12b)\begin{equation} \bar{\boldsymbol{V}} = [\bar{\boldsymbol{u}}(\boldsymbol{x},z,t), \bar{w}(\boldsymbol{x},z,t) ]\end{equation}

and

(2.12c) \begin{equation} \left[\begin{array}{@{}c@{}} B \\ \bar{\boldsymbol{u}}\\ \bar{w} \end{array}\right]={-}\mathrm{i}\int_{-\infty}^{\infty} \left[\begin{array}{@{}c@{}} 1 \\ \mathrm{i} (\boldsymbol{k}+\boldsymbol{k}_0)\\ |\boldsymbol{k}+\boldsymbol{k}_0| \end{array}\right] c_p(\boldsymbol{k}+\boldsymbol{k}_0) \hat{A}(\boldsymbol{k},t)\,\mathrm{e}^{(|\boldsymbol{k}+\boldsymbol{k}_0|-k_0)z} \,\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d} \boldsymbol{k}^2, \end{equation}

in which $\boldsymbol {V}$ denotes linear velocity vector and $\bar {\boldsymbol {V}}$ its magnitude in the envelope-type form; $\boldsymbol {u}^{(1)}(\boldsymbol {x},z,t)$ and $w^{(1)}(\boldsymbol {x},z,t)$ denote linear velocity vector in the horizontal plane and vertical component, respectively; $\bar {\boldsymbol {u}}$ and $\bar {w}$ are their corresponding magnitude in the envelope-type form; and $c_p(\boldsymbol {k}) = {\sqrt {g /|\boldsymbol {k}|}}$ denotes the phase velocity of wave $\boldsymbol {k}$. It is worth noting that (2.12c) is fundamental to the theory presented in this paper due to the following two aspects: it allows for a nonlinear evolution equation for $A$ being expressed in terms of $\bar {\boldsymbol {V}}$; and it also facilitates simple numerical implementations of the evolution equation for $A$ presented in § 4.

Using the definition of linear velocity potential $\varPhi ^{(1)}$ and following similar procedures that lead to the first-order evolution equation (2.8a), it is understood that the following identities hold (to the first order in wave steepness):

(2.13a,b) \begin{equation} \bar{\boldsymbol{V}} = \nabla_3 B +\boldsymbol{k}^{(3d)}_0 B\quad \text{and}\quad \partial_t \left[\begin{array}{@{}c@{}} B \\ \bar{\boldsymbol{u}}\\ \bar{w} \end{array}\right]+ \mathrm{i}\omega_0\mathcal{L} \left[\begin{array}{@{}c@{}} B \\ \bar{\boldsymbol{u}}\\ \bar{w} \end{array}\right] = 0, \end{equation}

in which $\boldsymbol {k}^{(3\text {d})}_0 = [\mathrm {i} \boldsymbol {k}_0, k_0]$ is introduced for convenience.

2.4. Quadratic property of linear waves: velocity head

Before presenting the wave solutions at second and third order in $\epsilon$, we introduce a leading-order approximation to the wave-velocity head, $H_v(\boldsymbol {x},z,t)$, defined as

(2.14)\begin{equation} H_v(\boldsymbol{x},z,t) = \frac{1}{2g}|\nabla_3\varPhi^{(1)}|^2, \end{equation}

which is correct to $O(\epsilon ^2)$. The wave-velocity head plays an essential role in the forcing of bound waves at second order and, thus, the nonlinear evolution equation for $A$ presented in § 4. Inserting the envelope-type expression for $\varPhi ^{(1)}$ into (2.14) and using (2.12b), we obtain

(2.15a)\begin{equation} H_v = (H^{(22)} +\text{c.c.})+ H^{(20)}, \end{equation}

with

(2.15b)\begin{equation} H_v^{(22)}=\frac{1}{8g} \bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}} \,\mathrm{e}^{2k_0z}\,\mathrm{e}^{2\mathrm{i} (\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t)}\quad \text{and}\quad H_v^{(20)} = \frac{1}{4g} \bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}^*\,\mathrm{e}^{2k_0z}, \end{equation}

in which the asterisk ‘*’ denotes the complex conjugates, and we note that $\bar {\boldsymbol {V}}\boldsymbol {\cdot }\bar {\boldsymbol {V}} = \bar {u}^2+ \bar {v}^2+\bar {w}^2$ with $\bar {u}$ and $\bar {v}$ denoting the velocity component of $\bar {\boldsymbol {u}}$ ($\bar {\boldsymbol {u}}= [\bar {u},\bar {v} ]$) in the $x$ and $y$ direction, respectively. It is noticeable that $H_v^{(22)}$ varies with the superharmonic factor $\exp (2\mathrm {i}\boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-2\mathrm {i}\omega _0t )$ whereas the mean velocity head $H_v^{(20)}$ is independent of the carrier wave phase $\boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\omega _0t$. Especially for a quasi-monochromatic wave group (or a monochromatic wave) that admits $\boldsymbol {\nabla } A \ll O(1)$ (or $\boldsymbol {\nabla } A = 0$), we obtain a leading-order approximation to the mean velocity head $H_v^{(20)}$ given by

(2.16)\begin{equation} H_v^{(20)}(\boldsymbol{x},z,t) = \frac{1}{g}U_{SD}c_{g0}\quad \text{with} \ U_{SD} = k_0\omega_0|A|^2\,\mathrm{e}^{2kz} \text{ and } c_{g0} = \frac{\omega_0}{2k_0}, \end{equation}

where $U_{SD}$ and $c_{g0}$ denote Stokes drift and the group velocity, respectively. We understand from (2.16) that the mean velocity head $H_v^{(20)}$ maintains the energy for the propagation of the Stokes drift at group velocity $c_{g0}$ for a quasi-monochromatic wave group. We show in § 3.3.2 that it is also responsible for the generation of the Eulerian return flow at second order.

In contrast, the rapidly varying velocity head, associated with the factor of second harmonic $\exp (2\mathrm {i}\boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-2\mathrm {i}\omega _0t)$, has been rarely investigated as it is zero for unidirectional waves. In particular, we note

(2.17)\begin{equation} \textstyle\frac{1}{4} \bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}\,\mathrm{e}^{2k_0z}\,\mathrm{e}^{2\mathrm{i} (\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t)} +\text{c.c.} = 0, \end{equation}

which holds for unidirectional waves of a broad bandwidth, as can also be inferred from earlier works, e.g. equations (23) and (25) in Dalzell (Reference Dalzell1999). Due to this, we examine the effects of the rapidly varying velocity head in §§ 3.3.1 and 6.3.1 for multidirectional waves.

3.1. Governing equation and boundary conditions

The velocity potential at second order in wave steepness, $\varPhi ^{(2)} = \varPhi ^{(22)} +\text {c.c.} +\varPhi ^{(20)}$, are described by the following boundary value problem (Dalzell Reference Dalzell1999; Li et al. Reference Li, Zheng, Lin, Adcock and van den Bremer2021c):

(3.1a)\begin{gather} \nabla_3 \varPhi^{(2)} = 0 \quad \text{for}\ -\infty< z<0, \end{gather}

(3.1b)\begin{gather}(\partial_{tt} + g\partial_z)\varPhi^{(2)} = \underbrace{- \zeta^{(1)}\varGamma_z\varPhi^{(1)} }_{{=}0}- \partial_t(|\nabla_3\varPhi^{(1)}|^2) \quad \text{for}\ z= 0, \end{gather}

(3.1c)\begin{gather}\partial_z\varPhi^{(2)} = 0\quad \text{for}\ z\to-\infty, \end{gather}

in which $\varGamma _z$ is defined as $\varGamma _z= \partial _z \varGamma$. As highlighted in (3.1b), the first term on the right-hand side does not contribute to the forcing on a still water surface as it vanishes and, thus, only the second term needs to be considered. Using the second term and (2.15b) we define a forcing term $S$ for $z=0$ as follows:

(3.2a)\begin{equation} S(\boldsymbol{x},t) = [- \partial_t(|\nabla_3\varPhi^{(1)}|^2)]_{z=0}, \end{equation}

which leads to

(3.2b)\begin{equation} S = (S^{(22)} \,\mathrm{e}^{2\mathrm{i} (\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t)}+ \text{c.c.} ) + S^{(20)}, \end{equation}

with

(3.3a)\begin{equation} S^{(22)} ={-} 2g \partial_t H^{(22)}_v(\boldsymbol{x},z,t) \,\mathrm{e}^{{-}2\mathrm{i} (\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t)} \end{equation}

and

(3.3b)\begin{equation} S^{(20)} ={-}2g \partial_t H_v^{(20)}(\boldsymbol{x},z,t) \quad \text{for}\ z=0. \end{equation}

The forcing term $S$ (or the velocity head on a still water surface) is responsible for the forcing of the second-order bound waves. The negative sign on the right-hand side of (3.3a,b) implies a flow generated below a still water surface, propagating in the opposite direction to the change of the velocity head in time at $z=0$ – a so-called return flow which is known for second-order subharmonic waves.

In order to obtain the solution to boundary value problem (3.1a)–(3.1c) for $\varPhi ^{(2)}$, two different and novel approaches are proposed in this paper and presented in §§ 3.2 and 3.3, respectively. In particular, the first approach is applicable to three-dimensional waves of a broad bandwidth, and the second is an approximate method based on the assumption of a narrow bandwidth with $\delta _{x,y}\ll 1$.

3.2. Approach I: semianalytical approach for $\varPhi ^{(2)}$

Approach I, referred to as a semianalytical approach, proposes solving boundary value problem (3.1a)–(3.1c) in the Fourier $\boldsymbol {k}$ plane by using a pseudospectral and a finite difference method. The boundary value problem for $\varPhi ^{(2j)}$ from (3.1a)–(3.1c) is given by

(3.4a)\begin{gather} \nabla^2{\varPhi}^{(2j)} =0\quad \text{for}\ -\infty< z<0, \end{gather}

(3.4b)\begin{gather}\varGamma{\varPhi}^{(2j)} = {S}^{(2j)}\quad \text{for}\ z=0 \end{gather}

and

(3.4c)\begin{equation} \partial_z{\varPhi}^{(2j)} = 0\quad \text{for}\ z\to-\infty, \end{equation}

in which $j=0$ and $j=2$ denote the subharmonic and superharmonic wave, respectively. The solutions of (3.4a,b,c) can be expressed in a form, as follows:

(3.5)\begin{equation} \varPhi^{(2j)} (\boldsymbol{x},z,t) = \mathrm{e}^{j\mathrm{i} (\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t)} \int_{-\infty}^{\infty} 2^{-j/2} \hat{B}^{(2j)} (\boldsymbol{k},t)\exp \mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}+|\boldsymbol{k} +j\boldsymbol{k}_0|z} \,\mathrm{d}\boldsymbol{k}^2, \end{equation}

for $j=0$ and $j=2$. Substituting (3.5) into (3.4b) leads to the second-order differential equation for $\hat {B}^{(2j)}$, as follows:

(3.6)\begin{equation} (\partial_t - \mathrm{i} j \omega_0)^2\hat{B}^{(2j)} + g|\boldsymbol{k}+j\boldsymbol{k}_0| \hat{B}^{(2j)} = 2^{j/2}\hat{S}^{(2j)}(\boldsymbol{k},t), \end{equation}

where $\hat {S}^{(2j)}$ denotes the Fourier transform for ${S}^{(2j)}$ for $j=0$ and $j=2$. It is understood that the second-order waves are bound, and which do not admit a homogeneous solution of (3.6) (Phillips Reference Phillips1960; Hasselmann Reference Hasselmann1962). As a result, (3.6) for $\hat {B}^{(2j)}$ can be solved numerically by a semianalytical approach that proposes evaluating $S^{(2j)}(\boldsymbol {x},t)$ by using a pseudospectral method and obtaining the numerical solution of (3.6) by available time-marching methods. The numerical procedures are explained in § 4.4. For numerical implementations, prescribed conditions at an initial instant are required for $\varPhi ^{(2)}(\boldsymbol {x},0,t_0)$ and $\partial _t\varPhi ^{(2)}(\boldsymbol {x},0,t_0)$. In practice, multiple choices are available to this end, depending on the purpose of wave predictions. For instance, a second-order wavemaker theory based on Schäffer (Reference Schäffer1996), periodic boundary conditions as in Dommermuth & Yue (Reference Dommermuth and Yue1987), the framework by Bonnefoy, Le Touzé & Ferrant (Reference Bonnefoy, Le Touzé and Ferrant2006) for the waves generation in a numerical wave tank, stationary waves based on Dalzell (Reference Dalzell1999) and approximate initial conditions as presented in § 3.3. If the waves generation is based on linear wave theory, $\varPhi ^{(2)}(\boldsymbol {x},0,t_0)=0$ and $\partial _t\varPhi ^{(2)}(\boldsymbol {x},0,t_0)=0$ are implied, leading to inconsistency and additional generation of spurious waves (Schäffer Reference Schäffer1996). The semianalytical approach leads to improved computational efficiency compared with the analytical method by Dalzell (Reference Dalzell1999) based on Fourier integrals. For $N$ Fourier modes, the semianalytical approach requires computational operations at $O(N\,\text {In}(N))$, whereas Dalzell (Reference Dalzell1999) requires $O(N^2)$. The validation of the semianalytical approach is presented in § 6.1.

3.3. Approach II: an approximate method for $\varPhi ^{(2)}$

With a narrow-banded assumption, this section presents an approximate method for $\varPhi ^{(22)}$ (§ 3.3.1) and $\varPhi ^{(20)}$ (§ 3.3.2). For identifying the correct order in bandwidth, we assume small and non-dimensional bandwidth scaling parameters $\delta _x$ and $\delta _y$ that denote the bandwidth in the $x$ and $y$ direction, respectively. The order of accuracy of the approximate solutions is indicated by the product of $\epsilon ^2$ and $\delta _{x,y}^j$ with $j$ non-zero integers.

With a narrow-banded assumption, we seek the solutions of (3.1a)–(3.1c) for unknown $B^{(22)}$ and $\varPhi ^{(20)}$ in a form of power series in bandwidth, as follows:

(3.7a,b)\begin{equation} B^{(22)} =\delta_yB^{(22)}_{{\approx},1} + \delta^2_{x,y} B^{(22)}_{{\approx},2} +\cdots\quad \text{and}\quad \varPhi^{(20)} = \delta_x \varPhi^{(20)}_{{\approx},1} + \delta_{x,y}^2 \varPhi^{(20)}_{{\approx},2} +\cdots, \end{equation}

in which the subscripts ‘$\approx , j$’ denote an approximation at $O(\epsilon ^2\delta _{x,y}^j)$, and we consider leading-order approximations up to $j=2$.

3.3.1. Forcing of second-order superharmonic waves by multidirectional waves

It is understood that the forcing term on the right-hand side of (3.1b) (see also (2.17)) for $\varPhi ^{(22)}$ for unidirectional waves is zero, which suggests $\varPhi ^{(22)}(\boldsymbol {x},z,t)=0$. Hence, the derivations for $\varPhi ^{(22)}$ are only needed for multidirectional waves. Based on the boundary value problem described by (3.1a)–(3.1c), the forcing equation of the superharmonic bound waves for $B^{(22)}$ on a still water surface reads

(3.8)\begin{equation} (\partial_{tt}-4\mathrm{i}\omega_0\partial_{t} -4\omega_0^2+ g\partial_z +2gk_0)B^{(22)} = \mathrm{i} \omega_0\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}} -\bar{\boldsymbol{V}}\boldsymbol{\cdot}\partial_t \bar{\boldsymbol{V}} \quad \text{for}\ z=0, \end{equation}

in which the first term on the right-hand side is of $O(\epsilon ^2\delta _y)$ and the second term of $O(\epsilon ^2\delta _y\delta _x)$. Inserting the expanded form (3.7a) for $B^{(22)}$ into (3.8) and the Laplace equation for $\varPhi ^{(22)}$ and following the conventional procedures for perturbed solutions, we obtain

(3.9a)\begin{equation} B^{(22)}_{{\approx},1} = \int_{-\infty}^{\infty} \frac{\mathrm{i}\omega_0\hat{V}_{sq}^{(22)}(\boldsymbol{k},t)}{g|\boldsymbol{k}+2\boldsymbol{k}_0|-4\omega_0^2} \,\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x} + (|\boldsymbol{k}+2\boldsymbol{k}_0| -2k_0)z}\,\mathrm{d} \boldsymbol{k}^2 \end{equation}

and

(3.9b)\begin{align} B^{(22)}_{{\approx},2} = \int_{-\infty}^{\infty} \frac{ -4\omega^2_0\partial_t\hat{V}_{sq}^{(22)}+( 4\omega_0^2- g|\boldsymbol{k}+2\boldsymbol{k}_0|) \hat{V}^{(22)}_{sq,t}(\boldsymbol{k},t)}{(g|\boldsymbol{k}+2\boldsymbol{k}_0|-4\omega_0^2)^2} \,\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x} + (|\boldsymbol{k}+2\boldsymbol{k}_0| -2k_0)z}\,\mathrm{d} \boldsymbol{k}^2, \end{align}

where $\hat {V}_{sq}^{(22)}$ and $\hat {V}_{sq,t}^{(22)}$ denote the Fourier transform for $\bar {\boldsymbol {V}}\boldsymbol {\cdot }\bar {\boldsymbol {V}}$ and $\bar {\boldsymbol {V}}\boldsymbol {\cdot }\partial _t\bar {\boldsymbol {V}}$ with respect to $\boldsymbol {x}$ for $z=0$, respectively.

3.3.2. Potential for the second-order mean flows in two and three dimensions

Similar to the second-order superharmonic waves, the forcing equation of the subharmonic bound waves at second order for $\varPhi ^{(20)}$ on a still water surface (cf. (3.1b)) reads

(3.10)\begin{equation} \partial_{tt}\varPhi^{(20)} + g\partial_z \varPhi^{(20)} = S^{(20)}\quad \text{for}\ z =0, \end{equation}

where $S^{(20)}$ is defined by (3.3b) and $S^{(20)} \sim O(\epsilon ^2\delta _x)$. With a narrow-banded assumption and following the same procedures as for $B^{(22)}$, we obtain

(3.11a,b)\begin{align} \varPhi^{(20)}_{{\approx},1} = \int_{-\infty}^{\infty} \frac{\hat{S}^{(20)}}{g|\boldsymbol{k}|} \,\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x} + |\boldsymbol{k}|z} \,\mathrm{d} \boldsymbol{k}^2\quad \text{and}\quad \varPhi^{(20)}_{{\approx},2} = \int_{-\infty}^{\infty} \frac{-2 \hat{S}_t^{(20)}}{(g|\boldsymbol{k}|)^2} \,\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x} + |\boldsymbol{k}|z} \,\mathrm{d} \boldsymbol{k}^2, \end{align}

in which $\hat {S}^{(20)}(\boldsymbol {k},t)$ and $\hat {S}_t^{(20)}(\boldsymbol {k},t)$ denote the Fourier transform for $S^{(20)}$ and $\partial _t S^{(20)}$ with respect to $\boldsymbol {x}$, respectively.

3.4. Wave elevation $\zeta ^{(2)}$ and velocity $\boldsymbol {V}^{(2)}$ at second order

With second-order potential $\varPhi ^{(2)}$ given by the semianalytical approach presented in § 3.2 or the approximate method presented in § 3.3, the surface elevation at second order is obtained from

(3.12)\begin{equation} \zeta^{(2)}(\boldsymbol{x},t) ={-}\frac{1}{g} \left[\partial_t\varPhi^{(2)} +\zeta^{(1)} \partial_{tz}\varPhi^{(1)} + \displaystyle\frac{1}{2} |\nabla_3\varPhi^{(1)}|^2\right]\quad \text{for}\ z = 0, \end{equation}

where $\varPhi ^{(2)} = \varPhi ^{(20)}$ for unidirectional waves and $\varPhi ^{(2)} = \varPhi ^{(20)} + \varPhi ^{(22)}+\text {c.c.}$ for multidirectional waves. Inserting the envelope-type expression for $\varPhi ^{(1)}$ and $\varPhi ^{(2)}$, we arrive at

(3.13a)\begin{equation} \zeta^{(2)} = \zeta^{(20)} + (\zeta^{(22)} +\text{c.c.}), \end{equation}

with

(3.13b)\begin{equation} \zeta^{(22)} ={-}\frac{1}{g} \left[\partial_t\varPhi^{(22)} + \frac{1}{4} \left(A(\partial_t-\mathrm{i}\omega_0)\bar{w} +\frac{1}{2}\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}\right)\mathrm{e}^{2\mathrm{i} (\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t)} \right] \end{equation}

and

(3.13c)\begin{equation} \zeta^{(20)} ={-}\frac{1}{g} \left[\partial_t\varPhi^{(20)} + \displaystyle\frac{1}{4} \bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}^*+\frac{1}{4}(A^*(\partial_t-\mathrm{i}\omega_0)\bar{w}+\text{c.c.}) \right]\quad \text{for}\ z=0. \end{equation}

With $\varPhi ^{(22)} =0$ in mind for unidirectional waves, (3.13b) gives an exact expression for the elevation for second-order superharmonic waves that depend only on linear modulational parameters, e.g. $A$ and $\bar {\boldsymbol {V}}$.

Similarly, the velocity at second order, $\boldsymbol {V}^{(2)}(\boldsymbol {x},z,t)$, is given by

(3.14a)\begin{equation} \boldsymbol{V}^{(2)} = \boldsymbol{V}^{(20)}(\boldsymbol{x},z,t) + [ \bar{\boldsymbol{V}}^{(22)}(\boldsymbol{x},z,t)\,\mathrm{e}^{2k_0z}\,\mathrm{e}^{2\mathrm{i} (\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t)} +\text{c.c.}], \end{equation}

where

(3.14b)\begin{equation} \boldsymbol{V}^{(20)} = \boldsymbol{\nabla} \varPhi^{(20)} \quad \text{and}\quad \bar{\boldsymbol{V}}^{(22)} =\textstyle\frac{1}{2} (2\boldsymbol{k}^{(3\text{d})}_0 +\nabla_3 )B^{(22)}. \end{equation}

4.1. An evolution equation correct to $O(\epsilon ^3)$ for waves of a broad bandwidth

Collecting the terms at third order in wave steepness in (2.2), we obtain

(4.1)\begin{equation} (\partial_{tt} + g\partial_z)\varPhi^{(3)} ={-} Q^{(3)}(\boldsymbol{x},z,t)\quad \text{for}\ z=0, \end{equation}

where the forcing term $Q^{(3)}$ is given by

(4.2)\begin{align} Q^{(3)}& = \underbrace{ \zeta^{(2)} \varGamma_z \varPhi^{(1)} + \textstyle\frac{1}{2}( \zeta^{(1)})^2 \varGamma_{zz}\varPhi^{(1)} }_{= 0} + \zeta^{(1)}\partial_z [ \varGamma \varPhi^{(2)} + \partial_{t}(|\nabla_3\varPhi^{(1)}|^2) ] \nonumber\\ &\quad + 2\partial_t (\nabla_3\varPhi^{(1)}\boldsymbol{\cdot} \nabla_3\varPhi^{(2)}) + \textstyle\frac{1}{2} \nabla_3\varPhi^{(1)} \boldsymbol{\cdot} \nabla_3 (|\nabla_3\varPhi^{(1)}|^2), \end{align}

with $\varGamma _{zz} = \partial _z \varGamma _z$. As indicated in (4.2), one would show that $\varGamma _z \varPhi ^{(1)}=0$ and $\varGamma _{zz}\varPhi ^{(1)}=0$, which are not necessarily so for a finite depth. The terms in the square bracket on the right-hand side of (4.2) can be simplified due to (3.1b). In particular, it is noticed that the following identity holds:

(4.3)\begin{equation} \partial_z [ \varGamma \varPhi^{(2)} + \partial_{t}(|\nabla_3\varPhi^{(1)}|^2)]= k_0 \partial_t(\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}^*) \quad \text{for}\ z=0. \end{equation}

Inserting the harmonic expansion for $\varPhi ^{(1)}$ and $\varPhi ^{(2)}$ into (4.2) leads to $Q^{(3)}$ being expressed in a form as follows:

(4.4)\begin{equation} Q^{(3)} (\boldsymbol{x},z,t) = Q^{(31)}(\boldsymbol{x},z,t) + Q^{(33)}(\boldsymbol{x},z,t), \end{equation}

where $Q^{(31)}$ denotes the forcing term with a factor $\exp {(\mathrm {i}\boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\mathrm {i}\omega _0t)}$ and $Q^{(33)}$ the term with a factor $\exp {(3\mathrm {i}\boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-3\mathrm {i}\omega _0t)}$. Leaving the expression for $Q^{(33)}$ in § 5 and focusing on $Q^{(31)}$ in this section, we obtain (see details in Appendix A)

(4.5a)\begin{equation} Q^{(31)} = \bar{Q}^{(31)}(\boldsymbol{x},z,t) \,\mathrm{e}^{\mathrm{i} (\boldsymbol{k}_0\cdot\boldsymbol{x}-\omega_0t)} +\text{c.c.}\quad \text{for}\ z=0, \end{equation}

with $\bar {Q}^{(31)} =\bar {Q}_{all}^{(31)} + \bar {Q}_{dir}^{(31)}$,

(4.5b)\begin{align} \bar{Q}_{all} ^{(31)} &= (-\mathrm{i}\omega_0 + \partial_t) [\nabla_3\varPhi^{(20)} \boldsymbol{\cdot}\bar{\boldsymbol{V}}]+ \textstyle\frac{1}{8} [\nabla_3(\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}^*)]\boldsymbol{\cdot} \bar{\boldsymbol{V}}\nonumber\\ &\quad +\textstyle\frac{1}{4} (k_0\bar{w} +2 k_0 A\partial_t) (\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}^*), \end{align}

(4.5c)\begin{align} \bar{Q}_{dir} ^{(31)} &= (-\mathrm{i}\omega_0 + \partial_t) \{ \textstyle\frac{1}{2} [(2\boldsymbol{k}_0^{(3d)} + \nabla_3) B^{(22)}]\boldsymbol{\cdot}\bar{\boldsymbol{V}}^* \}\nonumber\\ &\quad + \textstyle\frac{1}{16}[(2\boldsymbol{k}_0^{(3d)} + \nabla_3) (|\bar{\boldsymbol{V}}|^2)]\boldsymbol{\cdot} \bar{\boldsymbol{V}}^*, \end{align}

in which the subscripts ‘all’ and ‘dir’ are used to distinguish two different cases; the former means the parameter expressed in terms that are non-constant for both unidirectional and multi-unidirectional waves, and ‘dir’ means the parameter composed of terms that are non-zero only for directional waves. Hence, $\bar {Q} ^{(31)}$ admits a much simpler expression for waves considered in two dimensions or for long-crested waves due to $\bar {Q}^{(31)}_{dir}=0$.

Physically, $Q^{(31)}$ would lead to secular solutions that should be removed (see the discussion on page 376 in Madsen & Fuhrman (Reference Madsen and Fuhrman2006)). In order to achieve this, a conventional approach is to introduce a nonlinear amplitude-dependent frequency $\omega _2$ (${\sim }O(\epsilon ^2)$) for the first-order wave fields through replacing factor $\exp [\mathrm {i} (\boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x} -\omega _0 t )]$ in (2.5a) and (2.12) with $\exp [\mathrm {i} (\boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\omega _0 t ) - \mathrm {i}\epsilon ^2\omega _2t]$. Thereby, it leads to an additional contribution of the updated first-order potential $\varPhi ^{(1)}$ to (4.1) at $O(\epsilon ^3)$. Especially arising from $\partial _{tt}\varPhi ^{(1)}$, we obtain an additional term associated with first harmonic at $O(\epsilon ^3)$ expressed as $-\mathrm {i}\omega _2(\partial _t-\mathrm {i}\omega _0)B\exp (\mathrm {i} (\boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\omega _0 t ) - \mathrm {i}\epsilon ^2\omega _2t)+\text {c.c.}$ Mathematically, the secular solutions can be removed if the sum of this additional term and $\bar {Q} ^{(31)}$ equal zero. Hence, we arrive at

(4.6)\begin{equation} \mathrm{i}\omega_2(\partial_t-\mathrm{i}\omega_0)B = \bar{Q}^{(31)}\quad \text{for}\ z=0. \end{equation}

Noticing that $\zeta ^{(1)} = -\partial _t\varPhi ^{(1)}/g$ for $z=0$ which gives $(\partial _t -\mathrm {i}\omega _0 )B = -gA$, substituting this expression for $B$ into (4.6) leads to the equation for $A$, as follows:

(4.7)\begin{equation} {-}\mathrm{i} g\omega_2 A = \bar{Q}^{(31)}\quad \text{for}\ z=0. \end{equation}

The evolution equation for $A'$ with $A'=A\exp (-\mathrm {i}\omega _2t)$, correct to third order in wave steepness, can now be expressed as

(4.8)\begin{equation} \partial_t A' + \mathrm{i} \omega_0\mathcal{L} A' + \mathrm{i}\omega_2A' =0. \end{equation}

Inserting (4.7) into (4.8) and omitting the prime leads to

(4.9)\begin{equation} \partial_t A + \mathrm{i} \omega_0 \mathcal{L} A + \mathcal{N}(\boldsymbol{x},z,t)=0\quad \text{for}\ z=0, \end{equation}

where $\mathcal {N}=-\bar {Q}^{(31)} /g$ denotes the nonlinear term, given by

(4.10a)\begin{equation} \mathcal{N} = \mathcal{N}_{all} + \mathcal{N}_{dir}, \end{equation}

with

(4.10b)\begin{gather} \mathcal{N}_{all} = \frac{1}{g} (\mathrm{i}\omega_0 - \partial_t) [\boldsymbol{V}^{(20)} \boldsymbol{\cdot}\bar{\boldsymbol{V}}]- \frac{1}{8g} [2k_0\bar{w} + \bar{\boldsymbol{V}} \boldsymbol{\cdot} \nabla_3 +4k_0 A\partial_t](\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}^*), \end{gather}

(4.10c)\begin{gather}\mathcal{N}_{dir} = \frac{1}{g} (\mathrm{i}\omega_0 - \partial_t) (\bar{\boldsymbol{V}}^{(22)}\boldsymbol{\cdot}\bar{\boldsymbol{V}}^* ) - \frac{1}{16g}[(2\boldsymbol{k}_0^{(3d)} + \nabla_3) (|\bar{\boldsymbol{V}}|^2)]\boldsymbol{\cdot} \bar{\boldsymbol{V}}^*. \end{gather}

Equation (4.9) denotes the new NLSE that describes the evolution of linear wave envelope $A$ in time and space. In addition to the linear term associated with the time derivative $\partial _t(\ldots )$, (4.9) is composed of a term that denotes linear dispersion relation due to $\mathrm {i}\omega _0 \mathcal {L}$ and a nonlinear term denoted by $\mathcal {N}$ that is correct to $O(\epsilon ^3)$. For (4.9), no assumptions associated with the bandwidth have been introduced. Thereby, it has relaxed the limitation arising from a narrow-banded assumption on which most NLSEs are based, such as the NLSEs by Dysthe (Reference Dysthe1979), Trulsen & Dysthe (Reference Trulsen and Dysthe1996) and Trulsen et al. (Reference Trulsen, Kliakhandler, Dysthe and Velarde2000); i.e. the new NLSE should permit the prediction of waves with bandwidth $\delta _{x,y}>1$. The implementation of the new NLSE for the wave envelope, $A$, of the first-harmonic elevation requires the evaluation of linear term $\partial _tA$ in the Fourier plane based on (2.8c) and the evaluation of the linear envelope (vector) of both the linear and second-order velocity for $\mathcal {N}$, which will be explained in detail in § 4.4.

Extensions of the new NLSE (i.e. (4.9)) to more general cases are simple and straightforward. For example, the new NLSE can be extended to consider a finite water depth following the derivations in this paper, except that one would expect that the nonlinear term $\mathcal {N}$ may contain a few more terms that introduce a negligible additional computational cost. Following Dias et al. (Reference Dias, Dyachenko and Zakharov2008) and Kharif et al. (Reference Kharif, Kraenkel, Manna and Thomas2010), the new NLSE would allow for weakly damped/forced free surface flows through adding the viscous effects in the dynamic and kinematic boundary condition (cf. equations (35, 37) by Dias et al. (Reference Dias, Dyachenko and Zakharov2008)) at the water surface. Moreover, despite that the new NLSE is formulated in the framework of potential flow theory for simplicity, it is essentially an equation associated with linear envelope $A$, linear velocity $\bar {\boldsymbol {V}}$ and second-order velocity $\bar {\boldsymbol {V}}^{(2)}$. This suggests the weak effects due to viscosity and vorticity in a subsurface layer can be easily incorporated, allowing for waves interacting with ambient environments, e.g. turbulence and background rotational flows (cf. McWilliams, Restrepo & Lane Reference McWilliams, Restrepo and Lane2004). These aspects will be considered in future work.

4.2. Leading-order approximations to the new NLSE

Approximations to (4.9) correct to $O(\epsilon ^3\delta _{x,y}^j)$, with $j=0,1, 2, 3$, are presented in this section with a narrow-banded assumption, i.e. $O(\delta )\ll 1$. In particular, we assume that the bandwidths both in the longitudinal ($x$) direction and the transverse direction to the propagation of the carrier wave are small and that $O(\delta _y)\sim O(\delta ^2_x)$. Small directional spread is, hence, implied. Based on previous papers (e.g. Chu & Mei Reference Chu and Mei1971; Li et al. Reference Li, Zheng, Lin, Adcock and van den Bremer2021c), the order of magnitude associated with the relevant terms in (4.5) is assumed, as follows:

(4.11a–c)\begin{equation} [\bar{\boldsymbol{V}},\bar{\boldsymbol{V}}^*]\sim O(\epsilon),\quad [|\bar{\boldsymbol{V}}|^2, \bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}^*]\sim O(\epsilon^2\delta^2_y,\epsilon^2),\quad \varPhi^{(20)}\sim O(\epsilon^2\delta_x) \end{equation}

and

(4.11d)\begin{equation} B^{(22)} \sim O(\epsilon^2\delta^2_y). \end{equation}

Moreover, it is understood that the derivative operators $(\partial _t, \nabla _3 )$ on the terms in (4.11) would lead to an order of magnitude higher in bandwidth $\delta _{x,y}$. With the right orders of magnitude in mind, collecting the terms up to a particular order in $\epsilon$ and $\delta$ based on (4.5) and (4.9b), and, thus, leading-order approximations to (4.9) are obtained. In particular, they are expressed in a form, as follows:

(4.12)\begin{equation} \partial_{t} A + \mathrm{i}\omega_0\mathcal{L} A + \mathcal{N}^{(3,j)}_{{\approx}} = O(\epsilon^3\delta_{x,y}^{j+1}) \quad \text{for}\ z=0, \end{equation}

where $\mathcal {N}^{(3,j)}_{\approx }$ denotes an approximation to $\mathcal {N}$ correct to $O(\epsilon ^3\delta _{x,y}^j)$, with $j=0, 1, 2$ and $3$, given by

(4.13a)\begin{gather} \mathcal{N}^{(3,0)}_{{\approx}} ={-}\frac{1}{4g}(\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}^*)k_0\bar{w}, \end{gather}

(4.13b)\begin{gather}\mathcal{N}^{(3,1)}_{{\approx}} = \mathcal{N}^{(3,0)}_{{\approx}} -\frac{1}{8g} [\nabla_3(\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}^*)]\boldsymbol{\cdot} \bar{\boldsymbol{V}} - \frac{k_0}{2g} A\partial_t(\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}^*), \end{gather}

(4.13c)\begin{gather}\mathcal{N}^{(3,2)}_{{\approx}} = \mathcal{N}^{(3,1)}_{{\approx}} + \frac{\mathrm{i}\omega_0 }{g} (\nabla_3\varPhi_{{\approx},1}^{(20)} \boldsymbol{\cdot}\bar{\boldsymbol{V}}) + \mathcal{N}^{(3,2)}_{{\approx},{dir}}, \end{gather}

(4.13d)\begin{gather}\mathcal{N}^{(3,3)}_{{\approx}} = \mathcal{N}^{(3,2)}_{{\approx}} + \frac{\mathrm{i}\omega_0 }{g} (\nabla_3\varPhi_{{\approx},2}^{(20)} \boldsymbol{\cdot}\bar{\boldsymbol{V}}) -\frac{1}{g} \partial_t(\nabla_3\varPhi_{{\approx},1}^{(20)} \boldsymbol{\cdot} \bar{\boldsymbol{V}}) + \mathcal{N}^{(3,3)}_{{\approx},{dir}}, \end{gather}

(4.13e)\begin{gather}\mathcal{N}^{(3,4)} = \mathcal{N}, \end{gather}

where (4.13e) denotes that $\mathcal {N}^{(3,4)}$ has an identical expression to $\mathcal {N}$. As noted earlier, the subscript ‘dir’ means the parameter being composed of nonlinear terms that are non-zero only for multidirectional waves, given by

(4.13f)\begin{gather} \mathcal{N}^{(3,2)}_{{\approx},{dir}} = \frac{1}{g} \left( B_{{\approx},1}^{(22)} -\frac{1}{8}\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}} \right) ( \mathrm{i}\boldsymbol{k}_0\boldsymbol{\cdot}\bar{\boldsymbol{u}}^*+ k_0\bar{w}^*), \end{gather}

(4.13g)\begin{gather} \mathcal{N}^{(3,3)}_{{\approx},{dir}} ={-} \frac{1}{16 g}[\nabla_3(|\bar{\boldsymbol{V}}|^2)]\boldsymbol{\cdot} \bar{\boldsymbol{V}}^* - \frac{1}{g} \partial_t [( \mathrm{i}\boldsymbol{k}_0\boldsymbol{\cdot}\bar{\boldsymbol{u}}^*+ k_0\bar{w}^*) B^{(22)}_{{\approx},1}] \nonumber\\ \hspace{-5.2pc}+ \frac{1}{g}B_{{\approx},2}^{(22)} ( \mathrm{i}\boldsymbol{k}_0\boldsymbol{\cdot}\bar{\boldsymbol{u}}^*+ k_0\bar{w}^*). \end{gather}

Equation (4.13b) denotes that a leading-order approximation to $\mathcal {N}$, correct to $O(\epsilon ^3\delta )$, does not require an evaluation of the second-order mean potential. Equation (4.13c) suggests that the approximations to $\varPhi ^{(20)}$ and $B^{(22)}$ that are presented in § 3.3 are sufficient to give a leading-order approximation to (4.9), correct up to $O(\epsilon ^3\delta ^3)$. Any approximations correct to an order higher than $O(\epsilon ^3\delta ^2)$ would require the same computational cost as $\mathcal {N}$ because the approximations $\varPhi ^{(20)}_{\approx ,2}$ and $B^{(22)}_{\approx ,2}$ do not lead to enhanced computational efficiency, compared with solving for $\varPhi ^{(2)}$ by the semianalytical method.

4.3. Comparisons with two reduced-form equations for $A$

Two limiting cases are considered in this section in order to demonstrate that the new NLSE can recover two reduced-form equations for $A$, including the NLSE for the evolution of a train of Stokes wave (§ 4.3.1) and the NLSEs by Dysthe (Reference Dysthe1979), Trulsen & Dysthe (Reference Trulsen and Dysthe1996) and Trulsen et al. (Reference Trulsen, Kliakhandler, Dysthe and Velarde2000) with a narrow-banded assumption (§ 4.3.2).

4.3.1. Uniform Stokes waves

For uniform Stokes waves on an infinite depth, that denote $A$ and $\bar {\boldsymbol {V}}$ having no slow modulation in $\boldsymbol {x}$ and $z$, the new NLSE leads to

(4.14)\begin{equation} \mathcal{N} = \mathcal{N}^{(3,0)}\quad \text{and hence},\ \mathcal{N} = \tfrac{1}{2} \mathrm{i} k^2_0 \omega_0|A|^2A. \end{equation}

Introducing $\epsilon = k_0^2|A|^2$ and $\alpha _1 =\frac {1}{2}$ for Stokes waves, we obtain the evolution equation for $A$ as follows:

(4.15)\begin{equation} \partial_t A + \mathrm{i} \alpha_1\epsilon^2\omega_0 A =0, \end{equation}

which agrees with the classic NLSE for a train of Stokes waves as presented in many papers, e.g. Benjamin & Feir (Reference Benjamin and Feir1967) and Zakharov (Reference Zakharov1968). Equation (4.15) will be used for the analysis of sideband instability of Stokes waves presented in § 6.2.

4.3.2. Waves of a narrow bandwidth ($O(\epsilon ^3\delta _{x,y})$)

In order to recover the previous versions of NLSEs by Dysthe (Reference Dysthe1979), Trulsen & Dysthe (Reference Trulsen and Dysthe1996) and Trulsen et al. (Reference Trulsen, Kliakhandler, Dysthe and Velarde2000) from the NLSE (4.9), two aspects are addressed. First, the term associated with linear dispersion relation $\mathrm {i}\omega _0\mathcal {L}$ is identical to equations (12) and (13) in Trulsen et al. (Reference Trulsen, Kliakhandler, Dysthe and Velarde2000) using the dimensional versions. It has been demonstrated in Trulsen et al. (Reference Trulsen, Kliakhandler, Dysthe and Velarde2000) how this term can recover the terms associated with linear dispersion relation derived by Trulsen & Dysthe (Reference Trulsen and Dysthe1996) and Dysthe (Reference Dysthe1979). Second, the nonlinear term in the NLSEs by Trulsen et al. (Reference Trulsen, Kliakhandler, Dysthe and Velarde2000) and Trulsen & Dysthe (Reference Trulsen and Dysthe1996) are directly obtained from the Dysthe equation by Dysthe (Reference Dysthe1979). Hence, we only have to demonstrate that $\mathcal {N}$ can recover the nonlinear terms in the Dysthe equation. To this end, we start with (4.9b) for unidirectional waves as the Dysthe equation has neglected the contribution from $\varPhi ^{(22)}$ (which is non-zero only for multidirectional waves). Moreover, a narrow-banded assumption as in Dysthe (Reference Dysthe1979) is introduced, i.e. $O(\delta _{x,y})\sim O(\epsilon )$ is assumed in this section.

With a narrow-banded assumption, the velocity (magnitude) $\bar {\boldsymbol {V}}$ correct to $O(\epsilon \delta _{x,y})$ reads (cf. equations (13.2.21) and (13.2.30) by Mei et al. (Reference Mei, Stiassnie and Yue2005) for deep-water waves)

(4.16)\begin{equation} \bar{\boldsymbol{V}}(\boldsymbol{x},z,t) = \frac{-\mathrm{i} g}{\omega_0} \left[\boldsymbol{k}^{(3d)}_0 A + \left(\begin{array}{c} -\mathrm{i} (\boldsymbol{k}_0 z +\mathrm{i} \boldsymbol{e}_x)\boldsymbol{\nabla} A \\ -\mathrm{i}(k_0z +1)\partial_xA \end{array}\right)\right] + O(\epsilon\delta^2), \end{equation}

in which ${\boldsymbol {e}}_x = [1,0]$ denotes a unit vector along the $x$ direction in the horizontal plane. Inserting (4.16) into (4.9b) leads to a leading-order approximation to $g\mathcal {N}$, as follows (see Appendix A.2 for details):

(4.17)\begin{align} \varGamma_{Dysthe}C&\equiv{-}g\mathcal{N}=4k^4_0C|C|^2+ 8\mathrm{i} k_0^3C ( C\partial_xC^*-C^*\partial_xC)-4\mathrm{i} k_0C\partial_x(C^*C) \nonumber\\ &\quad +2k_0\omega_0C(\partial_x \varPhi^{(20)} -\partial_z \varPhi^{(20)} ), \end{align}

where $C$ is defined as $C \equiv B/2= - \mathrm {i}\omega _0 A/(2k_0)$. The subscript ‘Dysthe’ means that the parameter corresponds to the ‘$\varGamma$’ in Dysthe (Reference Dysthe1979); i.e. (4.17) recovers ‘$\varGamma$’ on the right-hand side of equation (2.17) of Dysthe (Reference Dysthe1979), while observing the following two aspects. Firstly, the notation $C$ defined here is denoted by ‘$A$’ in Dysthe (Reference Dysthe1979). Secondly, there is one misprint in the equation for ‘$\varGamma$’, as also pointed out by Brinch-Nielsen & Jonsson (Reference Brinch-Nielsen and Jonsson1986) (cf. p. 461) and Janssen (Reference Janssen1983) – i.e. the second term in the expression ‘$3\mathrm {i} k(\ldots )$’ for ‘$\varGamma$’ should read $8\mathrm {i} k(\ldots )$. For completeness, the detailed derivations for (4.17) are presented in Appendix A.2. The modified NLSE by Trulsen et al. (Reference Trulsen, Kliakhandler, Dysthe and Velarde2000) (i.e. equations (19–22) therein) is essentially a leading-order approximation to (4.9), correct to $O(\epsilon ^3\delta _{x,y})$.

4.4. Numerical implementation

Figure 1 shows the detailed procedures of the numerical implementation of the new NLSE (4.9) for $A$, including the solution of the second-order equation, (3.6), for $B^{(2j)}$ and the evaluation of $\mathcal {N}$ by a pseudospectral method and a finite-difference method. We remark here that it is conventional to solve an NLSE numerically for $A$ in a moving coordinate system through the transforms $\xi =x-c_{g0}t$ and $\eta = t$, which would also follow similar numerical procedures as follows, but for $A'(\xi ,y,\eta )$ with $A'(\xi ,y,\eta ) = A(\boldsymbol {x},t)$ instead.

Figure 1. Flow diagram of the numerical implementation of the new NLSE for $A$, and of the numerical solution of (3.6) for $B^{(2j)}(\boldsymbol {x},z,t)$. A dot denotes the derivative with respect to time.

At any temporal step $t_n =t_0 + n\Delta t$ with $t_0$ denoting the initial time, it is understood that the value of $A(\boldsymbol {x},t_n)$ for all $\boldsymbol {x}_{{s}}<\boldsymbol {x}<\boldsymbol {x}_e$, with $\boldsymbol {x}_s$ and $\boldsymbol {x}_e$ denoting the start and end position vector of a chosen spatial domain, respectively, is given from earlier computations. The spatial domain is discretized by $2N_x\times 2N_y$ points at which $A(\boldsymbol {x},t_n)$ is known, with $2N_x$ and $2N_y$ denoting the total number of points in $x$ and $y$ direction, respectively. Transforming $A$ to the discretized Fourier space using a fast Fourier transform (FFT) (Frigo & Johnson Reference Frigo and Johnson1998), we obtain

(4.18a)\begin{gather} \hat{A}(\boldsymbol{k},t_n) = \mathcal{F}\{A(\boldsymbol{x},t) \},\quad \text{with}\ \boldsymbol{k} = (k_x,k_y), \end{gather}

(4.18b)\begin{gather}k_x = (0, \pm 1,\pm 2,\ldots,\pm N_x)\times \mathrm{d} k_x,\quad k_y = (0, \pm 1,\pm 2,\ldots,\pm N_y)\times \mathrm{d} k_y, \end{gather}

where $\mathrm {d} k_x$ and $\mathrm {d} k_y$ denote the interval between two adjacent points in the Fourier $k_x$ and $k_y$ direction, respectively and $\mathcal {F}$ denotes a FFT with respect to $\boldsymbol {x}$.

As highlighted in grey in figure 1 with the known value of $\hat {A}(\boldsymbol {k},t_n)$, the new NLSE is solved numerically by using a split-step Fourier method (Tappert Reference Tappert1974). The main difference between the numerical procedures presented here and these presented in Lo & Mei (Reference Lo and Mei1985), lie in the fact that the evaluation of $\mathcal {N}$ per time step is obtained differently as $\mathcal {N}$ has different expressions. The numerical procedures for $\mathcal {N}$ at $t=t_n$ are decomposed into three consecutive steps, as explained in the following.

(Step I) The first step is to use the known value of $\hat {A}(\boldsymbol {k},t_n)$ to evaluate $\bar {\boldsymbol {V}}(\boldsymbol {x},0,t_n)$ and its derivative with respect to $t$, $x$, $y$ and $z$ using an inverse FFT, as highlighted in figure 1 in blue. Let $\widehat {\bar {\boldsymbol {V}}}$ be the transformed linear velocity, $\bar {\boldsymbol {V}}$, in the Fourier $\boldsymbol {k}$ space. We understand from (2.12) that

(4.19a)\begin{equation} \widehat{\bar{\boldsymbol{V}}}(\boldsymbol{k},z) =\left[\begin{array}{@{}c@{}} \boldsymbol{k}+\boldsymbol{k}_0\\ -\mathrm{i} |\boldsymbol{k}+\boldsymbol{k}_0| \end{array}\right] c_p(\boldsymbol{k}+\boldsymbol{k}_0)\hat{A} \,\mathrm{e}^{|\boldsymbol{k}_0+\boldsymbol{k}|z-k_0z}, \end{equation}

and thus

(4.19b)\begin{equation} \mathcal{F}\{ [\partial_t, \partial_x, \partial_y,\partial_z] \bar{\boldsymbol{V}} \} =\mathrm{i} [-\omega(\boldsymbol{k}+\boldsymbol{k}_0)-\omega_0,k_x+k_0,k_y, |\boldsymbol{k}_0+\boldsymbol{k}|-k_0 ] \hat{\bar{\boldsymbol{V}}}(\boldsymbol{k},z). \end{equation}

Hence, $\bar {\boldsymbol {V}}$ and also its derivatives $[\partial _t, \partial _x, \partial _y, \partial _z]\bar {\boldsymbol {V}}$ are evaluated for all spatial points through an inverse FFT using (4.19) for $z=0$.

(Step II) The second step (highlighted in green in figure 1) is to solve the second-order equations (3.6) for $B^{(2j)}$ with $j=0$ and $j=2$ by the semianalytical approach presented in §3.2. The forcing terms $S^{(2j)}$ for $j=0$ and $j=2$ are computed based on (3.3) in the physical plane using $\partial _t\bar {\boldsymbol {V}}$, $\bar {\boldsymbol {V}}$ and $\bar {\boldsymbol {V}}^*$ evaluated in step I. The second-order potentials $\hat {{B}}^{(2j)}$ and their derivatives with respect to time, $\partial _t{\hat {B}}^{(2j)}$ for $j=0$ and $j=2$, with prescribed initial conditions, can be evaluated based on (3.6) by using a time-marching method, e.g. the midpoint method as indicated in figure 1. With the potentials computed per time step, second-order velocity at different harmonics, $\bar {\boldsymbol {V}}^{(2j)}$ for $j=0$ and $j=2$, and their derivatives with respect to $t$, $x$, $y$ and $z$ can be readily obtained by an inverse FFT.

(Step III) The third step, as highlighted in purple in figure 1, is to evaluate $\mathcal {N}$ in the physical space based on (4.9b) using the terms evaluated in step I and step II.

Some remarks on the numerical efficiency of the new NLSE are explained here for two-dimensional waves, which allow for a fair comparison with the efficiency of existing NLSEs, e.g. the Dysthe equation as explained in Lo & Mei (Reference Lo and Mei1985). In the pseudospectral method, the nonlinear terms are evaluated through FFT at each $\boldsymbol {x}$ and $t_n$ before the integration in time is performed. The evaluation requires 20 FFTs (cf. the areas highlighted in blue and green in figure 1). The third-order nonlinear terms are evaluated in the physical domain, and an additional eight FFT computations arise from the terms associated with $\partial _t(\partial _x,\partial _z)B^{(2j)}$ and $(\partial _x,\partial _z)B^{(2j)}$ due to an inverse FFT (cf. the areas highlighted in green and purple in figure 1). Two FFTs are needed in the linear equation (i.e. the equation in the last line highlighted in grey in figure 1). If $N$ Fourier modes are used, a total of $30N\,\text {In}(N)$ operations are needed for advancing the solution by one step in time whereas $10N\,\text {In}(N)$ operations are required with NLSEs, e.g. the Dysthe equation and modified NLSE (Dysthe Reference Dysthe1979; Lo & Mei Reference Lo and Mei1985; Trulsen et al. Reference Trulsen, Kliakhandler, Dysthe and Velarde2000). Although a slightly increased computational cost relative to the previous versions is inevitable, as expected, the new NLSE offers much better computational efficiency compared with other methods, e.g. solving the Zakharov integral equation that requires roughly $O(N^3)$ operations when an explicit time-differencing scheme is used. Moreover, it is worth noting that the kinematics of the system, i.e. the elevation and velocities at the first and second orders, are obtained at almost no additional cost throughout the numerical implementations of the NLSE.

In contract to the HOS method which is also based on efficient FFTs, a pseudospectral method and numerical methods for time marching (cf. Dommermuth & Yue Reference Dommermuth and Yue1987; West et al. Reference West, Brueckner, Janda, Milder and Milton1987), the new NLSE has the following main features. First, as aforementioned, the new NLSE allows for a computational domain chosen based on the scaling of the wave envelope (or the maximum of the side bandwidth of a spectrum), in a manner similar to other NLSE-based models. This suggests a smaller number of discrete points than the HOS method can be chosen in a computational domain. It is mainly achieved by the semianalytical approach for second-order superharmonic bound waves, which are expressed in an envelope-type form. It is equivalent to shifting the spectrum of the second-order superharmonic bound waves by $2k_0$ towards the origin of the Fourier plane. In practice, the choice of the carrier wave (i.e. $k_0$) allows for flexibility and, therefore, can contribute to reducing the computational cost, in particular for the ocean spectra of a long tail. With a reasonable selection of $k_0$, the side bandwidth of an asymmetrical spectrum can be much reduced and, therefore, a less fine mesh would be allowed for the implementation of the new NLSE. For example, a numerical test for a JONSWAP spectrum (see Appendix C) was carried out with $k_0=4k_p$ chosen for computations, where $k_p$ denotes the spectrum peak. Secondly, the new NLSE separately obtains the wave fields due to free waves and second-order bound waves at a still water surface. Last, with proper separation of the perturbed solutions and multiple scales in time (in terms of wave steepness), the author conjectures that one would show that the solution of the third-order HOS method following Onorato, Osborne & Serio (Reference Onorato, Osborne and Serio2007) (i.e. equations (13) and (14) in the paper) can recover the new NLSE for the envelope of the first-harmonic wave elevation. This will be addressed in future work.

5.1. Potential $\varPhi ^{(33)}$ at third harmonic

Collecting the terms at third order in wave steepness in (2.1)–(2.3) and keeping the terms that have a factor $\exp (3\mathrm {i}\boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-3\mathrm {i}\omega _0t)$, we obtain the following equations for $\varPhi ^{(33)}$:

(5.1a,b)\begin{equation} \nabla^2_3\varPhi^{(33)} = 0\quad \text{for}\ -\infty< z<0,\quad \varGamma\varPhi^{(33)} = Q^{(33)}\quad \text{for}\ z=0 \end{equation}

and

(5.1c)\begin{equation} \partial_z\varPhi^{(33)} =0\quad \text{for}\ z\to-\infty, \end{equation}

where the third-order term $Q^{(33)}(\boldsymbol {x},t)$ denotes the forcing of bound waves at third harmonic on a still water surface, as noted in § 4, expressed as

(5.2a)\begin{equation} Q^{(33)} =\bar{Q}^{(33)}(\boldsymbol{x},z,t) \,\mathrm{e}^{3\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{x} - 3\mathrm{i}\omega_0t}\quad \text{for}\ z=0 \end{equation}

and

(5.2b)\begin{equation} \bar{Q}^{(33)} = \textstyle\frac{1}{2}(\partial_t-3\mathrm{i}\omega_0)[\bar{\boldsymbol{V}}\boldsymbol{\cdot} (\nabla_3 + 2\mathrm{i} \boldsymbol{k}_0^{(3d)} )B^{(22)} ] + \textstyle\frac{1}{16} \bar{\boldsymbol{V}}\boldsymbol{\cdot} [(\nabla_3 + 2\mathrm{i} \boldsymbol{k}_0^{(3d)} )(\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}) ]. \end{equation}

For unidirectional waves, it is understood that both $B^{(22)}$ and $(\bar {\boldsymbol {V}}\boldsymbol {\cdot }\bar {\boldsymbol {V}})$ have no contribution to the third-order nonlinear terms as they vanish, meaning $\bar {Q}^{(33)} =0$ and thus, $\varPhi ^{(33)} = 0$.

The solution of equations (5.1a,b,c) for $\varPhi ^{(33)}$ can be obtained in a way similar to the semianalytical approach for $\varPhi ^{(2)}$, i.e. the solution can be obtained by a finite-difference method in the Fourier $\boldsymbol {k}$ space based on the forcing equation on a still water surface, as follows:

(5.3)\begin{equation} (\partial_{t}-3\mathrm{i}\omega_0)^2 \hat{B}^{(33)}(\boldsymbol{k},z,t) + g|\boldsymbol{k}+3\boldsymbol{k}_0|\hat{B}^{(33)}(\boldsymbol{k},z,t) = 2\hat{Q}^{(33)}(\boldsymbol{k},z,t)\quad \text{for}\ z=0, \end{equation}

where, as defined in previous sections, ‘$\hat {(\ldots )}$’ denotes a Fourier transform to the $\boldsymbol {k}$ plane with respect to $\boldsymbol {x}$. Thereby, $\varPhi ^{(33)}$ is expressed as

(5.4a)\begin{equation} \varPhi^{(33)}= \textstyle\frac{1}{2} B^{(33)} \,\mathrm{e}^{3k_0z}\,\mathrm{e}^{3\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{x} - 3\mathrm{i}\omega_0t}, \end{equation}

with

(5.4b)\begin{equation} B^{(33)}=\int_{-\infty}^{\infty} \hat{B}^{(33)}(\boldsymbol{k},0,t) \,\mathrm{e}^{|\boldsymbol{k}+3\boldsymbol{k}_0|z-3k_0z}\,\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d} \boldsymbol{k}^2. \end{equation}

5.2. Elevation $\zeta ^{(3)}$ at third order

Based on the boundary conditions (2.2a,b) it is understood that the wave elevation at third order can be obtained from (cf. (13.2.3) by Mei et al. (Reference Mei, Stiassnie and Yue2005))

(5.5)\begin{align} \zeta^{(3)} &={-}\frac{1}{g}\left[ \partial_t \varPhi^{(3)} + \zeta^{(1)} \partial_{tz}\varPhi^{(2)} + \zeta^{(2)} \partial_{tz}\varPhi^{(1)}+ \nabla_3\varPhi^{(1)}\boldsymbol{\cdot} \nabla_3\varPhi^{(2)} \right.\nonumber\\ &\quad + \left.\frac{1}{2}(\zeta^{(1)})^2\partial_{tzz}\varPhi^{(1)} + \frac{1}{2}\zeta^{(1)}\partial_z(\nabla_3\varPhi^{(1)}\boldsymbol{\cdot} \nabla_3\varPhi^{(1)}) \right]. \end{align}

Based on (5.5), it is understood $\zeta ^{(3)}$ can be expressed in a form, as follows:

(5.6)\begin{equation} \zeta^{(3)} =( \textstyle\frac{1}{2} A^{(31)}(\boldsymbol{x},t) \,\mathrm{e}^{\mathrm{i} (\boldsymbol{k}_0\cdot\boldsymbol{x}-\omega_0t)} + \text{c.c.} )+ ( \textstyle\frac{1}{2} A^{(33)}(\boldsymbol{x},t) \,\mathrm{e}^{3\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{x} - 3\mathrm{i}\omega_0t} + \text{c.c.}), \end{equation}

where amplitudes $A^{(31)}(\boldsymbol {x},t)$ and $A^{(33)}(\boldsymbol {x},t)$ are in a form given by, respectively,

(5.7)\begin{align} A^{(31)}& ={-} \frac{1}{2g}[(\partial_t-2\mathrm{i}\omega_0)(2k_0+\partial_z)B^{(22)} ]A^* -\frac{1}{g} A \partial_{tz}\varPhi^{(20)} - \frac{1}{2g} A^{(22)}[(\partial_t+\mathrm{i}\omega_0)\bar{w}^*]\nonumber\\ &\quad - \frac{1}{g} \zeta^{(20)}[(\partial_t-\mathrm{i}\omega_0)\bar{w}] - \frac{1}{2g}\bar{\boldsymbol{V}}^*\boldsymbol{\cdot} [(2\boldsymbol{k}^{(3d)}_0+\nabla_3)B^{(22)} ] - \frac{A^2}{8g} (\partial_t+\mathrm{i}\omega_0)(k_0+\partial_{z})\bar{w}^*\nonumber\\ &\quad - \frac{1}{4g} AA^*(\partial_t-\mathrm{i}\omega_0)(k_0+\partial_{z})\bar{w} -\frac{1}{8g} A^* [(\partial_z+2k_0)(\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}})]\nonumber\\ &\quad - \frac{1}{4g} A(\partial_z+2k_0)(\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}^*) -\frac{1}{g}(\partial_{t}+\mathrm{i}\omega_0\mathcal{L} )B, \quad \text{for}\ z=0, \end{align}

(5.8)\begin{align} A^{(33)} &={-}\frac{1}{g}(\partial_{t} -3\mathrm{i}\omega_0)B^{(33)}- \frac{1}{2g}[(\partial_t-2\mathrm{i}\omega_0)(2k_0+\partial_z)B^{(22)} ]A - \frac{1}{2g} A^{(22)}[(\partial_t+\mathrm{i}\omega_0)\bar{w}]\nonumber\\ &\quad - \frac{1}{2g}\bar{\boldsymbol{V}}\boldsymbol{\cdot} [(2\boldsymbol{k}^{(3d)}_0+\nabla_3)B^{(22)} ] - \frac{A^2}{8g}(\partial_t-\mathrm{i}\omega_0)(k_0+\partial_{z})\bar{w}\nonumber\\ &\quad-\frac{1}{8g} A [(\partial_z+2k_0)(\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}})], \end{align}

in which the nonlinear terms are obtained in the physical plane using a pseudospectral method. As noted in §§ 2.4 and 3.3.1, both $B^{(22)}$ and $\bar {\boldsymbol {V}}\boldsymbol {\cdot }\bar {\boldsymbol {V}}$ vanish for unidirectional waves, under which circumstances the amplitudes at third order embrace a much simpler expressions given by

(5.9a)\begin{align} A^{(31)}_{uni} &={-}\frac{1}{g} A \partial_{tz}\varPhi^{(20)} -\frac{1}{2g} A^{(22)}(\partial_t+\mathrm{i}\omega_0)\bar{w}^*-\frac{1}{g} \zeta^{(20)} (\partial_t-\mathrm{i}\omega_0)\bar{w}\nonumber\\ &\quad - \frac{A^2}{8g} (\partial_t+\mathrm{i}\omega_0)(k_0+\partial_{z})\bar{w}^* - \frac{1}{4g} AA^*(\partial_t-\mathrm{i}\omega_0)(k_0+\partial_{z})\bar{w}\nonumber\\ &\quad - \frac{1}{4g} A[(\partial_z+2k_0)(\bar{\boldsymbol{V}}\boldsymbol{\cdot}\bar{\boldsymbol{V}}^*)]-\frac{1}{g}(\partial_{t}+\mathrm{i}\omega_0\mathcal{L} )B, \end{align}

(5.9b)\begin{align} A^{(33)}_{uni} &={-} \frac{A^2}{8g} (\partial_t-\mathrm{i}\omega_0)(k_0+\partial_{z})\bar{w}-\frac{1}{2g} A^{(22)}(\partial_t-\mathrm{i}\omega_0)\bar{w}, \end{align}

in which the subscript ‘uni’ refers to the cases for unidirectional waves. With a narrow-banded assumption $O(\delta )\ll 1$, we readily obtain from (5.9) (through omitting terms associated with the derivatives with respect to $t$ and $z$ and the second-order mean fields and noting $\bar {w}= -\mathrm {i} \omega _0 A + O(\epsilon \delta )$), as follows:

(5.10a,b)\begin{equation} A^{(31)}_{{uni},\approx} ={-}\textstyle\frac{3}{8} k_0^2|A|^2A + O(\epsilon^3\delta)\quad \text{and}\quad A^{(33)}_{{uni},\approx} = \textstyle\frac{3}{8} k_0^2A^3 + O(\epsilon^3\delta), \end{equation}

in which the subscript ‘$\approx$’ denotes an approximate expression due to a narrow-banded assumption. Inserting (5.10a,b) into (5.6) leads to a leading-order approximation for $\zeta ^{(3)}$ that agrees with (2.12) in Lo & Mei (Reference Lo and Mei1985).