# Stochastic particle transport by deep-water irregular breaking waves

D. Eeltink1,2,3  R. Calvert4,5  J.E. Swagemakers4  Qian Xiao2  T.S. van den
Bremer2,4<EMAIL_ADDRESS>1Department of Mechanical Engineering,
Massachusetts Institute of Technology, Cambridge, Massachusetts, USA
2Department of Engineering Science, University of Oxford, Parks Road, Oxford
OX1 3PJ, UK 3 Laboratory of Theoretical Physics of Nanosystems, EPFL, Ch-1015
Lausanne, Switzerland 4Faculty of Civil Engineering and Geosciences, Delft
University of Technology, 2628 CD, Delft, The Netherlands 5School of
Engineering, University of Edinburgh, Edinburgh, EH9 3FB, UK

###### Abstract

Correct prediction of particle transport by surface waves is crucial in many
practical applications such as search and rescue or salvage operations and
pollution tracking and clean-up efforts. Recent results by Deike et al. (2017)
and Pizzo et al. (2019) have indicated transport by deep-water breaking waves
is enhanced compared to non-breaking waves. To model particle transport in
irregular waves, some of which break, we develop a stochastic differential
equation describing both mean particle transport and its uncertainty. The
equation combines a Brownian motion, which captures non-breaking drift-
diffusion effects, and a compound Poisson process, which captures jumps in
particle positions due to breaking. From the corresponding Fokker–Planck
equation for the evolution of the probability density function for particle
position, we obtain closed-form expressions for its first three moments. We
corroborate these predictions with new experiments, in which we track large
numbers of particles in irregular breaking waves. For breaking and non-
breaking wave fields, our experiments confirm that the variance of the
particle position grows linearly with time, in accordance with Taylor’s
single-particle dispersion theory. For wave fields that include breaking, the
compound Poisson process increases the linear growth rate of the mean and
variance and introduces a finite skewness of the particle position
distribution.

##  1 Introduction

Correctly understanding and predicting the motion of objects on the ocean
surface is crucial for several applications. For example, floating plastic
marine litter has rapidly become an acute environmental problem (e.g., Eriksen
et al. (2014)). A mismatch exists between the estimated amount of land-
generated plastic entering coastal waters ($5$-$12$ million tonnes yr-1,
Jambeck et al. (2015)) and the estimated total amount of plastic floating at
sea (less than $0.3$ million tonnes, Cózar et al. (2014); Eriksen et al.
(2014); van Sebille et al. (2015)), which necessitates more accurate transport
and dispersion models (van Sebille et al., 2020). In addition, efforts to
clean up floating plastics rely on an accurate prediction of the distribution
and trajectories of particles to deploy clean-up devices at the right location
(e.g., Sainte-Rose et al. (2020)). Similarly, estimating environmental impact
for oil spills and open-sea rescue operations rely crucially on transport
models (e.g., Christensen et al. (2018)).

During the periodic motion of a (deep-water) surface gravity wave, a fluid
parcel does not follow a perfectly circular trajectory. Instead, it
experiences a net drift in the direction of wave propagation, known as the
Stokes drift (Stokes, 1847). Wave models such as WAM and WaveWatch III predict
wave fields properties averaged over longer timescales (typically 3 hours),
based on which estimates of the mean Stokes drift over that period can be made
(Webb & Fox-Kemper, 2011; Breivik et al., 2014). This Stokes drift prediction
is typically superimposed onto a Eulerian flow field, often given by ocean
general circulation models or measurements. A number of recent studies have
shown that the inclusion of the Stokes drift in this way can alter the
predicted direction of plastic pollution transport, shifting convergence
regions (Dobler et al., 2019) and pushing microplastics closer to the coast
(Delandmeter & van Sebille, 2019; Onink et al., 2019).

However, for several reasons, the Stokes drift alone should not be the
velocity with which waves actually transport objects on the ocean surface.
First and foremost, it is the wave-induced Lagrangian-mean velocity, made up
from the sum of the Stokes drift and the wave-induced Eulerian-mean velocity,
with which waves transport particles. On the rotating Earth, the Coriolis
force in combination with the Stokes drift drive an Eulerian-mean current in
the turbulent upper-ocean boundary layer, known as the Ekman–Stokes flow,
which includes the effect of boundary-layer streaming (Longuet-Higgins, 1953).
This Ekman–Stokes flow needs to be added to the Stokes drift in order to
predict the wave-induced Lagrangian-mean flow with which particles are
transported (Higgins et al., 2020a), or wave and ocean circulation models need
to be properly coupled (e.g., Staneva et al. (2021)). The Ekman–Stokes flow
can have significant consequences for global floating marine litter
accumulation patterns (Cunningham et al., 2022). Wave-induced Eulerian-mean
currents can also play a role in laboratory experiments (van den Bremer et
al., 2019; van den Bremer & Breivik, 2018). Second, the properties of the
object itself such as its shape, size and buoyancy can cause the object’s
trajectory to be different from that of an infinitesimally small Lagrangian
particle with a different mean transport as a result (Santamaria et al., 2013;
Huang et al., 2016; Alsina et al., 2020; Calvert et al., 2021; DiBenedetto et
al., 2022). Third, breaking waves are known to transport particles much faster
than non-breaking waves. For steep waves, particles may surf on the wave
(Pizzo, 2017) and be subject to transport faster than the Stokes drift due to
wave breaking (Deike et al., 2017; Pizzo et al., 2019). In this paper, we will
focus on (almost) perfect Lagrangian particles and consider the influence of
unidirectional irregular deep-water waves and wave breaking, thus ignoring the
effects of Coriolis forces, wind and non-wave-driven currents, which evidently
also determine the drift of an object in the real ocean.

For non-breaking waves, the Stokes drift can be estimated estimated for
various sea states (e.g., Webb & Fox-Kemper (2011); Breivik et al. (2014))
with an important but not always resolved role for the spectral tail (Lenain &
Pizzo, 2020). In an irregular or random wave field, the Stokes drift becomes a
stochastic process itself. That is, due to the random wave field, particles
will diffuse with respect to this mean velocity, yielding a distribution in
their predicted position (Herterich & Hasselmann, 1982). Estimates of the
variance of the Stokes drift can either be obtained from the spectrum
(Herterich & Hasselmann, 1982) or from the joint distribution of the
significant wave height and peak period (Longuet-Higgins, 1983; Myrhaug et
al., 2014). Alternatively, deterministic simulations of the particle
trajectories can be performed (e.g., Li & Li (2021); Farazmand & Sapsis
(2019)), which could then allow for calculation of statistics of the evolution
for different initial conditions using Monte Carlo methods.

The particle transport of breaking waves can be described at different scales.
At the scale of individual waves or wave groups, deterministic models for
particle trajectories can be formulated to include wave breaking. Direct
numerical simulations of breaking focused wave groups show a linear scaling of
the net particle transport with the theoretical steepness of the wave group at
the focus point $\mathcal{S}$ (Deike et al., 2017), in contrast to the
quadratic scaling of Stokes drift with steepness for non-breaking waves.
Restrepo & Ramirez (2019) confirm the linear scaling with $\mathcal{S}$ and
calculate the variance of the drift. Experimental confirmation of the enhanced
drift and the linear scaling with steepness for wave groups is provided in
Lenain et al. (2019) and Sinnis et al. (2021), where Sinnis et al. (2021) also
consider the effects of bandwidth. Considering much larger scales, relevant
for application to the real ocean, requires a stochastic approach. Pizzo et
al. (2019) extend the result obtained by Deike et al. (2017) for wave groups
to sea states. Based on the wave spectrum (peak wavenumber) and wind speed,
one can estimate the breaking statistic $\Lambda(c){\rm d}c$ (Banner et al.,
2000; Dawson et al., 1993; Ochi & Tsai, 1983; Holthuijsen & Herbers, 1986;
Sullivan et al., 2007), defined as the average length of breaking crests
moving with a velocity in the range $(c,c+{\rm d}c)$, where $c$ is the phase
speed (Phillips, 1985). Consequently, the breaking drift speed found in Pizzo
et al. (2019) is weighted by the percentage of broken sea surface per unit
area, which in turn can be computed as a function of peak wavenumber and wind
speed. Comparing their prediction of the drift speed to the Stokes drift for
non-breaking waves shows that, as the wind velocity and wave steepness
increase, the wave breaking component of drift becomes more important and can
be as large as 30% of the Stokes drift for non-breaking waves (Pizzo et al.,
2019).

A series of papers have examined stochastic Stokes drift (Jansons & Lythe,
1998; Bena et al., 2000), focusing on diffusion in the case of two opposing
waves, for which the mean drift should be zero, yet diffusion still causes
transport. Their insights cannot immediately be applied to realistic ocean
waves. More generally, several authors have examined the Taylor particle
diffusion of a random surface gravity wave field (Herterich & Hasselmann,
1982; Sanderson & Okubo, 1988; Weichman & Glazman, 2000; Balk, 2002, 2006).
Bühler & Holmes-Cerfon (2009) derive the Taylor single-particle diffusivity
for random waves in a shallow-water system under the influence of the Coriolis
force and corroborate their results with Monte Carlo simulations. A generally
applicable stochastic framework is provided in Bühler & Guo (2015), who derive
a stochastic differential equation (SDE) for the particle position and a
Kolmogorov backward equation for particles along quasi-horizontal
stratification surfaces induced by small-amplitude internal gravity waves that
are forced by white noise and dissipated by nonlinear damping designed to
model attenuation of internal waves.

In this paper, we propose a stochastic differential equation (SDE) for the
evolution of particle position in a unidirectional irregular deep-water wave
field with wave breaking and obtain the corresponding Fokker–Planck equation
for the evolution of the distribution. The SDE combines a Brownian motion,
which captures non-breaking drift-diffusion effects, and a compound Poisson
process, which captures jumps in particle positions due to breaking. We focus
on the short-time regime over which the properties of the sea state (i.e., its
spectrum) stay constant and corroborate the predictions of our SDE with new
laboratory experiments in which we track a large number of particles.

The paper is organized as follows. First, §2 introduces the Brownian drift-
diffusion process to model particle position evolution without breaking and
the Poisson process to model the surfing behavior experienced when a particle
encounters a breaking crest. The corresponding Fokker–Planck equation for the
evolution of the probability density function of particle position and its
mean, variance and kurtosis are also derived in §2. Then, §3 outlines the wave
basin experiments performed, where particles were tracked in irregular waves.
In §4, we corroborate our theoretical predictions with the experimental
results for different wave steepnesses and, consequently, different fractions
of breaking waves. Finally, we conclude in §5.

##  2 Stochastic model for particle transport

For a given initial particle position $X_{0}=X(t=0)$, we seek to determine the
(long-term) evolution of the particle position $X(t)$ and its distribution,
where we will only consider wave-averaged (Eulerian-mean or Lagrangian-mean)
quantities. Particles are transported with the Lagrangian-mean velocity, which
consists of the sum of the Eulerian-mean velocity and the Stokes drift (e.g.,
Bühler (2014)). In our case, we consider the Lagrangian-mean velocity:

$u_{\text{L}}=u_{\text{E,NB}}+\underbrace{u_{\text{S}}+u_{\text{B}}}_{u_{\text{D}}},$
(1)

which consists of the sum of the Eulerian-mean velocity that excludes the
effects of wave breaking $u_{\text{E,NB}}$ and a drift velocity
$u_{\text{D}}$, which in turn consists of the Stokes drift for non-breaking
waves $u_{\text{S}}$ and a breaking contribution $u_{\text{B}}$. For
simplicity, we subsequently ignore the wave-induced Eulerian-mean velocity
$u_{\text{E,NB}}$ in our model, as its contribution to drift on the surface of
the (non-rotating) ocean is negligible for deep-water waves (e.g., van den
Bremer & Taylor (2015); Higgins et al. (2020b)). To compare to basin
experiments (in §3), we will take the effect of Eulerian flow (wave-induced or
otherwise) into account.

We propose to model the (wave-averaged) particle position $X(t)$ as a jump-
diffusion process for which the stochastic differential equation (SDE) can be
written as:

$\displaystyle{\rm d}X$ $\displaystyle=$
$\displaystyle\underbrace{b(X(t),t){\rm d}t}_{\text{Mean drift}}\hskip
14.22636pt$ $\displaystyle+$ $\displaystyle\underbrace{\sigma(X(t),t){\rm
d}W(t)}_{\text{Diffusion process}}\hskip 14.22636pt$ $\displaystyle+$
$\displaystyle\underbrace{{\rm d}J(t),}_{\text{Compound jump process}}$ (2)
$\displaystyle=$ $\displaystyle\langle u_{\text{S}}\rangle{\rm d}t\hskip
14.22636pt$ $\displaystyle+$ $\displaystyle\sigma{\rm d}W(t)\hskip 14.22636pt$
$\displaystyle+$ $\displaystyle{\rm d}J(\Lambda(\epsilon),G(s;\epsilon);t),$
(3)

where $b(X(t))=\langle u_{\text{S}}\rangle$ is the mean Stokes drift of a
stochastic or irregular non-breaking wave field (angular brackets denote the
mean of a stochastic process); $\sigma(X(t),t)=\sigma$ the standard deviation
of the wave-averaged drift rate caused by stochastic nature of the individual
waves, with $D=\sigma^{2}/2$ the resulting diffusion coefficient; and $W(t)$
denotes a Wiener process (Brownian motion). The effect of wave breaking is
captured by the compound Poisson process $J(t)$, which represents the jumps in
particle position induced by breaking. Note this process has a non-zero mean,
$\langle u_{\text{B}}\rangle\neq 0$, reflecting the contribution of breaking
to mean drift.

Taking the terms in (2) in turn, we will proceed to outline how particle
displacement can be viewed as a drift-diffusion process in the absence of
breaking waves (§2.1) and then introduce a compound-jump process to account
for wave breaking (§2.2). Finally, in §2.3, we will propose a Fokker–Planck
equation for the evolution of the probability density function of particle
position $P(X,t)$ and give analytical solutions for the first three moments of
$P(X,t)$ (given a Dirac delta function for the initial particle distribution).

###  2.1 Stochastic Stokes drift in the absence of breaking

The first term in (2) pertains to the average drift experienced by a particle
in the absence of breaking, known as the Stokes drift. For a deep-water,
monochromatic, unidirectional wave, the Stokes drift is given by (Stokes,
1847):

$U_{\text{S}}(z)=a^{2}\omega ke^{2kz},$ (4)

where $\omega$ is the angular frequency of the wave, $k$ its wavenumber, $a$
its amplitude, and $z$ the vertical position measured upwards from the still-
water level. The value of the Stokes drift at the surface $u_{\text{S}}$ is
consistently approximated as $u_{\text{S}}=U_{\text{S}}(z=0)$. Written in
terms of the (commonly estimated) wave period $T$ and wave height for periodic
linear waves $H=2a$:

$u_{\text{S}}=(ak)^{2}c=\frac{a^{2}\omega^{3}}{g}=\frac{(2\pi)^{3}}{g}\frac{H^{2}}{T^{3}},$
(5)

where $c=\omega/k$, and we have used the linear deep-water dispersion
relationship $\omega^{2}=gk$ with $g$ the gravitational constant. The Stokes
drift is proportional the square of the steepness $ak$, as is evident from
(5).

An irregular or stochastic wave field consists of a distribution of different
wave periods and heights. The mean Stokes drift can be obtained by summing up
the Stokes drift contributions of the different spectral components (Kenyon,
1969; Webb & Fox-Kemper, 2011):

$\langle u_{\text{S}}\rangle=\frac{16\pi^{3}}{g}\int_{0}^{\infty}S(f)f^{3}{\rm
d}f,$ (6)

where the unidirectional frequency spectral density $S(f)$ is defined so that
$\int_{0}^{\infty}S(f){\rm d}f=\langle\eta^{2}\rangle$, with $\eta(t)$ the
surface elevation time series, and $f$ the wave frequency. Naturally, for a
monochromatic wave, equation (6) gives the same result as (5).

The second term in (2), ${\rm d}X=\sigma{\rm d}W$, models stochastic
deviations from the mean as a Wiener or normal diffusion process.
Conceptually, deviations from the mean arise because, on the wave-averaged
time scale, each different wavelength and wave height in an irregular sea
contributes a different Stokes drift and, therefore, a different displacement.
Regardless of the distribution of the Stokes drift, as long as its variance is
finite, the central limit theorem states that this will result in a normal
distribution for particle position $X(t)$. For a normal process, the second
central moment or variance $\langle(X(t)-\langle X(t)\rangle)^{2}=2Dt$, where
$D$ is the diffusion coefficient. Assuming a stationary underlying random
process, the statistics of fluctuations in the drift velocity (i.e.,
$\tilde{u}_{\text{S}}=u_{\text{S}}-\langle u_{\text{S}}\rangle$) do not depend
on time. One can write
$D=\langle\tilde{u}_{\text{S}}^{2}\rangle\tau=\sigma_{\text{S}}^{2}\tau$ ,
with $\tau\propto\Delta\omega^{-1}$ the integral correlation time of the drift
velocity, $\Delta\omega$ a measure of the width of the underlying spectrum,
and with $\sigma_{\text{S}}^{2}=\langle\tilde{u}_{\text{S}}^{2}\rangle$ the
variance of the Stokes drift (Taylor, 1922; Herterich & Hasselmann, 1982;
Farazmand & Sapsis, 2019). Therefore,
$\sigma=\sqrt{2D}=\sqrt{2\tau}\sigma_{\text{S}}$ in (2).

To estimate the diffusion coefficient $D$, it is therefore necessary to obtain
a distribution for the Stokes drift $P(u_{\text{S}})$ and derive from this its
standard deviation $\sigma_{\text{S}}$. In Herterich & Hasselmann (1982) the
spectrum is assumed to be narrow, implying a constant (non-stochastic) value
for the period $T=T_{p}$. Consequently, the distribution for the Stokes drift
can be derived directly from the Rayleigh distribution for the wave height
using (5) (Longuet-Higgins, 1952). Since $u_{\text{S}}\propto H_{s}^{2}$, the
Stokes drift follows an exponential distribution. Specifically, the
probability density function for $u_{\text{S}}$ reads:

$P(u_{\text{S}};H_{\text{s}},T_{p})=\frac{gT_{p}^{3}}{4\pi^{3}}\frac{1}{H_{\text{s}}^{2}}\exp{\left[-\frac{gT_{p}^{3}}{4\pi^{3}}\frac{u_{\text{S}}}{H_{\text{s}}^{2}}\right]}.$
(7)

Alternatively, Myrhaug et al. (2014) and Longuet-Higgins (1983) derive an
exponential distribution for the Pierson–Moskowitz spectrum, taking into
account the joint distribution of wave heights and wave periods. In both
cases, an exponential distribution is obtained, which has the property that
the mean is equal to the standard deviation. We therefore have
$\sigma_{\text{S}}=\langle u_{\text{S}}\rangle$, which we will use to predict
the diffusion coefficient $D$.

###  2.2 Breaking encounters as a jump process

Our experiments will show (cf. §3 and figure 2 in particular) that for waves
that are steep enough so that they start to break, in addition to the slow
gradual drift and diffusion of the particles predicted for non-breaking waves,
jumps in particle positions occur when waves break. Each ‘jump’ event
represents an encounter of a particle with (the crest of) a breaking wave. We
model these jump events by a compound Poisson process $J(t)$ in (2):

$J(t)=\sum_{i=0}^{N_{\Lambda}(t)}s_{i}.$ (8)

where $N_{\Lambda}(t)$ is a counting of a Poisson point process, and the
arrival rate $\Lambda$ corresponds to the expected number of jumps per unit
time for the compound Poisson process.

We expect the arrival rate $\Lambda$ to increase with the steepness
$\epsilon$. Assuming a gradual increase of the arrival rate with steepness
followed by followed by saturation at large enough steepness, we propose a
three-parameter sigmoid function (see figure 3a):

$\Lambda(\epsilon;\tau_{\Lambda},\phi_{\Lambda},\epsilon_{0,\Lambda})=\frac{1/\tau_{\Lambda}}{1+\exp{[-\phi_{\Lambda}(\epsilon-\epsilon_{0,\Lambda})}]},$
(9)

where the parameters $\tau_{\Lambda}$, $\phi_{\Lambda}$ and
$\epsilon_{0,\Lambda}$ will be estimated from our experimental data (see §3).
Furthermore, we assume that the amplitude of the jumps $s_{i}$ follows a two-
parameter Gamma distribution with probability density function
$G(s;\alpha(\epsilon),\beta(\epsilon))$ (see figure 3b):

$G(s;\alpha,\beta)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s},$ (10)

where $\alpha>0$ is the shape parameter, $\beta>0$ the rate parameter and
$\Gamma(\alpha)$ a Gamma function. Based on experimental data, we will show in
§3 that both $\alpha$ and $\beta$ can be effectively modelled as linear
functions of steepness $\epsilon$.

###  2.3 Fokker–Planck equation and analytical solutions

Two methods exist to integrate the stochastic differential equation (2) in
order to obtain a probability distribution for particle position $X(t)$. Monte
Carlo simulations of (2) can be numerically integrated (e.g., Milstein (1995);
Higham (2001)), and an empirical probability density function can be obtained
at each time step from the statistics of these trajectories. Alternatively,
the evolution equation for the probability density function $P(x,t)$ of the
random variable $X(t)$ corresponding to (2), the so-called Fokker–Planck
equation, can be solved directly. For (2), the Fokker–Planck equation is given
by (e.g., Gardiner (1983); Denisov et al. (2009); Gaviraghi et al. (2017)):

$\frac{\partial}{\partial t}P(x,t)=-\langle
u_{\text{S}}\rangle\frac{\partial}{\partial
x}P(x,t)+\frac{\sigma^{2}}{2}\frac{\partial^{2}}{\partial x^{2}}P(x,t)-\Lambda
P(x,t)+\Lambda\int_{-\infty}^{\infty}{\rm
d}x^{\prime}P(x^{\prime},t)J(x-x^{\prime}).$ (11)

This partial differential equation can be solved numerically (e.g., Gaviraghi
(2017)). Alternatively, by employing the characteristic functions of the
distribution,

$\hat{P}(l)=\int_{-\infty}^{\infty}{\rm d}xe^{ilx}P(x),$ (12)

we obtain the ordinary differential equation

$\frac{\partial}{\partial t}\hat{P}(l)=i\langle u_{\text{S}}\rangle
l\hat{P}(l)-\frac{\sigma^{2}}{2}l^{2}\hat{P}(l)-\Lambda\hat{P}(l)+\Lambda\hat{G}\hat{P}(l),$
(13)

where the characteristic function of the gamma distribution is

$\hat{G}(l)=\left(1-i\frac{l}{\beta}\right)^{-\alpha}.$ (14)

Equation 13 can be solved exactly, with solution:

$\hat{P}(l,t)=A(l)\exp{\left[\left(i\langle u_{\text{S}}\rangle
l-\frac{\sigma^{2}}{2}l^{2}-\Lambda+\Lambda\left(1-i\frac{l}{\beta}\right)^{-\alpha}\right)t\right]},$
(15)

where $A(l)$ is an unknown function that is determined by the initial
condition.

Raw moments of a probability density function can be readily evaluated using
its characteristic function:

$m_{n}=\langle X^{n}\rangle=\int_{-\infty}^{\infty}{\rm
d}xx^{n}e^{ilx}P(x)=\left.(-i)^{n}\frac{\partial^{n}\hat{P}(l)}{\partial
dl^{n}}\right|_{l=0}.$ (16)

For a Dirac delta function as the initial condition (i.e., $\hat{P}(t=0)=1$),
we have $A=1$ in (15) and we can obtain exact solutions for the first (three)
central moments of (11):

$\displaystyle\langle X(t)\rangle$ $\displaystyle=$
$\displaystyle\left.-i\frac{\partial\hat{P}(l,t)}{\partial l}\right|_{l=0}$
$\displaystyle=$ $\displaystyle\left(\langle
u_{\text{S}}\rangle+\frac{\alpha}{\beta}\Lambda\right)t,$ (17)
$\displaystyle\langle\tilde{X}(t)^{2}\rangle$ $\displaystyle=$
$\displaystyle\left.-\frac{\partial^{2}\hat{P}(l,t)}{\partial
l^{2}}\right|_{l=0}-m_{1}^{2}$ $\displaystyle=$
$\displaystyle\left(\sigma^{2}+\frac{\alpha(\alpha+1)}{\beta}\Lambda\right)t,$
(18) $\displaystyle\langle\tilde{X}(t)^{3}\rangle$ $\displaystyle=$
$\displaystyle\left.i\frac{\partial^{3}\hat{P}(l,t)}{\partial
l^{3}}\right|_{l=0}-3m_{1}m_{2}+2m_{1}^{3}$ $\displaystyle=$
$\displaystyle\frac{\alpha(1+\alpha)(2+\alpha)}{\beta^{2}}\Lambda t,$ (19)

where $\tilde{X}(t)=X(t)-\langle X(t)\rangle$. In (17)-(19), the expected
Stokes drift for non-breaking waves $\langle u_{\text{S}}\rangle(\epsilon)$,
the arrival rate $\Lambda(\epsilon)$, and the shape and scale factors
$\alpha(\epsilon)$ and $\beta(\epsilon)$ are all functions of steepness; their
dependence on $\epsilon$ will be estimated from experimental data in the next
section. Note from (17)-(19) that each central moment is linearly increasing
with time. Breaking increases both the mean drift (cf. (17)) and the variance
of particle position (cf. (18)); it also introduces a non-zero (positive)
skewness not predicted for non-breaking waves (cf. (19)).

##  3 Laboratory experiments

###  3.1 Laboratory set-up and input conditions

Figure 1: Experimental set-up, indicating the overhead particle-tracking
camera’s field of view (FOV), wave gauges (red dots), and EMS velocity
measurement probes (blue dots).

####  3.1.1 Wave conditions

Experiments were performed in the 8.7 m wide and 75 m long Atlantic Basin at
Deltares, the Netherlands. Figure 1 provides an overview of the set-up. The
basin was equipped with a segmented piston-type wavemaker, consisting of 20
wave paddles at one end and an absorbing beach at the other. The water depth
$d$ was 1 m. The experiments were carried out with irregular waves prescribed
by the JONSWAP spectral density $S(\omega)$,

$S(\omega)=\frac{Kg^{2}}{\omega^{5}}\exp\left[-\frac{5}{4}\left(\frac{\omega_{p}}{\omega}\right)^{4}\right]\gamma^{r},$
(20)

with $\omega$ the frequency of the waves (in rad/s), $g$ the gravitational
acceleration, $\omega_{p}$ the peak frequency, and
$r=\exp\left[-(\omega-\omega_{p})^{2}/(2\sigma_{\text{S}}^{2}\omega_{p}^{2})\right]$
with (non-dimensional) spectral width $\sigma_{\text{S}}=0.07$ for
$\omega\leq\omega_{p}$ and $\sigma_{\text{S}}=0.09$ for $\omega>\omega_{p}$.
The shape parameter $\gamma$ captures the ‘peakedness’ and thereby also the
bandwidth of the spectrum, and we set $\gamma=3.3$ for all experiments. The
bandwidth is important as it determines the correlation time $\tau$ and
thereby the variance of particle position predicted by (2) (i.e.,
$\sigma=\sqrt{2D}=\sqrt{2\tau}\sigma_{\text{S}}$). Using a second moment of
the spectrum, we calculate the spectral width as $\Delta\omega=1.39$ rad/s and
set $\tau=1/\Delta\omega$. The prefactor $K$ in (20) is adjusted to obtain the
desired significant wave height $H_{s}$. Phases of a discretised spectrum were
chosen randomly in order to create an irregular wave times series of 30 min
duration, which was used as (linear) forcing to the wavemakers. Reflections
were generally less than 5%.

$H_{s}$ | $H_{\text{s,exp}}$ | $\epsilon_{\text{exp}}$ | $\Delta t$ | $N_{\text{traj}}$ | $\Lambda$ | $\langle u_{\text{S}}\rangle$ | $\langle u_{\text{L,exp}}\rangle$ | $\langle u_{\text{E,fit}}\rangle$
---|---|---|---|---|---|---|---|---
(m) | (m) |  | (s) |  | (1/s) | (mm/s) | (mm/s) | (mm/s)
0.050 | 0.053 | 0.074 | 257 | 143 | $3.25\times 10^{-5}$ | 13.1 | 13.7 | +4.1
0.090 | 0.087 | 0.122 | 167 | 136 | $1.42\times 10^{-2}$ | 24.9 | 22.5 | +0.8
0.120 | 0.115 | 0.162 | 122 | 85 | $6.96\times 10^{-2}$ | 35.3 | 30.0 | -5.7
0.170 | 0.132 | 0.185 | 143 | 113 | $6.59\times 10^{-2}$ | 43.8 | 31.2 | -10.8
Table 1: Overview of experiments and parameter values. For all experiments
$T_{p}=1.2$ s. The subscript exp refers to the values measured in experiments,
$\Delta t$ is the time length of a (segmented) trajectory, $N_{\text{traj}}$
corresponds to the number of spheres tracked, and $\Lambda$ is the arrival
rate of jumps in particle position. The mean Stokes drift $\langle u_{\rm
S}\rangle$ is calculated using (6), the mean Lagrangian velocity $\langle
u_{\rm L,exp}\rangle$ is obtained from tracer particle positions in
experiments, and the Eulerian-mean velocity $\langle u_{\rm E,fit}\rangle$ is
obtained so that experiments and model predictions for the mean agree (see
§3.1.3).

Parameter values chosen for the experiments are given in table 1. The peak
period $T_{p}=1.2$ s, and four (input) significant wave heights are examined,
$H_{s}=$ 0.05, 0.09, 0.12, and 0.17 m (input), with an increasing fraction of
breaking waves. Defining a characteristic steepness as
$\epsilon=k_{p}H_{s}/2$, with $k_{p}$ the wavenumber corresponding to $T_{p}$,
we obtain $\epsilon=$ 0.070, 0.126, 0.168, and 0.238 (input). Experiments were
performed in deep water ($k_{p}d=2.8$). A total of 6 wave gauges were placed
in the basin. The time series of the surface elevation $\eta(t)$ at wave gauge
2 (co-located with the EMS probes) was used to calculate the values of the
parameters measured in experiments reported in table 1, where they are
labelled with the subscript expt to be distinguished from input values, and
the power spectrum $S(f)$ (used to estimate the Stokes drift according to
(6)). From table 1 it is evident that $H_{s}$ is under-produced in the
experiments for the larger values of $H_{s}$, which is in large part due to
wave breaking.

####  3.1.2 Lagrangian tracer particles

The tracer particles were 20 mm diameter yellow polypropylene spheres, which
were chosen to be as small as possible, but large enough to be tracked by the
overhead camera. The density of the particles was 920 kg/m3, so the particles
were as submerged as possible so that they would follow the motion of the
fluid, but remained detectable from above. The experiments were carried out in
fresh water (998 kg/m3). An automated device was used to drop a set of
approximately 20 spheres into the basin every 10 s, with a spacing of 15 cm
along a 3.0 m spanwise section of the basin ($y$-direction) a short distance
(in the $x$-direction) before entering the camera’s field of view. A Z-CAM
camera was mounted at a height of 11.9 m above the basin, allowing it to
capture an area of approximately 8 m along the length of the basin and 6 m of
its width (see figure 1).

To examine whether our tracer particles behave as idealized Lagrangian
particles, we compare the measured velocity spectrum of the tracer particles
to the velocity spectrum obtained from the measured surface elevation using
linear wave theory in appendix D. For frequencies up to a high-frequency cut-
off that lies much above the spectral peak, we find good agreement between the
two, confirming Lagrangian behaviour of the tracer particles. Based on Calvert
et al. (2021), we estimate that spherical particles with a diameter of less
than 10% of the wavelength typically behave as Lagrangian particles. For $20$
mm diameter spheres, this corresponds to a cut-off frequency of $2.8$ Hz (from
linear dispersion), which in turn agrees with what we find in the
aforementioned spectral comparison in appendix D. Since $D/\lambda_{p}=0.9\%$,
we do not expect our tracer particles to experience enhanced transport
compared to the Stokes drift in non-breaking waves due to mechanisms described
in Calvert et al. (2021).

####  3.1.3 Role of Eulerian-mean flow

While we do not take into account the Eulerian-mean flow $u_{\text{E}}$ in our
model (2), this flow is present in experiments and therefore has to be
accounted for to enable comparison (cf. (1)). Wave-induced Eulerian-mean flows
have often prevented observation of Stokes drift in laboratory wave flumes;
they are notoriously difficult to predict, specific to each laboratory basin
and experiment and, when observed in the laboratory, not representative of
wave-induced Eulerian-mean flows in the field (see reviews by van den Bremer &
Breivik (2017) and Monismith (2020)).

The Eulerian flow was measured using three EMS velocity measurement devices at
three different depths (always fully submerged) and one horizontal position
(($x,y=14.2,4.05$-$4.65)$ m). This allowed estimation of the basin-specific
Eulerian-mean flow for each experimental condition at fixed depths of $z=$
$-20$, $-30$, and $-50$ cm with $z=0$ the still-water level.

Although we have measured the Eulerian flow directly, we do not use the
Eulerian-mean flow that we obtain from these measurements (by wave averaging)
directly in the comparison between experiments and model predictions (cf.
(1)). The reason is that these Eulerian flow measurements, for reasons of
practicality, are only at a single point in space ($x$, $y$) and at a certain
distance below the free surface, making extrapolation to the surface a
potential source of error. Nevertheless, we can infer from the measurements
that the Eulerian-mean flow is non-stochastic, that is, it does not show
variability on the same time scale and with the same order of magnitude as the
measured Lagrangian-mean velocity.

We therefore correct $b$, the mean non-breaking drift in our model (2), with
an arbitrary (not measured) Eulerian-mean flow $\langle u_{\rm E,fit}\rangle$,
so that $b=\langle u_{\rm S}\rangle+\langle u_{\rm E,fit}\rangle$. The mean
Stokes drift $\langle u_{\rm S}\rangle$ is calculated using (6) from (a
JONSWAP spectrum fitted to) the measured surface elevation spectrum. The
arbitrary value of the Eulerian-mean flow $\langle u_{\rm E,fit}\rangle$ is
then chosen so that the mean Lagrangian drift predicted by the model $\langle
u_{\text{L,mod}}\rangle={\rm d}\langle X(t)\rangle/{\rm d}t$ (using $b=\langle
u_{\rm S}\rangle+\langle u_{\rm E,fit}\rangle$) is equal to the mean
Lagrangian drift in experiments $\langle u_{\text{L,exp}}\rangle$, that is
$\langle u_{\text{L,mod}}\rangle=\langle u_{\text{L,exp}}\rangle$. Our
interest in §4, where we compare model predictions to experiments, is
therefore in higher-order moments of particle position. The values $\langle
u_{\rm E,fit}\rangle$ thus obtained are within a reasonable margin of the
values measured at depth $\langle u_{\rm E,exp}\rangle$ (see table 3 in
appendix C). We note for completeness that for non-breaking waves, we expect
that $\langle u_{\text{L,mod}}\rangle=\langle u_{\text{S}}\rangle+\langle
u_{\text{E,fit}}\rangle$, while for breaking waves the jumps also make a
contribution to the mean, $\langle u_{\text{B}}\rangle$. The next section will
explain how this contribution is calculated, which must be done before
$\langle u_{\rm E,fit}\rangle$ can be found. Table 1 gives the values of the
different mean velocities for the different experiments.

###  3.2 Trajectory data processing

Our goal is to predict particle trajectories and properties of the probability
distribution of particle position based on information of the spectrum of the
waves. To this end, we first have to process the camera images to obtain
particle trajectories. We subsequently use these trajectories to obtain
properties of the jump process that is used to model breaking.

####  3.2.1 Image processing

The yellow spheres were tracked using OpenCV. The spheres were identified
using a Hue saturation filter and subsequently tracked using the CSRT
algorithm, obtaining the sub-pixel location ($x_{p},y_{p}$) of the centre of
each sphere at each point in time. The time step is determined by the frame
rate of the camera of 24 Hz. The trajectories were undistorted by calculating
the camera intrinsics using a checkerboard with 75 mm squares. The pixel
locations were then transformed into basin coordinates ($x,y$), assumed to be
in the plane of the still-water level, $z=0$, defined by a floating
checkerboard at a known position. See appendix A for further details.

####  3.2.2 Creating sample trajectories

To create samples that can be used to compare to predictions of our stochastic
model, particle trajectories $X_{i}(t)$ were segmented to trajectory lengths
of $\Delta t$ and all offset to have initial position $X_{i}(t=0)=0$, as
illustrated in figure 2. The resulting initial distribution is a Dirac delta
function, and we obtain $N_{\rm traj}$ trajectories of equal length $\Delta t$
(for each significant wave height).

The dashed lines in the insets in figure 2 show particle position $X(t)$ at
the time resolution of the camera (24 Hz), showing the oscillatory motion of
the particle with every wave. The stochastic model (2) is valid for a large
enough time scale, so that the oscillatory effects of the waves are averaged
out, but the stochastic effect of the irregular wave field on particle
transport is retained. We therefore sample the trajectories at time interval
$T_{p}$ to obtain stochastic wave-averaged trajectories. These are indicated
by the solid lines in the insets in figure 2 (and by all the lines in figure 2
itself).

Figure 2: Example particle trajectories in irregular waves, showing
experiments (light gray, and dark red) and Monte Carlo simulations of our
model (orange): (a) $H_{s}$ = 0.05 m, showing almost no breaking events and
trajectories evolving according to a Brownian motion, (b) $H_{s}$ = 0.17 m,
showing particles regularly surfing on a breaking wave, causing a jump in the
position. This is modelled by discrete jumps in the simulations. The insets
show a zoomed-in view of the dark red lines, where the red dashed lines show
an oscillating trajectory, corresponding to the effect of waves, corresponding
to the experimental data before wave averaging. All the other (solid) lines
show the wave-averaged position data (obtained from averaging using a period
corresponding to the peak period of the waves $T_{p}$).

####  3.2.3 Breaking detection

For the lowest significant wave height we have considered, breaking is
negligible, whereas for the highest, the particles have many encounters with
breaking crests, as indicated in table 1 by the jump arrival rate $\Lambda$,
which measures the number of jumps per unit time. Figure 2 illustrates the
evolution of particle position for both these extremes.

We classify particle motion as a ‘jump’, when the instantaneous horizontal
velocity (obtained from the trajectories before wave averaging) surpasses a
velocity threshold set as $u_{th}=0.3c$, where $c=\omega/k$ is the phase
velocity obtained from the linear dispersion relationship. Although this
threshold is arbitrary, lowering it results in a high number of jump events
for $H_{s}=0.05$ m, whilst breaking only rarely occurred for these
experiments. Physically, the threshold reflects the idea that particles during
breaking are transported with the crest at the phase velocity of the wave $c$
(they ‘surf’ the wave) instead of the much smaller Stokes drift velocity. The
distance covered at velocities higher than this threshold during one breaking
event is the jump amplitude $s$. See appendix B for an example of this
breaking detection procedure.

We use the jumps thus obtained to calibrate the jump process described by
(9)-(10). Figure 3a shows the jump arrival rate $\Lambda$ as a function of
steepness $\epsilon$ estimated from the experimental data for the four values
of steepness. Also shown is the sigmoid function for $\Lambda(\epsilon)$ given
by (9) with estimated values of the parameters in table 2. The variation in
jump amplitudes $s$ estimated from the experimental data in figure 3b is
captured well by the Gamma distribution (10). Finally, we estimate the
parameters of the Gamma distribution (10) as linear function of steepness:

$\alpha=a_{\alpha}+b_{\alpha}\epsilon,\quad\beta=a_{\beta}+b_{\beta}\epsilon,$
(21)

as shown in figure 3c with coefficients in table 2.

Figure 3: Calibration of the jump process for transport by breaking waves: (a) jump arrival rate $\Lambda$ as a function of steepness $\epsilon$, (b) probability density function of jump size $G(s)=\Gamma(s;\alpha,\beta)$ for $H_{s}$ = 0.09 (red), 0.12 (green) and 0.17 m (yellow) (the case of $H_{s}$ = 0.05 m is not shown, as no jumps are detected for this case), (c) scale and shape parameters $\alpha$ and $\beta$ for the probability density function of jump size as a function of $\epsilon$ as given by (21). $\Lambda(\epsilon)$ | $G(s;\alpha(\epsilon),\beta(\epsilon))$
---|---
$\tau_{\Lambda}$ [s] | $\phi_{\Lambda}$ | $\epsilon_{0,\Lambda}$ | $a_{\alpha}$ | $b_{\alpha}$ | $a_{\beta}$ | $b_{\beta}$
0.0677 | 503 | 0.128 | -6.49 | 3.46, | 0.305 | -0.003
Table 2: Calibrated values of the parameters describing the dependence of the
jump process for breaking waves on steepness $\epsilon$, distinguishing the
jump arrival rate $\Lambda(\epsilon)$ according to (9) and the jump size
probability density function $G(s)$ according to (10).

##  4 Results

In this section, we will compare Monte Carlo simulations of our model (2)
(using $10^{5}$ trajectories) to experiments. The Monte Carlo simulations
agree perfectly with the exact solutions for the first three moments
(17)-(19), so we will only show the former. We will examine non-breaking
(§4.1) and breaking waves (§4.2) in turn, followed by model predictions as a
function of steepness (§4.3).

###  4.1 Non-breaking waves: comparison between experiments and model
predictions

Figure 4: Comparison between model and experimental data for non-breaking
irregular waves ($H_{s}$ = 0.05 m): (a) normalized mean position $\langle
X(t)\rangle/\lambda_{p}$ (dashed-dotted) and $\pm 1$ standard deviation
(dashed) for experiments (red) and model simulations (orange), (b) normalized
mean position (dashed-dotted) for experiments (red), model simulations
(orange) and Stokes theory for non-breaking waves (dotted black), showing
proportionality to $t$. Normalized variance
$\langle|\tilde{X}(t)|^{2}\rangle/\lambda_{p}^{2}$ (dashed) for experiments
(red) and model simulations (orange) are approximately proportional to time
$t$. The normalized skewness
$\langle|\tilde{X}(t)|^{3}\rangle/\lambda_{p}^{3}$ (dotted) is negligible in
the experiments

Figure 4a shows the normalized mean $\langle X(t)\rangle/\lambda_{p}$ (dashed-
dotted), where $\lambda_{p}$ is the peak wavelength, for the experiments (red)
and the Monte Carlo simulations of (2) (orange) for the smallest steepness
waves ($H_{s}=$0.05 m, $\epsilon=0.074$). A negligible number of waves break
for this case such that $\Lambda\approx 0$, and the wave breaking term does
not contribute in (2). Therefore, the mean displacement by the Stokes drift
(dotted line) coincides with the mean drift predicted by the model, which in
turn is equal to that observed in the experiments (by definition here, as we
have used this agreement to estimate the Eulerian-mean flow, see §3.1.3). The
dashed lines show $\pm 1$ standard deviation.

More importantly, the normalized variance
$\langle|\tilde{X}(t)|^{2}\rangle/\lambda_{p}^{2}$, with
$\tilde{X}(t)=X(t)-\langle X(t)\rangle$, of the model simulations follows the
experiment closely in figure 4b. Indeed, extracting the power-law behavior in
figure 4b, the experimental particle position variance exhibits a linear $t$
dependence, in accordance with the single-particle Taylor diffusion prediction
by Herterich & Hasselmann (1982), which in turn agrees with our theoretical
prediction (18) when $\Lambda=0$ (no breaking). It is interesting to note that
a Wiener or normal diffusion process has a linear dependence on time, in
contrast to either sub- or superdiffusion, for which the conditions of the
central limit theorem are violated. Note that this result (the linear
dependence on time) is disagreement with Farazmand & Sapsis (2019), who
predict a superdiffusion $\langle|\tilde{X}(t)|^{2}\rangle\propto t^{4}$ for
$t>T_{p}$ based on the nonlinear John–Sclavounos equation. The particle
distribution we observe remains Gaussian with zero skewness (see figure 4), as
predicted by the analytical solutions for the third central moment in (19).

###  4.2 Breaking waves: comparison between experiments and model predictions

Figure 5: Comparison between model and experiment data for breaking irregular
waves ($H_{s}$ = 0.17 m): (a) normalized mean position $\langle
X(t)\rangle/\lambda_{p}$ (dashed-dotted) and $\pm 1$ standard deviation
(dashed lines) for experiments (red) and model simulations (orange), (b)
normalized mean position (dashed-dotted) for experiment (red), model
simulations (orange) and Stokes theory for non-breaking waves (dotted black),
showing proportionality to $t$. Normalized variance
$\langle|\tilde{X}(t)|^{2}\rangle/\lambda_{p}^{2}$ (dashed) for experiments
(red) and model simulations (orange) are approximately proportional to time
$t$. The normalized skewness
$\langle|\tilde{X}(t)|^{3}\rangle/\lambda_{p}^{3}$ (dotted) is finite for this
wave steepness with order-of-magnitude agreement between experiments (red) and
model simulations (orange).

Figure 5a shows the mean particle position (dashed-dotted) $\pm 1$ standard
deviation (dashed lines) in the experiments (red) and the Monte Carlo
simulations of (2) (orange) for the largest steepness waves ($H_{s}=0.17$ m
and $\epsilon=0.185$). For this significant wave height, many jumps in
position due to breaking occur, as illustrated earlier in figure 2b. Due to
the jump events, the mean drift (the slope of the line in figure 5a) is higher
than the theoretical prediction by the Stokes drift (dotted black line) and
instead follows that of (17) with $\Lambda\neq 0$. The power-law behavior in
figure 5b indicates that, in agreement with (18), the variance still scales
linearly with time. Here, the diffusion based on the Stokes drift alone (black
dashed line) underestimates the measured diffusion. Adding the effect of
breaking according to (18) gives good agreement with experimental data,
demonstrating that the contribution of breaking to particle diffusion can be
modelled effectively as a Poisson process. Finally, the finite skewness
predicted by the model (19) is also observed in experiments (orange and red
dotted lines in figure 5b, respectively).

###  4.3 Model predictions as a function of steepness

Figure 6: Model predictions with and without the jumps that transport by
breaking waves, showing (a) normalized mean velocity as a function of
characteristic steepness $\epsilon$, (b) variance of particle position, and
(c) skewness of particle position.

To quantify how the mean velocity and variance and skewness of particle
position are affected by the number of breaking encounters, we perform
simulations for the characteristic steepness range $\epsilon=[0.05,0.3]$
setting the Eulerian-mean flow to zero, as summarized in figure 6.

Figure 6a shows the mean drift $\langle u_{\rm L,mod}\rangle={\rm d}\langle
X(t)\rangle/{\rm d}t$ based on (2). If the breaking term is not taken into
account in (2) (i.e., $\Lambda=0$), this drift would be equal to the
theoretical Stokes drift value $\langle
u_{\text{S}}\rangle\propto\epsilon^{2}$ (dashed line). Taking into account the
jump-diffusion process of (2), $\langle u_{\rm L,mod}\rangle$ starts deviating
from $\langle u_{\text{S}}\rangle$ for steep waves in which many breaking
events occur. The nature of this deviation is given in (17), where $\langle
u_{\rm L,mod}\rangle$ is increased by a factor $(1+\alpha\Lambda/(\langle
u_{\text{S}}\rangle\beta))$. The second moment or variance deviates in a
similar fashion from the non-breaking case in figure 6b, deviating by a factor
$\left(1+\alpha(\alpha+1)\Lambda/(\sigma^{2}\beta)\right)$. The third moment
or skewness (figure 6c) becomes non-zero as the steepness increases, and the
particle position distribution only remains symmetric (zero skewness) for
small steepness.

Note that when studying drift induced by individual focused wave groups, Pizzo
et al. (2019) observed a sharp transition from a quadratic dependence of the
drift velocity on steepness below the breaking threshold to a linear
dependence above the breaking threshold. For individual waves or wave groups
there is a clear threshold for the steepness above which breaking occurs.
However, because we are considering many irregular waves only some of which
break, the number of breaking events and thus the mean drift velocity increase
continuously as a function of $\epsilon$ as shown in figure 3a.

##  5 Conclusions

In this paper we have developed a stochastic framework to describe particle
drift in irregular sea states using a jump-diffusion process to model the
enhanced drift due to breaking previously observed by Deike et al. (2017),
Pizzo et al. (2019) and Sinnis et al. (2021). The framework consisting of a
stochastic differential equation for particle transport and a corresponding
Fokker–Planck equation for the evolution of its probability distribution can
be used to predict mean drift and its higher-order statistical moments given
basic information describing the sea state, such as its spectrum or summary
parameters thereof (i.e., significant wave height and peak period). We compare
long-time laboratory experiments with a large number of particles with our
theoretical predictions and find good agreement, including specifically for
the contribution of the jump process to model enhanced transport by breaking
waves.

For an irregular wave field with negligible amount of breaking, we
experimentally verify that the variance of particle position or the single-
particle dispersion is proportional to time, confirming that the assumption of
normal diffusion is valid (in contrast to prediction by Farazmand & Sapsis
(2019)). Furthermore, we have described the evolution and quantified the
uncertainty of particle transport under the influence of wave breaking, by
modeling this as a compound Poisson process, where the amplitudes of the jumps
follow a Gamma distribution that is parameterized by the wave steepness. We
find that taking into account the jumps induced by breaking waves increases
the mean drift and the variance, and introduces a finite skewness into the
particle position distribution, whereas for the non-breaking diffusion
problem, the distribution remains normal and thus has zero skewness. Particle
tracking laboratory experiments corroborate this. Our results for the enhanced
mean drift due to breaking waves are in approximate quantitative agreement
with Pizzo et al. (2019), who show that the enhancement of drift due to
breaking can be up to 30% compared to the Stokes drift in certain cases.

In this paper we have considered long-crested (or unidirectional) waves. In
the ocean, sea state are almost always directionally spread (or short-
crested). We envisage our model can be readily extended to directional seas
using the results of Kenyon (1969) for non-breaking waves. Calibrating the
jump-diffusion process for directionally spread breaking waves will require
new laboratory experiments.

Looking forward, we believe that the simple stochastic framework we have
developed, which is based on stochastic differential equations used widely in
financial economics and climate physics, can be an important tool in
uncertainty quantification of prediction models in practical settings,
including search and rescue or salvage operations and pollution tracking and
clean-up efforts.

## Acknowledgements

We would like to thank Anton de Fockert and Wout Bakker at Deltares and Pieter
van der Gaag, Arie van der Vlies and Arno Doorn at Delft University of
Technology for their help setting up and conducting the experiment. This
research was performed in part through funding by the European Space Agency
(Grant no. 4000136626/21/NL/GLC/my). D.E. acknowledges financial support from
the Swiss National Science Foundation (P2GEP2-191480) and the ONR Grants
N00014-21-1-2357 and N00014-20-1-2366. T.S.vdB was supported by a Royal
Academy of Engineering Research Fellowship.

## Declaration of interests

The authors report no conflict of interest.

## References

  * Alsina et al. (2020) Alsina, José M., Jongedijk, Cleo E. & van Sebille, Erik 2020 Laboratory Measurements of the Wave‐Induced Motion of Plastic Particles: Influence of Wave Period, Plastic Size and Plastic Density. Journal of Geophysical Research: Oceans 125 (12).
  * Balk (2002) Balk, Alexander M. 2002 Anomalous behaviour of a passive tracer in wave turbulence. Journal of Fluid Mechanics 467, 163–203.
  * Balk (2006) Balk, A M 2006 Wave turbulent diffusion due to the Doppler shift. Journal of Statistical Mechanics: Theory and Experiment 2006 (08), P08018–P08018.
  * Banner et al. (2000) Banner, M. L., Babanin, A. V. & Young, I. R. 2000 Breaking probability for dominant waves on the sea surface. Journal of Physical Oceanography 30 (12), 3145–3160.
  * Bena et al. (2000) Bena, I., Copelli, M. & Van Den Broeck, C. 2000 Stokes’ drift: A rocking ratchet. Journal of Statistical Physics 101 (1-2), 415–424, arXiv: 9908338.
  * Breivik et al. (2014) Breivik, Øyvind, Janssen, Peter A.E.M. & Bidlot, Jean Raymond 2014 Approximate stokes drift profiles in deep water. Journal of Physical Oceanography 44 (9), 2433–2445.
  * van den Bremer & Breivik (2017) van den Bremer, T. S. & Breivik, Ø. 2017 Stokes drift. Phil. Trans. R. Soc. Lond. A 376, 20170104.
  * van den Bremer & Breivik (2018) van den Bremer, T. S. & Breivik, Ø 2018 Stokes drift. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 376 (2111), 20170104\.
  * van den Bremer & Taylor (2015) van den Bremer, T. S. & Taylor, P. H. 2015 Estimates of Lagrangian transport by wave groups: the effects of finite depth and directionality. J. Geophys. Res. 120 (4), 2701–2722.
  * van den Bremer et al. (2019) van den Bremer, T. S., Yassin, H. & Sutherland, B. R. 2019 Lagrangian transport by vertically confined internal gravity wavepackets. Journal of Fluid Mechanics 864, 348–380.
  * Bühler (2014) Bühler, O. 2014 Waves and Mean Flows, 2nd edn. Cambridge University Press, Cambridge, UK.
  * Bühler & Guo (2015) Bühler, Oliver & Guo, Yuan 2015 Particle dispersion by nonlinearly damped random waves. Journal of Fluid Mechanics 786, 332–347.
  * Bühler & Holmes-Cerfon (2009) Bühler, Oliver & Holmes-Cerfon, Miranda 2009 Particle dispersion by random waves in rotating shallow water. Journal of Fluid Mechanics 638, 5–26.
  * Calvert et al. (2021) Calvert, R., McAllister, M.L., Whittaker, C., Raby, A., Borthwick, A.G.L. & van den Bremer, T.S. 2021 A mechanism for the increased wave-induced drift of floating marine litter. Journal of Fluid Mechanics 915, A73.
  * Christensen et al. (2018) Christensen, Kai H., Øyvind Breivik, Dagestad, Knut-Frode, Röhrs, Johannes & Ward, Brian 2018 Short-term predictions of oceanic drift. Oceanography 31 (3), 59–67.
  * Cózar et al. (2014) Cózar, Andrés, Echevarría, Fidel, González-Gordillo, J. Ignacio, Irigoien, Xabier, Úbeda, Bárbara, Hernández-León, Santiago, Palma, Álvaro T., Navarro, Sandra, García-de Lomas, Juan, Ruiz, Andrea, Fernández-de Puelles, María L. & Duarte, Carlos M. 2014 Plastic debris in the open ocean. Proceedings of the National Academy of Sciences of the United States of America 111 (28), 10239–10244.
  * Cunningham et al. (2022) Cunningham, H. J., Higgins, C. & van den Bremer, T.S. 2022 The role of the unsteady surface wave-driven Ekman–Stokes flow in the accumulation of floating marine litter. J. Geophys. Res. Oceans 127 (6), e2021JC01810.
  * Dawson et al. (1993) Dawson, T. H., Kriebel, D. L. & Wallendorf, L. A. 1993 Breaking waves in laboratory-generated JONSWAP seas. Applied Ocean Research 15 (2), 85–93.
  * Deike et al. (2017) Deike, Luc, Pizzo, Nick & Melville, W. Kendall 2017 Lagrangian transport by breaking surface waves. Journal of Fluid Mechanics 829, 364–391.
  * Delandmeter & van Sebille (2019) Delandmeter, Philippe & van Sebille, Erik 2019 The Parcels v2.0 Lagrangian framework: new field interpolation schemes. Geoscientific Model Development 12 (8), 3571–3584.
  * Denisov et al. (2009) Denisov, S. I., Horsthemke, W. & Hänggi, P. 2009 Generalized Fokker-Planck equation: Derivation and exact solutions. European Physical Journal B 68 (4), 567–575, arXiv: 0808.0274.
  * DiBenedetto et al. (2022) DiBenedetto, M. H., Clark, L. K. & Pujara, N. 2022 Enhanced settling and dispersion of inertial particles in surface waves. Journal of Fluid Mechanics 936, A38.
  * Dobler et al. (2019) Dobler, Delphine, Huck, Thierry, Maes, Christophe, Grima, Nicolas, Blanke, Bruno, Martinez, Elodie & Ardhuin, Fabrice 2019 Large impact of Stokes drift on the fate of surface floating debris in the South Indian Basin. Marine Pollution Bulletin 148, 202–209.
  * Eriksen et al. (2014) Eriksen, Marcus, Lebreton, Laurent C.M., Carson, Henry S., Thiel, Martin, Moore, Charles J., Borerro, Jose C., Galgani, Francois, Ryan, Peter G. & Reisser, Julia 2014 Plastic Pollution in the World’s Oceans: More than 5 Trillion Plastic Pieces Weighing over 250,000 Tons Afloat at Sea. PLoS ONE 9 (12), 1–15.
  * Farazmand & Sapsis (2019) Farazmand, Mohammad & Sapsis, Themistoklis 2019 Surface Waves Enhance Particle Dispersion. Fluids 4 (1), 55, arXiv: 1902.04034.
  * Gardiner (1983) Gardiner, C. W. 1983 Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences. Springer Verlag, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/bbpc.19850890629.
  * Gaviraghi (2017) Gaviraghi, Beatrice 2017 Theoretical and numerical analysis of Fokker-Planck optimal control problems for jump-diffusion processes. PhD thesis, Julius-Maximilians-Universität Würzburg.
  * Gaviraghi et al. (2017) Gaviraghi, B., Annunziato, M. & Borzì, A. 2017 Analysis of splitting methods for solving a partial integro-differential fokker–planck equation 294, 1–17.
  * Herterich & Hasselmann (1982) Herterich, K. & Hasselmann, K. 1982 The horizontal diffusion of tracers by surface waves. J. Phys. Oceanogr. 12 (7 , Jul. 1982), 704–711.
  * Higgins et al. (2020a) Higgins, C., Vanneste, J. & Bremer, T. S. 2020a Unsteady Ekman‐Stokes Dynamics: Implications for Surface Wave‐Induced Drift of Floating Marine Litter. Geophysical Research Letters 47 (18).
  * Higgins et al. (2020b) Higgins, C., van den Bremer, T. & Vanneste, J. 2020b Lagrangian transport by deep-water surface gravity wavepackets: Effects of directional spreading and stratification. Journal of Fluid Mechanics 883, A42.
  * Higham (2001) Higham, D. J. 2001 An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Review 43 (3), 525–546.
  * Holthuijsen & Herbers (1986) Holthuijsen, L. H. & Herbers, T. H. C. 1986 Statistics of Breaking Waves Observed as Whitecaps in the Open Sea. Journal of Physical Oceanography 16 (2), 290–297.
  * Huang et al. (2016) Huang, Guoxing, Huang, Zhen Hua & Law, Adrian W. K. 2016 Analytical Study on Drift of Small Floating Objects under Regular Waves. Journal of Engineering Mechanics 142 (6), 06016002\.
  * Jambeck et al. (2015) Jambeck, Jenna R., Geyer, Roland, Wilcox, Chris, Siegler, Theodore R., Perryman, Miriam, Andrady, Anthony, Narayan, Ramani & Law, Kara Lavender 2015 Plastic waste inputs from land into the ocean. Science 347 (6223), 768–771.
  * Jansons & Lythe (1998) Jansons, Kalvis M. & Lythe, G. D. 1998 Stochastic Stokes Drift. Physical Review Letters 81 (15), 3136–3139, arXiv: 9808042.
  * Kenyon (1969) Kenyon, Kern 1969 Stokes Drift for Random Gravity Waves. J Geophys Res 74 (28), 6991–6994.
  * Lenain & Pizzo (2020) Lenain, Luc & Pizzo, Nick 2020 The contribution of high-frequency wind-generated surface waves to the Stokes drift. J. Phys. Oceanogr. 50 (12), 3455 – 3465.
  * Lenain et al. (2019) Lenain, Luc, Pizzo, Nick & Melville, W. Kendall 2019 Laboratory studies of Lagrangian transport by breaking surface waves. Journal of Fluid Mechanics 876, 1–12.
  * Li & Li (2021) Li, Yan & Li, Xin 2021 Weakly nonlinear broadband and multi-directional surface waves on an arbitrary depth: A framework, Stokes drift, and particle trajectories. Physics of Fluids 33 (7).
  * Longuet-Higgins (1952) Longuet-Higgins, M. S. 1952 On the Statistical Distribution of the Heights of Sea Waves. Journal of Marine Research 11, 245–266.
  * Longuet-Higgins (1953) Longuet-Higgins, M. S. 1953 Mass transport in water waves. Philos. T. R. Soc. A. 245 (903), 535–581.
  * Longuet-Higgins (1983) Longuet-Higgins, M. S. 1983 on the Joint Distribution of Wave Periods and Amplitudes in a Random Wave Field. Proceedings of The Royal Society of London, Series A: Mathematical and Physical Sciences 389 (1797), 241–258.
  * Milstein (1995) Milstein, G. N. 1995 Numerical Integration of Stochastic Differential Equations. Springer Netherlands.
  * Monismith (2020) Monismith, S. G. 2020 Stokes drift: Theory and experiments. J. Fluid Mech. 884, F1.
  * Myrhaug et al. (2014) Myrhaug, Dag, Wang, Hong & Holmedal, Lars Erik 2014 Stokes drift estimation for deep water waves based on short-term variation of wave conditions. Coastal Engineering 88, 27–32.
  * Ochi & Tsai (1983) Ochi, Michel K. & Tsai, Cheng-Han 1983 Prediction of Occurrence of Breaking Waves in Deep Water. Journal of Physical Oceanography 13 (11), 2008–2019.
  * Onink et al. (2019) Onink, Victor, Wichmann, David, Delandmeter, Philippe & Sebille, Erik 2019 The Role of Ekman Currents, Geostrophy, and Stokes Drift in the Accumulation of Floating Microplastic. Journal of Geophysical Research: Oceans 124 (3), 1474–1490.
  * Phillips (1985) Phillips, O. M. 1985 Spectral and statistical properties of the equilibrium range in wind-generated gravity waves. Journal of Fluid Mechanics 156, 505–531.
  * Pizzo et al. (2019) Pizzo, Nick, Melville, W. Kendall & Deike, Luc 2019 Lagrangian Transport by Nonbreaking and Breaking Deep-Water Waves at the Ocean Surface. Journal of Physical Oceanography 49 (4), 983–992.
  * Pizzo (2017) Pizzo, Nick E. 2017 Surfing surface gravity waves. Journal of Fluid Mechanics 823, 316–328.
  * Restrepo & Ramirez (2019) Restrepo, Juan M. & Ramirez, Jorge M. 2019 Transport due to Transient Progressive Waves. Journal of Physical Oceanography 49 (9), 2323–2336.
  * Sainte-Rose et al. (2020) Sainte-Rose, Bruno, Wrenger, Hendrik, Limburg, Hans, Fourny, Arthur & Tjallema, Arjen 2020 Monitoring and Performance Evaluation of Plastic Cleanup Systems: Part I — Description of the Experimental Campaign. International Conference on Offshore Mechanics and Arctic Engineering Volume 6B: Ocean Engineering.
  * Sanderson & Okubo (1988) Sanderson, Brian G. & Okubo, Akira 1988 Diffusion by internal waves. Journal of Geophysical Research 93 (C4), 3570\.
  * Santamaria et al. (2013) Santamaria, F., Boffetta, F., Martins Afonso, M., Mazzino, A., Onorato, M. & Pugliese, D. 2013 Stokes drift for inertial particles transported by water waves. Europhys. Lett. 102 (1), 14003.
  * van Sebille et al. (2020) van Sebille, Erik, Aliani, Stefano, Law, Kara Lavender, Maximenko, Nikolai, Alsina, José M, Bagaev, Andrei, Bergmann, Melanie, Chapron, Bertrand, Chubarenko, Irina, Cózar, Andrés, Delandmeter, Philippe, Egger, Matthias, Fox-Kemper, Baylor, Garaba, Shungudzemwoyo P, Goddijn-Murphy, Lonneke, Hardesty, Britta Denise, Hoffman, Matthew J, Isobe, Atsuhiko, Jongedijk, Cleo E, Kaandorp, Mikael L A, Khatmullina, Liliya, Koelmans, Albert A, Kukulka, Tobias, Laufkötter, Charlotte, Lebreton, Laurent, Lobelle, Delphine, Maes, Christophe, Martinez-Vicente, Victor, Morales Maqueda, Miguel Angel, Poulain-Zarcos, Marie, Rodríguez, Ernesto, Ryan, Peter G, Shanks, Alan L, Shim, Won Joon, Suaria, Giuseppe, Thiel, Martin, van den Bremer, Ton S & Wichmann, David 2020 The physical oceanography of the transport of floating marine debris. Environmental Research Letters 15 (2), 023003.
  * van Sebille et al. (2015) van Sebille, Erik, Wilcox, Chris, Lebreton, Laurent, Maximenko, Nikolai, Hardesty, Britta Denise, van Franeker, Jan A, Eriksen, Marcus, Siegel, David, Galgani, Francois & Law, Kara Lavender 2015 A global inventory of small floating plastic debris. Environmental Research Letters 10 (12), 124006\.
  * Sinnis et al. (2021) Sinnis, J. T., Grare, L., Lenain, L. & Pizzo, N. 2021 Laboratory studies of the role of bandwidth in surface transport and energy dissipation of deep-water breaking waves. J. Fluid Mech. 927, A5.
  * Staneva et al. (2021) Staneva, Joanna, Ricker, Marcel, Carrasco Alvarez, Ruben, Breivik, Øyvind & Schrum, Corinna 2021 Effects of wave-induced processes in a coupled wave–ocean model on particle transport simulations. Water 13 (4).
  * Stokes (1847) Stokes, G. G. 1847 On the theory of oscillatory waves. Trans. Camb. Phil. Soc. 8, 411–455.
  * Sullivan et al. (2007) Sullivan, Peter P., McWilliams, James C. & Melville, W. Kendall 2007 Surface gravity wave effects in the oceanic boundary layer: Large-eddy simulation with vortex force and stochastic breakers. Journal of Fluid Mechanics 593, 405–452.
  * Taylor (1922) Taylor, G. I. 1922 Diffusion by Continuous Movements. Proceedings of the London Mathematical Society s2-20 (1), 196–212.
  * Webb & Fox-Kemper (2011) Webb, A. & Fox-Kemper, B. 2011 Wave spectral moments and Stokes drift estimation. Ocean Modelling 40 (3-4), 273–288.
  * Weichman & Glazman (2000) Weichman, Peter B. & Glazman, Roman E. 2000 Passive scalar transport by travelling wave fields. Journal of Fluid Mechanics 420, 147–200.

##  Appendix A Particle tracking

A downward-looking camera mounted on a walkway above the basin was used to
track particles as they moved across the basin. The camera intriniscis were
found by detecting multiple images of a checkerboard in different orientations
across the entire field of view (FOV). Figure A.1(a) shows a top view of the
Atlantic Basin with more or less randomly distributed yellow plastic spheres
floating on the surface. The spheres are identified and tracked using OpenCV
in Python. The image is first filtered using a Hue filter, which produces
figure A.1(b). A Correlation Filter with Channel and Spatial Reliability
(CSRT) algorithm was used to track the objects between frames to create
trajectories in sub-pixel locations and time. The trajectories were then
transformed to the still-water plane, defined by an image of a floating
checkerboard, and tank coordinates ($x,y$) by detecting and inverting the
camera intrinsics and applying a measured translation from camera field of
view.

The trackers and trajectories were post-processed to eliminate any erroneous
tracking events, such as spheres colliding, loss of tracking, or jumps of the
particle tracked by the algorithm to a nearby particle, which sometimes
occurred when particles were lost momentarily under breaking waves. Finally,
the trajectories were manually inspected for quality control.

Figure A.1: (a) Footage of the Atlantic Basin from the downward-looking
camera. The small yellow spheres can be seen floating on the surface. (b)
Identification of the spheres after application of the Hue saturation filter
in OpenCV.

##  Appendix B Jump detection

Figure B.1 displays the steps of the jump detection and jump amplitude
(distance travelled in a jump) estimation process for two example
trajectories. For $H_{s}$ = 0.05 m, panel a shows the derivative of particle
position over a time-step $\delta t$ determined by the camera frame-rate:
$u_{\text{I}}=\delta X/\delta t$; we consider this to be the instantaneous
velocity (before wave averaging). The dashed line indicates the threshold
value $u_{th}=0.3c$, where $c=\omega/k$ is the phase velocity obtained from
the linear dispersion relationship. In panels c and d the Heaviside functions
$\mathcal{H}(u_{\text{I}}(t)-u_{th})$ mark the time segments where the
particle is classified as ‘jumping’. In panel e and f, based on the Heaviside
functions in panels c and d, the jump amplitude $s$ (the distance covered
during a jump) can be estimated.

Figure B.1: Jump detection and amplitude estimation process for (a,c,e) the
least steep waves with $H_{s}$ = 0.05 m and (b,d,f) the steepest waves with
$H_{s}$ = 0.17 m. (a,b) Instantaneous particle velocity (before wave
averaging) obtained from the derivative of particle position at the camera
frame-rate $u_{I}=\delta X/\delta t$. The dashed line indicates the threshold
value $u_{th}=0.3c$, with $c$ the phase velocity. (c,d) The Heaviside function
$\mathcal{H}(u_{\text{I}}(t)-u_{th})$ marks the time segments where the
particle is ‘jumping’. (e,f) Particle position $X(t)$ at the camera frame
rate, where the red crosses mark the positions for which the Heaviside
function is positive and behaviour is classified as ‘jumping’.

##  Appendix C Drift velocities

Figure C.1: Drift velocity estimations and measurements as a function of
steepness. Stokes drift based on a JONSWAP fit on the spectrum truncated at
$\langle u_{\text{S},f}\rangle$ ($\vartriangle$) and its quadratic fit (dashed
line), measured particle drift $\langle u_{\text{L,exp}}\rangle$ (thick
$\times$), measured particle drift corrected by model mean flow $\langle
u_{\text{L,exp}}\rangle-\langle u_{\text{E,fit}}\rangle$ (thin $\times$)

Figure C.1 shows various drift velocity estimations. For our experiments,the
Stokes drift is based on a JONSWAP spectrum fitted to the measured spectrum by
fitting a JONSWAR spectrum to the experimental spectrum $\langle
u_{\text{S},f}\rangle$, indicated by the blue triangles, with the dashed blue
line its quadratic fit.

We correct the experimentally measured mean drift velocity $\langle
u_{\text{L,exp}}\rangle$ by the mean flow in the model $\langle
u_{\text{E,fit}}\rangle$, resulting in the thin red crosses ($\times$). The
difference between these and dashed blue line, i.e. the Stokes drift, shows
that for low characteristic steepness this coincides with the stokes drift,
whereas for high steepness, the stokes drift underestimates the mean drift
velocity by a certain fraction. The difference is attributed to the jump
process in the model.

$H_{s}$ | $\langle u_{\text{E, exp}}\rangle$ | $\langle u_{\text{E, fit}}\rangle$
---|---|---
(m) | (mm/s) | (mm/s)
0.050 | -2 | +4.1
0.090 | -10 | +0.8
0.120 | -6 | -5.6
0.170 | -17 | -10.8
Table 3: Difference between the experimentally measured Eularian mean flow,
and its value used in the model.

##  Appendix D Velocity spectrum particles

Figure D.1 compares the velocity spectrum of the particle trajectories (dark
red) to that of the first order velocity spectrum calculated from the surface
elevation (pink). The later has much spectral power for higher frequencies
that is not present in the former. Therefore, if this first order velocity
does not contribute to the particle movement at these frequencies, the higher
order velocity (the Stokes drift) cannot contribute either.

Figure D.1: $H_{s}$ = 0.05 m. Spectrum of the velocity of the particles (red),
dashed lines indicate $\pm 1$ std. The linear velocity calculated from the
spectrum of the surface elevation $\eta$ is shown in pink. The gray dashed
line indicates $f=3f_{0}$. (a) Linear scale (b) Log scale.

Figure D.2: $H_{s}$ = 0.09 m. Lines same as in fig. D.1

Figure D.3: $H_{s}$ = 0.12 m. Lines same as in fig. D.1

Figure D.4: $H_{s}$ = 0.17 m. Lines same as in fig. D.1