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. Reference Eriksen, Lebreton, Carson, Thiel, Moore, Borerro, Galgani, Ryan and Reisser2014). A mismatch exists between the estimated amount of land-generated plastic entering coastal waters (
$5$–
$12$ million tonnes yr
$^{-1}$, Jambeck et al. Reference Jambeck, Geyer, Wilcox, Siegler, Perryman, Andrady, Narayan and Law2015) and the estimated total amount of plastic floating at sea (less than 
$0.3$ million tonnes, Cózar et al. Reference Cózar2014; Eriksen et al. Reference Eriksen, Lebreton, Carson, Thiel, Moore, Borerro, Galgani, Ryan and Reisser2014; van Sebille et al. Reference van Sebille, Wilcox, Lebreton, Maximenko, Hardesty, van Franeker, Eriksen, Siegel, Galgani and Law2015), which necessitates more accurate transport and dispersion models (van Sebille et al. Reference van Sebille2020) as well as better modelling of vertical concentrations (DiBenedetto et al. Reference DiBenedetto, Donohue, Tremblay, Edson and Law2023). 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. Reference Sainte-Rose, Wrenger, Limburg, Fourny and Tjallema2020). Similarly, estimating environmental impact for oil spills and open-sea rescue operations relies crucially on transport models (e.g. Christensen et al. Reference Christensen, Breivik, Dagestad, Röhrs and Ward2018).
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 Reference Stokes1847). Wave models such as WAM (The WAMDI Group 1988) and WaveWatch III (Tolman et al. Reference Tolman2009) predict wave field properties averaged over longer time scales (typically 3 hours), based on which estimates of the mean Stokes drift over that period can be made (Webb & Fox-Kemper Reference Webb and Fox-Kemper2011; Breivik, Janssen & Bidlot Reference Breivik, Janssen and Bidlot2014). 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. Reference Dobler, Huck, Maes, Grima, Blanke, Martinez and Ardhuin2019) and pushing microplastics closer to the coast (Delandmeter & van Sebille Reference Delandmeter and van Sebille2019; Onink et al. Reference Onink, Wichmann, Delandmeter and Sebille2019).
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 of 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 Reference Longuet-Higgins1953). 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, Vanneste & Bremer Reference Higgins, Vanneste and Bremer2020a), or wave and ocean circulation models need to be properly coupled (e.g. Staneva et al. Reference Staneva, Ricker, Carrasco Alvarez, Breivik and Schrum2021). The Ekman–Stokes flow can have significant consequences for global floating marine litter accumulation patterns (Cunningham, Higgins & van den Bremer Reference Cunningham, Higgins and van den Bremer2022). Wave-induced Eulerian-mean currents can also play a role in laboratory experiments (van den Bremer & Breivik Reference van den Bremer and Breivik2018; van den Bremer, Yassin & Sutherland Reference van den Bremer, Yassin and Sutherland2019). 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. Reference Santamaria, Boffetta, Martins Afonso, Mazzino, Onorato and Pugliese2013; Huang, Huang & Law Reference Huang, Huang and Law2016; Alsina, Jongedijk & van Sebille Reference Alsina, Jongedijk and van Sebille2020; Calvert et al. Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021; DiBenedetto, Clark & Pujara Reference DiBenedetto, Clark and Pujara2022). Third, breaking waves are known to transport particles much faster than non-breaking waves. For steep waves, particles may surf on the wave (Pizzo Reference Pizzo2017) and be subject to transport faster than the Stokes drift due to wave breaking (Deike, Pizzo & Melville Reference Deike, Pizzo and Melville2017; Pizzo, Melville & Deike Reference Pizzo, Melville and Deike2019). 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 for various sea states (e.g. Webb & Fox-Kemper Reference Webb and Fox-Kemper2011; Breivik et al. Reference Breivik, Janssen and Bidlot2014) with an important but not always resolved role for the spectral tail (Lenain & Pizzo Reference Lenain and Pizzo2020). 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 Reference Herterich and Hasselmann1982). Estimates of the variance of the Stokes drift can either be obtained from the spectrum (Herterich & Hasselmann Reference Herterich and Hasselmann1982) or from the joint distribution of the significant wave height and peak period (Longuet-Higgins Reference Longuet-Higgins1983; Myrhaug, Wang & Holmedal Reference Myrhaug, Wang and Holmedal2014). Alternatively, deterministic simulations of the particle trajectories can be performed (e.g. Farazmand & Sapsis Reference Farazmand and Sapsis2019; Li & Li Reference Li and Li2021), 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. Reference Deike, Pizzo and Melville2017), in contrast to the quadratic scaling of Stokes drift with steepness for non-breaking waves. Restrepo & Ramirez (Reference Restrepo and Ramirez2019) 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, Pizzo & Melville (Reference Lenain, Pizzo and Melville2019) and Sinnis et al. (Reference Sinnis, Grare, Lenain and Pizzo2021), where Sinnis et al. (Reference Sinnis, Grare, Lenain and Pizzo2021) also consider the effects of bandwidth. Considering much larger scales, relevant for application to the real ocean, requires a stochastic approach. Pizzo et al. (Reference Pizzo, Melville and Deike2019) extend the result obtained by Deike et al. (Reference Deike, Pizzo and Melville2017) for wave groups to sea states. Based on the wave spectrum (peak wavenumber) and wind speed, one can estimate the breaking statistic 
$\varLambda (c)\,{\rm d}c$ (Ochi & Tsai Reference Ochi and Tsai1983; Holthuijsen & Herbers Reference Holthuijsen and Herbers1986; Dawson, Kriebel & Wallendorf Reference Dawson, Kriebel and Wallendorf1993; Banner, Babanin & Young Reference Banner, Babanin and Young2000; Sullivan, McWilliams & Melville Reference Sullivan, McWilliams and Melville2007), 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 Reference Phillips1985). Consequently, the breaking drift speed found in Pizzo et al. (Reference Pizzo, Melville and Deike2019) 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 with 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. Reference Pizzo, Melville and Deike2019).
A series of papers have examined stochastic Stokes drift (Jansons & Lythe Reference Jansons and Lythe1998; Bena, Copelli & Van Den Broeck Reference Bena, Copelli and Van Den Broeck2000), 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 Reference Herterich and Hasselmann1982; Sanderson & Okubo Reference Sanderson and Okubo1988; Weichman & Glazman Reference Weichman and Glazman2000; Balk Reference Balk2002, Reference Balk2006). Bühler & Holmes-Cerfon (Reference Bühler and Holmes-Cerfon2009) 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 (Reference Bühler and Guo2015), 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 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. That is, given a wave spectrum with random phases, we examine if the drift-diffusion behaviour of particle position can be described by a SDE with a reasonable parameterization, which varies continuously with average steepness. We demonstrate that a stochastic framework, as opposed to a deterministic framework, is necessary, and that the model we propose consisting of a Brownian and jump process is a feasible option.
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 behaviour 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 are 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 Reference Bühler2014). In our case, we consider the Lagrangian-mean velocity
\begin{equation} u_{L} = u_{E,NB} + \underbrace{u_{S} + u_{B}}_{u_{D}},\end{equation}
which consists of the sum of the Eulerian-mean velocity that excludes the effects of wave breaking 
$u_{E,NB}$ and a drift velocity 
$u_{D}$, which in turn consists of the Stokes drift for non-breaking waves 
$u_{S}$ and a breaking contribution 
$u_{B}$. For simplicity, we subsequently ignore the wave-induced Eulerian-mean velocity 
$u_{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 Reference van den Bremer and Taylor2015; Higgins, van den Bremer & Vanneste Reference Higgins, van den Bremer and Vanneste2020a). To compare with 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 SDE can be written as
\begin{align} {\rm d}\,X &= \underbrace{b(X(t),t)\, {\rm d}t}_{{Mean\ drift}} + \underbrace{\sigma (X(t),t)\, {\rm d}W(t) }_{{Diffusion\ process}}+ \underbrace{ {\rm d}J(t), }_{{Compound\ jump\ process}} \end{align}
\begin{align} &= \langle u_{S}\rangle \,{\rm d}t + \sigma \,{\rm d}W(t) + {\rm d}J(\varLambda(\epsilon),G(s;\epsilon);t), \end{align}
where 
$b(X(t))=\langle u_{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_{B}\rangle \neq 0$, reflecting the contribution of breaking to mean drift.
Taking the terms in (2.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.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 Reference Stokes1847)
\begin{equation} U_{S}(z) = a^2 \omega k\, {\rm e}^{2k z}, \end{equation}
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_{S}$ is consistently approximated as 
$u_{S}=U_{S}(z=0)$. Written in terms of the (commonly estimated) wave period 
$T$ and wave height for periodic linear waves 
$H=2a$
\begin{equation} u_{S} = (ak)^2 c = \frac{a^2 \omega^3}{g} = \frac{(2 {\rm \pi})^3}{g}\frac{H^2}{T^3}, \end{equation}
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 (2.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 Reference Kenyon1969; Webb & Fox-Kemper Reference Webb and Fox-Kemper2011)
\begin{equation} \langle u_{S} \rangle = \frac{16{\rm \pi}^3}{g} \int_0^\infty S(f) f^3 \,{\rm d}f, \end{equation}
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, (2.6) gives the same result as (2.5).
The second term in (2.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}_{S} = u_{S}-\langle u_{S}\rangle$) do not depend on time. One can write 
$D=\langle \tilde {u}_{S}^2 \rangle \tau = \sigma _{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 _{S}^2=\langle \tilde {u}_{S}^2 \rangle$ the variance of the Stokes drift (Taylor Reference Taylor1922; Herterich & Hasselmann Reference Herterich and Hasselmann1982; Farazmand & Sapsis Reference Farazmand and Sapsis2019). Therefore, 
$\sigma = \sqrt {2D}=\sqrt {2 \tau }\sigma _{S}$ in (2.2).
To estimate the diffusion coefficient 
$D$, it is therefore necessary to obtain a distribution for the Stokes drift 
$P(u_{S})$ and derive from this its standard deviation 
$\sigma _{S}$. In Herterich & Hasselmann (Reference Herterich and Hasselmann1982) 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 (2.5) (Longuet-Higgins Reference Longuet-Higgins1952). Since 
$u_{S} \propto H_s^2$, the Stokes drift follows an exponential distribution. Specifically, the probability density function for 
$u_{S}$ reads
\begin{equation} P(u_{S};H_{s}, T_p) = \frac{g T_p^3}{4 {\rm \pi}^3} \frac{1}{H_{s}^2} \exp{\left[-\frac{gT_p^3}{ 4{\rm \pi}^3}\frac{u_{S}}{H_{s}^2}\right]}. \end{equation}
Alternatively, Myrhaug et al. (Reference Myrhaug, Wang and Holmedal2014) and Longuet-Higgins (Reference Longuet-Higgins1983) 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 _{S}=\langle u_{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, a significant increase in forward velocity can occur on the front of a breaking wave; see Pizzo (Reference Pizzo2017) and Deike et al. (Reference Deike, Pizzo and Melville2017) for detailed discussions of the dynamics of this surfing mechanism. We refer to this rapid and finite-time increase in forward velocity as a ‘jump’ in particle position. 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.2)
\begin{equation} J(t) = \sum_{i=0}^{N_\varLambda(t)}s_i, \end{equation}
where 
$N_\varLambda (t)$ is a counting of a Poisson point process, and the arrival rate 
$\varLambda$ corresponds to the expected number of jumps per unit time for the compound Poisson process.
We expect the arrival rate 
$\varLambda$ 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)
\begin{equation} \varLambda(\epsilon; \tau_\varLambda, \phi_\varLambda, \epsilon_{0,\varLambda}) = \frac{1/ \tau_\varLambda}{1+\exp{[-\phi_\varLambda(\epsilon-\epsilon_{0,\varLambda})}]}, \end{equation}
where the parameters 
$\tau _\varLambda$, 
$\phi _\varLambda$ and 
$\epsilon _{0,\varLambda }$ 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)
\begin{equation} G(s;\alpha,\beta) = \frac{\beta^\alpha}{\varGamma(\alpha)}s^{\alpha-1}\, {\rm e}^{-\beta s}, \end{equation}
where 
$\alpha >0$ is the shape parameter, 
$\beta >0$ the rate parameter and 
$\varGamma (\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.2) in order to obtain a probability distribution for particle position 
$X(t)$. Monte Carlo simulations of (2.2) can be numerically integrated (e.g. Milstein Reference Milstein1995; Higham Reference Higham2001), 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.2), the so-called Fokker–Planck equation, can be solved directly. For (2.2), the Fokker–Planck equation is given by (e.g. Gardiner Reference Gardiner1983; Denisov, Horsthemke & Hänggi Reference Denisov, Horsthemke and Hänggi2009; Gaviraghi, Annunziato & Borzì Reference Gaviraghi, Annunziato and Borzì2017)
\begin{equation} \frac{\partial}{\partial t}P(x,t) ={-}\langle u_{S} \rangle \frac{\partial}{\partial x}P(x,t) + \frac{\sigma^2}{2} \frac{\partial^2}{\partial x^2}P(x,t) - \varLambda P(x,t) + \varLambda \int_{-\infty}^{\infty}\,{\rm d}\kern0.06em x' P(x',t) J(x-x'). \end{equation}This partial differential equation can be solved numerically (e.g. Gaviraghi Reference Gaviraghi2017). Alternatively, by employing the characteristic functions of the distribution
\begin{equation} \hat{P}(l) = \int_{-\infty}^{\infty}\,{\rm d}\kern0.06em x\, {\rm e}^{{\rm i}lx} P(x), \end{equation}we obtain the ordinary differential equation
\begin{equation} \frac{\partial}{\partial t} \hat{P}(l) = i \langle u_{S} \rangle l \hat{P}(l) - \frac{\sigma^2}{2}l^2 \hat{P}(l) - \varLambda \hat{P}(l) + \varLambda \hat{G} \hat{P}(l), \end{equation}where the characteristic function of the gamma distribution is
\begin{equation} \hat{G}(l) = \left(1-i \frac{l}{\beta}\right)^{-\alpha}. \end{equation}Equation (2.13) can be solved exactly, with solution
\begin{equation} \hat{P}(l,t) = A(l) \exp{\left[\left(i \langle u_{S} \rangle l - \frac{\sigma^2}{2}l^2 - \varLambda + \varLambda \left(1-i \frac{l}{\beta}\right)^{-\alpha}\right)t\right]},\end{equation}
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
\begin{equation} m_n = \langle X^n\rangle = \int_{-\infty}^{\infty}\,{\rm d}\kern0.06em x x^n\, {\rm e}^{{\rm i}lx} P(x) =\left. ({-}i)^n \frac{\partial^n \hat{P}(l)}{\partial \,{\rm d}l^n}\right\vert_{l=0}. \end{equation}
For a Dirac delta function as the initial condition (i.e. 
$\hat {P}(t=0)=1$), we have 
$A=1$ in (2.15) and we can obtain exact solutions for the first (three) central moments of (2.11)
\begin{gather} \langle X(t) \rangle = \left.-{\rm i} \frac{\partial \hat{P}(l,t)}{\partial l} \right\vert_{l=0} = \left(\langle u_{S} \rangle + \frac{\alpha}{\beta}\varLambda\right)t, \end{gather}
\begin{gather}\langle \tilde{X}(t)^2 \rangle = \left.- \frac{\partial^2 \hat{P}(l,t)}{\partial l^2}\right\vert_{l=0} - m_1^2 = \left(\sigma^2 + \frac{\alpha (\alpha +1)}{\beta}\varLambda \right) t, \end{gather}
\begin{gather}\langle \tilde{X}(t)^3 \rangle = \left.i \frac{\partial^3 \hat{P}(l,t)}{\partial l^3} \right\vert_{l=0} - 3 m_1 m_2 + 2 m_1^3 = \frac{\alpha(1+\alpha)(2+\alpha)}{\beta^2}\varLambda t, \end{gather}
where 
$\tilde {X}(t)=X(t)-\langle X(t) \rangle$. In (2.17)–(2.19), the expected Stokes drift for non-breaking waves 
$\langle u_{S} \rangle (\epsilon )$, the arrival rate 
$\varLambda (\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 (2.17)–(2.19) that each central moment is linearly increasing with time. Breaking increases both the mean drift (cf. (2.17)) and the variance of particle position (cf. (2.18)); it also introduces a non-zero (positive) skewness not predicted for non-breaking waves (cf. (2.19)).
3. Laboratory experiments
3.1. Laboratory set-up and input conditions
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 Joint North Sea Wave Project (JONSWAP) spectral density 
$S(\omega )$
\begin{equation} S(\omega) = \frac{K g^2}{\omega^5}\exp\left[-\frac{5}{4}\left(\frac{\omega_p}{\omega} \right)^4\right]\gamma^r, \end{equation}
with 
$\omega$ the frequency of the waves (in rad s
$^{-1}$), 
$g$ the gravitational acceleration, 
$\omega _p$ the peak frequency and 
$r=\exp [-(\omega -\omega _p)^2/(2\sigma _{S}^2\omega _p^2)]$ with (non-dimensional) spectral width 
$\sigma _{S}= 0.07$ for 
$\omega \leq \omega _p$ and 
$\sigma _{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.2) (i.e. 
$\sigma = \sqrt {2D}=\sqrt {2 \tau }\sigma _{S}$). Using a second moment of the spectrum, we calculate the spectral width as 
$\Delta \omega = 1.39$ rad s
$^{-1}$ and set 
$\tau =1/\Delta \omega$. The prefactor 
$K$ in (3.1) is adjusted to obtain the desired significant wave height 
$H_s$. Phases of a discretized 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 %.
Figure 1. Experimental set-up, indicating the overhead particle-tracking camera's field of view (FOV), wave gauges (red dots) and velocity measurement probes (blue dots).
We note that the JONSWAP spectrum we use, which describes fetch-limited seas, is not representative of all realistic ocean conditions. In particular, the JONSWAP spectrum does not capture well the shape of the spectrum of seas that are not fetch limited. In particular, the transition from the equilibrium to the saturation range for high wavenumbers is not described well, for which parameterizations based on wind velocity have been proposed (Phillips Reference Phillips1985; Romero & Melville Reference Romero and Melville2010; Lenain & Melville Reference Lenain and Melville2017) and which could have significant consequences for the magnitude of the Stokes drift (Lenain & Pizzo Reference Lenain and Pizzo2020). However, as the laboratory (especially without directional spreading and wind generation) is not the right environment to study these effects, the JONSWAP spectrum suffices for our purposes.
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 velocity measurement probes) was used to calculate the values of the parameters measured in experiments reported in table 1, where they arelabelled with the subscript exp to be distinguished from input values, and the power spectrum 
$S(f)$ (used to estimate the Stokes drift according to (2.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.
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_{traj}$ corresponds to the number of spheres tracked and 
$\varLambda$ is the arrival rate of jumps in particle position. The mean Stokes drift 
$\langle u_{S} \rangle$ is calculated using (2.6), the mean Lagrangian velocity 
$\langle u_{L,exp}\rangle$ is obtained from tracer particle positions in experiments and the Eulerian-mean velocity 
$\langle u_{E,fit}\rangle$ is obtained so that experiments and model predictions for the mean agree (see § 3.1.3).
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 m
$^{-3}$, 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 m
$^{-3}$). 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 FOV. 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 with the velocity spectrum obtained from the measured surface elevation using linear wave theory in Appendix D. For frequencies up to a high-frequency cutoff 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. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021), 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 cutoff 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 with the Stokes drift in non-breaking waves due to mechanisms described in Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021).
For breaking waves, we cannot exclude the possibility of inertial (and therefore particle-property dependent) behaviour during the rapid accelerations (and decelerations) associated with breaking, although we have tried to minimize these effects by choosing particles that are as small and dense as possible. Using the definition of Stokes number in Santamaria et al. (Reference Santamaria, Boffetta, Martins Afonso, Mazzino, Onorato and Pugliese2013), 
$St =\omega \tau$, where 
$\tau =D^2/(12\beta \nu )$ and 
$\beta =3/(1+2\rho _w/\rho _p)$ with 
$\rho _w$ and 
$\rho _p$ the densities of water and the particle, respectively, we find a Stokes number that is a factor 3 larger than for the smaller (
$D=12$ mm vs 
$D=20$ mm) and denser (
$\rho _p=996$ kg m
$^{-3}$ vs 
$\rho _p=920$ kg m
$^{-3}$) used by Lenain et al. (Reference Lenain, Pizzo and Melville2019). If we use a Stokes number that does not take added mass into account, this factor becomes much larger. Lenain et al. (Reference Lenain, Pizzo and Melville2019) estimate there is a two-order-of-magnitude difference between the particle velocities during breaking and the particle rise velocity, justifying the assumption the particle behaves as a faithful tracer. The fact that we find similar ‘jump’ or ‘surf’ velocities as Lenain et al. (Reference Lenain, Pizzo and Melville2019) provides support that this assumption also applies here (see Appendix B).
3.1.3. Role of Eulerian-mean flow
While we do not take into account the Eulerian-mean flow 
$u_{E}$ in our model (2.2), this flow is present in experiments and therefore has to be accounted for to enable comparison (cf. (2.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 Reference van den Bremer and Breivik2017 and Monismith Reference Monismith2020).
The Eulerian flow was measured using three 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. (2.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.2), with an arbitrary (not measured) Eulerian-mean flow 
$\langle u_{E, fit}\rangle$, so that 
$b=\langle u_{S}\rangle +\langle u_{E, fit}\rangle$. The mean Stokes drift 
$\langle u_{S}\rangle$ is calculated using (2.6) from (a JONSWAP spectrum fitted to) the measured surface elevation spectrum. The arbitrary value of the Eulerian-mean flow 
$\langle u_{E, fit}\rangle$ is then chosen so that the mean Lagrangian drift predicted by the model 
$\langle u_{L,mod} \rangle ={\rm d}\langle X(t)\rangle /{\rm d}t$ (using 
$b=\langle u_{S}\rangle +\langle u_{E, fit}\rangle$) is equal to the mean Lagrangian drift in experiments 
$\langle u_{L,exp}\rangle$, that is 
$\langle u_{L,mod} \rangle =\langle u_{L,exp} \rangle$. Our interest in § 4, where we compare model predictions with experiments, is therefore in higher-order moments of particle position. The values 
$\langle u_{E, fit}\rangle$ thus obtained are within a reasonable margin of the values measured at depth 
$\langle u_{E, exp}\rangle$ (see Table 5 in Appendix C). We note for completeness that, for non-breaking waves, we expect that 
$\langle u_{L,mod} \rangle = \langle u_{S} \rangle + \langle u_{E,fit} \rangle$, while for breaking waves the jumps also make a contribution to the mean, 
$\langle u_{B}\rangle$. The next section will explain how this contribution is calculated, which must be done before 
$\langle u_{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 a discriminative correlation filter with channel and spatial reliability (Lukežič et al. Reference Lukežič, Vojíř, Čehovin Zajc, Matas and Kristan2018), 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 with 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_{traj}$ trajectories of equal length 
$\Delta t$ (for each significant wave height).
Figure 2. Example particle trajectories in irregular waves, showing experiments (light grey, blue) and Monte Carlo simulations of our model (orange): (a) 
$H_s = 0.05$ m, 
$\epsilon = 0.074$, showing almost no breaking events and trajectories evolving according to a Brownian motion; (b) 
$H_s = 0.17$ m, 
$\epsilon = 0.185$, 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 blue lines, where blue dashed lines show an oscillating trajectory, corresponding to the effect of waves, obtained from 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$).
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.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).
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 
$\varLambda$, 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 or 
$\epsilon = 0.074$, whilst breaking only rarely occurred for these experiments. The symbols 
$\blacktriangle$ and 
$\blacktriangledown$ respectively indicate the values of 
$\varLambda$ for 
$5\,\%$ lower and higher velocity thresholds to detect jumps, suggesting moderate sensitivity to the value of the threshold. 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 and for a comparison of the average jump velocity with the Stokes drift and with the results of Lenain et al. (Reference Lenain, Pizzo and Melville2019).
We use the jumps thus obtained to calibrate the jump process described by (2.9)–(2.10). Figure 3(a) shows the jump arrival rate 
$\varLambda$ as a function of steepness 
$\epsilon$ estimated from the experimental data for the four values of steepness. Also shown is the sigmoid function for 
$\varLambda (\epsilon )$ given by (2.9) with estimated values of the parameters in table 2. The variation in jump amplitudes 
$s$ estimated from the experimental data in figure 3(b) is captured well by the gamma distribution (2.10). Finally, we estimate the parameters of the gamma distribution (2.10) as linear function of steepness
\begin{equation} \alpha=a_\alpha+b_\alpha \epsilon, \quad \beta=a_\beta+b_\beta \epsilon,\end{equation}as shown in figure 3(c) with coefficients in table 2.
Figure 3. Calibration of the jump process for transport by breaking waves: (a) jump arrival rate 
$\varLambda$ as a function of steepness 
$\epsilon$. The symbols 
$\blacktriangle$ and 
$\blacktriangledown$ respectively indicate the values of 
$\varLambda$ for 
$5\,\%$ lower and higher velocity thresholds to detect jumps. (b) Probability density function of jump size 
$G(s) = \varGamma (s;\alpha,\beta )$ for 
$\epsilon = 0.122$ (red), 0.162 (green) and 0.185 m (yellow) (the case of 
$\epsilon =$ 0.074, 
$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 (3.2a,b).
4. Results
In this section, we will compare Monte Carlo simulations of our model (2.2) (using 
$10^5$ trajectories) with experiments. The Monte Carlo simulations agree perfectly with the exact solutions for the first three moments (2.17)–(2.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(a) 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.2) (orange) for the smallest steepness waves (
$\epsilon =0.074$, 
$H_s=0.05$ m). A negligible number of waves break for this case such that 
$\varLambda \approx 0$, and the wave breaking term does not contribute in (2.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.
Figure 4. Comparison between model (orange) and experimental data (blue) for non-breaking irregular waves (
$\epsilon = 0.074$, 
$H_s = 0.05$ m): (a) normalized mean position 
$\langle X(t) \rangle /\lambda _p$ (thick line) and 
$\pm 1$ standard deviation (dashed) for experiments (blue) and model simulations (orange), (b) normalized mean position (thick line) for experiments (solid blue), model simulations (dashed-dotted 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$ (medium-thickness line) for experiments (solid blue) and model simulations (dashed-dotted orange) are approximately proportional to time 
$t$. The normalized skewness 
$\langle |\tilde {X}(t)|^3\rangle /\lambda _p^3$ (thin line) is negligible in the experiments (solid blue).
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 4(b). Indeed, extracting the power-law behaviour in figure 4(b), the experimental particle position variance exhibits a linear 
$t$ dependence, in accordance with the single-particle Taylor diffusion prediction by Herterich & Hasselmann (Reference Herterich and Hasselmann1982), which in turn agrees with our theoretical prediction (2.18) when 
$\varLambda =0$ (no breaking). It is interesting to note that a Wiener or normal diffusion process results in a particle position variance that 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 predicted by Herterich & Hasselmann (Reference Herterich and Hasselmann1982) and observed in our experiments) is in disagreement with Farazmand & Sapsis (Reference Farazmand and Sapsis2019), 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 (2.19).
4.2. Breaking waves: comparison between experiments and model predictions
Figure 5(a) shows the mean particle position (dashed-dotted) 
$\pm 1$ standard deviation (dashed lines) in the experiments (red) and the Monte Carlo simulations of (2.2) (orange) for the largest steepness waves (
$\epsilon =0.185$, 
$H_s=0.17$ m). For this significant wave height, many jumps in position due to breaking occur, as illustrated earlier in figure 2(b). 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 (2.17) with 
$\varLambda \neq 0$. The power-law behaviour in figure 5(b) indicates that, in agreement with (2.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 (2.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 (2.19), thin orange dashed-dotted line in figure 5(b), is also observed in experiment (thin blue solid line).
Figure 5. Comparison between model (orange) and experiment data (blue) for breaking irregular waves (
$\epsilon = 0.074$, 
$H_s = 0.17$ m): (a) normalized mean position 
$\langle X(t) \rangle /\lambda _p$ (thick line) and 
$\pm 1$ standard deviation (dashed lines) for experiments (blue) and model simulations (orange), (b) normalized mean position (thick line) for experiment (solid, blue), model simulations (dashed-dotted, 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$ (medium thickness) for experiments (solid blue) and model simulations (dashed-dotted orange) are approximately proportional to time 
$t$. The normalized skewness 
$\langle |\tilde {X}(t)|^3\rangle /\lambda _p^3$ (thin line) is finite for this wave steepness with order-of-magnitude agreement between experiments (solid blue) and model simulations (dashed-dotted, orange).
Pizzo (Reference Pizzo2017) and Deike et al. (Reference Deike, Pizzo and Melville2017) have shown earlier that Lagrangian tracers on a sweet spot just below the crest ‘surf’ forward, having a forward velocity an order of magnitude higher than the Stokes drift, as confirmed by experiments in Lenain et al. (Reference Lenain, Pizzo and Melville2019) and Sinnis et al. (Reference Sinnis, Grare, Lenain and Pizzo2021). As Lenain et al. (Reference Lenain, Pizzo and Melville2019) and Sinnis et al. (Reference Sinnis, Grare, Lenain and Pizzo2021) used wave groups instead of random waves, a direct comparison cannot be made. Nevertheless, we observe comparable mean ‘jump’ or ‘surf’ velocities to those observed by Lenain et al. (Reference Lenain, Pizzo and Melville2019) (see figures B3–B4 and table 4 in Appendix B).
4.3. Model predictions as a function of steepness
As we describe particle transport as a stochastic process, we will examine how the first three moments, that is the mean velocity and variance and skewness of the particle position distribution, 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 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.
Figure 6(a) shows the mean drift 
$\langle u_{L,mod} \rangle = {\rm d}\langle X(t)\rangle /{\rm d}t$ based on (2.2). If the breaking term is not taken into account in (2.2) (i.e. 
$\varLambda =0$), this drift would be equal to the theoretical Stokes drift value 
$\langle u_{S}\rangle \propto \epsilon ^2$ (dashed line). Taking into account the jump-diffusion process of (2.2), 
$\langle u_{L,mod} \rangle$ starts deviating from 
$\langle u_{S} \rangle$ for steep waves in which many breaking events occur. The nature of this deviation is given in (2.17), where 
$\langle u_{L,mod} \rangle$ is increased by a factor 
$(1 + \alpha \varLambda /(\langle u_{S} \rangle \beta ))$. The second moment or variance deviates in a similar fashion from the non-breaking case in figure 6(b), deviating by a factor 
$(1 + \alpha (\alpha +1)\varLambda /(\sigma ^2 \beta ) )$. The jump process can therefore induce a significantly increased spread in particle position. 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. The finite skewness marks a deviation from the normal distribution and breaks the symmetry. Note that, physically, a Gaussian distribution can only evolve to a non-Gaussian distribution through a nonlinear process. A practical implication is that there will be more extremes as the tails of the distribution are larger relatively.
Note that when studying drift induced by individual focused wave groups, Pizzo et al. (Reference Pizzo, Melville and Deike2019) 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 3(a).
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. (Reference Deike, Pizzo and Melville2017), Pizzo et al. (Reference Pizzo, Melville and Deike2019), Lenain et al. (Reference Lenain, Pizzo and Melville2019) and Sinnis et al. (Reference Sinnis, Grare, Lenain and Pizzo2021). The framework consisting of a SDE 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. Furthermore, we have described the evolution and quantified the uncertainty of particle transport under the influence of wave breaking by modelling 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. We show that a jump-diffusion process is a feasible option for describing the experimentally obtained drift and diffusion of particle position by means of a parameterization that is a continuous function of mean wave steepness. Our results for the enhanced mean drift due to breaking waves are in approximate quantitative agreement with the experiments of Lenain et al. (Reference Lenain, Pizzo and Melville2019), who considered focused wave groups instead of irregular waves.
Looking forward, it will be desirable to develop a predictive rather than just a descriptive model for realistic ocean applications. That is, for given wave properties such as steepness (or even wind speed), one would like to be able to predict particle drift and diffusion. For our model to make such a prediction requires knowing the jump frequency (number of encounters with a breaking crest per unit time) and the jump amplitude distribution (modelled as a gamma distribution). Currently, we have fitted these parameters as a function of steepness based on four experiments. For these parametric relations with steepness to apply in more general settings, more experimental data need to be gathered. In this paper, we have considered long-crested (or unidirectional) waves. In the ocean, sea states are almost always directionally spread (or short crested). We envisage our model can be readily extended to directional seas using the results of Kenyon (Reference Kenyon1969) for non-breaking waves. Calibrating the jump-diffusion process for directionally spread breaking waves will be one example that requires new laboratory experiments. We also recommend future experiments consider the role of particle size and density to ensure particles act as faithful (i.e. Lagrangian) tracers during breaking and results are particle-property independent.
As an alternative to the Poisson process and the relationships with steepness in our model, the jump frequency and the jump amplitude could be related to the breaking statistic 
$\varLambda (c)$ (Phillips Reference Phillips1985). This would also allow existing parameterizations for 
$\varLambda (c)$ as a function of wind speed (Romero & Melville Reference Romero and Melville2010; Sutherland & Melville Reference Sutherland and Melville2013; Lenain & Melville Reference Lenain and Melville2017; Deike & Melville Reference Deike and Melville2018) to be used for particle transport predictions, such that this can be implemented based on spectral wave-model predictions (e.g. WaveWatchIII) of the sea state (Romero Reference Romero2019). New laboratory experiments would have to simultaneously measure 
$\varLambda (c)$, for example using infra-red cameras (see Deike & Melville Reference Deike and Melville2018).
Finally, we emphasize that any parameterization from observed quantities will come with an uncertainty. To be able to give a truly stochastic prediction, this uncertainty should be included by making the parameters in our model stochastic quantities too. 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 A. de Fockert and W. Bakker at Deltares and P. van der Gaag, A. van der Vlies and A. Doorn at Delft University of Technology for their help setting up and conducting the experiment. We thank L. Lenain and N. Pizzo of the Scripps Institution of Oceanography for fruitful discussions and for sharing their data to test our ideas.
Funding
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.
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 intrinsics were found by detecting multiple images of a checkerboard in different orientations across the entire FOV. Figure 7(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 7(b). A correlation filter with channel and spatial reliability 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 FOV.
Figure 7. (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.
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.
Appendix B. Jump detection
Figure 8 displays the steps of the jump detection and jump amplitude (distance travelled in a jump) estimation process for two example trajectories. For 
$\epsilon =0.074$, panel (a) shows the derivative of particle position over a time step 
$\delta t$ determined by the camera frame rate: 
$u_{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.3 c$, 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_{I} (t)-u_{th})$ mark the time segments where the particle is classified as ‘jumping’. In panels (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 8. Jump detection and amplitude estimation process for (a,c,e) the least steep waves with 
$\epsilon = 0.074$ and (b,d,f) the steepest waves with 
$\epsilon = 0.185$. (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.3 c$, with 
$c$ the phase velocity. (c,d) The Heaviside function 
$\mathcal {H}(u_{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’.
Figure 9 displays a histogram of 
$u_{I} = \delta X / \delta t$ at each time step 
$\delta t$. For 
$\epsilon = 0.074$ m, the distribution is symmetric, whereas for higher steepnesses, a heavier tail is observed on the right. Note that as the histogram shows (discretized) instantaneous velocities, this does not imply the mean velocity during a jump attains these higher values. In addition, we note that when the velocity exceeds the threshold for only one time step this is not considered to be a jump. Figure 10 then shows the mean velocity during each jump, calculated as 
$u_j=s/t_s$ with 
$s$ the distance travelled during the jump and 
$t_s$ the duration of the jump, estimated as illustrated in figure 8.
Figure 9. Histogram of the (discretized) instantaneous velocity 
$\delta X / \delta t$, normalized by 
$c$, for the four values of steepness considered.
Figure 10. Histogram of the mean jump velocity 
$u_j=s/t_s$, normalized by 
$c$, for the four values of steepness considered.
Table 3 shows a comparison between the normalized mean Stokes drift calculated from the measured wave spectrum (i.e. using (2.6)) and the normalized mean jump velocity
\begin{equation} \langle u_j\rangle / c = \frac{1}{N_j}\sum_{i=1}^{N_j} u_{j,i} /c, \end{equation}
where 
$N_j$ is the number of jumps. The normalized mean jump velocity is an order of magnitude higher.
Table 3. Comparison between the normalized average Stokes drift and the normalized average jump velocity.
Figure 11 shows that, as in Deike et al. (Reference Deike, Pizzo and Melville2017) and Lenain et al. (Reference Lenain, Pizzo and Melville2019), the scaling with steepness is also different for the two velocities. Performing a linear fit for the average jump velocity, 
$\langle u_{j}\rangle /c = a_2 \epsilon$ (black dashed line), we find a slope 
$a_2 = 2.02$. Performing a quadratic fit for the average Stokes drift 
$\langle u_{S} \rangle /c = a_1 \epsilon ^2$ gives 
$a_1 = 0.75$.
Figure 11. Normalized average jump velocity 
$\langle u_{j}\rangle /c$ and normalized average Stokes drift 
$\langle u_s\rangle /c$ as a function of steepness for the four steepness conditions. Dashed black line: linear fit 
$\langle u_{j}\rangle /c= 2.02\epsilon$, red dashed line: quadratic fit 
$\langle u_s \rangle /c=0.75\epsilon ^2$.
Deike et al. (Reference Deike, Pizzo and Melville2017) and Lenain et al. (Reference Lenain, Pizzo and Melville2019) calculate a similar metric, based on the displacement caused by a single focused wave group. In these papers, the surface drift velocity is not directly measured but calculated as 
$\langle \Delta x\rangle /T$, where 
$T = 1/f_c$, 
$f_c$ the central frequency of the packet and 
$\Delta x$ the total displacement. In addition, these papers calculate the velocities as a function of the maximum linear slope 
$S$ instead of the characteristic steepness 
$\epsilon$. That is, 
$\langle u_{j}\rangle /c = \chi _2 S$ and 
$\langle u_{S}\rangle /c = \chi _1 S^2$, and Deike et al. (Reference Deike, Pizzo and Melville2017) and Lenain et al. (Reference Lenain, Pizzo and Melville2019) estimate 
$\chi _2$ in the range 
$6.1$–
$9$ and 
$\chi _1$ in the range 
$0.82$–
$1.28$.
To make an approximate comparison with these values, we assume 
$S=a k_p$ and 
$H_s = 2 \sqrt {2}a$, yielding 
$\epsilon = \sqrt {2}S$ and, therefore, 
$a_2 = \chi _2/\sqrt {2}$ and 
$a_1= \chi _1/2$. Using this conversion, Table 4 shows the approximate values for 
$a_1$ and 
$a_2$ for Deike et al. (Reference Deike, Pizzo and Melville2017) and Lenain et al. (Reference Lenain, Pizzo and Melville2019). Our study gives comparable values. We ascribe the lower value of the average drift velocity to the fact that 
$\epsilon$ is an average steepness, and consequently there will be a higher number of lower-sloped breaking waves contributing to the average drift in B1.
Table 4. Fitted coefficients for 
$\langle u_{S}\rangle /c = a_1 \epsilon ^2$ and 
$\langle u_{j}\rangle /c = a_2 \epsilon$ for our study and Deike et al. (Reference Deike, Pizzo and Melville2017) and Lenain et al. (Reference Lenain, Pizzo and Melville2019).
Table 5. Difference between the experimentally measured Eulerian-mean flow and its value used in the model.
Appendix C. Drift velocities
Figure 12 shows various drift velocity estimations. For our experiments, the Stokes drift is based on a JONSWAP spectrum fitted to the measured spectrum, namely 
$\langle u_{{S}, f} \rangle$, indicated by the blue triangles, with the dashed blue line its quadratic fit.
Figure 12. Drift velocity estimations and measurements as a function of steepness. Stokes drift based on a JONSWAP fit on the spectrum truncated at 
$\langle u_{{S}, f} \rangle$ (
$\vartriangle$) and its quadratic fit (dashed line), measured particle drift 
$\langle u_{L,exp} \rangle$ (thick 
$\times$), measured particle drift corrected by model mean flow 
$\langle u_{L,exp} \rangle -\langle u_{E,fit} \rangle$ (thin 
$\times$).
We correct the experimentally measured mean drift velocity 
$\langle u_{L,exp} \rangle$ by the mean flow in the model 
$\langle u_{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.
Appendix D. Velocity spectrum particles
Figures 13–16 compare the velocity spectrum of the particle trajectories (dark red) with 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 13. Case 
$H_s = 0.05$ m, 
$\epsilon = 0.074$. Spectrum of the velocity of the particles (red), dashed lines indicate 
$\pm 1$ standard deviation. The linear velocity calculated from the spectrum of the surface elevation 
$\eta$ is shown in pink. The grey dashed line indicates 
$f=3f_0$. (a) Linear scale. (b) Log scale.
Figure 14. Case 
$H_s = 0.09$ m, 
$\epsilon = 0.122$. Lines same as in figure 13.
Figure 15. Case 
$H_s = 0.12$ m, 
$\epsilon = 0.162$. Lines same as in figure 13.
Figure 16. Case 
$H_s = 0.17$ m, 
$\epsilon = 0.185$. Lines same as in figure 13.
Open access
