2.1. Surface gravity waves

We consider a train of progressive, two-dimensional surface gravity waves travelling in the $x^{\prime }$-direction in water of depth $h$, as described by linear wave theory. In the $x'-z'$ plane, the fluid velocities are

(2.1a)$$\begin{gather} u_x^{\prime} = \frac{gka}{\omega} \frac{\cosh (k(z^{\prime}+h))}{\cosh kh} \cos (kx^{\prime} - \omega t^{\prime}), \end{gather}$$

(2.1b)$$\begin{gather}u_z^{\prime} = \frac{gka}{\omega} \frac{\sinh (k(z^{\prime}+h))}{\cosh kh} \sin (kx^{\prime} - \omega t^{\prime}). \end{gather}$$

Here, $t'$ is time, $a$ is the wave amplitude, $k$ is the wavenumber, $\omega$ is the angular frequency and $\boldsymbol {u}^{\prime } = [u_x^{\prime }, u_z^{\prime }]$ gives the fluid velocity field in $-h \leqslant z^{\prime } \leqslant 0$. The free-surface displacement is $\eta ^{\prime } =a \cos (kx^{\prime } - \omega t^{\prime })$, and the dispersion relationship is ${\omega }^{2} = gk \tanh kh$ where $g$ is gravitational acceleration. The phase velocity of the wave train is $c=\omega /k$. We make this system of equations dimensionless by defining the following variables: $t = \omega t^{\prime }, \boldsymbol {x} = k \boldsymbol {x}^{\prime }, \boldsymbol {u} = \boldsymbol {u}^{\prime } /c, \eta =k\eta '$. In dimensionless form, the free surface is defined as $\eta =ka\cos (x - t)$, and the velocity field is given by

(2.2a)$$\begin{gather} u_x = \varepsilon \frac{\cosh (z+kh)}{\cosh kh} \cos (x - t), \end{gather}$$

(2.2b)$$\begin{gather}u_z = \varepsilon \frac{\sinh (z+kh)}{\cosh kh} \sin (x - t), \end{gather}$$

where

(2.3)\begin{equation} \varepsilon = \frac{k a}{\tanh kh} \end{equation}

is a small parameter that we refer to as the wave nonlinearity parameter. In the deep-water limit ($kh \rightarrow \infty$), $\varepsilon \rightarrow ka$ and in the shallow-water limit ($kh \rightarrow 0$), $\varepsilon \rightarrow a/h$. The wave-induced flow field is completely characterised by the dimensionless variables $ka$ and $kh$.

2.2. Equations of particle motion

The equations of motion for a spherical particle of diameter $d_p$ moving in an unsteady and non-uniform flow are given by

(2.4a)$$\begin{gather} \frac{\mathrm{d}\kern0.06em \boldsymbol {x_p}^{\prime}}{\mathrm{d}t^{\prime}} =\boldsymbol{v}^{\prime} \end{gather}$$

(2.4b)$$\begin{gather}m_p \frac{\mathrm{d} \boldsymbol{v}^{\prime}}{\mathrm{d} t^{\prime}} = (m_p - m_f) \boldsymbol{g} + m_f \frac{\mathrm{D} \boldsymbol{u}^{\prime}}{\mathrm{D} t^{\prime}} - C_m m_f \left[\frac{\mathrm{d} \boldsymbol{v}^{\prime}}{\mathrm{d} t^{\prime}} - \frac{\mathrm{D} \boldsymbol{u}^{\prime}}{\mathrm{D} t^{\prime}} \right] + \boldsymbol{F}_{{drag}} , \end{gather}$$

where $\boldsymbol {x_p}^{\prime }$ is the particle position, $\boldsymbol {v}^{\prime }$ is the particle velocity, $\boldsymbol {u}^{\prime }$ is the fluid velocity at the particle's location and $\boldsymbol {g}$ is the gravitational acceleration vector. The particle mass is $m_p$, the mass of fluid displaced by the particle is $m_f$ and $C_m$ is the added mass coefficient for the particle. The added mass coefficient is defined as the added inertia associated with the fluid that has to accelerate around the particle normalised by the mass of fluid displaced by the particle. For spherical particles, $C_m = 0.5$ (Newman Reference Newman1977). The rate of change following the particle is $\mathrm {d}/\mathrm {d}t (= \partial /\partial t + \boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {\nabla })$ and the rate of change following a fluid parcel is $\mathrm {D}/\mathrm {D}t (= \partial /\partial t + \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla })$. The terms on the right-hand side of (2.4b) are the forces on the particle: buoyancy, fluid forcing, added mass and drag, respectively.

We make the equations of particle motion dimensionless in a similar manner to the fluid velocities but with the additional dimensionless variables $\boldsymbol {x_p} = k \boldsymbol {x_p}^{\prime }$, $\boldsymbol {v} = \boldsymbol {v}^{\prime } /c$, $\gamma = \rho _p/\rho$, and $C_D = F_{{drag}}/(\tfrac {1}{8}\rho (\boldsymbol {v}^{\prime } - \boldsymbol {u}^{\prime })^{2} {\rm \pi}d_p^{2})$, where $\rho _p$ is the particle density, $d_p$ is the particle diameter and $\rho$ is the fluid density. The dimensionless equations of motion are

(2.5a)$$\begin{gather} \frac{\mathrm{d}\kern0.06em\boldsymbol{x_p}}{\mathrm{d}t} =\boldsymbol{v} , \end{gather}$$

(2.5b)$$\begin{gather}\frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}t} ={-} \frac{\gamma - 1}{\gamma + C_m} \frac{1}{\tanh kh} \boldsymbol{e_z} + \frac{1 + C_m}{\gamma + C_m} \frac{\mathrm{D} \boldsymbol{u}}{\mathrm{D}t} - \frac{1}{\gamma + C_m} \frac{1}{kd_p} \frac{3}{4} C_D \lvert \boldsymbol{v} - \boldsymbol{u} \rvert (\boldsymbol{v} - \boldsymbol{u}), \end{gather}$$

where $\boldsymbol {e_z}$ is the unit vector antiparallel to gravity and where we have used the dispersion relation $\omega ^{2} = gk \tanh kh$.

From (2.5), we can see that four dimensionless parameters are required to compute particle motion: the density ratio $\gamma$, the dimensionless diameter $kd_p$, the added mass coefficient $C_m$ and the drag coefficient $C_D$. While $\gamma$ and $kd_p$ are fixed for a given particle in a given wave field, $C_D$ varies dynamically as a function of the particle Reynolds number based on the slip velocity

(2.6)\begin{equation} {{\textit{Re}}}_{p} = \frac{\lvert \boldsymbol{v}^{\prime} - \boldsymbol{u}^{\prime} \rvert d_p}{\nu}, \end{equation}

where $\nu$ is the kinematic viscosity of the fluid.

For ${{\textit {Re}}}_p \ll 1$, the drag coefficient is given by the Stokes drag model, $C_D = 24/{{\textit {Re}}}_{p}$, which results in a drag force that is linear with respect to the slip velocity. For finite values of ${{\textit {Re}}}_p$, empirical relationships for the drag coefficient show that the drag force is nonlinear with respect to the slip velocity. Here, we use the Schiller–Naumann drag model, $C_D = (24/{{\textit {Re}}}_p)(1 + 0.15 {{\textit {Re}}}_p^{0.687})$, which is $\pm 5\,\%$ accurate for ${{\textit {Re}}}_{p} \lesssim 10^{3}$ (Clift, Grace & Weber Reference Clift, Grace and Weber2005).

Linear drag: using the Stokes drag model, (2.5b) reduces to

(2.7)\begin{equation} \frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}t} = \frac{(\boldsymbol{u} - \boldsymbol{v})}{{St_{lin}}} + \beta \frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D} t} - \frac{v_{s,{lin}}}{{St_{lin}}} {\boldsymbol{e_{z}}} , \end{equation}

where ${St_{lin}} = \omega \tau _{p}$ is the particle Stokes number with $\tau _p = d_p^{2} (C_m + \gamma )/18\nu$ being the time scale over which the particle exponentially relaxes to match the flow velocity and $\beta = ({1 + C_m})/({\gamma + C_m})$. For spherical particles, the parameter $\beta \in [0, 3]$ with $\beta <1$ for negatively buoyant particles and $\beta >1$ for positively buoyant particles. The dimensionless particle terminal settling velocity is $v_{s,{lin}} = ({{St_{lin}} (1 - \beta )})/{\tanh kh}$, where the subscript ‘lin’ indicates the use of the linear drag model. Equation (2.7) is similar to Reference Maxey and RileyMaxey & Riley's (Reference Maxey and Riley1983) formulation (their (43)); the differences are that we have neglected the Basset history term and we have used $\mathrm {D}\boldsymbol {u}/\mathrm {D}t$ instead of $\mathrm {d}\boldsymbol {u}/\mathrm {d}t$ in the added mass term (Auton, Hunt & Prud'Homme Reference Auton, Hunt and Prud'Homme1988). We have also neglected the Faxén corrections, but they are irrelevant here since ${\nabla }^{2}\boldsymbol {u} \equiv 0$ for the wave-induced flow field.

Nonlinear drag: using the Schiller–Naumann drag model, (2.5b) reduces to

(2.8)\begin{equation} \frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}t} = \frac{(\boldsymbol{u} - \boldsymbol{v})}{{St_{SN}}} + \beta \frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D} t} - \frac{v_{s,{SN}}}{{St_{SN}}_{\!,t}} {\boldsymbol{e_{z}}} , \end{equation}

where ${St_{SN}} = \omega \tau _{p,SN}$ is the particle Stokes number for the Schiller–Naumann drag model with $\tau _{p,{SN}} = d_p^{2} (C_m + \gamma )/[18\nu (1 + 0.15 {{\textit {Re}}}_p^{0.687})]$ being the corresponding time scale. The subscript ‘SN’ indicates the use of the Schiller–Naumann drag model. In this model, the particle's terminal settling velocity must be found iteratively using the empirical drag relationship, which is then used to calculate the particle Stokes number at terminal velocity, ${St_{SN}}_{\!,t}$, using the particle Reynolds number at terminal velocity, ${{\textit {Re}}}_{p,t}=v_{s,{SN}}d_p/\nu$. While ${St_{SN}}_{\!,t}$ is constant for a given particle in a given wave field, ${St_{SN}}$ is a dynamic quantity that evolves as the particle travels through the flow. Under the simplification ${St_{SN}} \approx {St_{SN}}_{\!,t}$, which is equivalent to simplifying the drag coefficient as $C_D = (24/{{\textit {Re}}}_p)(1 + 0.15 {{\textit {Re}}}_{p,t}^{0.687})$, the equations for linear drag and nonlinear drag become identical except for the definition of the Stokes number and the terminal settling velocity. This simplification can be interpreted as assuming that the particle time scale is constant and independent of the flow, or equivalently, a ‘linearisation’ of the drag model such that drag force is linear with respect to the dynamic slip velocity, but contains the nonlinear aspects with respect to the particle's terminal settling velocity.

To show the applicability of linear and nonlinear drag models, the particle time scale and the particle Reynolds number based on the terminal velocity are plotted as functions of the particle density and particle diameter in figure 1. These quantities have been calculated using the Schiller–Naumann model. The range of particle densities and sizes cover small particles close to neutral buoyancy (e.g. diatoms, $\gamma \approx 1.01; d_p \approx 100\ \mathrm {\mu }{\rm m}$) to large dense particles (e.g. sediment pebbles, $\gamma \approx 2.65, d_p \approx 1\ {\rm cm}$). Over most of this parameter space, we see that $1 < {{\textit {Re}}}_p < 10^{3}$ (figure 1a), suggesting that the Schiller–Naumann nonlinear drag model is relevant for many particles of environmental relevance. We also see that the particle time scales are small compared with typical wave periods in the laboratory ($O(1$)s) and in the field ($O(10$)s) (figure 1b), suggesting that a small Stokes number approximation is appropriate.

Figure 1. Particle parameters calculated with the Schiller–Naumann drag model for $\rho = 10^{3}\ {\rm kg}\ {\rm m}^{-3}$ and $\nu = 10^{-6}\ {\rm m}^{2}\ {\rm s}^{-1}$: (a) particle Reynolds number at terminal velocity; (b) particle time scale in seconds.

2.3. Particle motion and drift

We further analyse the particle equation of motion under the approximation of a small Stokes number. We write the particle equation as

(2.9)\begin{equation} \frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}t} = \frac{(\boldsymbol{u} - \boldsymbol{v})}{{St}} + \beta \frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D} t} - \frac{v_{s}}{{St}} {\boldsymbol{e_{z}}} , \end{equation}

with the choice of drag model being implicit and affecting the values of $v_s$ and ${St}$. In the case of the Schiller–Naumann drag model, this also assumes that the Stokes number is constant. From an expansion of the particle velocity in ${St}$, we find kinematic solutions to the particle motion to be (as derived in Appendix A)

(2.10a)$$\begin{gather} {v_x} = \varepsilon A \frac{\cosh (z+kh)}{\cosh kh} \cos\left(x - t + \phi \right) + \varepsilon^{2} \frac{{St}(1 - \beta)}{2 \cosh^{2} kh} \sin (2(x - t)) , \end{gather}$$

(2.10b)$$\begin{gather}{v_z} ={-} v_s + \varepsilon A \frac{\sinh (z+kh)}{\cosh kh} \sin\left(x - t + \phi \right) - \varepsilon^{2} \frac{{St}(1 - \beta)}{2 \cosh^{2} kh} \sinh (2(z + kh)), \end{gather}$$

where $A$ and $\phi$ give the relative magnitude and phase shift for the leading-order particle orbital motion

(2.11a)$$\begin{gather} A = \sqrt{\left[1 - {St}^{2} \left( 1 - \beta \right) \right]^{2} + \left[{St}\left( 1 - \beta \right)\right]^{2}}, \end{gather}$$

(2.11b)$$\begin{gather}\phi = \arctan \left[ \frac{{St} \left( 1 - \beta \right)}{1 - {St}^{2} \left( 1 - \beta \right)} \right]. \end{gather}$$

In (2.10), the particle velocity includes terms of $O(\varepsilon )$ and $O(\varepsilon ^{2})$ as well as the constant terminal settling velocity $v_s$. The $O(\varepsilon )$ terms describe how inertia modifies particle motion relative to wave orbital motion experienced by fluid tracers. The tracer particle limit is described by $A=1$ and $\phi =0$, and deviations from this limit grow with increasing ${St}$ and increasing deviation from neutral buoyancy (where neutrally buoyant particles have $\gamma =1$ and $\beta =1$). The relative amplitude $A$ and phase shift $\phi$ are plotted for negatively buoyant particles as functions of $St$ and $\gamma$ in figure 2. For settling particles ($\gamma >1$), $A<1$ and $\phi >0$, showing that they trace out smaller wave orbitals than fluid parcels and lag the fluid velocity in time. These terms are equally valid for positively buoyant particles ($\gamma <1$); they show the opposite behaviour, tracing out larger wave orbitals than fluid parcels and leading the fluid velocity. The $O(\varepsilon ^{2})$ terms arise from the coupling between particle inertia and fluid inertia (specifically, the advective part of the fluid acceleration; (A6)). In the horizontal direction, this coupling produces a small oscillation in the particle motion at twice the wave frequency. In the vertical direction, it produces a depth-dependent enhancement of the particle vertical motion. For settling particles ($\gamma >1$), the settling velocity is enhanced, whereas for positively buoyant particles ($\gamma <1$), the rise velocity is enhanced.

Figure 2. Particle motion parameters as functions of Stokes number ${St}$ for different density ratios $\gamma$ calculated using (2.11). (a) Amplitude $A$ and (b) phase shift $\phi$ of particle velocity relative to tracer particles. Lines correspond to a different density ratios, $\gamma = 1, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7, 3$.

From the kinematic solution to particle motion (2.10), we can next derive the wave-averaged particle horizontal and vertical drift velocities using a two-time expansion in the wave nonlinearity parameter $\varepsilon$ (as shown in Appendix B). The resulting particle drift velocities can be written in terms of the classical Stokes drift velocity of tracer particles

(2.12a)$$\begin{gather} v_{x\text{-}drift} =\frac{A^{2}}{1+v_s^{2}} u_{{SD}} , \end{gather}$$

(2.12b)$$\begin{gather}v_{z\text{-}drift} ={-}v_s \left[1 + \frac{A^{2}}{1+v_s^{2}}u_{{SD}} + \frac{1}{2}\tanh kh\frac{\mathrm{d}u_{{SD}}}{\mathrm{d}z_{p0}}\right], \end{gather}$$

where $u_{{SD}}$ is the Stokes drift velocity, which is a function of the wave-averaged vertical position $z_{p0}$ (Stokes Reference Stokes1847; van den Bremer & Breivik Reference van den Bremer and Breivik2017)

(2.13)\begin{equation} u_{{SD}}(z_{p0}) = \varepsilon^{2} \frac{\cosh(2(z_{p0}+kh))}{2\cosh^{2} kh}. \end{equation}

Equations (2.12) show that the horizontal drift is reduced relative to the Stokes drift of tracers and that the settling velocity is enhanced relative to the particle's terminal settling velocity. (Although not the focus of this study, this result also applies to positively buoyant particles where the rise velocity is also enhanced.) In the tracer limit of ${St} \rightarrow 0$ and $\gamma \rightarrow 1$, we recover the classical Stokes drift velocity. In the deep-water limit of $kh \rightarrow \infty$, we recover results identical to Santamaria et al. (Reference Santamaria, Boffetta, Afonso, Mazzino, Onorato and Pugliese2013) up to $O({St})$, with slight differences occurring at higher orders due to differences in the dynamical equation (cf. (2.9) with their (4)).

The drift modifications are functions of wave and particle parameters. Figure 3 shows that the horizontal drift velocity of settling inertial particles is reduced relative to the Stokes drift of fluid tracers, with this effect becoming stronger in shallow-water waves (smaller $kh$), for denser particles (larger $\gamma$), and for particles with more inertia (larger ${St}$). This reduction in horizontal drift can be understood via the kinematic solution to particle motion (2.10): since the orbital motion amplitude of inertial particles is reduced relative to fluid tracers, the net displacement of these particles over a wave period is smaller than it is for fluid tracers.

Figure 3. Horizontal drift for settling particles relative to tracer particles as functions of Stokes number ${St}$ for different density ratios $\gamma$ at different relative depths $kh$, calculated using (2.12a). Lines correspond to different density ratios, $\gamma = 1, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7, 3$; (a) $kh=0.3$, (b) $kh=1.0$, (c) $kh=3.0$.

Similar to the horizontal drift, the vertical drift enhancement can be understood via the kinematic solution to particle motion (2.10). One part of the settling enhancement comes from settling particles sampling a stronger wave-induced flow during the downward part of their orbital motion than the subsequent upward part of their orbital motion since the particles settle lower into the water column during each wave cycle. Another part of the settling enhancement comes from the coupling between particle inertia and fluid inertia. In figure 4 we plot the settling enhancement across a range of particle and wave parameters. By varying $kh$ while keeping $\varepsilon$ constant, we find that the relative enhancement of settling velocity is larger in deeper water (larger $kh$), for less dense particles (smaller $\gamma$), and for particles with less inertia (smaller ${St}$). If instead $kh$ is varied while keeping $ka$ constant, the relative enhancement of settling velocity is larger in shallower water (smaller $kh$). Therefore, the vertical drift enhancement is a complex function of both $ka$ and $kh$. The settling enhancement also depends on the vertical position of the particle, with the effect always being strongest near the surface and decaying with increasing depth.

Figure 4. Vertical drift for settling particles relative to their terminal settling velocity as functions of Stokes number ${St}$ for different density ratios $\gamma$ at different relative depths $kh$, calculated using (2.12b) with $\varepsilon = 0.33$ and $z_{p0} = 0$. Lines correspond to different density ratios, $\gamma = 1.2, 1.5, 1.8, 2.1, 2.4, 2.7, 3$; note $\gamma = 1$ is not shown since particles do not settle in this case; (a) $kh=0.3$, (b) $kh=1.0$, (c) $kh=3.0$.

These results show that the modifications of horizontal and vertical drifts have different physical origins. The horizontal drift has the same structure and origin as classical Stokes drift: it arises from the correlation between the particle orbital motion and the Eulerian particle velocity field, and it is therefore a kinematic effect arising from the $O(\varepsilon )$ term in (2.10a). The reduction in the horizontal drift relative to classical Stokes drift is a dynamical effect because particles undergo smaller orbital motions because of inertial effects. The vertical drift has two components: a kinematic effect that is a vertical analogue to Stokes drift for settling particles and a dynamic effect from the coupling of particle inertia with the fluid inertia. The kinematic effect is the term proportional to $u_{{SD}}$ in (2.12b) and comes from the $O(\varepsilon )$ term in (2.10b), whereas the dynamic effect is the term proportional to $\mathrm {d}u_{{SD}}/\mathrm {d}z$ and comes from the $O(\varepsilon ^{2})$ term in (2.10b).

2.4. Particle dispersion

To quantify particle dispersion in the wave-averaged sense, we can compute the divergence of the particle drift velocity (2.12). We find that it is non-zero, where $\boldsymbol {\nabla }_{\boldsymbol {x}_{p0}} \boldsymbol {\cdot } \boldsymbol {v}_{{drift}} = {\partial v_{x\text {-}drift}} /{\partial x_{p0}} + {\partial v_{z\text {-}drift}} /{\partial z_{p0}}$, and

(2.14a)$$\begin{gather} \frac{\partial v_{x\text{-}drift}} {\partial x_{p0}} = 0 , \end{gather}$$

(2.14b)$$\begin{gather}\frac{\partial v_{z\text{-}drift}} {\partial z_{p0}} ={-} v_s \left[\frac{A^{2}}{1+v_s^{2}} \frac{\mathrm{d}u_{{SD}}}{\mathrm{d}z_{p0}} + \tanh kh \, u_{{SD}}\right]. \end{gather}$$

This shows that the wave-averaged inertial particle field has non-zero vertical divergence, but zero horizontal divergence. Negatively buoyant particles will converge in the vertical direction, whereas positively buoyant particles will diverge in the vertical direction. This vertical compressibility in the particle velocity field is associated with the same terms that lead to enhanced settling, but it is not the leading cause of particle dispersion.

To understand particle dispersion, we first consider the wave-averaged system. We initiate a cloud of settling particles near the surface and simulate their drift trajectories using (2.12). As shown in figure 5, a cloud of particles will experience dispersion in the horizontal direction in the wave-averaged sense despite the fact that the drift velocities do not show a horizontal divergence. The source of this effect is the influence of the initial particle positions. Since the drift velocities in (2.12) decay rapidly with depth in the same manner as the classical Stokes drift, particles with higher initial position will travel faster and further over time in the horizontal direction. In figure 5, the cloud of particles marked by their initial vertical position and tracked at regular time intervals expands in the horizontal direction because of differences in initial vertical positions while contracting slightly in the vertical direction as predicted by (2.14b). This mechanism for horizontal dispersion is stronger in deep-water waves (larger $kh$), for less dense particles (smaller $\gamma$), and for particles with less inertia (smaller ${St}$).

Figure 5. Wave-averaged particle positions computed using the wave-averaged drifts (2.12) for $ka = 0.1$, $kh = 1$, $\gamma = 1.1$ and ${St} = 0.1$. Particles are coloured according to their initial vertical position and particle positions are shown every 10 units of dimensionless time.

In the unaveraged system, figure 6 shows how particles can disperse when released at the same position continuously over time. The wave-averaged initial conditions depend on the wave phase at the initial particle position (Appendix B.3), and therefore particles released continuously over time will have different wave-averaged initial conditions, which will lead to horizontal dispersion. In this way, waves can disperse a continuous release of settling particles.

Figure 6. Particle trajectories computed with the linear drag model using (2.10) for particles released continuously at $(x,z) = (0,-0.2)$ (red asterisk) for $ka = 0.2$, $kh = 1$, $\gamma = 1.1$ and ${St} = 0.1$. Trajectories are coloured according to their initial starting time.

3.1. Numerical simulations

3.1.1. Particle trajectories

Linear drag: in the Stokes drag model, the entire problem is defined by four dimensionless parameters ($ka$, $kh$, $\gamma$ and ${St}$). Figure 7 compares particle trajectories obtained from the full dynamical equation with those obtained from the kinematic solution and wave-averaged drifts for different values of these four dimensionless parameters. In simulations using the full dynamical equations, the initial particle velocity is set to the sum of the local fluid velocity and the particle terminal velocity, though the solutions are not too sensitive to this choice. We can see that the kinematic solution tracks the full dynamical equation well in all cases, even when ${St} = O(1)$. The wave-averaged drifts also follow the unaveraged solutions faithfully in all cases, except when $ka$ is increased. This reflects the limitation of the drift velocities and initial condition corrections, which are correct to $O(\varepsilon ^{2})$ in the two-time expansion ( Appendix B). Therefore, discrepancies start to appear between the unaveraged and wave-averaged solutions as $\varepsilon$ increases, but remain small so long as $\varepsilon$ is not too large ($\varepsilon = 0.26$ in figure 7d).

Figure 7. Particle trajectories computed with the linear drag model using the full dynamical equation for particle motion ((2.7); solid blue line), the kinematic solution for particle motion ((2.10); dashed yellow line), and the wave-averaged drift velocities ((2.12); solid black line): (a) $ka=0.1$, $kh = 1$, $\gamma = 1.1$, ${St} = 0.1$; (b) $ka=0.1$, $kh = 3$, $\gamma = 1.1$, ${St} = 0.1$; (c) $ka=0.1$, $kh = 1$, $\gamma = 1.05$, ${St} = 0.75$; (d) $ka=0.2$, $kh = 1$, $\gamma = 1.1$, ${St} = 0.1$.

Nonlinear drag: in the Schiller–Naumann drag model, the terminal particle velocity must be calculated empirically and the drag coefficient becomes a nonlinear function of the slip velocity. This alters the particle dynamics significantly.

As an example, consider a typical laboratory experiment of a microplastic particle settling through waves: a particle of diameter 1 mm and density $1020\ {\rm kg}\ {\rm m}^{-3}$ settling through surface waves of amplitude 2.5 cm and period 1 s in water of depth 0.5 m. The fluid is of density $1000\ {\rm kg}\ {\rm m}^{-3}$ and kinematic viscosity $10^{-6}\ {\rm m}^{2}\ {\rm s}^{-1}$. This gives $\gamma = 1.02$, $ka = 0.10$ and $kh = 2.1$. Stokes drag predicts a terminal settling velocity of $1.1\ {\rm cm}\ {\rm s}^{-1}$ in quiescent water with an associated particle Reynolds number ${{\textit {Re}}}_{p,t} = 11$. The Schiller–Naumann drag model, which is expected to be more accurate in this case, predicts a terminal settling velocity of $0.7\ {\rm cm}\ {\rm s}^{-1}$ with an associated particle Reynolds number ${{\textit {Re}}}_{p,t} = 7$. Figure 8 shows the particle trajectories computed for this example using the full dynamical equation with linear drag (2.7), with nonlinear drag (2.8) and with nonlinear drag with constant Stokes number (2.8 with ${St_{SN}} = {St_{SN}}_{\!,t}$). The trajectories with nonlinear drag travel further in the horizontal direction as the particles settle slower through the water column because of a reduced settling velocity. The two nonlinear drag trajectories are similar to each other, suggesting that the slip velocity does not deviate too far from the terminal settling velocity. Figure 9 confirms this, showing that the particle slip velocity, the particle Reynolds number and the particle Stokes number all have small variations around their respective values based on the terminal velocity. In particular, ${St_{SN}}$ only fluctuates at most 2 % from its value computed with the terminal settling velocity ${St_{SN\!,t}}$.

Figure 8. Particle trajectories computed using the full dynamical equation with linear drag ((2.7); solid blue line), nonlinear drag ((2.8); solid red line) and nonlinear drag with constant Stokes number ((2.8) with ${St_{SN}} = {St_{SN}}_{\!,t}$; dashed green line) for $d_p = 1$ mm, $\rho _p = 1020\ {\rm kg}\ {\rm m}^{-3}$, $\rho = 1000\ {\rm kg}\ {\rm m}^{-3}$, $\nu = 10^{-6}\ {\rm m}^{2}\ {\rm s}^{-1}$, $a = 2.5\ {\rm cm}$, $\omega = 2{\rm \pi} \ {\rm s}^{-1}$, $h = 0.5$ m.

Figure 9. Slip velocity and derived quantities computed from a numerical solution of the full dynamical equation with nonlinear drag (2.8) for $d_p = 1$ mm, $\rho _p = 1020\ {\rm kg}\ {\rm m}^{-3}$, $\rho = 1000\ {\rm kg}\ {\rm m}^{-3}$, $\nu = 10^{-6} {\rm m}^{2}\ {\rm s}^{-1}$, $a = 2.5$ cm, $\omega = 2{\rm \pi} \ {\rm s}^{-1}$, $h = 0.5$ m. (a) Absolute slip velocity component normalised by $-v_s$, (b) ${{\textit {Re}}}_p$ and (c) ${St_{SN}}$ of the particle are plotted over dimensionless depth $z_p$. Corresponding terminal settling velocity values ${{\textit {Re}}}_{p,t}$ and ${St_{SN\!,t}}$ are marked with vertical grey lines.

If we consider the same particle settling in field conditions, longer wave periods in the field will result in lower Stokes numbers. Thus, the particle will come to local equilibrium with the flow faster in dimensionless terms and the fluctuations of ${St_{SN}}$ relative to the value of ${St_{SN\!,t}}$ will be even smaller in percentage terms. The approximation ${St_{SN}} = {St_{SN}}_{\!,t}$ is therefore even more justified in field conditions. Overall, this example serves to illustrate that the main effects of nonlinear drag are captured by how it alters terminal settling velocity and particle time scale to reach the terminal settling velocity in both laboratory and field conditions.

Given that a constant Stokes number approximation is reasonable, we now compare numerical simulations computed with the full dynamical equation (2.8) with the kinematic solution and wave-averaged particle drifts derived under the ${St_{SN}} = {St_{SN}}_{\!,t}$ approximation in figure 10. Similar to the linear drag cases, we see that the kinematic solution tracks the full dynamical equation for different values of the controlling dimensionless parameters and that the wave-averaged drifts also faithfully follow the unaveraged solutions, except when $ka$ becomes too large. Together, these results show that the kinematic solution and wave-averaged drifts are also applicable to scenarios with nonlinear drag, where the particle Reynolds number is significantly larger than unity.

Figure 10. Particle trajectories computed with the nonlinear drag model using the full dynamical equation for particle motion ((2.8); solid red line), the kinematic solution for particle motion ((2.10); dashed yellow line) and the wave-averaged drift velocities ((2.12); solid black line): (a) $ka=0.1$, $kh = 1$, $\gamma = 1.1$, ${St_{SN}}_{\!,t} = 0.1$, ${{\textit {Re}}}_{p,t} =212$; (b) $ka=0.1$, $kh = 3$, $\gamma = 1.1$, ${St_{SN}}_{\!,t} = 0.1$, ${{\textit {Re}}}_{p,t} = 55$; (c) $ka=0.1$, $kh = 1$, $\gamma = 1.05$, ${St_{SN}}_{\!,t} = 0.5$, ${{\textit {Re}}}_{p,t} = 1093$; (d) $ka=0.2$, $kh = 1$, $\gamma = 1.1$, ${St_{SN}}_{\!,t} = 0.1$, ${{\textit {Re}}}_{p,t} = 212$.

3.1.2. Settling velocity

In order to further compare the derived drift expressions with the full dynamical system, we consider the particle settling velocity directly. To achieve this comparison, the particle motion from numerical simulations needs to be wave averaged to remove the oscillations that are an order of magnitude larger than the drift velocities of interest. For settling particles, this averaging procedure is non-trivial for two reasons: (i) the oscillation frequency of a Lagrangian particle is different from that of the Eulerian flow because of a Doppler-shifting effect of the Stokes drift (Longuet-Higgins Reference Longuet-Higgins1986), and (ii) the oscillation frequency evolves as the particle settles through the wave field because the Stokes drift decays with depth. To account for this variation in frequency, we use a discrete Hilbert transform (Huang & Wu Reference Huang and Wu2008) to find the desired phase of the particle motion in each oscillation cycle. Figure 11(a) shows a particle trajectory with particle positions marked at the top and bottom of each oscillation using this procedure. The wave-averaged value for some quantity $f$ following the particle motion can then be found by time averaging between corresponding phases (which occur at times $t_a$ and $t_b$),

(3.1)\begin{equation} \bar{f}= \frac{1}{t_b-t_a}\int_{t_a}^{t_b} f(t) \,{\rm d}t. \end{equation}

Figure 11. Particle trajectory and wave-averaged drifts for $kh=1$, $ka=0.1$, $\gamma =1.05$, ${St}_{SN,t}=0.065$, $Re_{p,t}=6.5$. (a) Particle trajectories computed with the nonlinear drag model using the full dynamical equation for particle motion (2.8) with constant phase points marked with $\ast$. (b) Wave-averaged vertical velocity with starting phase points marked. (c) Twice wave-averaged vertical velocity of the trajectory $-\overline {\overline {v_z}}/v_s$.

Using this method, we average the particle's vertical velocity and position over its trajectory, and we plot the resulting wave-averaged quantities $\overline {v_z}$ and $\overline {z_p}$ and in figure 11(b). From inspection, we find an initial condition dependence exists: $\overline {v_z}$ is a function of the starting phase of the averaging. Its value oscillates around unity because the particle oscillatory motion does not settle uniformly; rather, the settling velocity is faster when averaged crest-to-crest compared with trough-to-trough. In other words, the tops of the particle orbitals move downward faster than the bottoms of the particle orbitals. This result suggests that waves can both enhance and reduce the particle settling velocity depending on where the wave-averaging process begins. The enhancement and reduction, however, are not symmetric. By wave averaging the signal once more, we obtain $\overline {\overline {v_z}}$ which is now only a function of $\overline {\overline {z_p}}$ (figure 11c). At this level of averaging, the initial condition dependence goes away, and we recover the enhanced vertical drift, which closely matches the analytical prediction of (2.12b). Only by wave averaging the particle velocities twice are we able to recover the vertical drift velocities and confirm that the particle settling velocity is enhanced by waves.

3.2. Laboratory experiments

Finally, we compare our theoretical results with the experimental data of settling spheres under surface waves reported in Clark et al. (Reference Clark, DiBenedetto, Ouellette and Koseff2020), where we have conducted a re-analysis of the previously reported data. The experiments investigated spheres with diameter 2.96 mm settling in three different wave conditions in water of depth 41.5 cm: $a = 3.5$ cm, $\omega = 2{\rm \pi} \ {\rm rad}\ {\rm s}^{-1}$ (‘shallow’; $ka = 0.15$, $kh = 1.78$), $a = 3.3$ cm, $\omega = 3{\rm \pi} \ {\rm rad}\ {\rm s}^{-1}$ (‘intermediate’; $ka=0.29$, $kh=3.7$) and $a = 2.3$ cm, $\omega = 4{\rm \pi} \ {\rm rad}\ {\rm s}^{-1}$ (‘deep’; $ka=0.32$, $kh = 5.9$). The flow was laminar and well approximated by linear wave theory (DiBenedetto et al. Reference DiBenedetto, Koseff and Ouellette2019). The terminal settling velocity of the spheres in quiescent fluid was $2.90\ {\rm cm}\ {\rm s}^{-1}$. In comparing the experimental data with our analytical results, we calculated a terminal settling velocity of $2.86\ {\rm cm}\ {\rm s}^{-1}$ using the Schiller–Naumann model, with specific gravity $\gamma =1.025$ and kinematic viscosity $\nu =10^{-6}\ {\rm m}^{2}\ {\rm s}^{-1}$. Thus, these particles exist in the nonlinear drag regime where ${{\textit {Re}}}_p=86$, ${St_{SN}}_{\!,t}=[1.1, 1.7,2.2]$, and dimensionless settling velocity $v_s=[0.019, 0.028, 0.037]$ for the shallow-, intermediate- and deep-water wave conditions, respectively.

The experimental data were measured with particle tracking velocimetry and showed enhanced settling under waves. The data were originally reported in terms of an ensemble mean of the instantaneous particle vertical velocities as a function of a dimensionless ‘velocity scale reflective of the local vertical flow velocities’ and the results appear to show an increase with this quantity (see their figure 3). This analysis allowed Clark et al. (Reference Clark, DiBenedetto, Ouellette and Koseff2020) to combine data across wave conditions, scaled by relative wave strength. To compare with our drift predictions, we average the instantaneous vertical velocities of the particles in a similar way, but now as a function of both wave conditions and dimensionless depth. We note that this ensemble average is distinct from the trajectory average that leads to the drift velocity predictions, and therefore we do not expect the laboratory results to be identical to the theory, but we do expect both to reflect the settling enhancement in a similar manner.

Figure 12 shows a comparison between the averaged experimental data and the theory, where we see that the theory predicts the correct overall trend with depth and wave conditions, but the experimental data show a larger velocity enhancement than predicted by (2.12b). These discrepancies between the theory and the data could be due to multiple reasons. First, the experimental results are averaged in a different way than the theory. Additionally, the experiments are performed at relatively large $ka$ values where the theory starts to break down, especially near the surface. Finally, other approximations made in our theoretical analysis, e.g. neglecting the Basset history and assuming that the drag coefficient in the unsteady system can be approximated as the steady drag coefficient at a particular ${{\textit {Re}}}_p$ value, could also play a role.

Figure 12. Experimental data (symbols) from Clark et al. (Reference Clark, DiBenedetto, Ouellette and Koseff2020) compared with predicted vertical drift (2.12b) (solid lines). Different colours represent different experimental cases: ‘shallow’: $ka = 0.15$, $kh = 1.77$; ‘intermediate’: $ka=0.30$, $kh=3.8$;‘deep’: $ka=0.37$, $kh = 6.7$. Error bars denote 95 % confidence intervals which were found by bootstrapping the experimental data.