On the role of turbulent large-scale streaks in generating sediment ridges

Abstract

The role of turbulent large-scale streaks or large-scale motions in forming subaqueous sediment ridges on an initially flat sediment bed is investigated with the aid of particle resolved direct numerical simulations of open channel flow at bulk Reynolds numbers up to 9500. The regular arrangement of quasi-streamwise ridges and troughs at a characteristic spanwise spacing between 1 and 1.5 times the mean fluid height is found to be a consequence of the spanwise organisation of turbulence in large-scale streamwise velocity streaks. Ridges predominantly appear in regions of weaker erosion below large-scale low-speed streaks and vice versa for troughs. The interaction between the dynamics of the large-scale streaks in the bulk flow and the evolution of sediment ridges on the sediment bed is best described as ‘top-down’ process, as the arrangement of the sediment bedforms is seen to adapt to changes in the outer flow with a time delay of several bulk time units. The observed ‘top-down’ interaction between the outer flow and the bed agrees fairly well with the conceptual model on causality in canonical channel flows proposed by Jiménez (J. Fluid Mech., vol. 842, 2018, P1, § 5.6). Mean secondary currents of Prandtl's second kind of comparable intensity and lateral spacing are found over developed sediment ridges and in single-phase smooth-wall channels alike in averages over ${O}(10)$ bulk time units. This indicates that the secondary flow commonly observed together with sediment ridges is the statistical footprint of the regularly organised large-scale streaks.

1. Introduction

Classically, secondary currents have been studied in flows with non-circular cross-section featuring smooth walls (Nikuradse 1926; Prandtl 1926), where they are generated in the region close to the sidewalls. It was also Prandtl who proposed the nowadays standard classification of secondary flows into two categories: secondary flows of Prandtl's first kind are due to the skewness of the mean flow axis as in meandering rivers or curved pipes, while secondary flows of the second kind are a pure turbulent phenomenon and originate in anisotropy and non-homogeneity of Reynolds stresses across the flow domain (Bradshaw 1987). Our current study focusses on the second family of turbulence-induced secondary currents exclusively and in the remainder of the text, we will always refer to this type when using the term secondary flow.

In hydraulic flows, secondary flow cells mutually interact with the mobile sediment bed, giving rise to long, essentially streamwise-aligned sediment ridges (cf. figure 1). Here, the transverse mean flow manifests itself as pairs of quasi-streamwise depth-spanning and counter-rotating flow cells on each side of a ridge, with a transverse spacing of $(1-2){H_f}$ (Nezu & Nakagawa 1993), where ${H_f}$ is the mean fluid height. The interplay between secondary currents and sediment ridges was first claimed in the context of early laboratory studies (Casey 1935; Vanoni 1946; Wolman & Brush 1961) and field observations of dried river beds (Culbertson 1967; Karcz 1967), but there was disagreement on the actual origin of the process. While it was first discussed that the evolution of sediment ridges could be a direct consequence of the secondary flow pattern induced by lateral sidewalls (Nezu & Nakagawa 1984), it was later discovered that ridges and secondary currents were simultaneously rising at different spanwise locations of wide experimental flumes (Nezu & Nakagawa 1989). These observations were in agreement with the conceptual idea of Ikeda (1981), who conjectured the underlying mechanism to be a flow–sediment bed instability that leads to the evolution of both sediment ridges and secondary flow cells, and that is generally independent of the presence of lateral sidewalls. These considerations were confirmed by Colombini (1993), who showed by means of linear stability analysis that a two-dimensional turbulent base flow in an infinitely wide channel is unstable with respect to infinitesimal spanwise disturbances of the sediment bed, with a most amplified wavelength of $1.3{H_f}$ for the considered parameter set. Later, Colombini & Parker (1995) extended the model by including the effect of laterally varying roughness and sand grain distributions as a second source of instability, in agreement with earlier experimental studies (Müller & Studerus 1979).

Figure 1. Conceptual sketch of sediment ridges in a turbulent channel flow with a mean lateral pattern spacing ${\lambda _{h,z}}$.

Outside the hydraulic community, interest in secondary currents has increased mainly during the past two decades, motivated by observations of secondary flows over complex three-dimensional industrial surfaces such as the replicated surface of a damaged turbine blade investigated by Mejia-Alvarez & Christensen (2010) and Barros & Christensen (2014). Secondary flows seem to represent a robust phenomenon that is ubiquitous in flow configurations that feature a significant lateral heterogeneity, including straight (Goldstein & Tuan 1998) and converging/diverging riblets (Kevin et al. 2017; Kevin, Monty & Hutchins 2019b), spanwise alternating roughness stripes (Hinze 1967,  1973; McLelland et al. 1999; Wang & Cheng 2005; Wangsawijaya et al. 2020), transverse alternating non-/hydrophobic surfaces (Türk et al. 2014; Stroh et al. 2016), streamwise-aligned artificial ridges of different cross-sectional shape (Wang & Cheng 2006; Vanderwel & Ganapathisubramani 2015; Vanderwel et al. 2019; Medjnoun, Vanderwel & Ganapathisubramani 2020; Stroh et al. 2020a; Zampiron, Cameron & Nikora 2020) or combinations of the above (Stroh et al. 2020b). The majority of the listed studies predominantly focus on the general structure of the mean secondary currents, their sense of rotation (Yang & Anderson 2018; Anderson 2019) or the influence of the lateral spacing of the heterogeneities on their intensity and arrangement (Chung, Monty & Hutchins 2018; Wangsawijaya et al. 2020).

Less attention has been given, on the other hand, to the relation between mean secondary flows and the dynamics of instantaneous coherent structures. In canonical laterally homogeneous channels, for instance, instantaneous quasi-streamwise vortices appear at all scales in connection with alternating streaks of high and low streamwise velocity (Jiménez 1998,  2018). The smallest vortices of this kind are those in the buffer layer that form a ‘self-sustaining process’ with the near-wall velocity streaks (Hamilton, Kim & Waleffe 1995; Jeong et al. 1997; Schoppa & Hussain 2002). With increasing distance to the wall, however, energy and enstrophy no longer reside in the same scales: velocity streaks are found to scale self-similarly across the logarithmic layer with the distance to the wall (Jiménez 2013b), whereas vorticity concentrates further away from the wall in structures of essentially the same scale as in the buffer layer (Jiménez 2013a). It is thus evident that larger-scale quasi-streamwise rotating motions that sustain the fully turbulent log-layer streaks cannot be attributed to a single large-scale vortex. Rather, they represent the collective effect of a large number of locally isotropic small-scale vortices that globally arrange in large-scale vortex clusters that are sufficiently anisotropic to generate an average upward motion in their centre (Del Álamo et al. 2006; Jiménez 2018). While these collective large-scale rotating motions are hardly detectable in instantaneous flow visualizations, they have been successfully uncovered in conditional averages (Del Álamo et al. 2006; Lozano-Durán, Flores & Jiménez 2012) or low-pass filtered fields (Motoori & Goto 2019), proving that the self-sustaining process between streaks and quasi-streamwise vortical motions is a self-similar cascade (Motoori & Goto 2021), just as the logarithmic layer is (Flores & Jiménez 2010; Kevin, Monty & Hutchins 2019a).

While there is a uniform terminology for the coherent structures in the buffer layer of canonical wall-bounded flows, different expressions exist for the large-scale coherent structures. In some studies, the term ‘large-scale streaks’ is preferred in order to highlight that they are the largest representatives of a cascade of self-similar streaks starting with the smooth buffer-layer streaks discussed by Kline et al. (1967) (cf., for instance, Jiménez 2018), whereas other authors refer to them as ‘large-scale motions’ (LSMs) (Smits, McKeon & Marusic 2011) to clearly express their different scaling behaviour in terms of the outer length scale, as opposed to the inner-scaling buffer-layer streaks. Also, the terms ‘high- and low-momentum regions’ (HMR/LMR) are established for large-scale spanwise alternating regions of relatively higher and lower streamwise velocity (Wu & Christensen 2010), often in the context of organized packets of hairpin vortices (Ganapathisubramani, Longmire & Marusic 2003). These HMRs and LMRs should not be confused with the terms ‘high- and low-momentum pathways’ (HMP/LMP) which are frequently used in the context of mean secondary currents for the laterally varying streamwise velocity in a time- and streamwise-averaged frame. Rather, one might think of these HMPs/LMPs as preferential pathways of the formerly discussed instantaneous three-dimensional large-scale streaks, or LMRs/HMRs (Mejia-Alvarez, Barros & Christensen 2013). In the remainder of the current study, we choose the first approach and refer to streamwise-elongated coherent motions whose dimensions are ${O}({H_f})$ as ‘large-scale streaks’, emphasizing, however, that these structures do not differ from the elsewhere discussed LSMs or HMRs/LMRs.

The clear connection between instantaneous coherent structures and the mean secondary flow was deduced by Uhlmann et al. (2007) and Pinelli et al. (2010) for the characteristic eight-vortex secondary flow pattern in a square duct. They showed that the preferential location of the buffer-layer streaks and quasi-streamwise vortices along the four solid sidewalls determines the statistics of the mean vorticity distribution, which is directly linked to the mean secondary flow (Pinelli et al. 2010; Pirozzoli et al. 2018). Sakai (2016) later elaborated that similar relations between coherent structures and the mean secondary flow also exist in open duct flows, in which the top wall is replaced by a free-slip boundary. The close relation between distinct turbulent coherent structures and the mean secondary flow was further supported by Uhlmann, Kawahara & Pinelli (2010), who presented a family of travelling wave solutions whose members produce essentially the same secondary flow pattern as the fully turbulent flow when averaged over the streamwise direction (Kawahara, Uhlmann & van Veen 2012). Also, it was recently concluded by Kevin et al. (2019a) that the mean secondary currents over spanwise heterogeneous bottom walls are the statistical footprint of the conditional or collective quasi-streamwise rollers observed in canonical wall-bounded flows (Lozano-Durán et al. 2012). The authors argue that the spanwise bottom heterogeneities cause a lateral locking of the large-scale structures and quasi-streamwise rollers in the vicinity of the roughness transition such that they get visible in the long-time average, in contrast to the flow over homogeneous walls where instantaneous flow structures will appear at any spanwise position with the same probability such that no mean secondary flow is detectable in the long-time average. Spatial locking should be understood in this context in a sense that the mean lateral position of the coherent structures does not vary significantly in time, while the log- and outer-layer streaks instantaneously meander in response to the quasi-streamwise rollers just as in canonical flows (Kevin et al. 2017,  2019a). Indeed, such lateral meandering has been detected for a variety of the above discussed spanwise heterogeneities, with the lateral spacing of the bottom wall heterogeneities serving as a control parameter for the organization of the flow, especially for the spacing and the location of the mean secondary currents (Kevin et al. 2017; Wangsawijaya et al. 2020; Zampiron et al. 2020). Interestingly, the findings reported in these latter studies suggest an optimal lateral spacing of the bottom roughness (optimal referring to a maximization of the secondary currents’ intensity) corresponds to the typical lateral dimension of the large-scale streaks in canonical wall-bounded shear flows of $(1-2){H_f}$ (Jiménez 2013b,  2018). For spacings clearly smaller or larger than the outer flow scale, on the other hand, the effect of the heterogeneities is confined to the near-wall region or the regions of roughness transition, respectively, while there is little effect on the outer flow. In the hydraulic context, sediment ridges similarly form at a lateral spacing of $(1-2){H_f}$ (Colombini 1993; Nezu & Nakagawa 1993), consistent with the lateral dimension of the turbulent large-scale structures with streamwise length ${O}(1-10){H_f}$ (Tamburrino & Gulliver 1999,  2007; Shvidchenko & Pender 2001; Roy et al. 2004), that have been detected as counterparts to the large-scale and very large-scale flow structures inherent to canonical turbulent flows (Adrian & Marusic 2012). A general relevance of these instantaneous structures for the formation of sediment ridges was formulated, for instance, by Nezu (2005) based on the smooth-wall duct flow experiments of Nezu & Rodi (1985). It was discovered that the mean secondary currents decay towards the duct centre, while instantaneous rotating motions remain observable. Nezu (2005) therefore claimed that the ‘instantaneous secondary currents may trigger the formation of smooth/rough striping of the bed and the development of sand ribbons.’ Shvidchenko & Pender (2001) further proposed that passing ‘macroturbulent structures’ with a streamwise spacing of $(4-5){H_f}$ induce laterally alternating regions of high and low sediment erosion and transport inside and outside their travelling path, respectively. The lateral variation of sediment transport is due to the fact that sediment erosion is predominantly driven through the action of large-scale high-speed sweeps reaching the bed, rather than by ejections lifting it up in the low-speed regions (Gyr & Schmid 1997; Cameron, Nikora & Witz 2020).

Even though it appears conclusive that organized recurrent large-scale structures are responsible for the formation of sediment ridges, to the best of our knowledge, clear evidence of the above elaborated mechanisms is still lacking, also because simultaneous measurements of three-dimensional flow structures and individual sediment motion remain challenging in experiments up to the present day (Wang & Cheng 2005). In the past decade, on the other hand, direct numerical simulations including fully resolved finite size particles have become a useful and feasible tool to overcome such limitations and have been proven to provide valuable insight into the morphodynamics of subaqueous sediment beds (Kidanemariam & Uhlmann 2014a,  2017; Vowinckel et al. 2017; Scherer, Kidanemariam & Uhlmann 2020).

In the current study, we will present a series of such fully resolved particle-laden simulations that provide numerical evidence that the organization of the streamwise velocity in large-scale high- and low-speed streaks directly causes a laterally varying erosion activity which, in turn, leads to the evolution of laterally organized sediment ridges. It will be shown that these large-scale flow structures that centre on the core of the channel interact in a kind of top–down mechanism with the flow in the region close to the bed as well as with the sediment bed itself, in good agreement with the conceptual model of Jiménez (2018, § 5.6 and references therein) on the organization of coherent structures in canonical channel flows. Eventually, we will close the loop by showing that the mean secondary currents usually attributed to sediment ridges are the statistical footprint of the aforementioned turbulent large-scale streaks and their associated Reynolds stress-carrying structures.

This manuscript is structured as follows: in the following § 2, we will first give a brief overview of the numerical scheme used in the current simulations, while details on the physical system under consideration will be given in § 3. Section 4 summarizes the applied procedure to extract the fluid–bed interface, which shall form the basis for the analysis of the simulation results in § 5. The observed ‘top–down mechanism’ between large-scale structures, the near-bed flow organization and sediment ridges will be discussed in § 6, before we close with a summary of the relevant outcomes in § 7.

2. Numerical method

The numerical scheme used in the current study was first used for the simulation of sediment transport in a carrying fluid in Kidanemariam & Uhlmann (2014b). Its main components are an immersed boundary technique based on the concept proposed by Uhlmann (2005) and a particle contact model whose features will be discussed below. The time evolution of the fluid phase is governed by the continuity and Navier–Stokes equations for incompressible Newtonian fluids, extended by an additional force term emanating from the immersed boundary formulation which enforces the no-slip boundary conditions at the phase boundaries at each particle's surface. This technique allows us then to discretize the physical domain by a fixed, uniform staggered finite difference grid, while for the numerical time integration, a mixed explicit–implicit framework is chosen, consisting of a Crank–Nicholson scheme for the diffusive terms and a three-step low storage Runge–Kutta scheme for the convective part of the equations. Within a fractional step method, the governing equations are first solved by disregarding the continuity equation, then projecting the velocity to the space of solenoidal velocity fields. Back-and-forth transformations of physical information between the global Eulerian fluid grid nodes and the Lagrangian marker points that represent the surface of the particles in the context of the immersed boundary method are performed based on regularized discrete delta functions (Peskin 2002).

The dynamics of the spherical particles follows the Newton–Euler equations for rigid body motion and are integrated in time alongside the fluid motion, using a sub-stepping procedure with ${O}(100)$ particle sub-steps per fluid time step to take care of the different characteristic time scales of fluid and particle motion (Kidanemariam 2015). The exchange of linear and angular momentum during particle–wall and inter-particle contacts, with the latter occurring frequently in bedload-dominated sedimentary flows, are modelled by a soft-sphere discrete element model that describes the interaction of solid objects by the analogy to simple mechanical mass–spring–damper systems. The individual force and torque components have a finite range of action, that is, two objects are in contact (i.e. they feel the contact force and torque) if and only if the minimal distance between their surfaces ${\varDelta }$ is below this force range ${\varDelta _c}$. In its current form, the collision force and torque include three components, namely a normal elastic, a normal damping and a tangential damping contribution. The normal elastic component is a linear function of the overlap length ${\delta _c}={\varDelta _c}-{\varDelta }$ with a constant stiffness coefficient ${k_n}$. The normal damping component is a linear function of the normal component of the relative particle velocity between the two particles at the contact point, with a constant proportionality coefficient ${c_n}$. In a similar fashion, the tangential damping component linearly depends on the relative tangential particle velocity at the contact point, premultiplied by a constant tangential friction coefficient ${c_t}$. Note that the latter force component has a natural upper traction limit in the Coulomb friction coefficient ${\mu _c}$. For a more detailed presentation of the described model and an extensive validation study, the interested reader is referred to Kidanemariam & Uhlmann (2014b).

The contact model thus comes with a set of five unknown parameters, that are, the four force parameters (${k_n}$, ${c_n}$, ${c_t}$, ${\mu _c}$) plus the force range ${\varDelta _c}$. As pointed out by Kidanemariam & Uhlmann (2014b), a quantity that is often more accessible than the force parameters is the dry restitution coefficient ${\varepsilon _d}$ which represents the absolute value of the ratio between the normal components of the relative velocity before and after a dry collision. It is therefore appropriate to use this quantity to eliminate the normal damping coefficient ${c_n}$ from the set of unknowns and to determine it as a function of ${\varepsilon _d}$ and ${k_n}$, such that the actually used set of model force parameters reads (${k_n}$, ${\varepsilon _d}$, ${c_t}$, ${\mu _c}$). A parameter sensitivity analysis has been recently presented in Scherer et al. (2020), outlining the influence of several of the parameters on the simulation results. For the simulations carried out in the present work, the normal stiffness coefficient is set in a range ${k_n}=17\,000-40\,000$ times the submerged weight of a single particle divided by the particle diameter $D$, the Coulomb friction limit is chosen in a range ${\mu _c}=0.4-0.5$ and the friction coefficient is set to ${c_t}={c_n}$. The force range ${\varDelta _c}$ is set equal to the uniform grid spacing ${\Delta x}$. The dry restitution coefficient is ${\varepsilon _d}=0.3$ except for case M850, where we have used a larger value ${\varepsilon _d}=0.9$. Note that the variation of ${\varepsilon _d}$ has no significant effect on the bed formation in the currently investigated bedload-dominated regime, as has been reported in the aforementioned sensitivity analysis.

3. Flow configuration and parameter values

In the course of the current study, we have conducted four individual simulations of turbulent open channel flow past a mobile sediment bed formed out of spherical sediment particles. The set of simulations is complemented by two additional reference simulations of single-phase smooth-wall open channel flow performed with a pseudo-spectral code using Fourier and Chebyshev expansions in the periodic and wall-normal directions, respectively (Kim, Moin & Moser 1987). The physical system simulated in the two-phase cases is sketched in figure 2. The Cartesian coordinate system has its origin on the bottom wall of the channel, such that the components of an arbitrary spatial position vector ${\boldsymbol {x}}=(x,y,z)^{\textrm {T}}$ are measured along the streamwise ($x$), wall-normal ($y$) and spanwise ($z$) directions, respectively. Accordingly, the fluid velocity vector field at position ${\boldsymbol {x}}$ has components ${\boldsymbol {u}_f}({\boldsymbol {x}})=({u_f},{v_f},{w_f})^\textrm {T}$. Subscripted letters ‘$f$’ and ‘$p$’ indicate fluid- and particle-related physical quantities. For the statistical analysis, we introduce the following averaging scheme

(3.1a)\begin{align} {\boldsymbol{u}_f} &= {{\langle {\boldsymbol{u}_f} \rangle_{t}}} + {\overline{\boldsymbol{u}}}_f, \end{align}

(3.1b)\begin{align} {\boldsymbol{u}_f} &= {{\langle {\boldsymbol{u}_f} \rangle_{xz}}} + {\boldsymbol{u}_f^{\prime}}, \end{align}

where (3.1a) and (3.1b) describe the decompositions with respect to time average and plane average, respectively. Averaging in the homogeneous directions and/or time is indicated by angular brackets with the respective indices, $\widebar {\bullet }$ and $\bullet ^{\prime }$ are the corresponding temporal and spatial fluctuations. We further introduce the following decomposition based upon streamwise averaging only:

(3.2)\begin{equation} {{\langle {\boldsymbol{u}_f} \rangle_{x}}} = {{\langle {\boldsymbol{u}_f} \rangle_{xz}}} + {\boldsymbol{u}_f^{\prime\prime}}, \end{equation}

where ${\boldsymbol {u}_f^{\prime \prime }}={\langle {{\boldsymbol {u}_f^{\prime }}} \rangle _{x}}$ are the spatial fluctuations with respect to the streamwise average.

Figure 2. Sketch of the physical system analysed in the multiphase simulations. Mean flow and gravity are pointing in positive $x$- and negative $y$-directions, respectively. The mean flow profile ${{\langle {\boldsymbol {u}_f} \rangle _{xzt}}}=({{\langle {u_f} \rangle _{xzt}}}(y),0,0)^\textrm {T}$ is shown in blue, while the green curve represents the wall-normal variation of the mean fluid shear stress $\tau _f(y)$; ${h_f}$ and ${h_b}$ are the local instantaneous fluid and bed height, respectively. Particles are coloured depending on their location: bed particles are coloured in black, interface particles in orange to yellow with increasing wall distance and transported particles are indicated by white colour.

The flow in both, single- and multi-phase systems, is driven by a time-dependent pressure gradient that maintains a constant fluid mass flow rate ${q_f}$ in the streamwise direction throughout the entire simulation interval. Gravity is acting in the negative $y$-direction with amplitude $g={|\boldsymbol {g}|}$. The simulation domain is periodically repeated in the wall-parallel $x$- and $z$-directions with fundamental periods ${L_x}$ and ${L_z}$, respectively, which allows an investigation of the self-forming sediment ridges in the absence of a sidewall-induced mean secondary flow. The lower bottom wall and the flat free surface are accounted for as no-slip and free-slip boundary conditions, respectively. In the wall-normal direction, the computational box has an extent ${L_y}$. Note that in the sediment-laden flows, the domain consists of a particle-dominated subdomain of mean height ${H_b}$, henceforth denoted as ‘the sediment bed’, and the upper fluid-dominated region of mean height ${H_f}={L_y}-{H_b}$. A rigorous definition of the two distinct sub-domains and their mean heights will be given in § 4. In the sediment-laden simulations, the fluid–bed interface takes the role of a virtual wall, and hence we will frequently refer to a shifted wall-normal coordinate ${\tilde {y}} = y - {H_b}$.

The mean wall-shear stress ${\tau _w}$ is evaluated at the wall-normal location of the virtual wall ${\tilde {y}}=0$ by interpolating the pure fluid stress

(3.3)\begin{equation} {\tau_f}({\tilde{y}})/{\rho_f}={\nu_f} \,\mathrm{d}{{\langle {\boldsymbol{u}_f} \rangle_{xzt}}}/\mathrm{d}{\tilde{y}} - \langle (u_{f}-\langle u_{f} \rangle_{xzt})\, v_{f}\rangle_{xzt} \end{equation}

from the bulk where it varies essentially linearly down to this location. Two additional characteristic length scales of the system are the particle diameter $D$ and the viscous length ${\delta _{\nu }}={\nu _f}/{u_\tau }$, respectively, where ${u_\tau }=\sqrt {{\tau _w}/{\rho _f}}$ is the friction velocity. Following the general conventions, we shall denote quantities that are non-dimensionalized with ${u_\tau }$ and/or ${\nu _f}$ with a superscript $(\bullet )^{+}$ and refer to them as scaled in inner or wall units.

From these three length scales, we derive the friction Reynolds number ${Re_\tau }={H_f}^{+}={H_f}/{\delta _{\nu }}$, the particle Reynolds number ${D^{+}}=D/{\delta _{\nu }}$ and the relative submergence ${H_f/D}$. Let us further introduce the bulk Reynolds number as ${Re_b}={q_f}/{\nu _f}={u_b}{H_f}/{\nu _f}$, where ${u_b}={q_f}/{H_f}$ is the bulk velocity. Note that in virtue of the increased parameter space that comes with the additional degrees of freedom of the mobile particles, it requires two additional non-dimensional numbers to fully describe the two-phase system. Here, we choose the density ratio of the two phases ${\rho _p}/{\rho _f}=2.5$, where the value of $2.5$ is close to that of glass beads or sand grains in water, and the Galileo number $Ga={u_{g}} D/{\nu _f}$, where ${u_{g}}=\sqrt {({\rho _p}/{\rho _f} - 1){|\boldsymbol {g}|} D)}$ is the gravitational velocity. The squared ratio of the Galileo and the particle Reynolds number ${\theta }=({D^{+}}/{Ga})^{2}=({u_\tau }/{u_{g}})^{2}$, the Shields number, can be understood as a ratio between destabilizing and stabilizing effects and is thus a measure for the ability of a flow to erode sediment. A conventionally applied critical value of ${\theta _c}= 0.03-0.05$ marks the onset of sediment erosion in a turbulent stream, that slightly depends on the Galileo number (Soulsby & Whitehouse 1997; Wong & Parker 2006; Franklin & Charru 2011).

Tables 1 and 2 summarize the relevant parameters of the simulations. As a shorthand notation, we will identify each case with a name combining the size of the underlying computational domain and its inherent friction Reynolds number in the remainder of this work. The smooth-wall single-phase flows are indicated by a respective subscript. The short (S), medium (M) and large (L) domains feature streamwise and lateral periods of approximately ${L_x}/{H_f}\in \{2,5,12\}$ and ${L_z}/{H_f}\in \{3,16\}$, respectively, and friction Reynolds numbers $250\lesssim {Re_\tau } \lesssim 850$. Let us remark that, for the smooth-wall reference simulations, the values of the friction Reynolds number are set to ${Re_\tau }=210$ and ${Re_\tau }=650$. These values match those of cases L250 and M850, respectively, but when the bed is still at rest. As the bed evolves, the values of ${Re_\tau }$ increase as a result of the sediment bed modulation. The box dimensions ${L_x}^{+}$ and ${L_z}^{+}$ are sufficiently larger than those of the Jiménez & Moin (1991) minimum flow unit and as such accommodate the full near-wall regeneration cycle. A comparison of the outer-scaled box width ${L_z}/{H_f}$ with the findings of Flores & Jiménez (2010) for closed channels, on the other hand, indicates that our narrow domains are just wide enough to host a single full regeneration cycle of the largest log-layer streaks, for which the authors have estimated a minimal box width ${L_z}/{H_f} \approx 3$. Surely, the comparison with our open channel is limited in the vicinity of the free surface, but should give a good estimate for the rest of the domain.

Table 1. Physical parameters of the simulations. Here, $Re_b$, $Re_{\tau }$ and $D^{+}$ are the bulk Reynolds number, the friction Reynolds number and the particle Reynolds number, respectively. The density ratio $\rho _p/\rho _f$ and the Galileo number $Ga$ are imposed in each simulation, whereas the Shields number $\theta$, the relative submergence $H_f/D$ and the relative sediment bed height $H_b/D$ are computed a posteriori (cf. table 2).

Table 2. Numerical parameters of the simulations. The computational domain has dimensions $L_i$ in the $i$th direction and is discretized using a uniform finite difference grid with mesh width $\Delta x=\Delta y=\Delta z$ for the multi-phase simulations, while the smooth-wall simulations were performed using a spectral method featuring a non-uniform distribution of the grid/collocation points in the three spatial directions. Here, $N_p$ is the total number of particles in the respective case and ${T_{obs}}$ is the total observation time of each simulation starting from the release of the moving particles at $t=0$. Time dependent physical and numerical parameters in tables 1 and 2 (${Re_\tau }$, ${D^{+}}$, ${H_f}$, ${H_b}$, $\theta$, ${\Delta y^{+}}$) are computed as an average over the entire simulation period.

The preparation of each individual simulation follows a careful procedure described in detail by Kidanemariam & Uhlmann (2014a), in which first a pseudo-randomly arranged sediment bed is created through settling of a large number of particles under the action of gravity in a quiescent fluid. The number of particles per simulation has been set such that all simulations feature a comparable relative submergence ${H_f/D}\approx 26-28$, while the Galileo number is adjusted a priori depending on the desired resulting Shields number. In all current simulations, the target Shields number is chosen to ensure a bedload-dominated state clearly above the onset of sediment erosion, in which there is an intense particle transport in the near-bed region while suspended sediment transport can be assumed to have negligible effects on the bedform evolution. First, a turbulent flow field is developed by fixing all particles in space. Then, the particles are released again at a time which we define as $t=0$. Let us remark that a small fraction of the particles, namely the layer closest to the bottom wall, remain fixed in space even for times $t>0$ to create a rough boundary. In the following, we will use the bulk time unit ${T_b}={H_f}/{u_b}$ as a reference time.

As discussed by Kidanemariam & Uhlmann (2017) and Scherer et al. (2020), ridges appear shortly after particles are released, but eventually, transverse bedforms appear and dominate the evolution process after $100-200 {T_b}$. For the investigation of ridge evolution at the current parameter point, consequently, there exist basically two options: either the observation time ${T_{obs}}$ in arbitrary long domains is limited to the initial simulation phase in which transverse bedform instabilities are still small enough to be of minor importance, or the streamwise domain length ${L_x}/D$ is reduced below the critical value of ${{L_x}}_{,c}/D \approx 80$ determined by Scherer et al. (2020) such that transverse bedform instabilities are effectively suppressed for arbitrary simulation times.

Again recalling that the main concern of the current study is the role of large-scale structures for the generation of sediment ridges, we prefer the former concept that allows the accommodation of longer large-scale streaks and reduces the domain size effects, accepting the limited observation interval as it is still sufficient to study all relevant processes due to the fast formation of the sediment ridges within a few bulk time units. In the single-phase simulations M650_smooth and L250_smooth, on the other hand, we have exploited the opportunity to gather statistics over clearly longer time intervals than in their multi-phase counterparts. For comparisons with longer time series over ridges, we have additionally performed simulation S250 that features a streamwise box length ${L_x}=51.2D$ to suppress the rise of transverse bed features. This concept is very similar to the streamwise-minimal channels of Toh & Itano (2005) in that sense that the large-scale streaks can be considered as infinitely long due to the missing spatial de-correlation, while the self-sustained regeneration cycle in the buffer layer (Hamilton et al. 1995) is captured by the domain without any restrictions.

Eventually, it should be stressed that case M250 is found at the same parameter point as case $H6$ of Kidanemariam & Uhlmann (2017), but represents an own independent simulation conducted in the context of this study. Due to technical reasons, data in the first roughly $5{T_b}$ of the simulation are not available for the study. Furthermore, to the best of the authors’ knowledge, the newly presented simulations L250 and M850 represent the largest number of fully resolved particles simulated and the highest ever obtained Reynolds number in direct numerical simulation (DNS) studies with fully resolved particles, respectively. The resolution which is required to correctly resolve all flow scales at comparably high Reynolds numbers naturally comes with immense computational costs, summing up to a total amount of approximately 9 million CPU hours consumed for the investigated simulation interval of around 60 bulk time units in case M850, including a number of around 16.8 billion grid nodes.

4. Definition of the fluid–bed interface

The determination of the two-dimensional interface that separates the domain in two distinct regions, namely the sediment bed and the fluid-dominated region above it, follows the procedure described in Scherer et al. (2020). For the sake of completeness, we will briefly recapitulate the concept in the following. For a more detailed presentation of the methodology, the interested reader is referred to the original study.

First, let us decompose the vertical domain length ${L_y}$ for each point of the $(x,z)$-plane into two contributions as

(4.1)\begin{equation} {L_y} = {h_b}(x,z,t) + {h_f}(x,z,t), \end{equation}

where ${h_b}$ and ${h_f}$ are the local bed height and thickness of the fluid-dominated region, respectively. In the special coordinate system that we have chosen, the local fluid–bed interface that separates the two distinct sub-domains is found at $y={h_b}(x,z,t)$. Implicitly contained in these considerations is that the sediment bed contour can be formulated in a continuous way, while the sediment bed is naturally only represented by a set of discrete particles. Several methods are known to approximate the presence of the sediment bed in a continuous way, for instance, by identifying the interface with the wall-normal location at which a threshold for the solid volume fraction is attained, as in Kidanemariam & Uhlmann (2014a,  2017).

The methodology proposed in Scherer et al. (2020) follows another approach that first identifies a set of ‘interface particles’ which form the uppermost sediment layer of the bed by an algorithm that will be outlined below. In a second step, a continuous two-dimensional manifold is obtained through triangulation between the interface particles and interpolation to a regular grid. The sorting algorithm is based on the conventional morphological classification of sediment as either belonging to the quasi-stationary bed or being transported within the bedload or suspended load layer (Bagnold 1956; van Rijn 1984). Potential candidates for ‘bed particles’ are only those that feature a negligible kinetic energy $(|{\boldsymbol {u}_p}|/{u_{g}})^{2}$ and a non-vanishing wall-normal contact force $F_{c,y}/F_{w}$, which a particle necessarily feels while being part of a densely packed bed. Here, $F_{w} = ({\rho _p}-{\rho _f}){|\boldsymbol {g}|} {\rm \pi}D^{3}/6$ is the submerged weight of a single spherical particle. All particles that fulfil these requirements are classified as bed particles. From this set, we determine in a second step the uppermost sediment layer (the interface particles) geometrically using the $\alpha$-shape algorithm of Edelsbrunner & Mücke (1994). While conceptually similar to a conventional convex hull around a set of points, the $\alpha$-shape allows non-convexity for length scales over some threshold radius $\alpha$ (here taken as $1.1$ times the particle diameter) while it is strictly convex for length scales smaller than this threshold. Under these restrictions, an enclosing surface can be generated by means of triangulation. Nodes that relate to interface particle centres are those that bound a triangle with an outward pointing normal with positive wall-normal component. The information about the local bed height is eventually transferred to a regular equidistant Eulerian grid in the $(x,z)$-plane with sampling width of $1D$ by means of linear interpolation. An exemplary visualization of bed, interface and mobile particles together with the generated fluid–bed interface is supplemented to figure 2.

5. Results

Classically, sediment ridges have been described as long, streamwise-aligned patterns (Casey 1935; Vanoni 1946), that are ‘parallel to each other, of little relief, and of a uniform transverse spacing’ (Allen 1968). The instantaneous snapshots of the sediment bed provided in figure 3 clearly reveal that the ridges observed in the current simulations share all these features. The sediment bed is covered by laterally alternating ridges and troughs of small amplitude ($1D-2D$) that are essentially parallel to each other and to the mean flow direction. The bedforms are seen to span over the entire streamwise domain for all cases, which is in complete agreement with experimental observations that ridges can easily reach streamwise extensions of ${O}(10{H_f})$ (Wolman & Brush 1961). Perhaps even more important is the fact that these observations confirm for the first time that sediment ridges can form due to mechanisms that are completely independent of sidewall-induced secondary currents (Ikeda 1981; Colombini 1993).

Figure 3. Instantaneous visualization of the evolved sediment ridges (a,c) compared with the instantaneous distribution of three-dimensional Reynolds stress-carrying $Q^{-}$ structures (b,d) characterized by connected regions fulfilling $|-{u_f^{\prime }}({\boldsymbol {x}},t){v_f^{\prime }}({\boldsymbol {x}},t)| > H u_{rms}(y)v_{rms}(y)$ with $H=1.75$ (Lozano-Durán et al. 2012), where $u_{i,rms}=\sqrt{\langle (u_{fi}-\langle u_{fi} \rangle_{xzt})^2\rangle_{xzt}}$ is the root-mean square of the i-th velocity component ($i=x,y,z$). Ejection structures are indicated by reddish colours, while sweeps are coloured in blue, with brighter colours indicating a larger distance to the bottom wall. Particles are coloured depending on their wall-normal location, ranging from dark to light brown with increasing coordinate $y$. For the sake of clarity, only bed and interface particles are shown (cf. § 4). In each panel, flow is from bottom left to top right; (a,b) M250 ($t/{T_b}=40$), (c,d) M850 ($t/{T_b}=59$).

To give a first impression on the comparable lateral organization of bed and large-scale structures, we have supplemented to figure 3 instantaneous visualizations of the Reynolds stress-carrying ejection (${u_f^{\prime }}<0$, ${v_f^{\prime }}>0$) and sweep structures (${u_f^{\prime }}>0$, ${v_f^{\prime }}<0$) (collectively termed as $Q^{-}$) introduced by Lozano-Durán et al. (2012) as a generalization of the classical quadrant analysis (Wallace, Eckelmann & Brodkey 1972) to three-dimensional objects. According to the former study, coherent ejection and sweep structures are connected sub-domains for which $|-{u_f^{\prime }}({\boldsymbol {x}},t){v_f^{\prime }}({\boldsymbol {x}},t)| > H u_{rms}(y) v_{rms}(y)$ with $H=1.75$ and $u_{rms}=\sqrt{\langle\bar{u}_{f}\bar{u}_{f}\rangle_{xzt}}$.

It is remarkable that regions of intense ejections seem to be reasonably well correlated with the lateral positions of preferential sediment deposition, i.e. ridges, while, vice versa, intense sweep events seem more likely to occur above the troughs where erosion is dominant. Our observations support those of Gyr & Schmid (1997) that sediment erosion is mainly due to local sweep events that are naturally directed towards the bed (Jiménez 2018). Note that, by definition, ejections live in the large-scale low-speed streaks whereas sweeps populate their high-speed counterparts, thus we can equivalently conclude that ridges are found below the large-scale low-speed streaks and troughs accordingly next to the large-scale high-speed streaks, which has been also verified by comparable visualizations (plots not shown).

5.1. Sediment ridge evolution

In the following, we provide more rigorous quantitative analysis of the relation between bedforms and the turbulent flow.

As mentioned in § 3, all simulations are initiated with an initially macroscopically flat sediment bed and ridges form solely under the action of the turbulent structures after a few bulk time units. Figure 4 shows the time evolution of the fluctuations of the streamwise-averaged sediment bed ${h_f^{\prime \prime }}(z,t)$. Note that the sediment ridges in the time period of interest are essentially parallel to and without relevant perturbations in the streamwise direction such that they can be considered as essentially statistically invariant in the mean flow direction.

Figure 4. Space–time plot of the streamwise-averaged sediment bed height fluctuations ${h_b^{\prime \prime }}(z,t)/D$. Blue and red regions refer to troughs and crests of the streamwise-averaged fluid–bed interface profiles, respectively. Cases (a) S250, (b) M250, (c) M850, (d) L250.

First, it is concluded that ridges form at different spanwise locations more or less simultaneously and hence independently from each other right after the release of the particles. It is interesting to see that in the early phase, say the first $20$ bulk time units of each simulation, the lateral spacing between adjacent crests and troughs is clearly lower than the conventionally found values of $(1-2){H_f}$ (Wolman & Brush 1961; Ikeda 1981; Colombini 1993). Advancing in time, however, a coalescence of the bedforms is observed that causes a reduction of the number of individual ridges. After approximately $40$ bulk time units, merging or splitting events between ridges are less frequent but still observable and the remaining bedforms now feature a lateral spacing approximately in the expected range and a more pronounced amplitude throughout all simulations. Note that an exact match between the observed/estimated spacings and our results is not to be expected, as the time scales at which the former are observed differ by several orders of magnitude and represent quasi-asymptotic states whereas the currently observed state is highly transient. Further recall that the majority of the available experimental datasets have been conducted under the influence of lateral sidewalls in low to moderate aspect ratio flumes and underlie the sidewall-induced secondary currents. The larger and further developed ridges appear to be rather immobile in the lateral direction over the observed time period, i.e. the mean spanwise position of these ridges is quite stable. This and the regular spacing of sediment ridges seem to be no effect of laterally narrow domain sizes, as ridge crests follow the same straight vertical space–time lines in the sufficiently wide channel of case L250.

We further conclude that the narrow boxes with a lateral domain period ${L_z}/{H_f} \approx 3$ are capable of accommodating one or two ridges, whereas the large domain of L250 exhibits nine to ten ridge units. We can therefore consider the small to medium domains as close to minimal in the context of the number of available ridges, which shall be favourable for the subsequent analysis in a sense that individual ridges and their relation to turbulent coherent structures can be investigated excluding possible merging or splitting effects between individual bedforms. The large domain simulation L250, on the other hand, contains a sufficient number of individual ridges to allow statements on the collective behaviour of the bedforms.

Figure 5 provides the time evolution of the mean bedform geometry and the particle transport rate. Figure 5(a,b) shows the development of the mean pattern height expressed by the root mean square of the sediment bed height fluctuations ${\sigma _{h,z}}$ and that of the lateral spacing in terms of the mean bedform wavelength ${\lambda _{h,z}}$ (cf. Appendix B.1 for definitions). In accordance with the discussed space–time plots, it is seen that the bedform amplitude globally increases with time during the initial approximately $40$ bulk time units. In the first approximately $10$ bulk time units, low-amplitude disturbances first form, before they merge in the subsequent phase between $t=10{T_b}$ and $t=20{T_b}$, leading to a net reduction of the number of ridges. The bedforms continue growing in amplitude predominantly during the time interval between $t=20{T_b}$ and $t=40{T_b}$. The lower growth rate in the first $5$ bulk time units of the three low Reynolds number cases S250, M250 and L250 is a consequence of the lower Shields numbers when compared with the high Reynolds number case M850, in which the high Shields number causes a stronger erosion and thus a faster pattern growth. For the sake of comparison, figure 5(a) additionally provides the evolution of the spanwise-averaged root mean square of the sediment bed height fluctuations ${\sigma _{h,x}}$, that is a measure for the amplitude of possibly evolving transverse ripple-like bed features. It is thereby verified that transverse bed features are indeed not relevant in this initial time frame, consistent with the instantaneous bed snapshots in figure 3. The variation of the number of individual patterns is revealed by the oscillations of the mean wavelength ${\lambda _{h,z}}$ during the first $40$ bulk time units. For $t>40{T_b}$, then, ${\lambda _{h,z}}$ settles without further strong oscillations in all cases, attaining final values of $1.47{H_f}$ (S250), $1.14{H_f}$ (M850) and $1.44{H_f}$ (L250) that are of comparable size to the dominant wavelength ${\lambda _{h,z}}=1.3{H_f}$ obtained by means of linear stability analysis (Colombini 1993). The clearly lower attained final value of $0.76{H_f}$ in case M250 is in agreement with the corresponding space–time plot in figure 4 which shows that the bed in the final phase indeed features three ridge patterns.

Figure 5. Time evolution of fluid–bed interface dynamics. (a) Root mean square of the fluctuation of the streamwise-averaged $\sigma _{h,z}/D$ (solid lines) and spanwise-averaged fluid–bed interface $\sigma _{h,x}/D$ (dashed lines). (b) Mean lateral ridge wavelength $\lambda _{h,z}/{H_f}$. The black solid line marks the most amplified wavelength determined in the linear stability analysis of Colombini (1993), while horizontal dashed lines of matching colour mark the lateral domain size ${L_z}/{H_f}$. (c) Streamwise particle flux ${{\langle {{q_{p,x}}} \rangle _{xz}}}/q_{ref}$ as a function of time. The reference particle flux $q_{ref}$ is computed based on the Wong & Parker (2006) version of the classical formula of Meyer-Peter & Müller (1948), that is, $q_{ref}/({u_{g}}\, D)=4.93({\theta }(t)-{\theta _c})^{1.6}$, with the critical Shields number ${\theta _c}=0.034$ (Soulsby & Whitehouse 1997). (d) Mean particle flux density ${\langle {\phi {u_p}} \rangle _{xzt}}/(q_{ref}/D)$ as a function of the bed/wall distance. The inset shows the same quantity in the near-bed region, with the wall distance scaled with the particle diameter $D$. S250 (blue), M250 (green), M850 (yellow), L250 (red).

Figure 5(c) contains the temporal evolution of the streamwise particle flux ${{\langle {{q_{p,x}}} \rangle _{xz}}}(t)$ (cf. (B4a) in Appendix B) showing that the particle flow rate asymptotically approaches a quasi-steady state that is well estimated based on the empirical formula of Meyer-Peter & Müller (1948) in the modified form of Wong & Parker (2006), viz.

(5.1)\begin{equation} q_{ref}/({u_{g}} D)=4.93({\theta}(t)-{\theta_c})^{1.6}, \end{equation}

where we have chosen a critical Shields number ${\theta _c}=0.034$ (Soulsby & Whitehouse 1997). This further supports earlier findings in the context of transverse patterns by Kidanemariam & Uhlmann (2017) that relation (5.1) reasonably well describes the mean particle flow in the presence of sediment bedforms, even though it was actually developed for flow over a flat sediment bed. It should be stressed that the apparent ’overshoot’ of ${{\langle {{q_{p,x}}} \rangle _{xz}}}(t)/q_{ref}$ during the initial 20 bulk time units is not a feature of the particle flux evolution ${{\langle {{q_{p,x}}} \rangle _{xz}}}(t)$ which increases in an essentially monotonic way in this time window, before it settles to a quasi-stationary plateau (not shown). Instead, it is observed that relation (5.1) underpredicts the actual particle flux in this initial transient regime owing to the different growth rates of the particle flux and the non-dimensional wall-shear stress, while it provides a good approximation in the quasi-steady interval in which both quantities vary only weakly.

All discussed simulations belong to the bedload-dominated regime, that is, particle transport focusses on a layer of thickness ${O}(D)$ above the bed in which particles are transported by rolling, sliding and small jumps (saltation) without losing the contact with the bed for longer time intervals (van Rijn 1984). This is verified in figure 5(d) which shows the wall-normal profile of the mean particle flux density ${\langle {\phi {u_p}} \rangle _{xzt}}(y)$ (cf. (B6)), highlighting that most of the transported particles indeed remain in the near bed region. It becomes evident from the wall-normal expansion of the flux density that the higher Shields number ${\theta }$ in case M850 leads to a thickening of the bedload layer by around $1D$ compared with the remaining cases. We conclude that in this latter case, a layer of $0.25{H_f}$ thickness above the interface is characterized by a non-negligible particle flux.

5.2. Turbulent mean flow

The presence of the mobile sediment severely modifies the mean fluid flow characteristics. Figure 6(a) shows the modification of the mean fluid velocity profile ${{\langle {\boldsymbol {u}_f} \rangle _{xzt}}}^{+}(y)$ compared with that of the single-phase cases. It is seen that the mean shear $\mathrm {d}{{\langle {\boldsymbol {u}_f} \rangle _{xzt}}}/\mathrm {d} y$ reduces to negligible value when approaching the mean fluid–bed interface at ${\tilde {y}}=0$, whereas in the smooth-wall case the mean shear peaks at the wall and is responsible for the entire wall-shear stress. In the multi-phase cases, on the other hand, the shear at the location of the virtual wall is clearly dominated by the contributions stemming from particle–fluid interactions (Kidanemariam & Uhlmann 2017). The mean velocity ${{\langle {\boldsymbol {u}_f} \rangle _{xzt}}}$ increases significantly slower with increasing distance to the bed in all particulate cases, leading to a velocity deficit compared with the smooth-wall cases, an essential feature arising in turbulent flows over rough walls due to an enhanced friction (Flores & Jiménez 2006; Chan-Braun, García-Villalba & Uhlmann 2011). The mean shear consequently increases further away from the bed compared with the smooth-wall cases in order to maintain the constant bulk flow rate ${q_f}$, resulting in a different slope of the velocity profile in the log and outer layers. For the low Reynolds number cases, the deviation is rather weak and most of the effect of the sediment bed can be compensated by the conventional shift of the log-layer velocity profile ${{\langle {u_f} \rangle _{xzt}}}^{+}=1/0.41 \ln ({\tilde {y}^{+}}) + 5.2 - {\Delta U^{+}}$ (Jiménez 2004) by an offset ${\Delta U^{+}}=3.59$ (M250) and ${\Delta U^{+}}=3.73$ (L250), respectively. For the high Reynolds number case M850, on the other hand, the deviation of the slope is more pronounced such that a simple vertical shift will not lead to a satisfactory match of the profiles.

Figure 6. (a) Wall-normal profiles of the mean velocity ${{\langle {u_f} \rangle _{xzt}}}^{+}$ as a function of ${\tilde {y}^{+}}$ (semi-logarithmic scaling). The dashed-dotted line in (a) shows the standard predictions for the linear viscous sublayer and logarithmic layer, respectively, ${{\langle {u_f} \rangle _{xzt}}}^{+}={\tilde {y}^{+}}$ and ${{\langle {u_f} \rangle _{xzt}}}^{+}=1/\kappa \ln ({\tilde {y}^{+}})+B$ ($\kappa =0.41$, $B=5.2$). (b) Wall-normal profiles of the mean turbulence intensities $u_{rms}$ (solid lines), $v_{rms}$ (dashed lines), $w_{rms}$ (dotted lines). Here, M250 (green), M850 (yellow), L250 (red), M650_smooth (grey), L250_smooth (black).

The modifications of the mean velocity profile have direct implications for the intensity and distribution of the Reynolds stresses, which are depicted in figure 6(b) exemplary for cases M250, M850 and M650_smooth, respectively. Recall that the presented statistics in the multi-phase cases cannot be taken as fully converged as the observation interval is quite short and in addition highly transient. While we can assume to have captured enough statistics for the short-living small scales in the near-wall region, deviations between single- and multi-phase simulations close to the free surface are expected to be a consequence of the restricted averaging interval ${T_{obs}}$ with respect to the lifetime of the largest scales. Figure 6(b) indicates that in the low Reynolds number particle-laden cases, the typical buffer-layer peak of the streamwise normal stress ${{\langle {{\overline {u}}_f{\overline {u}}_f} \rangle _{xzt}}}$ is visibly reduced but still detectable at a comparable wall-normal position (data for L250 not shown), whereas it is completely absent in the high Reynolds number case M850 as a consequence of the intense particle motion in the latter case which destroys the near wall regeneration cycle in a similar way as in the case of flows past fully rough walls (Jiménez 2004; Flores & Jiménez 2006; Mazzuoli & Uhlmann 2017). In the low Reynolds number cases, on the other hand, particles are smaller by a factor of three in terms of wall units ($D\approx 10{\delta _{\nu }}$), which has been observed to be sufficiently small for particles to be transported by the buffer-layer structures rather than to destroy the near-wall cycle (Kidanemariam et al. 2013). Deviations for the remaining normal stresses ${{\langle {{\overline {v}}_f{\overline {v}}_f} \rangle _{xzt}}}$ and ${{\langle {{\overline {w}}_f{\overline {w}}_f} \rangle _{xzt}}}$ are also observable, but are clearly weaker and visually restricted to a layer of thickness ${\tilde {y}}\approx 0.11{H_f}$ (${\tilde {y}}\approx 3D$) above the bed in the low Reynolds number cases and ${\tilde {y}}\approx 0.17{H_f}$ (${\tilde {y}}\approx 4.5D$) in the high Reynolds number case, respectively.

In the following, we scrutinize the impact of the mobile sediment bed on the distribution of the kinetic energy among the different scales away from the wall. To this end, the premultiplied streamwise energy spectra of cases M250 and M850 are analysed in figure 7 in the centre of the clear-fluid region (${\tilde {y}/{H_f}}=0.5$) for smooth-wall and particle-laden simulations. The instantaneous streamwise energy spectra at a given wall-normal location $y$ is introduced as

(5.2)\begin{equation} {\phi_{uu}}({k_x},y,{k_z},t)={\hat{u}_f}{\hat{u}_f^{*}}, \end{equation}

where ${\hat {u}_f}({k_x},y,{k_z},t)=\mathcal {F}({u_f^{\prime }}({\boldsymbol {x}},t))$ is the Fourier transform of the fluctuating field in the two wall-parallel directions, with the asterisk indicating complex conjugation. ${k_x}=2{\rm \pi} /{\lambda _x}$ and ${k_z}=2{\rm \pi} /{\lambda _z}$ are the streamwise and spanwise wavenumber–wavelength pairs. Velocity spectra for the remaining velocity components are defined accordingly.

Figure 7. Time-averaged premultiplied streamwise energy spectra ${k_x} {k_z} {{\langle {\phi _{uu}} \rangle _{t}}}({\lambda _x},y,{\lambda _z})\,{H_f}^{2}/{u_\tau}^{\!2}$ at the reference wall-normal location ${\tilde {y}/{H_f}} = 0.5$ for cases (a) M250, (b) M850. Coloured isolines are 0.2(0.2)0.6 times the maximum value of the respective energy spectrum, while grey-shaded areas represent the same quantity evaluated for the smooth-wall reference simulations (a) L250_smooth and (b) M650_smooth, respectively. The streamwise and spanwise domain periods ${L_x}$ and ${L_z}$ of the particle-laden simulations are highlighted by dashed lines in the respective colours.

Apart from some expected fluctuations due to the limited statistics, one may conclude a good agreement between the spectral patterns in both flow configurations. As expected, the medium domains cover only part of the full spectra, but this part is again observed to reasonably well collapse with that of the smooth-wall reference simulations. The good match between spectra in single- and multi-phase simulations is an indication that the evolving sediment bed does not significantly alter the streamwise energy spectra away from the bed and that the large-scale streaks over ridges are in fact essentially the same structures as those observed over smooth walls. While it is established for flow over fixed roughness elements that the roughness effect remains restricted to a layer of several multiples of their height (Jiménez 2004,  2018; Mazzuoli & Uhlmann 2017) and does not strongly affect the organization of the large scales outside the roughness, a similar statement was not evident for the current case.

The wall-normal variation of the energy spectra is clearly identifiable in figures 8 and 9, where the premultiplied streamwise-integrated energy spectra are presented as functions of the lateral wavelength ${\lambda _z}$ and the wall-normal distance ${\tilde {y}}$ in inner and outer scales, respectively. Note that we have scaled each spectrum by its maximum value over all wall-normal locations to highlight how the energy distribution over the wall-normal direction is modified through the mobile sediment. The presented energy spectra are consistent with our claim that the relative particle size ${D^{+}}$ together with the high particle transport rate in the bedload layer of case M850 cause the immediate breakdown of the buffer-layer regeneration cycle. We observe that the near-wall peak in the streamwise spectra at ${\lambda _z}^{+}\approx 100$ associated with the buffer-layer streaks is still present for the low Reynolds number case M250 and kinetic energy is concentrated around this region whereas it is absent in case M850, in which the kinetic energy peaks further away from the wall outside the buffer layer. As indicated by the dashed line, most of the energy agglomerates in this latter case in regions above the bedload layer with less intense particle transport. A direct comparison with the smooth-wall reference case shows that this region is connected to the already existing second energetic peak at ${\lambda _z}/{H_f}\approx 1-1.5$ (${\lambda _z}^{+}\approx 1000$) that reflects the large-scale streaks at higher Reynolds numbers.

Figure 8. Time-averaged and streamwise-integrated premultiplied energy spectra ${k_z} \int {{\langle {\phi _{uu}} \rangle _{t}}}({\lambda _x},y,{\lambda _z}) \,\mathrm {d}{k_x} {H_f}^{2}/{u_\tau }^{2}$ as a function of the inner-scaled spanwise wavelength ${\lambda _z}^{+}$ and wall distance ${\tilde {y}^{+}}$, respectively. Cases (a) M250, (b) M850. Coloured isolines are 0.2(0.2)0.6 times the maximum value of the respective energy spectra, while grey-shaded areas indicate the same quantities determined for the smooth-wall reference simulations (a) L250_smooth and (b) M650_smooth, respectively. The mean fluid height ${H_f}^{+}={Re_\tau }$ and the spanwise domain period ${L_z}^{+}$ of the particle simulations are marked by coloured dashed lines. The dashed black line refers to the wall-normal distance at which the mean particle flux density ${\langle {\phi {u_p}} \rangle _{xzt}}$ attains its maximum (cf. figure 5d).

Figure 9. Time-averaged and streamwise-integrated premultiplied energy spectra for cases M850 and M650_smooth as in figure 8, but with the lateral wavelength ${\lambda _z}/{H_f}$ and the wall-normal distance ${\tilde {y}/{H_f}}$ scaled in outer units. An additional reference height is added in form of the crest height of the mean fluid–bed interface averaged over $t/{T_b}\,\in \,[40,85]$ (black solid line); (a) streamwise ${{\langle {\phi _{uu}} \rangle _{t}}}$ and (b) wall-normal spectra ${{\langle {\phi _{vv}} \rangle _{t}}}$.

Turning our attention to the wall-normal energy spectra depicted in figure 9(b), we first conclude that the medium domains considered here are just wide enough to accommodate the major fraction of the wall-normal energy spectra and hence are at the limit of being minimal near the free surface (Jiménez 2013a), which further supports our observations in § 3. In the low Reynolds number case M250, the spectra of single-phase and particle laden simulations again almost collapse, highlighting the fact that the distribution of energy over the scales and the wall-normal layers is essentially the same as over the smooth wall (plot not shown). In the high Reynolds number case M850, on the other hand, the wall-normal kinetic energy is significantly reduced throughout the bedload layer and the region of intense vertical motion is found above the latter. However, the general ellipsoidal shape and the inclination angle are maintained. The contours of spectral energy density almost collapse in the vicinity of the free surface and for large wavelengths, whereas the intense energy region in the bulk is compressed, restricted by the bedload layer that damps the wall-normal energy even before the bed is reached, quite similar to the streamwise spectra.

5.3. Large-scale flow organization

Let us now turn to the dynamics of the large-scale flow structures in physical space. Figure 10 shows the space–time evolution of the streamwise velocity fluctuations with respect to the streamwise average, ${u_f^{\prime \prime }}(z,t)$. Note that streamwise averaging filters out lateral meandering of structures while retaining the footprint of the large streamwise-elongated structures.

Figure 10. Space–time plot of the streamwise-averaged fluctuations of the streamwise velocity component ${u_f^{\prime \prime }}/{u_\tau }$ extracted at ${\tilde {y}/{H_f}}=0.5$. Blue and red regions refer to streamwise-averaged zones of low and high streamwise velocity, respectively. Cases (a) M250, (b) M650_smooth, (c) M850, (d) L250.

In fact, the streamwise-averaged high and low-speed regions are found to develop comparably straight in the space–time plane, indicating that there is only weak lateral migration of these zones. It is noteworthy that there is no qualitative difference in this point between the ridge-featuring and smooth-wall cases during the current limited time interval, showing that even laterally homogeneous smooth channels are able to maintain a substantial spanwise modulation of the mean flow profile over an intermediate time scale as a consequence of the long-living large-scale streaks (Jiménez 2013a), which is expected to disappear for sufficiently long time averaging intervals leading to statistical homogeneity. We shall discuss the role of this intermediate time scale dynamics in the creation of ridges below. Note that a comparison between the medium domain simulations and case L250 in figure 10(d) verifies that the observations are not a consequence of the limited box width ${L_z}/{H_f}\approx 3$, but that comparable regular organization of the large-scale streaks over intermediate time intervals similarly occur in domains as wide as ${L_z}/{H_f}\approx 16$. Also, it shows that the large-scale streaks are of significant length, as they are clearly detectable even in streamwise averages over long streamwise extents ${L_x}/{H_f}\approx 12$.

Intermittently, the amplitude of the high- and low-speed regions is seen to reduce as, for instance, around $t=40{T_b}$ in case M650_smooth (figure 10b) or in a period $10\lesssim t/{T_b} \lesssim 20$ in case M850. It appears likely that at least some of these ‘events’ are related to regularly occurring bursting events in which the large-scale streaks bend and eventually break, just as their viscous counterparts in the buffer layer do (Flores & Jiménez 2010). Interestingly, there are events after which large-scale streaks recover at practically the same lateral position, while in other situations as after $t=45 {T_b}$ in case M250 (figure 10a), high- and low-speed regions reorganize and partly even change the number of large-scale streaks, as it is here the case. We will investigate the implications of such ‘events’ on sediment ridges in figure 12 below.

Throughout the previous sections, it has become clear that the large-scale streaks over ridges are no other structures than those found in smooth-wall turbulence, and we have in fact observed these structures throughout the entire ridge evolution interval. For homogeneous roughness, Flores, Jiménez & Del Álamo (2007) similarly report that the self-similar vortex clusters and the associated velocity streaks are not significantly modified outside the roughness layer. We have, however, not explicitly discussed the question of cause and effect: based on the results seen so far, it seems reasonable that the spanwise heterogeneity of the flow in form of the well-organized large-scale low and high-speed streaks generates a laterally varying bed-shear stress that leads accordingly to spanwise heterogeneous erosion from the very beginning, hence implying that the large-scale streaks trigger the first ridges and not vice versa. Note that this does not rule out the possibility of a feedback of developed ridges on the mobility or meandering tendency of the large-scale streaks in later stages of their lifetime.

Large-scale streaks are entirely turbulent structures and as such not as smooth as their viscous counterparts in the buffer layer (Jiménez 2013a). Also, large-scale instantaneous quasi-streamwise rotations are almost impossible to detect by the classical gradient-based vortex detection methods such as the $\lambda _2$-criterion of Jeong & Hussain (1995) as the velocity gradients are rather weak, as opposed to the near-wall quasi-streamwise vortices. Hence, for the understanding of the large-scale processes it is often customary to study fields in which small-scale motions are filtered out, either indirectly as in the context of ensemble averaging (see, for instance, Del Álamo et al. 2006; Lozano-Durán et al. 2012) and in some sense in our own streamwise-averaging procedure, or directly by means of instantaneous flow field filtering, for instance, with a Gaussian filter as in Motoori & Goto (2019,  2021). The advantage of direct filtering is that there is no information loss for all scales larger than the filter widths. In the current study, we therefore extend our analysis by applying a Gaussian cutoff filter to the flow field. The low-pass filtered field of the $i$th velocity component is obtained by the following convolution of the field with an anisotropic Gaussian kernel (Lozano-Durán, Holzner & Jiménez 2016):

(5.3)\begin{equation} {{\widetilde{u}}_{f,i}}({\boldsymbol{x}}) = \iiint_{V} G \cdot {u_{f,i}}({\boldsymbol{x}}-{\boldsymbol{x}}^{\prime}) \exp \left( - \left( \dfrac{{\rm \pi} x^{\prime}}{{\varDelta_x}} \right)^{2} - \left( \dfrac{{\rm \pi} y^{\prime}}{{\varDelta_y}} \right)^{2} - \left( \dfrac{{\rm \pi} z^{\prime}}{{\varDelta_z}} \right)^{2} \right)\mathrm{d} x^{\prime} \,\mathrm{d} y^{\prime} \,\mathrm{d}z^{\prime}. \end{equation}

Therein, $\varDelta _i$ ($i=x,y,z$) is the filter cutoff width in the three Cartesian directions and ${\boldsymbol {x}}^{\prime }$ is the inner-convolution coordinate. The boundary conditions are treated similarly as in Lozano-Durán et al. (2016), i.e. the flow field is periodically extended in the wall-parallel directions, while it is mirrored vertically at the bottom wall, reversing the sign of the wall-normal velocity component to take care of the fundamental symmetries. The free surface can be treated in a completely analogous way as the symmetries along this plane resulting from the free-slip boundary condition do not differ from that at a smooth wall. Here, $V$ is the filter volume and the constant coefficient $G$ is normalized to ensure that the integral of the kernel over $V$ is unity. Note that the filtered flow field in a layer of thickness ${\varDelta _y}/{H_f}$ will be affected by the velocity field inside the bed. In the following, however, we will study only filtering results in regions sufficiently away from this wall-normal region or for single-phase flow simulations, in which the above effect is naturally absent. If not otherwise stated, we conventionally use filter widths $[{\varDelta _x},{\varDelta _y},{\varDelta _z}]/{H_f}=[0.6,0.2,0.3]$ that obey the typical aspect ratio of the vortex clusters as reported by Del Álamo et al. (2006). In the following, we will also use a two-dimensional analogue of (5.3) in wall-parallel $(x,z)$-planes whenever we are mainly interested in the structure's wall-parallel extensions, using filter widths $[{\varDelta _x},{\varDelta _z}]/{H_f}= [0.6,0.3]$.

In the previous paragraphs, it has been discussed that the lateral organization of large-scale high- and low-speed velocity streaks and the associated ejection and sweep events may cause a laterally varying instantaneous shear stress along the bed or the lower wall, respectively, as recently discussed in the hydraulic context by Bagherimiyab & Lemmin (2018). Underlying this hypothesis is today's understanding that the near-wall region is not exclusively populated by short scales of order ${O}(100{\delta _{\nu }})$, but that it also features large scales of order ${O}({H_f})$ that scale in outer units and that are correlated with the structures away from the wall (Del Álamo & Jiménez 2003). Even though the near-wall regeneration cycle itself is capable of functioning autonomously (Jiménez & Pinelli 1999), large-scale influences have nonetheless been observed to modulate regions even close to the wall (Jiménez, Del Álamo & Flores 2004). Toh & Itano (2005) stated that there is a frequent interplay between small near-wall structures and large structures away from it, the former necessarily agglomerating below the large-scale ejections due to continuity (Jiménez 2018).

In figure 11, the organization of buffer-layer structures with dimensions ${O}(100{\delta _{\nu }})$ into larger essentially ${H_f}$-scaled streaks and the corresponding outer large-scale streaks at ${\tilde {y}/{H_f}}=0.3$ can be observed for the smooth-wall case M650_smooth in visualizations of unfiltered and corresponding spatially filtered fields. Note that qualitatively similar results are found for the flow organization over the sediment beds just before the particle release. To avoid confusion with the mean wall shear ${\tau _w}$, let us introduce the instantaneous bottom shear stress ${\tau _b}(x,z,t)={\nu _f} \,\mathrm {d}{u_f}/\mathrm {d}{\tilde {y}}|_{{\tilde {y}}=0}$ for the following analysis. In figure 11, the high-speed buffer-layer streaks whose position can be inferred by their induced zones of locally increased shear are seen to cluster in two streamwise-elongated streak-like structures, where one spans the entire box length in a slightly meandering way, while the other is roughly half as long and laterally shifted by an offset somewhat larger than ${H_f}$. These two large-scale structures are found to lay essentially below the corresponding large-scale high-speed streaks at the upper end of the logarithmic layer, ${\tilde {y}/{H_f}}=0.3$, and are of comparable lateral extension, even though the streamwise extension differs in particular for the patchy structure with centre at $z/{H_f}=0.5$. In order to quantify the correlation of the boundary shear with the flow field at each height of the channel in M650_smooth in the streamwise-averaged framework, let us define the instantaneous two-point correlation coefficient as follows:

(5.4)\begin{equation} {\rho_{u,\tau}}(y,t)={\langle {{\tau^{\prime\prime}_b}(z,t) {u_f^{\prime\prime}}(y,z,t)} \rangle_{z}}/[ {\langle {({\tau^{\prime\prime}_b}(z,t))^{2}} \rangle_{z}} \, {\langle {({u_f^{\prime\prime}}(y,z,t))^{2}} \rangle_{z}}]^{1/2}. \end{equation}

It turns out that, in accordance with the previous findings, regions of high and low momentum over the entire logarithmic layer are highly correlated with the corresponding wall shear regions along the bottom wall, with a mean value of the correlation in the channel centre of ${\langle {{\rho _{u,\tau }}} \rangle _{t}}({\tilde {y}/{H_f}}=0.5) \approx 0.5$. Note that the restricted periodic domain length makes it hard to visually detect a streamwise phase shift between the wall-shear object and the log-layer streaks which is, however, expected to exist as a consequence of the varying mean shear and thus propagation speed between the two wall-normal locations.

Figure 11. (a,b) Instantaneous wall-shear stress fluctuations ${\tau ^{\prime }_b}(x,z,t)/({\rho _f} {u_\tau }^{2})$ for the smooth-wall simulation M650_smooth (a) without applied filtering and (b) after a two-dimensional Gaussian cutoff filter has been applied in the two homogeneous directions with filter widths ${\varDelta _x}=3y_{ref}$ and ${\varDelta _z}=1.5y_{ref}$, respectively. The reference value has been set to the height $y_{ref}=0.3{H_f}$. (c,d) Corresponding streamwise velocity fluctuation field $({u_f^{\prime }})^{+}$ extracted at ${\tilde {y}/{H_f}}=0.3$ (c) non-filtered and (d) filtered with the same filter size as in (b). The time at which all snapshots have been extracted corresponds to the second marker point in figure 14(b), $t/{T_b}=83$.

The observed correlation between log-layer streaks and the organization of the bottom shear stress is of direct relevance for the current simulations, as it provides a mechanism that couples the erosion rate along the bed with the large-scale structures. While it appears inevitable that large scales drive the particle motion in case M850 due to the absence of the near-wall cycle, the observed link between large-scale structures and the local wall shear might explain why the same ridge spacing is found for the low Reynolds number cases in which small scales are, in general, still able to transport particles (Kidanemariam et al. 2013).

5.4. Large-scale streak–ridge interaction

The above arguments imply that particle erosion and hence bed topography are linked with the large-scale structures via their near-bed effects and the local wall shear. Indeed, it appears that the number of adjacent alternating high- and low-speed regions agrees reasonably well with the spanwise wavelength ${\lambda _z}/{H_f} \approx 1.1-1.3$ of the outer energetic peak determined from the streamwise energy spectra, which, in turn, is of striking similarity to the observed mean pattern wavelength ${\lambda _{h,z}}/{H_f}$. That, indeed, regions of high and low streamwise velocity are spatially correlated with the bed contour visualized in figure 12, in an exemplary manner for case M850. Here, the space–time evolution of the sediment bed fluctuation ${h_b^{\prime \prime }}(z,t)$ and the respective fluctuation of the streamwise velocity component at ${\tilde {y}/{H_f}}=0.5$ are repeated for convenience, supplemented by similarly presented data for the streamwise and spanwise particle flux components ${q_{p,x}^{\prime \prime }}$ and ${q_{p,z}^{\prime \prime }}$, respectively. The observed dependences can be summarized as follows: regions of high momentum which are the streamwise-averaged analogue of the large-scale streaks lead to spanwise alternating zones of locally increased erosion due to turbulent sweep structures that are located in the large-scale high-speed streaks. The wall-shear stress is locally increased and with it the local Shields number. Higher shear stress and erosion directly imply an increased streamwise particle flux ${q_{p,x}^{\prime \prime }}$ in these regions. The effect is further supported by the action of the lateral particle transport ${q_{p,z}^{\prime \prime }}$, which is predominantly directed in such sense that sediment is transported from the regions of high streamwise particle transport (${q_{p,x}^{\prime \prime }}>0$) to areas with lower streamwise particle flux, i.e. ${q_{p,x}^{\prime \prime }}<0$, as indicated by the red and blue lines in figure 12(d). The consequence is a local decrease of the sediment bed height in the form of a local trough. Corresponding opposite relations are found for low-speed regions, ridges and regions of low streamwise particle transport ${q_{p,x}^{\prime \prime }}<0$.

Figure 12. Space–time evolution in case M850 of (a) the sediment bed height fluctuation of the streamwise-averaged bed (identical to figure 4c); (b) the streamwise-averaged streamwise velocity fluctuation ${u_f^{\prime \prime }}/{u_\tau }$ at ${\tilde {y}/{H_f}}=0.5$ (identical to figure 10c); (c) the streamwise-averaged particle flux fluctuation ${q_{p,x}^{\prime \prime }}/q_{x,rms}$, with $q_{x,rms}(t)=({\langle {{q_{p,x}^{\prime \prime }}{q_{p,x}^{\prime \prime }}} \rangle _{z}})^{1/2}$. (d) Same data as in (c) are shown in the background as grey map, supplemented with red (blue) dots marking regions of lateral bed growth (decrease) that are regions with vanishing lateral particle flux ${q_{p,z}^{\prime \prime }}=0$ and negative (positive) lateral gradient $\partial _z {q_{p,x}^{\prime \prime }}(z,t)$.

The lateral particle flux is found to be two orders of magnitude smaller than the corresponding streamwise particle flux ${q_{p,x}^{\prime \prime }}$ in all simulations. It is therefore conceived that the initial sediment ridges readily appearing after the onset of particle motions are mostly a consequence of the laterally varying particle erosion in the streamwise direction, resulting in the comparably stronger streamwise particle flux. Lateral particle transport from the troughs to the crests which is due to the weaker lateral fluid motion then further supports the growth of these initial sediment ridges once the particle fluxes have reached a quasi-stationary level (cf. figure 5c).

The visually identified correlation between large-scale structures and the bed contour in the multi-phase simulations can be quantified with the aid of a two-point cross-correlation as defined in (5.4), replacing the wall-shear stress ${\tau ^{\prime \prime }_b}$ by the sediment bed height fluctuation ${h_b^{\prime \prime }}$. In agreement with the previous observations, it turns out that, over most of the simulation time, the bed contour shows a high instantaneous correlation with the flow field over most of the channel height reaching values of more than 0.6 at the bulk flow centreline ${\tilde {y}/{H_f}}=0.5$ (figure omitted). This is, however, not the case in the phase between $t=10{T_b}$ and $t=30{T_b}$, in which the bed contour is found to be essentially uncorrelated with the flow structures away from the bed. The missing correlation between the large-scale structures and the location of sediment ridges is also observed in figure 12(a,b); between $t=10{T_b}$ and $t=20{T_b}$, the amplitude of the streamwise-averaged velocity fluctuations reduces and the previously and subsequently observed intense regions temporarily ‘diffuse’ into several narrower regions. We have argued earlier that the large-scale streaks break up in this phase. Two large-scale streaks recover between $t=20{T_b}$ and $t=25{T_b}$, but at somewhat different lateral positions. Interestingly, the bed contour is seen to adapt to the new spanwise organization of the large-scale streaks, laterally migrating in time towards these spanwise positions, however, with a noticeable time delay.

Figure 13. Two-point cross-time correlations between ${u_f^{\prime \prime }}(y,z,t+{\delta t})$ and: (a) the sediment bed height fluctuations, $-{\rho _{u,h_b}^{t}}(y,{\delta t})$; (b) the velocity fluctuations at the reference height ${u_f^{\prime \prime }}({\tilde {y}/{H_f}}=0.5,z,t)$, ${\rho _{uu}^{t}}(y,{\delta t})$. The reference time and reference wall-normal position are indicated by a white cross. Red (blue) regions represent strong (weak) correlation of the compared quantities. Black dashed lines connect the maximum correlation values at each wall-normal distance ${\tilde {y}/{H_f}}$. Contours separating the coloured areas indicate values of 0(0.1)1 in all panels. Data are for case M850.

This time lag between the bed evolution and the large-scale streak organization is instructive as it indicates how information is propagating across the channel. In order to determine the time delay between the dynamics of the large-scale streaks and the deformation of the sediment bed contour, we introduce the two-time cross-correlation between the streamwise-averaged flow field and an arbitrary physical quantity $a$ as

(5.5)\begin{equation} {\rho_{ua}^{t}}(y,{\delta t})= {\langle {{u_f^{\prime\prime}}(y,z,t+{\delta t}) a^{\prime\prime}(z,t)} \rangle_{z}}/ [ {\langle {({u_f^{\prime\prime}}(y,z,t+{\delta t}))^{2}} \rangle_{z}} {\langle {(a^{\prime\prime}(z,t))^{2}} \rangle_{z}}]^{1/2}. \end{equation}

Figure 13(a) shows the cross-time correlation between the flow field and the bed contour, i.e. $a^{\prime \prime }={h_b^{\prime \prime }}(z,t)$, while figure 13(b) provides the cross-time correlation between the flow field at varying wall-normal positions and that at a reference height chosen as ${\tilde {y}/{H_f}}=0.5$, i.e. $a^{\prime \prime }={u_f^{\prime \prime }}({\tilde {y}/{H_f}}=0.5,z,t)$. Note that the correlation between bed and flow field is presented premultiplied with a global factor $-1$ to take account of the fact that a locally higher velocity implies a local decrease of the sediment bed height. In figures 13(a) and 13(b), a time lag is clearly identifiable. Figure 13(a) verifies our earlier observations that ridges align themselves with the location of the large-scale streaks with a time lag of around $10{T_b}$ (equivalent to approximately $0.9{H_f}/{u_\tau }$). A very similar time delay is observed in figure 13(b) between the dynamics of the flow field at the channel centre and that close to the bed. The observed time lag implies that the generation of the initial sediment ridges is indeed predominantly dictated by the organization of the large-scale streaks in the channel centre which interact in a ‘top–down mechanism’ with the fluid layers closer to the bed.

It also indicates that the existence of dominant large-scale velocity streaks in the core of the channel is of direct relevance to the stability and position of the underlying sediment ridges and troughs, respectively. In the following, we shall analyse such events during which a clear signature of high- and low-speed regions is missing, avoiding streamwise averaging for the moment to show that the modulation of ${u_f^{\prime \prime }}$ in such phases is indeed due to a break-up of the large-scale streaks. Figure 14 shows the time evolution of the most energetic Fourier modes of the streamwise velocity fluctuations, ${\phi _{uu}}({k_x},y,{k_z},t)$, in wall-parallel planes at ${\tilde {y}/{H_f}}=0.5$. For the sake of readability, let us refer to modes with wavelength pairs ${\lambda _x}={L_x}/i$ and ${\lambda _z}={L_z}/k$ ($i,k \in \mathbb {N}_0$) as the $(i,k)$-mode in the remainder. Note further that the shown curves are normalized by the total streamwise kinetic energy contained in the fluctuating field at that time and wall-normal position, respectively, such that all contributions sum up to unity. This is particularly noticeable regarding the transient nature of the currently studied systems: the insets show the absolute value of the streamwise ${\langle {{u_f^{\prime }}{u_f^{\prime }}} \rangle _{xz}}$ (black) and wall-normal energy ${\langle {{v_f^{\prime }}{v_f^{\prime }}} \rangle _{xz}}$ (red) contained in the same plane, indicating that in all sediment-laden cases the turbulent kinetic energy globally increases with time (at least initially) as a consequence of the increasing friction along the bed due to particle transport and bedform development.

Figure 14. Time evolution of the energy distribution between individual modes of the streamwise velocity spectrum ${\phi _{uu}}({k_x},{k_z},t)$ evaluated at a wall-normal position ${\tilde {y}/{H_f}}\approx 0.5$. Note that the spectrum at each time has been normalized such that the contributions of all individual modes sum up to unity. Note that the most dominant modes (those modes that carry more than $20\,\%$ of the total energy at least once during the observation time) are highlighted using the following colour scheme: blue line $(0,1)$-mode, red line $(0,2)$-mode, purple line $(1,2)$-mode. The smaller insets show the time evolution of the total streamwise and wall-normal fluctuation energy at the same wall-normal position, (black solid line) ${\langle {{u_f^{\prime }}{u_f^{\prime }}} \rangle _{xz}}$ and (red solid line) ${\langle {{v_f^{\prime }}{v_f^{\prime }}} \rangle _{xz}}$, normalized by their time-averaged mean. Cases (a) M250, (b) M650_smooth, (c) M850, (d) L250.

In the three medium-sized boxes, the range of wavelengths is limited by the box size and only a few large-scale modes can be accommodated, such that most of the time a large fraction of the total kinetic energy is contained in these few dominant modes that are highlighted in different colours. Here, we have classified those modes as dominant that contain more than $20\,\%$ of the total perturbation energy at least once during the observation interval. From the four dominant highly energetic mode pairs that are of relevance in the medium boxes, three feature the streamwise zero mode, i.e. the corresponding structures are of infinite streamwise length and do not appear in the premultiplied energy spectra.

These modes are in fact the alternating large-scale high- and low-speed streaks identified earlier, spanning the entire streamwise domain and are thus clearly detectable in the streamwise average. In the large box L250, on the other hand, there are no specially exposed modes, but energy is rather uniformly distributed among all modes. Comparing the dynamics of the individual modes in the medium boxes with the space–time plots in figure 10, the times during which the clear signature of high- and low-speed regions disappears correspond to instances in which there is temporarily no clear dominant mode in the spectra, in particular the energy in the infinitely long modes is markedly reduced. In case M250 (figure 14a), the earlier discussed change from two to only one pair of large-scale high- and low-speed streaks is clearly visible around $t\approx 50{T_b}$, accompanied by a reduction in the total energies ${\langle {{u_f^{\prime }}{u_f^{\prime }}} \rangle _{xz}}$ and ${\langle {{v_f^{\prime }}{v_f^{\prime }}} \rangle _{xz}}$, suggesting that at this time, the large-scale streaks indeed break up. Figure 14(b) shows a comparable behaviour for the smooth-wall case M650_smooth with the difference that after the breakdown of the infinitely long two-streak mode, the same mode recovers instead of being replaced by a single pair of large-scale streaks as in the previous case. The two states without and with a clear dominant mode indicated by the two symbols in figure 14(b) are visualized in physical space in terms of fluctuating streamwise velocity fields in figure 15(a,b), which also demonstrates that there is no strong lateral meandering. In the high Reynolds number sediment-laden case M850 provided in figure 14(c), however, the streamwise zero modes carry only slightly more energy than the remaining modes at the beginning of the simulation, but shortly after the particle release (between $t=10{T_b}$ and $t=20{T_b}$) even this slight plus in energy is lost and they cannot be differentiated from the remaining modes anymore. Again, this is in line with the space–time plots and the flow field visualization in figure 15(c) underlines the absence of a clear dominating large structure. Later, however, the mode that represents two pairs of elongated large-scale streaks increases first mildly, then more strongly until it settles at $t\approx 40{T_b}$ and dominates the spectra for the rest of the simulation. The dominance of this mode is also seen in figure 15(d).

Figure 15. Wall-parallel plane of the streamwise velocity fluctuation field $({u_f^{\prime }})^{+}$ extracted at ${\tilde {y}/{H_f}}\approx 0.5$ for cases (a,b) M650_smooth (c,d) M850 at two different times. The flow fields refer to the instances marked by filled circles in figures 14(b) and 14(c), respectively.

5.5. Large-scale streaks and mean secondary flow

Up to this point, we have seen that quasi-streamwise large-scale streaks characterize the appearance of the flow field in particular in the medium boxes where a considerable fraction of the energy resides in the infinitely long streamwise modes. Some authors have hypothesized that sediment ridges may ‘lock’ the spanwise position of these large-scale streaks such that they leave a footprint in the streamwise time average as up- and downwelling of the mean velocity profile together with mean secondary currents (Nezu 2005).

Having shown that, at least over intermediate time intervals, the large-scale streaks indeed propagate little in the lateral direction and observing that instantaneous spatial meandering appears to be rather weak, the assumption that mean secondary currents are the statistical footprints of large-scale streaks appears conclusive. Indeed, depth-spanning vortices are seen to populate the cross-section of all four periodic channels in figure 16, which provides mean flow field data averaged over short time intervals (between $20$ and $80$ bulk time units) and over the streamwise direction. The secondary mean flow is expressed in terms of the streamfunction ${\langle {\psi } \rangle _{xt}}(y,z)$, which is defined by the requirement that $\nabla _{2D} {\langle {\psi } \rangle _{xt}}(y,z)= (\partial _y {\langle {\psi } \rangle _{xt}},\partial _z {\langle {\psi } \rangle _{xt}})^{\textrm {T}} = (-{{\langle {w_f} \rangle _{xt}}},{{\langle {v_f} \rangle _{xt}}})^{\textrm {T}}$. For the particle-laden simulations, the conventionally discussed configuration of primary flow upwelling over ridge crests and downwelling over the troughs is observed with accordingly oriented secondary currents (Nezu & Nakagawa 1993), representing the short-time averaged equivalent of the LMPs/HMPs studied by, for instance, Mejia-Alvarez & Christensen (2013). A comparison between the cases over mobile sediment beds and the smooth-wall simulation M650_smooth in figure 16(b) shows that, over adequately chosen intermediate time intervals of the order of the large-scale streaks’ lifetime, even statistically spanwise homogeneous simulations do generate a mean secondary flow pattern that is similar to the secondary currents next to developed ridges in figure 16(a,c,d).

Figure 16. Mean primary and secondary flow patterns for cases (a) M250, (b) M650_smooth, (c) M850, (d) L250. Note that for the sake of comparison, in panel (d) we have arbitrarily chosen a sub-domain of the entire cross-section of case L250 with a lateral width similar to ${L_z}/{H_f}$ in the medium domain cases. Isolines of the primary flow ${{\langle {u_f} \rangle _{xt}}}/{u_b}$ are shown in intervals 0(0.1)1.2, while the secondary flow pattern $({{\langle {v_f} \rangle _{xt}}},{{\langle {w_f} \rangle _{xt}}})^{\textrm {T}}/{u_b}$ is shown in terms of the secondary mean flow streamfunction ${\langle {\psi } \rangle _{xt}}$. Clockwise (counter-clockwise) secondary flow rotation is indicated by red (blue) contours. The time-averaged fluid–bed interface profile is indicated by a red curve. The time window over which data have been accumulated can be seen in figure 17. For the shown contours, the range $[\min _{y,z}{\langle {\psi } \rangle _{xt}},\max _{y,z}{\langle {\psi } \rangle _{xt}}]$ is divided into 20 equally spaced subintervals.

Figure 17. Instantaneous secondary flow intensity ${u_{\perp }}/{u_b}$ as a function of time (solid lines) and the corresponding value determined for the time-averaged fields in figure 16 (dashed lines) following the definition in (5.6). The length of the dashed lines indicates the time window over which the fields shown in the respective panels of figure 16 have been averaged: M250 (green), M650_smooth (grey), M850 (yellow), L250 (red).

In particular, the secondary flow cells in all cases share roughly the same spanwise spacing, which is equivalent to the preferential spanwise wavelength of the instantaneous large-scale streaks, ${\lambda _z}/{H_f}\approx 1-1.5$. An exception is case M250 for which in the observed time interval intermittently only one pair of large-scale low- and high-speed streaks exists, as has been seen earlier. Local mean up- and downflows are the statistical footprints of ejections and sweeps, occurring in the respective large-scale low- and high-speed streaks.

Also, we do not observe any significant quantitative difference in the amplitude of the secondary currents from case to case. Figure 17 shows the time-dependent secondary flow intensity ${u_{\perp }}(t)$ defined from the cross-plane-averaged kinetic energy contained in the cross-flow field $({\langle {{v_f}} \rangle _{x}},{\langle {{w_f}} \rangle _{x}})^{\textrm {T}}$, viz.

(5.6)\begin{equation} {u_{{\perp}}}(t) = \left[ \dfrac{1}{{L_z} {\langle {h_f} \rangle_{xz}}(t)} \displaystyle \int_{0}^{{L_z}} \displaystyle \int_{{\langle {h_b} \rangle_{xz}}(t)}^{{L_y}} ( {\langle {{v_f}} \rangle_{x}}^{2} + {\langle {{w_f}} \rangle_{x}}^{2} )\, \mathrm{d} y \,\mathrm{d}z\right]^{1/2}. \end{equation}

Note that a comparable formulation has been originally presented by Sakai (2016) for quantification of sidewall-induced mean secondary flow in open ducts and has been adapted here to the multi-phase configuration. It is seen that ${u_{\perp }}(t)$ oscillates between $1.5\,\%$ and $2.5\,\%$ of the bulk velocity for all cases, with slightly lower values in case L250 as a consequence of the larger averaging ensemble in a sense that the wider domain allows accommodation of a larger number of secondary flows cells (by a factor of three to four). Note that the secondary flow intensity of the time-averaged cross-flow $({\langle {{v_f}} \rangle _{xt}},{\langle {{w_f}} \rangle _{xt}})^{\textrm {T}}$ in figure 16 is, as expected, systematically lower than their instantaneous counterparts. Interestingly, the latter values are of comparable size to the usually reported strength of sidewall-induced mean secondary flows in open ducts. For instance, Sakai (2016) determined the intensity of secondary flow within a distance $z=1{H_f}$ from the sidewalls of a smooth-wall open duct in a range $1.1\,\%$ to $1.3\,\%$ of the bulk velocity.

Note that the observed secondary flow pattern disappears for conventional time and streamwise averaging over sufficiently long time intervals in the smooth-wall case, and it is only recovered in the context of conditional averaging over individual log- or outer-layer streaks (Lozano-Durán et al. 2012; Kevin et al. 2019a). Long-time prognosis of the situation over ridge-covered beds is only possible for the shortest simulation S250 (cf. figure 18 in Appendix A) where the observation time is not restricted through the development of transverse bedforms. Here, observations over longer time intervals of roughly $600{T_b}$ suggest that the observed ridges undergo continuous variations in form of merging, splitting or short lateral propagations, probably reacting to changes in the structure of the large-scale streak organization as seen before in the initial phase. The global trend, however, shows that the lateral positions at which ridge crests are found do not vary strongly in this narrow computational domain. In figure 18 it can be seen that the dominant ridge recovers at essentially the same lateral position even after an intermediate disappearance during the period between $t=400{T_b}$ and $t=500{T_b}$.

6. Discussion

The previous analysis has revealed that large-scale velocity streaks in turbulent open channel flow are able to generate the characteristic pattern of spanwise alternating sediment troughs and ridges on an initially flat sediment bed. The observed generation mechanism can be summarized as follows: in the initial phase during which the bed is macroscopically flat and particles are still immobile, the streamwise velocity field is organized in streamwise-elongated streaks of laterally alternating high and low momentum just as in canonical smooth-wall channel flows. The largest representatives of this self-similar streak family are of dimensions ${O}({H_f})$. At the onset of particle motion, particles are predominantly eroded in regions below the large-scale high-speed streaks, where high-momentum fluid is brought towards the bed in form of sweep structures that have been identified as the main driver of sediment erosion (Gyr & Schmid 1997). The laterally non-homogeneous erosion of sediment, in turn, gives rise to a wave-like deformation of the initially flat sediment bed, with ridges found below large-scale low-speed streaks and troughs accordingly associated with regions of relatively higher streamwise velocity. This goes hand in hand with a lateral variation of the streamwise sediment flux. The described formation mechanism implies that, during the initial formation phase investigated in the current work, bedform evolution is predominantly controlled by the large-scale velocity structures in the channel core and their associated Reynolds stress-carrying structures. This hypothesis has been supported by the two-time cross-correlations in figure 13, which reveal a time lag between the lateral organization of large-scale streaks at ${\tilde {y}/{H_f}}=0.5$ and that of the flow structures in the vicinity of the bed on the one hand and that of the sediment ridges on the other hand. For the current case M850, the time lag has been estimated as approximately $10{T_b}$ or, equivalently in terms of the eddy turnover time, $0.9{H_f}/{u_\tau }$.

The observed mechanism conceptually differs from the instability studied in the classical linear stability analysis of Colombini (1993) in which cause and effect are reversed compared with the above discussed process: a strictly one-dimensional turbulent base profile is disturbed by an infinitesimal sinusoidal perturbation of the sediment bed across the spanwise direction that triggers the evolution of secondary flow cells. The fluid motion is captured by the Reynolds-averaged Navier–Stokes equations that have been closed with a nonlinear eddy-viscosity model (Speziale 1987) to allow for the necessary anisotropy of the Reynolds stress tensor. The fluid-related equations of the model are quasi-stationary and one-way coupled to those capturing the sediment bed dynamics such that the computation can be performed sequentially: all quantities feature a sinusoidal lateral perturbation that is ‘geometrically induced’ by the disturbances of the bed. The most amplified wavelength of the system is then directly obtained from the perturbed and linearized equations of fluid motion only and, as such, does not depend on the properties of the sediment bed and the particle flux. The sediment bed continuity equation (Exner 1925) is afterwards considered to determine the growth rate of the initial perturbations a posteriori, balancing the destabilizing lateral bed-shear stress exerted on the sediment bed by the perturbed flow field and a constant stabilizing gravitational term that is a function of the Shields number only. The latter is thus in particular independent of the wavenumber and, consequently, its effect reduces to a modification of the growth rate amplitude while the chosen bedform wavelength is entirely controlled by the bed-shear stress which features the most amplified wavelength determined from the linearized fluid equations.

We conclude that in the linear stability analysis, the initially deformed lower domain boundary is indispensable to trigger the spanwise disturbance of the flow field and, that way, to cause the secondary flow instability in the context of the time-independent RANS equations. By contrast, considering the full time-dependent Navier–Stokes equations as in the current study, lateral finite-amplitude modulations of the turbulent mean flow profile naturally exist in form of turbulent large-scale streaks, whose statistical footprints in the streamwise and time average are the secondary flow cells. In particular, large-scale streaks and the secondary currents exist over intermediate time intervals in open channels over flat smooth walls and macroscopically flat sediment beds likewise, highlighting that the currently observed mechanism does not require an initial spanwise disturbance of the sediment bed to trigger the lateral modulation of the flow field as in the linear stability model.

Despite the conceptual differences between the model of Colombini (1993) and the interaction of large-scale streaks and sediment ridges observed in the current study, interestingly, the organization of sediment ridges and secondary currents show remarkable similarities. The characteristic spacing of sediment ridges turns out to be in a very similar range ${\lambda _{h,z}}/{H_f}\approx 1-1.5$, which is comparable to values ${\lambda _{h,z}}/{H_f}\approx 1-3$ found in experiments (Wolman & Brush 1961; Ikeda 1981; McLelland et al. 1999). Also, a comparison of the secondary flow cells that are obtained in Colombini's analysis (recomputed, as plots are not provided in the article) with those presented in figure 16 shows a qualitatively similar structure concerning shape, wall-normal location of the centre of rotation and lateral spacing of the secondary currents. In particular, in both situations the lateral wavelength of the secondary currents and of the sediment ridges results from the momentum balance in the flow and is essentially independent of the bed properties. Therefore, we believe that a more detailed analysis of the implications of the present results for the linear stability of a sediment bed should be performed in the future.

The observed behaviour of the large-scale streaks and the interaction of flow structures at different wall distances agree fairly well with the conceptual model of Jiménez (2018, § 5.6 and references therein) in canonical smooth-wall-bounded flows according to which the main role of the lower domain boundary is to provide a mean shear that feeds turbulence with energy. The model turns away from the idea that turbulent structures are almost exclusively produced in the buffer layer where the mean shear peaks, and it assumes that structures can form at all wall distances. This is in contrast to hairpin-vortex-based models (Adrian, Meinhart & Tomkins 2000; Adrian 2007), in which near-wall generated structures are assumed to migrate outwards to form larger structures. It is therefore inherent to this class of models that the preferential ‘direction of causality’ is directed away from the wall. In the conceptual model of Jiménez, on the other hand, the lack of a single wall-parallel layer in which structures are preferentially generated reposes the question how scales at different wall-normal distances interact and whether there exists a net preferential direction of this interaction. In Jiménez (2018), this question is discussed after reanalysing the data obtained by Flores & Jiménez (2010) and it is concluded that, for the minimal simulation boxes considered therein, there is such a preferential direction in the logarithmic layer bringing information from outer regions towards the wall (i.e. ‘top–down’), rather than the other way around.

Support for the ‘top–down’ concept comes, amongst others, from direct numerical simulations which show that logarithmic- and outer-layer coherent structures are relatively persistent in situations in which the near-wall self-sustaining regeneration cycle is effectively suppressed, either artificially (Mizuno & Jiménez 2013) or through the presence of significant roughness (Flores et al. 2007). Our simulations support these findings, highlighting the fact that, even under severe perturbations of the near-bed flow through the presence of mobile particles and a consequently destroyed buffer-layer regeneration cycle, the regular organization of large-scale velocity streaks in the outer layer remains essentially unaffected. The large-scale structures in our simulations are seen to span the entire depth of the channel reaching down to the bed, and the discussed time lag of flow organization in the near-bed region compared with that of the large-scale streaks in the channel bulk agrees with the preferential direction of information propagation towards the lower boundary found by Jiménez (2018). Also, individual large-scale streaks are seen to intermittently break up (see, for instance, case M850 between $t=10{T_b}$ and $t=20{T_b}$) in a kind of bursting process comparable to that observed in Flores & Jiménez (2010). After such break-up events, the large-scale streaks are seen to first regenerate in the channel bulk essentially uncorrelated to the flow in other regions of the channel, while the flow near the bed later organizes in similar structures with the discussed time lag (cf. figure 13). This further strengthens the idea that coherent structures can indeed form at all distances from the bed/wall autonomously through the action of the local shear.

7. Conclusion

In the current study, we have investigated the role of turbulent large-scale velocity streaks in the formation cycle of sediment ridges in open channels and their connection to mean secondary currents. To this end, four direct numerical simulations with fully resolved spherical particles representing the immersed mobile sediment have been performed, supplemented with two reference simulations of single-phase smooth-wall open channel flow in matching domains and for comparable values of the friction Reynolds number ${Re_\tau }$. Simulations were performed in three different domain sizes with streamwise and spanwise periods ${L_x}/{H_f}=2-12$ and ${L_z}/{H_f}=2.5-16$, respectively. The variation of the domain size allows us to investigate pairs of ridges and large-scale streaks either isolated or in interacting groups. The friction Reynolds number has been varied in a range ${Re_\tau }=250-850$, from which case M850 features with ${Re_\tau }=850$, to the best of the authors’ knowledge, the highest Reynolds number ever reached in a direct numerical simulation with fully resolved particles.

The current simulations verify assumptions according to which sediment ridges can form out of an instability process between a mobile sediment bed and a turbulent stream that is devoid of sidewall-induced secondary currents (Ikeda 1981). The observed mean spanwise spacing ${{\lambda _{h,z}}/{H_f}=1-1.5}$ of sediment ridges in the current simulations is comparable to values ${\lambda _{h,z}}/{H_f}\approx 1-3$ found in experimental studies (Wolman & Brush 1961; Ikeda 1981; McLelland et al. 1999). The currently observed lateral sediment ridge wavelength and the overall structure of the secondary flow cells show similarly good agreement with predictions obtained by means of linear stability analysis (Colombini 1993), even though the exact physical mechanism underlying the therein proposed instability conceptually differs from the currently observed formation process. The qualitatively similar results therefore motivate a more detailed investigation of the instability found by Colombini (1993) in comparison with the results obtained in the current simulations in a future study.

The typical lateral spacing of sediment ridges in the present work is found to be a consequence of the preferential spanwise organization of turbulence in large-scale streamwise velocity streaks in the core of the fluid-dominated region at comparable spanwise wavelengths. Interestingly, the premultiplied energy spectra over developed sediment ridges do not significantly differ from their counterparts in single-phase smooth-wall channels in regions sufficiently away from the bed, in good agreement with data from rough-wall channel flows in which only a certain roughness layer in the vicinity of the wall is modified while large-scale structures remain essentially unaffected (Flores & Jiménez 2006; Flores et al. 2007).

Even though the bed in none of the current simulations can be classified as fully rough, the development of a bedload layer of thickness ${O}(D)$ above the bed in which intense particle transport takes places is seen to destroy the near-wall self-sustaining regeneration cycle just as in the context of fully rough boundaries (Jiménez 2004; Cameron, Nikora & Coleman 2008). Condition for this to occur is that particles are sufficiently large compared with the buffer-layer streak spacing, as in the current case M850 where ${D^{+}}\approx 30$. Irrespective of the existence of the characteristic buffer-layer cycle, large-scale streaks are seen to be correlated to those large-scale structures into which the wall-/bed-shear stress organizes, an effect that has been discussed earlier in the context of canonical closed channel flows (Del Álamo & Jiménez 2003; Jiménez et al. 2004).

The organization of the bed shear in laterally varying regions of high and low stresses causes a spanwise varying ability of the flow to erode sediment. It is consequently observed that at the onset of particle motion, significantly more particles are eroded in regions below large-scale high-speed streaks rapidly leading to the formation of a local trough, while local sediment ridges are preferentially located below large-scale low-speed streaks. The streamwise particle flux is accordingly higher in large-scale high-speed regions, where turbulent sweeps reside that are mainly responsible for particle erosion (Gyr & Schmid 1997).

An important implication of the formation process found in this study is that the initial appearance of sediment ridges is effectively controlled by large-scale streaks and their associated Reynolds stress-carrying structures. The preferred direction of information propagation has indeed been found to be oriented from the bulk flow towards the wall, by evaluating temporal cross-correlation data between the bulk flow and the bed. Sediment ridges and troughs adapt with a time delay to changes in the organization of the large-scale streaks, indicating that the coherent structures at different distances from the wall interact in a ‘top–down’ fashion. This agrees with the observations of Jiménez (2018) in minimal flow simulations of canonical closed channels that the ‘preferential direction of causality’ is directed towards the wall in the sense that Reynolds stresses near the wall are correlated to those away from the wall at earlier times.

The formation of new large-scale structures in the bulk of the channel has been observed to occur after intermittent break-ups of ‘older’ large-scale streaks, in a process that is consistent with the log-layer streak bursting reported by Flores & Jiménez (2010). Our simulations give thus further numerical evidence that coherent structures can form at significant distances to the lower domain boundary independently only due to the action of the local shear, in accordance with the model outlined in Jiménez (2018, §5.6 and references therein).

Finally, we have discussed how the instantaneous coherent structures correlate with the commonly observed mean secondary flow pattern in the presence of sediment ridges. Here, we have shown that the laterally alternating regions of high and low streamwise mean velocity over troughs and ridges, respectively, are the statistical footprints of the streamwise-elongated large-scale streaks when averaged over the streamwise direction and short time intervals of ${O}(10{T_b})$. The related upward and downward oriented secondary flow over ridges and troughs, respectively, are the collective effect of individual ejections and sweeps adjacent to the large-scale streaks. In that sense, the observed secondary flow cells are closely linked to the conditional quasi-streamwise rotational motions found between conditionally averaged log-layer streaks (Del Álamo et al. 2006; Lozano-Durán et al. 2012). This point is further strengthened by the observation that, over the limited time interval under consideration, a very similar secondary flow pattern with comparable intensity and lateral spacing is obtained also for the smooth-wall single-phase simulations. For sufficiently long averaging time intervals, the mean secondary flow pattern is only maintained for spanwise heterogeneous flow configurations (e.g. duct flow with sidewalls, Pinelli et al. 2010), while it is absent for smooth-wall channel flows. For fixed spanwise heterogeneities, such as alternating roughness stripes on the lower domain boundary, this effect is commonly explained by a ‘locking’ of the mean spanwise position of instantaneous large-scale flow structures due to the spanwise heterogeneities (Kevin et al. 2019a).

For self-formed sediment ridges, the situation is less clear as they are mobile themselves, and in the current study their position is seen to be relatively sensitive to rearrangements of the large-scale flow structures. In the current simulations, we cannot give a conclusive answer on the long-time evolution of sediment ridges as the observation time is limited in most of our cases due to the appearance of larger transverse ripple-like bedforms after ${O}(100)$ bulk time units. Experimentally observed sediment ridges, on the other hand, are usually described to be time persistent and to not significantly propagate laterally (Nezu & Nakagawa 1993). However, it should be kept in mind that the experimentally studied bedforms are usually subject to a stabilizing effect of lateral sidewalls in laboratory flumes, which are absent in the current channel configuration. To this end, it would be highly desirable to determine the effect of lateral domain boundaries on the formation and stability of sediment ridges in the future, which will then provide direct comparability between experiments and numerical simulations.