Interface-resolved direct numerical simulations of sediment transport in a turbulent oscillatory boundary layer

Abstract

The flow within an oscillatory boundary layer, which approximates the flow generated by propagating sea waves of small amplitude close to the bottom, is simulated numerically by integrating the Navier–Stokes and continuity equations. The bottom is made up of spherical particles, free to move, which mimic sediment grains. The approach allows one to fully resolve the flow around the particles and to evaluate the forces and torques that the fluid exerts on their surface. Then, the dynamics of sediments is explicitly computed by means of the Newton–Euler equations. For the smallest value of the flow Reynolds number presently simulated, the flow regime turns out to fall in the intermittently turbulent regime such that turbulence appears when the free-stream velocity is close to its largest value but the flow recovers a laminar-like behaviour during the remaining phases of the cycle. For the largest value of the Reynolds number, turbulence is significant during almost the whole flow cycle. The evaluation of the sediment transport rate allows one to estimate the reliability of the empirical predictors commonly used to estimate the amount of sediments transported by sea waves. For large values of the Shields parameter, the sediment flow rate during the accelerating phases does not differ from that observed during the decelerating phases. However, for relatively small values of the Shields parameter, the amount of moving particles depends not only on the bottom shear stress but also on flow acceleration. Moreover, the numerical results provide information on the role that turbulent eddies have on sediment dynamics.

1 Introduction

In nature, flows that involve the motion of solid particles coupled to that of a fluid are quite common, and different models have been developed to predict phenomena involving the motion of sediment particles in either air or water. The approaches employed to describe sediment and fluid motions are different depending on the spatial scale of interest, which can range from a few millimetres to hundreds of kilometres, i.e. from the scale of the sediment grains to the scale of the largest morphological patterns observed on the Earth’s surface (e.g. tidal sandbanks).

Depending on the problem under investigation, the interstitial fluid can play a minor role in the transport of momentum, and the rheology of the mixture is mainly controlled by phenomena occurring during direct grain–grain contacts. On the other hand, under different conditions, as happens in dilute suspensions, the motion of the fluid plays a primary role in the dynamics of the mixture. Finally, the hydrodynamic force acting on sediment grains and the force due to grain–grain contacts could be equally important, as happens at the bottom of water bodies (seas, lakes, rivers, estuaries, etc.), where flow drag can mobilise sediment grains arrested on the bed surface by gravity and frictional contacts.

The threshold conditions for the initiation of sediment transport and the sediment transport rate are usually determined by considering the average velocity field and neglecting the turbulent fluctuations (see e.g. Graf (1984), Fredsøe & Deigaard (1992), Soulsby (1997), Gyr & Hoyer (2006)). However, the vortex structures that characterise a turbulent flow might induce local high values of the fluid velocity and mobilise the sediment particles even when the average flow is relatively weak. Despite a lot of experimental studies having been devoted to investigate the mechanisms responsible for the initiation of sediment transport and the complex dynamics of sediment grains, a clear and detailed picture of the interaction of coherent vortex structures and sediment particles is still missing.

The flow generated by a monochromatic surface wave of small amplitude propagating over a flat sandy bottom provides a fair description of the actual flow that is observed in coastal environments seawards of the breaker zone. Close to the bottom, the surface wave induces oscillations of the pressure gradient and originates an oscillatory boundary layer (OBL). The OBL is characterised by (i) the amplitude $U_{0}^{\ast }$ of the irrotational velocity oscillations close to the bottom, (ii) the order of magnitude $\unicode[STIX]{x1D6FF}^{\ast }=\sqrt{2\unicode[STIX]{x1D708}^{\ast }/\unicode[STIX]{x1D714}^{\ast }}$ of the thickness of the viscous bottom boundary layer and (iii) the angular frequency $\unicode[STIX]{x1D714}^{\ast }=2\unicode[STIX]{x03C0}/T^{\ast }$ of the surface wave, where $T^{\ast }$ is the wave period. Hereinafter, an asterisk is used to denote a dimensional quantity, while the same symbol without the asterisk denotes its dimensionless counterpart. Moreover, we let the mechanical properties of sea water be assumed constant and represented by the density $\unicode[STIX]{x1D71A}^{\ast }$ and the kinematic viscosity  $\unicode[STIX]{x1D708}^{\ast }$ . The sediments are assumed to be cohesionless, monodisperse and characterised by the density $\unicode[STIX]{x1D71A}_{s}^{\ast }$ and the diameter $d^{\ast }$ of the grains.

The dynamics of the OBL over the seabed is rich because features typical of the laminar, transitional and turbulent regimes might coexist during a flow cycle depending on the values of the Reynolds number $R_{\unicode[STIX]{x1D6FF}}=U_{0}^{\ast }\unicode[STIX]{x1D6FF}^{\ast }/\unicode[STIX]{x1D708}^{\ast }$ and the dimensionless particle diameter $d=d^{\ast }/\unicode[STIX]{x1D6FF}^{\ast }$ .

The OBL over a smooth wall (i.e. for $d=0$ ) becomes turbulent if $R_{\unicode[STIX]{x1D6FF}}$ is larger than 550 (Costamagna, Vittori & Blondeaux 2003), while, in the presence of particles, this value progressively decreases, being approximately equal to 500 for $d=2.32$ , 400 for $d=2.80$ and 150 for $d=6.95$ (Ghodke & Apte 2016, 2018; Mazzuoli & Vittori 2016, 2019). Laboratory observations show that different flow regimes exist within the OBL over a smooth wall, namely, the laminar regime, the disturbed laminar regime, the intermittently turbulent regime and the fully developed turbulent regime. In the disturbed laminar regime, small perturbations of the Stokes flow appear but the average flow does not deviate significantly from that observed in the laminar regime. The intermittently turbulent regime is characterised by the appearance of turbulent bursts during the decelerating phases of the cycle but the flow recovers a laminar-like behaviour during the accelerating phases. Finally, in the fully developed turbulent regime, turbulence is present during the whole oscillation cycle (Hino, Sawamoto & Takasu 1976; Hino et al. 1983; Jensen, Sumer & Fredsøe 1989; Akhavan, Kamm & Shapiro 1991; Carstensen, Sumer & Fredsøe 2010). Later, the experimental observations found a theoretical interpretation by Wu (1992) and Blondeaux & Vittori (1994), who showed that the appearance of turbulence is due to both nonlinear three-dimensional effects and a receptivity mechanism. These theoretical findings were later supported by the results of direct numerical simulations (DNS) of Navier–Stokes and continuity equations (Akhavan et al. 1991; Verzicco & Vittori 1996; Costamagna et al. 2003; Ozdemir, Hsu & Balachandar 2014). It is worth pointing out that similar results were obtained by considering a rough wall, even though the roughness of the wall causes turbulence to appear for smaller values of the Reynolds number (Jensen et al. 1989; Carstensen, Sumer & Fredsøe 2012; Mazzuoli & Vittori 2019).

Even when, for small values of both $R_{\unicode[STIX]{x1D6FF}}$ and $d$ , the flow never becomes turbulent during the flow cycle, the prediction of the sediment transport rate is challenging because sediments are subject both to the viscous drag and to the effects of the wave-driven pressure gradient. Mazzuoli, Kidanemariam & Uhlmann (2019) investigated the formation of ripples in the OBL for $R_{\unicode[STIX]{x1D6FF}}=72$ and 128 over a bed of spherical mono-sized particles of dimensionless diameter $d$ equal to 0.25. By means of DNS, Mazzuoli et al. (2019) showed that the contribution of the pressure gradient to the particle dynamics can be significant. Indeed, they observed that a significant amount of sediment was mobilised also during phases characterised by small values of the bed shear stress. However, if the size $d^{\ast }$ of the sediment is noticeably smaller than $\unicode[STIX]{x1D6FF}^{\ast }$ , this contribution can be neglected and the sediment flow rate is well correlated with the bed shear stress, which is fairly well approximated by that of a Stokes boundary layer.

Laboratory experiments carried out by Lobkovsky et al. (2008), in the absence of turbulent fluctuations, revealed that sediment dynamics rapidly adapts to the slow changes of the driving flow, and the sediment flow rate could be estimated, also in unsteady conditions, by a power-law function of the excess of bed shear stress with respect to the critical value for incipient sediment motion. Mazzuoli et al. (2019) concluded that, as long as the appearance of turbulence is not triggered, both the viscous and pressure-gradient contributions were mainly controlled by the parameter $\unicode[STIX]{x1D6F9}/R_{\unicode[STIX]{x1D6FF}}$ , with $\unicode[STIX]{x1D6F9}=U_{0}^{\ast 2}/[((\unicode[STIX]{x1D71A}_{s}^{\ast }-\unicode[STIX]{x1D71A}^{\ast })/\unicode[STIX]{x1D71A}^{\ast })g^{\ast }d^{\ast }]$ denoting the mobility number and $g^{\ast }$ being the modulus of the gravitational acceleration.

For small values of $d$ , transition to turbulence occurs in the early stages of the decelerating phases, in a way apparently similar to that over a smooth wall (Mazzuoli & Vittori 2016). The effect of the transition to turbulence can be observed in figure 1, where the bottom shear stress measured by Jensen et al. (1989) at the bottom of an oscillatory boundary layer over a smooth bottom for $R_{\unicode[STIX]{x1D6FF}}=761$ is plotted versus the phase $\unicode[STIX]{x1D711}$ of the cycle. The phase variable, $\unicode[STIX]{x1D711}\in [\!0,2\unicode[STIX]{x03C0}\![$ , is expressed in radians and defined in order to be equal to zero when the maximum absolute value of the velocity far from the bottom is attained. Turbulence appears during the decelerating phases, but the flow recovers a laminar-like behaviour during the accelerating phases (intermittently turbulent regime). As discussed in Blondeaux, Vittori & Porcile (2018), it can be easily verified that the intermittently turbulent regime is present in a significant part of the coastal region, which shifts towards the shore during mild wave conditions whereas it shifts towards the offshore region during storms.

Then, by increasing either $R_{\unicode[STIX]{x1D6FF}}$ or $d$ , transition to turbulence occurs earlier and earlier, therefore pervading also the accelerating phase. It is noteworthy that, for $R_{\unicode[STIX]{x1D6FF}}>150$ and $d=6.95$ , the wall is hydrodynamically rough and turbulent fluctuations practically never disappear (Ghodke & Apte 2016; Mazzuoli & Vittori 2016) even though turbulence strength during the decelerating phases differs from that observed during the accelerating phases. In fact, one of the difficulties of modelling of the wave-averaged sediment transport induced by propagating surface waves lies in the fact that the sediment transport rate during the accelerating phases of the cycle differs from that observed during the decelerating phases even though the free-stream velocity has the same value. Although, in a large number of the empirical formulae used to quantify the sediment transport rate, the sediment flux is independent of the sign of the flow acceleration, figure 1 suggests that the bottom shear stress and the sediment transport rate observed during the decelerating phases are associated with levels of turbulence much larger than those observed during the accelerating phases. A further difficulty comes from the experimental evidence that the threshold conditions for the initiation of sediment motion differ from those leading particles to stop. In the former case, the probability that sediment particles are set into motion depends on the occurrence of a favourable particle–flow interaction, while in the latter situation, particles can stop moving depending on the likelihood that they find a stable configuration on the bed surface (Clark et al. 2017).

Figure 1. Dimensionless wall shear stress $|\unicode[STIX]{x1D70F}_{b}|=|\unicode[STIX]{x1D70F}_{b}^{\ast }|/(\frac{1}{2}\unicode[STIX]{x1D71A}^{\ast }U_{0}^{\ast }\unicode[STIX]{x1D714}^{\ast }\unicode[STIX]{x1D6FF}^{\ast })$ plotted versus the phase $\unicode[STIX]{x1D711}$ of the cycle for a smooth wall and $R_{\unicode[STIX]{x1D6FF}}=761$ . The dots are the experimental data by Jensen et al. (1989) and the continuous line is the dimensionless external velocity magnitude in the experiment, $|U_{e}|=|U_{e}^{\ast }|/U_{0}^{\ast }$ .

The results of the DNS, which are described in the following, are aimed at verifying that the picture we have drawn previously is realistic and significantly affects sediment dynamics. In particular, we want: (I) to verify whether the differences in the hydrodynamics of the boundary layer and the dynamics of sediment grains during the accelerating and decelerating phases of the wave-induced bottom flows are significant; and (II) to evaluate the dependence of the sediment flow rate $q_{s}^{\ast }$ on quantities characterising the flow properties, like the bottom shear stress $\unicode[STIX]{x1D70F}_{b}^{\ast }(t)$ or the turbulent kinetic energy, for values of the parameters such that the flow regime is intermittently turbulent (see table 1).

Table 1. Flow parameters for the present runs. From left to right: the Reynolds numbers $R_{\unicode[STIX]{x1D6FF}}=U_{0}^{\ast }\unicode[STIX]{x1D6FF}^{\ast }/\unicode[STIX]{x1D708}^{\ast }$ and $Re_{d}=U_{0}^{\ast }d^{\ast }/\unicode[STIX]{x1D708}^{\ast }$ ; the mobility number $\unicode[STIX]{x1D6F9}=U_{0}^{\ast 2}/v_{s}^{\ast 2}$ , with $v_{s}^{\ast }=\sqrt{(s-1)g^{\ast }d^{\ast }}$ indicating the sediment fall velocity; the Reynolds number of the sediment $Re_{p}=v_{s}^{\ast }d^{\ast }/\unicode[STIX]{x1D708}^{\ast }$ (also known as the Galileo number); the Keulegan–Carpenter number $K_{c}=U_{0}^{\ast }/(\unicode[STIX]{x1D714}^{\ast }d^{\ast })$ ; and the specific gravity $s=\unicode[STIX]{x1D71A}_{s}^{\ast }/\unicode[STIX]{x1D71A}^{\ast }$ . Note that the ratio $d^{\ast }/\unicode[STIX]{x1D6FF}^{\ast }=Re_{d}/R_{\unicode[STIX]{x1D6FF}}$ is equal to 0.335 for all the runs. Also shown are: the dimensionless displacement thickness $\unicode[STIX]{x1D6FF}_{dis}^{\ast }/\unicode[STIX]{x1D6FF}^{\ast }$ ; the grain Reynolds number $u_{\unicode[STIX]{x1D70F},max}^{\ast }d^{\ast }/\unicode[STIX]{x1D708}^{\ast }$ , with the maximum friction velocity defined as $u_{\unicode[STIX]{x1D70F},max}^{\ast }=\max _{t}\sqrt{|\unicode[STIX]{x1D70F}_{b}^{\ast }|/\unicode[STIX]{x1D71A}^{\ast }}$ ; and the maximum value of the Shields parameter $|\unicode[STIX]{x1D703}|$ . The last column gives the number of periods used in the post-processing, after the transient was removed.

The paper is structured as follows. In the next section we formulate the problem and we briefly describe the numerical approach used to evaluate the flow field and the sediment dynamics. In § 3, we describe the flow field and the sediment transport. Finally, § 4 is devoted to the conclusions.

2 Formulation of the problem and numerical approach

The flow within the boundary layer at the bottom of a sea wave is investigated by assuming that the wave steepness, i.e. the ratio between the amplitude and the length of the wave, is small and the linear Stokes theory describes the flow generated far from the bottom by wave propagation. Even though this approach neglects nonlinear effects and in particular the existence of a steady streaming, it provides a fair description of the oscillatory flow generated by propagating sea waves close to the bottom when their amplitude is small. Then, the flow within the bottom boundary layer can be determined by approximating it as the flow generated by an oscillating pressure gradient close to a fixed wall. Nonlinear effects, which become significant when the wave propagates into shallow water because of the increase of its amplitude, are neglected. In particular, the presence of steady streaming and wave asymmetry, which produce a skewness of the flow velocity and acceleration in the OBL (van der A et al. 2011; Scandura, Faraci & Foti 2016), are not presently considered. Hence the pressure gradient, which drives the flow, can be written in the form

(2.1a-c ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p^{\ast }}{\unicode[STIX]{x2202}x_{1}^{\ast }}=-\unicode[STIX]{x1D71A}^{\ast }U_{0}^{\ast }\unicode[STIX]{x1D714}^{\ast }\sin (\unicode[STIX]{x1D714}^{\ast }t^{\ast }),\quad \frac{\unicode[STIX]{x2202}p^{\ast }}{\unicode[STIX]{x2202}x_{2}^{\ast }}=0,\quad \frac{\unicode[STIX]{x2202}p^{\ast }}{\unicode[STIX]{x2202}x_{3}^{\ast }}=0,\end{eqnarray}$$

where $(x_{1}^{\ast },x_{2}^{\ast },x_{3}^{\ast })$ is a Cartesian coordinate system such that the $x_{1}^{\ast }$ -axis points in the direction of wave propagation and the $x_{2}^{\ast }$ -axis is vertical and points in the upward direction. The pressure gradient described by (2.1) drives the fluid motion as well as the motion of spherical particles of density $\unicode[STIX]{x1D71A}_{s}^{\ast }$ and diameter $d^{\ast }$ , which mimic actual sediment grains. The initial position of the spheres is obtained by simulating the settling of a large number $N_{p}$ of particles in the still fluid until the particles accumulate on the plane $x_{2}^{\ast }=0$ . Then, the particles in contact with the plane $x_{2}^{\ast }=0$ are kept fixed while the others are free to move. The thickness $x_{2\,bottom}^{(init)}$ of the particle layer, at the beginning of each run, is indicated in table 2. However, the reader should be aware that the fluid action is able to move only a few surficial layers of particles and many layers of particles practically do not move during the simulations.

Table 2. Size of the computational domain, initial bed elevation, number of grid points and number of sediment particles.

2.1 The fluid motion

The hydrodynamic problem is written in dimensionless form, introducing the following variables:

(2.2a,b ) $$\begin{eqnarray}t=t^{\ast }\unicode[STIX]{x1D714}^{\ast },\quad (x_{1},x_{2},x_{3})=\frac{(x_{1}^{\ast },x_{2}^{\ast },x_{3}^{\ast })}{\unicode[STIX]{x1D6FF}^{\ast }},\end{eqnarray}$$

(2.3a,b ) $$\begin{eqnarray}(u_{1},u_{2},u_{3})=\frac{(u_{1}^{\ast },u_{2}^{\ast },u_{3}^{\ast })}{U_{0}^{\ast }},\quad p=\frac{p^{\ast }}{\unicode[STIX]{x1D71A}^{\ast }U_{0}^{\ast 2}}.\end{eqnarray}$$

In (2.3), $t^{\ast }$ is time and $u_{1}^{\ast },u_{2}^{\ast },u_{3}^{\ast }$ are the fluid velocity components along the $x_{1}^{\ast }$ -, $x_{2}^{\ast }$ - and $x_{3}^{\ast }$ -directions, respectively.

Using (2.3), the continuity and momentum equations read

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{j}}=0, & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}t}+\frac{R_{\unicode[STIX]{x1D6FF}}}{2}u_{j}\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}=-\frac{R_{\unicode[STIX]{x1D6FF}}}{2}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x1D6FF}_{i1}\sin (t)+\frac{1}{2}\frac{\unicode[STIX]{x2202}^{2}u_{i}}{\unicode[STIX]{x2202}x_{k}\unicode[STIX]{x2202}x_{k}}+f_{i}, & \displaystyle\end{eqnarray}$$

where the pressure gradient is written as the sum of two terms. One term ( $-\text{sin}(t)$ ) is the imposed streamwise pressure gradient, which is uniform and drives the fluid oscillations. The other term ( $(R_{\unicode[STIX]{x1D6FF}}/2)\,\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}x_{i}$ ) is associated with the vortex structures shed by the sediment grains or with the turbulent eddies and is an output of the numerical simulations.

At the lower boundary of the fluid domain ( $x_{2}=0$ ), where a rigid wall is located, the no-slip condition is enforced,

(2.6) $$\begin{eqnarray}(u_{1},u_{2},u_{3})=(0,0,0),\end{eqnarray}$$

while at the upper boundary ( $x_{2}=L_{x2}$ ) the free-stream (free-slip) condition is enforced,

(2.7a,b ) $$\begin{eqnarray}\left({\displaystyle \frac{\unicode[STIX]{x2202}u_{1}}{\unicode[STIX]{x2202}x_{2}}},{\displaystyle \frac{\unicode[STIX]{x2202}u_{3}}{\unicode[STIX]{x2202}x_{2}}}\right)=(0,0),\quad u_{2}=0.\end{eqnarray}$$

Moreover, periodic boundary conditions are enforced in the homogeneous directions ( $x_{1},x_{3}$ ), because the computational box is chosen large enough to include the largest vortex structures of the flow.

The hydrodynamic problem is solved numerically by means of a finite difference approach in a computational domain of dimensions $L_{x1}$ , $L_{x2}$ and $L_{x3}$ in the streamwise, wall-normal and spanwise directions, respectively. A uniform grid is introduced, with $N_{x1},N_{x2},N_{x3}$ grid points along the three directions.

The numerical scheme is the same as that used by Kidanemariam & Uhlmann (2014a , 2017) and Mazzuoli et al. (2016, 2019). Standard centred second-order finite difference approximations are used to approximate the spatial derivatives, written using a uniform, staggered Cartesian grid, while the time advancement of the Navier–Stokes equations is made using a fractional-step method based upon the combination of explicit (three-step Runge–Kutta) and implicit (Crank–Nicolson) discretisations of the nonlinear and viscous terms, respectively.

The continuity and momentum equations are solved throughout the whole computational domain, including the space occupied by the solid particles, which are immersed in the fluid and move close to the bottom. The no-slip condition at the sediment–fluid interface is enforced, using the immersed-boundary technique (Uhlmann 2005), by means of the terms  $f_{i}$ , added to the right-hand side of (2.5). The numerical code has been widely tested (see e.g. Mazzuoli et al. (2016)).

2.2 The sediment motion

The sediment grains, which are modelled as spherical particles of uniform diameter  $d^{\ast }$ , are moved according to the Newton–Euler equations:

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle m_{p}^{\ast }\frac{\text{d}u_{i}^{(p)\ast }}{\text{d}t^{\ast }}=\int _{S^{\ast }}\unicode[STIX]{x1D70E}_{ij}^{(f)\ast }n_{j}\,\text{d}S^{\ast }+W_{i}^{\ast }+F_{i}^{(p)\ast }, & \displaystyle\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle I_{p}^{\ast }\frac{\text{d}\unicode[STIX]{x1D714}_{i}^{(p)\ast }}{\text{d}t^{\ast }}=\int _{S^{\ast }}\unicode[STIX]{x1D716}_{ijk}r_{j}^{\ast }\unicode[STIX]{x1D70E}_{km}^{(f)\ast }n_{m}\,\text{d}S^{\ast }+T_{i}^{(p)\ast }, & \displaystyle\end{eqnarray}$$

where $m_{p}^{\ast }$ is the mass of a single spherical particle, $I_{p}^{\ast }$ is its moment of inertia and $\unicode[STIX]{x1D716}_{ijk}$ denotes the Levi-Civita symbol. Moreover, $u_{i}^{(p)\ast }$ and $\unicode[STIX]{x1D714}_{i}^{(p)\ast }$ are the $i$ th components of the particle linear and angular velocity, respectively ( $i=1,2,3$ ). Finally, $\unicode[STIX]{x1D70E}_{ij}^{(f)\ast }$ is the fluid stress tensor, $r_{j}^{\ast }$ ( $j=1,2,3$ ) is the vector from the centre of the particle to the generic point on its surface, $n_{m}$ ( $m=1,2,3$ ) is a normal unit vector pointing outwards from the surface of the particle, $W_{i}^{\ast }$ is the weight of the particle and $F_{i}^{(p)\ast }$ and $T_{i}^{(p)\ast }$ indicate the force and torque due to inter-particle collisions. It follows that the phenomena associated with the grain size distribution and the irregular shape of the sand grains are not considered.

The motion of the sediment grains turns out to be controlled by their specific gravity $s=\unicode[STIX]{x1D71A}_{s}^{\ast }/\unicode[STIX]{x1D71A}^{\ast }$ and their dimensionless size $d=d^{\ast }/\unicode[STIX]{x1D6FF}^{\ast }$ even though it is common to use also the particle Reynolds number $R_{p}=\sqrt{(s-1)g^{\ast }d^{\ast 3}}/\unicode[STIX]{x1D708}^{\ast }$ (often known also as the Galileo number). The values of the parameters for the simulations presently considered are indicated in table 1, while the size of the computational domain and the number of grid points employed in each run are listed in table 2. In particular, one DNS was carried out, for the same values of the parameters as those of $\text{run}~2$ , by fixing the spheres at their resting positions. This run is indicated by ‘ $\text{run}~2~\text{(fix)}$ ’ in table 1.

The force and torque due to the grain–grain contacts are evaluated by means of a discrete-element model (DEM), which is based upon a linear mass–spring–damper model of particle interaction. More details on the evaluation of particle dynamics can be found in Kidanemariam & Uhlmann (2014b ).

Since the temporal scale of the grain collisions is $O(100)$ times smaller than the temporal scale of the oscillating flow, the position of colliding particles is evaluated by splitting each time step of the fluid solver into $O(100)$ substeps, during which the hydrodynamic force is assumed to be constant. The DEM model asks for the specification of the values of the following parameters: the ‘force range’, the normal stiffness, the Coulomb friction coefficient and the value of the restitution coefficient. These parameters are given values essentially equal to those of Mazzuoli et al. (2016) (see table 3).

Table 3. Values of the DEM parameters: dimensionless normal stiffness $k_{n}=(6/\unicode[STIX]{x03C0})\times k_{n}^{\ast }\unicode[STIX]{x0394}x^{\ast }/(d^{\ast 3}g^{\ast }\unicode[STIX]{x1D71A}_{s}^{\ast })$ , Coulomb friction coefficient $\unicode[STIX]{x1D707}_{cf}$ , restitution coefficient $\unicode[STIX]{x1D700}_{d}^{(part)}$ and the ‘force range’ $\unicode[STIX]{x1D6FF}_{c}^{\ast }/\unicode[STIX]{x0394}x^{\ast }$ .

2.3 Average operators

Since the bed surface preserves essentially a horizontal profile during each phase of the oscillation period for all the simulations presently considered, both the flow and the particle motion are assumed statistically homogeneous over horizontal planes. Thus plane averages are performed of quantities associated with either fluid or particle phases using the definitions provided by Kidanemariam & Uhlmann (2014b ). In particular, with reference to the sample box ${\mathcal{V}}(\boldsymbol{x})$ of size $L_{x1}\times \unicode[STIX]{x0394}x_{2}\times L_{x3}$ , centred in $\boldsymbol{x}$ , the plane average, $\langle \unicode[STIX]{x1D713}^{(f)}\rangle$ , of the generic fluid property $\unicode[STIX]{x1D713}^{(f)}$ is computed as

(2.10) $$\begin{eqnarray}\langle \unicode[STIX]{x1D713}^{(f)}\rangle (\boldsymbol{x},t)={\displaystyle \frac{\displaystyle \int _{{\mathcal{V}}}\unicode[STIX]{x1D713}^{(f)}(\boldsymbol{x},t)\displaystyle \mathop{\sum }_{\ell =1}^{N_{p}}\unicode[STIX]{x1D711}^{(f)}(\boldsymbol{x}-\boldsymbol{x}^{\ell }(t))\,\text{d}V}{\displaystyle \int _{{\mathcal{V}}}\displaystyle \mathop{\sum }_{\ell =1}^{N_{p}}\unicode[STIX]{x1D711}^{(f)}(\boldsymbol{x}-\boldsymbol{x}^{\ell }(t))\,\text{d}V}},\end{eqnarray}$$

while the plane average of a particle property, say $\unicode[STIX]{x1D713}^{(p)}$ , is defined as

(2.11) $$\begin{eqnarray}\langle \unicode[STIX]{x1D713}^{(p)}\rangle (\boldsymbol{x},t)={\displaystyle \frac{\displaystyle \int _{{\mathcal{V}}}\unicode[STIX]{x1D713}^{(p)}(\boldsymbol{x},t)\displaystyle \mathop{\sum }_{\ell =1}^{N_{p}}\unicode[STIX]{x1D711}^{(p)}(\boldsymbol{x}-\boldsymbol{x}^{\ell }(t))\,\text{d}V}{\displaystyle \int _{{\mathcal{V}}}\displaystyle \mathop{\sum }_{\ell =1}^{N_{p}}\unicode[STIX]{x1D711}^{(p)}(\boldsymbol{x}-\boldsymbol{x}^{\ell }(t))\,\text{d}V}},\end{eqnarray}$$

where $\boldsymbol{x}^{\ell }(t)$ is the position of the centre of the generic $\ell$ -particle. Moreover, the particle indicator function $\unicode[STIX]{x1D711}^{(p)}(\boldsymbol{x}-\boldsymbol{x}^{\ell }(t))$ is equal to 1 if $|\boldsymbol{x}-\boldsymbol{x}^{\ell }(t)|<d/2$ or to 0 otherwise, while $\unicode[STIX]{x1D711}^{(f)}(\boldsymbol{x}-\boldsymbol{x}^{\ell }(t))=1-\unicode[STIX]{x1D711}^{(p)}$ . In the following, where not explicitly indicated, plane-average quantities are shown omitting the angular brackets for the sake of clarity.

3 Discussion of the results

When the bottom is made up of moving sediment grains, the DNS of the Navier–Stokes and continuity equations within the bottom boundary layer generated by an oscillatory pressure gradient requires huge computational resources and ‘wall clock’ time. Hence, only a few cases are considered in the following and attention is focused on values of the parameters such that sediment particles are set into motion and the Reynolds number is large enough to trigger the appearance of turbulence. The values of the relevant dimensionless parameters of the present simulations are indicated in table 1. At this stage it is worth pointing out that neither ripples nor transverse bands of sediments, like those detected by Blondeaux, Vittori & Mazzuoli (2016), appear during the simulations. This is an effect of the size of the domains, which is too small for ripples to appear, and of the duration of the simulations, which is too short for ripples to develop. Indeed, we wanted to consider the plane bottom configuration.

For illustration, consider $\text{run}~2$ and $\text{run}~3$ , which are characterised by the dimensionless diameter $d$ equal to 0.335 and the Reynolds number $R_{\unicode[STIX]{x1D6FF}}$ equal to 750 and 1000, respectively. In order to relate these values of $R_{\unicode[STIX]{x1D6FF}}$ and $d$ to a field case, it can be easily verified that $R_{\unicode[STIX]{x1D6FF}}=750$ is the Reynolds number of the bottom boundary layer generated by a surface wave characterised by a period $T^{\ast }$ equal to 7 s and a height $H^{\ast }$ equal to approximately 1.4 m propagating in a coastal region where the water depth $h^{\ast }$ is equal to 10 m. It turns out that $\unicode[STIX]{x1D6FF}^{\ast }=1.49~\text{mm}$ and that $d=0.335$ implies $d^{\ast }=0.5~\text{mm}$ , i.e. a grain size that is coincident with the limit between medium and coarse sand.

3.1 Appearance of turbulence and turbulent kinetic energy

As discussed in Sleath (1988), on the basis of Kajiura’s (1968) criterion, the appearance of turbulence is certainly triggered during the oscillatory cycle for such values of the parameters ( $d=0.335,~R_{\unicode[STIX]{x1D6FF}}=750$ ). Indeed, the condition suggested by Kajiura (1968) for the initiation of the turbulent regime, i.e. $U_{0}^{\ast }d^{\ast }/\unicode[STIX]{x1D708}^{\ast }\geqslant 104$ , is widely satisfied. It is worth pointing out that, using the present notation, Kajiura’s criterion can be written in the form $R_{\unicode[STIX]{x1D6FF}}\geqslant 104/d$ . However, laboratory measurements (e.g. Sleath 1988) suggest that the random fluctuations of the velocity are not large and appear only during a part of the oscillatory cycle.

The flow and sediment dynamics were simulated within a computational box $24.50\unicode[STIX]{x1D6FF}^{\ast }$ long, $12.25\unicode[STIX]{x1D6FF}^{\ast }$ wide and $30.63\unicode[STIX]{x1D6FF}^{\ast }$ high (see table 2 summarising the main parameters of the numerical box and grid). The size of the box is similar to that used by Verzicco & Vittori (1996) and Vittori & Verzicco (1998) for their simulations of turbulence dynamics in an oscillatory boundary layer and turns out to be large enough for turbulence generation (minimal flow unit). Only the height of the box is significantly larger because, in the present simulation, a large number of spherical particles are deposited on the bottom. Before starting each simulation, particles were located randomly over the computational domain and allowed to settle in the resting fluid. Then, they were ‘shaken’ until a closely packed configuration was attained. Finally, the layer of particles in contact with the wall was fixed, while those whose centre was located above the elevation $x_{2\,bottom}^{\ast (init)}$ , indicated in table 2, were removed to guarantee the plane bottom configuration (cf. figure 2 a). The number of particles that remain is indicated in table 2. The grid size is equispaced and uniform along the three coordinates and such that 10 grid points are present per grain diameter.

Figure 2. Bed configuration (a) at the initial state and (b) at the early deceleration phase $\unicode[STIX]{x1D711}=0.02$ ( $\text{run}~3$ ). The broken red horizontal lines indicate the bottom surface elevation and the maximum elevation reached by particles. Light grey particles are resting while the other particles are moving and the black ones lay above $x_{2,bottom}$ .

Figure 3. Streamwise velocity component plotted versus the phase $\unicode[STIX]{x1D711}$ during the second cycle for $x_{1}=L_{x1}/2$ and $x_{3}=L_{x3}/2$ and different values of $x_{2}$ : (a)  $x_{2}^{\ast }=x_{2,bottom}^{\ast }+3\unicode[STIX]{x1D6FF}^{\ast }$ , (b)  $x_{2}^{\ast }=x_{2,bottom}^{\ast }+\unicode[STIX]{x1D6FF}^{\ast }$ , (c)  $x_{2}^{\ast }=x_{2,bottom}^{\ast }+0.5\unicode[STIX]{x1D6FF}^{\ast }$ , and (d)  $x_{2}^{\ast }=x_{2,bottom}^{\ast }$ . Continuous line, numerical results; broken line, Stokes solution. Here $R_{\unicode[STIX]{x1D6FF}}=750$ and $d=d^{\ast }/\unicode[STIX]{x1D6FF}^{\ast }=0.335$ ( $\text{run}~2$ ).

To give a qualitative idea of the flow field that is generated close to the bottom for these values of the parameters and to show how turbulence appearance and the moving grains affect the velocity field, figure 3 shows the streamwise velocity component plotted versus the phase $\unicode[STIX]{x1D711}$ during the second flow cycle, for different values of the distance $x_{2}-x_{2,bottom}$ from the bottom and for $x_{1}=L_{x1}/2$ and $x_{3}=L_{x3}/2$ . The Stokes solution is also plotted in figure 3 to allow an easy comparison of the numerical results with the laminar solution. The distance of the numerical velocity probes from the time-average bottom elevation $x_{2,bottom}$ is evaluated assuming that the instantaneous bottom surface elevation coincides with the plane where the plane-averaged volume fraction of the solid particles $\langle \unicode[STIX]{x1D719}_{s}\rangle$ reaches the value 0.1 (see Mazzuoli et al. (2019) for further details on the computation of the bed surface). It turns out that $x_{2,bottom}^{\ast }=6.60\unicode[STIX]{x1D6FF}^{\ast }$ , with fluctuations of the bottom elevation, during the flow cycle, ranging between $6.51\unicode[STIX]{x1D6FF}^{\ast }$ and $6.80\unicode[STIX]{x1D6FF}^{\ast }$ (i.e. of the order of $d^{\ast }$ ). This heuristic assumption might be modified taking into account that significant local instantaneous fluid velocities can be found even for $x_{2}$ smaller than $x_{2,bottom}$ when the Reynolds number is large enough to induce sediment motion and the spherical particles start to slide, roll and saltate on the resting particles. For example, the bottom position might be determined either by choosing a different threshold value of $\unicode[STIX]{x1D719}_{s}$ or by evaluating the value of $x_{2}$ at which either the average streamwise velocity component or the turbulent kinetic energy vanish. Figure 2(b) shows that saltating and floating particles can be present above $x_{2,bottom}$ up to the level $x_{2,free}$ , which delimits the particle-free region.

Figure 3 clearly shows that the velocity provided by the numerical simulation is characterised by large random fluctuations which appear when the velocity attains its largest values and at the beginning of the decelerating phases. However, these random velocity fluctuations are present only close to the bottom and they decrease moving far from it, till they assume negligible values when $x_{2}-x_{2,bottom}$ is larger than approximately 15 (not shown herein) where the velocity practically equals the free-stream velocity. In particular, figures 3(a) and 4, where the instantaneous velocity profile at $x_{1}=L_{x1}/2$ is plotted versus $x_{2}$ at different phases $\unicode[STIX]{x1D711}$ during the second cycle and for three different values of $x_{3}$ , show that the momentum transfer induced by turbulent fluctuations moves the overshooting of the velocity farther from the bottom and modifies its phase (cf. figure 4 b). Consequently, the displacement thickness of the boundary layer, defined as

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{dis}^{\ast }=\max _{t}\int _{x_{2\,bottom}^{\ast }}^{L_{x2}^{\ast }}\left(1-{\displaystyle \frac{U^{\ast }}{U_{e}^{\ast }}}\right)\text{d}x_{2}^{\ast },\end{eqnarray}$$

is 4.85 times larger than the displacement thickness computed for the Stokes boundary layer, which is constant and equal to $\unicode[STIX]{x1D6FF}^{\ast }$ (in (3.1) $U_{e}^{\ast }$ indicates the free-stream velocity). Such large increase of $\unicode[STIX]{x1D6FF}_{dis}^{\ast }$ is due to the presence of the particles saltating up to ${\sim}8d^{\ast }$ (i.e. $2.7\unicode[STIX]{x1D6FF}^{\ast }$ ) above the bed surface during the phases characterised by the maximum velocity (see the horizontal broken lines in figure 4 d), which, in turn, enhance turbulent fluctuations far from the bottom. In fact, for $\text{run}~2~\text{(fix)}$ ( $R_{\unicode[STIX]{x1D6FF}}=750$ , $d=0.335$ and the spheres fixed at their resting positions), $\unicode[STIX]{x1D6FF}_{dis}^{\ast }$ is equal to $2.81\unicode[STIX]{x1D6FF}^{\ast }$ (see also the values in table 1).

Figure 4. Streamwise velocity component at $x_{1}=L_{x1}/2$ and $x_{3}=2$ , 6 and 8, plotted versus the vertical coordinate $x_{2}-x_{2,bottom}$ in the near-bottom region: (a)  $\unicode[STIX]{x1D711}=1.25\unicode[STIX]{x03C0}$ , (b)  $\unicode[STIX]{x1D711}=1.50\unicode[STIX]{x03C0}$ , (c)  $\unicode[STIX]{x1D711}=1.75\unicode[STIX]{x03C0}$ and (d)  $\unicode[STIX]{x1D711}=2.00\unicode[STIX]{x03C0}$ , $\unicode[STIX]{x1D711}$ being the phase during the second cycle. Continuous line, local streamwise velocity; black broken line, plane-averaged streamwise velocity; red broken line, Stokes solution. The horizontal dash-dotted line indicates the elevation $x_{2,free}$ (cf. figure 2 b) above which the flow is free of particles. Here $R_{\unicode[STIX]{x1D6FF}}=750$ and $d=0.335$ ( $\text{run}~2$ ).

The intermittent appearance of turbulence and its vanishing value far from the bottom clearly appear in figure 5, where the dimensionless turbulent kinetic energy is plotted versus $x_{2}$ and time during the second cycle of $\text{run}~2$ and $\text{run}~3$ . The large computational costs do not allow a large number of oscillation cycles to be simulated and the turbulent kinetic energy is evaluated with respect to an average flow field that is not the phase-averaged value but the plane average of the velocity. This procedure makes the value of the turbulent kinetic energy at time $t$ slightly different from that at time $t+\unicode[STIX]{x03C0}$ and it makes the contour lines appearing in figure 5 not smooth because of the finite size of the computational domain. As already pointed out, turbulence is generated when the external velocity is maximum and turbulence generation takes place mainly close to the bottom. Then, turbulence spreads towards the irrotational region but meanwhile it decays and there are phases of the cycle such that a laminar-like flow is almost recovered. Even though the intensity of the dimensionless turbulent kinetic energy does not noticeably increase in $\text{run}~3$ with respect to that observed for $\text{run}~2$ (cf. figure 5), the larger value of the Reynolds number causes significant turbulent fluctuations to appear farther above the bed surface.

Figure 5. Turbulent kinetic energy per unit volume, normalised with $\unicode[STIX]{x1D71A}^{\ast }U_{0}^{\ast 2}$ , plotted versus the phase $\unicode[STIX]{x1D711}$ during the second cycle for $d=0.335$ : (a)  $R_{\unicode[STIX]{x1D6FF}}=750$ ( $\text{run}~2$ ) and (b)  $R_{\unicode[STIX]{x1D6FF}}=1000$ ( $\text{run}~3$ ). The red broken line indicates the instantaneous bottom surface elevation.

Figure 6. The dimensionless value of $\unicode[STIX]{x1D70F}_{b}$ plotted versus the phase $\unicode[STIX]{x1D711}$ during the second cycle for $R_{\unicode[STIX]{x1D6FF}}=750$ and $d=0.335$ ( $\text{run}~2$ ). The viscous contribution, the turbulent contribution and the contribution due to flow–particle interaction are also plotted along with the qualitative behaviour of the external velocity (broken line).

Figure 7. The dimensionless value of $\unicode[STIX]{x1D70F}_{b}$ and its contributions plotted versus the phase $\unicode[STIX]{x1D711}$ during the second cycle for $R_{\unicode[STIX]{x1D6FF}}=750$ and $d=0.335$ ( $\text{run}~2~\text{(fix)}$ ).

3.2 Evaluation of the bed shear stress

Since the velocity increases rapidly above the bed surface and the presence of particles above the bed surface is limited to a thin layer (cf. figure 4), it is reasonable to suppose that the sediment flow rate is closely related to the bed shear stress. Figures 6 and 7 show the time development of the dimensionless streamwise component $\unicode[STIX]{x1D70E}_{12}=\unicode[STIX]{x1D70E}_{12}^{\ast }/(\frac{1}{2}\unicode[STIX]{x1D70C}^{\ast }U_{0}^{\ast }\unicode[STIX]{x1D6FF}^{\ast }\unicode[STIX]{x1D714}^{\ast })$ of the averaged force per unit area exerted by the flow on the instantaneous bottom surface for $\text{run}~2$ and $\text{run}~2~\text{(fix)}$ . The reader should note that $\unicode[STIX]{x1D70E}_{12}^{\ast }(x_{2}^{\ast },t^{\ast })$ evaluated at the bed surface ( $x_{2}^{\ast }=x_{2\,bottom}^{\ast }$ ) coincides with what is commonly defined as the bottom shear stress $\unicode[STIX]{x1D70F}_{b}^{\ast }$ . Notwithstanding the fact that the force per unit area is averaged over the bottom surface, the value of $\unicode[STIX]{x1D70F}_{b}=\unicode[STIX]{x1D70F}_{b}^{\ast }/(\frac{1}{2}\unicode[STIX]{x1D70C}^{\ast }U_{0}^{\ast }\unicode[STIX]{x1D6FF}^{\ast }\unicode[STIX]{x1D714}^{\ast })$ is characterised by the presence of small random oscillations. To remove them, it would be necessary either to consider a much longer and wider computational box or to simulate a large number of cycles and to compute the phase-averaged value. The oscillations of the force per unit area are defined as ‘small’ when compared with the oscillations that are observed when the value of $\unicode[STIX]{x1D70F}_{b}$ is averaged over a much smaller horizontal surface. The value $\hat{\unicode[STIX]{x1D70F}}_{b}$ of $\unicode[STIX]{x1D70E}_{12}$ averaged over a portion of the instantaneous bottom surface that is $4\unicode[STIX]{x1D6FF}^{\ast }$ long, $2\unicode[STIX]{x1D6FF}^{\ast }$ wide and centred around the point $(x_{1}^{\ast },x_{3}^{\ast })=(L_{x1}/2,L_{x3}/2)\unicode[STIX]{x1D6FF}^{\ast }$ was computed to verify this point (the time development of $\hat{\unicode[STIX]{x1D70F}}_{b}$ is shown in the figure in the supplementary material, available at https://doi.org/10.1017/jfm.2019.1012).

Three contributions to the value of $\unicode[STIX]{x1D70F}_{b}$ appearing in figures 6 and 7 can be identified (Uhlmann 2008; Mazzuoli et al. 2018): (i) the contribution due to the viscous stress ( $\unicode[STIX]{x1D70F}_{visc}$ ), (ii) the contribution due to the turbulent stress ( $\unicode[STIX]{x1D70F}_{turb}$ ), and (iii) the contribution due to the flow–particle interactions ( $\unicode[STIX]{x1D70F}_{part}$ ). The procedure used to compute the different contributions is described in more detail by Mazzuoli et al. (2018, 2019). The time development of $\unicode[STIX]{x1D70F}_{b}$ for the $\text{run}~2~\text{(fix)}$ qualitatively agrees with that measured by Jensen et al. (1989) (see figure 1) even though the Reynolds number of the laboratory experiment is somewhat different from that of the numerical simulation but, more importantly, the bottom of the experimental apparatus was smooth. Of course, the sediment motion and, in particular, the saltating grains greatly affect turbulence dynamics. Indeed, in $\text{run}~2$ , the largest contributions to the bottom shear stress are those due to the turbulent stresses and the flow–particle interaction.

Figure 8. The dimensionless value of $\unicode[STIX]{x1D70F}_{b}$ plotted versus the phase $\unicode[STIX]{x1D711}$ during the second cycle of $\text{run}~2$ (mobile particle, solid line) and $\text{run}~2~\text{(fix)}$ (fixed particles, red dash-dotted line), i.e. for the values of the parameters $R_{\unicode[STIX]{x1D6FF}}=750$ and $d=0.335$ . Moreover, the black broken line indicates the values obtained in a run at $R_{\unicode[STIX]{x1D6FF}}=775$ over a smooth wall (results from Mazzuoli, Vittori & Blondeaux 2011b ).

To gain an idea of the role of the resting/moving sediment grains on turbulence dynamics, figure 8 shows the value of $\unicode[STIX]{x1D70F}_{b}$ computed for $\text{run}~2$ and $\text{run}~2~\text{(fix)}$ , and that computed by Mazzuoli, Vittori & Blondeaux (2011a ) and Mazzuoli et al. (2011b ), who made a DNS of the OBL over a smooth wall for $R_{\unicode[STIX]{x1D6FF}}=775$ . Even though the results of Mazzuoli et al. (2011a , b ) were obtained for a value of the Reynolds number slightly larger than that of the present simulations, figure 8 shows that the bottom shear stress of the smooth-bottom case is significantly smaller than that of $\text{run}~2~\text{(fix)}$ , where the particles were fixed and arranged in a plane-bed configuration. However, the phase when the inception of turbulence occurs and, in general, the time development of $\unicode[STIX]{x1D70F}_{b}$ are fairly comparable in these two cases. In particular, both in the simulation of Mazzuoli et al. (2011a , b ) and in $\text{run}~2~\text{(fix)}$ , turbulence appears when the free-stream velocity is decelerating.

On the other hand, the presence of mobile sediments enhances the effect of the sediment on the transition process. This fact is clearly shown by figure 9, where the total dimensionless shear stress  $\unicode[STIX]{x1D70E}_{12}$ , the Reynolds shear stress and the solid volume fraction for $\text{run}~2$ and $\text{run}~2~\text{(fix)}$ are plotted as functions of the wall-normal coordinate at flow reversal (figure 9 a) and at the phases when $\unicode[STIX]{x1D70F}_{b}$ is maximum in $\text{run}~2$ (figure 9 b) and in $\text{run}~2~\text{(fix)}$ (figure 9 c). The remarkable contribution of the Reynolds shear stress to $\unicode[STIX]{x1D70F}_{b}$ , which is associated with the presence of moving particles, can be observed in figure 9(b). Indeed, in semidilute and dense sheared suspensions (i.e. for $\unicode[STIX]{x1D719}_{s}>0.05$ ), the number of particles per unit volume that are exposed to the core flow is large enough to experience frequent contacts and generate significant overall drag. In order to understand the increase of the apparent roughness due to the particle motion, it is worth noting that, in a suspension of mono-sized spheres at $\unicode[STIX]{x1D719}_{s}=0.1$ , the distance between particles is approximately equal to their diameter and the flow resistance in steady conditions is a few times larger than that attained with the tightest arrangement of the spheres, i.e. in the ‘cannonball’ configuration. For instance, Schlichting (1936), who carried out experiments on channel flow over a plane layer of spheres of diameter $d^{\ast }$ arranged at the vertices of hexagons of side  $\ell ^{\ast }$ , found that the equivalent roughness obtained for $\ell ^{\ast }$ equal to $2d^{\ast }$ was approximately five times larger than that for closely packed spheres.

Figure 9. The dimensionless value of $\unicode[STIX]{x1D70F}_{b}$ (solid thick lines), of the Reynolds shear stress (solid thin lines) and the value of the solid volume fraction, $\unicode[STIX]{x1D719}_{s}$ (broken lines), plotted versus the wall-normal coordinate at the phases (a)  $\unicode[STIX]{x1D711}=0.5\unicode[STIX]{x03C0}$ , (b)  $\unicode[STIX]{x1D711}=0.9\unicode[STIX]{x03C0}$ and (c)  $\unicode[STIX]{x1D711}=1.035\unicode[STIX]{x03C0}$ of $\text{run}~2$ (mobile particles, black lines) and $\text{run}~2~\text{(fix)}$ (fixed particles, red lines), i.e. for the values of the parameters $R_{\unicode[STIX]{x1D6FF}}=750$ and $d=0.335$ . Dash-dotted horizontal lines indicate the bed surface elevation, $x_{2\,bottom}$ , where $\unicode[STIX]{x1D719}_{s}=0.1$ .

Because of turbulence growth, which takes place when the external flow is close to its largest values, and its subsequent damping, the values of $\unicode[STIX]{x1D70F}_{b}$ attained during flow acceleration differ from those observed during flow deceleration, even though the free-stream velocity is equal. Hence, if $\unicode[STIX]{x1D70F}_{b}$ is plotted versus the free-stream velocity $U_{e}$ , a hysteresis orbit can be observed (see figure 10 a). A similar orbit is observed even considering the curve $\unicode[STIX]{x1D70F}_{b}$ versus the dimensionless quantity $(U_{e}+\text{d}U_{e}/\text{d}t)$ , which in the laminar case should be a straight line passing through the origin of the axes, and the dimensionless quantity $U_{e}|U_{e}|$ . These results show that, even though empirical relationships do exist which allow a reliable estimate of the friction factor $f_{w}=\unicode[STIX]{x1D70F}_{b,max}^{\ast }/(\frac{1}{2}\unicode[STIX]{x1D71A}^{\ast }U_{0}^{\ast 2})$ and it is relatively easy to predict the maximum value $\unicode[STIX]{x1D70F}_{b,max}^{\ast }$ of the bottom shear stress, the time development of bottom shear stress during the flow cycle is much more difficult to predict. In order to fulfil the objectives (I) and (II) stated in § 1, let us see in the following section how the bed shear stress is related to the sediment flow rate.

Figure 10. The dimensionless value of $\unicode[STIX]{x1D70F}_{b}$ plotted versus (a) the free-stream velocity $U_{e}$ , (b) the quantity $\text{d}U_{e}/\text{d}t+U_{e}$ and (c) the quantity $U_{e}|U_{e}|$ for $R_{\unicode[STIX]{x1D6FF}}=750$ and $d=0.335$ ( $\text{run}~2$ ). In the laminar regime, the curve $\unicode[STIX]{x1D70F}_{b}$ versus $\text{d}U_{e}/\text{d}t+U_{e}$ is the straight line crossing the origin of the axes and plotted in panel (b).

3.3 Evolution of the sediment flow rate and dependence on the flow properties

Figure 11. Values of $q_{s}^{\ast }/\sqrt{(s-1)g^{\ast }d^{\ast 3}}$ (thick red line) and the signed values of $[|\unicode[STIX]{x1D70F}_{b}^{\ast }|/(\frac{1}{2}\unicode[STIX]{x1D71A}^{\ast }U_{0}^{\ast }\unicode[STIX]{x1D714}^{\ast }\unicode[STIX]{x1D6FF}^{\ast })]^{3/2}$ (thin black line) plotted versus the phase $\unicode[STIX]{x1D711}$ during the second cycle for $d=0.335$ : (a)  $R_{\unicode[STIX]{x1D6FF}}=750$ ( $\text{run}~2$ ) and (b)  $R_{\unicode[STIX]{x1D6FF}}=1000$ ( $\text{run}~3$ ).

Figure 11 shows the dimensionless sediment transport rate per unit width $q_{s}=q_{s}^{\ast }/\sqrt{(s-1)g^{\ast }d^{\ast 3}}$ (i.e. the volume of sediment grains that cross a plane $x_{1}=\text{const.}$ per unit time and unit width), averaged over the ‘homogeneous’ directions $x_{1}$ and $x_{3}$ , as a function of the phase $\unicode[STIX]{x1D711}$ during the second cycle. Since it is common practice to correlate the sediment flow rate to the power 3/2 of the bottom shear stress, in figure 11 the signed value of $|\unicode[STIX]{x1D70F}_{b}|^{3/2}$ is also plotted. This clearly shows that a fair correlation exists between $q_{s}$ and $|\unicode[STIX]{x1D70F}_{b}|^{3/2}$ for both $\text{run}~2$ (figure 11 a) and $\text{run}~3$ (figure 11 b), even though analogous differences can be observed in the two runs. Taking into account that (i) turbulence intensity during the decelerating phases of the cycle is different from that observed during the accelerating phases and (ii) the sediment particles are set into motion and transported more easily when turbulence intensity is high, it is reasonable to expect values of $q_{s}$ during the late accelerating phases to be larger than those observed during the decelerating phases, even for the same value of the bottom shear stress.

In figure 12, the value of $q_{s}$ is plotted versus the Shields parameter $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D70F}_{b}^{\ast }/((\unicode[STIX]{x1D71A}_{s}^{\ast }-\unicode[STIX]{x1D71A}^{\ast })g^{\ast }d^{\ast })$ , describing two nearly coincident orbits during one oscillation period. In fact, for values of the Shields parameter slightly larger than its critical value ( $10^{-2}\lesssim \unicode[STIX]{x1D703}\lesssim 10^{-1}$ ), the sediment transport rate depends not only on $\unicode[STIX]{x1D703}$ but also on the value of $\text{d}\unicode[STIX]{x1D703}/\text{d}t$ . A similar finding was obtained by Vittori (2003), who showed that, in an oscillatory flow, the amount of sediment grains picked up from the bed and carried into suspension correlates better with the production of turbulent kinetic energy than with the Shields parameter. The results plotted in figure 12 show that the sediment transport rate tends to vanish for a finite value of $\unicode[STIX]{x1D703}$ while it assumes a small but finite value for $\unicode[STIX]{x1D703}$ tending to zero. To understand this behaviour of $q_{s}$ , the reader should take into account that even for vanishing values of $\unicode[STIX]{x1D703}$ the sediments keep moving mainly because moving particles take some time to find a stable position on the bed surface (Clark et al. 2017) and also because of the effects of the imposed pressure gradient (Mazzuoli et al. 2019). Hence, the sediment transport tends to vanish for small values of $\unicode[STIX]{x1D703}$ when the effects of the pressure gradient and the viscous forces on the particles balance inertial effects. On the other hand, the results plotted in figure 12 show that fair predictions of the sediment transport rate can be obtained by assuming that $q_{s}$ depends only on $\unicode[STIX]{x1D703}$ when the Shields parameter assumes relatively large values.

Figure 12. Dimensionless sediment flow rate $q_{s}$ as a function of the Shields parameter $\unicode[STIX]{x1D703}$ for $R_{\unicode[STIX]{x1D6FF}}=750$ and $d=0.335$ ( $\text{run}~2$ ). Arrows indicate the orbital direction.

Some further insight into the effects of turbulence dynamics on sediment motion can be obtained looking at the time development of the vortex structures and at the associated dynamics of sediment grains. Figure 13 shows the streamwise velocity fluctuations in the horizontal plane $x_{2}-x_{2\,bottom}=0.9$ , along with sediment particles (white dots) that are picked up from the bottom and, during their motion, cross this plane. Figure 13(a) shows that low- and high-speed streaks characterise the flow field during the flow acceleration. Later, as also found by Costamagna et al. (2003), Mazzuoli et al. (2011b ) and Mazzuoli & Vittori (2016, 2019), who simulated an oscillatory boundary layer over a smooth bottom, the streaks oscillate, twist and interact (see figure 13 b,c) and then they break, generating small vortex structures and a fully turbulent flow (see figure 13 d–g). Eventually, the turbulent eddies decay because of viscous effects and the flow recovers a laminar-like behaviour (see figure 13 h). When the external velocity is maximum, the turbulence level is also high and a large amount of sediment particles are picked up from the bottom and transported by the flow. Then, at the late stages of flow deceleration, the turbulent eddies damp out, the bottom shear stress tends to vanish and the sediment particles settle down. The interested reader can look at movies 1 and 2 that are available in the supplementary material.

Figure 13. Streamwise velocity fluctuations at plane $x_{2}=7.5$ ( $x_{2}-x_{2,bottom}=0.9$ ) at different phases of the cycle: (a)  $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=1.73\unicode[STIX]{x03C0}$ ( $u_{\unicode[STIX]{x1D70F}}^{\ast }d^{\ast }/\unicode[STIX]{x1D708}^{\ast }=15$ , $\unicode[STIX]{x1D703}=0.12$ ); (b)  $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=1.75\unicode[STIX]{x03C0}$ ( $u_{\unicode[STIX]{x1D70F}}^{\ast }d^{\ast }/\unicode[STIX]{x1D708}^{\ast }=16$ , $\unicode[STIX]{x1D703}=0.13$ ); (c)  $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=1.78\unicode[STIX]{x03C0}$ ( $u_{\unicode[STIX]{x1D70F}}^{\ast }d^{\ast }/\unicode[STIX]{x1D708}^{\ast }=18$ , $\unicode[STIX]{x1D703}=0.16$ ); (d)  $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=1.80\unicode[STIX]{x03C0}$ ( $u_{\unicode[STIX]{x1D70F}}^{\ast }d^{\ast }/\unicode[STIX]{x1D708}^{\ast }=19$ , $\unicode[STIX]{x1D703}=0.18$ ); (e)  $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=1.85\unicode[STIX]{x03C0}$ ( $u_{\unicode[STIX]{x1D70F}}^{\ast }d^{\ast }/\unicode[STIX]{x1D708}^{\ast }=22$ , $\unicode[STIX]{x1D703}=0.24$ ); (f)  $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=1.91\unicode[STIX]{x03C0}$ ( $u_{\unicode[STIX]{x1D70F}}^{\ast }d^{\ast }/\unicode[STIX]{x1D708}^{\ast }=20$ , $\unicode[STIX]{x1D703}=0.20$ ); (g)  $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=2.05\unicode[STIX]{x03C0}$ ( $u_{\unicode[STIX]{x1D70F}}^{\ast }d^{\ast }/\unicode[STIX]{x1D708}^{\ast }=19$ , $\unicode[STIX]{x1D703}=0.17$ ); and (h)  $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=2.40\unicode[STIX]{x03C0}$ ( $u_{\unicode[STIX]{x1D70F}}^{\ast }d^{\ast }/\unicode[STIX]{x1D708}^{\ast }=6$ , $\unicode[STIX]{x1D703}=0.02$ ). Here $R_{\unicode[STIX]{x1D6FF}}=750$ and $d=0.335$ ( $\text{run}~2$ ). The full sequences of visualisations of streamwise velocity and spanwise vorticity fluctuations at plane $x_{2}=7.5$ can be found in movie 1 and movie 2 of the supplementary material, respectively.

We have discussed the flow field and sediment dynamics during the second cycle since, for $R_{\unicode[STIX]{x1D6FF}}=750$ , the average results obtained during this cycle are similar to those obtained during the first cycle, thus suggesting that the flow has attained its periodic status. The results of the numerical simulation carried out for $R_{\unicode[STIX]{x1D6FF}}=450$ and $d=0.335$ ( $\text{run}~1$ ) allow one to appreciate more easily the dynamics of both the vortex structures and the sediment particles because the process that leads to the appearance of turbulence is slower and similar to that observed experimentally over a smooth wall by Carstensen et al. (2010), who made flow visualisations for values of the Reynolds number close to its critical value. Indeed, Kajiura’s (1968) criterion suggests that, for $d=0.335$ , the critical value of the Reynolds number falls around 310. Hence, for $R_{\unicode[STIX]{x1D6FF}}=450$ , turbulence is expected to be weak and large fluctuations of turbulence intensity from cycle to cycle are expected to be present.

Figure 14. Dimensionless fluctuating spanwise component of vorticity $\unicode[STIX]{x1D6FA}_{3}^{\prime \ast }\unicode[STIX]{x1D6FF}^{\ast }/U_{0}^{\ast }=\unicode[STIX]{x1D6FA}_{3}^{\prime }=\unicode[STIX]{x1D6FA}_{3}-\langle \unicode[STIX]{x1D6FA}_{3}\rangle _{L_{x3},L_{x1}}$ at $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=5.7\unicode[STIX]{x03C0}$ and (a)  $x_{3}=2$ , (b)  $x_{3}=6$ and (c)  $x_{3}=8$ for $R_{\unicode[STIX]{x1D6FF}}=450$ and $d=0.335$ ( $\text{run}~1$ ).

Figure 14 shows the spanwise vorticity component at three streamwise–vertical planes defined by $x_{3}=2$ , $x_{3}=6$ and $x_{3}=8$ and at $t=5.7\unicode[STIX]{x03C0}$ , once the plane-averaged value is removed. The plots show that, close to the bottom, positive and negative regions alternate in the streamwise direction with a wavelength of approximately $12.5\unicode[STIX]{x1D6FF}^{\ast }$ . These coherent spanwise vortex structures appear when the free-stream velocity is almost maximum and then they are convected in the streamwise direction, generating almost regular oscillations of the velocity field as it appears in figure 15(a,b), where the streamwise velocity component is plotted versus $\unicode[STIX]{x1D714}^{\ast }t^{\ast }$ at two of the locations considered in figure 3. Later on, these vortex structures attain their maximum intensity and eventually decay because of viscous effects and disappear.

Figure 15. Streamwise velocity component plotted versus time for $x_{1}=L_{x1}/2$ , $x_{3}=L_{x3}/2$ and (a,c)  $x_{2}=x_{2,bottom}+1$ and (b,d)  $x_{2}=x_{2,bottom}+0.5$ . Continuous line, numerical results; broken line, Stokes solution. Here $R_{\unicode[STIX]{x1D6FF}}=450$ and $d=0.335$ ( $\text{run}~1$ ).

As already pointed out, similar spanwise vortices were visualised by Carstensen et al. (2010) during laboratory experiments and by Mazzuoli et al. (2011a , b ), who reproduced the experiments by Carstensen et al. (2010) by means of DNS and confirmed that the first vortex structures that appear during the transition process are two-dimensional spanwise vortices. In fact, during the first oscillation cycles, notwithstanding the presence of the sediment particles which make the bottom rough and certainly affect the transition process (see e.g. Blondeaux & Vittori 1991), the present numerical findings show vortex structures that are qualitatively in agreement with both the predictions of the linear stability analysis of Blondeaux & Seminara (1979) and the results of the DNS of Vittori & Verzicco (1998), Costamagna et al. (2003) and Bettencourt & Dias (2018), who considered a smooth bottom. Blondeaux & Seminara (1979), by using a momentary criterion of instability, showed that the laminar Stokes boundary layer is linearly and momentarily unstable when the Reynolds number $R_{\unicode[STIX]{x1D6FF}}$ is larger than 86 and the fastest-growing mode is two-dimensional and characterised by a streamwise wavelength of approximately $12.5\unicode[STIX]{x1D6FF}^{\ast }$ . For values of $R_{\unicode[STIX]{x1D6FF}}$ close to its critical value, the instability predicted by the linear analysis is restricted to a small part of the cycle, and during the remaining parts of the cycle, the amplification rate becomes negative and the flow recovers a laminar-like behaviour. Larger values of $R_{\unicode[STIX]{x1D6FF}}$ widen the unstable parts of the cycle and lead to larger perturbations. However, the flow is stable ‘on average’. Indeed the results of Blennerhassett & Bassom (2002) show that small perturbations of the Stokes solution are characterised by an average growth only when the Reynolds number is larger than 1416.

Transition to turbulence is triggered by nonlinear effects, which can no longer be neglected when the perturbations attain large values during the momentarily unstable phases, and it is also affected by wall imperfections (Blondeaux & Vittori 1994) and three-dimensional effects. Moreover, in the present simulations, the resting/moving sediment grains certainly affect the transition process. Figure 15(c,d), where the streamwise velocity component is plotted versus $\unicode[STIX]{x1D714}^{\ast }t^{\ast }$ for $R_{\unicode[STIX]{x1D6FF}}=450$ and $d=0.335$ , shows that turbulent velocity fluctuations appear during the fourth cycle, being triggered by the vortex structures generated during the previous cycles, which decayed but were still strong enough to trigger the growth of large perturbations of the laminar flow. This peculiar behaviour of turbulence dynamics clearly appears in figure 16, where the value of  $K^{\ast }$ , i.e. the plane-averaged turbulent kinetic energy per unit area of the bottom integrated over the whole computational domain, is plotted versus time. Indeed, large variations of $K^{\ast }$ from half-cycle to half-cycle can be observed in figure 16. During the phases characterised by the growth of $\int _{L_{x2}^{\ast }}K^{\ast }\,\text{d}x_{2}^{\ast }/\unicode[STIX]{x1D71A}^{\ast }U_{0}^{\ast 2}\unicode[STIX]{x1D6FF}^{\ast }$ , the sediment flow rate tends to increase but then it decreases even if the integral of $K^{\ast }$ still keeps growing. This result can be understood by observing that large values of $K^{\ast }$ far from the bottom are encountered later than close to the bottom, where turbulence is produced, as an effect of the diffusion of turbulent fluctuations (see e.g $\text{run}~2$ in figure 5). Therefore, the maxima of the integral of $K^{\ast }$ lag behind the maxima of $K^{\ast }$ close to the bed, which are associated with large values of the bed shear stress and of the sediment flow rate. The fact that the sporadic inception to turbulence in $\text{run}~1$ is due to the presence of particles appears reasonable, since Carstensen et al. (2010) could not observe the transition to turbulence for $R_{\unicode[STIX]{x1D6FF}}=450$ in the absence of particles. However, it is evident from the present results that the values of the parameters characterising $\text{run}~1$ lay on the edge of the intermittently turbulent region of the parameter space.

Figure 16. Time development of the dimensionless sediment flow rate $q_{s}$ (red line) and the turbulent kinetic energy $K$ (black line), integrated along the $x_{2}$ -direction from the bed surface to the top of the computational domain. Here $R_{\unicode[STIX]{x1D6FF}}=450$ and $d=0.335$ ( $\text{run}~1$ ).

Figure 17. Time development of the dimensionless sediment flow rate $q_{s}$ (red line) and the Shields parameter $\unicode[STIX]{x1D703}$ (black line). Here $R_{\unicode[STIX]{x1D6FF}}=450$ and $d=0.335$ ( $\text{run}~1$ ).

The value of the plane-averaged sediment flow rate $q_{s}$ per unit width is plotted in figure 17 along with the value of the Shields parameter as functions of time. First of all, it is worth pointing out that large fluctuations of the Shields parameter $\unicode[STIX]{x1D703}$ from cycle to cycle are present, because of the large fluctuations of turbulence intensity: the appearance of strong turbulent eddies close to the bottom gives rise to large velocity gradients at the bottom and large values of  $\unicode[STIX]{x1D70F}_{b}$ . Moreover, as already pointed out, the values of $q_{s}$ during flow acceleration differ from the values during flow deceleration, even though the Shields parameter assumes the same value. Hence, the results show that, for $R_{\unicode[STIX]{x1D6FF}}$ close to its critical value, a sediment transport rate predictor based on the assumption that $q_{s}$ depends only on $\unicode[STIX]{x1D703}$ cannot provide good predictions. If turbulence is strong, a large amount of sediment is picked up from the bed and easily transported by the flow in the saltating mode. If turbulence is weak, the moving grains roll and slide along the bottom, interacting with the resting particles, and the sediment transport rate is much smaller.

Figure 18. Top view of the computational domain. Green surfaces visualise the isocontour of $\unicode[STIX]{x1D706}_{2}^{\ast }\unicode[STIX]{x1D6FF}^{\ast 2}/U_{0}^{\ast 2}=-0.5$ at (a)  $t=3.70\unicode[STIX]{x03C0}$ , (b)  $t=3.75\unicode[STIX]{x03C0}$ , (c)  $t=3.80\unicode[STIX]{x03C0}$ and (d)  $t=4.00\unicode[STIX]{x03C0}$ . Spheres are coloured from blue to red on the basis of the magnitude of their instantaneous velocity ranging between 0 and $0.1U_{0}^{\ast }$ . White spheres are essentially at rest. Here $R_{\unicode[STIX]{x1D6FF}}=750$ and $d=0.335$ ( $\text{run}~2$ ). (e) The phase of the wave cycle to which panels (a–d) refer in terms of the velocity far from the bottom. A full sequence of similar top-view visualisations can be found in movie 3 in the supplementary material, which shows isosurfaces of spanwise-vorticity fluctuations.

To relate the sediment motion to the dynamics of the coherent vortex structures and the turbulent eddies generated by the transition from the laminar to the turbulent regime, let us consider again the numerical simulation carried out for $R_{\unicode[STIX]{x1D6FF}}=750$ and $d=0.335$ . Figure 18 shows a top view of the bed, where the sediment grains are coloured according to their velocity. Simultaneously, figure 18 shows the coherent vortex structures that characterise the turbulent flow and are visualised by the $\unicode[STIX]{x1D706}_{2}$ -criterion (Jeong & Hussain 1995). The surfaces that appear in figure 18 are characterised by a small negative value of the second eigenvalue $\unicode[STIX]{x1D706}_{2}$ of the matrix $\unicode[STIX]{x1D63F}^{2}+\unicode[STIX]{x1D734}^{2}$ , $\unicode[STIX]{x1D63F}$ and $\unicode[STIX]{x1D734}$ being the symmetric and antisymmetric parts of the gradient of the velocity field. The coherent vortex structures generated by turbulence appearance cause local high values of the bottom shear stress when they interact with the bottom and the sediment grains move with large velocities in the areas below the coherent vortices. In particular, it appears that patches of sediments randomly distributed over the bed are convected in the flow direction when coherent vortex structures are generated by turbulence appearance and move close to the bottom (see figure 18 a). Then, at $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=3.75\unicode[STIX]{x03C0}$ , the turbulent eddies give rise to a band aligned with the $x_{1}$ -direction and almost in the centre of the computational domain. Below the band of vortex structures, the sediment grains move with the largest velocities. Later, turbulence spreads over the whole domain and an intense sediment transport is observed over the entire bed. Eventually, the turbulence decays, thus the vortex structures that characterised the flow field weaken and the sediment transport decreases till it becomes negligible at approximately $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=4.3\unicode[STIX]{x03C0}$ .

Figure 19. Dimensionless fluctuating spanwise component of vorticity $\unicode[STIX]{x1D6FA}_{3}^{\prime \ast }\unicode[STIX]{x1D6FF}^{\ast }/U_{0}^{\ast }=\unicode[STIX]{x1D6FA}_{3}^{\prime }=\unicode[STIX]{x1D6FA}_{3}-\langle \unicode[STIX]{x1D6FA}_{3}\rangle _{L_{x3},L_{x1}}$ , in the plane $x_{3}=6\unicode[STIX]{x1D6FF}^{\ast }$ for $R_{\unicode[STIX]{x1D6FF}}=750$ and $d=0.335$ ( $\text{run}~2$ ) at (a)  $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=3.6\unicode[STIX]{x03C0}$ , (b)  $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=3.7\unicode[STIX]{x03C0}$ and (c)  $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=3.8\unicode[STIX]{x03C0}$ .

These findings are further supported by figure 19, which shows that the sediment transport rate is largely affected not only by the external flow, but also by the interaction of the turbulent eddies with the sediment grains. At $\unicode[STIX]{x1D714}^{\ast }t^{\ast }=3.6\unicode[STIX]{x03C0}$ , turbulence is weak and only the grains in unstable positions slowly roll and move to attain more stable positions being dragged by the external flow. Then, turbulent fluctuations become intense, in particular close to the bottom, and more sediments start to move because the fluctuating velocity components cause peaks of the hydrodynamic force acting on the sediment particles (see figure 19 b). Later, the external velocity becomes larger and the turbulent eddies become more intense as well, and the sediment grains not only move but also start to saltate, being picked up from the bed.

Figure 20. Turbulent kinetic energy per unit bottom area (black continuous line), integrated along the $x_{2}$ -direction (a) from the bed surface to the top of the computational domain and (b) within a layer $1\unicode[STIX]{x1D6FF}$ thick above the bed surface, plotted as a function of the phase $\unicode[STIX]{x1D711}$ during the second cycle for $R_{\unicode[STIX]{x1D6FF}}=750$ and $d=0.335$ ( $\text{run}~2$ ). The red broken lines indicate the dimensionless sediment flow rate.

The time development of $K^{\ast }$ is shown in figure 20(a), where it appears that turbulence grows around $\unicode[STIX]{x1D711}=0.67\unicode[STIX]{x03C0}$ and $1.67\unicode[STIX]{x03C0}$ , attains its maximum value and then slowly decreases later on, assuming relatively small values after flow inversion, thus being loosely related to sediment transport, which is also plotted in the same figure. Even though the growth of the turbulent kinetic energy and the growth of the sediment transport rate take place almost simultaneously, $K=K^{\ast }/(\unicode[STIX]{x1D70C}^{\ast }U_{0}^{\ast 2}\unicode[STIX]{x1D6FF}^{\ast })$ attains its maximum value later than $q_{s}$ and it keeps large values even when $q_{s}$ vanishes. This finding can be easily understood by taking into account that the sediment pick-up rate is mainly related to the upward velocity component generated by the turbulent eddies which are present close to the bottom. Figure 5 shows that strong turbulent vortex structures are generated close to the bed when transition to turbulence takes place. However, later on, turbulence diffuses far from the bottom and it no longer interacts with the moving sediments, which slow down and come to rest. Hence, the curve $q_{s}(t)$ follows more closely that obtained by considering the turbulent kinetic energy $K$ per unit area integrated from $x_{2}^{\ast }=0$ up to the horizontal plane located at distance $\unicode[STIX]{x1D6FF}^{\ast }$ from the instantaneous bottom surface (see figure 20 b).

Figure 21. Turbulent kinetic energy per unit bottom area (black continuous line) integrated along the $x_{2}$ -direction within a layer $1\unicode[STIX]{x1D6FF}$ thick above the bed surface, plotted as a function of time for $R_{\unicode[STIX]{x1D6FF}}=1000$ and $d=0.335$ ( $\text{run}~3$ ). The broken line indicates the dimensionless sediment flow rate.

Figure 22. (a,c) The p.d.f., in linear and semi-logarithmic scales, of the length of the jumps experienced by particles that cross the plane $x_{2}^{\ast }=\hat{x}_{2}^{\ast }=8.80\unicode[STIX]{x1D6FF}^{\ast }$ . (b,d) Similarly, the p.d.f. of the maximum height above the plane $x_{2}^{\ast }=\hat{x}_{2}^{\ast }$ reached by particles. Broken red lines, which approximate the p.d.f.s, correspond to the functions (c)  $\text{p.d.f.}\cdot \unicode[STIX]{x1D6FF}=0.1\exp (-0.07x_{range}^{\ast (tr)}/\unicode[STIX]{x1D6FF}^{\ast })$ and (d)  $\text{p.d.f.}\cdot \unicode[STIX]{x1D6FF}=2.72\exp (-1.8\max _{t}(x_{2}^{\ast (tr)}/\unicode[STIX]{x1D6FF}^{\ast })+15.84)$ . Here $R_{\unicode[STIX]{x1D6FF}}=1000$ and $d=0.335$ ( $\text{run}~3$ ).

If the Reynolds number $R_{\unicode[STIX]{x1D6FF}}$ is increased, turbulence strength increases as shown by figure 21, where $K$ is plotted versus the phase within the cycle for $R_{\unicode[STIX]{x1D6FF}}=1000$ ( $\text{run}~3$ ). Because of the high turbulence intensity, a great number of particles is picked up from the bed during the phases characterised by large values of the bottom shear stress. Once picked up, some of the particles cover large distances without interacting with the bottom, before coming to rest again. This result is clearly shown by figure 22(a,c), where the probability density function (p.d.f.) of the length of particle jumps is plotted for the particles that are moving above the plane $x_{2}^{\ast }-x_{2,bottom}^{\ast }=0.67\unicode[STIX]{x1D6FF}^{\ast }=2.01d^{\ast }$ , suitably chosen to distinguish the saltating particles from the particles that roll and slide on the resting particles. For $\text{run}~3$ , during the oscillation period, the bottom elevation fluctuates by an amount approximately equal to $1d^{\ast }$ above and $1d^{\ast }$ below the time-average bottom elevation $x_{2,bottom}^{\ast }=8.13\unicode[STIX]{x1D6FF}^{\ast }$ . Despite the fact that the average jump length is rather large, being approximately $36.2d^{\ast }$ , the median value is equal to $19.6d^{\ast }$ , since most of the suspended particles rapidly redeposit, as shown by figure 22(b,d), where the p.d.f. of the height of the particle jumps is plotted for the same particles considered in figure 22(a,c). Taking into account that the ratio $d^{\ast }/\unicode[STIX]{x1D6FF}^{\ast }$ is equal to 0.335, the average height of the particle jumps turns out to be much smaller than the thickness of the region above the bottom where turbulence is intense (see e.g. figure 5) and it can be assumed that no significant suspended load is present. Indeed, the displacement boundary layer thickness for $\text{run}~2$ , defined by (3.1), is $\unicode[STIX]{x1D6FF}_{dis}^{\ast }=7.15\unicode[STIX]{x1D6FF}^{\ast }$ . This numerical finding is consistent with the empirical criteria usually employed to determine the presence of sediment in suspension (e.g. Bagnold 1966; Sumer & Fredsøe 2002). The ratio $u_{\unicode[STIX]{x1D70F}}^{\ast }/v_{s}^{\ast }$ between the shear velocity and the fall velocity of the particles ( $v_{s}^{\ast }=\sqrt{(s-1)g^{\ast }d^{\ast }}$ ) is smaller than one and the Reynolds number, $u_{\unicode[STIX]{x1D70F},max}^{\ast }d^{\ast }/\unicode[STIX]{x1D708}^{\ast }$ , is larger than 20 when the maximum value of the bottom shear stress is considered (see table 1). Indeed, the maximum value of $u_{\unicode[STIX]{x1D70F}}^{\ast }/v_{s}^{\ast }$ for $\text{run}~3$ turns out to be approximately 0.6.

Notwithstanding the fact that the value of $K$ never vanishes and the random velocity fluctuations remain large during the whole flow cycle and tend to pick up the sediments from the bottom, there are phases such that $q_{s}$ vanishes as shown by figure 21. The vanishing of the sediment transport rate is related to the vanishing of the average bottom shear stress, which takes place twice during the cycle.

Figure 23. Dimensionless sediment flow rate $q_{s}$ as a function of the absolute value of the Shields parameter for runs at $R_{\unicode[STIX]{x1D6FF}}=450$ (turbulent cycles), 750 and 1000 ( $\text{run}~1$ , black line; $\text{run}~2$ , red line; and $\text{run}~3$ , green line). The blue dash-dotted line represents the formula by Wong & Parker (2006). Arrows indicate the orbital direction.

The results of the numerical simulations are summarised in figure 23, which shows the value of $q_{s}$ as a function of the Shields parameter $\unicode[STIX]{x1D703}$ for the three values of the Reynolds number presently considered. If small values of the Shields parameter are not considered, the present results suggest that fair predictions of the sediment transport rate can be obtained by looking for a correlation of $q_{s}$ with $\unicode[STIX]{x1D703}$ similar to that which describes $q_{s}$ versus $\unicode[STIX]{x1D703}$ in a steady flow. For example, figure 23 shows that the relationship $q=a(\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{cr})^{b}$ , proposed by Wong & Parker (2006) to evaluate $q_{s}$ in steady flows, with $a=4.93$ , $\unicode[STIX]{x1D703}_{cr}=0.047$ and $b=1.6$ , provides results that agree fairly well with those of the present simulations. When $\unicode[STIX]{x1D703}$ is close to $\unicode[STIX]{x1D703}_{cr}$ and $q_{s}$ is relatively small, it is necessary to take into account that the amount of sediment dragged by the fluid depends also on the time derivative of  $\unicode[STIX]{x1D703}$ . In fact, for the same value of  $\unicode[STIX]{x1D703}$ , $q_{s}$ is larger if $\text{d}\unicode[STIX]{x1D703}/\text{d}t$ is negative. However, this hysteresis is present for such small values of $\unicode[STIX]{x1D703}$ that it can be neglected for practical applications. The fact that there is a rather small but finite sediment transport rate when $\unicode[STIX]{x1D703}$ tends to zero is due to the particle inertia and to the small effects of the imposed pressure gradient. For the same reason, the sediment transport rate tends to vanish for finite values of $\unicode[STIX]{x1D703}$ when particle inertia and the imposed pressure gradient effects balance the viscous force acting on sediment particles.

4 Conclusions

New and interesting information of the sediment transport generated by sea waves is obtained by means of DNS, which allow one to evaluate the hydrodynamics within the oscillatory boundary layer generated by surface waves close to the bottom and to determine the dynamics of idealised sediment particles dragged by the flowing fluid.

The values of the flow Reynolds number fall in the intermittently turbulent regime, such that turbulence is significant only during part of the flow cycle. The other parameters are typical of medium sand. Hence, the results are useful to quantify the bedload sediment transport outside the breaking and surf regions where higher values of the Reynolds number are usually found such that the DNS of the turbulent flow field within the bottom boundary layer are presently unaffordable.

The main result of the investigation is the description of sediment dynamics under the action of the turbulent eddies that are generated within the boundary layer. The pressure fluctuations induced by the turbulent eddies penetrate within the porous bed and generate lift forces that are superimposed onto those due to the pressure difference between the bottom and the top of the sediment particles, which, in turn, is associated with the shear flow close to the bed surface. On average, the lift force due to the turbulent pressure fluctuations is directed upwards and the sediment grains tend to be picked up from the bed and then transported by the external flow in the saltation mode. On the other hand, when the flow relaminarises but the bed shear stress is large enough to induce sediment transport, the sediment grains tend to roll and slide one over the top of the others. This particle dynamics is typical of a laminar flow and it gives rise to sediment transport rates quite different from those observed when turbulence is present.

Figure 24. Dimensionless sediment flow rate $q_{s}$ as a function of the absolute value of the Shields parameter for simulations at $R_{\unicode[STIX]{x1D6FF}}=72,128$ and $d\simeq 0.25$ (solid blue and red lines) (non-turbulent, Mazzuoli et al. 2019) and for the non-turbulent cycles of $\text{run}~1$ at $R_{\unicode[STIX]{x1D6FF}}=450$ and $d=0.335$ (black line). The blue dash-dotted line represents the formula by Wong & Parker (2006). Arrows indicate the orbital direction.

The differences between the values of $q_{s}$ generated by a laminar and a turbulent oscillatory boundary layer can be easily appreciated if the results of figure 23 are compared with those obtained by Mazzuoli et al. (2019), which, for the reader’s convenience, are plotted in figure 24. Mazzuoli et al. (2019) investigated the formation of sea ripples by means of DNS and computed the sediment transport rate for the same values of the parameters as those of some of the laboratory experiments of Blondeaux, Sleath & Vittori (1988). In particular, the experiments characterised by $R_{\unicode[STIX]{x1D6FF}}=72$ and $R_{\unicode[STIX]{x1D6FF}}=128$ and by $d\simeq 0.25$ were considered by Mazzuoli et al. (2019). For such values of the Reynolds number, the flow regime is laminar. In figure 24, the results obtained for $R_{\unicode[STIX]{x1D6FF}}=450$ and $d=0.335$ ( $\text{run}~1$ ) during the half-cycles characterised by weak turbulence are also plotted. As already pointed out, in these cases, the values of $q_{s}$ during the accelerating phases are different from those computed during the decelerating phases, even if the Shields parameter $\unicode[STIX]{x1D703}$ is the same, because particle dynamics is affected not only by the bottom shear stress but also by the streamwise pressure gradient.

In the turbulent regime, the bedload sediment transport rate observed for large values of $\unicode[STIX]{x1D703}$ during the accelerating phases is practically equal to that observed during the decelerating phases because the pressure gradient plays a negligible role in particle dynamics. In fact, the magnitude of the sediment transport rate during the accelerating phases differs from that during the decelerating phases only when the Shields parameter is quite small, the differences being mainly due to the different values of the turbulence intensity observed during the accelerating and decelerating phases, even for the same value of the bottom shear stress. The reader should note that, when the flow regime is laminar, the sediment flow rate decreases if the Reynolds number is increased, while, if turbulence is present, the values of $q_{s}$ increase if the Reynolds number is increased.

Finally, it is worth pointing out that, in the turbulent regime, fair predictions of the bedload sediment flux can be obtained by means of empirical formulae obtained on the basis of experimental measurements carried out in steady flow, at least for the typical periods of sea waves. Indeed, for high values of the Reynolds number, the amplitude of the fluid displacement oscillations turns out to be much larger than the grain size. Hence, the Keulegan–Carpenter number $K_{c}$ of the flow around sediment grains is large and the sediment particles feel a succession of quasi-steady flows.