3.1. Reynolds and form-induced stresses

Figure 3 shows the bed-normal variation of all components of Reynolds and form-induced stresses for the PBL, PBM and PBH cases. The zero-displacement plane (see Appendix D for details) defines the SWI in this study and is used as a virtual origin while comparing and contrasting the primary and secondary statistics. The variables are normalized by $u_{\tau }^2$ (pressure by $\rho u_{\tau }^2$), and $y$ is shifted by $d$, the zero-displacement thickness, and then normalized by $\delta$, effectively making the virtual origin the same for all the cases. This implies that SWI planes for all cases align. It is observed that the magnitude of the form-induced stresses is generally smaller than the Reynolds stresses. The locations of the peaks are also different, with the form-induced stresses peaking below the SWI, whereas the Reynolds stresses peak close to the bed crest.

Figure 3. Comparison of Reynolds (lines) and form-induced (lines and symbols) stress profiles for PBL (green solid line, circles), PBM (blue dotted line, squares) and PBH (black dash-dotted line, triangles) cases. (a–c) Streamwise, bed-normal and spanwise components of the spatially averaged Reynolds stress tensor. (d) Spatially averaged shear Reynolds stress and shear form-induced stress. (e–g) Streamwise, bed-normal and spanwise components of form-induced stress. (h) Mean-square pressure fluctuations and form-induced pressure disturbances. Horizontal lines show the crest of sediment bed for PBL (gray solid line), PBM (gray dotted line) and PBH (gray dashed dotted line) cases. Pressure is normalized by $\rho u_{\tau }^2$.

The profiles of streamwise ($x$-direction) $\langle \overline {u^{\prime 2}}\rangle ^{+}$ and bed-normal ($y$-direction) $\langle \overline {v^{\prime 2}}\rangle ^{+}$ Reynolds stress (figures 3a,b) exhibit similarity for $(y+d)/\delta \geq 0.4$, substantiating the wall similarity hypothesis reported by Raupach, Antonia & Rajagopalan (Reference Raupach, Antonia and Rajagopalan1991), Breugem et al. (Reference Breugem, Boersma and Uittenbogaard2006) and Fang et al. (Reference Fang, Xu, He and Dey2018). Near the bed, the streamwise stress decreases, and bed-normal and spanwise ($z$-direction, figure 3c) stresses increase with increasing Reynolds number, suggesting weakening of the wall-blocking effect with increase in non-dimensional bed permeability. Moreover, not only does the flow penetrate deeper inside the bed (quantified in § 3.3) but the intensity also increases with increasing $Re_K$, as seen from the bed-normal component of the stress. However, this comes at the expense of loss in intensity in the streamwise direction. The peak values and their locations for the Reynolds stresses are shown in table 3. The peak values of $\langle \overline {u^{\prime 2}}\rangle ^{+}$ are larger than $\langle \overline {v^{\prime 2}}\rangle ^{+}$ for all Reynolds numbers, similar to turbulent flows over smooth or rough impermeable wall cases. The locations of peaks in $\langle \overline {u^{\prime 2}}\rangle ^{+}$ and $\langle \overline {w^{\prime 2}}\rangle ^{+}$ are close to the crest for all $Re_K$, whereas that in $\langle \overline {v^{\prime 2}}\rangle ^{+}$ shifts downwards, starting from above the crest level for the PBL case, and moving close to the crest for the PBH case. The Reynolds shear stress shown in figure 3(d) also peaks near the crest, increases with $Re_K$, and shows similarity in the outer region.

Table 3. The peak value and location ($(y+d)/\delta$), given in parentheses, of Reynolds stresses and form-induced stresses for the PBL, PBM and PBH cases. The peak values of the Reynolds and form-induced stresses are normalized by $u_{\tau }^2$ (pressure by $\rho u_{\tau }^2$).

The bed-normal variation of the form-induced stresses for the PBL, PBM and PBH cases are shown in figures 3(d–g). The influence of $Re_K$ is much more pronounced on the form-induced stresses. While the magnitudes for the streamwise stresses are comparable, the bed-normal and spanwise stresses are much smaller than the corresponding Reynolds stresses. More importantly, the peak values occur significantly below the sediment crest. Inside the bed, $\langle \tilde {\bar {u}}^2\rangle ^{+}$ decreases with $Re_K$, while both $\langle \tilde {\bar {v}}^2\rangle ^{+}$ and $\langle \tilde {\bar {w}}^2\rangle ^{+}$ increase with $Re_K$, mimicking the trend with Reynolds stresses. However, in contrast to $\langle \overline {v^{\prime 2}}\rangle ^{+}$, the peaks in $\langle \tilde {\bar {v}}^2\rangle ^{+}$ occur at a similar location (table 3) for all Reynolds numbers. This suggests that the penetration depth of $\langle \tilde {\bar {v}}^2\rangle ^{+}$ is independent of $Re_K$ and is dependent on the local porosity distribution in the top layer of the bed, which is similar for the three $Re_K$ cases studied in this work. It is interesting to note that the values of $\langle \tilde {\bar {w}}^2\rangle ^{+}$ are larger than $\langle \tilde {\bar {v}}^2\rangle ^{+}$ at lower $Re_K$, but are comparable and slightly smaller than $\langle \tilde {\bar {v}}^2\rangle ^{+}$ at higher $Re_K$. This may be attributed to the increased flow penetration, causing the bed-normal form-induced stresses peak further below the bed crest.

The turbulent fluctuations and form-induced disturbances in pressure, shown in figure 3(h), exhibit a very strong correlation with $Re_K$. Again, the form-induced pressure disturbances are generally much smaller than the turbulent fluctuations. There is more than a ten-fold increase in the peak value between the PBL and PBH cases for both the turbulent fluctuations and form-induced disturbances in pressure. The location of the peak in form-induced disturbances as well as turbulent fluctuations occurs just below the crest and is almost the same for all $Re_K$ studied. The fluctuations decay quickly for $(y+d)/\delta < -0.3$, indicating that a significant magnitude contribution comes from the top layer of the sediment bed. This suggests that the local protrusions of partially exposed sediment particles in the top layer and resultant stagnation flow are responsible for altering the flow structures that in turn produce larger form-induced disturbances and pressure fluctuations. The increased pressure disturbances due to the bed roughness elements can lead to enhanced mass transport rates at the SWIs, with potentially reduced residence times for pollutants and contaminants.

3.2. Turbulence structure and quadrant analysis

Distinct variations in the characteristics of primary turbulence structure are shown first in this subsection, followed by the quadrant analysis. Contours of instantaneous bed-normal vorticity, $\omega _{y}^{+}$ = $\omega _{y} \nu / u_{\tau }^2$, just above the crest ($y/\delta = 0.005$) are shown in figure 4. Results from a simulated smooth-wall case at $Re_{\tau } = 270$ (the same $Re_{\tau }$ as in the PBL case) are also shown in figure 4(a) for comparison. In the smooth-wall case, distinct long elongated streaky structures, which are a result of low- and high-speed streaks generated by quasi-streamwise vortices, are visible. The influence of a strong mean gradient and an impenetrable smooth wall results in these long streaky structures (Lee, Kim & Moin Reference Lee, Kim and Moin1990). In the low Reynolds number permeable bed case (PBL), shown in figure 4(b), the start of the breakdown of these structures can be seen. The roughness and permeability of the bed help in the breakdown. Although the long elongated streaks in the PBL case are shortened due to roughness, at this $Re_K$, the flow anisotropy is somewhat retained. With further increase in Reynolds number, the streaks are broken down even more, and in the PBH case (figure 4d), the streaky structures are significantly less pronounced, with the flow becoming more intermittent. As $Re_K$ increases, weakening of the wall-blocking effect due to strong bed-normal velocities prevents the formation of these long streaky structures and leads to reduction in flow anisotropy.

Figure 4. Contours of bed-normal vorticity $\omega _{y}^{+}$ normalized by $u_{\tau }^2/\nu$ at $y/\delta = 0.005$, for $0 < x/\delta <6.28$ and $0 < z/\delta <3.14$: (a) smooth-wall case simulated at $Re_{\tau } = 270$, (b) PBL case, (c) PBM case, and (d) PBH case.

Quadrant analysis is performed to understand the influence of near-bed flow structures on Reynolds stress. Joint p.d.f.s for turbulent velocity fluctuations $u^{\prime }$ and $v^{\prime }$ calculated at different elevations from the sediment bed are shown in figure 5. The product $u^{\prime } v^{\prime }$ is negative in the second and fourth quadrants, representing turbulence production. The second quadrant, where $u^{\prime } < 0$, $v^{\prime } > 0$, corresponds to an ejection event, whereas the fourth quadrant, where $u^{\prime } > 0$, $v^{\prime } < 0$, corresponds to a sweep event. Figures 5(a–l) show the correlations for different $Re_K$ at four bed-normal elevations within one particle diameter below and above the crest. Figures 5(a–c) are one diameter above the crest in the free-stream, figures 5(d–f) are at the bed crest, figures 5(g–i) are halfway into the top layer of the bed, and figures 5(j–l) are at the bottom of the top layer ($y/D_p = -1$).

Figure 5. Quadrant analysis of p.d.f.s of turbulent velocity fluctuations for (a,d,g,j) the PBL case, (b,e,h,k) the PBM case, and (c, f,i,l) the PBH case, at different elevations. Values of $y/\delta$ ($y/D_{p}$) are (a–c) 0.143 (0.5), (d–f) 0 (0), (g–i) $-$0.143 ($-0.5$), and (j–l) $-$0.286 ($-$1). The velocities have been normalized by $u_{\tau }$.

Inside the bed (figures 5j–l), at the bottom of the top layer, the probability for ejection and sweep events is small and nearly equal for all three $Re_K$ values. This suggests that the turbulent structures lose both their directional bias and strength, becoming nearly isotropic in nature below the top layer. At half the distance into the top layer (figures 5g–i), it can be seen that for the PBL and PBM cases, the ejection events are slightly more dominant, transporting fluid away from the bed; however, for the PBH case, both the ejection and sweep events are equally dominant, indicating that at higher Reynolds number, sweep events are enhanced, increasing transport of momentum towards the bed. At the bed crest (figures 5d–f), both ejection and sweep events become dominant over a greater range of velocity fluctuations, showcasing the interaction between fluid flow and permeable sediment bed. For lower $Re_K$, the joint p.d.f.s are more concentrated in a narrow band, indicating that large fluctuations in $u^{\prime }$ are associated with smaller excursions in $v^{\prime }$. With increase in Reynolds number, the p.d.f.s appear more diffused, highlighting increase in bed-normal velocity fluctuations, and reduction in anisotropy structures at the SWI, which is confirmed from the vorticity contours discussed earlier. Further away from the bed at $y/D_{p}=0.5$, as shown in figures 5(a–c), the ejection and sweep events again decrease in intensities compared to those at the SWI.

3.3. Penetration depths, mixing length, and similarity relations

Passage of turbulent structures over the sediment bed and resulting sweep events were shown to penetrate into the sediment bed in § 3.2. These sweep events induce momentum fluxes that can be of the order of the mean bed shear stress. The penetration of turbulence increases the flow resistance and effective roughness. The depth of turbulent shear penetration is associated with the characteristic size of the turbulent eddy across the SWI. Knowing how these scales are related to the permeability ($\sqrt {K}$) or the mean particle size ($D_p$) at different $Re_K$ is important for reduced-order models of turbulent momentum and mass transport across the interface.

The mean flow penetration depth $\delta _{b}$, known as the Brinkman layer thickness, is calculated from the mean velocity profiles below the sediment crest. Deep inside the bed, the mean velocity reaches a constant value (Darcy velocity), which is denoted as $U_{p}$. The Brinkman layer thickness is then calculated by measuring the the vertical distance from the SWI ($y = -d$) to a location inside the bed, where the difference between the local mean velocity ($\langle {\bar {u}}\rangle (y)$) and Darcy velocity ($U_{p}$) has decayed to 1 % of the velocity value at the SWI ($U_{i}$), i.e. $\langle {\bar {u}}\rangle (y)_{y+d=-\delta _{b}} = 0.01(U_{i} - U_{p}) + U_{p}$ (where subscript $i$ indicates the SWI location). The corresponding value of the Brinkman layer thickness measured from the crest of the bed is $\delta _b^{*} = \delta _b+d$. Figure 6(a) shows strong correlation of the normalized Brinkman layer thickness with permeability Reynolds number. Similar trends are observed in experimental data of Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017), despite the fact that the actual position and size of the particles in the random arrangement used in the present study are different from those in the experiments.

Figure 6. (a) Brinkman layer thickness $\delta _{b}$ normalized by $\sqrt {K}$; Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) indicated by blue diamonds, DNS indicated by black circles. (b) Turbulent shear stress penetration $\delta _{p}$, normalized by $\sqrt {K}$; Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) indicated by blue diamonds, DNS indicated by black circles. (c) Ratio of the flow penetration depth $\delta _{b}^*$ to the shear stress penetration depth $\delta _{p}^*$, measured relative to the crest; Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) indicated by blue diamonds, DNS indicated by black circles. (d) Left-hand axis: ratio of interfacial mixing length to the Kolmogorov length scale; Voermans et al. (Reference Voermans, Ghisalberti and Ivey2018) indicated by blue filled diamonds, DNS indicated by black circles. Right-hand axis: interfacial mixing length normalized by $\sqrt {K}$; DNS indicated by black squares.

The turbulent shear stress penetration is defined as the depth at which the Reynolds stress is 1 % of its value at the SWI, i.e. $\delta _p=\langle \overline {u^{\prime }v^{\prime }}\rangle _{y=-\delta _{p}} = 0.01\langle \overline {u^{\prime }v^{\prime }}\rangle _{i}$, and the value measured from the crest of the bed is $\delta _p^{*} = \delta _p+d$. Following the work of Ghisalberti (Reference Ghisalberti2009) on obstructed shear flows, Manes et al. (Reference Manes, Ridolfi and Katul2012) defined the penetration depth from the crest of the sediment ($\delta _p^{*}$), and showed that it is proportional to the drag length scale, i.e. $\delta _p^{*}\sim f(C_d a)^{-1}$, where $C_d$ is the drag coefficient of the medium, and $a$ is a length scale obtained based on the frontal area per unit volume of the solid medium. For monodispersed spheres, $a$ is proportional to the particle size, which in turn is related to the permeability. Thus, using a drag-force balance, Manes et al. (Reference Manes, Ridolfi and Katul2012) argued for a linear relation between $\delta _p^{*}$ and $\sqrt {K}$. However, figure 6(b) shows that the normalized $\delta _p$ is a function of $Re_K$. Both the mean flow ($\delta _b$) and turbulent shear stress penetration ($\delta _p$) show nonlinear correlation with the permeability, and increase with increasing $Re_K$. A deterministic relation is observed for the ratio $\delta _b^{*}/\delta _p^{*}$, with the permeability Reynolds number as the ratio approaches a constant value $1.1$, as shown in figure 6(c).

The dominant scale of the turbulent structures at the interface is affected by the interstitial spacing within the pore of which the permeability is a geometric measure, but it does not introduce any physical flow-dependent measure. Hence quantifying the dependence of the interfacial mixing length $\langle L_{m,i}\rangle = (\langle \overline {u^{\prime }v^{\prime }}\rangle _i / (\partial _y \langle \bar {u}\rangle )_i^2)^{1/2}$ on the permeability Reynolds number is important. The mixing length can be thought of as a representative length scale of the turbulent eddies at the SWI, responsible for turbulent transport of mass and momentum. It is of the order of the bed permeability and is greater than the Kolmogorov length scale ($L_k = (\nu ^3/\varepsilon )^{1/4}$, where $\varepsilon \approx u^2_{\tau }/\langle L_{m,i}\rangle$ is the dissipation rate of the turbulent kinetic energy; Tennekes & Lumley Reference Tennekes and Lumley1972). Figure 6(d) shows the interfacial mixing length $\langle L_{m,i}\rangle$ normalized by the Kolmogorov length scale ($L_{k}$) as well as normalized by the permeability ($\sqrt {K}$). The Brinkman layer thickness, the turbulent shear stress penetration and the mixing length at the interface show very similar dependence on $Re_K$ over the Reynolds numbers studied, suggesting that the mixing length is a relevant characteristic scale for transport of momentum and mass.

Flows over highly permeable boundaries share statistical similarities over different types of permeable geometries, as shown by Ghisalberti (Reference Ghisalberti2009). Asymptotic values predicted from present simulations for $\sigma _{v,c}/\sigma _{u,c} \sim 0.6$, $\sigma _{v,c}/u_{\tau } \sim 1.1$, $\sigma _{u,c}/u_{\tau } \sim 1.8$ and $U_c/u_{\tau } \sim 2.6$ match well with those observed by Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017). Here, $\sigma _{u,c}= \langle \overline {u^{\prime 2}}\rangle ^{1/2}_{c}$ and $\sigma _{v}= \langle \overline {v^{\prime 2}}\rangle ^{1/2}_{c}$, and $U_c$ is the mean velocity at the crest. The data are normalized by $-\langle \overline {u^{\prime }v^{\prime }}\rangle _c^{1/2}$, where the subscript $c$ indicates that these similarity relations are defined at the crest of the sediment bed, $y = 0$. The successful comparison of various turbulent quantities with the experimental work of Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017, Reference Voermans, Ghisalberti and Ivey2018) has important implications for flows over monodispersed sediment beds in general. Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) studied a range of $Re_K$ between 1 and 6.3 by varying the permeability, through the use of medium- and large-diameter particles, along with flow rates. In the present DNS, the particle sizes and permeability are kept constant, and $Re_K$ is varied by changing the bulk flow rate. Geometrical features of the sediment beds – such as size of particles, hence the permeability, local porosity variations, and tortuosity – between the experiments and the DNS are very different. The fact that the DNS predictions follow the experimental data closely suggests that for turbulent flow over randomly arranged monodispersed sediment beds, the actual location of sediment particles in the bed has little influence on the turbulence statistics. However, matching the $Re_K$ values is a necessary but not sufficient condition. A random as opposed to structured arrangement of monodispersed particles in the top layer is also important to achieve statistical similarity; this is investigated below.

3.4. Role of the top layer of the sediment bed

To quantify the influence of the top layer of the sediment bed on flow turbulence and statistics at the SWI, a rough impermeable wall case is investigated by matching the entire top layer of the permeable bed at medium Reynolds number (PBM). The top layer of the sediment is placed over a no-slip wall as shown in figure 2(d), corresponding to the case IWM-F in table 1.

Figures 7(a–h) compare the bed-normal variation of mean velocity, Reynolds and form-induced stresses, and pressure disturbances for the PBM and IWM-F cases. The mean velocity, Reynolds stress and turbulent pressure fluctuation profiles for both the cases overlap each other with no noticeable difference due to the presence of an underlying solid wall in IWM-F. The majority of high-magnitude bed-normal fluctuations are restricted to the top layer of the bed for both cases. The presence of a solid wall underneath the full layer of spherical roughness elements in IWM-F has minimal influence on both the magnitude and penetration of the mean flow and turbulent fluctuations. This is because the full layer of roughness elements creates pockets underneath, where the flow can penetrate. Since the turbulent kinetic energy within this layer is still small, the flow characteristics and momentum transport mechanisms resemble that of a permeable bed. Similar behaviour was observed and reported by Manes et al. (Reference Manes, Pokrajac, McEwan and Nikora2009) in their experimental work.

Figure 7. Comparison of the mean velocity, Reynolds (lines) and form-induced (lines and symbols) stress profiles for PBM (blue dotted line, squares) and IWM-F (red dash-dotted line, diamonds) cases. (a) Mean velocity. (b,c) Streamwise and bed-normal components of the spatially averaged Reynolds stress tensor. (d) Spatially averaged Reynolds stress and shear form-induced stress. (e–g) Streamwise, bed-normal and spanwise components of form-induced stresses. (h) Mean-square pressure fluctuations and form-induced pressure disturbances. The crest lines for the PBM and IWM-F cases overlap and are shown by the grey horizontal dotted and dashed lines, respectively. Pressure is normalized by $\rho u_{\tau }^2$.

As in the permeable bed cases, the form-induced stresses (figures 7e–g) have lower magnitudes compared to their corresponding Reynolds stresses, and the peak values occur significantly below the sediment crest even for IWM-F case. While the peak value for the streamwise component ($\langle \tilde {\bar {u}}^2\rangle ^{+}$) is well captured, the peaks in the bed-normal ($\langle \tilde {\bar {v}}^2\rangle ^{+}$) and spanwise ($\langle \tilde {\bar {w}}^2\rangle ^{+}$) components show differences between the two cases. Deeper penetration and higher magnitude of form-induced normal stresses are observed in PBM compared to the IWM-F case. The additional layers of sediment grains underneath the top layer in PBM provide connected pathways for the bed-induced flow disturbances to penetrate deeper with a gradual loss in intensity, while in IWM-F, the underlying solid wall blocks the flow penetration and redistributes the stresses tangentially into the spanwise direction, which is evident from the larger peak magnitude of $\langle \tilde {\bar {w}}^2\rangle ^{+}$ compared to $\langle \tilde {\bar {v}}^2\rangle ^{+}$.

The form-induced pressure disturbances ($\langle \tilde {\bar {p}}^2\rangle ^{+}$) shown in figure 7(h) penetrate deeper into the bed in PBM, resulting in a reduction of its peak value compared to the IWM-F case. The wall-blocking effect in IWM-F results in pressure disturbances extending above the crest level, up to about $({y+d})/{\delta } <0.3$. The turbulent pressure fluctuations, however, are well captured and are typically much larger than the form-induced pressure disturbances. The presence of the high-magnitude pressure fluctuations only in the top layer is an important observation provided by the present DNS data, as it showcases that inclusion of the effect of a single layer of roughness elements with a random arrangement for reach-scale hyporheic exchange models can potentially better capture the turbulent fluctuations.

3.5. Stress and force statistics on the particle bed

Direct measurements of shear stress or the drag and lift forces on sediment grains in a laboratory or in the field are challenging. The present pore-resolved simulations provide access to the spatio-temporal variations in these variables. Specifically, knowing higher-order statistics of bed shear stress as a function of Reynolds number can critically influence reduced-order models for mass and momentum transport that are based on the friction velocity ($u_{\tau }$). Furthermore, incipient motion, sliding, rolling and saltation, driving the bedload transport are modelled based on the bed shear stress exceeding a critical value. Higher-order statistics of bed shear stress are important in stochastic modelling of incipient motion (Ghodke & Apte Reference Ghodke and Apte2018b). Motion and rearrangement of sediments can alter local bed porosity and effective permeability, impacting hyporheic exchange directly. The DNS data are used to compute p.d.f.s and statistics of the local variation of net bed shear stress on particle surfaces as well as the net drag and lift forces on the particle bed at different $Re_K$.

The net stress ($\boldsymbol {\tau }^t = {\boldsymbol \tau }^v - p{\boldsymbol I}$) on the particle surface includes a contribution from the viscous stress (skin friction, ${\boldsymbol \tau }^v$) as well as the pressure stress (form, $p{\boldsymbol I}$). The normal and tangential components of the net stress can be evaluated directly on the particle surface from the velocity and pressure fields, and then transformed into the Cartesian frame using the surface normal vector to obtain the streamwise (${\tau }_x^t$), bed-normal (${\tau }_y^t$) and spanwise (${\tau }_z^t$) components, respectively. The local distribution of the net stress on the particle surface can be integrated over its surface area to obtain the net force on the particle,

(3.1)\begin{equation} {\boldsymbol F} = \int_{\varGamma} {\boldsymbol{\tau}}^t\boldsymbol{\cdot} {\boldsymbol n}\,{\rm d}\varGamma \equiv \int_{\varGamma} \left({\boldsymbol\tau}^v - p{\boldsymbol I}\right) \boldsymbol{\cdot} {\boldsymbol n}\,{\rm d}\varGamma, \end{equation}

where ${\rm d}\varGamma$ is the unit surface area of the particle, and $\boldsymbol n$ is the particle surface normal vector. Accordingly, statistics of the local distribution of the net stress as well as the drag and lift force components are evaluated.

3.5.1. Local distribution of net stress on the particle surface

The p.d.f.s of the streamwise ($\tau ^t_x$) and bed-normal ($\tau ^t_y$) net local stress normalized by the total stress in the $x$-direction integrated over all particles in the bed ($\tau _w^{t}$) are shown in figures 8(a,b). It is found that these distributions collapse nicely for all Reynolds numbers, suggesting that they are independent of $Re_K$. The dominance of positively skewed streamwise shear stress ($\tau _x^t$) events is clearly visible by this normalization. Positive skewness in the p.d.f.s is associated with the exposed protrusions of the spherical roughness elements to the free-stream, as instantaneous streamwise velocity generally increases in the bed-normal direction. For wall-bounded flows, the probability of negative wall shear stress is generally linked to the near-wall low- and high-speed regions; however, with protrusions from the rough sediment bed, these structures are destroyed, as was shown in § 3.2. The probability of negative $\tau _{x}^{t+}$ fluctuations is then highly associated with trough regions between the roughness elements, representative of reverse flow behind the exposed sediment particles. As the net bed stress increases with increase in $Re_K$, the probability of extreme streamwise stress events increases. The p.d.f.s of bed-normal ($\tau ^{t+}_y$) and spanwise ($\tau ^{t+}_z$, not shown) stress are more symmetric due to the absence of directional influence of a strong mean flow gradient.

Figure 8. The p.d.f.s of the net bed stress components: (a) streamwise ($\tau _{x}^t$), and (b) bed-normal ($\tau _{y}^t$). Superscript $+$ denotes a quantity normalized by the total bed stress in the $x$-direction ($\tau _w^{t}$). (c) Yaw angle $\psi _{\tau }$ based on the viscous components of the stresses. PBL indicated by green solid line, PBM indicated by blue dotted line, and PBH indicated by black dash-dotted line.

The relationship between the two viscous components of shear stress can be understood from their yaw angle $\psi _{\tau } = {\rm atan}(\tau _{z}^v(t)/\tau _{x}^v(t))$. Jeon et al. (Reference Jeon, Choi, Yoo and Moin1999) reported that the shear stress yaw angles for smooth walls are within the range $-45^{\circ }$ to $45^{\circ }$, indicating that events with large values of $\tau _{x}^v$ are associated with small $\tau _{z}^v$. However, in flow over rough permeable beds, the probability of yaw angles above $45^{\circ }$ is much higher, as seen in figure 8(c). This shows that comparable magnitudes of $\tau _{x}^v$ and $\tau _{z}^v$ occur more frequently in flows over sediment beds. However, the yaw angle is independent of $Re_K$, suggesting that it is influenced more by the roughness distribution of the bed rather than the flow. The roughness elements reduce the directional bias of near-bed vortex streaks, resulting in more isotropic vorticity fields, wherein the probability of large-scale fluctuations occurring simultaneously in both components of shear stress increases.

As the normalized p.d.f.s of the local distribution of the net bed stress collapse for different $Re_K$, it is conjectured that a simple deterministic fit to the p.d.f.s of fluctuations in stress is possible, which is investigated. The p.d.f.s of the local distribution of the net stress fluctuations on the particle surface normalized by their root-mean-square values are shown in figures 9(a–c). The higher-order statistics for the net stress are shown in table 4. The mean (not shown) and standard deviation of shear stresses increase with $Re_K$, suggesting higher probability of extreme events for larger $Re_K$. The fluctuation p.d.f.s are symmetric but non-Gaussian, and show peaky distribution with heavy tails and high kurtosis values given in table 4. A $t$-location model fit based on the variance, zero skewness and shape factor determined by excess kurtosis represents the distributions very well.

Figure 9. The p.d.f.s of net stress fluctuations normalized by their root-mean-square (rms) values for (a) streamwise ($\tau ^{t^\prime }_{x}$), (b) bed-normal ($\tau ^{t^\prime }_{y}$), and (c) spanwise ($\tau ^{t^\prime }_{z}$) components, for PBL (green solid line), PBM (blue dotted line), and PBH (black dash-dotted line) cases. Also shown is a non-Gaussian $t$-location scale fit (red circles).

Table 4. Higher-order statistics for streamwise ($\tau _x^{t^\prime }$), bed-normal ($\tau _y^{t^\prime }$) and spanwise ($\tau _z^{t^\prime }$) stress fluctuations, showing the standard deviation $\hat {\sigma }({\cdot })$, skewness $Sk({\cdot })$ and kurtosis $Ku({\cdot })$.

For smooth wall-bounded flows, the root-mean-square fluctuations of streamwise stress follow a logarithmic correlation as proposed by Örlü & Schlatter (Reference Örlü and Schlatter2011). Accordingly, a logarithmic correlation between the root-mean-square fluctuations and the friction Reynolds number is also assumed in the present permeable bed DNS:

(3.2)$$\begin{gather} \tau^{t+}_{x,rms} = \tau_{x,rms}^t/\tau_{w}^t =2.10\ln Re_{\tau} - 8.11, \end{gather}$$

(3.3)$$\begin{gather}\tau^{t+}_{y,rms} = \tau_{y,rms}^t/\tau_{w}^t =2.4\ln Re_{\tau} - 9.83, \end{gather}$$

(3.4)$$\begin{gather}\tau^{t+}_{z,rms} = \tau_{z,rms}^t/\tau_{w}^t = 2.4\ln Re_{\tau} - 10.10, \end{gather}$$

and is shown in figures 10(a–c). Here, the root-mean-square fluctuations are obtained by computing the net stress fluctuations on the particle surface over the entire bed, and then time averaging. For the present monodispersed bed, the friction ($Re_{\tau }$) and permeability ($Re_K$) Reynolds numbers are related to each other, thus the above relations can also be plotted in terms of the permeability Reynolds number as shown.

Figure 10. Variation of the root-mean-square fluctuations of the net bed stress normalized by the total bed stress ($\tau _w^t$) with the friction and permeability Reynolds numbers: (a) $\tau ^{+}_{x,rms}$, (b) $\tau ^{+}_{y,rms}$ and (c) $\tau ^{+}_{z,rms}$. Functions of $Re_{\tau }$: DNS data indicated by open black circles, fitted logarithmic correlation from DNS data for permeable sediment beds indicated by black dashed line. Functions of $Re_{K}$: DNS data indicated by solid black circles, fitted logarithmic correlation from DNS data for permeable sediment beds indicated by black dash-dotted line.

The logarithmic dependence of shear stress fluctuations, together with a symmetric, non-Gaussian distribution for local stress fluctuations, is an important result, as in the field measurements for natural stream or river bed studies, typically the friction velocity $u_\tau$, is measured (Jackson et al. Reference Jackson, Haggerty and Apte2013a, Reference Jackson, Apte, Haggerty and Budwig2015) to compute $Re_{\tau }$. Equations (3.2)–(3.4) can then be used to evaluate the bed stress variability for different $Re_K$. Together with the non-Gaussian model distribution for the p.d.f.s of stress fluctuations, a stochastic approach for mobilization and incipient motion of sediment grain can be developed.

3.5.2. Net drag and lift force distribution

The local stress distributions on the particle surface are integrated over each individual particle and then time averaged to obtain the mean and fluctuating drag and lift forces. Percentage contributions of the viscous and pressure stresses to the average drag and lift forces in the full bed as well as just the top layer of the bed are given in table 5. The majority of the contribution to the drag and lift forces comes from the pressure distribution. It was also found that the top layer of the bed results in average forces that are 3–4 times the average in the full bed, indicating that a significant contribution to the lift and drag force comes from the top layer of the bed.

Table 5. Drag and lift force statistics for the full bed and the top layer, showing the percentage contributions of viscous ($F_d^v$, $F_{\ell }^p$) and pressure ($F_d^p$, $F_{\ell }^p$) components to the mean force, and the standard deviation ($\sigma$), skewness ($Sk$) and kurtosis ($Ku$) of the fluctuations in the force.

The p.d.f.s of the fluctuations of drag ($F^{\prime }_d$, $x$-component) and lift ($F^{\prime }_{\ell }$, $y$-component) forces on all particles within the bed normalized by their respective standard deviation values are shown in figures 11(a,b). Similar to the local stress distributions, the drag and lift force fluctuations also exhibit a symmetric, non-Gaussian distribution with heavy tails. Higher-order statistics of the forces (see table 5) indicate minimal skewness and very high kurtosis, suggesting probability of extreme forces due to turbulence. A model $t$-location fit function based on the variance and excess kurtosis data fits well for all Reynolds numbers.

Figure 11. The p.d.f.s of (a,c) drag force and (b,d) lift force fluctuations, normalized by their respective standard deviation values. (a,b) Permeable bed results for PBL (green solid line), PBM (blue dotted line), and PBH (black dash-dotted line) cases. Also shown is a $t$-location scale model fit (red circles). (c,d) Comparison of the values for the full permeable bed (PBM, blue dotted line) and the top layer of the permeable wall (PBMTL, cyan solid line) for the medium $Re_K$.

Finally, to assess the contribution of the top layer to force statistics, figures 11(c,d) compare the p.d.f.s of the drag and lift forces for PBM full bed as well as only the top layer of the bed, referred to as the PBMTL data. Note that for the PBMTL data, force fluctuations are normalized by the total force only in the top layer. A close match suggests that the top layer of the bed contributes to the majority of the net drag and lift force fluctuations in the bed. This has implications for large-scale Reynolds-averaged Navier–Stokes modelling where a low-Reynolds-number model is used to estimate shear stress on the bottom solid wall of the domain. Including roughness effects through a single layer of sediments may help to improve prediction of reduced-order models.

3.6. Pressure distributions at the SWI

Pressure fluctuations at the SWI play a critical role in hyporheic transport even for a flat bed. Specifically, pressure fluctuations due to turbulence are conjectured to have significant impact on mass transport within the hyporheic zone as they can influence directly the residence times through turbulent advection. Reach-scale modelling of hyporheic zone transport typically uses a one-way coupling approach, wherein the pressure fields obtained from the streamflow calculations are used as boundary conditions for a separate mass-transport computation within the hyporheic zone using Darcy-flow-like models (Chen et al. Reference Chen2022). These studies have shown that better characterization of the pressure distributions at the SWI can have a significant impact on predicting transport. Specifically, using the present DNS data, the variation of pressure at the SWI with Reynolds number is quantified.

The p.d.f.s of pressure fluctuations and disturbances normalized by their respective standard deviations, and averaged over multiple flow through times for the PBL, PBM and PBH cases, are compared at their respective zero-displacement planes ($y=-d$) or the SWI. Figures 12(a,b) show the p.d.f.s for the normalized turbulent pressure fluctuations $p^{\prime }$ and normalized form-induced pressure disturbances $\tilde {\bar {p}}$. The turbulent fluctuations exhibit close to a normal distribution, whereas the form-induced pressure disturbances have skewed distributions with longer positive tails. This is attributed to the roughness protrusions that create positive pressure stagnation regions. Figure 12(c) shows the p.d.f.s of the sum of the normalized turbulent and form-induced ($p^{\prime } + \tilde {\bar {p}}$) pressure for the three cases. The pressure sum p.d.f.s for all cases are statistically similar and symmetric, slightly heavier-tailed than a Gaussian; however, the Gaussian distribution nearly captures the pressure data within $\pm 3 \hat {\sigma }$. This result suggests that the pressure behaviour inside the bed can be approximated with a simpler Gaussian distribution across a range of permeability Reynolds numbers typical of natural stream or river beds. Table 6 lists all higher-order statistics for turbulent fluctuations, form-induced disturbances and their sum. The skewness and kurtosis values for $\tilde {\bar {p}}$ are higher than for $p^{\prime }$, in alignment with the skewed distribution. The mean and variance values for the two are roughly of the same order at the SWI. However, the peak in turbulent pressure fluctuations is much larger than the peak in the form-induced disturbances, as seen from the bed-normal variations shown earlier, in figure 3(h).

Figure 12. The p.d.f.s of (a,d) turbulent pressure fluctuations ($p^{\prime }$), (b,e) form-induced pressure disturbances ($\tilde {\bar {p}}$), and (c, f) sum of $p^{\prime }$ and $\tilde {\bar {p}}$. (a–c) Data for the permeable bed cases PBL (green solid line), PBM (blue dotted line) and PBH (black dash-dotted line), and a Gaussian model fit (red circles). (d–f) Comparison of the p.d.f.s for the full permeable bed (PBM, blue dotted line) and the impermeable rough wall (IWF-M, red dash-dotted line).

Table 6. Higher-order statistics for turbulent fluctuations and form-induced pressure disturbances, showing the standard deviation $\hat {\sigma }({\cdot })$, skewness $Sk({\cdot })$ and kurtosis $Ku({\cdot })$.

The p.d.f.s of normalized pressure fluctuations and disturbances for the PBM and IWM-F cases are also compared at their respective zero-displacement planes ($y=-d$) or SWI in figure 12(d–f). The variance in turbulent pressure fluctuations and form-induced disturbances is generally larger in the IWM-F case than in the permeable bed case (PBM). The probability of higher negative $\tilde {\bar {p}}$ increases in the IWM-F case due to the blocking effect of the impermeable wall. However, the sum of the normalized distributions of $p^{\prime }$ and $\tilde {\bar {p}}$ shows a reasonable match between the PBM and IWM-F cases, as shown in figure 12( f). Therefore, the distribution of the sum of normalized pressure fluctuation and disturbances in the top layer of the sediment bed is found to be sufficient to capture the pressure variation at the SWI.