5.1 The minimum and the most amplified unstable pattern wavelengths

Our strategy to find the minimal box length $L_{x}$  which will accommodate the lower bound of the unstable wavelengths $\unicode[STIX]{x1D706}_{th}$  (at a given parameter point) is to perform a series of numerical experiments in which the streamwise domain length is successively reduced. The concept is similar to the minimal flow unit of Jiménez & Moin (Reference Jiménez and Moin1991). Below a threshold value of $L_{x}=\unicode[STIX]{x1D706}_{th}$ , pattern evolution will be hindered and a perturbed bed should in principle be stable even though it is in a regime where instability is expected. Note that determining the lower threshold of $\unicode[STIX]{x1D706}_{th}$  for bed instability does not necessarily mean determining the most amplified wavelength $\unicode[STIX]{x1D706}_{c}$ . The latter, if it uniquely exists, exhibits the maximum growth rate of the bed instability.

Figure 2 shows instantaneous snapshots of the particle positions in the different cases after approximately 200 bulk time units have elapsed. It can be seen that the sediment bed of case H3 does not feature dune-like patterns. However, one can observe that it exhibits distinct streamwise aligned alternating ridges and troughs, although the sediment bed is essentially flat when the spanwise average of the sediment-bed profile is observed. Such an organization of particles is linked to the streamwise-aligned near-wall turbulent structures (streaks) which are responsible for the non-homogeneous spatial distribution of particles. The highly regular ridge-trough pattern in this case is a result of the very small box length which hinders spatial de-correlation of the turbulent structures in the streamwise direction. On the other hand, the sediment bed of the remaining cases with $L_{x}\geqslant 4H_{f}$ is seen to be unstable, allowing for the formation of spanwise-oriented patterns. The streamwise-aligned patterns are also visible superimposed to the ripples (in the remainder of the present text we will refer to the observed bedforms as ripples). This indicates that the value of $\unicode[STIX]{x1D706}_{th}$  lies somewhere between the domain lengths of cases H3 and H4. Another observation is that the initial wavelength of the emerging patterns seems to be different depending on the chosen domain length. In particular, the case with the largest domain size, which has a streamwise box length 12 times the smallest domain length of the ripple-featuring cases, is observed to accommodate more than ten initial ripple units. The dimension of the ripples in the large domain case seems to have evolved less constrained by the box size when compared to those in the small domain cases. Let us provide a more quantitative analysis of the above observation by evaluating the r.m.s. sediment-bed-height fluctuation $\unicode[STIX]{x1D70E}_{h}$  of each case and using it as a criterion to determine whether a sediment bed in a given simulation is stable or not. That is, if $\unicode[STIX]{x1D70E}_{h}$  remains bounded within a certain threshold value, then the bed is said to be stable. On the other hand, for an unstable bed, $\unicode[STIX]{x1D70E}_{h}$  will initially exhibit a continuous increase as a function of time. Note that even a stable flat (mobile) sediment bed features a small but finite value of $\unicode[STIX]{x1D70E}_{h}$  as a result of the random uncorrelated bed undulations which stem from the discreteness of the bed at the grain scale. For instance, Coleman & Nikora (Reference Coleman and Nikora2009a ) report $\unicode[STIX]{x1D70E}_{h}/D\approx 0.17$ for their macroscopically ‘flat’ bed.

Figure 2.

Top view of the three-dimensional instantaneous particle positions and the corresponding spanwise-averaged fluid–bed interface (shown below each snapshot) at time $t/T_{b}\approx 200$ of cases (a) H3, (b) H4, (c) H6, (d) H7, (e) H12 and (f) H48. The wall-normal particle location is indicated by a global colour code, as shown in (b).

Figure 3.

(a) Time evolution of the r.m.s. sediment-bed-height fluctuation normalized by the particle diameter: — (purple), case H3; — (orange), case H4; — (green), case H6; — (blue), case H7; — (brown), case H12; —, case H48. The horizontal dashed line indicates the value $\unicode[STIX]{x1D70E}_{h}/D\approx 0.17$ which corresponds to the random bed fluctuations of featureless flat bed (Coleman & Nikora Reference Coleman and Nikora2009a ). (b) Same as in (a) but the data for cases H4 and H12 are ensemble averaged over the number of simulations performed. The dashed-dot line corresponds to an exponential fit curve $\unicode[STIX]{x1D70E}_{h}/D=0.0668\exp (0.0140t/T_{b})$ obtained from an average amplitude growth in the initial exponential interval. (c) Same as in (b) but plotted in log-linear scale to highlight the exponential growth regime during the first instant of the bed instability.

Figure 4.

(a) Time evolution of the mean wavelength $\unicode[STIX]{x1D706}_{h}$  of the sediment-bed height normalized by the particle diameter. Colour coding similar to figure 3. (b) The same as in (a) but the data for cases H4 and H12 are ensemble averaged over the number of simulations performed.

Figure 3 shows the time evolution of $\unicode[STIX]{x1D70E}_{h}$  for all the considered cases listed in table 1. It is seen that, starting from time $tu_{b}/H_{f}\approx 20$ , the value of $\unicode[STIX]{x1D70E}_{h}$  grows with time for all the cases with domain length $L_{x}\geqslant 102.4D$ ( $L_{x}\geqslant 4H_{f}$ ), indicating that the chosen streamwise box dimension of these cases is larger than the lower threshold of unstable wavelengths, such that the box can accommodate at least one of the unstable modes. On the contrary, the evolution of $\unicode[STIX]{x1D70E}_{h}$  in case H3 exhibits no growth with time except for small fluctuations; $\unicode[STIX]{x1D70E}_{h}$  in this case is effectively bounded by the threshold value $\unicode[STIX]{x1D70E}_{h}/D\approx 0.17$ as per Coleman & Nikora (Reference Coleman and Nikora2009a ). This indicates that the streamwise box length $L_{x}=76.8D$ ( $L_{x}=3H_{f}$ ) is not sufficient to accommodate the lower threshold of unstable modes. Thus, it can be concluded that a cutoff length scale for pattern formation (for the considered parameter point) lies in the range $77D<\unicode[STIX]{x1D706}_{th}<102D$ , or in terms of the mean fluid height, $3H_{f}<\unicode[STIX]{x1D706}_{th}<4H_{f}$ . Additional simulations with $L_{x}$  in the range $3H_{f}$ to $4H_{f}$ would be required in order to determine $\unicode[STIX]{x1D706}_{th}$  more accurately.

Furthermore, figure 3 highlights the fact that the evolution of the pattern amplitude exhibits distinct regimes of growth. During the first approximately 150 to 180 bulk time units (excluding the first of approximately 20 bulk time units during which the sediment bed initially dilates after the start of the simulations), $\unicode[STIX]{x1D70E}_{h}$  is observed to evolve exponentially for all the cases, with a growth rate which seems to be independent of the chosen domain size. Note that the time evolution of $\unicode[STIX]{x1D70E}_{h}$  for cases H4 and H12 presented in figure 3(b,c) is an ensemble average over the three separate simulations of the respective cases. An exponential curve of the form

with $A=0.0668$ and $B=0.0140$ , best fits the average growth over all cases in this interval. In the subsequent growth regime, the trend of the amplitude evolution is observed to be markedly different among the different cases, although all were carried out at the same imposed parameter values. After the initial exponential growth (up to $t\approx 200T_{b}$ ), and after a small transition interval of approximately $100T_{b}$ , cases H4, H6 and H7 exhibit a plateau of the pattern amplitude, showing no further growth with time. The attained final values are $\unicode[STIX]{x1D70E}_{h}/D=0.93$ , $2.08$ and $2.36$ for each of these cases, respectively. The evolution of $\unicode[STIX]{x1D70E}_{h}$  for case H12 also attains a plateau regime with a final value of $\unicode[STIX]{x1D70E}_{h}$  comparable to that of case H6. However, it settles to this value not immediately after the exponential regime, rather it gradually increases (approximately linearly) in the interval between $t=200T_{b}$ and $600T_{b}$ . On the other hand, since the observation interval of case H48 is shorter, it has covered only the exponential growth regime by the end of the simulation.

The different trend of the amplitude growth observed is closely related to the influence of the limited computational box size on the initially accommodated mean wavelength and its subsequent evolution. In cases H4, H6, H7 and H12, even though the streamwise length of the computational box is larger than the threshold for pattern formation, due to the fact that $L_{x}/\unicode[STIX]{x1D706}_{th}=O(1)$ , the discrete wavelengths from which the system has to select are very sparse, and it is not guaranteed that the initially selected wavelength will be the same as the one which an infinitely long system would select. On the other hand, case H48 has a streamwise box length $L_{x}$  which is between 12–16 multiples of the minimum unstable wavelength. Thus the initial wavelength selected in case H48 is expected to be sufficiently close to that of an infinitely long system.

Figure 4 shows the time evolution of the mean pattern wavelength. In close correlation with the evolution of the amplitude, one can see that the adopted computational box size strongly influences the selected mean wavelength and its subsequent evolution. In cases H4 and H6, $\unicode[STIX]{x1D706}_{h}$  jumps within the first 100–200 bulk time units to the maximum possible wavelength, i.e.  $\unicode[STIX]{x1D706}_{h}=L_{x}$ , and evolves constrained to this value. This indicates that in these boxes, the system is forced to select the most unstable wavelength from the available modes, which turns out to be $\unicode[STIX]{x1D706}_{h}\approx 102D$ for case H4 and $\unicode[STIX]{x1D706}_{h}\approx 154D$ for case H6. The value of the mean wavelength of case H7 is observed to oscillate between the harmonics $\unicode[STIX]{x1D706}_{1}=L_{x}$ and $\unicode[STIX]{x1D706}_{2}=L_{x}/2$ and finally settles at $\unicode[STIX]{x1D706}_{h}=\unicode[STIX]{x1D706}_{2}\approx 179D$ at $t/T_{b}\approx 180$ . On the other hand, from the ensemble averaged evolution of wavelength for case H12, it can be seen that the system initially selects a wavelength $\unicode[STIX]{x1D706}_{h}\approx \unicode[STIX]{x1D706}_{3}\approx 102D$ , which is very close to the wavelength selected by the large box size case H48. It then grows monotonically and settles at a wavelength $\unicode[STIX]{x1D706}_{h}\approx \unicode[STIX]{x1D706}_{2}\approx 154D$ at approximately $t=500T_{b}$ , subsequently evolving constrained close to this wavelength until the end of the simulation interval. Again, the constraint of the evolution of the wavelength is an indication of the influence of box size on the subsequent evolution processes, although the box could be considered marginally sufficient to capture the initial wavelength. The nonlinear nature of the evolution process is also highlighted by observing the difference in the evolution of the wavelength for the three independent simulations of case H12 (which only differ from one another by the respective initial condition) in the first interval up to $t\approx 500T_{b}$ . In case H48, the availability of relatively fine graded harmonics allows for a more steady growth of the mean wavelength. In the initial time $t\approx 100T_{b}$ up to $t\approx 200T_{b}$ , the selected wavelength grows from $\unicode[STIX]{x1D706}_{h}\approx \unicode[STIX]{x1D706}_{12}=L_{x}/12$ to $\unicode[STIX]{x1D706}_{h}\approx \unicode[STIX]{x1D706}_{11}=L_{x}/11$ . It exhibits a further growth and attains a value $\unicode[STIX]{x1D706}_{h}\approx 145D$ at the end of the simulation interval.

The above-discussed evolution of the amplitude and mean wavelength are integral representations of the individual modes which make up the resolved discrete spectrum. On the other hand, from the stability analysis point of view, it might be more interesting to analyse the time evolution of the individual modes, since for small amplitudes these can be expected to grow independently. To this end, we have computed the instantaneous single-sided amplitude spectrum $A_{j}(t)$ (for non-negative wavenumbers $j\geqslant 0$ ) which is defined as twice the absolute value of the coefficient ${\hat{h}}_{bj}$ , where ${\hat{h}}_{bj}(t)$ is the $j$ th harmonic of the Discrete Fourier Transform of the bed-height perturbation $h_{b}^{\prime }(x,t)$ . For future reference let us denote with $S_{h_{b}h_{b}}(\unicode[STIX]{x1D705}_{j})$ the power spectral-density corresponding to the $j$ th Fourier mode. In order to provide a dispersion relation, first we assume that, in the initial exponential growth regime of the r.m.s. bed-height fluctuation (cf. figure 3), the individual modes exhibit an exponential growth as well. Strictly speaking, this assumption is only true when the bed fluctuations are so small that linear instability holds. Nevertheless, it can be safely assumed that the nonlinear interaction among the modes is weak in the first approximately $150$ – $180$ bulk time units (in the complementary hypothetical experiment, this time duration amounts to 7–9 s). Following, we fit an exponential curve of the form

to the time evolution of the amplitude of each mode in the interval between $t\approx 20T_{b}$ up to $t\approx 180T_{b}$ and infer the growth rate $\unicode[STIX]{x1D70E}_{j}(\unicode[STIX]{x1D705}_{j})$ , where $\unicode[STIX]{x1D705}_{j}\equiv 2\unicode[STIX]{x03C0}j/L_{x}$ is the wavenumber of the $j$ th harmonic.

Figure 5.

Amplitude growth rate $\unicode[STIX]{x1D70E}_{j}$ (as per definition (5.2)) as a function of the wavelength in the first exponential growth regime: — (purple), case H3; — (orange), case H4; — (green), case H6; — (blue), case H7; — (brown), case H12; —, case H48. The vertical dashed line represents the value of the streamwise box length in case H3 which does not feature patterns.

Figure 5 shows the growth rate as a function of the wavelength $\unicode[STIX]{x1D706}_{j}\equiv 2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D705}_{j}$ . As is expected, none of the modes are seen to grow at any significant rate in case H3. Moreover, consistent with the evolution of the mean wavelength, the fastest growing mode in cases H4 and H6 is $\unicode[STIX]{x1D706}_{1}=L_{x}$ , with a growth rate which is comparable to that of the r.m.s. bed height $\unicode[STIX]{x1D70E}_{h}$  (coefficient $B$ in (5.1)). The next few modes in these two cases also exhibit a small but noticeable growth rate although their wavelength is smaller than the threshold. This could be attributed to nonlinear effects. A similar trend is observable for cases H7 and H12, except that in these cases, it is not only one harmonic which is growing the fastest. In case H7 both $\unicode[STIX]{x1D706}_{1}=L_{x}$ and $\unicode[STIX]{x1D706}_{2}=L_{x}/2$ grow at a comparable rate, indicating that both harmonics contribute substantially to the growth of $\unicode[STIX]{x1D70E}_{h}$ . On the other hand, in case H12, in addition to the expected fastest growing modes $\unicode[STIX]{x1D706}_{3}=L_{x}/3$ and $\unicode[STIX]{x1D706}_{2}=L_{x}/2$ , the first harmonic $\unicode[STIX]{x1D706}_{1}=L_{x}$ is also seen to contribute to the overall growth. It is not guaranteed that the box length of H48 is large enough to capture all the unstable modes (i.e. the upper threshold of pattern formation). Nevertheless, the fact that $\unicode[STIX]{x1D70E}_{j}$ exhibits a clear local maximum is an indication that the box is sufficient to capture the most amplified mode(s), which turn out to be in the range $\unicode[STIX]{x1D706}_{8}$ to $\unicode[STIX]{x1D706}_{12}$ .

Finally, it is worth mentioning that the computational box sizes are chosen such that a turbulent flow state is sustained. This is confirmed by monitoring the time evolution of the box-averaged turbulent kinetic energy (plots not shown). Note that, in smooth wall channel flows, the minimum spanwise and streamwise box dimensions required in order for turbulence to be self-sustained are $L_{z}^{+}\approx 100$ and $L_{x}^{+}\approx 350$ , respectively (Jiménez & Moin Reference Jiménez and Moin1991). The minimum box dimensions adopted in our study (corresponding to case H3) are substantially larger, i.e.  $L_{x}^{+}=L_{z}^{+}\approx 740$ .

5.2 Ripple morphology

An additional important advantage of adopting a relatively ‘small’ computational box length in cases H4, H6 and H7 is that, once the maximum possible wavelength is attained, further ripple dimension growth is effectively hindered. Thus, the accommodated ripple subsequently evolves, steadily maintaining its dimensions and propagation velocity. Similarly, in case H12, within the simulated interval, the ripples evolve steadily, constrained approximately at a mean wavelength $\unicode[STIX]{x1D706}_{h}=L_{x}/2$ . Such a numerical manipulation, which experimentally is not possible, allows us to address steady-state bedform characteristics at a particular wavelength. In the following, we analyse the two-dimensional shape and the propagation velocity of the ripples.

It is well know that bedforms of finite amplitude evolve towards an asymmetric shape even in the absence of flow separation (Best Reference Best2005; Colombini & Stocchino Reference Colombini and Stocchino2008; Seminara Reference Seminara2010). For instance, spatio-temporal plots provided by Ouriemi et al. (Reference Ouriemi, Aussillous and Guazzelli2009) show that the evolution of the bedforms towards their asymmetric shape is a phenomenon of both the laminar ‘small dunes’ which exhibit no flow separation in their leeside, as well as the turbulent ‘vortex dunes’ which are characterized by flow separation downstream of their crests. The degree of asymmetry is generally larger in the latter as a result of the flow recirculation. Moreover, the animations of particle motion available online at https://doi.org/10.1017/jfm.2014.284 and the space–time plots depicted in figure 4 of Kidanemariam & Uhlmann (Reference Kidanemariam and Uhlmann2014a ) have qualitatively shown the asymmetric evolution of sediment-bed shape in good agreement with what is observed experimentally.

In order to quantitatively characterize the statistically two-dimensional shape of the bedforms, in the steady ripple evolution interval of cases H4, H6, H7 and H12, we have evaluated the phase-averaged fluid–bed interface which is defined as follows:

where $\tilde{x}\equiv x-u_{D}t$ is the $x$ -coordinate in a frame of reference which is moving at the mean ripple propagation velocity $u_{D}$  (see § 5.3 for the definition of $u_{D}$ ). Note that, in cases H4, H6 and H7, $\widetilde{H}_{b}$  corresponds to the mean profile of a single ripple, while in case H12, $\widetilde{H}_{b}$  is an average profile over the two available ripples.

Figure 6.

(a) Mean two-dimensional profile of the patterns averaged over the last steady ripple propagation interval, and over the number of ripples when applicable, of each case. The location of the maximum of each profile is shown by the filled symbols. Note that the profiles are plotted not to scale and thus the aspect ratio is visually exaggerated. (b) The corresponding spectra of the sediment-bed-height fluctuation. In both plots: — (orange), case H4; — (green), case H6; — (blue), case H7; — (brown), case H12. Power laws with exponent $-3.3$ and $-4.3$ are shown by the red and blue dashed lines respectively. The vertical dashed line represents the value of the average fluid height $H_{f}\approx 25.27D$ (averaged over all cases).

As shown figure 6(a), it can be seen that in all cases, the ripple profile is clearly asymmetric with an upstream face (upstream of the crest) approximately three times longer in streamwise length and thus milder in bed slope than the downstream face. To within the statistical uncertainty, the bed shape of cases H6, H7 and H12, once scaled by the corresponding mean wavelength, seems to fairly collapse into a single shape, exhibiting small differences. The attained value of the aspect ratio, which is defined as the ratio between the height of the ripple $H_{D}=\max (\widetilde{H}_{b})-\min (\widetilde{H}_{b})$ (i.e. the wall-normal distance between the trough and the crest) and the mean wavelength $\unicode[STIX]{x1D706}_{h}$ , is approximately 0.037. Moreover, the degree of asymmetry, which can be defined as the ratio between the streamwise distance from the crest to the downstream trough and the mean wavelength, is approximately 0.28. On the other hand, the mean ripple shape of case H4 is seen to exhibit a noticeable smaller value of the aspect ratio (approximately 0.026) and a larger value of the degree of asymmetry (approximately 0.35) than the other cases. This difference is an indication of the ripple shape evolution process. Some experimental studies reported that, during the natural ripple coarsening process, the aspect ratio of the ripple geometry evolves with time and attains a value of approximately $1/15$ ( $0.067$ ) once the steady-state shape is attained (Charru et al. Reference Charru, Andreotti and Claudin2013). On the other hand, Fourrière, Claudin & Andreotti (Reference Fourrière, Claudin and Andreotti2010), based on field measurements, report a steady-state aspect ratio of approximately $0.045$ . In our configurations, since the final ‘natural’ steady-state regime is not yet reached, we cannot confidently say that the attained values of the aspect ratios are representative of those values which would naturally be attained. Nevertheless, the fact that the shape of the ripple with mean wavelength $\unicode[STIX]{x1D706}_{h}\approx 150D$ and that with $\unicode[STIX]{x1D706}_{h}\approx 180D$ do not exhibit substantial differences could be an indication that the ripple shape is not far from its final invariant shape. This aspect requires further investigation.

It is generally believed that the nonlinear interaction between the evolving sediment bed and the driving flow results in the above mentioned asymmetric ripple shape (Charru et al. Reference Charru, Andreotti and Claudin2013). That is, although the ripples have one well-defined length scale (the mean wavelength $\unicode[STIX]{x1D706}_{h}$  defined in (4.6)), a wide spectrum of modes interact nonlinearly to maintain a shape-invariant ripple which propagates downstream. There are several theoretical studies as well as field and experimental measurements which report that the spectra of fully developed bedforms exhibit a $-3$ power-law relationship with respect to the large wavenumbers, bounded by a threshold value (Hino Reference Hino1968; Jain & Kennedy Reference Jain and Kennedy1974; Nikora, Sukhodolov & Rowinski Reference Nikora, Sukhodolov and Rowinski1997; Coleman & Nikora Reference Coleman and Nikora2011). To gain further insight, we present in figure 6(b) the spectra $\langle S_{h_{b}h_{b}}\rangle _{t}$ of the sediment-bed-height fluctuation averaged over the steady ripple evolution interval of the above-mentioned cases. It can be seen that, for the dominant modes (with $\unicode[STIX]{x1D706}_{j}\gtrsim H_{f}$ ), $\langle S_{h_{b}h_{b}}\rangle _{t}$ of all the cases (only for the even modes of case H12) features a power-law variation as a function of $\unicode[STIX]{x1D705}_{j}$ . The value of the scaling exponent in cases H6, H7 and for the even modes in H12 is observed to be approximately $-3.3$ whereas that of case H4 is approximately $-4.3$ . The fact that case H4 has a larger value of the exponent than that of the other cases is again an indication of the evolution process. That is, the influence on the mean ripple shape of the modes with a smaller wavelength than the dominant one is relatively smaller compared to the other cases. Moreover, the fact that we have not attained a naturally developed mature ripple state could be a reason why the attained exponents are slightly larger than the $-3$ value reported in the literature. The reason for the odd–even separation of the spectra in case H12 could be explained by the fact that the domain length accommodates two ripples and thus modes with wavelength $L_{x}/(2i)$ mainly contribute to the attained shape.

5.3 Ripple migration velocity

The migration velocity of bedforms is an important quantity of interest in the study of morphodynamics. For instance, the rate at which bedforms propagate is closely related to the rate of sediment transport (Coleman & Melville Reference Coleman and Melville1994; Nikora et al. Reference Nikora, Sukhodolov and Rowinski1997; Betat et al. Reference Betat, Kruelle, Frette and Rehberg2002; Coleman et al. Reference Coleman, Fedele and Garca2003; Coleman & Nikora Reference Coleman and Nikora2009b ; Ouriemi et al. Reference Ouriemi, Aussillous and Guazzelli2009; Seminara Reference Seminara2010). From the sediment mass-conservation equation (Exner equation), it can be easily shown that, for a statistically shape-invariant two-dimensional ripple with height $H_{D}$ and which migrates downstream at a constant velocity $u_{D}$ ,

where $\langle q_{p}\rangle$  is the mean particle flow rate, which is assumed to be constant, while $\langle q_{p,min}\rangle$ is the corresponding minimum value at the trough of the ripple, $\unicode[STIX]{x1D6F7}_{b}$  is the mean solid volume fraction of the sediment bed and $\unicode[STIX]{x1D6FD}$ is a ripple shape parameter. From the approximation (5.4), it is evident that the ripple migration velocity is inversely proportional to the ripple height (and thus the ripple mean wavelength) whereas it is directly proportional to the mean particle flow rate (Charru et al. Reference Charru, Andreotti and Claudin2013). In the following, based on the DNS data, we evaluate the mean ripple migration velocity and assess its relation to the mean particle flow rate.

In the aforementioned steady evolution interval, an average migration velocity of the patterns $u_{D}$  can be determined from the shift of the maximum of space–time correlation function of the fluid–bed interface fluctuation $h_{b}^{\prime }$ (Nikora et al. Reference Nikora, Sukhodolov and Rowinski1997; Coleman & Nikora Reference Coleman and Nikora2011). First, the space–time correlation function is defined as follows:

where $t_{1}$ and $t_{2}$ are the start and end times of the considered interval. The streamwise location $x_{max}$ of the maximum of $R_{ht}$ , which is given by

varies approximately linearly with respect to $\unicode[STIX]{x1D6FF}t$ , and thus the slope of the line which is fitted to the data $r_{x}^{m}$ versus $\unicode[STIX]{x1D6FF}t$ gives an accurate estimate of $u_{D}$ .

Figure 7.

(a) Average propagation speed of the patterns, normalized in wall units, as a function of the mean wavelength (computed in the final steady ripple evolution interval as shown in table 2): ▴ (orange), case H4; ◂ (green), case H6; ▸ (blue), case H7; ● (brown), case H12. The corresponding cross-marks ( $\times$ ) represent the mean ripple migration velocity estimated by relation (5.4). The data points in the turbulent channel flow experiment of Cardona Florez & Franklin (Reference Cardona Florez and Franklin2016) for their particles with diameter $D=550~\unicode[STIX]{x03BC}\text{m}$ are shown by the open symbols ○ (initial phase of ripple evolution) and ▫ (final phase of the evolution). (b) Mean propagation speed of individual Fourier modes in the final steady ripple evolution interval: — (orange), case H4; — (green), case H6; — (blue), case H7; — (brown), case H12. The corresponding mean ripple migration velocity is plotted by the horizontal dashed lines. The vertical dashed line represents the value of the average fluid height $H_{f}\approx 25.27D$ (averaged over all cases). In both plots, ensemble averaging is performed for cases H4 and H12 over the available number of runs.

Figure 7(a) shows the migration velocity, normalized by $u_{\unicode[STIX]{x1D70F}}$ , as a function of the mean wavelength for cases H4, H6, H7 and H12. The variation of $u_{D}$  exhibits an inverse relationship with the ripple wavelength, varying from a value $u_{D}^{+}\approx 0.49$ for case H4 down to a value $u_{D}^{+}\approx 0.3$ for case H7. It should be noted, however that, as will be shown in § 5.5, as a consequence of an increase of the mean bottom shear stress, the value of the mean particle flow rate increases with increasing mean wavelength. Thus the observed trend of the ripple migration velocity is a net result of an increase of particle flow rate as well as an increase of ripple dimensions. In the following, we assess the mass-conservation-based relationship between the mean particle flow rate and the ripple migration velocity, by reverse-computing the latter from relation (5.4). To this end, a mean particle flow rate $\langle q_{p}\rangle$  is computed over the steady ripple migration interval (cf. § 5.5). The average minimum particle flow rate at the ripple troughs $\langle q_{p,min}\rangle$ was computed as well. It turns out that the value $\langle q_{p,min}\rangle$ is negligibly small in cases H6, H7 and H12, whereas in case H4 it attains a value of $\langle q_{p,min}\rangle \approx 0.2\langle q_{p}\rangle$ . Furthermore, the ripple shape parameter $\unicode[STIX]{x1D6FD}$ is defined as the ratio between the two-dimensional area of the ripple to that of a bounding rectangle of dimensions $H_{D}$ and $\unicode[STIX]{x1D706}_{h}$ , defined as

In most of the engineering-purpose bedload transport rate estimators, the value of $\unicode[STIX]{x1D6FD}$ is typically taken to be $0.5$ , assuming that the shape of the ripple is triangular, while some adopt values up to $\unicode[STIX]{x1D6FD}=0.6$ based on field and experimental measurements (see e.g. Gaeuman & Jacobson (Reference Gaeuman and Jacobson2007) and references therein). For the mean ripple geometries encountered in the present study, a value of $\unicode[STIX]{x1D6FD}$ in the range between $0.5$ and $0.52$ is recovered after evaluating (5.7) for all the cases. Subsequently, the mean ripple migration velocity is estimated based on relation (5.4) and the result is presented in figure 7(a). It can be seen that the estimation yields values which are in good agreement with the actual computed ripple migration velocity values.

Additionally, figure 7(a) shows recent experimental measurements of ripple migration velocity in a closed-conduit set-up (Cardona Florez & Franklin Reference Cardona Florez and Franklin2016). The experimental data shown correspond to particles with diameter $D=550~\unicode[STIX]{x03BC}\text{m}$ ( $Re_{b}\approx 7600{-}8900$ , $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=2.5$ , $H_{f}/D\approx 39$ , $\unicode[STIX]{x1D703}\approx 0.056{-}0.075$ , $D^{+}\approx 12{-}14$ , $Ga\approx 50$ ) which is not far from our parameter point. Although the experimental data exhibit relatively large scatter, it can be seen that the trend of the DNS data falls within the experimental data cloud.

Finally, let us turn to the celerity of the individual Fourier modes of the bedforms. The phase angle spectrum $\unicode[STIX]{x1D711}_{j}$ of the sediment-bed-height fluctuation is defined as

where ‘atan2’ is the four-quadrant inverse tangent. The instantaneous value of the phase speed $c_{j}$ is derived from (5.8) as:

In the steady ripple propagation interval of cases H4, H6, H7 and H12, we have evaluated an average mode-wise phase speed $\langle c_{j}\rangle$ which is presented in figure 7(b) as a function of the wavelength $\unicode[STIX]{x1D706}_{j}$ . In all the considered cases, the celerity data feature two distinct types of modes: dispersive and non-dispersive. All harmonics with a wavenumber value smaller than a threshold are observed to be non-dispersive, propagating at the corresponding mean ripple migration velocity. To within the statistical uncertainty of the data, the threshold value is the same in all cases and turns out to be $\unicode[STIX]{x1D706}_{j}\approx H_{f}$ . Note that the non-dispersive modes are those dominant ones which exhibit a power-law relationship in figure 6. On the other hand, the remaining harmonics with a wavelength smaller than the threshold are observed to be dispersive, exhibiting approximately the same variation of the phase speed as a function of the wavelength in all cases. The trend is seen to be well defined when scaled by the friction velocity, increasing from its value at the threshold location and attaining values of $\langle c_{j}\rangle ^{+}\approx 2$ at wavelengths of approximately $\unicode[STIX]{x1D706}_{j}\approx 10D$ . At first glance, the observations of dispersive modes seems incompatible with the fact that the ripples are migrating at a constant velocity maintaining their two-dimensional shape. Looking at figure 6, however, the contribution of the fast-moving modes to the mean ripple shape (and mean migration velocity) is orders of magnitude smaller than that of the most dominant one, thus negligible. These modes are footprints of the particle erosion/deposition cycle which takes place on the top of the ripples. Based on field measurements, Nikora et al. (Reference Nikora, Sukhodolov and Rowinski1997) similarly report that fully developed ripples are composed of non-dispersive modes with wavelength $\unicode[STIX]{x1D706}\gg H_{f}$ , as well as dispersive modes with $\unicode[STIX]{x1D706}\ll H_{f}$ , corroborating our findings.

5.4 Mean shear stress at the fluid–bed interface

For a given imposed flow rate, an increase of the mean interface shear stress with increasing ripple dimensions is expected since the ripples are effectively roughness elements of the sediment bed (cf. Jiménez Reference Jiménez2004, for instance). This in turn is expected to influence the mean rate of particle transport in a given channel, since it is directly proportional to the mean shear stress. The inter-dependency among the evolving ripple, the mean interface shear stress and the mean particle flow rate can be assessed by evaluating these quantities on the steady interval of each case (cf. table 2).

Let us recall that the driving volume force exerted by the imposed mean pressure gradient $\langle \text{d}p/\text{d}x\rangle$  in a particle-laden channel flow is balanced by the sum of resisting force of the fluid shear stress and stress contribution from the fluid–particle interaction. In the context of the immersed boundary method, the volume force, which imposes the no-slip condition at the fluid–particle interface, corresponds to the stress imposed on the system as a result of the fluid–solid interaction (see e.g. Uhlmann Reference Uhlmann2008; Kidanemariam et al. Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013). When the momentum equation is averaged over the entire domain comprising the fluid and particles, the streamwise momentum balance reduces to

where the last term in the right-hand side represents the fluid–solid interaction, $\langle u\rangle$ is the mean composite velocity, and $\langle u^{\prime }v^{\prime }\rangle$ is the covariance with respect to $\langle u\rangle$ . When (5.10) is applied to a statistically one-dimensional channel flow configuration (for instance case H3), $\langle u^{\prime }v^{\prime }\rangle$ represents the momentum flux due to turbulent fluctuations. On the other hand, the mean flow in those ripple-featuring cases is two-dimensional and $\langle u^{\prime }v^{\prime }\rangle$ represents not only the turbulent fluctuations, but also contribution from the deviation of the mean flow streamlines from the streamwise direction (see e.g. Yalin Reference Yalin1977; Nikora et al. Reference Nikora, McEwan, McLean, Coleman, Pokrajac and Walters2007). A detailed analysis of the flow field over bedforms is not part of the current study and will appear in a followup paper. Here, in order to define a mean interface shear stress, evaluating (5.10) suffices. In the steady ripple propagation interval of the current configuration (cf. table 2), $\langle \unicode[STIX]{x1D70F}_{tot}\rangle$  should vary linearly as a function of wall-normal distance. This is confirmed in figure 8(a), which shows the wall-normal variation of the different contributions to the total shear stress for case H7. The remaining cases exhibit similar trends (plots not shown). The figure highlights that, in the region above the highest crest of the ripples, the driving force is dominantly balanced by the fluid stress term while in the domain below the lowest trough of the ripples, the particle-related forcing entirely balances the pressure gradient forcing. The mean interface shear stress ( $\unicode[STIX]{x1D70C}_{f}u_{\unicode[STIX]{x1D70F}}^{2}$ ) is thus defined as the value of $\langle \unicode[STIX]{x1D70F}_{tot}\rangle$  at the location of the mean fluid–bed interface $y_{0}=H_{b}$ . It should be noted that, in the ripple-featuring cases, due to the roughness introduced as a result of the evolving bedforms, a horizontal fluid–bed interface does not exist. Nevertheless, a virtual wall could be defined to be located at a wall-normal location $y_{0}=H_{b}$ .

Figure 8.

(a) Wall-normal profiles of the different contributions to the total shear stress of case H7 (computed in the final steady ripple evolution interval as shown in table 2): — (magenta), $\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D708}\unicode[STIX]{x2202}_{y}\langle u\rangle$ ; — (blue), $-\unicode[STIX]{x1D70C}_{f}\langle u^{\prime }v^{\prime }\rangle$ ; — (red), $\int _{y}^{L_{y}}\langle f_{x}\rangle \,\text{d}y$ ; —, $\langle \unicode[STIX]{x1D70F}_{tot}\rangle$ . The horizontal dashed lines represent the wall-normal location of mean fluid–bed interface as well as its all-time maximum and minimum extents over the considered observation interval. (b) Friction velocity based Reynolds number as a function of the mean wavelength $\unicode[STIX]{x1D706}_{h}$ in the corresponding steady ripple evolution interval: ▴ (orange), case H4; ◂ (green), case H6; ▸ (blue), case H7; ● (brown), case H12. The horizontal dashed line corresponds to the value of $Re_{\unicode[STIX]{x1D70F}}$  in the featureless case H3.

Figure 8(b) shows the mean interface shear stress (expressed in terms of the friction velocity based Reynolds number) as a function of the mean ripple wavelength for the different cases (computed in the final steady ripple evolution interval as shown in table 2). It can be seen that the value of the friction Reynolds number in the featureless case H3 is $Re_{\unicode[STIX]{x1D70F}}\approx 245$ . This value is larger than the value of $Re_{\unicode[STIX]{x1D70F}}$  in a rough-wall channel flow at a comparable bulk Reynolds number and at roughly the same particle Reynolds number (cf. Chan-Braun, García-Villalba & Uhlmann Reference Chan-Braun, García-Villalba and Uhlmann2011). The main difference between the latter case and the present case H3 is the particle arrangement and particle mobility, i.e. the erosion, entrainment and re-deposition cycle. Furthermore, when comparing the friction Reynolds number of cases H4, H6, H7 and H12, it can be observed that there is a monotonic increase with increasing wavelength, attaining values $Re_{\unicode[STIX]{x1D70F}}\approx 265$ for $\unicode[STIX]{x1D706}_{h}\approx 102D$ , $Re_{\unicode[STIX]{x1D70F}}\approx 300$ for $\unicode[STIX]{x1D706}_{h}\approx 154D$ and $Re_{\unicode[STIX]{x1D70F}}\approx 310$ for $\unicode[STIX]{x1D706}_{h}\approx 180D$ . Recalling the fact that the imposed fluid flow rate and the mean fluid height are essentially the same in all cases, the increase in the shear stress among the ripple-featuring cases is entirely a consequence of the increase of the amplitude of the evolving sediment patterns (i.e. the macroscopic roughness height).

5.5 Mean particle flow rate

The instantaneous volumetric flow rate of the particle phase (per unit span), $q_{p}$ , is given by the sum (over all particles) of the streamwise particle velocity times the particle volume, divided by the product of the streamwise and spanwise extent of the domain (Kidanemariam & Uhlmann Reference Kidanemariam and Uhlmann2014 b ), viz

where $u_{p}^{(l)}(t)$ is the streamwise component of the instantaneous velocity of the $l$ th mobile particle at time $t$ . Averaging $q_{p}$ over a given stationary interval results in the mean particle flow rate $\langle q_{p}\rangle$  in that interval. We remark that, for the parameter point considered in the present study, particles are dominantly transported as ‘bedload’ material. Only a very small fraction of the mobile particles is suspended and entrained by the mean flow. Thus, the particle flow rate defined in (5.11) represents overwhelmingly the former mode of transport. In figure 9(a), as a consequence of the increased mean bottom friction, the mean particle flow rate is observed to increase with increasing ripple dimensions.

Of particular relevance to engineering applications is to express the particle flow rate as a function of the Shields number and to pursue scaling laws which relate the two quantities. In order to analyse the particle flow rate obtained in the present work in light of such scaling laws, the non-dimensional particle flow rate (normalized by the inertial scale $q_{ref}=U_{g}D$ ) is shown as a function of the excess Shields number $\tilde{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{c}$ in figure 9(b), using the value $\unicode[STIX]{x1D703}_{c}=0.034$ according to the empirical law proposed by Soulsby & Whitehouse (Reference Soulsby and Whitehouse1997). Since all the ripple-featuring cases have started from an initially flat sediment bed, the values of both the shear stress and the particle flow rate in these cases increase with time from their initial values to their final values at the end of the simulation interval. In order to capture this time evolution, the entire observation interval of each case is decomposed into smaller intervals of approximately 25 bulk time units. A mean particle flow rate and bottom shear stress is then computed for each interval, substantially increasing the data samples in the figure. The duration of the time intervals is chosen to be much smaller than the time scales of the ripple evolution. Figure 9(b) also shows the empirical power law of Wong & Parker (Reference Wong and Parker2006), which in turn is a modified version of the Meyer-Peter & Müller (Reference Meyer-Peter and Müller1948) formula for turbulent flows,

where $A=4.93$ and $\unicode[STIX]{x1D6FC}=1.6$ . Note that Wong & Parker’s formula is valid for a macroscopically flat sediment bed. As is observed in the figure, the DNS data points, which represent both plane sediment beds as well as pattern-featuring beds, are in very good agreement with the scaling law (5.12). Although the evolution of the ripples does increase the net particle transport rate, the net bottom shear stress simultaneously increases. As a net result, the formation of patterns does not strongly affect the particle transport scaling as a function of the excess shear stress (in the considered observation interval).

Figure 9.

(a) Volumetric particle flow rate $\langle q_{p}\rangle$  as a function of the mean wavelength $\unicode[STIX]{x1D706}_{h}$ (computed in the final steady ripple evolution interval as shown in table 2): ▴ (orange), case H4; ◂ (green), case H6; ▸ (blue), case H7; ● (brown), case H12. The horizontal dashed line corresponds to the value of $\langle q_{p}\rangle$  in the featureless case H3. (b) Particle flow rate $\langle q_{p}\rangle _{i}$ as a function of the excess Shields number $\tilde{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{c}$ . The entire observation interval of each case is decomposed into smaller intervals of approximately 25 bulk time units: ▴ (purple), case H3; ▴ (orange), case H4; ◂ (green), case H6; ▸ (blue), case H7; ● (brown), case H12; ▪, case H48. A mean particle flow rate and bottom shear stress is then computed for each interval (thus the notations $\langle q_{p}\rangle _{i}$ and $\unicode[STIX]{x1D703}_{i}$ ). Note that values for every fifth interval are plotted for clarity. The value $\unicode[STIX]{x1D703}_{c}=0.034$ is used in order to compute $\tilde{\unicode[STIX]{x1D703}}$ for the current DNS data, according to the empirical law proposed by Soulsby & Whitehouse (Reference Soulsby and Whitehouse1997). The dashed line is the Wong & Parker (Reference Wong and Parker2006) version of the Meyer-Peter & Müller (Reference Meyer-Peter and Müller1948) formula for turbulent flow which reads: $\langle q_{p}\rangle /q_{ref}=4.93\,\tilde{\unicode[STIX]{x1D703}}^{1.6}$ .