4.1. Spatial structure of the mean flow

The streamwise and wall-normal components of the phase-averaged fluid velocity vector $\tilde {\boldsymbol {U}}_f\equiv ({{\tilde {U}_f}},{{\tilde {V}_f}})$ and the corresponding particle velocity vector ${{\tilde {\boldsymbol {U}}_p}}\equiv ({{\tilde {U}_p}},{{\tilde {V}_p}})$ for case H7 are shown in figure 4 (cf. Appendix B for the definitions of $\tilde {\boldsymbol {U}}_f$ and ${{\tilde {\boldsymbol {U}}_p}}$). The figures highlight significant modulation of the flow by the ripple morphology. As is expected, below the fluid–bed interface, mean fluid and particle velocities are found to be negligibly small. Above the fluid–bed interface, on the other hand, velocities of both phases show large spatial variations that are seen to be strongly correlated with the ripple geometry. At and near the interface, ${{\tilde {U}_f}}$ is observed to be higher in the region above the crest of the ripple. The value of ${{\tilde {U}_f}}$ decreases downstream of the crest, attaining negative values in the re-circulation region highlighted by the negative ${{\tilde {U}_f}}$ contour shown in figure 4(a). Let us remark that a zero-valued contour level does not precisely demarcate the extent of the re-circulation region as there is very small but non-zero mean velocity inside the bed. The value of ${{\tilde {U}_f}}$ gradually increases along the ripple stoss side up to the crest. Similarly, in the outer flow region, ${{\tilde {U}_f}}$ exhibits higher velocity values in the flow contraction region above the crest and lower velocity values in the expansion region above the ripple stoss side. The wall-normal fluid velocity also exhibits an upward moving fluid region (${{\tilde {V}_f}}>0$) above the stoss side and a downward moving fluid region (${{\tilde {V}_f}}<0$) above the region downstream of the crest. A small region of positive ${{\tilde {V}_f}}$ is observed immediately downstream of the crest in the vicinity of the fluid–bed interface that corresponds to the up-lifting fluid of the re-circulation region. As $y$ tends to the upper boundary, ${{\tilde {V}_f}}$ tends to zero due to the imposed boundary condition.

Figure 4. (a) Phase-averaged fluid velocity components (${{\tilde {U}_f}}^{+},{{\tilde {V}_f}}^{+}$) of case H7. Contour levels of ${{\tilde {U}_f}}^{+},{{\tilde {V}_f}}^{+}= 0.1$ and $-0.1$ are marked by red and magenta curves. The phase-averaged sediment-bed height, ${\tilde {H}_b}$, is shown by the dashed curve. The vertical thin lines indicate the locations where wall-normal profiles of each quantity are extracted. (b) The corresponding phase-averaged particle velocity components (${{\tilde {U}_p}}^{+},{{\tilde {V}_p}}^{+}$). (c) Values of ${{\tilde {U}_f}}$ and ${{\tilde {U}_p}}$ extracted along the fluid–bed interface (shown by the curve at the top of the figure). The red dashed line represents the mean ripple migration velocity. (d) Same as panel (c) but for ${{\tilde {V}_f}}$ and ${{\tilde {V}_p}}$. (e) Wall-normal profiles ${{\tilde {U}_f}}$ and ${{\tilde {U}_p}}$ at selected streamwise locations. Profiles are, from left to right, at ${{\tilde {x}}}/D = 0$, $30$ and $112$ consecutively. The horizontal dashed lines indicate the location of the fluid–bed interface at the corresponding streamwise locations. ( f) Same as panel (e) but for ${{\tilde {V}_f}}$ and ${{\tilde {V}_p}}$.

The mean particle velocities also exhibit generally the same spatial structure as that of the fluid. From the intensity of the colour plots, it can be seen that the magnitude of streamwise particle velocity is on average smaller than that of the fluid in the region far above the fluid–bed interface, except in the re-circulation region. The opposite is true for the wall-normal component ${{\tilde {V}_p}}$, which is observed to be higher in magnitude than that of the fluid. More quantitative comparison of the fluid and particle velocities can be inferred from the profiles along the location of the fluid–bed interface shown in figure 4(c,d) as well as from the wall-normal profiles of the same quantities at selected streamwise locations ${{\tilde {x}}}/D = 0$ (near crest), ${{\tilde {x}}}/D = 30$ (in trough) and ${{\tilde {x}}}/D = 112$ (on the stoss side of the ripple) shown in figure 4(e, f). It can be seen that the fluid/particle velocities strongly vary along the fluid–bed interface. Note that the fluid–bed interface is permeable and is subject to active mass and momentum exchange of both phases as a result of erosion, transport and deposition of sediment particles as well as the fluid motion therein. At the ripple crest the streamwise velocity of both phases reaches local maximum values of up to ${{\tilde {U}_f}}^{+}\approx 4$ and ${{\tilde {U}_p}}^{+} \approx 3$, whereas in the ripple troughs, velocities decrease considerably, attaining negative values in the re-circulation region. Further observation is that particles are on average moving at a smaller streamwise velocity than that of the fluid along the fluid–bed interface, except in the re-circulation region, where the opposite trend is observed. On the stoss side of the ripples, the value of the velocity lag between the two phases ranges between $1u_{\tau }$ and $1.5u_{\tau }$, whereas in the re-circulation region, a negative velocity lag of up to $-0.35u_{\tau }$ is observed. The origin of the velocity lag could be explained by the combined effect of the preferential distribution of inertial particles with respect to the near-wall high- and low-speed fluid regions (Kidanemariam et al. Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013) and the fact that finite-sized particles feel the traction effect of neighbouring particles, which are located at a lower wall-normal location and are moving at a lower velocity, more than fluid particles do. Note that the particle velocity is much higher than the mean ripple migration velocity ${{u_{D}}}$ over the entire stoss side of the ripples, while the opposite is true at locations downstream of the crest. This observation is in line with the fact that ripples propagate as a result of erosion and deposition of sediment grains at the bed surface. ${{\tilde {V}_f}}$ and ${{\tilde {V}_p}}$ along the fluid–bed interface exhibit a notably different trend when compared with their streamwise counterparts. On the stoss side, ${{\tilde {V}_f}}$ and ${{\tilde {V}_p}}$ are positive while the difference between their values is not substantial. Moreover, the location of the local maximum of ${{\tilde {V}_f}}$ (${{\tilde {V}_p}}$) is observed to be located further upstream of the crest when compared with the streamwise components. Both phases exhibit an approximately zero value of the wall-normal velocity in the vicinity of the crest.

The wall-normal profiles of fluid and particle velocities presented in figure 4(e, f) further corroborate the above observations. Note that the apparent velocity lag is significant even at wall-normal locations well above the fluid–bed interface. Namely, positive velocity lag on the stoss side of the ripples (${{\tilde {U}_f}}>{{\tilde {U}_p}}$) and negative velocity lag in the re-circulation region (${{\tilde {U}_f}}<{{\tilde {U}_p}}$). In the region well above the fluid–bed interface, the preferential sampling of low-speed fluid regions by inertial particles (Kidanemariam et al. Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013) is expected to be the dominant reason for the observed positive velocity lag. The negative velocity lag in the re-circulation region is again attributed to particle inertia. Sediment particles that are suspended from the crest with high momentum into the re-circulation region do not immediately adapt to the new environment, thus resulting in a statistically higher mean particle velocity than that of the fluid. Turning to the wall-normal velocity component, the profiles further highlight the mean positive and negative velocities upstream and downstream of ripple crests, respectively. Furthermore, upstream of the crest and above the fluid–bed interface, the particles’ wall-normal velocity is observed to be larger in magnitude than that of the fluid whereas downstream of the crest and above the fluid–bed interface, the particle velocity is observed to be substantially smaller (larger in absolute value) than that of the fluid. This latter effect is in part attributed to the effect of gravity on the suspended sediment grains that settle in the re-circulation region.

Finally, we remark that the corresponding results for cases H4 and H6 are in general of similar spatial structure and plots are not included. In figure 5 we report the streamline plot of the mean flow for the three cases. The plots effectively show the development of the re-circulation region at different stages of the ripple evolution. Incidentally, the streamlines also show the structure of the minute flow inside the sediment bed. The overall trend of the spatial variation of the mean streamwise and wall-normal components of the fluid velocity is consistent with the experiments by Charru & Franklin (Reference Charru and Franklin2012) of the flow over barchan dunes. However, a rigorous analysis of the structure of the mean flow and turbulent statistics of the fluid and particle phases, as well as comparison with available experiments, is left for future work.

Figure 5. Streamlines of the phase-averaged fluid velocity (close-up in the region downstream of the ripple crests). The grey-scale plot in the background represents the mean solid volume fraction $\tilde {\phi }_p$. (a), H4; (b), H6; (c), H7.

4.2. Total shear stress and its spatial variation along the interface

As has been reported in KU2017, the driving mean pressure gradient in the flow is balanced by the sum of the fluid shear stress and stress contribution from the fluid–particle interaction resulting in a linearly varying plane-averaged total shear stress across the depth of the channel. This linearity can also be recovered from the phase-averaged statistics. In the two-dimensional flow evolving over the ripples, the total plane-averaged fluid shear stress comprises the viscous and turbulent Reynolds stresses as well as contribution from the ripple form-induced dispersive stress (Raupach & Shaw Reference Raupach and Shaw1982; Nikora et al. Reference Nikora, McEwan, Mclean, Coleman, Pokrajac and Walters2007a). Similarly, the Eulerian particle–fluid interaction contribution (the last term in 5.10 in Kidanemariam & Uhlmann Reference Kidanemariam and Uhlmann2017) could in principle be recovered from its phase-shifted counterpart $\langle \tilde {f}_x \rangle _{zt}$. However, the latter quantity is not explicitly available for evaluation because we have only accumulated time history of the Lagrangian hydrodynamic ($\boldsymbol {F}_p^{H(l)}$) and collision ($\boldsymbol {F}_p^{C(l)}$) forces for each particle $l$ by integrating over the particle surface. As has been demonstrated by Uhlmann (Reference Uhlmann2008), the mean Eulerian force density can be approximated from the total hydrodynamic force acting on the particles, viz.

(4.1)\begin{equation} \langle \tilde{\boldsymbol{f}} \rangle_{zt} \approx{-}\langle \tilde{\boldsymbol{f}}^{H}_p\rangle_{zt} \tilde{\phi}_p/V_p, \end{equation}

where the field $\boldsymbol {f}^{H}_p$ is computed from $\boldsymbol {F}^{H(l)}_p$ using relation (B5) in Appendix B and $V_p\equiv {\rm \pi}D^{3}/6$ is the volume of a single spherical particle. Here, $\tilde {\phi }_p$ is the mean solid volume fraction defined in (B12). The error of this approximation is expected to be small as long as there is no sharp local gradient in $\tilde {\phi }_p$. Thus, the total plane-averaged shear stress is given by

(4.2)\begin{equation} {{\tau_{tot}}} = \underbrace{{{\rho_f}}\nu\partial_y {\langle \tilde{u}_f \rangle_{xzt}} }_{\tau_{\rm vis} = {{\rho_f}}\nu\partial_y U_f} -\underbrace{ \left( \underbrace{{{\rho_f}}{{\langle\tilde{u}^{\prime}_f\tilde{v}^{\prime}_f\rangle_{xzt}}}}_{\tau_{Rey}} + \underbrace{{{\rho_f}}{{\langle\langle\tilde{u}^{\prime\prime}_f\tilde{v}^{\prime\prime}_f\rangle_{zt}\rangle_{x}}}}_{{\tau_{form}}} \right)}_{ = {{\rho_f}} {{\langle u^{\prime}_f v^{\prime}_f\rangle_{xzt}}}} - \underbrace{\frac{1}{V_p} \int_y^{L_y}\langle \langle \tilde{f}^{H}_{p,x} \rangle_{zt} \tilde{\phi}_p \rangle_x \,{\rm d} y }_{{{\tau_{part}}} = \int_y^{L_y}\langle f_x \rangle \,{\rm d} y}. \end{equation}

Definitions of the velocity fluctuation covariances ${{\langle \tilde {u}^{\prime }_f\tilde {v}^{\prime }_f\rangle _{xzt}}}$ and $\langle \tilde {u}^{\prime \prime }_f\tilde {v}^{\prime \prime }_f\rangle _{zt}$ are given in Appendix B. Figure 6 shows the wall-normal profiles of the different contributions in (4.2). It can be seen that the actual total shear stress is sufficiently recovered with this approach. A slight deviation from the linear variation is observable, in regions where the above introduced errors are expected to be non-negligible (in the vicinity of the fluid–bed interface location). We remark that for the purpose of verifying the convergence of statistics we actually evaluate the forcing term $\langle f_x\rangle _{xzt}$ during run time. As shown in Appendix C, ${{\tau _{tot}}}$ is indeed varying linearly throughout the channel height. Moreover, figure 6 shows that the form-induced stress contribution is of opposite sign to that of the turbulent Reynolds stress, except in the outer fluid region where it exhibits small positive values and subsequently tends to zero for increasing $y$. This means that there is high momentum fluid that is carried away from the bottom into the outer flow and vice versa as a result of the form-induced stress. However, the magnitude of the form-induced stress is small when compared with that of the Reynolds stress. The form-induced momentum transfer is explainable by observing the streamlines of the mean flow. Upstream of the crest, the mean flow increases in the positive streamwise direction and is characterised by a positive wall-normal velocity (first quadrant stress), while downstream of the crest, the flow decreases and is downward moving (third quadrant stress).

Figure 6. Wall-normal profiles of the different contributions to the plane-averaged total shear stress computed based on the phase-shifted statistics (4.2): (a) H4; (b) H6; (c) H7.

The analysis presented above provides an integral information on the total shear stress of the flow over a space and time varying sediment bed. However, it fails to provide information on the local boundary shear stress and its spatial variation with respect to the evolving bed. This quantity, and its relation to the local sediment flow rate, is highly relevant in sediment transport modelling and thus its accurate determination is important. In developed bedforms, the total boundary shear stress is the sum of the dune-form drag that results from the pressure difference between the regions upstream and downstream of the crest of a ripple, as well as the skin friction due to fluid viscous stress and fluid–particle interaction at the grain scale (see e.g. Yalin Reference Yalin1977; Nikora et al. Reference Nikora, McEwan, Mclean, Coleman, Pokrajac and Walters2007a). Usually in experiments, the form drag is estimated by integration of pressure measurements along a stationary bedform and the skin friction is estimated by assuming a logarithmic velocity profile below the lowest velocity measurement point and performing extrapolation in which a certain value for the hydrodynamic roughness height has to be prescribed (see e.g. McLean & Smith Reference McLean and Smith1986; Maddux, Mclean & Nelson Reference Maddux, Mclean and Nelson2003; Nikora et al. Reference Nikora, Mclean, Coleman, Pokrajac, Mcewan, Campbell, Aberle, Clunie and Koll2007b). Similarly, Charru & Franklin (Reference Charru and Franklin2012) estimate the boundary stress by extrapolation of the total shear-stress profile (corrected for streamline curvature) within the internal developing boundary layer to the surface of the bed. However, the degree of uncertainty inherent in such approaches is substantial, and in the present case, proved to be unreliable due to the strong degree of particle modulation therein. Instead, we have determined the local boundary shear stress by evaluating the momentum balance in a local control volume as described below.

First, we define quadrilateral control volumes $\mathcal {V}$ with a fixed shape and of unit spanwise extent, enclosed by a surface $\mathcal {S}$, an exemplary one of which is shown in figure 7. Here, $\mathcal {V}$ is moving at the constant mean ripple velocity ${{u_D}}$ in the positive streamwise direction and is demarcated by the locations of the points $A$, $B$, $C$ and $D$. The top boundary (face $BD$) coincides with the top boundary of the computational domain while faces $AB$ and $CD$ are aligned perpendicular to the streamwise direction. The streamwise extent of the control volume is chosen to be two particle diameters (i.e. $\varDelta _{V} = 2D$) so that the ripple geometry is sufficiently resolved. Points $A$ and $C$ coincide with the fluid–bed interface. The small interface curvature between $A$ and $C$ is ignored and thus face $AC$ is a straight line segment of constant slope $\tan (\alpha )$. Let us recall that the fluid–bed interface is defined based purely on a threshold of the solid volume fraction and in general is not exactly aligned with the streamlines of the mean flow. It is expected that there would be advective momentum transfer across face $AC$ as a result of the mean fluid/particle flux across it.

Figure 7. A control volume $\mathcal {V}$ (with unit spanwise extent), enclosed by a surface $\mathcal {S}$ over which the momentum equation is integrated for the determination of the boundary shear stress. $\mathcal {V}$ is moving at the constant mean ripple velocity ${{u_D}}$ in the positive streamwise direction. Note that $\mathcal {V}$ is a quadrilateral, i.e. the edge $AC$ is a straight line with a slope $\tan {(\alpha )}$. Here, ${{\boldsymbol {n}}}$ and $\boldsymbol {t}$ are the unit normal and tangential vectors at the edge $AC$, respectively and $\boldsymbol {i}$ and ${{\boldsymbol {j}}}$ are the unit vectors in the ${{\tilde {x}}}$ and $y$ directions, respectively.

The Reynolds phase-averaged momentum equation, integrated over $\mathcal {V}$ can be written as

(4.3)\begin{align} \oint_\mathcal{S} (\tilde{\boldsymbol{U}}_f \tilde{\boldsymbol{U}}_f)\boldsymbol{\cdot} {{\boldsymbol{n}}}\,{\rm d}S + \oint_{\mathcal{S}} \left(\langle{{\tilde{\boldsymbol{u}}}}^{\prime}_f {{\tilde{\boldsymbol{u}}}}^{\prime}_f \rangle_{zt} + \frac{1}{\rho_f}{{\langle \tilde{p} \rangle_{zt}}} \bar{\bar{\boldsymbol{I}}} - \langle \tilde{\bar{\bar{\boldsymbol{\tau}}}} \rangle_{zt} \right)\boldsymbol{\cdot} {{\boldsymbol{n}}}\,{\rm d}S - \int_{\mathcal{V}} \langle \tilde{{{\boldsymbol{f}}}} \rangle_{zt}\, {\rm d}V = \int_{\mathcal{V}} \langle \varPi\rangle_t \,{\rm d}V, \end{align}

where ${{\langle \tilde {p} \rangle _{zt}}} \bar {\bar {\boldsymbol {I}}}$ and $\langle \tilde {\bar {\bar {\boldsymbol {\tau }}}} \rangle _{zt}$ are the phase-averaged pressure and fluid shear-stress tensor, respectively. Also, $\langle \varPi \rangle _t$ is the driving pressure gradient term averaged in the steady ripple propagation interval. Note that the time rate-of-change term of the momentum equation vanishes in the moving frame of reference in the steady ripple propagation interval. We are interested in determining the total shear stress at the bottom face $AC$ by integrating the momentum balance (4.3) over $\mathcal {V}$. The total shear force acting on face $AC$ has contributions both from the fluid stresses as well as the particle–fluid interactions

(4.4)\begin{equation} \tilde{\boldsymbol{R}} = \int_{AC} (\langle \tilde{\bar{\bar{\boldsymbol{\tau}}}} \rangle_{zt} - \langle{{\tilde{\boldsymbol{u}}}}^{\prime}_f{{\tilde{\boldsymbol{u}}}}^{\prime}_f \rangle_{zt})\boldsymbol{\cdot} {{\boldsymbol{n}}}\,{\rm d}S + \int_\mathcal{V} \langle \tilde{{{\boldsymbol{f}}}} \rangle_{zt} \,{\rm d}V, \end{equation}

and the average total boundary shear stress $\tilde {\tau }_{b}$ (per unit spanwise width) is then equal to the component of $\tilde {\boldsymbol {R}}$ parallel to $AC$, divided by the length of the segment $AC$ ($\varDelta _V\cos {\alpha }$) viz.

(4.5)\begin{equation} \tilde{\tau}_{b} = \frac{1}{\varDelta_V\cos{\alpha}} (\tilde{\boldsymbol{R}} - (\tilde{\boldsymbol{R}}\boldsymbol{\cdot}{{\boldsymbol{n}}}){{\boldsymbol{n}}})\boldsymbol{\cdot}\boldsymbol{t}. \end{equation}

Here, ${{\boldsymbol {n}}} = \sin (\alpha ) \boldsymbol {i} - \cos (\alpha ) {{\boldsymbol {j}}}$ is the unit normal perpendicular to the interface while $\boldsymbol {t}$ is the corresponding tangential unit vector; $\boldsymbol {i}$ and ${{\boldsymbol {j}}}$ are the unit vectors aligned in the positive streamwise and wall-normal directions, respectively (cf. schematics in figure 7). The fluid–solid interaction term in (4.4) implicitly accounts for the total stress as a result of all particles within the control volume, whether moving as bedload or suspended.

The value of $\tilde {\tau }_{b}$ can be computed either by directly evaluating (4.5), or from the balance (4.3) where each term is integrated over the volume $\mathcal {V}$ and all bounding faces (assuming statistical stationarity). In practice, we have resorted to the latter approach as approximating the fluid–particle interaction term from the available integral Lagrangian particle forces introduces non-negligible errors in regions of sharp solid volume fraction gradients at the location of face $AC$ (see discussion above). Moreover, especially in experiments, the fluid–solid interaction term is not easily accessible. Thus the approach chosen here at the same time demonstrates a procedure for accurately retrieving the boundary shear stress from experimentally measurable quantities.

Figure 8(a) shows the streamwise variation of the boundary shear stress $\tilde {\tau }_{b}$ along the fluid–bed interface, computed with the above-described method. Profiles of $\tilde {\tau }_{b}$ for all cases are seen to exhibit a similar evolution along the ripple geometry: a sharp decrease downstream of the crest attaining negative values beneath the re-circulation region, with the location of the minimum value in the range ${{\tilde {x}}} \approx 0.14\lambda$–$0.16\lambda$, followed by a sharp increase in the following region up to ${{\tilde {x}}} \approx 0.45\lambda$. Subsequently, $\tilde {\tau }_{b}$ exhibits a quasi-plateau region with mild increase up to the location of the maximum at ${{\tilde {x}}} \approx 0.85\lambda$ followed by a mild decrease region in the interval $0.85\lambda$ up to the location of the crest. The location of the minimum $\tilde {\tau }_{b}$ is observed to be upstream of the bed's trough while that of the maximum of $\tilde {\tau }_{b}$ is upstream of the ripple crest (with a shift of approximately $15D$, $23D$ and $27D$ for cases H4, H6 and H7). The observed negative boundary shear stress downstream of the crest is expected due to the reversed flow in the re-circulation region. Note that, in correspondence to the almost non-existent re-circulation region of case H4 (cf. the streamline plots in figure 5), the region where $\tilde {\tau }_{b}$ of case H4 attains a negative value is very small. The observed shift of the maximum of $\tilde {\tau }_{b}$ with respect to the ripple crest, which is a consequence of the fluid inertia, is also in line with previous literature on flow over undulating boundaries (Benjamin Reference Benjamin1959; Jackson & Hunt Reference Jackson and Hunt1975; Zilker & Hanratty Reference Zilker and Hanratty1979; Charru et al. Reference Charru, Andreotti and Claudin2013). The variation of $\tilde {\tau }_{b}$ along the stoss side of the ripple is comparable to those reported by Charru & Franklin (Reference Charru and Franklin2012) for the flow over evolving barchan dunes, although those authors report the maximum of the shear stress to be located at the brink of the barchan, contrary to our findings. However, we remark that the variation between the maximum value of $\tilde {\tau }_{b}$ and the value at the crest is fairly small and is therefore difficult to capture subject to measurement and extrapolation technique uncertainties. We also remark that, in Charru & Franklin (Reference Charru and Franklin2012), the reported boundary shear stress is normalised by the unperturbed smooth wall friction velocity upstream of the barchan dune and thus the normalised values of $\tilde {\tau }_{b}$ are larger in magnitude when directly compared with our numerical data.

Figure 8. (a) Streamwise variation of the boundary shear stress $\tilde {\tau }_{b}/(\rho _fu_{\tau }^{2})$: blue line, H4; red line, H6; black line, H7. The global minima and maxima of $\tilde {\tau }_{b}$ are marked by the dot symbols. The mean ripple profile of each of the cases, along with the location of the crest (maximum of ${\tilde {H}_b}$) and trough (minimum of ${\tilde {H}_b}$) are shown at the top of the panel for reference. The colour code at the bottom of the panel and the corresponding vertical dashed lines demarcate different regions of the $\tilde {\tau }_{b}$ profile, and it is later used in colour coding the data points in figure 12 (for instance, the grey region shows the region where $\tilde {\tau }_{b}$ attains negative values). (b) Upstream phase shift, normalised by the particle diameter, of the first few dominant Fourier modes of $\tilde {\tau }_{b}$ relative to the phase of the corresponding modes of ${\tilde {H}_b}$ for the different cases (H4, —$\bullet$—, blue; H6, —$\bullet$—, red; H7, —$\bullet$—, black). The inset shows the phase shift in radians.

The shift of the extrema with respect to the ripple geometry discussed above does not provide a complete picture of the phase shift between the shear stress and topography, since the shapes of $\tilde {\tau }_{b}$ and ${\tilde {H}_b}$ curves differ visibly. Nevertheless, since $\tilde {\tau }_{b}$ and ${\tilde {H}_b}$ are periodic, a more robust ‘average’ phase shift can be estimated from the separation length at which the cross-correlation between the two has a maximum value. The average phase shift, computed by the latter method, reads approximately $16D$, $18D$ and $19D$ for cases H4, H6 and H7 (cf. table 3). That is, the shear stress is in advance of the topography by approximately $16$–$19$ particle diameters for all simulated cases, showing only few particle diameters of variation (in contrast to the observation regarding the locations of the maxima). Note that, due to the nonlinear nature of the system, the reported phase shift is an integral quantity and does not reflect the phase shift of the individual modes. To this end, we have computed the phase shift of the Fourier modes of $\tilde {\tau }_{b}$ relative to the phase of the corresponding modes of ${\tilde {H}_b}$. More precisely, the phases of $\tilde {\tau }_{b}$ and ${\tilde {H}_b}$ are computed based upon the real and imaginary parts of their respective Fourier coefficients and subsequently a relative phase shift $\Delta \varphi _j$ for mode $j$ is defined as the difference between the two. Figure 8(b) shows the mode-wise phase shift as a function of the wavenumber $\kappa _j$ of the first few significant Fourier modes. Interestingly, a dispersion among the different modes can be observed, $\Delta \varphi _j/\kappa _j$ exhibiting a decreasing trend with decreasing wavelength ($\Delta \varphi _j$ in radians shows the opposite trend, largest wavelength associated with smallest phase shift and vice versa). A further observation is that the phase-shift dispersion collapses fairly well for the three cases only being limited at the lower $\kappa$ end of the spectrum by the discreteness of the available numerical harmonics. This indicates that the phase shift depends purely on the imposed flow conditions and not on the state of the ripple evolution. The above findings corroborate previous works (Fourrière et al. Reference Fourrière, Claudin and Andreotti2010; Claudin, Charru & Andreotti Reference Claudin, Charru and Andreotti2011) that typically consider a single mode in their modelling efforts. The phase shift between the local boundary shear stress and the sediment-bed height is believed to be the main factor for destabilising a mobile sediment bed during the initial stages of ripple formation as well as during the subsequent nonlinear evolution process (Andreotti & Claudin Reference Andreotti and Claudin2013; Charru et al. Reference Charru, Andreotti and Claudin2013).

Table 3. Average phase advance in particle diameters of the boundary shear stress with respect to the ripple profile, $\bar {\varphi }_{\tilde {\tau }_d \tilde {H}_b}$, and with respect to the particle flow rate, $\bar {\varphi }_{\tilde {\tau }_d \tilde {q}_p}$. The values are computed from the separations at which the $\tilde {\tau }_{b}$–${\tilde {H}_b}$ and $\tilde {\tau }_{b}$–${\langle \tilde {q}_p\rangle _{zt}}$ cross-correlation functions are maximum, respectively.

Another way of looking at the phase advance of the shear stress is to distinguish between its component that is in quadrature and the one that is in phase with respect to the bed topology. In theoretical analysis, a sinusoidal bottom perturbation of small amplitude is usually considered. Then the flow over a bottom perturbation is analysed as a layered structure: an outer inertia dominated region in balance with the form-induced pressure gradient and where viscous effects are less important and an inner logarithmic region close to the boundary. The matching between these two in some intermediate region results in the phase advance of the shear stress with respect to the bottom undulation, the shear-stress maximum being located upstream of the crest (Charru et al. Reference Charru, Andreotti and Claudin2013). To understand the physical mechanisms, the Fourier coefficients $\hat {\tau }_{j}$ of the boundary shear stress $\tilde {\tau }_{b}$ are commonly expressed in terms of two hydrodynamical parameters $\mathcal {A}$ and $\mathcal {B}$ as (Charru et al. Reference Charru, Andreotti and Claudin2013; Duran Vinent et al. Reference Duran Vinent, Andreotti, Claudin and Winter2019)

(4.6)\begin{equation} \hat{\tau}_{j}^{+} \exp(-{\rm i}\varphi_{b,j}) = {\kappa}_j|\hat{A}_j|(\mathcal{A}_j + {\rm i}\mathcal{B}_j), \end{equation}

where $\hat {A}_j$ are the corresponding Fourier coefficients of the ripple morphology, $\varphi _{b,j}$ is the phase shift of mode $j$ and where ${\rm i}\equiv \sqrt {-1}$ is the imaginary unit. Note that, in the linear analysis, a single sinusoidal mode (with zero phase shift) is considered for the bottom perturbation. In our case, the interaction between the flow and the self-formed finite-size, asymmetric ripples is far from linear. The ripple morphology consists of a band of modes each with distinct phase shift $\varphi _{b,j}$. Fuller understanding requires rigorous analysis of the broadband modulation that is beyond the scope of the paper. Nevertheless, it is worthwhile to look at our data in the context of the linear theories that have been the basis for modelling sediment-bed evolution. Thus, $\hat {\tau }_{j}$ is multiplied by $\exp (-{\rm i}\varphi _{b,j})$. The parameters $\mathcal {A}$ and $\mathcal {B}$ exhibit different spatial structures and different physical meanings. In the linear regime, $\mathcal {A}$ is in-phase (symmetric) while $\mathcal {B}$ is in-quadrature (anti-symmetric) with respect to the bottom surface. A positive value of $\mathcal {B}$ is associated with particle erosion from the trough and deposition at crest while a positive value of $\mathcal {A}$ is associated with downstream migration of the perturbation (Charru et al. Reference Charru, Andreotti and Claudin2013). Figure 9 shows the shear-stress components $\mathcal {A}$ and $\mathcal {B}$ corresponding to the first few modes extracted from the DNS data as a function of the wavenumber scaled with the hydrodynamic roughness $k_0$ (cf. § 3). The data points collapse well for the three cases (except for the scatter in the data of case H4), corroborating our previous argument that the phase shift of individual modes is independent of the ripple evolution stage. Compared with theoretical linear predictions, the values of $\mathcal {A}$ and $\mathcal {B}$ fall in general within the same order of magnitude (Fourrière et al. Reference Fourrière, Claudin and Andreotti2010; Charru et al. Reference Charru, Andreotti and Claudin2013; Duran Vinent et al. Reference Duran Vinent, Andreotti, Claudin and Winter2019). However, we recall that the hydrodynamic roughness height $k_0$, which is an important ingredient of the above models, is inferred from the flow over stationary roughness. Moreover, the flow separation downstream of the ripple crests is known to alter the shear-stress phase advance (Charru et al. Reference Charru, Andreotti and Claudin2013). Therefore, a direct comparison of the data presented in figure 9 with the theoretical linear predictions is not meaningful.

Figure 9. (a) The symmetric real component $\mathcal {A}$ and (b) the anti-symmetric imaginary component $\mathcal {B}$ of the first few modes of the bottom shear stress (relative to the phase of the corresponding modes of the ripple morphology): blue circles, H4; red circles, H6; black circles, H7.

4.3. Relationship between the local sediment flux and the boundary shear stress

The phase-averaged total sediment flow rate per unit width can be evaluated from the phase-averaged particle statistics as

(4.7)\begin{equation} {\langle \tilde{q}_p\rangle_{zt}}(x) = \int_0^{L_y} \langle \tilde{q}^{2D}_{p,x}\rangle_{zt} \,{\rm d} y, \end{equation}

where $\langle \tilde {q}^{2D}_{p,x}\rangle _{zt} \equiv {{\tilde {U}_p}} \tilde {\phi }_p$ is the streamwise component of the phase-averaged two- dimensional particle flow rate per unit area. Note that ${\langle \tilde {q}_p\rangle _{zt}}$ incorporates all particle transport modes, bedload and suspended load. Figure 10(a) shows the streamwise variation of ${\langle \tilde {q}_p\rangle _{zt}}$. The values of ${\langle \tilde {q}_p\rangle _{zt}}$ are different among the cases, that of case H7 being the largest while that of case H4 being the smallest for most of the ripple profile, including the stoss side. This is expected due to the difference in the average bottom friction among the cases that is proportional to the ripple height (see the bottom friction evolution in figure 2). However, in the region $0.1\lambda \lesssim {{\tilde {x}}} \lesssim 0.35 \lambda$, the opposite trend is observed, i.e. ${\langle \tilde {q}_p\rangle _{zt}}$ of H4 is the largest. This can be explained by a combined effect of particle inertia and hindrance related to the smaller ripple amplitude and the very mild recirculation (the boundary layer of case H4 is not as disrupted as that of the other cases). For all cases, the particle flow rate exhibits a minimum value in the location very close to the trough and a maximum value very near the crest essentially in phase with the ripple geometry. Such an in-phase variation of the ripple morphology and the sediment flux is a consequence of mass conservation (under the condition that the bedforms are steadily migrating and are of constant shape) as described by the Exner equation (Charru et al. Reference Charru, Andreotti and Claudin2013). However, a closer look into the spectral decomposition of the phase shift (figure 10b) reveals that, although the dominant modes have essentially zero phase shift, there is a clear trend of (somewhat linear) phase-shift increase with increasing wavenumber, albeit only to within two particle diameters. The mode-wise resolved phase shift will be discussed in more detail below in relation to that of the bottom shear stress. In table 3 we report the average phase shift between $\tilde {\tau }_{b}$ and ${\langle \tilde {q}_p\rangle _{zt}}$ that measures approximately $18$–$19$ particle diameters, fairly independent of the ripple dimension.

Figure 10. Streamwise variation of the volumetric particle flow rate ${\langle \tilde {q}_p\rangle _{zt}}$, normalised by the inertial scale $q_{i} = U_{g} D$. The global minima and maxima are marked by the dot symbols. The mean ripple profile of each of the cases, along with the location of the crest (maximum of ${\tilde {H}_b}$) and trough (minimum of ${\tilde {H}_b}$) are shown at the top of the panel. (b) Upstream phase shift, normalised by particle diameter, of the first few dominant Fourier modes of ${\langle \tilde {q}_p\rangle _{zt}}$ relative to the phase of the corresponding modes of ${\tilde {H}_b}$ for the different cases. The inset shows the phase shift in radians. Colour coding is the same as for figure 8.

A key quantity in sediment transport modelling is the local relationship between $\tilde {\tau }_{b}$ and ${\langle \tilde {q}_p\rangle _{zt}}$. It has been a common practice to relate the two through algebraic semi-empirical formulae (e.g. Meyer-Peter & Müller Reference Meyer-Peter and Müller1948; Wong & Parker Reference Wong and Parker2006) that are popular in engineering applications. However, it is recognised that the sediment flux does not immediately adapt to local change in the shear stress (Sauermann et al. Reference Sauermann, Kroy and Herrmann2001; Charru Reference Charru2006; Fourrière et al. Reference Fourrière, Claudin and Andreotti2010; Claudin et al. Reference Claudin, Charru and Andreotti2011). Models such as that by Meyer-Peter & Müller (Reference Meyer-Peter and Müller1948) do not incorporate such a relaxation effect and thus have been shown to be not adequate to predict the spatial structure of the particle flow rate, in particular in the context of bedform formation and evolution. Comparing figures 8(a) and 10(a) it is evident that the two profiles do not follow the same trend, corroborating the relevance of a relaxation process. It is noteworthy that the location of the respective maxima/minima do not fall at the same location. Moreover, ${\langle \tilde {q}_p\rangle _{zt}}$ is observed to be entirely of positive value across the entire ripple profile, including in the re-circulation region where $\tilde {\tau }_{b}$ is negative, for all cases. This shows that the spatial structure of sediment transport cannot be entirely predicted by classical models where the particle flow rate is expressed in terms of local bottom shear stress alone. The non-local correlation between $\tilde {\tau }_{b}$ and ${\langle \tilde {q}_p\rangle _{zt}}$ can be scrutinised by analysing the amplitude and phase shift of the corresponding Fourier modes. Figure 11(a) presents the amplitude of the first few dominant Fourier coefficients of ${\langle \tilde {q}_p\rangle _{zt}}$ (normalised by the inertial scaling $q_{i} = U_gD$) as a function of the corresponding amplitude of the bottom shear-stress coefficients (expressed in terms of the non-dimensional local Shields number ${\tilde {\theta }_{b}} = \tilde {\tau }_{b}/((\rho _p -\rho _f)gD)$), while figure 11(b) shows the mode-wise phase shift between the two quantities (both in particle diameters and in radians). Interestingly, the correlation of the amplitudes for all cases sufficiently collapses and it is well represented by a power-law fit of the form

(4.8)\begin{equation} |\hat{q}_{p}|/q_{i} = \alpha |{{\hat{\theta}_{b}}}|^{\beta}, \end{equation}

with fit parameters $\alpha =4.63$ and $\beta = 1.5$. The phase shift in radians seems to exhibit no particular trend and is fairly independent of the wavenumber, in particular at the lower wavenumber end of the spectrum (subject to the scatter of the data). For our considered data point, the average phase shift reads $\Delta \varphi _{av} \approx 0.336{\rm \pi}$ (the dashed line in that figure). This translates to a wavenumber-dependent phase shift in particle diameters $\Delta \varphi _j/D = \Delta \varphi _{av}/(\kappa _jD)\approx 27/\kappa _j$. The above observation motivates us to write a simplified relationship between the bottom shear stress and particle flow rate, which incorporates the phase shift between the two, as follows:

(4.9)\begin{equation} q_p = \sum_{j=0}^{modes}\alpha|{{\hat{\theta}_{b}}}|^{\beta-1} {{\hat{\theta}_{b}}}\exp{(-{\rm i}\Delta \varphi_{\rm av})}\exp{({\rm i}\kappa_j{{\tilde{x}}})}, \end{equation}

where $\varphi _{av}$ is the average mode-wise phase shift for modes $j>0$ (there cannot be a phase shift in the zeroth-mode). We remark that (4.9) is obtained from a fit to the ${\langle \tilde {q}_p\rangle _{zt}}$–${\tilde {\theta }_{b}}$ relationship that has resulted from a nonlinear process. Moreover, an average phase shift (in radians) is considered that ignores possible mode-wise dispersion. However, it can be argued that inclusion of the phase shift between the shear stress and sediment flux that encodes the relaxation of the sediment bed (which will be discussed below) is superior when compared with models based on a local relationship. Let us mention that aeolian sand transport models that account for a similar phase shift between topology and shear stress (link with the sediment flux through the Exner equation) have been shown to provide better prediction of the formation and evolution of dunes (Andreotti et al. Reference Andreotti, Claudin and Douady2002; Kroy, Sauermann & Herrmann Reference Kroy, Sauermann and Herrmann2002; Hersen et al. Reference Hersen, Andersen, Elbelrhiti, Andreotti, Claudin and Douady2004).

Figure 11. (a) Amplitude of the Fourier coefficients of the first few modes of ${\langle \tilde {q}_p\rangle _{zt}}$ as a function of the corresponding amplitudes of non-dimensional local Shields number. The dashed line is a power-law function given by $|\hat {q}_{p}|/q_{i} = 4.63|{{\hat {\theta }_{b}}}|^{1.5}$ that best fits the data: blue circles, H4; red circles, H6; black circles, H7. (b) Upstream phase shift of the boundary shear stress relative to the phase of particle flow rate. The dashed line shows an average phase shift of the dominant modes (averaged over all cases). Colour coding is the same as for figure 8.

Figure 12 shows the variation of the non-dimensional sediment flow rate resolved along the ripple profile as a function of the non-dimensional local Shields number ${\tilde {\theta }_{b}}$. We note that, in plotting the data points in the figure, ${\langle \tilde {q}_p\rangle _{zt}}$ is averaged over bins of width $\varDelta _V$ that match those used to compute $\tilde {\tau }_{b}$. The above-discussed complexity of the relationship is evident, which classical models such as that by Meyer-Peter & Müller (Reference Meyer-Peter and Müller1948) fail to predict. It is only for the blue coloured symbols (corresponding to the blue demarcated region in figure 8 roughly upstream of the maximum shear-stress location) that the particle flow rate predicted by MPM formula approximately matches the data. Note that, as has been reported in KU2017 and is reproduced in figure 12(b), the space- and time-averaged particle flow rate is nevertheless well predicted by the Wong & Parker (Reference Wong and Parker2006) version of the Meyer-Peter & Müller (Reference Meyer-Peter and Müller1948) power-law formula. On the other hand, as demonstrated in the figure, the scale-resolved, phase-shift accounting model (4.9), evaluated using ten modes, provides a reasonable representation of the data.

Figure 12. (a) Particle flow rate ${\langle \tilde {q}_p\rangle _{zt}}$ as a function of the local phase-averaged Shields number $\tilde {\theta }$: $\vartriangle$, H4; $\circ$, H6; $\square$, H7. Marker colours indicate the different locations along the ripple profile (cf. colour coding in figure 8). The space- and time-averaged mean particle flow rate ${{\langle q_p\rangle _{zxt}}}$ for the three cases are shown by the symbols $\otimes$. The solid curves correspond to the model (4.9), evaluated using the first ten Fourier modes; blue line, H4; red line, H6; black line, H7. (b) The same data as in panel (a) but plotted logarithmically. The dashed line is the Wong & Parker (Reference Wong and Parker2006) version of the Meyer-Peter and Müller (MPM) formula for turbulent flows: $q_p/q_{i}=4.93\ (\theta -\theta _c)^{1.6}$.

For a closer inspection of the sediment-bed response to the spatial changes in the shear stress, it is essential to distinguish the particle flux within the mobile transport layer from possible contributions due to suspended particles, in particular in the re-circulation region downstream of the crest. Moreover, although minute in magnitude, the motion deep within the bulk sediment bed is not directly related to the boundary shear stress but rather to the form-induced pressure gradient from the outer flow (Mazzuoli et al. Reference Mazzuoli, Kidanemariam and Uhlmann2019). To this end, we have decomposed $\langle \tilde {q}^{2D}_{p,x}\rangle _{zt}$ into four transport mode regions (cf. figure 13). First we distinguish between the ‘resting’ sediment bed and the active mobile sediment layer – usually termed in the community as the ‘transport layer’ – by defining an interface where the phase-averaged solid volume fraction $\tilde {\phi }_p$ equals the mean value in the resting bed ${{\varPhi _{bed}}}$ (the interface between the orange and red regions in figure 13b). More precisely, for a streamwise position ${{\tilde {x}}}$, we fit a straight line to $\tilde {\phi }_p$ in the interval between the wall-normal locations where $\tilde {\phi }_p=0.1$ and $\tilde {\phi }_p=0.4$. Then, we pick the wall-normal location $y_b({{\tilde {x}}})$ where the fitted line crosses the value ${{\varPhi _{bed}}}$. The motivation for the above definition stems from the solid volume-fraction-dependent rheology of the sheared sediment layer in which $\tilde {\phi }_p$ tends to ${{\varPhi _{bed}}}$ as the local shear rate tends to zero (Boyer, Guazzelli & Pouliquen Reference Boyer, Guazzelli and Pouliquen2011). We subsequently demarcate the transport layer from the dilute suspension region by the interface located at $y_t({{\tilde {x}}})$ . In the latter, sediments are transported mainly in suspended mode (the blue region in figure 13b). The interface is defined based on a threshold value of the phase-averaged wall-normal inter-particle collision force ${{\langle \tilde {f}^{C}_{p,y}\rangle _{zt}}}$ (Scherer et al. Reference Scherer, Kidanemariam and Uhlmann2020). Additionally, we distinguish the sediment flux immediately downstream of the ripple crest that is dominated by fully suspended particles that have retained a substantial amount of their momentum from upstream. This ’inertial sediment flux’ region is demarcated based on a threshold value of the phase-averaged streamwise hydrodynamic force acting on the particles $\langle \tilde {f}^{H}_{p,x} \rangle _{zt}$ as particles are characterised by increased negative $\langle \tilde {f}^{H}_{p,x} \rangle _{zt}$ due to fluid-drag-induced deceleration. The inertial region is shown by the green region in figure 13(b). Defining the inertial region by an indicator function $I_{r}$ that has a value of unity in the region and zero elsewhere, the particle flow rate within the transport layer and excluding the inertial region can be defined as

(4.10)\begin{equation} {{{\langle \tilde{q}_p^{TL}\rangle_{zt}}}}({{\tilde{x}}}) = \int_{y_{b}({{\tilde{x}}})}^{y_{t}({{\tilde{x}}})}\langle \tilde{q}^{2D}_{p,x}\rangle_{zt}(1-I_{r})\,{\rm d} y. \end{equation}

The other contributions to ${\langle \tilde {q}_p\rangle _{zt}}$ are analogously defined. As can be seen from figure 13(c), ${\langle \tilde {q}_p\rangle _{zt}}$ is dominated by ${{{\langle \tilde {q}_p^{TL}\rangle _{zt}}}}$ in almost the entire ripple profile, except in the region immediately downstream of the crest. It is also noteworthy that ${{{\langle \tilde {q}_p^{TL}\rangle _{zt}}}}$ attains very small but negative values in the recirculation region.

Figure 13. (a) Phase-averaged two-dimensional particle flow rate (streamwise component) per unit area, $\langle \tilde {q}^{2D}_{p,x}\rangle _{zt}$, normalised by $U_gD^{2}$ for case H7. The pink dashed-dot curves indicate the boundaries of the bedload transport layer. The green solid curve demarcates the particle inertia dominated region. The fluid–bed interface is indicated by the black dashed curve. (b) Contour lines of $\langle \tilde {q}^{2D}_{p,x}\rangle _{zt}/U_gD^{2}$ in the region downstream of the ripple crest (the box region shown in panel (a)). The different sediment particle transport regions are colour coded as follows: orange, quasi-stationary sediment bed; maroon, bedload transport layer; blue, dilute suspension region; green, particle inertia dominated region. (c) ${\langle \tilde {q}_p\rangle _{zt}}/q_i$ of case H7 (thick black line) decomposed into the different sediment flow rate components. Line colour coding is the same as for panel (b).

In figure 14, the relationship between the local Shields number and local ${{{\langle \tilde {q}_p^{TL}\rangle _{zt}}}}$ is presented. Only data points corresponding to the stoss side of the ripple (defined as the region between the trough and crest) are shown. The relaxation of ${{{\langle \tilde {q}_p^{TL}\rangle _{zt}}}}$ to changes in ${\tilde {\theta }_{b}}$ is evident from the figure: in the ‘shear-stress recovery’ region, i.e. the region where the shear stress increases sharply from its almost zero value downstream of the trough (the red interval in figure 8), ${{{\langle \tilde {q}_p^{TL}\rangle _{zt}}}}$ is seen to increase modestly exhibiting a linear relationship with ${\tilde {\theta }_{b}}$, which is smaller when compared with the $3/2$ power law that is expected to hold in equilibrium. Subsequently, while the shear-stress evolution reaches a plateau (green and blue regions in figure 8), the particle flux sharply increases and then tends towards a $3/2$ power-law relationship in the blue region.

Figure 14. Particle flow rate in the bedload transport layer ${{{\langle \tilde {q}_p^{TL}\rangle _{zt}}}}$ as a function of the local Shields number $\tilde {\theta }$: $\blacktriangle$, blue, H4; $\bullet$, red, H6; $\blacksquare$, H7. Only data points on the stoss side of the ripple are shown. The black dashed curve is the MPM formula, while the red dashed line shows a linear scaling.

4.4. Relaxation of the sediment flux

Although it is established that the phase lag between $\tilde {\tau }_{b}$ and ${{{\langle \tilde {q}_p^{TL}\rangle _{zt}}}}$ is a consequence of the inertia of the subaqueous sediment bed, there has been no direct measurement of the relaxation effect. In the following, we evaluate this quantity from our DNS data, independent of the influence of the bedforms. Let us recall that we have initiated the simulations by first developing a turbulent flow over a macroscopically flat bed composed of stationary particles. After the coupled fluid–solid simulations were switched on by allowing the particles to move, the erosion–deposition process ramps up over a relatively short time in response to the shear imposed by the flow (see detailed description of the startup procedure in Kidanemariam & Uhlmann Reference Kidanemariam and Uhlmann2014a). Our startup procedure allows us to investigate the sediment flux relaxation behaviour during the first instants of the sediment grain motion. To this end, we have computed the space-averaged instantaneous particle flow rate within the transport layer given as

(4.11)\begin{equation} {{\langle q_p^{TL} \rangle_{xz}}}(t) = \int_{y_b(t)}^{y_t(t)} \langle u_p \phi_p \rangle_{xz}\,{\rm d} y, \end{equation}

for the time interval before any bed patterns have formed. Here, $y_b$ and $y_t$ are the lower and upper extents of the transport layer and are defined in § 4.2. Note that, during this stage, the transients of all the simulation cases reported in KU2017 are statistically identical. Therefore, in order to increase the sample size, we have analysed the data corresponding to case H48 in KU2017. That case possesses a relatively large domain size ($L_x = 48H$) with over a million particles and is suitable to scrutinise the transient evolution of the particle flow rate accurately. As can be seen in figure 15(a), ${{\langle q_p^{TL} \rangle _{xz}}}$ requires a finite time to attain a saturated value in this ’featureless bed’ interval (once bedforms emerge the particle flow rate increases again in accordance with the temporal increase of the bottom friction). In order to compute a characteristic time for the transient behaviour we fit an exponential of the form

(4.12)\begin{equation} q = q_s(1 - \exp({-}t/T_q)), \end{equation}

with fit parameters $q_s$ and $T_q$, to $\langle q_p \rangle _{xz}$ in the interval $0 < t/T_b < 80$. Our analysis follows that by Andreotti, Claudin & Pouliquen (Reference Andreotti, Claudin and Pouliquen2010) and Pähtz et al. (Reference Pähtz, Omeradžić, Carneiro, Araújo and Herrmann2015) who have carried out experiments and DEM simulations of aeolian sand transport saturation length, respectively. Subsequently, we define a transient time scale $T^{0.99}_q$ as the time required for $q$ to reach $0.99q_s$. It turns out that, for case H48, $q_s \approx 0.08q_i$, $T_q \approx 5.4T_b$ and $T^{0.99}_q\approx 25T_b$ (cf. figure 15a). Furthermore, $T^{0.99}_q$ can be translated into an equivalent length scale $L^{0.99}_q$ as

(4.13)\begin{equation} L^{0.99}_q = \int_0^{T^{0.99}_q} u_p^{TL}\, {\rm d}t, \end{equation}

where $u_p^{TL}$ is the average instantaneous particle velocity within the bedload transport layer, given as

(4.14)\begin{equation} u_p^{TL}(t) = \frac{{{\langle q_p^{TL} \rangle_{xz}}}}{ \varPhi_p^{TL}} . \end{equation}

Here, $\varPhi _p^{TL}$ is the volume of particles in the transport layer per unit horizontal area, viz.

(4.15)\begin{equation} \varPhi_p^{TL}(t) = \int_{y_b(t)}^{y_t(t)} \langle \phi_p\rangle_{xz} \,{\rm d} y . \end{equation}

The evolution of $u_p^{TL}$ shown in figure 15(b) exhibits a similar relaxation trend as ${{\langle q_p^{TL} \rangle _{xz}}}$ albeit at a shorter time scale. We emphasise that the relaxation of $u_p^{TL}$ is mainly controlled by the particle acceleration/deceleration due to fluid drag and inter-particle interaction within the transport layer while that of ${{\langle q_p^{TL} \rangle _{xz}}}$ is controlled both by the relaxation of $u_p^{TL}$ and that of the mobile particle number density, or equivalently $\varPhi _p^{TL}$ (Charru Reference Charru2006; Pähtz et al. Reference Pähtz, Kok, Parteli and Herrmann2013, Reference Pähtz, Parteli, Kok and Herrmann2014). As shown in the inset of figure 15(b), $L_q^{0.99} \approx 20D$ and agrees very well with the computed average phase shift between the boundary shear stress and the particle flow rate over the bedforms, thus suggesting the link between the two. The relaxation of $u_p^{TL}$ by itself is an important model ingredient (cf. Pähtz et al. Reference Pähtz, Kok, Parteli and Herrmann2013, Reference Pähtz, Parteli, Kok and Herrmann2014, Reference Pähtz, Omeradžić, Carneiro, Araújo and Herrmann2015) and analogous to the particle flow rate, we fit the same exponential function defined with parameters $u_{p,s}$ and $T_u$ (or equivalently $T^{0.99}_u$) to the evolution of $u_p^{TL}$. Incidentally, $u_{p,s} \approx 0.49u_{\tau }$, $T_u \approx 2.8T_b$ and $T^{0.99}_u\approx 13T_b$. Furthermore, $T^{\rm 0.99}_u$ can be translated into an equivalent length scale $L^{0.99}_u$, and, as shown in the inset of figure 15(b), it turns out that $L_u^{0.99} \approx 10D$. To the best of our knowledge, there is no direct experimental measurement of the relaxation length/time scales of subaqueous sediment to compare our results with. The quantity is usually reconstructed from measured dune shape and migration velocity (cf. Fourrière et al. Reference Fourrière, Claudin and Andreotti2010; Franklin & Charru Reference Franklin and Charru2011).

Figure 15. (a) Time evolution of the particle flow rate within the transport layer of case H48 in KU2017. The red curve is an exponential fit to the data of the form $q = q_{s}(1 - \exp (-t/T_q))$. (b) Corresponding evolution of the mean particle velocity within the transport layer along with the exponential fit $u_p = u_{p,s}(1 - \exp (-t/T_u))$. The inset shows the equivalent distance $s$ travelled by the mobile layer as a function of time defined as $s(t) = \int _0^{t} u_p^{TL} \,{\rm d}t$. The green and blue dashed lines in the inset show $T^{0.99}_u$ and $T^{0.99}_q$ respectively, while the ordinates at which the dashed lines cross the curve $s(t)$ correspond to $L^{0.99}_u$ and $L^{0.99}_q$, respectively.