Pore-resolved investigation of turbulent open channel flow over a randomly packed permeable sediment bed

Abstract

Pore-resolved direct numerical simulations are performed to investigate the interactions between streamflow turbulence and groundwater flow through a randomly packed porous sediment bed for three permeability Reynolds numbers, $Re_K=2.56$, 5.17 and 8.94, representative of natural stream or river systems. Time–space averaging is used to quantify the Reynolds stress, form-induced stress, mean flow and shear penetration depths, and mixing length at the sediment–water interface (SWI). The mean flow and shear penetration depths increase with $Re_K$ and are found to be nonlinear functions of non-dimensional permeability. The peaks and significant values of the Reynolds stresses, form-induced stresses, and pressure variations are shown to occur in the top layer of the bed, which is also confirmed by conducting simulations of just the top layer as roughness elements over an impermeable wall. The probability distribution functions (p.d.f.s) of normalized local bed stress are found to collapse for all Reynolds numbers, and their root-mean-square fluctuations are assumed to follow logarithmic correlations. The fluctuations in local bed stress and resultant drag and lift forces on sediment grains are mainly a result of the top layer; their p.d.f.s are symmetric with heavy tails, and can be well represented by a non-Gaussian model fit. The bed stress statistics and the pressure data at the SWI potentially can be used in providing better boundary conditions in modelling of incipient motion and reach-scale transport in the hyporheic zone.

1. Introduction

The interchange of mass and momentum between streamflow and groundwater occurs across the sediment–water interface (SWI) and into the porous bed underneath, termed the hyporheic zone. Hyporheic transient storage or retention, and transport of solutes such as chemicals and pollutants, dissolved oxygen, nutrients and heat across the SWI, is one of the most important concepts for stream ecology, and has enormous societal value in predicting sources of fresh drinking water, transport, biogeochemical processing of nutrients, and sustaining diverse aquatic ecosystems (Bencala et al. 1983; D'angelo et al. 1993; Harvey, Wagner & Bencala 1996; Valett et al. 1996; Anderson et al. 2008; Briggs et al. 2009; Grant, Gomez-Velez & Ghisalberti 2018).

A broad range of spatio-temporal scales corresponding to disparate physical and chemical processes contribute to the mixing within the hyporheic zone (Hester et al. 2017). Turbulent transport across the SWI, coherent flow structures, and non-Darcy flow within the sediment bed have been hypothesized as critical mechanisms impacting transient storage. The importance of penetration of turbulence within the bed and near-bed pressure fluctuations are crucial, and their impact on the hyporheic transient storage is poorly understood (Hester et al. 2017). Moreover, turbulence characteristics over a permeable bed are different compared to a rough, impermeable wall (Jiménez 2004), impacting long time scales of retention within the bed.

Mass and momentum transport in turbulent flow over a naturally occurring permeable sediment bed are characterized by bed permeability $K$ (which depends upon its porosity and average grain size according to the Carman–Kozeny relation), sediment bed arrangement (flat versus complex bedforms), and friction or shear velocity $u_{\tau }$ (based on (2.3)). Permeability Reynolds number $Re_K = u_{\tau } \sqrt {K}/\nu$ (where $\nu$ is the kinematic viscosity), representing the ratio between the permeability scale ($\sqrt {K}$) and the viscous scale (${\nu /u_{\tau }}$), is typically used for granular beds to identify different flow regimes based on the dominant transport mechanisms across the SWI. For flat beds, three flow regimes (see figure 1) characterized by Voermans, Ghisalberti & Ivey (2017, 2018) and Grant et al. (2018) have been identified: (i) the molecular regime, $Re_K<0.01$, where the bed is nearly impermeable and the transport is governed by molecular diffusion; (ii) the dispersive regime, $0.01< Re_K<1$, where dispersive transport associated with laminarization of stream turbulence is important; and (iii) the turbulent regime, $Re_K>1$, where turbulence is dominant near the highly permeable interface. Based on the data collected for flat beds from several local streams and rivers near Oregon State University by coworkers (Jackson, Haggerty & Apte 2013a; Jackson et al. 2015), it is found that the gravel grain sizes varied over the range 5–70 mm, with friction velocity $u_{\tau }=0.004\unicode{x2013}0.088\,{\rm m}\,{\rm s}^{-1}$, and $Re_K=2\unicode{x2013}70$ (figure 1). Free surface flow and waves typically do not affect the hyporheic exchange in natural streams and rivers under subcritical conditions with small Froude numbers.

Figure 1. Effective dispersion coefficient versus $Re_{K}$ (modified based on Voermans et al. 2017; Grant et al. 2018).

Few experimental studies have evaluated turbulence characteristics over flat permeable beds (Zagni & Smith 1976; Zippe & Graf 1983; Manes et al. 2009; Suga et al. 2010; Manes, Poggi & Ridolfi 2011; Voermans et al. 2017; Kim et al. 2020; Rousseau & Ancey 2022). Bed permeability was found to increase friction coefficient, reduce the wall-blocking effect due to impermeable rough walls, and reduce near-bed anisotropy in turbulence intensities. Manes et al. (2009) studied turbulent flow over uniform cubic patterns of single and multiple layers of spheres at $Re_K=31.2$ and 44.6. Permeability was shown to influence flow resistance dramatically, and the conventional assumption of a hydraulically rough regime, wherein the friction factor is dependent upon the relative submergence or the ratio of roughness size to flow thickness, does not apply to permeable beds. Friction factor was shown to increase progressively with increasing Reynolds number.

Voermans et al. (2017) studied the influence of different $Re_K$ on the interaction between surface and subsurface flows at the SWI of a synthetic sediment bed composed of randomly arranged monodispersed spheres using refractive-index-matched particle tracking velocimetry. Their experiments covered a wide range, $Re_K=0.36$–$6.3$, and varied the permeability of the beds by investigating three different sphere sizes. The results demonstrated a strong relationship between the structure of the mean and turbulent flow at the SWI and $Re_K$. Their data show that for $Re_K= {O}(1-10)$, the turbulence shear penetration depth, a measure of true roughness felt by the flow, normalized by the permeability scale ($\sqrt {K}$), is a nonlinear function of $Re_K$, as opposed to a commonly assumed linear relationship for $Re_K < {{O}(100)}$ (Ghisalberti 2009; Manes, Ridolfi & Katul 2012). Kim et al. (2020) also investigated – through experimental observations at $Re_K=50$ – the dynamic interplay between surface and subsurface flow in the presence of smooth and rough permeable walls, composed of a uniform cubic arrangement of packed spheres. They confirmed the existence of amplitude modulation, a phenomenon typically identified in impermeable boundaries, whereby the outer large scales modulate the intensity of the near-wall small-scale turbulence. They postulated that amplitude modulation of subsurface flow is driven by large-scale pressure fluctuations at the SWI and is generated by the passage of large-scale motions in the log-law region of the surface flow. However, detailed data on the pressure field at the SWI are needed to confirm these findings, a task that is difficult for experimental measurements, thus requiring pore-resolved direct simulations.

There have been a few pore-resolved direct numerical simulations (DNS) or large-eddy simulations of turbulent flow over permeable flat beds (Breugem & Boersma 2005; Breugem, Boersma & Uittenbogaard 2006; Kuwata & Suga 2016; Fang et al. 2018; Leonardi et al. 2018). However, all of these studies used structured, arranged packings of either non-touching particles or compactly packed particles. Although these numerical studies provide considerable insights into the fundamental aspects of turbulent flow over permeable beds, natural systems involve randomly packed beds with varying particle shapes and sizes. Furthermore, with structured packings such as simple cubic or body-centred cubic, there are open flow pathways that can lead to significant flow penetration and transport, which are not generally present in randomly packed natural systems (Finn 2013; Finn & Apte 2013; Fang et al. 2018). At the time of writing, there has been only one pore-resolved DNS study (Shen, Yuan & Phanikumar 2020), involving a randomly packed arrangement of monodispersed spheres with multiple layers at $Re_K=2.62$ and bed porosity 0.41. This provided significant insights into the flow physics of turbulence over randomly packed beds. To investigate the effect of bed roughness in regular versus randomly packed spheres, Shen et al. (2020) changed only the top layer of the bed to regularly arranged spheres. More intense mixing was observed near the random interface due to increased Reynolds and form-induced stresses, which resulted in a deeper penetration of turbulence (44 % higher) than the uniform, regularly arranged interface. Although this is one of the first studies of flat beds closely resembling natural systems, the investigation was carried out only at low $Re_K$ that falls in the marginally turbulent regime. In addition, since only the top layer of the sediment bed was changed from random to arranged, the open flow pathways that are characteristics of arranged packing were potentially absent in their study.

To the authors’ best knowledge, pore-resolved DNS of randomly packed flat beds over $Re_K\sim {{O}}(10)$ have not been conducted. This range is of critical importance to streamflows as measurements in local creeks show $Re_K$ values of this order (Jackson et al. 2013b, 2015) (see figure 1). The sweep–ejection events over permeable beds can cause significant spatio-temporal variability in the bed shear stresses and pressure forces on the sediment grains as well as at the sediment water interface. Direct measurements of these quantities in experiments pose significant challenges and hence have not been conducted. The present pore-resolved DNS studies aim to provide detailed data on local distribution of the bed stresses as well as drag and lift forces, which can be of importance for incipient motion models. Similarly, models for spatio-temporal pressure variation on the sediment bed can be used as boundary conditions in reach-scale modelling of hyporheic exchange (Chen, Cardenas & Chen 2018). Thus conducting a detailed analysis of the data for a range of $Re_K$ representative of stream and creek flows is of direct relevance in modelling transport across the SWI.

In the present study, pore-resolved DNS of flow over a bed of randomly packed, monodispersed, spherical particles for $Re_K= 2.56$, $5.17$ and $8.94$, representative of transitional to fully turbulent flows, are performed. The main goals of this study are to: (i) characterize the nature of the turbulent flow, Reynolds and form-induced stresses, turbulence penetration depths and sweep–ejection events as functions of $Re_K$; (ii) quantify the spatio-temporal variability of the bed stress and resultant forces on the sediment grains and pressure fluctuations at the SWI using higher-order statistics, and propose a model fit for the p.d.f.s that can be used in reduced-order models; and (iii) quantify the contribution of the top sediment layer on the turbulent flow characteristics over permeable beds.

The rest of the paper is arranged as follows. The methodology, flow domain and simulation parameters are described in § 2. To focus on main results and insights from this work, details on grid refinement and validation studies are presented in Appendices A–D. Details of turbulence structure, mean, turbulent and dispersive stresses, turbulence penetration depths, quadrant analysis, and the role of the top sediment layer, followed by detailed statistics of the bed stress, drag and lift forces on the sediment grains, and pressure variations at the SWI, are presented in § 3. The importance of the results to hyporheic exchange and transport across the SWI is summarized in § 4.

2. Simulation set-up and mathematical formulation

In this section, the flow domain and parameters for cases studied, numerical approach, grid resolution and averaging procedure for analysis are described.

2.1. Simulation domain and parameters

Various non-dimensional parameters relevant to the turbulent flow over a permeable bed made of monodispersed spherical particles are the permeability Reynolds number ($Re_K= u_{\tau } \sqrt {K}/\nu$), the friction or turbulent Reynolds number ($Re_{\tau }= u_{\tau } \delta /\nu$), the roughness Reynolds number ($D^{+}= D_p u_{\tau }/\nu$), the bulk Reynolds number ($Re_b = \delta U_b/\nu$), the ratio of sediment depth to the free-stream height ($H_s/\delta$), the ratio of the sediment grain diameter to the free-stream height ($D_p/\delta$), bed porosity ($\phi$), type of particle packing (random versus arranged), and the domain lengths in the streamwise and spanwise directions normalized by the free-stream height ($L_x/\delta$, $L_z/\delta$). Here, $u_{\tau }$ is the friction velocity, $U_b$ is the bulk velocity, $K$ is the bed permeability, and $\nu$ is the kinematic viscosity. For monodispersed, spherical particles, the size of the roughness element, $k_s$, scales with the permeability ($k_s/\sqrt {K}\approx 9$) (Wilson, Huettel & Klein 2008; Voermans et al. 2017). It should be noted that for a given monodispersed, randomly packed bed of certain porosity and flow rate, only one of the non-dimensional Reynolds numbers is an independent parameter.

Figure 2(a) shows a schematic of the sediment bed and the computational domain used in the present study. A doubly periodic domain in the streamwise ($x$) and spanwise ($z$) directions with four layers of randomly packed, monodispersed sediment grains of porosity $\phi = 0.41$ is shown in figure 2(c). The random packing of monodispersed, spherical particles is generated using the code developed by Dye et al. (2013). The bed porosity profile for low, medium and high permeability Reynolds number cases is shown in figure 2(b). In order to quantify the role and influence of the roughness due to only the top layer of the sediment bed on turbulence structure, and bed pressure and shear stresses, a rough impermeable wall configuration is generated as shown in figure 2(d). The roughness elements match exactly with the top layer of the porous sediment bed and are placed on top of a no-slip solid wall.

Figure 2. (a) Schematic of the permeable bed. (b) Porosity profile. (c) Permeable bed with randomly packed sediment particles (inset shows close-up view in the $xy$-plane). (d) Impermeable-wall medium Reynolds number case with full layer of particles. Periodicity is used in the streamwise ($x$) and spanwise ($z$) directions. The top and bottom surfaces are defined as slip boundaries for the permeable bed cases, whereas for the impermeable-wall full layer (IWM-F) case, the bottom surface is a no-slip wall.

Table 1 shows detailed simulation parameters for the cases used to investigate the structure and dynamics of turbulence over a porous sediment bed. Four permeable bed cases (VV, PBL, PBM and PBH) are studied by varying the permeability Reynolds number. Case VV is used to verify and validate the DNS simulations of turbulent flow over a sediment bed with experimental data from Voermans et al. (2017) as well as simulation data of Shen et al. (2020), and its details are given in Appendix C. The free-stream height for the VV case is based on the DNS study of Shen et al. (2020).

Table 1. Parameters used in the present pore-resolved DNS: $D_p$ is the sphere diameter, $\delta$ is the free-stream height, $H_s$ is the sediment depth, $\phi$ is the porosity, $L_x$ and $L_z$ are the streamwise and spanwise domain lengths, and $Re_K$, $Re_{\tau }$, $Re_b$ and $D^+$ are the permeability, friction, bulk and roughness Reynolds numbers, respectively. Superscript ${+}$ denotes wall units.

Cases PBL, PBM and PBH correspond to low (2.56), medium (5.17) and high (8.94) permeability Reynolds numbers, and are used to investigate the influence of $Re_K$ over the turbulent flow regime shown in figure 1. Finally, case IWM-F is an impermeable-wall case with only the top layer of the PBM sediment bed used as roughness elements over a no-slip surface. The free-stream height $\delta$ for the PBL, PBM, PBH and IWM-F cases is set to be 3.5$D_{p}$ and is similar to the experimental domains of Voermans et al. (2017) and Manes et al. (2009). The length of the streamwise and spanwise domains is based on the DNS of smooth channel turbulent flow (Moser, Kim & Mansour 1999). Moreover, since roughness is expected to break the long, elongated streamwise flow structures observed commonly in smooth channel flow, the domain size used is sufficient to impose the periodicity condition, which was also confirmed by evaluating integral length scales (see Appendix A).

2.2. Numerical method

The numerical approach is based on a fictitious domain method to handle arbitrary-shaped immersed objects without requiring the need for body-fitted grids (Apte, Martin & Patankar 2009). Cartesian grids are used in the entire simulation domain, including both fluid and solid phases. An additional body force is imposed on the solid part to enforce the rigidity constraint and satisfy the no-slip boundary condition. The absence of a highly skewed unstructured mesh at the bead surface has been shown to accelerate the convergence and lower the uncertainty (Finn & Apte 2013). The following governing equations are solved over the entire domain, including the region within the solid bed, and a rigidity constraint force $\boldsymbol f$ that is non-zero only in the solid region is applied to enforce the no-slip condition on the immersed object boundaries:

(2.1)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol u} = 0, \end{gather}$$

(2.2)$$\begin{gather}\rho \left[\frac{\partial {\boldsymbol u}}{\partial t} + \left({\boldsymbol u}\boldsymbol{\cdot} \boldsymbol{\nabla}\right) {\boldsymbol u} \right] ={-}\boldsymbol{\nabla} p + \mu\, \nabla^2{\boldsymbol u} + {\boldsymbol f} , \end{gather}$$

where $\boldsymbol u$ is the velocity vector (with components given by ${\boldsymbol u}=(u,v,w)$), $\rho$ is the fluid density, $\mu$ is the fluid dynamic viscosity, and $p$ is the pressure. A fully parallel, structured, collocated grid, finite volume solver has been developed and verified and validated thoroughly for a range of test cases, including flow over a cylinder and sphere for different Reynolds numbers, flow over touching spheres at different orientations, and flow developed by an oscillating cylinder. The details of the algorithm as well as very detailed verification and validation studies have been published elsewhere (Apte et al. 2009; Finn 2013; Finn & Apte 2013). The solver was used to perform direct one-to-one comparison with a body-fitted solver with known second-order accuracy for steady inertial, unsteady inertial and turbulent flow through porous media (Finn & Apte 2013) to show very good predictive capability. It has also been used recently for DNS of oscillatory, turbulent flow over a sediment layer (Ghodke & Apte 2016, 2018a), and pore-resolved simulations of turbulent flow within a porous unit cell with face-centred cubic packing (He et al. 2018, 2019).

2.3. Averaging

Since the flow properties are highly spatially heterogeneous near the rough sediment bed boundary, time averaging followed by spatial averaging is applied. Flow statistics such as Reynolds stresses, form-induced disturbances, shear stress and pressure fluctuations, among others, are computed using the time–space averaging. This consecutive time–space averaging involves Reynolds decomposition ($\psi = \bar {\psi } + \psi ^{\prime }$) accompanied by spatial decomposition of the time-averaged variable $\bar {\psi } = \langle {\bar {\psi }}\rangle + \tilde {\bar {\psi }}$ (Nikora et al. 2007, 2013). Here, $\psi$ represents an instantaneous flow variable, $\bar {\psi }$ is its temporal average, and $\psi ^{\prime } = \psi - \bar {\psi }$ is the instantaneous turbulent fluctuation. The angular brackets, $\langle \ \rangle$, denote the spatial averaging operator. The quantity $\tilde {\bar {\psi }}$, known as the form-induced or dispersive disturbance in space, is defined as $\tilde {\bar {\psi }} = \bar {\psi } - \langle \bar {\psi }\rangle$. This represents the deviation of the time-averaged variable $\bar {\psi }$ from its spatially averaged value $\langle \bar {\psi }\rangle$. Nikora et al. (2013) proposed to denote the form-induced disturbance quantity as $\tilde {\bar {\psi }}$, whereby a horizontal overbar is added to emphasize that time averaging has been done prior to spatial averaging. This modified notation is adopted in this work. The original notation of this quantity proposed in Nikora et al. (2007) was without the overbar, $\tilde {\psi }$, and has been used in the literature published in this field (Voermans et al. 2017; Fang et al. 2018; Shen et al. 2020).

Quantities are averaged over the fluid domain, giving the intrinsic spatial average of the time-averaged variable, $\langle {\overline \psi }\rangle = 1/V_f \int _{V_f} \overline \psi \,{\rm d}V$, where $V_f$ is the volume occupied by the fluid. In other words, while calculating the volume average of variables in the particle bed regions, only the portion of the volume occupied by the fluid is taken into account. The representative averaging volume lengths used for spatial averaging in the streamwise and spanwise directions are the same as the grid resolution in those directions. Since the grids are uniform in $x$ and $z$, this implies that the averaging volumes have the same lengths in the streamwise and spanwise direction. However, in the bed-normal direction, a variable volume averaging is used. In the boundary layer region, especially near the crest of the bed where steep gradients in flow quantities are present, thin volumes are used for averaging, whereas deeper inside the bed, thicker averaging volumes are used, as described in detail in Appendix B.

2.4. Grid resolution and flow set-up

Table 2 gives the details of grid resolution used for all cases. The grid resolutions required for these configurations are based on two main considerations: (i) minimum bed-normal grid resolution near the bed to capture the bed shear stress; and (ii) minimum resolution required to capture all details of the flow over spherical particles. Since the intensity of turbulence and mean flow penetration into the bed is expected to reduce further deep into the bed, finer grid resolutions are used in the bed-normal region in the top layer compared to other layers.

Table 2. Grid parameters used in the present DNS, where superscript ${+}$ denotes wall units. Details on the grid transition regions are given in Appendix B.

For DNS of channel flows, the bed-normal grid resolution in wall units should be $\Delta y^{+}<1$, in order to capture accurately the bed shear stress in the turbulent flow. The grid resolutions in the streamwise and spanwise directions are typically 3–4 times coarser, following the smooth channel flow simulations by Moser et al. (1999). Note that the roughness features and permeability are known to break the elongated flow structures along the streamwise direction in smooth walls, reducing the anisotropy in the near-bed region (Ghodke & Apte 2016). To capture the inertial flow features within the pore and around spherical particles, grid refinement studies were conducted on flow over a single sphere at different Reynolds numbers representative of the cases studied here (details in Appendix A). Accordingly, roughly 90 (PBL), 180 (PBM) and 548 (PBH) grid points are used in the bed-normal direction in a region covering the top layer and extending slightly into the free-stream. In the $x$- and $z$-directions, a uniform grid with a minimum of 26 (PBL), 38 (PBM) and 40 (PBH) grid points per sediment grain is used to resolve the bed geometry. The effect of uniform but non-cubic grids within the sediment bed was evaluated thoroughly by comparing drag coefficients to those obtained from cubic grids over a single layer of particles to show no discernible differences over the range of Reynolds numbers studied here (see Appendix A).

Below the top layer, the bed-normal resolution is reduced slowly and a grid that is nearly uniform in all directions is used deep inside the bed, as the frictional Reynolds number in the bottom layers of sediment decreases significantly. From the crest of the top sediment layer, the grid is stretched, coarsening gradually towards the top of the channel using a standard hyperbolic tangent function (Moser et al. 1999). Based on these grid resolutions, the total grid count for the PBL case is $\sim$232 million cells, for the PBM case is $\sim$200 million cells, and for the PBH case is $\sim$428 million cells, as given in table 2.

The flow in the simulations is driven by a constant mass flow rate. A target mass flow rate is adjusted until the friction velocity $u_{\tau }$ that results in the required $Re_K$ is obtained. Then $Re_{\tau }$ is calculated based on the free-stream height $\delta$. Pokrajac et al. (2006) noted that there is lack of a general definition of $u_{\tau }$ applicable to the boundary layers with variable shear stress where the roughness height is comparable with the boundary layer thickness. Accordingly, $u_{\tau }$ can be specified based on (i) bed shear stress, (ii) total fluid shear stress based on the roughness crest, (iii) fluid shear stress extrapolated to the zero-displacement plane, or (iv) shear stress obtained by fitting data to log law. Pokrajac et al. (2006) proposed to use $u_{\tau }$ based on the the fluid shear stress at the roughness crest, to obtain the least ambiguous definition. Since the present work builds upon the experimental data of Voermans et al. (2017), the friction velocity $u_{\tau }$ is calculated from the maximum value of the time–space-averaged total fluid stress, which happens to be very close to the sediment crest. The friction velocity is then based on the sum of the viscous, turbulent and form-induced shear stresses (Nikora et al. 2004; Voermans et al. 2017) as

(2.3)\begin{equation} \tau(y) = \rho\nu\,\partial(\phi\langle\bar{u}\rangle)/\partial{y} - \rho\phi\,\langle\overline{u^{\prime}v^{\prime}}\rangle - \rho\phi\,\langle\tilde{\bar{u}} \tilde{\bar{v}}\rangle. \end{equation}

Following smooth-wall DNS studies by Moser et al. (1999), between 20–25 flow-through times (computed as the length over the average bulk velocity $L_x/U_{b}$) are needed for the turbulent flow to reach stationary state. Once a stationary flow field is obtained, computations were performed for an additional time period $T = 25 \delta /u_{\tau }$ to collect single-point and two-point statistics. For the PBL case, where the domain in the streamwise and spanwise directions is twice that used in the PBM and PBH cases, the statistics were collected for over a time period $T = 13 \delta /u_{\tau }$. The flow statistics were monitored to obtain statistically stationary values over the above averaging periods.

3. Results

The main results for the different cases studied in this work are discussed. The Reynolds and form-induced stresses are first compared for different Reynolds numbers (§ 3.1). The structure of turbulence is visualized using vorticity contours, followed by a detailed quadrant analysis describing the sweep and ejection events (§ 3.2). Variation of the turbulence penetration depth, interfacial mixing length and similarity relations is investigated as a function of $Re_K$ (§ 3.3). Next, the role of the top layer of the sediment is quantified by comparing the permeable bed statistics with an impermeable rough wall with roughness equivalent to the permeable bed (§ 3.4). Finally, p.d.f.s of the viscous and pressure components of the normalized bed shear stress are presented in § 3.5, followed by the pressure fluctuations at the SWI in § 3.6.

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 (1991), Breugem et al. (2006) and Fang et al. (2018). 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 1990). 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. (2017), 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. (2017) indicated by blue diamonds, DNS indicated by black circles. (b) Turbulent shear stress penetration $\delta _{p}$, normalized by $\sqrt {K}$; Voermans et al. (2017) 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. (2017) 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. (2018) 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 (2009) on obstructed shear flows, Manes et al. (2012) 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. (2012) 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 1972). 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 (2009). 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. (2017). 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. (2017, 2018) has important implications for flows over monodispersed sediment beds in general. Voermans et al. (2017) 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. (2009) 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 2018b). 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. (1999) 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 (2011). 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. 2013a, 2015) 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. 2022). 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.

4. Conclusions

Pore-resolved direct numerical simulations (DNS) of turbulent flow over a randomly packed, monodispersed bed of spherical particles are performed for three permeability Reynolds numbers, $Re_K= 2.56$, 5.17 and 8.94 ($Re_{\tau } = 270$, 545 and $943$), in the hydrodynamically fully rough regime representative of natural stream systems. A thoroughly validated fictitious domain based numerical approach (Apte et al. 2009; Finn & Apte 2013; Ghodke & Apte 2016, 2018a; He et al. 2019) is used to conduct these simulations. The numerical computations are first validated against the experimental data by Voermans et al. (2017) at $Re_K \sim 2.56$, and are then used to investigate different Reynolds numbers. The time–space averaging methodology is used to compute the mean velocity, Reynolds stresses and form-induced stresses. Differences in the near-bed turbulence structure, statistics of the local distribution of the net bed stress on the sediment grains and the resultant drag and lift forces, pressure distributions at the sediment–water interface (SWI), and the contribution of the top layer of sediment grains to turbulence statistics were quantified in detail. The key findings of this work are summarized below.

1. (i) The peak and significant values of the Reynolds stresses occur in the top layer of the bed for all three $Re_K$ cases, decreasing quickly below one grain diameter from the sediment crest. While peak values in streamwise stress decrease, those in bed-normal and spanwise stresses increase with increasing Reynolds number. Streamwise, bed-normal and shear Reynolds stresses exhibit similarity in the free-stream region, substantiating the wall similarity hypothesis.

2. (ii) Form-induced stresses are typically lower in magnitude than their respective Reynolds stress counterparts, with the locations of their peak values occurring further below the crest. For low $Re_K$, the spanwise form-induced stress is typically larger than the bed-normal values, a result similar to a rough, impermeable wall. However, at higher $Re_K$, the bed-normal stresses are comparable to the spanwise stresses due to increased flow penetration.

3. (iii) Mean flow penetration (Brinkman layer thickness) and shear penetration show nonlinear, increasing correlation with $Re_K$, their ratio approaching a constant deterministic value. The length scale of dominant turbulent eddies at the SWI is better represented by the mixing length obtained from the Reynolds stress and mean flow gradient. The normalized interfacial mixing length at the SWI increases with $Re_K$ and shows behaviour similar to that of the Brinkman layer thickness and shear penetration depth, suggesting that the mixing length is relevant as the characteristic length scale for transport of momentum and mass across the SWI. Quadrant analysis for turbulent fluctuations showcases domination of ejection and sweep events at the SWI. Within one particle diameter inside the bed, the turbulent structures lose both their directional bias and strength, becoming more isotropic in nature.

4. (iv) To quantify the role of the top layer of the sediment bed on the flow, an impermeable rough wall with the same roughness elements as the top layer was investigated. The mean and Reynolds stress profiles show very little difference between the permeable bed and impermeable rough wall. Form-induced stresses are, however, influenced by the impermeability, redistributing stresses tangentially along the wall.

5. (v) The form-induced disturbances and the turbulent fluctuations in pressure are strongly dependent on Reynolds number, with a ten-fold increase in the peak value between the lowest and highest $Re_K$ cases studied. This increase is attributed to the nature of the flow over the exposed particles in the top layer. A majority of the high-magnitude pressure fluctuations are restricted to the top layer of the bed for all Reynolds numbers. The standardized p.d.f.s of the sum of the pressure fluctuations and form-induced pressure disturbances at the SWI are statistically similar and symmetric, and collapse for different $Re_K$.

6. (vi) The p.d.f.s of local distribution of the net bed stress computed directly on the sediment grains, and normalized by the total bed stress in the streamwise direction, collapse for all Reynolds numbers. The p.d.f.s of fluctuations in the bed stress normalized by their root-mean-square values are symmetric and exhibit a peaky, non-Gaussian distribution with heavy tails. A logarithmic correlation between the root mean-square stress fluctuations and the friction Reynolds number (as well as $Re_K$) is observed, which together with the non-Gaussian distribution for fluctuations in stress can be used to develop mechanistic force balance models for incipient motion of sediment grains.

7. (vii) The mean and fluctuations in drag and lift forces on the particle are computed by integrating the local bed stress on the particle surface. The majority of the contribution to the drag and lift forces comes from the pressure distribution for all $Re_K$. In addition, the top layer of the bed results in average forces that are 3–4 times the average value in the full bed, indicating that a significant contribution to the lift and drag force comes from the top layer of the bed. Fluctuations in drag and lift forces have minimal skewness and high kurtosis, indicative of a symmetric, non-Gaussian distribution with heavy tails. A $t$-location shape function model based on the skewness and excess kurtosis data fits well for all $Re_K$. Since the local distributions of the net bed stress and the drag and lift forces on particles are influenced mainly by the top layer of the sediment grain, including the roughness effects through a single layer of randomly arranged sediments potentially can improve reach-scale predictions based on reduced-order models.