2.1 The fluid motion

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

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

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,

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

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 (Reference Kidanemariam and Uhlmann2014a , Reference Kidanemariam and Uhlmann2017) and Mazzuoli et al. (Reference Mazzuoli, Kidanemariam, Blondeaux, Vittori and Uhlmann2016, Reference Mazzuoli, Kidanemariam and Uhlmann2019). 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 Reference Uhlmann2005), 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. (Reference Mazzuoli, Kidanemariam, Blondeaux, Vittori and Uhlmann2016)).

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:

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 (Reference Kidanemariam and Uhlmann2014b ).

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. (Reference Mazzuoli, Kidanemariam, Blondeaux, Vittori and Uhlmann2016) (see table 3).

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 (Reference Kidanemariam and Uhlmann2014b ). 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

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

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.1 Appearance of turbulence and turbulent kinetic energy

As discussed in Sleath (Reference Sleath1988), on the basis of Kajiura’s (Reference Kajiura1968) 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 (Reference Kajiura1968) 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 Reference Sleath1988) 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 (Reference Verzicco and Vittori1996) and Vittori & Verzicco (Reference Vittori and Verzicco1998) 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. (Reference Mazzuoli, Kidanemariam and Uhlmann2019) 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

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 Reference Uhlmann2008; Mazzuoli et al. Reference Mazzuoli, Blondeaux, Simeonov and Calantoni2018): (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. (Reference Mazzuoli, Blondeaux, Simeonov and Calantoni2018, Reference Mazzuoli, Kidanemariam and Uhlmann2019). The time development of $\unicode[STIX]{x1D70F}_{b}$ for the $\text{run}~2~\text{(fix)}$ qualitatively agrees with that measured by Jensen et al. (Reference Jensen, Sumer and Fredsøe1989) (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 Reference Mazzuoli, Vittori and Blondeaux2011b ).

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 (Reference Mazzuoli, Vittori and Blondeaux2011a ) and Mazzuoli et al. (Reference Mazzuoli, Vittori and Blondeaux2011b ), 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. (Reference Mazzuoli, Vittori and Blondeaux2011a ,Reference Mazzuoli, Vittori and Blondeaux 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. (Reference Mazzuoli, Vittori and Blondeaux2011a ,Reference Mazzuoli, Vittori and Blondeaux 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 (Reference Schlichting1936), 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 (Reference Vittori2003), 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. Reference Clark, Shattuck, Ouellette and O’Hern2017) and also because of the effects of the imposed pressure gradient (Mazzuoli et al. Reference Mazzuoli, Kidanemariam and Uhlmann2019). 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. (Reference Costamagna, Vittori and Blondeaux2003), Mazzuoli et al. (Reference Mazzuoli, Vittori and Blondeaux2011b ) and Mazzuoli & Vittori (Reference Mazzuoli and Vittori2016, Reference Mazzuoli and Vittori2019), 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. (Reference Carstensen, Sumer and Fredsøe2010), who made flow visualisations for values of the Reynolds number close to its critical value. Indeed, Kajiura’s (Reference Kajiura1968) 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. (Reference Carstensen, Sumer and Fredsøe2010) during laboratory experiments and by Mazzuoli et al. (Reference Mazzuoli, Vittori and Blondeaux2011a ,Reference Mazzuoli, Vittori and Blondeaux b ), who reproduced the experiments by Carstensen et al. (Reference Carstensen, Sumer and Fredsøe2010) 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 Reference Blondeaux and Vittori1991), the present numerical findings show vortex structures that are qualitatively in agreement with both the predictions of the linear stability analysis of Blondeaux & Seminara (Reference Blondeaux and Seminara1979) and the results of the DNS of Vittori & Verzicco (Reference Vittori and Verzicco1998), Costamagna et al. (Reference Costamagna, Vittori and Blondeaux2003) and Bettencourt & Dias (Reference Bettencourt and Dias2018), who considered a smooth bottom. Blondeaux & Seminara (Reference Blondeaux and Seminara1979), 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 (Reference Blennerhassett and Bassom2002) 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 Reference Blondeaux and Vittori1994) 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. (Reference Carstensen, Sumer and Fredsøe2010) 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 Reference Jeong and Hussain1995). 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 Reference Bagnold1966; Sumer & Fredsøe Reference Sumer and Fredsøe2002). 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 (Reference Wong and Parker2006). 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 (Reference Wong and Parker2006) 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.