3.1.1. Solution structure for a finite uniform depth

We proceed to explain how the velocity potential and vertical velocity can be expressed in the form of functions of both $\zeta$ and $\psi$. Following Dommermuth & Yue (Reference Dommermuth and Yue1987) and West et al. (Reference West, Brueckner, Janda, Milder and Milton1987), the HOS method proposes to letting

(3.1)\begin{equation} \varPhi^{(1)}(\boldsymbol{x},0,t) = \psi(\boldsymbol{x},t) \equiv \varPhi(\boldsymbol{x},\zeta(\boldsymbol{x},t),t), \end{equation}

which, due to the perturbation expansion (2.6), leads to the velocity potential on the free water surface given by

(3.2)\begin{equation} \varPhi(\boldsymbol{x},\zeta,t) = \sum_{m=1}^{M} \varPhi^{(m)}(\boldsymbol{x},z,t) \quad \text{for}\ {z=\zeta}. \end{equation}

An expression for $\varPhi ^{(m)}$ for $m>1$ can be obtained through the subsequent procedures: Taylor expanding the terms on the right-hand side of (3.2) about $z=0$, inserting (3.1) and collecting the same orders in wave steepness. Hence, $\varPhi ^{(m)}$ is expressed as functions of the lower-order parameters as

(3.3)\begin{equation} \varPhi^{(m)}_0(\boldsymbol{x},t) =- \sum_{k=1}^{m-1} \dfrac{1}{k!}\zeta^k\partial_z^k \varPhi^{(m-k)}(\boldsymbol{x},z,t)\quad\text{for}~z=0\quad \text{and}\quad m \in \{2,3,\ldots M\}, \end{equation}

where, as noted, the subscript ‘0’ denotes the parameters evaluated at a still water surface, $z=0$; i.e. $\varPhi ^{(m)}_0 \equiv \varPhi ^{(m)}(\boldsymbol {x},0,t)$. The expression (3.3) suggests that if $\psi$ and $\zeta$ are given, $\varPhi ^{(m)}$ for $m>1$ will be obtained in sequence from the lowest to higher orders. Similarly, the vertical velocity ($W$) on the unknown free water surface can be obtained from Taylor expanding its definition given by (2.4b) about $z=0$ and inserting the expression for $\varPhi ^{(m)}_0$ to give

(3.4a)$$\begin{gather} W(\boldsymbol{x},t) = \sum_{m=1}^{M} \epsilon^m W^{(m)}(\boldsymbol{x},t), \quad \text{with} \end{gather}$$

(3.4b)$$\begin{gather}W^{(m)}(\boldsymbol{x},t) = \sum_{k=0}^{m-1} \dfrac{\zeta^k}{k!}\partial_z^{k+1} \varPhi^{(m-k)}(\boldsymbol{x},z,t)\quad \text{for}\ z=0\quad \text{and}\quad m\in \{1,2,\ldots,M \}. \end{gather}$$

Therefore, $W^{(m)}$ can be obtained in sequence from the lowest order $m=1$ to higher orders due to (2.7) for $\varPhi ^{(m)}(\boldsymbol {x},z,t)$, (3.3) for $\varPhi ^{(m)}_0(\boldsymbol {x},t)$ and (3.4), with $\psi$ and $\zeta$ as input.

3.1.2. The $M$th order accurate equations in the HOS method

Inserting the perturbed solution (3.1) for $\psi$ and (3.4a,b) for $W$ into the boundary conditions (2.5a,b) and keeping the terms up to order $M$ gives rise to

(3.5a,b)\begin{equation} \partial_t \zeta - \partial_z\varPhi^{(1)} = \mathcal{W}_{M}(\boldsymbol{x},t) \quad \text{and}\quad \partial_t\psi + g\zeta = \mathcal{T}_{M}(\boldsymbol{x},t), \quad \text{for}\ z=0, \end{equation}

where the subscript ‘$M$’ denotes the truncated order of accuracy in wave steepness,

(3.6a,b)\begin{equation} \mathcal{W}_{M}(\boldsymbol{x},t) = \sum_{m=1}^M \mathcal{W}^{(m)}(\boldsymbol{x},t)\quad \text{and}\quad \mathcal{T}_{M}(\boldsymbol{x},t) = \sum_{m=1}^M \mathcal{T}^{(m)}(\boldsymbol{x},t), \end{equation}

where $\mathcal {W}_{1}\equiv \mathcal {W}^{(1)} =0$ and $\mathcal {T}_{1}\equiv \mathcal {W}^{(1)}=0$, and the nonlinear forcing terms, $\mathcal {W}^{(m)}(\boldsymbol {x},t)$ with $m\geq 2$, are non-vanishing and given by

(3.7a)$$\begin{gather} \mathcal{W}^{(2)}= W^{(2)}- \boldsymbol{\nabla}\psi\boldsymbol{\cdot}\boldsymbol{\nabla}\zeta, \end{gather}$$

(3.7b)$$\begin{gather}\mathcal{W}^{(3)} = W^{(3)}+ W^{(1)} (\boldsymbol{\nabla}\zeta)^2, \end{gather}$$

(3.7c)$$\begin{gather}\mathcal{W}^{(4)}= W^{(4)} + W^{(2)}(\boldsymbol{\nabla}\zeta)^2\quad \text{and} \end{gather}$$

(3.7d)$$\begin{gather}\mathcal{W}^{(5)} = W^{(5)} + W^{(3)}(\boldsymbol{\nabla}\zeta)^2, \end{gather}$$

which are explicitly expressed up to the fifth order in wave steepness; similarly, $\mathcal {T}^{(m)}(\boldsymbol {x},t)$ are

(3.8a)$$\begin{gather} \mathcal{T}^{(2)} = \tfrac{1}{2} ( W^{(1)})^2-\tfrac{1}{2} (\boldsymbol{\nabla}\psi)^2, \end{gather}$$

(3.8b)$$\begin{gather}\mathcal{T}^{(3)} = W^{(1)} W^{(2)} , \end{gather}$$

(3.8c)$$\begin{gather}\mathcal{T}^{(4)}= \tfrac{1}{2} (W^{(2)})^2+ W^{(1)} W^{(3)}+ \tfrac{1}{2} (W^{(1)})^2(\boldsymbol{\nabla}\zeta)^2, \end{gather}$$

(3.8d)$$\begin{gather}\mathcal{T}^{(5)} = W^{(2)} W^{(3)} + W^{(1)} W^{(4)} + (W^{(1)} W^{(2)})(\boldsymbol{\nabla}\zeta)^2. \end{gather}$$

We can summarize that the HOS method has derived the $M$th order accurate equations (3.5a,b) that can be numerically solved for $\zeta$ and $\psi$, with $\varPhi ^{(m)}$ and $W^{(m)}$ obtained from (3.3) and (3.4), respectively, in sequence from the lowest order to higher orders (West et al. Reference West, Brueckner, Janda, Milder and Milton1987). The nonlinear terms on the right-hand side of (3.5a,b) are given by (3.7) and (3.8).

4.3.1. Methodology illustration

We proceed to explain the fundamental methodology of the new envelope framework using the velocity potential in the first- and second-order approximations as examples. For later reference, we define

(4.15a,b)\begin{equation} B_0(\boldsymbol{x},t) = B_s(\boldsymbol{x},t)\quad \text{and}\quad B^{(11)}(\boldsymbol{x},z,t) = B(\boldsymbol{x},z,t). \end{equation}

Due to the definitions $\varPhi _0^{(1)}=\psi$ and also (2.12b) for the potential on the free water surface, we propose to let $B(\boldsymbol {x},z,t)=[\varPhi ^{(1)}(\boldsymbol {x},z,t) ]^{[1]}_+$ and, therefore,

(4.16a)\begin{equation} \varPhi^{(1)}(\boldsymbol{x},z,t) = \dfrac{1}{2} B(\boldsymbol{x},z,t)\exp\left({\mathrm{i} (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \text{c.c.}, \end{equation}

where the relation $\hat {B}_0(\boldsymbol {k},t) = \hat {B}_s(\boldsymbol {k},t)$ holds by definition, and the Laplace equation and seabed boundary condition require

(4.16b)\begin{equation} B(\boldsymbol{x},z,t) = \int_{-\infty}^{\infty} \hat{B}_0(\boldsymbol{k},t) \dfrac{\cosh |\boldsymbol{k}+\alpha \boldsymbol{k}_0|(z+h)}{\cosh |\boldsymbol{k}+\alpha\boldsymbol{k}_0|h} \mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}} \,\mathrm{d}\boldsymbol{k}. \end{equation}

The expression (4.16a) for $\varPhi ^{(1)}$ is in a new form we intend to obtain at the first order in wave steepness. We next proceed to $m=2$ for the velocity potential $\varPhi ^{(2)}$. Inserting (2.12a,b) and (4.16a) into (3.3) for $\varPhi ^{(2)}$, we obtain

(4.17a,b)\begin{align} \varPhi_0^{(2)} = \varPhi^{(20)}_0(\boldsymbol{x},t) + \varPhi_0^{(22)}(\boldsymbol{x},t)\quad \text{and, thus,}\quad \varPhi^{(2)} = \varPhi^{(20)}(\boldsymbol{x},z,t) + \varPhi^{(22)}(\boldsymbol{x},z,t), \end{align}

where the superscript in a form of ‘($mj$)’ denotes $O(\epsilon ^m)$ and the $j$th wave harmonic. Especially the superscripts ‘$20$’ and ‘$22$’ denote the potential for the second-order subharmonic and superharmonic bound waves, respectively, which are obtained through the separation of harmonics given by

(4.18a,b)\begin{align} {\varPhi}^{(20)}_0=-\tfrac{1}{4} A^*\partial_zB+\text{c.c.} \quad \text{and}\quad \varPhi^{(22)}_0 =- \tfrac{1}{4} A\partial_zB\mathrm{e}^{2\mathrm{i} (\alpha\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)} +\text{c.c.}, \quad \text{for}\ z=0. \end{align}

We next define the envelope of the two potentials as

(4.19a,b)\begin{align} B^{(2j)}(\boldsymbol{x},z,t) =\left[ \varPhi^{(2j)}\right]^{[j]}_{+} \quad \text{and}\quad B^{(2j)}_0(\boldsymbol{x},t) = \left[ \varPhi^{(2j)}_0\right]^{[j]}_{+},\quad \text{for}\ j=0 \quad \text{and}\quad j=2, \end{align}

which show the relations between the envelope and potential due to the second-order superharmonic ($\,j=2$) and subharmonic ($\,j=0$) waves. As the second-order subharmonic envelopes based on (4.19a,b) and (4.18a) depend only on slowly varying envelopes $A$ and $B$, they are used in the new framework. Nevertheless, the envelopes of the waves with a second harmonic given by (4.19a,b) and (4.18b) have the same temporal and spatial variation as the second-order velocity potential used in the HOS method, and thereby do not introduce merits in the efficiency in numerical implementations compared with using the HOS method. To make a difference, the second-order superharmonic envelope in the form

(4.20)\begin{equation} B^{(22)}_0 =-\tfrac{1}{2} A\partial_zB \quad \text{for}\ z=0 \end{equation}

is used instead, which depends only on the slowly varying envelopes. Due to the definition of the second-order envelopes $B^{(2j)}$, the second-order potentials can also be given by

(4.21)\begin{equation} \varPhi^{(2j)} = \tfrac{1}{2} B^{(2j)}(\boldsymbol{x},z,t)\exp({\mathrm{i} j(\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t)}) +\text{c.c.}, \quad \text{for}\ j=0\quad \text{and}\quad j=2. \end{equation}

The Laplace equation for $\varPhi ^{(2j)}$ and the seabed boundary condition lead to

(4.22)\begin{equation} B^{(2j)}(\boldsymbol{x},z,t) = \int_{-\infty}^{\infty} \hat{B}^{(2j)}_0(\boldsymbol{k},t) \dfrac{\cosh |\boldsymbol{k}+j\alpha\boldsymbol{k}_0|(z+h)}{\cosh |\boldsymbol{k}+j\alpha\boldsymbol{k}_0|h} \mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}} \,\mathrm{d}\boldsymbol{k}, \end{equation}

where $j=0$ and $j=2$ for the subharmonic and superharmonic envelopes, respectively; as noted, the subscript ‘0’ denotes the evaluation at $z=0$ and the hat added denotes the Fourier transform.

4.3.2. Velocity potential and vertical velocity

At second order, we have obtained an equating form for the second-order potential that is in a form of the linear superposition of different wave harmonics and envelopes. All second-order envelopes are functions of slowly varying envelopes in the lower order in wave steepness. Following the methodology presented in § 4.3.1, the envelope of an individual field can be obtained to arbitrary order in wave steepness, which has been given explicitly here up to the fourth order and the general expressions up to arbitrary order are derived in Appendix A. In particular, we propose to obtain a new expression for the potential at different orders in wave steepness based on the Laplace equation for an individual potential, the seabed condition, and the perturbation expansion (3.3). They have a general form as follows:

(4.23a,b)\begin{equation} \varPhi^{(m)} = \sum_{j=0}^{j=m} \varPhi^{(mj)}(\boldsymbol{x},z,t),\quad \text{with}\ \varPhi^{(mj)} \equiv \dfrac{1}{2} B^{(mj)}(\boldsymbol{x},z,t)\exp\left({\mathrm{i} j (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \text{c.c.}. \end{equation}

Therefore,

(4.24a)\begin{align} \varPhi^{(m)}_0 &= \sum_{j=0}^{j=m} \varPhi^{(mj)}_0(\boldsymbol{x},t),\quad \text{with} \end{align}

(4.24b)\begin{align} \varPhi^{(mj)}_0 &\equiv \tfrac{1}{2} B^{(mj)}_0(\boldsymbol{x},t)\exp\left({\mathrm{i} j (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \text{c.c.} \end{align}

(4.24c)\begin{align} &\equiv \tfrac{1}{2} \bar{\varPhi}^{(mj)}_0(\boldsymbol{x},t)\exp\left({\mathrm{i} j (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \text{c.c.}, \end{align}

where the velocity potential of the $j$th harmonic in the $m$th order in wave steepness in the form of (4.24b) is used in the new framework. The differences between the envelope $B_0^{(mj)}$ and $\bar {\varPhi }_0^{(mj)}$ lie in the fact that the latter is obtained based on (3.3) that corresponds to the factor in front of the $j$th harmonic due to $\exp [\mathrm {i} j(\alpha \boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\beta \omega _0t)]$, and thereby $B_0^{(mj)}=\bar {\varPhi }_0^{(mj)}$ only if $m=j$. The envelopes are

(4.25a,b)\begin{equation} B^{(mj)}(\boldsymbol{x},z,t) = \left(\varPhi^{(mj)}\right)_+^{[j]} \quad \text{and}\quad B_0^{(mj)}(\boldsymbol{x},t) = \left(\varPhi^{(mj)}_0\right)_+^{[j]}, \end{equation}

which hold by definition. Due to the Laplace equation and the seabed condition, we arrive at

(4.26)\begin{equation} B^{(mj)}(\boldsymbol{x},z,t) = \int_{-\infty}^{\infty} \hat{B}^{(mj)}_0(\boldsymbol{k},t) \dfrac{\cosh |\boldsymbol{k}+j\alpha\boldsymbol{k}_0|(z+h)}{\cosh |\boldsymbol{k}+j\alpha\boldsymbol{k}_0|h} \mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}} \,\mathrm{d}\boldsymbol{k}. \end{equation}

Thereby, if $B_0^{(mj)}(\boldsymbol {x},t)$ is given, $\varPhi ^{(m)}$ and $\varPhi ^{(mj)}$ will be explicitly obtained from (4.23a,b), respectively. Next we only have to explain how to obtain the envelope on the still water surface, $B^{(mj)}_0(\boldsymbol {x},t)$, and their Fourier transform having appeared in the integrand of the integral given by (4.26) in practice, following the same procedures as for $m=2$ presented in § 4.3.1. Inserting (2.12a) and (2.12b) for the surface elevation and the potential on the free water surface, respectively, into (3.3), we obtain in sequence up to the fourth order in wave steepness

(4.27a)$$\begin{gather} \varPhi^{(2)}_0(\boldsymbol{x},t) = \varPhi^{(22)}_0(\boldsymbol{x},t) + \varPhi^{(20)}_0(\boldsymbol{x},t) , \end{gather}$$

(4.27b)$$\begin{gather}\varPhi^{(3)}_0(\boldsymbol{x},t) = \varPhi^{(31)}_0 (\boldsymbol{x},t) + \varPhi^{(33)}_0 (\boldsymbol{x},t), \end{gather}$$

(4.27c)$$\begin{gather}\varPhi^{(4)}_0 (\boldsymbol{x},t) = \varPhi^{(40)}_0 (\boldsymbol{x},t) + \varPhi^{(42)}_0 (\boldsymbol{x},t) + \varPhi_0 ^{(44)}(\boldsymbol{x},t) , \end{gather}$$

where $\varPhi ^{(mj)}_0$ is the (real) potential of the $j$th harmonic at $O(\epsilon ^m)$ and

(4.28)\begin{equation} {\varPhi}_0^{(mj)} = 0,\quad \text{for}\ mj\in \{21,30,32,41,43 \}. \end{equation}

The new framework aims to express the non-vanishing potentials $\varPhi ^{(mj)}_0$ in the form of (4.24b), relying on the middle step for $\varPhi ^{(mj)}_0$ given by (4.24c) that depends on the explicit expression for $\bar {\varPhi }_0^{(mj)}$. Thereby, as noted, $\bar {\varPhi }_0^{(mj)}$ are obtained from (3.3) through collecting the $j$th harmonics at an individual order in wave steepness from the lowest to higher orders in sequence; explicitly, for $z=0$,

(4.29a)$$\begin{gather} \bar{\varPhi}^{(20)}_0=- \tfrac{1}{2} A^*\partial_zB, \end{gather}$$

(4.29b)$$\begin{gather}\bar{\varPhi}_0^{(22)} =-\tfrac{1}{2} A\partial_zB, \end{gather}$$

(4.29c)$$\begin{gather}\bar{\varPhi}^{(31)}_0=- \tfrac{1}{2}\left( \partial_zB^{(22)} A^* + 2\mathcal{R}(\partial_zB^{(20)}) A + \tfrac{1}{2} |A|^2 \partial_{zz}B+\tfrac{1}{4} A^2\partial_{zz} B ^* \right), \end{gather}$$

(4.29d)$$\begin{gather}\bar{\varPhi}_0^{(33)}=-\tfrac{1}{8} A^2 \partial_{zz}B- \tfrac{1}{2} A\partial_z B^{(22)}, \end{gather}$$

(4.29e)$$\begin{gather}\bar{\varPhi}^{(40)}_0=- \tfrac{1}{2} \left( \partial_zB^{(31)} A^* + {\tfrac{1}{2}}\left(\partial_{zz}B^{(20)}\right) |A|^2 + \tfrac{1}{4}\partial_{zz}B^{(22)} (A^2)^* + \tfrac{1}{8} |A|^2A^*\partial_{zzz}B \right),\\ \bar{\varPhi}^{(42)}_0=- \tfrac{1}{2}\partial_zB^{(31)} A - \tfrac{1}{2} \partial_zB^{(33)} A^* - \tfrac{1}{4}\partial_{zz}B^{(22)} |A|^2- \tfrac{1}{4}\mathcal{R}\left(\partial_{zz}B^{(20)}\right) A^2 \nonumber\\ - \tfrac{1}{16} |A|^2A\partial_{zzz}B - \tfrac{1}{48}A^3\partial_{zzz}B^*, \notag \end{gather}$$

(4.29f)$$\begin{gather}\bar{\varPhi}_0^{(44)} =-\tfrac{1}{48} A^3 \partial_{zzz}B-\tfrac{1}{8} A^2\partial_{zz} B^{(22)} -\tfrac{1}{2} A^2\partial_{z}B^{(33)}, \end{gather}$$

where it is clear that $\bar {\varPhi }^{(mj)}_0$ depends only on the slowly varying envelopes. The envelopes of the velocity potential $B^{(mj)}$ rely on their values at the still water surface, $B^{(mj)}_0$, due to their explicit form given by (4.25). To this end, the envelopes at the still water surface $B^{(mj)}_0$ are obtained through their relation with $\bar {\varPhi }_0^{(mj)}$ due to (4.24b,c). As noted, $B_0^{(mj)}=\bar {\varPhi }_0^{(mj)}$ for $m=j$ leads to

(4.30a)$$\begin{gather} B_0^{(22)} = \bar{\varPhi}_0^{(22)}\equiv-\tfrac{1}{2} A\partial_zB, \end{gather}$$

(4.30b)$$\begin{gather}B_0^{(33)}= \bar{\varPhi}_0^{(33)}\equiv-\tfrac{1}{8} A^2 \partial_{zz}B- \tfrac{1}{2} A\partial_z B^{(22)} , \end{gather}$$

(4.30c)$$\begin{gather}B_0^{(44)} =\bar{\varPhi}_0^{(44)}\equiv-\tfrac{1}{48} A^3 \partial_{zzz}B-\tfrac{1}{8} A^2\partial_{zz} B^{(22)} -\tfrac{1}{2} A^2\partial_{z}B^{(33)}. \end{gather}$$

The other non-vanishing $B_0^{(mj)}$ are obtained through combining the relation with $\bar {\varPhi }_0^{(mj)}$ given by (4.24b,c) and their definitions by (4.25b). Using in addition the properties of Fourier transforms, they can be especially obtained through an inverse Fourier transform as follows:

(4.31)$$\begin{gather} \hat{B}^{(mj)}_0(\boldsymbol{k}+j\alpha\boldsymbol{k}_0,t)\exp\left({-\mathrm{i} j\beta\omega_0t}\right) = \varTheta[(\boldsymbol{k}+j\alpha\boldsymbol{k}_0)\boldsymbol{\cdot}\boldsymbol{k}_0] \left\{ \hat{\bar{\varPhi}}^{(mj)}(\boldsymbol{k}+j \alpha\boldsymbol{k}_0,t)\exp\left({-\mathrm{i} j\beta\omega_0t}\right) \right.\nonumber\\ +\left.\left[\hat{\bar{\varPhi}}^{(mj)}(-\boldsymbol{k}-j\alpha\boldsymbol{k}_0,t)\exp\left({-\mathrm{i} j\beta\omega_0t}\right) \right]^* \right\}. \end{gather}$$

Here $\hat {\bar {\varPhi }}_0^{(mj)}(\boldsymbol {k},t)$ is the Fourier transform of $\bar {\varPhi }_0^{(mj)}(\boldsymbol {x},t)$.

It should be noted that the use of (4.30a–c) and (4.31) for the envelopes, $B^{(mm)}$ and $\hat {B}^{(mj)}$ with $m\neq j$, respectively, contribute to the improvement of the computational efficiency, compared with using the original definition of the envelopes given by (4.25b). It is by virtue of the fact that $\bar {\varPhi }^{(mj)}$ always have the same spatial (long) scale as envelope $A$ and $B_s$. However, due to the linear translation operator indicated by the independent variable $\boldsymbol {k}+j\alpha \boldsymbol {k}_0$ in (4.31), great care is needed in the numerical implementation when using a (inverse) fast Fourier transform (FFT).

It is now understood that potential $\varPhi _0^{(m)}$ at a nonlinear order in wave steepness admits three equating forms, which are in a form given by (3.3) and

(4.32)\begin{align} \varPhi^{(m)}_0 &\equiv \sum_{j=0}^{m} \left[ \tfrac{1}{2} B^{(mj)}_0 \exp\left({\mathrm{i} j (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \text{c.c.} \right]\nonumber\\ &\equiv \sum_{j=0}^{m} \left[ \tfrac{1}{2} \bar{\varPhi}^{(mj)}_0 \exp\left({\mathrm{i} j (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \text{c.c.} \right]. \end{align}

Inserting the expression of $\varPhi ^{(m)}$ and $\varPhi _0^{(m)}$ given by (4.23a) and (4.24a), respectively, into (3.4) leads to the vertical velocity given by

(4.33a)\begin{equation} w^{(1)} = \tfrac{1}{2} \partial_zB\exp\left({\mathrm{i} (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) +\text{c.c.},\quad W^{(1)} = w^{(1)}_0, \end{equation}

and

(4.33b)$$\begin{gather} W^{(m)} = \sum_{j=0}^{j=m} \left[\tfrac{1}{2}\bar{W}^{(mj)}(\boldsymbol{x},t) \exp\left({\mathrm{i} j (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \text{c.c.} \right]\quad \text{for}\ m=2,3,\ldots \end{gather}$$

(4.33c)$$\begin{gather}w^{(m)} = \sum_{j=0}^{j=m} \tfrac{1}{2} \partial_z{B}^{(mj)}(\boldsymbol{x},z,t) \exp\left({\mathrm{i} j (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \text{c.c.} \quad \text{for}\ m=2,3,\ldots \end{gather}$$

where the non-vanishing terms are expressed as

(4.33d)$$\begin{gather} \bar{W}^{(20)} = \partial_zB^{(20)} + \tfrac{1}{2} A^*\partial_{zz}B, \end{gather}$$

(4.33e)$$\begin{gather}\bar{W}^{(22)} = \partial_zB^{(22)} + \tfrac{1}{2} A\partial_{zz}B, \end{gather}$$

(4.33f)$$\begin{gather}\bar{W}^{(31)} = \partial_zB^{(31)} + \tfrac{1}{2} A^*\partial_{zz}B^{(22)} + A{\mathcal{R}\left(\partial_{zz}B^{(20)}\right)} + \tfrac{1}{8} A^2\partial_{zzz}B^* + \tfrac{1}{4} |A|^2\partial_{zzz}B , \end{gather}$$

(4.33g)$$\begin{gather}\bar{W}^{(33)} = \partial_zB^{(33)} + \tfrac{1}{2} A\partial_{zz}B^{(22)} + \tfrac{1}{8} A^2\partial_{zzz}B, \end{gather}$$

(4.33h)$$\begin{gather}\bar{W}^{(40)} = \partial_zB^{(40)} + \tfrac{1}{2} A^*\partial_{zz}B^{(31)} + \tfrac{1}{{4}} |A|^2 \left(\partial_{zzz}B^{(20)}\right) +\tfrac{1}{8}(A^*)^2\partial_{zzz}B^{(22)}\nonumber\\ + \tfrac{1}{16} |A|^2 A^*\partial_{zzzz}B,\end{gather}$$

(4.33i)$$\begin{gather}\bar{W}^{(42)} = \partial_zB^{(42)} + \tfrac{1}{2} \partial_{zz}B^{(31)} A +\tfrac{1}{2} \partial_{zz}B^{(33)} A^* + \tfrac{1}{4}|A|^2\partial_{zzz}B^{(22)} + \tfrac{1}{4}A^2\mathcal{R}\left(\partial_{zzz}B^{(20)}\right)\nonumber\\ + \tfrac{1}{16} |A|^2A\partial_{zzzz}B+ \tfrac{1}{{48}} A^3\partial_{zzzz}B^*, \end{gather}$$

(4.33j)$$\begin{gather}\bar{W}^{(44)} = \partial_zB^{(44)} + \tfrac{1}{2} A\partial_{zz}B^{(33)} + \tfrac{1}{8} A^2\partial_{zzz}B^{(22)} + \tfrac{1}{48}A^3\partial_{zzzz}B, \end{gather}$$

for $z=0$, which depend only on the slowly varying envelopes.

4.3.3. The nonlinear forcing terms on a still water surface

The forcing terms at different orders in wave steepness, described by (3.7a–e) and (3.8a–e), can also be expressed in the form of a function of envelope $A$ and $B_s$ (bearing in mind that $B_s=B_0$). Based on (3.7a–e) and (3.8a–e), we obtain

(4.34a)$$\begin{gather} \mathcal{W}^{(m)}(\boldsymbol{x},t) = \sum_{j=0}^{j=m} \underbrace{ \left[ \tfrac{1}{2} \bar{\mathcal{W}}^{(mj)}(\boldsymbol{x},t) \exp\left({\mathrm{i} j (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) +\text{c.c.} \right] }_{{\equiv} \mathcal{W}^{(mj)}(\boldsymbol{x},t)} \quad \text{and} \end{gather}$$

(4.34b)$$\begin{gather}\mathcal{T}^{(m)}(\boldsymbol{x},t) = \sum_{j=0}^{j=m} \underbrace{ \left[ \tfrac{1}{2} \bar{\mathcal{T}}^{(mj)}(\boldsymbol{x},t) \exp\left({\mathrm{i} j (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) +\text{c.c.} \right] }_{{\equiv} \mathcal{T}^{(mj)}(\boldsymbol{x},t)}, \end{gather}$$

where both $\mathcal {W}^{(mj)}(\boldsymbol {x},t)$ and $\mathcal {T}^{(mj)}(\boldsymbol {x},t)$ are introduced by definition and they are real functions, with their non-zero complex envelopes $\bar {\mathcal {W}}^{(mj)}$ given by

(4.35a)$$\begin{gather} \bar{\mathcal{W}}^{(20)} = \bar{W}^{(20)}- \tfrac{1}{2} (\boldsymbol{\nabla} B_s+\mathrm{i}\boldsymbol{k}_0B_s)^*\boldsymbol{\cdot}(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)A, \end{gather}$$

(4.35b)$$\begin{gather}\bar{\mathcal{W}}^{(22)} = \bar{W}^{(22)}- \tfrac{1}{2}(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)B_s\boldsymbol{\cdot}(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)A, \end{gather}$$

(4.35c)$$\begin{gather}\bar{\mathcal{W}}^{(31)} = \bar{W}^{(31)} + \tfrac{1}{4}\partial_zB^* \left[(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)A\right]^2 + \tfrac{1}{2} \partial_zB \left|(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)A\right|^2, \end{gather}$$

(4.35d)$$\begin{gather}\bar{\mathcal{W}}^{(33)} = \bar{W}^{(33)}+\tfrac{1}{4}\partial_zB \left[(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)A\right]^2, \end{gather}$$

(4.35e)$$\begin{gather}\bar{\mathcal{W}}^{(40)} = \bar{W}^{(40)} + \tfrac{1}{2} \bar{W}^{(20)} \left|(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)A\right|^2 + \tfrac{1}{4} \bar{W}^{(22)} \left[(\boldsymbol{\nabla} A+\mathrm{i}\boldsymbol{k}_0 A)^*\right]^2 , \end{gather}$$

(4.35f)$$\begin{gather}\bar{\mathcal{W}}^{(42)} = \bar{W}^{(42)} + \tfrac{1}{2} \mathcal{R}\left[\bar{W}^{(20)} \right](\boldsymbol{\nabla} A+\mathrm{i}\boldsymbol{k}_0 A)^2 + \tfrac{1}{2} \bar{W}^{(22)} \left|(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)A\right|^2 , \end{gather}$$

(4.35g)$$\begin{gather}\bar{\mathcal{W}}^{(44)} = \bar{W}^{(44)} + \tfrac{1}{4} \bar{W}^{(22)} (\boldsymbol{\nabla} A+\mathrm{i}\boldsymbol{k}_0 A)^2, \end{gather}$$

for $z=0$, and the non-zero $\bar {\mathcal {T}}^{(mj)}$ given by

(4.36a)$$\begin{gather} \bar{\mathcal{T}}^{(20)} =-\tfrac{1}{4} |(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)B|^2 +\tfrac{1}{4} |\partial_zB|^2, \end{gather}$$

(4.36b)$$\begin{gather}\bar{\mathcal{T}}^{(22)} =-\tfrac{1}{4}\left[(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)B\right]^2 +\tfrac{1}{4} (\partial_zB)^2, \end{gather}$$

(4.36c)$$\begin{gather}\bar{\mathcal{T}}^{(31)} = \tfrac{1}{2} \bar{W}^{(22)}\partial_zB^* + \mathcal{R}\left(W^{(20)} \right)\partial_zB , \end{gather}$$

(4.36d)$$\begin{gather}\bar{\mathcal{T}}^{(33)} = \tfrac{1}{2} \bar{W}^{(22)}\partial_zB ,\end{gather}$$

(4.36e)$$\begin{gather}\bar{\mathcal{T}}^{(40)} = \tfrac{1}{2} \mathcal{R}(\bar{W}^{(20)} )^2+ \tfrac{1}{4} |\bar{W}^{(22)}|^2 + \tfrac{1}{2} \bar{W}^{(1)}(\bar{W}^{(31)})^* + \tfrac{1}{8}|\bar{W}^{(1)}|^2|(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)A|^2\nonumber\\ + \tfrac{1}{16}(\bar{W}^{(1)})^2[(\boldsymbol{\nabla} A+\mathrm{i}\boldsymbol{k}_0 A)^*]^2 , \end{gather}$$

(4.36f)$$\begin{gather}\bar{\mathcal{T}}^{(42)} = \mathcal{R}\left[\bar{W}^{(20)} \right] \bar{W}^{(22)} + \tfrac{1}{2} \bar{W}^{(1)}\bar{W}^{(31)} +\tfrac{1}{2} (\bar{W}^{(1)})^*\bar{W}^{(33)} \nonumber\\ + \tfrac{1}{8} |\bar{W}^{(1)}|^2\left[(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)A\right]^2 + \tfrac{1}{8} (\bar{W}^{(1)})^2|\boldsymbol{\nabla} A+\mathrm{i}\boldsymbol{k}_0A|^2 , \end{gather}$$

(4.36g)$$\begin{gather}\bar{\mathcal{T}}^{(44)} =\tfrac{1}{4}( \bar{W}^{(22)} )^2 + \tfrac{1}{2} \bar{W}^{(1)}\bar{W}^{(33)} +\tfrac{1}{16} (\bar{W}^{(1)})^2 (\boldsymbol{\nabla} A+\mathrm{i}\boldsymbol{k}_0A)^2, \end{gather}$$

for $z=0$. Using envelope $A$ and $B_s$ as input, envelopes $B_0^{(mj)}$, $\bar {W}^{(mj)}$, $\bar {\mathcal {W}}^{(mj)}$ and $\bar {\mathcal {T}}^{(mj)}$ are obtained in sequence from the lowest to higher orders in wave steepness, which will be directly used in the CEEEs derived in the following section.

6.1.1. A train of Stokes waves

A train of Stokes waves is considered to have a wave vector of $\boldsymbol {k}_0$ and phase of $\theta _0$ and we choose $\alpha =1$ and $\beta =1$ for the implementation of the CEEEs. In a traditional perturbation method as noted in § 3.2, it is typical of solving for the wave-perturbed parameters on a still water surface, in contrast to the CEEEs method where the primary unknowns are those defined on the free water surface. In the CEEEs, envelope $A=A(\boldsymbol {x},t)$ and $B_s=B_s(\boldsymbol {x},t)$ depend on both time and the horizontal position vector. They can be expressed in the form of a power series of the wave steepness $\epsilon _0$, up to the second order,

(6.1a,b) \begin{align} A &= \epsilon_0 A_1 + \epsilon_0^2 A_2\exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right) \quad \text{and}\notag\\ B_s &= \epsilon_0 B_{s,1} + \epsilon_0^2 B_{s,2}\exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right), \end{align}

where the constant mean of both the envelope and the potential are neglected for simplicity but can be additionally considered; $\epsilon _0$ denotes the small dimensionless wave steepness in a traditional perturbation method as noted in § 3.2; subscript ‘1’ and ‘2’ denotes the envelopes at first and second order in wave steepness, $\epsilon _0$, respectively; $A_1=A_1(\epsilon _0^2t)$ and $B_{s,1}=B_{s,1}(\epsilon _0^2t)$ are a real and imaginary time-dependent amplitude of the linear elevation and the potential on a still water surface in a traditional perturbation method, respectively (see, e.g. Fenton Reference Fenton1985). Similarly, a leading-order approximation to envelope $\bar {W}^{(1)}$ is, to $O(\epsilon ^2_0)$,

(6.2)\begin{equation} \bar{W}^{(1)} = k_0\left(\epsilon_0\tanh (k_0h) B_{s,1} + 2\epsilon_0^2\tanh 2k_0h B_{s,2}\exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right) \right), \end{equation}

where the factor of $2$ arises from the derivative with respect to the vertical axis due to the second-order superharmonic waves. Inserting (6.1a,b) and (6.2) into the elevation, the potential on the free water surface and the velocity potential at an arbitrary depth in the framework of the CEEEs leads to, up to second order in wave steepness $\epsilon _0$,

(6.3a)$$\begin{gather} \zeta(\boldsymbol{x},t) = \tfrac{1}{2}\epsilon_0 A_1\exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right) +\tfrac{1}{2} \epsilon_0^2 A_2\exp\left({2\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t+\theta_0)}\right) + \text{c.c.}, \end{gather}$$

(6.3b)$$\begin{gather}\psi(\boldsymbol{x},t) = \tfrac{1}{2} \epsilon_0 B_{s,1}\exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right) + \tfrac{1}{2} \epsilon_0^2 B_{s,2}\exp\left({2\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t+\theta_0)}\right) +\text{c.c.}, \end{gather}$$

(6.3c)$$\begin{gather} \varPhi(\boldsymbol{x},z,t) = \dfrac{1}{2} \epsilon_0 B_{s,1} \dfrac{\cosh k_0(z+h)}{\cosh k_0h}\exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right) \nonumber\\ +\dfrac{1}{2} \epsilon_0^2\left(B_{s,2}- \dfrac{1}{2} k_0\tanh k_0h A_1 B_{s,1}\right)\dfrac{\cosh 2k_0(z+h)}{\cosh 2k_0h} \exp\left({2\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t+\theta_0)}\right) . \end{gather}$$

Therefore, the CEEEs for the evolution of a train of Stokes waves can be much simplified to

(6.4a,b)\begin{equation} (\partial_t-\mathrm{i}\omega_0) A - \bar{W}^{(1)} = \mathcal{N}_{A,2} \quad \text{and}\quad (\partial_t-\mathrm{i}\omega_0) B_s + g A = \mathcal{N}_{B,2}, \end{equation}

where

(6.5a,b)\begin{equation} \mathcal{N}_{A,2} = \mathcal{N}_A^{(22)} \exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right) \quad \text{and}\quad \mathcal{N}_{B,2} = \mathcal{N}_B^{(22)} \exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right) , \end{equation}

and

(6.6a)$$\begin{gather} \mathcal{N}_A^{(22)} = \tfrac{1}{2} \left( 1- 2\tanh k_0h \tanh 2k_0h \right)k_0^2 A_1B_{s,1}, \end{gather}$$

(6.6b)$$\begin{gather}\mathcal{N}_B^{(22)} = \tfrac{1}{4} k_0^2(1+\tanh^2k_0h ) B_{s,1}^2. \end{gather}$$

Inserting $A_1=A_1(\epsilon _0^2t)$ and $B_{s,1}=B_{s,1}(\epsilon _0^2t)$ into (6.6a,b) and performing the analysis of scales gives rise to

(6.7a)$$\begin{gather} O\left( \mathcal{N}^{(22)}_A, \mathcal{N}^{(22)}_B \right) \sim O(A_1B_{s,1}) \sim \epsilon_0^2 \quad \text{and} \end{gather}$$

(6.7b)$$\begin{gather}O\left( \partial_t\mathcal{N}^{(22)}_A, \partial_t\mathcal{N}^{(22)}_B \right) \sim O(\epsilon_0\partial_t A_1) \sim O(\epsilon_0\partial_t B_{s,1}) \sim \epsilon_0^4, \end{gather}$$

and therefore, $\mathcal {N}_{A}^{(mj)} = \mathcal {N}_{A}^{(mj)}(\epsilon _0^2t)$ and $\mathcal {N}_{B}^{(mj)}=\mathcal {N}_{B}^{(mj)}(\epsilon _0^2t)$, which suggests $\varepsilon _\text {f}=\epsilon _0^2$ as is well known for the temporal evolution of the amplitude of a train of Stokes waves (Zakharov Reference Zakharov1968; Fenton Reference Fenton1985), e.g. $A_1\equiv A_1(\epsilon ^2_0t)=A_1(\varepsilon _\text {f} t)$ as noted. A traditional perturbation method proposes to solve the equations in sequence from the first to higher orders in wave steepness $\epsilon _0$ and thereby we obtain

(6.8a,b)\begin{equation} O(\epsilon_0): -\mathrm{i}\omega_0 A_1 - k_0\tanh k_0h B_{s,1} = 0 \quad \text{and}\quad -\mathrm{i}\omega_0 B_{s,1} - gA_1 = 0, \end{equation}

which are well known for the linear evolution of monochromatic waves; at $O(\epsilon _0^2)$,

(6.9a)$$\begin{gather} A_2\exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right) -\dfrac{\mathrm{i} k_0\tanh 2k_0h}{\omega_0} B_{s,2}\exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right) \nonumber\\ = \dfrac{\mathrm{i}}{2\omega_0} \mathcal{N}_A^{(22)}\exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right) , \end{gather}$$

(6.9b)$$\begin{gather} 2\mathrm{i}\omega_0 B_{s,2}\exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right) + gA_2\exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right) \nonumber\\ =- \mathcal{N}^{(22)}_B\exp\left({\mathrm{i} ( \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\omega_0t +\theta_0)}\right) , \end{gather}$$

where, despite the fact that it can be eliminated, the factor $\exp \left ({\mathrm {i} ( \boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\omega _0t +\theta _0)}\right )$ was kept with the intention to demonstrate clearly that the numerical solution of (6.9a,b) for the unknown parts of envelope $A$ and $B_s$, i.e. $A_2\exp \left ({\mathrm {i} ( \boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\omega _0t +\theta _0)}\right )$ and $B_{s,2}\exp \left ({\mathrm {i} ( \boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\omega _0t +\theta _0)}\right )$, depends only on the temporal scale of the forcing terms $\mathcal {N}^{(22)}_A$ and $\mathcal {N}^{(22)}_B$. Especially in the limiting case of second-order Stokes waves, the solution of (6.9a,b) obtained numerically predicts no mathematical truncation error (Atkinson Reference Atkinson2008, § 1.3) as the unknown envelopes are approximately constants, and, therefore, the numerical solution can return the same results as the analytical solution from trivial algebras. Comparing the expressions given by (6.3a–c) in the CEEEs and those by (3.12a,b) for the same parameter, the following relations hold:

(6.10a–c)\begin{equation} A'_{mm} = A_{m}, \quad B'_{11} = B_{s,1}, \quad \text{and}\quad B'_{22} = B_{s,2} -\tfrac{1}{2} k_0\tanh k_0h A_{1}B_{s,1}. \end{equation}

Here $m=1$ and $m=2$. Replacing the parameters in (6.9a,b) with those used in a traditional perturbation method through the relations (6.10a–c), eliminating the factor $\exp \left ({\mathrm {i} ( \boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\omega _0t +\theta _0)}\right )$ and using the relation $(\partial _{t}-\mathrm {i}\omega _0)A_{11}' =\partial _{z}B'_{11}$ readily leads to the envelope equations given by (3.13). Hence, the relation between the CEEEs and a traditional perturbation method for the evolution of a train of weakly nonlinear Stokes waves has been established.

6.1.2. The semi-analytical framework for directionally spread broadband waves

Similar to § 6.1.1, the main focus of this section is to show how the CEEEs lead to the framework by Li & Li (Reference Li and Li2021) for the evolution of three-dimensional broadband waves with large directional spreading. Again, we start from envelope $A$ and $B_s$ in the CEEEs in the form of a power series in $\epsilon _0$ up to the second order in $\epsilon _0$,

(6.11a)$$\begin{gather} A = \epsilon_0 A_{11} +\epsilon^2_0 A_{20}\exp\left({-\mathrm{i}(\alpha\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right)+ \epsilon_0^2 A_{22}\exp\left({\mathrm{i} (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) , \end{gather}$$

(6.11b)$$\begin{gather}B_s = \epsilon_0 B_{s,11} + \epsilon^2_0 B_{s,20}\exp\left({-\mathrm{i}(\alpha\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right)+ \epsilon_0^2 B_{s,22}\exp\left({\mathrm{i} (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) , \end{gather}$$

where $A_{mj}=A_{mj}(\boldsymbol {x},t)$ and $B_{s,mj}=B_{s,mj}(\boldsymbol {x},t)$. The relationship between (6.11b) and the envelope of the potential and vertical velocity on the free water surface leads to

(6.12a)$$\begin{gather} B(\boldsymbol{x},z,t) =\epsilon_0 B_{11}(\boldsymbol{x},z,t)+ \epsilon_0^2 B_{20} \exp\left({-\mathrm{i}(\alpha\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right)\notag\\ + \epsilon_0^2 B_{22}\exp\left({\mathrm{i} (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) \end{gather}$$

(6.12b)$$\begin{gather} \bar{W}^{(1)}(\boldsymbol{x},t) = \epsilon_0 \bar{W}^{(1)}_{11}(\boldsymbol{x},0,t)+ \epsilon_0^2 \bar{W}^{(1)}_{20}(\boldsymbol{x},0,t) \mathrm{e}^{-\mathrm{i}(\alpha\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)} \nonumber\\ +\epsilon_0^2 \bar{W}^{(1)}_{22}(\boldsymbol{x},0,t)\exp\left({\mathrm{i} (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) , \end{gather}$$

with $B_{s,mj}(\boldsymbol {x},t) = B_{mj}(\boldsymbol {x},0,t)$ and

(6.12c)$$\begin{gather} B_{mj} = \int_{-\infty}^{\infty} \hat{B}_{s,mj} \dfrac{\cosh |\boldsymbol{k}+j\alpha \boldsymbol{k}_0|(z+h)}{\cosh |\boldsymbol{k}+j\alpha \boldsymbol{k}_0|h}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{k}, \end{gather}$$

(6.12d)$$\begin{gather}\bar{W}^{(1)}_{mj}(\boldsymbol{x},z,t) = \partial_zB_{mj} \equiv \int_{-\infty}^{\infty} |\boldsymbol{k}+j\alpha \boldsymbol{k}_0|\hat{B}_{s,mj}\dfrac{\sinh |\boldsymbol{k}+j\alpha \boldsymbol{k}_0|(z+h)}{\cosh |\boldsymbol{k}+j\alpha \boldsymbol{k}_0|h}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{k}. \end{gather}$$

We next derive the parameters needed at second order in wave steepness in the CEEEs. Substituting the envelopes on the still water surface given by (6.11a,b) and (6.12a,b) into the envelopes of the approximate wave fields at the second order in wave steepness presented in § 4.3.2, and keep the terms up to $O(\epsilon _0^2)$ gives rise to

(6.13a)$$\begin{gather} B_0^{(20)} =- \epsilon_0^2\left[\mathcal{R}\left(\tfrac{1}{2} A^*_{11}B_{s,11}\right) \right]^{[0]}_+, \end{gather}$$

(6.13b)$$\begin{gather}B_0^{(22)} =-\tfrac{1}{2} \epsilon^2_0 A_{11}B_{s,11}, \end{gather}$$

(6.13c)$$\begin{gather}\bar{W}^{(22)} = \epsilon_0^2\int_{-\infty}^{\infty} |2(\boldsymbol{k}_0+\boldsymbol{k})| \hat{B}_0^{(22)} \dfrac{\sinh |(2\alpha\boldsymbol{k}_0+\boldsymbol{k})(z+h)|}{\cosh |2(\boldsymbol{k}_0+\boldsymbol{k})h|}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{k} + \dfrac{1}{2} \epsilon_0^2 A_{11} \partial_z\bar{W}^{(1)}_{11}, \end{gather}$$

(6.13d)$$\begin{gather}\bar{W}^{(20)} = \epsilon_0^2 \int_{-\infty}^{\infty} |\boldsymbol{k}| \hat{B}_0^{(20)} \dfrac{\sinh |\boldsymbol{k}|(z+h)}{\cosh |\boldsymbol{k}| h}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{k} + \epsilon_0^2\mathcal{R} \left(\dfrac{1}{2} A^*_{11} \partial_z\bar{W}^{(1)}_{11} \right). \end{gather}$$

Similarly, inserting the resulting expressions for the envelopes in the first and second orders in this section into the nonlinear forcing terms in the second order in the CEEEs leads to

(6.14a)\begin{align} \mathcal{N}_A^{(22)} &= \epsilon_0^2 \int_{-\infty}^{\infty} |2(\boldsymbol{k}_0+\boldsymbol{k})| \hat{B}_0^{(22)} \dfrac{\sinh |(2\alpha\boldsymbol{k}_0+\boldsymbol{k})(z+h)|}{\cosh |2(\boldsymbol{k}_0+\boldsymbol{k})h|}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{k} \nonumber\\ &\quad +\epsilon_0^2 \left[\dfrac{1}{2} A_{11} \partial_z\bar{W}^{(1)}_{11} -\dfrac{1}{2} (\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)A_{11}\boldsymbol{\cdot}(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)B_{s,11} \right],\end{align}

(6.14b)\begin{align} \mathcal{N}_A^{(20)} &= \epsilon_0^2 \int_{-\infty}^{\infty} |\boldsymbol{k}| \hat{B}_0^{(20)} \dfrac{\sinh |\boldsymbol{k}|(z+h)}{\cosh |\boldsymbol{k}| h}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{k},\nonumber\\ &\quad + \epsilon_0^2 \left[ \mathcal{R} \left(\dfrac{1}{2} A^*_{11} \partial_z\bar{W}^{(1)}_{11} -\dfrac{1}{2} (\boldsymbol{\nabla}-\mathrm{i}\boldsymbol{k}_0)A^*_{11}\boldsymbol{\cdot}(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)B_{s,11} \right)\right]_+^{[0]}, \end{align}

(6.14c)\begin{align} \mathcal{N}_B^{(22)} &= \dfrac{1}{4} \left[-(\boldsymbol{\nabla} B_{s,1}+\mathrm{i}\boldsymbol{k}_0B_{s,1})^2 +(\bar{W}_{11}^{(1)})^2 \right], \end{align}

(6.14d)\begin{align} \mathcal{N}_B^{(20)} &= \dfrac{1}{4} \left[-|\boldsymbol{\nabla} B_{s,1}+\mathrm{i}\boldsymbol{k}_0B_{s,1}|^2 +|\bar{W}_{11}^{(1)}|^2 \right]_+^{[0]}. \end{align}

Inserting the parameters in the form of a power series of $\epsilon _0$ back to the CEEEs, collecting the terms in the second order in $\epsilon _0$ and separating the wave harmonics leads to the evolution equations for the second-order super-harmonic waves

(6.15a,b)\begin{align} (\partial_t-2\mathrm{i}\beta\omega_0) A_{22} - \bar{W}^{(1)}_{22} = \mathcal{N}^{(22)}_A \quad \text{and}\quad (\partial_t-2\mathrm{i}\beta\omega_0) B_{s, 22} + g A_{22} = \mathcal{N}^{(22)}_B , \end{align}

with the factor $\exp \left ({\mathrm {i} (\alpha \boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\beta \omega _0t)}\right )$ eliminated on both sides of the equation, and the equations for the subharmonic waves

(6.16a,b)\begin{equation} \partial_t A_{20} - \bar{W}^{(1)}_{20} = \mathcal{N}^{(20)}_A \quad \text{and}\quad \partial_t B_{s, 20} + g A_{20} = \mathcal{N}^{(20)}_B, \end{equation}

with the factor $\exp {[-\mathrm {i}(\alpha \boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\beta \omega _0t)]}$ eliminated. The next step is to establish the relationship between $B_{s,22}$ and $\bar {W}^{(1)}_{22}$ and the corresponding parameter given by Li & Li (Reference Li and Li2021), respectively. Comparing (6.11a) for the envelope of the elevation in the CEEEs and that given by (3.12c) in the framework by Li & Li (Reference Li and Li2021), it is evidently understood that $A_{mj}= A_{mj}'$ and

(6.17a)$$\begin{gather} B'_{11}(\boldsymbol{x},0,t)= B_{11}(\boldsymbol{x},0,t)=B_{s,11}, \end{gather}$$

(6.17b)$$\begin{gather}B'_{22}(\boldsymbol{x},0,t) = B_{s,22} - \tfrac{1}{2} A'_{11}\partial_zB'_{11}, \end{gather}$$

(6.17c)$$\begin{gather}B'_{20}(\boldsymbol{x},0,t) = B_{s,20} -\left[\mathcal{R}\left(\tfrac{1}{2} A^{'*}_{11}B'_{11}\right) \right]^{[0]}_+, \end{gather}$$

(6.17d)$$\begin{gather}\partial_zB'_{22}(\boldsymbol{x},0,t) = \bar{W}^{(1)}_{22} + \partial_zB^{(22)}(\boldsymbol{x},0,t), \end{gather}$$

(6.17e)$$\begin{gather}\partial_zB'_{20}(\boldsymbol{x},0,t) = \bar{W}^{(1)}_{20} +\partial_zB^{(20)}(\boldsymbol{x},0,t). \end{gather}$$

Replacing the terms in the CEEEs with those by Li & Li (Reference Li and Li2021) through the connections given by (6.17a–e) leads to the envelope equations of Li & Li (Reference Li and Li2021) given by (3.12), for which the identity $(\partial _t-\mathrm {i}\omega _0)A'_{11}= \partial _zB'_{11}$ from the linearized kinematic condition on a still water surface was used.

Based on (6.17a–e) and the relations

(6.18a)$$\begin{gather} \zeta'_2(\boldsymbol{x},t) = \tfrac{1}{2} A'_{22}\exp\left({2\mathrm{i} (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \tfrac{1}{2} A'_{20} + \text{c.c.}, \end{gather}$$

(6.18b)$$\begin{gather}\varPhi'_2(\boldsymbol{x},z,t) = \tfrac{1}{2} B'_{22}\exp\left({2\mathrm{i} (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \tfrac{1}{2} B'_{20} + \text{c.c.}, \end{gather}$$

the evolution equations for the second-order elevation $\zeta _2'$ and potential $\varPhi _2'(\boldsymbol {x},0,t)$ can be derived, which conform with (3.11).

Again, (6.15a,b) and (6.16a,b) in the framework of CEEEs suggest that the computational accuracy and efficiency for numerical solutions relies only on the spatial and temporal scales of the nonlinear forcing terms $\mathcal {N}_A^{(mj)}$ and $\mathcal {N}_B^{(mj)}$, in a way similar to the numerical implementation of (3.13a,b) and (3.14a,b) by Li & Li (Reference Li and Li2021). In particular, Li & Li (Reference Li and Li2021) have demonstrated that the (second-order accurate) envelope framework due to (3.13a,b) and (3.14a,b) permits a significant improvement in the numerical efficiency at no cost to the accuracy for weakly nonlinear waves, compared with the numerical simulations based on the HOS method (Dommermuth & Yue Reference Dommermuth and Yue1987; West et al. Reference West, Brueckner, Janda, Milder and Milton1987) and the Fourier kernels (Hasselmann Reference Hasselmann1962; Zakharov Reference Zakharov1968; Dalzell Reference Dalzell1999). Therefore, it can be conjectured that the CEEEs should have the same numerical advantage, at least up to second order in wave steepness due to the same numerical features from comparing (6.15a,b) and (6.16a,b) from the CEEEs with (3.13a,b) and (3.14a,b) by Li & Li (Reference Li and Li2021), respectively.

6.2.1. Comparisons of the nonlinear forcing terms

As noted above, the differences between the CEEEs and the HOS method mainly arise from the different choices of the main unknowns that are solved for numerically. They should in principle return exactly the same numerical results if the initial conditions are the same. As the detailed derivations in § 4 have analytically demonstrated this point, we choose to keep the numerical validation to a minimum to avoid unnecessary repetition. In particular, we focus on the fact that the nonlinear forcing terms given by (3.7) and (3.8) should be identical to (4.34a) and (4.34b) at different orders (i.e. different values of $m$) in wave steepness. The HOS method and the CEEEs shall require different parameters as input, the former of which uses $\zeta (\boldsymbol {x},t)$ and $\psi (\boldsymbol {x},t)$ whereas the latter $A$ and $B_s$. The relation given by (2.11a) between $\zeta$ (or $\psi$) and $A$ (or $B_s$) should hold for all the times if both methods are used for the computation of the same case.

A case of a right-propagating (i.e. in the positive $x$ direction) focused wave group in two dimensions was chosen to show the comparisons between the two methods. In particular, the parameters for the HOS method at a time instant were prescribed according to

(6.19a,b)\begin{align} \zeta = \mathcal{R}\left[A_p\int \dfrac{|\hat{\zeta}(\omega_g)|}{\sqrt{m_0}} E(x,t)\,\mathrm{d} \omega_g\right] \quad \text{and}\quad \psi = \mathcal{R}\left[ A_p\int \dfrac{-\mathrm{i} g|\hat{\zeta}(\omega_g)|}{\omega_g\sqrt{m_0}} E(x,t) \,\mathrm{d}\omega\right], \end{align}

where the factor $E(x,t)= \exp \{ \mathrm {i} [k (x-x_\text {f})-\omega _g(t-t_\text {f})+\theta _\text {f}]\}$, with $A_p$, $x_\text {f}$, $t_\text {f}$ and $\theta _\text {f}$ the prescribed peak (real) amplitude, position, time and phase for the wave group at linear focus, respectively; $k$ and $\omega _g$ denote the wavenumber and angular frequency of a train of a right-propagating monochromatic wave, respectively, obeying the linear relationship $\omega _g=\omega (k,h)$ as defined in § 2.4; $m_0$ denotes the zeroth moment of a JONSWAP power energy spectrum $S(\omega _g)$ used with the enhancement peak factor $\gamma$ of $3.3$, and $|\hat {\zeta }(\omega _g)|=\sqrt {2S(\omega _g)\Delta \omega _g}$ denotes the amplitude evaluated based on the JONSWAP with the interval between two adjacent frequencies prescribed on a numerical grid. It is highlighted that the initial conditions given by (6.19a,b) for the initial time instant $t=t_0$ are not an exact solution to the fully nonlinear potential flow boundary value problem but are the exact solution of the linearised problem. These conditions will not affect the comparisons of the nonlinear forcing terms $\mathcal {W}_M$ and $\mathcal {T}_M$ presented in this section between the HOS method and CEEEs as long as the initial conditions are consistent. A detailed procedure for the initialization of nonlinear waves for the solutions of an initial-value problem can be found in works of, for example, Dommermuth (Reference Dommermuth2000) and Slunyaev, Sergeeva & Didenkulova (Reference Slunyaev, Sergeeva and Didenkulova2016) among others.

Similarly, envelope $A$ and $B_s$ are expressed as, respectively,

(6.20a,b)\begin{equation} A = A_p \int \dfrac{|\hat{\zeta}(\omega_g)|}{\sqrt{m_0}}E_c(x,t)\,\mathrm{d}\omega_g \quad \text{and}\quad B_s = A_p \int \dfrac{-\mathrm{i} g|\hat{\zeta}(\omega_g)|}{\omega\sqrt{m_0}} E_c(x,t)\,\mathrm{d}\omega_g, \end{equation}

with $E_c(x,t)=\exp \{ \mathrm {i} [(k-k_p) (x-x_\text {f})-\omega _g(t-t_\text {f})+\omega _pt+\theta _\text {f}]\}$, which suggest that the following values were chosen for different parameters in the implementation of the CEEEs: $k_0=k_p$ and $\omega _0=\omega _p$ with $k_p=0.045$ m$^{-1}$ and $\omega _p=\omega (k_ph)$ the peak wavenumber and frequency of the JONSWAP spectrum, respectively, $\alpha =\beta =1$ and $k_ph=1.395$. We note that the elevation $\zeta$ (potential $\psi$) given by (6.19a) (6.19b) and the elevation envelope $A$ (the potential envelope $B_s$) given by (6.20a) (6.20b) should obey the equation (2.11a) (2.11b) as it indeed does. A comparison of the nonlinear forcing terms at second, third and fourth order in wave steepness is shown in figures 2 and 3 for two different time instants using two different values for the dimensional depth, $k_ph$. The computations from the HOS method were obtained from inserting (6.19a,b) for the elevation and potential into (3.7) and (3.8) and the results based on the CEEEs were through substituting (6.20a,b) for the envelopes into (4.34a) and (4.34b). In the case implemented in figures 2 and 3, an unrealistically high value of $0.8$ for the dimensionless wave steepness $\epsilon =k_pA_p$ was used with the intention of demonstrating the possible differences between the two methods, if any. Figures 2 and 3 show evidently good agreement of the nonlinear terms in all orders (in wave steepness) between the HOS method and CEEEs. In particular, the panels in the lowest two bottom rows of figures 2 and 3 suggest that the absolute differences of the predictions between the two methods are a few orders of magnitude smaller than the individual predictions in the highest order in wave steepness (e.g. compared with the nonlinear forcing terms in the fourth order in wave steepness). The results shown in figures 2 and 3 are in accordance with the theoretical derivations presented in § 4.

Figure 2. Comparisons of the nonlinear forcing terms (panels (a–f, i–vi)) at different orders in wave steepness, and the differences of the nonlinear forcing terms $\mathcal {W}_{M}$ and $\mathcal {T}_{M}$ for $M=4$ (panels (g,h, vii, viii)), between the HOS method (blue dashed) and the CEEEs-based model (red dot-dashed) based on the equations presented in §§ 3.1.1 and 4.3.3, respectively. Results are shown for (a–h) $t=-15\times T_p$ and (i–viii) $t=0\times T_p$ for the wavepacket at the linear focus, with $T_p$, $c_{g,p}$ and $\lambda _p$ the period, group velocity and wavelength of the spectral peak wave of a JONSWAP spectrum, respectively; $k_pA_\text {f}=0.8$ and $k_ph=1.5$ were used, where $k_p=2{\rm \pi} /\lambda _p$, $h$ is the water depth and $A_\text {f}$ is the amplitude of the focus wave at linear focus; $R(\chi )$ denotes the absolute differences of an arbitrary field $\chi$ obtained from CEEEs and the HOS method.

Figure 3. Caption is the same as figure 2 but $k_ph=3$.

6.2.2. Computational complexity

We compare the numerical performances between the HOS method and CEEEs for which both are based on FFTs for the evaluation of the nonlinear terms in the evolution equations. It should be noted that the comparison is not to seek the optimal numerical methods for computations. Instead, it aims to examine three main aspects essential to numerical efficiency for achieving the same level of accuracy: (i) the number of FFTs (including inverse FFTs) needed for advancing one time step; (ii) the criterion for the choice of the domain spacing and size and the time interval for numerical convergence and instability; (iii) the mathematical truncation error; see, e.g. Atkinson (Reference Atkinson2008, § 1.3), which indicates the numerical accuracy that can be achieved. The comparisons of these aspects are shown in table 2.

Table 2. Comparison of the numerical performances for convergent and stable computations between the HOS method and the CEEEs. The columns below the order (i.e. $M$) of accuracy in wave steepness ($\epsilon$) show the total number of FFTs (including inverse FFTs) needed for advancing one time step. In the table, $f_s$ and $k_s$ denote the frequency and wavenumber of the shortest wave that can be represented numerically to a sufficiently good level of accuracy, respectively; $N_s$ denotes the selected number of points used per characteristic wavelength; $\epsilon$, $\varepsilon _\text {f} = (\,f_s-M\beta f_0)/f_s$ ($0\lesssim \varepsilon _\text {f}\leq 1$) and $\varepsilon _k=(k_s-M\alpha k_0)/k_s$ ($0\lesssim \varepsilon _k\leq 1$) denote the dimensionless wave steepness and characteristic dimensionless bandwidth in frequency and wavenumber, respectively, with $f_s\gtrsim \beta Mf_0$, $k_s\gtrsim \alpha Mk_0$, $k_s=M\alpha k_0/(1-\varepsilon _k)$, and $\alpha$ and $\beta$ arbitrarily chosen non-negative parameters. Different choices of approximate methods were made for the time integration in the numerical implementation of the HOS method; ‘FE1’, ‘MP2’ and ‘RK4’ denote the first-, second- and fourth-order accurate forward Euler, mid-point and Runge–Kutta method for the time integration, respectively, and ‘ExpInt1’, ‘ExpInt2’ and ‘ExpInt3’ denote a first-, second- and third-order accurate approximate method chosen for the time integration based on an exponential integrator, respectively.

A total number of FFTs (including inverse FFTs) needed for updating one time step depends on a specific numerical algorithm for time integration. Therefore, in order to make the comparisons as fair as possible for the HOS method, three widely known numerical algorithms are chosen to this end, including the first-, second- and fourth-order accurate forward Euler (‘FE1’), mid-point (‘MP2’) and Runge–Kutta (explicit, ‘RK4’) methods, respectively. The main difference between the use of an exponential integrator and the previous three algorithms lie in the fact that a numerical algorithm for the former and the latter for time integration is carried out in the Fourier $\boldsymbol {k}$ and physical plane, respectively, due to which the former leads to an accurate evaluation of the temporal evolution arising from the linear terms in equations. How to apply an exponential integrator for implementing the CEEEs is presented in § 5.1. Similar procedures can also be taken for the HOS method. An exponential integrator is used in table 2 for both the HOS method and the CEEEs. An exponential integrator together with a first-order accurate numerical algorithm to approximate the time integral in the Fourier $\boldsymbol {k}$ space (‘ExpInt1’) is sufficient to demonstrate the numerical features of the CEEEs, as shown in table 2, whereas additional second- (‘ExpInt2’) and third-order (‘ExpInt3’) accurate algorithms are listed for the HOS method. How to optimize the numerical implementation of the CEEEs is not the main focus of this paper and is open for studies in future works.

We introduce a few parameters for the discussion about numerical performances. Let $L_s$, $f_s$ and $k_s$ be the wavelength, frequency and wavenumber of the shortest wave that can be represented numerically to a sufficiently good level of accuracy with $k_sL_s=2{\rm \pi}$. The use of FFTs and inverse FFTs requires an evenly spaced computational domain that has a characteristic length of $L$ and the spacing of $|\Delta \boldsymbol {x}_n|$ between two adjacent discrete points. The following dimensionless parameters are defined:

(6.21a–c)\begin{equation} \varepsilon_t = f_s\Delta t_n, \quad \varepsilon_\text{f} = \dfrac{f_s-\beta Mf_0}{f_s} , \quad \text{and}\quad \varepsilon_k =\dfrac{k_s-\alpha Mk_0}{k_s}. \end{equation}

Here $\varepsilon _t$ denotes the dimensionless time interval for indicating the mathematical truncation error for computing the time integration, $f_0=\omega _0/(2{\rm \pi} )$ and $\Delta t_n$ the time interval. It should be noted that $f_s\gtrsim \beta Mf_0$ and $k_s\gtrsim \alpha Mk_0$ are assumed and should hold in most practical situations as the range $0\leq \alpha \lesssim 2$ and $0\leq \beta \lesssim 2$ are recommended. Thus, the inequalities $0\lesssim \varepsilon _\text {f}\leq 1$ and $0\lesssim \varepsilon _k\leq 1$ hold.

Examining the algorithms that implement the HOS method in table 2, we find that an exponential integrator together with a numerical algorithm for the time integration seems to have advantages in the numerical performances against using the other three methods. It is shown in table 2 that the first-order accurate forward Euler algorithm requires the time interval to be extremely small for achieving a sufficient level of accuracy. As a result, the additional computational cost needed due to a larger number of time steps would not be made up by the computational efficiency saved from a smaller number of FFTs required per step, compared with the other three methods. For weakly nonlinear waves where the wave steepness $\epsilon \to 0^+$ and $\varepsilon _t\ll 1$ for convergent numerical results, and, therefore, the assumption of the scales $O(\epsilon )\sim O(\varepsilon _t)$ can be made, the computational efficiency of an exponential integrator (i.e. ‘ExpInt1’ and ‘ExpInt2’) is superior to the mid-point method (‘MP2’) due to either a lesser number of FFTs or a higher accuracy when the other aspects are kept the same. For steeper waves where the wave steepness can be one order of magnitude larger compared with weakly nonlinear waves, and, thus, $O(\epsilon )\sim O(\sqrt {\varepsilon _t})$ can be assumed, the same conclusion can be drawn for a second-order accurate exponential integrator (‘ExpInt2’) as it leads to obviously a higher level of accuracy despite a slightly increased cost arising from a larger number of FFTs. Compared with the Runge–Kutta algorithm (‘RK4’), the third-order accurate exponential integrator (‘ExpInt3’) shows a slightly better performance due to a smaller number of FFTs required for achieving the same level of accuracy with $O(\epsilon )\sim O(\varepsilon _t)$.

Due to the above discussion about the HOS method examined in table 2, the comparison between the CEEEs and the HOS method will focus only on their implementation through using an exponential integrator that has been demonstrated to have favoured features for the HOS method. The numerical performances of the CEEEs shown in table 2 especially depend on the choice of the values for $\alpha$ and $\beta$. Thereby, discussions are made based on the categories of the choices listed in the following, starting from the least advantageous category for the CEEEs.

1. (a) $\alpha =\beta =0$ and, therefore, $\varepsilon _\text {f}=\varepsilon _k\equiv 1$. It suggests that the CEEEs are simply the HOS but with the newly introduced envelope transform for the evaluation of the nonlinear forcing terms and wave parameters. The newly proposed evaluation through the envelope transform is at the expense of an increased number of FFTs and, therefore, an additional computational cost with all the other parameters the same as the HOS method (‘ExpInt1’), as clearly seen in table 2. As a result, this category of the value for $\alpha$ and $\beta$ should be dropped out if the CEEEs will be implemented for numerical computations as it does not introduce additional advantages compared with the HOS method.

2. (b) $\alpha =0$ and $0< \beta$ and, thus, $\varepsilon _k=1$ and $0\lesssim \varepsilon _\text {f}< 1$, respectively. This suggests that a (small) dimensionless bandwidth parameter (i.e. $\varepsilon _\text {f}$) in wave frequency has been introduced. The characteristic wave frequency, $\omega _0$, can be chosen as the angular frequency of either a carrier wave or the peak wave of a wave spectrum. A specific value for $\beta$ can be selected that permits $\varepsilon _\text {f}\sim \varepsilon _t$ in practice. For instance, $\beta =1$ has been used in § 6.1.1 that has led to $\varepsilon _\text {f}\sim \epsilon ^2$. Table 2 indicates that the first-order accurate exponential integrator (‘ExpInt1’) for the CEEEs can reach the same accuracy as a second-order accurate exponential integrator (‘ExpInt2’) for the HOS method but with a slightly smaller number of FFTs required. Moreover, the CEEEs can allow for a much larger time interval, e.g. $\varepsilon _t\sim O(1)$, for numerical convergence as long as the inequality $\varepsilon _\text {f}\varepsilon _t\ll 1$ holds. Physically, a larger time interval (e.g. $1/f_s\ll \Delta t_n$) achieved by the CEEEs attributes to the fact that it depends on the rate of change of a wave spectrum that has a much slower temporal scale relative to that (i.e. $1/f_s$) of the phase of the fastest wave. This feature is especially similar to a NLS equation-based model by using a split-step method as explained in Lo & Mei (Reference Lo and Mei1985).

3. (c) $\alpha >0$ and $\beta >0$ and thereby $0\lesssim \varepsilon _k< 1$ and $0\lesssim \varepsilon _\text {f}< 1$, respectively. Compared with category (b), additional advantageous features are introduced due to a small value permitted for $\varepsilon _k$ that denotes the bandwidth in wavenumber. In contrast to $\varepsilon _k=1$, the parameter $\varepsilon _k$ of a small value for the CEEEs suggests that a much larger domain with a length of $L\sim N_k |\Delta \boldsymbol {x}_n|/\varepsilon _k$ can be achieved at no expense to the computational efficiency and accuracy if the same number, $N_k$, of Fourier modes are used. This is regardless of the choice of the time-dependent parameters (e.g. $\varepsilon _\text {f}$ and $\varepsilon _t$). Compared with a second-order accurate exponential integrator (‘ExpInt2’) for the HOS method, a first-order accurate exponential integrator (‘ExpInt1’) implementing the CEEEs has the following three features. Firstly, it has a slightly smaller computational cost examining the number of FFTs needed. Secondly, it permits a much larger temporal scale for computational instability and convergence. Thirdly, with the same Fourier modes chosen for computations, a much larger domain can be allowed. It should be highlighted that the introduction of bandwidth parameters, $\varepsilon _k$ and $\varepsilon _\text {f}$, is similar to a NLS equation-based model where evolution of a slowly varying envelope is described and where the linear terms of the equation can be accurately solved with a split step together with a pseudo-spectral method (Li Reference Li2021). Both the CEEEs and a NLS equation-based model permit a relatively coarse spatial domain and larger time interval for numerics while allowing for resolving of the wave phase on a prescribed computational domain.

Recall that the primary objective of this paper is to propose a new framework that can combine the advantages of both the HOS method and a NLS equation-based model, in the sense that the new framework can reach the same accuracy as the HOS method and, similar to a NLS equation-based method, it permits for both a larger spatial domain and slower temporal scale but at no expense to the computational efficiency. The choice of values for $\alpha$ and $\beta$ indicated by category (c) has demonstrated that the newly derived CEEEs can indeed reach this objective.