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

The role of turbulent large-scale streaks 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 preferential spanwise organization of turbulence in large-scale streamwise velocity streaks. Sediment 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 a `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 Jimenez (J. Fluid Mech., vol. 842, 2018, P1, section 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 organized large-scale streaks.


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 side-walls. 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. Hinze (1967Hinze ( , 1973 later found that secondary flow cells appeared also in flows over bottom-walls covered with laterally alternating smooth and rough regions, centered over the roughness transition. He proposed an interaction between secondary momentum transport and the local production-dissipation balance of the turbulent kinetic energy, through which high-turbulent fluid is carried out of regions with an excess in turbulence production while, simultaneously, low-turbulent fluid is brought into it and vice versa for regions with dominating dissipation. In hydraulic flows, secondary flow cells mutually interact with the mobile sediment bed, giving rise to long, essentially streamwise-aligned sediment ridges. 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 and 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 and Brush, 1961) and field observations of dried river beds (Culbertson, 1967;Karcz, 1967), but there was disagreement on the actual origin of the process. Nezu and Nakagawa (1984) proposed a cascade-like process, according to which the secondary flow cell induced by lateral domain boundaries such as river banks or the side-walls of a laboratory flume give rise to a first ridge in its close proximity, that, in turn, will then trigger a next ridge further apart with the corresponding secondary current and so on. In a later study, however, Nezu and Nakagawa (1989) could observe ridges and secondary currents simultaneously rising at different spanwise locations all across the domain. Ikeda (1981) followed a different approach and conjectured the underlying mechanism to be a flow-sediment bed instability that is able to generate sediment ridges and the according secondary flow cells. Although side-wall induced secondary currents may still influence the formation of sediment ridges, they are no necessary condition for their evolution and ridges can accordingly form in absence of sidewalls likewise. These considerations were later confirmed by Colombini (1993), who showed by means of a linear stability analysis that a two-dimensional turbulent base flow in an infinitely wide channel is, indeed, linearly unstable with respect to infinitesimal spanwise disturbances of the sediment bed. For his model, he closed the Reynolds-averaged Navier-Stokes (RANS) equations with a non-linear eddy-viscosity model originally proposed by Speziale (1987) which allows, in contrast to simpler linear closures, the required anisotropy of the Reynolds stress tensor (Speziale, 1982). For the set of investigated parameters, a most amplified wavelength of λ z = 1.3H f could be determined. Later, Colombini and 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 and Studerus, 1979). In fact, the two different instability mechanisms allow a useful classification of secondary currents in strip type (e. g. alternating roughness stripes) and ridge type roughness (e. g. topographical variation of the lower wall) or a combination of both (Colombini and Parker, 1995;Hwang and Lee, 2018).
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 and Christensen (2010) and Barros and Christensen (2014). Secondary flows seem to represent a robust phenomenon that is ubiquitous in flow configurations that feature a significant laterally heterogeneity. For example, secondary currents have been observed in the flow over straight riblets (Goldstein and Tuan, 1998), converging/diverging riblets (Kevin et al., 2017(Kevin et al., , 2019a, spanwise alternating roughness stripes (McLelland et al., 1999;Wang and Cheng, 2005;Wangsawijaya et al., 2020), flow over transverse alternating non-/hydrophobic roughness (Stroh et al., 2016;Türk et al., 2014), streamwise aligned artificial ridges of different cross-sectional shape (Medjnoun et al., 2020;Stroh et al., 2020a;Vanderwel and Ganapathisubramani, 2015;Vanderwel et al., 2019;Wang and Cheng, 2006;Zampiron et al., 2020) or combinations of the above as in Stroh et al. (2020b). The majority of the listed studies focusses on the general structure of mean secondary motion over the respective heterogeneities, their sense of rotation (Anderson, 2019;Yang and Anderson, 2018) or the influence of the lateral spacing of the heterogeneities on their intensity and arrangement (Chung et al., 2018;Wangsawijaya et al., 2020).
Less attention has been given to the transient nature of secondary flows. Recently however, fullyresolved direct numerical simulations (DNS) and experiments have allowed to shed some light on instantaneous coherent structures and their relation to the mean secondary flow. In this spirit, Uhlmann et al. (2007) and  were the first to show that the characteristic eight-vortex secondary flow pattern in a square duct is indeed a direct consequence of instantaneous large-scale structures, while the preferential location of the buffer-layer streaks and quasi-streamwise vortices along the four solid side-walls determines the statistics of the mean vorticity distribution. The close correlation between distinct turbulent coherent structures and the mean secondary flow was further supported by , who presented a family of traveling wave solutions whose members produce essentially the same secondary flow patterns when averaged over the streamwise direction (Kawahara et al., 2012). Sakai (2016) later elaborated that similar relations between coherent structures and the statistical mean secondary flow exist for open ducts just as much.
In canonical laterally homogeneous channels, instantaneous quasi-streamwise vortices appear at all scales in connection with alternating streaks of high and low streamwise velocity (Jiménez, 1998(Jiménez, , 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 et al., 1995;Jeong et al., 1997;Schoppa and Hussain, 2002). With increasing wall-normal distance, 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, 2013a), whereas vorticity concentrates further away from the wall in structures of essentially the same scale as in the buffer layer (Jiménez, 2013b). It is thus evident that the role of larger-scale quasistreamwise rotations that sustain the fully turbulent log-layer streaks cannot be overtaken by a single large-scale vortex, but by the collective effect of a large number of locally isotropic small-scale vortices that globally arrange in large-scale vortex clusters, that are anisotropic enough to generate an average upward motion in their middle (DelÁlamo et al., 2006;Jiménez, 2018). While these collective large-scale rotations are hardly detectable in instantaneous flow visualizations, they have been successfully uncovered in conditional averages (DelÁlamo et al., 2006;Lozano-Durán et al., 2012) or in low-pass filtered fields (Motoori and Goto, 2019), proving that the self-sustaining process between streaks and quasi-streamwise rotations is a self-similar cascade (Motoori and Goto, 2021), just as the logarithmic layer is (Flores and Jiménez, 2010;Kevin et al., 2019b).
Recently, it has been argued by Kevin et al. (2019b) that the mean secondary currents conventionally found over heterogeneous bottom-walls are nothing else than the statistical footprint of these conditional or collective quasi-streamwise rollers. The latter authors argue that the role of the lateral bottom heterogeneities is to laterally lock the large scales structures and quasi-streamwise rollers in the vicinity of the roughness transition, such that they get visible in the long-time average. Structures over homogeneous walls, on the other hand, will appear at any spanwise position with the same probability, explaining the absence of their statistical footprint in form of mean secondary currents when averaged over a sufficiently long time interval. Note that spatial locking should be understood in a sense that the mean lateral position does not vary significantly in time, while the streaks instantaneously meander in response to the quasi-streamwise rollers just as in canonical flows (Kevin et al., 2017(Kevin et al., , 2019b. Indeed, such lateral meandering has been detected for a variety of spanwise heterogeneities such as converging/diverging riblets (Kevin et al., 2017), alternating smooth/rough patches (Wangsawijaya et al., 2020) and streamwisealigned artificial ridges (Zampiron et al., 2020). In contrast to homogeneous smooth walls, the lateral spacing of bottom wall heterogeneities allows to control the organization of the flow, especially the spacing and location of mean secondary currents. The above studies show that an optimal lateral spacing of the bottom roughness elements, which enhances the intensity of the secondary currents, corresponds to the typical lateral dimension of the large-scale motions (1-2)H f (Jiménez, 2013a(Jiménez, , 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 and 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 and Nakagawa, 1993), consistent with the lateral dimension of the turbulent large-scale structures with streamwise length O(1-10)H f (Roy et al., 2004;Shvidchenko and Pender, 2001;Gulliver, 1999, 2007), that have been detected as counterparts to large-and very large-scale motions inherent to canonical turbulent flows (Adrian and Marusic, 2012). A general relevance of those 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 and Rodi (1985). It was discovered that the mean secondary currents decay towards the duct center, 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 and 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 traveling 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 (Cameron et al., 2020;Gyr and Schmid, 1997).
Even though it appears conclusive that organized recurrent large-scale structures are responsible for the formation of sediment ridges, it lacks, to the best of our knowledge, clear experimental and numerical evidence of the above elaborated mechanisms. While simultaneous measurements of three-dimensional flow structures and individual sediment motion remain challenging (Wang and Cheng, 2005), direct numerical simulations featuring fully-resolved particles have recently proven useful to study complex interactions between mobile sediment and the different scales of wall-bounded turbulence. So, Chan-Braun et al. (2011) and Mazzuoli and Uhlmann (2017) investigated the flow over fixed arrays of fullyresolved spherical roughness elements from the hydraulically smooth to rough regime and quantified the pressure and forces exerted on the particles. Kidanemariam et al. (2013) showed for moderate Reynolds numbers that the lateral migration of particles in a smooth-walled channel into the nearwall low speed streaks can be attributed to the quasi-streamwise streaks in the buffer-layer. Vowinckel et al. (2017) studied the dynamics of a limited number of particles over a fixed pattern of roughness elements and reported particle clustering in streamwise aligned chains appearing together with mean secondary currents. Recently, Kidanemariam et al. investigated the evolution of transverse bedforms over a thick mobile sediment bed under laminar and turbulent conditions in a series of simulations with fully-resolved sediment grains (Kidanemariam andUhlmann, 2014a, 2017;Scherer et al., 2020). To investigate the minimal length of those patterns, Kidanemariam and Uhlmann (2017) performed a sequence of simulations successively reducing the streamwise box dimension in a similar fashion as the minimal flow units of Jiménez and Moin (1991). That way, it was possible to show that there exists a minimal streamwise domain length below which the domain cannot accommodate any unstable pattern wavelength and consequently no ripple-like bedforms evolve. The correct scaling of this minimal unstable wavelength with the particle diameter D rather than with the fluid height H f was shown by Scherer et al. (2020) for the investigated parameter range, and it was possible to quantify the minimal box length as L x,min ≈ 80D for any ripple-like pattern to rise.
In most of the simulations conducted by Kidanemariam and Uhlmann (2017) and Scherer et al. (2020), sediment ridges were observed to be the first bedforms to evolve rapidly after the onset of particle motion, while the actually investigated transverse patterns required roughly two orders of magnitude longer to develop and eventually take over the smaller-amplitude ridges. The current study aims to provide numerical evidence for the ridge formation process by, first, clarifying the role of large-scale turbulent structures in the generation process of sediment ridges and, second, by investigating the relation between coherent structures and mean secondary currents in the vicinity of the bedforms. 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. § 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.

Numerical method
The numerical scheme used in the current study has been first used for the simulation of sediment transport in a carrying fluid in Kidanemariam and 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 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 disregarding the continuity equation, then projecting the velocity to the space of solenoidal velocity fields. Back-and forthtransformations 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 follow 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, from which the latter occur frequently in bedload dominated sedimentary flows, are modeled 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 ∆ is below this force range ∆ 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 δ c = ∆ c − ∆ 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, Figure 1: Sketch of the physical system analyzed in the multiphase simulations. Mean flow and gravity are pointing in positive x-and negative y-direction, respectively. The mean flow profile u f xzt = ( u f xzt (y), 0, 0) T is shown in blue, while the green curve represents the wall-normal variation of the mean fluid shear stress τ f (y). h f and h b are the local instantaneous fluid and bed height, respectively. Particles are colored depending on their location: bed particles are colored in black, interface particles in orange to yellow with increasing wall distance and transported particles are indicated by white color. For the sake of comparison, the bed at the downstream end of the domain is overlain by the two-dimensional surface that represents the fluid-bed interface as defined in § 4, the color ranging from dark to bright blue with increasing height.
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 µ c . For a more detailed presentation of the described model and an extensive validation study, the interested reader is referred to Kidanemariam and 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 ,µ c ) plus the force range ∆ c . As pointed out by Kidanemariam and Uhlmann (2014b), a quantity that is often more accessible than the force parameters is the dry restitution coefficient ε 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 ε d and k n , such that the actually used set of model force parameters reads (k n ,ε d ,c t ,µ 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 µ c = 0.4-0.5 and the friction coefficient is set to c t = c n . The force range ∆ c is set equal to the uniform grid spacing ∆x. The dry restitution coefficient is ε d = 0.3 except for case M 850, where we have used a larger value ε d = 0.9. Note that the variation of ε 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.

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 wallnormal directions, respectively (Kim et al., 1987). The physical system simulated in the two-phase cases is sketched in figure 1. The Cartesian coordinate system has its origin on the bottom wall of the channel, such that the components of an arbitrary spatial position vector x = (x, y, z) T are measured along the streamwise (x), wall-normal (y) and spanwise (z) direction, respectively. Accordingly, the fluid velocity vector field at position x has components u f (x) = (u f , v f , w f ) T . Subscripted letters 'f ' and 'p' indicate fluid-and particle-related physical quantities. For the statistical analysis, we introduce the following   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 ∆x = ∆y = ∆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. 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 τ , D + , H f , H b , θ, ∆y + ) are computed as an average over the entire simulation period.
averaging scheme where (1a) and (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, • and • are the corresponding temporal and spatial fluctuations. We further introduce the following decomposition based upon streamwise averaging only: where u f = u f x are the spatial fluctuations with respect to the streamwise average. 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 = |g|. The simulation domain is periodically repeated in the wall-parallel x-and z-directions with fundamental periods L x and L z , respectively. 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ỹ = y − H b .
The mean wall-shear stress τ w is evaluated at the wall-normal location of the virtual wallỹ = 0 by interpolating the pure fluid stress 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 δ ν = ν f /u τ , respectively, where u τ = τ w /ρ f is the friction velocity. Following the general conventions, we shall denote quantities that are non-dimensionalized with u τ and/or ν f with a superscript (•) + and refer to them as scaled in inner or wall units. From these three length scales, we derive the friction Reynolds number Re τ = H + f = H f /δ ν , the particle Reynolds number D + = D/δ ν and the relative submergence H f /D, respectively. Let us further introduce the bulk Reynolds number as 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 twophase system. Here, we choose the density ratio of the two phases ρ p /ρ 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/ν f , where u g = (ρ p /ρ f − 1)|g|D) is the gravitational velocity. The squared ratio of the Galileo and the particle Reynolds number θ = (D + /Ga) 2 = (u τ /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 θ c = 0.03-0.05 marks the onset of sediment erosion in a turbulent stream, that slightly depends on the Galileo number (Franklin and Charru, 2011;Soulsby et al., 1997;Wong and Parker, 2006).
Tables 1 and 2 summarize the relevant parameters of the simulations. As a short-hand 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 ∈ {2, 5, 12} and L z /H f ∈ {3, 16}, respectively, and friction Reynolds numbers 250 Re τ 850. Let us remark that for the smooth-wall reference simulations, the values of the friction Reynolds number are set to Re τ = 210 and Re τ = 650. These values match those of cases L250 and M 850, respectively, but when the bed is still at rest. As the bed evolves, the values of Re τ increase as a result of the sediment bed modulation.
The box dimensions L + x and L + z are sufficiently larger than those of Jiménez and Moin (1991)'s 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 and 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 ≈ 3. Surely, the comparison to our open channel is limited in the vicinity of the free-surface, but should give a good estimate for the rest of the domain (Bauer et al., 2021).
The preparation of each individual simulation follows a careful procedure described in detail by Kidanemariam and Uhlmann (2014a), in which first a pseudo-randomly arranged sediment bed is created through settling of a large amount 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 ≈ 26-28, while the Galileo number is chosen a priori such that the a resulting Shields number is sufficiently larger than its critical value to allow for particle erosion. 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 has been stated earlier, ridges appear shortly after particles are released, but eventually, transverse bedforms appear and dominate the evolution process after 100-200T 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 ≈ 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 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 M 650 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 and 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 M 250 is found at the same parameter point as case H6 of Kidanemariam and Uhlmann (2017), but represents an own independent simulation conducted in the context of this study. Due to technical reasons, data in the first roughly 5T b of the simulation are not available for the study. Furthermore, to the best of the authors' knowledge, the newly presented simulations L250 and M 850 represent the largest number of fully-resolved particles simulated and the highest ever obtained Reynolds number in 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 M 850, including a number of around 16.8 billion grid nodes.

Definition of the fluid-bed interface
The determination of the two-dimensional interface that separates the domain in two distinct regions, that are, 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 in two contributions as 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 andUhlmann (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 bases 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 (|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 = (ρ p − ρ f )|g|πD 3 /6 is the submerged weight of a single spherical particle. All particles that fulfill 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 α-shape algorithm of Edelsbrunner and Mücke (1994). While conceptually similar to a conventional convex hull around a set of points, the α-shape allows non-convexity for length scales over some threshold radius α (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 centers 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 1.

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 2 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(10H f ) (Wolman and 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 side-wall induced secondary currents (Colombini, 1993;Ikeda, 1981).
To give a first impression on the comparable lateral organization of bed and large-scale structures, we have supplemented to figure 2 instantaneous visualizations of the Reynolds stress carrying ejection (u f < 0,v f > 0) and sweep structures (u f > 0,v f < 0) (collectively termed as Q − 's) introduced by Lozano-Durán et al. (2012) as a generalization of the classical quadrant analysis (Wallace et al., 1972) to three dimensional objects. According to the former study, coherent ejection and sweep structures are connected sub-domains for which It is remarkable that regions of intense ejections seem to reasonably well correlate 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 and 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, large-scale ejections live in the low-speed streaks whereas sweeps populate the high-speed streaks of the logarithmic layer, thus we can equivalently 0 1 2 0 20 40 60 80 conclude that ridges are found below the large low-speed streaks and troughs accordingly next to the high-speed streaks, which has been also verified by comparable visualizations (plots not shown). To avoid confusion, the term 'streak' will always refer to the velocity structures of the logarithmic layer in the remainder of this work if not otherwise stated, whereas the nearest-wall structures will be always explicitly termed as buffer layer structures.

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 3 shows the time evolution of the fluctuations of the streamwise-averaged sediment bed h f (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. 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 (Colombini, 1993;Ikeda, 1981;Wolman and Brush, 1961). 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 those 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 side-walls in low to moderate aspect ratio flumes and underlie the side-wall 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 ≈ 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 4 provides the time evolution of the mean bedform geometry and the particle transport rate. Panels 4(a,b) show the development of the mean pattern height expressed by the root mean square of the sediment bed height fluctuations σ h,z and that of the lateral spacing in terms of the mean bedform wavelength λ 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 = 10T b and t = 20T b , leading to a net reduction of the number of ridges. The bedforms continue growing in amplitude predominantly during the time interval between t = 20T b and t = 40T b . The lower growth rate in the first 5 bulk time units of the three low Reynolds number cases S250, M 250 and L250 is a consequence of the lower Shields numbers when compared to the high Reynolds number case M 850, in which the high Shields number causes a stronger erosion and thus a faster pattern growth. For the sake of comparison, panel 4(a) additionally provides the evolution of the spanwise-averaged root mean square of the sediment bed height fluctuations σ 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  (Colombini, 1993). The clearly lower attained final value of 0.76H f in case M 250 is in agreement with the corresponding space-time plot in figure 3 which shows that the bed in the final phase indeed features three ridge patterns.
Panel 4(c) contains the temporal evolution of the streamwise particle flux q p,x xz (t) (cf. equation (13a) 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 and Müller (1948) in the modified form of Wong and Parker (2006), viz.
where we have chosen a critical Shields number θ c = 0.034 (Soulsby et al., 1997). This further supports earlier findings in the context of transverse patterns by Kidanemariam and Uhlmann (2017) that relation (3) 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 q p,x xz (t)/q ref during the initial 20 bulk time units is not a feature of the particle flux evolution q p,x xz (t) which increases in an essentially monotonic way in this time window, before it settles to a quasi-stationary plateau (cf., for instance, Kidanemariam et al., 2021). Instead, it is observed that relation (3) 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 in a layer of thickness O(D) above the bed in which particles are transported by rolling, sliding and small jumps (saltation) without loosing the contact with the bed for longer time intervals (van Rijn, 1984). This is verified in panel 4(d ) which shows the wall-normal profile of the mean particle flux density φu p xzt (y) (cf. equation (15)), 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 θ in case M 850 leads to a thickening of the bedload layer by around 1D compared to the remaining cases. We conclude that in this latter case, a layer of 0.25H f thickness above the interface is characterized by a non-negligible particle flux.

Turbulent mean flow
The presence of the mobile sediment severely modifies the mean fluid flow characteristics. Figure 5 shows the mean fluid velocity profile u f xzt (y) compared to that of the single-phase cases. First, the sediment-laden cases feature inflection points in the near-wall region as a consequence of the underlying porous beds, which are susceptible to a Kelvin-Helmholtz instability (Jiménez et al., 2001). Generally, the near wall region is strongly modified in the particle-laden simulations in the sense that the mean shear d u f xzt /dy reduces to negligible value when approaching the mean fluid-bed interface atỹ = 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 and Uhlmann, 2017). The mean velocity u f xzt increases significantly slower with increasing distance to the bed in all particulate cases, leading to a velocity deficit compared to the smooth-wall cases, an essential feature arising in turbulent flows over rough-walls due to an enhanced friction (Chan-Braun et al., 2011;Flores and Jiménez, 2006). Farther away from the wall, in particular forỹ/H f > 0.5, the opposite is true, i.e. the particle-laden flows feature a higher mean velocity than their smooth wall counterparts. This is expected as the reduced mass flow rate in the vicinity of the sediment bed as a consequence of the fluid-particle interactions has to be compensated by a higher flow rate in the outer flow to maintain the constant bulk flow rate q f .
As a consequence, the slope of the velocity profile in the semi-logarithmic visualization deviates from that in the logarithmic layer of the smooth wall channels. For the low Reynolds number simulations, the deviation is rather weak due to the comparably lower particle flux, such that the effect of the sediment bed can be mostly compensated by the conventional shift of the log-layer velocity profile by an offset ∆U + known from analysis of rough wall turbulence Using a standard choice of von Kármán's constant κ = 0.41 and the coefficient B = 5.2, we obtain ∆U + = 3.59 (M 250) and ∆U + = 3.73 (L250) as values for the roughness function. For the high Reynolds number case M 850, on the other hand, it is observed that the slope in the semi-logarithmic representation differs markedly compared to the smooth-wall reference solutions, such that a simple vertical shift will not lead to a satisfactory match of the profiles. 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 exemplary for cases M 250 and M 850, respectively, together with the corresponding single-phase smooth-wall data. 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.
In all simulations, the typical buffer-layer peak of the streamwise normal stress u f u f xzt is visibly reduced when compared to the smooth-wall cases. This is not unexpected, as the mean shear which represents the source of energy for the fluctuating field has been observed to be lower in the vicinity of the bed (Chan-Braun et al., 2011). For the low Reynolds number cases M 250 and L250 (plot of the latter case not shown), the modulation w.r.t. the smooth-wall reference case is moderate in the sense that the peak is, even though damped, clearly detectable and found at essentially the same wall-normal position. In the high Reynolds number case M 850, on the contrary, the differences are more pronounced: the intense particle motion in the bedload layer of particles with significant size results in the complete absence of the characteristic buffer-layer peak of u f u f xzt . Note that the essentially same phenomenon appears in flows past fully-rough walls where the roughness height scaled in wall-units is also large enough to destroy the self-sustained near-wall cycle (Flores and Jiménez, 2006;Mazzuoli and 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 ≈ 10δ ν ), which has been observed to be sufficiently small for particles to be transported by the buffer-layer structures rather then to destroy the near-wall cycle (Kidanemariam et al., 2013). Deviations for the remaining normal stresses v f v f xzt and w f w f xzt are also observable, but are clearly weaker and visually restricted to a layer of thicknessỹ ≈ 0.11H f (ỹ ≈ 3D) above the bed in the low Reynolds number cases andỹ ≈ 0.17H f (ỹ ≈ 4.5D) in the high Reynolds number case, respectively. The insets in figure 6 show the Reynolds stress decay in the near-bed region in double-logarithmic scaling. It is seen that all Reynolds stresses are higher for the multi-phase simulations as there is no strict boundary condition that enforces them to vanish atỹ = 0 and thus they are observed to settle at small but finite levels when entering the bed whereas their counterparts over smooth walls decay at the theoretically predicted decay rates when approaching the physical wall (cf. Pope, 2000, p.284).
Our conclusions from the afore-discussed global flow statistics highlight the fact that a general mechanism for the creation of sediment ridges cannot be due to the dynamics of the usual buffer-layer structures, as these readily disappear once particles are released in case M 850. Also, natural flows over fully-rough beds with D + 1 lack a viscous near-wall cycle just as much. It is thus conclusive to focus in the following on the dynamics of the outer flow structures as the main candidates to drive ridge formation. This idea is further strengthened by our observation that ridges feature a very similar lateral spacing of λ h,z /H f ≈ 1-2 in the current small to intermediate Reynolds number simulations and in the high Reynolds number experiments.
In the following, we first compare whether the mobile sediment bed has a detectable influence 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 M 250 and M 850 are analyzed in figure 7 in the center of the clear-fluid region (ỹ/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 (c) (d ) Figure 8: Time-averaged and streamwise-integrated premultiplied energy spectra k z φ αα t (λ x , y, λ z )dk x H 2 f /u 2 τ as a function of the inner-scaled spanwise wavelength λ + z and wall distanceỹ + , respectively. (a,c) streamwise φ uu t and (b,d ) wall-normal φ vv t energy spectra for cases (a,b) M 250, (c,d ) M 850. Colored isolines are 0.2(0.2)0.6 times the maximum value of the respective energy spectra, while gray-shaded areas indicate the same quantities determined for the smooth-wall reference simulations (a,b) L250 smooth and (c,d ) M 650 smooth , respectively. The mean fluid height H + f = Re τ and the spanwise domain period L + z of the particle simulations are marked by colored dashed lines. The dashed black line refers to the wall-normal distance at which the mean particle flux density φu p xzt attains its maximum (cf. figure 4(d )).
whereû f (k x , y, k z , t) = F(u f (x, t)) is the Fourier transform of the fluctuating field in the two wallparallel directions, with the asterisk indicating complex conjugation. k x = 2π/λ x and k z = 2π/λ z are the streamwise and spanwise wavenumber-wavelength pairs. Velocity spectra for the remaining velocity components are defined accordingly.
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 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, , 2018Mazzuoli and 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 streamwise (left column) and wall-normal energy spectra (right column) is clearly identifiable in figures 8 and 9, where the premultiplied streamwise-integrated energy spectra are presented as functions of the lateral wavelength λ z and the wall-normal distanceỹ 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. Let us first conclude that the medium domains 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, 2013b), which further supports our observations in § 3. Also, it is consistent with our claim that the relative particle size D + together with the high particle transport rate in the bedload layer of case M 850 cause the immediate breakdown of the buffer layer regeneration cycle. We observe that the near-wall peak in the streamwise spectra at λ + z ≈ 100 associated with the buffer-layer streaks is still present for the low Reynolds number case M 250 and kinetic energy is concentrated around this region whereas it is absent in case M 850, 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 λ z /H f ≈ 1-1.5 (λ + z ≈ 1000) that reflects the large-scale motions at higher Reynolds numbers.
Naturally, the wall-normal energy spectra in smooth-wall channels has to concentrate in regions away from both, the lower wall and the upper free surface, as the boundary conditions prevent intense wallnormal motion in both regions. Again, the spectra of the low Reynolds number case M 250 almost collapse with that of the smooth-wall dataset with slight deviations only due to the slight difference in Re τ , highlighting that the distribution of energy over the scales and the wall-normal layers is essentially the same as over the smooth wall. In the high Reynolds number case M 850, 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. As can be seen in the outer-scaled figure 9(b), the spectra contours 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.

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 w.r.t. the streamwise average, u f (z, t). Note that streamwise averaging filters out lateral meandering of structures while retaining the footprint of the large streamwise elongated structures. In fact, the streamwise-averaged low-momentum and highmomentum 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 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  (Jiménez, 2013b), which is expected to disappear for sufficiently long time averaging intervals leading to statistical homogeneity. We shall discuss the role of these 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 ≈ 3, but that comparable regular organization of the streaks over intermediate time intervals similarly occur in domains as wide as L z /H f ≈ 16. Also, it shows that the streaks are of significant length, as they are clearly detectable even in streamwise averages over long streamwise extents L x /H f ≈ 12.
Intermittently, the amplitude of the high and low-momentum regions is seen to reduce as, for instance, around t = 40T b in case M 650 smooth (panel 10(b)) or in a period 10 t/T b 20 in case M 850. It appears likely that at least part of these 'events' are related to regularly occurring bursting events in which the large streaks bend and eventually break, just as their viscous counterparts in the buffer layer do (Flores and Jiménez, 2010). Interestingly, there are events after which streaks recover at practically the same lateral position, while in other situations as after t = 45T b in case M 250 (panel 10(a)), high and low-momentum regions reorganize and partly even change the number of 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 largest 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 et al. (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 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 streaks trigger the first ridges and not vice versa. Note that this does not rule out the possibility of a feed-back of developed ridges on the mobility or meandering tendency of the 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, 2013b). Also, large-scale instantaneous quasi-streamwise rotations are almost impossible to detect by the classical gradient-based vortex detection methods such as the λ 2 -criterion of Jeong and 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 Goto (2019, 2021). The advantage x 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 cut-off 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 et al., 2016): Therein, ∆ i (i = x, y, z) is the filter cut-off width in the three Cartesian directions and x is the innerconvolution 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. 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 ∆ 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 [∆ x , ∆ y , ∆ 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 equation (6) in wall-parallel (x, z)-planes whenever we are mainly interested in the structure's wall-parallel extensions, using filter widths 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 and 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δ ν ), 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 and Jiménez, 2003). Even though the near-wall regeneration cycle itself is capable of functioning autonomously (Jiménez and Pinelli, 1999), large-scale influences have nonetheless been observed to modulate regions even close to the wall . Toh and 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δ ν ) into larger essentially H f -scaled streaks and the corresponding outer large-scale streaks atỹ/H f = 0.3 can be observed exemplary for the smooth-wall case M 650 smooth in visualizations of unfiltered and corresponding spatiallyfiltered 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 τ w , let us introduce the instantaneous bottom-shear stress τ b (x, z, t) = ν f du f /dỹ|ỹ =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, from which 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 high-speed streaks at the upper end of the logarithmic layer,ỹ/H f = 0.3, and are of comparable lateral extension, even though the streamwise extension differs in particular for the patchy structure with center around 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 M 650 smooth in the streamwise-averaged framework, let us define the instantaneous two-point correlation coefficient as follows: 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 center of ρ u,τ t (ỹ/H f = 0.5) ≈ 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. 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 M 850 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).

Streak-ridge interaction
The above arguments imply that particle erosion and hence bed topography are linked with the largescale 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 λ z /H f ≈ 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 λ h,z /H f . That, indeed, regions of high and low streamwise velocity are spatially correlated with the bed contour is visualized in figure 12, exemplary for case M 850. Here, the space-time evolution of the sediment bed fluctuation h b (z, t) and the respective fluctuation of the streamwise velocity component atỹ/H f = 0.5 are repeated for convenience, supplemented by similarly presented data for the streamwise and spanwise particle flux components q p,x and q p,z , respectively. The observed dependences can be summarized as follows: regions of high momentum which are the streamwise-averaged analogue to the large-scale streaks lead to spanwise alternating zones of locally increased erosion due to turbulent sweep structures that are located in the 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 in these regions. The effect is further supported by the action of the lateral particle transport q p,z , which is predominantly directed in such sense that sediment is transported from the regions of high streamwise particle transport (q p,x > 0) to areas with lower streamwise particle flux, i.e. q p,x < 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 form of a local trough. Corresponding opposite relations are found for low-momentum regions, ridges and regions of low streamwise particle transport q p,x < 0.
The lateral particle flux is found to be two orders of magnitude smaller than the corresponding streamwise particle flux q p,x in all simulations. It is therefore conceived that the initial sediment ridges  (c)); (c) The streamwise-averaged particle flux fluctuation q p,x /q x,rms , with q x,rms (t) = ( q p,x q p,x z ) 1/2 . (d ) Same data as in (c) is shown in the background as gray map, supplemented with red (blue) dots marking regions of lateral bed growth (decrease), that are, regions with vanishing lateral particle flux q p,z = 0 and negative (positive) lateral gradient ∂ z q p,x (z, t).
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 4(c)).
The visually identified correlation between large-scale structures and the bed contour in the multiphase simulations can be quantified with the aid of a two-point cross-correlation as defined in equation (7), replacing the wall-shear stress τ b by the sediment bed height fluctuation h b . 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 centerlineỹ/H f = 0.5 (figure omitted). This holds, however, not true in the phase between t = 10T b and t = 30T 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 figures 12(a,b): between t = 10T b and t = 20T 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 streaks are believed to break up in this phase. Two large-scale streaks recover between t = 20T b and t = 25T b , but at somewhat different lateral positions. Interestingly, the bed contour is seen to adapt to the new spanwise organization of the streaks, laterally migrating in time towards these spanwise positions, however with a noticeable time delay.
This time lag between the bed evolution and the 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 motions 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 ρ t ua (y, δt) = u f (y, z, t + δt) a (z, t) z / (u f (y, z, t + δt)) 2 z (a (z, t)) 2 z 1/2 .
Panel 13(a) shows the cross-time correlation between the flow field and the bed contour, i.e. a = h b (z, t), while panel 13(b) provides the cross-time correlation between the flow field at varying wallnormal positions and that at a reference height chosen asỹ/H f = 0.5, i.e. a = u f (ỹ/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 both panels, figure 13(a) and 13(b), a time lag is clearly identifiable. Figure 13 our earlier observations that ridges align themselves with the location of the large-scale streaks with a time lag of around 10T b (equivalent to approximately 0.9H f /u τ ). A very similar time delay is observed in figure 13(b) between the dynamics of the flow field in the channel center 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 center which interact in a 'top-down mechanism' with the fluid layers closer to the bed.
It also indicates that the existence or absence of a dominant large-scale velocity streaks in the core of the channel is of direct relevance for the stability and position of the underlying sediment ridges and troughs, respectively. In the following, we shall analyze such events during which a clear signature of high-and low-momentum regions is missing, avoiding streamwise averaging for the moment to show that the modulation of u f in such phases is indeed due to a break-up of the streaks. Figure 14 shows the time evolution of the most-energetic Fourier modes of the streamwise velocity fluctuations, φ uu (k x , y, k z , t), in wall-parallel planes atỹ/H f = 0.5. For the sake of readability, let us refer to modes with wavelength pairs λ x = L x /i and λ z = L z /k (i, k ∈ N 0 ) as (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 u f u f xz (black) and wall-normal energy v f v f 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.
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 colors. 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 large alternating 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 spacetime plots in figure 10, the times during which the clear signature of high and low-momentum 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 M 250 (panel 14(a)), the earlier discussed change from two to only one pair of high and low-speed streaks is clearly visible around t ≈ 50T b , accompanied by a reduction in the total energies u f u f xz and v f v f xz , suggesting that at this time, the streaks indeed break up. Panel 14(b) shows a comparable behavior for the smooth-wall x case M 650 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 streaks as in the previous case. The two states without and with 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 sedimentladen case M 850 provided in panel 14(c), though, 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 = 10T b and t = 20T b ) even this slight plus in energy is lost and they cannot be differentiated from the remaining modes any more. 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 streaks increases first mildly, then stronger until it settles at t ≈ 40T b and dominates the spectra for the rest of the simulation. The dominance of this mode is also seen in figure 15(d ).

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 streaks such that they leave a footprint in the streamwise time-average as upand 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 stream function ψ xt (y, z), which is defined by the requirement that ∇ 2D ψ xt (y, z) = (∂ y ψ xt , ∂ z ψ xt ) T = (− w f xt , v f xt ) T . For the particle-laden simulations, the con- ventionally discussed configuration of primary flow upwelling over ridge crests and downwelling over the troughs is observed with accordingly oriented secondary currents (Nezu and Nakagawa, 1993). A comparison between the cases over mobile sediment beds and the smooth-wall simulation M 650 smooth in panel 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 panels 16 (a,c,d ).
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 streaks, λ z /H f ≈ 1-1.5. An exception is case M 250 for which in the observed time interval intermittently only one pair of 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 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 ⊥ (t) defined from the cross-plane averaged kinetic energy contained in the cross-flow field ( v f x , w f x ) T , viz.
Note that a comparable formulation has been originally presented by Sakai (2016) for quantification of side-wall induced mean secondary flow in open ducts and has been adapted here to the multi-phase configuration. It is seen that u ⊥ (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 ( v f xt , w f xt ) 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 side-wall induced mean secondary flows in open ducts. For instance, Sakai (2016) determined the intensity of secondary flow within a distance z = 1H f from the side-walls 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 streaks (Kevin et al., 2019b;Lozano-Durán et al., 2012). 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 600T 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 = 400T b and t = 500T b .

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 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 and 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ỹ/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 M 850, the time lag has been estimated as approximately 10T b or, equivalently in terms of the eddy turnover time, 0.9H f /u τ .
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 to 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 RANS 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, 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 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 turn out to be in a very similar range λ h,z /H f ≈ 1-1.5, which is comparable to values λ h,z /H f ≈ 1-3 found in experiments (Ikeda, 1981;McLelland et al., 1999;Wolman and Brush, 1961). 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 center 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, 2007;Adrian et al., 2000), 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 reanalyzing the data obtained by Flores and 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 and Jiménez, 2013) or through the presence of significant roughness (Flores et al., 2007). Our simulations support these findings, highlighting 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 largescale 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 to that of the streaks in the channel bulk agrees with the preferential direction of information propagation towards the lower boundary found by Jiménez (2018). Also, individual streaks are seen to intermittently break up (see, for instance, case M 850 between t = 10T b and t = 20T b ) in a kind of bursting process comparable to that observed in Flores and 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.

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 smoothwall open channel flow in matching domains and for comparable values of the friction Reynolds number Re τ . 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 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 τ = 250-850, from which case M 850 features with Re τ = 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 side-wall induced secondary currents (Ikeda, 1981). The observed mean spanwise spacing λ h,z /H f = 1-1.5 of sediment ridges in the current simulations is comparable to values λ h,z /H f ≈ 1-3 found in experimental studies (Ikeda, 1981;McLelland et al., 1999;Wolman and Brush, 1961). 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 roughwall 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 and 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 (Cameron et al., 2008;. Condition for this to occur is that particles are sufficiently large compared to the buffer-layer streak spacing, as in the current case M 850 where D + ≈ 30. Irrespective of the existence or absence 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 and Jiménez, 2003;. 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 low-speed streaks. The streamwise particle flux is accordingly higher in high-speed regions, where turbulent sweeps reside that are mainly responsible for particle erosion (Gyr and 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' streaks, in a process that is consistent with the log-layer streak bursting reported by Flores and 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 characteristic upwelling and downwelling of the mean primary flow over ridges and troughs are the statistical footprints of the streamwise elongated streaks when averaged over the streamwise direction and short time intervals of O(10T b ), while upward and downward secondary motions are the collective effect of individual sweeps and ejections adjacent to the large-scale streaks. In that sense, the observed secondary flow cells are closely linked to the conditional quasi-streamwise rotations found between conditionally averaged 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 singlephase 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 side-walls, , 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., 2019b).
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 and Nakagawa, 1993). However, it should be kept in mind that the experimentally studied bedforms are usually subject to a stabilizing effect of lateral side-walls 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 simulation.

Declaration of Interests
The authors report no conflict of interest. A long-time evolution of case S250 Figure 18 provides the space-time evolution of the streamwise-averaged sediment bed height fluctuation h b (z, t) and streamwise particle flux fluctuation q p,x (z, t) for case S250, respectively, for the entire observation interval T obs . As such, figure 3(a) represents a close-up of the initial 85T b of figure 18(a).

B.1 Quantification of bedform dimensions
In an abstract sense, we can interpret the evolution of sediment bedforms as wave-like deformations of the defined two-dimensional surface that represents the fluid-bed interface. In this context, it is natural to quantify the wall-parallel and wall-normal bedform dimensions as wave amplitude and wavelengths, respectively. To distinguish between the effects of streamwise and spanwise bedforms, we investigate the length scales for the streamwise-and spanwise-averaged interface h b x (z, t) and h b z (x, t) separately. As a measure for the ridge height, we choose a statistical approach based on the root mean square of the sediment bed height fluctuations (Kidanemariam and Uhlmann, 2017;Langlois and Valance, 2007), viz. Space-time evolution in case S250 of (a) the sediment-bed height fluctuation of the streamwise-averaged bed and (b) the streamwise-averaged particle flux fluctuation q p,x /q x,rms , with q x,rms (t) = ( q p,x q p,x z ) 1/2 . σ h,x is accordingly defined based on the fluctuations of the spanwise-averaged interface. Alternative approaches to determine bedform dimensions are reviewed in Coleman and Nikora (2011). The spanwise mean wavelength of the sediment bed, λ h,z , is defined below based on the instantaneous two-point correlation coefficient where δz denotes the spanwise separation between two positions. For a given time t, we identify the mean wavelength as twice the separation length for which ρ hh attains its first and global minimum, ρ hh (δz min , t) ≤ ρ hh (δz, t) ∀ δz ∈ [0, L z /2] λ h,z (t) = 2δz min .
B.2 Quantification of particle transport In order to transform the discretely spaced information of Lagrangian particle properties to the Eulerian observation framework, we apply a binning technique similar to that of Kidanemariam (2015), but generalized to the two-dimensional case. Here we choose to discretize the two periodic directions in bins of width ∆x bin = ∆z bin ≈ 1.5D spanning over the entire wall-normal box length L y , resulting in a number of N x,bin and N z,bin bins in the streamwise and spanwise direction, respectively. The local particle flux in the (i, k)-th bin (1 ≤ i ≤ N x,bin , 1 ≤ k ≤ N z,bin ) is then defined as the volumetric particle flow rate averaged over that bin, that is, the sum of the particle velocity of all particles centered in the bin times the particle volume divided by the bin base area ∆x bin ∆z bin , viz.
where q p,x , q p,z are the streamwise and spanwise components of the particle flux vector q p . V p = (π/6)D 3 denotes the volume of a single spherical particle and u (l) p , w p are the streamwise and spanwise component of the particle velocity vector u (l) p of particle l, respectively. The bin centers are located at x i = ((i − 1) + i)∆x bin /2 and z k = ((k − 1) + k)∆z bin /2. I (i,k) (l) (t) is an indicator function that is unity if the l-th particle is centered in the (i, k)-th bin, i.e.
Note that the particle flux is an integral measure and as such does not reflect the inhomogeneous distribution of particle transport intensity over the wall-normal direction. Let us therefore additionally introduce the streamwise particle flux density as with the wall-normal bin centers y j = ((j − 1) + j)∆y bin /2 (1 ≤ j ≤ N y,bin ) and the indicator function I In a similar way as in the previous definitions, we have subdivided the wall-normal box length into N y,bin bins with however smaller bin height ∆y bin ≈ 0.5D, this time spanning over the entire box length L x and width L z . By definition, the following relation holds (Chiodi et al., 2014;Lobkovsky et al., 2008): q p,x xzt = Ly y=0 φu p xzt dy.