2. Nonlinear waves with non-zero mean

The surface water-wave problem is generally described by the Euler equations with boundary conditions at the bottom and free surface. Since the fluid is incompressible and the flow is assumed irrotational, the problem may be formulated in terms of the velocity potential $\phi (x,z,t)$ and the surface elevation $\eta (x,t)$ . In two dimensions, the variables $x$ and $z$ denote the horizontal and vertical coordinates, respectively, whereas the variable $t$ measures the time evolution. The pressure is eliminated with help of Bernoulli’s equation and the potential $\phi$ satisfies Laplace’s equation

(2.1) \begin{equation} {\nabla} ^2 \phi = 0, \quad \textrm {for} \ -h \lt z \lt \eta (x,t). \end{equation}

At the free surface $z = \eta (x,t)$ , the kinematic and dynamic boundary conditions are formulated in terms of the velocity potential and the surface excursion by

(2.2) \begin{align} \eta _t+\phi _x\eta _x-\phi _z & = 0, \end{align}

(2.3) \begin{align} \phi _t+\frac {1}{2} \big ( \phi ^2_x+\phi ^2_z \big ) + g\eta & = 0. \end{align}

At the bottom $z = -h$ (set to be uniform), the no-flow boundary condition amounts to

(2.4) \begin{equation} \phi _z = 0. \end{equation}

Hereafter, subscripts denote partial derivatives, unless stated otherwise.

For steadily progressive waves propagating at speed $c$ in the $x$ -direction, (2.2)–(2.3) are rewritten as

(2.5) \begin{align} \left . \begin{aligned} -c \eta _x +\phi _x\eta _x-\phi _z & = 0, \\ -c \phi _x+\frac {1}{2} \big ( \phi ^2_x+\phi ^2_z \big ) + g \eta & = 0, \end{aligned} \right \}\ \textrm {on} \ z=\eta (x). \end{align}

Now suppose that we impose a solution by shifting $\eta$ by a small (indeed minute) constant amount $s$ , i.e. a setup if $s \gt 0$ (increase of the water level) or a setdown if $s \lt 0$ (decrease of the water level). This parameter $s$ may be viewed as the mean water level with respect to the coordinate system $(x,z)$ recognizing that, if $\tilde {\eta }(x) = \eta (x) + s$ , then

(2.6) \begin{align} \frac {1}{L} \int _0^L \tilde {\eta } \, {\rm d}x = s + \frac {1}{L} \int _0^L \eta \, {\rm d}x, \end{align}

where $L$ is the length of the domain under consideration or, in the context of the present paper, the wavelength $\lambda$ of a travelling wave. If the original free surface $\eta$ has zero mean, i.e. $\int _0^L \eta \, {\rm d}x = 0$ , then $s$ represents the shifted mean water level. The equations then become

(2.7) \begin{align} \left . \begin{aligned} -C \tilde {\eta }_x + \tilde {\phi }_x \tilde {\eta }_x - \tilde {\phi }_z & = 0, \\ -C \tilde {\phi }_x + \frac {1}{2} \big ( \tilde {\phi }^2_x + \tilde {\phi }^2_z \big ) + g \tilde {\eta } & = 0, \end{aligned} \right \}\ \textrm {on} \ z = \tilde {\eta }(x) = \eta (x) + s. \end{align}

The functions $\tilde {\eta }$ and $\tilde {\phi }$ are solutions of this slightly altered problem, now with a propagation speed $C$ . Let us assume that the original propagation speed $c$ and the new propagation speed $C$ are related by $C=c+V$ . Let us also assume that the only change in the velocity potential is an augmented horizontal speed, i.e. $\tilde {\phi } = \phi + V x$ . We can then use the second equation to write

(2.8) \begin{align} -C (\phi _x + V ) + \frac {1}{2} \phi _x^2 + \phi _x V + \frac {1}{2} V^2 + \frac {1}{2} \phi _z^2 + g \eta + g s = 0. \end{align}

Rearranging gives

(2.9) \begin{align} -( C - V ) \phi _x - CV + \frac {1}{2} \phi _x^2 + \frac {1}{2} V^2 + \frac {1}{2} \phi _z^2 + g \eta + g s = 0. \end{align}

We see that some terms in the last equation simplify if we assume that the original dynamic boundary condition is still approximately valid, and we are left with

(2.10) \begin{align} \frac {1}{2} V^2 -CV + gs = 0. \end{align}

Remembering that $C=c+V$ and reducing further leads to

(2.11) \begin{align} V^2 + 2cV - 2gs = 0. \end{align}

Solving this quadratic equation yields two roots

(2.12) \begin{align} V = -c \pm \sqrt {c^2 + 2gs}. \end{align}

We are only interested in the solution

(2.13) \begin{equation} V = -c + \sqrt {c^2 + 2gs}, \end{equation}

such that $V = 0$ when $s = 0$ in the unperturbed configuration ( $\tilde {\eta } = \eta$ , $\tilde {\phi } = \phi$ ). Asymptotically, for small $s$ (i.e. $g s /c \ll 1$ ), we have $V \sim g s/c$ . In subsequent sections, we will propose (2.13) as a correction to linear or nonlinear potential flow for estimation of the Stokes drift velocity and for simulation of particle dynamics in the presence of a wave setup or setdown. As it will appear, the corrected flow speed given by (2.13) captures the effect that a setup or setdown has on the background flow of a surface wave. If this background current $V$ is introduced together with a setup or setdown $s$ in a periodic travelling wave, the resulting combination yields a new approximate steady wave solution which cannot be readily obtained by other means. We also point out that a background current can be added to this nonlinear problem (e.g. Constantin & Strauss (Reference Constantin and Strauss2010), Hsu (Reference Hsu2013) and many others). However, here the current is strictly connected with the altered mean level of the surface wave profile.

Figure 1. This schematic elucidates the geometric setting of the problem. The undisturbed water level is indicated by a dashed line. A surface wave is indicated in light grey (initial time) and black (final time). Since this is a travelling wave, the two waveforms coincide exactly after one period. This wave features a setup, i.e. a positive mean level $s$ which is indicated by a red line. The wave height is $H$ , and it can be seen that a fluid particle located at the free surface stays there, but instead of cycling through the wave as in a typical Stokes wave, the particle experiences a significant forward drift.

3. Particle paths for linear waves and Stokes drift

As can be found in many authoritative works on fluid mechanics (e.g. Debnath Reference Debnath1994, Kundu & Cohen Reference Kundu and Cohen2008), the Stokes drift can be explained in the linear case by invoking the exponential decrease with depth of the wave-induced horizontal velocity field. As fluid particles trace out a nearly circular path in the fluid column, the upper part of the path experiences a stronger horizontal velocity than the lower part. This slight discrepancy leads to a small offset as the particle approaches the original position and results in a small drift in the direction of wave propagation (figure 1).

Consider a sinusoidal wave propagating on a fluid of depth $h$ . The deflection of the free surface for a wave of amplitude $a$ (wave height $H = 2a$ ) is assumed to be $\eta (x,t)=a \cos (kx-\omega t)$ . Solving the linearized Euler equations with bottom and surface boundary conditions shows that the velocity potential is

(3.1) \begin{equation} \phi (x,z,t) = \frac {\omega a}{k} \frac {\cosh [k(h+z)]}{\sinh (kh)} \sin (kx-\omega t). \end{equation}

Here, $k = 2 \pi / \lambda$ is the wavenumber, $\lambda$ is the wavelength, $\omega = 2 \pi / \tau$ is the angular frequency and $\tau$ is the wave period. The velocity field $\boldsymbol {u} = (u,v)^\top = \boldsymbol{\nabla }\phi$ is obtained by taking spatial derivatives of $\phi$ , as follows:

(3.2) \begin{equation} \begin{aligned} \phi _x(x,z,t) & = \omega a \frac {\cosh [k(h+z)]}{\sinh (kh)} \cos (kx-\omega t), \\ \phi _z(x,z,t) & = \omega a \frac {\sinh [k(h+z)]}{\sinh (kh)} \sin (kx-\omega t). \end{aligned} \end{equation}

Particle paths given parametrically by $(x,z) = (\xi (t), \zeta (t))$ can now be computed by solving the coupled system of ordinary differential equations (ODEs)

(3.3) \begin{equation} \begin{aligned} \frac {{\rm d}\xi }{{\rm d}t} & = \omega a \frac {\cosh [k (h+\zeta )]}{\sinh (kh)} \cos (k \xi - \omega t), \\ \frac {{\rm d}\zeta }{{\rm d}t} & = \omega a\frac { \sinh [k(h+\zeta )]}{\sinh (kh)} \sin (k \xi -\omega t). \end{aligned} \end{equation}

Figure 2. (a) Wave profile near the wave crest of a linear periodic wave for the original wave at time $t=0$ and after one period at $t=\tau$ . A particle located initially at the wave crest, and constantly staying at the free surface has drifted slightly forward. This is the Stokes drift. However, the particle has not yet reached quite the same height as in the beginning. (b) The wave profile after time $t=\tau _L$ is shown together with the path of the same particle. Here $\tau _L$ is based on the particle shown, and it is found by specifying that the particle reach the same height as where it started. Since $\tau _L \gt \tau$ , the wave profile has completed slightly more than a full cycle.

Figure 3. (a) Zoom-in of figure 2(a). (b) Zoom-in of figure 2(b). The dashed curve continuing the particle trajectory indicates that the particle is at the highest position at time $\tau _L$ .

Assuming that the particle position stays close to a centre $(x_0, z_0)$ allows for replacement of the position $(\xi (t), \zeta (t))$ on the right-hand side by the centre position, leading to closed elliptic orbits of the form

(3.4) \begin{equation} \begin{aligned} \xi (t) & = x_0 - a \frac {\cosh [k (h+z_0)]}{\sinh (kh)} \sin (k x_0 - \omega t), \\ \zeta (t) & = z_0 + a \frac {\sinh [k (h+z_0)]}{\sinh (kh)} \cos (k x_0 - \omega t). \end{aligned} \end{equation}

If the velocity field is not evaluated at the centre point $(x_0, z_0)$ , but it is assumed that the difference between the particle position and the centre point stays small, i.e. $|x-x_0|$ and $|z - z_0|$ are small during one wave cycle, then the second-order approximation of the paths yields a net movement of fluid particles in the direction of wave propagation. This net transport of a particle is the Stokes transport

(3.5) \begin{equation} x_L = a^2 \omega k \tau \frac {\cosh [2k (h+z_0)]}{2 \sinh ^2(kh)}, \end{equation}

which can be obtained analytically by integrating over one period $\tau = 2\pi /\omega$ (Debnath Reference Debnath1994). A mathematically rigorous treatment of the Stokes drift for linear water waves is given by Constantin & Villari (Reference Constantin and Villari2008). Another quantity which is often used is the Stokes drift velocity

(3.6) \begin{align} u_L = \frac {x_L}{\tau }. \end{align}

Furthermore, it is well known that the time for a fluid particle to complete one cycle is not exactly equal to the wave period, but (in the inviscid theory) a bit longer (see figures 2 and 3). This period is known as the Lagrangian period $\tau _L$ and varies with depth, as can also be understood if the linear surface water-wave problem is solved in the Lagrangian framework (Lamb Reference Lamb1932). If particle trajectories are obtained numerically from solving (3.3) without any further asymptotic approximation, the Lagrangian period $\tau _L$ may be found and the Stokes drift velocity then takes the form

(3.7) \begin{equation} u_L = \frac {x_L}{\tau _L}. \end{equation}

However, the difference between $\tau _L$ and $\tau$ is usually minute so that, for most purposes, we can write the Stokes drift velocity as

(3.8) \begin{equation} u_L = a^2 \omega k \frac {\cosh [2k (h+z_0)]}{2\sinh ^2(kh)}, \end{equation}

after dividing (3.5) by $\tau$ . Equation (3.8) for the Stokes drift velocity $u_L$ is compared with the background current (2.13) induced by a wave setup or setdown $s$ in figure 4. As can be seen in this figure, the setup-induced velocity dominates the Stokes drift unless the relative wave height $H/h$ is significantly larger than $1$ which is already outside of the linear regime. For particles located near the fluid bed, this discrepancy is even more extreme. It is also noticed that the relative position of a fluid particle along the water column affects the Stokes drift and this influence is more pronounced in the case of shorter waves, which is an expected result.

Figure 4. Comparison of magnitude of the current $V/c_0$ (dashed line) due to a setup or setdown $s/h$ with the Stokes drift velocity $u_L/c_0$ (solid line) as a function of the wave height $H/h$ , where $c_0 = \sqrt {g h}$ . The drift is compared at the undisturbed free surface (a) $z/h = 0$ , at the upper third of the fluid column (b) $z/h = -0.3$ , and at the fluid bed (c) $z/h = -1$ . The red and blue lines represent $u_L$ for $kh = 2/5$ and $8/5$ , respectively.

Alternatively, if a particle path via the second-order approximation is integrated over one Lagrangian period $\tau _L$ rather than one wave period $\tau$ , then a more explicit expression for (3.7) can be obtained. Recall that (3.5) for $x_L$ results from integrating over $\tau$ , not $\tau _L$ . For this purpose, we directly consider the second-order correction $u'$ to the horizontal velocity $u = \phi _x$ . Taylor expanding (3.2) for $\phi _x(\xi (t),\zeta (t),t)$ about $(x_0,z_0)$ and combining with the Taylor expansions (3.4) for $(\xi (t),\zeta (t))$ yields

(3.9) \begin{align} u'(t) = \frac {a^2 \omega k}{\sinh ^2(kh)} \big[ \cosh ^2[k (h+z_0)] \sin ^2(k x_0 - \omega t) + \sinh ^2[k (h+z_0)] \cos ^2(k x_0 - \omega t) \big]. \end{align}

Then integrating this expression over one Lagrangian period gives

(3.10) \begin{align} u_L & = \frac {1}{\tau _L} \int _0^{\tau _L} u'(t) {\rm d}t \nonumber \\ & = a^2 \omega k \frac {\cosh [2k (h+z_0)]}{2 \sinh ^2(kh)} - \frac {a^2 k}{4 \tau _L \sinh ^2(kh)} \big[\sin (2 k x_0) - \sin (2 k x_0 - 2 \omega \tau _L) \big], \end{align}

which follows from

(3.11) \begin{align} \int _0^{\tau _L} \sin ^2(k x_0 - \omega t) {\rm d}t & = \frac {1}{2\omega } \int _{k x_0 - \omega \tau _L}^{k x_0} \big [ 1 - \cos (2t) \big ] {\rm d}t, \nonumber \\ \int _0^{\tau _L} \cos ^2(k x_0 - \omega t) {\rm d}t & = \frac {1}{2\omega } \int _{k x_0 - \omega \tau _L}^{k x_0} \big [ 1 + \cos (2t) \big ] {\rm d}t. \end{align}

The first term in (3.10) coincides with (3.8) while the remainder measures the deviation caused by a phase shift because $\tau _L \neq 2\pi /\omega = \tau$ . This deviation vanishes exactly if $\tau _L = \tau$ since $\sin (2 k x_0)$ and $\sin (2 k x_0 - 2 \omega \tau _L)$ cancel out. Hence, if $\tau _L \simeq \tau$ , the distinction between (3.8) and (3.10) will be insignificant. From numerical simulations of (3.3), we have indeed not observed any discernible discrepancies between (3.8) and (3.10) in this linear regime. Consequently, we will choose (3.8) as the reference analytical estimate for $u_L$ in our subsequent tests.

Note that (3.7) is a general definition of the Stokes drift velocity which also applies to fluid particles under nonlinear waves, in which case $x_L$ and $\tau _L$ would be calculated numerically.

4. Particle paths and linear wave setup

If a wave setup or setdown $s$ is present, then the components of the linear fluid velocity change to

(4.1) \begin{equation} \begin{aligned} \phi _x(x,z,t) & = \omega a \frac {\cosh [k(h+z)]}{\sinh (kh)} \cos (kx-\omega t) + \frac {s g }{ c_0}, \\ \phi _z(x,z,t) & = \omega a \frac {\sinh [k(h+z)]}{\sinh (kh)} \sin (kx-\omega t). \end{aligned} \end{equation}

Integrating over one wave period, the Stokes drift velocity becomes

(4.2) \begin{equation} u_L = a^2 \omega k \frac {\cosh [2k (h+z_0)]}{2\sinh ^2(kh)} + \frac {s g }{c_0}. \end{equation}

This modification of the velocity field can be understood by interpreting the departure of the mean water level from the overall average as an infinitely long wave of the form

(4.3) \begin{align} s = \lim _{k \rightarrow 0} s \cos {(kx - \omega t)}. \end{align}

The associated velocity field is

(4.4) \begin{align} \lim _{k \rightarrow 0} s \sqrt {gk \tanh (k h)} \, \, \, \frac {\cosh [k(h+z)]}{\sinh (kh)} \cos (kx-\omega t), \end{align}

which reduces to

(4.5) \begin{equation} s \sqrt {\frac {g}{h}} = \frac {s g }{ \sqrt {g h} } = \frac {s g }{ c_0}, \end{equation}

as a function of the linear wave speed $c_0 = \sqrt {gh}$ in the long-wave shallow-water regime. Alternatively, this correction

(4.6) \begin{equation} V \sim \frac {s g}{c_0}, \end{equation}

may be interpreted as the limit of (2.13) for small $s$ and for $c \sim c_0$ (i.e. linear waves in the shallow-water limit). Equation (4.6) also gives a clue about the size of $s$ relative to the depth $h$ as it can be rewritten in the form

(4.7) \begin{align} \frac {V}{c_0} \sim \frac {s}{h}. \end{align}

Hence, the relative sizes of the speed correction $V$ and the setup/setdown $s$ are comparable in the long-wave limit.

In the presence of a non-zero mean level $s$ , (3.3) become

(4.8) \begin{equation} \begin{aligned} \frac {{\rm d}\xi }{{\rm d}t} & = \omega a \frac {\cosh [k (h+\zeta )]}{\sinh (kh)} \cos (k \xi - \omega t) + \frac {s g }{ c_0}, \\ \frac {{\rm d}\zeta }{{\rm d}t} & = \omega a\frac { \sinh [k(h+\zeta )]}{\sinh (kh)} \sin (k \xi -\omega t). \end{aligned} \end{equation}

In general, (3.3) or (4.8) must be integrated numerically via an ODE solver to determine the particle trajectories, even for a linear velocity field because of the nonlinear dependence on $\xi (t)$ and $\zeta (t)$ .

5. The KdV theory for particle paths

Extending the previous linear analysis to a weakly nonlinear regime, we now outline the modelling of particle paths in the KdV theory. For the purpose of this study, we write the KdV equation in the form

(5.1) \begin{equation} \eta _t + c_0 \eta _x + \frac {3 c_0}{2h} \eta \eta _x + \frac {c_0 h^2}{6} \eta _{xxx} = 0. \end{equation}

A close evaluation of the physical derivation of the KdV equation reveals that it follows the same idea as the derivation of standard Boussinesq systems for shallow water waves (Whitham Reference Whitham1974; Ali & Kalisch Reference Ali and Kalisch2014). In particular, the derivation gives rise to an approximate velocity potential $\phi (x,z,t) = f - (1/2)(z+h)^2 f_{xx}$ given in terms of a function $f(x,t)$ that represents the trace of the potential at the flat bottom boundary ( $z = -h$ ). The function $f$ can be related to the horizontal fluid velocity $u^{(b)}$ at the bottom (remember this is an inviscid theory) by $f_x = u^{(b)}$ . More specifically, the free-surface elevation $\eta$ and this velocity component $u^{(b)}$ are related by

(5.2) \begin{align} u^{(b)}(x,t) = \frac {c_0}{h} \eta - \frac {c_0}{4h^2} \eta ^2 + \frac {c_0 h}{3}\eta _{xx}. \end{align}

More generally, given a solution $\eta (x,t)$ of the KdV equation, we can find the horizontal velocity field across the water column in the form

(5.3) \begin{equation} \phi _x(x,z,t) = \frac {c_0}{h} \left [ \eta - \frac {1}{4h} \eta ^2 + \left ( \frac {h^2}{3} - \frac {z^2}{2} \right )\eta _{xx} \right ]. \end{equation}

Similarly, the vertical velocity field can be evaluated by

(5.4) \begin{equation} \phi _z(x,z,t)= -\frac {c_0}{h} (z + h) \eta _x. \end{equation}

As shown in Borluk & Kalisch (Reference Borluk and Kalisch2012), the velocity field (5.3)–(5.4) can be integrated to describe the Lagrangian motion of fluid particles

(5.5) \begin{equation} \frac {{\rm d}\xi }{{\rm d}t} = \phi _x(\xi (t),\zeta (t),t), \quad \frac {{\rm d}\zeta }{{\rm d}t} = \phi _z(\xi (t),\zeta (t),t). \end{equation}

These equations define a particle trajectory $(\xi (t),\zeta (t))$ emanating from any initial point $(\xi (0),\zeta (0))$ in the fluid domain, provided that the velocity field is known. If a solution of the KdV equation has been obtained, either numerically or in exact form such as a solitary or cnoidal wave, then the particle trajectories can be approximated numerically. In this study, we will use a fourth-order Runge–Kutta scheme to solve (5.5) from the KdV solution. As suggested by Borluk & Kalisch (Reference Borluk and Kalisch2012), it can be beneficial to include a higher-order correction term in these equations, although this is not strictly required. Constantin (Reference Constantin2007) proposed an alternate method for constructing very accurate velocity fields associated with solitary wave solutions of the KdV equation but, in the present work, we follow the approach of Borluk & Kalisch (Reference Borluk and Kalisch2012). As seen in Bjørnestad et al. (Reference Bjørnestad, Buckley, Kalisch, Streßer, Horstmann, Frøysa, Ige, Cysewski and Carrasco-Alvarez2021), the KdV particle paths are able to capture particle trajectories measured in field experiments.

The cnoidal solution of the KdV equation is given by

(5.6) \begin{equation} \eta (x,t) = f_2 - (f_2 - f_1) \, \textrm{cn}^2 \left [ \frac {2K(m)}{\lambda } ( x - c \, t ); m \right ], \end{equation}

where the function ${cn}(F;m)$ is the Jacobi elliptic cosine function defined in terms of the inverse of the incomplete elliptic integral of first kind (Lawden Reference Lawden2013). In other words, given a number $0 \lt m \lt 1$ , if we define

(5.7) \begin{align} F(\gamma ;m) = \int _0^{\gamma } \frac {{\rm d} \theta }{\sqrt {1 - m \sin ^2 \theta }}, \end{align}

and let $\text{am}(F;m) = \gamma$ be the Jacobi amplitude, then ${cn}(F;m) = \cos \text{am}(F;m)$ . The number $m$ is the elliptic modulus defined in terms of the constants $f_1$ , $f_2$ and $f_3$ by $m=(f_{1}-f_{2})/(f_{1}-f_{3})$ , and $K(m)$ is the complete elliptic integral of first kind. The wavelength $\lambda$ and wave speed $c$ are given by

(5.8) \begin{equation} \lambda = K(m) \sqrt {\frac {16 h^3}{3(f_1 - f_3)}}, \end{equation}

and

(5.9) \begin{equation} c = c_0 \Big {(} 1 + \frac {f_1 + f_2 + f_3}{2h}\Big {)}, \end{equation}

respectively.

When using these formulae in practice, it is convenient to take the wave height $H$ , the mean surface level $s$ and the elliptic modulus $m$ as given parameters. The constants $f_1$ , $f_2$ and $f_3$ can then be computed explicitly as follows:

(5.10) \begin{equation} f_3 = s - \frac {H E(m)}{mK(m)}, \quad f_1 = f_3 + \frac {H}{m}, \quad f_2 = f_1 - H, \end{equation}

where $E(m)$ is the complete elliptic integral of second kind (Lawden Reference Lawden2013).

Since the function ${cn}$ oscillates between $0$ and $1$ , it is apparent that $f_1$ denotes the $z$ -value of the wave crest, and $f_2$ represents the wave trough. In the non-dimensionalized situation where the undisturbed depth $h$ is taken as a unit length, a cnoidal wave is completely determined if the wavelength $\lambda$ , the wave height $H$ and the mean surface level $s$ are given. This implies in particular that a cnoidal wave can be readily computed in the presence of a setup or setdown, with $s \gtrless 0$ in (5.10) while the formulae (5.3)–(5.4) and (5.6) for the velocity field and surface elevation remain unchanged. The elliptic modulus $m$ must be evaluated numerically by minimizing the difference between (5.8) and a prescribed wavelength $\lambda$ (with fixed $H$ and $s$ ). The expression (5.6) is then substituted into (5.3) and (5.4) which are needed in the ODEs (5.5).

Note that (5.5) is a general definition of the dynamical equations governing the Lagrangian motion of any fluid particle, which will also be used for the full nonlinear problem. However, a difficulty in this case is the calculation of the fluid velocity field $\boldsymbol{u} = \boldsymbol{\nabla }\phi$ for which a numerical procedure will be devised as reported in § 9.

6. Hamiltonian formulation for nonlinear waves

After examining the linear and KdV regimes, we now come back to the full nonlinear problem for a more detailed investigation of flow kinematics. We first present a mathematical reformulation of the water wave equations (2.1)–(2.4) which forms the basis for our nonlinear potential-flow solver, before describing our numerical approach to reconstruct the velocity field and to simulate particle motions. An adjustment will be proposed to accommodate the effects of a wave setup or setdown in this nonlinear setting.

Zakharov (Reference Zakharov1968) and Craig & Sulem (Reference Craig and Sulem1993) showed that the Laplace problem (2.1)–(2.4) is equivalent to a Hamiltonian system in canonical form

(6.1) \begin{equation} \begin{pmatrix} \eta _t \\ \psi _t \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} \mathbb{H}_\eta \\ \mathbb{H}_\psi \end{pmatrix}, \end{equation}

where the two conjugate variables are the surface elevation $\eta (x,t)$ and

(6.2) \begin{equation} \psi (x,t) = \phi (x,\eta (x,t),t), \end{equation}

which is the trace of the velocity potential evaluated on $z = \eta (x,t)$ . The corresponding Hamiltonian is given by

(6.3) \begin{align} \mathbb{H} = \frac {1}{2} \int _{-\infty }^\infty \big [ \psi G(\eta ) \psi + g \eta ^2 \big ] {\rm d}x, \end{align}

in terms of the DNO

(6.4) \begin{equation} G(\eta ) : \psi \mapsto \sqrt {1 + \eta _x^2} \, \phi _n \big |_{z=\eta } , \end{equation}

which by definition takes Dirichlet data $\psi$ on the free surface, solves Laplace’s equation (2.1) subject to (2.4) and returns the associated Neumann data.

More explicitly, (6.1) read

(6.5) \begin{align} \eta _t & = G(\eta ) \psi , \end{align}

(6.6) \begin{align} \psi _t & = -g \eta + \frac {1}{2} \psi _x^2 + \frac {1}{2} \frac {\big [ G(\eta ) \psi + \eta _x \psi _x \big ]^2}{1 + \eta _x^2}, \end{align}

which constitute a closed system of evolution equations in terms of surface variables alone. Note that the function $\eta$ is assumed to be a graph of $x$ , so overturning waves with a multivalued profile are not permitted here. More details on this Hamiltonian formulation of (2.1)–(2.4) can be found in Craig et al. (Reference Craig, Guyenne and Sulem2021). Because the DNO is analytic for sufficiently smooth $\eta$ , it can be expressed via a convergent Taylor series expansion

(6.7) \begin{equation} G(\eta ) = \sum _{j=0}^\infty G_j(\eta ), \end{equation}

about the quiescent state $\eta = 0$ , where each term $G_j$ is homogeneous of degree $j$ in $\eta$ and can be determined by a recursion formula (Craig & Sulem Reference Craig and Sulem1993; Xu & Guyenne Reference Xu and Guyenne2009; Li et al. Reference Li, Xu, Guyenne and Yu2010).

To introduce our reconstruction method for the internal velocity field, we find it suitable to recall the derivation of (6.5)–(6.6) as these are linked to the required surface data. First, by differentiating (6.2) and using the chain rule, we note that

(6.8) \begin{equation} \psi _t = \phi _t + \eta _t \phi _z \big |_{z=\eta } , \quad \psi _x = \phi _x + \eta _x \phi _z \big |_{z=\eta }. \end{equation}

The definition (6.4) of the DNO implies

(6.9) \begin{equation} G(\eta ) \psi = \phi _z - \eta _x \phi _x \big |_{z=\eta }, \end{equation}

which yields

(6.10) \begin{align} \phi _z \big |_{z=\eta } = G(\eta ) \psi + \eta _x \phi _x \big |_{z=\eta } = G(\eta ) \psi + \eta _x (\psi _x - \eta _x \phi _z) \big |_{z=\eta }, \end{align}

hence

(6.11) \begin{equation} v^{(s)}(x,t) = \phi _z \big |_{z=\eta } = \frac {G(\eta ) \psi + \eta _x \psi _x}{1 + \eta _x^2}, \end{equation}

after inserting the second equation from (6.8). Then substituting (6.11) back into (6.8) leads to

(6.12) \begin{equation} u^{(s)}(x,t) = \phi _x \big |_{z=\eta } = \psi _x - \frac {\eta _x}{1 + \eta _x^2} \big ( G(\eta ) \psi + \eta _x \psi _x \big ) = \frac {\psi _x - \eta _x G(\eta ) \psi }{1 + \eta _x^2}. \end{equation}

Adding up the squares of (6.11) and (6.12) results in

(6.13) \begin{equation} \phi _x^2 + \phi _z^2 \big |_{z=\eta } = \frac {\psi _x^2 + \big ( G(\eta ) \psi \big )^2}{1 + \eta _x^2}. \end{equation}

It is easy to see that the kinematic boundary condition (2.2) coincides with (6.5) by virtue of (6.9).

On the other hand, substituting (6.5) and (6.11) into the first equation of (6.8) gives

(6.14) \begin{align} \phi _t \big |_{z=\eta } = \psi _t - \frac {G(\eta ) \psi }{1 + \eta _x^2} \big ( G(\eta ) \psi + \eta _x \psi _x \big ). \end{align}

Then combining this expression with (6.13) shows the equivalence between the dynamic boundary condition (2.3) and the Hamiltonian equation (6.6).

7. Nonlinear surface kinematics

For potential flow, (6.11) and (6.12) yield, respectively, the vertical and horizontal components of the fluid velocity at the free surface (denoted by $v^{(s)}$ and $u^{(s)}$ ), and thus can be used to compute these quantities explicitly given the DNO and the solution $(\eta ,\psi )^\top$ of (6.5)–(6.6) at any instant $t$ .

Furthermore, given the surface data $(\eta ,\psi )^\top$ , the Lagrangian motion of any fluid particle labelled $\ell$ and located at $(x,z) = (\xi _\ell (t),\zeta _\ell (t))$ on the free surface is governed by

(7.1) \begin{equation} \frac {{\rm d} \xi _\ell }{{\rm d}t} = u^{(s)}(\xi _\ell ,t) = \frac {\psi _x - \eta _x G(\eta ) \psi }{1 + \eta _x^2} \Big |_{x = \xi _\ell }, \end{equation}

according to (6.12), where $\zeta _\ell (t) = \eta (\xi _\ell (t),t)$ . An evolution equation to determine the vertical position $\zeta _\ell$ is not required because the fluid particle is constrained to lie on the free surface.

Equation (7.1) may be viewed as an alternate version of the JS equations to describe fluid particles on a free surface (Fedele et al. Reference Fedele, Chandre and Farazmand2016). In the two-dimensional setting, the JS system reduces to a single ODE for $\xi _\ell$ , which is second-order in time and depends explicitly on the surface elevation $\eta$ along with its derivatives. By contrast, (7.1) is a first-order ODE that depends on surface variables $\eta$ and $\psi$ . This dual dependence is not an issue since the pair $(\eta ,\psi )^\top$ is typically given by the solution of the water-wave equations. We point out again that (7.1) is deduced from an exact nonlinear representation of the horizontal fluid velocity as derived in the previous section. Moreover, the dual dependence of (7.1) on $\eta$ and $\psi$ is explicit in the present formulation, via the DNO as defined by (6.7).

Equation (7.1) may thus have advantages when it comes to simulating the dynamics of surface particles, in combination with the numerical solution of (6.5)–(6.6) for surface waves. However, this equation will not be solved here because, as detailed in a subsequent section, we will propose a more general (albeit less explicit) approach to evaluate the velocity field on or under the free surface, which enables the description of Lagrangian particles anywhere in the fluid domain. Instead of considering (7.1), this general approach incorporates $u^{(s)}$ as a constraint in the calculation of the internal velocity field to match it with surface values.

8. Numerical methods

Assuming periodic boundary conditions in $x$ , we use a pseudospectral method to discretize in space the equations of motion (6.5)–(6.6) together with the DNO. Both functions $\eta$ and $\psi$ are represented by truncated Fourier series

(8.1) \begin{equation} \begin{pmatrix} \eta _j \\[3pt] \psi _j \end{pmatrix} = \sum _{n=-N/2+1}^{N/2} \begin{pmatrix} \widehat \eta _n \\[3pt] \widehat \psi _n \end{pmatrix} e^{i k_n x_j}, \quad k_n = \frac {2\pi }{L} n, \quad x_j = \frac {L}{N} j, \quad j = \{ 0, \ldots , N-1 \}, \end{equation}

where $\eta _j \simeq \eta (x_j)$ and $\psi _j \simeq \psi (x_j)$ in the periodic domain $0 \leqslant x_j \lt L$ . The time dependence is omitted here for convenience. Spatial derivatives and Fourier multipliers are evaluated in the Fourier space, while the nonlinear products are calculated in the physical space on a regular grid with $N$ collocation points. All operations going back and forth between the Fourier and physical spaces are performed via the fast Fourier transform.

The DNO is approximated by a truncated version of (6.7), namely

(8.2) \begin{align} G(\eta ) \simeq \sum _{j=0}^M G_j(\eta ), \end{align}

for which a small number $M$ of terms is sufficient to achieve highly accurate results. The value $M = 6$ is selected based on previous extensive tests (Craig et al. Reference Craig, Guyenne, Hammack, Henderson and Sulem2006; Xu & Guyenne Reference Xu and Guyenne2009). Time integration of (6.5) and (6.6) is implemented in the Fourier space so that the linear terms can be solved exactly by the integrating factor technique. The associated nonlinear system is marched in time by using the fourth-order Runge–Kutta scheme with constant step $\Delta t$ . More details on this numerical model can be found in Craig & Sulem (Reference Craig and Sulem1993) and Guyenne & Nicholls (Reference Guyenne and Nicholls2007). We also point out that artificial filtering was not employed in any of the cases that we have examined here. Filtering was not applied to either simulations of the surface waves ((6.5)–(6.6)) or simulations of the fluid particles (see (11.1) in § 11).

9. Calculation of the internal velocity field

Recall that (6.5)–(6.6) form a closed system of evolution equations that leads to an efficient and accurate numerical solver for the nonlinear surface wave problem. The question now is how to infer the internal flow kinematics (i.e. the internal fluid velocity) from surface data, which in turn governs the dynamics of Lagrangian particles in the fluid domain.

The starting point is the general solution

(9.1) \begin{equation} \phi (x,z,t) = \sum _{n=-N/2+1}^{N/2} q_n(t) \frac {\cosh \big [ k_n (z + h) \big ]}{\cosh (k_n h)} e^{i k_n x}, \quad 0 \leqslant x \lt L, \quad -h \lt z \lt \eta (x,t), \end{equation}

to Laplace’s equation (2.1) subject to (2.4) with periodic boundary conditions in $x$ , consistent with the pseudospectral method being employed to solve (6.5)–(6.6) for $\eta$ and $\psi$ . Equation (9.1) generalizes (3.1) to the prospect of nonlinear multimodal waves. From (9.1), the components of the fluid velocity are given by

(9.2) \begin{align} u(x,z,t) = \phi _x(x,z,t) & = \sum _{n=-N/2+1}^{N/2} i k_n q_n(t) \frac {\cosh \big [ k_n (z + h) \big ]}{\cosh (k_n h)} e^{i k_n x}, \end{align}

(9.3) \begin{align} v(x,z,t) = \phi _z(x,z,t) & = \sum _{n=-N/2+1}^{N/2} k_n q_n(t) \frac {\sinh \big [ k_n (z + h) \big ]}{\cosh (k_n h)} e^{i k_n x}, \end{align}

at any point $(x,z)$ on or under the free surface. The complex-valued coefficients $\{ q_n(t) \}$ are to be found whenever data on $(u,v)^\top$ are needed. Note that they are only functions of $k_n$ and $t$ . For this purpose, at every such instant $t$ , we match

(9.4) \begin{equation} u^{(s)}(x_j,t) = u(x_j,\eta (x_j,t),t), \end{equation}

at all nodes $\{ x_j \}$ on the pseudospectral grid. The left- and right-hand sides of (9.4) are given by (6.12) and (9.2), respectively. This produces an algebraic linear system

(9.5) \begin{equation} A X = B, \quad B = \{ u^{(s)}(x_j,t) \}, \end{equation}

to be solved for $X = \{ k_n q_n(t) \}$ where

(9.6) \begin{equation} A = \{ A_{jn} \} = i \big [ \cosh (k_n \eta (x_j,t)) + \tanh (k_n h) \sinh (k_n \eta (x_j,t)) \big ] e^{i k_n x_j}. \end{equation}

The coefficient matrix $A$ is a square $N \times N$ matrix, consistent with the pseudospectral discretization (8.1) for the surface wave dynamics. In practice, a lower resolution $\lt N$ may be prescribed when assembling and solving the linear system (9.5), to allow for more efficient computations of the internal velocity field. It is also preferable to absorb $k_n$ into the vector $X$ of unknowns in order to avoid multiplying the entries of the matrix $A$ by a factor that may amplify numerical errors at high wavenumbers during the solution process of (9.5).

Some general comments are in order.

1. (i) The coefficients $\{ q_n \}$ are in general not directly related to the Fourier coefficients of $u^{(s)}$ (except in the linear case where quantities are evaluated at $z = 0$ ) because the condition (9.4) requires that (9.2) be evaluated at $z = \eta (x_j,t)$ through the $\cosh$ function, hence the necessity to solve a non-local problem via the linear system (9.5).

2. (ii) This linear system involves a dense coefficient matrix and must be solved repeatedly during the time evolution. In the present two-dimensional case, we use a direct Gaussian elimination solver which still leads to simulations at a reasonable cost. It is expected that a more efficient solver (e.g. an iterative scheme like GMRES (generalized minimal residual method)) would be needed in the three-dimensional setting to deal with larger systems (9.5).

3. (iii) At any instant $t$ , after the coefficients $\{ q_n(t) \}$ have been determined, the velocity components $(u,v)$ can be reconstructed at any location $(x,z)$ in the fluid domain via (9.2)–(9.3). We point out that this direct computation is not perturbative and thus is not sensitive to the smoothness of the free surface, in contrast to the perturbative approach proposed by Bateman et al. (Reference Bateman, Swan and Taylor2003).

4. (iv) To solve for $\{ q_n \}$ , we choose to impose a matching condition (9.4) on the horizontal velocity $u$ rather than on the velocity potential $\phi$ , because $\phi$ is only known up to an additive constant by definition. Using (9.4) allows us to accommodate effects from a mean flow or background current which may be relevant, especially in shallow water (e.g. nearshore areas) or in a confined environment (e.g. wave-tank experiments). By absorbing $k_n$ into $X$ , a non-trivial zeroth mode may be deduced for the internal velocity field, given the surface data $u^{(s)}$ .

5. (v) The basic idea following (9.4) to compute the coefficients in the harmonic decomposition of the fluid velocity field was actually suggested by Bateman et al. (Reference Bateman, Swan and Taylor2003). However, instead of tackling this calculation in a direct manner, these authors devised an iterative scheme via a Taylor series representation of the flow kinematics about the free surface. Due to convergence issues for locations far beneath the free surface or associated with the surface smoothness, their method requires the sequential application of two iterative schemes (their so-called $H$ and $H_2$ operators).

Figure 5. Profiles of horizontal velocity $u/c_0$ along the water column $z/h$ for nonlinear travelling waves with varying wave steepness $k a$ and varying water depth $k h$ . The reconstruction results (red dots) are compared with steady-state predictions by Fenton’s method (blue line). (a) $\textit{ka} = 0.005, \textit{kh} = 0.3$ , (b) $\textit{ka} = 0.05, \textit{kh} = 1$ , (c) $\textit{ka} = 0.05, \textit{kh} = 2\pi$ , (d) $\textit{ka} = 0.01, \textit{kh} = 0.3$ , (e) $\textit{ka} = 0.1, \textit{kh} = 1$ , (f) $\textit{ka} = 0.1, \textit{kh} = 2\pi$ , (g) $\textit{ka} = 0.02, \textit{kh} = 0.3$ , (h) $\textit{ka} = 0.15, \textit{kh} = 1$ , (i) $\textit{ka} = 0.2, \textit{kh} = 2\pi$ , (j) $\textit{ka} = 0.03, \textit{kh} = 0.3$ , (k) $\textit{ka} = 0.2, \textit{kh} = 1$ , (l) $\textit{ka} = 0.3, \textit{kh} = 2\pi$ .

Figure 6. Reciprocal condition number $\kappa$ of the coefficient matrix (9.6) with Fenton’s solution for water depth $k h = 0.3$ (blue) and $k h = 1$ (red) as a function of (a) wave steepness $k a$ and (b) number $N$ of Fourier modes.

10. Numerical tests

We assess the performance of this reconstruction algorithm by testing it against Fenton’s method for periodic travelling waves (Fenton Reference Fenton1988). The latter method solves a steady version of (2.1)–(2.4) in a reference frame moving at constant speed $c$ for solutions of the form (8.1)–(9.3) with $x$ replaced by $x - c \, t$ . The coefficients in these expansions (including the wave speed $c$ ) are computed numerically by solving the corresponding nonlinear system iteratively up to an arbitrary order of truncation (see also Vanden-Broeck (Reference Vanden-Broeck2010) and references therein).

In the present tests, we specify Fenton’s solution as an initial condition for our numerical solver of (6.5)–(6.6). All equations are non-dimensionalized according to shallow-water theory, i.e. lengths are divided by $h_0$ and times are divided by $\sqrt {h_0/g}$ (where $h_0$ is a characteristic water depth) so that $g = 1$ in dimensionless units. We simulate the travelling wave over some time interval (up to $t = T$ ) at which point we reconstruct the internal velocity field using the simulated surface data and compare our results with the steady-state predictions by Fenton’s method.

Figure 5 shows profiles of the horizontal velocity $u$ along the water column for varying wave steepness $k a$ and varying water depth $k h$ . The values of $u$ are normalized by the linear phase speed $c_0 = \sqrt {g \tanh (kh)/k}$ , while the values of $z$ are normalized by $h$ . In all these computations, the domain length and the carrier wavenumber are set to $L = 2\pi$ and $k = 1$ , respectively. Varying the wave steepness and water depth affects the properties of wave nonlinearity and dispersion, with $k h = 0.3$ , $1$ , $2\pi$ corresponding to a regime of shallow, intermediate and deep water, respectively. In each case, the internal points are distributed vertically from near the bottom all the way up to the free surface, at an abscissa $x_0$ near the central wave crest. This location $x_0$ need not be selected among the grid points $\{ x_j \}$ . We see an excellent agreement between these two predictions in all the graphs of figure 5. The numerical solution of (6.5)–(6.6) is computed up to $t = T = 10$ with resolutions $\Delta t = 0.01$ and $\Delta x = L/N = 0.05$ ( $N = 128$ ). As expected, the vertical structure of the flow is more uniform in shallower water, while it exhibits more variation in deeper water, decreasing quickly in magnitude with the depth. For the highest point in each plot (which is taken at the free surface), we have also estimated the horizontal velocity by interpolating the simulated surface values (6.12) at $x = x_0$ (using linear interpolation) and found that it coincides with the reconstructed value at $(x,z) = (x_0,\eta (x_0,t))$ based on (9.2). The surface elevation $\eta (x_0,t)$ is also estimated by linear interpolation of the simulated values (8.1). This confirms that the matching condition (9.4) is properly enforced by our reconstruction method, and it is reassuring to see in figure 5 that the velocity profiles vary smoothly along the entire water column.

It is well known that such systems as defined by (9.5)–(9.6) suffer from ill-conditioning, which may lead to inaccurate solutions (Wilkening & Vasan Reference Wilkening and Vasan2015). We investigate this issue by computing the reciprocal condition number $\kappa$ of the coefficient matrix $A$ , where $\eta$ is prescribed by Fenton’s solution. Typically, if $A$ is well conditioned, $\kappa$ is near $1$ . If $A$ is badly conditioned, $\kappa$ is near $\varepsilon _m$ where $\varepsilon _m \sim 10^{-16}$ represents machine epsilon in double precision (Golub & Van Loan Reference Golub and Van Loan1996). Estimates of $\kappa$ are reported in figure 6 for various values of water depth $k h$ and wave steepness $k a$ as illustrated in figure 5. We also examine the dependence of $\kappa$ on the number $N$ of Fourier modes (which determines the size of $A$ ). These tests confirm that the conditioning of $A$ deteriorates in cases of deeper water, steeper waves or for a larger matrix (i.e. a finer spatial resolution). As will be discussed in § 12, the present study considers the problem of wave setup and setdown with nonlinear waves of moderate steepness ( $k a \leqslant 0.16$ ) on shallow or intermediate water depth ( $k h \leqslant 1.6$ ), for which ill-conditioning is not severe in accordance with figure 6.

11. Particle paths for nonlinear waves

After deducing the fluid velocity, we can address the simulation of particle dynamics. Similar to (5.5), we compute the Lagrangian trajectory of any fluid particle labelled $\ell$ on or under the free surface by solving the equations of motion

(11.1) \begin{equation} \frac {{\rm d}\xi _\ell }{{\rm d}t} = u(\xi _\ell ,\zeta _\ell ,t), \quad \frac {{\rm d}\zeta _\ell }{{\rm d}t} = v(\xi _\ell ,\zeta _\ell ,t), \end{equation}

where $(x,z) = (\xi _\ell (t),\zeta _\ell (t))$ denote the particle coordinates at time $t$ along its trajectory. Equations (11.1) mean that a fluid particle is advected by the flow according to the velocity field (9.2)–(9.3).

Note that (11.1) are evaluated at any time $t$ after solving (9.5) for $\{ q_n(t) \}$ given $(\eta ,\psi )$ . Therefore, we integrate (11.1) by employing the fourth-order Adams–Bashforth scheme

(11.2) \begin{equation} \boldsymbol{x}_\ell ^{(m+1)} = \boldsymbol{x}_\ell ^{(m)} + \frac {\Delta t}{24} \big[ 55 \boldsymbol{u}_\ell ^{(m)} - 59 \boldsymbol{u}_\ell ^{(m-1)} + 37 \boldsymbol{u}_\ell ^{(m-2)} - 9 \boldsymbol{u}_\ell ^{(m-3)} \big], \quad m \geqslant 3, \end{equation}

to update the particle position from time $t_m$ to time $t_{m+1} = t_m + \Delta t$ , where $\boldsymbol{x}_\ell ^{(m)} = (\xi _\ell ^{(m)},\zeta _\ell ^{(m)})$ , $\boldsymbol{u}_\ell ^{(m)} = ( u(\xi _\ell ^{(m)},\zeta _\ell ^{(m)},t_m),v(\xi _\ell ^{(m)},\zeta _\ell ^{(m)},t_m) )$ with $\xi _\ell ^{(m)} \simeq \xi _\ell (t_m)$ , $\zeta _\ell ^{(m)} \simeq \zeta _\ell (t_m)$ .

Advantages of (11.2) include its explicit character and fourth-order accuracy, consistent with the time integrator for the water-wave system (6.5)–(6.6) and for the particle equations (5.5) in the KdV approximation. However, unlike a Runge–Kutta scheme, (11.2) does not involve intermediate steps between $t_m$ and $t_{m+1}$ , which would require additional calculations of (9.5) and thus would entail a higher computational cost. For the first three time steps, the (first-order) forward Euler scheme is adopted to compute $\{ \boldsymbol{x}_\ell ^{(1)}, \ldots , \boldsymbol{x}_\ell ^{(3)} \}$ from the initial position $\boldsymbol{x}_\ell ^{(0)}$ at $t = t_0$ . Because the water-wave system (6.5)–(6.6) and associated particle equations (11.1) are solved simultaneously, the same time step $\Delta t$ is specified in their respective numerical schemes.

For a linear travelling wave of amplitude $a$ , wavenumber $k$ and frequency $\omega$ , (11.1) describing a particle trajectory labelled $\ell$ reduce to (3.3) and can be rewritten as

(11.3) \begin{align} \frac {{\rm d} \widetilde \xi _\ell }{{\rm d}t} & = \omega a \frac {\cosh [k (h+\zeta _\ell )]}{\sinh (kh)} \cos (k \widetilde \xi _\ell ) - c, \end{align}

(11.4) \begin{align} \frac {{\rm d}\zeta _\ell }{{\rm d}t} & = \omega a\frac { \sinh [k(h+\zeta _\ell )]}{\sinh (kh)} \sin (k \widetilde \xi _\ell ), \end{align}

for $\widetilde \xi _\ell = \xi _\ell - c \, t$ where $c = \omega /k$ is the phase speed. In this form, (11.3)–(11.4) possess a canonical Hamiltonian structure

(11.5) \begin{equation} \frac {{\rm d} \widetilde \xi _\ell }{{\rm d}t} = \frac {\partial {\mathcal H}_\ell }{\partial \zeta _\ell }, \quad \frac {{\rm d} \zeta _\ell }{{\rm d}t} = -\frac {\partial {\mathcal H}_\ell }{\partial \widetilde \xi _\ell }, \end{equation}

where the Hamiltonian

(11.6) \begin{equation} {\mathcal H}_\ell = c \, a \frac {\sinh [k (h+\zeta _\ell )]}{\sinh (kh)} \cos (k \widetilde \xi _\ell ) - c \, \zeta _\ell , \end{equation}

is conserved over time. The change of variables $\xi _\ell \to \widetilde \xi _\ell$ helps accommodate the explicit dependence on $t$ so that

(11.7) \begin{align} \frac {{\rm d}{\mathcal H}_\ell }{{\rm d}t} = \frac {\partial {\mathcal H}_\ell }{\partial t} + \frac {{\rm d}\widetilde \xi _\ell }{{\rm d}t} \frac {\partial {\mathcal H}_\ell }{\partial \widetilde \xi _\ell } + \frac {{\rm d}\zeta _\ell }{{\rm d}t} \frac {\partial {\mathcal H}_\ell }{\partial \zeta _\ell } = 0. \end{align}

This Hamiltonian structure readily extends to nonlinear travelling waves by exploiting their special functional dependence on $x - c \, t$ . With the solution form (9.1) in mind, this property implies that $q_n(t) = Q_n e^{-i k_n c \, t}$ where $Q_n$ are constants (Fenton Reference Fenton1988; Chang et al. Reference Chang, Chen and Liou2009). As a consequence, the Hamiltonian formulation (11.5) of (11.1) still holds in this nonlinear case, modulo the following series expansion for the Hamiltonian:

(11.8) \begin{equation} {\mathcal H}_\ell = \sum _{n=-N/2+1}^{N/2} i Q_n \frac {\sinh \big [ k_n (h + \zeta _\ell ) \big ]}{\cosh (k_n h)} e^{i k_n \widetilde \xi _\ell } - c \, \zeta _\ell . \end{equation}

By construction, this Hamiltonian is also a conserved quantity but now $c$ is meant to be the actual speed of the nonlinear travelling wave (rather than the linear phase speed $c = \omega /k$ for a monochromatic wave).

Figure 7. Particle trajectories near the bottom (c,d) and at the free surface (a,b), for a linear travelling wave of height $H = 2 a = 0.02$ (a,c) and $0.04$ (b,d) with wavenumber $k = 2/5$ in the domain $0 \leqslant x \lt 10 \pi$ , $h = 1$ . The blue line represents the numerical solution with initial (dot) and final (triangle) locations at $t = 0$ and $40$ , respectively. The red line represents the closed orbit as predicted by the first-order approximation. (a) $H = 0.02$ , at the free surface, (b) $H = 0.04$ , at the free surface, (c) $H = 0.02$ , near the bottom, (d) $H = 0.04$ , near the bottom.

However, a Hamiltonian formulation is unclear in the more general (time-evolving) nonlinear setting due to the temporal dependence of the coefficients $\{ q_n(t) \}$ via the matching condition (9.4). Supposedly, such a Hamiltonian formulation (if it exists) might also involve $q_n$ among the conjugate variables in addition to $\xi$ and $\zeta$ .

We test the performance of the numerical integrator (11.2) for particle trajectories by applying it to the case of a linear travelling wave where the equations of motion are given explicitly by (3.3). Note that, even for a linear velocity field, these equations are nonlinear in the particle path parameterization and thus must be solved numerically in general. Simulations of particle trajectories up to $t = T = 40$ are illustrated in figure 7 for wave heights $H = 2 a = 0.01$ and $0.04$ with wavenumber $k = 2/5$ (two wavelengths over the interval $0 \leqslant x \lt L$ ). The domain dimensions are chosen to be $h = 1$ and $L = 10\pi$ with resolution $\Delta x = 0.12$ ( $N = 256$ ). Low wave steepness is selected in order to comply with the linear wave regime being considered here. Accordingly, a small time step $\Delta t = 0.0001$ is used to accurately resolve variations in the short particle paths. As a reference, the closed orbit predicted by the first-order analytical approximation (3.4) is also depicted in this figure, where the centre position $(x_0,z_0)$ is taken to be the midpoint between the minimum and maximum $(x,z)$ -coordinates along the numerically calculated path, namely

(11.9) \begin{align} x_0 = \frac {\min \big ( \xi _\ell (t) \big ) + \max \big ( \xi _\ell (t) \big )}{2}, \quad z_0 = \frac {\min \big ( \zeta _\ell (t) \big ) + \max \big ( \zeta _\ell (t) \big )}{2}, \quad 0 \leqslant t \leqslant T, \end{align}

for any particle $\ell$ . Overall, these two solutions compare well together. However, we can discern that the numerical orbit is not closed exhibiting a tiny drift, which is consistent with the mathematical proof by Constantin (Reference Constantin2006).

Figure 8. Stokes drift velocity $u_L$ for $16$ Lagrangian particles along the water column $z$ under a linear travelling wave of height (a) $H = 0.02$ and (b) $H = 0.04$ with wavenumber $k = 2/5$ . Numerical predictions based on the exact (blue line) and reconstructed (red circles) velocity field are compared with an analytical second-order approximation (black dashed line). (a) $H = 0.02$ , (b) $H = 0.04$ .

Figure 8 portrays the Stokes drift velocity $u_L$ of 16 fluid particles distributed along the water column (from near the bottom up to the free surface) under such a linear wave. Initially (at $t = 0$ ), these 16 particles are evenly distributed along the vertical axis $x = L/2$ under the central wave crest. However, over time, they move individually according to the internal velocity field. As explained in § 3, for each particle trajectory, the Stokes drift velocity (3.7) is computed by evaluating the Stokes transport $x_L$ covered during a Lagrangian period $\tau _L$ , which corresponds to the horizontal distance and time interval between two consecutive maxima of the particle elevation, and then by dividing $x_L$ by $\tau _L$ . Both quantities $x_L$ and $\tau _L$ are determined numerically from the particle simulations. Here the Lagrangian period $\tau _L$ essentially coincides with the linear wave period $\tau = 2\pi /\omega$ . These numerical predictions using a linear velocity field (3.2) are compared with the analytical result based on a second-order approximation (3.8) for the Stokes drift, where the central elevation $z_0$ of each particle path is chosen as indicated above. An excellent agreement is again found considering the very low magnitude of this Stokes drift velocity. The smaller the wave amplitude, the smaller the drift, and the smaller the difference between the analytical and numerical predictions. A perfect match is not expected because the analytical prediction is an asymptotic approximation.

For further assessment, numerical estimates of $u_L$ via reconstruction of the internal velocity field from the numerical solution of (6.5)–(6.6) are also included in figure 8. To accommodate this linear wave regime, the nonlinear terms are omitted from (6.5)–(6.6) with the DNO reducing to $G(\eta ) = G_0 = D \tanh (h D)$ for finite water depth $h$ , where $D = -i \, \partial _x$ . The initial conditions for (6.5)–(6.6) are given by

(11.10) \begin{align} \eta (x,0) = a \cos (kx), \quad \psi (x,0) = \frac {a g}{\omega } \sin (kx). \end{align}

The same time integrator (11.2) with the same time step $\Delta t$ is employed to compute particle trajectories from the reconstructed velocity field. On one hand, the fact that the two numerical predictions (based on the exact form (3.2) and reconstructed form (9.2)–(9.3) of $(u,v)$ ) look indistinguishable in this figure further validates our reconstruction method. On the other hand, the close match between these numerical estimates and the analytical approximation of $u_L$ attests to the high accuracy of the time integrator (11.2) that describes the particle dynamics. We have examined cases with a lower wave amplitude and observed an even closer agreement. Because discrepancies become less noticeable, these results are not shown for convenience. Additionally, we have also checked the conservation of the total Hamiltonian

(11.11) \begin{align} {\mathcal H} = \sum _{\ell =1}^{16} {\mathcal H}_\ell , \end{align}

for all 16 particles in the two previous cases $H = 0.02$ and $0.04$ , where ${\mathcal H}_\ell$ represents each individual contribution. Figure 9 plots the relative error

(11.12) \begin{align} \frac {\Delta {\mathcal H}}{{\mathcal H}_0} = \frac {|{\mathcal H} - {\mathcal H}_0|}{{\mathcal H}_0}, \end{align}

on $\mathcal H$ as a function of $t$ , where ${\mathcal H}_0$ denotes the initial value of $\mathcal H$ at $t = 0$ . Again, numerical estimates based on the exact velocity field (3.2), with ${\mathcal H}_\ell$ given by (11.6), and on the reconstructed internal flow (9.2)–(9.3) from (6.5)–(6.6), with ${\mathcal H}_\ell$ given by (11.8), are displayed in both cases $H = 0.02$ , $0.04$ . These two sets of values are found to be comparable in magnitude and to evolve similarly in time. Overall, the very low levels attained by these errors (around $10^{-3}$ %) as well as their stable evolution confirm that $\mathcal H$ is very well conserved. Not only does this test verify the existence of such a conserved quantity for Lagrangian particles under a travelling wave, but it also further attests to the good performance of our numerical solver (flow field reconstruction and time integration) for the particle dynamics.

Figure 9. Relative error on Hamiltonian $\mathcal H$ as a function of time $t$ for $16$ Lagrangian particles along the water column under a linear travelling wave of height (a) $H = 0.02$ and (b) $H = 0.04$ with wavenumber $k = 2/5$ . Numerical estimates based on the exact (blue line) and reconstructed (red line) velocity field are shown. (a) $H = 0.02$ , (b) $H = 0.04$ .

12. Nonlinear wave setup and comparison with KdV theory

We now perform direct simulations of (6.5)–(6.6) in the presence of a wave setup or setdown and compare them with weakly nonlinear predictions by the KdV equation. The main goal is to highlight differences in the flow dynamics between the setup and setdown configurations. We will see that, although the surface wave dynamics in these two cases look similar, their respective internal flow dynamics can be quite distinct, even for a small difference in the mean water level. We focus on the shallow water regime where conditions are more suitable for a wave setup or setdown to occur, e.g. due to effects from currents as commonly observed in nearshore areas.

For this comparison, we use exact cnoidal wave solutions of the KdV equation, which are periodic travelling waves with analytical expressions for the surface elevation and velocity field. Such explicit formulae can easily accommodate perturbations due to a wave setup or setdown, as discussed in § 5.

On the other hand, to tackle this problem using the full nonlinear model, the equations of motion for the surface wave dynamics as well as the reconstruction algorithm for the internal velocity field need to be adjusted. The point of this adjustment is to include additional effects from a uniform horizontal flow $U$ as explained in § 2. Accordingly, the extended form of (6.5)–(6.6) is

(12.1) \begin{align} \eta _t & = G(\eta ) \psi - U \eta _x, \end{align}

(12.2) \begin{align} \psi _t & = -g \eta - U \psi _x + \frac {1}{2} \psi _x^2 + \frac {1}{2} \frac {\big [ G(\eta ) \psi + \eta _x \psi _x \big ]^2}{1 + \eta _x^2}, \end{align}

where $U$ contributes to linear advection, as presented in Guyenne (Reference Guyenne2019). This is still a canonical Hamiltonian system (6.1) for $\eta$ and $\psi$ with Hamiltonian

(12.3) \begin{align} \mathbb{H} = \frac {1}{2} \int _{-\infty }^\infty \big [ \psi G(\eta ) \psi + 2 U \eta \psi _x + g \eta ^2 \big ] {\rm d}x. \end{align}

Figure 10. Stokes drift velocity $u_L$ for $16$ Lagrangian particles along the water column $z$ under an unperturbed nonlinear travelling wave ( $s = 0$ ). Different cases of $(H,k)$ are shown: (a) $(0.1,2/5)$ , (b) $(0.2,2/5)$ , (c) $(0.3,2/5)$ , (d) $(0.2,8/5)$ . Numerical predictions from the Euler system (blue), KdV equation (red) and linear approximation (black) are compared.

Figure 11. Stokes drift velocity $u_L$ for $16$ Lagrangian particles along the water column $z$ under a nonlinear travelling wave with setup ( $s = +0.025$ ). Different cases of $(H,k)$ are shown: (a) $(0.1,2/5)$ , (b) $(0.2,2/5)$ , (c) $(0.3,2/5)$ , (d) $(0.2,8/5)$ . Numerical predictions from the Euler system (blue), KdV equation (red) and linear approximation (black) are compared.

Figure 12. Stokes drift velocity $u_L$ for $16$ Lagrangian particles along the water column $z$ under a nonlinear travelling wave with setdown ( $s = -0.025$ ). Different cases of $(H,k)$ are shown: (a) $(0.1,2/5)$ , (b) $(0.2,2/5)$ , (c) $(0.3,2/5)$ , (d) $(0.2,8/5)$ . Numerical predictions from the Euler system (blue), KdV equation (red) and linear approximation (black) are compared.

Figure 13. Surface profile $\eta$ for an unperturbed nonlinear travelling wave with (a) $(H,k) = (0.3,2/5)$ and (b) $(H,k) = (0.2,8/5)$ at $t = 0$ . Initial conditions for the Euler system (blue) and KdV equation (red) are compared.

The reconstruction method for $(u,v)$ is adjusted via the surface data in the matching condition (9.4), so that the linear system (9.5) to be solved for $X = \{ k_n q_n(t) \}$ now becomes

(12.4) \begin{equation} A X = B, \quad B = \{ u^{(s)}(x_j,t) + U \}. \end{equation}

The initial conditions are Stokes waves $(\eta _0(x),\psi _0(x))$ computed by Fenton’s method and perturbed by a small constant shift $s$ of the water level, so that

(12.5) \begin{align} \eta (x,0) = \eta _0(x) + s, \quad \psi (x,0) = \psi _0(x), \end{align}

with $s \gt 0$ (respectively, $s \lt 0$ ) for a wave setup (respectively, setdown). The induced mean flow is given by $U = -c + \sqrt {c^2 + 2 g s}$ as derived in § 2, where $c$ denotes the speed of the unperturbed Stokes wave, and it is prescribed in (12.1)–(12.4). Such a perturbation of the Stokes wave is not an exact travelling wave solution of (12.1)–(12.2) and produces some minor unsteadiness, especially at early stages of the simulation while the numerical solution transitions to a more steady state. In contrast, cnoidal waves with setup or setdown remain exact travelling wave solutions of the KdV equation, by exploiting the associated explicit expressions for the surface elevation and velocity field.

Considering the same domain as in the previous section ( $h = 1$ , $L = 10\pi$ ), figures 10–12 plot the Stokes drift velocity $u_L$ for 16 fluid particles initially located under the central wave crest at $x = L/2$ (from near the bottom up to the free surface). We examine $s = 0$ and $s = \pm 0.025$ in four different cases: $H = 0.1$ , $0.2$ , $0.3$ with $k = 2/5$ (two wavelengths over $0 \leqslant x \lt L$ ), and $H = 0.2$ with $k = 8/5$ (eight wavelengths over $0 \leqslant x \lt L$ ). Simulations are run until $t = T = 40$ (respectively, $T = 20$ ) in cases with $k = 2/5$ (respectively, $k = 8/5$ ). For convenience, we will focus on $s = \pm 0.025$ to illustrate our results about a wave setup or setdown, but we have made similar observations with other small values of $s$ .

For these computations, the spatial and temporal resolutions are set to $\Delta x = 0.12$ and $\Delta t = 0.01$ . Unlike the previous tests, a very small time step is not needed to accurately describe the surface wave and internal flow in the present cases with larger wave steepnesses. For each configuration $(H,k,s)$ , we compare estimates of $u_L$ from the Euler system (12.1)–(12.2), the KdV equation (5.1) and the linear approximation (4.8). Recall that the KdV equation is a weakly nonlinear and weakly dispersive model for water waves, hence $(H,k) = (0.3,2/5)$ and $(0.2,8/5)$ represent more discriminating cases for this model because the former is a more nonlinear regime whereas the latter is a more dispersive regime. The Euler and KdV wave profiles at $t = 0$ for these parameter values are depicted in figure 13 and small discrepancies can indeed be discerned around the wave crests or troughs.

In each configuration determined by $s$ , all curves of $u_L$ are displayed over the same range of values to highlight differences between the various cases of wave parameters. As expected, results from these three models tend to coincide when the wave steepness is decreased, and they also tend to be more uniform along the water column. As the wave steepness is increased (with increasing $H$ and fixed $k$ ), the Euler and KdV predictions remain close together while deviating from the linear approximation. For higher wavenumber (with increasing $k$ and fixed $h$ ), the Stokes drift velocity exhibits more variation in the vertical direction and is larger in magnitude. Additionally, notable differences are observed in this limit depending on the value of $s$ .

For $kh = 8/5$ and $s = 0$ , the Euler and linear predictions on $u_L$ look very much alike, while the KdV graph is more concave. This higher concavity reflects the fact that the KdV equation becomes less accurate in resolving the vertical structure of the flow in a more dispersive wave regime (as $kh$ is increased). This feature is clearly revealed by an inflection of the KdV graph near the free surface. Indeed, we find that integration of (5.5) for the KdV equation leads to a particle trajectory that (while lying initially on the free surface) gradually drops below the cnoidal wave elevation over time. This issue is only noticeable for the KdV equation and arises in all three configurations $s = 0$ , $\pm 0.025$ . To highlight this decline of the surface particle trajectories, we substitute the value of the vertical position $\zeta$ obtained from (5.5) for the highest particle by its interpolation along the cnoidal wave profile $\eta$ , hence the slight inflection of the KdV curves near the free surface as shown in figures 10–12(d).

Compared with $s = 0$ , the picture is reversed for $s = \pm 0.025$ ( $kh = 8/5$ ) with the KdV and linear estimates being closer together while the Euler estimates are larger in magnitude. Such a contrast is due to the higher levels of wave dispersion and nonlinearity in this case ( $H = 0.2$ , $k = 8/5$ ). The fact that the KdV and linear datasets are close together may be explained by the definition (4.5) of the mean-flow correction to the linear velocity field (4.1) in the presence of a setup or setdown, which is given in terms of the long-wave approximation $c_0 = \sqrt {gh}$ as discussed in § 4. Note that the setup $s = +0.025$ induces a positive (i.e. forward) drift with $x_L \gt 0$ and $u_L \gt 0$ in the direction of wave propagation, whereas the setdown $s = -0.025$ generates a negative (i.e. backward) drift with $x_L \lt 0$ and $u_L \lt 0$ . This drift is strongest at the free surface (respectively, near the bottom) for a wave setup (respectively, setdown). Overall, these values of $u_L$ are found to be comparable in magnitude between the two configurations $s = \pm 0.025$ and are consistent with the correction $\sim g s/c_0$ as appearing in the linear estimate (4.2).

Figure 14. Stokes drift velocity $u_L$ for $16$ Lagrangian particles along the water column $z$ under a nonlinear travelling wave for different values of $s$ . Two cases of $(H,k)$ are shown: (a) $(0.3,2/5)$ , (b) $(0.2,8/5)$ . Numerical predictions for $s = +0.025$ (blue), $s = -0.025$ (red) and $s = 0$ (black) are compared. (a) $H = 0.3, k = 2/5$ , (b) $H = 0.2, k = 8/5$ .

Figures 10–12 also reveal a small discrepancy in the wave elevation between the linear approximation and the two nonlinear solutions, which is indicated by a small gap in the position $z_0$ of their surface particles. This feature is especially pronounced for small $kh$ and is expected considering that nonlinear waves in this long-wave limit produce an appreciable mean flow, leading to a higher water level than in the linear case. We summarize these findings on $u_L$ for $(H,k) = (0.3,2/5)$ and $(0.2,8/5)$ in figure 14 where all three datasets $s = 0, \pm 0.025$ are plotted together in each of these wave regimes. The forward versus backward dichotomy of the Stokes drift when $s \lessgtr 0$ is evident from this figure. For each choice of $(H,k)$ , all three profiles of $u_L$ look similar along the water column, and the magnitude of the deviation from $s = 0$ is seen to be essentially the same for either $s = \pm 0.025$ , which correlates with $|s|$ as can be expected.

Our results on the Stokes drift velocity show that there is a significant qualitative difference in the internal flow dynamics between a wave setup and setdown, even for a relatively small shift of the water level and despite the fact that the surface wave dynamics look similar in both settings. Moreover, we find that the KdV predictions are robust over a wide range of parameter values in this shallow-water regime, agreeing well with direct computations based on the Euler system. To further illustrate these points, figures 15 and 16 depict the trajectories of the 16 Lagrangian particles for $s = 0$ , $\pm 0.025$ in the more discriminating cases $(H,k) = (0.3,2/5)$ and $(0.2,8/5)$ . We now restrict ourselves to analysing and comparing the KdV and Euler predictions. We see that these particles typically follow a trochoidal path with cyclic arches and loops. For $(H,k) = (0.3,2/5)$ , the particle paths from both models look quite similar whereas, for $(H,k) = (0.2,8/5)$ , discrepancies are more discernible, which is consistent with our previous observations when inspecting $u_L$ . In this case, over the same duration $0 \leqslant t \leqslant T$ and for a given $z_0$ , the Euler particle paths yield a longer horizontal drift $x_T$ with an additional arch, as compared with the KdV ones. On the other hand, each KdV trochoidal cycle seems to be longer in arclength than its Euler counterpart. Here, $x_T$ is defined as the total horizontal drift

(12.6) \begin{align} x_T = \xi _\ell (T) - \xi _\ell (0), \end{align}

over $0 \leqslant t \leqslant T$ for any particle $\ell$ .

Figure 15. Trajectories (black) of $16$ Lagrangian particles along the water column under a nonlinear travelling wave $(H,k) = (0.3,2/5)$ from $t = 0$ (red) to $t = 40$ (blue). Configuration: $s = 0$ (a,b), $s = +0.025$ (c,d), $s = -0.025$ (e,f). Numerical predictions from the KdV equation (a,c,e) and Euler system (b,d,f) are shown. (a) $\textrm{KdV}, s = 0$ , (b) $\textrm{Euler}, s = 0$ , (c) $\textrm{KdV}, s = +0.025$ , (d) $\textrm{Euler}, s = +0.025$ , (e) $\textrm{KdV}, s = -0.025$ , (f) $\textrm{Euler}, s = -0.025$ .

Figure 16. Trajectories (black) of $16$ Lagrangian particles along the water column under a nonlinear travelling wave $(H,k) = (0.2,8/5)$ from $t = 0$ (red) to $t = 20$ (blue). Configuration: $s = 0$ (a,b), $s = +0.025$ (c,d), $s = -0.025$ (e,f). Numerical predictions from the KdV equation (a,c,e) and Euler system (b,d,f) are shown.

Figures 15 and 16 confirm that this total drift $x_T$ is positive (respectively, negative) for a wave setup (respectively, setdown). In both settings, the drift magnitude is found to be appreciable along the entire water column, even near the bottom where it turns out to be largest for a wave setdown. This phenomenon is even more apparent in the more dispersive wave regime ( $H = 0.2$ , $k = 8/5$ ) where the distinction between the unperturbed ( $s = 0$ ) and perturbed ( $s \neq 0$ ) configurations is particularly striking. Indeed, the particle paths for $s = 0$ exhibit the typical pattern of the deep-water regime where the fluid motion rapidly decays with the depth. By contrast, the particle paths for $s \neq 0$ bear closer resemblance to a shallow-water flow. Such perturbations of the internal flow dynamics in the presence of a wave setup or setdown have important implications for mass and momentum transfer (both horizontally and vertically) below the water surface, including effects on sediment transport at the seafloor.

Figure 17. Zoom-in on trajectories of Lagrangian particles at the free surface (a,b) and near the bottom (c,d) under a nonlinear travelling wave $(H,k) = (0.3,2/5)$ from $t = 0$ (dot) to $t = 40$ (triangle). Configuration: $s = +0.025$ (a,c), $s = -0.025$ (b,d). Numerical predictions from the KdV equation (red) and Euler system (blue) are compared. (a) $s = +0.025$ , at the free surface, (b) $s = -0.025$ , at the free surface, (c) $s = +0.025$ , near the bottom, (d) $s = -0.025$ , near the bottom.

Figure 18. Zoom-in on trajectories of Lagrangian particles at the free surface (a,b) and near the bottom (c,d) under a nonlinear travelling wave $(H,k) = (0.2,8/5)$ from $t = 0$ (dot) to $t = 20$ (triangle). Configuration: $s = +0.025$ (a,c), $s = -0.025$ (b,d). Numerical predictions from the KdV equation (red) and Euler system (blue) are compared. (a) $s = +0.025$ , at the free surface, (b) $s = -0.025$ , at the free surface, (c) $s = +0.025$ , near the bottom, (d) $s = -0.025$ , near the bottom.

The setup and setdown configurations are further distinguished by noting that the particle paths are trochoids of prolate type with concave-down arches for $s \gt 0$ , whereas they are of prolate type with concave-up arches for $s \lt 0$ . This can be seen in figures 17 and 18 that show close-ups with a direct comparison of particle paths predicted by the KdV equation and the Euler system for $s = \pm 0.025$ in the two previous cases of $(H,k)$ . Discrepancies between these two approaches for $(H,k) = (0.2,8/5)$ are particularly apparent in these close-ups, and are more significant near the bottom where the Euler particle trajectories display more arches but with lower amplitude than their KdV counterparts. In the setdown configuration, especially for $(H,k) = (0.3,2/5)$ , the return loop after each cycle turns out to be quite prominent, thus contributing to a shorter drift $x_T$ than in the setup configuration, as demonstrated by figures 11 and 12. Consistent with our previous comments on figures 10–12, the good agreement between Euler and KdV particle dynamics for $(H,k) = (0.3,2/5)$ in figure 17 may serve to provide a cross-validation of these two strategies, assessing on one hand the regime of validity for the KdV equation and on the other hand the reconstruction method for the fluid velocity with a mean flow via the Euler system.

Again, such differences in the geometry of particle orbits between a wave setup and setdown may have important implications for mass and momentum transfer under ocean surface waves, especially in nearshore areas. As noted earlier, the fact that the Euler solution yields a longer drift $x_T$ than its KdV counterpart for $(H,k) = (0.2,8/5)$ with both $s = \pm 0.025$ is clearly revealed in the close-ups of figure 18. Overall, these observations on $x_T$ over the interval $0 \leqslant t \leqslant T$ correlate well with our previous results on $u_L$ ( $x_L$ over $\tau _L$ ) as presented in figures 11 and 12.