On coupled envelope evolution equations in the Hamiltonian theory of nonlinear surface gravity waves

Abstract

This paper presents a novel theoretical framework in the Hamiltonian theory of nonlinear surface gravity waves. The envelope of surface elevation and the velocity potential on the free water surface are introduced in the framework, which are shown to be a new pair of canonical variables. Using the two envelopes as the main unknowns, coupled envelope evolution equations (CEEEs) are derived based on a perturbation expansion. Similar to the high-order spectral method, the CEEEs can be derived up to arbitrary order in wave steepness. In contrast, they have a temporal scale as slow as the rate of change of a wave spectrum and allow for the wave fields prescribed on a computational (spatial) domain with a much larger size and with spacing longer than the characteristic wavelength at no expense of accuracy and numerical efficiency. The energy balance equation is derived based on the CEEEs. The nonlinear terms in the CEEEs are in a form of the separation of wave harmonics, due to which an individual term is shown to have clear physical meanings in terms of whether or not it is able to force free waves that obey the dispersion relation. Both the nonlinear terms that can only lead to the forcing of bound waves and those that are capable of forcing free waves are demonstrated, in the case of the latter through the analysis of the quartet and quintet resonant interactions of linear waves. The relations between the CEEEs and two other existing theoretical frameworks are established, including the theory for a train of Stokes waves up to second order in wave steepness (Fenton, ASCE J. Waterway Port Coastal Ocean Engng, vol. 111, issue 2, 1985, pp. 216–234) and a semi-analytical framework for three-dimensional weakly nonlinear surface waves with arbitrary bandwidth and large directional spreading by Li & Li (Phys. Fluids, vol. 33, issue 7, 2021, 076609).

1. Introduction

For the description of non-breaking surface water waves in the open ocean and coastal waters, there are many available approaches in the framework of potential flow. The high-order spectral (HOS) method (Dommermuth & Yue 1987; West et al. 1987), Boussinesq-type formulations (Wei et al. 1995; Agnon, Madsen & Schäffer 1999), volumetric methods (Engsig-Karup, Bingham & Lindberg 2009; Bihs et al. 2020) and the fast computational method developed by Clamond & Grue (2001) are a few examples of these that can account for waves up to arbitrary order in wave steepness. In terms of the numerical efficiency for a given level of accuracy, it is without doubt that the HOS method is the preferred choice compared with the aforementioned alternatives (Klahn, Madsen & Fuhrman 2020). Two aspects have especially contributed to its high efficiency: (i) it permits an explicit method for the time integration and vertical velocity on the free water surface, and (ii) it takes advantage of spectral methods for numerical computations (Dommermuth & Yue 1987; Ducrozet et al. 2016).

In cases for the evolution of weakly nonlinear waves, a nonlinear Schrödinger (NLS) equation-based model (Benney & Newell 1967; Zakharov 1968; Chu & Mei 1971; Davey & Stewartson 1974; Dysthe 1979; Trulsen & Dysthe 1996; Trulsen et al. 2000; Li 2021) as well as the Hasselmann/Zakharov integral equation (Hasselmann 1962; Zakharov 1968; Janssen 1983; Stiassnie & Shemer 1984; Krasitskii 1994) have been widely known as a powerful analytical tool. The NLS equation-based model is especially superior to the HOS method in the sense of the computational efficiency, arising from the fact that it describes the evolution of an envelope that varies slowly in time and depends on a long length scale compared with the rapidly varying wave phase with a short length scale. A NLS equation, e.g. those by Trulsen et al. (2000), Gramstad & Trulsen (2011) and Li (2021), can efficiently resolve the phase of free waves in a large computational domain while it well captures the wave energy transfers due to quartet interactions of waves within a narrow bandwidth. However, approximations in addition to a perturbation expansion are necessary throughout the derivations of a NLS equation, having limited its wide validity (see, e.g. p.202 by Janssen 2004). This especially indicates that their capability of accounting for the physics due to moderately nonlinear and steeper waves is likely compromised.

Similar to a NLS-based equation, the reduced Zakharov equation has been widely used for elucidating the nonlinear physical properties of water waves; see, e.g. Crawford, Saffman & Yuen (1980), Janssen & Onorato (2007), Stiassnie & Gramstad (2009) and Gramstad (2014) among others. It produces explicit expressions for the interaction of a number of up to five waves and its numerical and theoretical potential has been extensively explored in recent years (Annenkov & Shrira 2001, 2006, 2009; Dyachenko & Zakharov 2011; Dyachenko, Kachulin & Zakharov 2017). One distinctive feature of the reduced Zakharov equation and its compact form is that they describe the evolution of a complex function that is defined using the Hamiltonian structure of physical variables for eliminating non-resonant interaction terms. This means other wave fields such as the surface elevation and velocity are evaluated at an additional cost. Nevertheless, existing wave-averaged equations of ocean mean flows commonly rely on the input of the wave-induced (Eulerian or Lagrangian) velocity and surface elevation (see, e.g. Sullivan & McWilliams (2010); Suzuki & Fox-Kemper (2016) and references therein). In such a context, the use of the new complex function may not necessarily introduce the advantage of the possible elimination of the nonlinear non-resonant terms for surface waves in a non-conservative system as they are expected to play a role. Thus, the Zakharov equation and its compact and numerical versions are not seemingly ideal as a wave-phase resolved coupled model in large-scale physical processes in the open ocean, either.

The superior features of a NLS equation-based model to the HOS method are especially important in the study of the roles of surface waves in the dynamics of the upper ocean, e.g. vertical mixing and the circulation of submesoscale currents, attributing to two aspects. Firstly, various important physical processes, e.g. the exchange of the momentum and energy flux between surface waves and a submesoscale flow, occur in a temporal and length scale at a magnitude significantly larger than that of surface waves, as has been especially found in wave-averaged equations for the dynamics of the upper ocean flows (McWilliams, Restrepo & Lane 2004; Sullivan & McWilliams 2010; Suzuki & Fox-Kemper 2016). This suggests the need for capturing the long-term energy evolution of waves described on an extremely large domain. Secondly, recent studies find that wave phases play an important role in distorting turbulence with a scale smaller than the surface waves (Teixeira & Belcher 2002; Thorpe et al. 2004; D'asaro 2014) and the generation of turbulence by non-breaking surface waves (Babanin 2006; Benilov 2012), addressing the additional need for resolving wave phases for important physical mechanisms with a small scale.

Indeed, a NLS equation-based model has been included in the framework of a regional ocean modelling system (ROMS) based on McWilliams et al. (2004), despite the fact that a few important nonlinear wave physics such as the Benjamin–Feir instability (Benjamin & Feir 1967; Longuet-Higgins 1978; Janssen & Herbers 2009) and quartet resonant interaction of waves (Phillips 1960) have not been considered yet. Moreover, due to the multiple scales involving a few orders of magnitude, the understanding of the coupled effects between small-scale turbulence, middle-scale surface waves and large-scale submesoscale currents have been extremely limited. To make a difference, the author believes that it relies on an accurate and efficient model of surface waves, which should have the potential of bridging the connections between the smaller and larger scale physical processes with the scale in the middle being characterized by surface waves. To this end, neither the HOS method nor a NLS equation-based model is seemingly ideal, the former of which due to the relatively low numerical efficiency for wave parameters required on an extremely large domain and the latter of which due to the restricted accuracy and validity.

Following the above discussion, an obvious question is whether it is possible to derive a framework that combines the advantages of both the HOS method and a NLS equation-based model, with the potential of being applied in more general works that directly bridge the coupled physical processes of ocean surface waves with small-scale turbulence and submesoscale current. It means that such a framework should be as accurate as the HOS method while permitting the main numerical features of a NLS equation-based model. Specifically, it is desired to include the computational efficiency that allows for a large and coarse computation domain on which both the amplitude and wave phase can be well resolved. Addressing this question defines the primary objective of this paper. It aims to present a new framework originally inspired by a NLS equation-based model in the manner that envelopes are introduced as the starting point. The coupled envelope evolution equations that can be derived accurate to arbitrary order in wave steepness are presented in the Hamiltonian theory. It should be noted that the idea of the envelope equations, which takes the advantages of both Fourier transforms and a newly defined linear operator, is different from that of the localized Zakharov equation (LZE). The idea of LZE deals with wave field dynamics in a manner of multiple interacting wave packets, posing challenges in its mathematical implementations (Rasmussen & Stiassnie 1999; Gramstad, Agnon & Stiassnie 2011). With a need for extension, the newly derived framework is expected to have especially wide applicability in terms of coupling the surface wave driven processes with regional oceanic dynamics in the upper ocean.

This paper is laid out as follows. The statement of the problem is presented in § 2, followed by a review of the HOS method (§ 3.1) and a traditional perturbation method (§ 3.2) in § 3. A new proposed framework is presented in § 4. In § 4.1 the (slowly spatial-temporal varying) envelope of both surface elevation and velocity potential on the free water surface – which are firstly introduced in this paper – are shown to be a new pair of canonical variables. The detailed derivations for the coupled envelope evolution equations (CEEEs) are presented in §§ 4.2–4.4. A few important features of the CEEEs are explored in § 5 and § 6. Three aspects of the newly derived CEEEs are discussed in § 5. How to apply an exponential integrator with the CEEEs is presented in § 5.1, together with the numerical implementation of the CEEEs that leads to the accurate description of linear waves. It is shown in § 5.2 that the CEEEs can lead to the energy balance equation. In the limiting cases where wave nonlinearity can be neglected, the energy balance equation is shown to naturally conserve the energy, implying no exchange of energy between linear waves as it should be. The nonlinear forcing terms in the CEEEs have clear physical meanings that are discussed in § 5.3, including those that can only contribute to the forcing of bound waves and those that are capable of forcing free waves arising from the quartet and quintet resonant interactions of waves. The CEEEs are compared with a traditional perturbation method in § 6.1, where the relations between the two methods are analytically shown for the evolution of a train of both Stokes waves (Fenton 1985) and three-dimensional waves with arbitrary bandwidth and with large directional spreading (Li & Li 2021). The comparisons between the CEEEs and the HOS method are discussed in § 6.2 in numerical computations for a limiting case and the numerical performances illustrated through a few numerical algorithms used for numerical implementations. The main conclusions are presented in § 7.

2. Mathematical formulation

2.1. Problem definition

We consider ocean surface waves propagating on waters of a finite depth in the framework of potential flow theory, thereby assuming incompressible inviscid flows and irrotational fluid motions, and negligible effects of surface tension. A Cartesian coordinate system is chosen with the undisturbed water surface located at $z = 0$. A list of the main symbol notations used in this work is given in table 1. 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 \text{for}\ -h< 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, $h$ is the water depth assumed to be constant, $t$ is the time and $\boldsymbol {\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) on the free water surface $z= \zeta (\boldsymbol {x},t)$, respectively,

(2.2a,b)\begin{equation} \partial_t\zeta +\boldsymbol{\nabla} \varPhi\boldsymbol{\cdot} \boldsymbol{\nabla}\zeta - \partial_z\varPhi = 0 \quad \text{and}\quad \partial_t\varPhi +g\zeta+\tfrac{1}{2}\left(\boldsymbol{\nabla}_3\varPhi\right)^2=0, \end{equation}

where $g$ denotes the gravitational acceleration; a seabed boundary condition

(2.3)\begin{equation} \partial_z \varPhi = 0\quad \text{for}\ z =-h, \end{equation}

where a constant uniform water depth $h$, is assumed. It should be noted that the extension of the new envelope equations presented in § 4 to permit a slowly varying water depth with the horizontal position would be straightforward following Dommermuth & Yue (1987) and the detailed derivations in this paper.

Table 1. Nomenclature.

2.2. Boundary conditions on the free water surface

Following Zakharov (1968) and Krasitskii (1994), we introduce the potential ($\psi$) and vertical velocity ($W$) defined on the unknown free water surface $z=\zeta (\boldsymbol {x},t)$ as

(2.4a,b)\begin{equation} \psi(\boldsymbol{x},t) = \varPhi(\boldsymbol{x}, \zeta(\boldsymbol{x},t), t) \quad \text{and}\quad W(\boldsymbol{x},t) = \partial_z\varPhi(\boldsymbol{x}, z, t), \quad \text{for}\ z= \zeta(\boldsymbol{x},t). \end{equation}

Inserting the definition of $\psi$ and $W$ into (2.2a,b) leads to the boundary conditions on the free water surface expressed as equations for unknowns $\psi (\boldsymbol {x},t)$, $W(\boldsymbol {x},t)$ and $\zeta (\boldsymbol {x},t)$, given by

(2.5a,b)\begin{align} \partial_t\zeta -W =-\boldsymbol{\nabla} \psi\boldsymbol{\cdot} \boldsymbol{\nabla}\zeta +W(\boldsymbol{\nabla}\zeta)^2\quad \text{and}\quad \partial_t\psi + g\zeta =- \dfrac{1}{2}(\boldsymbol{\nabla}\psi)^2+\dfrac{1}{2} W ^2\left[1+(\boldsymbol{\nabla}\zeta)^2\right]. \end{align}

The system described by the equations composed of (2.1), (2.3) and (2.5a,b) is known as the fully nonlinear (potential flow) boundary value problem in a Hamiltonian theory (see, e.g. West et al. 1987). An approximate solution to this problem can be obtained by using various methods as noted in the introduction, e.g. a HOS method (Dommermuth & Yue 1987; West et al. 1987), Hasselmmann/Zakharov integral equation (Hasselmann 1962; Zakharov 1968; Krasitskii 1994) and the CEEEs that are derived for the first time in § 4.

2.3. Velocity potential

We seek the solution for the unknown potential ($\varPhi (\boldsymbol {x},z,t)$) of the Laplace equation and the seabed boundary condition given by (2.1) and (2.3), respectively, in the form of a power series in wave steepness denoted by $\epsilon$ that stands for a small non-dimensional scaling parameter

(2.6)\begin{equation} \varPhi(\boldsymbol{x},z,t) = \sum_{m=1}^{M} \epsilon^m\varPhi^{(m)}(\boldsymbol{x},z,t), \end{equation}

where the terms are kept up to the $M$th order in wave steepness and the superscript ‘($m$)’ denotes $O(\epsilon ^m)$, and the unknown (real) potential at the $m$th order in wave steepness is given by (see, e.g. Dommermuth & Yue 1987)

(2.7)\begin{equation} \varPhi^{(m)}(\boldsymbol{x},z,t) = \int_{-\infty}^{\infty} \hat{\varPhi}_0^{(m)}(\boldsymbol{k},t) \dfrac{\cosh |\boldsymbol{k}|(z+h)}{\cosh |\boldsymbol{k}|h}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{k}, \end{equation}

where $\hat {\varPhi }^{(m)}_0(\boldsymbol {k},t)=\hat {\varPhi }^{(m)}(\boldsymbol {k},0,t)$ denotes the $m$th order velocity potential evaluated at $z=0$ in the Fourier $\boldsymbol {k}$ plane; the subscript ’0’ is used to denote the evaluation at a still water surface $z=0$.

2.4. Definition of two operators

For later reference and simplicity, we introduce a characteristic wave vector and angular wave frequency, denoted by $\boldsymbol {k}_0=(k_0,0)$ and $\omega _0$, respectively, with the magnitude $k_0=|\boldsymbol {k}_0|$. It is highlighted that the positive $\boldsymbol{x}$ direction is chosen along the direction of the wave vector $\boldsymbol {k}_0$, which also accords to the main direction of wave propagation. The characteristic wave vector and frequency obey the linear dispersion relation $\omega _0=\omega (\boldsymbol {k}_0,h)$ with $\omega (\boldsymbol {k},h) =\sqrt {g|\boldsymbol {k}|\tanh |\boldsymbol {k}|h}$.

We introduce two operators that are demonstrated by an arbitrary temporal-spatial function, $\chi (\boldsymbol {x},t)$, including a Fourier transform with respect to the horizontal position vector $\boldsymbol {x}$ and a new operator referred to as the envelope transform, given by, respectively,

(2.8a)$$\begin{gather} \hat{\chi}(\boldsymbol{k},t) = \dfrac{1}{4{\rm \pi}^2}\int_{-\infty}^{\infty} \chi(\boldsymbol{x},t)\mathrm{e}^{-\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d}\kern0.06em \boldsymbol{x}, \end{gather}$$

(2.8b)$$\begin{gather}\chi^{[j]}_+(\boldsymbol{x},t; \alpha, \beta) = \mathrm{e}^{\mathrm{i} j\beta \omega_0t}\int_{-\infty}^{\infty} 2\varTheta[(\boldsymbol{k}+j\alpha \boldsymbol{k}_0)\boldsymbol{\cdot}\boldsymbol{k}_0]\hat{\chi}(\boldsymbol{k}+\alpha\boldsymbol{k}_0,t) \mathrm{e}^{\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}} \,\mathrm{d}\boldsymbol{k}, \end{gather}$$

in which the hat denotes a parameter transformed to the Fourier $\boldsymbol {k}$ plane; $\alpha$, $\beta$ and $j$ are arbitrary non-negative constants that can be freely chosen to facilitate the numerical computations as will be shown in the following; a combination of the subscript ‘+’ and superscript ‘[$j$]’ are used to denote the envelop transform; $\varTheta$ denotes the Heaviside step function. The relationship between the two operators defined in (2.8b) is shown in figure 1. It is seen from figure 1 that the envelope transform is composed of three consecutive procedures. Firstly, it collects the components of $2\hat {\chi }(\boldsymbol {k})$ in the wavenumber region for $\boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {k}>0$. Subsequently, it employs the translation operator defined as $\exp [-\mathrm {i} j (\alpha \boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\beta \omega _0 t)]$ in the Fourier plane. Thirdly, it operates an inverse Fourier transform given by (2.8a). For $\alpha =0=\beta$ or $j=0$, it is understood that the envelope transform simply recovers the inverse Fourier transform with respect to $2\hat {\chi }$ in the positive wavenumber plane where $\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {k}_0>0$.

Figure 1. Diagram of the operators in the Fourier wavenumber space in two dimensions; for the envelope transform, (a) $\alpha = 1$ and (b) $\alpha =1.5$.

Due to the symmetrical properties of a Fourier transform and the definition of the envelope transform, we readily obtain

(2.9)\begin{equation} \chi(\boldsymbol{x},t) = \tfrac{1}{2} \chi^{[j]}_+\exp\left({\mathrm{i} j (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \text{c.c.}, \end{equation}

where $\text {c.c.}$ denotes the complex conjugates, and

(2.10)\begin{equation} \hat{\chi}^{[j]}_+(\boldsymbol{k},t) = 2\varTheta[(\boldsymbol{k}+j\alpha \boldsymbol{k}_0)\boldsymbol{\cdot}\boldsymbol{k}_0]\hat{\chi}(\boldsymbol{k}+j\alpha\boldsymbol{k}_0,t)\mathrm{e}^{\mathrm{i} j\beta\omega_0t}. \end{equation}

It is understood that $\hat {\chi }_+^{[j]}$ varies with the different combinations of the constants $(\alpha, \beta, j)$, as shown in figure 1 where two special cases with $\alpha =1$ and $\alpha =1.5$ are shown. With $j=1$, the envelope transform leads to the definition of the envelope of the wave elevation and potential at the free water surface, which will be used in § 4,

(2.11a,b)\begin{equation} A(\boldsymbol{x},t; \alpha,\beta) = \zeta^{[1]}_+\quad \text{and}\quad B_{s}(\boldsymbol{x},t; \alpha,\beta) = \psi^{[1]}_+, \end{equation}

thereby

(2.12a,b)\begin{equation} \zeta = \tfrac{1}{2} A \exp\left({\mathrm{i} (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \text{c.c.} \quad \text{and}\quad \psi = \tfrac{1}{2} B_s \exp\left({\mathrm{i} (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \text{c.c.} \end{equation}

The main intention of allowing for an arbitrary choice of ($\alpha,\beta,j$) is to facilitate numerical implementations and, therefore, improve the numerical efficiency. For example, with the choice of $\alpha =0$ and $\beta =1$, we will show in § 6 that the computational efficiency of a HOS method can be improved. With $\alpha =\beta =1$ and additional assumptions required, the CEEEs would be reduced to a third-order accurate NLS equation-based model that has been demonstrated with an excellent performance.

3. A review of two methods

In this section we review two methods for the description of non-breaking surface waves, which are the HOS method presented in § 3.1 and a so-called traditional perturbation method in § 3.2. Both methods rely on an unknown velocity potential being expressed in the form of a perturbation expansion and Fourier transform, as explained in § 2.3. Their distinctive difference lies in the fact that they seek different approaches for the unknowns: surface elevation, potential and vertical velocity on the free water surface.

3.1. The HOS method

The HOS method proposes to solve the fully nonlinear (potential flow) boundary value problem in a Hamiltonian theory, as introduced in § 2.2, for two main unknowns that are the surface elevation ($\zeta$) and potential ($\psi$) on the free water surface. It consists of two procedures. It firstly seeks to express the unknowns including the velocity potential ($\varPhi$) and the vertical velocity ($W(\boldsymbol {x},t)$) on the free water surface in the form of a function of $\zeta$ and $\psi$, which will be presented in § 3.1.1. Secondly, through using the boundary conditions (2.5a,b) it leads to the evolution equations that can be numerically solved for the two main unknowns, as presented in § 3.1.2.

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 (1987) and West et al. (1987), 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. 1987). The nonlinear terms on the right-hand side of (3.5a,b) are given by (3.7) and (3.8).

3.2. Traditional perturbation method

Note that the derivations presented in this section will only be used in § 6.1 for comparison and completeness. Different from the methods based on the boundary conditions given by (2.5a,b) in § 2.2, reference is made to one of these methods that have the following features as a traditional perturbation method. It primarily seeks the approximate solution to the boundary value problem described by (2.1), (2.2a,b) and (2.3). The boundary conditions (2.2a,b) on the free water surface are especially expanded about the still water surface $z=0$ for all $z$-dependent parameters, in contrast to the boundary conditions given by (3.5a,b) where only the vertical velocity is expanded about $z=0$. The primary unknowns are the velocity potential on a still water surface and the surface elevation that are expressed in the form of a power series in wave steepness. Substituting these approximate expressions into the boundary value problem given by (2.1), (2.2a,b) expanded about $z=0$, and (2.3), and collecting the same orders in wave steepness will lead to the boundary value problems at different order in wave steepness. These boundary value problems are solved in sequence from the first to the $M$th order in wave steepness. A few examples that are based on a traditional perturbation method are Chu & Mei (1971), § 13 by Mei, Stiassnie & Yue (2005); a NLS equation-based model like Davey & Stewartson (1974), Dysthe (1979), Trulsen et al. (2000), Slunyaev (2005) and Li (2021); the fifth-order Stokes waves by Fenton (1985) and the second-order broadband framework by Li & Li (2021).

The main derivations needed in a traditional perturbation method are explained in the following, for which the leading-order approximations are kept only to the second order in wave steepness for simplicity. In order to indicate the differences with the main results presented in most chapters of this paper, a prime is added to denote the parameters used in a traditional expansion method and subscripts are used to denote the different orders in wave steepness and wave harmonics. An approximate form for both the unknown velocity potential and elevation are assumed, to the second order in wave steepness $\epsilon _0$,

(3.9a,b)\begin{equation} \varPhi = \epsilon_0 \varPhi'_{11}+ \epsilon^2_0 \underbrace{(\varPhi'_{22} + \varPhi'_{20})}_{{\equiv} \varPhi'_2(\boldsymbol{x},z,t)} \quad \text{and}\quad \zeta= \epsilon_0 \zeta'_{11}+ \epsilon^2_0 \underbrace{(\zeta'_{22} + \zeta'_{20})}_{{\equiv} \zeta'_2(\boldsymbol{x},t)}, \end{equation}

where $\epsilon _0$ denotes the dimensionless steepness of linear waves that obeys $O(\epsilon _0)\sim O(k_0\zeta' _{11})$ to primarily distinguish it from $\epsilon$ defined in § 2.3; the subscript ‘$mj$’ denotes $O(\epsilon _0^m)$ and the $j$th wave harmonic; $\zeta '_{mj}=\zeta '_{mj}(\boldsymbol {x},t)$ and $\varPhi _{mj}'= \varPhi '_{mj}(\boldsymbol {x},z,t)$. The potential and vertical velocity on the free water surface are given in the form of a Taylor expansion about $z=0$, to the second order in wave steepness $\epsilon _0$,

(3.10a,b)\begin{equation} \psi = \epsilon_0 \varPhi'_{11} +\epsilon^2_0 \zeta'_{11}\partial_z\varPhi'_{11} \quad \text{and}\quad W= \epsilon_0 \partial_z\varPhi'_{11} + \epsilon^2_0\zeta'_{11}\partial_{zz}\varPhi'_{11}, \quad \text{for}\ z=0. \end{equation}

Inserting (3.9a,b) for the unknown potential and elevation, respectively, into the surface boundary conditions (2.2a,b), expanding the equations about $z=0$, collecting the terms at second order in $\epsilon _0$ leads to

(3.11a)$$\begin{gather} \partial_t \zeta'_{2} -\partial_z\varPhi'_{2} = \zeta'_{11}\partial_{zz}\varPhi'_{11} - \boldsymbol{\nabla}\zeta'_{11}\boldsymbol{\cdot}\boldsymbol{\nabla}\varPhi'_{11}, \end{gather}$$

(3.11b)$$\begin{gather}\partial_t\varPhi'_{2} + g\zeta'_{2} = \zeta'_{11}\partial_{tz}\varPhi'_{11}-\tfrac{1}{2} (\boldsymbol{\nabla}_3\varPhi'_{11})^2, \end{gather}$$

which are used to solve for the unknowns (i.e. $\zeta _{2}$ and $\varPhi _2$) at second order with the linear parameters obtained from the linearized equations of (2.1), (2.2a,b) and (2.3); see, e.g. § 13 by Mei et al. (2005). With the forcing terms on the right-hand side of (3.11) being separated according to the wave harmonics, the unknown fields with the subscript ‘$mj=22$’ and ‘$mj=20$’ can be obtained due to the second-order superharmonic and subharmonic waves, respectively (Li et al. 2021).

A different but equivalent framework to (3.11) has been proposed by Li & Li (2021) where envelopes have been introduced, which are the primary unknowns for the waves of different harmonics up to the second order in wave steepness, $\epsilon _0$. The elevation and potential are in particular given by

(3.12a) $$\begin{align} \zeta &= \tfrac{1}{2} \epsilon_0 A_{11}'\exp({\mathrm{i}(\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\mathrm{i}\omega_0t)}) +\text{c.c.}\notag\\ &\quad + \epsilon^2_0 \left( \tfrac{1}{2} A'_{22}\exp({2\mathrm{i}(\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\mathrm{i}\omega_0t)}) +\tfrac{1}{2} A'_{20} +\text{c.c.} \right) , \end{align}$$

(3.12b) $$\begin{align}\varPhi &= \tfrac{1}{2} \epsilon_0 B_{11}'\exp({\mathrm{i}(\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\mathrm{i}\omega_0t)}) +\text{c.c.}\notag\\ &\quad + \epsilon^2_0 \left( \tfrac{1}{2} B'_{22}\exp({2\mathrm{i}(\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\mathrm{i}\omega_0t)})+\tfrac{1}{2} B'_{20} +\text{c.c.} \right) , \end{align}$$

where the (complex) envelopes $A'_{mj}(\boldsymbol {x},t)$ and $B'_{mj}(\boldsymbol {x},z,t)$ are the main unknowns. With the second-order elevation and potential, the envelopes obey

(3.12c)\begin{equation} A'_{mj}(\boldsymbol{x},t)=\left[\zeta'_{mj}\right]_+^{[j]} \quad \text{and}\quad B'_{mj}(\boldsymbol{x},z,t)=\left[\varPhi'_{mj}\right]_+^{[j]}, \end{equation}

for ‘$mj=11$’, ‘$mj=20$’ and ‘$mj=22$’. Due to the Laplace equation and the seabed boundary condition, the vertical structure of the envelope $B'_{mj}$ is obtained and given by

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

where $\hat {B}'_{mj}(\boldsymbol {k},t)$ denotes the envelope $B'_{mj}(\boldsymbol {x},z,t)$ at $z=0$ transformed to the Fourier space. Similarly, inserting (3.12a,b) into (2.2a,b), expanding the $z$-dependent wave parameters about $z=0$, collecting the terms at the second order and separating the wave harmonics leads to the boundary conditions for the second-order superharmonic waves on a still water surface, $z=0$,

(3.13a)$$\begin{gather} (\partial_{t}-2\mathrm{i}\omega_0)A'_{22}-\partial_{z}B'_{22} = \tfrac{1}{2} A'_{11}\partial_{zz}B'_{11} - \tfrac{1}{2} (\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)A'_{11}\boldsymbol{\cdot}(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)B'_{11}, \end{gather}$$

(3.13b)$$\begin{gather}(\partial_{t}-2\mathrm{i}\omega_0)B'_{22}+gA'_{22} = \tfrac{1}{2} A'_{11}\partial_{tz}B'_{11} - \tfrac{1}{2} [(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)B'_{11}]^2 - \tfrac{1}{2} (\partial_zB'_{11})^2; \end{gather}$$

the boundary conditions for the second-order subharmonic waves on a still water surface are

(3.14a)$$\begin{gather} \partial_tA'_{20}-\partial_{z}B'_{20} = \left[ \mathcal{R}\left( \tfrac{1}{2} A^{'*}_{11}\partial_{zz}B'_{22} - \tfrac{1}{2} (\boldsymbol{\nabla}-\mathrm{i}\boldsymbol{k}_0)A^{'*}_{11}\boldsymbol{\cdot}(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)B'_{11}\right) \right]^{(0)}_{+}, \end{gather}$$

(3.14b)$$\begin{gather}\partial_t B'_{20}+gA'_{20} = \left[ \mathcal{R}\left( \tfrac{1}{2} A^{'*}_{11}\partial_{tz}B'_{22} - \tfrac{1}{2} |(\boldsymbol{\nabla}+\mathrm{i}\boldsymbol{k}_0)B'_{11}|^2 - \tfrac{1}{2} |\partial_zB'_{11}|^2\right) \right]^{(0)}_{+}, \end{gather}$$

where $\mathcal {R}$ denotes the real component. With the linear envelope $A'_{11}$ and $B'_{11}$ solved from the linearized problem, the second-order envelopes can be obtained from (3.13) and (3.14) for the superharmonic and subharmonic waves, respectively. Both the envelope equations, (3.13) and (3.14), for the envelopes as well as the equations, (3.11), for the second-order elevation and potential will be used in § 6.1 for a comparison with the CEEEs derived in § 4, showing how the CEEEs can recover to these equations if waves are considered only up to the second order in wave steepness, $\epsilon _0$.

4. Coupled envelope evolution equations in a Hamiltonian theory

In this section we derive the new equations, referred to as the CEEEs, in a Hamiltonian theory for numerical solutions, which are the main results of this paper. To this end, the starting point is a new pair of canonical variables shown in § 4.1. In contrast to the HOS, the new canonical variables become our primary unknowns that can be numerically solved for. In particular, the wave parameters (velocity potential and vertical velocity) at different orders in wave steepness will be obtained through harmonic separations, as shown in §§ 4.2 and 4.3. The CEEEs describing the evolution of the new canonical variables are presented in § 4.4. Similar to the HOS method, the CEEEs can be derived up to an arbitrary order in wave steepness, with the general expressions presented in Appendix A.

4.1. A new pair of canonical variables

It is understood that $\zeta$ and $\psi$ is a pair of canonical variables (Zakharov 1968; Krasitskii 1994). Due to the definition of the envelopes, we proceed to show that $A$ and $B_s$ are a new pair of canonical variables. It is understood that, due to the properties of a Fourier transform, the following identities hold:

(4.1a–d)\begin{align} \hat{\zeta}(\boldsymbol{k}) = \hat{\zeta}^*(-\boldsymbol{k}),\quad \hat{\zeta}(-\boldsymbol{k}) = \hat{\zeta}^*(\boldsymbol{k}),\quad \hat{\psi}(\boldsymbol{k}) = \hat{\psi}^*(-\boldsymbol{k}),\quad \text{and}\quad \hat{\psi}(-\boldsymbol{k}) = \hat{\psi}^*(\boldsymbol{k}). \end{align}

Here the asterisk ‘*’ denotes the complex conjugates. Let $H'$ be the Hamiltonian for $\hat {\zeta }$ and $\hat {\psi }$, suggesting that

(4.2a,b)\begin{equation} \partial_t\hat{\zeta} = \dfrac{\delta H'}{\delta \hat{\psi}^*} \quad \text{and} \quad \partial_t\hat{\psi}=- \dfrac{\delta H'}{\delta \hat{\zeta}^*}, \end{equation}

where the Hamiltonian $H'=H'(\hat {\zeta },\hat {\zeta }^*,\hat {\psi }, \hat {\psi }^*)$ and $\delta$ denotes the functional derivative. By definition we obtain $\hat {\zeta }(\boldsymbol {k}+\alpha \boldsymbol {k}_0) = \hat {A}(\boldsymbol {k},t)\exp {(-\mathrm {i}\beta \omega _0t)}$ and $\hat {\psi }(\boldsymbol {k}+\alpha \boldsymbol {k}_0) = \hat {B}(\boldsymbol {k},t)\exp {(-\mathrm {i}\beta \omega _0t)}$, and by inserting these into the Hamiltonian $H'$, we obtain

(4.3)\begin{equation} H'(\hat{\zeta},\hat{\psi},\hat{\zeta}^*, \hat{\psi}^*) \equiv H(\hat{A},\hat{B}_s,\hat{A}^*,\hat{B}_s^*), \end{equation}

where $H$ denotes the new Hamiltonian obtained through replacing the elevation and potential with their envelopes in $H'$. Next, we perform the following functional derivatives based on (4.3) and obtain

(4.4a,b)\begin{equation} \delta_{\hat{\psi}^*} H' = \exp{(-\mathrm{i}\beta\omega_0t)} \delta_{B^*} H \quad \text{and}\quad \delta_{\hat{\zeta}^*} H' = \exp{(-\mathrm{i}\beta\omega_0t)} \delta_{A^*} H. \end{equation}

Inserting (4.4a,b) into the right-hand sides of (4.2a,b), replacing the elevation and potential with their envelopes on the left-hand sides of (4.2a,b) and eliminating the factor $\exp {(-\mathrm {i}\beta \omega _0t)}$, we obtain

(4.5a,b)\begin{equation} \partial_t \hat{A} -\mathrm{i}\omega_0 \hat{A} = \delta_{B_s^*} H \quad \text{and}\quad -\partial_t \hat{B}_s +\mathrm{i}\omega_0 \hat{B}_s = \delta_{A^*} H. \end{equation}

Multiplying (4.5a,b) by $\hat {B}_s^*$ and $\hat {A}^*$, respectively, leads to

(4.6a,b)\begin{equation} \hat{B}^*\partial_t \hat{A} -\mathrm{i}\omega_0 \hat{B}_s^*\hat{A} = \hat{B}_s^*\delta_{B^*_s} H' \quad \text{and}\quad - \hat{A} ^*\partial_t \hat{B}_s +\mathrm{i}\omega_0 \hat{A} ^* \hat{B}_s =\hat{A} ^*\delta_{A^*} H', \end{equation}

based on which we introduce a new Hamiltonian defined as

(4.7)\begin{equation} H_{AB}(\hat{A},\hat{B}_s,\hat{A}^*,\hat{B}_s^*)= \int H' + \mathrm{i} \omega_0 \left[ \hat{B}_s^*\hat{A} -\hat{A} ^* \hat{B}_s\right]\,\mathrm{d} \boldsymbol{k}. \end{equation}

Performing the following functional derivatives on the new Hamiltonian $H_{AB}$ leads to

(4.8a,b)\begin{equation} \delta_{B_s^*} H_{AB} = \delta_{\hat{B}_s^*} H' + \mathrm{i} \omega_0 \hat{A} \quad \text{and}\quad \delta_{A^*} H_{AB} = \delta_{\hat{A}^*} H' - \mathrm{i} \omega_0 \hat{B}_s. \end{equation}

Inserting (4.8a,b) for $\delta _{B^*} H'$ and $\delta _{A^*} H'$ into the right-hand sides of (4.5a,b), respectively, leads to

(4.9a,b)\begin{equation} \partial_t\hat{A} = \dfrac{\delta H_{AB}'}{\delta \hat{B}_s^*} \quad \text{and}\quad \partial_t\hat{B}_s =-\dfrac{\delta H_{AB}'}{\delta \hat{A}^*}, \end{equation}

meaning that $\hat {A}$ and $\hat {B}_s$ are a pair of canonical variables. As noted in Krasitskii (1994), due to the fact that the inverse Fourier transform is a canonical one, (4.9a,b) therefore also imply that envelope $A$ and $B_s$ are a pair of canonical variables.

4.2. Three different methods for the evaluation of a quadratic term

In this paper we take advantage of the symmetrical properties of a Fourier transform. The separation of wave harmonics presented in § 4.3 builds upon two key features. Firstly, it is understood that nonlinear terms (e.g. $\mathcal {W}^{(3)}$ and $\mathcal {T}^{(4)}$) at orders higher than the second can always be written in a form of the linear superposition of the product of two parameters. We use $\mathcal {T}^{(4)}$ described by (3.8c) as an example. We define

(4.10a,b)\begin{equation} W_{sq}^{(1)} = (W^{(2)})^2 \quad \text{and}\quad \zeta_x^{(2)} = (\boldsymbol{\nabla}\zeta)^2. \end{equation}

Inserting (4.10a,b) into (3.8c) leads to

(4.11)\begin{equation} \mathcal{T}_0^{(4)}= \tfrac{1}{2} (W^{(2)})^2+ W^{(1)} W^{(3)}+ \tfrac{1}{2} W_{sq}^{(1)} \zeta_x^{(2)}, \end{equation}

which is in a form of the linear superposition of quadratic terms but it is obvious that it is not at second order in wave steepness. Second, by virtue of the symmetrical properties of a Fourier transform, we next consider a function of two arbitrary real parameters $\chi (\boldsymbol {x})$ and $\xi (\boldsymbol {x})$ defined as $\boldsymbol {f}(\boldsymbol {x}) = \boldsymbol {\nabla }\chi (\boldsymbol {x})\xi (\boldsymbol {x})$. Assuming the Fourier transform of both $\chi (\boldsymbol {x})$ and $\xi (\boldsymbol {x})$ exist, $\boldsymbol {f}(\boldsymbol {x})$ can also be expressed in the form of an inverse Fourier transform

(4.12a)$$\begin{gather} \boldsymbol{f}(\boldsymbol{x})= \int_{-\infty}^{\infty} \mathrm{i}\boldsymbol{k}_1\hat{\chi}_1\hat{\xi}_2\mathrm{e}^{\mathrm{i}(\boldsymbol{k}_1+\boldsymbol{k}_2)\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d} \boldsymbol{k}_1\,\mathrm{d} \boldsymbol{k}_2 \quad \text{or} \end{gather}$$

(4.12b)$$\begin{gather}\boldsymbol{f}(\boldsymbol{x}) = \left[ \int_{-\infty}^{\infty} \mathrm{i}\boldsymbol{k}_1 \hat{\chi}_1\mathrm{e}^{\mathrm{i}\boldsymbol{k}_2\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d} \boldsymbol{k}_1 \right]\times \left[ \int_{-\infty}^{\infty} \hat{\xi}_2\mathrm{e}^{\mathrm{i}\boldsymbol{k}_2\boldsymbol{\cdot}\boldsymbol{x}}\,\mathrm{d} \boldsymbol{k}_2 \right], \end{gather}$$

where $\hat {\chi }_1 = \hat {\chi }(\boldsymbol {k}_1)$ and $\hat {\xi }_2=\hat {\xi }(\boldsymbol {k}_2)$. Due to the symmetrical property of a Fourier transform (i.e. $\hat {\chi }(-\boldsymbol {k})= \hat {\chi }^* (\boldsymbol {k})$) and decomposing the entire integral region of the integrals in (4.12a) into four equal quarters leads to

(4.13)\begin{align} &\boldsymbol{f}(\boldsymbol{x})= \left[ \exp\left({2\mathrm{i} (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) \int_{\varGamma_{1+}} \int_{\varGamma_{2+}} \mathrm{i}\boldsymbol{k}_1\hat{\chi}_1\hat{\xi}_2\right.\nonumber\\ &\quad\qquad\left. \exp\left({\mathrm{i}(\boldsymbol{k}_1-\alpha\boldsymbol{k}_0+\boldsymbol{k}_2-\alpha\boldsymbol{k}_0)\boldsymbol{\cdot}\boldsymbol{x} +2\mathrm{i}\beta\omega_0t}\right)\,\mathrm{d}\boldsymbol{k}_1\,\mathrm{d} \boldsymbol{k}_2 + \text{c.c.} \vphantom{\int_{\varGamma_{1+}}}\right]\nonumber\\ &\qquad + \left[ \int_{\varGamma_{1+}} \int_{\varGamma_{2+}} \mathrm{i}\boldsymbol{k}_1 \hat{\chi}_1\hat{\xi}^*_2 \exp\left({\mathrm{i}(\boldsymbol{k}_1-\alpha\boldsymbol{k}_0)\boldsymbol{\cdot}\boldsymbol{x} -\mathrm{i} (\boldsymbol{k}_2-\alpha\boldsymbol{k}_0)\boldsymbol{\cdot}\boldsymbol{x}} \right)\,\mathrm{d}\boldsymbol{k}_1\,\mathrm{d} \boldsymbol{k}_2 + \text{c.c.} \right], \end{align}

in which $\varGamma _{j+}$ for $j=1$ and $j=2$ defines the region where $\boldsymbol {k}_j\boldsymbol {\cdot }\boldsymbol {k}_0> 0$. Replacing the terms that correspond to the definition of the envelope transform in (4.13) leads to

(4.14)$$\begin{gather} \boldsymbol{f}(\boldsymbol{x})= \left[\dfrac{1}{4} (\boldsymbol{\nabla}+\mathrm{i}\alpha\boldsymbol{k}_0) \chi^{[1]}_+{\xi}^{[1]}_+ \exp\left({2\mathrm{i} (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) + \text{c.c.} \right]\nonumber\\ + \left[ \dfrac{1}{4} (\boldsymbol{\nabla}+\mathrm{i}\alpha\boldsymbol{k}_0)\chi^{[1]}_+ \left(\xi^{[1]}_+\right)^*+ \text{c.c.}\right]. \end{gather}$$

The above discussion suggests that the evaluation of $\boldsymbol {f}(\boldsymbol {x})$ admits at least three different forms, which lead to the main differences between different methods for the description of water waves as explained in the following. The HOS method relies on (4.12b) that is typically used when the derivatives with respect to $x$ or $y$ of a parameter are involved based on a pseudo-spectral method. The Hasselmann/Zakharov integral equation (Hasselmann 1962; Zakharov 1968; Stiassnie & Shemer 1984; Krasitskii 1994) uses (4.12a) for the evaluation of $\boldsymbol {f}(\boldsymbol {x})$. The CEEEs proposed in this paper rely on (4.14) instead, which is in principle the so-called separation of wave harmonics and does not rely on the narrowband assumption. In contrast to the HOS and Zakharov equations, we will obtain the equations for unknown envelopes in the following sections for the CEEEs in a manner similar to a NLS-based model but at no cost of the accuracy.

4.3. Separation of wave harmonics

We proceed to seek a different but equating expression for the description of wave fields, especially the potential and vertical velocity, in the form of functions of unknown envelope $A(\boldsymbol {x},t)$ and $B_s(\boldsymbol {x},t)$ through the separation of harmonics presented in § 4.2. Doing so will permit us to take advantage of a pseudo-spectral Fourier method in a more coarse and larger grid at a little additional cost to the numerical computations but not at the expense of the accuracy, as compared with the HOS method. Based on the solution structure for both the potential and velocity presented in § 3.1.1, we start with § 4.3.1 for the velocity potential ($\varPhi ^{(m)}$) in the form of a function of $A$ and $B_s$ up to second order in wave steepness. For simplicity, an example of $M$ up to $M=4$ for both the potential and vertical velocity is presented in § 4.3.2 with the forcing terms derived in § 4.3.3. The general procedures for the derivations up to an arbitrary order are presented in Appendix A. The new framework presented in this section can be made numerically feasible with a varying parameter of $M$, similar to the HOS method.

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.

4.4. The CEEEs

The $M$th order accurate (in wave steepness) CEEEs are obtained through the following sequential procedures: (i) inserting wave parameters and forcing terms in a form of the separation of wave harmonics presented in § 4.3 into the evolution equations (3.5a,b), (ii) keeping the components in the Fourier wavenumber region where $\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {k}_0> 0$, and (iii) multiplying all terms by a factor of $\exp (-\mathrm {i}\alpha \boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}+\mathrm {i}\beta \omega _0t)$. Hence, the CEEEs obtained are

(4.37a,b)\begin{equation} ( \partial_t-\mathrm{i}\beta\omega_0) A - \partial_zB = \mathcal{N}_{A,M} \quad \text{and}\quad ( \partial_t-\mathrm{i}\beta\omega_0) B_s + g A=~ \mathcal{N}_{B,M}, \end{equation}

with the terms for the complex conjugates removed, and the nonlinear forcing terms on the right-hand side of the equations given by

(4.38a)$$\begin{gather} \mathcal{N}_{A,M}(\boldsymbol{x},t) = \sum_{m=1}^{m=M} \mathcal{N}_A^{(m) } \equiv \sum_{m=1}^{m=M} \sum_{j=0}^{j=m} \mathcal{N}_{A}^{(mj)} \exp({\mathrm{i}(\,j-1)(\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta \omega_0t)}) \quad \text{and} \end{gather}$$

(4.38b)$$\begin{gather}\mathcal{N}_{B,M}(\boldsymbol{x},t) = \sum_{m=1}^{m=M} \mathcal{N}_B^{m) } \equiv \sum_{m=1}^{m=M} \sum_{j=0}^{j=m} \mathcal{N}_{B}^{(ij)} \exp({\mathrm{i}(\,j-1)(\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta \omega_0t)}), \end{gather}$$

where $\mathcal {N}^{(mj)}_A =0$ and $\mathcal {N}^{(ij)}_B=0$ for $ij=10$ and $ij=11$, and the other non-vanishing components are expressed as

(4.39a,b)\begin{equation} \mathcal{N}_{A}^{(mj)} = \left[ \mathcal{W}^{(mj)}\right]^{[j]}_+ \quad \text{and}\quad \mathcal{N}_{B}^{(mj)} = \left[ \mathcal{T}^{(mj)}\right]^{[j]}_+, \end{equation}

where, with $m=j$, the relations

(4.40a,b)\begin{equation} \mathcal{N}_{A}^{(mm)} = \bar{\mathcal{W}}^{(mm)} \quad \text{and}\quad \mathcal{N}_{B}^{(mm)} = \bar{\mathcal{T}}^{(mm)} \end{equation}

hold, attributing to the fact that both $\bar {\mathcal {W}}^{(mj)}$ and $\bar {\mathcal {T}}^{(mj)}$ with $m=j$ are the product of the terms that are only non-zero in the half wavenumber plane where $(\boldsymbol {k}+j\boldsymbol {k}_0)\boldsymbol {\cdot }\boldsymbol {k}_0>0$ by definition, e.g. $\hat {B}^{(22)}(\boldsymbol {k},t)$ and $\hat {B}^{(33)}(\boldsymbol {k},t)$. We will see in § 5.3.2 that the waves forced by the nonlinear terms with $m=j$ can not be free. The exploration of the newly derived CEEEs given by (4.37a,b) in this paper are presented in §§ 5 and 6.

5. Discussion of the CEEEs

The CEEEs given by (4.37a,b) have a few key features which cannot be fully explored in this paper. In this section, three aspects are especially highlighted. The first is associated with their potential of high numerical efficiency through using an exponential integrator, as examined in § 5.1 and also § 6.2.2. Secondly, due to that the two main unknowns in the CEEEs are a pair of canonical variables, it is shown in § 5.2 that the CEEEs can lead to the nonlinear evolution equation of the wave action. The third illustrates the clear physical meanings of the nonlinear forcing terms (i.e. $\mathcal {N}_{A}^{(mj)}$ and $\mathcal {N}_{B}^{(mj)}$ ) of different harmonics in § 5.3 in their capability in the nonlinear forcing of different waves.

5.1. Analytical solution and numerical implementation using an exponential integrator

For the solution of the CEEEs, there are many applicable time integration methods. In this paper we propose to using an exponential integrator; see, e.g. Hochbruck & Ostermann (2010) for details. This choice is made due to two aspects: (i) the terms involved in the CEEEs can be highly oscillatory, and (ii) we can easily identify the terms with a highly oscillatory nature from those that are slowly varying. Preforming a Fourier transform on both sides of the CEEEs gives rise to

(5.1) \begin{equation} \left[\begin{array}{@{}c@{}} \dot{\hat{A}} \\ \dot{\hat{B}}_s \end{array} \right] = \left[ \begin{array}{@{}cc@{}} \mathrm{i}\beta\omega_0 & |\boldsymbol{k}+\boldsymbol{k}_0|\tanh |\boldsymbol{k}+\boldsymbol{k}_0|h \\ - g & \mathrm{i}\beta\omega_0 \end{array} \right] \left[ \begin{array}{@{}c@{}} \hat{A} \\ \hat{B}_s \end{array} \right] + \left[ \begin{array}{@{}c@{}} \widehat{\mathcal{N} }_{A,M}(\boldsymbol{k},\tau)\\ \widehat{\mathcal{N} }_{B,M}(\boldsymbol{k},\tau) \end{array} \right], \end{equation}

in which the dot denotes the derivative with respect to time. Following an exponential integrator, the analytical solution of (5.1) can be expressed in the form

(5.2a) \begin{equation} \left[ \begin{array}{@{}c@{}} {\hat{A}}(\boldsymbol{k},t) \\ {\hat{B}}_s(\boldsymbol{k},t) \end{array} \right] = \pmb{\mathcal{E}}((t-t_0)\varOmega) \left[ \begin{array}{@{}c@{}} {\hat{A}}(\boldsymbol{k},t_0) \\ {\hat{B}}_s(\boldsymbol{k},t_0) \end{array} \right] + \sum_{m=1}^{m=M} \sum_{j=0}^{j=m} \left[ \begin{array}{@{}c@{}} \mathcal{I}^{(mj)}_A(\boldsymbol{k},t)\\ \mathcal{I}^{(mj)}_B(\boldsymbol{k},t) \end{array} \right], \end{equation}

where the initial value problem is considered at the initial time instant $t_0$ when the envelopes are given, $\pmb {\mathcal {E}}$ denotes a matrix exponential given by

(5.2b)\begin{align} \pmb{\mathcal{E}}((t-t_0)\varOmega) = \exp({\mathrm{i}\beta\omega_0(t-t_0)}) \left[ \begin{array}{@{}cc@{}} \cos \left( (t-t_0)\varOmega \right) & \dfrac{\varOmega}{g} \sin \left((t-t_0)\varOmega \right)\\ - \dfrac{g}{ \varOmega} \sin \left((t-t_0)\varOmega \right) & \cos \left((t-t_0)\varOmega\right) \end{array} \right], \end{align}

with $\varOmega = \omega (\boldsymbol {k}+\alpha \boldsymbol {k}_0, h)$, and

(5.2c) \begin{align} \left[ \begin{array}{@{}c@{}} \mathcal{I}^{(mj)}_A(\boldsymbol{k},t)\\ \mathcal{I}^{(mj)}_B(\boldsymbol{k},t) \end{array} \right] = \int_{t_0}^t \pmb{\mathcal{E}}((t-\tau)\varOmega) \exp({-\mathrm{i}(\,j-1)\beta \omega_0\tau}) \left[ \begin{array}{@{}c@{}} \widehat{\mathcal{N} }^{(mj)}_{A}(\boldsymbol{k}+(\,j-1)\alpha\boldsymbol{k}_0,\tau)\\ \widehat{\mathcal{N} }^{(mj)}_{B}(\boldsymbol{k}+(\,j-1)\alpha\boldsymbol{k}_0,\tau) \end{array} \right] \,\mathrm{d} \tau. \end{align}

Similar to a NLS equation-based model, the analytical solution of the envelope evolution equations in the form of (5.2a) consists of two terms: a linear and nonlinear term corresponding to the first term and the term of double sums on the right-hand side of (5.2a), respectively. Equation (5.2a) is obviously accurate for the evolution of linear surface waves. Owing to the fact that the integrand of the integrals in (5.2c) depends on the time-dependent envelopes, i.e. in an implicit form, the computation of (5.2a) requires a time integration method for the temporal-spatial evolution of nonlinear waves. To this end, there are many available approaches, e.g. the mid-point or the fourth-order Runge–Kutta method. We note that a leading-order scale of both $\widehat {\mathcal {N}}^{(mj)}_{A}$ and $\widehat {\mathcal {N}}^{(mj)}_{B}$ is $O(\epsilon ^2)$ and their time derivative $\partial _t\widehat {\mathcal {N}}^{(mj)}_{A}$ and $\partial _t\widehat {\mathcal {N}}^{(mj)}_{B}$ have the scale $O(\epsilon ^2\varepsilon _\text {f})$, where $\varepsilon _\text {f}$ denotes the dimensionless frequency bandwidth. The aforementioned scale of analysis can be clearly demonstrated in § 6.1.1. Therefore, for updating (5.2c) for one time step in the time interval $[t_n, t_n+\Delta t]$ with $t_n$ a time instant when $\hat {A}(\boldsymbol {k},t_n)$ and $\hat {B}(\boldsymbol {k},t_n)$ were computed and $\Delta t$ a small time interval, a numerical algorithm for time integration is required. The forward Euler method that is first-order accurate in the temporal interval leads to

(5.3) \begin{equation} \left[ \begin{array}{@{}c@{}} \mathcal{I}^{(mj)}_A(\boldsymbol{k},t)\\ \mathcal{I}^{(mj)}_B(\boldsymbol{k},t) \end{array} \right] = \boldsymbol{I}_n(t_{n+1}) \left[ \begin{array}{@{}c@{}} \widehat{\mathcal{N} }^{(mj)}_{A}({\boldsymbol{k}+(\,j-1)\alpha\boldsymbol{k}_0},t_n)\\ \widehat{\mathcal{N} }^{(mj)}_{B}({\boldsymbol{k}+(\,j-1)\alpha\boldsymbol{k}_0},t_n) \end{array} \right] + O(\epsilon^2{\varepsilon_\text{f}}\varepsilon_t \omega_0\Delta t), \end{equation}

where $t_{n+1}=t_{n}+\Delta t$, $\varepsilon _t\sim 1/(\omega _0\Delta t)\ll 1$ such that the approximation to the time integration given by (5.3) numerically converges, and the rapidly varying integral matrix $\boldsymbol {I}_n$ is defined as, and thereby given by,

(5.4a)\begin{align} &\boldsymbol{I}_n(t_{n+1}) \equiv \int_{t_n}^{t_{n+1}} \pmb{\mathcal{E}}((t_{n+1}-\tau)\varOmega) \exp({-\mathrm{i}(\,j-1)\beta \omega_0\tau}) \,\mathrm{d} \tau \end{align}

(5.4b)\begin{align} &\quad =-\dfrac{1}{2} \mathrm{i} \exp({-\mathrm{i}(\,j-1)\beta \omega_0t_{n+1}}) \left\{ \dfrac{1-\exp({\mathrm{i} [\omega_0+ (\,j-1)\beta\omega_0+\varOmega]\Delta t})}{\omega_0+ (\,j-1)\beta\omega_0+\varOmega } \left[ \begin{array}{@{}cc@{}} 1 & \dfrac{\varOmega}{g} \\ - \dfrac{g} { \varOmega} & 1 \end{array} \right]\right. \nonumber\\ & \qquad +\left. \dfrac{1-\exp({\mathrm{i} [\omega_0+ (\,j-1)\beta\omega_0-\varOmega]\Delta t})}{\omega_0+ (\,j-1)\beta\omega_0-\varOmega } \left[ \begin{array}{@{}cc@{}} 1 & -\mathrm{i} \dfrac{\varOmega}{g}\\ \mathrm{i} \dfrac{g}{ \varOmega} & 1 \end{array} \right] \right\}. \end{align}

It is worth noting that a numerical algorithm more accurate than the forward Euler in the sense of time integration can be used to evaluate (5.2c). An exponential integrator can also be used for the HOS method following similar procedures in this section.

Therefore, through an exponential integrator and the forward Euler method, we obtain

(5.5)$$\begin{gather} \left[ \begin{array}{c} {\hat{A}}(\boldsymbol{k},t_{n+1}) \\ {\hat{B}}_s(\boldsymbol{k},t_{n+1}) \end{array} \right] = \pmb{\mathcal{E}}((t_{n+1}-t_0)\varOmega) \left[ \begin{array}{c} {\hat{A}}(\boldsymbol{k},t_0) \\ {\hat{B}}_s(\boldsymbol{k},t_0) \end{array}\right] \nonumber\\ + \sum_{m=1}^{m=M} \sum_{j=0}^{j=m} \boldsymbol{I}_n(t_{n+1}) \left[ \begin{array}{@{}c@{}} \widehat{\mathcal{N} }^{(mj)}_{A}(\boldsymbol{k}+(\,j-1)\alpha\boldsymbol{k}_0,t_n)\\ \widehat{\mathcal{N} }^{(mj)}_{B}(\boldsymbol{k}+(\,j-1)\alpha\boldsymbol{k}_0,t_n) \end{array} \right] + O(\epsilon^2\varepsilon_t \omega_0\Delta t), \end{gather}$$

which needs to be updated step by step from the initial instant $t=t_0$ with given initial conditions. Compared with the HOS method implemented by Ducrozet et al. (2016) where the time interval depends on the shortest wave period, (5.5) permits a time interval that admits $\omega _0\Delta t\sim O(1/(\epsilon ^2\varepsilon _\text {f}))$ for numerical stability and convergence and thereby a much larger value without compromising the numerical efficiency with a careful choice of the value for $\beta$, as explained in § 6.2.2.

5.2. The energy balance equation

Based on the CEEEs in the Fourier plane expressed as (5.1), it is straightforward to derive the energy balance equation. To this end, the parameter

(5.6)\begin{equation} \hat{a} = \sqrt{\dfrac{g}{2\varOmega}}\hat{A} + \mathrm{i} \sqrt{\dfrac{\varOmega}{2g}} \hat{B} \end{equation}

is introduced, which is a function of the same (slow) time as envelope $\hat {A}$ and $\hat {B}_s$. The following sequential procedures are used: (i) multiplying the first and second (sub-) equation of (5.1) by $\sqrt {g/(2\varOmega )}$ and $\mathrm {i} \sqrt {\varOmega /(2g)}$, respectively, (ii) adding up the two resulting equations, (iii) replacing the terms that correspond to the definition of the new variable $\hat {a}$. Thereby, we obtain

(5.7)\begin{equation} \partial_t\hat{a} + \mathrm{i}(\varOmega-\beta\omega_0) \hat{a} = \sqrt{\dfrac{g}{2\varOmega}} \widehat{\mathcal{N}}_A + \sqrt{\dfrac{\varOmega}{2g}}\widehat{\mathcal{N}}_B. \end{equation}

Multiplying (5.7) and its complex conjugates by $\hat {a}^*$ and $\hat {a}$, respectively, and adding up the resulting equations leads to

(5.8a)\begin{equation} \partial_{t} ( \hat{a}\hat{a}^* )= \dfrac{g}{2\varOmega} (\widehat{\mathcal{N}}_{B,M}^*\hat{a} + \widehat{\mathcal{N}}_{B,M}\hat{a}^*) + \dfrac{\varOmega}{2g} (\widehat{\mathcal{N}}_{B,M}^*\hat{a} + \widehat{\mathcal{N}}_{B,M}\hat{a}^*), \end{equation}

in which $\hat {a}\hat {a}^*$ denotes the wave action similar to that defined in an extensive body of literature, e.g. Zakharov (1968), Stiassnie & Shemer (1984), Krasitskii (1994), Annenkov & Shrira (2009) and Gramstad (2014). Equation (5.8a) is known as the energy balance equation or the wave action equation. With the nonlinear effects neglected, the conservation of linear wave actions is evident,

(5.8b)\begin{equation} \partial_{t} ( \hat{a}\hat{a}^* )= 0, \end{equation}

which suggests that the transfer of wave actions does not occur between linear waves, as it should be.

5.3. Nonlinear forcing of waves

Due to the nonlinear effects that lead to the forcing terms on the right-hand side of the CEEEs described by (4.37a,b), it is understood that both free and bound (locked) waves, which do and do not obey the dispersion relation, respectively, can be forced. The former can arise from resonant and instability conditions that are the result of a combination of bandwidth and nonlinearity, as studied by numerous works, noticeably Phillips (1960), Hasselmann (1962), Benjamin & Feir (1967), Zakharov (1968), Longuet-Higgins (1978) and McLean (1982a,b). We show in this section how the nonlinear terms of different wave harmonics on the right-hand side of CEEEs lead to the forcing of bound waves and resonant free waves.

For later reference, we introduce the dimensionless wave vector $\boldsymbol {e}_n=\boldsymbol {k}_n/k_0$ and wave frequency $\sigma _n=\omega (|\boldsymbol {k}_n|h)/\sqrt {gk_0}$, which are given by, respectively,

(5.9a,b)\begin{equation} \boldsymbol{e}_{n} = (1+p_n, q_n)\quad \text{and}\quad \sigma_n = \sqrt{\tanh (k_0|\boldsymbol{e}_n|h)} \quad\text{for}~n\in\{1,2,3,\ldots\}, \end{equation}

where $n$ denotes the $n$th free wave and $\boldsymbol {e}_n\boldsymbol {\cdot }\boldsymbol {e}_0>0$ with $\boldsymbol {e}_0 =(1,0)$ and $p_n>-1$ and $q_n$ an arbitrarily chosen parameter. For an infinitesimal wave steepness, the resulting dimensionless wave vector, $\boldsymbol {e}_N$, and angular frequency, $\sigma _N$, due to the nonlinear forcing terms $\mathcal {N}_{A}^{(mj)}\exp ({\mathrm {i}(\,j-1)(\alpha \boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\beta \omega _0t)})$ and $\mathcal {N}_{B}^{(mj)}\exp ({\mathrm {i}(\,j-1)(\alpha \boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\beta \omega _0t)})$, can be obtained through the analysis in the Fourier plane and the superposition of linear waves. Especially, they are obtained by inserting the linear approximations of the unknown envelopes

(5.10a)$$\begin{gather} \hat{A}(\boldsymbol{k},t)\approx \hat{A}(\boldsymbol{k},t_0)\exp[\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}-\mathrm{i}(\varOmega-\beta\omega_0)(t-t_0)] \quad \text{and} \end{gather}$$

(5.10b)$$\begin{gather}\hat{B}_s(\boldsymbol{k},t)\approx \hat{B}_s(\boldsymbol{k},t_0)\exp[\mathrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}-\mathrm{i}(\varOmega-\beta\omega_0)(t-t_0)], \end{gather}$$

into the envelopes of the potential and vertical velocity and thereafter the nonlinear forcing term of the $j$th harmonic at $O(\epsilon ^m)$; explicitly, we arrive at

(5.11a)$$\begin{gather} \boldsymbol{e}_{N} -\alpha \boldsymbol{e}_0= (\,j-1)\alpha\boldsymbol{e}_0 + \sum_{n=1}^m\pm( \boldsymbol{e}_{n}- \alpha \boldsymbol{e}_0 ), \end{gather}$$

(5.11b)$$\begin{gather}\sigma_{N} -\beta \sigma_0= (\,j-1)\beta\sigma_0 + \sum_{n=1}^m \pm( \sigma_n-\beta \sigma_0 ), \end{gather}$$

where $\sigma _0=\sqrt {\tanh k_0h}$, $\boldsymbol {e}_N\boldsymbol {\cdot }\boldsymbol {e}_0>0$ holds for non-vanishing $\hat {A}$ and $\hat {B}_s$ by definition, and $N=m+1$ denotes the number of waves involved in the interaction. As to which sign to choose between ‘$\pm$’ depends on $j$; for $j=m$, the ‘+’ sign needs to be taken for all $n$ values whereas it is not permitted to choose the ‘$-$’ sign for all $n$ values as it will lead to the inequality $\boldsymbol {e}_N\boldsymbol {\cdot }\boldsymbol {e}_0<0$ where $\hat {A}$ and $\hat {B}_s$ vanish. We note that $N=3$, $N=4$ and $N=5$ correspond to triad, quartet and quintet wave interactions, which occur at the second, third and fourth order in wave steepness, respectively. With $\alpha =0$ and $\beta =0$, one can readily see that (5.11a,b) are simply the kernels arising from the $N$-wave interaction based on the Zakharov integral equation (Stiassnie & Shemer 1984; Shrira, Badulin & Kharif 1996; Janssen 2003), as it should be.

Physically, the resonant condition for $N$ waves interaction means that the resulting wave vector and frequency obey the dimensionless linear dispersion relation

(5.12)\begin{equation} \sigma_N - \sqrt{\tanh (k_0|\boldsymbol{e}_N|h)} = 0, \end{equation}

which can be used in the analysis of the nonlinear forcing of free waves arising from the interaction between linear waves. Mathematically, it corresponds to the particular terms (i.e. the nonlinear forcing terms) of the non-homogeneous CEEEs that have components that coincide with the eigenvalues of the homogeneous CEEEs in frequency, leading to a linear growth in time similar to the discussion by Hasselmann (1962). This point in principle determines whether the nonlinear forcing of waves are free or bound, which are discussed in §§ 5.3.1 and 5.3.2, respectively.

5.3.1. Forcing of free waves due to resonant effects

The nonlinear forcing of free waves due to the class I resonant condition occurs at third order in wave steepness for waves of the first harmonic and corresponds to the effects of the nonlinear forcing terms with $m=3$ and $j=1$ in the CEEEs. Following the previous works, e.g. Stiassnie & Shemer (1984), the class I resonant condition due to quartet (linear) wave interaction in the CEEEs is given by

(5.13a,b)\begin{equation} \boldsymbol{e}_4 = \boldsymbol{e}_1+\boldsymbol{e}_2-\boldsymbol{e}_3\quad \text{and}\quad \sigma_4 = \sigma_1+\sigma_2-\sigma_3, \end{equation}

where $\boldsymbol {e}_1=(1+p, q) , \boldsymbol {e}_2= (|1-p|,-q), \boldsymbol {e}_3=\boldsymbol {e}_4\equiv (1,1)$ and $\sigma _4=\sigma _0$. Similarly, the class II resonant condition due to quintet wave interaction occurs due to the nonlinear forcing terms with $m=4$ and $j=2$ in the CEEEs where

(5.14a,b)\begin{equation} \boldsymbol{e}_5 = \boldsymbol{e}_1+\boldsymbol{e}_2+\boldsymbol{e}_3-\boldsymbol{e}_4\quad \text{and}\quad \sigma_5 = \sigma_1+\sigma_2+\sigma_3-\sigma_4, \end{equation}

with $\boldsymbol {e}_5=(1,0)$, $\sigma _5=\sigma _0$. The effect of wave nonlinearity on the dispersion relation is neglected in (5.14a,b). In order to account for this, the ‘Stokes-corrected’ nonlinear dispersion relation, which can be obtained from (2.21c) by Stiassnie & Shemer (1984), should be used for $\sigma _n$ with $n\in \{1,2,3,4,5\}$.

5.3.2. Bound waves

With $j=m$ for $m\geq 2$, the inequality $\boldsymbol {e}_n\boldsymbol {\cdot }\boldsymbol {e}_0>0$ holds for all wave vectors in (5.11a,b) due to the definition of the envelope transform and the expression of $B^{(nn)}$ and $A$. This can be demonstrated by a simple example. To this end, we choose the third term in (4.33g). It is understood that

(5.15)$$\begin{gather} \tfrac{1}{8}A^2\partial_{zzz}B \exp({2\mathrm{i}(\alpha\boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}) \equiv \int \hat{A}(\boldsymbol{k}_1,t)\hat{A}(\boldsymbol{k}_2,t)|\boldsymbol{k}_3+\boldsymbol{k}_0|^3 \tanh(|\boldsymbol{k}_3+\boldsymbol{k}_0|h) \nonumber\\ \times\hat{B}_s(\boldsymbol{k}_3,t) \exp({\mathrm{i}(\boldsymbol{k}_1+\boldsymbol{k}_2+\boldsymbol{k}_3+2\alpha\boldsymbol{k}_0-2\beta\omega_0t)\boldsymbol{\cdot}\boldsymbol{x}})\,\mathrm{d} \boldsymbol{k}_1\,\mathrm{d}\boldsymbol{k}_2\,\mathrm{d}\boldsymbol{k}_3, \end{gather}$$

where $\boldsymbol {k}_n+\boldsymbol {k}_0 =k_0\boldsymbol {e}_n$ and $(\boldsymbol {k}_n+\boldsymbol {k}_0)\boldsymbol {\cdot }\boldsymbol {k}_0>0$ for non-vanishing $\hat {A}(\boldsymbol {k}_n,t)$ and $\hat {B}_s(\boldsymbol {k}_n,t)$ by definition for $n=1,2$ and $n=3$. The frequency superposition can be obtained through the linear approximation to the CEEEs,

(5.16)\begin{align} [\hat{A}(\boldsymbol{k}_n,t) , \hat{B}_s(\boldsymbol{k}_n,t) ]=[\hat{A}(\boldsymbol{k}_n,t_0) , \hat{B}_s(\boldsymbol{k}_n,t_0) ]\exp({-\mathrm{i}(\omega_n-\beta\omega_0 )(t-t_0)}) + O(\epsilon^2), \end{align}

where $\omega _n= \omega (\boldsymbol {k}_n+\alpha \boldsymbol {k}_0,h)$. Inserting (5.16) into (5.15) and taking the frequency superposition readily leads to the frequency combination on the right-hand side of (5.11b). The last step is to multiply the factor $\exp [\mathrm {i}(\alpha \boldsymbol {k}_0\boldsymbol {\cdot }\boldsymbol {x}-\beta \omega _0t)]$ on both sides of the resulting equation and, hence, we arrive at

(5.17)$$\begin{gather} \tfrac{1}{8}A^2\partial_{zzz}B \exp\left({3\mathrm{i} (\alpha \boldsymbol{k}_0\boldsymbol{\cdot}\boldsymbol{x}-\beta\omega_0t)}\right) \equiv \int \hat{A}(\pmb{\kappa}_1-\alpha\boldsymbol{k}_0,t_0)\hat{A}(\pmb{\kappa}_2-\alpha\boldsymbol{k}_0,t_0)\hat{B}_s(\pmb{\kappa}_3-\alpha\boldsymbol{k}_0,t_0) \nonumber\\ \times|\pmb{\kappa}_3|^3 \tanh(|\pmb{\kappa}_3|h) \exp\left({ \mathrm{i} \sum_{n=1}^3\left[\pmb{\kappa}_n\boldsymbol{\cdot}\boldsymbol{x}-\sqrt{g|\pmb{\kappa}_n|\tanh(|\pmb{\kappa}_n|h)}(t-t_0) \right] }\right)\,\mathrm{d} \pmb{\kappa}_1\,\mathrm{d}\pmb{\kappa}_2\,\mathrm{d}\pmb{\kappa}_3, \end{gather}$$

where the relation $\pmb {\kappa }_n=\boldsymbol {k}_n+\alpha \boldsymbol {k}_0$ was used for the change of the integral variables. It becomes evident that the non-vanishing integrand in (5.17) leads to the resulting waves whose dimensionless wave vector and frequency are given by

(5.18)\begin{equation} \boldsymbol{e}_4 = \sum_{n=1}^{3}\pmb{\kappa}_n/k_0\quad \text{and}\quad \sigma_4 = \sum_{n=1}^{3} \sqrt{g|\pmb{\kappa}_n|\tanh(|\pmb{\kappa}_n|h)} /\omega_0, \quad \text{with}\ \pmb{\kappa}_n\boldsymbol{\cdot}\boldsymbol{k}_0>0, \end{equation}

which can be written in a form similar to (5.13a,b), except that the third wave on the right-hand side needs to take the positive sign instead. Due to this, the resonant quartet wave condition cannot be satisfied. Thereby, the nonlinear term chosen as an example can only lead to the forcing of bound waves that do not obey the linear dispersion relation. Similar analysis can be carried out for the other remaining components of the nonlinear forcing terms with $m=j$. The third and fourth wave given by the quintet and quartet resonant condition in (5.13) and (5.14), respectively, needs to take a negative sign. In contrast, all waves in the nonlinear terms on the right-hand of the CEEEs with $m=j$ can only take the positive sign as the inequality $\boldsymbol {e}_n\boldsymbol {\cdot }\boldsymbol {e}_0>0$ should hold as noted. Therefore, it can be readily inferred that only bound waves can be forced by the nonlinear terms in the CEEEs with $m=j$.

The discussions above did not cover the forcing of either free or bound waves arising from the nonlinear forcing terms in the CEEEs for $j=0$ and $m\geq 2$ that only appear in the even orders in wave steepness and that are often responsible for the forcing of mean flows (or waves with a low or vanishing frequency). These nonlinear forcing terms $j=0$ cannot force free waves since their resulting wave vectors and wave frequencies do not obey the resonant conditions presented in § 5.3.1 for the interaction of infinitesimal waves. They are also associated with the singular terms in the Zakharov's kernel functions, as clearly stated in the introduction of Gramstad (2014). Evidently, nonlinearity effects on the dispersion relation have been neglected in this section, and thereby how they affect the forcing of mean flows has not been explored. At the second order with ‘$mj=20$’, we understand that subharmonic bound waves can be forced, as is well known; see, e.g. Phillips (1960) and Hasselmann (1962). The author conjectures that novel physics may be elucidated through the exploration of these terms (which can force nonlinear mean flows) with the consideration of higher-order nonlinearity. As it is not the main focus of this paper, this aspect will be left for future explorations.

6. Comparisons with two other methods

For further demonstrating the potential, the CEEEs are firstly compared with a traditional perturbation expansion for the evolution of a train of Stokes waves (§ 6.1.1) and irregular waves with an arbitrary bandwidth and directional spreading (§ 6.1.2). Next, we proceed to comparisons with the HOS method in § 6.2 of the nonlinear forcing terms in a limiting case (§ 6.2.1) and of the computational complexity (§ 6.2.2) to especially demonstrate the numerical efficiency of the CEEEs.

6.1. Relation with a traditional perturbation expansion

Using an example of both a train of Stokes waves in § 6.1.1 and the more general weakly nonlinear three-dimensional waves in § 6.1.2, we show in this section how to establish the relation between a traditional perturbation method and the CEEEs. The analytical analysis in the section has a twofold sub-goal. It firstly demonstrates that the CEEEs are correctly derived. Secondly, it shows the first few steps that are essential to more general derivations for bridging the relations between the CEEEs and other higher-order frameworks, e.g. the different versions of third-order NLS equations. For example, if both the third orders in $\epsilon _0$ and a narrow bandwidth are additionally considered in § 6.1.1, the classic third-order accurate NLS equation would be recovered based on the CEEEs. Or if the Stokes waves are considered up to the fifth order in $\epsilon _0$ in § 6.1.1, one would deduce the framework by Fenton (1985) starting from the CEEEs.

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 1985). 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 1968; Fenton 1985), 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 2008, § 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 (2021) 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 (2021), 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 (2021), 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 (2021) through the connections given by (6.17a–e) leads to the envelope equations of Li & Li (2021) 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 (2021). In particular, Li & Li (2021) 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 1987; West et al. 1987) and the Fourier kernels (Hasselmann 1962; Zakharov 1968; Dalzell 1999). 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 (2021), respectively.

6.2. Comparisons with the HOS method

The CEEEs are next compared with the HOS method in the perspective of the numerical implementation of the nonlinear forcing terms (§ 6.2.1) in a limiting case and of the computational efficiency for the same level of accuracy (§ 6.2.2). It is especially noted that the relation between the HOS method and the Zakharov integral equation has been shown in Onorato, Osborne & Serio (2007) and thereby the relation between the CEEEs and the Zakharov integral equation can be established through the HOS method.

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 (2000) and Slunyaev, Sergeeva & Didenkulova (2016) 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 (2008, § 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 (1985).

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 2021). 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.

7. Conclusions

This paper deals with the description of surface gravity waves on a finite water depth in the framework of potential flow theory. The main objective of this paper is to propose a framework that combines the merits of both the HOS method (Dommermuth & Yue 1987; West et al. 1987) and a NLS equation-based model (Zakharov 1968; Davey & Stewartson 1974; Dysthe 1979). In particular, it can be as accurate as the HOS method at no additional cost of numerical efficiency on the one hand. Similar to a NLS equation-based model, it shall be capable of describing the slow temporal evolution of a wave spectrum but at no expense of accuracy, permitting a coarse and large computational domain compared with the characteristic length of the wave phase, on the other hand. To this end, a novel theoretical framework has been presented in the Hamiltonian theory based on a perturbation expansion, which leads to the CEEEs given by (4.37a,b). The CEEEs can be derived up to arbitrary order in wave steepness and have a much slower temporal and larger spatial scale compared with the characteristic time and wavelength of the wave phase, respectively. In the new theoretical framework, the envelope of both surface elevation and the potential on the free water surface are introduced, which are shown to be a pair of canonical variables for the first time. The two envelopes are used for expressing wave fields and thereby are the main unknowns solved for numerically from the CEEEs.

A few main features of the CEEEs are explored in this paper. Firstly, based on the CEEEs, the energy balance equation for the evolution of wave actions is derived. Secondly, due to the fact that the CEEEs are composed of both linear and nonlinear terms in wave steepness, it is proposed to solve the CEEEs by using an exponential integrator that leads to the analytical description of linear wave fields. Furthermore, the nonlinear terms are expressed in a form of the separation of different wave harmonics, due to which they can be especially split into two categories: one which can only force bound waves that do not obey the dispersion relation and the other that is capable of the forcing of free waves if particular conditions are met. Analytical derivations are presented showing how the forcing of free waves can arise from quartet and quintet resonant interactions of linear waves. Much more physical implications remain to be explored in future works.

The newly derived framework has been compared with three different methods. Analytical relations between the CEEEs and two traditional perturbation methods are established, including the theory for the evolution of a train of Stokes waves up to second order by Fenton (1985) and the second-order semi-analytical framework for three-dimensional surface waves with arbitrary bandwidth and large directional spreading by Li & Li (2021). One would find that the relations established can be extended for more general cases. For example, proceeding to an order higher the CEEEs would be shown to recover a third-order accurate NLS equation in the limiting case of narrow-band waves (Trulsen et al. 2000) or of the removal of secular terms at the third order in wave steepness (Li 2021). From examining numerically the case of the evolution of a nonlinear focus wave group, the nonlinear terms from the second to fourth order based on the CEEEs are shown to be identical to those based on the HOS method. Compared with the HOS method for the same level of accuracy, the CEEEs do not require a larger number of FFTs for computations but allow for a much larger time interval and computational domain with both a larger size and spacing.

This new framework has potential for applicability in more general cases that account for the interaction between surface waves and their ambient environments with different scales, e.g. sub-mesoscale currents and small-scale turbulence. Despite the fact that surface tension and the slowly varying water depth are neglected in the framework, an extension to additionally consider these features would be straightforward. The new framework is in principle a spectral method and thereby it would suffer from the drawbacks due to the use of a FFT in a way similar to the HOS method.