1. Introduction
The dynamics of particulate matter in surface gravity waves determine the transport and dispersion of natural particles, such as sediment, ice crystals and plankton, as well as synthetic particles, such as microplastics. For small, neutrally buoyant particles that can be treated as tracers of the fluid motion, there is a net horizontal transport in the wave-averaged sense known as Stokes drift (Stokes Reference Stokes1847; van den Bremer & Breivik Reference van den Bremer and Breivik2017). However, the dynamics of most particles in the environment deviate from the tracer limit because of their inertia, stemming from finite size and density difference with the fluid. These deviations alter particle transport in surface waves from the classical Stokes drift.
Recent work in inertial particle transport in surface waves, much of which has been motivated by marine litter, has shown that altered forms of Stokes drift can arise from several mechanisms. For positively buoyant particles that float on the surface, Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021) showed that the horizontal drift velocity is enhanced for larger particles because of dynamic buoyancy forces. For weakly inertial particles suspended in the water column, Alsina, Jongedijk & van Sebille (Reference Alsina, Jongedijk and van Sebille2020) found that horizontal particle drift was reasonably well predicted by the classical Stokes drift velocity in laboratory experiments.
For negatively buoyant particles, several works have examined how particle inertia alters the horizontal and vertical drifts (Eames Reference Eames2008; Santamaria et al. Reference Santamaria, Boffetta, Afonso, Mazzino, Onorato and Pugliese2013; Bakhoday-Paskyabi Reference Bakhoday-Paskyabi2015). In particular, Santamaria et al. (Reference Santamaria, Boffetta, Afonso, Mazzino, Onorato and Pugliese2013) predicted a vertical drift that enhances particle settling velocity due to the combined effects of particle inertia and the flow field. This effect can be thought of as a vertical Stokes drift. Laboratory experiments have similarly measured enhanced settling (Clark et al. Reference Clark, DiBenedetto, Ouellette and Koseff2020; De Leo et al. Reference De Leo, Cutroneo, Sous and Stocchino2021), although these experiments have been outside the regime of the low particle Reynolds number assumed in theoretical works, making it difficult to perform a direct comparison. Additionally, non-spherical particles have been shown to experience an angular analogue to Stokes drift that causes particles to adopt preferred orientations (DiBenedetto & Ouellette Reference DiBenedetto and Ouellette2018; DiBenedetto, Koseff & Ouellette Reference DiBenedetto, Koseff and Ouellette2019) that alter their trajectories in shape-dependent ways (DiBenedetto, Ouellette & Koseff Reference DiBenedetto, Ouellette and Koseff2018; Clark et al. Reference Clark, DiBenedetto, Ouellette and Koseff2020). The combined findings of these studies show the value of analysing particle dynamics from a Stokes drift framework.
Previous theory on inertial particle transport has generally adopted the well-known equations of motion for small, rigid spheres whose motion relative to the fluid is characterised by small Reynolds number (Maxey & Riley Reference Maxey and Riley1983). If the assumption of small particle Reynolds number is relaxed, the drag force deviates from the Stokes drag law and becomes a nonlinear function of the slip velocity. The effects of this nonlinear drag on particle settling have been explored in vertically oscillating fluid columns, which partially mimic the oscillatory flow of surface waves. In the Newtonian drag regime, where the drag coefficient is a constant, it was found that particles have a reduced settling velocity relative to their terminal settling velocity (Ho Reference Ho1964; Nielsen Reference Nielsen1984; Hwang Reference Hwang1985; Ikeda & Yamasaka Reference Ikeda and Yamasaka1989). However, it is unclear whether these findings from vertically oscillating columns translate directly to surface waves, and if they do, how this reduction in settling velocity in the Newtonian drag regime competes with the mechanism that produces enhanced settling in the linear drag regime.
Here, we reconsider the settling of inertial spheres in surface waves. While previous analytical work by Eames (Reference Eames2008) and Santamaria et al. (Reference Santamaria, Boffetta, Afonso, Mazzino, Onorato and Pugliese2013) examined this problem in deep-water waves in the limit of small particle Reynolds number, we expand the scope by considering waves in arbitrary depth and particles whose motion relative to the fluid need not be in the small Reynolds number limit. Starting from the dynamical equation of particle motion, we first derive a kinematic solution to particle motion that gives the particle velocity as a function of position within the flow and then use a two-time expansion to isolate the wave-induced particle oscillatory motion from the wave-averaged particle drift. We also describe a mechanism for horizontal dispersal of particles that are settling vertically. We compare our theoretical results with numerical simulations and laboratory experiments, being careful to describe procedures to accurately translate between the wave-resolving and wave-averaged systems. Our results show that inertial particles that are denser than the fluid undergo oscillatory motion that has a smaller amplitude and a phase lag relative to tracer particles. We also find that such particles experience an enhancement in the settling velocity, which stems from a dynamical part related to particle inertia and a kinematic part related to how the particles sample the flow field, akin to a vertical Stokes drift. The enhanced settling velocity can be expressed in terms of the classical Stokes drift. For larger and denser particles, we use a well-known model of the drag coefficient in the transitional regime of Reynolds number and make a simplification that allows analytical solutions of the particle motion to be extended to particles with nonlinear drag. The theory compares well with numerical simulations and reasonably well with laboratory experiments.
The remainder of this paper is organised as follows. We begin with the equations of motion and derive the theoretical results in § 2, with the comparisons with simulations and experiments presented in § 3 and conclusions in § 4.
2. Particle motion in surface gravity waves
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
$$\begin{gather} u_x^{\prime} = \frac{gka}{\omega} \frac{\cosh (k(z^{\prime}+h))}{\cosh kh} \cos (kx^{\prime} - \omega t^{\prime}), \end{gather}$$
$$\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
$$\begin{gather} u_x = \varepsilon \frac{\cosh (z+kh)}{\cosh kh} \cos (x - t), \end{gather}$$
$$\begin{gather}u_z = \varepsilon \frac{\sinh (z+kh)}{\cosh kh} \sin (x - t), \end{gather}$$where
\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
$$\begin{gather} \frac{\mathrm{d}\kern0.06em \boldsymbol {x_p}^{\prime}}{\mathrm{d}t^{\prime}} =\boldsymbol{v}^{\prime} \end{gather}$$
$$\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
$$\begin{gather} \frac{\mathrm{d}\kern0.06em\boldsymbol{x_p}}{\mathrm{d}t} =\boldsymbol{v} , \end{gather}$$
$$\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
\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
\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
\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
\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)
$$\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}$$
$$\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
$$\begin{gather} A = \sqrt{\left[1 - {St}^{2} \left( 1 - \beta \right) \right]^{2} + \left[{St}\left( 1 - \beta \right)\right]^{2}}, \end{gather}$$
$$\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
$$\begin{gather} v_{x\text{-}drift} =\frac{A^{2}}{1+v_s^{2}} u_{{SD}} , \end{gather}$$
$$\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)
\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
$$\begin{gather} \frac{\partial v_{x\text{-}drift}} {\partial x_{p0}} = 0 , \end{gather}$$
$$\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. Comparisons with numerical simulations and laboratory experiments
Having derived a kinematic solution for particle motion (2.10) and the associated wave-averaged particle drift velocities (2.12), we now compare them with numerical simulations of the full dynamical equation of particle motion (2.7), (2.8) and laboratory experiments.
When comparing particle trajectories between the unaveraged system and the wave-averaged system, it is important to project the initial conditions from the unaveraged system to the averaged system (Appendix B.3). Even a change as small as 
$O(\varepsilon ^{2})$ in the initial vertical position of the particle in the averaged system can lead to a large difference in particle trajectories as shown in figures 5 and 6.
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$),
\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.
4. Conclusions
From the dynamical equation of inertial spheres in surface waves, we have found a kinematic solution that gives the particle velocity as a function of its position within the background flow field (2.10) and a solution for the particle drift velocity that gives the horizontal and vertical particle motion in the wave-averaged sense (2.12). While the former can be used in wave-resolving models, the latter can be used in wave-averaged or spectral models to compute inertial particle motion. These solutions are derived using expansions that assume small Stokes number and small wave amplitude, with no restrictions on the relative water depth. Comparisons with numerical simulations of the full dynamical equation show that the solutions accurately predict particle motion up to 
${St} = O(1)$ and 
$\varepsilon = O(10^{-1})$, suggesting that they are applicable within the confines of linear wave theory even when the particle inertia is significant.
The kinematic solution shows that negatively buoyant inertial particles move in orbitals that are smaller and lag behind fluid parcels in time (figure 2). This effect is especially noticeable for 
${St} > 0.1$. The particle drift velocity solution shows that the horizontal drift is proportional to the Stokes drift for fluid tracers, but reduced in magnitude for negatively buoyant inertial particles (figure 3). The vertical drift velocity solution shows that inertial particles experience enhanced settling in surface waves (figure 4). This enhancement relative to the terminal settling velocity is both a kinematic and a dynamic effect. While the kinematics create a bias for the particle to sample downward wave-induced velocities, akin to a vertical Stokes drift, the dynamics enhance settling via a coupling between the particle inertia and fluid inertia. These conclusions are supported by theory (§ 2.3), numerical simulations (§ 3.1) and experimental evidence (§ 3.2).
Additionally, we find that a cloud of particles released together or a continuous release of particles at a single location will undergo horizontal dispersion because of the differences in the particles’ wave-averaged initial vertical positions (§ 2.3). Particles with a higher initial wave-averaged position travel faster and further than those with a lower initial position, which leads to particle trajectories diverging in the horizontal direction.
Finally, we have examined different regimes of drag, including linear (Stokes) drag and nonlinear drag (as parameterised by the Schiller–Naumann model). A simplifying assumption in the nonlinear drag formulation allows us to extend the kinematic solution and wave-averaged drift velocity solution to also include the nonlinear drag regime where the particle Reynolds number can be up to 
${{\textit {Re}}}_p = O(10^{3})$. This simplification can be interpreted as assuming that the particle relaxation time scale is independent of the flow such that 
$\tau _{p} = d_p^{2} (C_m + \gamma )/[18\nu \, f({{\textit {Re}}}_p)] \approx d_p^{2} (C_m + \gamma )/[18\nu \, f({{\textit {Re}}}_{p,t})]$ or that the drag model is ‘linearised’ such that 
$C_D = (24/{{\textit {Re}}}_p)\, f({{\textit {Re}}}_p) \approx (24/{{\textit {Re}}}_p)\, f({{\textit {Re}}}_{p,t})$, where 
${{\textit {Re}}}_{p,t}$ is the particle Reynolds number at terminal velocity in quiescent fluid determined from the nonlinear drag model. Comparisons with numerical simulations show that this simplification is effective and that the kinematic solution and wave-averaged drifts apply equally well to the nonlinear drag regime as long as the terminal settling velocity and Stokes number are based on the nonlinear drag model.
Acknowledgements
We acknowledge helpful discussions with K. Ma, J.-L. Thiffeault, N. Ouellette and J. Koseff, as well as constructive comments from three anonymous referees.
Funding
N.P. acknowledges support from the Freshwater Collaborative of Wisconsin (grant number T2-5/20-06) and the U.S. National Science Foundation (grant number: OCE-2048676 ); L.K.C. acknowledges support from the U.S. National Science Foundation (grant number: CBET-1706586).
Declaration of interests
The authors report no conflict of interest.
Author contributions
N.P. and M.H.D. derived the theory and performed the simulations. L.K.C. performed the experiments and L.K.C. and M.H.D. analysed the data. N.P., L.K.C and M.H.D. wrote the paper.
Appendix A. Particle motion for small 
${St}$ number
In this appendix, we find the kinematic solution for particle motion in the limit of small Stokes number. Expanding the particle velocity in the Stokes number
\begin{equation} \boldsymbol{v} = \boldsymbol{v}_0 + {St} \, \boldsymbol{v}_1 + {St}^{2} \, \boldsymbol{v}_2 + \dots, \end{equation}and substituting it into (2.9), which is re-written as
\begin{equation} \boldsymbol{v} = \boldsymbol{u} + {\boldsymbol{v}_s} - {St} \, \frac{{\mathrm{d}\boldsymbol{v}}}{{\mathrm{d}t}} + {St} \,\beta \frac{{\mathrm{D} \boldsymbol{u}}}{{\mathrm{D} t}}, \end{equation}and collecting terms of different orders gives
\begin{equation} \boldsymbol{v} = \boldsymbol{u} + \boldsymbol{v}_s + {St} \left[ - \frac{\mathrm{d} \boldsymbol{u}}{\mathrm{d} t} + \beta \frac{\mathrm{D} \boldsymbol{u}}{\mathrm{D} t} \right] + {St}^{2} \left[ \frac{\mathrm{d}^{2} \boldsymbol{u}}{\mathrm{d} t^{2}} - \beta \frac{\mathrm{d} }{\mathrm{d} t} \left( \frac{\mathrm{D} \boldsymbol{u}}{\mathrm{D} t} \right) \right] + O({St}^{3}). \end{equation} By expanding 
${\mathrm {d} \boldsymbol {u}}/{\mathrm {d} t}$
\begin{equation} \frac{\mathrm{d} \boldsymbol{u}}{\mathrm{d} t} = \frac{{\partial} \boldsymbol{u}}{{\partial} t} + \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u} = \frac{\mathrm{D} \boldsymbol{u}}{\mathrm{D} t} + \boldsymbol{v}_s \boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u} + O(\varepsilon^{2} {St}), \end{equation}
and substituting it into (A3) we can find explicit expressions for the particle velocity. The expressions depend on the relative magnitudes of 
${St}$, 
$\varepsilon$ and 
$v_s$. We take 
$\varepsilon$ and 
${St}$ to be of the same order and 
$v_s \ll St$. This last choice is driven by an analysis of the ratio of the dimensionless settling velocity and the particle Stokes number, given by
\begin{equation} \frac{v_s}{{St}} = \frac{1}{\tanh kh} \frac{v_s'}{g \tau_p}, \end{equation}
where 
$v_s'$ is the dimensional terminal settling velocity and 
$\tau _p$ is the particle time scale. Figure 13 shows that 
${v_s'}/({g \tau _p})$ is typically small, and hence, 
${v_s} < {{St}}$ for most applications as long as the waves are not in extremely shallow water.
Figure 13. Ratio of dimensionless terminal velocity to Stokes number for deep-water waves 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}$.
From (A3), the particle velocity is given by
\begin{equation} \boldsymbol{v} = \boldsymbol{u} + \boldsymbol{v}_s - {St}(1 - \beta)\frac{\mathrm{D} \boldsymbol{u}}{\mathrm{D} t} + {St}^{2} (1 - \beta) \frac{\partial^{2} \boldsymbol{u}}{\partial t^{2}} + O(\varepsilon^{2}{St}^{2},\varepsilon {St}^{3}), \end{equation}which can be evaluated using (2.2) to give
$$\begin{gather} {v_x} = \varepsilon A \frac{\cosh (z+kh)}{\cosh kh} \cos\left(x - t + \phi \right) + \varepsilon^{2} {St}(1 - \beta) \frac{\sin (2(x - t))}{2 \cosh^{2} kh} , \end{gather}$$
$$\begin{gather}{v_z} ={-} v_s + \varepsilon A \frac{\sinh (z+kh)}{\cosh kh} \sin\left(x - t + \phi \right) - \varepsilon^{2} {St}(1 - \beta) \frac{\sinh (2(z + kh))}{2 \cosh^{2} kh}. \end{gather}$$The particle orbital motion at leading order has modified amplitude and phase shift relative to the motion of fluid parcels given by
$$\begin{gather} A = \sqrt{\left[1 - {St}^{2} \left( 1 - \beta \right) \right]^{2} + \left[{St}\left( 1 - \beta \right)\right]^{2}}, \end{gather}$$
$$\begin{gather}\phi = \arctan \left[ \frac{{St} \left( 1 - \beta \right)}{1 - {St}^{2} \left( 1 - \beta \right)} \right]. \end{gather}$$Appendix B. Wave-averaged particle motion
In this appendix, we employ a two-time expansion on the particle motion found in Appendix A to remove the wave-induced oscillations, find the wave-averaged particle drift velocities and find suitable initial conditions for the wave-averaged particle motion.
B.1. Two-time expansion
To perform a two-time expansion to separate the fast wave-induced oscillations from the slow wave-averaged behaviour, we look for a solution of the form
\begin{equation} \boldsymbol{x}_p (t; \varepsilon) = \boldsymbol{x}_p (\tau, T; \varepsilon), \quad \text{with} \ \tau = t, \ T = \varepsilon^{2} t. \end{equation} The slow time scale is taken to be order 
$\varepsilon ^{2}$ since the Stokes drift and the dynamical settling enhancement for inertial settling particles both arise at the second order. With this, the particle velocity (A7) is given by
\begin{align} \partial_{\tau}{x_p} + \varepsilon^{2} \partial_{T}{x_p} &= \varepsilon A \frac{\cosh (z_p+kh)}{\cosh kh} \cos (x_p - \tau + \phi) \nonumber\\ & \quad + \varepsilon^{2} {St}(1 - \beta) \frac{\sin (2(x_p - \tau))}{2 \cosh^{2} kh} \end{align}
\begin{align} \partial_{\tau}{z_p} + \varepsilon^{2} \partial_{T}{z_p} &={-} v_s + \varepsilon A \frac{\sinh (z_p+kh)}{\cosh kh} \sin(x_p - \tau + \phi) \nonumber\\ & \quad - \varepsilon^{2} {St}(1 - \beta) \frac{\sinh (2(z_p + kh))}{2 \cosh^{2} kh}. \end{align} We expand the particle position about the small parameter 
$\varepsilon$ that characterises the dimensionless wave amplitude,
\begin{equation} \boldsymbol{x}_p (\tau, T) = \boldsymbol{x}_{p0} (\tau, T) + \varepsilon \boldsymbol{x}_{p1} (\tau, T) + \dots, \end{equation}
treating 
$v_s$ and 
${St}$ as constants. The order at which the terminal settling velocity 
$v_s$ enters the problem depends on the magnitude of the settling number
\begin{equation} {Sv} = \frac{v_s}{\varepsilon}, \end{equation}
which is the ratio of the terminal particle settling velocity to the fluid velocity scale. We will assume 
$v_s = O(1)$ (
${Sv} \gg 1$) so that settling enters the problem at leading order. Although the settling number is not necessarily always large, following this assumption ensures that we capture how particle settling kinematics affect wave-averaged drift velocities.
Substituting this expansion into (B2) and collecting terms of the same order gives the following.
At 
$\varepsilon ^{0}$
$$\begin{gather} \partial_{\tau} x_{p0} = 0 , \end{gather}$$
$$\begin{gather}\partial_{\tau} z_{p0} ={-}v_s. \end{gather}$$
The leading-order solution is a constant settling superimposed on functions of the slow time scale (
$X_p (T), Z_p (T$))
$$\begin{gather} x_{p0} = X_p (T) , \end{gather}$$
$$\begin{gather}z_{p0} = Z_p(T)-v_s\tau. \end{gather}$$ At 
$\varepsilon ^{1}$
$$\begin{gather} \partial_{\tau} x_{p1} = A\frac{\cosh (z_{p0} +kh)}{\cosh kh} \cos(x_{p0} - \tau + \phi) , \end{gather}$$
$$\begin{gather}\partial_{\tau} z_{p1} = A \frac{\sinh (z_{p0} +kh)}{\cosh kh} \sin(x_{p0} - \tau + \phi) , \end{gather}$$whose solution is given by
\begin{align} x_{p1} &={-}\frac{A}{\cosh kh (1 + v_s^{2})} \left[\cosh(z_{p0}+ kh) \sin(x_{p0} - \tau + \phi)\right.\nonumber\\ &\quad \left. + v_s \sinh(z_{p0} + kh) \cos(x_{p0} - \tau + \phi) \right] \end{align}
\begin{align} z_{p1} &= \frac{A}{\cosh kh (1 + v_s^{2})} \left[ \sinh(z_{p0} + kh) \cos(x_{p0} - \tau + \phi) \right.\nonumber\\ &\left. \quad - v_s \cosh(z_{p0} + kh) \sin(x_{p0} - \tau + \phi) \right]. \end{align} At 
$\varepsilon ^{2}$
\begin{align} \partial_{\tau} x_{p2} + \partial_{T} X_p &= \frac{A^{2}}{\cosh^{2} kh(1 + v_s^{2})} \left[ \sinh^{2}( z_{p0} + kh) +\sin^{2}(x_{p0} - \tau + \phi)\right] \nonumber\\ & \quad + {St}(1-\beta) \frac{\sin (2(x_{p0} - \tau))}{2\cosh^{2} kh} \end{align}
\begin{align} \partial_{\tau} z_{p2} + \partial_{T} Z_p &={-}\frac{v_sA^{2}}{\cosh^{2} kh(1 + v_s^{2})} \left[ \sinh^{2}(z_{p0} + kh) +\sin^{2}(x_{p0} - \tau + \phi)\right] \nonumber\\ & \quad - {St}(1-\beta) \frac{\sinh (2(z_{p0} + kh))}{2\cosh^{2} kh}. \end{align}The solvability condition at this order is given by averaging this equation over a wave period
$$\begin{gather} \partial_{T} X_p =\frac{A^{2}}{2\cosh^{2} kh(1 + v_s^{2})} \cosh(2( z_{p0} + kh)) , \end{gather}$$
$$\begin{gather}\partial_{T} Z_p ={-}\frac{v_sA^{2}}{2\cosh^{2} kh(1 + v_s^{2})} \cosh (2( z_{p0} + kh)) - \frac{{St}(1-\beta)}{2\cosh^{2} kh} \sinh (2(z_{p0} + kh)). \end{gather}$$
Note that we have made a simplifying approximation to obtain explicit expressions: we have assumed that the constant settling is small enough such that changes in 
$z_{p0} = Z_p - v_s \tau$ over a wave period can be neglected. The validity of this simplification is investigated by comparisons with numerical solutions.
To find the solution at the next order, we substitute (B10) into the second-order equation to get
$$\begin{gather} \partial_{\tau} x_{p2} ={-} \frac{A^{2}}{2\cosh^{2} kh(1 + v_s^{2})}\left[\cos(2(x_{p0} - \tau + \phi )) + {St}(1-\beta)\sin (2(x_{p0} - \tau)) \right], \end{gather}$$
$$\begin{gather}\partial_{\tau} z_{p2} = \frac{v_s A^{2}}{2\cosh^{2} kh(1 + v_s^{2})} \cos(2(x_{p0} - \tau + \phi )). \end{gather}$$The solution is given by
$$\begin{gather} x_{p2} = \frac{A^{2}}{4\cosh^{2} kh(1 + v_s^{2})} \left[\sin(2(x_{p0} - \tau + \phi )) + {St}(1-\beta) \cos (2(x_{p0} - \tau))\right], \end{gather}$$
$$\begin{gather}z_{p2} ={-}\frac{v_s A^{2}}{4\cosh^{2} kh(1+v_s^{2})} \sin(2(x_{p0} - \tau + \phi)). \end{gather}$$B.2. Wave-averaged particle drift velocities
From (B10), we can obtain the leading-order wave-averaged particle velocities that include the constant settling component by using 
$\partial _t = \partial _{\tau } + \varepsilon ^{2} \partial _T$
$$\begin{gather} \partial_{t} x_{p0} =\varepsilon^{2} \frac{A^{2}}{2\cosh^{2} kh(1+v_s^{2})}\cosh(2( z_{p0} + kh)) , \end{gather}$$
$$\begin{gather}\partial_{t} z_{p0} ={-}v_s - \varepsilon^{2} \frac{v_s A^{2}}{2\cosh^{2} kh(1 + v_s^{2})}\cosh (2( z_{p0} + kh)) - \varepsilon^{2} \frac{{St}(1-\beta)}{2\cosh^{2} kh} \sinh (2(z_{p0} + kh)). \end{gather}$$These are the particle drift velocities and are related to the classical Stokes drift velocity of tracer particles
\begin{equation} u_{{SD}} = \varepsilon^{2} \frac{\cosh(2(z_{p0} + kh))}{2\cosh^{2} kh}. \end{equation}
Hence, we refer to them as 
$v_{x\text {-}drift}$ and 
$v_{z\text {-}drift}$, and we re-write them in terms of the Stokes drift
$$\begin{gather} v_{x\text{-}drift} =\frac{A^{2}}{1+v_s^{2}} u_{{SD}} , \end{gather}$$
$$\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}$$B.3. Initial conditions for wave-averaged particle motion
Equation (B15) can be used to compute the wave-averaged trajectories of particles, but this requires that the initial conditions first be transformed into the wave-averaged system. To find this transformation, we expand the wave-averaged initial condition as
\begin{equation} \boldsymbol{x}_{p0} (0) = \boldsymbol{x}_{p0,0} + \varepsilon \boldsymbol{x}_{p0,1} + \varepsilon^{2} \boldsymbol{x}_{p0,2} + \dots, \end{equation}
and substitute it into the expansion (B3) at the initial time 
$t=0$
\begin{equation} \boldsymbol{x}_p (0) = \boldsymbol{x}_{p0} (0) + \varepsilon \boldsymbol{x}_{p1} (0) + \varepsilon^{2} \boldsymbol{x}_{p2} (0) + \dots, \end{equation}and collect terms of the same order.
At 
$\varepsilon ^{0}$, we get 
$\boldsymbol {x}_{p0,0} = \boldsymbol {x}_p (0)$. At 
$\varepsilon ^{1}$, we get
\begin{align} {x}_{p0,1} &= \frac{A}{\cosh kh (1 + v_s^{2})} \cosh (z_{p0,0} + kh ) \sin (x_{p0,0} + \phi ) \nonumber\\ & \quad + \frac{v_s A }{\cosh kh (1 + v_s^{2})} \sinh (z_{p0,0} + kh ) \cos (x_{p0,0} + \phi ) \end{align}
\begin{align} {z}_{p0,1} &={-}\frac{A}{\cosh kh (1 + v_s^{2})} \sinh (z_{p0,0} + kh) \cos (x_{p0,0} + \phi ) \nonumber\\ & \quad + \frac{v_s A}{\cosh kh (1 + v_s^{2})} \cosh (z_{p0,0} + kh ) \sin (x_{p0,0} + \phi ). \end{align} At 
$\varepsilon ^{2}$, we get
$$\begin{gather} {x}_{p0,2} = \frac{A^{2}}{4\cosh^{2} kh(1+v_s^{2})}\sin(2(x_{p0,0} + \phi )) - \frac{{St}(1- \beta)}{4\cosh^{2} kh} \cos (2 x_{p0,0}) , \end{gather}$$
$$\begin{gather}{z}_{p0,2} = \frac{A^{2}}{2\cosh^{2} kh(1+v_s^{2})}\sinh(2(z_{p0,0} + kh )) + \frac{v_s A^{2}}{4 \cosh^{2} kh(1+v_s^{2})}\sin(2(x_{p0,0} + \phi)). \end{gather}$$ Thus, for a given set of initial conditions 
$\boldsymbol {x}_p (0)$, the initial conditions for the wave-averaged system correct to 
$O(\varepsilon ^{2})$ are
\begin{align} {x}_{p0}(0) &= {x}_{p}(0) + \varepsilon \frac{A}{\cosh kh(1 + v_s^{2})} \cosh (z_{p} (0) + kh) \sin (x_{p} (0) + \phi ) \nonumber\\ & \quad + \varepsilon \frac{v_s A }{\cosh kh(1 + v_s^{2})}\sinh (z_{p} (0) + kh) \cos (x_{p} (0) + \phi) \nonumber\\ & \quad + \varepsilon^{2} \frac{A^{2}}{4\cosh^{2} kh(1+v_s^{2})} \sin(2(x_{p}(0) + \phi)) \nonumber\\ & \quad - \varepsilon^{2} \frac{{St}(1-\beta)}{4\cosh^{2} kh} \cos (2 x_{p}(0)) \end{align}
\begin{align} {z}_{p0}(0) &= {z}_{p}(0) - \varepsilon \frac{A}{\cosh kh(1 + v_s^{2})}\sinh (z_{p} (0) + kh) \cos (x_{p} (0) + \phi ) \nonumber\\ & \quad + \varepsilon \frac{v_s A }{\cosh kh(1 + v_s^{2})}\cosh (z_{p} (0) + kh ) \sin (x_{p} (0) + \phi ) \nonumber\\ & \quad + \varepsilon^{2} \frac{A^{2}}{2\cosh^{2} kh(1+v_s^{2})}\sinh(2(z_{p}(0) + kh )) \nonumber\\ & \quad + \varepsilon^{2} \frac{v_s A^{2}}{4\cosh^{2} kh(1+v_s^{2})}\sin(2(x_{p}(0) + \phi)). \end{align}
Open access
