1. Introduction
In the last half-century, large concentrations of plastic have polluted the oceans, with harmful effects on marine wildlife and potentially on human health (Cole et al. Reference Cole, Lindeque, Halsband and Galloway2011; Cózar et al. Reference Cózar2014; Ostle et al. Reference Ostle, Thompson, Broughton, Gregory, Wootton and Johns2019). Plastic pollution may have lasting impact, noting that it has been estimated that plastic may take hundreds or thousands of years for plastic to decay in the ocean (Cole et al. Reference Cole, Lindeque, Halsband and Galloway2011), although such estimates are subject to considerable uncertainty (Ward & Reddy Reference Ward and Reddy2020). Floating plastic debris is transported and dispersed by three key mechanisms: currents, wind and waves (van Sebille et al. Reference van Sebille2020). This paper will investigate wave-induced transport.
To leading order and in deep water, the Lagrangian motion induced by waves takes the form of circular orbits with Lagrangian particles following these orbits in a periodic fashion. The imbalance between the forward orbital velocity when under the crest and backward orbital velocity when under the trough, caused by the decay in velocity with depth, and the fact that particles spend more time under the forward-moving crest than under the backward-moving trough results in orbits that do not close, i.e. a Lagrangian-mean drift, known as Stokes drift (Stokes Reference Stokes1847). Stokes drift in deep water is proportional to the square of wave steepness and decays with depth at twice the rate of the oscillatory water particle velocity (see e.g. the review by van den Bremer & Breivik Reference van den Bremer and Breivik2017). Ocean surface gravity waves are driven by wind, and thus Stokes drift has often been assumed to be locally proportional to the wind forcing (Weber Reference Weber1983). However, waves are slow to build and, once established as swell, waves can travel long distances with little dispersion (Hanley, Belcher & Sullivan Reference Hanley, Belcher and Sullivan2010; Ardhuin et al. Reference Ardhuin2019), and so their magnitude is not always proportional to the local wind forcing. Wave models, such as WaveWatch III (The WaveWatch III Development Group 2016), can be used to predict Stokes drift (Webb & Fox-Kemper Reference Webb and Fox-Kemper2011, Reference Webb and Fox-Kemper2015).
Several authors have considered the effect of Stokes drift on the transport of floating marine litter. In an early study, Kubota (Reference Kubota1994) found that Stokes drift derived from local wind fields did not make a significant contribution towards debris transport. However, more recent studies that included the entire wave field showed that Stokes drift could play an important role. For example, Iwasaki et al. (Reference Iwasaki, Isobe, Kako, Uchida and Tokai2017) found that Stokes drift transported plastic towards the coast in the Sea of Japan during winter, and Delandmeter & Van Sebille (Reference Delandmeter and Van Sebille2019) reported similar behaviour in the Norwegian Sea. Stokes drift could enable debris to leak out of the Indian Ocean (Dobler et al. Reference Dobler, Huck, Maes, Grima, Blanke, Martinez and Ardhuin2019), cause drifting debris to cross the strong circumpolar winds and currents to reach the Antarctic coast (Fraser et al. Reference Fraser, Morrison, Hogg, Macaya, van Sebille, Ryan, Padovan, Jack, Valdivia and Waters2018), and thus promote increased transport to polar regions (Onink et al. Reference Onink, Wichmann, Delandmeter and Van Sebille2019). Isobe et al. (Reference Isobe2014) modelled the plastic beaching process by including Stokes drift and sinking velocity and observed that larger plastic debris was selectively moved onshore. All the foregoing studies have simply assumed that floating marine litter objects are transported with the Stokes drift; in other words, that they are perfect Lagrangian tracers.
If a particle is infinitesimally small and has the same density as water, it will behave purely as a Lagrangian tracer and will be transported with the Stokes drift. This is not necessarily true for an object of finite size or of a density different to that of water. As the inertia of such an object becomes important, the fluid will exert a drag on the object owing to the relative velocity between the object and fluid. Furthermore, the object may rise, sink, or float depending on the density difference. The literature distinguishes between fully submerged and floating objects, discussed separately below.
The motion of a fully submerged sphere in unsteady flow with viscous drag can be described by the Maxey–Riley equations (Maxey & Riley Reference Maxey and Riley1983). Based on this pioneering work, Eames (Reference Eames2008) and Santamaria et al. (Reference Santamaria, Boffetta, Martins Afonso, Mazzino, Onorato and Pugliese2013) examined how far slightly positively or negatively buoyant objects would be transported by regular waves. They defined the distance transported as either the horizontal distance transported whilst a negatively buoyant object sinks from the free surface to the sea floor or the horizontal distance transported whilst a positively buoyant object rises from the sea floor to the free surface. Eames (Reference Eames2008) and Santamaria et al. (Reference Santamaria, Boffetta, Martins Afonso, Mazzino, Onorato and Pugliese2013) used an expansion in wave steepness and Stokes number to arrive at analytical solutions for small objects. To leading order and for negatively buoyant objects, Eames (Reference Eames2008) showed such small objects are transported with a mean horizontal Stokes drift velocity and sediment with their terminal fall velocity. Santamaria et al. (Reference Santamaria, Boffetta, Martins Afonso, Mazzino, Onorato and Pugliese2013) predicted that positively buoyant objects would experience an increase in drift owing to their inertia. Although Eames (Reference Eames2008) and Santamaria et al. (Reference Santamaria, Boffetta, Martins Afonso, Mazzino, Onorato and Pugliese2013) considered the object's inertia when examining transport by waves, both considered completely submerged objects.
Also considering fully submerged objects, DiBenedetto & Ouellette (Reference DiBenedetto and Ouellette2018) first showed that non-spherical objects have a preferential orientation under waves, confirming this result numerically (DiBenedetto & Ouellette Reference DiBenedetto and Ouellette2018) and experimentally (DiBenedetto, Koseff & Ouellette Reference DiBenedetto, Koseff and Ouellette2019) but not examining the effect of the object's inertia. The orientation changes the drag on slightly negatively buoyant objects, which results in objects of different shapes being transported different distances before ‘raining out’ (DiBenedetto, Ouellette & Koseff Reference DiBenedetto, Ouellette and Koseff2018).
Analysis of the motion of floating objects commences with the extension of Maxey–Riley equation (Maxey & Riley Reference Maxey and Riley1983) to include a free surface, as undertaken by Rumer, Crissman & Wake (Reference Rumer, Crissman and Wake1979). These authors considered the free surface to be an oscillating slope with a vertical force balance between gravity and buoyancy, whilst the horizontal part of the buoyancy force induces object motion in what Rumer et al. (Reference Rumer, Crissman and Wake1979) termed the slope-sliding effect. Shen & Zhong (Reference Shen and Zhong2001) further extended the slope-sliding model, proceeding to find analytical solutions of the object motion in limit of no added mass or no resistance. Huang, Huang & Law (Reference Huang, Huang and Law2016) found that the drift of relatively large floating discs, used to model floating ice sheets, increased beyond the Stokes drift in physical experiments. This could be explained by numerical solutions to an equation of motion based on a rotating coordinate system which aligned with the free surface, leaving the physical mechanism at work unclear.
Although not focusing on waves, Beron-Vera, Olascoaga & Lumpkin (Reference Beron-Vera, Olascoaga and Lumpkin2016) showed that the inertia of an undrogued drifter is important for their accumulation in subtropic gyres. The study integrated a Maxey–Riley equation that modelled the variable submergence of surface drifters and included forcing from current and wind velocities, by varying the relative effect of each with the submerged volume of the drifter. The drag formulation assumed linear dependence of force on the density ratio between the object and water, as has been experimentally validated by Miron et al. (Reference Miron, Medina, Olascoaga and Beron-Vera2020). The Maxey–Riley equation has been extended to model floating Sargassum rafts (Beron-Vera & Miron Reference Beron-Vera and Miron2020).
Surface tension can be important in the response of small inertial particles under wave action, as shown by Falkovich et al. (Reference Falkovich, Weinberg, Denissenko and Lukaschuk2005), who found that hydrophobic and hydrophilic particles concentrate in antinodes and nodes of a standing wave, respectively. Denissenko, Falkovich & Lukaschuk (Reference Denissenko, Falkovich and Lukaschuk2006) demonstrated the importance of surface tension when predicting time scales of small particle clusters in standing waves. In this paper, we do not examine the effect of surface tension, which places a lower limit on the size of particles for which our model is valid.
This paper examines the transport of inertial, finite-size floating marine litter under the influence of non-breaking waves. Our derivation starts from Newton's second law, with buoyancy, gravity and drag force components. Using a transformed coordinate system, similar but not equivalent to Huang et al. (Reference Huang, Huang and Law2016), that vertically translates and is oriented orthogonally to the time-varying free surface, we ensure that the dynamic buoyancy term is directed normal to the free surface. In this model, the drag force changes with submergence of the object, and we formulate a drag coefficient that is valid across a range of Reynolds numbers. We use perturbation methods to derive a closed-form solution for the transport of inertial, finite-size floating spherical objects, which is then used to interpret the physical mechanism for their enhanced transport compared to the Stokes drift. Numerical and analytical solutions are compared for viscous drag. In order to observe the predicted response, we perform experiments in a laboratory wave flume.
This paper is laid out as follows. Section 2 presents the theoretical model. Section 3 describes solutions obtained using perturbation methods for viscous drag. Section 4 compares the analytical solutions thus obtained against numerical solutions of the model. The numerical solutions are also used to compare model predictions of viscous and non-viscous drag. Conclusions are drawn in § 5.
2. Mathematical model
2.1. Equation of motion of a floating object
The motion of a floating inertial object is described by Newton's second law
\begin{equation} m\dot{\boldsymbol{v}}= \boldsymbol{F}\equiv\boldsymbol{B}+\boldsymbol{M}+\boldsymbol{G}+\boldsymbol{R} , \end{equation}
where 
$m$ is the mass of the object and 
$\boldsymbol {v}$ its velocity with the dot denoting a derivative with respect to time. The total force on the object 
$\boldsymbol {F}$ can be decomposed into a buoyancy force 
$\boldsymbol {B}$, an added-mass force 
$\boldsymbol {M}$, a gravity force 
$\boldsymbol {G}$ and a resistance force 
$\boldsymbol {R}$, which are formulated below. The buoyancy and added-mass forces arise from the integral of pressure around the object. For simplicity, we will assume the object is spherical with diameter 
$D$. Throughout, it is assumed that the object is small relative to the wavelength, such that 
$D/\lambda _0\ll 1$, with 
$D$ the diameter of the object and 
$\lambda _0$ the wavelength. This has four important consequences. First, the wave field is unaffected by the presence of the object; in other words, there is no diffraction. Second, the free surface can be approximated as an (inclined) straight line on the scale of the object. Third, we can approximate the (relative) velocity field between the liquid and object, which determines the drag on the object, as the velocity at a point. Fourth, the buoyancy force can be computed from the submergence measured relative to the free surface. Nevertheless, the model neglects surface tension. This assumption is reasonable for floating objects provided the following threshold criterion (e.g. Falkovich et al. Reference Falkovich, Weinberg, Denissenko and Lukaschuk2005) is met: 
$D/2 > \sqrt {\gamma /(\rho g)}$, where 
$\gamma$ is surface tension, 
$\rho$ is density of water and 
$g$ is gravitational acceleration. For water, the criterion is satisfied for objects of diameter exceeding 5.4 mm, resulting in the findings being invalid for microplastic. However, such small plastics are likely to behave as purely Lagrangian tracers.
We first adopt a stationary two-dimensional laboratory coordinate system (
$x$, 
$z$) with the vertical coordinate 
$z$ measured upwards from the undisturbed free surface. To define the forces on the object, a second, moving coordinate system (
$\tau$, 
$n$) is established that moves vertically with the free surface 
$z=\eta (x,t)$ and aligns locally with the 
$\tau$-axis tangential to the free surface at the position of the object 
$x_p$ and the 
$n$-axis normal to it, as shown in figure 1. The coordinate transformation takes the form of a vertical translation followed by a clockwise rotation through the angle 
$\theta =\arctan (\partial \eta /\partial x)$, both at the horizontal position of the object 
$x_p$
\begin{gather} \begin{bmatrix} \tau \\ n \end{bmatrix}= \begin{bmatrix} 1 & \partial_x \eta(x,t)|_{x_p}\\ -\partial_x \eta(x,t)|_{x_p} & 1 \end{bmatrix} \begin{bmatrix} x \\ z-\eta(x_p,t) \end{bmatrix} \varXi(x_p,t)\nonumber\\ \textrm{with}\ \varXi(x_p,t)\equiv (1+(\partial_x \eta(x,t)|_{x_p})^2)^{{-}1/2},\end{gather}
where a small-angle approximation on 
$\theta$ has been used and 
$\varXi$ is required for the determinant of the transformation matrix to be unity and thus conserve area. The quantities 
$\partial _x\eta (x,t)$, 
$\eta (x,t)$ and 
$\varXi (x,t)$ are evaluated at the object position 
$x_p(t)$ and are thus solely functions of time 
$t$. The coordinate system (
$\tau$, 
$n$) does not translate in the horizontal direction, enabling direct estimation of the object's horizontal drift 
$\bar {v}_x=\bar {\dot {x}}_p$, where the overbar denotes an average over the wave cycle. The time-dependent unit normal vectors are
\begin{equation} \boldsymbol{e}_{\tau} = [1, \partial_x \eta(x,t)|_{x_p}] \varXi(x_p,t) \quad \text{and} \quad \boldsymbol{e}_n=[-\partial_x \eta(x,t)|_{x_p}, 1] \varXi(x_p,t). \end{equation}
It should be emphasised that (
$\tau$, 
$n$) is an accelerating coordinate system, both in terms of rotation and vertical translation. Inverting (2.3)
\begin{align} \boldsymbol{e}_x= (\boldsymbol{e}_{\tau}(t) -\partial_x \eta(x,t)|_{x_p} \boldsymbol{e}_n(t) )\varXi(x_p,t) \quad \text{and} \quad \boldsymbol{e}_z= (\partial_x \eta(x,t)|_{x_p}\boldsymbol{e}_{\tau}(t) + \boldsymbol{e}_n(t))\varXi(x_p,t). \end{align}
For the time-dependent unit normal vectors 
$\boldsymbol {e}_\tau (t)$ and 
$\boldsymbol {e}_n(t)$
\begin{equation} \frac{\text{d} \boldsymbol{e}_{\tau}(t)}{\text{d} t}=\dot{\theta}_p \boldsymbol{e}_{n}(t) \quad \text{and}\quad \frac{\text{d} \boldsymbol{e}_n(t)}{\text{d} t}={-}\dot{\theta}_p \boldsymbol{e}_{\tau} (t) \quad \text{with}\ \dot{\theta}_p=\text{d}_{t}(\partial_{x}\eta(x,t)|_{x_p}) \varXi_p^2, \end{equation}
in which 
$\theta _p(t)\equiv \theta (x_p(t),t)$, 
$\varXi _p(t)=\varXi (x_p(t),t)$ and 
$\textrm {d}_t\equiv \textrm {d}/\textrm {d}t$.
Figure 1. Diagram of the two coordinate systems used to describe a floating object of diameter 
$D$: a stationary laboratory coordinate system (
$x$, 
$z$) and a vertically translating and rotating coordinate system (
$\tau$, 
$n$) with its origin at the vertical position of the free surface 
$z=\eta _p$ and the 
$\tau$-axis aligned tangential to the free surface. The vector 
$\boldsymbol {x}_p$ locates the centre of the object relative to the origin of the stationary coordinate system, 
$\tan \theta =\partial \eta /\partial x$ is the slope of the free surface and 
$s$ is the (variable) submergence.
Denoting the position of the object as 
$\boldsymbol {x}_p=x_p\boldsymbol {e}_x+z_p\boldsymbol {e}_z=\eta _p\boldsymbol {e}_z+\tau _p\boldsymbol {e}_\tau +n_p\boldsymbol {e}_n$ with 
$\eta _p(t)\equiv \eta (x_p(t),t)$, its velocity may be written as
\begin{equation} \boldsymbol{v}=v_x\boldsymbol{e}_x +v_z \boldsymbol{e}_z = (\dot{\tau}_p-\dot{\theta}_p n_p+\dot{\eta}_p\partial_x \eta|_{x_p}\varXi_p)\boldsymbol{e}_{\tau}+(\dot{n}_p+\dot{\theta}_p \tau_p +\dot{\eta}_p\varXi_p)\boldsymbol{e}_n, \end{equation}
where we have used (2.5) for the time derivatives of the unit vectors 
$\boldsymbol {e}_\tau$ and 
$\boldsymbol {e}_n$, and 
$\boldsymbol {e}_z$ was substituted for from (2.4b). The velocity in the translating reference frame 
$\boldsymbol {v}^*$ is related to the velocity in the stationary reference frame 
$\boldsymbol {v}$ by 
$\boldsymbol {v}^*=\boldsymbol {v}-\dot {\eta }_p\boldsymbol {e}_z$, where both vectors can be expressed in any arbitrary set of orthogonal components, such as 
$\boldsymbol {e}_x$ and 
$\boldsymbol {e}_z$ or 
$\boldsymbol {e}_\tau$ and 
$\boldsymbol {e}_n$. Accordingly, the acceleration of the object can be written as
\begin{gather} \dot{\boldsymbol{v}}=\dot{v}_x\boldsymbol{e}_x +\dot{v}_z \boldsymbol{e}_z= (\ddot{\tau}_p-\ddot{\theta}_p n_p-2\dot{\theta}_p \dot{n}_p-(\dot{\theta}_p)^2\tau_p+\ddot{\eta}_p\partial_x \eta|_{x_p}\varXi_p)\boldsymbol{e}_\tau\nonumber\\ +\,(\ddot{n}_p+\ddot{\theta}_p \tau_p+2\dot{\theta}_p \dot{\tau}_p-(\dot{\theta}_p)^2 n_p +\ddot{\eta}_p\varXi_p)\boldsymbol{e}_n. \end{gather}
To evaluate (2.7), the double time derivatives 
$\ddot {\theta }_p$ and 
$\ddot {\eta }_p$ must be evaluated explicitly. The double time derivative 
$\ddot {\theta }_p$ can be obtained by differentiating with respect to time twice using the relationship 
$\theta _p=\arctan (\partial \eta /\partial x |_{x_p})$, noting that 
$x_p$ is a function of time requiring the chain rule, to obtain
\begin{gather} \ddot{\theta}_p=(\partial_{ttx}\eta|_{x_p}+2\dot{x}_p \partial_{txx} \eta|_{x_p} + (\dot{x}_p)^2\partial_{xxx}\eta|_{x_p}+\ddot{x}_p\partial_{xx}\eta|_{x_p})\varXi_p^2 \nonumber\\ +\,(\partial_{tx}\eta|_{x_p} + \dot{x}_p \partial_{xx}\eta|_{x_p}) 2 \varXi_p \dot{\varXi}_p. \end{gather}
Similarly, the double time derivative 
$\ddot {\eta }_p$ takes into account the dependence of the free surface 
$\eta _p(x_p(t),t)$ on time 
$t$ and the time-dependent horizontal position 
$x_p(t)$, which gives through the chain rule after differentiating twice
\begin{equation} \ddot{\eta}_p=\partial_{tt}\eta|_{x_p} + 2\dot{x}_p \partial_{tx}\eta|_{x_p}+ (\dot{x}_p)^2\partial_{xx}\eta|_{x_p} + \ddot{x}_p \partial_x \eta|_{x_p}. \end{equation}
Substituting (2.8) and (2.9) into (2.7) and (2.7) thence into (2.1) results in two second-order differential equations in the (
$n,\tau$) coordinate system, which are explicitly given by (A1) and (A2) in Appendix A. These two equations contain three second-order time derivatives, and so a third (kinematic) equation relating the second-order derivatives is required to solve the system. Such an equation can for example be found by taking the dot product of (2.7) and 
$\boldsymbol {e}_x$ (see (A3) in Appendix A).
For convenience, we express the normal coordinate of the centre of the object 
$n_p$ in terms of the submergence depth 
$s$ (see figure 1). To do so, we assume that 
$D/\lambda _0\ll 1$ so that the free surface is a locally straight line with 
$n$-coordinate 
$n_s=-\partial _x \eta |_{x_p}x_p\varXi _p$ (using (2.2), setting 
$x=x_p$ and 
$z=\eta _p$). The submergence depth is then given by 
$s=D/2-(n_p-n_s)=D/2-n_p-x_p\partial _x\eta |_{x_p}\varXi _p$, where 
$D$ is the diameter of the object. From (2.6), the following expression is obtained for the horizontal velocity of the object:
\begin{equation} \dot{x}_p=(\dot{\tau}_p-\dot{\theta}_p(n_p+\tau_p\partial_x\eta|_{x_p}) -\dot{n}_p\partial_x\eta|_{x_p})\varXi_p. \end{equation}
It should be noted that 
$\dot {s}=-\dot {n}_p-\textrm {d}_t(x_p\partial _x\eta |_{x_p}\varXi _p)$.
2.1.1. Buoyancy and added mass
We decompose total pressure 
$p$ into an undisturbed component 
$p_{undisturbed}$ and a disturbed component 
$p_{disturbed}$ owing to the presence of the object. Assuming an object that is small relative to the wavelength (
$D/\lambda _0\ll 1$), the undisturbed pressure varies as 
$p_{undisturbed}=\rho _f g(\eta (x,t)-z)$ on the scale of the object with 
$\rho _f$ the density of the fluid, so that the dynamic free surface boundary condition 
$p_{undisturbed}(z=\eta )=0$ is satisfied, the variation with depth is hydrostatic, and any depth-dependent variation owing the waves (cf. 
$\exp (k_0 z)$ with 
$k_0$ the wavenumber) is ignored.
The undisturbed pressure integrated around the wetted surface results in a buoyancy force acting in the normal direction to the free surface,
\begin{equation} B_n(t)=\frac{g m }{\beta} \frac{V_{{s}}}{V}\varXi^{{-}1}_p = \frac{g m }{\beta}\left(3\left(\frac{s(t)}{D}\right)^2-2\left(\frac{s(t)}{D}\right)^3\right)\varXi^{{-}1}_p , \end{equation}
where 
$g$ is the gravitational constant, 
$V_s$ is the submerged, 
$V$ the total volume of the sphere and 
$\beta \equiv \rho _o/ \rho _f$ is the ratio of object to fluid density. By including 
$\rho _f g\eta (x,t)$ in the undisturbed pressure, we have included the Froude–Krylov force resulting from the waves.
The disturbed component of pressure leads to added-mass terms, as derived by Maxey & Riley (Reference Maxey and Riley1983)
\begin{equation} M_\tau = \frac{C_{m,\tau}(s) m}{\beta} (\dot{u}_\tau(\boldsymbol{\tilde{x}}_p,t)-\dot{v}_\tau)\quad \text{and}\quad M_n = \frac{C_{m,n}(s)m}{\beta} (\dot{ u}_n(\boldsymbol{\tilde{x}}_p,t)-\dot{v}_n), \end{equation}
where 
$\boldsymbol {C}_m=(C_{m,\tau },C_{m,n})$ is the added-mass coefficient, which is deliberately left as an unspecified function of submergence 
$s(t)$ at this stage of the derivation.
The small-diameter assumption leaves the vertical location, where we should evaluate the velocity of the surrounding fluid in (2.12), unspecified. We set this location to be at the free surface, 
$\tilde {\boldsymbol {x}}_p=(x_p,\eta _p)$.
2.1.2. Gravity forces
The gravity force acts in the vertical direction, and has the following components in the moving coordinate system:
\begin{equation} G_\tau(t)={-}m g\partial_x \eta|_{x_p} \varXi_p(t) \quad \text{and}\quad G_n(t)={-}m g\varXi_p(t). \end{equation}2.1.3. Resistance forces
The resistance terms are caused by drag on the object when it has a velocity relative to that of the surrounding liquid. To begin, we assume viscous drag. We assume this drag depends on the submergence of the object and, specifically, we assume the drag is proportional to the submerged projected area of the sphere in the tangential and normal directions (see figure 2). Other drag formulations are discussed and examined in § 4. The resistance force in the tangential direction is,
\begin{equation} R_{\tau}= 3{\rm \pi} \rho_f \nu D \hat{A}_{{s},\tau}(u^*_{\tau}-v^*_{\tau}), \end{equation}
where 
$u^*_{\tau }$ and 
$v^*_{\tau }$ are the velocity components in the 
$\tau$-direction of the surrounding fluid and the object velocity respectively (in the moving reference frame). The normalised area in the tangential direction 
$\hat {A}_{{s},\tau }$ is the projected area of the submerged sphere,
\begin{equation} A_{{s},\tau}=\frac{D^2}{8}\left(\zeta-\sin(\zeta)\right) \quad \text{with}\ \zeta\equiv 2\cos^{{-}1}\left(1-2s/D\right), \end{equation}
normalised by the maximum projected area 
$A= {\rm \pi} D^2/4$, so that 
$\hat {A}_{{s}, \tau } = A_{{s},\tau }/A$. Assuming the drag is proportional to the submerged projected area following Beron-Vera et al. (Reference Beron-Vera, Olascoaga and Lumpkin2016), which has been validated for steady flows (Miron et al. Reference Miron, Medina, Olascoaga and Beron-Vera2020; Olascoaga et al. Reference Olascoaga, Beron-Vera, Miron, Triñanes, Putman, Lumpkin and Goni2020), we evaluate the fluid velocity 
$u^*_{\tau }$ at the free surface, 
$\tilde {\boldsymbol {x}}_p=(x_p,\eta _p)$.
Figure 2. Diagrams of (a) the submerged volume 
$V_{{s}}$ as a function of the variable submergence 
$s(t)$; (b) the projected area of a submerged sphere moving in the normal direction (
$\boldsymbol {e}_n$); and (c) the projected area of a submerged sphere moving in the tangential direction (
$\boldsymbol {e}_\tau$). All diagrams are shown in the (
$\tau$, 
$n$) coordinate system.
Similar to the 
$\tau$-direction, we have for the 
$n$-direction,
\begin{equation} R_{n}=3{\rm \pi} \rho_f \nu D \hat{A}_{{s},n}(u^*_{n}-v^*_{n}), \end{equation}
where we have evaluated the velocity of the surrounding fluid at the same location 
$\tilde {\boldsymbol {x}}_p$ as for the tangential resistance force. The submerged projected area of a sphere in the normal direction is given by (see figure 2)
\begin{equation} A_{{s},n}={\rm \pi} s(t)\left(D-s(t)\right), \end{equation}
which again, is normalised by the maximum projected area of a sphere 
$A= {\rm \pi} D^2/4$, so that 
$\hat {A}_{{s},n}= A_{{s},n}/A$. Later, in § 4, other drag formulations are considered to examine the robustness of the model's predictions.
2.2. Fluid velocity for surface gravity waves
We consider unidirectional deep-water surface gravity waves propagating over a horizontal bed in the 
$(x,z)$-coordinate system, with 
$z$ measured vertically upwards from still water level, and the free surface located at 
$z=\eta$. For irrotational flow of inviscid, incompressible fluid, the governing (Laplace) equation is,
\begin{equation} \nabla^2 \phi=0, \quad \text{for}\ -d \leq z \leq \eta, \end{equation}
where 
$\phi$ is the velocity potential and 
$d$ depth. Equation (2.18) is solved subject to the no-flow bottom boundary condition,
\begin{equation} \partial_z \phi=0, \quad \text{for} \ z={-}d, \end{equation}and the kinematic and dynamic linear free surface boundary conditions,
\begin{equation} u_z-\partial_t \eta- u \partial_x \eta=0 \quad \text{and} \quad g \eta +\partial_t \phi +\tfrac{1}{2}(\boldsymbol{\nabla} \phi)^2=0 \quad \text{at} \ z=\eta, \end{equation}
where the velocity components are 
$u_x=\partial _x \phi$ and 
$u_z=\partial _z \phi$.
3. Perturbation theory for viscous drag
To interpret the physical mechanism behind the drift predicted by the model derived in § 2, we use perturbation theory to establish an analytical solution. We do so here for the case of viscous drag, as this allows inclusion of drag at first order in our expansion. We will discuss limitations of viscous drag in § 3.4 and consider numerical solutions of our model in § 4 in which the assumption of viscous drag is relaxed. We consider only periodic, weakly nonlinear, deep-water surface gravity waves, so that 
$k_0d\gg 1$ with 
$k_0$ the wavenumber. We perturb the object position 
$\boldsymbol {x}_p$ in a Stokes-type expansion in wave steepness (
$\alpha =k_0 a_0$, where 
$a_0$ the wave amplitude), giving
\begin{equation} \boldsymbol{x}_p(t)=\boldsymbol{x}_p^{(0)} +\alpha \left.\boldsymbol{x}_p^{(1)}(t)\right|_{\boldsymbol{x}_{p}^{(0)}} + \alpha^2 \left.\boldsymbol{x}^{(2)}_p(t)\right|_{\boldsymbol{x}_p^{(0)}}+{O}(\alpha^3), \end{equation}
where the superscript corresponds to the order in 
$\alpha$, and 
$\boldsymbol {x}_p^{(0)}$ is the object label and thus not a function of time. As we are interested in wave-induced drift, which arises at second order, we only pursue those terms necessary to obtain this drift.
Applying a perturbation expansion in the same small parameter 
$\alpha$ to the governing equation of the fluid (2.18) and its boundary conditions (2.19) and (2.20) allows the free surface 
$\eta$ and the velocity potential 
$\phi$ to be determined, and we do so up to second order.
Although the perturbation theory solutions in this section are for regular waves, the experiments introduced in Appendix B make use of long (or narrow-bandwidth) wave packets for practical reasons. We assume that inertial effects do not arise on the scale of the packets, as justified in Appendix C, so that we can correct for the presence of a wave packet simply by accounting for its Eulerian mean flow. Table 1 lists the resulting solutions, whose derivation and laboratory validation is given in more detail by van den Bremer et al. (Reference van den Bremer, Whittaker, Calvert, Raby and Taylor2019) for deep water and Calvert et al. (Reference Calvert, Whittaker, Raby, Taylor, Borthwick and van den Bremer2019) for intermediate depth. We consider only deep-water waves here (
$k_0d\gg 1$). The solutions for the Eulerian return flow and the second-order surface elevation are based on wave packets with envelope 
$|A_0|$. It is assumed that the wave packets are narrow banded and that the Eulerian return flow is shallow, corresponding to a depth that is small relative to the packet length (Calvert et al. (Reference Calvert, Whittaker, Raby, Taylor, Borthwick and van den Bremer2019) establish the Eulerian return flow without the shallow return flow assumption). In practice, inclusion of the effect of the return flow merely leads to a small correction of less than 
$2\,\%$ for our laboratory experiments.
Table 1. First- and second-order solutions for the kinematic properties of deep-water surface gravity waves, with 
$A_0= a_0\hat {A}_0$ the wave amplitude envelope, 
$a_0$ its amplitude, 
$\hat {A}_0$ a non-dimensional envelope, 
$\omega _0$ the carrier wave frequency and 
$k_0$ the carrier wavenumber. Where complex fields are given, the real part is understood, and 
$\varphi = k_0 x- \omega _0 t$. The first three rows are first-order solutions, valid for regular waves or wave packets. The remaining rows comprise second-order solutions for the wave-averaged Eulerian and Stokes velocities and the set-down. The second-order wave-averaged Eulerian velocity only arises for wave packets, considered in the experiments in Appendix B.
3.1. Zeroth order in wave steepness: 
${O}(\alpha ^0)$
At zeroth order in wave steepness, wave forcing evidently does not play a role. Only the normal direction of (2.1) has any forcing at zeroth order, where the following leading-order static balance is achieved between buoyancy force and gravity,
\begin{equation} F^{(0)}_n= \frac{g m}{\beta}\left[ 3 \left(\frac{s^{(0)}}{D} \right)^2- 2\left(\frac{s^{(0)}}{D} \right)^3\right] - gm =0. \end{equation}
We have used the fact that 
$\varXi _p=1$ at zeroth order and note that (3.2) is only valid for a floating sphere, i.e. 
$|D/2-s^{(0)}| \leq D/2$. Equation (3.2) is a cubic equation, which can be readily solved numerically for the depth of submergence of a floating sphere in the absence of waves 
$s^{(0)}$.
3.2. First order in wave steepness: ${O}(\alpha ^1)$
We begin by expressing the projected areas of the sphere required to calculate the tangential and normal resistance forces as series expansions around 
$s^{(0)}$. The submerged projected area of a sphere in the tangential direction (2.15) can be approximated by
\begin{equation} A_{{s},\tau}(s)=A_{{s},\tau}(s^{(0)}) +2D\sqrt{\frac{s^{(0)}}{D} -\left(\frac{s^{(0)}}{D}\right)^2}s^{(1)}+{O}(\alpha^2), \end{equation}
where we have obtained 
$\partial _s(A_{{s},\tau })$ from (2.15) by implicit differentiation. For the submerged projected area of a sphere in the normal direction, it is sufficient for our purposes to evaluate 
$A_{{s},n}(s)$ at zeroth order, i.e. 
$A_{{s},n}(s)=A_{{s},n}(s^{(0)})+{O}(\alpha ^1)$.
3.2.1. The tangential direction
To first order of approximation, the velocity and acceleration in the horizontal coordinate 
$x$ and the tangential coordinate 
$\tau$ are equal, i.e. 
$\dot {x}_p^{(1)}=v_x^{(1)}=v_{\tau }^{(1)}$ and 
$\ddot {x}_p^{(1)}=\dot {v}_x^{(1)}=\dot {v}_{\tau }^{(1)}$. The only forces that play a role are the tangential components of the added mass, gravity and the resistance force. The first-order added-mass terms in the tangential direction are
\begin{equation} M_\tau^{(1)}=\frac{C_m m}{\beta}( \dot{u}_x^{(1)} -\ddot{x}_p^{(1)}), \end{equation}
where we now assume for simplicity that the added-mass coefficient 
$C_m$ is a constant and independent of direction. Other added-mass formulations are discussed and examined in § 4.
In a potential flow, a fully submerged sphere has an added-mass coefficient of 
$1/2$. Instead of deriving the complicated dependence of 
$C_m$ on the object's density, we interpolate linearly between the values for a sphere that is fully submerged (
$\beta =1$, 
$C_m=1/2$) and a sphere that is entirely out of the water (
$\beta =0$, 
$C_m=0$) and set 
$C_m= \beta /2$. The robustness of this assumption is investigated numerically in § 4.
The resistance force (2.14) can be approximated as
\begin{equation} R_\tau^{(1)}=\varGamma_R m\omega_0 \hat{A}_{{s},\tau}^{(0)}(u_x^{(1)}|_{\tilde{\boldsymbol{x}}_{p}^{(0)}}-\dot{x}_p^{(1)}) \quad \text{with}\ \varGamma_R\equiv \frac{3 {\rm \pi} \nu D}{\beta V \omega_0},\end{equation}
where the non-dimensional coefficient 
$\varGamma _R$ measures the importance of the resistance force.
From the object's equation of motion (2.1) we thus obtain
\begin{equation} \left(1+\frac{C_m}{\beta}\right) \ddot{x}_p^{(1)}=\frac{C_m}{\beta} \dot{u}_x^{(1)}|_{\tilde{x}_p^{(0)}} - g \partial_x \eta^{(1)}|_{x_{p}^{(0)}}+\varGamma_R \hat{A}^{(0)}_{{s},\tau}\omega_0 (u_x^{(1)}|_{\tilde{\boldsymbol{x}}_{p}^{(0)}}-\dot{x}_p^{(1)}). \end{equation}
We seek a solution to the forced second-order ordinary differential equation (3.6) of the form 
$x_p^{(1)}=\mathcal {R}(i {X}^{(1)}a_0 \exp (\textrm {i} \varphi _{p}^{(0)}))$ with 
$\varphi _{p}^{(0)}=k_0x^{(0)}_p-\omega _0 t+\varphi _0$ and 
$\varphi _0=\arg (A_0)$, ignoring initial transients. The complex coefficient 
$X^{(1)}$ represents the amplitude and phase change of the horizontal motion of the object relative to those of an idealised Lagrangian object under the influence of waves at the same order, 
$x_{L}^{(1)}=\mathcal {R}(i a_0 \exp (\textrm {i} \varphi _{p}^{(0)}))$. We obtain 
$X^{(1)}=1$, i.e. there is no horizontal motion amplification compared to that of a Lagrangian particle.
3.2.2. The normal direction
Expressing the submergence depth 
$s$ in terms of the vertical coordinate 
$z_p$, we have, without approximation, that 
$s=D/2-(z_p-\eta _p)\varXi _p$. Therefore, the velocity and acceleration in the vertical coordinate 
$z$ and the normal coordinate 
$n$ are related to first order by
\begin{equation} \dot{z}_p^{(1)}=v_z^{(1)}={-}\dot{s}^{(1)}+\dot{\eta}_p^{(1)} \quad \text{and}\quad \ddot{z}_p^{(1)}=\dot{v}_z^{(1)}={-}\ddot{s}^{(1)}+\ddot{\eta}_p^{(1)}. \end{equation}We first approximate the buoyancy force (2.11) by
\begin{equation} B_n^{(1)}= \varGamma_B m\omega_0^2 s^{(1)} \quad \text{with} \ \varGamma_B\equiv \frac{6}{\beta k_0 D}\left(\frac{s^{(0)}}{D}-\left(\frac{s^{(0)}}{D}\right)^2\right), \end{equation}the added-mass terms by
\begin{equation} M_n^{(1)}=\frac{C_m m}{\beta} \ddot{s}^{(1)}, \end{equation}and the resistance force (2.16) by
\begin{equation} R_n^{(1)}=\varGamma_R m\omega_0 \hat{A}_{{s},n}^{(0)}\dot{s}^{(1)}, \end{equation}
where we have used 
$u_z^{(1)}(z=0)=\dot {\eta }_p^{(1)}$ from the linearised kinematic free surface boundary condition and 
$v_n^{(1)}=\dot {z}_p^{(1)}$. The new non-dimensional coefficient 
$\varGamma _B$ measures the strength of dynamic buoyancy, and 
$\varGamma _R$ measures the strength of the resistance force, as for the tangential resistance force in (3.5). From the object's equation of motion (2.1) we thus obtain
\begin{equation} \left(1+\frac{C_m}{\beta}\right)(\ddot{\eta}^{(1)}_p-\ddot{s}^{(1)}) =\frac{C_m}{\beta}\dot{u}_z^{(1)}|_{\tilde{x}_p^{(0)}} +\varGamma_B \omega_0^2s^{(1)} +\varGamma_R \hat{A}^{(0)}_{{s},n}\omega_0 \dot{s}^{(1)}, \end{equation}
where we note gravity only enters at zeroth order. As for the tangential direction, we seek a solution to the forced second-order ordinary differential equation (3.11) of the form 
$s^{(1)}=\mathcal {R}(\mathcal {S}^{(1)}a_0 \exp (\textrm {i} \varphi _{p}^{(0)}))$ with 
$\varphi _{p}^{(0)}=k_0x^{(0)}_p-\omega _0 t+\varphi _0$ and 
$\varphi _0=\arg (A_0)$, ignoring initial transients. We find for the non-dimensional submergence at first order 
$\mathcal {S}^{(1)}$
\begin{equation} \mathcal{S}^{(1)}=\frac{1+\displaystyle \frac{C_m}{\beta}-\varGamma_B-\textrm{i}\varGamma_R \hat{A}^{(0)}_{{s},n}}{\left(1+\displaystyle \frac{C_m}{\beta}-\varGamma_B\right)^2+(\varGamma_R \hat{A}^{(0)}_{{s},n})^2}. \end{equation} Figures 3 and 4 respectively show the magnitudes and arguments of the first-order solutions for the horizontal motion amplification 
$X^{(1)}$ and the variable submergence 
$\mathcal {S}^{(1)}$. In these figures, the purely Lagrangian limit, in which the object is simply transported with the Stokes drift and floats on the moving surface, corresponds to 
$X^{(1)}=1$, 
$\mathcal {S}^{(1)}=0$. This limit is obtained as the object size tends to zero. Note that the phase of variable submergence in this limit is non-zero, 
$\text {arg}(\mathcal {S}^{(1)})\rightarrow {\rm \pi} /2$. This is because both imaginary and real parts of the variable submergence tend to zero, with the imaginary part approaching zero at a faster rate. As our model is only valid for objects that are small relative to the wavelength, we truncate the 
$x$-axis at 
$D/\lambda _0=6\,\%$. Diffraction of the wave field typically only becomes important for 
$D/\lambda _0>20\,\%$.
Figure 3. For viscous drag, magnitudes of the first-order horizontal motion amplification 
$X^{(1)}$ (a) and the variable submergence 
$\mathcal {S}^{(1)}$ (b) as functions of dimensionless object size 
$D/\lambda _0$ for different density ratios 
$\beta =\rho _o/\rho _f$, where the density ratio for each colour is shown in the legend. We have set 
$C_m= \beta /2$. Numerical and analytical solutions from perturbation theory are denoted by crosses and solid lines, respectively.
Figure 4. For viscous drag, arguments of the first-order horizontal motion amplification 
$X^{(1)}$ (a) and the variable submergence 
$\mathcal {S}^{(1)}$ (b) as functions of dimensionless object size 
$D/\lambda _0$ for viscous drag and for different density ratios 
$\beta =\rho _o/\rho _f$, as shown in the legend. We have set 
$C_m= \beta /2$. Numerical and analytical solutions from perturbation theory are denoted by crosses and solid lines, respectively.
As confirmed in figure 3(a), the magnitude of the horizontal motion 
$|X^{(1)}|$ is equivalent to that of a purely Lagrangian tracer. Turning to figure 4(a), the argument of the horizontal motion 
$\arg (X^{(1)})$ is evidently also zero. As shown in figure 3(b), the magnitude of the variable submergence 
$|\mathcal {S}^{(1)}|$ increases monotonically with object size and does so at a larger rate for density ratios closer to unity. Variable submergence is driven by the free surface elevation and governed by drag, dynamic buoyancy and (added) mass, which are respectively the resistance, spring and inertia terms of a forced spring–mass–damper system (cf. (3.11)). The larger the object, the more dominant is the acceleration of the free surface, which acts as an apparent force in the moving reference frame in which the variable submergence is defined, thus increasing the ‘bobbing’ of the object. The lower the density ratio, the stronger the buoyancy force and the stiffer the ‘spring’. The response in variable submergence for a stiffer ‘spring’ is smaller. The argument of variable submergence 
$\arg (\mathcal {S}^{(1)})$ decreases monotonically with object size and growing importance of inertia but is dependent on the density ratio, as shown in figure 4(b).
At first order in steepness the tangential and normal directions are independent, and so it is possible for there to be a significant change in first-order variable submergence whilst the first-order horizontal motion remains unchanged. As can be seen in the next section, a change in first-order variable submergence results in a change in horizontal motion at second order.
3.3. Second order in wave steepness: ${O}(\alpha ^2)$
The equation of motion (2.1) resolved in the horizontal direction and at second order of approximation gives
\begin{equation} \ddot{x}^{(2)}_p= \frac{1}{m}(F^{(2)}_{\tau} - \partial_x \eta^{(1)}|_{x_p^{(0)}} F^{(1)}_n). \end{equation}
In order to examine the wave-induced drift of a floating object in periodic waves, we consider the steady wave-averaged transport and set 
$\bar {\ddot {x}}^{(2)}_p=0$, so that the resultant force must be zero. We will now consider the tangential and normal force contributions to (3.13) in turn.
3.3.1. Tangential and normal directions
In the tangential direction, the added-mass terms at second order can be obtained from the combination of an expansion in the horizontal and vertical displacements of the object, a coordinate transformation and evaluation of the advective derivative, respectively
\begin{align} M_\tau^{(2)} &= \frac{C_m m}{\beta} (\dot{u}_x^{(2)} +x_p^{(1)}\partial_x \dot{u}_x^{(1)}|_{\tilde{x}_p^{(0)}} +\eta_p^{(1)} \partial_z \dot{u}_x^{(1)}|_{\tilde{x}_p^{(0)}}\nonumber\\ &\quad +\dot{u}_z^{(1)}|_{\tilde{x}_p^{(0)}} \partial_x \eta^{(1)}|_{x_p^{(0)}} +\dot{x}_p^{(1)} \partial_x u_x^{(1)}|_{\tilde{x}_p^{(0)}} +\dot{\eta}_p^{(1)} \partial_z u_x^{(1)}|_{\tilde{x}_p^{(0)}} -\dot{v}_{\tau}^{(2)}). \end{align}In addition to the added-mass terms, the tangential force consists of a correction to the tangential component of gravity due to the object's horizontal displacement,
\begin{equation} G_{\tau}^{(2)}=\left.-mg\partial_{xx}\eta^{(1)}\right|_{x_p^{(0)}}x_{p}^{(1)}, \end{equation}and a tangential resistance force,
\begin{equation} R_{\tau}^{(2)}=3{\rm \pi}\rho_f\nu D(\hat{A}_{{s},\tau}^{(1)}(u_{\tau,p}^{(1)}-v_{\tau}^{(1)}) +\hat{A}_{{s},\tau}^{(0)}(u_{\tau,p}^{(2)}-v_{\tau}^{(2)})). \end{equation}
For the first-order velocity components, we have 
$u_{\tau ,p}^{(1)}=u_x^{(1)}|_{\tilde {\boldsymbol {x}}_p^{(0)}}$ and 
$v_{\tau }^{(1)}=\dot {x}_{p}^{(1)}$. Noting from the coordinate transformation that 
$u_\tau =u_x+\partial _x\eta |_{x_p}u_z+{O}(\alpha ^3)$, we obtain for the second-order accurate horizontal fluid velocity component at the object position
\begin{equation} u_{\tau,p}^{(2)}= u_{x}^{(2)}|_{\tilde{\boldsymbol{x}}_p^{(0)}} +\partial_x u_{x}^{(1)}|_{\tilde{\boldsymbol{x}}_p^{(0)}}x_p^{(1)} +\partial_z u_{x}^{(1)}|_{\tilde{\boldsymbol{x}}_p^{(0)}}\tilde{z}_p^{(1)} +\partial_x\eta^{(1)}|_{x_p^{(0)}}u_z^{(1)}|_{\tilde{\boldsymbol{x}}_p^{(0)}}. \end{equation}
We set the second-order Eulerian wave-induced velocity 
$u_{x}^{(2)}$ to zero for the regular waves considered here. The object's horizontal velocity component at second order is
\begin{equation} v_{\tau}^{(2)}=\dot{x}_{p}^{(2)}+\partial_x \eta^{(1)}|_{x_p^{(0)}}\dot{z}_{p}^{(1)}, \end{equation}
where 
$\dot {x}_{p}^{(2)}$ is the quantity that is ultimately of interest. Combining (3.17) and (3.18) and substituting into (3.16) gives
\begin{align} R_{\tau}^{(2)}&=3{\rm \pi}\rho_f\nu D( \hat{A}_{{s},\tau}^{(1)}(u_x^{(1)}|_{\tilde{\boldsymbol{x}}_p^{(0)}}-\dot{x}_{p}^{(1)})\nonumber\\ &\quad +\hat{A}_{{s},\tau}^{(0)}( \partial_x u_{x}^{(1)}|_{\tilde{\boldsymbol{x}}_p^{(0)}}x_p^{(1)} +\partial_z u_{x}^{(1)}|_{\tilde{\boldsymbol{x}}_p^{(0)}}\eta_p^{(1)} -\dot{x}_{p}^{(2)}+\partial_x \eta^{(1)}|_{x_p^{(0)}}\dot{s}^{(1)})),\end{align}
where we have substituted 
$u_x^{(2)}=0$, 
$\dot {z}_{p}^{(1)}=\dot {\eta }_p^{(1)}-\dot {s}^{(1)}$ and 
$u_z^{(1)}|_{\tilde {\boldsymbol {x}}_p^{(0)}}=\dot {\eta }_p^{(1)}$ from the linearised kinematic free surface boundary condition. We use the notation 
$\hat {A}^{(1)}_{{s},\tau }=\hat {A}^{\prime (0)}_{{s},\tau }(s^{(1)}/D)$ with 
$\hat {A}^{\prime (0)}_{{s},\tau }\equiv \partial _{\hat {s}} \hat {A}_{{s},\tau }(\hat {s})|_{\hat {s}^{(0)}}$ and 
$\hat {s}\equiv s/D$ according to (3.3).
In the normal direction, the total force at first order consists of a buoyancy force, an added mass and a resistance force already evaluated in (3.8), (3.9) and (3.10), respectively.
3.3.2. The wave-induced drift
Substituting the first-order solutions for 
$x_p^{(1)}$ (i.e. 
$X^{(1)}=1$) and for 
$s^{(1)}$ from (3.12) and for the wave quantities from table 1 and averaging over the waves, we obtain the following expression from (3.13) for the wave-induced drift of the object 
$\bar {v}_x={\bar {\dot {x}_p^{(2)}}}$:
\begin{align} \bar{v}_x=\frac{u_{S}}{2}\left[\overbrace{2-\underbrace{\mathcal{R}(\mathcal{S}^{(1)})}_{\substack{\text{Increases}\\\text{drift}}}}^{\text{Adjusted Stokes drift}} +\frac{1}{\hat{A}^{(0)}_{{s},\tau}\varGamma_{R}}\left(\overbrace{\underbrace{-\varGamma_B \mathcal{I}(\mathcal{S}^{(1)})}_{\text{Increases drift}}}^{\substack{\text{Buoyancy} \\ \text{resolved into} \\ \text{the } x\text{-direction}}} + \overbrace{\underbrace{\frac{C_m \mathcal{I}(\mathcal{S}^{(1)})}{\beta }}_{\text{Negligible effect}}}^{\text{Added mass}}\right) +\overbrace{\underbrace{\frac{\hat{A}^{(0)}_{{s},n}}{\hat{A}^{(0)}_{{s},\tau}}\mathcal{R}(\mathcal{S}^{(1)})}_{\text{Reduces drift}}}^{\text{Normal drag}}\right], \end{align}
where 
$u_S=k_0\omega _0a_0^2$ is the Stokes drift. We define the drift amplification factor 
$X^{(2)}\equiv {\bar {v_x}}/u_S$, so that 
$X^{(2)}$ corresponds to the terms inside the square brackets in (3.20) divided by 
$2$. Equation (3.20) is the main result of this paper, and we will interpret it below. The text above the terms explains their physical origins, and the text below their effect on the wave-induced drift of the object compared to the Stokes drift.
We begin by examining the wave-induced drift amplification factor 
$X^{(2)}$ as a function of object size and for different density ratios in figure 5. It is evident that the drift is enhanced and increasingly so for larger and heavier objects. Figure 6 examines the contributions to 
$X^{(2)}$ of the four components in (3.20): the adjusted Stokes drift, buoyancy resolved in the 
$x$-direction, normal drag and added mass, which we will discuss in turn. In (3.20) and figure 6, 
$X^{(2)}=1$ corresponds to objects that do not experience an increase in drift and are simply transported with the Stokes drift (i.e. 
$\bar {v}_x=u_{S}$).
Figure 5. For viscous drag, wave-induced drift amplification 
$X^{(2)}$ as a function of dimensionless object size 
$D/\lambda _0$ for different density ratios 
$\beta =\rho _o/\rho _f$ (see legend). We have set 
$C_m= \beta /2$. Numerical and analytical solutions from perturbation theory are denoted by crosses and solid lines, respectively.
Figure 6. For viscous drag, contributions to the wave-induced drift amplification 
$X^{(2)}$ from the five components in (3.20) as a function of non-dimensional object size 
$D/\lambda _0$ for density ratio 
$\beta =0.8$ and 
$C_m= \beta /2$.
3.3.3. Adjusted Stokes drift
The adjusted Stokes drift terms in (3.20) reflect change in linear object trajectory. For unmodified horizontal motion (
$X^{(1)}=1$) and zero variable submergence (
$\mathcal {S}^{(1)}=0$), we obtain 
$X^{(2)}=1$ from the adjusted Stokes drift terms alone. For larger objects, the increase in the vertical motion due to ‘bobbing’ of the object effectively enhances the Stokes drift, as shown in figure 6. This mechanism occurs because the linear variable submergence changes the object's orbit and hence its velocity and time spent under trough and crest. Integration of the linear velocity component along the linear orbit results in Stokes drift. Hence, changes to velocity and orbit result in an adjusted Stokes drift.
3.3.4. Buoyancy resolved in the $x$-direction
The mechanism through which buoyancy, when resolved in the 
$x$-direction and averaged over the wave cycle, can increase the drift of an object is illustrated in figure 7. Without variable submergence (left column), the dynamic buoyancy force is simply zero. With variable submergence but without drag in the normal direction (middle column), the first-order buoyancy force resolved in the 
$x$-direction does not result in a net force on the object, as the first-order buoyancy force and the first-order slope required to resolve this force into the 
$x$-direction are out of phase. It is only in the presence of a drag component in the normal direction (right column) that a phase lag in the submergence depth arises and a net force results. As shown in figure 6, the buoyancy force thus makes a relatively large contribution to the object's drift.
Figure 7. Schematics of the object trajectory (red) and free surface (blue) for three cases: no variable submergence, variable submergence with no normal drag and variable submergence with normal drag. The schematics illustrate the physical mechanism for increased drift arising from variable submergence 
$s^{(1)}$, where variable submergence and drag are in the 
$n$-direction, and a mean motion in the 
$x$-direction is created due to the slope of the free surface 
$\partial _x\eta ^{(1)}$. For this illustration, we have chosen a density ratio 
$\beta =1/2$.
3.3.5. Normal drag
Although normal drag is required to create the phase difference that leads to the net buoyancy force resolved in the 
$x$-direction, normal drag also acts to reduce the magnitude of the ‘bobbing’ mechanism and thus reduces the drift motion, as shown in figure 3. The horizontal direction component of normal drag opposes the horizontal direction component of buoyancy force, with the balance resulting in a drift that is greater than the adjusted Stokes drift discussed above. Tangential drag, through the inverse dependence of 
$X^{(2)}$ on the projected area 
$\hat {A}_{s,\tau }^{(0)}$ and the effective drag coefficient 
$\varGamma _R$ in (3.20), acts to reduce the increase in object drift, by effectively ‘anchoring’ the object to the fluid and its Stokes drift.
3.3.6. Added mass
At first order, the object accelerates in the normal direction, experiencing an inertia force in addition to the buoyancy force and the normal drag discussed above, and so an added-mass term has to be take into account. As shown in figure 3, the contribution by added mass is relatively small and acts to reduce drift.
3.4. Limitation on validity of viscous drag
Although the preceding analysis has demonstrated how enhanced drift of non-infinitesimal objects may arise, the underlying assumption of viscous drag places an upper limit on object size. The maximum Reynolds number that arises from the linear motion in the normal direction is estimated from
\begin{equation} Re_{max} = \frac{a_0 \omega_0\ | \mathcal{S}^{(1)}| D }{\nu} \leq 2, \end{equation}
where we take 
$2$ to be the maximum Reynolds number for drag to be considered viscous. Noting that 
$\mathcal {S}^{(1)}(D/\lambda _0,\beta )$ and taking 
$\beta =0.8$, we obtain from (3.21) for the maximum diameter that
\begin{equation} \mathcal{S}^{(1)}\left(D_{max}/\lambda_0,\beta=0.8\right)(D/\lambda_0) =\frac{k_0^2\nu }{\alpha\omega_0{\rm \pi}}. \end{equation}
For a typical laboratory water wave of steepness 
$\alpha =0.1$ and frequency 
$f_0=1.25\ \textrm {Hz}$, the right-hand side of (3.22) becomes equal to 
$1.6\times 10^{-5}$. Fitting a linear curve 
$S^{(1)}=5.8 D/\lambda _0$ to figure 3(b), we can solve the quadratic (3.22) in 
$D/\lambda _0$ and obtain a maximum diameter to wavelength ratio of 
$0.2\,\%$ corresponding to 
$Re_{max}=2$. Examining figure 5, we can conclude that drift enhancement is negligible for such small objects. We will therefore have to use a realistic, non-viscous drag formulation, as discussed in the next section.
4. Numerical solutions
To validate the perturbation theory for viscous drag in § 3 and to explore the predictions of our model for realistic, non-viscous drag, we set out to obtain numerical solutions of our model. Specifically, we solved the set of differential equations ((A1)–(A3)) with the forces described in detail in § 2 using a numerical ordinary differential equation solver. The fluid velocity and free surface elevation from table 1 were used as input. We first consider viscous drag in § 4.1 and then non-viscous drag in § 4.2, distinguishing conditions (notably Reynolds numbers) that are representative of laboratory (§ 4.2.1; see Appendix B for further details) and field scale (§ 4.2.2). Appendix D discusses the small-object limit of the numerical solutions. Alternative drag and added-mass formulations are examined in Appendix E.
The numerical solutions commenced from an initial condition in the absence of waves with the object depth set at the static submergence given by numerical solution of (3.2). Numerical integration in time was carried out using an explicit Runge–Kutta method with variable time step based on Dormand & Prince's (Reference Dormand and Prince1980) formulation, which is fifth order in time and fourth order in accuracy. Avoiding initial transients, wave forcing was ramped up using half of a Gaussian envelope to steady state. A convergence study showed that a Gaussian half-width set to 20 wavelengths was sufficient to avoid initial transients, whilst the spatial and temporal convergences were in part resolved by the variable time step method and checked explicitly for the largest objects. Once the object motion reached steady state, its motion components in the 
$x$ and 
$z$ directions were effectively linearised using a band-pass filter between 
$0.8 f_0$ and 
$1.2 f_0$. The linear phase was determined using the cross-correlation of the linearised object motion and the linearised Eulerian velocity evaluated at the object position in both directions. The object drift velocity was calculated as the gradient of a straight line fitted to the sub-harmonic 
$x(t)$ motion obtained by low-pass filtering at 
$0.5f_0$.
4.1. Viscous drag
The crosses in figures 3–5 display the numerical solutions of the model with a viscous drag formulation for a (small) steepness 
$\alpha =0.02$. Near perfect agreement is evident with the perturbation theory solutions shown as continuous lines for both the first-order amplitudes (figure 3) and phases (figures 4) and the second-order drift (figure 5). Tiny discrepancies between perturbation theory and numerical simulations in these figures are due to the inherent inclusion of higher-order terms (beyond second order) in steepness in the numerical simulations. The comparison verifies both the numerical model and the second-order perturbation theory.
4.2. Non-viscous drag
To overcome the maximum-Reynolds-number limit of the viscous drag formulation (of 
$Re\equiv |\boldsymbol {u}-\boldsymbol {v}|D / \nu = 2.5 \times 10^4$), we also consider the following non-viscous drag formulation:
\begin{equation} R_{j}(t)=\tfrac{1}{2}C_{d}\left(Re\right) \rho_f A_{{s},j} |u^*_{j}(\tilde{\boldsymbol{x}}_p,t)-v^*_{j}(t)|(u^*_{j} (\tilde{\boldsymbol{x}}_p,t)-v^*_{j}(t)),\end{equation}
where the indices 
$j =n,\tau$ represent the tangential and normal directions; and drag is determined using an experimentally fitted, non-viscous drag coefficient 
$C_{d}$. We choose a formulation of the drag coefficient 
$C_d(Re)$ that captures both viscous drag at small Reynolds number, which is linear in velocity difference, and form drag at high Reynolds number. Specifically, we use the fit to experimental data for drag on a sphere obtained by Morrison (Reference Morrison2013, p. 625), which is accurate for 
$Re<1\times 10^6$
\begin{align} C_d(Re)&=\frac{24}{Re}+ 2.6\frac{Re/5}{(1+Re/5)^{1.52}}+0.411 \frac{(Re/(2.63\times 10^5))^{{-}7.94}}{(1+Re/(2.63\times 10^5))^{{-}8}} \nonumber\\ &\quad +0.25\frac{Re/(1\times 10^6)}{1+Re/(1\times 10^6)}, \end{align}
where (4.2) is the same in both directions because the Reynolds number is independent of direction (
$Re \equiv |\boldsymbol {u}-\boldsymbol {v}|D / \nu$). Taking the small-object and thus the small-Reynolds-number limit of the drag force in (4.1) we can recover the viscous drag on a partially submerged sphere (2.14) and (2.16).
4.2.1. Laboratory-scale results
At laboratory scale, we set 
$f_0=1.25\ \textrm {Hz}$, corresponding to 
$\lambda _0=1.0\ \textrm {m}$ and 
$\alpha =0.1$. With object diameters up to 
$D=60\ \textrm {mm}$, we obtain 
$D/\lambda _0=6\,\%$, where the limit of validity for viscous drag is 
$D/\lambda _0=0.2\,\%$ (see § 3.4). At laboratory scale, figure 8 compares the analytically predicted linear motion using viscous drag with the corresponding numerical results using non-viscous drag. The response in the normal direction is unchanged because the forcing is inertial with little effect from drag. As the object size increases, inertia increasingly dominates over drag. A small decrease in horizontal linear motion is evident reaching a few per cent for larger objects. The results for small objects are the same because the non-viscous drag recovers viscous drag in the small-object limit.
Figure 8. Laboratory-scale numerical simulation results using non-viscous drag for magnitudes of the first-order horizontal motion amplification 
$X^{(1)}$ (a) and the variable submergence 
$\mathcal {S}^{(1)}$ (b) as functions of dimensionless object size 
$D/\lambda _0$ for different density ratios 
$\beta =\rho _o/\rho _f$, where the density ratio corresponding to each colour is listed in the legend. Here, 
$C_m= \beta /2$. Numerical and analytical solutions from perturbation theory are denoted by crosses and solid lines, respectively.
The drift amplification increases slightly when using non-viscous drag for larger objects, as seen in figure 9. This is because the (tangential) drag force for larger objects is lower for non-viscous drag than for viscous drag, resulting in reduced resistance to increased drift compared to the Stokes drift. The maximum Reynolds number reached in the numerical solutions at laboratory scale was 
$Re_{{max}}=3.1 \times 10^4$.
Figure 9. Laboratory-scale numerical simulation results using non-viscous drag for wave-induced drift amplification 
$X^{(2)}$ as a function of dimensionless object size 
$D/\lambda _0$ and for different density ratios 
$\beta =\rho _o/\rho _f$ (see legend). Here, 
$C_m=\beta /2$. Analytical solutions using viscous drag from perturbation theory are denoted by solid lines.
4.2.2. Field-scale results
We set a wave frequency of 
$f_0=0.2\ \textrm {Hz}$ and a steepness of 
$\alpha =0.05$ to represent a typical wind wave at field scale. The frequency of 
$0.2\ \textrm {Hz}$ corresponds to the peak in the spectrum with 
$\alpha =0.05$ at the upper end of the steepness range for wind waves in the ocean (Toffoli & Bitner-Gregersen Reference Toffoli and Bitner-Gregersen2017). This steepness corresponds to a dimensional wave amplitude of 
$a_0=0.3\ \textrm {m}$. The difference between viscous and non-viscous drag results will be larger at field scale owing to the higher value of Reynolds numbers, which reached a maximum of 
$Re_{{max}}= 7.3\times 10^5$ in the numerical simulations.
Figure 10(a) shows the linear horizontal motion, which is mostly unchanged from the perturbation theory result. The magnitude of variable submergence is inertia driven and thus very similar to the viscous analytical result shown in figure 10(b).
Figure 10. Field-scale numerical simulation results using non-viscous drag for magnitudes of the first-order horizontal motion amplification 
$X^{(1)}$ (a) and variable submergence 
$\mathcal {S}^{(1)}$ (b) as functions of dimensionless object size 
$D/\lambda _0$ for different density ratios 
$\beta =\rho _o/\rho _f$, where the density ratio corresponding to each colour is shown in the legend. Field scale here denotes a 
$0.2\ \textrm {Hz}$ wave with a steepness of 
$\alpha =0.05$. Here, 
$C_m= \beta /2$. Numerical and analytical solutions from perturbation theory are denoted by crosses and solid lines, respectively.
The drift amplification for field-scale simulations using non-viscous drag shown in figure 11 is greater than the perturbation theory result based on viscous drag, and even more so than at laboratory scale. This is because the non-viscous drag force is now considerably smaller than its viscous equivalent (taken outside the range of Reynolds numbers for which it is valid). The (tangential) drag force obtained for larger objects is lower for non-viscous drag than for a viscous drag formulation, resulting in reduced resistance to increased drift compared to the Stokes drift.
Figure 11. Field-scale numerical simulation results using non-viscous drag for the wave-induced drift amplification 
$X^{(2)}$ as a function of dimensionless object size 
$D/\lambda _0$ for different density ratios 
$\beta =\rho _o/\rho _f$ (see legend). Field scale is modelled by a 
$0.2\ \textrm {Hz}$ wave with a steepness of 
$\alpha =0.05$. Here, 
$C_m=\beta /2$. Analytical solutions using viscous drag from perturbation theory are denoted by solid lines.
Using the results from field-scale numerical simulations for non-viscous drag, a 
$1\ \textrm {m}$ diameter object of density 
$\rho _p = 0.9\ \text {g}\,\text {cm}^{-3}$ leads to a 
$50\,\%$ increase in drift (
$X^{(2)}=1.5$). This is a significant increase compared to the Stokes drift infinitesimal objects would experience. By comparison, a 
$0.1\ \textrm {m}$ diameter object in the same wave field does not experience any drift amplification (
$X^{(2)}=1$) and behaves as a perfectly Lagrangian tracer.
5. Conclusions
In this paper, we have developed a model for the transport of spherical, finite-size, floating marine debris by deep-water waves. Using a Stokes-like expansion in wave steepness, we have derived closed-form solutions for the linear response and the wave-induced drift of an object forced by regular waves and experiencing viscous drag. These closed-form solutions match numerical solutions of our model in the case of viscous drag. Our model recovers the Lagrangian limit as object size tends to zero, meaning that small objects are simply transported with the Stokes drift of surface gravity waves.
Through our perturbation solutions, we have identified two mechanisms for increased drift. The first arises from the change in magnitude of the linear orbits, especially its vertical component. The second arises when an out-of-phase variable submergence is resolved in the horizontal direction by the slope of the free surface. The second mechanism requires buoyancy and drag to be acting normal to the free surface, where the drag is required to create the phase difference that gives rise to the drift when averaged over the wave cycle. In any realistic oceanographic scenario, an non-viscous drag is required in order for the drift amplification to be significant. To observe the predicted effect, we have carried out laboratory wave flume experiments for a range of object sizes and densities (see Appendix B). The experiments show that an increase in wave-induced drift occurs. However, due to large experimental error, the present results have not been used to validate the theoretical model or choice of physics contained within.
The main driver for an increased drift is predicted to be an object's size relative to the wavelength. Thus, in the real ocean, where wavelengths range from 
$10$ to 
$10^3\ \textrm {m}$, increased drift will likely only be observed where shorter wavelengths are present, such as in gulfs or smaller seas. Modelling an object with a diameter of 
$1\ \textrm {m}$ and density of 
$0.9\ \text {g}\,\text {cm}^{-3}$ floating on a wave with a 5 s period and a steepness of 
$\alpha \equiv k_0 a_0=0.05$, typical of a moderately steep wind wave, results in a 
$50\,\%$ increase in wave-induced drift compared to the Stokes drift for such a wave. In the same wave field, an object with a diameter of 
$0.1\ \textrm {m}$ would not experience an increase in drift at all. High-quality experiments are recommended at larger scale, covering a wider range of object sizes and considering the effect of object shape. Insights from the present work should be useful in the development of more sophisticated models for tracking floating marine litter.
Acknowledgements
T.S.v.d.B. acknowledges a Royal Academy of Engineering Research Fellowship.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Equations of motion
Substituting (2.8) and (2.9) into (2.7), and (2.7) into (2.1) results in two second-order differential equations in the (
$n,\tau$) coordinate system
\begin{align} &\ddot{\tau}_p -(-(\partial_x\eta|_{x_p})^2 \varXi_p + \partial_{xx}\eta|_{x_p} \varXi_p^2 n_p) \ddot{x}_p= \frac{1}{m\left(1+\dfrac{C_{m,\tau}}{\beta}\right)} F_\tau +\{2 \dot{\theta}_p\dot{n}_p +(\dot{\theta}_p)^2 \tau_p \nonumber\\ &\quad -\partial_x \eta|_{x_p} \varXi_p(\partial_{tt} \eta|_{x_p} +2 \dot{x}_p \partial_{tx}\eta|_{x_p} +(\dot{x}_p)^2 \partial_{xx}\eta|_{x_p}) +n_p[(\partial_{tx}\eta|_{x_p} +\dot{x}_p \partial_{xx} \eta|_{x_p}) 2 \varXi_p \dot{\varXi}_p\nonumber\\ &\quad +(\partial_{ttx}\eta|_{x_p} + 2\dot{x}_p \partial_{txx} \eta|_{x_p} +(\dot{x}_p)^2\partial_{xxx}\eta|_{x_p}) \varXi_p^2]\}, \end{align}
\begin{align} &\ddot{n}_p +(\partial_x\eta|_{x_p} \varXi_p + \partial_{xx}\eta|_{x_p} \varXi_p^2 \tau_p)\ddot{x}_p= \frac{1}{m\left(1+\dfrac{C_{m,n}}{\beta}\right)}F_n -\{2 \dot{\theta}_p\dot{\tau}_p -(\dot{\theta}_p)^2 n_p \nonumber\\ &\quad + \varXi_p(\partial_{tt} \eta|_{x_p} +2 \dot{x}_p \partial_{tx}\eta|_{x_p} +(\dot{x}_p)^2 \partial_{xx}\eta|_{x_p}) +\tau_p[(\partial_{tx}\eta|_{x_p} +\dot{x}_p \partial_{xx} \eta|_{x_p}) 2 \varXi_p \dot{\varXi}_p\nonumber\\ &\quad -(\partial_{ttx}\eta|_{x_p} + 2\dot{x}_p \partial_{txx} \eta|_{x_p} +(\dot{x}_p)^2\partial_{xxx}\eta|_{x_p}) \varXi_p^2]\}, \end{align}
where we have kept all the second-order time derivatives on the left-hand side. We now have two equations in terms of three second-order time derivatives, namely 
$\ddot {\tau }_p$, 
$\ddot {n}_p$ and 
$\ddot {x}_p$, and require a third equation to solve the system. We obtain this third (kinematic) equation by taking the dot product of (2.7), in which we have substituted for 
$\ddot {\theta }_p$ and 
$\ddot {\eta }_p$ from (2.8) and (2.9), and 
$\boldsymbol {e}_x$, giving
\begin{align} &\ddot{x}_p [ 1
+\partial_{xx}\eta|_{x_p} \varXi_p^3 ( n_p+\partial_x
\eta|_{x_p} \tau_p ) ] - \ddot{\tau}_p \varXi_p
+ \ddot{n}_p \partial_x \eta|_{x_p} \varXi_p\nonumber\\
&\quad = \varXi_p\{{-}n_p[ (\partial_{tx}\eta|_{x_p}
+\dot{x}_p \partial_{xx} \eta|_{x_p}) 2 \varXi_p
\dot{\varXi}_p +( \partial_{ttx}\eta|_{x_p} +
2\dot{x}_p \partial_{txx}
\eta|_{x_p}\nonumber\\
&\qquad +(\dot{x}_p)^2\partial_{xxx}\eta|_{x_p} )
\varXi_p^2 ] - 2\dot{\theta}_p\dot{n}_p
-(\dot{\theta}_p)^2 \tau_p -\partial_x \eta|_{x_p}[
\tau_p[ (\partial_{tx}\eta|_{x_p} +\dot{x}_p
\partial_{xx} \eta|_{x_p}) 2 \varXi_p
\dot{\varXi}_p\nonumber\\
&\qquad+( \partial_{ttx}\eta|_{x_p} +
2\dot{x}_p \partial_{txx} \eta|_{x_p}
+(\dot{x}_p)^2\partial_{xxx}\eta|_{x_p} ) \varXi_p^2
] +2\dot{\theta}_p \dot{\tau}_p -(\dot{\theta}_p)^2
n_p] \}. \end{align}
Appendix B. Wave flume experiments
B.1. Set-up and data acquisition
A series of object tracking experiments were conducted in the Sediment Wave Flume in the Coastal, Ocean and Sediment Transport (COAST) Laboratory at the University of Plymouth, UK. The flume has length 
$35\ \text {m}$, width 
$0.60\ \text {m}$ and was filled with water to 
$0.50\ \text {m}$ depth, as shown in figure 12. A double-element piston-type wavemaker supplied by Edinburgh Designs Ltd was used to generate a wave packet with a spectral shape that linearly focuses to a Gaussian packet, 
$A_0=a_0\exp (-(x_f-c_{g,0}t)^2/2\sigma ^2)$, at a measurement zone centred 
$x_f=9.75\ \text {m}$ from the rest position of the wavemaker. The wave packet was made as long as possible to make it quasi-monochromatic whilst avoiding reflection (
$\epsilon = 1/(k_0 \sigma )=0.04$) with a steepness 
$\alpha = a_0 k_0=0.1$ and peak frequency 
$f_0=1.25\ \textrm {Hz}$.
Figure 12. Experimental set-up used to track the motion of floating objects under wave motion generated by a double-element piston-type wave maker at the COAST Laboratory, University of Plymouth, UK.
Despite our perturbation theory solutions being for periodic waves, we used quasi-monochromatic wave packets in our laboratory experiments because wave-induced transport is much easier to measure experimentally for wave packets (see van den Bremer, Yassin & Sutherland Reference van den Bremer, Yassin and Sutherland2019; Calvert et al. Reference Calvert, Whittaker, Raby, Taylor, Borthwick and van den Bremer2019 and the discussion in Monismith Reference Monismith2020). In Appendix C, we confirm that the slow modulation associated with the wave packet does not result in any additional non-inertial behaviour of the object. As a result, our model predictions for periodic waves and the wave packets considered in our experiments are equivalent.
We controlled the wavemaker using linear wave theory. Although sub-harmonic error waves at second order generated for wave packets (e.g. Nielsen & Baldock Reference Nielsen and Baldock2010; Orszaghova et al. Reference Orszaghova, Taylor, Borthwick and Raby2014) can lead to spurious wave-induced displacements (Calvert et al. Reference Calvert, Whittaker, Raby, Taylor, Borthwick and van den Bremer2019), these displacements are negligibly small for the deep-water waves we consider (van den Bremer et al. Reference van den Bremer, Yassin and Sutherland2019).
Seven resistance-type wave gauges provided 
$128\ \text {Hz}$ free surface elevation measurements. Five gauges were located close to the focus location at 15 cm intervals, as shown in figure 12. Two gauges were located significant distances before and after the focus location. After propagating through the measurement zone, the dispersed wave packets were absorbed by mesh-filled wedges within an absorption zone located at the downstream end of the wave flume. To ensure near-quiescent initial conditions for each experiment, the water surface was allowed to settle for 
$10\ \text {min}$ between experiments. A Photron SA4 high-speed camera captured the object motions at 
$125\ \text {frames}\,\text {s}^{-1}$, resolution of 
$1024$ by 
$1024\ \text {pixels}$ and shutter speed of 
$1/125\ \text {s}$. Optical distortion was removed using 
$35\ \text {mm}$ chequerboard images and MATLAB's inbuilt image processing package.
B.2. Matrix of experiments
In the experiments, we selected a peak frequency of 
$f_0=1.25\ \textrm {Hz}$, corresponding to a wavelength of 
$\lambda _0=1.0\ \textrm {m}$ and non-dimensional water depth 
$k_0d=3.1$. We then varied systematically the diameter 
$D$ and the density 
$\rho _o$ of the spherical floating object, with values for the 16 experiments listed in table 2. Object size was limited by camera resolution and the MATLAB tracking algorithm. Density was varied by filling hollow spheres with different ratios of epoxy to glass micro-ball filler. Each experiment was repeated five times.
Table 2. Matrix of experiments listing dimensional object diameter 
$D$, object density 
$\rho _o$, non-dimensional object diameter 
$D/\lambda _0$ and density ratio 
$\beta = \rho _o /\rho _f$.
B.3. Data processing
B.3.1. Free surface elevation
Wave packets were created from narrow-banded spectra to allow frequency filtering to separate the linear and second-order sub-harmonic components in the wave gauge signal. A band-pass filter between 
$0.8f_0$ and 
$1.2f_0$ was used to extract the linear free surface elevation. The measured envelope 
$A_0$ was calculated using the Hilbert transform of the linear free surface elevation. Use of the measured envelope at the location where the trajectories were measured, to calculate purely Lagrangian displacement, accounted for any dissipation or nonlinear dispersion between the wavemaker and the zone of interest.
B.3.2. Object tracking
Profile images of the floating white spheres were illuminated from various angles and captured by the Photron camera. The trajectories of the floating objects were tracked by identifying their position in each frame using a circle finding algorithm. The apparent size of the circle in the image was used to calibrate the pixel scale against the known size of the sphere. This also reduced any errors from out-of-plane motion not captured by the single camera. The horizontal components of the raw trajectories, repeated five times, are shown in figure 13.
Figure 13. Time histories of object horizontal position for each experiment. Each panel shows the five repeated experiments in different colours.
Every effort was made to settle the sphere at the start of each experiment in order to give it a zero initial velocity. This was not completely possible due to air flows over the water surface and slight disturbance from human touch. A linear fit in the time domain, assuming a constant pre-existing drift velocity, was used to remove motion before the arrival wave packet from the raw orbits in figure 13. The focus location was determined as coinciding with the position of the maximum of the linearised vertical motion envelope of the object. The difference in object location and exact focus location in the flume had negligible effect because of the very long wave packets used.
The magnitudes of the linear response were determined by filtering the horizontal and vertical motion components with a band-pass filter of 
$0.8$–
$1.2f_0$, followed by a Hilbert transform to obtain the envelope 
$A_0$. Note that frequency filtering was only applied to velocities, and numerical integration was used to calculate displacements. The maximum magnitude of the envelope was then normalised by wave amplitude 
$a_0$ to obtain 
$X^{(1)}$ and unity subtracted from the normalised vertical motion to give 
$\mathcal {S}^{(1)}$ (the normal and vertical directions equivalent up to first-order accuracy). We were not able to extract the linear phase from the experiments because exact spatial and temporal matching of Eulerian wave gauge data and Lagrangian object positions could not be achieved. A low-pass filter at 
$0.5 f_0$ was used to extract the sub-harmonic horizontal velocity component. The drift value 
$X^{(2)}$ was then determined by subtracting the Eulerian return flow from the maximum value of the sub-harmonic horizontal velocity component flow and dividing by the Stokes drift.
B.4. Comparison between theory and experiments
B.4.1. First order in wave steepness: ${O}(\alpha )$
Figure 14 presents the first-order magnitudes 
$|X^{(1)}|$ and 
$|\mathcal {S}^{(1)}|$ as functions of dimensionless diameter (
$D/\lambda _0$) for each experiment, with colour corresponding to density ratio. Comparison is made with numerical solutions of our model for non-viscous drag and analytical solutions using viscous drag. Overall, the horizontal motion in figure 14(a) is of similar magnitude to what is theoretically predicted (
$X^{(1)}$) with some variability, as quantified by the error bars. We note that a decrease of a few per cent in the numerical simulation solutions to 
$|X^{(1)}|$ is equivalent to a (small) dimensional decrease in the horizontal motion less than 
$1\ \textrm {mm}$. The first-order variable submergence 
$|\mathcal {S}^{(1)}|$ in figure 14(b) increases monotonically with dimensionless diameter (
$D/\lambda _0$), as predicted by theory.
Figure 14. Magnitude of the first-order motion as a function of non-dimensional object size 
$D/\lambda _0$ for different density ratios (see legend): analytical solution with viscous drag (solid lines) and experiments (circles). The density ratios for the numerical solutions are listed in the legend; density ratios for the experiments are labelled using the same colour scale. The error bars are obtained from repeated experiments and correspond to two standard deviations.
The experiments do not show a consistent trend with density for either linear motion component. We note that the densities are not equally spaced or the same for each size sphere owing to practical constraints on filling the spheres with different ratios of epoxy to glass micro-ball filler (see table 2 for the experimental matrix). The error bars shown for each experiment, which are twice the standard deviation of the five repeats, are large enough to mask any trend in density. Although we could measure the overall density of the spheres accurately, we emphasise that we were not able to measure its uniformity within the sphere.
Errors could have arisen from various physical sources that can account for the relatively large standard deviations. The initial motion of the object was hard to eliminate. Air conditioning was switched off, but there were occasional air flows over the flume. The method of taking the value of sub-harmonic velocity at the peak of the wave packet has been shown numerically to match regular waves in Appendix C. However, inertia at packet scale can be seen in figure 15 as the velocity does not go to zero after the packet passes. Although a 10 minute delay was prescribed between experiments to allow water in the flume to settle, there may have been residual currents still present. The theoretical model also has uncertainty, as can be seen in the sensitivity analysis in Appendix E, which arises from the choice of drag and added-mass formulations, and the exclusion of certain physics from the model, such as surface tension.
Figure 15. Sub-harmonic horizontal object velocity relative to the Eulerian mean flow, normalised by the Stokes drift at the centre of the wave packet: 
$X^{(2)}_{{exp}}=(v_x^{(2)|_{t=0}}-u_x^{(2)}|_{t=0})/(u_{{s}}|_{t=0})$ where 
$u_{{s}}|_{t=0}=\omega _0 k_0 a_0^2$. The mean of the five repeated experiments is shown as a continuous red line, and the confidence band corresponding to two standard deviations is shaded in grey, with five lines overlaid for each individual experiment.
B.4.2. Second order in wave steepness: ${O}(\alpha )$
Figure 15 presents time histories of the normalised sub-harmonic horizontal object velocity component for all 16 experiments, having first removed motion ahead of the wave packet and the Eulerian mean flow associated with the wave packet. In all cases, the non-dimensional sub-harmonic horizontal object velocity exceeds or is very close to unity near focus, and has a Gaussian-like profile, reducing close to zero within approximately 25 s either side of focus. The distributions are slightly skewed, with a faster rising limb than falling. There is more variability after focus than before. Using the peak values from figure 15, figure 16 shows the dimensionless drift factor 
$X^{(2)}$ for each experiment as a function of dimensionless diameter, with colour indicating density ratio. Drift increases with non-dimensional diameter and, as for the first-order results, the trend with density is unclear from the experiments and masked by substantial variability. We note that the density of floating plastic in the ocean typically has a small range between 
$800$ and 
$1000\ \text {kg}\,\text {m}^{-3}$ and may thus be a less important variable than object size. The trend with object size is consistent between experiments and theory, both presenting a similar increase with size.
Figure 16. Second-order drift amplification factor 
$X^{(2)}$ as a function of non-dimensional object size for different density ratios (see legend): analytical solution with viscous drag (solid lines) and experiments (circles). The density ratios for the numerical solutions are listed in the legend; density ratios for the experiments are labelled using the same colour scale. The error bars are obtained from repeated experiments and correspond to two standard deviations.
The experiments show that sufficiently large floating objects experience an increase in wave-induced drift. However, the experimental results are not sufficiently accurate to validate the theoretical model. In future work, it is therefore intended to carry out more experiments aimed at validating the model.
Appendix C. Wave packets vs. periodic waves
We use numerical solutions (see § 4) to the model developed in § 2 to examine the difference in predictions for objects subject to the quasi-monochromatic wave packets we use in our experiments and periodic waves. The processing of the trajectory data from the numerical simulations using wave packets was the same as for the experiments described in Appendix B. Figure 17 shows the almost identical first-order response as a function of non-dimensional object diameter at different density ratios for periodic waves (crosses) versus wave packets of the same bandwidth as in experiments (circles). Figure 18 shows the corresponding second-order drift amplification factors. Very slight differences are only predicted for larger object sizes for which the role of inertia is more dominant. For wave packets, a slightly smaller drift motion is predicted, because the time required for inertial objects to reach steady state is longer for larger objects.
Figure 17. Numerical predictions of the magnitude of the first-order horizontal motion amplification 
$X^{(1)}$ (a) and the variable submergence 
$\mathcal {S}^{(1)}$ (b) as functions of dimensionless object size 
$D/\lambda _0$ for non-viscous drag and for different density ratios 
$\beta =\rho _o/\rho _f$ (see legend). In the figure, periodic waves are denoted by crosses and wave packets of the same bandwidth as in the experiments by circles.
Figure 18. Numerical predictions of the magnitude of the second-order horizontal motion amplification 
$X^{(2)}$ as functions of dimensionless object size 
$D/\lambda _0$ for non-viscous drag and for different density ratios 
$\beta =\rho _o/\rho _f$ (see legend). In the figure, periodic waves are denoted by crosses and wave packets of the same bandwidth as in the experiments by circles.
Appendix D. Limiting behaviour of the numerical solutions
To confirm the model developed in § 2 is correct, including its cumbersome coordinate transforms, we examine the perfectly Lagrangian limit (§ D.1) and the small-object limit (§ D.2) of its numerical solutions obtained using MATLAB's ODE15s solver.
D.1. The Lagrangian limit
To obtain the Lagrangian limit, we replace the forces on the object by the accelerations a Lagrangian particle would experience under linear periodic waves
\begin{equation} \ddot{x}_p=a_0 \omega_0^2 \sin(\varphi) \exp\left(k_0 z_p\right), \quad \ddot{z}_p={-}a_0 \omega_0^2 \cos(\varphi) \exp\left(k_0 z_p\right), \end{equation}
where 
$\varphi = k_0 x_p - \omega _0 t+\varphi _0$. The accelerations are then mapped to the translating coordinate system and expressed in the (
$n,\tau$)-directions. The system is then solved numerically in (
$n,\tau$)-coordinates and the results mapped back onto (
$x,z$)-coordinates, providing confirmation our transformations are correct. As shown in figure 19, we obtain the correct amplitude of the vertical and horizontal linear motion and the correct Stokes drift.
Figure 19. Trajectory of a perfectly Lagrangian tracer obtained using a numerical solution of the present model with forcing provided by (D1). Panels (a,b) display the horizontal and vertical motions 
$x_p(t)$ and 
$z_p(t)$, with the blue dashed line showing the theoretical Stokes drift displacement and the red lines the superimposed wave amplitudes. Panels (c,d) show the tracer particle positions in the (
$n,\tau$)-coordinate system.
D.2. Small-object limit
As object size tends to zero, 
$D \rightarrow 0$, the solution should recover the behaviour of a perfectly Lagrangian tracer. This has been explicitly checked by numerically solving for an object of non-dimensional diameter 
$D/\lambda _0= 1 \times 10^{-6}$, which results in 
$X^{(1)}=1.00$, 
$\mathcal {S}^{(1)}=0.00$ and 
$X^{(2)}=1.00$.
Appendix E. Alternative drag and added-mass formulations
This appendix examines several alternative approaches to modelling the drag (§ E.1) and added-mass (§ E.2) forces on a floating object. Results are obtained from numerical solutions at laboratory-scale conditions as in § 4.
E.1. Drag
Although drag on a fully submerged sphere away from a free surface and in steady flow is well defined across a large range of Reynolds numbers (e.g. Morrison Reference Morrison2013), the drag force on a partially submerged, floating object in the unsteady flow field arising from surface waves is not. To understand the implications for our model's predictions, we consider the following drag formulations: viscous drag with 
$C_d=24/Re$, non-viscous drag with 
$C_d=C_d(Re)$ based on Morrison (Reference Morrison2013) and turbulent drag with 
$C_d=1/2$.
E.1.1. Viscous drag: $C_d=24/Re$
For the viscous drag coefficient 
$C_d=24/Re$, we consider three cases: a case based on submergence-dependent and thus time-varying projected area 
$\boldsymbol {A}_{PA}(t)=(A_{{s},n}(t),A_{{s},\tau }(t))$, as in the paper, a case that ignores the time dependence and sets 
$\boldsymbol {A}_{PA}=\boldsymbol {A}_{PA}^{(0)}\equiv \boldsymbol {A}_{PA}(s^{(0)})$, and a case that is based on the time-varying, direction-independent submerged surface area 
$A_{SA}(t)$. To compute the drag force, we use (2.14) and (2.16). For a sphere, the submerged surface area 
$A_{{SA}}(t)= {\rm \pi} D s(t)$. We normalise this by the surface area of a sphere 
$A_{{FS}}= {\rm \pi} D^2$, so that 
$\hat {A}_{{SA}}(t)=s(t)/D$ and replace both 
$\hat {A}_{s,\tau }$ in (2.14) and 
$\hat {A}_{s,n}$ in (2.16) by 
$\hat {A}_{{SA}}$. As a result of this normalisation, the drag forces on a fully submerged sphere based on projected area and based on submerged area are equal.
The first-order horizontal motion remains unchanged and so is not presented here. Variable submergence and second-order drift solutions are shown in figure 20. It is evident that inclusion of time-varying submergence in the projected area and replacing projected by submerged area has a negligible effect on the first-order submergence and only a very minor effect on the drift.
Figure 20. The effect of alternative drag formulations on the numerical predictions of first-order variable submergence 
$S^{(1)}$ (a) and second-order drift 
$X^{(2)}$ (b) as a function of non-dimensional object size 
$D/\lambda _0$ for a density ratio 
$\beta =0.8$ at laboratory-scale conditions. The lines correspond to different drag formulations, labelled in the legend, using either viscous drag (
$C_d= 24/Re$, solid lines), non-viscous drag (
$C_d=C_d(Re)$, dashed lines) or turbulent drag (
$C_d=1/2$, dotted lines), which either vary with the time-varying projected area in the respective directions (
$\boldsymbol {A}(t)=\boldsymbol {A}_{{PA}}(t)$), with the constant projected area in the respective directions (
$\boldsymbol {A}=\boldsymbol {A}_{{PA}}^{(0)}$), or with the submerged surface area (
${A}(t)={A}_{{SA}}(t)$).
E.1.2. Non-viscous drag: $C_d=C_d(Re)$
For the non-viscous drag coefficient, which is based on a fit to experimental data for a fully submerged sphere (4.2) (from Morrison Reference Morrison2013), we consider two cases. First, we set the drag to be proportional to the submergence-dependent, time-varying projected area 
$\boldsymbol {A}_{{PA}}(t)$, which is the approach used in the paper. Second, we ignore the time dependence and use the projected area of the sphere without waves 
$\boldsymbol {A}_{PA}=\boldsymbol {A}_{PA}^{(0)}\equiv \boldsymbol {A}_{PA}(s^{(0)})$.
Again, the first-order horizontal motion is unchanged and not presented here. The magnitude of the variable submergence and the drift are presented in figure 20. The variable submergence responses in these two cases are very similar to each other and to the viscous drag cases discussed above. The solutions for drift are similar to the viscous solution for small objects, diverging as the object size increases. For larger objects, the drift is significantly larger than when modelled with viscous drag. This is caused by the relative reduction in the drag force. There is a slight increase in drift when the projected area is time dependent.
E.1.3. Turbulent drag: $C_d=1/2$
We capture the turbulent drag limit by setting 
$C_d=1/2$, which we consider to be the practical large-object limit of (4.2). We consider two cases; similar to non-viscous drag, we have used the time-dependent projected areas 
$\boldsymbol {A}_{PA}$ and also consider time-independent projected areas of a sphere in the absence of waves 
$\boldsymbol {A}_{PA}^{(0)}$.
Again, the linear horizontal motion is unchanged and so not presented. The variable submergence is slightly decreased when compared with the viscous and non-viscous cases for larger object sizes, which results in a smaller adjusted Stokes drift. The increase in drift is larger than the viscous cases because of the relative reduction in drag, but smaller than the non-viscous cases. The comparative increase observed when using time-dependent submerged projected area, seen for non-viscous drag, can also be observed with turbulent drag.
E.2. Added mass
Maxey & Riley (Reference Maxey and Riley1983) derived the added mass for a fully submerged sphere in a low Reynolds regime and found the added-mass coefficient to be 
$C_m=1/2$. Hulme (Reference Hulme1982) studied a floating hemisphere under wave forcing and derived independent surge and heave added-mass coefficients as functions of non-dimensional object size 
$k_0 D/2$. The range of non-dimensional object sizes in the present study is 
$0 < k_0D/2\le 0.16$, which corresponds to added-mass coefficients in the range 
$0.83 \le C_{m,n}\le 0.86$ in heave and 
$0.5 \le C_{m,\tau }\le 0.53$ in surge (Hulme Reference Hulme1982).
We consider two categories of added-mass formulations: direction independent and dependent. In the first category, we consider 
$C_m=0$, 
$C_m=0.5$ representative of a submerged sphere in a low Reynolds regime, and 
$C_m= 0.5 \beta$ for an added mass that increases linearly with depth of submergence in the absence of waves but remains time independent. In the second category, we consider constant 
$C_m= (0.53,0.83)$ representative of a hemisphere (Hulme (Reference Hulme1982)), 
$C_m= 2 \beta (0.53,0.83)$ so that the added mass recovers Hulme's (Reference Hulme1982) result for a hemisphere and is zero for an entirely unsubmerged sphere. Finally, we extend this to a submergence and time-dependent added mass: 
$C_m= 2(0.53,0.83) s(t) /D$.
As for the different drag formulations, the first-order horizontal motion is insensitive to our choice of added-mass formulation. Figure 21 shows the first-order variable submergence and drift responses obtained for the different added-mass formulations considered. Figure 21(a) shows the relative insensitivity of the variable submergence response to the different added-mass formulations. The variable submergence exhibits a slight increase when the added mass is directionally dependent and a function of submergence. Drift, shown in figure 21(b), is more sensitive to the choice of added-mass formulation. Direction-independent formulations result in a smaller increase in drift compared to their direction-dependent counterparts. The smallest increase in wave-induced transport (excluding the special case of zero added mass 
$C_m=0$) is 
$C_m =0.5\beta$ which is used to generate the analytical and numerical solutions presented in the paper.
Figure 21. The effect of alternative added-mass formulations on the numerical predictions of first-order variable submergence 
$S^{(1)}$ (a) and drift 
$X^{(2)}$ (b) as a function of non-dimensional object size 
$D/\lambda _0$ for a density ratio 
$\beta =0.8$ at laboratory-scale conditions for non-viscous drag. The lines correspond to different added-mass formulations, described in the legend, with solid lines for directionally independent added-mass formulations, and dashed lines for added-mass formulations decomposed into normal and tangential directions.
