Particle based modelling and simulation of natural sand dynamics in the wave bottom boundary layer

Abstract

Sand transport and morphological change occur in the wave bottom boundary layer due to sand particle interactions with an oscillatory flow and granular interactions between particles. Although these interactions depend strongly on the characteristics of the particle population, i.e. size and shape, little is known about how natural sand particles behave under oscillatory conditions and how variations in particle size influence transport behaviour. To enable this to be studied numerically, an Euler–Lagrange point-particle model is developed which can capture the individual and collective dynamics of subaqueous natural sand grains. Special treatments for particle collision, friction and hydrodynamic interactions are included to take into account the wide size and shape variations in natural sands. The model is used to simulate sand particle dynamics in two asymmetric oscillatory flow conditions corresponding to the vortex ripple experiments of Van der Werf et al. (J. Geophys. Res., vol. 112, 2007, F02012) and the sheet-flow experiments of O’Donoghue & Wright (Coast. Engng, vol. 50, 2004, pp. 117–138). A comparison of the phase resolved velocity and concentration fields shows overall excellent agreement between simulation and experiments. The particle based datasets are used to investigate the spatio-temporal dynamics of the particle-size distribution and the influence of three-dimensional vortical features on particle entrainment and suspension processes. For the first time, it is demonstrated that even for the relatively well-sorted medium-size sands considered here, the characteristics of the local grain size population exhibit significant space–time variation. Both conditions demonstrate a distinct coarse-over-fine armouring at the bed surface during low-velocity phases, which restricts the vertical mobility of finer fractions in the bed, and also results in strong pickup events involving disproportionately coarse fractions. The near-bed layer composition is seen to be very dynamic in the sheet-flow condition, while it remains coarse through most of the cycle in the vortex ripple condition. Particles in suspension spend more time sampling the upward directed parts of these flows, especially the smaller fractions, which delays particle settling and enhances the vertical size sorting of grains in suspension. For the submillimetre grain sizes considered, most particle–particle collisions occur at low impact Stokes numbers and can be expected to have low rebound velocities. The results are considered to be fundamental to the understanding of particle sorting and transport mechanisms, and can help to improve current transport modelling approaches.

1 Introduction

In a coastal environment, sandy particles on the seabed surface will start to move once the shear force exerted by surface waves and currents exceeds critical values. Such motion involves complex fluid–particle and particle–particle interactions and can result in various visible features such as vortex ripples in the low-flow regime, or a flat bed with intensive transport under large stormy waves, i.e. sheet flows (van Rijn 1993). Increasing evidence from theoretical, experimental and field studies suggests that the physical processes involved in these grain scale interactions often underpin the mechanisms for transport observed on natural beaches. In particular, the size and shape variations of natural sand particles have long been recognized as important factors influencing local sand transport under oscillatory flows. For example, Black & Oldman (1999) observed wave induced grain size sorting and subsequent effects on sand ripple development on the continental shelf, and Vincent, Stolk & Porter (1998) and Roos et al. (2007) identified spatial variations in the mean grain sizes on the Middelkerke Bank that have clear influences on the sand banks and seabed features along the Belgian coast.

Broadly, grain size influences on transport processes can be classified into two effects. First, for the suspended particles, vertical size sorting and selective transport can lead to very different transport capacity for different size fractions. Especially under oscillatory flows due to surface waves, the strong onshore peak flow under the wave crest can lift a large amount of sand from the bed surface. Once in suspension, the fine fractions take a significant time to settle back to the bed, during which the flow already changes its direction, resulting in enhanced offshore suspended transport, i.e. phase lag effects (Dohmen-Janssen et al. 2002). The overall net transport direction for fine sand therefore is often found against the wave propagation direction. In contrast, medium and coarse sands are transported close to the bed by the mean flow, settling during low-velocity phases, and their wave-period averaged net transport tends to be in line with the wave propagation direction (Camenen & Larson 2006; Hassan & Ribberink 2010; Kranenburg, Hsu & Ribberink 2014). Second, within the bed, agitation by the oscillatory shear force together with the bed level changes due to bedform migration or bed scour can lead to inverse size gradation (Bagnold 1954; Legros 2002), where coarse fractions rise to the surface, leaving the finer fractions below with very limited mobility. The most noticeable consequence of these armouring and sheltering effects is the shift in transport regime, i.e. a switch from offshore transport towards onshore transport or vice versa (Egiazaroff 1965).

Under oscillatory flows, particle exchange between the bedload and the suspended load is fairly dynamic through competing settling and entrainment effects (Hassan & Ribberink 2005). Therefore, in order to understand the grain size effects, the transport processes within both the suspension layer and the high-concentration bed region need to be examined hand in hand. Yet, most existing model concepts decouple these two regions in some way and largely rely on a single-particle representation of the population of natural sand sizes (van Rijn et al. 2013). This means that the coupled selective transport in the water column and the sorting/segregation phenomena in the bed are yet to be fully understood (Blondeaux 2012). Several series of oscillatory flow experiments have provided detailed measurements of net transport rates, flow hydrodynamics and sediment concentrations for both ‘uniform’ and graded sands in sheet-flow conditions (Ribberink & Al-Salem 1994; Ahmed & Sato 2003; O’Donoghue & Wright 2004a , b ; Hassan & Ribberink 2005) as well as above vortex ripples (Ribberink & Al-Salem 1994; Fredsøe, Andersen & Sumer 1999; Van der Werf et al. 2007). Due to difficulties involved in measuring individual particle accelerations and grain size while transport is taking place, information about particle dynamics, grain size sorting and the impact on transport is very limited. The recent laboratory measurements of Revil-Baudard et al. (2015) are perhaps the most detailed datasets to date of near-bed particle scale transport processes. Their time resolved velocity and concentration measurements for a steady sheet flow show a number of features compatible with existing modelling concepts, e.g. the Rouse (1937) profile, and strongly damped turbulence (Ribberink & Al-Salem 1994), but a strong intermittency observed in the near-bed transport emphasizes the importance of grain scale interactions with coherent flow structures.

These factors motivate the development of numerical modelling approaches to identify the mechanisms of particle–fluid interaction, grain sorting and selective transport with relation to the bottom boundary layer hydrodynamics. Balachandar (2009) and Balachandar & Eaton (2010) reviewed simulation techniques for turbulent multiphase flows, and identified the range of applicability for several methods in terms of particle length and time scales, $d_{p}$ and ${\it\tau}_{p}$ , non-dimensionalized by the Kolmogorov scales, ${\it\eta}$ and ${\it\tau}_{{\it\eta}}$ . When the particle Stokes number, $St={\it\tau}_{p}/{\it\tau}_{{\it\eta}}$ , is small ( $St\lesssim 0.2$ ), particles respond rapidly to the smallest eddies in the flow, and the relative particle phase velocity may be expressed using a first-order expansion in terms of $St$ (Ferry & Balachandar 2001). The resulting equilibrium Eulerian or single mixture equations then become a very effective tool for simulating transport of small particles, as demonstrated by Penko et al. (2013) in the vortex ripple regime and Ozdemir, Hsu & Balachandar (2010) for oscillatory sheet flows. For larger Stokes numbers, a continuum description of the particle phase can be retained by solving two distinct sets of conservation equations within a classical Eulerian–Eulerian approach (Druzhinin & Elghobashi 1998). Several one-dimensional two-fluid frameworks have been proposed for simulation of sheet flows, typically in conjunction with modified two-equation Reynolds averaged turbulence closures to account for particle–fluid interactions in various $St$ regimes (Hsu, Jenkins & Liu 2004; Amoudry, Hsu & Liu 2008; Chen et al. 2011; Revil-Baudard & Chauchat 2013; Kranenburg et al. 2014). Recently, Chen & Yu (2015) have extended their model to 2D, and simulated sediment motion over fixed vortex ripples. Their simulations, as well as the discrete-vortex model calculations of Van der Werf et al. (2008) and Malarkey, Magar & Davies (2015), reveal the effectiveness of sediment trapping by near-bed vortices, and the correlation of the observed suspended sediment concentration peaks with the vortex formation and ejection events. Polydisperse particles and particles with larger inertia ( $St\gtrsim 1$ ) are problematic to model with continuum approaches, although adaptations based on the method of moments (Fox & Vedula 2009) are promising. In such cases a Lagrangian representation of the sediment phase, where hydrodynamic interactions between the fluid and particles are computed on a per-particle basis, becomes an attractive option. In the context of multiphase turbulent flow simulation, Lagrangian particle models can roughly be classified at either fully resolved or point-particle approaches.

In fully resolved simulations (FRS), the no-slip boundary condition on individual particle surfaces is enforced directly, and all scales of fluid motion, including ones introduced by the particles, are resolved by the continuous phase grid. Impressive simulations have been performed by several groups on spherical particle bedload transport in steady currents (Derksen 2011; Ji et al. 2013; Kidanemariam & Uhlmann 2014b ), flows of dense suspensions (Picano, Breugem & Brandt 2015) and the formation and propagation of small rolling grain and vortex dunes (Kidanemariam & Uhlmann 2014a ). It is expected that FRS can contribute significantly to knowledge of microscale sediment transport processes, and lead to improved closures for the averaged equations (Prosperetti 2015). However, for simulating sediment transport at length and time scales relevant to the wave bottom boundary layer, FRS are likely to be limited by computational resources for some time to come.

In the point-particle approach, the particle–fluid interface is not resolved by the continuous phase grid. Instead, particles are treated as point sources of mass and momentum and are coupled to the flow through the appropriately averaged equations (i.e. Anderson & Jackson 1967) and closure relationships. This approach is now well established as a research tool in a number of fields (see Van der Hoef et al. (2008), Subramaniam (2013) for recent reviews). Early efforts at point-particle modelling in the wave bottom boundary layer concentrated on bedload transport of massive particles, where collisional forces are dominant, and one-dimensional hydrodynamic models could be justified for the near-bed flow (Jiang & Haff 1993; Drake & Calantoni 2001). More recent work has coupled point-particle models to three-dimensional direct and large-eddy simulations in order to simulate both bedload and suspended load particle transport in steady currents (Apte, Mahesh & Lundgren 2008; Capecelatro & Desjardins 2013b ; Schmeeckle 2014). For many geophysical flows, including full-scale coastal boundary layers, large Reynolds numbers result in prohibitively small grid sizes for FRS. At the same time, large $St$ complicates existing continuum based models. To simulate natural sand dynamics in relevant conditions, the point-particle approach may then be the only viable option, provided that the important particle–fluid interactions can be captured or modelled.

Within the sand bed and in the dense near-bed region of the bottom boundary layer, energy is exchanged and dissipated between particles during collisions, and a model is needed for these particle–particle interactions. Existing models based on the soft-sphere discrete element method (DEM) (Cundall & Strack 1979) typically assume a common spherical shape for all of the particles. When dealing with natural sands, several fundamental mechanisms at the particle scale are absent as a result of this simplification.

1. (i) The packing fraction of a natural sand bed at rest and its maximum angle of repose both depend strongly on the particle-size distribution (Soulsby 1997).

2. (ii) The effective drag coefficient on individual angular sand particles is known to be stronger than that of equivalent sized spheres (Fredsøe & Deigaard 1992), and there is evidence that the hindered settling velocity of a particle suspension may be influenced by particle shape (Baldock et al. 2004).

3. (iii) Immersed collisions of subaqueous sand particles appear to be more influenced by lubrication forces than collisions of spherical particles (Schmeeckle et al. 2001), resulting in lower coefficients of restitution for similarly energetic collisions.

Without resolving these critical challenges, application of the DEM so far is limited to particle-laden flows with fairly uniform particle size and shape. The present study therefore aims to tackle these difficulties in order to extend the capability of the DEM approach to simulate natural sediment transport processes. The inclusion of realistic size distributions is straightforward within a point-particle model, and lognormal size distributions characterized by a median diameter, $d_{50}$ , and geometric standard deviation, ${\it\sigma}_{g}$ , are used here. Explicit simulation of irregular shaped particles has been performed (e.g. Calantoni, Holland & Drake 2004), but in the present model it is preferred to retain the computational simplicity of spherical particles, and instead irregular shape effects are represented implicitly through several mechanisms in the model equations. The present model builds upon the point-particle approach described by Finn, Shams & Apte (2011), Shams, Finn & Apte (2011) and Cihonski, Finn & Apte (2013), extending that model to capture these phenomena in densely laden flows with natural sand grains. The coupled multiphase equations for Eulerian fluid motion and Lagrangian particle motion are described in § 2, and some unique features of their numerical implementation are described in § 3. In § 4, two mobile bed oscillatory flow simulations are described corresponding to the vortex ripple experiments of Van der Werf et al. (2007) and the sheet-flow experiments of O’Donoghue & Wright (2004a , b ). Simulation results are compared directly with measurements before exploring the size-dependent particle motion and spatio-temporal size sorting phenomena. Additionally, a systematic validation of the model’s ability to predict the individual and collective behaviour of sand particles, and a model sensitivity analysis are provided in appendices A and B.

2 Mathematical model

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 (1983) 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 1954). 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. (2015) for a recent review). For this work, the correlation of Tenneti, Garg & Subramaniam (2011) 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 1992) rather than the typical relation for smooth spherical particles (Schiller & Naumann 1935). 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 (2007a , 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 (1965) 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 (2001) 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 (2008) and Zastawny et al. (2012) 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 (1979) 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 2001). 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 2013a ),

(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 1998). 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. (1999b ),

(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 1967; Cihonski et al. 2013):

(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 (2004) 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 (1941) equation (Ferrini et al. 1979) 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 2005). 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. 1979) 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. (2011) 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 Numerical implementation

The two systems of equations for the particles (equations (2.1)–(2.3)) and fluid (equations (2.20)–(2.21)) are coupled and solved in a structured Cartesian grid framework using a finite volume discretization and a pressure based second-order fractional time-stepping scheme based on the work of Finn et al. (2011), Shams et al. (2011) and Cihonski et al. (2013). Extensive validation of the method for dilute bubble- and particle-laden flows has been reported in these works. In appendix A, several new tests are presented to validate and calibrate the model for simulating subaqueous natural sand behaviour, namely (i) immersed collision dynamics, (ii) macroscale packing fraction and angle of repose, and (iii) hindered settling behaviour at both dense and dilute concentrations. In appendix B, model sensitivity to the parameters described in this section is examined for the case of unidirectional flow.

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 2001) 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. 2001; Yang & Hunt 2006) and FRS (Simeonov & Calantoni 2012) 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 2006). A much more limited characterization of immersed natural sand grain collisions exists, with the observations of Schmeeckle et al. (2001) 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 2015). 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 (2001).

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 2005), and the behaviour of non-spherical particle suspensions is compatible with increased values of $\left[{\it\mu}\right]$ (Ferrini et al. 1979). Following the recommendations of Penko, Calantoni & Slinn (2009), $\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 2013a ; Simeonov, Bateman & Calantoni 2015), 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. 2008; Shams et al. 2011), 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 2005) 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 Model application: sand transport in oscillatory boundary layers

The remainder of this paper considers two simulations of sediment transport in oscillatory boundary layer flow, corresponding to condition Mr5b63 from the vortex ripple (VR) experiments of Van der Werf et al. (2007) and condition MA5010 from the sheet-flow (SF) experiments of O’Donoghue & Wright (2004a , b ). Both experiments were conducted in a closed-loop piston-driven oscillatory flow tunnel and resulted in phase resolved velocity and concentration measurements.

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 (2007) 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. 2007). 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. (2008). 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. 1996). 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 (1974) 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 (2009) 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 2009). 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 2009) 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 (2009),

(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. (2007) (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 1986; Nielsen 1992; Ribberink & Al-Salem 1994), 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 (1994). Van der Werf et al. (2008) and Malarkey et al. (2015) 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 (2004b ) 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 (1994),

(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.3.3 Net sediment transport

In the experiments, the net sediment transport rate was estimated by accounting for bed profile changes as well as the amounts of the sediment collected in traps at both ends of the tunnel. For the vortex ripple condition, Van der Werf et al. (2007) found a small offshore directed net transport, $q_{vr}=-3.69\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ , while the ripples were observed to migrate onshore with a velocity of $u_{mig}\approx 18~\text{mm}~\text{s}^{-1}$ . The simulated ripple migration speed, $u_{mig}\approx 14~\text{mm}~\text{s}^{-1}$ , is in very good agreement. However, the net sediment transport rate in the simulations is found to be onshore directed, with a magnitude of $q_{vr}=8.70\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ , suggesting that some of the offshore directed suspended sediment fluxes are underpredicted by the model. For the sheet-flow condition, O’Donoghue & Wright (2004b ) reported $q_{sf}=53\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ , which is in reasonable agreement with the simulated value of $q_{sf}=42.3\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ . Ribberink & Al-Salem (1994) found a strong correlation between the net transport and the third velocity moment for both rippled bed and sheet-flow conditions, $q=A\langle U^{3}\rangle$ . Using $A=150\times 10^{-6}$ , which corresponds to their flows with similar periods to the ones considered here, gives $q_{vr}=3.12$ and $q_{sf}=49.4$ , quite close to and in the same direction as both simulation results.

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. 2010; Penko et al. 2013). 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. (1999a ) 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 1998; Scandura, Vittori & Blondeaux 2000; Zedler & Street 2006; Schmeeckle 2014). 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. (2007). 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. 2002). 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. (2008) 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 1966; Wei & Willmarth 1991) combined with a loitering tendency (Tooby, Wick & Isaacs 1977; Nielsen 1992), 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 (2005). 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. (2001) 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. (2001).

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.

5 Conclusions

An Euler–Lagrange point-particle model has been developed to simulate the dynamics of subaqueous natural sand particles. Several characteristics of the model allow it to faithfully reproduce the well-known individual and collective behaviour of natural sand grains, specifically the following.

1. (i) A variable model for the coefficient of normal restitution is used to model immersed particle–particle impacts. This allows the model to capture the dynamics of lubricated collisions at low and transitional impact Stokes numbers, as well as particle rebound of more massive particles involved in energetic collisions (§ A.1).

2. (ii) A rolling friction model is combined with the typical sliding friction approach to capture the enhanced friction of angular natural sand grains. The frictional coefficients were calibrated to match the close-packing fraction and angle of repose of natural quartz sands (§ A.2).

3. (iii) The particle–fluid interactions, and crucially the drag law and the mixture viscosity model, implicitly take into account the effects of particle shape, particle concentration and the local particle-size distribution. This allows the model to reproduce the settling velocity behaviour of both uniform spheres and angular sands for solid fractions ranging from dilute to close packed (§ A.3).

The model was applied to simulate sand particle motion in asymmetric oscillatory flow conditions typical of the wave bottom boundary layer in both the vortex ripple and sheet-flow regimes. Conditions were matched with the medium sand flow tunnel experiments of O’Donoghue & Wright (2004a , b ) and Van der Werf et al. (2007), and the predictions of phase resolved velocity and concentration fields are in overall excellent agreement with the experimental measurements. The three-dimensional phase resolved data were used to understand particle motion in these two very different conditions, and to characterize the flow induced particle-size sorting. From these detailed particle based datasets the following conclusions can be drawn.

1. (1) Despite these flows being nominally one- or two-dimensional, the near-bed vortical structures responsible for sediment pickup, suspension and mixing are highly three-dimensional. In the vortex ripple condition, the dominant flow features, detected with the ${\it\lambda}_{ci}$ criterion, involve a pseudo-two-dimensional core and a mass of entangled three-dimensional vortex filaments, usually accompanied by large sediment clouds. In the sheet-flow condition, near-bed vortical streaks are created during each half-cycle, and bring significant amounts of sediment into suspension when they are lifted from the bed at flow reversal.

2. (2) Even for the fairly narrow particle-size distributions considered here, strong sorting of the suspended sediment by particle size is significant under both sheet-flow and vortex ripple conditions. The size sorting demonstrates significant space and time dependence due to the influence of specific suspension/settling events that lead to segregation and mixing of the grain size population.

3. (3) Both conditions show a strong inverse size gradation of particles in the surface layers, with a layer of the coarsest particles resting above the finer fractions at moments of low velocity. This armouring of the fine fractions restricts their vertical mobility and means that only the most energetic suspension events are able to suspend significant amounts of sediment due to the low mobility of the surface particles. This was quantified by examining the composition of the near-bed layer, which is very different under these two conditions. In the vortex ripple case, this layer remains predominantly coarse throughout most of the cycle. The strong near-bed vortical streaks and major disturbance of the bed by the lee vortex result in a mixing of all size fractions around maximum onshore velocity. In the sheet-flow condition, the near-bed layer composition is much more dynamic, and all fractions participate in the pickup events during different phases of the flow.

4. (4) There is a preference for suspended particles to sample upward velocity regions of the flow generated by three-dimensional vortices. This is due to both a positively skewed velocity probability in the flow and a tendency of the particles to loiter in regions of upward directed flow created by 3D vortices. In the vortex ripple conditions, the velocity of maximum probability is sufficient to keep the smaller fractions permanently suspended, enhancing the vertical size sorting of grains in suspension.

5. (5) For the medium-size sands considered, most particle–particle impacts occur at low impact Stokes numbers, where lubrication forces can be expected to result in near or complete damping of rebound velocities. Some more energetic collisions do occur, especially in high-concentration regions of the ripple trough, but it is expected that a low constant coefficient of restitution will yield satisfactory results when simulating fine to medium-size sand with a spring–damper type model for normal collision force.