2.1. Mean Lagrangian drift of general flows

While directly solving the Euler (2.7)–(2.8) subject to ${\mathcal{J}}=1$ , $q=0$ and the boundary conditions (2.11)–(2.12) will yield particle trajectories that explicitly contain the mean Lagrangian drift, this offers little physical insight into its origin. Previous studies connected the mean Lagrangian drift, or equivalently the mean Lagrangian momentum density, to other physical quantities such as vorticity and energy (Pizzo et al. Reference Pizzo, Murray, Smith and Lenain2023; Blaser et al. Reference Blaser, Benamran, Bôas, Lenain and Pizzo2024), but these results necessarily assumed waves that were steady and monochromatic. In this section, we introduce a new method of constraining the mean Lagrangian drift for completely general flows. To do so, we start by considering the circulation of a material loop $\mathcal{C}$ , which is defined as

(2.13) \begin{equation} \varGamma \equiv \oint _{\mathcal{C}} \boldsymbol{u} \boldsymbol{\cdot }{\textrm {d}} \boldsymbol{\ell } = \iint _{\varOmega } (\boldsymbol{\nabla } \times \boldsymbol{u}) \boldsymbol{\cdot } \boldsymbol{\hat {n}} \, {\textrm {d}} x {\textrm {d}} y , \end{equation}

with the last relation due to Stokes’ theorem for the area $\varOmega$ enclosed by the contour where $\boldsymbol{\hat {n}}$ is the unit outward normal. We simplify here to two-dimensional flow, but the following results may be readily extended to three dimensions (Salmon Reference Salmon1988). Just as with the vorticity, we can rewrite the circulation in Lagrangian coordinates via the chain rule

(2.14) \begin{equation} \varGamma = \oint _{\mathcal{C}_0} \big (x_\tau \boldsymbol{\nabla _\alpha }x + y_\tau \boldsymbol{\nabla _\alpha }y\big ) \, \boldsymbol{\cdot } {\textrm {d}} \boldsymbol{\alpha } = \iint _{\varOmega _0} {\mathcal{J}} q \, {\textrm {d}} \alpha \, {\textrm {d}} \beta , \end{equation}

where $\boldsymbol{\nabla _\alpha } = ({\partial }_\alpha ,{\partial }_\beta )$ is the gradient operator in label space. The contour $\mathcal{C}_0$ is now a contour in label space and is therefore fixed in time by definition; the same goes for the enclosed area $\varOmega _0$ . If we decompose the Lagrangian trajectories into deviations from a quiescent reference state $(x,y) = (\alpha ,\beta )$

(2.15) \begin{equation} x = \alpha + \xi (\alpha ,\beta ,\tau ) , \quad y = \beta + \zeta (\alpha ,\beta ,\tau ), \end{equation}

so that $\alpha$ and $\beta$ can be seen as ‘horizontal’ and ‘vertical’ labels respectively, we can rewrite the circulation as

(2.16) \begin{equation} \varGamma = \oint _{\mathcal{C}_0} (\boldsymbol{x}_\tau - \boldsymbol{p}) \boldsymbol{\cdot } {\textrm {d}} \boldsymbol{\alpha } , \end{equation}

where $\boldsymbol{x}_\tau = (\xi _\tau ,\zeta _\tau )$ is the Lagrangian velocity, and $\boldsymbol{p} =-(\xi _\tau \xi _\alpha + \zeta _\tau \zeta _\alpha , \, \xi _\tau \xi _\beta + \zeta _\tau \zeta _\beta )$ is identified as the Lagrangian pseudomomentum. While its form looks identical to pseudomomentum as defined in GLM theory (Andrews & McIntyre Reference Andrews and McIntyre1978; Bühler Reference Bühler2014), they are still distinct since the displacement vector in GLM is a function of the Lagrangian-mean trajectory, not Lagrangian particle labels. For irrotational flows where $\varGamma =0$ for all closed loops, (2.16) implies that the label space curl of the velocity must be everywhere equal to the label space curl of the pseudomomentum

(2.17) \begin{equation} \boldsymbol{\nabla }_{\boldsymbol{\alpha }} \times \boldsymbol{x}_\tau = \boldsymbol{\nabla }_{\boldsymbol{\alpha }}\times \boldsymbol{p}, \end{equation}

analogous to the celebrated result in GLM (Bühler Reference Bühler2014, Ch. 10). Since what we are interested in is the mean component of the velocity, we can take an average of (2.17) to get

(2.18) \begin{equation} \boldsymbol{\nabla }_{\boldsymbol{\alpha }} \times \langle \boldsymbol{x}_\tau \rangle = \boldsymbol{\nabla }_{\boldsymbol{\alpha }} \times \langle \boldsymbol{p} \rangle , \end{equation}

where the angle brackets represent any general averaging operator in the Lagrangian frame that commutes with the curl in label space, such as a time mean or convolutional average. It is worth pausing here for a moment to unpack this result, which states that, for irrotational flow, the curl of the mean Lagrangian drift is exactly set by the curl of the mean pseudomomentum so that any modification to one immediately affects the other. Viewing the mean Lagrangian drift in relation to the dynamic mean pseudomomentum highlights its role as not simply a passive byproduct of the waves, but as a dynamic mean flow in its own right. This view will be especially helpful when we turn to the mean Lagrangian drift of narrow-banded wave packets. However, for completeness, we will use this new general framework to compute the mean Lagrangian drift for linear waves in the following subsections.

2.2. Monochromatic waves

We start with the classical example of a linear deep-water monochromatic wave with wavenumber $k$ and constant amplitude $A$ where the non-dimensional steepness $Ak$ is assumed to be small. Following the method of (Salmon Reference Salmon2020, § 1), we assume a wavelike solution for $x$ and $y$ after expanding about a hydrostatic state of rest ( $x=\alpha$ , $y=\beta$ , $p = -g\beta$ )

(2.19) \begin{align} x(\alpha ,\beta ,\tau ) &= \alpha - A e^{k \beta } \sin (k\alpha - \omega \tau ), \\[-12pt]\nonumber \end{align}

(2.20) \begin{align} y(\alpha ,\beta ,\tau ) &= \beta + Ae^{k \beta }\cos (k\alpha - \omega \tau ), \\[-12pt]\nonumber \end{align}

(2.21) \begin{align} p(\beta ) &= -g \beta , \end{align}

where $\omega ^2 = g k$ is the linear deep-water dispersion relation determined by substituting (2.19)–(2.21) into the Euler (2.7)–(2.8). These simple circular trajectories are in fact exact solutions to the Euler equations known as Gerstner (Reference Gerstner1802) waves. However, these waves are not irrotational, which can be seen by computing their vorticity using (2.9). From (2.18), we see that irrotational flow requires that the curl of the mean Lagrangian drift be equal to the curl of the mean pseudomomentum. Computing the pseudomomentum of (2.19)–(2.20) yields only a horizontal component

(2.22) \begin{equation} \mathrm{p}^{(x)} = -(\xi _\tau \xi _\alpha + \zeta _\tau \zeta _\alpha ) = A^2 k \omega e^{2 k \beta }, \end{equation}

that varies only with depth. Taking the mean to be a long time average following a fixed particle, from (2.18) we require

(2.23) \begin{equation} \frac {{\partial } \langle y_\tau \rangle }{{\partial } \alpha } - \frac {{\partial } \langle x_\tau \rangle }{{\partial } \beta } = -\frac {{\partial } \langle \mathrm{p}^{(x)}\rangle }{{\partial } \beta } = -2 A^2 k^2 \omega e^{2 k \beta }. \end{equation}

On physical grounds we can assume there is no mean vertical velocity, so that the solution to (2.23) is

(2.24) \begin{equation} \langle x_\tau \rangle = A^2 k \omega e^{2 k \beta }, \end{equation}

where the arbitrary constant of integration can be removed in the frame where the velocity of fluid at depth vanishes. This the classical Stokes drift. We reproduce it here as an example of our general method but also because it shows how the second-order mean flow is constrained by first-order orbital motion, due to the pseudomomentum being a quadratic quantity. This carries to higher-order corrections as well; since the particle displacements in surface gravity wave fields $(\xi ,\zeta )$ are always first-order quantities or higher, one needs only to constrain trajectories valid to order $n$ to constrain the drift to order $n+1$ .

2.3. Multiple waves – lowest-order theory

Following Pierson (Reference Pierson1961), if we instead consider a discrete spectrum of $N$ deep-water plane waves travelling in the same direction, to first order we have

(2.25) \begin{align} x &= \alpha - \sum _{n=1}^N A_n e^{k_n \beta } \sin (k_n\alpha - \omega _n \tau + \varphi _n), \\[-12pt]\nonumber \end{align}

(2.26) \begin{align} y &= \beta + \sum _{n=1}^N A_n e^{k_n \beta } \cos (k_n \alpha - \omega _n \tau + \varphi _n), \end{align}

where $A_n$ , $k_n$ , $\omega _n$ and $\varphi _n$ are the amplitude, wavenumber, frequency and initial phase of each wave component, respectively. It is assumed that each wave’s steepness $A_n k_n$ is small and, importantly for this analysis, constant. The horizontal component of the pseudomomentum is given by products of sums, but taking the mean to be a long time average following a fixed particle, we have

(2.27) \begin{equation} \langle \mathrm{p}^{(x)}\rangle =-\langle \xi _\tau \xi _\alpha + \zeta _\tau \zeta _\alpha \rangle = \sum _{n=1}^N A_n^2 k_n \omega _n e^{2 k_n \beta }, \end{equation}

where any cross-terms vanish in the time mean due to to the fact that $A_n$ , $k_n$ , $\omega _n$ and $\varphi _n$ are constant. From (2.18),

(2.28) \begin{equation} \frac {{\partial } \langle y_\tau \rangle }{{\partial } \alpha } - \frac {{\partial } \langle x_\tau \rangle }{{\partial } \beta } = -\frac {{\partial } \langle \mathrm{p}^{(x)} \rangle }{{\partial } \beta } \ = -\sum _{n=1}^N 2A_n^2 k_n^2 \omega _n e^{2 k_n \beta }. \end{equation}

Once again, the only physically valid solution is balanced by $\langle x_\tau \rangle$ , so that

(2.29) \begin{equation} \langle x_\tau \rangle = \sum _{n=1}^N A_n^2 k_n \omega _n e^{2 k_n \beta }, \end{equation}

where the constant of integration vanishes in the frame where the fluid interior is at rest. We see that according to lowest-order theory, the total mean Lagrangian drift for a wave field is a simple linear sum of the individual drifts of each wave component. While the full second-order particle trajectory solutions contain bounded second-order harmonics (Pierson Reference Pierson1961; Nouguier, Chapron & Guérin Reference Nouguier, Chapron and Guérin2015), which can be interpreted as local fluctuations to the mean Lagrangian drift, these terms are fully oscillatory and do not contribute to the long time transport regardless of how the initial phases $\varphi _n$ are tuned. We will henceforth refer to this as the linear theory, on account of its additive property, although this should not be confused with linearity with respect to wave steepness. From this theory, the effect of local steepness fluctuations to the mean Lagrangian drift is symmetric; any local increases during constructive interference are cancelled by local decreases during destructive interference. We would therefore expect the total transport of a passing wave packet, expressed as a sum of plane waves, to be similarly invariant to local wave focusing. In the following section, we investigate the Lagrangian transport of focusing wave packets, presenting numerically simulated particle trajectories alongside laboratory data.

3.1. Packet initialisation

We define our packets as in Rapp & Melville (Reference Rapp and Melville1990), Drazen, Melville & Lenain (Reference Drazen, Melville and Lenain2008) and Sinnis et al. (Reference Seitz, Freilich and Pizzo2021) to focus according to linear theory at a prescribed space and time

(3.1) \begin{equation} \eta (x,t) = \sum _{n=1}^N A_n \cos (k_n (x - x_{\!f}) - \omega _n (t - t_{\!f})), \end{equation}

where $\eta (x,t)$ is the Eulerian free surface displacement, $A_n$ is the amplitude of each discrete wave, $k_n$ and $\omega _n$ represent the respective wavenumber and frequency of each component, both positive as all wave components travel to the right, and $x_{\!f}$ and $t_{\!f}$ denote the focusing location and time respectively according to linear theory. We consider a uniformly distributed spectrum in frequency space, so that our frequencies can be expressed as

(3.2) \begin{equation} \omega _n = \omega _c \left(1 + \varDelta \left(\frac {n}{N} - \frac {1}{2}\right)\right)\!, \end{equation}

where $\omega _c$ is the central frequency (so that $k_c = \omega _c^2/g$ is the central wavenumber) and $\varDelta$ is the non-dimensional bandwidth equivalent to $(\omega _N - \omega _1)/\omega _c$ which sets the time and space scales of the focusing event and must be less than $2$ to ensure positive frequencies. In addition, as the slope of waves is an indicator of their nonlinearity, we wish to define the amplitudes $A_n$ such that at focusing, the linear prediction of the maximum slope equals some prescribed value $S$ . Therefore, we define

(3.3) \begin{equation} S = \sum _{n=1}^N A_n k_n . \end{equation}

Thus, given the linear deep-water dispersion relationship $\omega _n^2 = g k_n$ , we can determine the values of $A_n$ , assuming the slope of each mode is equal following Drazen et al. (Reference Drazen, Melville and Lenain2008) (i.e. $A_n = S/(Nk_n)$ ). Placed in this formulation, the wave packets we consider are primarily functions of two independent variables, $S$ and $\varDelta$ , which will be used as our parameter space.

The linear prediction of the surface mean Lagrangian drift is given by a simple sum of the lowest-order contributions (2.29)

(3.4) \begin{equation} \langle x_\tau \rangle \big |_{\beta =0} = \sum _n^N A_n^2 k_n \omega _n = \sum _n^N \frac {S^2}{N^2} \frac {\omega _n}{k_n} . \end{equation}

Based solely on (3.4), the surface mean Lagrangian drift scales as $ {1}/{N}$ , which implies that a packet with more waves experiences less drift, despite the fact that $N$ simply represents the spectral resolution of the packet, whose form converges as $N \rightarrow \infty$ . This is due to the fact that the temporal periodicity of (3.1) is given by

(3.5) \begin{equation} T_{\!p} = \frac {2\pi N}{\omega _c \varDelta }, \end{equation}

so that as $N$ increases, the time between subsequent packets also increases. Since what we are after is not the mean Lagrangian drift itself, which according to (2.29) is the same for all particles at all times since it treats the packet as a sum of monochromatic plane waves, we instead compute the total linear surface Lagrangian transport $\delta x_{\textit{lin}}$ after a single packet has passed. This is done by integrating (3.4) in time over the temporal periodicity of the packet (3.5)

(3.6) \begin{equation} \delta x_{\textit{lin}} = \int _0^{T_{\!p}} \langle x_\tau \rangle \, {\textrm {d}}\tau \big |_{\beta =0}= \langle x_\tau \rangle T_{\!p} \big |_{\beta =0}= \frac {2\pi S^2}{\omega _c \varDelta } \left(\frac {1}{N}\sum _n^N \frac {\omega _n}{k_n} \right)\!. \end{equation}

We see that the linear prediction of the total surface Lagrangian transport should scale with $S^2$ as one should expect for the lowest-order theory. The transport scaling inversely with $\varDelta$ should also be expected as the packet width in physical space is inversely proportional to its width in wavenumber space (Sinnis et al. Reference Seitz, Freilich and Pizzo2021). The last term represents a spectrally weighted phase speed. For $N$ large, (3.6) can be approximated in closed form from (3.2) as

(3.7) \begin{equation} \delta x_{\textit{lin}} \approx \frac {2 \pi S^2}{k_c} f(\varDelta ), \quad f(\varDelta ) = \frac {1}{\varDelta ^2}\ln \left(\frac {1 + \varDelta /2}{1 - \varDelta /2} \right)\!, \end{equation}

where $f(\Delta )$ represents the linear bandwidth dependence on the total transport, found by approximating the sum in (3.6) as an integral. While this full expression is slightly more complicated than the heuristic argument given above, $f(\Delta )$ is well approximated by $\Delta ^{-1}$ when $\varDelta$ is small.

3.2. Numerical simulations of Lagrangian trajectories

To simulate the Lagrangian trajectories of surface particles within these packets, we employ a fully nonlinear mixed Eulerian–Lagrangian potential flow solver (Dold Reference Dold1992). Originally developed by Longuet-Higgins & Cokelet (Reference Longuet-Higgins and Cokelet1976), this method takes advantage of the fact that, at a fixed time, the Eulerian and Lagrangian velocities are equal since a particle occupies a single fixed location at a fixed time. Because solutions to Laplace’s (2.1) are uniquely determined by the boundary conditions, only the surface needs to be simulated, assuming a constant or infinite depth and a periodic domain. By initialising Lagrangian particles with initial positions $(x_0,y_0)$ and velocity potential $\phi _0$ , this solver computes the gradient of $\phi$ at the surface given its value via Cauchy’s integral theorem at each time step. This allows for the particle positions $(x,y)$ to evolve via the pathline (2.4), with the velocity potential evolving according to Bernoulli’s equation at the free surface (2.3). Because Lagrangian particles naturally cluster at wave crests where the spatial curvature is strongest, the resolution of this method is naturally adaptive, and numerous studies (Dommermuth et al. Reference Dommermuth, Yue, Lin, Rapp, Chan and Melville1988; Skyner Reference Sinnis, Grare, Lenain and Pizzo1996) have validated the accuracy and validity of this numerical method.

For our simulations, we chose a central wave frequency of $1$ Hz, so that $\omega _c = 2\pi$ rad s⁻¹ and $k_c= (2\pi )^2/g$ rad m⁻¹. The domain length was chosen to be $L = 150$ m, long enough so that the entire packet could fully pass over a large enough collection of particles to obtain an unambiguous measure of total Lagrangian transport before any signal wrapped around due to the periodicity of the domain. To fully resolve the free surface, we systematically increased the number of Lagrangian particles used until convergence was reached at 2048 particles, or around 20 per central wavelength. The depth of the water is taken to be infinitely deep, and the packets were initialised to start $20$ m away from the prescribed focusing location so that there were sufficient particles within the focusing region that both started and ended at rest. We defined the linear prediction of the focusing time as $t_{\!f} = x_{\!f}/c_g$ , where $c_g = ({1}/{2})\sqrt {g/k_c}$ is the central group velocity according to linear theory. Lastly, we chose to use $N=1000$ wave modes so that the spectral resolution is sufficiently high to converge the physical shape of the packet.

The procedure for simulating these packets is as follows. First, we initialise the horizontal positions $x_0$ to be evenly spaced along the domain. Then, using (3.1), the vertical initial positions $y_0$ are found for prescribed values of $S$ and $\varDelta$ . To ensure only one packet is used and that the domain is totally periodic, a windowing function is applied to $y_0$ with minimal energy loss (less than one part in 100). The initial velocity potential is found according to linear theory by performing a Fourier transform on the windowed $y_0$ , multiplying each Fourier amplitude by $\omega _n/k_n$ , and performing an inverse transform. The simulations were repeated for a parameter space spanning $\varDelta = 0.3$ , incremented by $0.1$ until $\varDelta = 1.3$ , and $S = 0.05$ , incremented by $0.02$ until the packet broke, which agreed well with the results of Pizzo & Melville (Reference Pizzo and Wagner2019); see also Pizzo et al. (Reference Pizzo and Melville2021), who numerically investigated the breaking threshold $S_*$ of these same packets as a function of bandwidth. To improve parameter space resolution near $S_*$ , we ran additional simulations near this threshold.

Figure 1 shows a typical output of the surface particle trajectories during one such focusing event with bandwidth $\varDelta = 0.8$ and linear prediction of maximum slope at focusing $S = 0.27$ . For particles far downstream and upstream of focusing, represented by the blue and green curves respectively, their trajectories evolve gradually as the packet passes over. Their measured total transport $\delta x$ , represented by the difference from their final and initial positions, mostly follows linear theory (3.6). In contrast, for the particles at or near the focusing location, highlighted by the red curve, the transport occurs in one short burst as the focused packet passes over, well surpassing the predictions of linear theory and violating the supposed spatial invariance of the transport. For each particle, the total Lagrangian transport $\delta x$ is computed by taking the horizontal position averaged over the final two seconds of the simulation, roughly two central wave periods, and subtracting from it the particle’s fixed initial position $x_0$ . Both here and for the rest of this paper, we only show results for particles that began and ended at rest (i.e. they experienced the full packet passing) so that total transport $\delta x$ is unambiguous.

Figure 1.

Surface particle trajectories in a focusing wave packet with $\varDelta = 0.8$ and $S = 0.27$ . In panel $(a)$ , the vertical elevation of individual fluid particles is plotted as a function of time. Each curve represents a different particle, labelled by its initial horizontal distance from the focusing region $(x_0 - x_{\!f})$ , normalised by the central wavenumber $k_c$ , as shown on the vertical axis. The coloured lines represent particles downstream of focusing (blue), at focusing (red) and upstream of focusing (green). Likewise, on the right, panels (b,c,d) show the physical particle trajectories of these downstream, at focusing, and upstream particles respectively, normalised by $k_c$ . Note that the total transport during focusing (red) is much greater than that away from focusing, contrary to linear theory (3.6) (dashed line) which states that all particles should experience the same transport.

Figure 2.

The total Lagrangian transport $\delta x$ of surface particles as a function of their initial distance from the linear prediction of maximum focusing $(x_0 - x_{\!f})$ , normalised by the central wavenumber $k_c$ for the same simulation as in figure 1, $\varDelta = 0.8$ and $S = 0.27$ . The normalised linear prediction of the total transport $k_c \delta x_{\textit{lin}}$ (3.6), constant in space, is shown in red.

Figure 2 shows the computed total surface Lagrangian transport $\delta x$ as a function of its initial position relative to the linear focusing prediction, both normalised by the central wavenumber $k_c$ for the same simulation as in figure 1. Plotted also is $\delta x_{\textit{lin}}$ , normalised by $k_c$ , constant for each particle. From figure 2, we see a strong spatial dependence of the total Lagrangian transport, with a maximum transport 75 % higher than linear theory predicts. At larger values of $\varDelta$ , where the physical packet width at focusing is smaller, the maximum transport was even found to be up to double that of linear theory. These are surprising results as it might be expected that any higher-order corrections to linear theory would be necessarily small. Here, we show that these corrections are comparable in magnitude to the linear prediction and exhibit a strong spatial dependence. In addition to the general increase around the focusing location, all simulations have oscillations in their transport curve near focusing with a spatial periodicity that matches the wavelength of the central wave. To compare these results with laboratory experiments, we also introduce a measure of the mean surface transport over the focusing region following Sinnis et al. (Reference Seitz, Freilich and Pizzo2021),

(3.8) \begin{equation} \langle \delta x \rangle = \frac {1}{x_{02}-x_{01}}\int _{x_{01}}^{x_{02}} \delta x \, {\textrm {d}} x_{0}, \end{equation}

where in our study we define $x_{01}$ and $x_{02}$ as the first and last points, respectively, where the deviation of the transport from linear theory $(\delta x - \delta x_{\textit{lin}})$ exceeds 10 % of its maximum value. In this case the mean transport is $23.6\,\%$ higher than linear theory. Although this particular definition of the mean transport is by no means unique, it was chosen to best match the approach used in Sinnis et al. (Reference Seitz, Freilich and Pizzo2021).

While these simulations provide the first detailed account of the increased transport of steep non-breaking focusing wave packets, this study was motivated by earlier laboratory experiments (Lenain et al. Reference Lenain, Pizzo and Melville2019; Sinnis et al. Reference Seitz, Freilich and Pizzo2021). These wave tank experiments measured the spatially varying surface transport in primarily breaking focusing wave packets described by (3.1), with several steep non-breaking cases included for comparison. They found that wave breaking produces a large local increase to the surface transport. Wave breaking, in this case, breaks both the translational symmetry of the system and the transport in an obvious way. However, this symmetry breaking is also present for non-breaking focusing waves, allowing for a spatially dependent non-breaking transport which can be seen for example in Sinnis et al. (Reference Seitz, Freilich and Pizzo2021) figure 5.

Figure 3 $(a)$ compares the mean surface transport from laboratory experiments and simulations as a function of $S$ , normalised by the central wavenumber $k_c$ and linear bandwidth dependence $f(\Delta )$ so that the linear theory (3.7) coincides for both cases. For the laboratory experiments, conducted in a finite-depth tank of mean depth $\ell = 0.5$ m, this requires using the full dispersion relationship $\omega ^2 = g k \tanh (k \ell )$ to numerically compute $f(\Delta )$ . Despite these differences, a clear trend emerges: the mean surface transport exceeds the predictions of linear theory as $S$ increases. To better visualise these increases, figure 3(b) plots the same data instead as a percentage deviation from linear theory, revealing that these enhancements are of comparable magnitude to the linear theory itself.

Figure 3.

Mean surface transport $\langle \delta x \rangle$ as a function of the linear prediction of maximum wave slope $S$ . Panel $(a)$ shows the mean transport normalised by the central wavenumber $k_c$ and linear bandwidth dependence $f(\Delta )$ so that the prediction of linear theory (3.6) (red) collapses to a single curve for both the simulation and laboratory parameters. A polynomial fit of the discrete simulation points is shown in green. Panel $(b)$ shows the same data plotted as a percentage increase from linear theory.

Figure 4.

Percentage increases of the maximum $(a)$ and mean $(b)$ surface Lagrangian transport relative to linear theory (3.6) for numerically simulated focusing wave packets as a function of parameter space $(S,\varDelta )$ . Discrete simulation runs are shown via coloured markers, with interpolated values in between. Note the two distinct colour bar scalings for panels (a,b). The red line outlining the parameter space represents the breaking slope threshold numerically determined by Pizzo et al. (Reference Pizzo and Melville2021) which we found to be consistent with our simulations.

The enhancements of the maximum and mean surface transport relative to linear theory for all simulations are shown in figure 4 as discrete points, with interpolated values in between. The dashed red line indicates the breaking slope threshold $S_*$ numerically determined by Pizzo et al. (Reference Pizzo and Melville2021) for equivalently defined deep-water packets. Individual particles, shown in figure 4 $(a)$ , can be transported up to twice as far as linear theory predicts, with the mean transport over the focusing region surpassing linear theory by up to $30\,\%$ as shown in figure 4 $(b)$ . Although the enhancements to the surface Lagrangian transport primarily scale with increasing $S$ , there is also a noticeable $\varDelta$ dependence close to the breaking threshold. To investigate why these spatially varying enhancements occur when waves steepen, we next turn to a theoretical derivation of the local mean Lagrangian drift in deep-water narrow-banded waves.

4.1. The mean Lagrangian drift of steep narrow-banded waves

Because we enforced a scale separation between the fast orbital motion and the slow envelope evolution through the small parameter $\gamma$ , the averaging operator $\langle \boldsymbol{\cdot }\rangle$ is defined as a spatial convolution over an intermediate scale – large enough to remove the fast oscillations but small enough to retain slow envelope modulations. Spatial convolutions commute with the curl, so that from (2.18) we know that the curl of the mean Lagrangian drift is exactly equal to the curl of the mean pseudomomentum. The pseudomomentum itself is computed from products of the Lagrangian particle displacements $(\xi ,\zeta )$ whose leading-order contributions are $O(\epsilon )$ ; consequently, the fourth-order pseudomomentum is determined entirely by third-order terms. It is easy to check that the second-order mean terms in (4.9)–(4.10) only contribute to the curl of the mean pseudomomentum beginning at fifth order, so that only the oscillatory terms are required. Direct evaluation of this quantity from (4.13)–(4.14) shows that the mean Lagrangian drift must satisfy

(4.21) \begin{align} &\frac {{\partial } \langle y_\tau \rangle }{{\partial } \alpha } - \frac {{\partial } \langle x_\tau \rangle }{{\partial } \beta }= -2 \epsilon ^2 |{\mathcal{F}}|^2 e^{2\beta } + \epsilon ^2\gamma \Big(2i (1 + \beta ) \big ({\mathcal{F}}_{\mathcal{A}} {\mathcal{F}}^* - {\mathcal{F}} {\mathcal{F}}_{\mathcal{A}}^*\big ) - i \big ({\mathcal{F}}_T{\mathcal{F}}^* - {\mathcal{F}} {\mathcal{F}}_T^*\big )\Big )e^{2\beta } \nonumber \\ &\quad - 8 \epsilon ^4 |{\mathcal{F}}|^4 e^{4\beta } + \epsilon ^2 \gamma ^2 \left(\left(\frac {1}{2} + \frac {5}{2}\beta + \beta ^2 \right)|{\mathcal{F}}|^2_{{\mathcal{A}}{\mathcal{A}}} e^{2\beta } - \big (4 + 10\beta + 4\beta ^2\big )|{\mathcal{F}}_{\mathcal{A}}|^2 e^{2\beta } \right)\!, \end{align}

where we have used the identity

(4.22) \begin{equation} \big ({\mathcal{F}}_{{\mathcal{A}}{\mathcal{A}}}{\mathcal{F}}^* + {\mathcal{F}} {\mathcal{F}}_{{\mathcal{A}}{\mathcal{A}}}^*\big )= |{\mathcal{F}}|^2_{{\mathcal{A}}{\mathcal{A}}} - 2|{\mathcal{F}}_{\mathcal{A}}|^2 , \end{equation}

to consolidate various terms. In its current form (4.21) is underdetermined since we cannot ignore ${\partial }_\alpha \langle y_\tau \rangle$ at higher orders. Although the fluid is incompressible, the mean Lagrangian velocity need not be divergence free (2.6). As pointed out by Vanneste & Young (Reference Sulem and Sulem2022), the fact that waves modify the potential energy of a fluid implies a changing centre of mass, which, when paired with the bottom boundary condition, requires a divergent Lagrangian-mean flow. We can account for this, however, in the Lagrangian-mean water level, whose value can be computed from the Jacobian (2.5) independently of the mean Lagrangian drift to fourth order from (4.13)–(4.14)

(4.23) \begin{align} \langle y\rangle _{\textit{MWL}} = \frac {1}{2}\epsilon ^2|{\mathcal{F}}|^2 e^{2\beta } - \frac {1}{2}\epsilon ^2 \gamma i \bigg (\beta + \frac {1}{2}\bigg )\big ({\mathcal{F}}_{{\mathcal{A}}}{\mathcal{F}}^* - {\mathcal{F}} {\mathcal{F}}_{\mathcal{A}}^*\big ) e^{2\beta } + \frac {3}{2}\epsilon ^4 |{\mathcal{F}}|^4 e^{4\beta }\nonumber \\ + \epsilon ^2 \gamma ^2 \bigg ( \big(\beta + \beta ^2\big)|{\mathcal{F}}_{\mathcal{A}}|^2 + \frac {1}{8}\big(1 - 2\beta -2\beta ^2 \big) |{\mathcal{F}}|_{{\mathcal{A}}{\mathcal{A}}}^2 \bigg )e^{2\beta } . \end{align}

The Lagrangian-mean water level (4.23), or more aptly changes thereof, fully constrain the divergent part of the total mean Lagrangian velocity. This can be seen by taking the average of the incompressibility condition (2.6), which from (4.13), (4.14) and (4.23) imply

(4.24) \begin{equation} \langle x_\tau \rangle _\alpha + \langle y_\tau \rangle _\beta = O(\epsilon ^4 \gamma ), \end{equation}

so that, to fourth order, the mean Lagrangian drift is divergence free (at higher orders products of the mean Lagrangian drift enter (4.24)). This allows us to define a streamfunction for the mean Lagrangian drift $\psi$ such that

(4.25) \begin{equation} \langle x_\tau \rangle = -\psi _\beta, \quad \langle y_\tau \rangle = \psi _\alpha . \end{equation}

Equation (4.21) does not distinguish between the mean Lagrangian drift or the Lagrangian-mean water level, so to only constrain the drift we must account for the mean water level’s contribution to the curl, which only emerges at fourth order

(4.26) \begin{equation} {\partial }_\alpha {\partial }_\tau (\langle y \rangle _{\textit{MWL}}) = \frac {1}{2}\epsilon ^2 \gamma ^2|{\mathcal{F}}|^2_{{\mathcal{A}} T} \approx -\frac {1}{4}\epsilon ^2 \gamma ^2 |{\mathcal{F}}|^2_{{\mathcal{A}}{\mathcal{A}}}, \end{equation}

where for simplicity we have used the fact that, to lowest order, ${\mathcal{F}}_T = -({1}/{2}){\mathcal{F}}_{\mathcal{A}}$ . Combining (4.21) with (4.25) therefore results in

(4.27) \begin{align} &{\nabla} ^2_{\boldsymbol{\alpha }} \psi = -2 \epsilon ^2 |{\mathcal{F}}|^2 e^{2\beta } + \epsilon ^2\gamma \bigg (2i (1 + \beta ) \big ({\mathcal{F}}_{\mathcal{A}} {\mathcal{F}}^* - {\mathcal{F}} {\mathcal{F}}_{\mathcal{A}}^*\big ) - i \big ({\mathcal{F}}_T{\mathcal{F}}^* - {\mathcal{F}} {\mathcal{F}}_T^*\big )\bigg )e^{2\beta } \nonumber \\ &\quad - 8 \epsilon ^4 |{\mathcal{F}}|^4 e^{4\beta } + \epsilon ^2 \gamma ^2 \left(\left(\frac {3}{4} + \frac {5}{2}\beta + \beta ^2\right)|{\mathcal{F}}|^2_{{\mathcal{A}}{\mathcal{A}}} e^{2\beta } - \big (4 + 10\beta + 4\beta ^2\big )|{\mathcal{F}}_{\mathcal{A}}|^2 e^{2\beta }\right)\!, \end{align}

which is just the linear Poisson equation. One can interpret (4.27) as equating the vorticity of the mean Lagrangian drift to the curl of the mean pseudomomentum, which acts here as a wave-induced source of vorticity (Salmon Reference Salmon2020). It is important to reiterate that, despite this interpretation, the vorticity of the fluid is still exactly zero everywhere in the fluid. The monochromatic mean Lagrangian drift shows how a mean flow that is sheared in the Lagrangian reference frame can still describe perfectly irrotational flow.

To close the system, we require boundary conditions on $\psi$ . The first one is simple: the flow vanishing at infinite depth implies $\psi \rightarrow 0$ as $\beta \rightarrow -\infty$ . The surface boundary condition is more subtle, as it was implicitly introduced in § 2 from the fact that the surface of the fluid is always defined by particles with $\beta =0$ . This is the Lagrangian equivalent of the kinematic boundary condition, (2.2) in the Eulerian frame, which in plain language states that particles which start at the surface always remain at the surface. From the perspective of the mean Lagrangian drift, there can therefore be no mass flux through the surface, making it a streamline which we can without loss of generality set to $\psi = 0$ at $\beta =0$ . This is not to say that there can be no mean vertical motion, only that any such motion at the surface must correspond with changes to the surface geometry, which is already set by the Lagrangian-mean water level. Slow changes in the Lagrangian-mean water level, found by taking a time derivative of (4.23), can be interpreted as a ‘vertical drift’ which several studies investigate (e.g. Vanneste & Young Reference Sulem and Sulem2022). We choose to separate these effects due to their distinct dynamical origins; the mean Lagrangian drift is fundamentally set by the vorticity (or lack thereof), whereas the Lagrangian-mean water level is constrained by the geometry of material curves and would be present even in the absence of the mean Lagrangian drift, such as in the rotational Gerstner (Reference Gerstner1802) wave.

The full solution to (4.27) can be found in Appendix A, although its general character is determined solely by considering the pseudomomentum forcing term. To leading order, this forcing is negative and concentrated near the surface. Assuming the wave packet has a finite width, so that the streamlines must be closed, a clockwise circulation will develop which moves with the packet. This circulation presents itself as a strong jet near the surface, as determined in (4.9) as the classical mean Lagrangian drift, but also includes a slow deep return flow in a direction opposite to that of wave propagation that is well known in the literature (Longuet-Higgins & Stewart Reference Longuet-Higgins and Stewart1962; McIntyre Reference McIntyre1981; Salmon Reference Salmon2020; Pizzo & Wagner Reference Pizzo and Melville2025). These studies, however, only constrain the lowest-order mean flow response.

At higher orders, additional forcing terms arise that modify the structure of the mean Lagrangian drift. The third-order forcing terms in (4.27) depend on the quantities $i({\mathcal{F}}_{\mathcal{A}}{\mathcal{F}}^* - {\mathcal{F}} {\mathcal{F}}_{\mathcal{A}}^*)$ and $i({\mathcal{F}_T}{\mathcal{F}}^* - {\mathcal{F}} {\mathcal{F}}_T^*)$ , which represent local fluctuations to the wavenumber and frequency from their characteristic values $k_0$ and $\omega _0$ , respectively. These terms therefore account for modifications to the mean Lagrangian drift associated with local variations in the phase speed.

At fourth order, there is a forcing term proportional to the fourth power of the envelope magnitude, which will always act to strengthen the near-surface mean Lagrangian drift, particularly during focusing when its magnitude is most pronounced. The quantity $|{\mathcal{F}}|^2_{{\mathcal{A}}{\mathcal{A}}}$ is the curvature of the squared envelope magnitude; it enhances the forcing where the curvature of the envelope is most negative (near the packet centre), and reduces it at the edges where the curvature is positive. Owing to the non-trivial vertical dependence, this effect is reversed at depth. The final contribution, proportional to $|{\mathcal{F}}_{\mathcal{A}}|^2$ , similarly enhances the forcing near the surface. Together, these fourth-order corrections provide enhancements to the forcing of the mean Lagrangian drift near the surface in regions where wave envelopes are both steep and concave, precisely the conditions present during focusing which led to the greatest observed transport.

Despite the complexity of the full solution, which can be found in Appendix A, the mean Lagrangian drift at the surface can be written explicitly as

(4.28) \begin{align} \langle x_\tau \rangle \big |_{\beta =0} =\epsilon ^2 |{\mathcal{F}}|^2 + \epsilon ^2\gamma \bigg (\frac {1}{2}i \big ({\mathcal{F}}_T {\mathcal{F}}^* - {\mathcal{F}} {\mathcal{F}}_T^*\big ) - \frac {1}{2} i\big ({\mathcal{F}}_{\mathcal{A}}{\mathcal{F}}^* - {\mathcal{F}} {\mathcal{F}}^*_{\mathcal{A}}\big ) +\frac {1}{2}\mathcal{H}\big (|{\mathcal{F}}|^2_{\mathcal{A}}\big )\bigg ) \nonumber \\ + 2\epsilon ^4 |{\mathcal{F}}|^4 + \epsilon ^2 \gamma ^2 \bigg ( \frac {1}{2}|{\mathcal{F}}_{\mathcal{A}}|^2 - \frac {1}{4}|{\mathcal{F}}|^2_{{\mathcal{A}}{\mathcal{A}}} + \frac {1}{4}i {\partial }_{\mathcal{A}}\mathcal{H}\big ({\mathcal{F}}_{T}{\mathcal{F}}^* - {\mathcal{F}}{\mathcal{F}}_{T}^*\big )\Big )\bigg ), \end{align}

where $\mathcal{H}$ represents the spatial Hilbert transform.

To check that our solutions reduce to known results for steep monochromatic waves, we see what happens when the envelope $\mathcal{F}$ is constant in space. From (4.16), we have

(4.29) \begin{equation} i {\mathcal{F}}_T = \frac {\epsilon ^2}{\gamma }\frac {1}{2}|{\mathcal{F}}|^2 {\mathcal{F}}, \end{equation}

which admits the exact solution (recalling $T = \gamma \tau$ )

(4.30) \begin{equation} {\mathcal{F}} = {\mathcal{F}}_0 e^{-\frac {1}{2}i \epsilon ^2 |{\mathcal{F}}_0|^2 \tau } , \end{equation}

where ${\mathcal{F}}_0$ is a complex constant, still $O(1)$ . This slow time modulation to $\mathcal{F}$ is just the classical Stokes correction to the phase speed for finite amplitude waves (Stokes Reference Skyner1847). Assuming without loss of generality that ${\mathcal{F}}_0 = 1$ (the physical dimensions can be added later), inserting (4.30) into the trajectories (4.13)–(4.15) yields

(4.31) \begin{align} x &= \alpha - \big(\epsilon + 2\epsilon ^3 e^{2\beta }\big)e^\beta \sin (\alpha - (c-U(\beta ))\tau ) + U(\beta )\tau , \\[-12pt]\nonumber \end{align}

(4.32) \begin{align} y &= \beta + \big(\epsilon + \epsilon ^3 e^{2\beta }\big) e^\beta \cos (\alpha - (c-U(\beta ))\tau ) + \frac {1}{2}\epsilon ^2 e^{2\beta }, \\[-12pt]\nonumber \end{align}

(4.33) \begin{align} p &= -\beta + \epsilon ^3\big(e^\beta - e^{3\beta }\big)\cos (\alpha - (c-U(\beta ))\tau ), \end{align}

where $c = ( 1 + ({1}/{2}) \epsilon ^2 )$ is the nonlinear, non-dimensionalised phase speed, and $U(\beta )$ is the mean Lagrangian drift, governed by

(4.34) \begin{equation} {\nabla} _{\boldsymbol{\alpha }}^2 \psi = -2\epsilon ^2 e^{2\beta } - 8 \epsilon ^4 e^{4\beta } - \epsilon ^4 e^{2\beta }, \end{equation}

which only depends on $\beta$ . Therefore, (4.34) reduces to an ordinary differential equation, and by simple integration the solution becomes

(4.35) \begin{equation} U(\beta ) = -\psi _\beta = \epsilon ^2 \underbrace {\Big (1 + \frac {1}{2}\epsilon ^2\Big )}_ce^{2\beta } + 2\epsilon ^4 e^{4\beta }. \end{equation}

While these are uniformly valid solutions that match those of Clamond (Reference Clamond2007), he defines his small steepness parameter $\epsilon '$ to be equal to $k_0 h/2$ , where $h$ is the crest to trough distance at the surface. The surface crest to trough distance from our solution (4.32) implies the relationship between these small parameters is

(4.36) \begin{equation} \epsilon ' = \epsilon + \epsilon ^3, \end{equation}

so that using the standard definition of wave steepness $\epsilon '$ , the drift at the surface can be expressed as

(4.37) \begin{equation} U(0) = \epsilon '^2 c', \quad c' = \left(1 + \frac {1}{2}\epsilon '^2 \right), \end{equation}

to fourth order in $\epsilon '$ , which matches the results in the literature (Longuet-Higgins Reference Longuet-Higgins1987). Because unsteady narrow-banded waves do not have an unambiguous geometric reference such as crest to trough height, we must instead settle for its more basic definition above based on the steepness of the first-order Lagrangian expansions.

4.2. Comparison with simulation

As a test of the above theory, we directly apply our above expression for the local mean Lagrangian drift of steep, narrow-banded waves (4.28) to estimate the surface transport of focusing wave packets from measurements of the wave envelope, comparing these predictions with the results from our fully nonlinear simulations. Note that the independent variables of the simulation, the initial positions $(x_0,y_0)$ , only coincide with the Lagrangian labels $(\alpha ,\beta )$ for particles which begin at rest (i.e. where $\mathcal{F}$ and its derivatives vanish). Since the Lagrangian displacement is only computed for particles which both begin and end at rest, for these particles we may treat both sets of independent variables identically. At the surface the transport is given by integrating (4.28) in time

(4.38) \begin{align} \delta x = \int \epsilon ^2 |{\mathcal{F}}|^2 + \epsilon ^2\gamma \bigg (\frac {1}{2}i \big ({\mathcal{F}}_T {\mathcal{F}}^* - {\mathcal{F}} {\mathcal{F}}_T^*\big ) - \frac {1}{2} i\big ({\mathcal{F}}_{\mathcal{A}}{\mathcal{F}}^* - {\mathcal{F}} {\mathcal{F}}^*_{\mathcal{A}}\big ) +\frac {1}{2}\mathcal{H}\big (|{\mathcal{F}}|^2_{\mathcal{A}}\big )\bigg ) \nonumber \\ + 2\epsilon ^4 |{\mathcal{F}}|^4 + \epsilon ^2 \gamma ^2 \bigg ( \frac {1}{2}|{\mathcal{F}}_{\mathcal{A}}|^2 - \frac {1}{4}|{\mathcal{F}}|^2_{{\mathcal{A}}{\mathcal{A}}} + \frac {1}{4}i {\partial }_{\mathcal{A}}\mathcal{H}\big ({\mathcal{F}}_{T}{\mathcal{F}}^* - {\mathcal{F}}{\mathcal{F}}_{T}^*\big )\bigg ) \, {\textrm {d}} \tau . \end{align}

To compare this prediction with the simulations, $\mathcal{F}$ is estimated by using the spatial Hilbert transform of the vertical Lagrangian positions at each time to produce an analytic envelope signal. While this correctly estimates the magnitude of the envelope $\mathcal{F}$ , its phase must be corrected by removing the phase of the carrier wave, which is given by (4.12). Because the theory outlined above is non-dimensional, all quantities must be dimensionalised by a characteristic wavenumber $k_0$ and frequency $\omega _0 = \sqrt {g k_0}$ . Note that $k_0$ and $\omega _0$ need not be equivalent to the initially specified $k_c$ and $\omega _c$ . Because the phase of $\mathcal{F}$ represents narrow-banded deviations from the phase of the carrier wave, we choose $k_0$ such that these deviations have zero mean when integrated in time. Once the envelope $\mathcal{F}$ is computed, all that remains is to compute the surface mean Lagrangian transport (4.38) (where $\epsilon$ and $\gamma$ are implicitly included in measurements of $\mathcal{F}$ and its slow derivatives) for each particle and integrate in time over the duration of the packet passing to estimate the total transport predicted by this theory. From this, we can compute the maximum and mean transport predicted by narrow-banded theory to compare directly with the simulations.

Figure 5.

The mean surface transport $\langle \delta x \rangle$ within the focusing region as a function of $S$ computed both directly from simulation (circles) and using our higher-order theory (4.38) (lines) for each simulation. The mean surface transport for all laboratory experiments (Lenain et al. Reference Lenain, Pizzo and Melville2019; Sinnis et al. Reference Seitz, Freilich and Pizzo2021) is additionally shown (triangles). In $(a)$ , $\langle \delta x \rangle$ is normalised by the central wavenumber $k_c$ , and each line represents the theoretical prediction of mean transport for each bandwidth value. In $(b)$ , $\langle \delta x \rangle$ is also normalised by the linear bandwidth dependence $f(\Delta )$ (3.7) which collapses the results. The theory performs best at lower values of $\varDelta$ where the narrow-banded envelope assumption is most valid.

Figure 5 presents the mean surface transport obtained directly from the simulations, together with the corresponding theoretical prediction evaluated for each case using (4.38). Also included is the mean surface transport from all available laboratory experiments (Lenain et al. Reference Lenain, Pizzo and Melville2019; Sinnis et al. Reference Seitz, Freilich and Pizzo2021). In figure 5(a), the mean transport is scaled only by the central wavenumber $k_c$ , and good agreement between simulation and theory is found, especially for lower values of $\varDelta$ and $S$ . This is expected since (4.38) is valid to fourth order in the small parameters $(\epsilon ,\gamma )$ which are proportional to $(S,\Delta )$ respectively. The theory still performs well across a wide range of parameter space, indicating that the narrow-banded wave approximation is able to capture the local enhancements to the mean Lagrangian drift. Figure 5(b) collapses these results, normalising by both $k_c$ and the linear bandwidth dependence $f(\Delta )$ , showing that the bandwidth dependence of the surface transport, even in steep packets, is still reasonably approximated by linear theory. This collapse demonstrates good qualitative agreement between the simulations, theory and experiments, even though neither the theory nor the simulations account for finite-depth effects.