2.1. Surface waves

In small-amplitude (linear) wave theory, the dimensionless free-surface displacement for monochromatic progressive waves is given by $\eta = ka \cos (x-t)$ , and the dimensionless velocity field is given by

(2.1a) \begin{align} u_x=\varepsilon _w \frac {\cosh (z+kh)}{\cosh {kh}}\cos (x-t), \end{align}

(2.1b) \begin{align} u_z=\varepsilon _w \frac {\sinh (z+kh)}{\cosh {kh}}\sin (x-t), \end{align}

where $\varepsilon _w = ka/\tanh {kh}$ is the wave nonlinearity parameter. Here, $a$ is the wave amplitude, $k$ is the wavenumber, $\omega$ is the wave angular frequency, and $h$ is the undisturbed water depth. The dimensionless coordinates $x$ and $z$ point in the direction of wave propagation and upwards against gravity, respectively, with the origin of $z$ being at the undisturbed free surface. All variables are made dimensionless using the wavenumber $k$ and angular frequency $\omega$ .

2.2. Dimensionless parameters

The parameters that characterise the settling of small, weakly negatively buoyant spheres in surface gravity waves are as follows. The flow field in small-amplitude surface gravity waves is characterised by three length scales, namely the wavelength $L$ (or wavenumber $k = 2\unicode{x03C0} /L$ ), the wave amplitude $a$ , and the water depth $h$ , and by a single time scale, the wave period $T$ (or angular frequency $\omega = 2\unicode{x03C0} /T$ ). The spherical particle is characterised by its density $\rho _p$ , diameter $d_p$ , and terminal settling velocity in quiescent fluid $v^{\prime }_s$ . Together with the fluid properties (namely, density $\rho$ and kinematic viscosity $\nu$ ), and the gravitational acceleration $g$ , there are a total of ten dimensional variables and three fundamental dimensions, which yields seven independent dimensionless parameters. Note that the need to include $v^{\prime}_s$ as an additional parameter arises since its value cannot be determined from other variables when the particle Reynolds number is not small such that the settling velocity is no longer given by the classical Stokes settling velocity. Additionally, the gravitational acceleration is implicitly included in the variables describing the waves since it links the wave period and water depth to the wavelength via the dispersion relation ( $\omega ^2 = gk \tanh kh$ ), but it needs to be included separately when considering particle settling in waves due to its role in the vertical dynamics of particle motion.

We choose the seven dimensionless parameters as follows. The wave-induced flow is characterised by the wave steepness $ka$ and the relative water depth $kh$ . The particle settling in quiescent fluid is characterised by the specific gravity $\gamma = \rho _p/\rho$ and the terminal settling velocity particle Reynolds number ${\textit {Re}}_{p,t} = v^{\prime}_s d_p/\nu$ , which quantifies the importance of fluid inertia. Note that since the terminal settling velocity is not known a priori, a better indicator of the importance of fluid inertia, instead of $Re_{p,t}$ , is the Galileo number ${\textit{Ga}} = v^{\prime}_g d_p/\nu$ , where the characteristic settling velocity is the inertia-gravitational velocity scale $v^{\prime}_g = {((\gamma -1) gd_p)^{1/2}}$ . We will use both dimensionless numbers as appropriate, but only one of them should be included in the list of seven dimensionless parameters.

For the remaining three dimensionless parameters, we choose to characterise the interaction between the particle and wave-induced flow by the ratios of length, time and velocity scales between the particle and the fluid, given respectively by

(2.2) \begin{equation} {{Sz}} = kd_p, \quad {{St}} = \frac {\omega d_p^2 (\gamma +C_m)}{18\nu (1+0.15\, Re_{p,t}^{0.687})} (=\omega \tau _{p,SN}), \quad {{Sv}} = \frac {v^{\prime}_s}{gka/\omega }. \end{equation}

Here, ${\textit{Sz}}$ is the size number, which quantifies the ratio between the particle diameter ( $d_p$ ) and the wavelength ( $L$ ); ${\textit{St}}$ is the Stokes number, which quantifies the particle inertia via the ratio between the time scale for an inertial particle to relax to the flow using the Schiller–Naumann drag correlation ( $\tau _{p,SN}$ ) and the wave period ( $T$ ); and ${\textit{S}v}$ is the settling number (sometimes called the Rouse number), which quantifies the importance of gravitational settling via the ratio between the particle terminal settling velocity ( $v^{\prime}_s$ ) and the characteristic velocity scale of fluid motion in waves ( $u_w=gka/\omega$ ). In the particle time scale, the added mass coefficient for a sphere is taken to be $C_m = {1}/{2}$ , and the particle Reynolds number is determined based on the terminal settling velocity.

2.3. Particle motion under the MRG framework

A widely used simplified version of the MRG equation (without the history force and Faxén corrections, e.g. see Maxey & Riley Reference Maxey and Riley1983; Bergougnoux et al. Reference Bergougnoux, Bouchet, Lopez and Guazzelli2014; DiBenedetto et al. Reference DiBenedetto, Clark and Pujara2022) is given in dimensionless form by

(2.3) \begin{equation} \frac {{\rm d}\boldsymbol{v}}{{\rm d}t} = \frac {(\boldsymbol{u} - \boldsymbol{v} - v_s \boldsymbol{e}_z)}{{\textit{St}}_{{lin}}} + \beta \frac {{\rm D}\boldsymbol{u}}{{\rm D}t}, \end{equation}

where all variables are made dimensionless using $k$ and $\omega$ . Here, $\boldsymbol{v}$ is the particle velocity, $\boldsymbol{u}$ is the undisturbed fluid velocity, $v_s$ is the dimensionless terminal settling velocity, $\boldsymbol{e}_z$ is the unit vector antiparallel to gravity, and $\beta = ({1 + C_m})/({\gamma + C_m})$ . The term on the left-hand side signifies the particle acceleration, the first term on the right-hand side signifies the combined forces of drag and buoyancy, and the second term on the right-hand side signifies the forces due to added mass and fluid acceleration (Bergougnoux et al. Reference Bergougnoux, Bouchet, Lopez and Guazzelli2014). To emphasise that the particle motion is in the small particle Reynolds number regime, the particle Stokes number is written as ${\textit{St}}_{{lin}}=\omega \tau _{p,\textit{lin}}$ where ‘lin’ refers to the fact that the (Stokes) drag is linear with respect to the slip velocity. The particle relaxation time scale is given by $\tau _{p,\textit{lin}} = d_p^2 (\gamma +C_m)/18\nu$ , and the dimensional terminal settling velocity is given by the classical Stokes settling velocity $v^{\prime}_{s} = \tau _{p,\textit{lin}}(1-\beta )g$ .

In the limit of small and weakly negatively buoyant particles ( $(1 - \beta ) \rightarrow 0,\ {\textit{St}} \rightarrow 0$ ), the equation reduces to the simple kinematic model

(2.4) \begin{equation} \boldsymbol{v} = \boldsymbol{u} - v_s \boldsymbol{e}_z, \end{equation}

where the particle velocity is the sum of the undisturbed fluid velocity at the particle position and the particle terminal settling velocity. In this state, the particle drag force and buoyancy are the only forces in play, and they are in balance. When ${\textit{St}}$ is small but finite, a kinematic solution for particle motion can be obtained from an ${\textit{St}}$ expansion. However, for larger or denser particles, an increase in ${\textit{St}}$ (quantifying particle inertia) is also accompanied by an increase in particle Reynolds number (quantifying fluid inertia) since larger particle size and density also result in faster relative motion between the fluid and the particle. This can render the problem outside the bounds of validity of the MRG formulation since the relative motion between the particle and the fluid is no longer characterised by a low Reynolds number.

Bergougnoux et al. (Reference Bergougnoux, Bouchet, Lopez and Guazzelli2014) investigated the applicability of (2.3) and (2.4) for particle settling based on laboratory experiments of spheres settling in a two-dimensional, steady cellular flow. Comparing particle trajectories from the full MRG equation ((2.3) with the history force and Faxén corrections included) with experimental data, they found that the motion of settling particles at low Stokes numbers ( ${\textit{St}} \lesssim 0.1$ ) and finite particle Reynolds numbers ( ${\textit {Re}}_{p,t} \lesssim 10$ ) was well predicted by the simple kinematic model (2.4). This was true as long as one used the measured settling velocity rather than the classical Stokes settling velocity. However, as Bergougnoux et al. (Reference Bergougnoux, Bouchet, Lopez and Guazzelli2014) noted, the applicability of (2.3) and (2.4) for larger particle diameters and densities is questionable since the larger Stokes numbers (quantifying particle inertia) would also be accompanied by larger particle Reynolds numbers (quantifying fluid inertia).

DiBenedetto et al. (Reference DiBenedetto, Clark and Pujara2022) used an ad hoc form of the MRG equation to study particle settling in waves that takes the dimensionless form

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

where the Stokes number based on the Schiller–Naumann drag correlation is given by ${\textit{St}}_{\textit{SN}} = {\omega d_p^2(C_m+\gamma )}/{[18\nu (1+0.15\, {\textit{Re}}_p^{0.687})]}$ and ${\textit{St}}_{{\textit{SN},t}}$ is the same quantity replacing ${\textit {Re}}_p$ with ${\textit {Re}}_{p,t}$ . Using the Schiller–Naumann drag correlation to approximate the drag force for ${\textit {Re}}_p \lesssim 10^3$ yields a modified definition of the particle Stokes number, denoted as ${\textit{St}}_{\textit{SN}}$ , where the subscript SN indicates the use of the Schiller–Naumann drag model. While $\textit{St}_{\textit{SN}}$ is a dynamic quantity that evolves as the particle travels through the flow, ${\textit{St}_{\textit{SN}}}_{,t} = \omega d_p^2 (C_m + \gamma )/[18\nu (1 + 0.15\, {\textit {Re}}_{p,t}^{0.687})]$ is constant. Here, ${\textit{S}v}$ appears implicitly in the dimensionless particle terminal settling velocity since it can be written as $v_s = \varepsilon _w \, {\textit{S}v}$ .

DiBenedetto et al. (Reference DiBenedetto, Clark and Pujara2022) derived an explicit expression for the particle velocity based on an expansion in ${\textit{St}}$ to quantify effects of particle inertia, and used a multiple time scale analysis to obtain the wave-averaged particle drift velocity. They also discussed the mechanisms for enhanced settling and dispersion of inertial particles in waves, extending the Santamaria et al. (Reference Santamaria, Boffetta, Martins Afonso, Mazzino, Onorato and Pugliese2013) results to arbitrary depth ( $kh$ value) and drag force that is nonlinear with respect to the particle slip velocity using the Schiller–Naumann drag correlation.

3.1. Waves and wave-induced flow

The waves and wave-induced flow were characterised as follows. A wave gauge (Toughsonic 3, Senix Corporation) was used to measure the free surface displacement time series $\eta ^{\prime}(t^{\prime})$ . The wave angular frequency $\omega$ was found from the location of the peak in the spectrum of the wave gauge data. The wave amplitude $a$ was estimated by fitting the data to $\eta^{\prime} = a \cos \omega t^{\prime}$ , allowing for an arbitrary phase shift. Finally, the wavenumber $k$ was calculated from the dispersion relationship $\omega ^2 = gk\tanh {kh}$ . Table 1 summarises the wave characteristics and dimensionless wave parameters for each wave case, including the wave speed $c_w$ , characteristic wave-induced fluid velocity $u_w = gka/\omega$ , wave steepness $ka$ , and relative depth $kh$ .

Table 1.

Parameters of each wave case, where C and T in the wave case column indicate that waves are generated by the cylindrical and triangular wavemakers, respectively.

We measured the wave-induced flow using particle image velocimetry (PIV). The PIV field of view spanned 16.5 cm in the $x^{\prime}$ direction, and 26 cm in the $z^{\prime}$ direction, covering 80 % of the water depth over $-26\ \text{cm}\lt z^{\prime}\lt 0$ . No measurements were made in the bottom 5.5 cm of the water column. It was illuminated with an LED line light with average light sheet thickness 1 cm, and seeded with tracer particles ( $d_{50} = 10 \unicode{x03BC}$ m, $\gamma = 1.1$ ; 110P8 Potter Industries). We took images at 20 Hz with a camera (LaVision Imager MX 2M-160, 1936 $\times$ 1216 px) fitted with a 16 mm lens (Tamron) to get magnification factor 7.14 px $\text{mm}^{-1}$ . To obtain the flow velocity field, we used a PIV algorithm with an iterative multi-pass cross-correlation method (DaVis v10.2; LaVision), where the last pass used 64 $\times$ 64 px with 50 % overlap and a Gaussian fit for subpixel accuracy. The final velocity vectors had spatial resolution approximately 4.5 mm, and were checked for quality by examining the relative heights of the first and second correlation peaks. Low-quality vectors ( $\lt 2\, \%$ ) were removed without interpolation. For wave cases W1 and W2, the reflection of the moving water surface on the flume side wall contaminated the image data near the surface. For wave case W3, the vertical gradients in the flow were large, and the flow near the moving surface resulted in tracer streaks. Despite these challenges, we were able to obtain velocity data for all wave cases below depth $z^{\prime} \lesssim -3a$ (figure 2).

Figure 2.

Wave-induced flow spanning $-26\ \text{cm}\lesssim z^{\prime} \lesssim -3a$ . (a,b) Velocity amplitude profiles for oscillatory flow; data with 95 % confidence intervals (circles) and linear wave theory (solid lines). (c) Horizontal Eulerian-mean flow (W1), $\overline {U}^{\prime}_x$ . (d) Vertical Eulerian-mean flow (W1), $\overline {U}^{\prime}_z$ .

To evaluate whether the wave-induced flows conform to linear wave theory (2.1), we compared the PIV velocity measurements with the dimensional velocity field given by

(3.1a) \begin{align} u^{\prime}_x = U^{\prime}_x(z^{\prime}) \cos (kx^{\prime}-\omega t^{\prime}) + \overline {U}^{\prime}_x, \end{align}

(3.1b) \begin{align} u^{\prime}_z = U^{\prime}_z(z^{\prime}) \sin (kx^{\prime}-\omega t^{\prime}) + \overline {U}^{\prime}_z, \end{align}

where $U^{\prime}_x(z^{\prime})$ and $U^{\prime}_z(z^{\prime})$ denote the oscillatory velocity amplitudes in the horizontal and vertical directions, respectively, and $\overline {U}^{\prime}_x$ and $\overline {U}^{\prime}_z$ represent the Eulerian-mean velocities in the horizontal and vertical directions, respectively. For pure monochromatic progressive waves, $U^{\prime}_x(z^{\prime}) = (gka/\omega )\cosh (k(z^{\prime}+h))/\cosh kh$ and $U^{\prime}_z(z^{\prime}) = (gka/\omega )\sinh (k(z^{\prime}+h))/\cosh kh$ , and $\overline {U}^{\prime}_x=\overline {U}^{\prime}_z=0$ . Figure 2 displays the wave-induced velocity amplitudes for wave cases W1–W3, and the Eulerian-mean flow field for wave case W1, obtained from the experimental data by fitting the velocities at each PIV grid point to (3.1). The velocity amplitudes are averaged in the horizontal direction, with the 95 % confidence intervals calculated using the bootstrap technique. The velocity amplitudes show good agreement with linear wave theory, except for a small deviation near the free surface. Furthermore, the magnitudes of Eulerian-mean flow are small compared to the amplitudes of the oscillatory flow.

3.2. Particle settling

We measured the settling dynamics of five solid spherical particles whose properties are summarised in table 2. The nylon (N1, N2, N3) and Torlon (T1) particles were purchased (US Plastics Corp.), whereas the nylon 12 particles (N4) were 3-D printed (selective laser sintering, Fuse 1, Formlabs). We checked that the manufacturer’s reported particle diameters ( $d_p$ ) were accurate to the resolution of standard callipers. Since particle mass (and density) are expected to be altered by water absorption, we stored the particles in water before and in between experiments. The particle mass was then measured using a precision balance scale, allowing us to calculate the particle density ( $\rho _p$ ) and the relative density ( $\gamma$ ).

Table 2.

Dimensional particle parameters. Uncertainty bounds show the 95 $\%$ confidence intervals using the bootstrap method.

We measured the terminal settling velocity ( $v^{\prime}_s$ ) of each particle with 10 repetitions by particle tracking (described below) in quiescent water in the same tank. We also estimated the particle relaxation time scale using the Schiller–Naumann drag correlation $\tau _{p,SN}=d^2_p(\gamma +({1}/{2}))/[18\nu (1+0.15\, Re_{p,t}^{0.687})]$ . Data from particles settling in quiescent fluid showed that the particle vertical velocities reached the terminal settling velocities over a time that was $O(\tau _{p,SN})$ , similar to the exponential approach to the terminal settling velocity in the Stokes drag regime. Thus, we concluded that $\tau _{p,SN}$ provides a reliable estimate for the particle relaxation time scale of our inertial particles.

For N4 particles, the mass measurements were less reliable due to variability in water adherence. So we estimated its relative density from the measured terminal settling velocity, solving for $\gamma$ in $v^{\prime}_s = \sqrt {(4/3) (\gamma - 1)d_p g/C_D}$ , where $C_D=(24/Re_{p,t})(1+0.15\,Re_{p,t}^{0.687})$ is the Schiller–Naumann correlation for the drag coefficient.

Table 3 lists all experiments for particles settling in waves with the relevant dimensionless parameters. In all cases, the particles are small compared to the wavelengths so that the size number is small ( ${\textit{Sz}} \ll 1$ ), and the particles are weakly negatively buoyant so that the specific gravity is just above unity, $(\gamma - 1) \ll 1$ , with the possible exception of T1 particles, where $(\gamma - 1) = 0.4$ . However, the particle relaxation time scale is not small compared to the wave period, and the Stokes number is ${\textit{St}} =O(1)$ . The settling number covers the range $0.8 \leqslant Sv \leqslant 2.4$ , with values below, at and above unity. The $Re_{p,t}$ and $\textit{Ga}$ values span a range that includes particles both below the threshold for path instabilities (N1) and above the threshold for path instabilities (N2–N4, T1).

Table 3.

Dimensionless parameters for each experiment. Ratios of the particle relaxation time to the wave period ( $\tau _{p,SN}/T$ ) are also presented, alongside ${\textit{St}}$ .

3.2.1. The N1, N2, N3, N4 and T1 particles in wave case W1

Experiments with N1, N2, N3, N4 and T1 particles in wave case W1 were used to investigate the particle settling dynamics at different Galileo numbers ( $126 \leqslant {\textit{Ga}} \leqslant 549$ ) at ${\textit{St}} = O(1)$ ( $0.78 \leqslant {\textit{St}} \leqslant 3.0$ ) and ${\textit{S}v} = O(1)$ ( $0.84 \leqslant {\textit{S}v} \leqslant 2.4$ ). To gather data for particle settling dynamics, we simultaneously tracked the spheres and measured the velocity field using a combined particle tracking and PIV approach.

In each experiment, a single sphere was released with a tweezer at approximately 1 cm below the wave trough once the waves were established. We followed this procedure after preliminary experiments suggested that if a particle is dropped from above the water surface, then its initial velocity in the water may exceed its terminal settling velocity, and its trajectory may be altered by air bubbles trapped underneath the particle. The release points were confined to within a few centimetres of the origin in the horizontal direction.

We conducted 10 repetitions for N1, N3, N4 and T1 particles, and 25 repetitions for N2 particles (as detailed in § 3.2.2), varying the initial phase of the wave when the particle was dropped, and excluding drops where the particle travelled outside the light sheet in the spanwise direction. We waited at least 10 minutes between drops to minimise the influence of the previous drop. Occasionally, we had to remove settled particles from the tank before the next release, and this procedure likely altered the structure of the Eulerian-mean flow even though we waited approximately 30 minutes between releases in such cases.

To obtain particle trajectories from the image data, we first removed the background image over the set of images for each trajectory by subtracting the minimum count for each pixel. We determined the particle centroid in each image via the imfindcircles function in MATLAB, which employs a circular Hough transform, setting the radius range based on the expected particle size in pixels. This procedure had an uncertainty of 0.5 pixel ( $\approx 0.07$ mm), which was estimated by checking that circles with centroids shifted from the estimated centroid did not cover the entire particle in the image. Finally, a Savitzky–Golay filter (sgolayfilt in MATLAB) with a second-order polynomial was applied to the data to reduce noise in particle centroid positions. Particle velocities were then calculated using a first-order finite difference scheme using the filtered data. Using the finite difference scheme ensures that the particle velocity data are synchronised with the fluid velocity PIV data.

Since the filter performance depends on the fitting window $L_f$ , we compared the trajectories of the raw and filtered particle centroids for different $L_f$ for the N1 particle (figure 3 a). The filtered data with $L_f = 5$ help to reduce the noise in the particle centroid positions, but larger $L_f$ values change the particle positions more than the measurement uncertainty (0.5 pixel or 0.07 mm) (inset in figure 3 a). Additionally, the effect of filtering on particle velocity shows that using $L_f=5$ reduces noise, but further increases in $L_f$ result in excessive smoothing (figure 3 b). This is evident from the reduction in the magnitude of the velocity spectra at the peak (figure 3 c). Thus, we use $L_f = 5$ points across all data.

Figure 3.

Representative raw and filtered data for N1 particles: (a) positions, (b) velocities, (c) velocity spectra. Insets offer close-up views.

The fluid velocity data were obtained using PIV (as described above), with low-quality vectors ( $\lt 2\, \%$ ) removed without interpolation, which also included vectors where the interrogation window overlapped with any portion of the settling particle. The undisturbed fluid velocity at the particle position was calculated with linear interpolation of the velocity vectors in space, using the particle position at the midpoint between successive images. Since the fluid velocity data are well resolved in space and time, the linear interpolation method was found to be sufficiently accurate. In the interpolation, we excluded fluid velocity vectors within a radius $3d_p$ from each particle centroid and from near the free surface ( $z^{\prime} \gtrsim -3a$ ).

We examined the effect of the exclusion radius on the interpolation; see the example of flow velocities at the particle position for an N3 particle (largest size and $\textit{Ga}$ ) in figure 4. The interpolated fluid velocity at the particle centroid is not sensitive to the size of the exclusion region at $3d_p$ or above. The maximum velocity difference between using exclusion radii $3d_p$ and $5d_p$ is approximately $10^{-3}$ , which is much smaller than the magnitude of the slip velocities (see figure 7(e) below). Thus, we concluded that an exclusion radius of $3d_p$ was sufficient to obtain the background fluid velocities at the particle centroid positions.

It is worth noting that while our experimental set-up allows us to obtain well-resolved simultaneous measurements of the particle trajectories and background fluid velocities over most of the water column, we are not able to capture the particle wake within the fluid since the light is delivered from below the tank, and the wake flow immediately above the particle is blocked by the particle. Some part of the wake structure may influence the fluid velocity vectors in the vicinity of the particle, but since the wake flow is not measurable, these vectors are excluded from the analysis as explained above.

Figure 4.

Representative dimensionless fluid velocities interpolated to the particle centroid for varying exclusion radius; N3 particle in wave case W1.

4.1. The N1 particle

4.1.1. Approach to the kinematic model

As mentioned above, the MRG equation (2.3), which is valid for small particle Reynolds number ( ${\textit {Re}}_p \ll 1$ ), reduces to the kinematic model (2.4) in the limit of small Stokes number ( ${\textit{St}} \ll 1$ ) for weakly negatively buoyant particles ( $(1 - \beta) \ll 1$ ) because the first term on the right-hand side of the equation dominates. Here, we have found that N1 particles in wave case W1, which have Stokes number ${\textit{St}} = O(1)$ and particle Reynolds number ${\textit {Re}}_p = O(100)$ , also follow the kinematic model. In the following, we show how this finding can be rationalised in the context of the ad hoc extension of the MRG equation for large particle Reynolds number (2.5).

Starting with (2.5), we consider the scenario where the fluctuations in the particle slip velocities are small relative to the slip velocity: $\lvert \Delta v_{\textit{slip}} \rvert /v_{\textit{slip}} \ll 1$ . We observe this in the N1 particle data when the slip approaches the terminal settling velocity. As ${\textit {Re}}_p\ (= |\boldsymbol{v}^{\prime}-\boldsymbol{u}^{\prime}|\,d_p/\nu )$ remains nearly constant in these condition, ${\textit{St}}_{\textit{SN}}$ can be approximated by a constant, which in this case is ${\textit{St}}_{\textit{SN},t}$ . This gives

(4.1) \begin{equation} \frac {\mathrm{d}\boldsymbol{v}}{\mathrm{d}t} = \frac {(\boldsymbol{u} - \boldsymbol{v} - v_s \boldsymbol{e}_z)}{\textit{St}_{{\textit{SN},t}}} + \beta \frac {\mathrm{D}\boldsymbol{u}}{\mathrm{D} t} . \end{equation}

We now simplify this equation further using the following scalings, where $\varepsilon _w$ is the wave nonlinearity parameter: $\boldsymbol{u} = O(\varepsilon _w)$ , $\boldsymbol{v} = O(\varepsilon _w)$ , ${\textit{St}} = O(1)$ , $(1-\beta) =O(\varepsilon _w)$ and ${\textit{S}v} = O(1)$ , which implies $v_s = O(\varepsilon _w)$ since $ v_s=\varepsilon _w\,{\textit{S}v}$ . We can write the last term as

(4.2) \begin{align} \beta \frac {\mathrm{D}\boldsymbol{u}}{\mathrm{D}t} & = \frac {\mathrm{D}\boldsymbol{u}}{\mathrm{D}t} - (1-\beta )\frac {\mathrm{D}\boldsymbol{u}}{\mathrm{D}t} \nonumber \\[6pt] & = \frac {\mathrm{d}\boldsymbol{u}}{\mathrm{d}t} + (\boldsymbol{u}-\boldsymbol{v})\boldsymbol{\cdot} \nabla \boldsymbol{u} - (1-\beta ) \frac {\mathrm{D}\boldsymbol{u}}{\mathrm{D}t} \nonumber \\[6pt] & = \frac {\mathrm{d}\boldsymbol{u}}{\mathrm{d}t} + O(\varepsilon _w^2). \end{align}

By dropping terms of $O(\varepsilon _w^2)$ , (4.1) can be written as

(4.3) \begin{equation} \frac {\mathrm{d} (\boldsymbol{v}-\boldsymbol{u}+v_{s}\boldsymbol{e}_z)}{\mathrm{d}t} + \frac {(\boldsymbol{v}-\boldsymbol{u}+v_{s}\boldsymbol{e}_z)}{{\textit{St}}_{\textit{SN},t}}=0, \end{equation}

where we have used $\mathrm{d}v_{s}/\mathrm{d}t =0$ . The solution is given by

(4.4) \begin{equation} \frac {(\boldsymbol{v}-\boldsymbol{u}+v_{s}\boldsymbol{e}_z)}{(\boldsymbol{v}-\boldsymbol{u}+v_{s}\boldsymbol{e}_z)_{t_0}} = \exp \left (-\frac {(t-t_0)}{{\textit{St}}_{\textit{SN},t}}\right ), \end{equation}

which shows that the particle slip velocities approach the kinematic model with an exponential decay with a dimensionless time scale given by ${\textit{St}}_{\textit{SN},t}$ , with an ‘initial condition’ $(\boldsymbol{v}-\boldsymbol{u}+v_{s}\boldsymbol{e}_z)_{t_0}$ at $t=t_0$ . At $(t-t_0) \gg {\textit{St}}_{\textit{SN},t}$ , we recover the kinematic model $\boldsymbol{v} = \boldsymbol{u} - v_{s}\boldsymbol{e}_z$ .

Figure 8 shows how the data compare against (4.4) during the particles’ approach to the kinematic model. We have to select the ‘initial time’ $t_0$ in (4.4) with care since the noise in the data is amplified when $(\boldsymbol{v}-\boldsymbol{u}+v_{s}\boldsymbol{e}_z)_{t_0}$ is small, and the particle release mechanism will impart some unknown initial motion to the particle. In the horizontal direction, we select $t_{x0}$ as the time when the horizontal slip velocity reaches its peak (see figure 7 a), after which it decays towards zero. For trajectories where the horizontal slip velocity does not have a clear peak, $t_{x0}$ is chosen as the time when $|(v-u)_x|$ begins to be smaller than $0.15 v_s$ . In the vertical direction, $t_{z0}$ is defined as the time when $|(v-u)_z|$ exceeds $0.7 v_s$ , by which time ${\textit{St}}_{\textit{SN}}$ is also within 20 % of ${\textit{St}}_{\textit{SN},t}$ .

Overall, the data show a reasonable agreement with the prediction in the horizontal direction, while it appears that the approach to the kinematic model in the vertical direction takes longer than predicted. The slower approach to the kinematic model in the vertical could be due to the history force (ignored here), which has been shown in previous experiments to have this effect (Mordant & Pinton Reference Mordant and Pinton2000). Additionally, we see that forces related to drag, buoyancy, added mass and fluid acceleration are all important in the approach to the kinematic model. This is different to the case for particles with ${\textit{St}} \ll 1$ , where the kinematic model results from the drag–buoyancy forces being dominant. Finally, the Stokes number is important only to describe the transient phase in approaching the kinematic model. After that, only ${\textit{S}v}$ influences the particle motion.

Figure 8.

Evolution of $\boldsymbol{v}-\boldsymbol{u}+v_s \boldsymbol{e}_z$ for N1 particle data normalised by the initial condition in different directions: (a) horizontal, (b) vertical. Black solid lines are the analytical model (4.4).

4.1.2. Enhanced settling and dispersion

Since we found that ${\textit{St}}$ does not influence the particle motion after an initial transient phase, we can use the experimentally validated kinematic model (2.4) for particles below the $\textit{Ga}$ threshold for path instabilities to perform the same multiple time scale analysis as in DiBenedetto et al. (Reference DiBenedetto, Clark and Pujara2022) to obtain the wave-averaged particle drift velocities. We find these (in dimensionless form) to be

(4.5) \begin{equation} v_{x\text{-}drift} = \frac {1}{1 + v_s^2} u_{{SD}}, \quad v_{z\text{-}drift} = -v_s\left [1 + \frac {1}{1 + v_s^2} u_{{SD}} \right ], \end{equation}

where $u_{{SD}} = \varepsilon _{w}^2 \cosh {2(Z_p +kh)}/2\cosh ^2kh$ is the classical dimensionless Stokes drift velocity with wave-averaged vertical particle position $Z_p$ , and $v_s$ is the dimensionless particle settling velocity, with all velocities made dimensionless using the wave speed $\omega /k$ . Equation (4.5) predicts that the horizontal drift is reduced while the vertical drift is enhanced, with both effects related to how the particle samples the fluid flow, and quantified by the term ${u_{{SD}}}(1/({1 + v_s^2}))$ . The settling velocity enhancement increases in percentage terms as the wave steepness increases and as the dimensionless particle settling velocity decreases.

To check how well (4.5) predicts the settling velocity enhancement, we compute the net settling velocity using a single particle trajectory, following the double wave-averaging approach in DiBenedetto et al. (Reference DiBenedetto, Clark and Pujara2022). Defining the wave average of quantity $s$ over wave period $T$ as $\overline {s} = ({1}/{T}) \int _{t}^{t+T} s(\tau )\,{\rm d}\tau$ , we wave average the particle vertical position and velocity once to remove the wave-induced oscillations, and once more to remove the influence of the initial wave phase at the release point. The double wave-averaged vertical velocity $\overline {\overline {v_z}}$ is then only a function of the double wave-averaged vertical position $\overline {\overline {z^{\prime}_p}}$ .

Figure 9.

(a) Double wave-averaged velocity profiles of the N1 particles made dimensionless with the terminal settling velocity. Markers show data with 95 % confidence intervals, dashed lines show predicted enhanced settling (4.5), and vertical black dash-dotted lines indicate 95 % confidence interval of the terminal settling velocity (table 2). (b) Vertical profile of the standard deviation of the wave-averaged horizontal particle position for N1 particles.

Figure 9(a) shows profiles of $\overline {\overline {v_z}}$ as a function of $\overline {\overline {z^{\prime}_p}}$ for N1 particle data binned in vertical increments $0.1(z^{\prime}/h)$ and averaged across different particle releases. For $\overline {\overline {z^{\prime}_p}}/h\gt -0.3$ , the vertical drift is affected by the initial relaxation of the particle to the wave-induced flow, where its vertical slip velocity is smaller than the terminal settling velocity. Below this region, where the particle motion follows the kinematic model, the data appear to show a very small ( $1\,\%$ ) settling enhancement, which in turn is much greater than the enhancement predicted by (4.5) at these depths.

However, the settling enhancement suggested by the data should be interpreted with care since it falls within the uncertainty in the terminal settling velocity (dash-dotted lines). While (4.5) does not correctly predict the settling enhancement, it also does not account for the presence of Eulerian-mean flow. As seen in figure 2, although the Eulerian-mean flow is small compared to the wave-induced oscillatory flow, it is approximately $1\,\%$ of $v_s$ , and could be the cause of the settling enhancement in the data. We conclude that although a wave-induced settling enhancement cannot be ruled out, it is challenging to conclusively demonstrate it given the unavoidable sources of uncertainty and the smallness of the enhancement predicted from theory.

Apart from enhanced settling, DiBenedetto et al. (Reference DiBenedetto, Clark and Pujara2022) also provided a discussion of the mechanisms by which inertial particles disperse in the horizontal direction while settling in waves. While no horizontal dispersion is predicted by divergence of the drift velocities (4.5), there is dispersion in the wave-resolved sense due to variations in the wave phase at the point of initial particle release. The variation in the initial phase can be understood as vertical variations in the effective wave-averaged initial positions, which are acted upon by the shear in the Stokes drift near the surface.

To quantitatively inspect the horizontal dispersion of the N1 particles, we compute the standard deviation of their wave-averaged horizontal positions. As the N1 particles undergo five wave periods while settling in the measurement area, we obtain their wave-averaged particle positions ( $\overline {x^{\prime}_p}$ and $\overline {z^{\prime}_p}$ ) for each wave period. The standard deviations of the horizontal particle positions with corresponding mean vertical positions are then determined by bootstrapping using data for each wave period. Note that since the particles are subjected to similar wave-induced flows, their vertical displacements are nearly uniform across the particles for each wave period, leading to negligible confidence intervals in $\overline {z^{\prime}_p}$ .

Figure 9(b) shows the vertical profile of the dimensionless standard deviations of the wave-averaged horizontal particle positions. While the standard deviation increases near the surface as particles with different initial wave phases spread out, it remains relatively constant in $\overline {z^{\prime}_p}/h\lt -0.2$ , which is consistent with the mechanism discussed in DiBenedetto et al. (Reference DiBenedetto, Clark and Pujara2022). Note that the small changes in the standard deviation in $\overline {z^{\prime}_p}/h\lt -0.2$ could be due to the weak Eulerian-mean flows, and they are not investigated further.

4.2. Path instabilities

To understand why certain particles show large deviations from the kinematic model, we compared how spheres settle in waves with the corresponding data in quiescent water. Figure 10 shows images of particles N1 and N3 (smallest and largest $\textit{Ga}$ values, respectively) with black markings on them to visualise particle rotation. The markings on the N1 particle do not move, showing that the particle had little to no rotation, as expected for a particle settling in a quiescent fluid and an irrotational flow. However, the markings on the N3 particle change consistently, indicating that the particle rotates while settling in both quiescent fluid and in irrotational waves. Since there is no mechanism by which irrotational wavy flow can cause rotations of spherical particles, we conclude that the oscillations in horizontal slip velocity in wavy flow for our larger $\textit{Ga}$ particles must be due to path instabilities. In other words, the deviations of particle motion from the kinematic model are due to fluid inertia rather than particle inertia.

Figure 10.

The N1 and N3 particles with surface markings settling in quiescent water and waves.

To explore the path instabilities more quantitatively, we calculated spectra of the particle slip velocity $P(f)$ , where $f$ is the frequency. Specifically, we used the dimensional slip velocity in waves $(\boldsymbol {v}^{\prime }-\boldsymbol {u}^{\prime })$ , and the dimensional particle velocity in quiescent conditions $\boldsymbol v^{\prime }$ . Note that by using the particle slip velocity, we avoid the spectral peak produced by the wave-induced oscillations in particle motion. We used the data over the vertical region $-0.8 \lt z^{\prime}/h \lt -0.3$ , where the vertical slip has already relaxed to the particle’s terminal settling velocity, but the particle motion is still affected by the wave-induced flow. To compute the spectra, we first matched the length of the slip velocity data time series for all data of a given particle in both waves and quiescent conditions by excluding the last few data points in order to ensure consistency in the frequency resolution. Typically, no more than three data points were excluded. We also subtracted the mean from the time series of each particle release to obtain time series of zero-mean velocity fluctuations (slip velocity in waves $(\boldsymbol {v}^{\prime }-\boldsymbol {u}^{\prime })^*$ and particle velocity in quiescent conditions $(\boldsymbol{v}^{\prime })^*$ ). The spectrum for each particle release was then computed using the fast Fourier transform on the zero-mean time series before being ensemble averaged across different particle releases (as recommended in Bendat & Piersol Reference Bendat and Piersol2011). Finally, the ensemble mean spectrum and its confidence interval at each frequency were determined via bootstrapping across the spectra.

Figure 11.

Normalised ensemble mean spectra of particle slip velocity fluctuations in waves and quiescent fluid for (a) N1, (b) N2, (c) N3, (d) N4, (e) T1. Vertical dashed lines indicate the dimensionless wave frequency ( $d_p/Tv^{\prime}_g$ ), with $d_p/Tv^{\prime}_g \approx 0.22$ in (d) not shown.

Figure 11 shows the ensemble mean spectra for each particle type using 10 particle releases in both quiescent fluid and waves. While the spectra for the N1 particle ( ${\textit{Ga}} = 126$ , below the ${\textit{Ga}} \approx 155$ threshold for path instabilities) show no peaks, the spectra of the other particles show clear peaks with magnitudes 1–2 decades larger in the horizontal direction compared to the vertical. The normalised frequencies at which the spectral peaks occur in our data agree with previously reported values for similar $\textit{Ga}$ - and $\gamma$ -valued particles settling in quiescent conditions (Veldhuis & Biesheuvel Reference Veldhuis and Biesheuvel2007; Raaghav et al. Reference Raaghav, Poelma and Breugem2022; Cabrera-Booman et al. Reference Cabrera-Booman, Plihon and Bourgoin2024). Since the wake instabilities of spheres settling in quiescent fluid are determined solely by particle and fluid properties, the spectral peaks occur at the same frequencies regardless of whether the particles are settling in quiescent fluid or in waves, and these frequencies are distinct from the wave frequency, we conclude that the deviations from the wave-induced flow in the horizontal direction (figure 7) are indeed due to asymmetric wakes and path instabilities, and that these path instabilities are unaffected by the waves. Note that the lack of influence of the oscillatory wave-induced flow on the path instabilities caused by particle wakes is different to, for example, the case where spheres or cylinders are forced to oscillate transversely across a free stream flow (Williamson & Roshko Reference Williamson and Roshko1988; Peter & De Reference Peter and De2016), in which the wake vortex frequency is influenced by the forced oscillation frequency.

Figure 12.

The N2 particles settling in quiescent fluid: (a) example trajectories; (b) normalised horizontal particle velocities; (c) horizontal particle velocity autocorrelation function; (d) normalised spectra of horizontal particle velocities.

While the data in figure 11 suggest that the wave-induced flow has little influence on the particle path instabilities, the spectral peaks are broad and the frequency resolution is coarse. Furthermore, the data are from a single wave case, leaving open the question of whether the waves might affect the path instabilities at different wave frequencies. We address these points using N2 particle data spanning wave cases W1, W2 and W3 (§ 3.2.2) and further analysis.

Figure 12 shows further analysis of N2 particles settling in quiescent fluid. In the example trajectories in figure 12(a), we observe two distinct regimes of path instabilities: a low-frequency oscillating oblique regime (circles), and a chaotic regime (asterisks). This is not surprising since N2 particles are at the boundary between the two different regimes in $\textit{Ga}$ – $\gamma$ space (see figure 5 in Raaghav et al. Reference Raaghav, Poelma and Breugem2022), and the particle properties have small variations across different particles of each type. In the oscillating oblique regime, there are multiple oscillations, but the sign of the horizontal velocity remains unchanged, resulting in a non-negligible mean (figure 12 b). In contrast, in the chaotic regime, there is one or fewer oscillations, but they occur around a near-zero mean horizontal velocity with a change in sign (figure 12 b). These differences are reflected in the horizontal velocity autocorrelations ( $R_{xx} ({\text{lag}}) = \langle v^{\prime *}_x(t^{\prime})\, v^{\prime *}_x(t^{\prime}+{\text{lag}})\rangle /\langle (v^{\prime *}_x)^2\rangle$ , where $\langle \boldsymbol{\cdot} \rangle$ denotes a mean over the time series; see figure 12 c), which show a distinct peak for the oscillating oblique regime, and no distinct peak for the chaotic regime. The horizontal velocity spectra (figure 12 d) show that the peaks for the two regimes are at different frequencies. While the spectral peaks for the oscillating oblique regime occur at frequencies corresponding to the lags associated with the peaks in the autocorrelation functions ( $\approx 0.6$ and $0.8$ s), the spectral peaks for the chaotic regime occur at the smallest resolved frequency (i.e. the first non-zero bin) since the trajectories are too short to observe complete oscillations. We also note that the variance of the horizontal particle velocity for the oscillating oblique regime is much smaller than that for the chaotic regime, which results in a smaller spectral peak.

Figure 13.

The N2 particles settling in waves. Examples of normalised horizontal slip velocities and the corresponding horizontal slip velocity autocorrelation functions for different wave cases: (a,d) W1, (b,e) W2, (c,f) W3.

Figure 13 shows the corresponding analysis for N2 particles settling in waves. The wave-induced oscillatory motion of particles makes it difficult to identify different path instability regimes from particle trajectories, but we find that the horizontal slip velocities show features that resemble the same two different regimes: a low-frequency oscillating oblique regime (circles), and a chaotic regime (asterisks) (figure 13 a–c). As in the data from quiescent fluid, there are multiple oscillations with the same sign for the oscillating oblique regime, whereas there are one or fewer oscillations with a change in sign in the chaotic regime. The corresponding slip velocity autocorrelation functions (figure 13 d–f) show distinct peaks for the former, while no distinct peak is observed for the latter, just as in the quiescent fluid case.

Figure 14.

Normalised ensemble mean spectra of N2 particle slip velocity fluctuations in waves and quiescent fluid: (a) oscillating oblique regime; (b) chaotic regime. Vertical dashed lines indicate the dimensionless wave frequency ( $d_p/Tv^{\prime}_g$ ), with different colours corresponding to different wave cases.

With the knowledge of the two different path instability regimes, we can return to the comparison between the spectra in quiescent fluid and in waves presented in figure 11, and use the N2 particle data to confirm that the magnitude and width of the spectral peaks are unaffected by the waves. Figure 14 shows the ensemble mean spectra split into the low-frequency oscillating oblique regime and the chaotic regime. Each spectrum comes from the mean of approximately 12 particle releases (from a total of 25 releases, with both regimes being equally likely), and the particle in each release experiences a significant wave-induced flow for at least one wave period while settling through the region of analysis ( $-0.8 \lt z^{\prime}/h \lt -0.3$ ). We note that while the coarse frequency resolution is a fundamental limitation of the short particle trajectories (see also Veldhuis & Biesheuvel Reference Veldhuis and Biesheuvel2007; Cabrera-Booman et al. Reference Cabrera-Booman, Plihon and Bourgoin2024), which themselves are limited by the finite depth of the wave-influenced region, the time domain analysis in figures 12 and 13 confirms the robustness of the spectral peaks.

We observe that in both path instability regimes, the presence of the waves does not affect either the magnitude or the width of the spectral peak in figure 14. This is true across the three different wave cases with different frequencies, where $T/T_{\textit{oscillating oblique}} = 0.80, 1.4, 1.0$ and $T/T_{\textit{chaotic}} = 0.40, 0.70, 0.51$ , with $T_{\textit{oscillating oblique}}$ and $T_{\textit{chaotic}}$ being the periods associated with spectral peaks for the two path instability regimes, and $T$ being the wave period. Moreover, the spectral peaks occur at dimensionless frequencies that are distinct from those associated with the wave period. These results demonstrate that path instabilities of settling particles observed in quiescent fluid persist in the unsteady wave-induced flows, and account for the deviation of particle velocities from the kinematic model. We suspect that the lack of interaction between the particle wakes and the wavy flow is due to the smallness of ${\textit{Sz}}$ (table 3), which implies a large-scale separation between the particles and the flow length scale.