On the scaling of the instability of a flat sediment bed with respect to ripple-like patterns

We investigate the formation of subaqueous transverse bedforms in turbulent open channel flow by means of direct numerical simulations with fully-resolved particles. The main goal of the present analysis is to address the question whether the initial pattern wavelength scales with the particle diameter or with the mean fluid height. A previous study (Kidanemariam and Uhlmann, J. Fluid Mech., vol. 818, 2017, pp. 716-743) has observed a lower bound for the most unstable pattern wavelength in the range 75-100 times the particle diameter, which was equivalent to 3-4 times the mean fluid height. In the current paper, we vary the streamwise box length in terms of the particle diameter and of the mean fluid height independently in order to distinguish between the two possible scaling relations. For the chosen parameter range, the obtained results clearly exhibit a scaling of the initial pattern wavelength with the particle diameter, with a lower bound around a streamwise extent of approximately 80 particle diameters. In longer domains, on the other hand, patterns are observed at initial wavelengths in the range 150-180 times the particle diameter, which is in good agreement with experimental measurements. Variations of the mean fluid height, on the other hand, seem to have no significant influence on the most unstable initial pattern wavelength. Furthermore, for the cases with the largest relative submergence, we observe spanwise and streamwise sediment waves of similar amplitude to evolve and superimpose, leading to three-dimensional sediment patterns.


Introduction
Sediment beds under streaming water bodies often exhibit characteristic bedforms which can significantly affect the transport properties of the flow. Among these various bedforms, transverse bedforms are usually classified as ripples, dunes or antidunes (Yalin, 1977). Ripples and dunes both exhibit a roughly triangular shape with a gentle slope on the upstream surface and a steeper downstream face, the latter approximately inclined at the angle of repose (Best, 2005). However, these two bedforms differ in their scaling properties: while ripples are small compared to the fluid height and thus their wavelength is believed to scale with the particle diameter, the height of dunes is large enough to distinctly modify the flow field over the entire flow depth. As a consequence, the amplitude of dunes should scale with the fluid height (Engelund and Fredsoe, 1982). In contrast to the remaining bedform types, antidunes exhibit a symmetric cross-section and they can, depending on the hydraulic conditions, travel in the upstream or the downstream direction. Antidunes can only form in free-surface flows, whereas ripples and dunes can form in configurations without a free surface as well, such as pipe flows (see e.g. Ouriemi et al., 2009) or closed-conduit flows (see e.g. Cardona Florez and Franklin, 2016;Coleman et al., 2003).
The origin of these transverse bedforms has been found to be an instability process in which an initially flat sediment bed deforms after an initial perturbation in the flow system leading to the formation of transverse patterns (Yalin, 1977). Inglis (1949) proposed that inhomogeneities of the sediment bed such as small particle agglomerations could act as initial perturbations that trigger subsequent bedform evolution. In a more recent paper, Coleman and Melville (1996) proposed an instability process originating in an isolated random pileup that grows with time and, once it has reached a critical height of approximately 3 − 4 times the particle diameter D, induces the generation of further pileups downstream. Nikora (2009, 2011) describe the formation of the initial bedforms as a two-stage formation process. In a first stage, random sediment patches with a typical length of 7 − 15D interact on an active (but still planar) sediment bed due to sediment transport events that are believed to be related to coherent turbulent structures. Once the height of these initial random patches exceeds some threshold height, the initial disturbance stabilizes by agglomerating sediment particles, with the consequence that further downstream regular patterns of sand-wavelets form. Venditti et al. (2005) observed different modes of bedform initiation in their experiments depending on the flow rate and, consequently, on the Reynolds and Froude number. At lower flow rates, local defects are seen to initiate a similar formation cycle as that described by Coleman and Melville (1996), whereas at higher flow rates, the authors observed patterns to spontaneously evolve over the entire bed, leading to a regular 'cross-hatch pattern'.
Unfortunately, up to the present date, accurate measurements of the initial bedform dimensions remain challenging. On the one hand, this is due to the very short time window in which the initial bedform evolution can be observed, and, on the other hand, due to the small height of the initial patterns of only a few particle diameters which makes them hard to detect. It is for this reason that most of the experimental studies in the last decades have focused on the description of fully-developped patterns in the equilibrium state (see e.g. Yalin, 1985), where the wavelength of the initial bedforms is typically much higher than in the initial phase (Langlois and Valance, 2007). On the other hand, only a small number of experimentalists were able to quantitatively describe the wavelength of the initial bedforms. In turbulent open channel flow, Coleman and Melville (1996) and Coleman and Nikora (2009) observed the initial wavelength to scale mainly with the sediment size and to be rather unaffected by the flow conditions. Similar observations were made by Coleman et al. (2003) in turbulent closed conduit flows as well as by Coleman and Eling (2000) in laminar open channel flows, the latter indicating in particular that the formation of initial bedforms is not restricted to the turbulent regime and the accompanying turbulent bursts, as earlier suggested by Raudkivi (1997). Similarly, Langlois and Valance (2007) and Cardona Florez and Franklin (2016) both report under turbulent channel flow conditions that the initial wavelength depends mainly on the particle diameter, while the influence of the shear velocity and the flow conditions is rather weak. Ouriemi et al. (2009), on the contrary, measured initial wavelengths of the order of the fluid height in pipe flows, while Franklin (2008) observed a dependence of the initial wavelength on both the particle diameter and the shear velocity in experimental measurements of turbulent closed conduit flows.
In theoretical studies based on linear stability analysis, the flow system is usually described by a simplified model for the driving flow (such as RANS models or potential flow solutions) combined with a sediment bed continuity equation for the evolution of the bed (see e.g. Charru, 2006;Kennedy, 1963Kennedy, , 1969. In order to close the system of equations, a formulation for the particle flux is required. In most stability analysis, it is assumed that the particle flux and the local bed shear stress are in phase, which allows then to express the particle flux as a function of the local bed shear stress (Charru et al., 2013). In recent studies, this assumption has been removed and additional relaxation equations are used to take into account a possible phase lag between both quantities (Charru, 2006). A linear stability analysis is then performed in order to determine regions of instability in the parameter space as well as the most amplified wavelength for a given flow configuration. During the past decades, a large number of stability analysis for different flow configurations including different stabilizing and destabilizing effects has been presented for turbulent (see e.g. Colombini, 2004;Colombini and Stocchino, 2011;Fourriere et al., 2010;Richards, 1980;Sumer and Bakioglu, 1984) as well as for laminar flows (see e.g. Charru and Hinch, 2006;Charru and Mouilleron-Arnould, 2002). A detailed overview of the different approaches used in linear stability analysis can be found in the reviews of Engelund and Fredsoe (1982), Seminara (2010) and Charru et al. (2013). Recently, Zgheib and Balachandar (2019) have presented a combined numericaltheoretical stability analysis, in which the bed shear stress for the linear stability analysis is computed by means of direct numerical simulations (DNS) in which the sediment bed is represented in a continuous and impermeable fashion.
However, despite an increasing complexity of the used models, the wavelengths predicted by most linear stability analysis still differs by more than one order of magnitude from the wavelengths observed in experiments (Langlois and Valance, 2007;Ouriemi et al., 2009). Furthermore, the observations concerning the scaling of the initial wavelengths differ markedly between the different studies. For instance, Fourriere et al. (2010) found a single region with unstable wavelengths independent of the fluid height, leading them to the conclusion that only ripples can form directly from a flat bed, while dunes develop by a ripple coarsening processes only. Colombini and Stocchino (2011), in contrast, observed two separate regions of instability with most amplified wavelengths of the size of the fluid height and of the particle diameter, respectively, which suggests that either ripples or dunes may form out of the same instability mechanism.
It should be kept in mind that the validity of linear stability investigations is limited to the very first instances of sediment bed evolution (Coleman and Melville, 1996), i.e. the theory is not able to predict the bedform dimensions correctly, once non-linear effects become dominant. In more recent studies, thus, weakly non-linear analysis have been presented (e.g. Colombini and Stocchino, 2008) to take into account also non-linear effects.
In recent years an alternative method for the investigation of sediment transport in its early stages has become available in form of DNS featuring fully-resolved particles, which allow to resolve all relevant flow scales even below the particle length scale (e.g. Derksen, 2015;Uhlmann, 2014a,b, 2017;Vowinckel et al., 2014Vowinckel et al., , 2017. In a set of numerical experiments, Kidanemariam and Uhlmann (2017) (henceforth denoted KU2017) have recently determined a lower bound for the minimal unstable wavelength, λ th , by reducing the streamwise domain length L x successively in a comparable concept as the minimal flow unit of Jiménez and Moin (1991). It was observed that below a threshold L x = λ th , the formation of transverse bedforms is effectively hindered and a perturbed bed remains stable as a consequence of the limited domain size, even though the conditions would otherwise allow the bed to become unstable. For the considered parameter values, KU2017 found λ th /D to be in the range 75−100, which is equivalent to a range of λ th /H f = 3 − 4 in their cases. Since, however, the relative fluid height H f /D has been kept constant over all simulations, it was not possible to further distinguish between the two alternative scaling relations.
The purpose of the present work is to investigate the scaling of the lower threshold for the minimal unstable wavelength. In order to be able to distinguish between the two alternative scalings, i.e. either with the particle diameter D or the fluid height H f , two new series of numerical experiments are performed in which the relative streamwise domain lengths L x /D and L x /H f are varied independently. In the subsequent analysis, we investigate the influence of both length scales on the stability of the subaqueous bedforms using the newly computed DNS data. The present manuscript is organized as follows. In § 2, we briefly describe the numerical method which we use for the simulation of subaqueous sediment transport in this work. An overview over the relevant physical and numerical parameters as well as the chosen flow configurations is given in § 3. Subsequently, in § 4, we shortly present two different ways of defining the fluid-bed interface, depending on whether the sediment bed is analysed in a spanwise-averaged framework or in its full extent. In § 5, we analyse the data obtained in the context of the two simulation series separately, focussing on the bedform geometry and its temporal evolution. A discussion of the results as well as a comparison with experimental and theoretical studies follows in § 6. We conclude the study with a summary of the main findings in § 7.

Numerical method
We use the same numerical method used in Uhlmann (2014a,b, 2017) to solve the coupled fluid-solid problem. For the simulation of the fluid phase, the incompressible Navier-Stokes equations are solved numerically in the entire computational domain using a second order finite difference scheme together with a fractional step algorithm on a uniform Cartesian grid. Time integration of the governing equations is done in a semi-implicit way, including a Crank-Nicholson scheme for the viscous terms and a low-storage three-step Runge-Kutta scheme for the non-linear terms. The immersed boundary formulation of Uhlmann (2005) is then used to couple the flow field with the solid phase: localized force terms are introduced into the Navier-Stokes equations which impose the no-slip condition at the interface between fluid and solid phase. The motion of the particles is obtained by time integration of the Newton-Euler equations for rigid-body motion. The driving force and torque comprises hydrodynamic and gravitational contributions as well as those resulting from particle-particle and particle-wall contact. Due to the fact that the characteristic time scale of particle collisions is typically several orders of magnitude smaller than those of the turbulent fluid motion, a sub-stepping method is used for the numerical time integration of the Newton-Euler equations (Kidanemariam and Uhlmann, 2014b). Momentum exchange due to particle collisions is computed using a soft-sphere discreteelement model (DEM). In the frame of the chosen approach, two particles are defined as 'being in contact', if the minimal distance between their surfaces ∆ falls below a force range ∆ c . In this case, a contact force and torque act on both particles which is defined as the sum of three individual contributions, i.e. an elastic normal force, a normal damping force and a tangential frictional force. The elastic normal force is a linear function of the penetration length δ c = ∆ c − ∆ with a constant stiffness coefficient k n . The normal damping force is a linear function of the normal component of the relative particle velocity between the two particles at the contact point, with a constant normal damping coefficient c n . Similarly, the tangential frictional force is defined as a linear function of the tangential component of the relative particle velocity at the contact point, with a constant tangential friction coefficient c t . Note that the magnitude of the tangential frictional force has an upper traction limit in the form of the Coulomb friction limit with a friction coefficient µ c . A detailed description of the collision model and extensive validation can be found in Kidanemariam and Uhlmann (2014b).
For each simulation, the model thus requires the choice of the four force parameters (k n ,c n ,c t ,µ c ) as well as the force range ∆ c . Introducing a dry restitution coefficient ε d , which is defined as the absolute value of the ratio between the normal components of the relative velocity before and after a dry collision, allows to relate the normal stiffness coefficient k n and the normal damping coefficient c n , and to formulate the model depending on the alternative set of force parameters (k n ,ε d ,c t ,µ c ). Due to the varying particle size and submerged weight in the present simulations, the set of parameters used in Kidanemariam and Uhlmann (2014a) and Kidanemariam and Uhlmann (2017) has been adapted to the respective cases. The force range ∆ c is set equal to the uniform grid spacing ∆x for all cases with a particle diameter D ≤ 15∆x and equal to 2∆x for all cases with larger particles. The stiffness parameter k n has a value between approximately 8400 and 17000 times the submerged weight of a single particle, divided by the particle diameter in the respective case. The parameters are chosen such that the maximum overlap δ c is within a few percent of ∆ c . The dry restitution coefficient is set to ε d = 0.3, which determines, together with the chosen value of k n , the constant normal damping coefficient c n . The constant tangential friction coefficient is set to the same value c t = c n . Finally, the Coulomb friction coefficient was fixed at µ c = 0.4 except for cases H2D052 and H6D102, in which a slightly higher value µ c = 0.5 was unintentionally used. In Appendix A we show, however, that this difference in the limiting Coulomb friction has only a minor influence on the eventually developed bedform and that it does not affect the stability or instability of the sediment bed.

Flow configuration and parameter values
In the course of the current study, we have performed seven new independent direct numerical simulations of bedform evolution over an erodible subaqueous sediment bed in a turbulent open channel. Two additional cases from KU2017 have been included in our analysis. The studied flow configuration is shown in figure 1. As can be seen, a Cartesian coordinate system is centered at the lower boundary of the open channel such that the x-, y-and z-direction are the streamwise, wall-normal and spanwise direction, respectively. Mean flow is directed in positive x-direction and gravity in the negative y-direction. In the streamwise and spanwise directions, periodic boundary conditions are imposed, whereas a no-slip-and a free slip-condition are imposed at the bottom and at the top plane of the channel, respectively.
As characteristic length scales, we define the mean fluid height H f and the mean sediment bed height H b as an average over both spatial directions (streamwise and spanwise) as well as over time, i.e. H f = h f zxt and H b = h b zxt , respectively. A precise definition of h f (x, z, t) and h b (x, z, t) will be given in section 4.1. Before the simulations with mobile particles are started, sediment bed and flow field have completed a start-up procedure, which is described in detail in Kidanemariam and Uhlmann (2014a). The names of the individual simulations are chosen according to their streamwise domain length in terms of the mean fluid height and in terms of the particle diameter. For instance, in case H4D052, the relative streamwise box length is L x /H f ≈ 4 and L x /D = 51.2, respectively (cf . table 1 and table 2).
In all cases, the flow is driven by a time-dependent streamwise pressure gradient, that is adjusted at each time step to ensure a constant flow rate q f . Therefore, the bulk Reynolds number can be computed a priori as where the friction velocity u τ is computed a posteriori from the streamwise and spanwise averaged total mean shear stress, which is composed of viscous stresses, turbulent  Table 1: Physical parameters of the simulations. Re b , Re τ and D + are the bulk Reynolds number, the friction Reynolds number and the particle Reynolds number, respectively. The density ratio ρ p /ρ f and the Galileo number Ga are imposed in each simulation, whereas the Shields number θ, the relative submergence H f /D and the relative sediment bed height H b /D are computed a posteriori (cf. table 2). The last column provides information about the source of the listed cases, dinstinguishing between simulations that have been computed in the course of the current study (present) and cases that are from Kidanemariam and Uhlmann (2017)(KU2017). It should be further mentioned that the physical parameters presented for case H4D102 1,2,3 have been averaged over the three individual simulations. A list of the physical parameters for each individual run can be found in KU2017.  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 grid with mesh width ∆x, ∆x + denoting the grid width in wall-units. N p is the total number of particles in the respective case. The time is scaled in bulk time units T obs is the total observation time of each simulation starting from the release of the moving particles. Time dependent physical and numerical parameters in tables 1 and 2 (Re τ , D + , H f , H b , θ, ∆x + ) are computed as an average over a final time interval T s obs .
Reynolds stresses and the stresses resulting from the fluid-particle interaction. Due to the absence of a horizontal bottom wall, we evaluate the mean friction velocity u τ at the location of the mean fluid-bed interface y = H b , which can be interpreted as a virtual wall. For a more detailed description of the determination of the shear stress distribution, the reader is referred to KU2017. From dimensional analysis, it can be concluded that by adding sediment to a turbulent flow, the parameter space increases and in total, four non-dimensional numbers are required to fully describe a given system. In addition to Re b , we choose the density ratio ρ p /ρ f , the non-dimensional length scale H f /D as well as the Galileo number Ga = u g D/ν, which expresses the ratio between gravity and viscous forces, where the gravitational velocity scale is u g = (ρ p /ρ f − 1) |g|D.
In the present simulations, we set the density ratio at ρ p /ρ f = 2.5 which is comparable to the values reported for glass in water. To allow pattern formation, the non-dimensional boundary shear stress, expressed by the Shields number θ = u 2 τ /u 2 g = (D + /Ga) 2 , has to be larger than the critical value for incipient sediment motion. In turbulent flows, the critical Shields number has been observed to be in a range θ c = 0.03 − 0.05, slightly depending on the Galileo number (Franklin and Charru, 2011;Soulsby et al., 1997;Wong and Parker, 2006).
The non-dimensional mean fluid height H f /D is varied in the different simulations to elucidate the relevant length scales that dominate the scaling of subaqueous bedforms by either increasing the particle diameter or the mean fluid height while keeping the same dimensions for the domain. As a consequence, the number of fully-resolved particles lies in a range between O(10 3 ) in the largest particle case and up to O(10 5 ) in the case with the smallest particles.
Note that in all cases, the dimensions of the computational domain are sufficiently large to allow self-sustained turbulence. In particular, the case with the shortest streamwise and spanwise dimensions (scaled in viscous lengths) is case H2D102 2 with L + x ≈ 726 and L + z ≈ 555. By comparison, Jiménez and Moin (1991) report the dimensions of their minimal flow-unit as L + x ≈ 250 − 350 and L + x ≈ 100, respectively. In section 5, we will analyse the two simulation series in which either the particle diameter or the mean fluid height have been varied separately. It should be noted that we have performed two additional simulations H2D052 and H4D052, which do not fit into either of these two series. Therefore, we discuss the results obtained in these latter cases in section 6 only.

Extraction of bedform dimensions
For the following definition of the fluid-bed interface, we will consider the domain as composed of two distinct regions, that are, a lower particle-dominated region, hereafter termed as the sediment bed, and the overlaying fluid-dominated region. Hence, the wallnormal dimension of the channel height L y can be written as the sum of the instantaneous local height of the fluid phase h f (x, z, t) and that of the sediment bed h b (x, z, t), viz.
In the chosen Cartesian coordinate system, the wall-normal location of the fluid-bed interface is identical to the value of the function h b (x, z, t). In the following, we will present two different approaches to define the instantaneous fluid-bed interface. The first method defines the fluid-bed interface through a spanwise averaging, whereas the second extracts the sediment bed as a two-dimensional surface.

Definition and analysis of the spanwise-averaged fluid-bed interface
Assuming statistical homogeneity of the observed system, equation (1) can be averaged in the spanwise direction as The spanwise averaged sediment bed height h b z is determined depending on a threshold for the solid volume fraction. Since the method is presented in detail in our previous works Uhlmann, 2014a, 2017), we restrict ourself to a short summary of the most important points. First, a solid-phase indicator function φ p (x, t) is defined, which attains the value of unity for Eulerian grid points being located inside the particle domain Ω p and zero elsewhere. Second, the spanwise-averaged sediment bed height h b z is defined as the wall-normal location, at which the spanwise averaged solid indicator function φ p z attains a threshold of φ p thresh z = 0.1 (Kidanemariam and Uhlmann, 2014b) The spanwise-averaged fluid height h f z is computed using relation (2). Eventually, we perform another decomposition in the streamwise direction, similar to the spanwise decomposition used in equation (2). The resulting expression divides the spanwise-averaged interface in an instantaneous mean height h b zx (t) and a fluctuation h b z (x, t) with respect to the former. In the following, the fluctuation will form the basis for the analysis of the bedform evolution. The size and shape of twodimensional transverse patterns is usually quantified by the pattern length and height. Over the last decades, a variety of different approaches has been developed to define these length scales. An overview over some of these methods is presented in Coleman and Nikora (2011). In the current work, we choose a statistical definition for the pattern height as well as for the wavelength of the bedforms (Langlois and Valance, 2007). As a measure for the pattern height, we use the root mean square of the sediment bed height fluctuation In order to determine the average pattern wavelength, let us define the instantaneous two-point correlation coefficient of the sediment bed height fluctuation where δx expresses the streamwise separation length between two points. With increasing δx, the evolution of R h (δx, t) at a given time t behaves similar to a damped oscillation curve around a zero mean. The global (and first) minimum of the curve occurs at a separation δx min , i.e.
The mean wavelength is finally obtained as twice this separation length, i.e. λ h (t) = 2δx min .

Definition and analysis of the two-dimensional fluid-bed interface
The spanwise-averaged definition and analysis, as described in the previous section, is restricted to a fluid-bed interface that is statistically homogeneous in the spanwise direction. However, in the general case, this procedure may be not adequate to describe possible three-dimensional sediment patterning.
Here, we use a criterion based on determining the top-layer particles of the fluid-bed interface. In classical morphodynamics, the sediment bed is distinguished from bedload and suspended load (Bagnold, 1956;van Rijn, 1984). Following this classification, the developed algorithm sorts out the latter two categories, leaving only the particles that are part of the sediment bed itself. First, single suspended particles are detected based on the wall-normal collision force component F c,y . While particles inside the bed are exposed to a permanent wall-normal contact force of the order of their submerged weight F w = (ρ p − ρ f ) πD 3 /6 |g|, single suspended particles are typically not in contact with surrounding particles and, accordingly, only a small contact force acts on them. Thus, particles that are only exposed to a negligible wall-normal contact force (|F c,y |/F w < 10 −5 ) will not be considered as 'bed particles'.
The remaining particles include both the actual sediment bed and the bedload layer. Particles inside the latter are in motion close to the bed, such that they lose the contact to the bed only for short time intervals (Yalin, 1977). In order to separate the bedload transport particles from the bed particles, a criterion based upon the particle speed is employed, eliminating all particles which exceed a threshold value. The threshold value can be chosen by defining a non-dimensional particle Shields number θ p = (|u p |/u g ) 2 with the norm of the particle velocity |u p | and the gravitational velocity scale u g . The sediment bed particles are then found as the set of all particles for which θ p ≤ θ c . Here we choose as the critical value θ c = 0.05 (Wong and Parker, 2006). Note that the particle velocity criterion also eliminates suspended particle pairs which might instantaneously experience collision forces above the chosen threshold of the wall-normal contact force.
The fluid-bed interface is then defined through a three-dimensional α-shape (Edelsbrunner and Mücke, 1994) enclosing the set of sediment bed particles. This enclosing surface can be thought of as a generalization of a convex hull, with an imposed radius α (here taken as 1.1 times the particle diameter) that defines the length scale above which non-convexity is allowed. The fluid-bed interface is then made up of those triangular faces of the α-shape which have an outward-pointing face-normal with a positive y-component. The data consists of a function h(x, z) which is sampled at the projection of the vertex points of the surface triangulation upon the (x, z)-plane. In a final step, this function h(x, z) is interpolated to a uniform Eulerian grid with a mesh width equal to one particle diameter and the result is smoothed using a two-dimensional box-filter with a width of 5 particle diameters.
For the analysis of this two-dimensional interface, let us extend our definition of the root mean square of the sediment bed height fluctuation, as defined in equation (5), to the multidimensional case in the following form:

Variation of the particle diameter
In a first series, the influence of the particle diameter D on the minimal unstable wavelength λ th at a given parameter point is investigated. To this end, a series of four simulations with different particle size is performed, in which the particle diameter D is successively increased from D/H f ≈ 0.04 in case H6D154 to D/H f ≈ 0.13 in case H6D052. It should be stressed that in all four simulations, the streamwise domain length L x is almost constant in the range 6.04 − 6.52H f and thus clearly above the critical value observed in KU2017. Note that case H6D154 with the smallest particle diameter is identical to case H6 presented in KU2017, in which a single bedform has been observed to evolve and to eventually reach a quasi-equilibrium state. Figure 2 shows instantaneous snapshots of the sediment bed in the different cases after a simulation period of approximately 300 bulk time units. It is seen that in case H6D052, no transverse pattern has formed. Instead, the bed has remained essentially flat and eroded sediment seems to distribute over the entire channel length and width, without accumulating in specific locations. In particular, case H6D052 does not exhibit a similar regular pattern of streamwise aligned alternating ridges and troughs as those KU2017 found in the stable sediment bed of their case H3. On the other hand, the remaining new cases with a domain length L x ≥ 76.8D are unstable, each featuring one single transverse bedform with an initial wavelength of the order of the streamwise box length. This indicates that at the given parameter point, there is indeed a lower limit for λ th , which depends on the particle diameter and which is found in the range 51.2−76.8D. As can already be inferred from the top-view visualizations of the sediment bed (cf. figure 2), the shape of those bedforms which do emerge differs, in particular in the vicinity of the threshold.
In order to provide a quantitative description of these patterns, we will now discuss various geometrical measures. First, we will focus on the pattern height evolution by studying the root mean square sediment bed height fluctuation σ h , which can be seen as a measure for the inhomogeneity of the sediment bed height Uhlmann, 2014a, 2017;Zgheib et al., 2018). In a case, in which no transverse bedforms evolve, σ h will show some fluctuations due to random uncorrelated bed undulations, but it will remain bounded by a small value in the course of the simulation. A sediment bed that shows such evolution will be classified as stable. In unstable systems, on the other hand, σ h will grow more or less monotonically during the initial phase of the simulation and in particular, it will exceed the aforementioned threshold that bounds the stable cases. Figure 3 shows the time evolution of σ h . After a short initial transient of a few bulk time units during which the chaotic particle motion leads to small finite values of σ h , we observe that in all cases with L x ≥ 76.8D, σ h increases with time starting at t ≈ 20T b , indicating that the chosen relative box length is sufficiently long to cover at least one unstable mode and thus to allow the bed to evolve transverse patterns. On the contrary, σ h in case H6D052 does not exhibit any substantial growth period. Instead, it remains around a value of σ h ≈ 0.3D, only featuring fluctuations of small amplitude, bounded by an upper value of approximately 0.47D. This value is in good agreement with the results of Coleman and Nikora (2009), who observed "mobile sediments on planar but active bed" for σ h in the range 0.40 − 0.50D. Unstable bedforms are usually observed to run through different phases of bedform evolution (KU2017). During an initial growth period, the bed increases exponentially as predicted by linear stability theory. After some time, non-linear contributions become relevant and let the bed height tend to its quasi-steady equilibrium value (Charru et al., 2013). KU2017 observed that all their unstable cases exhibited an exponential growth at a very similar growth rate, apparently independent of the streamwise box length of the respective case. They found an exponential function with amplitude A = 0.0668 and growth rate B = 0.0140 to best fit the initial increase of σ h with time over an interval of approximately 150 bulk time units. In contrast, the temporal evolution in the first instances of the current unstable cases strongly varies from case to case. While the initial growth of case H6D154 follows the above described exponential function of KU2017, the remaining two cases H6D102 and H6D077 grow more slowly. In the subsequent phase, case H6D154 and H6D102 attain what could be called as asymptotic state after approximately 280 and 350 bulk time units, respectively. Case H6D077 exhibits stronger oscillations with higher amplitude compared to the previous two cases, a behaviour similar to that of case H4D102 3 from KU2017. Averaging σ h over the final time interval T s obs leads to values of σ h t /D ≈ 2.08, σ h t /D ≈ 0.89 and σ h t /D ≈ 0.67 for cases H6D154, H6D102 and H6D077, respectively. It is remarkable that the mean values in the final interval differ by more than a factor of two between cases H6D154 and H6D102, which possess a comparable relative domain length of L x /H f ≈ 6. On the other hand, cases H6D102 and H4D102 1,2,3 attain very similar values of σ h t /D ≈ 0.89 and σ h t /D ≈ 0.96, respectively, although case H4D102 1,2,3 has a smaller relative box length L x /H f ≈ 4. This indicates that the attained mean value of σ h in the final interval mainly depends on the chosen L x /D ratio, whereas it is not very sensitive to a variation of the mean fluid height H f . Figure 4 shows the time evolution of the mean pattern wavelength λ h . It can be observed that the chosen mean wavelength for all cases with L x ≥ 76.8D settle, after some fluctuation in the first 100−200 bulk time units, at the maximum possible wavelength, i.e. λ h = L x , and maintains this value until the end of the simulation. KU2017 observed a similar evolution of the mean wavelength settling at λ h = L x for their shorter cases up to a box length L x = 179.2D. The current observations further support their findings that these systems cannot freely choose their initial wavelength, but that they are constrained by the limitation of the streamwise domain size. As a consequence, the system chooses the maximum possible unstable wavelength as the dominant one, which, however, is not necessary the same as the one it would choose in a system without spatial limitations (cf. the very long box with L x /H f = 48 simulated by KU2017). In contrast to the observed unstable cases, λ h in case H6D052 jumps between the first three harmonics λ 1 = L x , λ 2 = L x /2 and λ 3 = L x /3 for the entire observation interval. This behaviour is caused by random uncorrelated bed disturbances and it indicates that the system is not able to develop a pattern at a finite wavelength.
In the following, we will investigate the fully-developed spanwise-averaged profile of the sediment pattern along the streamwise direction. To this end, we compute the phaseaveraged bed heightH b (x) as defined in KU2017. For this purpose we use a constant pattern migration velocity u D which is determined from a linear fit of the space-time correlation. Phase-averaging is performed over the final time interval T s obs (as indicated in table 2). We further introduce the aspect ratio AR and the degree of asymmetry LR of the pattern as (KU2017) where from the crest to the neighbouring downstream trough and λ h is the mean wavelength. Figure 5 shows the mean pattern profile averaged over the final time interval of all simulations as well as the attained aspect ratio of the bedforms. It can be seen, that the profiles of cases with L x /D ≥ 153.6 almost collapse upon a self-similar profile exhibiting a strong asymmetry (LR ≈ 0.28) and an aspect ratio of AR ≈ 0.037 (KU2017). Some experimental studies report higher values for the steady-state aspect ratio, i.e. AR ≈ 0.045 (Fourriere et al., 2010), AR ≈ 0.05  and AR ≈ 0.067 (Charru et al., 2013). KU2017 explain their lower value with the fact, that in their cases, "the 'natural' steady-state regime is not yet reached" and thus, a larger aspect ratio might be attained in a later phase, whereas the experimental data reflects the long-time state of the system which in those cases has evolved for at least one order of magnitude longer times. On the other hand, the pattern profiles in cases H6D102 and H4D102 1,2,3 with L x = 102.4D are characterized by a smaller aspect ratio of AR ≈ 0.024 − 0.028 and a more symmetric shape (LR ≈ 0.33 − 0.036). Important to note is, however, that all four profiles with L x = 102.4D are again nearly self-similar, despite the fact that L x /H f varies by up to a factor of 1.5. This is another indication that the pattern shape is mainly controlled by the particle diameter D, but almost unaffected by the mean fluid height H f . Further reducing the relative domain length to a value L x /D = 76.8 in case H6D077 leads to an even more symmetric profile (LR ≈ 0.41), whereas the aspect ratio still attains a value AR ≈ 0.024, comparable to the cases with L x = 102.4D. An additional observation in case H6D077 is that, in contrast to the previous cases, the curve representing the pattern profile is not smooth, but exhibits wavy disturbances. This is a direct consequence of the lower amount of particles in case H6D077, which leads to stronger disturbances of the spanwise-averaged pattern profile.

Variation of the mean fluid height
In order to study the dependence of the minimal unstable wavelength λ th on the mean fluid height H f , a second set of three simulations will be analysed in the following, including case H4D102 1,2,3 of KU2017 as well as two new simulations H2D102 1 and H2D102 2 . All three cases have the same initial sediment bed configuration and the streamwise and spanwise box dimensions match. However, in case H2D102 1 and H2D102 2 , H f has been increased by a factor of two compared to case H4D102 1,2,3 . Consequently, the bulk and friction Reynolds number are higher in the former two cases (cf. table 1). In cases H2D102 1 and H2D102 2 , two different values have been chosen for the Shields number, i.e. θ = 0.18 in the former and θ = 0.13 in the latter, by adjusting the Galileo number, while the remaining parameters are identical. The dimensions of the channel configuration in the new simulations have been chosen in such a way, that their relative domain length (L x /H f ≈ 2) lies below the lower bound for the most unstable wavelength as initially reported by KU2017. As a consequence, none of the two cases should allow the sediment bed to become unstable and transverse patterns to evolve, if the minimal unstable wavelength scales with the mean fluid height. Note, however, that these authors did not vary L x /D and L x /H f independently, and, therefore, could not distinguish between the two alternative scaling relations. Here, we are now in a position to do this. Figure 6 and figure 7 show selected instantaneous snapshots of the particle bed evolution for the new cases H2D102 1 and H2D102 2 , respectively. The sediment bed surface of case H2D102 1 is deformed showing streamwise and spanwise elongated crestlines. In some phases, the system alters between patterns at either of the two orientations (cf. figure 6(b) and 6(c)), whereas in figure 6(a), for instance, streamwise and spanwise oriented crestlines appear at the same time. Similarly, the sediment bed of case H2D102 2 exhibits one crest line parallel and one perpendicular to the mean flow direction in figure 7(a). In contrast, a later snapshot in figure 7(b) shows a three-dimensional bedform with a horseshoe like shape with horns on both sides pointing downstream. This pattern resembles in its geometry a barchan dune (Franklin and Charru, 2011). In the last snapshot provided in figure 7(c) a diagonal crest line, crossing the whole extent of the x-z-plane from the lower left to the upper right corner, has evolved. These observations represent a marked difference between the systems with larger clear fluid height (H f /D ≈ 50) and those with smaller relative submergences (H f /D ≈ 25). In the latter case (simulations H4D102 1,2,3 of KU2017 -images not shown), regular streamwise-aligned ridges are either displaced by the higher, dominating transverse bedforms or do exist superimposed to the former ones.
An important consequence of the observed behaviour in cases H2D102 1 and H2D102 2 is that the bed height is clearly two-dimensional. This means that the spanwise-averaged sediment bed and flow field will not correctly describe the physical processes that lead to the formation of these type of bedforms. In the remainder of the current work, we will therefore analyse the sediment bed in its entire streamwise and spanwise dimension and omit the spanwise averaging (cf. the definition in section 4.2). Figure 8 shows the time evolution of the root mean square of the two-dimensional sediment bed height fluctuation σ 2D h . All cases exhibit an exponential growth interval, but with clearly different growth rates, as indicated by the fitted exponential curves presented in figure 8. It is important to keep in mind that in contrast to the spanwise-averaged case, the two-dimensional measure σ 2D h takes into account streamwise and spanwise variations of the fluid-bed interface, which means that a change in σ 2D h with time can be a sign of evolving ridges or transverse patterns likewise. Indeed, the detailed analysis of a two-dimensional Fourier decomposition of the bed height perturbation which is presented below shows that the amplitude of variations in both directions is of similar order. In order to determine the scaling of the initial growth rate of the sediment bed height fluctuation, however, one would require a larger database and a larger number of individual runs at each parameter point to allow for ensemble averaging similar to case H4D102 1,2,3 . Nevertheless, it would be of high interest to elucidate the influence of the relevant parameters on the initial growth rate of σ 2D h such as the Reynolds number Re b (Re τ , respectively) and the relative submergence H f /D in a future study. After roughly  Figure 8: (a) Time evolution of the two dimensional root mean square of the bedform amplitude normalized by the particle diameter σ 2D h /D. Time is scaled in bulk time units T b . Cases H4D102 1,2,3 ( ), H2D102 1 ( ), H2D102 2 ( ). The data of case H4D102 1,2,3 is presented as ensemble average over the three simulations ( , thick line). Additionaly, the individual evolution of run H4D102 3 is presented ( , thin line). The dashed-dotted lines represent exponential curves of the form σ h /D = A exp(Bt/T b ) which have been found to best fit the evolution of σ 2D h during the initial growth period of the respective case. The coefficients are as follows: H4D102 1,2,3 : A = 0.2215, B = 0.0070; H2D102 1 : A = 0.2441, B = 0.0266; H2D102 2 : A = 0.1926, B = 0.0130. Note that in case H4D102 1,2,3 , the exponential curve best fits the ensemble average over the three runs. (b) Same data as (a), but represented in semi-logarithmical scale. 300 bulk time units, all cases eventually reach what could be called an asymptotic state. We observe that the pattern amplitude attains fully-developed averaged values of similar order comparable to one particle diameter for all cases, which suggests that at the given parameter point shown in figure 8, the averaged value of σ 2D h in the final interval does not strongly depend on the value of the mean fluid height H f . The analysis of σ 2D h has shown that the sediment bed becomes unstable in all three cases. However, as σ 2D h is an integral measure for the bed evolution, it cannot provide further information about the evolution and interaction of single unstable modes. In particular, it does not allow to determine the role of streamwise and transverse modes separately. For this reason, the evolution of the bed height perturbation will be further analysed in Fourier space for the different streamwise and spanwise wavenumbers. To this end, we compute the single-sided amplitude spectraÂ (k,l) (for the physically relevant non-negative wavenumbers k, l ≥ 0) as twice the absolute value of the coefficient h b(k,l) (κ 1 k , κ 3 l , t), which is the Discrete Fourier Transform (DFT) of the sediment bed height perturbation in the physical space, h b (x, z, t). In the above expression, κ d k is the wavenumber of the k-th mode in spatial direction d (where d = 1, 3 corresponds to the x-and z-direction, respectively). Figure 9 shows the time evolution of the single-sided amplitude spectra for the most dominant modes in case H4D102 3 , H2D102 1 and H2D102 2 , respectively. In case t/(H f /u b ) Figure 9: Time evolution of the single-sided amplitude spectra for the most dominant modes normalized with the particle diameterÂ (k,l) /D for cases (a) H4D102 3 , (b) H2D102 1 and (c) H2D102 2 . Note that only those modes are defined as dominant, that exceed a value ofÂ (k,l) = 0.30D during the simulation. The notationÂ (k,l) indicates that the amplitude corresponds to the mode with streamwise and spanwise wavenumber κ 1 k and κ 3 l , respectively. Colouring of the individual dominant modes is as follows: , A (1,0) ; ,Â (2,0) ; ,Â (0,1) ; ,Â (0,2) ; ,Â (0,3) ; ,Â (0,4) ; ,Â (1,1) ; ,Â (1,2) ; , H4D102 3 , the wave with modes (1, 0) (first harmonic in the streamwise, constant in the spanwise direction) attains the highest amplitudes and thus clearly dominates the spectra. In good agreement with the evolution of σ 2D h , this mode first increases exponentially during the initial 200 − 250 bulk time units, before it settles at values which fluctuate somewhat around an asymptotic mean. Note that only the first two harmonics of the pure streamwise waves λ (1,0) and λ (2,0) exceed a value ofÂ (k,l) = 0.30D during the observation interval, highlighting that in these cases with a comparably short domain, a very sparse range of wavenumbers is amplified, and is then available to form the sediment pattern. On the other hand, waves of the first three pure spanwise harmonics λ (0,1...3) are found to be significant. In the first approximately 100 − 150 bulk time units, these pure spanwise modes dominate the amplitude spectra, which reflects the formation of initial streamwise aligned ridges that has been observed by KU2017. In the subsequent quasi-steady phase of the bedform, pure spanwise modes are still present with finite but smaller amplitudes compared toÂ (1,0) , reaching maximum valuesÂ (0,1) = 0.80D. This is in good agreement with the observation of the aforementioned authors, that streamwise elongated patterns are still visible on the upstream face of the transverse bedforms, once these latter have reached a quasi-steady state. In addition, several diagonal waves including modes larger than zero in both directions evolve, but they do not reach amplitudes higher thanÂ = 0.5D. The time evolution of the amplitude spectra in the new simulations H2D102 1 and H2D102 2 (cf. figures 9(b),(c), respectively), on the other hand, is dominated by the first pure spanwise oscillating wave mode (0, 1) during almost the whole observation interval, attaining maximum values ofÂ (0,1) /D above unity. It is worth to note that this latter is the only pure streamwise wave which is significant. In general, the amplitude of the single modes exhibit higher fluctuations than those in case H4D102 3 , leading to differences of approximately 1D between the highest and the lowest attained values of a single mode. During the entire simulation, none of the modes reach a plateau regime which would indicate a quasi-steady state of the bed. These findings match our observations, that both systems H2D102 1 and H2D102 2 show a sequence of different alternating bedform configurations composed of a number of streamwise and spanwise unstable modes without eventually reaching a quasi-steady state. In particular, transverse sediment bed waves evolve and, even though they are observed to be of smaller amplitude than their spanwise directed counterparts, they are involved in the formation of three dimensional patterns. Eventually, this reveals that, indeed, a sediment bed can become unstable for a domain length L x /H f < 3, indicating that the lower bound for the most-unstable wavelength as reported by KU2017 does not scale with the mean fluid height H f . The three-dimensional pattern evolution, however, requires further investigation.

Discussion
Our observations in the previous section suggest that at the investigated parameter point, the minimal unstable wavelength depends on the particle diameter, whereas it seems to be rather unaffected by variations of the mean fluid height. In the following, we will compare our results with measurements as well as proposed scaling relations from experimental and theoretical studies.
In figure 10, we present an overview of initial mean pattern wavelengths from our current and previous direct numerical simulation studies together with values determined in laboratory experiments. The observed wavelengths are shown as functions of the particle diameter D and of the mean fluid height H f , respectively. Note that for stable sediment beds (indicated by open symbols in figure 10), the definition of an initial pattern wavelength is meaningless. Here, we give instead the streamwise domain length L x , which indicates the maximum possible wavelength that did not evolve during the simulation. It should be noted that for stable systems as well as for cases H2D102 1 and H2D102 2 , in which three-dimensional patterns evolve, the streamwise boxlength is given instead of the spanwise-averaged mean wavelength (cf. the detailed discussion in the text). Wavelengths measured in experiments are presented with the following symbols: Coleman et al. (2003) (closed-conduit, ), Langlois and Valance (2007) (channel, ), Cardona Florez and Franklin (2016) (closed-conduit, for the first and for the last measured wavelength). Note that in the experimental studies, no free surface is present. Thus, the experimentally determined wavelengths are normalized with the half mean fluid height. Based on the results of our analysis, we expect the minimal unstable wavelength in the vertical grey region around λ/D = 80.
In cases H2D102 1 and H2D102 2 , on the other hand, we have observed the sediment bed to evolve both in the streamwise and spanwise direction. In these two cases we equally associate the value of the domain length L x with the most unstable wavelength, since the study of the single-sided amplitude spectra of the sediment bed height fluctuation shows that the most dominant streamwise mode is the one with indices (1, 0), i.e. the first streamwise harmonic which is constant in the spanwise direction. It should be further mentioned that with cases H2D052 and H4D052, two additional new cases have been added to the parameter plane, which were not part of the analysis in the previous section (cf . table 1 and table 2 for the physical and numerical parameters). In both cases, however, the bed is observed to remain stable, i.e. no transverse patterns form during the simulation. The experimental data shown in figure 10 comes from measurements in channel flow by Langlois and Valance (2007) as well as in closed-conduit flows by Coleman et al. (2003) andCardona Florez andFranklin (2016).
In our direct numerical simulations, we have observed transverse pattern formation only in cases with a relative streamwise domain length L x /D ≥ 76.8, while the sediment bed remained stable for all simulations below this limit, irrespective of the fluid height. At a domain length L x /D = 76.8, we have observed case H6D077 to become unstable, but the shape and asymmetry of the evolving pattern is seen to substantially deviate from the observed values in sufficiently long domains such as the cases of KU2017 with a box length of L x /D = O(10 3 ). A possible reason for this observation is that in H6D077, the system can only choose between 8 harmonics with a wavelength higher than O(D). The interaction of this small number of discrete modes may not lead to the same interface that would form in a sufficiently long domain. Interestingly, KU2017 report a case with the same relative domain length L x /D = 76.8 in which the sediment bed remained flat. We therefore expect that a box length L x /D = 76.8 is in the vicinity of the sought threshold wavelength, such that at this parameter point, already small differences in the configuration might tip the system in either way. If the domain length is even further decreased (as in cases H2D052, H4D052 and H6D052) the bed remains flat.
The existence of such a lower threshold for the most amplified wavelength originates in the interaction of several stabilizing and destabilizing effects in the flow system, which allow only a certain range of wavelengths to evolve. In order to show that this is indeed the case and that the absence of patterns in cases H2D052, H4D052 and H6D052 is not due to purely geometrical constraints, we have estimated the height of a pattern with wavelength λ/D = L x /D = 51.2 and a downstream face inclined at the angle of repose, which represents the highest possible pattern for a given domain length. On the other hand, a pattern at the same wavelength in the asymptotic state (i.e. with similar aspect ratio and degree of asymmetry as observed in the cases with L x /D ≥ 150D) would attain a height which is only half this maximum possible height. This observation thus verifies that all domains which are investigated in the current study are sufficiently long to accommodate, in principle, a fully-developped pattern without exceeding the limitations given in form of the angle of repose.
In the considered experimental studies, the shortest measured wavelength λ/D ≈ 105 is clearly above the critical wavelength determined in the present work. In our simulations, the shortest box which allowed a single bedform to reach an asymptotic state with an asymmetric triangular shape was observed for a very similar box length L x /D = 102.4. Nevertheless, the aspect ratio is still smaller than in the longer cases H6, H7 and H12 of KU2017, which might be an indication that the patterns did not yet reach their final state and would further evolve if their growth was not hindered by the limited domain size. Indeed, the majority of the experimental data points concentrates in a range of λ/D = 120 − 250 (λ/H f = 0.5 − 4.0). In good agreement, all numerical simulations with 150 ≤ L x /D ≤ 300 develop patterns with a mean wavelength in the range 150 ≤ λ/D ≤ 180 and a self-similar profile. These observations further support the suggestion of KU2017, that the bedforms finally settle at a wavelength of approximately 150 ≤ λ/D ≤ 180. Similarly, Nikora (2009, 2011) report a lower bound for the preferred initial wavelength of λ/D = 130. Eventually, the occurrence of patterns with a wavelength λ/H f < 3 in numerical simulations and experiments which can be seen in figure 10 excludes a possible scaling of the threshold with the mean fluid height which KU2017 could not exclude based on their limited data. This observation is further supported by the results of the analysis in the preceding section: cases with identical relative streamwise domain length L x /D but varying L x /H f ratio attain similar final values for all studied geometrical parameters and spanwise-averaged bedform profiles.
Compared to the wavelengths observed in experiments, the initial wavelengths predicted by most linear stability analysis are smaller by at least one order of magnitude (Langlois and Valance, 2007;Ouriemi et al., 2009). However, in some recent studies, the prediction of the most unstable wavelength could be improved taking additional physical effects into account. For instance, initial unstable pattern wavelengths similar to those determined in experiments were found after introducing an additional stabilizing effect in form of a phase-lag between the boundary shear stress and the particle flow rate into the linear stability analysis by several authors (Charru, 2006;Charru and Hinch, 2006;Fourriere et al., 2010). This phase shift, which is usually termed as characteristic saturation length L sat , has been reported to be a function of the particle diameter, attaining values of the order of 10 times the particle diameter for pure bedload transport (Claudin et al., 2011). For the special case of sand grains, Fourriere et al. (2010) present saturation lengths between 7 and 15 grain diameters. In their subsequent stability analysis, these authors predict for sufficiently large fluid heights most amplified wavelengths in a range of λ/L sat = 15 − 20, which is equivalent to a range λ/D = 105 − 300, when assuming that L sat = 7 − 15. This range of most amplified wavelengths overlaps with the values found in simulations and the mentioned experimental studies.
The findings of Colombini and Stocchino (2011), on the other hand, differ from our observations and the results of Fourriere et al. (2010). For the range of Galileo numbers considered in the present study, the authors report only one unstable region that depends on the fluid height exclusively. In contrast, unstable wavelengths scaling with the particle diameter, which could be related to a ripple instability, appear in their model for Ga ≤ 14 only. It should be however noted that the model of Colombini and Stocchino (2011) strongly depends on the roughness height and, consequently, on the relative submergence H f /D, which is in their case at least one order of magnitude higher than in our numerical simulations. Therefore, a direct comparison of the predicted critical wavelengths with our results is not possible.

Conclusion
In the current study, we have investigated the initial stage of subaqueous pattern formation in a turbulent open channel flow by means of direct numerical simulations with fully-resolved particles. The main contribution is to address the question of scaling of the initial bedform wavelength with either the particle diameter D or with the main fluid height H f . Both scaling relations have been proposed by different authors in the past decades, but up to the present day, there is no clear consensus about the correct scaling length. In a recent paper, Kidanemariam and Uhlmann (2017) have observed a lower bound for the most unstable wavelength λ th in a range 75 − 100D (equivalent to 3 − 4H f , respectively). However, it was not possible in that study to further tackle the problem of the correct scaling, since the relative submergence H f /D was constant in all of their simulations. In the present work, we have therefore performed and analysed two sets of simulations, in which we have varied the relative streamwise boxlengths L x /D and L x /H f independently, which allowed us to investigate the influence of both length scales on the initial pattern wavelength exclusively. Consequently, the relative submergence has been varied in the range H f /D = 7.86 − 50.59.
Our findings imply a scaling of the initial wavelength with the particle diameter, while it seems to be rather unaffected by variations of the mean fluid height. In particular, a lower bound for the most unstable wavelength has been found around a streamwise domain length of approximately L x /D = 80. In cases with a shorter relative box length L x /D, transverse pattern formation was effectively suppressed, indicating that the domain length is not sufficient to accommodate the minimal unstable wavelength. However, the observed pattern with a wavelength close to the proposed threshold exhibit remarkable differences compared to patterns in sufficiently long domains (cf. for instance the cases with L x /D = O(10 3 ) of Kidanemariam and Uhlmann (2017)), i.e. in very marginal boxes it does not reach an asymptotic state during the observation time and shows a more symmetric spanwise-averaged profile than in the longer cases.
The predicted range of initial pattern wavelengths is in good agreement with values measured in laboratory channel and closed-conduits experiments, which concentrate predominantly in a range λ/D = 120 − 250. Furthermore, good agreement is also observed with wavelengths predicted as λ/D = 105 − 300 in the stability analysis of Fourriere et al. (2010). Note that in their study, the dependence of the most amplified wavelength on the particle diameter appears indirectly in form of a scaling with a saturation length L sat , which in turn is a function of the particle diameter.
In contrast, the linear stability analysis of Colombini and Stocchino (2011) predicts only unstable wavelengths that scale with the fluid height for the range of Galileo numbers considered in the current study. It should be noted that since the relative submergence H f /D, which is a crucial parameter in their study, is at least one order of magnitude larger than in our simulations, a direct comparison of predicted wavelengths from their stability analysis and the values observed in our simulations was not possible. Generally, we could not observe such a scaling with the fluid height in our numerical simulations. For simulations with the same L x /D ratio, all relevant geometrical parameters attained comparable averaged values and profiles in the final phase irrespective of the varying L x /H f ratio. In addition, initial wavelengths λ/H f < 3 have been observed in experiments and simulations, excluding a scaling of the threshold which scales with the mean fluid height.
Another remarkable phenomenon has been observed in the cases with the highest relative submergence H f /D ≈ 50 in form of streamwise and spanwise sediment features, that are seen to either concur or interact with each other, allowing for the formation of three-dimensional patterns. The subsequent analysis of the single-sided amplitude spectra of the sediment bed height perturbation has shown that the evolution of the sediment bed is dominated by several streamwise and spanwise oriented waves of comparable amplitude. In the remaining simulations, streamwise aligned patterns appeared, if at all, mainly as a predecessor of the transverse patterns in the first few bulk time units of the simulations, but at clearly lower amplitude than their transverse oriented counterparts.
Since the focus of the current study was on the scaling of transverse patterns, it was out of the scope to investigate in detail the reasons for this different sediment bed evolution. However, it would be of high interest to determine the parameters that are responsible for the amplification of the spanwise waves. In this context, it should be investigated how initial subaqueous pattern formation changes when further increasing the Reynolds number.