Three-dimensional surface gravity waves of a broad bandwidth on deep water

Abstract

A new nonlinear Schrödinger equation (NLSE) is presented for ocean surface waves. Earlier derivations of NLSEs that describe the evolution of deep-water waves have been limited to a narrow bandwidth, for which the bound waves at second order in wave steepness are described in leading-order approximations. This work generalizes these earlier works to allow for deep-water waves of a broad bandwidth with large directional spreading. The new NLSE permits simple numerical implementations and can be extended in a straightforward manner in order to account for waves on water of finite depth. For the description of second-order waves, this paper proposes a semianalytical approach that can provide accurate and computationally efficient predictions. With a leading-order approximation to the new NLSE, the instability region and energy growth rate of Stokes waves are investigated. Compared with the exact results based on McLean (J. Fluid Mech., vol. 511, 1982, p. 135), predictions by the new NLSE show better agreement than by Trulsen et al. (Phys. Fluids, vol. 12, 2000, pp. 2432–2437). With numerical implementations of the new NLSE, the effects of wave directionality are investigated by examining the evolution of a directionally spread focused wave group. A downward shift of the spectral peak is observed, owing to the asymmetry in the change rate of energy in a more complex manner than that for uniform Stokes waves. Rapid oblique energy transfers near the group at linear focus are observed, likely arising from the instability of uniform Stokes waves appearing in a narrow spectrum subject to oblique sideband disturbances.

1. Introduction

An important task in the study of surface gravity waves is the development of the theoretical description of flow fields. Theoretical models of surface gravity waves are essential in a wide range of fields, such as in engineering practices and for research purposes. A theoretical study has obvious advantages in elucidating the underlying physics, thereby advancing the understanding of realistic wave problems, compared with other approaches such as experiments, field studies and direct numerical computations. Theoretical findings have contributed to providing possible explanations of the formation mechanism of extremely large waves that appear suddenly with much larger amplitude than their surroundings, known also as ‘rogue’ or ‘freak’ waves. A few examples of the possible mechanisms proposed are modulational or Benjamin–Feir instability in deep water (Onorato et al. 2009), refraction by ambient currents or bathymetry (Janssen & Herbers 2009; Onorato, Proment & Toffoli 2011) and effects on weakly nonlinear waves of a depth change in shallow or intermediate water depth (Trulsen et al. 2020; Li et al. 2021a,b).

The modulational instability of weakly nonlinear Stokes waves subject to sideband disturbances was discovered by Benjamin & Feir (1967). Theoretical advances made in the late 1960s in this perspective include contributions by Lighthill (1965), Benjamin (1967), Whitham (1967), Zakharov (1968) and Benney & Roskes (1969). Since the discovery of the instability of a train of Stokes waves, progress has been made both numerically (Janssen 1983; Lo & Mei 1985; Trulsen & Dysthe 1997; Dysthe et al. 2003) and experimentally (Lake et al. 1977; Melville 1982; Su 1982; Chabchoub, Hoffmann & Akhmediev 2011) in understanding the evolution properties of surface waves on deep water. For smaller values of wave steepness, symmetrical upper and lower sidebands tend to grow equally in the initial stage. This is due to the degenerate resonant interaction in the prestigious ‘figure of eight’ quartet resonance loop (Phillips 1960, 1967; Longuet-Higgins 1976; Lake et al. 1977; McLean 1982), which corresponds to the special cases where two out of the four wavenumbers obeying the ‘figure of eight’ loop are identical. As nonlinearity increases, the sidebands appear to grow unequally, with the lower sideband growing faster than the upper sideband (Lake et al. 1977; Lo & Mei 1985), reaching a maximum larger than the minimum to which the spectrum peak drops. This unequal growth in energy, as a result of a combination of dissipation, wave breaking and nonlinear wave evolution, likely causes the downward shift of the spectral peak of wind waves (Lake et al. 1977; Lo & Mei 1985; Trulsen & Dysthe 1997).

The so-called nonlinear Schrödinger equation (NLSE) has been widely known as a convenient approach for analysing the instability of Stokes waves. For deep-water waves, Zakharov (1968) found that the wave envelope satisfies the NLSE which is cubic in wave slope $\epsilon$ (or wave steepness). The NLSE is derived based on the assumption of a small wave slope and the so-called narrow-banded assumption that restricts the modulation of the wave envelope to be slow in both space and time relative to a rapidly varying wave phase.

Many attempts have been made to expand the applicability scope of the NLSE for deep-water waves through adding higher-order terms, with a primary focus on reducing the restrictions due to bandwidth, as shown in table 1. Table 1 indicates the lower-order wave fields considered in a NLSE for the nonlinear and dynamic evolution of the (potential or elevation) envelope of the first harmonic. Specifically, the second-order wave fields do not appear alone in a NLSE but in combination with the first-order wave fields. Let $\delta _{x,y}$ be a small non-dimensional parameter as a measure of the bandwidth in the main propagation ($\delta _x$) and the transverse ($\delta _y$) directions, respectively. With one further step, Dysthe (1979) has obtained a NLSE that is correct to $O(\epsilon ^3\delta _{x,y}, \epsilon \delta ^n_{x}\delta _y^{3-n})$ ($n=0,1,2,3$), known as Dysthe's equation, or the fourth-order NLSE, as $O(\delta _{x,y})\sim O({\epsilon })$ is implied. With the fourth-order NLSE, Dysthe (1979) found that second-order mean flows play a significant role in the instability of Stokes waves. It is shown by Stiassnie (1984) that Dysthe's equation can be derived from Zakharov's integral equation (Zakharov 1968) through a narrow-banded assumption. Higher-order extensions in bandwidth are derived by Trulsen & Dysthe (1996) with a modified NLSE correct to $O(\epsilon \delta ^n_{x}\delta _y^{5-n},\epsilon ^3\delta _{x,y})$ ($n=0,1,\ldots ,5$) and by Trulsen et al. (2000) in which the linear evolution of the envelope of linear wave is described by the exact linear dispersion relation, giving an NLSE correct to $O(\epsilon \delta _{x,y}^\infty ,\epsilon ^3\delta _{x,y})$. Using some of the earlier versions of NLSEs, Trulsen & Dysthe (1997) and Dysthe et al. (2003) found that the NLSEs can recover a downward shift of the spectral peak in a three- and two-dimensional narrow spectrum, respectively. The latter occurs after the spectrum has reached a quasi-steady state and hence, is unlikely due to the modulational instability which results in symmetrical growth of sidebands. Strong frequency dependence developed in the temporal evolution of directional waves initially of no frequency dependence is observed in NLSE-based numerical simulations and field data (Simanesew et al. 2016). Simanesew et al. (2017) found the modified NLSE of Trulsen et al. (2000) gives less accurate predictions for short crest waves with larger directional spread.

Table 1. A summary of NLSEs for waves on deep water and the regimes of applicability which refer to the accuracy of lower-order wave fields considered in an NLSE for the dynamic evolution of the (potential or elevation) envelope of the first harmonic, in addition to the NLSEs for an arbitrary depth. Here $\epsilon$ denotes the dimensionless wave steepness; $\delta _x$ and $\delta _y$ denote the dimensionless bandwidth in the main propagation direction and in the direction normal to the main, respectively; $\delta ^\infty _{x,y}$ denotes arbitrary order in bandwidth.

Recent progress has also been made in generalizing an NLSE in order to take into account the interaction of waves with ambient environments, e.g. dissipation and forcing effects due to wind actions and turbulence (cf. Wu, Liu & Yue 2006; Dias, Dyachenko & Zakharov 2008; Kharif et al. 2010; Slunyaev, Sergeeva & Pelinovsky 2015) and wave current interaction (Dysthe & Das 1981; Stocker & Peregrine 1999; Hjelmervik & Trulsen 2009; Curtis, Carter & Kalisch 2018). In addition to the NLSEs for deep-water waves, coupled NLSEs are derived to consider crossing sea states (see Gramstad & Trulsen (2011), Trulsen et al. (2015) and references therein). A generalization to water of finite uniform depth is done by Davey & Stewartson (1974) and higher-order corrections are added by Johnson (1977), Brinch-Nielsen & Jonsson (1986), Slunyaev (2005) and Gramstad (2014), among others. Modified cubic equations are obtained that allow for waves in water of a finite, varying depth due to a mildly sloping seabed (Djordjevié & Redekopp 1978; Kirby & Dalrymple 1983).

Many other approaches have been developed as alternative tools to the NLSE for understanding the properties of surface deep-water waves, such as the Zakharov integral equation and its leading-order approximations (see Zakharov 1968; Crawford, Saffman & Yuen 1980; Crawford et al. 1981; Krasitskii 1994), numerical wave tanks (e.g. Bateman, Swan & Taylor 2001; Ducrozet et al. 2012) based on the high-order spectral (HOS) method (Dommermuth & Yue 1987; West et al. 1987) and direct numerical solutions of fully nonlinear potential-flow equations (Engsig-Karup, Bingham & Lindberg 2009; Bihs et al. 2020; Zheng et al. 2020). Through studying the temporal evolution of three-dimensional steep focused wave groups, rapid energy changes in a wave spectrum are reported by earlier works e.g. Gibbs & Taylor (2005) and Gibson & Swan (2007). They attribute this phenomenon to the third-order resonant interactions, which are considered as the cause for the formation of the so-called ‘wing waves’ that appear in the transverse direction after the group at nonlinear focus reported in Barratt, Bingham & Adcock (2020). Direct numerical solutions were used by Barratt et al. (2020) for the evolution of a steep focused wave group with directional spread. In three dimensions, oblique energy transfer in a wave spectrum and a downward shift of the spectral peak are observed in previous works, e.g. Trulsen & Dysthe (1997) and Barratt et al. (2021).

There is still room and the need for further development of the models based on an NLSE, with two aspects addressed here. Firstly, ocean surface waves are generally three-dimensional and broad banded. Nevertheless, it should be noted that the NLSEs mentioned above, e.g. Dysthe (1979) and Trulsen et al. (2000), are only correct to order $(\epsilon ^3\delta _{x,y})$, as shown in table 1. This suggests the need for extending the bandwidth at third order in wave steepness. Secondly, surface waves interact with ambient environments such as subsurface currents, turbulence and wind actions in the atmosphere, which often require explicit and accurate vertical structures of flow fields below the surface to allow for coupling. In contrast, it is known that a NLSE is based on the description of the envelope of the surface displacement, accompanied by mostly lowest-order approximations to the vertical profiles of flow fields at second order in wave steepness, as shown in table 1.

This paper aims to make attempts to fill in gaps in the aforementioned two aspects, through the development of a new framework that would allow for the study of three-dimensional waves of a broad bandwidth and would take into account the structures of flow fields below the surface without compromising much the computational efficiency of an NLSE-based model. The objective of this paper is threefold. First, a new NLSE for the evolution of three-dimensional deep-water surface waves is derived in §§ 2–4 that does not rely on the assumption of a narrow bandwidth as aforementioned NLSEs. Specially, the newly developed NLSE should allow for $\delta _{x,y}>1$ since it has relaxed the narrowband assumption for wave fields at second order, as shown in table 1. Simple numerical implementations of the new NLSE as performed with earlier versions of NLSEs are explained. Second, a semianalytical approach for the description of wave fields for bound waves at second order in wave steepness is proposed in § 3 which allows for simpler and more efficient numerical implementations than by Dalzell (1999). Finally, a study of the instability region and energy growth rate of Stokes waves and the temporal evolution of a directionally spread focused wave group are presented in § 6, where comparisons with Longuet-Higgins (1978), Dysthe (1979), Crawford et al. (1981), McLean (1982) and Trulsen et al. (2000) are made.

2. Mathematical formulation and methodology

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 1974) 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 1971; Davey & Stewartson 1974).

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. (2000) 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 (1979) and Stiassnie (1984). 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 2005), 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 (1999). 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. Second-order solutions $O(\epsilon ^2)$

In this section, the solutions for the second-order superharmonic and subharmonic waves forced by linear surface waves are presented. Compared with Dalzell (1999), the solutions in this section are as exact but expressed in an envelope-type form. In contrast to Chu & Mei (1971), Davey & Stewartson (1974) and Dysthe (1979), the derivations do not rely on a narrow-banded assumption except for § 3.3 where leading-order approximations are presented.

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 1999; Li et al. 2021c):

(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 1960; Hasselmann 1962). 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 (1996), periodic boundary conditions as in Dommermuth & Yue (1987), the framework by Bonnefoy, Le Touzé & Ferrant (2006) for the waves generation in a numerical wave tank, stationary waves based on Dalzell (1999) 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 1996). The semianalytical approach leads to improved computational efficiency compared with the analytical method by Dalzell (1999) based on Fourier integrals. For $N$ Fourier modes, the semianalytical approach requires computational operations at $O(N\,\text {In}(N))$, whereas Dalzell (1999) 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. A new NLSE for $A$ ($O(\epsilon ^3)$)

In this section, a new NLSE for linear wave envelope $A$ is presented in § 4.1 for multidirectional waves of a broad bandwidth. Based on a narrow-banded assumption, leading-order approximations to the new NLSE are derived in § 4.2, and how this new evolution equation recovers existing NLSEs is explained in § 4.3. The numerical implementations of the new NLSE are presented in § 4.4.

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 (2006)). 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 (1979), Trulsen & Dysthe (1996) and Trulsen et al. (2000); 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. (2008) and Kharif et al. (2010), 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. (2008)) 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 2004). 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 1971; Li et al. 2021c), 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 (1979), Trulsen & Dysthe (1996) and Trulsen et al. (2000) 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 (1967) and Zakharov (1968). 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 (1979), Trulsen & Dysthe (1996) and Trulsen et al. (2000) 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. (2000) using the dimensional versions. It has been demonstrated in Trulsen et al. (2000) how this term can recover the terms associated with linear dispersion relation derived by Trulsen & Dysthe (1996) and Dysthe (1979). Second, the nonlinear term in the NLSEs by Trulsen et al. (2000) and Trulsen & Dysthe (1996) are directly obtained from the Dysthe equation by Dysthe (1979). 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 (1979) 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. (2005) 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 (1979); i.e. (4.17) recovers ‘$\varGamma$’ on the right-hand side of equation (2.17) of Dysthe (1979), while observing the following two aspects. Firstly, the notation $C$ defined here is denoted by ‘$A$’ in Dysthe (1979). Secondly, there is one misprint in the equation for ‘$\varGamma$’, as also pointed out by Brinch-Nielsen & Jonsson (1986) (cf. p. 461) and Janssen (1983) – 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. (2000) (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 1998), 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 1974). The main difference between the numerical procedures presented here and these presented in Lo & Mei (1985), 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 (1985). 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 1979; Lo & Mei 1985; Trulsen et al. 2000). 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 1987; West et al. 1987), 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 (2007) (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. Solutions at third order $O(\epsilon ^3)$

For completeness, this section presents the solutions at third order for both potential $\varPhi ^{(33)}$ given in §5.1 and elevation $\zeta ^{(3)}$ in § 5.2.

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. (2005))

(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 (1985).

6. Results

This section focuses on three aspects in order to explore and validate the theory presented, including the second-order solutions presented in § 3, sideband instability of Stokes waves and the roles of wave directionality. They are examined in §§ 6.1, 6.2 and 6.3, respectively.

For numerical implementations, a Gaussian amplitude spectrum ($\varXi (k)$) in wavenumber and a (unnormalized) directional distribution ($D(\theta )$) were chosen to generate a short (focused) wave group, given by, respectively,

(6.1a,b)\begin{equation} {\varXi}(k) = \exp\left[- \frac{(k-k_p)^2}{2k^2_w}\right] \quad \text{and}\quad D(\theta) = \exp\left[- \frac{(\theta-\theta_p)^2}{2\theta^2_w}\right], \end{equation}

where $k_p$ denotes the peak wavenumber and $\theta _p$ denotes the propagation direction of the peak wave at linear focus, $k_w$ and $\theta _w$ denote the standard deviation of the spectrum in wavenumber and direction, respectively. An asymmetrical Gaussian spectrum is obtained from (6.1a) through using two different values of $k_w$ for the upper sideband $k>k_p$ and lower sideband $k\leq k_p$, respectively. For the relevant cases computed in this section, we used $\theta _p=0$ and $k_p= 0.02769$ m$^{-1}$ which is typical of the North Sea (Barratt et al. 2021). Using the spectra described by (6.1a,b), linear elevation $\zeta ^{(1)}$ at $t=t_0$ is given as input by

(6.2a)\begin{equation} \zeta^{(1)}(x,t) = \frac{\epsilon_0 }{k_0 {\varXi}_{k}}\sum^{N_k}_{n=1} {\varXi}(k_n)\cos[k_n (x-x_{{f}}) - \sqrt{gk_n}(t_0-t_{f})+\psi_{{f}}], \end{equation}

with

(6.2b)\begin{equation} {\varXi}_{k} = \sum^{N_k}_{n=1} {\varXi}(k_n), \end{equation}

for unidirectional waves, where $\epsilon _0$ denotes the dimensionless steepness of a wave group at linear focus; $k_n$ and $N_k$ denote the evenly spaced discrete wavenumbers and the total number chosen for numerical implementations, respectively; the subscript ‘$f$’ denotes the position or time for the wave group at linear focus; $\psi _{{f}}$ denotes the wave phase at linear focus; for multidirectional waves,

(6.3a)\begin{equation} \zeta^{(1)}(\boldsymbol{x},t) = \frac{\epsilon_0 }{k_0 \varXi_{k,\theta}}\sum^{N_k}_{n=1} \sum ^{M_\theta}_{m= 1} \varXi(k_n)D(\theta_m) \cos[\boldsymbol{k}_{n,m}\boldsymbol{\cdot}(\boldsymbol{x}-\boldsymbol{x}_{{f}}) - \omega({\boldsymbol{k}_{n,m}})(t_0-t_{{f}})+\psi_{{f}}],\end{equation}

with

(6.3b)\begin{equation} \varXi_{k,\theta} = \sum^{N_k}_{n=1} \sum ^{M_\theta}_{m= 1} \varXi(k_n)D(\theta_m) \quad \text{and}\quad \boldsymbol{k}_{n,m} = [k_n\cos\theta_m, k_n\sin\theta_m], \end{equation}

where $\theta _m$ and $M_\theta$ denote the evenly spaced discrete wave directions in the range of $(-{\rm \pi} , {\rm \pi})$ and the total number chosen for numerical implementations, respectively. Despite that the directional distribution described by (6.1b) is in an unnormalized form, (6.3a,b) suggest that the integral unity of the directions of the linear multidirectional waves generated has implicitly been satisfied due to the additional scaling factor, $\varXi _{k,\theta }$, defined in (6.3b).

6.1. Second-order solutions

We validate the second-order solutions presented in § 3 in this section through comparisons with the analytical method by Dalzell (1999) and with existing approximate methods. We focus on the temporal–spatial evolution of a wave group. For simplicity, the effects at third order are neglected in this section as they do not affect, qualitatively, the discussion presented here. Let a wave group start to propagate at $t=t_0$ at which the linear elevation $\zeta ^{(1)}(\boldsymbol {x},t_0)$ of a group is given in space by (6.2) for a unidirectional group and by (6.3a) for a directionally spread focused wave group. Inserting $\zeta ^{(1)}(\boldsymbol {x},t_0)$ into (2.6) leads to an expression for $A(\boldsymbol {x},t)$. For a fair comparison, the input at $t=t_0$ for second-order wave fields is based on the results by Dalzell (1999). Following the numerical procedures explained in § 4.4, $\zeta ^{(22)}$ and $\zeta ^{(20)}$ are readily obtained from (3.13b) and (3.13c), respectively.

6.1.1. Unidirectional and multidirectional wave group

Figure 2 shows a comparison of the second-order superharmonic wave elevation for a unidirectional wave group at two different times (noting that $\varPhi ^{(22)}(x,t)=0$), between the results predicted by (3.13b), Dalzell (1999) and a leading-order approximation based on a second-order Stokes theory given by

(6.4)\begin{equation} \zeta^{(22)}_{{\text{Stokes}}} = \textstyle\frac{1}{2} k_0A^2(\boldsymbol{x},t)\cos[k_0(x-x_{{f}})-\omega_0(t-t_{{f}})+\psi_0]. \end{equation}

Due to an asymmetrical Gaussian amplitude spectrum chosen in figure 2, a good estimate of the side bandwidth can be given by $\delta = 3k_w/k_p$ for $k>k_p$, which measures from the wavenumber (i.e. $k=k_p+3k_w$) dropping by $99\,\%$ in magnitude relative to the spectrum's peak. This gives $\delta = 4.86$ for the case examined in figures 2 and 3. The good agreement between (3.13b) and Dalzell (1999) is evident in figure 2, whereas the less satisfactory agreement between the second-order superharmonic elevation (6.4) and the other two methods is also shown. For completeness, the expressions for both $\zeta ^{(22)}$ and $\zeta ^{(20)}$ based on Dalzell (1999) are presented in Appendix B. It is understood Dalzell (1999) is capable of providing exact predictions of second-order wave fields based on Fourier integrals. The agreement with Dalzell (1999) in figure 2 clearly demonstrates that the semianalytical approach presented in § 3.2 can also work for this purpose. Similar observations are also shown in figure 4 for $\zeta ^{(22)}$ predicted for the evolution of a directionally spread focused wave group at different times.

Figure 2. A comparison of second-order superharmonic wave elevation $\zeta ^{(22)}$ for a unidirectional wave group at different times between the predictions by (3.13b) (black solid line), Dalzell (1999) (red dashed line) and (6.4) (blue dot–dashed line) based on a second-order Stokes wave theory. The wave group started to propagate at an initial time $t_0 = -15 T_0$ with $T_0$ the period of the peak wave; it focuses linearly at $x_{f}=-25\lambda _0$ with $\lambda _0= 2{\rm \pi} /k_p$, $k_0=k_p$ and $t_{{f}} = -10 T_0$. At $t=t_0$, an asymmetrical Gaussian spectrum, with $k_w = 0.27k_0$ for $k< k_0$ and $k_w=1.62 k_0$ for $k>k_0$ and $\epsilon _0=0.3$, was chosen for $\zeta ^{(1)}$ based on (6.2). Panels (a,b) show for $t=-10T_0$, (c,d) for $t=20 T_p$ and (e,f) for $t=40T_0$. Panels (b,d,f) show the area highlighted by the thick blue box in (a,c,e), respectively.

Figure 3. A comparison of surface wave elevation $\zeta ^{(20)}$ for the second-order subharmonic waves forced by a unidirectional wave group at different times between the prediction by (3.13c) (black solid line), Dalzell (1999) (red dashed line) and Trulsen et al. (2000) (blue dot–dashed line). In the figure, panels (a,b) show for $t= -10T_0$, panels (c,d) for $t = 20T_0$ and panels (e,f) $t = 40T_0$; the parameters were chosen the same as in figure 2.

Figure 4. Second-order superharmonic wave elevation $\zeta ^{(22)}$ for a directionally spread focused wave group at different times; the left-hand subpanels show the results predicted by (3.13b); the right-hand subpanels show a comparison of the prediction by (3.13b) (blue dot–dashed line) and Dalzell (1999) (red solid line) for $y=0$. The wave group started to propagate at an initial time $t_0 = -15 T_0$ with $T_0$ denoting the period of the peak wave $k_p=k_0$; it focuses linearly at $x_{f}=0$ and $t_{{f}} = 0$. At $t=t_0$, a symmetrical Gaussian spectrum with $k_w = 0.3k_0$, $\epsilon _0=0.3$, $\theta _w=15^\circ$ was chosen for $\zeta ^{(1)}$ based on (6.3a). In the figure, $\lambda _0=2{\rm \pi} /k_0$ denotes the wavelength of the carrier wave of the group. Here (a) $t=-15T_0$, (b) $t=-1.74T_0$, (c) $t=11.5T_0$ and (d) $t=14.8T_0$.

Figure 3 shows a comparison of the elevation for second-order subharmonic waves forced by a unidirectional linear wave group on a still water surface between the results predicted by (3.13c), Dalzell (1999) (cf. Appendix B) and Trulsen et al. (2000) (cf. Appendix D). It is clearly seen in figure 3 that (3.13c) can predict $\zeta ^{(20)}$ as well as Dalzell (1999), whereas Trulsen et al. (2000) leads to less satisfactory agreement with Dalzell (1999), owing to a narrow banded assumption.

6.2. Sideband instability of Stokes waves

In this section, we examine the capabilities of the new NLSE by investigating the sideband instability of Stokes waves through comparisons with earlier works, e.g. Dysthe (1979), Crawford et al. (1981), McLean (1982) and Trulsen et al. (2000). The approximations to $\mathcal {N}$ presented in § 4.2 are used in this section as they are sufficient to this end. The contribution due to $\mathcal {N}_{\approx , {dir}}$ is neglected for simplicity. Based on (4.15), the stability of Stokes waves to sideband perturbations can be investigated by assuming small disturbances in amplitude and phase for $A$ in a form given by (cf. Trulsen et al. 2000)

(6.5)\begin{equation} A = a_0(1+a+\mathrm{i} \theta)\,\mathrm{e}^{-\mathrm{i}\alpha_1\epsilon_0^2\omega_0t}\quad \text{with}\ [a,\theta] =2 [\hat{a},\hat{\theta}]\cos\psi(\boldsymbol{x},t), \end{equation}

in which $a_0$ denotes the amplitude of a train of Stokes wave, $\epsilon _0=k_0a_0$ denotes the dimensionless steepness of the Stokes wave, $\hat {a}$ and $\hat {\theta }$ are infinitesimal real parameters (i.e. $\ll O(1)$) that denote small perturbations in amplitude and phase, respectively, $\psi (\boldsymbol {x},t)=k_0[(1+\delta _x)x + \delta _y y] -\omega _0\delta _\varOmega t$ denotes the modulated phase, with $k_0$ and $\omega _0$ the wavenumber and angular frequency of the Stokes wave, respectively, $\delta _x$ and $\delta _y$ the dimensionless sideband in wavenumber in the longitudinal ($k_x$) and transverse/spanwise ($k_y$) direction, respectively, and $\delta _\varOmega$ the dimensionless sideband in frequency; as defined in § 4.3.1, $\alpha _1$ is defined as $\alpha _1 =1/2$.

Treating the expression (6.5) for $A$ as the modulated amplitude of a plane wave $\boldsymbol {k}_0$ and inserting (6.5) into (2.12) leads to $\bar {\boldsymbol {V}}(\boldsymbol {x},z,t)$ given by

(6.6) \begin{equation} \bar{\boldsymbol{V}}={-}\mathrm{i}\frac{\omega_0 }{k_0}\epsilon_0 \left\{\left[\begin{array}{@{}c@{}} \mathrm{i} \\ 0 \\ 1 \end{array}\right] +\left[\begin{array}{@{}c@{}} (\mathrm{i} \hat{a} -\hat{\theta} ) V_{x,+}(z)\cos\psi +( \hat{a}+\mathrm{i}\hat{\theta}) V_{x,-}(z)\sin\psi \\ (\mathrm{i} \hat{a} -\hat{\theta} ) V_{y,+}(z)\cos\psi +( \hat{a}+\mathrm{i}\hat{\theta}) V_{y,-}(z)\sin\psi \\ (\hat{a} +\mathrm{i}\hat{\theta} ) V_{z,+}(z)\cos\psi + (\mathrm{i}\hat{a} -\hat{\theta} ) V_{z,-}(z)\sin\psi \end{array}\right]\right\} \mathrm{e}^{-\mathrm{i}\alpha_1|a_0|^2t}, \end{equation}

where the $z$ dependent coefficients $V_{i,\pm }(z)$ are

(6.7a)\begin{gather} V_{x,\pm}(z) = \frac{\sqrt{1+\delta_{z,+}}}{\delta_{z,+}}(1+\delta_x)\,\mathrm{e}^{k_0z\delta_{z,+}} \pm \frac{\sqrt{1+\delta_{z,-}}}{\delta_{z,-}}(1-\delta_x)\,\mathrm{e}^{k_0z\delta_{z,-}}, \end{gather}

(6.7b)\begin{gather}V_{y,\pm}(z) = \frac{\sqrt{1+\delta_{z,+}}}{\delta_{z,+}}\delta_y\,\mathrm{e}^{k_0z\delta_{z,+}} \pm \frac{\sqrt{1+\delta_{z,-}}}{\delta_{z,-}}(-\delta_y)\,\mathrm{e}^{k_0z\delta_{z,-}} \end{gather}

and

(6.7c)\begin{equation} V_{y,\pm}(z) = ~ \sqrt{1+\delta_{z,+}}\,\mathrm{e}^{k_0z\delta_{z,+}} \pm \sqrt{1+\delta_{z,-}}\,\mathrm{e}^{k_0z\delta_{z,-}}, \end{equation}

with

(6.7d)\begin{equation} \delta_{z,\pm} =~ \sqrt{(1\pm \delta_x)^2+\delta_y^2} -1. \end{equation}

The full solution for $\bar {\boldsymbol {V}}$ in a form as (6.6) allows for $\min (|\delta _x|, |\delta _y|)>1$. It should be noted that, time dependant disturbances (small) for $\bar {\boldsymbol {V}}$ except for these due to linear dispersion relation are neglected for simplicity. Otherwise, the derivations following Benjamin (1967) (e.g. starting from equations (11) and (21) in the paper) are needed that would lead to numerical analysis for the sideband instability. The insertion of (6.6) into (4.12) for $j=3$ leads to an eigenvalues equation for $\delta _\varOmega$, with the terms at $O(\epsilon _0^2\delta _y^2)$ and higher neglected. In particular, the eigenvalues equation for $\delta _\varOmega$ reads

(6.8)\begin{equation} (-\delta_\varOmega+\mathcal{L}_I)^2 + (-\delta_\varOmega+\mathcal{L}_I)(\alpha_{11}+\alpha_{22}) +\alpha_{12}\alpha_{21}= 0, \end{equation}

where $\mathcal {L}_r = \frac {1}{2} [\mathcal {L}(\mathrm {i} k_0\delta _x, \mathrm {i} k_0\delta _x) + \mathcal {L}(-\mathrm {i} k_0\delta _x, -\mathrm {i} k_0\delta _x)]$ and $\mathcal {L}_I = \frac {1}{2} [\mathcal {L}(\mathrm {i} k_0\delta _x, \mathrm {i} k_0\delta _y) -\mathcal {L} (-\mathrm {i} k_0\delta _x , -\mathrm {i} k_0\delta _x)]$ denote the contribution of the linear dispersion, and coefficients $\alpha _{ij}$ are given by

(6.9a)\begin{equation} \alpha_{11} =\frac{1}{4} V_{z,-}\epsilon_0^2 + (V_{x,+} + V_{z,+}) \left[ -\frac{1}{8}\delta_x +\left(\frac{1}{2}\mathcal{L}_2 +\frac{\mathcal{L}_2(1+\epsilon_0^2)\delta_{20} - \delta_x\mathcal{L}_2^2} {\delta_{20} -\mathcal{L}_2^2} \right)\right]\epsilon_0^2, \end{equation}

(6.9b)\begin{align} \alpha_{12} &= \alpha_1\epsilon_0^2 -\mathcal{L}_r -\frac{1}{4} V_{z,+}\epsilon_0^2 - (V_{x,-}-V_{z,-})\nonumber\\ &\quad \times \left[ \frac{1}{8}\delta_x -\frac{1}{2} \left(\mathcal{L}_2 +\frac{2\mathcal{L}_2(1+\epsilon_0^2)\delta_{20} - 2\delta_x(\mathcal{L}_2)^2} {\delta_{20} -(\mathcal{L}_2)^2} \right) \right]\epsilon_0^2, \end{align}

(6.9c)\begin{align} \alpha_{22} &=\left[ - \frac{1}{4} ( V_{x,-} -2V_{z,-} ) - \frac{1}{8}\partial_z( V_{x,-} -V_{z,-} )+ ( V_{x,-} -V_{z,-} ) \right.\nonumber\\ &\quad \left. \times \frac{\delta_x(1+\epsilon_0^2)\mathcal{L}_2-\delta_{20}\mathcal{L}_2^2} {\delta_{20} -\mathcal{L}_2^2} \right]\epsilon_0^2, \end{align}

(6.9d)\begin{align} \alpha_{21} &=\mathcal{L}_r -\alpha_1\epsilon_0^2 + \frac{1}{4} (V_{x,+} + 2V_{z,+})\epsilon_0^2 + \frac{1}{8}\partial_z(V_{x,+} + V_{z,+})\epsilon_0^2 \nonumber\\ &\quad - \epsilon_0^2(V_{x,+} + V_{z,+}) [\delta_x(1+\epsilon_0^2)\mathcal{L}_2-\delta_{20}\mathcal{L}_2^2]/ ({\delta_{20} -\mathcal{L}_2^2}), \end{align}

in which the $z$ dependent variables associated with $V_{i,\pm }(z)$ are evaluated for $z=0$, $\delta _{20} = \sqrt {\delta _x^2 + \delta _y^2}$ denotes the dimensionless magnitude of the sidebands away from the wavenumber of the Stokes wave, and the approximation $\mathcal {L}_2 \approx \frac {1}{2}(\mathcal {L}_I + \mathcal {L}_r)$ (cf. equations (11) and (15) by Crawford et al. (1981) for similar discussions) was used due to the terms associated with $\partial _t(\bar {\boldsymbol {V}}\boldsymbol {\cdot }\bar {\boldsymbol {V}}^*)$ and, therefore, $\boldsymbol {V}^{(20)}$, arising from the symmetry of the disturbances to the Stokes wave.

It is understood that the instability occurs if (6.8), a quadratic equation in $\delta _\varOmega$, has a positive imaginary root for $\delta _\varOmega$. Thereby, the instability occurs if the following inequality holds:

(6.10)\begin{equation} R_\varOmega \equiv \textstyle\frac{1}{4} (\alpha_{11}+\alpha_{22})^2 - \alpha_{12}\alpha_{21}<0, \end{equation}

in which $R_\varOmega$ is defined for later reference.

Figure 5 shows a comparison of the instability region predicted by (6.10) and its lower-order expression, equations (18–22) by Trulsen et al. (2000), and the exact computations by McLean (1982) for three different values of wave steepness. The equations of Trulsen et al. (2000) are presented in Appendix D for reference. In the cases chosen by McLean (1982) the value of a different steepness, defined as half the crest-to-trough height of Stokes waves and referred to as $\epsilon _M$ here, is given. For a direct comparison, the identity derived by Trulsen & Dysthe (1996) (i.e. equation (30) therein) was used to evaluate $\epsilon _0$, as follows:

(6.11)\begin{equation} \epsilon_0 = \epsilon_M - \textstyle\frac{1}{2} \epsilon_M^3. \end{equation}

Figure 5. A comparison of the instability region predicted by Trulsen et al. (2000), the exact computations based on McLean (1982) and the present theory based on (6.10) for three different values of wave steepness defined as $\epsilon _0=k_0a_0$: (a,d) $\epsilon _0 = 0.0995$, (b,e) $\epsilon _0=0.197$ and (c,f) $\epsilon _0 = 0.327$. (a–f) The exact results based on McLean (1982) (indicated by black dashed lines) are computed; (a–c) computations are based on equation (18) by Trulsen et al. (2000) (blue dots) and a leading-order approximation based on (4.12) and (4.13a) (red solid) – both approximations are correct to $O(\epsilon ^3_0\delta ^0)$ with the exact linear dispersion relation; (d–f) the approximations based on equations (19–22) by Trulsen et al. (2000) (blue dot–dashed), correct to $O(\epsilon _0^2\delta _x^2, \epsilon _0^2\delta ^0_y)$ and (6.10) derived in this paper (red dot–dashed).

It is seen from figure 5(a–c) that, compared with McLean (1982), both approximations behave similarly; good agreement is shown for small $\delta _x$ and small steepness $\epsilon _0$, whereas the agreement becomes less satisfactory for larger $\epsilon _0$. As approximations of higher-order accuracy are used, better agreement with McLean (1982) is clearly seen in figure 5(d–f) compared with figure 5(a–c). In particular, as seen in figure 5(d–f), equation (6.10) agrees with McLean (1982) in a nearly perfect manner for $\epsilon _0\approx 0.1$ and $\epsilon _0\approx 0.2$ in the range $\delta _x\lesssim 1$, while the predictions by (6.10) are also satisfactory for $\epsilon _0\approx 0.33$. Comparing the results by Trulsen et al. (2000) and (6.10) with McLean (1982) shown in figure 5(d–f), it is clear that (6.10) provides a better prediction for all three cases. In particular, (6.10) is capable of proving more accurate predictions for the instability of Stokes waves subject to sidebands in the transverse direction, as clearly shown in figure 5(e,f). Focusing on the instability due to transverse sidebands ($\delta _y$), figure 5 shows that equations (19–22) in Trulsen et al. (2000) fail to improve the accuracy of the lowest-order approximation (i.e. equation (18) by Trulsen et al. (2000)) for small $\delta _x$ for all three cases. The reason to this is explained as follows. Equations (19–22) of Trulsen et al. (2000), with the notation in this paper, lead to the instability region described by the inequality, as follows:

(6.12)\begin{equation} R_{T} \equiv \mathcal{L}_r \left[\mathcal{L}_r +\epsilon_0^2 \left(1- \frac{\delta_x^2}{\delta_{20} }\right) \right]+\frac{\epsilon^4_0\delta_x^2}{16}<0, \end{equation}

where $R_T$ is introduced for later reference. Equation (6.12) gives the instability region based on equation (18) by Trulsen et al. (2000) by setting $\delta _x$ to zero. For small $\epsilon _0$ ($\epsilon _0\lesssim 0.5$), the last term in $R_T$ is negligible, which gives two conditions for instability to occur,

(6.13)\begin{equation} \mathcal{L}_r <0\quad \text{or}\quad \mathcal{L}_r +\epsilon_0^2 \left(1- \frac{\delta_x^2}{\delta_{20} }\right)<0. \end{equation}

Noting that $\delta _x^2 /\delta _{20} <1$ for $\delta _x\lesssim 1$, the first condition in (6.13) dominates for small $\delta _x$. Using a leading-order approximation to $\mathcal {L}_r$, i.e. $\mathcal {L}_r\approx (2\delta _y^2-\delta ^2_x)/8$, the first condition is approximately

(6.14)\begin{equation} 2\delta_y^2 -\delta_ x^2 <0, \end{equation}

which corresponds to a (neutral) instability slope near the origin defined by $\delta _y/\delta _x = \pm 1/\sqrt {2}$ (corresponding to $\sim \pm 35.26^\circ$ to the carrier wave $\boldsymbol {k}_0$), as pointed out by Longuet-Higgins (1976) for the analysis of resonant energy transfer in a narrow spectrum and clearly shown by the predictions based on Trulsen et al. (2000) in figure 5 (see also § 6.3.2).

In contrast, (6.10) is capable of providing an accurate prediction of the instability region of Stokes waves subject to spanwise ($y$ direction) sideband perturbations for moderate values of wave steepness, $\epsilon _0\lesssim 0.2$, as can be clearly seen in figure 5(d,e). This is due to that (6.6) has taken into account the transverse sidebands $\delta _y$ in $\delta _{z,\pm }$ which denotes the effects on wave kinematics (i.e. $\bar {\boldsymbol {V}}$) in addition to the wave envelope $A$. For $\epsilon _0\approx 0.33$, the difference between (6.10) and McLean (1982) is likely due to that the terms described by $\mathcal {N}_{\approx ,{dir}}$ presented in § 4.2 at $O(\epsilon _0^2\delta _y^2)$ were neglected in the computations.

We now proceed to examining the energy growth rate of Stokes waves, which can be defined as the positive imaginary component of the roots of (6.8), as follows:

(6.15)\begin{equation} \text{Im}(\delta_\varOmega)=\text{Im}(\sqrt{R_\varOmega}), \end{equation}

where Im denotes the imaginary component. When instability occurs, the rate of amplitude growth/decline behaves like $\exp [\frac {1}{2}\text {Im}(\delta _\varOmega )\omega _0t]$ based on (6.5). This will be further examined in § 6.3.2.

Figure 6 shows the energy growth rate, $\textrm {Im}(\delta _\varOmega )$, varying with the sideband $\delta _x$ for five values of $\epsilon _0$ and varying with $\epsilon _0$ for two values of $\delta _x$. In addition to the experimental measurements and the predictions of (6.15), figure 6 also presents the results due to Crawford et al. (1981), for which small corrections following Krasitskii (1994) and Janssen (2004) (cf. § 4.11) were used. Crawford et al. (1981) is based on the Zakharov integral equation for the study of a uniform wave train and the results are correct to third order in wave steepness with all higher-order dispersion effects included. It is seen in figure 6 that (6.15) shows better agreement with Crawford et al. (1981) than Trulsen et al. (2000) (i.e. predicted by Im$\sqrt {R_T}$), which is more so for $0.15 \lesssim \epsilon _0$ and $0.4\lesssim \delta _x$.

Figure 6. The growth rate $\textrm {Im}(\delta _\varOmega )$ varying with the longitudinal bandwidth $\delta _x$ (a) and $\epsilon _0=k_0a_0$ (b) when the instability of Stokes waves occurs. In the figure, the predictions are computed based on Crawford et al. (1981), (6.15), and equations (19–22) by Trulsen et al. (2000) for three different values of steepness (a) and two different values of $\delta _x$ (b). Panel (b) shows, in addition, experimental measurements by Benjamin & Feir (1967) and Lake et al. (1977).

Figure 7 shows the instability boundary for growth of unstable perturbations for two-dimensional wave trains, predicted due to Benjamin & Feir (1967), Longuet-Higgins (1978), Crawford et al. (1981) and by this paper based on (6.15). Compared with the exact numerical results due to Longuet-Higgins (1976), (6.15) is in qualitative good agreement whereas the quantitative discrepancy is $\sim$23 % for the restabilization of steeper waves perturbed by long waves with $\delta _x \to 0$. Although the results based on (6.15) in figure 7 show slightly worse performance than the Zakharov equation, it can obviously yield restabilization for larger wave steepness, being an improvement over earlier versions of NLSEs.

Figure 7. Stability boundary for growth of unstable perturbations for two-dimensional wave trains predicted based on (6.15), compared with the results of Benjamin & Feir (1967), Longuet-Higgins (1978) and the Zakharov equation due to Crawford et al. (1981).

For large values of wave steepness, i.e. $0.2\lesssim \epsilon _0$, as shown in figures 6(a,b) and 7, the disagreement between (6.15) and Crawford et al. (1981) can be observed, which likely arises from that the additional time-dependant (small) disturbances are neglected in (6.15) to the linear velocity envelope (vector), $\bar {\boldsymbol {V}}$. The analysis in this section shows the results due to (6.15) are satisfactory.

Figures 5–7 suggest that the new NLSE described by (4.9) should work for $\epsilon _0\lesssim 0.25$ and $\delta _{x,y}> 1$ for the evolution of surface waves with directional spread.

6.3. Effects of wave directionality

The presented theory permits us to investigate the effects of wave directionality based on an order in wave steepness. Thereby, we focus the wave-directionality effects on waves at second order (§ 6.3.1) and on nonlinear energy transfer in a narrow wave spectrum by implementing the new NLSE numerically in § 6.3.2.

6.3.1. Waves at second order

As explained in § 2.4, one of the major effects of wave directionality arises from the fact that the superharmonic velocity head $H_v^{(22)}$ throughout the water columns becomes non-zero for directionally spread waves. This leads to both non-zero $\varPhi ^{(22)}$, additional contributions to elevation $\zeta ^{(22)}$ due to both $\varPhi ^{(22)}$ and $H_v^{(22)}$. These second-order effects on the spatial–temporal evolution of potential or elevation envelope of the first harmonic are considered negligible in most models based on an earlier version of NLSE, e.g. the NLSE by Dysthe (1979), Stiassnie (1984) and Trulsen et al. (2000). We focus on examining the velocity head $H_v^{(22)}$ and $H^{(20)}_v$ by studying the evolution of a directionally spread focused wave group at different times. For simplicity, third-order effects were neglected in the computations as they do not affect the discussions here in a qualitative manner.

Figures 8 and 9 show the dimensionless superharmonic and subharmonic velocity head defined in § 2.4, respectively. It is seen from figure 8 that, except for the group at linear focus (figure 8b,e,h), larger values of directional spread lead to smaller magnitudes in superharmonic velocity head, owing to the fact that, as the total wave energy is kept the same, the wave energy distributes in larger areas for larger values of directional spread. The magnitude of the superharmonic velocity head is shown to be $\lesssim 10 \sim 30\,\% \times \epsilon _0$ in figure 8. Similar observations can be found in figure 9 for the subharmonic velocity head which is of a similar order of magnitude to the superharmonic velocity head shown in figure 8. This indicates both components can be equally important for waves with directional spread.

Figure 8. The non-dimensional superharmonic velocity head $k_0(H_v^{(22)}+\text {c.c.})$ for the evolution of a focused wave group at different times for $z=0$. The parameters are the same as figure 4 except that $\theta _w$ differs: (a–c) $\theta _w = 15$ deg; (d–f) $\theta _w = 30$ deg; (g–i) $\theta _w = 45$ deg; (a,d,g) $t=-15 \times T_0$; (b,e,h) $t=0$ for the groups at linear focus with $\epsilon _0=0.3$; (c,f,i) $t=15\times T_0$ after focus.

Figure 9. The dimensionless velocity head $k_0H_v^{(20)}$ at different times for $z=0$. See the caption of figure 8 for the specific parameters.

6.3.2. Nonlinear energy transfer in a narrow wave spectrum

In order to examine the nonlinear energy transfer in a narrow wave spectrum, the temporal evolution of a ‘steep’ wave group is investigated numerically with the implementation of the new NLSE and equations (19–22) by Trulsen et al. (2000). Figure 10 shows the evolution of the discrete amplitude spectra, $\hat {A}(\boldsymbol {k},t)$, in the first quadrant of the Fourier space for a wave group at different times with $\epsilon _0 = 0.2$. It is seen in figure 10 that the results based on (4.9) and Trulsen et al. (2000) do not show differences in a qualitative manner and the energy conservation described by

(6.16)\begin{equation} \int^{\boldsymbol{x}_e}_{\boldsymbol{x}_s} |A(\boldsymbol{x},t)|^2\,\mathrm{d}\kern0.06em \boldsymbol{x} = \text{constant}, \end{equation}

were verified (not shown) to be within $0.01\,\%$ for both models. The dimensionless sidebands in the longitudinal and transverse direction are defined as, respectively,

(6.17a,b)\begin{equation} \delta_x = \frac{k_x-k_0}{k_0} \quad \text{and}\quad \delta_y = \frac{k_y}{k_0}, \end{equation}

in which $k_0$ was chosen to be the peak wavenumber of the initial wave spectrum, $\delta _x>0$ and $\delta _x<0$ denote the region of upper and lower sidebands in the longitudinal direction, respectively.

Figure 10. Discrete dimensionless amplitude spectra, $k_0|\hat {A}(\boldsymbol {k},t)|$, in the Fourier space at different times in the first quadrant for the spatial evolution of a wave group with $\epsilon _0=0.2$. The group starts to propagate at $t= -15 T_0$ and focuses linearly at $\boldsymbol {x}_{f} =0$ and $t_{f} = 0$. (a–c) The approximate NLSE with (4.13b) was used for computations; (d–f) (19–22) by Trulsen et al. (2000) were implemented. Contour levels (black solid line) are evenly distributed between 0.005 and 0.12 m in intervals of 0.005 m. (a) The red circles denote the wavenumbers chosen to investigate the growth rate of $|\hat {A}(\boldsymbol {k},t)|$ in figure 11. The longitudinal and transverse sidebands are defined as $\delta _x = (k_x-k_0)/k_0$ and $\delta _y= k_y/k_0$, respectively.

Compared with the initial time shown in figure 10(a,d), figure 10(b,e) shows broadened spectral amplitude in both the lower and higher sidebands in the longitudinal direction but narrowed in the region of larger sidebands in the transverse ($k_y$) direction, suggesting energy transfers from the latter to the former.

From $t=0$ to $t = 24\times T_0$ with $T_0$ denoting the wave period of the spectral peak, two main features can be observed in figure 10. Firstly, oblique energy transfers are clear and the change of the spectral amplitude in the regions of the lower ($k_x/k_0<1$ or $\delta _x<0$) and upper ($k_x/k_0>1$ or $\delta _x<0$) sidebands behaves differently. This indicates asymmetrical change rate of the spectral amplitude in sidebands. Particularly, in the region of the lower sidebands described by $0<\delta _x = (k_x-k_0)/k_0<-\sqrt {2}k_y/k_0$, a decline in the amplitude is indicated whereas an energy expansion is obvious in the region of upper sidebands, especially for $\delta _y>\delta _x/\sqrt {2}$. Secondly, relative to the original spectral peak at $k_x/k_0 =1$, a small downward shift of the spectral peak can be observed in figure 10(c,f), as is clearly shown by the nearest contour in the vicinity of $(1,0)$ that is no longer symmetrical relative to $k_x/k_0=1$.

In order to examine the change rate of the spectral amplitude in time, $|\hat {A}(\boldsymbol {k},t)|$ were calculated based on (4.9) at different times, and $|\hat {A}(\boldsymbol {k},t)|$ at 35 different wavenumbers were chosen and shown in figure 11. Figure 11 shows the non-dimensional spectral amplitude as a function of the dimensionless time $\frac {1}{2}\epsilon _0^2\omega _0t$. As is seen in figure 11, the change rate of $|\hat {A}(\boldsymbol {k},t)|$ behaves differently at different wavenumbers. The most interesting range is obviously for $\frac {1}{2}\epsilon _0^2\omega _0(t-t_0)$ in $\sim [1, 3]$ in which the amplitude changes rapidly. This range is in the vicinity of $t= t_{{f}}$ at which the wave group focuses linearly. We refer to $t\lesssim t_{{f}}$ and $t_{{f}}\lesssim t$ in this time range as the prefocus and postfocus phase, respectively.

Figure 11. Amplitude spectra $k_0|\hat {A}(\boldsymbol {k},t)|$ at 35 wavenumbers in the first quadrant of the Fourier $\boldsymbol {k}$ space as a function of the dimensionless time $\frac {1}{2}\epsilon _0^2\omega _0(t-t_0)$. The results are based on the new NLSE and the parameters were chosen the same as figure 10. In the figure, $\delta _x = (k_x-k_0)/k_0$ and $\delta _y=k_y/k_0$ denote the dimensionless sidebands in the longitudinal and transverse/spanwise direction, respectively. Short vertical lines near $\frac {1}{2}\epsilon _0^2\omega _0(t-t_0)= 2$ indicates for $t = t_{{f}}$ at which the wave group focuses linearly.

Focusing on the prefocus phase first, figure 11 shows the amplitude grows for the wavenumbers in the range $\delta _y\lesssim 0.1$ and wavenumbers in the upper sidebands with $\delta _x=0.36$ for $\delta _y\lesssim 0.35$. In contrast, the amplitudes for larger sidebands in the transverse direction, i.e. for $0.18\lesssim \delta _y$, experience a decline except for $\delta _x=0.36$. In the postfocus phase, the amplitudes that have experienced a growth (decline) in the prefocus phase tend to decrease (grow) exponentially. Especially, for the amplitudes at the wavenumbers in small transverse sidebands, i.e. $\delta _y\lesssim 0.09$, the rate of decline is larger than that of increase. The asymmetry in the rate of amplitude change in different regions is likely the cause for the downward shift of the spectral peak shown in figure 10(c,f).

The change rate of In$(k_0|\hat {A}|)$ in the range for $\frac {1}{2}\epsilon _0^2\omega _0(t-t_0)$ in $\sim [1, 3]$ in figure 11 shows that $\textrm {In}(k_0|\hat {A}|)$ decreases/grows with $\frac {1}{2}\epsilon _0^2\omega _0(t-t_0)$ (approximately) linearly with absolute slopes ${\lesssim }1$. This suggests that $|\hat {A}(\boldsymbol {k},t)|$ changes exponentially with the dimensionless time $\frac {1}{2}\epsilon _0^2\omega _0(t-t_0)$, unlikely arising from resonant energy transfers as resonance would lead to the change rate to be linear in time. The rapid change of the spectral amplitudes shown in figure 11 may be due to instability. This paper argues that it is likely due to a mechanism similar to oblique sideband instability of Stokes waves. To illustrate this, the amplitude change rate in the instability region of Stokes waves is shown in figure 12 for five different values of wave steepness. The results in figure 12 were computed based on (6.15) and equations (19–22) in Trulsen et al. (2000) (i.e. $\textrm {Im}\sqrt {R_T}$). As noted in § 6.2, it is understood that, when instability occurs, the amplitude behaves like

(6.18a,b)\begin{equation} |A|\sim a_0 \,\mathrm{e}^{\text{Im}(\delta_\varOmega )\omega_0(t-t_0)} \quad \text{and thus,}\ R_c \equiv \frac{\text{In}(k_0|A|)}{\frac{1}{2} \omega_0\epsilon_0^2 (t-t_0) } \sim \text{In} (\epsilon_0) \frac{\text{Im}(\delta_\varOmega )}{\frac{1}{2}\epsilon_0^2}, \end{equation}

where $R_c$ denotes the change rate of amplitude in the logarithmic scale, introduced for convenience. It also corresponds to the gradient of $\text {In} (k_0|A|)$ with respect to $\frac {1}{2} \omega _0\epsilon _0^2(t-t_0)$ shown in figure 11.

Figure 12. Instability growth rate of Stokes waves, defined as $\textrm {Im}(\delta _\varOmega )/(\epsilon _0^2/2)$, subject to oblique sideband perturbations for five values of wave steepness $\epsilon _0=k_0a_0$, obtained using the prediction by (6.8) (a–e) and equations (19–22), (f–j) by Trulsen et al. (2000) for five different values of steepness.

Based on (6.18b), it is clear that $R_c$ depends on both $\epsilon _0$ and the magnitude of $\text {Im}(\delta _\varOmega )$. Due to $\epsilon _0 < 0.5$ in general, which means $\text {In}(\epsilon _0)<0$, whether or not the amplitude grows in time depends on the sign of $\text {Im}(\delta _\varOmega )$; a positive sign leads to a decline whereas a negative sign a growth. In practice $0.1\lesssim \epsilon _0\lesssim 0.3$ would give a good estimate of the steepness range, yielding $-2.3\lesssim \text {In}(\epsilon _0) \lesssim -1.2$. Figure 12 shows $\delta _\varOmega /(\frac {1}{2}\epsilon _0^2)$ varying with transverse and longitudinal sidebands at five different values of $\epsilon _0$ varying from $0.1$ to $0.3$. The contours of large growth rates shown in figure 12 suggest that oblique energy transfers would be favoured when instability occurs, although the maximum $\delta _\varOmega$ always lies on the $\delta _x$ axis for deep-water Stokes waves. The magnitude of $\delta _\varOmega /(\frac {1}{2}\epsilon _0^2)$ predicted by both the present theory and Trulsen et al. (2000) is $\lesssim 1$. Obvious ‘coincidences’ between figures 11 and 12 are that both indicate the same magnitude for $R_c$, i.e. $\lesssim 1$, and that the wave amplitude in both cases changes exponentially. Due to the two points, this paper argues the instability of Stokes waves subject to sidebands disturbances in oblique directions would be a possible cause for the rapid change of wave amplitude shown in figure 11 in the neighbourhood of the linear focused time $t=t_{f}$.

Examining figures 10, 11 and 12 together, the following features are highlighted.

1. (i) The rapid energy transfers between upper and lower sidebands shown in figure 10 are likely due to the instability of Stokes waves in addition to the degenerate resonant interaction in the ‘figure of eight’ quartet resonant loop (Phillips 1967), as a result of that, the change rate of the amplitude shown in figure 11 is of the same order of magnitude as those in figure 12 when instability occurs. It is well understood that the instability of Stokes waves is mostly accompanied by the degenerate resonant interaction, as pointed out by Phillips (1967) where it is shown that the neutral instability modes (i.e. $\textrm {Im}(\delta _\varOmega )=0$) of Stokes waves can be directly derived from the degenerate resonant interaction (Hammack & Henderson 1993) and, most importantly, the instability of Stokes waves subject to oblique wave-train disturbances is likely to occur in a narrow spectrum.

2. (ii) A downward shift of the spectral peak arises probably from the asymmetrical change rates of energy between the different regions of sidebands. This shift is consistent with the finding of Trulsen & Dysthe (1997). When a downward shift of spectral peak is observed, rapid energy transfers are shown between lower and upper sidebands in the longitudinal direction, and between the transverse and longitudinal sidebands.

3. (iii) Oblique energy transfers towards large upper sidebands ($0.3\lesssim \delta _x$ or $0.27\lesssim \delta _y$) are obviously seen in figure 10(c,f), probably due to the different change rates within the instability region of Stokes waves; the contour levels of large change rates of amplitude have demonstrated that the change of energy favours more in oblique directions. This observation is consistent with Trulsen & Dysthe (1997) and the rapid energy transfer reported in Barratt et al. (2021). Barratt et al. (2021) focus on a numerical study of the propagation of a steep wave group based on direct numerical solutions of fully nonlinear potential flow equations. Barratt et al. (2021) reports an angle along $\text {ang}(\delta _y/\delta _x )=\pm 55^\circ$ to the spectral peak in the region of upper sidebands. This paper argues that the angle as such, $\pm 55^\circ$, may not exist. It is likely a coincidence of the effect shown in figure 11(d,e) that larger upper sidebands ($0.3\lesssim \delta _x$ or $0.27\lesssim \delta _y$) tend to grow exponentially at a rate among the largest. Moreover, Barratt et al. (2021) attribute the rapid oblique energy transfers, similar to those shown in the region (i.e. $\delta _y>\delta _x/\sqrt {2}$ for $0.3\lesssim \delta _x$) of the upper sidebands in figure 10(c,f), to non-degenerate resonant interaction. This argument is different from what is shown in figures 10, 11 and 12, i.e. that the amplitude changes exponentially is unlikely a result of resonant interactions.

7. Concluding remarks

Using a perturbation expansion, this paper has derived a new NLSE for the linear envelope evolution of three-dimensional surface gravity waves without the assumption of a narrow bandwidth. The new NLSE can be extended for more general cases, e.g. waves on water of finite uniform or slowly varying depth. Simple numerical implementations of the NLSE for the wave envelope using a pseudospectral, a split-step and a finite-difference method are demonstrated. The computational cost is of the same order as the numeral implementations of an existing NLSE, e.g. the Dysthe equation of Dysthe (1979) and the modified NLSE of Trulsen et al. (2000). With leading-order approximations to the NLSE, its capabilities for studying the instability region and energy growth rate of Stokes waves are demonstrated. The approximation is shown to be capable of providing satisfactory predictions of the instability region, compared with the exact computations by McLean (1982) for dimensionless wave steepness equal to as large as $0.33$. It can also provide good estimates of the energy growth rate and stability boundary of unstable perturbations for two-dimensional wave trains, compared with the results due to Crawford et al. (1981) and the exact numerical results based on Longuet-Higgins (1978). In these comparisons, the approximation has been demonstrated to provide a better performance than the modified NLSE by Trulsen et al. (2000), particularly for sideband disturbances of waves in an oblique direction to a train of Stokes waves.

This paper has also proposed a semianalytical approach for the description of wave fields at second order in wave steepness, based on a pseudospectral and a finite-difference method. It is demonstrated that the semianalytical approach is capable of providing exact predictions as that of Dalzell (1999). Compared with Dalzell (1999), the semianalytical approach can significantly improve the computational efficiency; the numerical implementations of Dalzell (1999) require computational operations at $O(N^2)$ and the semianalytical approach requires $O(N\,\text {In}(N))$ if $N$ Fourier modes are used. In the analysis of a velocity head at second order in wave steepness, it is suggested that the superharmonic and subharmonic waves may be equally important in directional sea states due to the same order of magnitude.

Through a study of the temporal evolution of a steep focused wave group, energy transfers in a narrow spectrum have been investigated. For the time range in the vicinity of the group at linear focus, rapid oblique energy transfers between different regions of sidebands are observed, likely arising from degenerate resonant interactions in the so-called ‘figure of eight’ quartet resonant loop (Phillips 1967) and the instability of Stokes waves subject to oblique sideband perturbations. The change rate of wave amplitude is found to be of the same order of magnitude as that of oblique instability of Stokes waves. A downward shift of the spectral peak is also observed in this paper, consistent with the finding in Trulsen & Dysthe (1997). This shift is observed in the course of rapid energy transfers between different sideband regions when oblique modulational instability occurs, in contrast to the temporal evolution of a unidirectional spectrum which is found to lead to a permanent shift of the spectrum peak only in a quasi-steady state in Dysthe et al. (2003). Rapid energy transfers towards the angle, ${\sim } \pm 55^\circ$, to the spectral peak in the region of upper sideband, that are reported in Barratt et al. (2021), are not supported by this study. Instead, the upper sideband range of large oblique directions (${>}45^\circ$) to the spectral peak is found to be favoured the most for the energy growth, coinciding with the sideband region where the instability growth rate of three-dimensional Stokes waves is among the largest.