1. Introduction
Floating marine objects, moored, propelled or freely floating, are all exposed to and influenced by the ocean environment. These objects vary greatly in size, shape and density. The assessment of wave-induced drift of floating objects in the ocean is of importance for environmental and offshore engineering alike (Arikainen Reference Arikainen1972; Wilson Reference Wilson1982; Perrie & Hu Reference Perrie and Hu1997; Law & Huang Reference Law and Huang2007; Webb & Fox-Kemper Reference Webb and Fox-Kemper2011; Meylan et al. Reference Meylan, Yiew, Bennetts, French and Thomas2015; van den Bremer et al. Reference van den Bremer, Whittaker, Calvert, Raby and Taylor2019; Monismith Reference Monismith2020). Recently, there has been much interest in the topic because of concerns about marine plastic pollution (e.g. Law et al. Reference Law, Morét-Ferguson, Maximenko, Proskurowski, Peacock, Hafner and Reddy2010; van Sebille et al. Reference van Sebille2020).
An unrestrained object floating in a surface gravity wave field will normally experience a net drift in the direction of wave propagation, known as the Stokes drift (Stokes Reference Stokes1847), in addition to the oscillatory motion associated with the waves. This net drift typically only becomes relevant over long time scales due to its small magnitude (typically, of a few 
${\rm cm}\ {\rm s}^{-1}$ in the ocean). Unlike a perfectly Lagrangian tracer, whose drift is equal to the Stokes drift in the absence of Eulerian-mean flows, an object of finite size may display a different behaviour, and a velocity difference between the object and an idealized (i.e. Lagrangian) fluid parcel may emerge (Santamaria et al. Reference Santamaria, Boffetta, Afonso, Mazzino, Onorato and Pugliese2013; Meylan et al. Reference Meylan, Yiew, Bennetts, French and Thomas2015; Calvert et al. Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021; DiBenedetto, Clark & Pujara Reference DiBenedetto, Clark and Pujara2022).
The drift of small floating objects in periodic waves was investigated experimentally by Nath (Reference Nath1978). For small wave amplitudes, Lagrangian drift behaviour was found for very small objects, while for a spar-type drifting buoy with a deep draft, enhanced drift compared with the Stokes drift was reported. Huang, Law & Huang (Reference Huang, Law and Huang2011) explored the drift motion of objects of different shapes with two different submergence depths. Enhanced drift was found for all shapes, and objects with a larger submergence depth experienced a greater increase in drift regardless of shape. The studies of Tanizawa, Minami & Imoto (Reference Tanizawa, Minami and Imoto2002) and He, Ren & Qiu (Reference He, Ren and Qiu2016) showed that small objects behave like Lagrangian particles, following the Stokes drift, while large objects drift faster than Lagrangian particles with wave reflection off the object evident.
Theoretical models developed for wave-induced loads can be grouped into two main categories: models that take the object to be part of the boundary of the fluid domain allowing for calculation of diffraction effects based on potential-flow theory (Haskind Reference Haskind1946; Faltinsen & Løken Reference Faltinsen and Løken1979; Chen Reference Chen1994; Stansberg & Kristiansen Reference Stansberg and Kristiansen2011; Pessoa & Fonseca Reference Pessoa and Fonseca2015), and a second class of models that express loads in terms of the velocity field in the absence of the object considering both viscosity (drag) and fluid inertia, for example, using Morison's equation (Morison, Johnson & Schaaf Reference Morison, Johnson and Schaaf1950; Shen & Zhong Reference Shen and Zhong2001; Grotmaack & Meylan Reference Grotmaack and Meylan2006; Huang, Huang & Law Reference Huang, Huang and Law2016).
For objects of small size relative to the incident wavelength, the disturbance of the wave field by the object can be neglected, and thus, Morison's equation can provide an acceptable approximation. Morison's equation is applied by considering the motion of the object at its centre of mass and calculating the total force due to the waves as the sum of the inertial force, including the effect of added mass, and the drag force (caused by the velocity difference between the object and surrounding water). Rumer, Crissman & Wake (Reference Rumer, Crissman and Wake1979) conducted pioneering work by extending Morison's equation to study wave-induced drift of small floating objects including inertia, buoyancy, added mass and drag effects. The concept underlying their approach is to regard the free surface as an oscillating slope. A (dynamic) force balance normal to the free surface is achieved through the combined effect of a gravity force component and buoyancy, while the tangential component of gravity causes the drift motion of the object, and this is termed the slope-sliding effect (Rumer et al. Reference Rumer, Crissman and Wake1979). The slope-sliding concept has been applied and developed to study wave-induced motions of various objects by Shen & Ackley (Reference Shen and Ackley1991) and Huang et al. (Reference Huang, Huang and Law2016). They showed that a model that includes the slope-sliding term predicts enhanced drift but tends to underestimate the enhancement of the wave-induced drift of small floating objects compared with experiments. Also making use of a slope-sliding term, Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021) used a transformed coordinate system and employed perturbation methods to derive a closed-form solution for the drift of spherical floating objects. Enhanced drift motion is explained by two mechanisms in Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021). First, the magnitude of the linear motion (normal to the free surface) of a floating particle is enhanced compared with a Lagrangian particle, and, second, the dynamic buoyancy force has a net effect when averaged over the wave cycle in a similar fashion to the slope-sliding term of Rumer et al. (Reference Rumer, Crissman and Wake1979).
To accurately predict the drift when the object is large relative to the wavelength it is essential to account for the disturbance in the fluid field caused by the presence of the object. For models based on potential-flow theory, the fluid can be described by a velocity potential, which satisfies the Laplace equation subject to boundary conditions on the wetted body's surface as well as on the free surface, bottom boundary conditions and a radiation condition. When exposed to an incident wave field, objects experience forces and moments due to the waves. These encompass both unsteady forces, leading to oscillatory motion, and steady (or wave-averaged) forces arising from nonlinear effects. The steady forces, often referred to as drift forces, affect the magnitude and direction of objects’ drift, resulting in a slow and steady drift motion unless the object is moored (Suyehiro Reference Suyehiro1924; Watanabe Reference Watanabe1938; Havelock Reference Havelock1942). Two approaches to calculate second-order forces are highlighted here: the first solves for the far-field velocity potential and the second solves for the near-field velocity potential. Newman (Reference Newman1967) utilized conservation of momentum to relate the drift forces to the far-field potential and derived the horizontal steady second-order forces on a freely floating body in regular waves. The drift forces are found to differ considerably both in magnitude and sign depending on the wavelength and direction relative to the object. Pinkster & Hooft (Reference Pinkster and Hooft1976), Pinkster & Van Oortmerssen (Reference Pinkster and Van Oortmerssen1977) and Pinkster & Huijsmans (Reference Pinkster and Huijsmans1982) calculated the mean (or low-frequency) forces for different directions in regular and irregular waves by directly integrating the pressure distribution on the object. Their results show that the mean horizontal force due to the relative elevation between the object and surrounding waves can be significant in certain cases. Importantly, viscosity is not included in these calculations, but may have strong effects (Huse Reference Huse1977).
This paper examines the net drift of floating objects under the influence of unidirectional regular waves in deep water for objects that are sufficiently large to diffract the wave field. To do so, we use a hybrid numerical model that employs a fully nonlinear potential-flow model to capture the incident wave field and a Navier–Stokes (NS) model to calculate the detailed flow pattern near the object. Both viscous effects and (nonlinear) wave–body interactions are modelled. Objects with different sizes, drafts (submergence depths) and shapes in waves with different wave steepness are investigated in the presence and absence of viscosity with the objective of understanding the effect of these variables on drift and the mechanisms involved. To help explain our results, we propose a diffraction-modified Stokes drift. In this case, we use a simplified linear boundary element method (BEM) to generate the linear wave fields solving the (linear) wave–body interaction problem based on potential-flow theory, which in turn forms the basis of the diffraction-modified Stokes drift.
This paper is laid out as follows. We commence with the description of our two models, the hybrid numerical model qaleFOAM and the diffraction-modified Stokes drift model in § 2. The validation and verification of the hybrid numerical model are presented in § 3 (convergence tests are reported in Appendix B). In § 4, drift estimates based on the hybrid numerical model are presented for objects with a range of sizes, different steepnesses and two different shapes: rectangular boxes with round corners and round objects. The effects of viscosity are explored in each case. For the largest box, different submergence depths (drafts) and different radii of the round corners (corresponding to different submerged shapes) are simulated. In § 5 the results from non-viscous simulation of the hybrid numerical model are compared with the diffraction-modified Stokes drift model. Finally, conclusions are drawn in § 6.
2. Models
2.1. Problem formulation
We examine the wave-induced motion and drift of objects with different size, submergence depth and shape. We consider two shapes: a round-cornered box (RCB) and a round object (RO) with dimensions shown in figure 1. We define the size of the object as its length 
$l$ in the direction of wave propagation (for ROs, 
$l=D$ with 
$D$ the diameter), the submergence depth is 
$h_d$, the height of the RCB is 
$h$, and the radius of the rounded corner is 
$r$. Objects are placed in a regular incident wave field with wave amplitude 
$a_w$ and angular frequency 
$\omega$ in deep water (i.e. 
$kd>3$, where 
$k$ is the wavenumber and 
$d$ the depth of the fluid).
Figure 1. Shapes and dimensions of the two objects considered: RCBs and ROs.
Two models are used: the hybrid model qaleFOAM and a diffraction-modified Stokes drift model, which is solved based on the linearized potential-flow BEM. We use both models to conduct two-dimensional (2-D) simulations. In the hybrid qaleFOAM model, an inertial coordinate system (
$X,Z$) is chosen with its origin 
$O$ located at the bottom left corner of the fluid domain, with waves propagating from left to right, the 
$X$-axis positive in the direction of wave propagation, and the 
$Z$-axis positive upwards, as shown in figure 2(a). In the diffraction-modified Stokes drift model, we establish a Cartesian coordinate system (
$x,z$) with its origin 
$o$ located on the still-water level at the horizontal centre of the object, the 
$x$ axis in the direction of wave propagation, and the 
$z$-axis positive upwards, as shown in figure 3. Both coordinate systems, (
$X,Z$) and (
$x,z$), are inertial, earth-fixed coordinate systems and do not move with the objects. The only difference between these two coordinate systems is the position of the origin.
Figure 2. Domains, domain boundaries and coordinate system used in the hybrid numerical model qaleFOAM. (a) Schematic of the hybrid computational domain. (b) Boundaries of the NS domain.
Figure 3. Domain and coordinate system in the diffraction-modified Stokes drift model for the two objects considered, also showing object dimensions. Results are shown for a (a) rectangular box (RB) and (b) RO.
2.2. Hybrid numerical model: qaleFOAM
The hybrid numerical model qaleFOAM is used in this paper. The model is based on the domain-decomposition method, which couples the quasi arbitrary Lagrangian–Eulerian finite element method (QALE-FEM) potential-flow model with the two-phase incompressible NS model InterDyMFOAM in OpenFOAM. For details, see Ma & Yan (Reference Ma and Yan2010), Jacobsen, Fuhrman & Fredsøe (Reference Jacobsen, Fuhrman and Fredsøe2012), Li et al. (Reference Li, Wang, Yan, Gong and Ma2018), Gong et al. (Reference Gong, Yan, Ma and Li2020), Yan et al. (Reference Yan, Wang, Wang, Ma and Xie2020) and references therein. QaleFOAM has been applied to study various wave-structure interaction problems (Li et al. Reference Li, Wang, Yan, Gong and Ma2018; Yan et al. Reference Yan, Ma, Wang and Wang2019; Gong et al. Reference Gong, Yan, Ma and Li2020). The structures in these studies are either moored or self-propelled, and their sizes are at least 0.2 times the characteristic wavelength. Thus, the application of this model to smaller and unmoored objects (down to 0.01 times the wavelength) is new.
In the hybrid numerical model, the larger outer domain is solved by QALE-FEM to capture the incident waves; a smaller inner domain surrounding the object uses OpenFOAM to solve the NS equations, as shown in figure 2(a). In the NS model, both the air and water phases are assumed incompressible, and the volume-of-fluid method is used to identify the phases and their interface. The coupling approach employed in this paper is a one-way coupling, which means that at the interfaces, the NS model only takes the solutions of the QALE-FEM solver but does not feed its solutions back. The wave diffraction problem is thus solved in the NS domain, and we have to ensure that this domain is large enough so its finite size does not affect the solution, while not too large to become computationally prohibitive. By performing simulations with different domain lengths (in the range of 1–4 wavelengths), we demonstrate that our results (notably for object drift) are independent of the length of the NS domain. The left and right interfaces of the NS domain are equipped with passive wave absorbers (Yan et al. Reference Yan, Wang, Wang, Ma and Xie2020). The left interface of the NS domain is coupled with the QALE-FEM solver, where the waves generated in the QALE-FEM domain using a flap-type wavemaker are transferred into the NS domain. The boundaries of the NS domain are shown in figure 2(b). We note that to ensure the two-dimensionality of the simulations performed in this paper (using a numerical model that is, in principle, three dimensional), the front and back interfaces in the NS domain are not used. Furthermore, a laminar viscosity model is employed (i.e. we do not use a turbulence model).
Waves are generated on the left boundary in the QALE-FEM domain and absorbed on the right boundary. It takes time for waves generated by the wavemaker to propagate to the NS domain. In order to save computational cost, a reference time period 
${t_R}$ is set, during which the NS model is turned off. The tank length 
${L_x}$ for all simulations in this paper is chosen to be sufficiently long so that simulations finish before the reflected waves reach the object's location. The distance between the NS domain and the wavemaker is chosen to be at least 3 wavelengths in order to minimise the effects of evanescent waves from the wavemaker.
2.3. Diffraction-modified Stokes drift model
To provide an estimate of how object drift is affected by the diffraction of the wave field (without viscosity), we first use a BEM to solve the linearized potential-flow problem. From these linearized potential-flow solutions we obtain an estimate of the object drift (second order in steepness) in a fashion akin to Stokes (Reference Stokes1847) but by taking into account the modified wave field and object motion. To do so, we need to define all the boundaries of the fluid domain (see figure 3): 
$d$ is the depth of the fluid (at 
$z=-d$ a no-flow bottom boundary condition must be satisfied) and 
$x=\pm L_{BEM}/2$ correspond to the left and right boundaries of the fluid domain, where the radiation condition must be satisfied. We choose a value of 
$L_{BEM}/2$ that is large enough for the far-field truncation of the radiation condition not to affect our results. The boundaries 
$C_1$, 
$C_2$ and 
$C_3$ in figure 3(a) and 
$C_r$ in figure 3(b) require kinematic agreement of the fluid velocity and the (rigid) object's motion for the rectangular box (RB) and the RO, respectively. For the diffraction-modified Stokes drift model, we only consider a RB with square corners (i.e. 
$r=0$), whereas for the hybrid numerical model, we explore the effect of the radius of the rounded corner for the RCB.
At first order in steepness (i.e. for linear waves) the flow is described by a velocity potential 
$\varPhi$, which can be further divided into an incident potential 
${\varPhi _I}$, diffraction potential 
${\varPhi _D}$ and radiation potential 
${\varPhi _R}$, i.e.
\begin{align} \varPhi (x,z,t) &= {\varPhi _I}(x,z,t) + {\varPhi _D}(x,z,t) + {\varPhi _R}(x,z,t)\nonumber\\ &={{\rm{Re}}} \{ {\phi _I(x,z)}{{\rm e}^{ - \iota \omega t}}\} +{{\rm Re}} \left\{ {{\phi _D}(x,z){{\rm e}^{ - \iota\omega t}}} \right\}+ {{\rm Re}} \left\{ {{\phi _R}(x,z){{\rm e}^{ - \iota\omega t}}} \right\}, \end{align}
where all three components oscillate with the same angular frequency 
$\omega$. We denote the incident wave amplitude as 
${a_w}$, and the wavenumber 
$k$ is obtained from the linear dispersion relationship 
$\omega ^2=gk\tanh (kd)$, where 
$g$ is the gravitational acceleration. Although we will consider deep-water waves in this paper (i.e. 
$kd>3$ so that 
$\tanh (kd)\approx 1$), our diffraction-modified Stokes drift model is valid for general water depth. The time-invariant part of the incident wave potential 
${\phi _I}$ can be expressed as
\begin{equation} {\phi _I} = \frac{{{a _w}g}}{{\iota\omega }}\frac{{\cosh k(z + d)}}{{\cosh kd}}{{\rm e}^{\iota kx}}. \end{equation}
The boundary value problems for 
${\phi _R}$ and 
${\phi _D}$ are governed by the Laplace equation and solved using the Green's function method; the corresponding forces and equations of motions can then be found using standard methods (e.g. Newman Reference Newman2018). We use the implementation of these standard methods by Chen et al. (Reference Chen, Zhu, Zhao, Zhou and Fan2018), Yang, Zhu & Hong (Reference Yang, Zhu and Hong2019a) and Yang et al. (Reference Yang, Zhu, Hong and Huang2019b) (see Appendix A for details).
2.3.1. Estimating wave-induced object drift velocity
To obtain a leading-order (in steepness) estimate of object drift, we perform the same perturbation expansion, up to second order in wave steepness as Stokes (Reference Stokes1847) originally used to estimate the Stokes drift (see also van den Bremer & Breivik Reference van den Bremer and Breivik2018). Instead of only the linear incident wave field, we use the total linear wave field (cf. (2.1)) to estimate the ‘diffraction-modified Stokes drift’ for objects that are large enough to diffract the wave field:
\begin{align} u_{S,O} &= \overline{{\xi _x}\frac{{{\partial ^2}\varPhi }}{{\partial {x^2}}} + {\xi _z} \frac{{{\partial ^2}\varPhi }}{{\partial x\partial z}}} \nonumber\\ &=\overline{ {\xi _x}\left( {\frac{{{\partial ^2}{\varPhi _I}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\varPhi _R}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\varPhi _D}}}{{\partial {x^2}}}} \right) + {\xi _z}\left( {\frac{{{\partial ^2}{\varPhi _I}}}{{\partial x\partial z}} + \frac{{{\partial ^2}{\varPhi _R}}}{{\partial x\partial z}} + \frac{{{\partial ^2}{\varPhi _D}}}{{\partial x\partial z}}} \right)}. \end{align}
Here 
${\xi _x}={{\rm Re}} \{ {A_x}{{\rm e}^{ - \iota \omega t}}\}$ and 
${\xi _z}={{\rm Re}} \{ {A_z}{{\rm e}^{ - \iota \omega t}}\}$ are the linear horizontal and vertical harmonic oscillatory motions of the object, and the overline denotes averaging over the wave period. We term our estimate of the object drift in (2.3) the ‘diffraction-modified Stokes drift’, as it is based on the (linear) diffracted wave field. The diffraction-modified Stokes drift 
$u_{S,O}$ in (2.3) is calculated conceptually by (a leading-order estimate of) the difference between the Lagrangian-mean object speed and the Eulerian-mean fluid speed (i.e. 
$u_{S,O}=\bar {u}_{L,O}-\bar {u}_{E}$, cf. Bühler Reference Bühler2014). We do not compute the second-order Eulerian-mean flow, as this is generally very small for deep-water waves (e.g. van den Bremer & Taylor Reference van den Bremer and Taylor2015; van den Bremer et al. Reference van den Bremer, Whittaker, Calvert, Raby and Taylor2019) in the absence of the Earth's rotation, which is not considered here (cf. Higgins, Vanneste & van den Bremer Reference Higgins, Vanneste and van den Bremer2020), and thus, set 
$\bar {u}_{E}=0$. Below we show that we have negligible Eulerian flow in our simulations.
The different contributions to the diffraction-modified Stokes drift (2.3) can be made more explicit. We denote
$$\begin{gather} {\xi _x} = {A_x}{{\rm e}^{\iota{\theta _x}}},\quad {\xi _z} = {A_z}{{\rm e}^{\iota{\theta_z}}}; \end{gather}$$
$$\begin{gather}\frac{{{\partial ^2}{\phi _I}}}{{\partial {x^2}}} = k \omega {A_{xx}^I}{{\rm e}^{\iota{\theta _{xx}^I}}},\quad \frac{{{\partial ^2}{\phi _I}}}{{\partial x\partial z}} = k \omega {A_{xz}^I}{{\rm e}^{\iota{\theta _{xz}^I}}}; \end{gather}$$
$$\begin{gather}\frac{{{\partial ^2}{\phi _R}}}{{\partial {x^2}}} =k \omega {A_{xx}^R}{{\rm e}^{\iota{\theta _{xx}^R}}},\quad \frac{{{\partial ^2}{\phi _R}}}{{\partial x\partial z}} = k \omega {A_{xz}^R}{{\rm e}^{\iota{\theta _{xz}^R}}}; \end{gather}$$
$$\begin{gather}\frac{{{\partial ^2}{\phi _D}}}{{\partial {x^2}}} =k \omega {A_{xx}^D}{{\rm e}^{\iota{\theta _{xx}^D}}},\quad \frac{{{\partial ^2}{\phi _D}}}{{\partial x\partial z}} = k \omega {A_{xz}^D}{{\rm e}^{\iota{\theta _{xz}^D}}}, \end{gather}$$
where symbols 
$A$ denote the (potentially spatially dependent) magnitudes of the terms (given as amplitudes, in metres), and the phase and (oscillatory) spatial dependencies are captured by symbols 
$\theta$ with sub- and super-scripts on both 
$A$ and 
$\theta$ used as to indicate the different terms (and not derivatives). These amplitudes and phases can be obtained from the linear BEM model, which includes the equation of motion of the object. Now, diffraction-modified Stokes drift (2.3) can be rewritten as
\begin{align} \frac{u_{S,O}}{u_S}
&=\overline{{\rm Re} \{ {\hat{A}_x}{{\rm e}^{\iota({\theta
_x}- \omega t)}}\} [ {{\rm Re}} \{
{\hat{A}_{xx}^I}{{\rm e}^{\iota{(\theta _{xx}^I-\omega
t)}}}\} + {{\rm Re}} \{{\hat{A}_{xx}^R}{{\rm
e}^{\iota{(\theta _{xx}^R-\omega t)}}}\} + {{\rm Re}}
\{{\hat{A}_{xx}^D}{{\rm e}^{\iota{(\theta _{xx}^D-\omega
t)}}} \} ] + \dots }\nonumber\\ &\qquad \overline{
{{{\rm Re}} \{ {\hat{A}_z}{{\rm e}^{\iota{(\theta_z-\omega
t)}}}\} [ {{{\rm Re}} \{ {\hat{A}_{xz}^I}{{\rm
e}^{\iota{(\theta _{xz}^I-\omega t)}}}\} + {{\rm Re}} \{
{\hat{A}_{xz}^R}{{\rm e}^{\iota{(\theta _{xz}^R-\omega
t)}}}\} + {{\rm Re}} \{ {\hat{A}_{xz}^D}{{\rm
e}^{\iota{(\theta _{xz}^D-\omega t)}}}\} } ] }},
\nonumber\\ &=\overline{{\hat{A}_x}\cos ({\theta _x} -
\omega t)[ { {\hat{A}_{xx}^I}\cos ({\theta _{xx}^I} -
\omega t) + {\hat{A}_{xx}^R}\cos ({\theta _{xx}^R} - \omega
t) + {\hat{A}_{xx}^D}\cos ({\theta _{xx}^D} - \omega t)}
] + \dots}\nonumber\\ &\qquad
\overline{{\hat{A}_z}\cos ({\theta_z} - \omega t)[
{{\hat{A}_{xz}^I} \cos ({\theta _{xz}^I} - \omega t) +
{\hat{A}_{xz}^R}\cos ({\theta _{xz}^R} - \omega t) +
{\hat{A}_{xz}^D}\cos ({\theta _{xz}^D} - \omega t)}
]},\nonumber\\ &=\frac{{{\hat{A}_x}}}{2}[ {
{\hat{A}_{xx}^I}\cos ({\theta _x} - {\theta _{xx}^I}) +
{\hat{A}_{xx}^R} \cos ({\theta _x} - {\theta_{xx}^R}) +
{\hat{A}_{xx}^D}\cos ({\theta _x} - {\theta _{xx}^D})}
] \nonumber\\ &\quad +\dots
\frac{{{\hat{A}_z}}}{2}[ {{\hat{A}_{xz}^I}\cos
({\theta_z} - {\theta _{xz}^I}) + {\hat{A}_{xz}^R}\cos
({\theta_z} - {\theta _{xz}^R}) + {\hat{A}_{xz}^D} \cos
({\theta_z} - {\theta_{xz}^D})} ],
\end{align}
where the overline denotes wave averaging in time, upon which all the super-harmonic terms (of the form 
$\cos (2\omega t+\beta )$ with 
$\beta$ an arbitrary phase) disappear, and the symbol 
$\hat {A}$ denotes the normalization of the corresponding magnitude 
$A$ by the incoming wave amplitude 
$a_w$ (i.e. 
$\hat {A}=A/a_w$). The theoretical Stokes drift 
$u_S$ is given by (Stokes Reference Stokes1847)
\begin{equation} u_S = {{a_w}^2}\omega k{{\rm e}^{2kz}} = {\varepsilon ^2}\displaystyle\frac{\omega}{k}{{\rm e}^{2kz}}, \end{equation}
in which 
$\epsilon =ka_w$ is the incident wave steepness, where we have used 
$z=0$ in the normalization in (2.8).
3. Validation and verification of the hybrid numerical model (qaleFOAM)
In this section, validation and verification are conducted for the qaleFOAM hybrid numerical model (see Appendix A for validation and verification of the BEM model). To do so, we first examine the (Stokes) drift of a Lagrangian particle (i.e. a fluid parcel) through analysis of the Lagrangian-mean velocity (the Eulerian-mean velocity field as well as grid convergence are examined in Appendix B). We then examine the Lagrangian drift behaviour of small floating objects.
3.1. Drift of a Lagrangian particle
First, we consider, in turn, regular waves in deep water with two different frequencies. For each frequency, a series of waves are simulated with different wave amplitudes, and the horizontal drift velocities of fluid particles are calculated to confirm these are equal to the theoretical Stokes drift based on (2.9). The mean drift velocity of a fluid particle in quasi-steady state is obtained by applying the best linear fit to its horizontal trajectory and determining the slope of the linear fit line. The trajectories themselves are obtained from solving the ordinary differential equation 
${\rm d}\kern 0.06em {\boldsymbol {x}}_L/{\rm d}t={\boldsymbol {u}}({\boldsymbol {x}}_L(t),t)$, where 
${\boldsymbol {x}}_L(t)$ is the position of a Lagrangian particle.
The properties of the waves and the numerical parameters of the simulations are given in table 1, where 
${T_{dur}}$ and 
${T}=2{\rm \pi} /\omega$ refer to the total time duration of the simulations and the wave period, respectively, and 
$L_x$ is the (horizontal) length of the total domain. The parameters 
$\Delta x$, 
$\Delta z$ and 
$\Delta t$ denote the horizontal and vertical grid sizes and time steps, respectively. The horizontal positions 
${x_L}$ and 
${x_R}$ represent the left and right boundary locations of the NS domain, respectively, and 
${z_A}$ denotes the vertical location of the top of the air phase relative to the still-water level in the NS domain. Horizontal grid size is given as a fraction of the wavelength 
$\lambda$ and vertical grid size as a fraction of the linear wave amplitude 
$a_w$. Finally, the maximum Courant number 
${Co} = \Delta t| {\boldsymbol {u}} |/\Delta x=0.25$, in which 
$| {\boldsymbol {u}} |$ refers to the maximum absolute velocity.
Table 1. Wave properties and numerical parameter values for the hybrid numerical model (qaleFOAM) simulations, distinguishing the potential-flow (QALE-FEM) and the NS domains.
We use a Crank–Nicolson scheme for time integration and a non-uniform mesh with finer resolution close to the free surface in the 
$z$ direction in both the QALE-FEM and NS domains. Specifically, the grid density in the 
$z$ direction in QALE-FEM increases exponentially with distance to the free surface, and the vertical grid size is defined by the number of layers in the vertical direction, for which 20 are typically enough for deep-water simulations. The mesh sizes in table 1 all refer to those in the region near the free surface.
For InterFoam and InterDyMFoam solvers, Larsen, Fuhrman & Roenby (Reference Larsen, Fuhrman and Roenby2019) provide a detailed analysis of different combinations of discretization schemes, mesh sizes and Courant numbers for surface waves to maintain stable amplitudes over long times. For these solvers, Devolder et al. (Reference Devolder, Schmitt, Rauwoens, Elsaesser and Troch2015) reported instability of the added mass term and suggested how to choose the initial values of the added mass relaxation factor in order to obtain fast convergence and stable motion. Moradi, Zhou & Cheng (Reference Moradi, Zhou and Cheng2015), Palm et al. (Reference Palm, Eskilsson, Paredes and Bergdahl2016), Mohseni, Esperanca & Sphaier (Reference Mohseni, Esperanca and Sphaier2018) and Palm et al. (Reference Palm, Eskilsson, Bergdahl and Bensow2018) have investigated wave–body interaction. Based on the above, we choose the PISO algorithm to solve the pressure–velocity coupling, a limited second-order Crank–Nicolson scheme (implicit) with a blending factor of 0.9 is used for time integration (ddtSchemes), a minimally diffusive gradient limiter (cellMDLimited Gauss linear 1, which is second order and bounded) is used for gradients (gradSchemes) to avoid over and undershooting. To compute the divergence term (divSchemes), a second-order total variation diminishing (TVD) scheme (Gauss MUSCL) is used for the momentum convection term. A second-order and bounded TVD scheme (Gauss vanLeer01 with Gauss interfaceCompression) is used to compute the volume fraction.
We vary the steepness of the simulated waves of the two different frequencies from 0.03 to 0.13, and let the relaxation zone vary in length from 1 to 1.5 wavelengths as the wave steepness increases (Yan et al. Reference Yan, Ma, Wang and Wang2019). Correspondingly, the location of the right-hand side of the NS domain is adapted so that the length of the NS domain is equal to the relaxation zone length plus the necessary length for particles to move during the proposed time duration of the simulation. The initial position of the tracked particles is chosen to avoid disturbance by transition through the relaxation zone. The initial horizontal position 
${x}_{L,0}$ is chosen some distance to the right of the left relaxation zone, and the initial vertical position 
${z}_{L,0}$ is chosen immediately below the trough of the wave.
Figure 4(a) displays the horizontal trajectory 
${x_L}(t)$ of the tracked particles (
${x_{L0}}, z_{L0}$) for 
$\epsilon =0.034$, 
$\omega =4.09\ {\rm rad}\ {\rm s}^{-1}$ as an example. It is evident that in addition to the oscillatory motion of the waves, the Lagrangian particle undergoes a mean drift that agrees well with the theoretical Stokes drift according to (2.9), in which we set 
$z$ equal to the initial vertical position of the tracked particle 
${z_{L0}}$. A comparison between the numerical prediction of the mean drift and the theoretical Stokes drift for different steepnesses is shown in figure 4(b). Excellent agreement is achieved for both higher and lower frequencies and for a range of steepnesses. Finally, we confirm that the Eulerian-mean velocity is negligibly small everywhere in our domain in the case without viscosity (shown in Appendix B, so that 
$\bar {u}_L=u_S$). Together, this validates our model and confirms its ability to predict the drift velocity of an infinitely small object. The deviation from perfectly quadratic behaviour in figure 4(b) results from the initial vertical position 
$z_{L0}$ being chosen just below the trough for each steepness; this vertical position is also used to evaluate the theoretical Stokes drift according to (2.9), hence, the very good agreement.
Figure 4. Drift of Lagrangian particles in simulations of the hybrid numerical model (qaleFOAM) without viscosity (
$\nu =0$): (a) time history of the horizontal motion of a Lagrangian particle (
${x_{L_0}}=22.5$ m, 
${z_{L_0}}=-25.0$ mm) for 
$\epsilon =0.034$ for the lower-frequency case in table 1, where the theoretical Stokes drift is evaluated using (2.9). The wave amplitude here is 
$a_w=20.0$ mm; thus, the vertical position of the particle is always below the trough of the wave. The black dashed line denotes the time at which a quasi-steady state has been achieved and the drift speed has become constant (
$t_s=32$ s). The drift speed in (b) has been obtained from the average speed for 
$t>t_s$. (b) Comparison of non-dimensional drift velocities of a Lagrangian particle 
$\bar {u}_L/c$ between numerical solutions (red and blue squares) and theoretical Stokes drift (red and blue lines) as a function of wave steepness for higher and lower frequencies, where 
$c=\omega /k$ is the phase speed of the waves.
3.2. Drift of very small objects
To further validate our model, we examine whether it correctly predicts the drift of very small but finite-size objects, which should be equal to the Stokes drift (in the absence of Eulerian-mean flows). Clear experimental evidence exists that when the size of an object is very small, its behaviour is purely Lagrangian (Nath Reference Nath1978; van den Bremer et al. Reference van den Bremer, Whittaker, Calvert, Raby and Taylor2019). Two different object shapes are examined here: a RB and a RO. We choose the lower-frequency (
$\omega =4.09\ {\rm rad}\ {\rm s}^{-1}$) wave condition from table 1 and a small steepness of 
$\epsilon =0.034$. The round-cornered RB has a length (in the direction of wave propagation) of 
$l=0.036$ m (
$l/\lambda =0.97\,\%$), a draft 
${h_d}=0.025$ m and the radius of the rounded corner is 
$r=0.006$ m. The diameter of the RO is 
$D=0.05$ m 
$(D/\lambda =1.3\,\%)$. Both of the objects have a density of 
$\rho =500\ {\rm kg}\ {\rm m}^{-3}$. The height of the box and the radius of the RO are chosen to make sure that the object will not be submerged by the waves.
Table 2 gives the object drift velocities of both objects with and without viscous effects modelled in the simulation (i.e. 
$\nu =1.00\times 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$ and 
$\nu =0$). The total number of cells of the discretization mesh 
${N_c}$ reported in table 2 is the one for the non-viscous simulation, based on which six vertical layers are added and one level of refinement is applied near the object for the corresponding viscous simulation. The results confirm the drifts of both very small objects (
${\sim }1\,\%$ of the wavelength) are approximately equal to the theoretical Stokes drift in both the viscous and the non-viscous simulation. By comparing a coarser grid to a finer grid, table 2 also shows the convergence of several physical quantities (wave amplitude, horizontal viscous forces on the particle and the object drift velocity) for both the viscous and non-viscous simulations. We note the horizontal resolution changes with distance to the object in order to make the grid's aspect ratio approximately unity in the object's near field.
Table 2. Numerically predicted wave amplitude as a fraction of the input amplitude 
$a/a_w$ and object drift as a fraction of theoretical Stokes drift 
$u_{O}/u_{S}$ for very small objects of two different shapes and for different mesh sizes, with 
$l$ the length of the object. The horizontal grid size 
$\Delta x$ varies with distance to the object from the farthest location where 
$\Delta x=0.02$ m to the nearest location where 
$\Delta x=0.0005\unicode{x2013}0.001$ m. The total number of cells of the mesh is denoted by 
$N_c$. The amplitude 
$a$ is the time-averaged wave amplitude at the location of more than one wavelength downstream away from the object (near the outlet, defined in figure 2b) scaled by the input wave amplitude 
${a_w}$. The force 
$F'_{x,vis}=F_{x,vis}/(\rho g{a_w}l^2/4)$ is the non-dimensional magnitude of the horizontal viscous force on the object.
4. Results from the hybrid numerical model (qaleFOAM)
In this section we explore the role played by a variety of factors, namely size, shape, submergence depth, viscosity and wave steepness, in determining an object's drift behaviour. The wave field is equivalent to that considered before in § 3 (see table 1). The two object shapes considered along with the definitions of object dimensions are given in figure 1. We use RCBs instead of boxes with sharp corners to avoid generation of undesirable vorticity that would complicate the analysis (cf. Moradi et al. Reference Moradi, Zhou and Cheng2015). To organize our findings, we define the following dimensionless parameters: relative object size is described by 
$l/\lambda$ (for ROs 
$l=D$) and the radius of the corners of the RCB by 
$r/{h_d}$. The total duration 
${T_{dur}}$ is around 
$60$–
$80$ s for all simulations; it takes around 
$25$–
$30$ s for object drift to achieve a steady state, and another 
$25$–
$35{T}$ is sufficient to estimate the drift velocity. The spatial resolution of the mesh for all cases without viscosity lies in the range from 
$\Delta x=\lambda /200$, 
$\Delta z={a_w}/20$ at a location 
$1\unicode{x2013}1.5\lambda$ away from the object to 
$\Delta x={a_w}/20$, 
$\Delta z={a_w}/20$ within a distance of 2–4
$l$ surrounding the object. For the cases with viscosity, we add six vertical layers near the surface and apply one level of refinement near the object. The maximum Courant number is 
${Co}=0.25$.
We conduct three categories of simulations. In the first category (category I), we consider the effect of size for a RCB and a RO (§ 4.1). For RCBs, we keep the aspect ratio 
$h_d/l$ and submerged shape 
$r/h_d$ constant, and we consider only RCBs to avoid the effect of undesirable vorticity from sharp corners. To begin our analysis with a case of the simplest possible geometry, we set the object density 
$\rho =500\ {\rm kg}\ {\rm m}^{-3}$, so that both objects are exactly half-submerged. As the RCBs are hydrodynamically unstable for this density, we constrain the object rotation to be zero in the category I simulations, the consequences of which we discuss in § 4.1. In the second category (category II), we keep the size of the object constant but vary its submergence depth and submerged shape by changing the radius of the round corner (§ 4.2). Instead of changing the aspect ratio of the objects, we vary the submergence and roundness of the objects, which is intended to examine the effect of ‘streamlining’ of objects of constant size. The density of the objects in category II is different from category I, as we no longer wish to constrain the object's rotation; specifically, we choose a density (
$\rho =781\ {\rm kg}\ {\rm m}^{-3}$) that is high enough for RCBs to become hydrodynamically stable while maintaining the same size and aspect ratio 
$h_d/l$ as the objects in category I. Unlike the first and second categories, which are all conducted on low-steepness waves, in the third category (category III), we simulate drift behaviour for a range of wave steepnesses (§ 4.3). We examine the role of viscosity explicitly in § 4.4 for the experiments in categories I and II.
4.1. Effect of size (category I)
To study the effect of size, we vary the size of the two objects measured relative to the wavelength from 
$1\,\%$ to 
$10\,\%$ in 
$1\,\%$-point steps. Detailed object dimensions are given in table 3. For the RCB, we set 
$r/{h_d}=0.24$. The differences between a RCB and a RO of equivalent relative size are the submerged shape and submergence depth of the object. The simulations are performed with and without the inclusion of viscosity. Simulation results are shown in figure 5. The non-dimensional magnitude of oscillatory motion in the horizontal and vertical directions are shown in figures 5(a) and 5(b), respectively. The amplitudes of the oscillatory motions 
${A_x}$ and 
${A_z}$ are obtained by filtering out sub- and super-harmonic components, before obtaining amplitudes after the quasi-steady state has been achieved. Figure 5(c) shows the influence of sizes and shape on drift. Finally, figure 5(d) depicts the local wave amplitude as a function of its horizontal distance to the centre of mass of the objects for RCBs of relative sizes 
$l/\lambda =1\,\%$, 
$l/\lambda =9\,\%$ and 
$l/\lambda =10\,\%$ and ROs of relative size 
$l/\lambda =10\,\%$.
Table 3. For category I, simulations exploring the effect of size, object dimensions of the two different objects considered: RCBs and ROs (
$\rho =500\ {\rm kg}\ {\rm m}^{-3}$ for both).
Figure 5. Effect of object size on object motion and drift (category I simulations): (a) horizontal oscillatory motion amplitude 
${A_x}$, normalized by input wave amplitude 
${a_w}$, as a function of relative object size; (b) vertical oscillatory motion amplitude 
${A_z}$, normalized by wave amplitude 
${a_w}$, as a function of relative object size; (c) celerity-normalized object drift 
$u_{O}/c$ as a function of relative object size; (d) non-dimensional amplitude of the local surface elevation 
$a(x)/{a_w}$ as a function of horizontal distance (scaled by wavenumber) from the centre of mass 
${x_c}$ without viscosity (
$\nu =0$). The gap between the two vertical red lines in (d) represents the position of the object for 
$l/\lambda =10\,\%$ and corresponds to its left and right sides, respectively. The percentage (
$\%$) in (d) represents the relative size of the object 
$l/\lambda$. The results for Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021) in panels (a–c) are their results for spheres with an equivalent diameter to our ROs.
We start by examining the oscillatory motion for the RCB, shown in figures 5(a) and 5(b). When the object is very small, the magnitudes of oscillatory motion in both directions are very close to the incident wave amplitude 
${a_w}$, suggesting the object behaves as a Lagrangian particle. As object size increases, the horizontal oscillatory motion is reduced as an approximately linear function of relative size, while the vertical oscillatory motion increases nonlinearly at an increasing rate. For the object drift velocity in figure 5(c), we observe that when the box is very small, its drift rate is equal to the Stokes drift, while, as the box becomes larger, the drift speed is enhanced significantly. The enhanced drift is minimal for RCBs with a relative size less than approximately 
$7\,\%$, but becomes more evident for larger boxes. Significant increases in the drift for RCBs only become evident at greater relative size compared with increases in the vertical oscillatory motion. For completeness, we note the drift is slightly reduced compared with the theoretical Stokes drift for intermediate-size RCBs with a relative size of 
$3\,\% \leq l/\lambda \leq 7\,\%$ and ROs with 
$5\,\% \leq l/\lambda \leq 8\,\%$.
From the wave-field analysis in figure 5(d), a standing wave pattern emerges in the case of a large RCB (with large submergence depth). On the upstream side, the time-averaged wave amplitudes show a pattern of partial nodes and anti-nodes with smaller and larger amplitudes locally compared with the incident undisturbed wave amplitude (and compared with the 
$1\,\%$ relative size object, for which a standing wave pattern is not discernible). Wave amplitudes on the far downstream side for large objects are unaffected, while for locations near the object on the downstream side, smaller amplitudes are found. Differences in surface elevation between the two sides of the object become most evident for the larger RCBs. We note (numerical) wave gauges are set at a fixed spatial interval; thus, there may be small errors in determining maxima and minima. The (partial) standing wave pattern becomes difficult to discern for RCBs with a relative size smaller than 
$7\,\%$ (not presented here for brevity). All the above trends are similar whether viscosity is included or not; the drift of RCBs is enhanced further by viscous effects and even more so for larger objects, by up to 
$20\,\%$, as shown in figure 5(c). For practical computational reasons, the local wave amplitude is obtained by analysing the surface elevation in a stationary reference frame and not in the referencing frame moving with the object, in which the standing wave pattern most likely forms. As the object drift is always small relative to the phase speed (i.e. 
$u_{O}/c\leq 2.5\times 10^{-3}$), we anticipate the standing wave patterns in both reference frames to be similar albeit likely smaller and more diffused in the stationary reference frame shown here.
Comparing the RO and the RCB, both display a similar linear decrease in horizontal oscillatory motion with relative size (figure 5a), whereas the vertical motion of ROs is only enhanced by a very small amount compared with RCBs (figure 5b). Furthermore, for ROs, as depicted in figure 5(c), no obvious enhancement of the drift speed is found in the absence of viscosity even when the relative size is as large as 
$10\,\%$. In the presence of viscosity, a small amount of enhanced drift is observed as the RO becomes larger (relative size larger than 
$8\,\%$). The motion of ROs is thus distinctly different from that of RCBs, especially their vertical oscillatory motion and enhanced drift. To explain this, we note that the standing wave pattern in figure 5(d) for the largest RO (
$10\,\%$) is even smaller than for the 
$8\,\%$ relative size RCB (not shown in the figure). The standing wave pattern is an indication of the presence of a diffracted wave field; the extent to which diffraction occurs depends on the streamlining of the object. In § 6 the effect of shape is examined further.
Figure 6. Effect of the shape of the submerged part of a RCB as determined by the radius of the rounded corner 
$r$, scaled by submergence depth 
${h_d}$, on object motion and drift (category II simulations, group 1): (a,b) non-dimensional amplitudes of oscillatory motion in the horizontal and vertical directions, respectively; (c) celerity-normalized object drift velocity. Blue square markers in (a–c) represent the results for RCBs. Red lines in (a,b) denote the oscillatory motions of a Lagrangian particle, these in (c) denote the theoretical Stokes drift. (d) Non-dimensional amplitude of the local surface elevation 
$a(x)/{a_w}$ as a function of horizontal distance (scaled by wavenumber) from the centre of mass 
${x_c}$. The gap between the two vertical red lines in (d) represents the position of the object and corresponds to its left and right sides, respectively. The line denoted by FEM in (d) corresponds to simulations of the incident wave field only.
We note that in the above simulations (category I), we have constrained the object's rotation. This is necessary, as in keeping the object density constant at 
$\rho =500\ {\rm kg}\ {\rm m}^{-3}$, we have considered a RCB where height 
$h$ exceeds length 
$l$ (i.e. 
$h/l=1.33$). This is hydrodynamically unstable, and would normally start to rotate upon small perturbations from its vertically upright position. In presenting the results here, we have thus assumed the interaction between the motions in the different degrees of freedom (translation and rotation) is small. In Appendix C we keep the submergence depth and submerged shape of the RCB the same as in table 3 but change its density to properly explore the effects of rotation. We show that allowing rotation does not affect the conclusions presented in this section. In the following sections, we will allow rotation.
4.2. Effect of shape and submergence depth (category II)
Motivated by the difference in drift enhancement between box-shaped and round objects of equivalent, relatively large size in category I above, we proceed to examine how the shape and size of the submerged part of a RCB affect the standing wave pattern and the drift enhancement (category II). Unlike in category I, the objects in category II are allowed to rotate but, to simplify the analysis, we do not include viscosity. In all cases, the object size and density are kept constant at a relative size of 
$10\,\%$ and a density of 
$\rho =781\ {\rm kg}\ {\rm m}^{-3}$, respectively. Two groups of RCBs, one where each object has a different round-corner radius 
$r$ (group 1) and the other where each object has a different height submergence depth 
$h_d$ achieved through varying the height of the box 
$h$ (group 2), are simulated. For group 1, the boxes have the same height, and we vary the radius of the round corners to change their submerged shape. Object dimensions for group 1 are given in table 4. For group 2, we vary the submergence depth of the box by varying its height, keeping the radius of the round corner constant. Object dimensions for group 2 are given in table 5.
Table 4. For category II, simulations exploring the effect of object shape with fixed height (group 1): object dimensions for different shapes (
$l=0.37$ m, 
$h=0.32$ m, 
$\rho =781\ {\rm kg}\ {\rm m}^{-3}$).
Table 5. For category II, simulations exploring the effect of submergence depth with fixed radius of the round corners (group 2): object dimensions for different submergence depths (
$l=0.37$ m, 
$\rho =781\ {\rm kg}\ {\rm m}^{-3}$).
For group 1, the amplitudes of the oscillatory part of the motion in both the horizontal and vertical directions are given as a function of normalized round-corner radius 
$r/{h_d}$ in figures 6(a) and 6(b), respectively. Amplitudes are obtained in the same way as for category I simulations. Figure 6(c) presents the normalized time-averaged object drift speed 
$u_{O} / c$ as a function of normalized round-corner radius 
$r/{h_d}$ for group 1. Finally, figure 6(d) depicts the spatial wave amplitude distribution as a function of the wavenumber-normalized distance from the centre of mass of the object for group 1. Figure 7 gives analogous results for group 2.
Figure 7. Effect of the submergence depth 
${h_d}$ of a RCB on object motion and drift (category II, group 2): (a,b) non-dimensional amplitudes of oscillatory motion in the horizontal and vertical directions, respectively; (c) celerity-normalized object drift velocity. Blue square markers in (a–c) represent the results for RCBs. Red lines in (a,b) denote the oscillatory motions of a Lagrangian particle, these in (c) denote the theoretical Stokes drift. (d) Non-dimensional amplitude of the local surface elevation 
$a(x)/{a_w}$ as a function of horizontal distance (scaled by wavenumber) from the centre of mass 
${x_c}$. The gap between the two vertical red lines in (d) represents the position of the object and corresponds to its left and right sides, respectively. The line denoted by FEM in (d) corresponds to simulations of the incident wave field only.
We begin examining the influence of shape by returning to the results for category I. As shown in figure 5(a), RCBs and ROs experience a similar linear decrease of the horizontal oscillatory motion amplitude with relative size. The amplitude of the vertical oscillatory motion of ROs increases much less with relative size compared with RCBs (figure 5b). The difference in object drift becomes large when the relative size is larger than approximately 
$7\,\%$ (figure 5c). We note that increased drift is always accompanied by an increase in amplitude of the vertical motion.
We now turn to the simulations in category II, group 1, in which we vary the radius of the rounded corners. Figure 6(c) shows that as the radius of the rounded corner becomes larger, which corresponds to a more rounded shape, the drift speed decreases. So does the amplitude of the vertical motion (figure 6b). The amplitude of the horizontal oscillatory motion increases by only a small amount with increasing radius (figure 6a). Furthermore, the standing wave pattern becomes less apparent with increasing radius (figure 6d). Accordingly, the wave amplitude difference between the two sides of the object decreases. To sum up, figure 6 provides strong support for the idea that the increase in object drift compared with the Stokes drift is related to the standing wave pattern and is determined by how ‘streamlined’ the object is. We note that even for the RCB with the largest rounded-cornered radius, the enhanced drift is still significant, which is due to its large submergence depth, as we will examine next.
For RCBs with different submergence depths (category II, group 2), it is evident from figure 7 that as the submergence depth increases, the object drift increases significantly, as does the amplitude of the oscillating vertical motion. The horizontal oscillatory motions increase by only a small amount with increasing submergence depth. Figure 7(d) reveals that the increase in object drift is accompanied by an increase in the standing wave pattern. The largest wave amplitude on the upstream side and the relative difference in wave amplitudes between both sides of the object both increase as the submergence depth increases, further supporting our finding that the standing wave pattern drives enhanced drift.
Taking the above analysis of the effects of size, shape (of the submerged part of the object) and submergence depth together, the role of the standing wave pattern generated by the object and the relative wave amplitude difference between the two sides of the object stands out. All these effects described above can be understood in terms of the ability of the object to ‘hinder’ the flow pattern associated with the incident wave field. The larger the submergence depth and the more box-like the submerged shape, the more the objects hinder the flow. Enhanced drift is always accompanied by a (small) reduction in horizontal oscillatory motion and a (large) increase in vertical oscillatory motion.
4.3. Effect of wave steepness (category III)
Simulations in category I and II have all been conducted for low-wave steepness 
$k{a_w}=0.034$ 
$({a_w}=0.02$ m). In category III we examine the dependence of drift on wave steepness. We select three different relative sizes (
$l/\lambda =5.1\,\%, 8.0\,\%, 10.0\,\%$; see table 8 in Appendix C for all object properties) and perform simulations with steepness in the range 
$k{a_w}=0.02$ to 
$k{a_w}=0.13$. Cases with and without viscosity are considered.
The dimensionless amplitudes of the horizontal oscillatory motion 
${A_x}/{a_w}$ and the vertical oscillatory motion 
${A_z}/{a_w}$ are shown as a function of wave steepness 
$k{a_w}$ in figure 8(a) and 8(b), respectively. The wave celerity-normalized object drift velocities of the objects of different sizes are shown as a function of wave steepness 
$k{a_w}$ in figure 8(c); the object drift is shown as a ratio of the Stokes drift, namely 
$u_{O}/u_{S}$, in figure 8(d), noting 
$u_{S}\sim (k{a_w})^2$. Finally, the local wave amplitude distribution of a RCB with 
$l/\lambda =10\,\%$ is shown as a function of horizontal distance from the centre of mass of the object in figures 8(e) and 8(f) for three values of wave steepness. In (e) the wave amplitude is given in dimensional form as a difference between the local wave amplitude 
$a(x)$ and the input wave amplitude 
${a_w}$. In (f) the local wave amplitude is scaled by its corresponding input wave amplitude, which is different for each steepness.
Figure 8. Effect of wave steepness on object motion and drift (category III): (a) non-dimensional horizontal oscillatory motion amplitude; (b) non-dimensional vertical oscillatory motion amplitude; (c) celerity-normalized object drift velocity; (d) Stokes drift-normalized object drift velocity; (e) difference between the local wave amplitude distribution 
$a(x)$ and the input wave amplitude 
${a_w}$ for three different values of wave steepness for the object size of 
$l/\lambda =10\,\%$; (f) normalized local wave amplitude distribution 
$a(x)/{a_w}$ for three different values of wave steepness for the object size of 
$l/\lambda =10\,\%$. Diamond, star and square markers represent the results for RCBs of 
$l/\lambda =5\,\%$, 
$8\,\%$, 
$10\,\%$ relative sizes, respectively. Red lines in (a,b) denote input wave amplitudes, while in (c,d) they denote theoretical Stokes drift. The local wave amplitude distributions 
$a(x)$ in (e,f) are given as functions of horizontal distance (scaled by wavenumber) from the centre of mass 
${x_c}$, and the gaps between the two vertical red lines in (e) and (f) represent the position of the object (
$l/\lambda =10\,\%$) and correspond to its left and right sides, respectively.
We commence our analysis with the oscillatory motions of the objects. For each box, the horizontal oscillatory motion amplitude, scaled by 
${a_w}$, does not show any obvious variation with steepness (figure 8a), while the vertical motion amplitude, scaled by 
${a_w}$ shows a small decrease as the wave steepness increases (figure 8b). We note that this is consistent with the reduction in the heave response amplitude operator with increasing wave height reported for wave energy devices of similar 2-D shape (Palm et al. Reference Palm, Eskilsson, Bergdahl and Bensow2018). This is attributed therein to the induced drag and nonlinearity of the force the waves exert on the object. For the object drift, figures 8(c) and 8(d) show that an object of relative size 
$l/\lambda =5\,\%$ continues to follow the Stokes drift without notable enhancement (2 %) despite the increase in wave steepness (for 
$\nu =0$). As the object becomes large enough to induce a drift enhancement at low steepness (i.e. 
$l/\lambda =8\,\%$ and 
$10\,\%$; cf. § 4.1), the drift is further enhanced as the waves become steeper. The amplification factors 
$u_{O}/u_S$ of these large boxes initially decrease somewhat with increasing steepness for low values of steepness, but then reach constant values as steepness increases. This is consistent with what has been found in the experiments conducted by Huang et al. (Reference Huang, Law and Huang2011), Huang & Law (Reference Huang and Law2013), He et al. (Reference He, Ren and Qiu2016) and Tanizawa et al. (Reference Tanizawa, Minami and Imoto2002). To be more specific, in the experiments of Huang et al. (Reference Huang, Law and Huang2011) for ‘small’ floating objects (
$l/\lambda =13\,\%\unicode{x2013}16\,\%$), these authors found that object drift is enhanced more as wave steepness increases and that the amplification factor 
$u_{O}/u_{S}$ behaved in a similar fashion as presented here. Due to the large computational cost, we do not increase the relative size of our object beyond 
$10\,\%$.
Focusing on the standing wave pattern, identified in §§ 4.1 and 4.2 as intimately related to drift enhancement, figures 8(e) and 8(f) show the wave amplitude distribution for a round-cornered box with 
$l/\lambda =10\,\%$ for three different values of steepness 
$ka_w$. The absolute magnitude of the standing wave pattern increases with steepness, which is consistent with the increase in object drift shown for this object in figure 8(c). However, the normalized wave amplitude distribution in figure 8(f) shows a modest decrease in the amplitude of the wave pattern for the higher-steepness cases (
$k{a_w}=0.07, 0.09$), which is consistent with the behaviour of the amplification of object drift relative to Stokes drift in figure 8(d).
4.4. Effects of viscosity (categories I and III)
Finally, we examine the role played by viscosity, re-examining the category I and III simulations. We do so by comparing results from our hybrid numerical model qaleFOAM that are based on the inviscid Euler equations (
$\nu =0$) and those that are based on the viscous NS equations (
$\nu =1\times 10^{-6}\ {\rm m}^{2}\ {\rm s}^{-1}$). We define and estimate Reynolds and Keulegan–Carpenter numbers in Appendix D, where we also present simulations using the Reynolds-averaged NS equations to examine the potential role of turbulence. These results show that the flow is laminar in our cases (at laboratory rather than field scale) with low Reynolds numbers and the inclusion of a turbulence model to ensure convergence is not necessary.
We begin by re-examining the simulations in category I (§ 4.1). As shown in figures 5(a) and 5(b), the inclusion of viscosity induces negligible change to the oscillatory motion in the horizontal direction and a small increase in the vertical direction. This is more evident for RCBs. Turning to the object drift velocity (figure 5c), we start by examining ROs because no obvious change to the standing wave pattern arises from the inclusion of viscosity (not shown). In the absence of viscosity no significant drift enhancement is found for ROs of all sizes considered, whereas enhanced drift becomes evident for ROs larger than 
$8\,\%$ when viscosity is considered.
Next, we consider RCBs, for which the standing wave pattern comes into play for large enough relative sizes. When the standing wave pattern is small, which is the case for objects with a relative size smaller than 
$7\,\%$, the presence of viscosity contributes to drift enhancement in a way consistent with the behaviour of ROs. As a RCB becomes larger, the draft (submergence depth) of the box becomes larger and the standing wave pattern starts to drive drift enhancement. Viscosity works to promote enhanced drift and yields a larger drift increase compared with the case without viscosity included. For the largest RCB, we observe a nearly 20 % increase as a result of the inclusion of viscosity. We note the effect of the standing wave pattern and the effect of viscosity appear independent, with viscosity generally not affecting the standing wave pattern (not shown explicitly). From the simulations in category III, we observe from figure 8 that the inclusion of viscosity enhances the drift for all boxes. However, as a ratio of the Stokes drift, the enhanced drift speed reduces with wave steepness for low steepness, then reaches a constant value (figure 8d).
The fact that drift increases with relative size when viscosity is considered (in the form of viscous drag) is consistent with the findings of Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021), who do not consider diffraction of the wave field and who examine three-dimensional (3-D) spheres instead of the 2-D ROs considered here. Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021) propose two mechanisms to explain enhanced drift motion. First, they note that the linear motion (normal to the free surface) of a floating particle has a larger magnitude compared with that of a Lagrangian particle, leading to an increased drift. Second, the dynamic buoyancy force has a net effect when averaged over the wave cycle in a similar fashion to the slope-sliding term of Rumer et al. (Reference Rumer, Crissman and Wake1979). This net effect arises after averaging over the wave cycle because of a phase change that is introduced by the effect of a (viscous) drag force acting in the direction normal to the free surface. A comparison between our results for 2-D ROs and their results for 3-D spheres is made in figures 5(a), 5(b) and 5(c). To evaluate the model of Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021), we have taken the diameter of our (2-D) ROs to be equal to the diameter of the spheres in Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021). Due to the difference in geometry (two dimensional vs three dimensional), we emphasise that we do not expect quantitative agreement. As shown in figures 5(a) and 5(b), Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021) found that the horizontal linear motion remains unaffected by viscous drag, but the magnitude of vertical linear motion increases with relative size. We observe similar results for the vertical linear motion, although the vertical motion we observe is much smaller for the same relative size. Unlike Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021), we predict the horizontal linear is reduced. Figure 5(c) shows a reasonable level of agreement for the drift between our results and Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021) when the object is small, but neither theory predicts significant drift enhancement for such small objects. For larger objects, we observe less enhanced drift than predicted by Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021). The discrepancies in linear motion and drift may be due to the fact that the theoretical model in Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021), based on the slope-sliding concept, does not consider two-way coupling between the waves and the object and assumes that the wave field is unaffected by the presence of the object, thus causing the object in Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021) to lose less energy, as only viscous and no wave damping is considered therein. Wave damping could lead to reduced horizontal and vertical linear motions (Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021) predict larger linear motion responses in both directions), which in turn results in less enhanced drift. Wave damping could also contribute to the phase difference and potentially enhance the second mechanism. Nevertheless, both mechanisms in Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021) could play a role in explaining the phenomenon that the inclusion of viscosity for relatively large objects enhances the drift. For completeness, we note that we cannot rule out the occurrence of some boundary-layer streaming in the near-surface wave-driven boundary layer (see, e.g. Grue & Kolaas Reference Grue and Kolaas2017), which would also enhance drift and also only arises in the presence of viscosity, although this boundary layer only has a very short distance to develop (namely, only in the QALE-FEM domain). We examine the effects of viscosity in more detail in Appendix C.
4.5. Relationship between local standing wave pattern and the object drift
Thus far, we have investigated the effects of various factors on the drift behaviour of finite-size floating objects. All the results indicate that drift enhancement is closely related to the diffraction of the wave field. To gain a more quantitative understanding of this relationship, figure 9 shows the drift speed as a function of the maximum local wave amplitude 
$a_{max}$ obtained from the standing wave pattern. All the results shown in this figure are based on the same input incident wave and, thus, identical theoretical Stokes drift. We note that the ‘wave gauges’ we have used to output information about the free surface elevation from the code are set at fixed intervals, making it challenging to precisely predict the local wave amplitude maxima. According to figure 9, there is a positive correlation between the local maximum amplitude 
$a_{max}$ and the object drift, which is most clearly observable when the local maximum wave amplitude 
$a_{max}$ is relatively large (i.e. 
$a_{max}/a_w \ge 1.05$). Figure 9 shows that the correlation between the maximum local wave amplitude (as a measure of how much diffraction takes place) and object drift is similar, regardless of differences in shape, submergence or size of the object, as long as the object is large enough to diffract the wave field.
Figure 9. The celerity-normalized object drift 
$u_{O}/c$ as a function of the local maximum wave amplitude normalized by the input wave amplitude 
$a_{max}/a_w$: (a) the results for the RCBs and ROs of different sizes (category I simulations); (b) the results for the objects of 
$l/\lambda =10\,\%$ with different round-corner radii and depths of submergence (category II simulations). All results shown are for the same input wave condition: 
$a_w=0.02$ m and 
$u_S/c=0.0012$.
It is instructive to note that an increase in the magnitude of the vertical oscillatory motion is always accompanied by a more distinct standing wave pattern, as particularly evident in the case of large RCBs. This is because the formation of the standing wave pattern results from the disturbance to the fluid field caused by the presence of the object, which would be largest if the object were stationary (in which case we have only diffraction, no radiation). However, when the object is free to move, its motion serves as a response to the waves, mitigating the effects of diffraction. The radiated wave field arises due to the object's oscillation (as if it were in calm water, in the linear approximation), generating a wave field that weakens the standing wave pattern present in the combined diffracted and incident wave field.
5. Comparison between the hybrid numerical model and the diffraction-modified Stokes drift model
To develop the hypothesis developed in § 4 that drift enhancement is a result of diffraction of the wave field by the object and gain a better understanding of the underlying mechanism (in the absence of the viscosity), we compare the predictions of the hybrid numerical model qaleFOAM presented in § 4 with the diffraction-modified Stokes drift model introduced in § 2. In particular, the diffraction-modified Stokes drift model can distinguish the contributions to the object drift of the incident, diffracted and radiated parts of the wave field. The objects considered are the same as for the hybrid numerical model, although we do not include rounded corners for the RB in the diffraction-modified Stokes drift model (we set 
$r/{h_d}=0.24$ in the simulations of the hybrid numerical model shown in this section). The grid sizes of all BEM simulations are chosen to be 
$\Delta x/{h_d}=\Delta z/{h_d}=0.01$ for the RB and 
$\Delta l/{h_d}=0.01$ for the RO (following BEM model verification in Appendix A and a convergence study not shown here). We examine oscillatory motion (§ 5.1) and object drift (§ 5.2) in turn.
5.1. Oscillatory motion
5.1.1. Amplitude of oscillatory motion
Figure 10 provides the non-dimensional amplitudes of oscillatory motion as a function of relative size for the different objects defined in table 3 predicted by both the BEM model and the qaleFOAM model for low-wave steepness 
$k{a_w}=0.034$. Specifically, for the BEM model, these take the form
\begin{equation} {\xi_x}=A_{x}\cos(\theta _{x}-\omega t) \quad \textrm{and}\quad {\xi_z}=A_{z}\cos(\theta_{z}-\omega t). \end{equation}
Figure 10. Comparison of amplitudes of oscillatory motion for objects of different sizes and shapes predicted by the qaleFOAM (QF) and BEM models for 
$k{a_w}=0.034$: (a) horizontal oscillatory motion; (b) vertical oscillatory motion of objects. Motion amplitudes in (a,b) are normalized by input wave amplitude 
${a_w}$, with the red lines then corresponding to the behaviour of a Lagrangian particle. Lines correspond to the results from BEM simulations, while markers are those from qaleFOAM (QF) simulations. Square markers represent results for boxes, RBs for BEM simulations and RCBs for qaleFOAM simulations, and circle markers denote ROs.
Figure 10 demonstrates that the amplitude of motion between the two models is in agreement. The effects of size and shape on the oscillatory motions in the BEM model are consistent with those in the qaleFOAM model discussed in §§ 4.1 and 4.2. The BEM and the qaleFOAM models agree well for the full range of steepnesses studied in this paper (
$k a_w=0.02\unicode{x2013}0.11$) (not shown explicitly).
5.1.2. Phase of oscillatory motion
Having demonstrated the ability of the BEM model to capture the amplitude of the oscillatory motion, we now examine its phase. The phases of the vertical and horizontal oscillatory motions predicted by the (linear) BEM model are defined in (5.1a,b), and we compare these to the phases of the linear incident wave field. For a linear incident wave of the form 
$\eta =a_w \cos (kx-\omega t+\theta _w)$, (linearized) horizontal and vertical components of the motion of a Lagrangian particle have the form
$$\begin{gather} \xi_{I,x}(t)=a_{w}\sin(\omega t-kx-\theta_{w})_{x=0}=a_{w}\cos(\theta_{w}+{\rm \pi} /2-\omega t), \end{gather}$$
$$\begin{gather}\xi_{I,z}=a_{w}\cos(\omega t-kx-\theta_{w})_{x=0}=a_{w}\cos(\theta_{w}-\omega t), \end{gather}$$
which have been evaluated at the particle's initial location in the BEM model (
$x=0$). Figure 11 shows the phase difference between a finite-size object and a Lagrangian particle as a function of relative size of the object for both shapes for the horizontal (i.e. 
$\theta _{x}-\theta _{w}-{\rm \pi} /2$) and vertical (i.e. 
$\theta _{z}-\theta _{w}$) oscillatory motions predicted by the BEM model.
Figure 11. Phase difference of the oscillatory motion between a finite-size object and a Lagrangian particle as a function of relative object size from BEM simulations: (a) horizontal direction, (b) vertical direction. Blue lines correspond to RBs and black lines to ROs.
Figure 11 shows that when objects are very small, the phase difference is zero, confirming that the objects behave as Lagrangian particles. As the objects become larger, the magnitude of the phase differences of both horizontal and vertical motions become larger, and this relationship is nonlinear. Specifically, the larger the object is, the more phase lag it gains vertically while the more phase lead it shows horizontally. The magnitudes of the phase difference of the vertical motion are much larger than that of the horizontal motion, which are negligibly small.
We note that RBs have larger phase lags vertically but smaller phase leads horizontally compared with ROs. Furthermore, for objects with greater submergence depths, we find greater phase lags in the vertical and smaller phase leads in the horizontal. The phase differences are found to be independent of steepness. These results are not shown here explicitly in the interest of brevity.
5.2. Object drift
We now compare predictions of object drift by the qaleFOAM model already examined in § 4 to predictions of object drift based on the diffraction-modified Stokes drift model (i.e. using (2.8)), i.e. an estimate of the drift accounting for the radiated and diffracted as well as the incident waves. We consider waves with a low input steepness of 
$k{a_w}=0.034$, and the dimensions of the RBs and the ROs we consider are given in table 3 (and table 8 for boxes larger than 
$l/\lambda =4\,\%$). Figure 12 makes the comparison between drift predictions by the qaleFOAM and the diffraction-modified Stokes drift model.
Figure 12. Diffraction-modified Stokes drift velocity 
$u_{S,O}$ as a function of object size: (a) comparison between BEM and qaleFOAM (QF) models with (
$\nu =1.00\times 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$) and without (
$\nu =0$) viscosity for ROs and rectangular boxes (RBs, RCBs); (b) different components of the diffraction-modified Stokes drift velocity predicted by the diffraction-modified Stokes drift model, where 
$u_{S,O,I}$, 
$u_{S,O,R}$ and 
$u_{S,O,D}$, respectively, represent the incident, radiation and diffraction components of the diffraction-modified Stokes drift velocity.
It is clear from figure 12 that the diffraction-modified Stokes drift model captures the main trend well, predicting a significant increase of object drift with increasing relative size for RBs. The diffraction-modified Stokes drift model and the qaleFOAM model without viscosity (
$\nu =0$) show good agreement except for the largest box with 
$l/\lambda =10\,\%$, for which the diffraction-modified Stokes drift model underestimates the drift velocity compared with the prediction of the qaleFOAM model (for 
$\nu =0$). This discrepancy can be attributed to two factors: the slight differences in shape between the boxes in the diffraction-modified Stokes drift model (RB with sharp corners) and in the qaleFOAM model (RB with round corners) and the approximate nature of the diffraction-modified Stokes drift model (cf. (2.3)), whose underlying assumptions become less valid as the object becomes larger.
To analyse the physical mechanism underlying the drift enhancement, we decompose the object drift predicted by the diffraction-modified Stokes drift model (i.e. (2.8)) into contributions from the incident, radiated and diffracted waves, i.e.
\begin{equation} u_{S,O}=u_{S,O,I}+u_{S,O,R}+u_{S,O,D}, \end{equation}where
$$\begin{gather} u_{S,O,I}=\frac{{{A_x}}}{2}k \omega {A_{xx}^I}\cos ({\theta _x}- {\theta _{xx}^I})+\frac{{{A_z}}}{2}k \omega{A_{xz}^I}\cos ({\theta_z} - {\theta _{xz}^I}); \end{gather}$$
$$\begin{gather}u_{S,O,R}= \frac{{{A_x}}}{2}k \omega {A_{xx}^R}\cos ({\theta _x}- {\theta _{xx}^R})+\frac{{{A_z}}}{2}k \omega {A_{xz}^R}\cos ({\theta_z} - {\theta _{xz}^R}); \end{gather}$$
$$\begin{gather}u_{S,O,D}= \frac{{{A_x}}}{2}k \omega {A_{xx}^D}\cos ({\theta _x}- {\theta _{xx}^D})+\frac{{{A_z}}}{2}k \omega {A_{xz}^D}\cos ({\theta_z} - {\theta _{xz}^D}). \end{gather}$$
Each term (i.e. (5.5), (5.6) and (5.7)) includes contributions from the horizontal (first term) and vertical (second term) motion of the object. Figures 12(b) and 13 show these contributions as a function of relative object size. We start by examining the contribution of the incident wave field to the diffraction-modified Stokes drift, 
$u_{S,O,I}$. For very small objects, we have 
${A_x}={a_w}$, 
${\theta _x}-\theta _w={\rm \pi} /2$, 
${A_z}={a_w}$ and 
${\theta _z}-\theta _w=0$ (cf. figures 10 and 11). Combined with 
${A_{xx}^I}={a _w}gk/\omega ^2={a _w}$, 
${\theta _{I,x}}-\theta _w={\rm \pi} /2$, 
${A_{xz}^I}={a _w}gk/\omega ^2={a _w}$ and 
${\theta _{I,z}}-\theta _w=0$ (from (2.5)) for the deep-water incident wave field we consider, it is readily evident from (5.5) that
\begin{equation} u_{S,O,I}=\frac{{{a_w}}}{2} k \omega {a_w} \cos ({\rm \pi}/2- {\rm \pi}/2)+\frac{{{a_w}}}{2}k \omega {a_w} \cos (0 - 0)={a_w}^2\omega k=u_{S}. \end{equation}Any contributions to the diffraction-modified Stokes drift from the radiated and the diffracted wave field for very small objects must therefore be equal and opposite (cf. (5.4)), which is consistent with what is presented in figure 12(b).
Figure 13. Decomposition of the diffraction-modified Stokes drift predicted by the diffraction-modified Stokes drift model for RBs according to (5.4) and (5.5): (a) contributions of the incident (
$u_{S,O,I}$) and the sum of the radiated and diffracted waves (
$u_{S,O,R}+u_{S,O,D}$); and (b) further decomposition into contributions due to horizontal (
$x$) and vertical (
$z$) components of object motion. The diffraction-modified Stokes drift velocity components are normalized by the Stokes drift (of an infinitesimally small object).
As the object becomes larger, the amplitude of horizontal oscillatory motion 
${A_x}$ becomes smaller (figure 10a), its phase difference does not change significantly (figure 11a), while the amplitude of vertical oscillatory motion 
${A_z}$ becomes larger (figure 10b), but the phase difference of the vertical motion increases (figure 11a), diminishing the effect of the enhanced amplitude of vertical motion (cf. 
$A_z\cos ({\theta _{I,z}}-\theta _w)$).
A careful reader may observe that the incident component of drift 
$u_{S,O,I}$ in figures 12(b) and 13(a) experiences a slight decrease before undergoing a more significant increase. To explain this, we note that as object size increases, the amplitude of vertical oscillatory motion increases (cf. figure 10b) but its positive effect on drift is diminished by the increasing phase difference (cf. figure 11b), while the amplitude of horizontal oscillatory motion decreases (cf. figure 10a), acting to reduce drift. When the negative effects resulting from reduced horizontal oscillatory motion compete over the positive contribution of the enhanced vertical motion, the incident drift component 
$u_{S,O,I}$ is reduced. This is evident in figures 12(b) and 13(a) for objects with a relative size between 
$3\,\%$ and 
$7\,\%$.
Noting that 
${A_x}$ decreases with relative size in a linear fashion (figure 10a), whereas 
${A_z}$ increases nonlinearly at an increasing rate (figure 10b), figures 13(a) and 13(b) show that the effect of the increased vertical oscillatory motion generally dominates and the contribution of the incident waves acts to increase the object drift for larger objects. However, the increased vertical oscillatory motion cannot explain the majority of the large drift enhancement observed for large objects.
Turning to the contributions of the radiated and diffracted waves to the diffraction-modified Stokes drift, 
$u_{S,O,R}$ and 
$u_{S,O,D}$, we observe from figure 12(b) that both terms decrease rapidly with increasing relative size. Since these two terms have opposing signs, the fact that 
$u_{S,O,D}$ decreases more rapidly leads to a net positive contribution to the diffraction-modified Stokes drift that increases with relative size, as illustrated in figure 13(a). This explains most of the enhanced drift for large objects found in this paper.
However, if the small reduction of the incident drift component 
$u_{S,O,I}$ for objects with intermediate size, described above, cannot be compensated for by the net positive contribution from the imbalance of diffraction and radiation components, the overall drift will be reduced. This helps to explain the slight reduction in drift compared with the theoretical Stokes drift for objects with a relative size of 
$2\,\% \leq l/\lambda \leq 6\,\%$ in figure 12(a).
From figure 13(b) it is further evident that the horizontal object motion is responsible for the increase in 
$u_{S,O,R}$ and 
$u_{S,O,D}$ and, thus, the total diffraction-modified Stokes drift for large objects. Note that the motions evaluated in § 4 are linear oscillatory motions, while the object motion evaluated in figure 13 is the drift motion. According to (5.5)–(5.7), the drift components are dependent not only on oscillatory motions (amplitudes and phases) but also on derivatives of the velocity field. Only their combined effect determines the drift. Physically, such large objects are less able to follow the horizontal motion the waves would induce for a Lagrangian particle (cf. figure 10a) and are therefore transported at a faster speed.
It may seem somewhat counter-intuitive that the smallest object has the largest diffraction/radiation drift components (i.e. 
$u_{S,O,R}$ and 
$u_{S,O,D}$) and that the absolute values of these components decrease as the object size increases, as shown in figure 12(b). However, what really matters here is the combined contribution of the diffraction and radiation components, as shown in figure 13(a), as they do not contribute to drift independently.
For ROs, the decrease in 
$u_{S,O,I}$ is not balanced by a sufficiently large increase in the sum of 
$u_{S,O,R}$ and 
$u_{S,O,D}$, and these objects do not experience an increase in diffraction-modified Stokes drift for large size (in the absence of viscosity). For increased submergence depth, the effects discussed above for a RB are only enhanced. Furthermore, the diffraction-modified Stokes drift model (cf. (2.3)) and the qaleFOAM model agree well for the full range of steepnesses studied in this paper (
$k a_w=0.02\unicode{x2013}0.11$). These results are not shown here explicitly in the interest of brevity.
6. Conclusions
In this paper we have investigated the fluid mechanics that can lead to enhancements in the drift of floating objects under the influence of gravity waves beyond that of the well-understood Stokes drift. We restrict our analysis to unidirectional waves on deep water and where the object is less than 
$10\,\%$ the size of the wavelength. Based on numerical modelling we have identified two mechanisms that can explain increased drift compared with the Stokes drift: a mechanism that relies on viscosity and a mechanism that is related to the diffraction of the wave field by the object and the standing wave pattern that arises. Both these mechanisms come into play when the size of the object is larger than a few percent of the wavelength. When the object is smaller than this, the inertial (i.e. non-Lagrangian) behaviour of the object becomes less evident and the difference in velocity between the object and the surrounding fluid can be ignored. There is no obvious increase in the amplitude of the vertical motion, and the drift motion becomes that of a Lagrangian particle. As the object becomes larger, the amplitude of the motion normal to the free surface increases, which creates a drag force because of the difference between the object motion and the fluid surrounding it. This effect can cause enhanced horizontal drift (Calvert et al. Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021). In addition, as the object size becomes larger, the draft of the object (submergence depth) becomes larger, and the submerged part of the object acts to impede the fluid motion associated with the waves and thereby changes the waves themselves. That is, the object diffracts the wave field. Dependent on how large the submerged part of the object is and on its shape, the impeding effect is different, and thus, the drift is enhanced to a different degree. The larger the submergence depth is and the closer the shape of the object is to a box (i.e. not streamlined), the stronger is the impeding effect, yielding a larger increase in horizontal drift.
We consider objects of up to 
$10\,\%$ of the wavelength and, for the largest of these drift enhancements over that of a Lagrangian particle, can be as large 
$92\,\%$ of the Stokes drift for simulations without viscosity or 
$113\,\%$ with viscosity. Most of the increased drift results from diffraction for RCBs with viscosity typically contributing a further 20 %. For ROs, diffraction only has a small effect, and the much smaller enhanced drift arises because of the effects of viscosity.
To enable quantitative predictions to be made about the contribution of diffraction to drift, we have derived a diffraction-modified Stokes drift akin to Stokes (Reference Stokes1847), but accounting for the combination of the incident, diffracted and radiated wave fields rather than simply the first of these. To calculate the necessary diffracted and radiated fields, a linear BEM model based on potential-flow theory is used. Two effects become clear. First, the increased vertical oscillatory motion of the object causes a greater contribution from the incident wave field to the diffraction-modified Stokes drift. Second, the combination of diffracted and radiated waves makes an additional contribution to the diffraction-modified Stokes drift that is not present when the object is small. Although we have not analysed the force and momentum balance resulting in the object's (enhanced) drift motion, we foresee this will give valuable insight into the mechanism and recommend it as future work.
Various authors have found evidence for enhanced drift in different circumstances. Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021) (based on previous work by Rumer et al. (Reference Rumer, Crissman and Wake1979), Shen & Ackley (Reference Shen and Ackley1991) and Huang et al. (Reference Huang, Huang and Law2016)) explored the influence of viscosity described above using a theoretical model that considers viscous forces but ignores the diffraction of the wave field caused by the presence of the object (required for the second mechanism). Future work should quantitatively compare the findings of the present work on the effect of viscosity to the predictions of Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021), taking account of Reynolds number and wave steepness, but most importantly for the same geometry, that is, by extending the 3-D model of Calvert et al. (Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021) to a 2-D model or our 2-D numerical simulations to 3-D numerical simulations. According to the experiments and theoretical analysis by Longuet-Higgins (Reference Longuet-Higgins1953, Reference Longuet-Higgins1960), Grue & Kolaas (Reference Grue and Kolaas2017), the presence of viscosity should also be accompanied by an additional mechanism for (Eulerian-mean) wave-induced drift, namely boundary-layer streaming. In principle, boundary-layer streaming of the free surface is included in the viscous simulations performed in this paper but it is not explicitly investigated and likely only small, as the boundary layer only has a short distance to develop in the NS part of the hybrid numerical model. As boundary-layer streaming is associated with strong vertical shear, its differential effect on objects of different submergence depths is foreseeable and should be investigated for inertial objects in future work.
Enhanced object drift due to diffraction has probably been observed in previous experiments, although it has not been identified as such. Harms (Reference Harms1987) showed using experiments that, for ice floes (box shaped) with large submergence depths, the drift velocity increased with the draft of the object, keeping size constant (for relative sizes smaller than approximately 
$25\,\%$). Enhanced drift was also found for larger submergence depths in experiments conducted by Huang et al. (Reference Huang, Law and Huang2011), and a similar trend of the drift scaled by Stokes drift as a function of wave steepness can be identified in their results. In the experiments carried out by He et al. (Reference He, Ren and Qiu2016) for regular waves in finite depth, drift enhancement is seen to increase with wave steepness for boxes with 
$l/\lambda =9\,\%$ and 
$10\,\%$. Future experiments should focus explicitly on identifying the standing wave pattern associated with diffraction and should explore the roles of both length and width (i.e. 3-D effects) of objects relative to the wavelength.
Funding
Q.X. is supported by the China Scholarship Council-PAG Oxford Scholarship. R.C. was supported by funding from the European Space Agency (grant no. 4000136626/21/NL/GLC/my). T.S.vdB was supported by a Royal Academy of Engineering Research Fellowship.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Mathematical formulation and verification of the BEM model
In the BEM model we deal with the problem based on linearized potential-flow theory, which assumes the fluid is homogeneous, incompressible, inviscid and the flow is irrotational. The flow is described by a velocity potential 
$\varPhi$, which, under the linear approximation, can be further divided into an incident potential 
${\varPhi _I}$, a diffraction potential 
${\varPhi _D}$ and a radiation potential 
${\varPhi _R}$ (cf. (2.1)), where the incident potential neglects the influence caused by the presence of the object, the diffraction potential considers the disturbance induced by the presence of the object on the incident wave field and ignores the oscillatory motion of the object, and the radiation potential describes the wave field forced by the oscillatory motion of the object in still water. The incident potential is given in (2.2).
A.1. Radiation potential
Based on potential-flow theory (e.g. Newman (Reference Newman2018), specifically, §§ 6.15–6.19), the radiation potential 
${\phi _R}$ can be further decomposed into components representing radiation of waves due to horizontal, vertical and rolling object motions, i.e.
\begin{equation} {\phi _R}(x,z) = \sum_{j=1}^{3} - \iota \omega {A _j}{\hat{\phi} _j}(x,z), \end{equation}
in which 
${A_j}$ (for 
$j=x,z,\alpha$) represent the amplitudes of the corresponding three modes of motion of the object, and 
$\widehat {\phi _j}$ (
$\kern 1.5pt j=x,z,\alpha$) represents the radiation velocity potential due to unit body motion in mode 
$j$. Three motion modes are included: horizontal (
$\kern 1.5pt j=x$), vertical (
$\kern 1.5pt j=z$) and rolling motion (
$\kern 1.5pt j=\alpha$, where 
$\alpha$ is defined as the angle of rotation around the 
$y$ axis in the (
$x,z$) plane). The potential 
${\widehat {\phi _j}}$ is governed by Laplace's equation [
$L$] and subject to a linearized free surface boundary 
${C_F}$ condition [
$F$], body surface 
${C_B}$ condition [
$B$], bottom 
${C_D}$ condition [
$D$] as well as a far-field radiation 
${C_R}$ condition [
$R$], shown as
\begin{equation} \left.\begin{array}{lll} [L] & \dfrac{{{\partial ^{\rm{2}}}{\widehat{\phi _j}}}}{{\partial {x^2}}} + \dfrac{{{\partial ^{{2}}}{\widehat{\phi _j}}}}{{\partial {z^2}}} =0, & \\[10pt] [F] & \dfrac{{\partial {\widehat{\phi _j}}}}{{\partial z}} - \dfrac{{{\omega ^2}}}{g}{\widehat{\phi _j}} = 0 & \textrm{for } z=0,{C_F}, \\[10pt] [B] & \dfrac{{\partial {\widehat{\phi _j}}}}{{\partial N}} = {N_j} & \textrm{for }{C_B},\\[10pt] [D] & \dfrac{{\partial {\widehat{\phi _j}}}}{{\partial z}} = 0 & \textrm{for } z=-d, {C_D},\\[10pt] [R] & \mathop{\lim}\limits_{x \to \pm \infty } \left( {\dfrac{{\partial {\widehat{\phi _j}}}}{{\partial x}} \mp {\rm i}k{\widehat{\phi _j}}} \right) = 0 & \textrm{for } {C_R}, \end{array}\right\} \end{equation}
in which 
$C_B$ refers to boundaries of 
$C_1$, 
$C_2$ and 
$C_3$ in figure 3(a) and 
$C_r$ in figure 3(b), 
${N}$ is the unit normal vector of the object surface and 
${N_j}$ represents the projection of the unit normal vector of the relevant boundary in the 
$j$th direction.
In numerical simulations, we truncate the domain at 
$x=\pm L_{BEM}/2$ for both horizontal sides (shown in figure 3) and rewrite the radiation condition in a uniform expression for both ends as
\begin{equation} \frac{{\partial {\widehat{\phi _j}}}}{{\partial N}} = {\rm i}k{\widehat{\phi _j}}. \end{equation}A.2. Diffraction potential
The governing equation and boundary conditions of the diffraction potential 
${\phi _D}$ are
\begin{equation} \left.\begin{array}{lll} [L] & \dfrac{{{\partial ^{\rm{2}}}{\phi _D}}}{{\partial {x^2}}} + \dfrac{{{\partial ^{\rm{2}}}{\phi _D}}}{{\partial {z^2}}} =0, & \\[10pt] [F] & \dfrac{{\partial {\phi _D}}}{{\partial z}} - \dfrac{{{\omega ^2}}}{g}{\phi _D} = 0 & \textrm{for } z=0,{C_F}, \\[10pt] [B] & \dfrac{{\partial {\phi _D}}}{{\partial N}} =-\dfrac{{\partial {\phi _I}}}{{\partial N}} & \textrm{for }{C_B},\\[10pt] [D] & \dfrac{{\partial {\phi _D}}}{{\partial z}} = 0 & \textrm{for } z=-d, {C_D,}\\[10pt] [R] & \dfrac{{\partial {\phi _D}}}{{\partial N}} = {\rm i}k{\phi _D} & \textrm{for }{C_R}. \end{array}\right\} \end{equation}From (A4) and (A2), it can be seen that the description of the diffraction potential is very similar to that of the radiation potential; the only difference is in the body surface condition.
A.3. Equations of motion
The equations of motion of the object are given by
\begin{equation} \boldsymbol{Ma}=\boldsymbol{F}_{R}+\boldsymbol{F}_{C}+\boldsymbol{F}_{K}+ \boldsymbol{F}_{D}, \end{equation}
in which 
$\boldsymbol {M}$ is the general mass matrix, 
$\boldsymbol {a}$ is the acceleration, 
$\boldsymbol {F}_{R}$, 
$\boldsymbol {F}_{C}$, 
$\boldsymbol {F}_{K}$, 
$\boldsymbol {F}_{D}$ are the radiation, restoring, incident and diffraction forces, respectively. The radiation forces 
$F_{R,kj}$ can be expressed as
\begin{equation} {F_{R,kj}} =- {{\rm Re}}\{ ( - {\omega ^2}{\mu _{kj}} - \iota \omega {\lambda _{kj}}){A_j}{{\rm e}^{ - \iota \omega t}}\}, \end{equation}
where 
${\lambda _{kj}}$, 
${\mu _{kj}}$ are the added mass coefficient and damping coefficient, respectively. They can be calculated by
$$\begin{gather} {\mu _{kj}} = {{\rm Re}} \left( {\rho_f \int\limits_{{C_B}} {{\widehat{\phi _j}}{N_k}\,{\rm d}l} } \right), \end{gather}$$
$$\begin{gather}{\lambda _{kj}} = \omega \,{{\rm Im}} \left( {\rho_f \int\limits_{{C_B}} {{\widehat{\phi _j}}{N_k}\,{\rm d}l} } \right), \end{gather}$$
where Re and Im represent the real and imaginary parts of the complex number, and 
$\rho _f$ is the fluid density, 
${N_k}$ represents the projection of the unit normal vector of the boundary in the 
$k$th direction. The restoring forces 
$F_{C,k}$ can be calculated by
\begin{equation} {F_{C,k}} = {{\rm Re}} \{ - {C_{kj}} \cdot {A_j}{{\rm e}^{ - \iota \omega t}}\}, \end{equation}
where 
${C_{kj}}$ is the matrix of restoring force coefficient. Einstein notation is used here to imply the summation over a set of 
$j=1,2,3$. The incident wave forces (Froude–Krylov) 
${F_k^I}$ can be calculated by
\begin{equation} F_k^I = {{\rm Re}} \{\,{f_k^I {{\rm e}^{ - \iota \omega t}}} \}={{\rm Re}} \left\{\rho_f \iota \omega \int\limits_{{C_B}} {\phi _I}{N_k}\,{\rm d}l {{\rm e}^{ - \iota \omega t}} \right\} . \end{equation}
The diffraction wave forces 
${F_k^D}$ of the 
$k$th mode in the BEM model can be calculated in two ways:
$$\begin{gather} {F_k^D}={{\rm Re}} \{\,{f_k^D {{\rm e}^{ - \iota \omega t}}} \}={{\rm Re}} \left\{ {\rho_f} \iota\omega \int\limits_{{C_B}} {{\phi _D}{N_k}\,{\rm d}l} {{\rm e}^{ - \iota \omega t}} \right\}, \end{gather}$$
$$\begin{gather}{F^D_k}={{\rm Re}} \{\,{f_k^D {{\rm e}^{ - \iota \omega t}}} \}={{\rm Re}} \left\{ -{\rho_f} \iota\omega \int\limits_{{C_B}} {\widehat{\phi _k}\frac{{\partial {\phi _I}}}{{\partial N}}\,{\rm d}l} {{\rm e}^{ - \iota \omega t}} \right\}. \end{gather}$$
Here (A11) calculates the diffraction force by directly integrating the diffraction potential, while (A12) calculates the diffraction force using the Haskind formula. If we substitute (A6), (A9), (A10), (A11) into (A5) and take the time factor 
$\textrm {e}^{-\iota \omega t}$ out, we obtain the equations of motions in the frequency domain for the 
$k$th motion mode (
$k=1,2,3$):
\begin{equation} \sum_{j = 1}^3 {[ { - {\omega ^2}({{m_{kj}} + {\mu _{kj}}})-\iota \omega {\lambda _{kj}} + {C_{kj}}} ]} {\kern 1pt} {\kern 1pt} {A_j} = f_k^{I} + f_k^d. \end{equation}The equations of motion for the object contain the (added) mass, hydrodynamic damping and restoring forces on its left side and forces due to the incident and diffracted wave field on the right side. We thus take the object's inertia and wave–body interaction into account.
A.4. Two-dimensional Green's function method
The potentials 
$\widehat {\phi _j}$ and 
$\phi _D$ are harmonic functions and are governed by the Laplace equation. Assuming 
$P(x,z)$ is a field point in the fluid domain and 
$Q(x', z')$ is the source point in the field, we choose a Green's function that satisfies 
$\nabla ^2G(P,Q)=\delta (P,Q)$, then according to Green's second formula, the value of the potential 
$\widehat {\phi _j}$ and 
$\phi _D$ can be determined uniquely by giving its value and normal derivative over all boundaries (Newman Reference Newman2018). We have
$$\begin{gather} \alpha (P){\widehat{\phi _j}}(P) =- \int\limits_C \left[ {{\widehat{\phi _j}}(Q)\frac{{\partial G(P,Q)}}{{\partial N}} - G\frac{{\partial {\widehat{\phi _j}}(Q)}}{{\partial N}}} \right] {\rm d}l, \end{gather}$$
$$\begin{gather}\alpha (P){\phi _D}(P) =- \int\limits_C \left[ {{\phi _D}(Q)\frac{{\partial G(P,Q)}}{{\partial N}} - G\frac{{\partial {\phi _D}(Q)}}{{\partial N}}} \right] {\rm d}l, \end{gather}$$
in which 
$C$ represents the boundaries of the fluid domain, including free surface boundaries 
${C_F}$, body surface boundaries 
${C_B}$, far-field radiation boundaries 
${C_R}$ and bottom boundaries 
${C_D}$. For a smooth boundary, 
$\alpha (P)$ is valued as
\begin{equation} \alpha (P)=\left\{\begin{array}{ll} 2{\rm \pi}, & P \in\varOmega, \\ {\rm \pi}, & P \in S, \\ 0, & P \notin S\cup \varOmega, \end{array}\right. \end{equation}
where 
$\varOmega$ represents the region inside of the fluid domain and 
$S$ denotes the boundary of the fluid domain. Here, the simple Green's function 
$G(P,Q)=\ln ({1}/{r(P,Q)})$ is used. As 
$\ln ({1}/{r(P,Q)})$ is the general solution of the governing equation and does not satisfy any boundary conditions, this requires the source to be distributed over all boundaries.
A.5. Second derivatives of the velocity potential
Calculation of the diffraction-modified Stokes drift based on (2.8) requires the evaluation of second derivatives of velocity potentials in (2.6) and (2.7). Due to the singularity of the Green's function method employed here, direct numerical evaluation of these second derivatives based on finite differences is challenging as it may cause a loss of accuracy (Zhang, Bandyk & Beck Reference Zhang, Bandyk and Beck2010; Chen et al. Reference Chen, Zhu, Zhao, Zhou and Fan2018). We follow Zhang et al. (Reference Zhang, Bandyk and Beck2010) and use the so-called desingularized source distribution method. Different from the standard source distribution method (Newman Reference Newman2018), the desingularized method enforces the boundary conditions that are satisfied exactly on the boundary (denoted by 
$P$) but distributes the source points 
$Q$ slightly outside the boundary, so that the singularities on the boundaries are removed (see also Raven (Reference Raven1988) and Kim & Kim (Reference Kim and Kim2007) for more details). We set the distance of the source points to the boundaries to be 1–2 times the grid size in the direction normal to the boundary. Second derivatives at point 
$P$ are thus calculated as
\begin{equation} \frac{\partial ^2\phi(P)}{\partial x^2} =\int_{{C_l}'}\sigma (Q)\frac{\partial ^2G(P,Q)}{\partial x^2}{\rm d}l; \quad \frac{\partial ^2\phi(P)}{\partial x \partial z} =\int_{{C_l}'}\sigma (Q)\frac{\partial ^2G(P,Q)}{\partial x \partial z}{\rm d}l, \end{equation}
in which 
${C_l}'$ denotes the boundaries where the source points 
$Q$ are located, obtained by moving a certain distance from the original fluid boundaries (i.e. 
${C_1}$, 
${C_2}$ and 
${C_3}$ in figure 3a) according to the desingularized source distribution method. In (A17), 
$\sigma (Q)$ is the source strength at source point 
$Q$. Based on (A17), we can obtain the second derivatives in (2.8) once the source strength 
$\sigma (Q)$ is known. The source strength 
$\sigma (Q)$ is solved following the standard source distribution method. Derivatives of the incident potential are directly calculated from (2.2).
It is worth noting that, theoretically, based on (A17), the diffraction-modified Stokes drift could be evaluated on the surface of the body or at the object's centre of mass. For large objects, we use the second derivatives obtained by (A17) evaluated at the object's centre of mass. However, when the size of the object is very small, the second derivatives evaluated at the centre of mass or at the boundaries of the object become very sensitive to small changes in position in the direction normal to the boundary due to the use of the desingularized source distribution method. Instead, we make use of the boundary conditions on the object boundary (
$C_1$, 
$C_2$, 
$C_3$ or 
$C_r$) to calculate second derivatives.
To improve the accuracy of second derivatives evaluated on object boundaries in the BEM model (i.e. in (2.6) and (2.7)) when the object is small, we take advantage of the boundary conditions, which are themselves given in the form of normal derivatives on the boundary. For the RB defined in figure 3(a), there are three object boundaries: 
${C_1}$, 
${C_2}$ and 
${C_3}$. We will examine the general principle of our method using 
${C_2}$ as an example. The boundary conditions on 
${C_2}$ for the radiation and the diffraction problem are
$$\begin{gather} \displaystyle \frac{{\partial \hat{\phi }_j}}{{\partial z}} = {N_j} = [0,1, - x]\quad \text{for}\ j=x,z,\alpha \quad \text{on}\ {C_2}; \end{gather}$$
$$\begin{gather}\displaystyle\frac{{\partial {\phi _D}}}{{\partial z}} =-\frac{\partial {\phi_I}}{\partial z} \quad \text{on}\ {C_2}, \end{gather}$$
where 
${N_j}$ represents the projection of the unit normal vector of the relevant boundary in the 
$j$th direction. For the object boundary 
${C_2}$ (see figure 3a), the normal vector is in the vertical direction and the normal derivative of its velocity potential in the form of a Green's function is continuous in the 
$x$ direction but not continuous in the 
$z$ direction. We can therefore evaluate horizontal derivatives directly along the boundary and we have
$$\begin{gather} \frac{{{\partial ^2}\widehat{\phi _j}}}{{\partial x\partial z}} = [0,0, - 1]; \end{gather}$$
$$\begin{gather}\frac{{{\partial ^2}{\phi _D}}}{{\partial x\partial z}} =- \frac{{{\partial ^2}{\phi _I}}}{{\partial x\partial z}}. \end{gather}$$
The second derivative 
${{\partial ^2}{\phi _R}}/{\partial x\partial z}$ in (2.6) can now be calculated directly based on (A20) and (A1). The second derivatives 
${{\partial ^2}{\phi _D}}/{\partial x\partial z}$ in (2.7) can be calculated directly based on (A21). The second derivatives 
${{\partial ^2}\widehat {\phi _j}}/{\partial x^2}$, 
${{\partial ^2}{\phi _D}}/{\partial x^2}$ can be calculated numerically directly from the potential as the latter is continuous over the boundary 
${C_2}$ in the horizontal direction. Finally, to obtain a single value to be used to estimate the diffraction-modified Stokes drift, we evaluate the second derivative at the centre of the boundary 
${C_2}$.
A.6. Verification of the BEM model
To verify the BEM model we use in this paper, we evaluate the radiation and diffraction solutions for three specific examples involving rectangular objects in regular waves and compare these numerical solutions to the theoretical solutions of Zheng, You & Shen (Reference Zheng, You and Shen2004). In their theoretical model, the added mass coefficient 
$\mu _{kj}$ and radiation damping 
$\lambda _{kj}$ are calculated based on (A7) and (A8) based on an analytical solution for 
$\widehat {\phi _j}$. The wave excitation forces in their paper are
\begin{equation} {f_k}=f_k^I+f_k^D={\rho_f} \iota\omega \int\limits_{{C_B}} {{(\phi _I+\phi _D)}{N_j}\,{\rm d}l}, \end{equation} In example 1 the object's size and density are chosen so that 
$d/{h_d}=3$ and 
$l/{h_d}=1$. The (truncation) length of the domain 
$L_{BEM}/2=10 h_d$, and the grid size is chosen to be 
$\Delta x/{h_d}=\Delta z/{h_d}=0.01$. Figure 14 compares the normalized added mass and hydrodynamic damping coefficients 
$\mu$ and 
$\lambda$ predicted by our BEM model to their theoretical counterparts by Zheng et al. (Reference Zheng, You and Shen2004). Good agreement is achieved for both added mass and hydrodynamic damping coefficients for a broad range of water depths 
$kd$, including the deep-water values we examine in the paper.
Figure 14. Non-dimensional added mass (
$\mu _{ii}$) and damping (
$\lambda _{ii}$) coefficients for example 1 (
$d/{h_d}=3$, 
$l/{h_d}=1$) of the BEM model verification. The black squares are the predictions by the BEM model, the red lines correspond to the theoretical solutions of Zheng et al. (Reference Zheng, You and Shen2004) based on (A7) and (A8), 
${\rho _0}$ is the density of water and 
$d$ the water depth.
For example 2 and 3, we consider objects with 
$d/{h_d}=2$, 
$l/{h_d}=2$ and 
$d/{h_d}=2$, 
$l/{h_d}=6$, respectively, and we compare the wave-induced forces predicted by our BEM model to their theoretical counterparts by Zheng et al. (Reference Zheng, You and Shen2004). We choose the (truncation) domain length to be 
$L_{BEM}/2=15 l$ and 
$\Delta x/{h_d}=\Delta z/{h_d}=0.01$ for both cases. The diffraction wave forces 
${F_j^D}$ of the 
$j$th mode in the BEM model can be calculated in two ways based on (A11) and (A12). Given the accuracy with which our BEM model solves the radiation problem, as verified in example 1 (figure 14), consistency between the two approaches (i.e. (A11) and (A12)) confirms the diffraction potential is solved correctly. The results of this comparison and the comparison to the theoretical solutions of Zheng et al. (Reference Zheng, You and Shen2004) are given in figure 15 for examples 2 and 3. The BEM model performs well in predicting the wave forces for a range of water depths 
$kd$, and the two different approaches agree well for both examples, further verifying the model.
Figure 15. Non-dimensional wave-induced forces for example 2 (
$d/{h_d}=2$, 
$l/{h_d}=2$) and example 3 (
$d/{h_d}=2$, 
$l/{h_d}=6$) of the BEM model verification. The black and blue squares are the predictions by the BEM model using (A11) and (A12), respectively. The red lines correspond to the theoretical solutions of Zheng et al. (Reference Zheng, You and Shen2004) based on (A22).
Appendix B. Convergence of the hybrid numerical model (qaleFOAM)
Our convergence tests focus on the NS domain, as corresponding tests for the QALE-FEM domain used to simulate the incident wave field have been performed extensively and are well documented in the literature (e.g. Ma & Yan Reference Ma and Yan2009; Li et al. Reference Li, Wang, Yan, Gong and Ma2018). To ensure optimal relaxation zone lengths, we have conducted a series of simulations with different lengths and draw similar conclusions to Li et al. (Reference Li, Wang, Yan, Gong and Ma2018) and Yan et al. (Reference Yan, Ma, Wang and Wang2019), namely that for the high-wave steepness cases, 1.5 wavelengths are required, while a single wavelength is sufficient for the low-wave steepness cases. In the interest of brevity, these results are not shown here.
We note that in previous studies the surface elevation is typically considered in a convergence test, whereas in our simulations the focus is on the velocity field. Our targets for the convergence tests are surface elevation (wave amplitude), Eulerian-mean velocity and (Lagrangian-mean) drift rates. Here, we report results for the lower-frequency waves (from table 1) of the lowest steepness 
$k{a_w}=0.034$ and the highest steepness 
$k{a_w}=0.126$ examined in § 3.1. In each case, four sets of grids have been tested, which are defined by their spatial resolution, and the three target quantities are examined and compared.
For the lowest-steepness case (
$k{a_w}=0.034$), figure 16 shows the spatial Eulerian-mean (time-averaged) velocity distribution covering the region where our object is placed. The Eulerian-mean velocity 
$\bar {u}_E$ is obtained by time averaging the Eulerian velocity after a quasi-steady state has been achieved, in which the drift speed is constant. The figure demonstrates that, as the spatial resolution becomes higher, the Eulerian-mean velocity becomes very small (at most 
$1\,\%$ of the Stokes drift for the highest steepness), which confirms the (near) absence of Eulerian currents in our numerical wave tank, so that the Lagrangian velocity becomes equal to the Stokes drift (as already shown in § 3.1).
Figure 16. Eulerian-mean velocity scaled by the theoretical Stokes drift 
${u}_{S}$ a small distance below the wave trough 
$a_w=20.0$ mm (
$z=-25.0$ mm) for different spatial resolutions as a function of horizontal position; 
${x_l}$ indicates the location of the left boundary of the NS domain.
Tables 6 and 7 outline the values of our three target quantities obtained for four sets of grids for the lowest-steepness (
$k{a_w}=0.034$) and the highest-steepness (
$k{a_w}=0.126$) cases, respectively. Results are given for wave amplitudes, Eulerian-mean and Lagrangian-mean velocities after a quasi-steady state has been achieved at the location 
$x=22.5$ m, 
$z=-25.0$ mm for 
$k{a_w}=0.034, a_w=20.0$ mm and 
$x=25.2$ m, 
$z=-85.0$ mm for 
$k{a_w}=0.126, a_w=74.0$ mm. The Lagrangian-mean velocities are obtained in the same way as in § 3.1. We find that as the wave steepness is increased, a finer spatial resolution is required for sufficient convergence. Eulerian-mean flows remain small even for the highest-steepness case (at most 
$1\,\%$ of the Stokes drift).
Table 6. Values of the three target quantities of the convergence tests: wave amplitude 
${a_w}$, Eulerian-mean velocity 
$\bar {u}_E$ and drift rate of a Lagrangian particle 
$\bar {u}_L$. Results are shown for four different spatial resolutions for the low-wave steepness case, where 
$\Delta x$ and 
$\Delta z$ represent the grid size in 
$x$ and 
$z$ directions, respectively, 
${N_c}$ is the total number of cells in the NS domain, and 
${a_{in}}$ is the input wave amplitude.
Table 7. Values of the three target quantities of the convergence tests: wave amplitude 
${a_w}$, Eulerian-mean velocity 
$\bar {u}_E$ and drift rate of a Lagrangian particle 
$\bar {u}_L$. Results are shown for four different spatial resolutions for the high-wave steepness case, where 
$\Delta x$ and 
$\Delta z$ represent the grid size in 
$x$ and 
$z$ directions, respectively, 
${N_c}$ is the total number of cells in the NS domain, and 
${a_{in}}$ is the input wave amplitude.
Table 8. Object dimensions of RCBs of different sizes with 
$\rho =781\ \textrm {kg}\ \textrm {m}^{-3}$.
Appendix C. Effects of object rotation and density
The boxes in § 4.1 (table 3) were constrained to prevent rotation during the simulations as these boxes (
$\rho =500\ \textrm {kg}\ \textrm {m}^{-3}, h/l=1.33$) are hydrodynamically unstable. To examine the implications of this assumption and to investigate the effects of object density and rotation, we simulate the drift of a new set of boxes that are hydrodynamically stable. We keep the submergence depth and submerged shape of the new set of RCBs the same as for the previous set but change the height and the density of the boxes. Dimensions of the new set of boxes are given in table 8. The density of the new boxes is 
$\rho =781\ \textrm {kg}\ \textrm {m}^{-3}$. For this density, we do not consider relative sizes smaller than 
$l/\lambda =5\,\%$ to ensure the freeboard of the object (i.e. 
$h-h_d$) is larger than the incident wave amplitude, thus avoiding green-water/plunging impact of the wave onto the object. For this new set of RCBs, we consider three scenarios: without rotation (i.e. with the rotational degree of freedom constrained as in § 4.1) and without viscosity; with rotation but without viscosity; and with rotation and with viscosity. The corresponding object drift velocities are given in figure 17(a). Figure 17(b) provides a comparison of object drift velocities between the RCBs with 
$\rho =500\ \textrm {kg}\ \textrm {m}^{-3}$ (defined in table 3 and considered in § 4.1) and the new set of RCBs with 
$\rho =781\ \textrm {kg}\ \textrm {m}^{-3}$ (defined in table 8) without rotation and without viscosity.
Figure 17. Effects of rotation and density on object drift: (a) celerity-normalized drift velocity as a function of relative size of the objects in table 8 for three different scenarios; (b) celerity-normalized drift as a function of relative size for objects of density 
$\rho =500\ \textrm {kg}\ \textrm {m}^{-3}$ (cf. table 3) and 
$\rho =781\ \textrm {kg}\ \textrm {m}^{-3}$ (cf. table 8) without rotation and without viscosity. The red lines correspond to the theoretical Stokes drift.
It is clear from figure 17 that restraining the rotational degree of freedom causes only a very small change to the object drift. Object rotation also does not significantly change the role that the viscosity plays in altering the drift of the objects. Furthermore, figure 17(b) shows that as long as the submergence depth and submerged shape of the objects are unchanged, the density of the object only has a very minor effect on the object drift.
Appendix D. Effects of viscosity and turbulence
D.1. Reynolds and Keulegan–Carpenter numbers
The problem we consider is one of flow around a 2-D bluff body (a RCB) or a circular object that is freely floating in a surface wave field without a current. We use the characteristic length of the object 
$l$ and the velocity difference between the object and the fluid to define a Reynolds number: 
${{Re}_x}=|{u_{o,x}}-{u_{x}}|l/\nu$, where 
${u_{o,x}}$ and 
${u_{x}}$ represent the magnitude of the horizontal velocity of the object and the fluid, respectively, and 
$l=D$ for ROs. The magnitude of the horizontal velocity of the object 
${u_{o,x}}$ is obtained after a quasi-steady state has been achieved, in which the object oscillates harmonically and drifts at a constant speed. We estimate the magnitude of the horizontal fluid velocity as 
$u_{x}={a_w} \omega \exp (-k h_d)$ for boxes and 
$u_{x}=\int _{-D/2}^{D/2}{a_w} \omega \exp (- k\sqrt {{r^2} - {x^2}})\,\textrm {d}\kern0.06em x/D$ for ROs. The Reynolds numbers for all simulations in category I (see § 4.1) and category III (see § 4.3) are given in tables 9 and 10, respectively. These tables also report the grid size near the moving object boundary: 
${\varDelta _{min}}={\Delta x_{min}}={\Delta z_{min}}$ (the aspect ratio of the mesh near the object is 1). Because OpenFOAM uses collocated grids, which means all of the flow variables are calculated and stored at the cell centroids and these variables vary linearly within a cell, we report 
${\varDelta _{min}}/2$. In order to evaluate whether the mesh density in the vicinity of the boundary is sufficient, we estimate the normal-wall distance 
${y_d}$. We estimate 
${y_d}$ from 
${y_d} = \nu y^+ /u_*$, where the shear velocity is estimated as 
$u_*= \sqrt {(0.058/2)\textrm {Re}^{-0.2}| u_{x,o} - u_{x,f}|^2}$, and the non-dimensional wall distance 
${y^+}$ is set to 1 (Schlichting & Kestin Reference Schlichting and Kestin1961). By comparing 
$y_d$ to 
${\varDelta _{min}}/2$, which is much smaller, we can conclude from tables 9 and 10 that the mesh used in our simulations is fine enough to capture the detailed boundary-layer flows around the object.
Table 9. For category I simulations, exploring the effect of size, Reynolds numbers, smallest mesh sizes near the boundary 
$\varDelta _{min}$, estimates of normal-wall distance 
${y_d}$ and total numbers of cells in the mesh 
$N_c$, and Keulegan–Carpenter numbers 
$K_c$ for RCBs and ROs.
Table 10. For category III simulations, exploring the effect of steepness, Reynolds numbers, smallest mesh sizes near the boundary 
$\varDelta _{min}$, estimates of normal-wall distance 
${y_d}$ and total numbers of cells in the mesh 
$N_c$, and Keulegan–Carpenter numbers 
$K_c$ for RCBs and ROs.
To determine the relative importance of drag and inertial forces and thereby determine the likelihood of boundary-layer separation, we also estimate the Keulegan–Carpenter number: 
${K_c}=|{u_{o,z}}-{u_{z}}|T /l$, 
$T=2 {\rm \pi}/\omega$ is the wave period, and we use the size of the object as the characteristic length scale. In our problem, separation can occur in both the horizontal and the vertical boundary layers and we thus estimate the Keulegan–Carpenter number in both directions. We find the Keulegan–Carpenter number in the vertical directions to always be larger and we therefore report this number in tables 9 and 10. We estimate the magnitude of the vertical fluid velocity as 
${u_{z}}={a_w} \omega \exp (-k{h_d})$ for both boxes and ROs. The Keulegan–Carpenter number 
$K_c$ can be interpreted as the ratio of the magnitude of the oscillatory motion of fluid particles to the length of the object. When the 
$K_c$ number is small, fluid moves only a small distance along the object's boundaries without flow separation, and inertial or diffraction forces will be dominant. For large 
${K_c}$, the fluid particle travels a large distance relative to the size of the object, leading to flow separation and vortex formation. When 
${K_c}<3$, the flow is inertia dominated, and the effects of boundary-layer separation and vorticity are small (e.g. Sumer Reference Sumer2006; Yoon et al. Reference Yoon, Kim, Sadat-Hosseini, Yang and Stern2016; Mohseni et al. Reference Mohseni, Esperanca and Sphaier2018). All of our simulations are in this regime (cf. tables 9 and 10). Furthermore, we do not observe vortex formation and boundary-layer separation in the streamlines and in the velocity and vorticity (using the 
$Q$ criterion) fields.
D.2. Turbulent simulations
Our maximum Reynolds numbers in tables 9 and 10 are 
${Re}>3000$; these numbers are in the sub-critical Reynolds number regime for typical flow around a cylinder, which suggests the boundary layer is laminar but the wake becomes turbulent (Sumer Reference Sumer2006). Although our analysis shows that there is no distinct wake in our simulations, given the Reynolds number of the problem, the flow around the object may become turbulent. To investigate whether the effects of turbulence need to be taken into account (following Yu & Li Reference Yu and Li2013; Li et al. Reference Li, Wang, Yan, Gong and Ma2018), we implement an unsteady Reynolds-averaged NS (URANS) model by introducing a shear stress transport 
$k$-
$\omega$ turbulence model. We consider the cases with the largest Reynolds number in category I (the RCB with 
$l/\lambda =10\,\%, k{a_w}=0.034$) and category III (the RCBs with 
$l/\lambda =8\,\%$, 
$k{a_w}=0.08$ and 
$l/\lambda =10\,\%$, 
$k{a_w}=0.09$). The results are given in tables 11 and 12. Our mesh is fine near the object boundary, as an accurate prediction of viscous forces (wall shear stress) on the object is important. Therefore, in terms of the near-wall treatment, we choose a wall-resolving approach and compare these results with a low-Reynolds-number wall function approach. The boundary conditions for the object boundary are given in table 11.
Table 11. Near-wall treatment and effect of including a turbulence model for a RCB with 
$l/\lambda =10\,\%$ and 
$k{a_w}=0.034$ for two different boundary conditions (BC1 and BC2) in the URANS simulation. Here, 
$k$ is the turbulent kinetic energy (its initial value on both BCs is set to a uniform value of 
$k={1.00\times 10^{-12}}\ \textrm {m}^2\ \textrm {s}^{-2}$), 
${\nu _t}$ is the turbulent viscosity (its initial value on both BCs is set to a uniform value of 
${\nu _t}=0\ \textrm {m}^2\ \textrm {s}^{-1}$), 
$\omega$ is the specific turbulence dissipation rate (its initial value is calculated as 
$\omega =60\nu /(0.075y^2)$, where 
$y$ is the normal distance from the boundary to the first cell centre), 
$\bar {y}^+$ is the time-averaged mean value of 
$y^+$ on the object boundary.
Table 12. Effect of including a turbulence model for RCBs with 
$l/\lambda =8\,\%$, 
$k{a_w}=0.08$ and 
$l/\lambda =10\,\%$, 
$k{a_w}=0.09$ (category III). The boundary conditions BC1 and BC2 are those in table 11, 
$\bar {y}^+$ is the time-averaged mean value of 
$y^+$ on the object boundary.
The mesh used for simulations with and without the turbulence model (for both boundary conditions) is the same. It is clear from table 11 that there is no difference between the two boundary condition (BC) settings, which also confirms that our mesh is fine enough for a wall-resolving approach. Tables 11 and 12 that, compared with the results of the laminar model, a URANS model predicts a similar albeit very slightly lower value of the object drift (
$u_O/c$) along with a similar albeit very slightly larger horizontal motion (
$A_x/a_w$) and a similar albeit slightly smaller vertical motion (
$A_z/a_w$). Sensitivity to the initial value of the specific turbulence dissipation rate 
$\omega$ is small.
Open access
