2.1 Sediment motion

The motion of the sediment phase is computed by evaluating the net force and torque acting on each Lagrangian particle and integrating the following equations for the particle positions, $\boldsymbol{x}_{p}$ , velocities, $\boldsymbol{u}_{p}$ , and angular velocities, ${\bf\omega}_{p}$ ,

(2.1a-c ) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}(\boldsymbol{x}_{p})=\boldsymbol{u}_{p},\quad \frac{\text{d}}{\text{d}t}(\boldsymbol{u}_{p})=\frac{1}{m_{p}}\sum \boldsymbol{F}_{p},\quad \frac{\text{d}}{\text{d}t}({\bf\omega}_{p})=\frac{1}{i_{p}}\sum \boldsymbol{T}_{p},\end{eqnarray}$$

where $m_{p}$ and $i_{p}$ are the particle mass and moment of inertia respectively. Throughout this work the particle diameter, $d_{p}$ , is assumed to be of the order of or smaller than the Eulerian grid spacing, ${\it\Delta}$ . The particle–fluid interface is therefore unresolved, and a point-particle formulation is used to model the net force, $\boldsymbol{F}_{p}$ , and torque, $\boldsymbol{T}_{p}$ , as a summation of specific hydrodynamic and interparticle contributions:

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{p}=\underbrace{\boldsymbol{F}_{g}}_{\text{Gravity}}+\underbrace{\boldsymbol{F}_{pr}}_{\text{Pressure}}+\underbrace{\boldsymbol{F}_{d}}_{\text{Drag}}+\underbrace{\boldsymbol{F}_{l}}_{\text{Lift}}+\underbrace{\boldsymbol{F}_{am}}_{\text{Added mass}}+\underbrace{\boldsymbol{F}_{c}}_{\text{Collision}}, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{T}_{p}=\underbrace{\boldsymbol{T}_{h}}_{\text{Hydrodynamic}}+\underbrace{\boldsymbol{T}_{c}}_{\text{Collision}}+\underbrace{\boldsymbol{T}_{r}}_{\text{Rolling}}. & \displaystyle\end{eqnarray}$$

The formulation of each term is discussed below.

2.1.1 Hydrodynamic forces

In this work, the total hydrodynamic force acting on each particle contains contributions from gravity, pressure, drag, lift and added mass effects along the lines of a modified Maxey & Riley (Reference Maxey and Riley1983) equation for inertial particle motion. When developing a system of Lagrangian equations, inclusion or exclusion of individual terms can be justified on a case by case basis, with model closures selected for the desired application. The relatively low specific gravity of marine sediments and the oscillatory nature of the wave bottom boundary layer have motivated the inclusion of lift and added mass forces, in addition to the more commonly used drag, gravity and pressure contributions.

Table 1. Models for hydrodynamic forces and torques acting on Lagrangian point particles. The subscript $_{|p}$ denotes a fluid property evaluated at the particle centre.

The expressions used to compute the hydrodynamic contributions to (2.2) and (2.3) are given in table 1, including closure models where applicable. The equations for hydrodynamic forces and torques require several fluid properties to be interpolated to individual particle locations. The interpolation of fluid property ${\it\phi}_{f}$ from the control volume (CV) centres, $\boldsymbol{x}$ , to the particle coordinates, $\boldsymbol{x}_{p}$ , is defined as

(2.4) $$\begin{eqnarray}{\it\phi}_{f|p}=\displaystyle \mathop{\sum }_{cv}\mathscr{I}(\boldsymbol{x},\boldsymbol{x}_{p}){\it\phi}_{f},\end{eqnarray}$$

where the interpolant function, $\mathscr{I}(\boldsymbol{x},\boldsymbol{x}_{p})$ , is a trilinear function which uses the eight nearest grid points on the three-dimensional grid. This compact second-order kernel is sufficient to obtain smooth particle trajectories for the sand particles of interest when small time steps are used to resolve particle collisions, and requires a lower computational overhead compared with higher-order interpolants.

The gravity force, $\boldsymbol{F}_{g}$ , is the weight of the particle, with $g_{y}$ taken to be $9.81~\text{m}~\text{s}^{-2}$ . When combined with the pressure force, $\boldsymbol{F}_{pr}$ , which involves the gradient of the total fluid pressure, the particles respond to both buoyancy and local variations in dynamic pressure.

The drag force, $\boldsymbol{F}_{d}$ , is computed based on the relative particle velocity, $|\boldsymbol{u}_{p}-\boldsymbol{u}_{f|p}|$ , the Stokes flow relaxation time, ${\it\tau}_{st}={\it\rho}_{p}d_{p}^{2}/18{\it\mu}_{0}{\it\Theta}_{f}$ , and the drag coefficient, $C_{d}(Re_{p},{\it\Theta}_{p})$ . Here, ${\it\Theta}_{p|p}$ is the local particle fraction ( ${\it\Theta}_{f}=1-{\it\Theta}_{p}$ corresponds to the fluid fraction) interpolated to the particle centre, and the particle Reynolds number based on the superficial velocity is

(2.5) $$\begin{eqnarray}Re_{p}=\frac{{\it\Theta}_{f|p}{\it\rho}_{f}d_{p}|\boldsymbol{u}_{p}-\boldsymbol{u}_{f|p}|}{{\it\mu}_{0}},\end{eqnarray}$$

where ${\it\mu}_{0}$ is the dynamic viscosity of the fluid phase without considering the effects of the particle suspension (§ 2.2). The functional form of $C_{d}(Re_{p},{\it\Theta}_{p})$ can significantly influence the collective behaviour of particles and is especially important for simulating coastal sediments, where settling velocities can vary by an order of magnitude in the boundary layer due to concentration gradients (Richardson & Zaki Reference Richardson and Zaki1954). Fully resolved simulations of flow through fixed arrangements of spherical particles have allowed for parametrization of the effective drag coefficient for a wide range of ${\it\Theta}_{p}$ and $Re_{p}$ , and recently several relationships for $C_{d}(Re_{p},{\it\Theta}_{p})$ have been proposed (see Tang et al. (Reference Tang, Peters, Kuipers, Kriebitzsch and Hoef2015) for a recent review). For this work, the correlation of Tenneti, Garg & Subramaniam (Reference Tenneti, Garg and Subramaniam2011) is chosen because (i) it is valid up to $Re_{p}=300$ and (ii) it has been developed as a correction to the dilute limit, $C_{d}(Re_{p},0)$ , rather than the low-Reynolds-number limit, $C_{d}(0,{\it\Theta}_{p})$ . This attribute allows particle shape effects to be more easily incorporated into the drag law by, for example, specifying the $C_{d}(Re_{p},0)$ relationship for angular natural sand particles (Fredsøe & Deigaard Reference Fredsøe and Deigaard1992) rather than the typical relation for smooth spherical particles (Schiller & Naumann Reference Schiller and Naumann1935). To account for the wide distribution of particle sizes found in naturally graded sands, the polydisperse correction to the drag force proposed by Beetstra, Van der Hoef & Kuipers (Reference Beetstra, Van der Hoef and Kuipers2007a ,Reference Beetstra, Van der Hoef and Kuipers b ) is used. They found that computing $Re_{p}$ based on the Sauter mean diameter of the local particle population, $\langle d\rangle$ , and scaling the individual particle drag force by a factor $d_{p}/\langle d\rangle$ provided good results for both bi- and polydisperse systems with wide size variations.

Only shear based contributions to the lift force, $\boldsymbol{F}_{l}$ are considered, and the closure of Saffman (Reference Saffman1965) is used for the lift coefficient, $C_{l}$ . As with the drag closures, new experimental techniques and FRS will hopefully lead to improved relations for $C_{l}$ in conditions relevant to coastal boundary layers, and future developments can easily be incorporated into the present model. Similarly, investigations of virtual mass effects on particles in non-dilute suspensions should provide improved understanding on how to effectively model the added mass force, $\boldsymbol{F}_{am}$ , in these conditions. In the present work, the expression for an isolated rigid sphere in non-uniform flow is used, with added mass coefficient $C_{am}=0.5$ .

A particle will also experience a hydrodynamic torque, $\boldsymbol{T}_{h}$ , due to its relative rate of rotation in a viscous fluid, ${\bf\omega}_{rel}=(1/2)(\boldsymbol{{\rm\nabla}}\times \boldsymbol{u}_{f})_{|p}-{\bf\omega}_{p}$ . The expression for $\boldsymbol{T}_{h}$ contains a torque coefficient, $C_{t}$ , which is modelled as a function of the particle rotational Reynolds number,

(2.6) $$\begin{eqnarray}Re_{r}={\displaystyle \frac{{\it\rho}_{f}d_{p}^{2}|{\bf\omega}_{rel}|}{4{\it\mu}_{0}}}.\end{eqnarray}$$

Values of $C_{t}$ were determined by Pan, Tanaka & Tsuji (Reference Pan, Tanaka and Tsuji2001) by matching the Stokes solution at low $Re_{r}$ to experimental data at higher $Re_{r}$ .

Figure 1. Schematic of two particles colliding and the variables used in the soft-sphere collision treatment.

As our fundamental understanding of multiphase flow physics improves through new experimental and FRS databases, it is expected that the predictive capability of point-particle models will also improve, especially where irregularly shaped natural sediments are concerned. For example, Hölzer & Sommerfeld (Reference Hölzer and Sommerfeld2008) and Zastawny et al. (Reference Zastawny, Mallouppas, Zhao and Van Wachem2012) have developed orientation-dependent closures for shape-dependent lift, drag and torque coefficients obtained from FRS of irregular shaped particles (fibres and ellipsoids). Book-keeping of individual particle orientations within a point-particle framework represents a significant barrier to implementing such relationships, but such an approach may be a promising next step for simulating natural sand particles.

2.1.2 Particle collisions

In regions of high particle concentration encountered near the sediment bed, particle motion is dominated by collisions and enduring contacts between particles. To include these collisions in a physically realistic way, a soft-sphere DEM based on the work of Cundall & Strack (Reference Cundall and Strack1979) is employed. Figure 1 schematically illustrates a pair of particles, denoted $i$ and $j$ , undergoing collision and the variables used to compute the collision force, $\boldsymbol{F}_{c}$ , collision torque, $\boldsymbol{T}_{c}$ , and rolling torque, $\boldsymbol{T}_{r}$ . During collision a unit normal vector, $\boldsymbol{n}_{ij}$ , points from particle $i$ to particle $j$ ,

(2.7) $$\begin{eqnarray}\boldsymbol{n}_{ij}=\frac{\boldsymbol{x}_{p,j}-\boldsymbol{x}_{p,i}}{|\boldsymbol{x}_{p,j}-\boldsymbol{x}_{p,i}|}.\end{eqnarray}$$

In order to mimic the elastic deformation that occurs during collision, the two colliding particles are allowed to overlap slightly in the normal direction, by an amount

(2.8) $$\begin{eqnarray}{\it\delta}_{ij}^{n}=0.5(d_{p,i}+d_{p,j})-|\boldsymbol{x}_{p,i}-\boldsymbol{x}_{p,j}|.\end{eqnarray}$$

The total relative velocity of the two spheres at the contact point can be written as

(2.9) $$\begin{eqnarray}\boldsymbol{u}_{ij}=\boldsymbol{u}_{p,i}-\boldsymbol{u}_{p,j}+\left({\textstyle \frac{1}{2}}d_{p,i}{\bf\omega}_{p,i}+{\textstyle \frac{1}{2}}d_{p,j}{\bf\omega}_{p,j}\right)\times \boldsymbol{n}_{ij},\end{eqnarray}$$

and can be decomposed into the normal and tangential components,

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}_{ij}^{n}=(\boldsymbol{u}_{ij}\boldsymbol{\cdot }\boldsymbol{n}_{ij})\boldsymbol{n}_{ij}, & \displaystyle\end{eqnarray}$$

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}_{ij}^{t}=\boldsymbol{u}_{ij}-\boldsymbol{u}_{ij}^{n}. & \displaystyle\end{eqnarray}$$

The normal force generated by the collision is modelled by considering the overlapping particles as a linear spring–damper system, with spring constant $k_{n}$ and damping constant ${\it\xi}_{n}$ :

(2.12) $$\begin{eqnarray}\boldsymbol{F}_{c,ij}^{n}=\left\{\begin{array}{@{}ll@{}}-k_{n}{\it\delta}_{ij}\boldsymbol{n}_{ij}-{\it\xi}_{n}\boldsymbol{u}_{ij}^{n}, & \text{ for }|\boldsymbol{x}_{p,i}-\boldsymbol{x}_{p,j}|<0.5(d_{p,i}+d_{p,j})+{\it\alpha},\\ 0, & \text{ otherwise}.\end{array}\right.\end{eqnarray}$$

In (2.12) the parameter ${\it\alpha}$ is the so-called radius of influence (Patankar & Joseph Reference Patankar and Joseph2001). Values of ${\it\alpha}$ greater than zero will initiate a repulsive force slightly before the particles overlap. This allows the model to be robust for high-speed collisions, but will result in close-packed states that are less dense than in real systems. To achieve correct close-packing densities, ${\it\alpha}$ is adjusted linearly as a function of the collision CFL number (Capecelatro & Desjardins Reference Capecelatro and Desjardins2013a ),

(2.13) $$\begin{eqnarray}\displaystyle {\it\alpha}=\frac{{\it\alpha}_{0}}{2}(d_{p,i}+d_{p,j})\left(\frac{\text{CFL}_{ij}^{c}}{\text{CFL}_{max}^{c}}\right),\quad \text{where }\text{CFL}_{ij}^{c}=\frac{2|\boldsymbol{u}_{ij,n}|{\rm\Delta}t_{p}}{(d_{p,i}+d_{p,j})}. & & \displaystyle\end{eqnarray}$$

Here, ${\rm\Delta}t_{p}$ is the particle time step and ${\it\alpha}_{0}$ is the maximum radius of influence, which is used when the collision CFL is equal to the maximum permitted collision CFL number, $\text{CFL}_{max}^{c}$ .

A similar linear spring–damper analogy with the addition of a frictional limiter is used to model the frictional forces generated by relative motion in the tangential direction as

(2.14) $$\begin{eqnarray}\boldsymbol{F}_{c,ij}^{t}=\left\{\begin{array}{@{}ll@{}}-k_{t}{\bf\delta}_{ij}^{t}-{\it\xi}_{t}\boldsymbol{u}_{ij}^{t}, & \text{for }\left|k_{t}{\bf\delta}_{ij}^{t}\right|\leqslant {\it\vartheta}_{s}\left|\boldsymbol{F}_{\!c,ij}^{n}\right|,\left|\boldsymbol{F}_{\!c,ij}^{n}\right|>0\text{ (sticking)},\\ -{\it\vartheta}_{s}\left|\boldsymbol{F}_{\!c,ij}^{n}\right|\boldsymbol{t}_{ij}, & \text{for }\left|k_{t}{\bf\delta}_{ij}^{t}\right|>{\it\vartheta}_{s}\left|\boldsymbol{F}_{\!c,ij}^{n}\right|,\left|\boldsymbol{F}_{\!c,ij}^{n}\right|>0\text{ (sliding)},\\ 0, & \text{for }|\boldsymbol{F}_{\!c,ij}^{n}|=0\text{ (no contact).}\end{array}\right.\end{eqnarray}$$

Here, ${\it\vartheta}_{s}$ is the coefficient of static/sliding friction, and the tangential unit vector is defined from the extension of the tangential spring, $\boldsymbol{t}_{ij}={\bf\delta}_{ij}^{t}/|{\bf\delta}_{ij}^{t}|$ . The model distinguishes between the sticking and sliding regimes based on the tangential displacement history, a feature crucial to obtaining stable heap formation (Brendel & Dippel Reference Brendel and Dippel1998). This requires storing and updating the length of the tangential spring, ${\bf\delta}_{ij}^{t}$ , during each time step according to

(2.15) $$\begin{eqnarray}{\bf\delta}_{ij}^{t}=\left\{\begin{array}{@{}ll@{}}{\bf\delta}_{ij}^{t}+\boldsymbol{u}_{ij}^{t}{\rm\Delta}t_{p} & \text{ (sticking)},\\ -{\displaystyle \frac{1}{k_{t}}}\left({\it\vartheta}_{s}|\boldsymbol{F}_{ij}^{n}|+{\it\xi}_{t}\boldsymbol{u}_{ij}^{t}\right) & \text{ (sliding)},\\ 0 & \text{ (no contact)}.\end{array}\right.\end{eqnarray}$$

Because the tangential plane may change during the duration of the contact, the length of the spring is projected back into the current tangential plane after every step as

(2.16) $$\begin{eqnarray}{\bf\delta}_{ij}^{t}={\bf\delta}_{ij}^{t}-({\bf\delta}_{ij}^{t}\boldsymbol{\cdot }\boldsymbol{n}_{ij})\boldsymbol{n}_{ij}.\end{eqnarray}$$

The total force acting on particle $i$ due to collision with all other particles $j$ is then

(2.17) $$\begin{eqnarray}\boldsymbol{F}_{c}=\mathop{\sum }_{j/=i}\boldsymbol{F}_{\!c,ij}^{n}+\boldsymbol{F}_{\!c,ij}^{t},\end{eqnarray}$$

and from Newton’s third law it follows that $\boldsymbol{F}_{\!c,ji}=-\boldsymbol{F}_{\!c,ij}$ . The corresponding collisional torque on particle $i$ is due to the tangential collision forces described above,

(2.18) $$\begin{eqnarray}\boldsymbol{T}_{c}=\mathop{\sum }_{j/=i}\frac{d_{i}}{2}\boldsymbol{F}_{\!c,ij}^{t}\times \boldsymbol{n}_{ij}.\end{eqnarray}$$

We also include a rolling torque to account for the increased rolling resistance of irregularly shaped particles, and use the directional constant torque model of Zhou et al. (Reference Zhou, Wright, Yang, Xu and Yu1999b ),

(2.19) $$\begin{eqnarray}\boldsymbol{T}_{r}=\mathop{\sum }_{j/=i}-{\it\vartheta}_{r}|\boldsymbol{F}_{\!c,ij}^{n}|r_{ij}\frac{{\bf\omega}_{rel}}{|{\bf\omega}_{rel}|}.\end{eqnarray}$$

Here, ${\it\vartheta}_{r}$ is a dimensionless coefficient of rolling friction, $r_{ij}=((d_{i}d_{j})/(d_{i}+d_{j}))/2$ is the reduced radius and ${\bf\omega}_{rel}={\bf\omega}_{p,i}-{\bf\omega}_{p,j}$ is used to form a unit vector in the direction of rolling resistance for particle $i$ .

2.2 Fluid motion

The discrete Lagrangian sediment motion is coupled to the continuous Eulerian fluid motion by solving the volume filtered Navier–Stokes equations (Anderson & Jackson Reference Anderson and Jackson1967; Cihonski et al. Reference Cihonski, Finn and Apte2013):

(2.20) $$\begin{eqnarray}\displaystyle \frac{\partial }{\partial t}({\it\rho}_{f}{\it\Theta}_{f})+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}_{f}{\it\Theta}_{f}\boldsymbol{u}_{f})=0, & & \displaystyle\end{eqnarray}$$

(2.21) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\partial }{\partial t}({\it\rho}_{f}{\it\Theta}_{f}\boldsymbol{u}_{f})+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}_{f}{\it\Theta}_{f}\boldsymbol{u}_{f}\boldsymbol{u}_{f})\nonumber\\ \displaystyle & & \displaystyle \qquad =-\boldsymbol{{\rm\nabla}}p+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left({\it\mu}_{eff}\left(\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{f}+\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{f}^{T}-\frac{2}{3}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}_{f}\right)\right)-{\it\Theta}_{f}{\it\rho}_{f}\boldsymbol{g}+\boldsymbol{f}_{p\rightarrow f}.\end{eqnarray}$$

In the above form the conservation equations account for the volume of fluid that is locally displaced by the motion of the sediment through the fluid fraction, ${\it\Theta}_{f}$ . In addition to volume exclusion effects, (2.21) also contains the typical interphase momentum transfer term, $\boldsymbol{f}_{p\rightarrow f}$ . This term includes the equal and opposite reaction from the particle surface forces back to the flow. To define these quantities, it is necessary to introduce the filter operation which distributes a property of the particles, ${\it\phi}_{p}$ , to the continuous field, ${\it\phi}_{f}$ , located at the Eulerian grid points,

(2.22) $$\begin{eqnarray}{\it\phi}_{f}(\boldsymbol{x})=\displaystyle \mathop{\sum }_{i_{p}=1}^{n_{p}}\mathscr{F}(\boldsymbol{x},\boldsymbol{x}_{p}){\it\phi}_{p}.\end{eqnarray}$$

The truncated polynomial filter function, $\mathscr{F}$ , suggested by Deen, van Sint Annaland & Kuipers (Reference Deen, van Sint Annaland and Kuipers2004) is used,

(2.23) $$\begin{eqnarray}\mathscr{F}(\boldsymbol{x},\boldsymbol{x}_{p})=\left\{\begin{array}{@{}ll@{}}{\displaystyle \frac{15}{16}}\left[{\displaystyle \frac{\left|\boldsymbol{x}-\boldsymbol{x}_{p}\right|^{4}}{h^{5}}}-2{\displaystyle \frac{\left|\boldsymbol{x}-\boldsymbol{x}_{p}\right|^{2}}{h^{3}}}+{\displaystyle \frac{1}{h}}\right], & \left|\boldsymbol{x}-\boldsymbol{x}_{p}\right|<h,\\ 0, & \text{otherwise},\end{array}\right.\end{eqnarray}$$

where $h$ is the filter half-width which controls the region of influence of each particle.

Using (2.22) and (2.23), the fluid fraction of a control volume with volume $\mathscr{V}_{cv}$ is defined from the particle positions as

(2.24) $$\begin{eqnarray}{\it\Theta}_{f}(\boldsymbol{x})=1-\frac{1}{\mathscr{V}_{cv}}\left[\displaystyle \mathop{\sum }_{i_{p}=1}^{n_{p}}\mathscr{F}(\boldsymbol{x},\boldsymbol{x}_{p})\frac{{\rm\pi}}{6}d_{p}^{3}\right].\end{eqnarray}$$

Similarly, the interphase momentum transfer term is

(2.25) $$\begin{eqnarray}\boldsymbol{f}_{p\rightarrow f}\left(\boldsymbol{x}\right)=-\mathop{\sum }_{i_{p}=1}^{n_{p}}\mathscr{F}\left(\boldsymbol{x},\boldsymbol{x}_{p}\right)(\boldsymbol{F}_{pr}+\boldsymbol{F}_{d}+\boldsymbol{F}_{l}+\boldsymbol{F}_{am}).\end{eqnarray}$$

The variation of the effective viscosity of the sediment–fluid suspension is modelled using a modified version of the Eilers (Reference Eilers1941) equation (Ferrini et al. Reference Ferrini, Ercolani, De Cindio, Nicodemo, Nicolais and Ranaudo1979) as

(2.26) $$\begin{eqnarray}{\it\mu}_{eff}={\it\mu}_{0}\left[1+\frac{0.5\left[{\it\mu}\right]{\it\Theta}_{p}}{\left(1-{\it\Theta}_{p}\right)/{\it\Theta}_{cp}}\right]^{2}.\end{eqnarray}$$

Here, ${\it\Theta}_{cp}$ is the packing fraction under close-packed conditions and $\left[{\it\mu}\right]$ is the intrinsic viscosity which accounts for the effect of particle shape on the rheology of the mixture. Aside from particle shape, factors such as microstructural changes with shear rate, and particle polydispersivity are known to influence ${\it\mu}_{eff}$ in complex ways (Stickel & Powell Reference Stickel and Powell2005). Equation (2.26) does not explicitly account for these factors, but has been shown to closely reproduce the behaviour of a variety of particle–fluid mixtures at high shear rates (Ferrini et al. Reference Ferrini, Ercolani, De Cindio, Nicodemo, Nicolais and Ranaudo1979) relevant to the present conditions.

The unresolved subgrid-scale turbulent stresses can be modelled using a large-eddy simulation approach, thereby adding a turbulent contribution to ${\it\mu}_{eff}$ . For the high- $Re$ oscillatory flow conditions considered here, a variable density dynamic Smagorinsky model described in Shams et al. (Reference Shams, Finn and Apte2011) is used. Under the present conditions, where the flow is resolved at near particle scale, the turbulent contribution to the effective viscosity is found to be minor relative to the contribution from the suspension viscosity (equation (2.26)).

The particles therefore have a direct influence on the resolved continuous phase stress tensor via the fluid fraction (volume displacement) and the interphase momentum transfer. These resolved fields in turn influence the modelled subgrid dissipation via the dynamic Smagorinsky and Eilers equation closures for ${\it\mu}_{eff}$ . We have not attempted to further model any subparticle fluid motions (wakes) introduced by the particles, or any direct effect that the particles have on the subgrid scales. Similarly, the effects of subgrid-scale fluid motions on particle dynamics are neglected, and a scaling analysis has indicated (§ 4.2) that the particles in the present study interact primarily with resolved fluid length scales.

3.1 Model parameters

The particle collision model contains several parameters, which can potentially influence both the micro- and the macroscale particle behaviour. In addition, the modified Eilers equation for the effective mixture viscosity contains the parameters ${\it\Theta}_{cp}$ and $\left[{\it\mu}\right]$ , and the transfer of information from particles to the flow relies on the filter scale parameter, $h$ . The parameter values used for input to the natural sand simulations in § 4 are provided in table 2 and are discussed below. Model sensitivity to reasonable variations in these parameters is explored in appendices A and B.

Table 2. Baseline parameters chosen for simulation of natural sands.

From a numerical standpoint, the main challenge is to select an appropriate value for the normal particle stiffness, $k_{n}$ . It is possible to derive values for $k_{n}$ and $k_{t}$ directly from the Young’s and shear moduli of a material, but such values can result in prohibitively small contact durations, which must be resolved by the simulation particle time step, ${\rm\Delta}t_{p}$ . Instead, $k_{n}=10^{3}~\text{N}~\text{m}^{-1}$ is chosen on the basis that it is large enough to avoid significant particle overlap ( ${>}1\,\%$ ). Arguably, the bulk particle motion will be rather insensitive to variations in stiffness in the presence of large dissipative forces (i.e. drag), and this has been observed to be true in prior numerical studies (i.e. Drake & Calantoni Reference Drake and Calantoni2001) as well as in appendix B. The ratio of tangential stiffness to normal stiffness is fixed to $k_{t}/k_{n}=0.4$ .

For most materials the dry coefficient of normal restitution, $e_{n,dry}$ , is well characterized (approximately 0.65 for quartz sands), and can be used to relate the normal damping parameter, ${\it\xi}_{n}$ , to the spring constant, $k_{n}$ , through an analytical solution to the linear mass–spring–damper problem,

(3.1) $$\begin{eqnarray}{\it\xi}_{n}=\frac{-2\ln e_{n}\sqrt{m_{ij}k_{n}}}{\sqrt{{\rm\pi}^{2}+\ln ^{2}e_{n}}},\end{eqnarray}$$

where $e_{n}$ is taken to be $e_{n,dry}$ and $m_{ij}=m_{i}m_{j}/(m_{i}+m_{j})$ is the reduced mass. During the rebound phase of an immersed (wet) collision, the flow in the gap between two particles may result in a significantly suppressed rebound velocity and lower effective coefficient of restitution, $e_{n,wet}$ . Point-particle models that do not resolve this flow should model the effects of these unresolved lubrication forces. Controlled experiments (Schmeeckle et al. Reference Schmeeckle, Nelson, Pitlick and Bennett2001; Yang & Hunt Reference Yang and Hunt2006) and FRS (Simeonov & Calantoni Reference Simeonov and Calantoni2012) have demonstrated the strong dependence of $e_{n,wet}$ on the impact Stokes number,

(3.2) $$\begin{eqnarray}St_{i}=\frac{m_{ij}\left|\boldsymbol{u}_{ij}^{n}\right|}{6{\rm\pi}{\it\mu}_{0}r_{ij}^{2}},\end{eqnarray}$$

where $St_{i}$ is the ratio of the inertial force required to stop the two approaching particles to the pressure force generated by the flow in the small gap between them. Results show there is a limit, $St_{c}$ , below which no rebound occurs ( $e_{n,wet}=0$ ) and an elastic limit, $St_{e}$ , above which the rebound is unaffected by lubrication forces ( $e_{n,wet}=e_{n,dry}$ ). Experiments with perfect spheres have shown that $St_{c}\approx 10{-}15$ , and that $e_{wet}$ asymptotes to $e_{dry}$ around a value of $St_{e}\approx 100{-}500$ (Yang & Hunt Reference Yang and Hunt2006). A much more limited characterization of immersed natural sand grain collisions exists, with the observations of Schmeeckle et al. (Reference Schmeeckle, Nelson, Pitlick and Bennett2001) being an exception. Although there is significant scatter in this dataset, they estimate a significantly smaller transition region with $St_{c}\approx 39$ and $St_{e}\approx 105$ due to both hydrodynamic and grain shape effects. These values are used in the present work, and it is assumed that $e_{n,wet}$ varies linearly with $St$ for impacts in the transition range ( $St_{c}<St<St_{e}$ ),

(3.3) $$\begin{eqnarray}e_{n,wet}=\min \left[\max \left[\frac{St_{i}-St_{c}}{St_{e}-St_{c}},0\right],1\right]e_{n,dry}.\end{eqnarray}$$

Similarly to the tangential spring history, $e_{n,wet}$ is computed during the initiation of contact between two particles and is stored as a property of the collision pair, $ij$ , for the duration of each contact. It is then used in (3.1) to compute ${\it\xi}_{n}$ .

Tangential damping associated with ${\it\xi}_{t}$ is set to zero for the present simulations. Lubrication influences on the tangential collision forces may be included, for example by including Stokesian lubrication theory for spherical particles (Marzougui, Chareyre & Chauchat Reference Marzougui, Chareyre and Chauchat2015). However, such effects are poorly quantified for irregular shaped particles, and in the present model the lubrication effects on the mixture shear stresses are included implicitly in the modified Eilers equation.

Values for the friction coefficients ${\it\vartheta}_{s}$ and ${\it\vartheta}_{r}$ may be determined experimentally for spherical particles, but the frictional interactions between angular sand grains are more difficult to characterize. Ideally, the model should be able to match the collective behaviour of angular sand particles, especially the close-packing fraction and angle of repose, while retaining a spherical representation for numerical convenience. To accomplish this, the influence of ${\it\vartheta}_{s}$ and ${\it\vartheta}_{r}$ on these two parameters is studied in a calibration test detailed in § A.2. This leads to the choice of ${\it\vartheta}_{s}=0.4$ and ${\it\vartheta}_{r}=0.06$ .

To compute the radius of influence, ${\it\alpha}_{0}=0.075$ and $\text{CFL}_{MAX}^{c}=0.1$ are chosen, in line with Patankar & Joseph (Reference Patankar and Joseph2001).

The maximum packing fraction used in the modified Eilers equation (2.26) is set to ${\it\Theta}_{cp}=0.64$ . For monodisperse spherical particles the intrinsic viscosity, $\left[{\it\mu}\right]$ , in such a model is accepted to be 2.5 (Stickel & Powell Reference Stickel and Powell2005), and the behaviour of non-spherical particle suspensions is compatible with increased values of $\left[{\it\mu}\right]$ (Ferrini et al. Reference Ferrini, Ercolani, De Cindio, Nicodemo, Nicolais and Ranaudo1979). Following the recommendations of Penko, Calantoni & Slinn (Reference Penko, Calantoni and Slinn2009), $\left[{\it\mu}\right]=3.0$ is chosen for natural sands in this paper.

The choice of the filter half-width, $h$ , is not straightforward (Capecelatro & Desjardins Reference Capecelatro and Desjardins2013a ; Simeonov, Bateman & Calantoni Reference Simeonov, Bateman and Calantoni2015), especially for polydisperse systems. It should be large enough that the volume filtered description of the flow is valid. On the other hand, it is desirable to keep it small enough that particle interactions with meso-scale flow features can be captured. A sensitivity analysis of the effects of $h$ (appendix B) led to the selection of $h=3d_{50}$ for the present study.

Figure 2. Strategy for isolating mobile grains on the surface of bedforms during the simulation.

3.2 Particle time advancement

For the submillimetre particle sizes considered here, collision durations can be quite small relative to the required fluid time step, and subcycling of the particle equations (2.1)–(2.3) is performed with ${\rm\Delta}t_{p}\ll {\it\tau}_{c}$ , where

(3.4) $$\begin{eqnarray}{\it\tau}_{c}=\sqrt{\frac{{\rm\pi}^{2}+\left(\ln e_{dry}\right)^{2}}{k_{n}/m_{min}}}\end{eqnarray}$$

is the minimum collision duration, with $m_{min}$ corresponding to the mass of the smallest particle in the system. Simulations of purely granular systems often employ high-order schemes for accurate particle time advancement, but in the presence of large dissipative forces (drag), first-order schemes perform adequately with low computational overhead (Van der Hoef et al. Reference Van der Hoef, van Sint Annaland, Deen and Kuipers2008; Shams et al. Reference Shams, Finn and Apte2011), and first-order Euler time advancement is used here.

Simulation of bedforms with a particle based approach can be very computationally intensive due to the severe separation of space and time scales: sand vortex ripples are typically tens of centimetres in length, meaning that several million particles can be needed to simulate both bedload and suspended load over a single ripple. When combined with the short collision durations, solution of the particle equations can quickly become a simulation bottleneck. To address this issue and enable the model to simulate full-scale bedforms, the model performs an extra step to identify which grains may be mobilized by the flow at any given instant. The concept is similar in spirit to the ‘bottom to top reconstruction’ methods sometimes used to simulate static granular assemblies (Pöschel & Schwager Reference Pöschel and Schwager2005) and is illustrated in figure 2. At regular intervals, the location of the three-dimensional bedform surface is computed by finding the level set ${\it\Theta}_{p}=0.5$ . Each particle is then classified based on its depth, $\mathscr{D}$ , below this interface as mobile ( $\mathscr{D}<9d_{50}$ ), fixed ( $9d_{50}<\mathscr{D}<18d_{50}$ ) or dormant ( $18d_{50}<\mathscr{D}$ ). Particle motion is then updated according to this classification as follows.

1. (i) Mobile: the particle is mobile, and will be advanced according to system (2.1), using the time step ${\rm\Delta}t_{p}\ll {\it\tau}_{c}$ .

2. (ii) Fixed: all particle forces are computed with the same ${\rm\Delta}t_{p}$ as for mobile particles, but $\boldsymbol{u}_{p}$ and ${\bf\omega}_{p}$ are set to zero.

3. (iii) Dormant: $\boldsymbol{u}_{p}$ and ${\bf\omega}_{p}$ are set to zero, and only hydrodynamic forces are computed, with a larger time step equal to the fluid phase time step, ${\rm\Delta}t_{f}$ (no subcycling).

In this approach, the fluid phase solution remains unchanged, and small velocities in and out of the fixed and dormant portions of the ripple bed may result, where the solid fraction is roughly ${\it\Theta}_{p}\approx 0.6$ . The purpose of the fixed particle layer is to support the mobile layer and isolate these particles from the dormant collisionless region below. The thickness of the mobile layer is chosen to be large enough that particle motion at the interface with the fixed particles is effectively zero, and the fixed layer thickness is chosen to ensure isolation of the mobile and dormant particles. If the number of dormant particles at any instant is large, the strategy described here becomes very efficient without sacrificing the ability to capture particle–fluid coupling and associated morphological change (i.e. ripple migration).

4.1 Simulation set-up

Table 3 provides an overview of the parameters used to set up these two simulations, and the numerical configuration is shown in figure 3. In both experiments, the piston velocity, $U(t)$ , corresponded to the near-bed flow under a second-order Stokes wave,

(4.1) $$\begin{eqnarray}U(t)=u_{1}\cos ({\it\omega}t-{\it\gamma})+u_{2}\cos (2{\it\omega}t-2{\it\gamma}).\end{eqnarray}$$

Here, ${\it\omega}=2{\rm\pi}/T$ is the angular frequency, $T$ is the period, $U_{1}$ and $U_{2}$ are the first and second harmonic velocity amplitudes, and ${\it\gamma}$ is the phase shift such that $U(0)=0$ ,

(4.2) $$\begin{eqnarray}{\it\gamma}=\arccos \left(\frac{\sqrt{U_{1}^{2}+8U_{2}^{2}}-U_{1}}{4U_{2}}\right).\end{eqnarray}$$

The piston velocity $U(t)$ is shown in figure 3(a) for both cases, and also above several phase-dependent results in this section. The values of $U_{1}$ and $U_{2}$ were set so that the maximum free-stream velocity, $U_{max}$ , was equal to $0.63~\text{m}~\text{s}^{-1}$ in the vortex ripple condition and $1.50~\text{m}~\text{s}^{-1}$ in the sheet-flow condition. Both cases have the same period, $T$ , similar velocity asymmetry, $a=(U_{max}/U_{max}-U_{min})\approx 0.6$ , and similar moments of maximum onshore velocity ( $t/T\approx 0.22$ ), maximum offshore velocity ( $t/T\approx 0.72$ ) and flow reversal ( $t/T\approx 0.44$ ). The root-mean-square orbital velocity, $U_{rms}=\sqrt{0.5U_{1}^{2}+0.5U_{2}^{2}}$ , provides orbital amplitudes, $A=\sqrt{2}U_{rms}/{\it\omega}$ , of 0.44 m and 1.0 m for vortex ripple and sheet-flow conditions respectively. In the following results, positive values of $U(t)$ are considered to be ‘onshore’ directed and negative values ‘offshore’ directed.

Figure 3. The set-up of the mobile bed oscillatory flow simulations. (a) The free-stream velocity used to force the flow. (c,d) The domain and initial particle configurations for the sheet-flow and vortex ripple conditions respectively. (b) The particle-size distributions used.

Table 3. The parameters used for the oscillatory flow simulations.

The domains used in both simulations (shown in figure 3 c,d) were periodic in both the streamwise ( $X$ ) and spanwise ( $Z$ ) directions. For the sheet-flow case, the domain should be large enough in these directions that the velocity autocorrelation decays within half the domain length, while in the wall normal ( $Y$ ) direction the domain must be sufficiently tall to eliminate damping or generation of turbulence by the upper boundary. Salon, Armenio & Crise (Reference Salon, Armenio and Crise2007) found that $L_{x}\gtrsim 50{\it\delta}_{s}$ , $L_{x}\gtrsim 40{\it\delta}_{s}$ and $L_{z}\gtrsim 25{\it\delta}_{s}$ were sufficient to satisfy these conditions for smooth-wall oscillatory flow simulations, where ${\it\delta}_{s}=\sqrt{{\it\nu}T/{\rm\pi}}$ is the Stokes layer thickness. Taking into account the added thickness and roughness of our mobile bed, the domain for the sheet-flow condition was chosen to be $L_{x}=72{\it\delta}_{s},L_{y}=72{\it\delta}_{s},L_{z}=36{\it\delta}_{s}$ . The domain for the rippled bed was made to be $L_{x}=1{\it\lambda}_{r},L_{y}=6h_{r},L_{z}=30{\it\delta}_{s}$ , where ${\it\lambda}_{r}=0.41~\text{m}$ and $h_{r}=0.076~\text{m}$ are the experimentally measured ripple wavelength and amplitude. For any periodic rippled bed simulation, $L_{x}$ is constrained to be an integer multiple of the ripple wavelength, and in the present paper only a single ripple was simulated under the assumption that ripple-to-ripple variations in the flow do not significantly influence the sediment behaviour. The spanwise thickness, $L_{z}=30{\it\delta}_{s}$ , was chosen to capture the near-bed three-dimensional hydrodynamics, but because of the periodic boundary condition, forced the ripple to remain approximately uniform in the spanwise direction, and restricted any morphological change to the $X$ – $Y$ plane. These are seen as reasonable first approximations for applying the model to rippled beds in light of the fact that the experimentally measured ripples were very uniform, repeatable and two-dimensional in the flow tunnel (Van der Werf et al. Reference Van der Werf, Doucette, O’Donoghue and Ribberink2007). The upper slip boundary was located $L_{y}\approx 6h_{r}$ above the ripple crest, corresponding to roughly the top of the flow tunnel. Both domains had an impenetrable bottom wall, made rough by fixing any particles that came into contact with it. The $Y=0$ level was considered to be the initial bed level for the sheet-flow condition and the ripple crest for the vortex ripple condition.

The sand used in both experiments was well sorted, with medium grain sizes. A lognormal particle-size distribution has been assumed, with geometric standard deviation of $\langle {\it\sigma}_{g}\rangle =1.46$ for both conditions and median grain sizes, $\langle d_{50}\rangle$ , equal to 0.28 mm for the sheet-flow condition and 0.44 mm for the vortex ripple condition. The $\langle ~\rangle$ brackets here are used to distinguish these global properties from the local space–time-dependent particle-size distributions examined later on. The terminal settling velocity of the median particle size is then $v_{t,50}=58~\text{mm}~\text{s}^{-1}$ for the vortex ripple condition and $34~\text{mm}~\text{s}^{-1}$ for the sheet-flow condition. Initial bed conditions were generated using the following procedure. First, the lognormal cumulative distribution functions (CDFs) shown in figure 3(b) were sampled to obtain the diameters of individual particles. Then, particles with these diameters were seeded in random non-overlapping positions in the lower portion of each domain. Totals of $N_{p}=3.8M$ and $16.8M$ particles were employed for the sheet-flow and vortex ripple simulations respectively. It was confirmed that the $\langle d_{10}\rangle$ , $\langle d_{30}\rangle$ , $\langle d_{50}\rangle$ , $\langle d_{70}\rangle$ and $\langle d_{90}\rangle$ values of the particle populations closely matched the values reported in the experiments. The sand particles were allowed to settle under the action of gravity in water for several seconds until particle motion ceased. Finally, in the vortex ripple condition, particles above the experimentally measured ripple surface were trimmed away to match the periodic ripple shape provided by Van der Werf et al. (Reference Van der Werf, Magar, Malarkey, Guizien and ODonoghue2008). Due to the asymmetry of the flow, this initial profile, shown in figure 3(d), is also asymmetric, with a steeper ‘lee’ slope facing onshore and a shallower ‘stoss’ slope facing offshore.

The continuous phase grid was uniform and cubic in the near-bed region, with a spacing of ${\it\Delta}=0.71~\text{mm}$ ( $1.6d_{50}$ ) for the vortex ripple condition and ${\it\Delta}=0.31~\text{mm}$ ( $1.1d_{50}$ ) for the sheet-flow condition. To reduce the total number of control volumes somewhat, the vortex ripple grid was stretched in the vertical direction for $Y>2.5h_{r}$ , to a maximum spacing of ${\it\Delta}_{y,max}=2.1~\text{mm}$ at the upper slip boundary. The resulting grids contain a total of 12M (sheet-flow) and 16M (vortex ripple) control volumes.

After the initial beds were created, the simulation was started from rest, and a time-dependent body force, representing the streamwise oscillatory pressure gradient, was applied to both the fluid and the particles,

(4.3) $$\begin{eqnarray}F_{x}(t)={\it\rho}_{f}\frac{\text{d}U(t)}{\text{d}t}.\end{eqnarray}$$

A total of four cycles were simulated for the sheet-flow condition, and nine cycles for the vortex ripple condition. To minimize any morphological change due to startup effects in the vortex ripple case, all particles were held fixed for the first two cycles while the flow developed.

4.2 Characteristic parameters

The transport regime for these two flows can be characterized using two non-dimensional parameters, namely the Shields parameter, ${\it\theta}$ , and the suspension number, $\mathscr{S}$ . For ${\it\theta}\gtrsim 0.5$ , bedforms become washed out and the transport occurs in the sheet-flow regime, while the value of $\mathscr{S}$ indicates the relative importance of bedload and suspended load to the overall transport. For conditions in the sheet-flow regime, a transition occurs from bedload dominated transport to suspension dominated transport when $\mathscr{S}$ becomes less than approximately 0.8–1.0 (Sumer et al. Reference Sumer, Kozakiewicz, Fredsøe and Deigaard1996). In oscillatory flows, it is typical to compute the Shields parameter as

(4.4) $$\begin{eqnarray}{\it\theta}=\frac{{\it\sigma}_{max}}{\left({\it\rho}_{p}-{\it\rho}_{f}\right)gd_{50}},\end{eqnarray}$$

where ${\it\sigma}_{max}$ is the maximum bed shear stress and can be estimated from the grain roughness friction factor, i.e. ${\it\sigma}_{max}=f_{w}U_{max}^{2}/2$ . The corresponding suspension number,

(4.5) $$\begin{eqnarray}\mathscr{S}=\frac{v_{t,50}}{u_{\star }},\end{eqnarray}$$

is the ratio of the particle settling velocity (of the mean size fraction) to the bed friction velocity, $u_{\star }=\sqrt{{\it\sigma}_{max}/{\it\rho}_{f}}$ . Using the expression of Swart (Reference Swart1974) for $f_{w}$ with a roughness height of $2.5d_{50}$ results in $[{\it\theta}_{VR},{\it\theta}_{SF}]=\left[0.37,2.24\right]$ and $[\mathscr{S}_{VR},\mathscr{S}_{SF}]=[1.12,0.34]$ . The transport is therefore expected to be suspension dominated for the sheet-flow condition, while a suspension number near 1 and relatively large Shields parameter indicate an energetic near-bed flow with significant suspended sediment for the vortex ripple condition.

Figure 4. Normalized particle length and time scales for the vortex ripple condition (—) and the sheet-flow condition (○). In (a), the particle scales are normalized by the Kolmogorov scales (top) and the integral scales (bottom). In (b), ${\it\tau}_{p}/{\it\tau}_{{\it\Delta}}$ is the ratio of the particle time scale to the time scale of the largest resolved eddy, and $d_{p}/l_{i}$ is the ratio of the particle size to the eddy size that has the same time scale as the particle, and therefore has a dominant effect on the relative particle velocity.

To understand the nature of the particle–turbulence interactions, the scaling analysis of Balachandar (Reference Balachandar2009) has been adapted to the present conditions. In figure 4(a), the ranges of the particle length and time scales ( $d_{p},{\it\tau}_{p}$ ) have been made non-dimensional by the approximate Kolmogorov scales ( ${\it\eta},{\it\tau}_{{\it\eta}}$ ) and integral scales ( $L,{\it\tau}_{L}$ ). Here, the particle time scale is computed as

(4.6) $$\begin{eqnarray}{\it\tau}_{p}=\frac{2{\it\rho}_{p}/{\it\rho}_{f}+1}{36}\frac{d_{p}^{2}}{{\it\nu}^{2}}\frac{1}{C_{d}(Re_{p},0)}.\end{eqnarray}$$

The Kolmogorov scales can be related to the Shields parameter by assuming that ${\it\eta}\approx {\it\nu}/u_{\star }$ , where ${\it\nu}={\it\mu}_{0}/{\it\rho}_{f}$ is the kinematic viscosity, resulting in

(4.7a,b ) $$\begin{eqnarray}{\it\eta}=\frac{{\it\nu}}{\sqrt{{\it\theta}\left({\displaystyle \frac{{\it\rho}_{p}}{{\it\rho}_{f}}}-1\right)gd_{50}}},\quad {\it\tau}_{{\it\eta}}=\frac{{\it\nu}}{{\it\theta}\left({\displaystyle \frac{{\it\rho}_{p}}{{\it\rho}_{f}}}-1\right)gd_{50}}.\end{eqnarray}$$

The integral time scale is estimated as ${\it\tau}_{L}=L/U_{rms}$ , where the integral length scale, $L$ , is taken to be either the ripple height, $h_{r}=0.076~\text{m}$ , or the maximum sheet-flow boundary layer thickness, $H_{bl}\approx 0.02~\text{m}$ (cf. figure 7). From figure 4(a), it can be seen that ${\it\eta}\ll d_{p}\ll L$ and also that ${\it\tau}_{{\it\eta}}\ll {\it\tau}_{p}\ll {\it\tau}_{L}$ for the range of particle sizes used in both conditions. The particles are between 10 and 50 times larger than ${\it\eta}$ , with corresponding Stokes numbers in the range $9\lesssim St\lesssim 34$ .

The relative velocity of the particles and the resulting momentum exchange between phases are therefore controlled primarily by an intermediate-size eddy in the inertial range, which has the same time scale as the $d_{p}$ size particle, in addition to the gravitational settling. The size of this eddy is $l_{i}={\it\tau}_{p}^{3/2}{\it\epsilon}^{1/2}$ , where ${\it\epsilon}\approx {\it\nu}^{3}/{\it\eta}^{4}$ is the energy dissipation rate. When adopting a point-particle LES approach, it is important to capture the effect of the $l_{i}$ sized eddies on the particles, meaning that their time scale should be larger than the time scale of the largest resolved eddy, ${\it\tau}_{{\it\Delta}}={\it\tau}_{{\it\eta}}({\it\Delta}/{\it\eta})^{2/3}$ (Balachandar Reference Balachandar2009). In figure 4(b), the ratio ${\it\tau}_{p}/{\it\tau}_{{\it\Delta}}$ is plotted versus the range of the non-dimensional particle size, $d_{p}/l_{i}$ . For the present conditions, the $l_{i}$ size eddies are between approximately 2.5 and 5 times larger than the particles, and these eddies are resolved by the computational grids ( ${\it\tau}_{p}/{\it\tau}_{{\it\Delta}}\gtrsim 1$ ). However, we note that especially for the smallest particle fractions the $l_{i}$ scale eddies may not be resolved by more than approximately two grid points, and the contribution of subgrid-scale eddies to the particle relative velocity may not be negligible. Our future efforts will involve incorporation of a stochastic model (i.e. Pozorski & Apte Reference Pozorski and Apte2009) to account for and investigate these effects.

The ratio $g/a_{k}$ , where $a_{k}={\it\epsilon}^{2/3}/{\it\eta}^{1/3}$ is the Kolmogorov acceleration, is much less than 1 for the present conditions, indicating that turbulence influences the relative particle velocity much more than gravitational settling for the present conditions. The particle Reynolds numbers can then be estimated from the expression provided by Balachandar (Reference Balachandar2009),

(4.8) $$\begin{eqnarray}Re_{p}=\sqrt{\frac{1}{2{\it\rho}_{p}/{\it\rho}_{f}+1}}\frac{{\it\rho}_{p}/{\it\rho}_{f}-1}{3\sqrt{C_{d}(Re_{p},0)}}\left(\frac{d_{p}}{{\it\eta}}\right)^{2},\end{eqnarray}$$

which results in $17\leqslant Re_{p}\leqslant 124$ for the vortex ripple condition and $25\leqslant Re_{p}\leqslant 175$ for the sheet-flow condition.

4.3 Experimental validation

To compare with the phase resolved experimental measurements, the simulation results have been averaged over the homogeneous directions ( $X$ and $Z$ for the sheet-flow condition, $Z$ for the vortex ripple condition) and have been ensemble averaged over the final three cycles of each simulation.

4.3.1 Near-bed velocities

Figure 5 shows the spanwise and ensemble averaged near-ripple fluid velocity vector fields at eight phases of the flow for the vortex ripple condition. The major features of the flow are consistent with the particle image velocimetry (PIV) results reported by Van der Werf et al. (Reference Van der Werf, Doucette, O’Donoghue and Ribberink2007) (not shown), which will briefly be summarized here. At the moments of flow reversal (a,e), a large vortex is ejected from the ripple trough. The vortex ejected over the lee slope at onshore reversal (e) is significantly stronger than the one ejected over the stoss slope during offshore reversal (a). This, combined with strong erosion of the ripple crest during the onshore half-cycle (b,c,d), contributes to the asymmetric shape and onshore migration of the ripples. The overshoot of the near-bed flow velocity can be seen clearly during the acceleration and deceleration phases (b and f) above the ripple crest.

Figure 5. Spanwise and ensemble averaged 2D vector fields shown at selected phases for the vortex ripple case.

A direct comparison with near-ripple experimental measurements (PIV) is made in figure 6 for eight points along the ripple. The data were extracted from a region $9~\text{mm}$ above the simulated ripple surface (determined as the surface where ${\it\Theta}_{p}=0.6$ ), corresponding as closely as possible to the experiments. Overall, quite good agreement is achieved at most positions along the ripple. The maximum on- and offshore velocities are very well predicted by the model, both above the crest and in the trough, indicating that vortex creation and ejection is well captured. The main disagreement can be found over the lee slope, shown in (a,c), during the creation of the stoss side vortex ( $0.6\lesssim t/T\lesssim 0.9$ ), which is possibly due to slight differences in the simulated and experimental ripple profiles. Some of the higher-frequency variation observed in the simulation results, but not the experimental measurements, could be due to the periodic simulation domain. Periodicity in the streamwise and spanwise directions is likely to promote a more identical cycle-to-cycle flow than in the laboratory, and this will preserve the signature of smaller-time-scale motions in the mean flow. Extending the ensemble average to the last five cycles (to match with the experiments) did not change the character of this result.

Figure 6. Comparison of predicted and measured velocities for the vortex ripple condition: —, simulation; ○, experiment (PIV); - - -, free-stream velocity; $x/{\it\lambda}=-0.09$ (a), 0.04 (b), $-0.22$ (c), 0.16 (d), $-0.35$ (e), 0.29 (f), $-0.48$ (g), 0.42 (h).

Figure 7 compares the measured and predicted velocity profiles at several phases for the sheet-flow condition. The simulations capture the velocity skewness of the boundary layer, and the agreement with the experimental measurements is quite good for all phases with the exception of the onshore acceleration ( $t/T=0.12$ ) phase, which lags behind the measured profile. It is worth noting that due to three-dimensional effects in the flow tunnel, the measured free-stream velocities are somewhat lower than the piston velocity used to drive the flow, and therefore exact hydrodynamic equivalence with the experiments is not achieved by forcing the flow using (4.1). However, the computed flow velocity profile follows the measured data very well from the fully packed bed up to the top of the suspension layer ( $Y\approx 40~\text{mm}$ ), which indicates the good performance of the model across regions with a wide range of sediment concentrations.

Figure 7. Comparison of predicted and measured velocities for the sheet-flow condition: —, simulation; ○, experiment.

Figure 8. Suspended sediment concentrations for the vortex ripple condition. The figure on the left-hand side shows the time and ripple averaged suspended sediment concentration as a function of height: —, simulation; ○, experiment (acoustic); ▾, experiment (optical); - - -, a fit to (4.9). The contour plots on the right-hand side show the ensemble averaged suspended sediment concentration at select phases ( $\log _{10}$ scale).

4.3.2 Suspended sediment concentrations

The time and ripple-length ( $X$ ) averaged suspended sediment concentrations predicted by the model are compared with the measurements on the left-hand side of figure 8 for the vortex ripple condition. The simulated concentration results generally lie between the acoustic and the optical measurements and illustrate the well-known negative exponential decay above the ripple crest (Bosman & Steetzel Reference Bosman and Steetzel1986; Nielsen Reference Nielsen1992; Ribberink & Al-Salem Reference Ribberink and Al-Salem1994), of the form

(4.9) $$\begin{eqnarray}C(z)=C_{0}\exp \left(\frac{-Yv_{t,50}}{{\it\epsilon}_{s}}\right).\end{eqnarray}$$

Here, $C_{0}$ is the concentration just above the ripple crest and ${\it\epsilon}_{s}$ is the sediment diffusivity. A fit to (4.9) with $C_{0}=2.37~\text{g}~\text{l}^{-1}$ and ${\it\epsilon}_{s}=0.0034~\text{m}^{2}~\text{s}^{-1}$ provides a good match to the data and is shown as a dashed line in figure 8. This suggests a vertical length scale of suspended sediment mixing of $r_{c}={\it\epsilon}_{s}/v_{t,50}=0.059~\text{m}$ , or approximately $0.8h_{r}$ , which agrees well with the data of Ribberink & Al-Salem (Reference Ribberink and Al-Salem1994). Van der Werf et al. (Reference Van der Werf, Magar, Malarkey, Guizien and ODonoghue2008) and Malarkey et al. (Reference Malarkey, Magar and Davies2015) show a similar comparison based on two-dimensional advection–diffusion based models. The much larger discrepancies in their results were attributed to an inadequate representation of turbulent sediment diffusion. A strength of the present model is that the time-dependent concentration field depends only on Lagrangian particle motion, with no uncertainties introduced to the concentration field through an advection–diffusion solution, reference concentration or pickup function. The results also suggest that in addition to the turbulence, particle-size effects may also play an important role in determining the cycle averaged concentration, as examined in the following section.

In figure 8(a–h), the spanwise and ensemble averaged suspended sediment concentration fields are shown at several phases of the cycle: at off–onshore flow reversal (a), a suspended sediment cloud with moderate concentrations can be seen over the stoss side slope, while a smaller region of high concentration is seen in the trough. This sediment is carried over one or two ripple lengths as the flow accelerates onshore, with a certain amount of sediment settling back to the bed during this time (b,c). At the same time, the ripple crest is eroded by the strong onshore flow (b,c,d). As the flow begins to decelerate (d), there is abundant sediment in suspension along the lee side slope due to the lee vortex generation. The lee vortex is therefore very rich in suspended sediment when it is ejected around the time of on–offshore flow reversal, resulting in a large suspended sediment cloud over the ripple crest (e,f). This cloud passes through the periodic domain and over the crest three times at roughly $t/T=0.5$ , $t/T=0.7$ and $t/T=0.9$ , consistent with experimental observations.

Figure 9. Sediment concentrations for the sheet-flow case. On the left-hand side, concentrations are shown at select phases in the sheet-flow layer: —, simulation; ○, experiment (CCM); $t/T=0$ (a), 0.12 (b), 0.22 (c), 0.34 (d), 0.44 (e), 0.54 (f), 0.74 (g), 0.92 (h). On the right-hand side, the period averaged concentration is shown with additional suction system measurements (▴) and the fit to (4.10) for $Y>1~\text{cm}$ (- - -) with ${\it\alpha}=2.1$ .

The concentration profiles in the lower sheet-flow layer are shown in figure 9 for the sheet-flow condition, and are compared with experimental measurements made using a conductivity concentration meter (CCM). The concentration measurements have been normalized by the undisturbed bed concentration ( $c_{bed}\approx 1600~\text{g}~\text{l}^{-1}$ ). The computed sediment concentration follows the CCM measurements very well in both the pickup layer ( $Y<-2~\text{mm}$ ) and the upper sheet-flow layer ( $-2~\text{mm}<Y<2~\text{mm}$ ). The main discrepancy seems to be around the initial deceleration towards on–offshore flow reversal (c,d), which suggests that the model underpredicts the pickup of sediment from the lower part of the bed at the maximum flow. It is possible that the effects of unresolved subgrid-scale fluid–particle interactions are important during this phase of high shear stress, which warrants future investigation of stochastic models for particle–fluid interaction under such conditions. The overall agreement, however, is considered to be good, which indicates the strength of the model in simulating particle–particle and fluid–particle interactions in such high-concentration regions.

In the plot on the right-hand side of figure 9, the period averaged concentration profile is compared with measurements made using a CCM and a suction sampling device by O’Donoghue & Wright (Reference O’Donoghue and Wright2004b ) over a larger range. The simulation is in reasonable agreement with the lower suction measurements and the upper CCM measurements. Also shown is the fit to the power law equation proposed by Ribberink & Al-Salem (Reference Ribberink and Al-Salem1994),

(4.10) $$\begin{eqnarray}C(z)=C_{0}\left(\frac{Y_{ref}}{Y}\right)^{{\it\alpha}},\end{eqnarray}$$

where $C_{0}$ is taken to be the reference concentration at 1 cm above the undisturbed bed level, and the exponent ${\it\alpha}\approx 2.1$ provided a good fit to their data for slightly finer sand than considered here ( $d_{50}=0.21~\text{mm}$ ). In the simulation, less sediment is retained in suspension for $Y>1~\text{cm}$ than was measured by the suction system, and a stronger decay is observed, with ${\it\alpha}\approx 3.3$ .

4.4 Three-dimensional description

Up to this point, streamwise and spanwise averaging have been used to reduce these flows to either one or two spatial dimensions (in addition to time). However, at the Reynolds numbers involved in the experiments, the turbulent structure of the flow and its interaction with the sediment are known to be highly three-dimensional (Ozdemir et al. Reference Ozdemir, Hsu and Balachandar2010; Penko et al. Reference Penko, Calantoni, Rodriguez-Abudo, Foster and Slinn2013). Of particular importance are the three-dimensional vortical structures which provide the regions of upward directed flow necessary for sediment to travel away from the bed and remain in suspension. The swirling strength criterion ( ${\it\lambda}_{ci}$ ) of Zhou et al. (Reference Zhou, Adrian, Balachandar and Kendall1999a ) is used to extract these features and examine their interaction with the particle phase by plotting isosurfaces of ${\it\lambda}_{ci}$ along with instantaneous particle positions in figures 10 and 11. To remove very-small-scale features from the visualizations, a low-pass filter has been applied to the instantaneous velocity fields before computing $\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{f}$ and its eigenvalues. To allow spatio-temporal size sorting to be identified, individual particles have been coloured and rendered according to their size.

Figure 10. Three-dimensional snapshots of vortex cores visualized as isosurfaces of ${\it\lambda}_{ci}=20~\text{s}^{-1}$ and particle positions at eight phases of the oscillatory flow for the vortex ripple condition. Individual particles have been rendered at twice their actual size to aid in visualization; $t/T=0$ (a), 0.12 (b), 0.22 (c), 0.34 (d), 0.44 (e), 0.54 (f), 0.74 (g), 0.92 (h).

Figure 11. As in figure 10, but for the sheet-flow condition, with an isovalue of ${\it\lambda}_{ci}=70~\text{s}^{-1}$ . Particles are rendered at their actual size; $t/T=0$ (a), 0.12 (b), 0.22 (c), 0.34 (d), 0.44 (e), 0.54 (f), 0.74 (g), 0.92 (h).

Figure 10 shows these results for the vortex ripple condition at eight phases of the wave cycle as in the previous section. The largest vortex cores can be seen to be pseudo-two-dimensional, spanning the entire width of the domain in the $Z$ direction. However, significant three-dimensional small-scale structures are also evident throughout the cycle in the form of entangled vortex filaments around the larger two-dimensional vortex cores and streaks along the bed surface, similar to observations from other numerical studies (e.g. Vittori & Verzicco Reference Vittori and Verzicco1998; Scandura, Vittori & Blondeaux Reference Scandura, Vittori and Blondeaux2000; Zedler & Street Reference Zedler and Street2006; Schmeeckle Reference Schmeeckle2014). These coherent vortical features have a clear influence on particle motion throughout the cycle: at off–onshore flow reversal (figure 10 a), a small vortex is detected at the bottom of the stoss slope in an area of high suspended sediment concentration (see also figure 8 a), also observed by Van der Werf et al. (Reference Van der Werf, Doucette, O’Donoghue and Ribberink2007). At the same time, remnants of several vortices ejected from the ripple trough during previous cycles can be seen higher in the water column, maintaining large numbers of particles around them. As the flow speed increases in the onshore direction (b), near-bed vortical streaks appear around the crest, and as the flow separates over the lee slope, a number of particles are picked up by these near-bed features. As the flow speed peaks (c), the 3D streaks cover the entire surface of the ripple, and some have rolled up into larger pseudo-two-dimensional vortex cores. Once the flow starts to reduce speed and reverse direction (d,e), these streaks are lifted from the bed and in the trough they roll up around a large two-dimensional vortex core, lifting a substantial number of particles from the bed. This particle cloud extends as high as $100~\text{mm}$ above crest level as the lee vortex is ejected during flow reversal (e) and the flow accelerates offshore (f). During the offshore half-cycle (e–h), a similar process takes place, but the near-bed streaks have a longer time to develop and become more elongated along the ripple surface (g,h) due to the asymmetry of the wave.

Comparatively, the same results for the sheet-flow condition, shown in figure 11, exhibit a lack of large-scale vortex shedding, and more uniform near-bed streaks that are significantly stronger (note that the magnitude of the ${\it\lambda}_{ci}$ isosurfaces is 3.5 times larger than in figure 10). Similarly to the vortex ripple condition, the flow is almost clear of vortical structures at the beginning of the cycle, apart from remnants of broken vortex cores from the last wave cycle. During the onshore flow acceleration (figure 11 b,c), densely located energetic vortical structures are generated close to the bed surface. These streaks remain close to the bed until maximum onshore flow (c), and as the flow decelerates (d), they are lifted 10–20 mm above the initial bed level, carrying large quantities of particles with them. At on–offshore flow reversal (e), most of the near-bed streaks have been dissipated, although many of them remain higher in the water column. During the subsequent offshore half-cycle (f–h), a similar process follows and, similarly to the vortex ripple condition, the near-bed streaks become better developed and elongated due to the asymmetry of the mean flow.

4.5 Behaviour of different size fractions

By rendering the particles according to size in figures 10 and 11, strong spatio-temporal particle-size sorting patterns are revealed, both near the bed and in suspension, which appear to be strongly influenced by the three-dimensional flow field.

4.5.1 Near the bed

Throughout most of the wave cycle in the vortex ripple condition, the bed surface is covered by large particles, and a similarly coarse surface layer can be seen in the low-velocity phases under sheet-flow conditions. The composition of this near-bed surface layer is important, as it directly influences which particles are brought into suspension and at what phase of the cycle, thus influencing the net transport behaviour (Dohmen-Janssen et al. Reference Dohmen-Janssen, Kroekenstoel, Hassan and Ribberink2002). This feature is explored in figure 12, where the relative abundance of different particle sizes in the near-bed layer, $Y^{\prime }=Y_{bed}+2.5d_{50}$ , has been plotted versus phase. To extract this result from the particle data, the instantaneous bed level, $Y_{bed}(x,z)$ , is taken to be the local elevation where ${\it\Theta}_{p}=0.6$ , and particles are binned into five fractions with mean diameters $d_{10}$ , $d_{30}$ , $d_{50}$ , $d_{70}$ and $d_{90}$ . In a perfectly mixed condition, each of these fractions would be expected to contribute 20 % to the total composition. The departure of the results in figure 12 from this behaviour demonstrates the effectiveness of different flow-driven particle-size sorting mechanisms.

Figure 12. Relative abundance of different particle-size fractions in the near-bed layer, $Y^{\prime }=Y_{bed}+2.5d_{50}$ : (- $\cdot$ -)  $d_{10}$ , (- - -) $d_{30}$ , (—) $d_{50}$ , (—▴—) $d_{70}$ , (—▪—) $d_{90}$ ; vortex ripple (a,c), sheet flow (b,d).

Under vortex ripple conditions, the composition of the near-bed layer is nearly constant for the low-velocity offshore half-cycle, and is disproportionately coarse with roughly 2/3 of the solid volume being coarser than $d_{60}$ . During onshore acceleration, when the near-bed vortical streaks touch the bed surface (figure 10 b,c), high-speed flow near the bed is able to lift the coarse fractions into suspension. This exposes the finer fractions below, especially along the stoss side of the ripple, which experiences some of the largest near-bed velocities (see figure 6). Consequently, the particle sizes in the near-bed layer are better mixed during the onshore half-cycle and during formation of the lee vortex. When the lee vortex is ejected and passes over the crest (figure 10 f,g), the largest fractions quickly settle down, coarsening the surface layer. Being exposed to the high velocities and vortical streaks for most of the cycle, the crest tends to lose coarse particles, while the trough area is less affected by the streaks and hence accumulates the coarse fractions. Such a sorting process agrees with time averaged measurements of the settling velocity distribution (see Van der Werf et al. (Reference Van der Werf, Magar, Malarkey, Guizien and ODonoghue2008) figure 4). The near-constant composition in the near-bed layer reflects the fact that mixing and particle-size sorting largely occur in the suspension under vortex ripple conditions, driven by large-scale vortex shedding.

Under sheet-flow conditions, the overall composition of the near-bed layer is much more dynamic. A strong coarse-over-fine layering of particles is seen in low-velocity phases (figure 11 a,e), which results in a very coarse near-bed composition. As the flow accelerates onshore, these coarse particles are lifted to a low level by the near-bed streaks (figure 11 b), and the layer composition changes from very coarse to very fine as the finer fractions below are mobilized around maximum onshore velocity (figure 11 c). Once these fine fractions become exposed to the flow, they are rapidly brought into suspension by the near-bed structures being lifted from the bed as the flow decelerates and reverses (figure 11 d,e). At the point of flow reversal, the coarsest fractions have settled back to the bed and the near-bed layer is again predominantly coarse. while many fine particles remain in suspension. During the offshore half cycle (figure 11 f–h), the flow strength is low enough that the coarse particles remain close to the bed and only a small number of fine particles are lifted into near-bed suspension.

4.5.2 In suspension

Once sediment is lifted from the bed, the three-dimensional vortical structures provide the vertical velocities needed to retain particles in suspension. A net upward flux of particles can be achieved by a skewed vertical velocity probability (Bagnold Reference Bagnold1966; Wei & Willmarth Reference Wei and Willmarth1991) combined with a loitering tendency (Tooby, Wick & Isaacs Reference Tooby, Wick and Isaacs1977; Nielsen Reference Nielsen1992), which causes settling particles to spend more time sampling the upward parts of the flow against their downward settling preference. In figure 13, the probability density function of the vertical fluid velocity experienced by suspended sediment particles ( ${\it\Theta}_{p|p}<0.08$ ) is plotted for both conditions. The vertical fluid velocities are normalized by the clear water terminal velocity for five different size fractions, so that the ratio $v_{f|p}/v_{t}>1$ indicates fluid velocities resulting in upward particle suspension for each fraction, while $v_{f|p}/v_{t}<1$ results in downward particle settling (assuming that drag is the dominant hydrodynamic force on the particles). In both conditions, the probability distributions are positively skewed, and have modes shifted to positive values of $v_{f|p}$ . This shift in the mode decreases with particle size, and reflects the decreased potential of massive particles to become trapped in the vortical regions. In the vortex ripple condition, the mode of the vertical velocity distribution experienced by the smaller fractions corresponds to $v_{f|p}/v_{t}\approx 1$ , underscoring the effectiveness of vortex trapping as a mechanism to retain sediment in suspension. This leads to the visible gradient in suspended particle sizes that will be explored in more detail in the next section.

Figure 13. Probability density function of the vertical fluid velocity experienced by suspended particles ( ${\it\Theta}_{p|p}<0.08$ ) normalized by the clear water terminal settling velocity for five different size fractions in (a) the ripple bed condition and (b) the sheet-flow condition: (- $\cdot$ -)  $d_{10}$ , (- - -) $d_{30}$ , (—) $d_{50}$ , (—▴—) $d_{70}$ , (—▪—) $d_{90}$ .

Figure 14. Space–time visualization of particle-size distributions in the oscillatory boundary layer simulations. (a,b) The time dependence of the layer averaged concentration ( $\log _{10}$ scale). (c,d) The time dependence of the layer averaged $d_{50}$ normalized by $\langle d_{50}\rangle$ . (e,f) The time dependence of the layer averaged ${\it\sigma}_{g}$ . (g,h) The time dependence of the layer averaged settling velocity, $v_{s}$ , normalized by $v_{t,50}$ , the terminal velocity of an isolated particle with $d=d_{50}$ . The dashed line corresponds to the ripple crest (vortex ripple (a,c,e,g)) or the initial bed level (sheet flow (b,d,f,h)).

4.6 Spatio-temporal dynamics of the particle-size distribution

To better understand the grain size specific dynamics seen in §§ 4.4 and 4.5, the spatio-temporal dynamics of the particle-size distribution for these conditions are analysed in figure 14 as functions of the vertical elevation, $Y$ , and phase, $t/T$ : (a) and (b) present the sediment concentrations; (c) and (d) show the ratio $d_{50}/\langle d_{50}\rangle$ , which indicates the local mean particle size relative to the bed mean value, $\langle d_{50}\rangle$ ; (e) and (f) show the local ${\it\sigma}_{g}$ value, which indicates the degree of mixing of the grain sizes; (g) and (h) show the mean volumetric settling velocity, $v_{s}=u_{p,y}-u_{f|p,y}$ , normalized by $v_{t,50}$ , the clear water terminal velocity of a particle with $d_{p}=\langle d_{50}\rangle$ . The results were obtained by streamwise and spanwise averaging over all particles at a given $Y$ elevation. However, for the vortex ripple condition, only particles that have become entrained above the instantaneous bed level ( ${\it\Theta}_{p}<0.08$ ) are included in the averages, thus the results are representative of the suspended particles in the vortex ripple condition. In all plots, the thick dashed line corresponds to the ripple crest in the vortex ripple condition or the initial bed level in the sheet-flow condition.

Many of the suspension phenomena already discussed have clear signatures in the space–time concentration plots (figure 14 a,b). Significant sediment concentration, $1~\text{g}~\text{l}^{-1}$ or more, can be found under vortex ripple conditions up to as high as $120~\text{mm}$ ( $1.5h_{r}$ ) above crest level throughout most of the cycle, with concentrations of $0.1~\text{g}~\text{l}^{-1}$ persisting up to 250 mm above crest level ( $3h_{r}$ ). In contrast, the majority of suspended sediment in the sheet-flow condition is contained in the $30~\text{mm}$ above the initial bed level, and there is significant settling of suspended particles during the low-velocity phases. This again highlights the differences in sediment suspension mechanisms, and the effectiveness of the ripple generated turbulence at retaining sediment in suspension.

The parameters describing the particle-size distribution demonstrate strong space–time dependence in these flows. In the vortex ripple condition, particles suspended in the high-concentration regions below crest level are on average coarser than $\langle d_{50}\rangle$ , due to the continuous pickup and settling of the coarse fractions from the ripple surface. The mean suspended particle size above the crest line bears a resemblance to the concentration pattern: a significant rise in both concentration and $d_{50}$ immediately after flow reversal $(t/T=0.4)$ at levels up to $100~\text{mm}$ is maintained until before the peak offshore flow $(t/T=0.7)$ , due to the passage of the lee vortex sediment cloud. The lee vortex is capable of lifting all particle sizes from the bed into this region, where ${\it\sigma}_{g}\approx 1.35$ approaches the initial bed value (1.46), indicating that the sediment cloud is well mixed. At elevations more than $100~\text{mm}$ above the crest, the mean particle size is more or less constant at a given level throughout the cycle, with its value decreasing from $0.8\langle d_{50}\rangle$ at 100 mm to $0.6\langle d_{50}\rangle$ at $300~\text{mm}$ . In this range, ${\it\sigma}_{g}$ is also seen to decrease significantly, indicating that the sorting of particles by the settling velocity is very effective. This is confirmed by figure 14(g), which shows a nearly constant vertical gradient of the settling velocity above $100~\text{mm}$ .

The corresponding plots for the sheet-flow condition reveal a very pronounced size sorting within the sheet-flow layer. The most noticeable feature, which persists for the entire cycle and is also clearly visible in figure 11, is a well-established layer of finer particles below the initial bed level. This layer develops very quickly (it is visible after just one flow cycle) and is armoured by a layer of coarser particles above, consistent with the laboratory findings of Hassan & Ribberink (Reference Hassan and Ribberink2005). During onshore acceleration, when the bed is rapidly mobilized, the high-concentration sheet-flow layer is very coarse up to the point of maximum velocity. It subsequently becomes progressively finer as the coarse particles are peeled away and the finer fractions become exposed to the flow. By the time flow reverses, the mean diameter in suspension is less than $0.85\langle d_{50}\rangle$ , almost all the way down to the initial bed level. During the offshore half-cycle, the sheet-flow layer remains coarse because, as discussed in §§4.4 and 4.5, the near-bed flow is not energetic enough to expose the fine fractions below. During this time, a particle mixing layer (high ${\it\sigma}_{g}$ ) exists at the interface between the high- and low-concentration regions where coarse particles rising from the bed meet the fine particles left in suspension from the onshore half-cycle. The variations in the fall velocity under sheet-flow conditions reflect the functional dependence on particle size as well as local concentration and therefore vary significantly throughout the water column.

Figure 15. Vertical profiles of normalized local $d_{50}$ at four different instances for the vortex ripple case. Each profile corresponds to one of the five zones shown above: (—▫—) stoss trough, (—♢—) stoss slope, (—○—) crest, (—▿—) lee slope, (—▵—) lee trough; off-onshore flow reversal (a), max on-shore velocity (b), on-offshore flow reversal (c), max offshore velocity (d).

Within $100~\text{mm}$ of the crest in the vortex ripple condition, the temporal variation of $d_{50}$ and ${\it\sigma}_{g}$ is due to specific pickup, suspension and transport events occurring in different lateral zones along the ripple. To better understand these dynamics, the vertical variation of $d_{50}$ is plotted in figure 15 at the moments of off–onshore flow reversal (a), peak onshore velocity (b), on–offshore flow reversal (c) and peak offshore velocity (d) for each of five lateral zones: the stoss side trough (ST), the stoss side slope (SS), the crest (CR), the lee side slope (LS) and the lee side trough (LT). The variability of $d_{50}$ among different lateral zones is quite strong at the moments of flow reversal (a,c), when the energetic vortex ejection events entrain coarse particles from one side of the ripple surface, but on the opposite slope smaller particles are able to settle due to the low velocities. At these moments, the mean suspended particle size at a given vertical level in the ripple trough varies by almost 50 %. There is comparably less lateral variation during the peak velocity phases (b,d), when turbulent advection/diffusion processes are dominant and can efficiently mix particles of different sizes. At maximum onshore velocity (b), all vertical size profiles are similar, while at maximum offshore velocity, the signature of the coarse sediment cloud passing over the crest and lee slope can be seen between 50 and 100 mm above the crest.

4.7 Collisional dynamics

To model immersed particle impacts in the wave bottom boundary layer, it is common for both continuum and particle based approaches to use a single coefficient of normal restitution, $e_{n}$ . For irregular particles with non-uniform sizes and conditions that result in spatio-temporal grain size sorting, an obvious choice for this coefficient may not exist, which is why the impact Stokes number is used to compute $e_{n}$ on a per-collision basis. It is therefore of interest to examine the impact Stokes number statistics from the simulation data for these two cases. In figure 16, the CDF of the impact Stokes number (3.2) is plotted for both conditions. The different symbols correspond to collisions occurring in regions with different solid fractions. Also plotted as a dashed line is the $St=39$ threshold for viscously damped rebound observed by Schmeeckle et al. (Reference Schmeeckle, Nelson, Pitlick and Bennett2001) for quartz sands. For both cases, at very low and high solid fractions, nearly all collisions occur at very low Stokes numbers, ${\lesssim}20$ . This is for different reasons. Low-concentration regions ( ${\it\Theta}_{p}<0.02$ ) in these flows are higher in the water column and contain finer fractions in suspension. These fractions will have faster response times (low fall velocities), making the likelihood of an energetic collision due to large relative velocity between particles very small. On the other hand, at very high solid fractions in and near the bed ( ${\it\Theta}_{p}>0.3$ ), relative particle motion is inhibited by much lower settling velocities, even for large particles. At intermediate solid fractions in the range $0.02\leqslant {\it\Theta}_{p}\leqslant 0.3$ , occurring in the sheet-flow layer and the ripple trough, a broader range of $St$ is observed for both conditions. However, most impacts occur well below the critical value of $St=39$ . In the vortex ripple case, which has a larger $d_{50}$ , roughly 10 % of all collisions occur above this threshold.

Figure 16. Time averaged cumulative distribution functions of the collisional Stokes number (3.2) for (a) vortex ripple conditions and (b) sheet-flow conditions. Individual CDFs are computed for collisions occurring in regions of different solid fractions: (—▫—) ${\it\Theta}_{p}<0.02$ , (—♢—) $0.02\leqslant {\it\Theta}_{p}<0.08$ , (—○—) $0.08\leqslant {\it\Theta}_{p}<0.2$ , (—▹—) $0.2\leqslant {\it\Theta}_{p}<0.3$ , (—▵—) $0.3\leqslant {\it\Theta}_{p}<0.4$ . The dashed line corresponds to $St=39$ , the threshold of zero rebound for natural sand grains found by Schmeeckle et al. (Reference Schmeeckle, Nelson, Pitlick and Bennett2001).

These statistics indicate that the use of a small constant coefficient of normal restitution to mimic the viscous damping of low-Stokes-number impacts may be a safe assumption for small- to medium-size sand. For much larger particles, however, caution should be taken because the impact Stokes numbers will grow nonlinearly with particle size. Moreover, in light of the strong particle-size sorting effects observed in these flows, unphysically low restitution coefficients may result in a sink of energy for the larger particle fractions, or during energetic suspension events at particular phases of the flow.