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 )$

(3.1)\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.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

(3.2a,b)\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).

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 $\varLambda (\epsilon )$ according to (2.9) and the jump size probability density function $G(s)$ according to (2.10).