2.1 Linear monochromatic waves: the classical Stokes drift

For completeness, we begin by considering a deep-water monochromatic periodic surface gravity wave of constant amplitude $a$ , with angular frequency $\unicode[STIX]{x1D714}$ and wavenumber given by the linear dispersion relationship $\unicode[STIX]{x1D714}^{2}=gk$ where $g$ is the acceleration due to gravity. To first order then,

In the Lagrangian reference frame, we label a particle $\boldsymbol{x}=(\unicode[STIX]{x1D709},\unicode[STIX]{x1D701})$ . The Eulerian and Lagrangian reference frames are then connected by

Now, for a particle originally at $\boldsymbol{x}_{0}=(\unicode[STIX]{x1D709}_{0},\unicode[STIX]{x1D701}_{0})$ this implies

The Lagrangian velocity at a point is found by Taylor expanding the Eulerian velocity about $(\unicode[STIX]{x1D709}_{0},\unicode[STIX]{x1D701}_{0})$ , that is

The averaged Lagrangian velocity, or drift, is then given by

where $T=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}$ and the phase velocity is given by $c=g/\unicode[STIX]{x1D714}$ . This is the classical Stokes drift for periodic linear deep-water surface gravity waves, and by definition $u_{s}\equiv u_{L}$ .

2.2 Linear wave packets

In the context of the ocean, a slightly more realistic example comes from considering a compact wave packet. This follows the work of Kenyon (Reference Kenyon1969). Let

where $\hat{a}(k)$ is the amplitude density and the reality condition implies

Furthermore, setting $\unicode[STIX]{x1D70C}=1$ , we define the energy spectral density $F(k)$ by

where the integral is over the (compact) wave packet and the last relation is via Parseval’s theorem.

The Lagrangian motion induced by this packet, found via equation (2.6), is

and

The Stokes drift then follows from (2.9), from which we find

In the limit $F=ga_{0}^{2}\unicode[STIX]{x1D6FF}(k)$ , we return equation (2.10). At the surface, the drift reads

This prediction is linear and does not give information about whether focusing or defocusing packets lead to greater convergence. To gain intuition about this phenomenon, in appendix B, we examine the drift associated with weakly nonlinear wave packets and find that focusing induces a positive correction to the Stokes drift while defocusing induces a negative correction. Furthermore, steep and breaking waves are not linear, so this analysis does not yield information about what happens in those cases. To gain intuition about this behaviour, we next turn to a scaling argument.

2.3 Scaling argument for the added drift due to breaking

In a recent paper, Pizzo (Reference Pizzo2017) examined the equation of John (Reference John1953), which is a relationship between the horizontal kinematics of a particle on the fluid surface and the free surface geometry, to derive a criterion for particles to surf an underlying wave. There, it is shown that in a localized region on the forward face of the wave, near the crest, particles accelerate and surf the underlying wave. The John equation (see also Sclavounos (Reference Sclavounos2005)) is

where again $\unicode[STIX]{x1D709}$ is a the horizontal position of the particle.

Now, following Drazen et al. (Reference Drazen, Melville and Lenain2008), Loewen & Melville (Reference Loewen and Melville1991), we assume the free surface displacement $\unicode[STIX]{x1D702}$ scales with $h$ while horizontal distances scale with the typical wavenumber $k$ and time with the typical angular frequency $\unicode[STIX]{x1D714}^{-1}$ . Furthermore, we assume that the particle velocity increases with increasing slope, consistent with numerical (Vinje & Brevig Reference Vinje and Brevig1981) and laboratory (Drazen et al. Reference Drazen, Melville and Lenain2008) observations. Therefore, the maximum acceleration, which we denote as $A_{m}$ , is given by

where $\unicode[STIX]{x1D712}$ is a $O(1)$ non-dimensional constant. Assuming this acceleration holds approximately over the duration of breaking, $T$ , which we take to be $2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}$ (Rapp & Melville Reference Rapp and Melville1990), we conclude

where $\unicode[STIX]{x1D712}^{\prime }$ is an order $O(10)$ constant and $\unicode[STIX]{x0394}x$ is the total horizontal displacement. We define the slope $S=hk$ , and following earlier work introduce a threshold for breaking $S_{0}$ , leading to

where $\unicode[STIX]{x1D712}^{\prime \prime }$ is an order $O(10)$ constant. This can be written as a drift velocity due to breaking $u_{L,s}=\unicode[STIX]{x0394}x/T$ , with $c=\unicode[STIX]{x1D714}/k$ the characteristic phase velocity, that is

As we will show later in the paper, this result is observed in our numerical results. Ultimately, this scaling should be tested in laboratory experiments. We now investigate the Lagrangian drift for non-breaking and breaking wave focusing packet by performing direct numerical simulations.

3.1 The Gerris flow solver

We solve the two-dimensional two-phase incompressible Navier–Stokes equations accounting for surface tension and viscous effects using the open-source solver Gerris (Popinet Reference Popinet2003, Reference Popinet2009), based on a quad/octree adaptive spatial discretization, multilevel Poisson solver. The interface between the high density liquid (water) and the low density gas (air) is reconstructed by a geometric volume of fluid method.

The multifluid interface is traced by a function ${\mathcal{T}}(\boldsymbol{x},t)$ , defined as the volume fraction of a given fluid in each cell of the computational mesh. The density and viscosity can thus be written as $\unicode[STIX]{x1D70C}({\mathcal{T}})=T\unicode[STIX]{x1D70C}_{w}+(1-{\mathcal{T}})\unicode[STIX]{x1D70C}_{a}$ , $\unicode[STIX]{x1D707}({\mathcal{T}})={\mathcal{T}}\unicode[STIX]{x1D707}_{w}+(1-{\mathcal{T}})\unicode[STIX]{x1D707}_{a}$ , with $\unicode[STIX]{x1D70C}_{w}$ , $\unicode[STIX]{x1D70C}_{a}$ and $\unicode[STIX]{x1D707}_{w}$ , $\unicode[STIX]{x1D707}_{a}$ the density and viscosity of the two fluids (water and air), respectively. The incompressible, variable density, Navier–Stokes equations with surface tension can be written as

with $\boldsymbol{u}=(u,w)$ the fluid velocity, $\unicode[STIX]{x1D70C}\equiv \unicode[STIX]{x1D70C}(\boldsymbol{x},t)$ the fluid density, $\unicode[STIX]{x1D707}\equiv \unicode[STIX]{x1D707}(\boldsymbol{x},t)$ the dynamic viscosity and $\unicode[STIX]{x1D63F}$ the deformation tensor defined as $\unicode[STIX]{x1D60B}_{ij}\equiv (\unicode[STIX]{x2202}_{i}u_{j}+\unicode[STIX]{x2202}_{j}u_{i})/2$ . The Dirac delta, $\unicode[STIX]{x1D6FF}_{s}$ applied on the surface only, expresses the fact that the surface tension term is concentrated on the interface, where $\unicode[STIX]{x1D6FE}$ is the surface tension coefficient, $\unicode[STIX]{x1D705}$ and $\boldsymbol{n}$ the curvature and normal to the interface, respectively.

The density and viscosity ratios of the two phases are those of air and water. The surface tension is $\unicode[STIX]{x1D6FE}=0.07~\text{N}~\text{m}^{-1}$ , corresponding to the surface tension between air and water. The Reynolds number in the liquid is defined by $Re=c\unicode[STIX]{x1D706}/\unicode[STIX]{x1D708}$ , with $c$ the linear deep-water gravity wave central phase speed of the wave packet (see below) and $\unicode[STIX]{x1D708}$ the kinematic viscosity of the liquid (water). Due to computational limitations, related to spatial resolution constraints, we limit ourselves to $Re=40\,000$ . This resolves the viscous boundary layer at the surface, and is still a large Reynolds number. However, it is smaller than the $O(10^{6})$ Reynolds number corresponding to waves generated in the laboratory. Nevertheless, we have shown in previous work on energy dissipation, air entrainment and the breaking-induced currents that this lower Reynolds number still leads to consistent results with laboratory studies, implying that the Reynolds number is sufficient to describe the phenomena in question (Deike et al. Reference Deike, Popinet and Melville2015, Reference Deike, Melville and Popinet2016; Pizzo et al. Reference Pizzo, Deike and Melville2016). This gives us confidence in studying the Lagrangian transport in the breaking-induced flow in this paper.

Note that we are performing direct numerical simulations (DNS) of the Navier–Stokes equations to solve for breaking waves, therefore no model or assumption is used in our computations. There is simply a limitation in the Reynolds number we can achieve with the molecular viscosity. This DNS approach should not be confused with other approaches to model breaking wave dissipation using a turbulent model, such as an eddy viscosity, which rely on laboratory experiment to calibrate the value of the eddy viscosity. Tian et al. (Reference Tian, Perlin and Choi2010), Tian, Perlin & Choi (Reference Tian, Perlin and Choi2012) used an eddy viscosity between $0.5\times 10^{-3}$ and $1.7\times 10^{-3}~\text{m}^{2}~\text{s}^{-1}$ to reproduce the dissipation rate observed in their breaking wave experiments. This would lead to a Reynolds number between $2\times 10^{3}$ and $8\times 10^{3}$ , between five to twenty times smaller than the Reynolds number we consider using the molecular diffusivity.

The present solver has been successfully used in multiphase problems like atomization (Fuster et al. Reference Fuster, Agbaglah, Josserand, Popinet and Zaleski2009; Tomar et al. Reference Tomar, Fuster, Zaleski and Popinet2010), the growth of instabilities at the interface (Fuster et al. Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013) and the nonlinear dynamics of wave breaking in two (Deike et al. Reference Deike, Popinet and Melville2015; Pizzo et al. Reference Pizzo, Deike and Melville2016) and three dimensions (Deike et al. Reference Deike, Melville and Popinet2016). These previous studies have validated these methods and shown the ability of the numerical simulations to reproduce properties observed in the laboratory regarding capillary effects on wave breaking (Deike et al. Reference Deike, Popinet and Melville2015), wave dissipation due to breaking (Deike et al. Reference Deike, Popinet and Melville2015, Reference Deike, Melville and Popinet2016), air entrainment and bubble statistics in breaking waves (Deike et al. Reference Deike, Melville and Popinet2016), as well as current generation (Pizzo et al. Reference Pizzo, Deike and Melville2016).

3.2 The focusing wave packet configuration

The numerical wave tank is $L=32~\text{m}$ long and $l=2~\text{m}$ high, with a water depth of $d=1~\text{m}$ . Top and bottom boundary conditions are no slip. The left-hand side of the tank corresponds to the numerical wave maker and consists of a forcing boundary condition on the velocity $(u(x=0,t),w(x=0,t))$ and on the interface $\unicode[STIX]{x1D702}(x=0,t)$ . Following laboratory experiments, we consider a focusing wave packet to create breaking waves (Rapp & Melville Reference Rapp and Melville1990; Drazen et al. Reference Drazen, Melville and Lenain2008) and the imposed boundary conditions at $x=0$ are

The linear prediction of the focusing time and location are $t_{b}=25~\text{s}$ and $x_{b}=12~\text{m}$ , in the middle of the tank. Following Rapp & Melville (Reference Rapp and Melville1990), Drazen et al. (Reference Drazen, Melville and Lenain2008), we consider packets with $N=32$ components, a central frequency $f_{c}=0.89~\text{Hz}$ (corresponding to a central wavelength $\unicode[STIX]{x1D706}_{c}=2\unicode[STIX]{x03C0}/k_{c}=1.9~\text{m}$ which is found using the linear dispersion relation for gravity waves, $\unicode[STIX]{x1D714}^{2}=gk\tanh (kd)$ ) and a bandwidth $\unicode[STIX]{x1D6E5}=\unicode[STIX]{x0394}f/f_{c}=0.75$ . We use the constant slope formulation (Drazen et al. Reference Drazen, Melville and Lenain2008) with $s=a_{i}k_{i}$ constant for all components. The maximum linear slope at focusing is then given by

The central phase speed of the wave packet is $c=\unicode[STIX]{x1D714}_{c}/k_{c}$ . For a non-breaking packet, the theoretical Stokes drift of the packet can be computed combining equation (2.16) and the initial shape of the packet given by (3.2).

The last 12 m of the tank is a numerical sponge layer, based on the one developed by Clement (Reference Clement1995) for a potential solver. This consists of an added empirical friction term which acts to attenuate the velocity. This leads to almost complete suppression of waves propagating towards the right end of the numerical tank. The transition between the physical part of the tank (the first 24 m) and the sponge layer (the last 12 m) is smooth to avoid wave reflections. Finally the end of the tank on the right corresponds to a wave radiation boundary condition, in order to remove any remaining small perturbations.

Adaptive mesh refinement is used to accurately solve for the interface and the vortical structures. We tested several resolutions and adaptive grid configurations and report in this paper on high and medium resolution configurations. The equivalent grid resolution for the high resolution case is a mesh size of $\unicode[STIX]{x1D6FF}_{h,i}=1~\text{mm}$ on the interface gradients and $\unicode[STIX]{x1D6FF}_{h,v}=4~\text{mm}$ on the vortical structure in the air or water, and the medium resolution has a mesh of $\unicode[STIX]{x1D6FF}_{m,i}=2~\text{mm}$ on the interface gradients and $\unicode[STIX]{x1D6FF}_{m,v}=8$ mm on the vortical structure. As will be discussed below, we observed convergence of the results for the total dissipation and total drift between these two resolutions. Note, although we are solving for some effects due to surface tension (the capillary length scale for air–water is $l_{gc}=\sqrt{\unicode[STIX]{x1D6FE}/(\unicode[STIX]{x1D70C}g)}\approx 3~\text{mm}$ , and the corresponding wavelength is $\unicode[STIX]{x1D706}_{gc}=2\unicode[STIX]{x03C0}l_{gc}\approx 1~\text{cm}$ ), we do not fully resolve the dynamics of micro-breakers and parasitic capillaries, for which a higher resolution is needed (around 0.2 mm was successfully used for parasitic capillaries in Deike et al. (Reference Deike, Popinet and Melville2015) for example).

The flow is seeded by Lagrangian particles that are systematically tracked. Each Lagrangian particle $i$ is described by its Lagrangian trajectory $\{x_{i}(t),z_{i}(t)\}$ , or velocity $\{u_{i}(t),w_{i}(t)\}$ . The initial location of the particles is denoted by $\{x_{i}(t=0),z_{i}(t=0)\}=\{x_{i,0},z_{i,0}\}$ . The flow is seeded with one particle every 10 cm horizontally, i.e. $x_{i+1,0}-x_{i,0}=10~\text{cm}$ ; and 1 particle every 1 cm vertically, i.e. $z_{i+1,0}-z_{i,0}=1~\text{cm}$ , leading to a total of over 30 000 Lagrangian particles over the tank. We mainly focus our analysis on the Lagrangian kinematics in an area around the focusing point of the wave packet. The accuracy of the tracking scheme has been tested in Tomar et al. (Reference Tomar, Fuster, Zaleski and Popinet2010) and has been verified in our case for propagating waves by comparing the Lagrangian drift for a monochromatic small amplitude wave against the theoretical Stokes drift. The test case is presented in the appendix A and we found very good agreement between the simulation and the theory.

3.3 Wave focusing, breaking and energy dissipation due to breaking

Figure 1(a) shows the space–time evolution of the free surface displacement $\unicode[STIX]{x1D702}(x,t)$ for a typical numerical experiment of a focusing wave packet leading to breaking. The packet is generated at the left side of the tank and focuses around $x\approx 12~\text{m}$ and $t\approx 25~\text{s}$ , as expected by linear theory and in agreement with laboratory experiments (Rapp & Melville Reference Rapp and Melville1990; Drazen et al. Reference Drazen, Melville and Lenain2008). The wave packet propagates at the weighted group velocity as observed experimentally (Drazen et al. Reference Drazen, Melville and Lenain2008; Tian et al. Reference Tian, Perlin and Choi2010, see also Pizzo Reference Pizzo2015, Pizzo & Melville Reference Pizzo and Melville2016), $c_{gs}=(\sum _{i=1}^{N}a_{i}^{2}c_{gi})/(\sum _{i=1}^{N}a_{i}^{2})$ , where $c_{gi}=\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{i}/\unicode[STIX]{x2202}k_{i}$ is the group velocity of the wave component. The wave breaks close to the focusing point, with some waves being radiated back towards the left boundary of the tank. The packet then continues to propagate and defocuses, before being damped when it reaches the numerical beach ( $x>24~\text{m}$ ), and no wave signal is visible for $x>30~\text{m}$ . No reflection is visible on the right-hand side of the tank.

Figure 1. (a) Space–time diagram of the wave height $\unicode[STIX]{x1D702}(x,t)$ , showing the propagation and focusing of the wave packet, for $S=0.38$ , $f_{c}=0.89~\text{Hz}$ , $\unicode[STIX]{x0394}f/f_{c}=0.75$ . Breaking occurs at $t\approx 25~\text{s}$ and at $x\approx 12~\text{m}$ . Colour scale is $\unicode[STIX]{x1D702}(x,t)$ in metres. The dashed line is the energy weighted group velocity $c_{gs}$ defined by Drazen et al. (Reference Drazen, Melville and Lenain2008), that describes the propagation of the wave packet. No significant reflection can be seen at the right end of the numerical tank due to the numerical sponge layer. (b) Breaking parameter $b$ as a function of wave slope $S$ , for the high resolution (full circle) and medium resolution (cross) runs. The present simulations are in agreement with earlier experimental and numerical work as well as the inertial scaling argument from Drazen et al. (Reference Drazen, Melville and Lenain2008), Romero et al. (Reference Romero, Melville and Kleiss2012). Laboratory data are from Melville (Reference Melville1994) (M), Banner & Peirson (Reference Banner and Peirson2007) (BP), Drazen et al. (Reference Drazen, Melville and Lenain2008) (DML), Grare et al. (Reference Grare, Peirson, Branger, Walker, Giovanangeli and Makin2013) and numerical data from Deike et al. (Reference Deike, Popinet and Melville2015, Reference Deike, Melville and Popinet2016).

Note, convergence of the dynamical and integrated quantities such as the time and location of focusing, the total wave dissipation, and average drift of patches of particles have been verified between the medium and high resolutions.

Figure 1(b) shows the breaking parameter $b$ as a function of the wave slope $S$ . The breaking parameter $b$ is related to the dissipation rate per unit length of breaking crest $\unicode[STIX]{x1D716}_{l}$ by (Duncan Reference Duncan1981; Phillips Reference Phillips1985)

and is computed from the dissipation due to breaking as in Drazen et al. (Reference Drazen, Melville and Lenain2008), by computing the wave energy budget before and after breaking. Excellent agreement is found between our numerical results in the present configuration and earlier laboratory work (Drazen et al. Reference Drazen, Melville and Lenain2008; Romero et al. Reference Romero, Melville and Kleiss2012; Grare et al. Reference Grare, Peirson, Branger, Walker, Giovanangeli and Makin2013), the inertial argument from Drazen et al. (Reference Drazen, Melville and Lenain2008) and recent numerical simulations (Deike et al. Reference Deike, Popinet and Melville2015, Reference Deike, Melville and Popinet2016; Derakhti & Kirby Reference Derakhti and Kirby2016; Pizzo et al. Reference Pizzo, Deike and Melville2016). The breaking threshold observed in our present configuration is $S_{\ast }=0.31$ , where breaking is defined by the appearance of a vertical interface, which is a classic geometric criterion (Perlin et al. Reference Perlin, Choi and Tian2013; Deike et al. Reference Deike, Popinet and Melville2015). This relatively high value is probably related to the limited mesh size which does not allow the proper resolution for parasitic capillaries and micro-breaking processes.

The mass in each fluids (air and water) is conserved to better than 0.2 % as shown in appendix A.

The results in the remainder of the paper are normalized by the characteristic wavenumber $k_{c}$ , angular frequency $\unicode[STIX]{x1D714}_{c}$ and phase velocity $c$ , all taken as the central component of the wave packets. For dynamical quantities, such as the group velocity or the wave energy, the weighted quantities (e.g. the weighted group velocity) arise naturally (Pizzo Reference Pizzo2015) and provide a better non-dimensionsalization of the variables (Drazen et al. Reference Drazen, Melville and Lenain2008; Tian et al. Reference Tian, Perlin and Choi2010). However, for the present study, the kinematics of the wave is considered and the weighted quantities do not arise naturally in the theoretical analysis which justifies our choice to work with the central wave packet quantities. Note that all results on the drift presented in § 4 are not modified if we were to use weighted quantities, as the characteristic scales are only used for normalization.

4.1 Particle trajectories

Figures 2 and 3 show the spatial wave surface $\unicode[STIX]{x1D702}(x,t)$ at different times $t$ before and after breaking. The propagation and focusing of the wave packet close to the focusing point is visible, in time and space, for a non-breaking case (figure 2, $S=0.26$ ) and a breaking case (figure 3, slope $S=0.38$ ). These results are representative of our breaking and non-breaking cases. The horizontal axis $x$ is centred by the linear focusing location $x_{b}$ and normalized by the central wavenumber, $(x-x_{b})k_{c}$ , the vertical axis is normalized to $zk_{c}$ , while the time is centred by the linear focusing time $t_{b}$ and normalized by the central angular frequency $\unicode[STIX]{x1D714}_{c}$ to be $(t-t_{b})\unicode[STIX]{x1D714}_{c}$ .

Focusing occurs within one wave period $1/f_{c}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{c}$ of the linear focusing time $t_{b}$ , and close to the linear focusing location $x_{b}$ , and is consistent with observations in laboratory experiments (Rapp & Melville Reference Rapp and Melville1990; Drazen et al. Reference Drazen, Melville and Lenain2008). Together with the free surface, these figures also show the displacement of particles $\{x_{i}(t),z_{i}(t)\}$ , at different initial depths ( $z_{i,0}=0k_{c},-0.64k_{c}$ and $-1.28k_{c}$ ), and initial horizontal locations upstream ( $(x_{i,0}-x_{b})k_{c}=-2.9k_{c}$ ), near ( $(x_{i,0}-x_{b})k_{c}=1.3k_{c}$ ) and downstream ( $(x_{i,0}-x_{b})k_{c}=8k_{c}$ ) of the focusing point.

Figure 2. Interface evolution $\unicode[STIX]{x1D702}(x,t)$ and particle trajectories during wave propagation for a non-breaking focusing packet, $f_{c}=0.89~\text{Hz}$ ( $k_{c}=5.58~\text{rad}~\text{m}^{-1}$ ),  $\unicode[STIX]{x0394}f/f_{c}=0.76$ and maximum linear slope $S=0.26$ . The interface $\unicode[STIX]{x1D702}(x,t)$ is shown in black, while the particle trajectories are colour coded by their horizontal velocity $u/c$ . Times between $(t-t_{b})\unicode[STIX]{x1D714}_{c}=-5.59$ before focusing (a) and $(t-t_{b})\unicode[STIX]{x1D714}_{c}=16.8$ after focusing (f) are shown. Focusing occurs at $t\approx t_{b}$ , with particles presenting a larger drift and velocity at that time. The drift of all particles is similar for all particles at the same depth and decreases with depth. Note that particle trajectories shown here above the interface are in fact particles at the interface but for earlier times.

Figure 3. Interface evolution and particle trajectories during the wave propagation and wave breaking event for a focusing packet, $f_{c}=0.89$ Hz ( $k_{c}=5.58~\text{rad}~\text{m}^{-1}$ ),  $\unicode[STIX]{x0394}f/f_{c}=0.76$ and maximum linear slope $S=0.38$ . The interface $\unicode[STIX]{x1D702}(x,t)$ is in black, while the particle trajectories are colour coded by their horizontal velocity $u/c$ . Times between $(t-t_{b})\unicode[STIX]{x1D714}_{c}=-5.59$ before breaking (a) and $(t-t_{b})\unicode[STIX]{x1D714}_{c}=16.8$ after breaking (f) are shown. Breaking occurs at $t\approx t_{b}$ , with particles presenting a much larger drift and velocity. The drift decreases with depth. Note that particle trajectories shown here above the interface are in fact particles at the interface but for earlier times.

Figure 2 shows a typical non-breaking focusing case (central frequency $f_{c}=0.89~\text{Hz}$ , bandwidth $\unicode[STIX]{x0394}f/f_{c}=0.76$ and maximum linear slope $S=0.26$ ). The wave reaches its largest slope near the focusing point. The trajectories of the particles at the surface are similar in all non-breaking cases and can be described as follows: small deviations to almost closed ellipses before the focusing time and the arrival of the main part of the wave packet. These small amplitude oscillations correspond to waves of small amplitude arriving before focusing. During focusing, a significant horizontal jump, or drift in the horizontal $x$ -direction is observed, and corresponds to the passing of the wave group. The total horizontal displacement during the whole simulation is almost entirely due to this event and it corresponds to the horizontal transport, similar to the Stokes drift described for a single frequency propagating wave. The horizontal velocity of the particles $u$ at the surface at the focusing time are $u/c\approx 0.2$ , where $c$ is the central phase speed of the wave packet. After focusing, small oscillations are visible again, with much smaller velocities. Next, the depth dependence is characterized by the same behaviour, with an amplitude that is attenuated by an exponential decay with e-folding scale corresponding to twice the wavenumber. Finally the particles close to the bottom of the tank experience a negative drift (a displacement opposite to the wave propagation), corresponding to the return flow, consistent with mass conservation (Longuet-Higgins Reference Longuet-Higgins1953).

Figure 3 shows a typical focusing packet leading to breaking (central frequency $f_{c}=0.89~\text{Hz}$ , bandwidth $\unicode[STIX]{x0394}f/f_{c}=0.76$ and maximum linear slope $S=0.38$ ). The focusing of wave energy leads to breaking and the impact of the jet occurs at $t\approx t_{b}$ . The displacement of the particles at the surface that are located far enough downstream of the focusing breaking point is similar to that discussed above for a non-breaking focusing packet: small oscillations before the focusing event, with small velocities ( $u/c<0.1$ ) and a small jump when the main packet passes with a velocity $u/c\approx 0.4$ . The net drift is therefore larger in the horizontal direction than for the non-breaking case, even out of the breaking point, due to the increase in the wave amplitude (or slope). A similar picture is visible for particles far enough upstream of the breaking point. Again, almost all the drift of the particles happens during a short period around focusing (typically one central wave period $1/f_{c}$ ); the displacement before and after being small compared to the jump when the focused wave passes by. The amplitude of the displacement decreases with depth for the particles away from the breaking area in a similar way to what we described for the non-breaking packet.

By comparison, the particles located close to the breaking point experience a large jump, or drift in the propagation direction during the breaking event, with a displacement that can be larger than one meter, i.e. of the order of $\unicode[STIX]{x1D706}_{c}$ . The velocity of the particles at breaking can be above the central phase velocity $c$ (Vinje & Brevig Reference Vinje and Brevig1981; Pizzo Reference Pizzo2017), with a value of $1.1c$ when the particle catches the breaker, and up to $1.5c$ during splashes for the particular example displayed here. This displacement is therefore an order of magnitude larger than for non-breaking particles at the focusing point, or for particles away from the breaking point. This corresponds to a Lagrangian transport due to the breaking wave which is the central point of this paper. Note also, particles can be transported vertically (downward) in the water when located around the breaking area corresponding to the vertical transport of horizontal momentum (Rapp & Melville Reference Rapp and Melville1990; Drazen & Melville Reference Drazen and Melville2009).

Particles located deeper in the water experience a smaller displacement in the propagation direction when located in the upper half of the tank and a negative displacement in the lower half of the tank, corresponding to the return flow, again consistent with mass conservation.

4.2 Drift at the surface

We now discuss the displacement of particles at the surface. Figure 4 shows the trajectories of such particles located just before, at, and just after the focusing point for a non-breaking and a breaking focusing packet, with the time colour coded.

Figure 4. Trajectories of particles for focusing packets, close to focusing (a) no breaking and close to breaking (b) breaking. $f_{c}=0.89~\text{Hz}$ , bandwidth $\unicode[STIX]{x0394}f/f_{c}=0.76$ . $S=0.26$ for the non-breaking packet and $S=0.38$ for the breaking packet. Particles just before the break point do not catch the wave, and experience a net drift smaller than the ones that are catching the breaking wave. Finally, particles further upstream or downstream of the breaking event are not significantly influenced by breaking.

As already described above, in the non-breaking case, all trajectories are similar, with small oscillations before the focusing time, and a small jump at focusing in the propagation direction. The net drift in the propagation direction is between $0.05\unicode[STIX]{x1D706}_{c}$ and $0.1\unicode[STIX]{x1D706}_{c}$ , for all particles at the surface. This observation applies over a wide area around the focusing point, where the horizontal extent of the initial particle positions is given by $-2<(x_{0}-x_{b})k_{c}<10$ , i.e. over approximately $2\unicode[STIX]{x1D706}_{c}$ . The particles initially near the focus point experience a slightly larger drift due to focusing which is explained by the weakly nonlinear theory presented in appendix B, where the focusing effect induces a positive correction to the classical Stokes drift. This effect increases with the wave slope $S$ , as suggested theoretically.

In the breaking case, particles initially in a narrow window around the effective breaking point (in this case close to $(x-x_{b})k_{c}\approx 1$ ) experience a much larger displacement (drift) in the propagation direction, than those located further upstream or downstream. The particle located at $(x_{0}-x_{b})k_{c}=-1.9$ moves only $0.06\unicode[STIX]{x1D706}_{c}$ , while that at $(x_{0}-x_{b})k_{c}=-0.3$ moves around $0.2\unicode[STIX]{x1D706}_{c}$ . The one at $(x_{0}-x_{b})k_{c}=1.3$ is caught by the breaking wave and moves around $0.36\unicode[STIX]{x1D706}_{c}$ , and finally the one slightly further at $(x_{0}-x_{b})k_{c}=1.6$ moves almost $1\unicode[STIX]{x1D706}_{c}$ , surfing the wave. Particles much further downstream have a smaller drift. This breaking case shows differences of an order of magnitude for the drift of particles located in a window of length less than one $\unicode[STIX]{x1D706}_{c}$ ; illustrating the importance of the local effect of breaking on the drift. This effect has been recently discussed theoretically by Pizzo (Reference Pizzo2017). The drift due to breaking is found to be an order of magnitude larger than the drift for a non-breaking packet, consistent with the scaling argument presented in § 2 for the drift due to the breaking event.

As discussed previously, most of the displacement occurs within a space–time window related to the passing of the focused wave, and scaled on the characteristic wavelength and period of the wave packet. This is even more dramatic in the breaking case, where the acceleration is very strong during breaking.

4.3 Spatial dependence of the total drift

We now focus on the total horizontal drift of the particles. We define the total drift as $\unicode[STIX]{x0394}x(x_{0},z_{0})=(x_{i}(t_{f})-x_{i}(t_{0}))$ , where $t_{f}$ is a time long after breaking, $(t_{f}-t_{b})\unicode[STIX]{x1D714}_{c}\gg 1$ (typically the final time of the simulation), and $t_{0}$ is the initial time, long before breaking, $(t_{0}-t_{b})\unicode[STIX]{x1D714}_{c}\ll -1$ (typically the initial time of the simulation). The total drift is not sensitive to the choice of $t_{f}$ , $t_{0}$ , as long as the conditions $(t_{0}-t_{b})\unicode[STIX]{x1D714}_{c}\ll -1$ and $(t_{f}-t_{b})\unicode[STIX]{x1D714}_{c}\gg 1$ are fulfilled. For example, the drift obtained for $(t_{f}-t_{b})\unicode[STIX]{x1D714}_{c}=56$ and $(t_{f}-t_{b})\unicode[STIX]{x1D714}_{c}=27$ (with $(t_{0}-t_{b})\unicode[STIX]{x1D714}_{c}=-140$ ) varies by less than 5 %.

Figure 5. (a) Total horizontal drift for particles initially at the surface, $z_{0}=0$ , normalized by the central wavelength $\unicode[STIX]{x1D706}_{c}$ , for a breaking case ( $S=0.38$ ) as a function of the normalized initial horizontal position of the particles $(x_{0}-x_{b})k_{c}$ . Triangles are the raw data $\unicode[STIX]{x0394}x(z_{0})$ , while the solid line, $\overline{\unicode[STIX]{x0394}x(z_{0})}$ , is obtained after a running mean averaged over $\unicode[STIX]{x1D706}_{c}$ . The breaking or focusing point $x_{b,n}$ is taken as the maximum drift location for each case, and the horizontal dashed line indicates the breaking area, $[x_{b,n}-\unicode[STIX]{x1D706}_{c}/2:x_{b,n}+\unicode[STIX]{x1D706}_{c}/2]$ . The area of breaking is where particles experience a stronger drift due to breaking and will be used to compute vertical profiles and drift at the surface. (b)  $\overline{\unicode[STIX]{x0394}x(z_{0})}/\unicode[STIX]{x1D706}_{c}$ as a function of $(x_{0}-x_{b})k_{c}$ and wave slope, including breaking and non-breaking packets. Non-breaking cases are indicated by dashed lines while breaking cases are solid lines. Breaking threshold in these simulations is $S=0.31$ , clearly separating the horizontal profiles. For each case, the horizontal solid line indicates the breaking area, $[x_{b,n}-\unicode[STIX]{x1D706}_{c}/2:x_{b,n}+\unicode[STIX]{x1D706}_{c}/2]$ . The width of the region of enhanced drift increases with the slope, while $x_{b,n}$ decreases as $S$ increases.

Figure 5(a) shows the total drift at the surface normalized by the central wavelength, $\unicode[STIX]{x0394}x(x_{0},z_{0}=0)/\unicode[STIX]{x1D706}_{c}$ , as a function of the normalized horizontal initial position $(x_{0}-x_{b})k_{c}$ for a typical breaking event ( $S=0.38$ ), both the raw data $\unicode[STIX]{x0394}x$ , each data point corresponding to one particle, as well as a running average drift $\overline{\unicode[STIX]{x0394}x}(x_{0})$ as a solid line (with the horizontal averaging window taken to be the characteristic wavelength $\unicode[STIX]{x1D706}_{c}$ ). The drift is small far upstream and downstream of the breaking area and increases strongly near the focusing point, before decreasing again after the focusing point. The extent of the breaking area is approximately given by the central wavelength of the packet $\unicode[STIX]{x1D706}_{c}$ . The location $x_{b,n}$ is defined as the peak location of the drift (found to be slightly after $x_{b}$ in this case), with the horizontal dashed line on the figure showing the area $[x_{b,n}-\unicode[STIX]{x1D706}_{c}/2:x_{b,n}+\unicode[STIX]{x1D706}_{c}/2]$ .

Figure 5(b) shows the normalized drift $\overline{\unicode[STIX]{x0394}x}/\unicode[STIX]{x1D706}_{c}$ , as a function of $(x_{0}-x_{b})k_{c}$ , for increasing slopes, for both non-breaking and breaking focusing wave packets. For each case, the maximum location is extracted, $x_{b,n}(S)$ , and the breaking area, $[x_{b,n}-\unicode[STIX]{x1D706}_{c}/2:x_{b,n}+\unicode[STIX]{x1D706}_{c}/2]$ , is indicated by horizontal dashed lines. For smaller slopes, the drift at the surface is almost constant as a function of the horizontal position ( $S<0.2$ ). For intermediate slopes of non-breaking wave packets ( $0.2<S<0.31$ ), a significant increase of the drift around the focusing location is observed. This is again in agreement with the argument discussed earlier for weakly nonlinear packets (appendix B). The breaking cases exhibit a strong increase in the horizontal drift around the breaking area. The breaking location defined as the location of the maximum drift $x_{b,n}$ generally moves upstream as the slope $S$ increases, which is consistent with earlier reports on breaking waves by focusing packets (Rapp & Melville Reference Rapp and Melville1990). The width of the region of enhanced drift increases with the slope, i.e. enhanced breaking strength leads to a larger area of enhanced drift.

4.4 Vertical profiles at the focusing point

We now look at the vertical profile of the horizontal drift. We compute the horizontally averaged displacement $\langle \overline{\unicode[STIX]{x0394}x}(z_{0})\rangle$ , for particles initially located at $z_{i,0}=z_{0}$ , and within the horizontal breaking area defined previously, i.e. particles initially located between $x_{b,n}-\unicode[STIX]{x1D706}_{c}/2$ and $x_{b,n}+\unicode[STIX]{x1D706}_{c}/2$ .

To facilitate the comparison with the Stokes drift velocity defined in § 2 and in the literature, we now define the drift velocity $u_{L}(z_{0})$ , with $T$ the typical time of the wave packet, taken as $T=1/f_{c}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{c}$ ,

Note that $u_{L}(z_{0})/c=\langle \overline{\unicode[STIX]{x0394}x}(z_{0})\rangle /(\unicode[STIX]{x1D706}_{c})$ .

Figure 6. Vertical profiles of the normalized horizontal particle drift $u_{L}/c=\langle \overline{\unicode[STIX]{x0394}x}\rangle /\unicode[STIX]{x1D706}_{c}$ , as a function of the initial depth of the particles $z_{0}k_{c}$ , obtained for each depth by averaging over a horizontal window $x\in [x_{b,n}-\unicode[STIX]{x1D706}_{c}/2:x_{b,n}+\unicode[STIX]{x1D706}_{c}/2]$ around the maximum drift location $x_{b,n}$ , for increasing slope $S$ (from blue to red). The drift at the surface increases with the slope, as does the return flow in the bottom half of the tank, ensuring mass conservation. Non-breaking cases are indicated by dashed lines while breaking cases are solid lines.

Figure 6 shows the normalized vertical profiles of the horizontal drift $u_{L}(z_{0})/c$ , for increasing wave slope $S$ , from small amplitude non-breaking packets (in blue) to stronger breaking waves (in red). The drift near the surface in the propagation direction increases significantly with the wave slope. The increase in the drift at the surface is further enhanced until we reach the breaking threshold ( $S>0.31$ here) where the drift shows a strong increase. It corresponds to the added Lagrangian transport due to breaking. This added drift due to breaking increases with the breaking strength. Moreover, a return flow is observed in the bottom part of the tank (extending from the bottom a distance $0.2\unicode[STIX]{x1D706}_{c}$ , or $0.4d$ ), with an intensity increasing with the slope. This return flow balances the downstream drift near the surface assuring mass conservation. Note that we have verified that $\int \langle \overline{\unicode[STIX]{x0394}x}(z_{0})\rangle \,\text{d}z_{0}=0$ , and that mass is indeed conserved, up to 0.2 % (as discussed in the appendix A for the total mass budget).

The general shape of the profile does not change significantly if the running averaging window is increased or reduced by a factor of $2$ ; however, the magnitudes of both the drift at the surface and the return flow as the bottom are increased and reduced, respectively, since the averaging process gives more or less weight to the enhancement at the breaking point. The choice of a window of $\unicode[STIX]{x1D706}_{c}$ appears to be the most consistent choice since it corresponds to the phase averaging used to compute the classical Stokes drift.

Figure 7. Rescaled vertical profiles of the particle drift velocity, averaged over a horizontal window $x\in [x_{b,n}-\unicode[STIX]{x1D706}_{c}/2:x_{b,n}+\unicode[STIX]{x1D706}_{c}/2]$ , for increasing slope $S$ . Slope is colour coded (same as in figure 6). (a,b) Rescaling of the non-breaking wave packet, showing $z_{0}k_{c}$ in the $y$ -axis and $u_{L}/c$ (a) and $u_{L}/(cS^{2})$ (b) in the $x$ -axis. The $S^{2}$ scaling for the drift at the surface is observed. Dashed line is $u_{L}/c=\unicode[STIX]{x1D712}_{1}S^{2}\text{e}^{-2k_{c}z_{0}}$ with $\unicode[STIX]{x1D712}_{1}=0.82$ obtained from (2.17). Solid line is the linear theory given by (2.16). A reasonable collapse of all data onto the theoretical curve is observed, the e-folding scale is then $1/2k_{c}$ . (c,d) Rescaling of the breaking packets. (c)  $u_{L}/c$ in $x$ -axis as a function of $z_{0}k_{c}$ in $y$ -axis. (d)  $u_{L}/cS$ in $x$ -axis as a function of $z_{0}k_{c}/(S-S_{\ast })$ in $y$ -axis. Better rescaling of the data is observed in (d), showing that the vertical scale is given by $(S-S_{\ast })/k_{c}$ , i.e. the height above the breaking threshold, while the drift at the surface also scales with $S$ . Dashed line is an exponential fit, $u_{L}/cS=\unicode[STIX]{x1D704}_{1}\exp [\unicode[STIX]{x1D704}_{2}(z_{0}k_{c})/(S-S_{\ast })]$ , with $\unicode[STIX]{x1D704}_{1}=2$ , $\unicode[STIX]{x1D704}_{2}=2$ fitted constant, and dotted line is an algebraic decay $u_{L}/cS=\unicode[STIX]{x1D704}_{3}/[1-(z_{0}k_{c})/(S-S_{\ast })]^{\unicode[STIX]{x1D704}_{4}}$ , with $\unicode[STIX]{x1D704}_{3}=2$ and $\unicode[STIX]{x1D704}_{4}=3$ fitted constants.

Figure 7 shows the rescaled vertical profiles of the horizontal velocity $u_{L}/c$ , with $c$ the characteristic phase velocity of the packet, as a function of the rescaled depth $z_{0}k_{c}$ .

Figure 7(a) shows the normalized Lagrangian drift $u_{L}/c$ as a function of $z_{0}k_{c}$ . The drift near the surface increases with the square of the slope $S^{2}$ and is well described by an exponential decay with depth. Figure 7(b) shows the drift velocity rescaled by the slope squared, as suggested by the Stokes drift $u_{L}/(cS^{2})$ . All data in the upper part ( $z_{0}k_{c}>-2$ ) of the tank collapse reasonably well on a single exponential profile. The solid line corresponds to the theoretical prediction given by (2.16), where the amplitude of each wave component was computed using the initial wave packet given by (3.2). Good agreement is observed between the numerical data and the theoretical profile, with a slight deviation observed for the largest non-breaking slope considered here, which might be related to nonlinear effects. The dashed line corresponds to an exponential decay $\unicode[STIX]{x1D712}_{1}S^{2}c\text{e}^{-2k_{c}z_{0}}$ , where $\unicode[STIX]{x1D712}_{1}=0.82$ is obtained from the Stokes drift at the surface, equation (2.17), where the amplitude of each wave component was computed using the initial wave packet given by (3.2). This simplified equation characterizes the bulk scale behaviour of the data, showing that the e-folding scale is given by the inverse of twice the characteristic wavenumber $1/2k_{c}$ for the non-breaking packets.

The dependency of the drift with $S^{2}$ at the surface and the e-folding scale $1/2k_{c}$ are no longer valid once the wave breaks, and the data are no longer described by this model. Figure 7(c,d) shows the breaking cases. In figure 7(c) $u_{L}/c$ is shown as a function of $k_{c}z_{0}$ , with the drift at the surface increasing with the slope $S$ . Figure 7(d) shows $u_{L}/cS$ as a function of $k_{c}z_{0}/(S-S_{\ast })$ , where $S_{\ast }=0.31$ is the breaking threshold in the present configuration. The typical vertical scale of the profile is now $k_{c}/(S-S_{\ast })$ , i.e. the height above the breaking threshold, while the drift at the surface scales (and increases) with $S$ . All data can be reasonably described by a single curve. The scaled profile can be described by an exponential fit, with $(S-S_{\ast })/k_{c}$ the e-folding scale, $u_{L}/cS=\unicode[STIX]{x1D704}_{1}\exp [\unicode[STIX]{x1D704}_{2}(z_{0}k_{c})/(S-S_{\ast })]$ or by an algebraic decay with depth, with the same characteristic vertical length scale, $u_{L}/cS=u_{L}/cS=\unicode[STIX]{x1D704}_{3}/[1-(z_{0}k_{c})/(S-S_{\ast })]^{\unicode[STIX]{x1D704}_{4}}$ , where $\unicode[STIX]{x1D704}_{1,2,3,4}$ are fitted constant. The latter gives a slightly more accurate vertical profile. Note that in the vicinity of the surface, the drift can be described by a linear decay with depth, which arises naturally from a Taylor expansion of the previous formula (i.e. for $(z_{0}k_{c})/(S-S_{\ast })\ll 1$ ), $u_{L}/cS=\unicode[STIX]{x1D704}_{3}/[1-(z_{0}k_{c})/(S-S_{\ast })]^{\unicode[STIX]{x1D704}_{4}}\approx \unicode[STIX]{x1D704}_{3}[1+\unicode[STIX]{x1D704}_{4}(z_{0}k_{c})/(S-S_{\ast })]$ . Furthermore this leads to $u_{L}(z_{0}\approx 0)/c\approx \unicode[STIX]{x1D704}_{4}S$ , at first order, which is compatible with the scaling argument developed in § 2.

Thus, the characteristic vertical length scale in the breaking cases is given by the height of the breaking wave (above the threshold for breaking). This is consistent with the laboratory work from Rapp & Melville (Reference Rapp and Melville1990) that shows that the depth of the broken fluid scales with the wave height. This further shows that the nature of the Lagrangian kinematics and transport is changed when we move from non-breaking to breaking waves. The fact that the drift due to breaking waves can be described by an algebraic decay with depth is a striking result as it represents a large departure from the evanescent models used almost ubiquitously in models of surface gravity waves.

4.5 Lagrangian drift at the surface

We now focus on the drift velocity at the surface,

Figure 8. Surface drift velocity $u_{L,s}/c$ , as a function of the wave slope $S$ . $u_{L,s}$ increases with the slope. For non-breaking waves, $u_{L,s}\propto S^{2}$ and can be described by the linear theory (solid red line), equation (2.17), or by the simplified equation, $u_{L,s}/c=\unicode[STIX]{x1D712}_{1}S^{2}$ (black dotted line), with $\unicode[STIX]{x1D712}_{1}=0.82$ a fitted parameter. The drift experiences a sharp increase for slopes above $S_{\ast }=0.31$ , which correspond to the appearance of breaking waves. For breaking waves we have $u_{L,s}\propto S$ , as described by our scaling argument in § 2. The drift is then described, within the scatter of the data, by $u_{L,s}/c=\unicode[STIX]{x1D712}_{2}(S-S_{\ast })$ , with $\unicode[STIX]{x1D712}_{2}=9$ a fitted parameter. Circle and diamond symbols are numerical data for different runs, for two different mesh resolutions: open symbols are medium resolution and full symbols are high resolution, while diamonds are for $f_{c}=0.87~\text{Hz}$ , $\unicode[STIX]{x0394}f/f_{c}=0.74$ , and circles are for $f_{c}=0.89~\text{Hz}$ , $\unicode[STIX]{x0394}f/f_{c}=0.75$ . The yellow star is estimated from the experiments of Grue & Kolaas (Reference Grue and Kolaas2017) obtained in intermediate depth water.

Figure 8 shows the surface Lagrangian drift velocity $u_{L,s}$ , scaled by $c$ , the characteristic phase velocity of the packet. For non-breaking packets, the drift increases with $S^{2}$ as expected by linear and weakly nonlinear theory, and as already discussed before. We also show on figure 8 the data from Grue & Kolaas (Reference Grue and Kolaas2017) obtained for periodic trains of steep non-breaking waves in intermediate depth water (their wave frequency is 0.8 Hz, and water depth is 0.2 m). A sudden increase of the Lagrangian drift is observed for breaking packets, where the drift increases linearly with the wave slope above threshold $S-S_{\ast }$ .

The drift velocity at the surface for the non-breaking packets is well described by (2.17), with the amplitude of the wave component being evaluated using the initial shape of the wave packet, equation (3.2), as shown by the solid line on figure 8. Note that since our packet has a constant slope, the non-breaking packets can also be described by this simplified equation,

which is the dotted line in figure 8, where $\unicode[STIX]{x1D712}_{1}=0.82$ is again obtained from the Stokes drift at the surface, equation (2.17), where the amplitude of each wave component was computed using the initial wave packet given by (3.2). This was already observed when rescaling the vertical profiles. These relations are valid for non-breaking packets, i.e. with slope $S<0.31$ in the present configuration.

When the waves start to break ( $S>0.31$ ), the Lagrangian drift increases linearly with the slope above threshold $S-S_{\ast }$ and can be modelled by

which is the dashed line in figure 8. The linear increase of the drift with $S$ is suggested by the scaling argument presented in § 2 and $\unicode[STIX]{x1D712}_{2}$ is fitted to the data and is indeed $O(10)$ as expected by our scaling argument for the particle acceleration. $S_{\ast }$ is an empirical breaking threshold determined for these focusing experiments.