2.1. Asymptotic analysis and the water-wave problem

The water-wave problem, whether written as a system of partial differential equations and boundary conditions or using the Hamiltonian reformulation, remains mathematically very complex. While the mathematical existence of periodic and solitary wave solutions has been established in ever more general settings (for example, the books by Constantin Reference Constantin2012 or Lannes Reference Lannes2013 offer a relatively recent overview of mathematical progress), much of what is known about qualitative and quantitative properties of ocean waves rests on an analysis of simplified equations.

These simplified equations are usually obtained through so-called asymptotic analysis or perturbation theory – a procedure of non-dimensionalisation and scaling the equations, followed by the identification of small dimensionless parameters, and an asymptotic expansion of the variables in terms of the small parameters selected. In the case of water waves, the most common small parameter is the wave steepness $\epsilon,$ written either as the wave amplitude divided by the wavelength or multiplied by the wavenumber. A very thorough discussion of asymptotic analysis as applied to water waves can be found in the textbook by Johnson (Reference Johnson1997).

Instead of expanding the water-wave equations (2.1a )–(2.1d ) in this manner (see Madsen & Fuhrman Reference Madsen and Fuhrman2012 for a detailed, relevant account), it is also possible to expand the Hamiltonian directly. This procedure was carried out by Zakharov (Reference Zakharov1968) and, with more detail, Krasitskii (Reference Krasitskii1994), and it has the advantage of enforcing structure that might otherwise be lost, such as that resulting equations conserve energy.

The first step in this procedure is to write the canonical equations in terms of the Fourier transformed variables $\hat {\zeta }(k)$ and $\hat {\phi }(k),$ as

(2.4) \begin{align} \frac {\partial \hat {\zeta }}{\partial t} = \frac {\delta H}{\delta \hat {{\phi ^{s}}}^*}, \quad \frac {\partial \hat {{\phi ^{s}}}}{\partial t} = - \frac {\delta H}{\delta \hat {\zeta }^*}, \end{align}

with $*$ denoting the complex conjugate. Then, one defines a new pair of canonical variables $a(k)$ and $ia^*(k)$ related to $\phi ^{s}$ and $\zeta$ via Fourier transforms as

(2.5) \begin{align} \zeta (x) = \frac {1}{2\pi } \int \left ( \frac {q(k)}{2 \omega (k)}\right )^{1/2} \left [ a(k) + a^*(-k)\right ] e^{ikx} \,{\rm d}k, \\[-12pt] \nonumber \end{align}

(2.6) \begin{align} {\phi ^{s}}(x) = \frac {-i}{2\pi } \int \left ( \frac {\omega (k)}{2 q(k)}\right )^{1/2} \left [ a(k) - a^*(-k)\right ] e^{ikx} \,{\rm d}k,\\[9pt] \nonumber \end{align}

where

(2.7) \begin{align} q(k) = |k| \tanh (|k|h), \end{align}

$\omega = \sqrt {gq(k)}$ is the frequency (in rad s⁻¹) and $g$ is the constant acceleration of gravity (taken to be 9.8 m s⁻ $^2$ in computations). As we will be interested in deep-water waves ( $h\rightarrow \infty )$ , we note that the dispersion relation simplifies to

(2.8) \begin{align} \omega ({k})^2 = g | {k}|. \end{align}

We will use this in all the computations that follow.

The subsequent steps in the expansion are somewhat lengthy and are given by Krasitskii (Reference Krasitskii1994). At this stage, it is important to emphasise only one key point: at lowest order, the problem is linear, so the superposition principle holds. This means that the Fourier amplitudes – which correspond essentially to the magnitudes of the terms $a(k)$ introduced in (2.5)–(2.6) – do not change with time. If gravity is the sole restoring force, the surprising fact is that the Fourier amplitudes do not change with time at the next order either. Only at third-order in nonlinearity, where resonance between quartets of gravity waves becomes possible, do we obtain an evolution equation. This equation gives the slow evolution of the complex amplitude $b_0=b(k_0,t),$

(2.9) \begin{align} i \frac {\partial b_0}{\partial t} = \omega _0 b_0 + \int \tilde {V}^{(2)}_{0123} b_1^* b_2 b_3 \delta _{0+1-2-3}\, {\rm d}k_{123}, \end{align}

and is known as the Zakharov equation, with a kernel term $\tilde {V}^{(2)}$ given by Krasitskii (Reference Krasitskii1994, (3.5)). We have introduced subscript notation to denote wavenumbers, so that, for example, the delta-function $\delta _{0+1-2-3}$ is a shorthand for $\delta (k_0 + k_1 - k_2 - k_3).$ In principle, $k_i$ can be either a scalar or a vector wavenumber, although our computations will deal only with the scalar case for simplicity.

If $b$ is known, then the canonical variable $a$ can be recovered via

(2.10) \begin{align} a_0 = b_0 &+ \int A^{(1)}_{012} b_1 b_2 \delta _{0-1-2}\, {\rm d}k_{12} + \int A^{(2)}_{012} b_1^* b_2 \delta _{0+1-2} \,{\rm d}k_{12} + \int A^{(3)}_{012} b_1^* b_2^* \delta _{0+1+2} \,{\rm d}k_{12} \nonumber \\ &+ \int B^{(1)}_{0123} b_1 b_2 b_3 \delta _{0-1-2-3}\, {\rm d}k_{123} + \int B^{(2)}_{0123} b_1^* b_2 b_3 \delta _{0+1-2-3} \,{\rm d}k_{123} \nonumber \\ &+ \int B^{(3)}_{0123} b_1^* b_2^* b_3 \delta _{0+1+2-3}\, {\rm d}k_{123} + \int B^{(4)}_{0123} b_1^* b_2^* b_3^* \delta _{0+1+2+3}\, {\rm d}k_{123} + \cdots, \end{align}

where we have written terms up to third order only. This expression is given by Krasitskii (Reference Krasitskii1994, (2.17)), as are the relevant interaction kernels $A^{(i)}, \, B^{(i)}$ . Finally, once $a$ is known, it is possible to recover $\zeta$ and $\phi ^{s}$ up to the desired order by means of (2.5)–(2.6).

2.2. Constant magnitude approximation

The procedure outlined previously would seem to suggest that we first need to solve the complicated integro-differential Zakharov equation (2.9). In fact, for certain situations of interest, we can circumvent this step. The key insight is that the evolution described by the Zakharov equation – being associated with cubic nonlinearities – occurs on the slow time scale $\epsilon ^{-2} T$ , where $\epsilon$ is a typical wave steepness and $T$ is a typical wave period. For processes taking place on faster time scales, it may be appropriate to consider the complex amplitudes to be constant.

This is indeed the approach we will adopt. There is, however, one important addition: for certain wavenumber combinations, the Zakharov equation (2.9) can be solved explicitly. In particular, for all symmetric quartets which satisfy $k_0 + k_0 = k_0 + k_0$ or $k_0 + k_1 = k_0 + k_1$ , the equation reduces to a system of one or two ordinary differential equations with constant-amplitude solutions, as shown by Mei, Stiassnie & Yue (Reference Mei, Stiassnie and Yue2018, Ch. 14). Such solutions represent nonlinear corrections to the frequency of the waves and were first found by Stokes for a single wave mode. For arbitrarily many modes, or a continuous spectrum, it is possible to write these corrections in a compact form (Stuhlmeier & Stiassnie Reference Stuhlmeier and Stiassnie2019).

Incorporating these corrections means that the complex amplitudes $b_i(t)$ have constant magnitudes but time-dependent phases,

(2.11) \begin{align} b_i(t) = |b_i(0)| \exp \bigg ({-}it \bigg [ \tilde {V}^{(2)}_{\textit{iiii}} |b_i(0)|^2 + 2 \sum _{j\neq i} \tilde {V}^{(2)}_{\textit{ijij}} |b_j(0)|^2 \bigg ] \bigg ). \end{align}

The terms in the argument of the exponential reduce to $\omega a^2 k^2/2$ , the well-known Stokes’ correction, when only a single wave-mode is present, or the mutual frequency correction found by Longuet-Higgins & Phillips (Reference Longuet-Higgins and Phillips1962) for two modes. The form (2.11) now contains all the mutual dispersion corrections of the cubically nonlinear water-wave theory. It has been demonstrated that taking these corrections into account allows one to predict the short-time propagation of laboratory waves with high accuracy without the need to solve for the evolution of the slow-time variables (Galvagno, Eeltink & Stuhlmeier Reference Galvagno, Eeltink and Stuhlmeier2021; Stuhlmeier & Stiassnie Reference Stuhlmeier and Stiassnie2021; Meisner et al. Reference Meisner, Galvagno, Andrade, Liberzon and Stuhlmeier2023), showing that the constant magnitude approximation is well justified for such cases. Of course, under specific circumstances – for example, with a sufficiently narrow spectrum or close to modulation instability – we may have significant energy transfer (Andrade & Stuhlmeier Reference Andrade and Stuhlmeier2023). However, given our focus on particle paths (with a characteristic time of the order of one period), the amplitude evolution will be considered to be negligible (see also the discussion by Stuhlmeier & Stiassnie (Reference Stuhlmeier and Stiassnie2019) and Appendix A therein).

If we sum up the consequences of weak-nonlinearity in the Fourier description of surface gravity waves in deep water as:

1. (1) dispersion corrections;

2. (2) changes to wave shape (bound harmonics);

3. (3) energy exchange between modes,

we will account for items (1) and (2) while neglecting the slow effects of item (3). Only the amplitude-dependence of the dispersion relation is a general property of nonlinear waves. The other effects arise from the viewpoint of Fourier analysis. Solitary waves, for example, do not have a counterpart in linear theory whose ‘shape’ is modified by the inclusion of higher-order effects. Appendix C provides some further discussion of the neglect of amplitude evolution, along with several examples.

2.3. Recovery of the third-order solution

The Hamiltonian formulation makes clear how the free surface and potential can be recovered from the canonical variables $a(k)$ and $a^*(k)$ up to a given order. In practice, discrete wavenumbers are considered, which means that the integrals in (2.10) can be simplified by the substitution

(2.12) \begin{align} b(k) = \sum _i B_i \delta (k-k_i). \end{align}

Here $B_i=B(k_i,t)$ is unknown in principle and, in practice, comes from a constant magnitude approximation to (2.9) as discussed in § 2.2.

Because the free surface equation (2.5) is recovered directly from substitution of (2.10), it is clear that we can distinguish linear, quadratic and cubic terms in $\zeta$ based on the respective terms in (2.10). The same is not true of $\phi$ , as we shall see later. Taking the inverse Fourier transform yields

(2.13) \begin{align} \zeta _1 = \frac {1}{2 \pi } \sum _i \sqrt {\frac {\omega _i}{2g}} A_i \big ( e^{i \xi _i} + e^{-i \xi _i} \big ). \end{align}

Here, we have substituted $B_i = A_i e^{-i\varphi _i(t)}$ with amplitude $A_i \in \mathbb{R}$ and phase $\varphi _i(t),$ and employ the compact notation $\xi _i = k_i x - \varOmega _i t + \theta _i.$ Clearly, this yields the expected sinusoidal free surface elevation upon defining the amplitude

(2.14) \begin{align} \frac {1}{\pi } \sqrt {\frac {\omega _i}{2g}} A_i = a_i. \end{align}

The nonlinear corrected frequency is given by

(2.15) \begin{align} \varOmega _i = \omega _i + \tilde {V}^{(2)}_{\textit{iiii}} |A_i|^2 + 2 \sum _{j\neq i} \tilde {V}^{(2)}_{\textit{ijij}} |A_j|^2, \end{align}

and for unidirectional waves in infinite depth – which will be our focus – we can simplify the kernels to

(2.16) \begin{align} \tilde {V}^{(2)}_{\textit{iiii}} = \frac {1}{4 \pi ^2} k_i^3 \text{ and } \tilde {V}^{(2)}_{\textit{ijij}} = \frac {1}{4 \pi ^2} k_i k_j \min (k_i,k_j). \end{align}

The quadratic components of the free surface elevation can be obtained in analogous manner, with slightly more in the way of algebra required to resolve all delta-functions. We find

(2.17) \begin{align} \zeta _2(x,t) = & \frac {1}{2\pi } \sum _{i,j} A_i A_j \left \lbrace \sqrt {\frac {\omega _{i+j}}{2g}} \left ( A^{(1)}_{i+j,i,j} + A^{(3)}_{-i-j,i,j} \right ) \big [ e^{i(\xi _i + \xi _j)} + e^{-i(\xi _i + \xi _j)} \big ] \right . \nonumber \\ & \left . + \sqrt {\frac {\omega _{i-j}}{2g}} A^{(2)}_{j-i,i,j} \big [ e^{i(\xi _j-\xi _i)} + e^{i(\xi _i-\xi _j)} \big ] \right \rbrace . \end{align}

The shorthand notation $\omega _{i\pm j}$ is used to denote $\omega (k_i \pm k_j) = \sqrt {g |k_i \pm k_j|}.$

We can write the cubic components of the free surface elevation as

(2.18) \begin{align} \zeta _3(x,t) = & \frac {1}{2 \pi } \sum _{i,j,k} A_i A_j A_k \left \lbrace \sqrt {\frac {\omega _{i+j+k}}{2g}} \left ( B^{(1)}_{i+j+k,i,j,k} + B^{(4)}_{-i-j-k,i,j,k} \right ) \big [ e^{i(\xi _i + \xi _j + \xi _k)} \big . \right . \nonumber \\& \left. \big . + e^{-i(\xi _i + \xi _j + \xi _k)} \big ] + \sqrt {\frac {\omega _{i-j-k}}{2g}} B^{(2)}_{j+k-i,i,j,k} \big [ e^{i(\xi _j+\xi _k-\xi _i)} + e^{i(\xi _i -\xi _j - \xi _k)} \big ] \right . \nonumber \\& \left . + \sqrt {\frac {\omega _{i+j-k}}{2g}} B^{(3)}_{k-i-j,i,j,k} \big [ e^{i(\xi _k-\xi _i-\xi _j)} + e^{i(\xi _i+\xi _j-\xi _k)} \big ] \right \rbrace . \end{align}

The velocity potential at the free surface, $\phi ^{s}$ , is obtained from (2.6) via an identical procedure.

The contributions of the second-order and third-order free surface elevation change only the wave shape. The resulting harmonics of the form $\exp (i(\xi _i \pm \xi _j))$ behave as waves with wavenumber $k_i \pm k_j$ and frequency $\varOmega _i \pm \varOmega _j,$ but because these sums and differences do not satisfy the dispersion relation, they are referred to as ‘bound waves’, in contrast to the free waves of (2.13). The impact of these bound waves on sharpening the crests and flattening the troughs is shown in figure 1 for a bichromatic sea, together with the effect of the third-order frequency correction (2.15) on the wave phases. We also note the important effect of the difference harmonic terms $k_i-k_j$ in creating a ‘set-down’ below the crest, see van den Bremer et al. (Reference van den Bremer, Whittaker, Calvert, Raby and Taylor2019).

2.3.1. Recovery of the bulk potential

The reformulation of the problem in surface variables is very convenient from many perspectives, but is something of a disadvantage when we wish to recover information about the bulk fluid below the free surface. To obtain this, we start with the general solution of the Laplace equation which decays for infinite depth

(2.19) \begin{align} \phi (x,z) = \frac {1}{2\pi } \int \phi (k) e^{kz} e^{ikx}\, {\rm d}k, \quad \phi (k) = \phi ^*(-k), \end{align}

the exact analogue of Krasitskii (Reference Krasitskii1994, (4.1)).

Figure 1.

Free surface of a bichromatic wave train with linear frequencies $\omega _1=1$ rad s⁻¹ and $\omega _2=1.25$ rad s⁻¹, and wave slopes $\epsilon _1=\epsilon _2=0.15,$ showing first-order (linear), second-order and third-order solutions and their bound wave (B.W.) constituents. Note the phase shift appearing at third order due to the frequency correction (4.4)–(4.5).

We relate the bulk potential $\phi$ and its surface trace $\phi ^{s}$ to one another in terms of an expansion. Taking the problem in infinite depth, we need to solve

(2.20) \begin{align} &\Delta \phi = 0, \\[-12pt] \nonumber \end{align}

(2.21) \begin{align} &\phi = {\phi ^{s}}\quad \text{ on } z = \zeta , \\[-12pt] \nonumber \end{align}

(2.22) \begin{align} &\phi _z \rightarrow 0 \quad\text{ as } z \rightarrow -\infty , \\[9pt] \nonumber \end{align}

which we do by assuming the free surface elevation is small ( $\zeta \rightarrow \epsilon \zeta$ ) and then inserting this expansion into the boundary condition (2.21), so that

(2.23) \begin{align} \phi (\zeta ) = \phi (0) + \epsilon \zeta \phi _z(0) + \frac {\epsilon ^2 \zeta ^2}{2} \phi _zz(0) + \cdots . \end{align}

This gives rise to a hierarchy of problems in physical space, whose successive solution yields

(2.24) \begin{align} \hat {\phi }(k,z,t) = e^{|k|z} \left [ \hat {\phi }^s - \frac {1}{2 \pi } \int \hat {\zeta }_1 \hat {\phi }^s_2 |k_2| \delta _{0-1-2}\, {\rm d}k_{12} - \int D_{0123} \hat {\phi }^s_1 \hat {\zeta }_2 \hat {\zeta }_3 \delta _{0-1-2-3} \,{\rm d}k_{123} + \cdots \right ], \end{align}

with

(2.25) \begin{align} D_{0123} = -\frac {1}{4} |k_1| \left ( 2 |k_1| - |k_1 + k_2| - |k_1 + k_3| - |k_0-k_3| - |k_0-k_2| \right ). \end{align}

With $\hat {{\phi ^{s}}}$ and $\hat {\zeta }$ known from the formulation in (2.5)–(2.6), we can now recover the bulk potential at each order desired. This procedure is also detailed by Janssen (Reference Janssen2004, below (4.7)), Krasitskii (Reference Krasitskii1994, p.15) and Zakharov (Reference Zakharov1968, (1.8)), although care should be taken in the symmetrisation of (2.25).

The linear contribution consists simply of the first term in (2.24),

(2.26) \begin{align} \phi _1(x,z,t) = \frac {1}{\pi } \sum _i A_i e^{|k_i|z} \sqrt {\frac {g}{2\omega _i}} \sin (\xi _i). \end{align}

The quadratic contribution is

(2.27) \begin{align} \phi _2(x,z,t) & = \sum _{i,j} A_i A_j \left [ \mathcal{C}^{(2)}_{i+j} \sin (\xi _i + \xi _j) e^{|k_i+k_j|z} + \mathcal{C}^{(2)}_{i-j} \sin (\xi _i - \xi _j) e^{|k_i-k_j|z} \right ], \end{align}

where we need to employ the linear and quadratic parts of $\hat {\phi }^s$ and $\hat {\zeta }$ (and products thereof) in the formulation. The coefficients $\mathcal{C}^{(2)}$ are given in Appendix B. Up to second order, for unidirectional waves in deep water, the potential reduces to the simple expression

(2.28) \begin{align} \phi = \sum _j \frac {a_j g}{\omega _j} e^{|k_j|z} \sin ( k_j x - \omega _j t) - \sum _{i\gt j} \omega _{i} a_i a_j \sin (\xi _{i} - \xi _{j} ) e^{|k_i - k_j| z}, \end{align}

which follows also from Dalzell (Reference Dalzell1999).

The cubic contribution to the bulk potential is

(2.29) \begin{align} \phi _3(x,z,t) &= \sum _{ijk} A_i A_j A_k \left \lbrace \mathcal{C}^{(3)}_{i+j+k} \sin (\xi _i+\xi _j+\xi _k) e^{|k_i+k_j+k_k|z} \right . \nonumber \\ &+ \mathcal{C}^{(3)}_{i-j-k} \sin (\xi _i-\xi _j-\xi _k) e^{|k_i-k_j-k_k|z} \left . + \mathcal{C}^{(3)}_{i+j-k} \sin (\xi _i + \xi _j - \xi _k) e^{|k_i+k_j-k_k|z} \right . \nonumber \\ &+ \left . \mathcal{C}^{(3)}_{i-j+k} \sin (\xi _i-\xi _j+\xi _k) e^{|k_i-k_j+k_k|z} \right \rbrace . \end{align}

The coefficients $\mathcal{C}^{(3)}$ are given in Appendix B. Unfortunately, efforts to find a compact simplification for the third order have failed except in the case of a single monochromatic wave, where the analytical expression

(2.30) \begin{align} \phi _3 = -\frac {a^3 \omega k}{4} e^{kz} \sin \xi \end{align}

is recovered (see Gao et al. Reference Gao, Sun and Liang2021 and Appendix A).

It should be emphasised that the bulk potential thus obtained is equivalent (except for the non-uniqueness at a given order mentioned in Appendix A) to the potential found through direct perturbation expansion of the governing equations (2.1a )–(2.1d^′ ). In particular, as we shall see later, this means that the fluid domain at each order is the half-space $\{ (x,z) \mid x\in \mathbb{R}, z\leq 0 \},$ owing to the transfer of the boundary conditions.

3.1. Lagrangian drift with depth

The depth-dependent Lagrangian drift can be obtained from integrating the particle trajectory mapping, provided particles do not cross the still-water level $z=0$ – otherwise, as mentioned previously, extrapolating the terms $\exp (kz)$ can lead to spuriously high velocities or even divergence of the trajectories. Because the forward drift of a particle is phase-dependent, this means averaging over the phases of particles lying at an initial level $z_0.$ In addition to the linear solution (3.1), we will employ Eulerian solutions

(3.10) \begin{align} \phi = \frac {a \omega }{k}e^{kz} \sin \xi + \frac {a^4 k^2 \omega }{2} e^{2kz} \sin 2 \xi + \frac {a^5 k^3 \omega }{12}e^{3 kz} \sin 3 \xi + \cdots \end{align}

up to fifth order (second and third order may be obtained from § 2.3, see also Appendix A, fourth and fifth order are given for Stokes waves by Zhao & Liu (Reference Zhao and Liu2022)), which can be readily inserted into the right-hand side of the particle trajectory mapping (3.4). From (3.10), it can be noted that $u_E$ is zero at all orders, so that the Lagrangian drift $u_L$ is equal to the Stokes drift $u_S,$ whose leading order constituent is given in (3.9).

In figure 3, 50 equally spaced initial conditions covering the wave phase are used at each depth to compute the Lagrangian drift. We compare results using first-, third- and fifth-order velocity fields from (3.10), as well as a fourth-order Stokes drift from the Lagrangian approach of Blaser, Lenain & Pizzo (Reference Blaser, Lenain and Pizzo2025) (see also Clamond Reference Clamond2007), and which we shall return to for unsteady cases later. The vertical axes show depth below $z=0$ , while the horizontal axes show drift velocity $u$ ; for labels 1, 2, 3 and L4, this is formally $u_L$ (although $u_E=0$ means $u_L=u_S$ ) while label SD shows the leading order approximation to $u_S$ given by (3.9).

Figure 3.

Comparison of horizontal drift velocities with depth beneath monochromatic waves with $k=1$ m⁻¹ and (a) $H=0.3$ , (b) $H=0.45$ and (c) $H=0.6$ . Particle trajectories are obtained at first, third and fifth order from integration of the particle trajectory ODEs and yield Lagrangian drift velocities (labels 1, 3, 5). These are compared with the Stokes drift approximation (3.9) (SD) and the fourth order Lagrangian solution (Blaser et al. Reference Blaser, Lenain and Pizzo2025) (L4).

It is noteworthy that the integration of the first-order velocity field (red curves, label 1) gives a significant overestimate of the Lagrangian drift with depth, while the approximation (3.9) based on that same velocity field (blue curves, label SD) gives much better agreement with both the fourth-order Lagrangian result and higher-order Eulerian velocity fields. Only for high steepness $kH/2=0.225$ and $0.3$ , shown in panels (b) and (c), is the approximate Stokes drift noticeably different from the higher-order solutions, although its asymptotics at large depth hew close to the first-order theory, while the higher-order solutions have a somewhat different behaviour at large depth (see inset in panel c). Nevertheless, these results bear out the fact that – for steady, monochromatic waves – the Stokes drift is very well approximated by the leading order (quadratic) contribution (3.9).

3.2. Lagrangian drift at the surface

It appears that the most natural way to obtain surface drift values from an approximate Eulerian theory is by solving for the horizontal particle displacement from

(3.11) \begin{align} x'(t)=u(x,\zeta, t) \end{align}

(see Grue & Kolaas Reference Grue and Kolaas2020, who employ this to second order to compute the Lagrangian period) where the right-hand side can be approximated as

(3.12) \begin{align} u(x,\zeta, t) = u(x,0,t) + \zeta u_z(x,0,t) + \frac {\zeta ^2}{2} u_{zz}(x,0,t) + \cdots , \end{align}

making use of the expansion of the free surface

(3.13) \begin{align} \zeta & = a \cos \xi \left [ 1 + \frac {1}{8} a^2 k^2 + \frac {121}{192} a^4 k^4 \right ] + a \cos 2 \xi \left [ \frac {1}{2}ak + \frac {5}{6}a^3 k^3 \right ] \nonumber \\[5pt]& \quad + a \cos 3 \xi \left [ \frac {3}{8}a^2 k^2 + \frac {171}{128}a^4 k^4 \right ] + a \cos 4 \xi \left [\frac {1}{3} a^3 k^3 \right ]+a \cos 5 \xi \left [ \frac {125}{384}a^4 k^4 \right ] + \cdots \end{align}

together with the potential given in (3.10).

Results are shown in figure 4, which compares (3.9) evaluated at $z_0=0$ (label SD) with solutions of (3.11) up to fifth order, as well as the exact surface drift obtained by Longuet-Higgins (Reference Longuet-Higgins1987, figure 2) (label LH). Taking $k=1$ m⁻¹, the horizontal axis denotes half the actual crest-to-trough height – not the leading-order amplitude $a$ that is used in the perturbation expansion – which must be adjusted every odd order to obtain a desired value of $H$ (see discussion in Appendix A). Without such an adjustment, the wave height grows with the order of the expansion, as does the surface drift.

Figure 4.

Plots of Lagrangian surface drift velocity $u_L$ for monochromatic waves of varying steepness $Hk/2$ , using the surface velocity from second to fifth order (2–5), compared with the exact solution of Longuet-Higgins (LH) and the approximate Stokes drift $u_S$ (3.9) (SD).

Figure 4 shows that for waves of low steepness (below $Hk/2=0.2$ ) the differences in the formulations are essentially negligible. The Stokes drift formula (3.9) gives a drift that is slightly too small (as shown by Longuet-Higgins (Reference Longuet-Higgins1987) and complementary work using the Lagrangian formulation, such as Clamond (Reference Clamond2007)), but higher-order theories provide good agreement up to very high steepness (though the trend clearly does not continue to Longuet-Higgins’ steepest wave with $kH/2=0.44316$ ). We note that dispersion corrections appearing at the third and fifth order have the effect of decreasing surface drift slightly, which we comment upon later. The seemingly excellent agreement between the second-order approximation and the exact solution at slopes above $Hk/2=0.35$ should be considered accidental.

3.3. Remarks on the higher-order contributions

In the formulation of the velocity potential used here (corresponding to (A1a ) in Appendix A), the only difference between the first and third order is only the inclusion of nonlinear dispersion, so that third-order waves travel with a characteristic coordinate $\xi = kx - \varOmega t$ for

(3.14) \begin{align} \varOmega = \omega (k) \left (1 + \frac {a^2 k^2}{2} \right ) \!, \end{align}

( $a$ here is related to the wave height $H_1$ , as described in (A3)). No additional harmonics appear until fourth order, so the difference between curves 1 and 3 in figure 3 is entirely a consequence of this dispersion correction. Indeed, the procedure used to derive (3.9) can be applied without alteration to the third-order potential $\phi$ i.e. by calculating $\overline {\Delta \boldsymbol{x}^\intercal \boldsymbol{\nabla }\phi }$ , yielding

(3.15) \begin{align} {u_S} = \frac {a^2 \omega ^2 k e^{2kz_0}}{\varOmega } = \frac {2 a^2 \omega k e^{2kz_0}}{2+a^2 k^2}. \end{align}

Because $\varOmega \gt \omega$ , it is clear that the drift velocity (3.15) is generally smaller than (3.9), as observed also from integrating the particle trajectory mapping.

We can summarise the results for monochromatic waves as follows: the approximate Stokes drift formula (3.9) gives results that are somewhat too small at the surface and somewhat too large at depth. The differences for waves of moderate steepness, however, are not sizeable, as borne out also by experimental work (Paprota & Sulisz Reference Paprota and Sulisz2018). In our description, the only effects of nonlinearity are in adjusting the frequency (which has consequences for the entire flow) and in the appearance of higher-harmonics (whose effect is confined close to the surface due to the $\exp (nkz)$ -terms, see (3.10)). We shall now see how these effects appear in unsteady (bichromatic) wave trains.

4.1. A Stokes drift approximation from second-order wave theory

The second-order solution (4.1) lends itself to a calculation of a modified Stokes drift in a rather straightforward manner: calculate $(u,w)$ from (4.1), Taylor expand about a point $(x_0,z_0)$ , keep only the constant terms in the Taylor expansion and calculate approximate trajectories $(x(t),z(t))$ , subsequently substitute these trajectories into the linear terms in the Taylor expansion and integrate the resulting particle trajectory equations in time. The result of this procedure is

(4.3) \begin{align} u_S = \underbrace {a_1^2 k_1 \omega _1 e^{2 k_1 z_0} + a_2^2 k_2 \omega _2 e^{2 k_2 z_0}}_{(\text{I})} + \underbrace {a_1^2 a_2^2 \omega _1^2 \frac {(k_1-k_2)^3}{\omega _1-\omega _2} e^{2(k_1-k_2)z_0}}_{(\text{II})}, \end{align}

which consists of the sum of the Stokes drifts of each mode separately, term (I), as well as a new term (denoted term (II)) that arises from the difference harmonics (recall that we assume $k_1\gt k_2$ ). The author is grateful to one of the reviewers for pointing out that, while this paper was under review, an article by Liao & Zou (Reference Liao and Zou2025) appeared which also derives this term, § 1 in their Appendix C. Just as the Stokes drift derived from linear theory includes quadratic terms, so the modified Stokes drift derived from second-order theory now includes (formally) quartic terms.

Figure 6.

Illustration of the bichromatic Stokes drift approximation (4.3) for $k_2 = 1$ m⁻¹, $\epsilon _1=\epsilon _2=0.075$ and variable $k_1\gt k_2$ , shown at a depth (a) $z_0=0$ m, (b) $z_0=-4$ m and (c) $z_0=-5$ m. Circles in panels (b) and (c) show the results of computing the drift by integration of the particle paths using the first-order and second-order velocity field for $k_1=1.2$ m⁻¹ at each depth.

While overall particle motion decreases as we descend into the fluid, there are regions where the second-order contribution to Stokes’ drift has a sizeable effect, as we show in figure 6, which plots the terms (I) and (II) of (4.3) at three depths: (a) $z_0=0$ , (b) $z_0=-4$ and (c) $z_0=-5$ for a fixed wave $k_2=1$ m⁻¹. In connection with this, it is important to recall that the formulation (4.3) formally assumes deep water also for the difference harmonic $k_1-k_2$ , i.e. that $(k_1-k_2)h$ is large.

It is simplest to analyse the relative importance of the terms of (4.3) by setting $k_1=1+\Delta$ and varying $\Delta \gt 0.$ We note that the contribution (II) at the top of the fluid domain $z_0=0$ grows monotonically with $\varDelta$ , as can be observed in figure 6(a) (note that the steepness of both modes is fixed at $\epsilon =0.075$ ). By contrast, for $z_0\lt 0$ , the quartic term (II) has a single maximum and subsequently decays, as shown in figure 6(b–c). As we descend into the fluid, this maximum moves towards smaller values of $k_1.$ At certain depths and for certain bichromatic waves, the second-order contribution (II) dominates the first-order contribution (I), as shown around $k_1=1.2$ in panel (c). These contributions at depth remain several orders of magnitude smaller than the surface Stokes drift, although their significance (or lack or significance) should be assessed based on the time scales of interest and over the entire water column (see § 5.3).

4.2. Lagrangian drift for higher-order bichromatic waves

In attempting to compare the formula (4.3) with other methods of obtaining the wave-induced drift from wave theory, the principal difficulty we encounter – for the first time with bichromatic waves – is connected to the unsteadiness of the flow field and the lack of a unique solution towards which different expansions are known to converge. Unlike the monochromatic waves of § 3, we cannot simply compare waves of the same length and height, which differ geometrically only in the curvature of the profile between crest and trough. An additional complication arises at third order with the appearance of an asymmetric frequency correction (for $k_1\gt k_2$ ) first found by Longuet-Higgins & Phillips (Reference Longuet-Higgins and Phillips1962),

(4.4) \begin{align} \varOmega _1 &= \omega _1 \left [ 1+ \frac {1}{2} \epsilon _1^2 + \epsilon _2^2 \left ( \frac {k_1}{k_2} \right )^{1/2} \right ] \!, \\[-12pt] \nonumber \end{align}

(4.5) \begin{align} \varOmega _2 &= \omega _2 \left [ 1+ \frac {1}{2} \epsilon _2^2 + \epsilon _1^2 \left ( \frac {k_2}{k_1} \right )^{3/2} \right ]\!. \\[9pt] \nonumber \end{align}

This means that bichromatic waves at third-order (and fifth, seventh and so on) will be out of phase with their lower-order counterparts, as seen in figure 1.

This is compounded by a paucity of results (theoretical or experimental) on drift associated with bichromatic waves, making comparisons difficult. However, we can employ the methodology tested for monochromatic waves in § 3 to obtain results for bichromatic waves at the surface and below the trough level, and compare these with the approximate Stokes drift derived in (4.3). In this connection, it is important to emphasise that the Stokes drift velocity $u_S$ is still defined as the difference between Eulerian and Lagrangian mean velocities. In dealing with periodic, bichromatic wave groups, our Eulerian mean velocity (with the mean taken over the period of the group) vanishes exactly as for monochromatic waves. The terms of $u_S$ in (4.3) are therefore an approximation to the Lagrangian mean velocity which takes difference harmonics into account.

4.2.1. Lagrangian drift with depth

Below the surface, exactly as shown for monochromatic waves, it is possible to use the particle trajectories obtained by integrating the Eulerian velocity field to obtain the Lagrangian drift. This is put to the test in figure 7, where we select a bichromatic wave train with significant subharmonic drift as per (4.3), namely $k_1=1.2$ and $k_2=1$ m⁻¹. This corresponds to the circles shown in figure 6.

Figure 7.

Comparison of horizontal drift velocities with depth for a bichromatic wave train with $k_1=1.2$ and $k_2=1$ m⁻¹ with $\epsilon _1=\epsilon _2=0.075,$ as in figure 1.

Figure 7 compares the Stokes drift formulations (I) and (I)+(II) of (4.3) with the fourth-order Lagrangian drift of Blaser et al. (Reference Blaser, Lenain and Pizzo2025), and the Lagrangian drift obtained from integration of the particle trajectory mapping (3.4) with first-, second- and third-order velocity fields. The second-order horizontal velocity can immediately be read off from (4.1), while the third-order velocity, for which a simple algebraic expression has not been found, is obtained from the Hamiltonian expansion of § 2, (2.29).

The Lagrangian drift obtained from integrating the first-order theory overlaps completely with the classical Stokes drift (SD I). However, there is very good agreement among second- and third-order theories, the improved Stokes drift formulation (SD I + II), and the fourth-order Lagrangian theory (L4). This points to the fact that the approximation (4.3) captures important features of the bichromatic Lagrangian drift. The same comparison for a pair of more widely separated wavenumbers (such as $k_1=2$ and $k_2=1$ m⁻¹, not shown) would demonstrate that all six curves essentially coincide, as might be predicted from figure 6 (noting that there are no notable flows induced by the difference harmonic in this case).

4.2.2. Lagrangian drift at the surface

To obtain drift at the surface, we integrate the equation for the horizontal particle position $x'(t)=u(x,\zeta, t),$ where the velocity field is given by Taylor expansion about $z=0$ to a given order, exactly as in § 3.2 for monochromatic waves. In addition to the velocity fields to third order, we also require the second-order bichromatic free surface, which is given in (2.17).

Figure 8.

Comparison of surface drift formulations for a bichromatic wave train with (a) $\omega _1=1.25, \omega _2=1$ rad s⁻¹ and (b) $\omega _1=2, \omega _2=1$ rad s⁻¹ and identical (linear) steepness $\epsilon _1=\epsilon _2$ . SD I and SD I + II denote the respective terms in (4.3). These are compared with the second-order and third-order bichromatic solution, as well as the fourth-order Lagrangian drift (L4) of Blaser et al. (Reference Blaser, Lenain and Pizzo2025).

Figure 8 compares the Lagrangian surface drift of second- and third-order theory with the classical Stokes drift (SD I) and the modified Stokes drift (SD I + II) evaluated at $z_0=0.$ In panel (a), where the two waves are close in frequency, we expect the contributions of term (II) to be negligible. Instead, the Lagrangian drift at the surface is influenced by sum-harmonics appearing in the second- and third-order theory, and which are not accounted for by the approximation (4.3).

When the constituent harmonics are more widely separated, such as in the case of $\omega _1=2$ and $\omega _2=1$ rad s⁻¹ shown in panel (b), the difference harmonic contributions of term (II) are more prominent. In both cases, the second- and third-order theories give an enhanced drift, which is consistent with the results obtained for monochromatic waves and shown in figure 4.

Figure 8 also compares with the fourth-order Lagrangian drift obtained from the third-order bichromatic solution in Lagrangian variables obtained in Appendix B of Blaser et al. (Reference Blaser, Lenain and Pizzo2025), which generalises Pierson’s (Reference Pierson1961) classical second-order result. Unfortunately, bichromatic waves in Eulerian and Lagrangian coordinates are not the same; rather, the wave profiles are geometrically distinct, with a phase-shift that occurs between successive crests of both solutions, as already noted by Fouques & Stansberg (Reference Fouques and Stansberg2008). (This phase shift is not due to differences in dispersion relation, as might be hypothesised from Pierson (Reference Pierson1961), since Blaser et al. ( Reference Blaser, Lenain and Pizzo2025) recover the Eulerian dispersion correction (4.4)–(4.5) at third order.) Whether it is possible to adapt the Lagrangian perturbation expansion to match the Eulerian solution (or vice versa) remains an open question. For well-separated frequencies (panel b), the Lagrangian drift obtained from our Eulerian particle trajectory mapping (3.11) and the Lagrangian solution of Blaser et al. (Reference Blaser, Lenain and Pizzo2025) remain remarkably close, while for two nearby frequencies (panel a), the notable difference grows with wave slope.

4.3. Difference harmonics and particle drift

In general, we see the importance of leading-order effects near the surface of bichromatic waves, while second-order effects become more prominent at depth. This is no surprise and follows from the structure of the difference harmonic terms themselves. As a prelude to our work on wave groups, it is expedient to investigate their role in bichromatic seas.

We first consider the velocity field which drives the particle trajectory mapping and which can be calculated from the second-order (4.1) or third-order potential (2.29). In fact, figure 5 makes clear that the third-order velocity field – though difficult to calculate – makes only a small contribution. In large part, this is due to the third-order dispersion correction (4.4)–(4.5), as can be seen by comparing panels (b) and (d) with panels (c) and (e) in figure 5. These dispersion corrections do not affect the instantaneous velocity field, only its time evolution.

As an example, we consider the waves studied in figure 6 with $k_1=1.2$ and $k_2=1$ m⁻¹, and focus only on horizontal velocities. When only the first-order contributions $u_1$ are retained, we find the familiar pattern of forward velocity below the crests and backward velocity below the troughs, with an exponential fall-off in depth as depicted in figure 9(a). If instead we focus only on the second-order contributions $u_2$ , the situation looks quite different. The exponential decay in velocities is slower because $\exp ((k_1-k_2)z)$ is larger than either $\exp (k_1 z)$ or $\exp (k_2 z)$ for $k_1 \sim k_2,$ and we find a negative horizontal velocity under the centre of the bichromatic wave train, as shown in figure 9(c). The full second-order picture is a sum of the first- and second-order terms, shown in panel (b), and it is clear that the first order dominates near the surface, while the second-order difference harmonic dominates at depth. This can be observed clearly by considering the zero contour lines of the horizontal velocity, shown in black for each case.

Figure 9.

Horizontal velocity for a bichromatic wave train with $k_1=1.2$ and $k_2 = 1$ m⁻¹, shown indicatively in the top panel. (a) First-order velocity only. (b) Sum of first- and second-order horizontal velocities. (c) Second-order horizontal velocity only. Black curves are contours of vanishing horizontal velocity.

The generic situation for particle paths in a bichromatic wave train is thus as follows: close to the surface, the first-order contributions dominate, giving the curlicue particle paths observed in figure 5(b). The effect of difference harmonics is to retard slightly the forward motion of the particles, due to the slow moving (in deep water, the group velocity is half the phase velocity) region of backward (negative) velocity under the highest crests observed in figures 5(e) and 9.

Deeper into the fluid column, the effect of the first-order free-modes becomes negligible. The fluid motion is dominated by the difference harmonic term $k_1-k_2$ , and the particle paths behave essentially as those of a monochromatic wave with that wavenumber. This is clearly illustrated in figure 5(e), where we see that the modes $k_1$ and $k_2$ that make up the first-order theory are only a small perturbation on the displacement due to the difference harmonic term.

It is worth pointing out the existence of a region beneath the bichromatic wave train where the particle displacement is largely against the direction of wave propagation, as seen between $t\approx [-5.2,5.2]$ s in figure 5(e) (see markers ‘+’). This should not be surprising, since it simply reflects the effective motion of the difference harmonic term, as shown in panel (d). Despite the existence of this ‘return flow’ – which we shall encounter again in our discussion of wave groups – the average particle drift over the period of the bichromatic group is in the direction of wave propagation at every depth in the fluid, as shown in figures 6 and 7. Because the solution is periodic, over many periods, the net effect of difference harmonics is to enhance forward drift, as shown in figure 7.

5.1. Focused linear wave groups

In treating waves with an increasing number of linear harmonic constituents, we also have a plethora of parameters that can be tuned. One case of interest, from both an experimental and a theoretical standpoint, is that of a (focused) wave group. Linear focusing is quite straightforward: if our free surface is a superposition of sinusoids as per (3.1), then writing

(5.1) \begin{align} \zeta _1 = \sum _i a_i \cos (k_i(x-x_0) - \omega _i(t-t_0)) \end{align}

readily gives a wave train that focuses – i.e. all components come into phase – at $(x_0,t_0)$ . Pulse-like wave groups can be created by suitable tuning of the amplitude spectrum $a_i.$ However, for any finite number of Fourier modes, the group so created will never be completely localised; depending on the choice of frequencies $\omega _i$ or wavenumbers $k_i$ , it may or may not be periodic in time or space, and the extent to which it appears to have a well-formed envelope depends on the choice of harmonics and amplitude spectrum.

Figure 10.

A crest-focused wave group with $k_p=1.5, \sigma =0.18$ and ten Fourier modes, with $S=0.24.$ Panel (a) shows the free surface in first-, second- and third-order theory. Panels (b) and (d) show the respective particle trajectories with initial position $(x_0,z_0)=(0,-0.2)$ and $(0,-3),$ beginning at time $t=-20$ prior to the passage of the group. Panels (c) and (e) show the horizontal particle position of panels (b) and (d) with time $t,$ illustrating the negative drift at depth also seen in figure 5.

One possible choice is a Gaussian amplitude spectrum of the form

(5.2) \begin{align} a_i = A \exp \left ( -\frac {(k_i-k_p)^2}{2 \sigma ^2}\right )\!. \end{align}

An example of a wave train with such an amplitude spectrum is shown in figure 10(a), where we have chosen $10$ equidistant modes between $k_1=1$ and $k_{10}=2$ m⁻¹, $k_p=1.5$ , $A=0.04$ , and $\sigma =0.18.$ The focusing location is chosen to be $(x_0,t_0)=(0,0)$ , and the focus occurs at a crest. Following Blaser et al. (Reference Blaser, Lenain and Pizzo2025), we define

(5.3) \begin{align} S = \sum _{n=1}^N a_n k_n \end{align}

as the maximum slope for the linear focused wave and use this as a measure of nonlinearity going forward.

The particle paths in figure 10(b,d) are in many ways similar to those in the simpler, bichromatic wave group shown in figure 5. Towards the top of the fluid column, first-order theory dominates, while at depth, the effect of the in-phase subharmonics associated with higher-order theories clearly accounts for the bulk of the particle motion. To further illustrate this point, panels (c) and (e) show only the horizontal component of the particle position as a function of time $t$ , as the wave group passes over the initial position. Towards the surface in panel (c), there is nearly no motion until $t=-10$ s, while at depth, we observe a steady forward drift in the positive $x$ -direction. Under the centre of the group, as expected, we see a strong forward drift towards the surface (panel c) and backward (return) flow at depth (panel e).

This effect can be seen even more clearly in figure 11, which illustrates (for a single particle) the depth dependence of the horizontal motion. The (leading order) free surface of a focused Gaussian wave train with $k_p=2.5$ m⁻¹ is shown in panel (a), in both the crest focused (blue, solid curve) and trough focused (red, dash-dotted curve) case. Panel (b) shows the horizontal displacement $\delta x$ of a particle as the wave group passes it. Higgins, van den Bremer & Vanneste (Reference Higgins, van Den Bremer and Vanneste2020) suggest that the net (long-time) displacement of the return flow to second order can be calculated by integrating the second-order velocity $u_2$ in time (see their (2.10)), and we see (panel b, black dotted curve) that this indeed matches our second-order displacement (shown in red) quite well at depth.

Figure 11.

Indicative horizontal drift $\delta x$ beneath a crest/trough focused wave train. Panel (a) shows a crest focused (blue) and trough focused (red, dash-dotted) Gaussian wave group with $k_p=2.5$ m⁻¹, $\sigma =0.55$ and 10 equispaced modes between $k_1=1$ and $k_{10}=4$ m⁻¹, with steepness $S=0.3.$ Panel (b) shows the horizontal (Lagrangian) displacement $\delta x$ of a particle as the wave group passes over it (from $t=-8.7$ to $t=8.7$ s, as shown in panel a). Solid lines show the drift for the crest-focused wave, dash-dotted lines the corresponding trough-focused values. Black dots refer to the second-order net displacement of the return flow calculated from (2.10) of Higgins et al. (Reference Higgins, van Den Bremer and Vanneste2020).

Interestingly, there is a marked difference in the return flow between crest- and trough-focused cases, with even first-order theory showing some return flow for trough-focused waves (dash-dotted curves). Such (localised) return flows are not captured by Stokes drift approximations, whether based simply on the linear theory or incorporating second-order effects via the extension of (4.3),

(5.4) \begin{align} u_S = \underbrace {\sum _j a_j^2 \omega _j k_j e^{2k_j z_0}}_{\text{(I)}} + \underbrace {\sum _{k_i\gt k_j} \frac {\omega _i ^2 a_i^2 a_j^2 (k_i-k_j)^3}{\omega _i-\omega _j} e^{2(k_i-k_j)z_0}}_{\text{(II)}}. \end{align}

Indeed, the effects shown in figure 11 are highly phase-dependent and rely on suitable choice of position in the spatio-temporal evolution of the group, while the approximate formulae (5.4) can only capture an averaged drift based on the Fourier amplitudes, but independent of the phases. This contrasts strongly with the theoretical result of Higgins et al. (Reference Higgins, van Den Bremer and Vanneste2020) for a spatially localised packet.

As we have done for monochromatic and bichromatic waves, we can apply second- and third-order, deep-water theory directly to calculate the Lagrangian surface drift of focused wave groups, using the particle trajectory mapping at the surface. The results of such a computation are reported in figure 12, which is the counterpart of figures 4 and 8. In this case, we have initialised our wave group with a central frequency $\omega _c=2 \pi$ , and using $N=12$ modes, have defined the amplitudes $a_n=S(Nk_n)^{-1},$ where the frequencies are chosen via $\omega _n=\omega _c(1+\Delta (n-N/2)/N)$ for bandwidth $\Delta =0.77$ and $k_n$ are the wavenumbers obtained from the linear dispersion relation. This is precisely the formulation used by Blaser et al. (Reference Blaser, Lenain and Pizzo2025), and allows us to compare our results with their experimental and numerical values.

Figure 12.

Comparison of surface Lagrangian drift beneath a wave group using second- and third-order theories, as in figures 4 and 8, compared with the leading-order approximation to the Stokes drift (SD I). The dimensionless bandwidth parameter $\Delta =0.77$ and steepness $S$ is varied from 0 to 0.3 to compare with figure 3(a) of Blaser et al. (Reference Blaser, Lenain and Pizzo2025). Circles ( $\bullet$ ) denote their fully nonlinear numerical simulations (with $\Delta =0.8$ ) and triangles ( $\blacktriangle$ ) denote the mean experimental result.

Following the trends for monochromatic and bichromatic waves, we note that the lowest order Stokes drift evaluated at $z_0=0$ (SD I) and which treats the Stokes drift as a sum of the drifts of each individual mode (3.9) gives the smallest value. The second- and third-order formulations show good agreement with both laboratory experiments ( $\blacktriangle$ ) and fully nonlinear simulations ( $\bullet$ ) up to rather high steepness $S$ , even though our formulation neglects amplitude evolution which is included in the numerical results (see also the discussion in Appendix C).

A further issue which we encounter in these calculations is the lack of joint spatio-temporal localisation of our wave groups. This fact obscures the scale over which a ‘mean’ quantity such as $u_E$ should be calculated and implies that we cannot always numerically enforce $u_E \equiv 0$ by suitable construction of the group (see the discussion in § 4). Thus, while it is possible to compare the Lagrangian drift obtained from integration of particle paths with the Lagrangian drift obtained computationally and experimentally (as in figure 12), the difference between Lagrangian drift and (a priori unknown) Stokes drift is an Eulerian mean flow which is sensitive to the domain over which the averaging is taken.

5.2. Random multichromatic waves

Random waves are the closest mathematical idealisation to the waves we might typically find in the ocean: the phases of one Fourier component are not correlated with the next and the resulting pattern – while occasionally exhibiting a distinct group – is highly irregular. Random phases at lowest order also means that sum and difference phases at higher orders are random, due to a neglect of slow phase-coherence effects associated with nonlinear coupling (Andrade & Stuhlmeier Reference Andrade and Stuhlmeier2023). In such cases, the averaging necessary to unambiguously define an Eulerian mean velocity is not achievable. At best, we can average over a suitably large number of individual waves in the expectation that the Eulerian mean velocity over such a large time will be small.

Where focused wave groups exhibit localised features, such as a clear division between flows in the direction of propagation at the surface and return flows at depth – which have been well described by using wave-envelope formulations by van den Bremer and co-workers (van den Bremer & Taylor Reference van den Bremer and Taylor2016; van den Bremer & Breivik Reference van den Bremer and Breivik2018; van den Bremer et al. Reference van den Bremer, Whittaker, Calvert, Raby and Taylor2019; Calvert et al. Reference Calvert, Whittaker, Raby, Taylor, Borthwick and Van Den Bremer2019) – the lack of clear ‘groupiness’ in random seas means that we should expect these localised effects to be averaged out. The advantage is that we can test whether (on average) the Stokes drift formulation (5.4) is suitable, since it clearly cannot capture spatio-temporally localised flows around a group itself (all terms are positive).

Figure 13.

Drift velocity with depth below a random wave field, calculated from the average of 200 realisations (shown as lighter curves) of an irregular sea surface with random phases. The horizontal Lagrangian drift velocities are calculated from the velocity fields as in figure 7 using first-order, second-order and third-order theory, and compared with the Stokes drift approximations of (5.4) (SD I) and (SD I + II).

Figure 13 shows the results of computing the wave-induced particle displacement for 200 realisations of a random sea surface. A Gaussian amplitude spectrum (5.2) with $A=0.05$ and the parameters $k_p=2.5$ m⁻¹ and $\sigma =0.7$ are used, but rather than assigning the phases so as to achieve linear focusing at a specific location, these are randomly generated for each realisation. Blue curves denote first-order theory, red curves second order, and yellow curves third order, with the prominent, darker curves representing the average over all realisations for a given order. In each case, the integration is for approximately 80 peak periods, which is long enough to clearly distinguish the particle drift.

With a peak wavenumber $k_p=2.5$ m⁻¹, most of the motion due to the free-wave constituents is filtered out at a depth of $z=1.25$ m. The linear theory gives a very narrow range of displacements, with the realisations in red barely showing deviation from the average and agreeing well with the linear Stokes drift (SD I). Second- and third-order displacements are quite large in some realisations, depending on the phase relationships involved. However, as expected, there is no systematic return flow at any depth for random waves, because such (spatially and temporally) localised flows are averaged out to yield the net forward drift which we observe in figure 13. Moreover, all theories approach zero from above with depth, exhibiting the same ordering we have hitherto observed: linear theory predicting the smallest drift followed by second and third order. The correction (II) in (5.4) is found to give a useful improvement over using only the classical Stokes drift term (I).

5.3. Parametric wave spectra

The comparisons of the previous sections, based on explicit computation at different orders, support the inclusion of quartic difference harmonic terms in calculations of the Stokes drift throughout the water column. In the absence of information about Eulerian currents in the open ocean, the Stokes drift is the constituent of the Lagrangian drift which can be estimated directly from the wave field. As we have seen, the difference harmonic constituents provide a useful approximation to this Stokes drift. For many applications in ocean science and engineering, it is advantageous to present such a result in terms of energy rather than amplitude spectra, and to explore its consequences for some common parametric spectral shapes.

For a continuous wavenumber spectrum $E(k)$ , we make the identification with the modal amplitudes ${a_n^2}/{2} =: E_n = E(k) \,{\rm d}k$ , where the uniform grid spacing ${\rm d}k$ is assumed to tend to zero. With this assumption, it is straightforward to rewrite the discrete multi-chromatic Stokes drift (5.4) with second-order contributions as

(5.5) \begin{align} u_S = 2 \int _0^\infty E(k) \omega (k) k e^{2kz} \,{\rm d}k + 4 \int _0^\infty \int _{k'}^\infty \frac {\omega ^2(k) E(k) E(k') (k-k')^3}{\omega (k)-\omega (k')} e^{2(k-k')z} \,{\rm d}k \,{\rm d}k', \end{align}

where the first term is identical to that found by Kenyon (Reference Kenyon1969) and the second term comes from the inclusion of second-order terms in deep water.

Numerous critical discussions have focused on the key role of the high-frequency cut-off on the Stokes drift near the surface, including Breivik, Janssen & Bidlot (Reference Breivik, Janssen and Bidlot2014) and more recent work by Clarke & Van Gorder (Reference Clarke and Van Gorder2018) and Lenain & Pizzo (Reference Lenain and Pizzo2020). As Clarke & van Gorder (Reference Clarke and Van Gorder2018) point out, if the Stokes drift is rewritten in terms of frequency $u_S \sim \int E(\omega ) \omega ^3 \,{\rm d} \omega,$ and the high-frequency tail of the spectrum is such that $E(\omega )\sim \omega ^{-n}$ , then we have a divergent expression for $n\leq 4$ . Their proposed resolution, introducing a wave breaking frequency as an upper limit for spectral Stokes drift calculations, was found to be reasonable for the subsurface Stokes drift by Lenain & Pizzo (Reference Lenain and Pizzo2020), but for the following calculations, we will employ the maximum wavenumber $k_M$ suggested by Lenain & Melville (Reference Lenain and Melville2017) and used by Lenain & Pizzo (Reference Lenain and Pizzo2020).

As a first comparison, we consider the Phillips spectrum for the equilibrium range of wind-generated waves, given in terms of radian frequency as

(5.6) \begin{align} E(\omega ) = \begin{cases} \alpha g^2 \omega ^{-5} &\text{ for } \omega \gt \omega _p,\\ 0 &\text{ for } \omega \leq \omega _p, \end{cases} \end{align}

where we take $\alpha =0.0083$ as in Breivik et al. (Reference Breivik, Janssen and Bidlot2014). We relate the peak frequency $\omega _p$ to a friction velocity $u_*$ (see Holthuijsen Reference Holthuijsen2007) to establish $k_M$ , and show the results in figure 14.

Figure 14.

Comparison of Stokes drift formulations for a Phillips spectrum with $T_p=10$ s. Panel (a) compares the Stokes drift near the surface up to the e-folding depth of the peak wavenumber ( $k_p=0.04$ m⁻¹). Panel (b) compares the Stokes drift formulations at depths below one peak wavelength.

Panel (a) shows the Stokes drift near the surface using either the first term in (5.5) only (SD I) or both terms in the same equation (SD I + II). The difference between the two is sizable, but strongly dependent on the choice of high-frequency cut-off $k_M.$ This tracks with the very pertinent discussion of Lenain & Pizzo (Reference Lenain and Pizzo2020) (therein, they report an underestimation of the Stokes drift of up to 50 % from failure to account for high frequencies, but they do not consider the effect of difference harmonics). This difference persists as we descend into the fluid and a noticeable difference on the order of 1 % is visible at depths below one peak wavelength in panel (b). Defining the Stokes transport $V_S$ as

(5.7) \begin{align} V_S = \int _{-\infty }^0 u_S(z) \,{\rm d}z, \end{align}

we find nearly a 9 % change in total Stokes transport. If we calculate the Stokes transport only below of the $e$ -folding depth of the peak frequency, we still observe a change of approximately 2 % owing to the inclusion of difference harmonic terms (II).

The Phillips spectrum describes only the high-frequency tail, so it is of interest to consider a more realistic parametric spectrum: we next consider a Pierson–Moskowitz (PM) spectrum, corresponding to a fully developed sea under a constant wind over unlimited fetch. The case shown in figure 15 corresponds to a wind speeds of $U_{10}=$ 10 m s⁻¹ (with corresponding $U_{19.5}\approx 1.026 U_{10}$ ) and peak frequency $f_p$ of 0.17 Hz. Blue curves show the formula of Kenyon (Reference Kenyon1969) (corresponding to the first term in (5.5), labelled SD I) and red curves show both terms of (5.5) (labelled SD I + II). For reference, we note that the high-frequency cut-off $f_M$ obtained from $k_M$ is 6.3 $f_p.$ We note also that it is possible to cut off the low frequencies and evaluate the outer integrals in (5.5) starting from $k_p$ instead of 0 without materially affecting the results, as suggested by Breivik et al. (Reference Breivik, Bidlot and Janssen2016) and Lenain & Pizzo (Reference Lenain and Pizzo2020).

Figure 15.

Comparison of Stokes drift formulations for a PM spectrum with $U_{10}=10$ m s⁻¹. Panel (a) compares the Stokes drift near the surface up to the $e$ -folding depth of the peak wavenumber ( $k_p=0.1$ m⁻¹). Panel (b) compares the Stokes drift formulations at depths below one peak wavelength.

The differences in the formulations are clearly visible, both near the surface and at depth. At the surface, we find increases of surface drift in excess of 40 % when difference harmonic terms are accounted for, reducing to just less than 2 % at half the $e$ -folding depth of the peak wavenumber $k_p.$ Interestingly, starker differences reappear at greater depths, as shown in panel (b), which highlights the different asymptotics of the two formulations at depths of more than one peak wavelength. We find that keeping both terms in (5.5) leads to an increase in Stokes transport of nearly 3.3 % for $U_{10}=10$ m s⁻¹. The majority of this difference comes from the near-surface transport, with transport below the peak $e$ -folding depth accounting for a difference of only 0.1 %. The same analysis can be carried out for JONSWAP spectral shapes, with qualitatively very similar results to those found in figure 15.

We can also compare the Stokes drift formulations of (5.5) for other spectral shapes, such as swell-wave spectra described by Ochi & Hubble (Reference Ochi and Hubble1976). In figure 16, we show the most probable Ochi–Hubble shape corresponding to a significant wave-height $H_s=3$ m and the two Stokes drift calculations. The surface drift shows a change in excess of 20 % when subharmonic terms (II) are included and the Stokes transport $V_S$ is slightly more than 1 % larger over the entire water column.

Figure 16.

Comparison of Stokes drift formulations for an Ochi–Hubble spectrum with $H_s=3$ m. Panel (a) shows the spectral shape as a function of wavenumber $k$ , while panels (b) and (c) compare the use of the first term of (5.5) (SD I) and both terms of (5.5) (SD I + II).