Skip to main content Accessibility help
×
Home
Hostname: page-component-78dcdb465f-nbrzn Total loading time: 1.394 Render date: 2021-04-18T21:45:07.715Z Has data issue: true Feature Flags: { "shouldUseShareProductTool": true, "shouldUseHypothesis": true, "isUnsiloEnabled": true, "metricsAbstractViews": false, "figures": false, "newCiteModal": false, "newCitedByModal": true }

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

Published online by Cambridge University Press:  03 January 2020

Marco Mazzuoli
Affiliation:
Department of Civil, Chemical and Environmental Engineering, University of Genoa, Via Montallegro 1, 16145 Genova, Italy
Paolo Blondeaux
Affiliation:
Department of Civil, Chemical and Environmental Engineering, University of Genoa, Via Montallegro 1, 16145 Genova, Italy
Giovanna Vittori
Affiliation:
Department of Civil, Chemical and Environmental Engineering, University of Genoa, Via Montallegro 1, 16145 Genova, Italy
Markus Uhlmann
Affiliation:
Institute for Hydromechanics, Karlsruhe Institute of Technology, 76131 Karlsruhe, Germany
Julian Simeonov
Affiliation:
Marine Geosciences Division, Naval Research Laboratory, Stennis Space Center, MS 39529, USA
Joseph Calantoni
Affiliation:
Marine Geosciences Division, Naval Research Laboratory, Stennis Space Center, MS 39529, USA
Corresponding
E-mail address:
Rights & Permissions[Opens in a new window]

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.

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© The Author(s), 2019. Published by Cambridge University Press

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 (Reference Graf1984), Fredsøe & Deigaard (Reference Fredsøe and Deigaard1992), Soulsby (Reference Soulsby1997), Gyr & Hoyer (Reference Gyr and Hoyer2006)). 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 Reference Costamagna, Vittori and Blondeaux2003), 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 Reference Ghodke and Apte2016, Reference Ghodke and Apte2018; Mazzuoli & Vittori Reference Mazzuoli and Vittori2016, Reference Mazzuoli and Vittori2019). 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 Reference Hino, Sawamoto and Takasu1976; Hino et al. Reference Hino, Kashiwayanagi, Nakayama and Hara1983; Jensen, Sumer & Fredsøe Reference Jensen, Sumer and Fredsøe1989; Akhavan, Kamm & Shapiro Reference Akhavan, Kamm and Shapiro1991; Carstensen, Sumer & Fredsøe Reference Carstensen, Sumer and Fredsøe2010). Later, the experimental observations found a theoretical interpretation by Wu (Reference Wu1992) and Blondeaux & Vittori (Reference Blondeaux and Vittori1994), 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. Reference Akhavan, Kamm and Shapiro1991; Verzicco & Vittori Reference Verzicco and Vittori1996; Costamagna et al. Reference Costamagna, Vittori and Blondeaux2003; Ozdemir, Hsu & Balachandar Reference Ozdemir, Hsu and Balachandar2014). 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. Reference Jensen, Sumer and Fredsøe1989; Carstensen, Sumer & Fredsøe Reference Carstensen, Sumer and Fredsøe2012; Mazzuoli & Vittori Reference Mazzuoli and Vittori2019).

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 (Reference Mazzuoli, Kidanemariam and Uhlmann2019) 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. (Reference Mazzuoli, Kidanemariam and Uhlmann2019) 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. (Reference Lobkovsky, Orpe, Molloy, Kudrolli and Rothman2008), 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. (Reference Mazzuoli, Kidanemariam and Uhlmann2019) 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 Reference Mazzuoli and Vittori2016). The effect of the transition to turbulence can be observed in figure 1, where the bottom shear stress measured by Jensen et al. (Reference Jensen, Sumer and Fredsøe1989) 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 (Reference Blondeaux, Vittori and Porcile2018), 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 Reference Ghodke and Apte2016; Mazzuoli & Vittori Reference Mazzuoli and Vittori2016) 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. Reference Clark, Shattuck, Ouellette and O’Hern2017).

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. (Reference Jensen, Sumer and Fredsøe1989) 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. Reference van der A, O’Donoghue, Davies and Ribberink2011; Scandura, Faraci & Foti Reference Scandura, Faraci and Foti2016), 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 (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:

(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 (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).

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 (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

(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 (Reference Blondeaux, Vittori and Mazzuoli2016), 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 (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

(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 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 dg). 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 (ad) 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.

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. Reference Mazzuoli, Kidanemariam and Uhlmann2019) 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 (Reference Wong and Parker2006). 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. (Reference Mazzuoli, Kidanemariam and Uhlmann2019), which, for the reader’s convenience, are plotted in figure 24. Mazzuoli et al. (Reference Mazzuoli, Kidanemariam and Uhlmann2019) 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 (Reference Blondeaux, Sleath and Vittori1988). 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. (Reference Mazzuoli, Kidanemariam and Uhlmann2019). 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.

Acknowledgements

This study has been funded by the Office of Naval Research (USA) (under the research project no. N62909-17-1-2144). The support of CINECA, which provided computational resources on Marconi under the PRACE project MOST SEA (Proposal ID: 2017174199), and a grant of computer time from the DoD High Performance Computing Modernization Program at the ERDC DSRC are also acknowledged. We also acknowledge that this research was also supported by the MIUR, in the framework of PRIN 2017, with grant no. 20172B7MY9, project ‘FUNdamentals of BREAKing wave-induced boundary dynamics’.

Declaration of interests

The authors report no conflict of interest.

Supplementary material and movies

Supplementary material and movies are available at https://doi.org/10.1017/jfm.2019.1012.

References

van der A, D. A., O’Donoghue, T., Davies, A. G. & Ribberink, J. S. 2011 Experimental study of the turbulent boundary layer in acceleration-skewed oscillatory flow. J. Fluid Mech. 684, 251283.CrossRefGoogle Scholar
Akhavan, R., Kamm, R. D. & Shapiro, A. H. 1991 An investigation of transition to turbulence in bounded oscillatory Stokes flows. Part 1. Experiments. J. Fluid Mech. 225, 395422.CrossRefGoogle Scholar
Bagnold, R. A.1966 An approach to the sediment transport for general physics. Geological Survey Professional Papers 422-I.CrossRefGoogle Scholar
Bettencourt, J. H. & Dias, F. 2018 Wall pressure and vorticity in the intermittently turbulent regime of the Stokes boundary layer. J. Fluid Mech. 851, 479506.CrossRefGoogle Scholar
Blennerhassett, P. J. & Bassom, A. P. 2002 The linear stability of flat Stokes layers. J. Fluid Mech. 464, 393410.CrossRefGoogle Scholar
Blondeaux, P. & Seminara, G. 1979 Transizione incipiente al fondo di unonda di gravitá. Acc. Naz. Lincei 67, 408411.Google Scholar
Blondeaux, P., Sleath, J. F. A. & Vittori, G. 1988 Experimental data on sand ripples in an oscillatory flow. Rep. 1, 88.Google Scholar
Blondeaux, P. & Vittori, G. 1991 A route to chaos in an oscillatory flow: Feigenbaum scenario. Phys. Fluids 3 (11), 24922495.CrossRefGoogle Scholar
Blondeaux, P. & Vittori, G. 1994 Wall imperfections as a triggering mechanism for Stokes layer transition. J. Fluid Mech. 264, 107135.CrossRefGoogle Scholar
Blondeaux, P., Vittori, G. & Mazzuoli, M. 2016 Pattern formation in a thin layer of sediment. Mar. Geol. 376, 3950.CrossRefGoogle Scholar
Blondeaux, P., Vittori, G. & Porcile, G. 2018 Modeling the turbulent boundary layer at the bottom of sea wave. Coast. Engng 141, 1223.CrossRefGoogle Scholar
Carstensen, S., Sumer, B. M. & Fredsøe, J. 2010 Coherent structures in wave boundary layers. Part 1. Oscillatory motion. J. Fluid Mech. 646, 169206.CrossRefGoogle Scholar
Carstensen, S., Sumer, B. M. & Fredsøe, J. 2012 A note on turbulent spots over a rough bed in wave boundary layers. Phys. Fluids 24 (11), 115104.CrossRefGoogle Scholar
Clark, A. H., Shattuck, M. D., Ouellette, N. T. & O’Hern, C. S. 2017 Role of grain dynamics in determining the onset of sediment transport. Phys. Rev. Fluids 2 (3), 034305.CrossRefGoogle Scholar
Costamagna, P., Vittori, G. & Blondeaux, P. 2003 Coherent structures in oscillatory boundary layers. J. Fluid Mech. 474, 133.CrossRefGoogle Scholar
Fredsøe, J. & Deigaard, R. 1992 Mechanics of Coastal Sediment Transport. World Scientific.CrossRefGoogle Scholar
Ghodke, C. D. & Apte, S. V. 2016 Dns study of particle-bed–turbulence interactions in an oscillatory wall-bounded flow. J. Fluid Mech. 792, 232251.CrossRefGoogle Scholar
Ghodke, C. D. & Apte, S. V. 2018 Roughness effects on the second-order turbulence statistics in oscillatory flows. Comput. Fluids 162, 160170.CrossRefGoogle Scholar
Graf, W. H. 1984 Hydraulics of Sediment Transport. Water Resources Publications, LLC.Google Scholar
Gyr, A. & Hoyer, K. 2006 Sediment Transport. Springer.Google Scholar
Hino, M., Kashiwayanagi, M., Nakayama, A. & Hara, T. 1983 Experiments on the turbulence statistics and the structure of a reciprocating oscillatory flow. J. Fluid Mech. 131, 363400.CrossRefGoogle Scholar
Hino, M., Sawamoto, M. & Takasu, S. 1976 Experiments on transition to turbulence in an oscillatory pipe flow. J. Fluid Mech. 75 (2), 193207.CrossRefGoogle Scholar
Jensen, B. L., Sumer, B. M. & Fredsøe, J. 1989 Turbulent oscillatory boundary layers at high Reynolds numbers. J. Fluid Mech. 206, 265297.CrossRefGoogle Scholar
Jeong, J. & Hussain, F. 1995 On the identification of a vortex. J. Fluid Mech. 285, 6994.CrossRefGoogle Scholar
Kajiura, K. 1968 A model of the bottom boundary layer in water waves. Bull. Earthq. Res. Inst. Univ. Tokyo 46, 75123.Google Scholar
Kidanemariam, A. G. & Uhlmann, M. 2014a Direct numerical simulation of pattern formation in subaqueous sediment. J. Fluid Mech. 750, R2.CrossRefGoogle Scholar
Kidanemariam, A. G. & Uhlmann, M. 2014b Interface-resolved direct numerical simulation of the erosion of a sediment bed sheared by laminar channel flow. Intl J. Multiphase Flow 67, 174188.CrossRefGoogle Scholar
Kidanemariam, A. G. & Uhlmann, M. 2017 Formation of sediment patterns in channel flow: minimal unstable systems and their temporal evolution. J. Fluid Mech. 818, 716743.CrossRefGoogle Scholar
Lobkovsky, A. E., Orpe, A. V., Molloy, R., Kudrolli, A. & Rothman, D. H. 2008 Erosion of a granular bed driven by laminar fluid flow. J. Fluid Mech. 605, 4758.CrossRefGoogle Scholar
Mazzuoli, M., Blondeaux, P., Simeonov, J. & Calantoni, J. 2018 Direct numerical simulation of oscillatory flow over a wavy, rough, and permeable bottom. J. Geophys. Res. Oceans 123 (3), 15951611.CrossRefGoogle Scholar
Mazzuoli, M., Kidanemariam, A. G., Blondeaux, P., Vittori, G. & Uhlmann, M. 2016 On the formation of sediment chains in an oscillatory boundary layer. J. Fluid Mech. 789, 461480.CrossRefGoogle Scholar
Mazzuoli, M., Kidanemariam, A. G. & Uhlmann, M. 2019 Direct numerical simulations of ripples in an oscillatory flow. J. Fluid Mech. 863, 572600.CrossRefGoogle Scholar
Mazzuoli, M. & Vittori, G. 2016 Transition to turbulence in an oscillatory flow over a rough wall. J. Fluid Mech. 792, 6797.CrossRefGoogle Scholar
Mazzuoli, M. & Vittori, G. 2019 Turbulent spots in an oscillatory flow over a rough wall. Eur. J. Mech. (B/Fluids) 78, 161168.CrossRefGoogle Scholar
Mazzuoli, M., Vittori, G. & Blondeaux, P. 2011a Turbulent spots in a Stokes boundary layer. J. Phys.: Conf. Ser. 318, 032032.Google Scholar
Mazzuoli, M., Vittori, G. & Blondeaux, P. 2011b Turbulent spots in oscillatory boundary layers. J. Fluid Mech. 685, 365376.CrossRefGoogle Scholar
Ozdemir, C. E., Hsu, T.-J. & Balachandar, S. 2014 Direct numerical simulations of transition and turbulence in smooth-walled Stokes boundary layer. Phys. Fluids 26 (4), 045108.CrossRefGoogle Scholar
Scandura, P., Faraci, C. & Foti, E. 2016 A numerical investigation of acceleration-skewed oscillatory flows. J. Fluid Mech. 808, 576613.CrossRefGoogle Scholar
Schlichting, H. 1936 Experimentelle untersuchungen zum rauhigkeitsproblem. Arch. Appl. Mech. 7 (1), 134.Google Scholar
Sleath, J. F. A. 1988 Transition in oscillatory flow over rough beds. J. Waterways Port Coast. Ocean Engng 114 (1), 1833.CrossRefGoogle Scholar
Soulsby, R. 1997 Dynamics of Marine Sands. Thomas Telford.Google Scholar
Sumer, B. M. & Fredsøe, J. 2002 The Mechanics of Scour in the Marine Environment. World Scientific.CrossRefGoogle Scholar
Uhlmann, M. 2005 An immersed boundary method with direct forcing for the simulation of particulate flows. J. Comput. Phys. 209 (2), 448476.CrossRefGoogle Scholar
Uhlmann, M. 2008 Interface-resolved direct numerical simulation of vertical particulate channel flow in the turbulent regime. Phys. Fluids 20 (5), 053305.CrossRefGoogle Scholar
Verzicco, R. & Vittori, G. 1996 Direct simulation of transition in Stokes boundary layers. Phys. Fluids 8 (6), 13411343.CrossRefGoogle Scholar
Vittori, G. 2003 Sediment suspension due to waves. J. Geophys. Res. 108 (C6), 31733189.CrossRefGoogle Scholar
Vittori, G. & Verzicco, R. 1998 Direct simulation of transition in an oscillatory boundary layer. J. Fluid Mech. 371, 207232.CrossRefGoogle Scholar
Wong, M. & Parker, G. 2006 Reanalysis and correction of bed-load relation of Meyer-Peter and Müller using their own database. ASCE J. Hydraul. Engng 132 (11), 11591168.CrossRefGoogle Scholar
Wu, X. 1992 The nonlinear evolution of high-frequency resonant-triad waves in an oscillatory Stokes layer at high Reynolds number. J. Fluid Mech. 245, 553597.CrossRefGoogle Scholar

Mazzuoli et al. supplementary movie 1

Streamwise velocity fluctuations at a plane above the mean bottom elevation. Run 2 (Rδ=750, d/δ=0.335).

Video 11 MB

Mazzuoli et al. supplementary movie 2

Spanwise vorticity fluctuations at a plane above the mean bottom elevation. Run 2 (Rδ=750, d/δ=0.335).