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

Shashank K. Karra1  Sourabh V. Apte1<EMAIL_ADDRESS>
Xiaoliang He2  Timothy D. Scheibe2
1School of Mechanical, Industrial and Manufacturing Engineering, Oregon State
University, Corvallis, OR, 97331 2Pacific Northwest National Laboratory,
Richland, WA, 99354

###### Abstract

Pore-resolved direct numerical simulations (DNS) 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}$, of 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 (PDFs) of normalized local bed stress are
found to collapse for all Reynolds numbers and their root mean-squared
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 PDFs 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 can potentially be used in
providing better boundary conditions in modeling of incipient motion and
reach-scale transport in the hyporheic zone.

##  1 Introduction

The interchange of mass and momentum between streamflow and ground water
occurs across the sediment-water interface (SWI) and into the porous bed
underneath, termed as 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
source of fresh drinking water, transport, biogeochemical processing of
nutrients, and sustaining diverse aquatic ecosystems (Bencala et al., 1983;
D’angelo et al., 1993; Valett et al., 1996; Harvey et al., 1996; Anderson et
al., 2008; Briggs et al., 2009; Grant et al., 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 equation 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}$)
to 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 et al. (2017, 2018); 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 et al., 2013a, 2015), it is found that
the gravel grain sizes varied over the range 5–70mm, with friction velocity,
$u_{\tau}=0.004$–$0.088~{}\rm m/s$, and $Re_{K}=2$–70 (figure 1). Free surface
flow and waves typically do not affect the hyporheic exchange in natural
stream 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,
2011; Suga et al., 2010; 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 pattern of single and multiple layers of
spheres at $Re_{K}$ of 31.2 and 44.6. Permeability was shown to influence flow
resistance dramatically, and the conventional assumption of hydraulically-
rough regime, wherein the friction factor is dependent upon the relative
submergence or the ratio of the roughness size to flow thickness, does not
apply to permeable beds. Friction factor was shown to progressively increase
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 of $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 shows that for $Re_{K}=\mathcal{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 non-linear
function of $Re_{K}$, as opposed to a commonly assumed linear relationship for
$Re_{K}<{\mathcal{O}(100)}$ (Ghisalberti, 2009; Manes et al., 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 are 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 is needed to confirm
these findings, a task difficult for experimental measurements, and thus
requiring a need for pore-resolved direct simulations.

There have been a few pore-resolved DNS or large-eddy simulations (LES) of
turbulent flow over permeable flat beds (Breugem & Boersma, 2005; Breugem et
al., 2006; Kuwata & Suga, 2016; Leonardi et al., 2018; Fang 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-centered 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 only been
one pore-resolved DNS study (Shen et al., 2020) involving randomly packed
arrangement of monodispersed spheres with multiple layers at $Re_{K}=2.62$ and
bed porosity of 0.41. They 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, they 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 only carried out 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{\mathcal{O}}(10)$ have not been conducted. This range is of
critical importance to stream flows as measurements in local creeks show
$Re_{K}$ values on 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 modeling of hyporheic
exchange (Chen et al., 2018). Thus, conducting a detailed analysis of the data
for a range of $Re_{K}$ representative of stream and creek flows are of direct
relevance in modeling transport across the sediment-water interface.

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) first characterize the nature of the
turbulent flow, Reynolds and form-induced stresses, turbulence penetration
depths, and sweep-ejection events as a function 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 probability distribution functions
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 section 2. To focus on main results
and insights from this work, details on grid refinement and validation studies
are presented in the Appendix. 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 section 3 . Importance of the
results to hyporheic exchange and transport across the SWI is summarized in
section 4.

##  2 Simulation Setup 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 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 of size of the
roughness element, $k_{s}$, scales with the permeability
($k_{s}/\sqrt{K}\approx 9$) (Wilson et al., 2008; Voermans et al., 2017)). It
should be noted that for a given mono-dispersed, randomly packed bed of
certain porosity and flow rate, only one of the non-dimensional Reynolds
numbers is an independent parameter.

Figure 2a shows the schematic of the sediment bed and the computational domain
used in the present study. A doubly periodic domain in streamwise ($x-$) and
spanwise ($z-$) directions with four layers of randomly packed, monodispersed
sediment grains of porosity $\phi=0.41$ is shown in figure 2c. 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 2b. 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 2d. 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.

(a) Domain schematic (b) Bed porosity (c) Permeable bed (d) IWM-F

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 xy-plane), and (e) impermeable-wall medium Reynolds number case with full layer of particles. Periodicity is used in the streamwise ($x$) and spanwise ($z$) directions. Top and bottom surface 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. Case | Domain | $Re_{K}$ | $Re_{\tau}$ | $Re_{b}$ | $D^{+}$ | $\phi$ | $H_{s}/{\delta}$ | $D_{p}/{\delta}$ | ${(L_{x},L_{z})/\delta}$
---|---|---|---|---|---|---|---|---|---
VV | permeable | 2.56 | 180 | 1,886 | 77 | 0.41 | 1.71 | 0.43 | ($4\pi$,$2\pi$)
PBL | permeable | 2.56 | 270 | 2,823 | 77 | 0.41 | 1.14 | 0.29 | ($4\pi$,$2\pi$)
PBM | permeable | 5.17 | 545 | 5,681 | 156 | 0.41 | 1.14 | 0.29 | ($2\pi$,$\pi$)
PBH | permeable | 8.94 | 943 | 9,965 | 270 | 0.41 | 1.14 | 0.29 | ($2\pi$,$\pi$)
IWM-F | impermeable | - | 545 | 5,683 | 156 | - | 0.29 | 0.29 | ($2\pi$,$\pi$)
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, $Re_{K}$, $Re_{\tau}$, $Re_{b}$ and $D^{+}$ are the
permeability, friction, bulk and roughness Reynolds numbers, respectively.
$(\hphantom{0})^{+}$ denotes wall units.

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).

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 be
3.5$D_{p}$ and is similar to the experimental domains of Voermans et al.
(2017); Manes et al. (2009). The length of the streamwise and spanwise domains
is based on the DNS of smooth channel turbulent flow (Moser et al., 1999).
Moreover, since roughness is expected to break the long, elongated streamwise
flow structures commonly observed 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 et al., 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 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, $\bf f$, that is non-zero only in the solid region
is applied to enforce the no-slip condition on the immersed object boundaries.
The governing equations are

$\displaystyle\nabla\cdot{\bf u}$ $\displaystyle=0,$ (1)
$\displaystyle\rho\bigg{[}\frac{\partial{\bf u}}{\partial t}+\left({\bf
u}\cdot\nabla\right){\bf u}\bigg{]}$ $\displaystyle=-\nabla
p+\mu\nabla^{2}{\bf u}+{\mathbf{f}}\>,$ (2)

where $\bf u$ is the velocity vector (with components given by ${\bf
u}=(u,v,w)$, $\rho$ the fluid density, $\mu$ the fluid dynamic viscosity, and
$p$ the pressure. A fully parallel, structured, collocated grid, finite volume
solver has been developed and thoroughly verified and validated for a range of
test cases including flow over a cylinder and sphere for different Reynolds
numbers, flow over touching spheres at different orientations, flow developed
by an oscillating cylinder, among others. 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 recently used for direct numerical
simulations 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-centered cubic packing (He et al., 2018, 2019).

###  2.3 Averaging

Since the flow properties are highly spatially heterogeneous near 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=\overline{\psi}+\psi^{\prime}$) accompanied by spatial
decomposition of the time averaged variable
$\overline{\psi}=\langle{\overline{\psi}}\rangle+\widetilde{\overline{\psi}}$
(Nikora et al., 2007, 2013). Here $\psi$ represents an instantaneous flow
variable, $\overline{\psi}$ is its temporal average,
$\psi^{\prime}=\psi-\overline{\psi}$ is the instantaneous turbulent
fluctuation. The angular brackets, $\langle\rangle$, denote spatial averaging
operator. The quantity $\widetilde{\overline{\psi}}$ known as the form-induced
or dispersive disturbance in space is defined as
$\widetilde{\overline{\psi}}=\overline{\psi}-\langle\overline{\psi}\rangle$.
This represents the deviation of the time-averaged variable,
$\overline{\psi}$, from its spatially averaged value,
$\langle\overline{\psi}\rangle$. Nikora et al. (2013) proposed to denote the
form-induced disturbance quantity as $\widetilde{\overline{\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, $\widetilde{\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}dV$, 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 setup

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.

For DNS of channel flows, the bed-normal grid resolution in wall units should
be $\Delta y^{+}<1$, in order to accurately capture 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 are used to resolve the bed
geometry. Effect of uniform, but non-cubic grids within the sediment bed was
thoroughly evaluated by comparing drag coefficients to those obtained from
cubic grids over a single and a 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 slowly reduced and nearly
uniform grid 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 it
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.

Case | $N_{x}\times N_{y}\times N_{z}$ | Bed-Normal Grid Distribution | ($\Delta x^{+},\Delta y^{+},\Delta z^{+})$
---|---|---|---
|  | Channel region | Top layer | Bottom layers |
VV | $768\times 288\times 384$ | 96 | 86 | 106 | (2.94, 0.95, 2.94)
PBL | $1152\times 350\times 576$ | 150 | 90 | 110 | (2.94, 0.95, 2.94)
PBM | $846\times 530\times 448$ | 184 | 180 | 166 | (4.01, 0.95, 3.8)
PBH | $882\times 1082\times 448$ | 342 | 548 | 192 | (6.74, 0.55, 6.63)
IWM-F | $846\times 364\times 448$ | 184 | 180 | - | (4.01, 0.95, 3.8)
Table 2: Grid parameters used in the present DNS. $(\hphantom{0})^{+}$ denotes
wall units. Details on the grid transition regions are given in appendix B.

The flow in the simulations is driven by constant mass flow rate. A target
mass flow rate is adjusted until the friction velocity, $u_{\tau}$, which
results in the required $Re_{K}$ is obtained. $Re_{\tau}$ is then 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 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 the form-induced shear stresses (Nikora et al., 2004; Voermans
et al., 2017),

$\displaystyle\tau(y)$
$\displaystyle=\rho\nu\partial(\phi\langle\overline{u}\rangle)/\partial{y}-\rho\phi\langle\overline{u^{\prime}v^{\prime}}\rangle-\rho\phi\langle\widetilde{{}\mkern
3.0mu\overline{\mkern-3.0muu}}\widetilde{{}\mkern
3.0mu\overline{\mkern-3.0muv}}\rangle.$ (3)

Following smooth wall DNS studies by Moser et al. (1999), between 20-25 flow-
through times (computed as the length over 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 of $T=25\delta/u_{\tau}$ to collect single-point and
two-points statistics. For the PBL case, where the domain in the streamwise
and spanwise directions is twice as that used in the PBM and PBH cases, the
stats 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 (section 3.1). The structure of turbulence is visualized
using vorticity contours followed by a detailed quadrant analysis describing
the sweep and ejection events (section 3.2). Variation of the turbulence
penetration depth, interfacial mixing length, and similarity relations are
investigated as a function of $Re_{K}$ (section 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 (section 3.4). Finally, probability distribution function (PDF)
of the viscous and pressure components of the normalized bed shear stress are
presented in section 3.5 followed by the pressure fluctuations at the
sediment-water interface in section 3.6.

###  3.1 Reynolds and form-induced stresses

Figure 3: Comparison of Reynolds (lines) and form-induced (lines and symbols)
stress profiles for PBL ( ,  ), PBM ( ,  ), and PBH ( ,  ) cases: (a-c)
streamwise, bed-normal, and spanwise components of 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, and (h) mean-square pressure fluctuations and form-induced
pressure disturbances. Horizontal lines show the crest of sediment bed for PBL
( ), PBM ( ), and PBH ( ) cases. Pressure is normalized by $\rho
u_{\tau}^{2}$.

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 sediment-water interface (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 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 location of the
peaks are also different, with the form induced stresses peaking below the
sediment-water interface, whereas the Reynolds stresses peak close to the bed
crest.

The profiles of streamwise ($x$–direction), $\langle\overline{u^{\prime
2}}\rangle^{+}$, and bed-normal ($y$–direction) Reynolds stress,
$\langle\overline{v^{\prime 2}}\rangle^{+}$, (figure 3a,b) exhibit similarity
for $(y+d)/\delta>=0.4$, substantiating the wall similarity hypothesis
reported by Raupach et al. (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 section 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 location of peak in $\langle\overline{u^{\prime
2}}\rangle^{+}$ and $\langle\overline{w^{\prime 2}}\rangle^{+}$ is close to
the crest for all $Re_{K}$, whereas that in $\langle\overline{v^{\prime
2}}\rangle^{+}$ shifts downward 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 3d also peaks near the crest, increases with $Re_{K}$,
and shows similarity in the outer region.

Case | $\langle\overline{u^{\prime 2}}\rangle^{+}$ | $\langle\overline{v^{\prime 2}}\rangle^{+}$ | $\langle\overline{w^{\prime 2}}\rangle^{+}$ | $\langle\overline{u^{\prime}v^{\prime}}\rangle^{+}$ | $\langle\overline{p^{\prime 2}}\rangle^{+}$
---|---|---|---|---|---
PBL | 3.46 (0.15) | 1.03 (0.25) | 1.79 (0.18) | 0.80 (0.16) | 9.41 (0.094)
PBM | 3.19 (0.15) | 1.18 (0.19) | 1.88 (0.16) | 0.84 (0.15) | 40.25 (0.095)
PBH | 2.97 (0.15) | 1.23 (0.16) | 2.05 (0.16) | 0.88 (0.15) | 120.88 (0.097)
Case | $\langle\widetilde{\overline{u}}^{2}\rangle^{+}$ | $\langle\widetilde{\overline{v}}^{2}\rangle^{+}$ | $\langle\widetilde{\overline{w}}^{2}\rangle^{+}$ | $\langle\widetilde{{}\mkern 3.0mu\overline{\mkern-3.0muu}}\widetilde{{}\mkern 3.0mu\overline{\mkern-3.0muv}}\rangle^{+}$ | $\langle\widetilde{\overline{p}}^{2}\rangle^{+}$
PBL | 2.71 (0.08) | 0.25 (-0.04) | 0.38 (0.014) | 0.21 (0.05) | 3.7 (0.07)
PBM | 2.45 (0.07) | 0.47 (-0.05) | 0.42 (0.01) | 0.29 (0.04) | 16.5 (0.063)
PBH | 2.33 (0.07) | 0.53 (-0.05) | 0.46 (0.009) | 0.31 (0.04) | 49.9 (0.06)
Table 3: The peak value and location [$(y+d)/\delta$] given in brackets, of
Reynolds stresses and form-induced stresses for the PBL, PBM, and PBH cases.
The peak values of 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 figure 3d–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 sediment crest. Inside the bed,
$\langle\widetilde{\overline{u}}^{2}\rangle^{+}$ decreases with $Re_{K}$,
while both $\langle\widetilde{\overline{v}}^{2}\rangle^{+}$ and
$\langle\widetilde{\overline{w}}^{2}\rangle^{+}$ values increase with $Re_{K}$
mimicking the trend with Reynolds stresses. However, in contrast to
$\langle\overline{v^{\prime 2}}\rangle^{+}$, the peaks in
$\langle\widetilde{\overline{v}}^{2}\rangle^{+}$ occur at a similar location
(table 3) for all Reynolds numbers. This suggests that the penetration depth
of $\langle\widetilde{\overline{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\widetilde{\overline{w}}^{2}\rangle^{+}$ are larger than the
$\langle\widetilde{\overline{v}}^{2}\rangle^{+}$ at lower $Re_{K}$, but are
comparable and slightly smaller than
$\langle\widetilde{\overline{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 3h, 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 formed induced
disturbances in pressure. The location of the peak in form-induced
disturbances as well as turbulent fluctuations occur just below the crest and
is almost the same for all $Re_{K}$ studied. The fluctuations quickly decay
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 in 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 sediment-water interfaces 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
first shown in this section 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$ (same $Re_{\tau}$ as the PBL case) are also shown in figure 4a
for comparison. In the smooth wall case, distinct long elongated streaky
structures, which are a result of low and high speed streaks generated quasi-
streamwise vortices, are visible. The influence of strong mean gradient and an
impenetrable smooth wall results in these long streaky structures (Lee et al.,
1990). In the low Reynolds number permeable bed case (PBL), shown in figure
4b, the starting 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 strong bed-normal velocities prevents the formation of
these long streaky structures and lead 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$, between $0<x/\delta<6.28$ and
$0<z/\delta<3.14$ for (a) smooth-wall case simulated at $Re_{\tau}=270$, (b)
PBL case, (c) PBM case and (d) PBH case.

Figure 5: Quadrant analysis of probability distribution functions of turbulent
velocity fluctuations for the PBL (left panel), PBM (middle panel), and PBH
(right panel) cases at different elevations: $y/\delta$ ($y/D_{p}$) of (a–c)
0.143 (0.5), (d–f) 0 (0), (g–i) -0.143 (-0.5), (j–l) -0.286 (-1). The
velocities have been normalized by $u_{\tau}$.

Quadrant analysis is performed to understand the influence of near bed flow
structures on Reynolds stress. Joint probability distribution functions (PDFs)
for turbulent velocity fluctuations, $u^{\prime}$ and $v^{\prime}$ calculated
at different elevations from the sediment bed are shown in figure 5. The
product of $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 correlation for different $Re_{K}$ at four bed-
normal elevations within one particle diameter below and above the crest.
Figures 5a–c are one diameter above the crest in the free-stream, figures 5d–f
are at the bed crest, figures 5g–i are half-way into the top layer of the bed,
and figures 5j–l are at the bottom of the top layer ($y/D_{p}=-1$),
respectively.

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}$. 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 PDFs are more concentrated in a narrow band,
indicating large fluctuations in $u^{\prime}$ are associated with smaller
excursions in $v^{\prime}$. With increase in Reynolds number, the PDFs 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 5a–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 section 3.2. These
sweep events induce momentum fluxes that can be on 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 to
the characteristic size of the turbulent eddy across the sediment-water
interface (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{\overline{u}}\rangle(y)$) and Darcy velocity ($U_{p}$) has decayed
to 1% of the velocity value at the SWI ($U_{i}$), i.e
$\langle{\overline{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 6a 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 than those in the
experiments.

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 6b
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 non-linear 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 of $1.1$ as shown in figure 6c.

Figure 6: (a) Brinkman layer thickness, $\delta_{b}$, normalised by
$\sqrt{K}$; Voermans et al. (2017) ( ), DNS ( ), (b) turbulent shear stress
penetration, $\delta_{p}$, normalised by $\sqrt{K}$; , Voermans et al. (2017)
( ), DNS ( ), (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) ( ), DNS ( ), and (d) Left axis: Ratio of
interfacial mixing-length to the Kolmogorov length scale; Voermans et al.
(2018) ( ), DNS( ), and Right axis: interfacial mixing-length normalised by
$\sqrt{K}$; DNS ( ).

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 interfacial mixing length,
$\langle
L_{m,i}\rangle=(\langle\overline{u^{\prime}v^{\prime}}\rangle_{i}/(\partial_{y}\langle\overline{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 on 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 dissipation rate of the turbulent
kinetic energy (Tennekes et al., 1972). Figure 6d 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 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-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 closely with the experimental data 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}$ is a necessary but not a
sufficient condition. A random as opposed to structured arrangement of
monodispersed particles in the top layer is also important to achieve
statistical similarity and 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 2d corresponding to the case IWM-F in table 1.

Figure 7: Comparison of the mean velocity, Reynolds (lines) and form-induced
(lines and symbols) stress profiles for PBM ( ,  ) and IWM-F ( , ) cases: (a)
mean velocity, (b-c) streamwise, and bed-normal components of 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, and (h) mean-square pressure fluctuations
and form-induced pressure disturbances. The crest lines for PBM and IWM-F
cases overlap and are shown by the horizontal lines ( ,  ). Pressure is
normalized by $\rho u_{\tau}^{2}$.

Figures 7a–h compare the bed-normal variation of mean velocity, Reynolds and
form-induced stresses, and pressure disturbances for the PBM and impermeable
wall 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 behavior was
observed and reported by Manes et al. (2009) in their experimental work.

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 streamwise component
($\langle\widetilde{\overline{u}}^{2}\rangle^{+}$) is well captured, the peaks
in bed-normal ($\langle\widetilde{\overline{v}}^{2}\rangle^{+}$) and spanwise
($\langle\widetilde{\overline{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 by larger peak magnitude of
$\langle\widetilde{\overline{w}}^{2}\rangle^{+}$ compared to
$\langle\widetilde{\overline{v}}^{2}\rangle^{+}$.

The form-induced pressure disturbances
($\langle\widetilde{\overline{p}}^{2}\rangle^{+}$), shown in figure 7h,
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
$\left({y+d}\right)/{\delta}<0.3$. The turbulent pressure fluctuations,
however, are well captured and are typically much larger than the form-induced
pressure disturbances. 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 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 modeled based on the bed shear stress
exceeding a critical value. Higher-order statistics of bed shear stress are
important in stochastic modeling of incipient motion (Ghodke & Apte, 2018b).
Motion and rearrangement of sediments can alter local bed porosity and
effective permeability, directly impacting hyporheic exchange. The DNS data
are used to compute probability distribution functions (pdfs) 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 ($\bm{\tau}^{t}={\bm{\tau}}^{v}-p{\mathbf{I}}$) on the particle
surface includes contribution from the viscous (skin-friction,
${\bm{\tau}}^{v}$) as well as pressure (form, $p{\mathbf{I}}$) stresses. The
normal and tangential components of the net stress can be directly evaluated
on the particle surface from the velocity and the pressure fields and then
transformed into the Cartesian frame using the surface normal vector to obtain
the streamwise (${\tau}_{x}^{t}$), the bed-normal (${\tau}_{y}^{t}$), and the
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,

${\mathbf{F}}=\int_{\Gamma}{\bm{\tau}}^{t}\cdot{\mathbf{n}}~{}d\Gamma\equiv\int_{\Gamma}\left({\bm{\tau}}^{v}-p{\mathbf{I}}\right)\cdot{\mathbf{n}}~{}d\Gamma,$
(4)

where $d\Gamma$ is the unit surface area of the particle, $\mathbf{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 probability distribution functions 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 8a-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 PDFs 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 section
3.2. The probability of negative $\tau_{x}^{t+}$ fluctuations is then highly
associated with trough regions between the roughness elements and are
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 PDFs 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.

The relationship between the two viscous components of shear stress can be
understood by 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 of range of
$-45$ to $45$ degrees, indicating that events with large values of
$\tau_{x}^{v}$ are associated with small $\tau_{z}^{v}$. However, in flow over
rough permeable beds, probability of yaw angles above $45$ degrees is much
higher as seen in figure 8c. 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.

Figure 8: PDFs of the net bed stress components (a) streamwise
($\tau_{x}^{t}$) and (b) bed-normal ($\tau_{y}^{t}$). Superscript $(.)^{+}$
denotes normalized quantity by the total bed stress in the $x$ direction
($\tau_{w}^{t}$). Also shown is the yaw angle $\psi_{\tau}$ based on the
viscous components of the stresses (c). Legend: PBL ( ), PBM ( ), and PBH ( )
cases.

Figure 9: PDFs of net stress fluctuations normalized by their root-mean-square
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 ( ), PBM ( ), and PBH ( ) cases. Also shown is a non-Gaussian $t$-location
scale fit ( ).

As the normalized PDFs of the local distribution of the net bed stress
collapse for different $Re_{K}$, it is conjectured that a simple deterministic
fit to the PDFs of fluctuations in stress is possible and is investigated. The
PDFs of the local distribution of the net stress fluctuations on the particle
surface normalized by their rms values are shown in figure 9a–c. The higher
order statistics for the net stress are shown in the 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 PDFs 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.

| $\tau^{t^{\prime}}_{x}$ | $\tau^{t^{\prime}}_{y}$ | $\tau^{t^{\prime}}_{z}$
---|---|---|---
Case | $\hat{\sigma}(\cdot)$ | $Sk(\cdot)$ | $Ku(\cdot)$ | $\hat{\sigma}(\cdot)$ | $Sk(\cdot)$ | $Ku(\cdot)$ | $\hat{\sigma}(\cdot)$ | $Sk(\cdot)$ | $Ku(\cdot)$
PBL | 3.3e-1 | -1.28e-1 | 26.51 | 3.2e-1 | 7.86e-2 | 26.67 | 3.0e-1 | 8.0e-3 | 22.75
PBM | 1.38 | -5.16e-1 | 19.92 | 1.41 | 1.57e-1 | 20.95 | 1.29 | -7.89e-2 | 19.5
PBH | 2.7 | 1.65e-1 | 14.86 | 2.79 | 1.3e-1 | 15.85 | 2.71 | -8.1e-2 | 15.2
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-squared fluctuations of
streamwise stress follow a logarithmic correlation as proposed by Örlü &
Schlatter (2011). Accordingly, a logarithmic correlation between the root
mean-squared fluctuations and the friction Reynolds number is also assumed in
the present permeable bed DNS,

$\displaystyle\tau^{t+}_{x,rms}=\tau_{x,rms}^{t}/\tau_{w}^{t}$
$\displaystyle=$ $\displaystyle 2.10\ln Re_{\tau}-8.11,$ (5)
$\displaystyle\tau^{t+}_{y,rms}=\tau_{y,rms}^{t}/\tau_{w}^{t}$
$\displaystyle=$ $\displaystyle 2.4\ln Re_{\tau}-9.83,$ (6)
$\displaystyle\tau^{t+}_{z,rms}=\tau_{z,rms}^{t}/\tau_{w}^{t}$
$\displaystyle=$ $\displaystyle 2.4\ln Re_{\tau}-10.10,$ (7)

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

Figure 10: Variation of the root-mean squared 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}$. Function of $Re_{\tau}$: DNS
data ( ), fitted logarithmic correlation from DNS data for permeable sediment
beds  , and Function of $Re_{K}$: DNS data ( ), fitted logarithmic correlation
from DNS data for permeable sediment beds  .

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 5–7 can then be used to
evaluate the bed stress variability for different $Re_{K}$. Together with the
non-Gaussian model distribution for the PDFs 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 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. Majority of the contribution to
the drag and lift forces comes from 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.

|  | Drag force | Lift force
---|---|---|---
| Case | %$F_{d}^{v}$ | %$F_{d}^{p}$ | ${\hat{\sigma}}(\cdot)$ | $Sk(\cdot)$ | $Ku(\cdot)$ | %$F_{\ell}^{v}$ | %$F_{\ell}^{p}$ | ${\hat{\sigma}}(\cdot)$ | $Sk(\cdot)$ | $Ku(\cdot)$
| PBL | 39.4 | 60.6 | 0.36 | -0.44 | 14.84 | 0.9 | 99.1 | 0.31 | 0.35 | 12.09
Full Bed | PBM | 23.7 | 76.3 | 1.59 | -0.88 | 14.49 | 9.9 | 90.1 | 1.26 | -0.15 | 11.15
| PBH | 16.7 | 83.25 | 2.82 | 0.46 | 11.20 | 3.2 | 96.8 | 2.83 | 0.58 | 8.46
| PBL | 36.8 | 63.2 | 0.6 | -0.35 | 6.79 | 20.1 | 79.9 | 0.62 | -3.0e-4 | 4.10
Top Layer | PBM | 22.15 | 77.85 | 2.53 | -0.63 | 5.35 | 19.3 | 80.7 | 2.5 | 0.11 | 3.47
| PBH | 13.4 | 86.6 | 4.50 | -0.18 | 4.42 | 4.6 | 95.4 | 4.99 | -0.08 | 3.81
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.

Figure 11: PDFs of drag and lift force fluctuations (a,c) drag force and (b,d)
lift force normalized by their respective standard deviation values. Top panel
shows permeable bed results for PBL ( ), PBM ( ), and PBH ( ) cases. Also
shown is a $t$-location scale model fit ( ). Bottom panel compares the values
for the full permeable bed (PBM,  ) and the top layer of the permeable wall
(PBMTL,  ) for the medium $Re_{K}$.

The probability distribution functions of the fluctuations of drag
($F^{\prime}_{d}$, $x-$component) and lift forces ($F^{\prime}_{\ell}$,
$y-$component) on all particles within the bed normalized by their respective
standard deviation values are shown in figures 11a,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.

Finally, to assess the contribution of the top layer on force statistics,
figures 11c,d compare the PDFs 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
modeling 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 improve prediction of reduced-order models.

###  3.6 Pressure distributions at 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 it can directly influence the residence times through turbulent
advection. Reach-scale modeling of hyporheic zone transport typically use 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., 2020). 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.

PDFs 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 sediment-water interface. Figures 12(a,b) show the pdf
distributions for the normalized turbulent pressure fluctuations,
$p^{\prime}$, and normalized form-induced pressure disturbances,
$\widetilde{\overline{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
12c shows the pdfs of the sum of the normalized turbulent and form-induced,
$p^{\prime}+\widetilde{\overline{p}}$, pressure PDFs for the three cases. The
pressure sum PDFs for all cases are statistically similar and symmetric,
slightly heavy-tailed than a Gaussian, however, the Gaussian distribution
nearly captures the pressure data within $\pm 3\hat{\sigma}$. This result
suggests that the pressure behavior 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 for $\widetilde{\overline{p}}$ is higher than
$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 3h.

Figure 12: PDFs of (a,d) turbulent pressure fluctuations ($p^{\prime}$), (b,e) form-induced pressure disturbances ($\widetilde{\overline{p}}$) (c,f) sum of $p^{\prime}$ and $\widetilde{\overline{p}}$. Top panel shows data for the permeable bed cases PBL ( ), PBM ( ), and PBH ( ), and a Gaussian model fit ( ). Bottom panel compares the PDFs for the full permeable bed (PBM) and the impermeable rough wall (IWF-M,  ).  | $p^{\prime}$ | $\widetilde{\overline{p}}$ | $\widetilde{\overline{p}}+p^{\prime}$
---|---|---|---
Case | $\hat{\sigma}(\cdot)$ | $Sk(\cdot)$ | $Ku(\cdot)$ | $\hat{\sigma}(\cdot)$ | $Sk(\cdot)$ | $Ku(\cdot)$ | $\hat{\sigma}(\cdot)$ | $Sk(\cdot)$ | $Ku(\cdot)$
PBL | 1.12 | -1.84e-1 | 4.22 | 1.60 | 1.70 | 8.04 | 1.86 | 1.09 | 6.76
PBM | 2.98 | -5.15e-1 | 5.68 | 2.38 | 2.07 | 9.76 | 3.68 | 0.48 | 6.65
PBH | 6.11 | -5.66e-1 | 7.52 | 4.27 | 1.87 | 8.77 | 7.54 | 0.68 | 8.52
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 PDFs 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 figures 12d–f. The variance in turbulent pressure
fluctuations and form-induced disturbances is generally larger in the IWM-F
case than those in the permeable bed case (PBM). The probability of higher
negative $\widetilde{\overline{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 $\widetilde{\overline{p}}$ show a reasonable
match between the PBM and IWM-F cases as shown in figure 12f. 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 sediment-water interface.

##  4 Conclusions

Pore-resolved direct numerical simulations 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, 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.

(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.

(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.

(iii) Mean flow penetration (Brinkman layer thickness) and shear penetration
show non-linear, increasing correlation with $Re_{K}$, their ratio approaching
a constant deterministic value. The length scale of dominant turbulent eddies
at the sediment-water interface is better represented by the mixing length
obtained from the Reynolds stress and mean flow gradient. The normalized
interfacial mixing length at the sediment-water interface increases with
$Re_{K}$ and shows similar behavior as 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.

(v) 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
differences between the permeable bed and impermeable rough wall. Form-induced
stresses are, however, influenced by the impermeability, redistributing
stresses tangential along the wall.

(vi) 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 PDFs of the sum of the pressure fluctuations and form-induced
pressure disturbances at the sediment-water interface are statistically
similar, symmetric, and collapse for different $Re_{K}$.

(vii) The PDFs 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 PDFs of
fluctuations in the bed stress normalized by their root-mean-squared values
are symmetric and exhibit a peaky, non-Gaussian distribution with heavy tails.
A logarithmic correlation between the root mean-squared 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.

(viii) The mean and fluctuations in drag and lift forces on the particle are
computed by integrating the local bed stress on the particle surface. 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 distribution of the net bed stress and the drag and lift forces on
particles are mainly influenced by the top layer of the sediment grain,
including the roughness effects through a single layer of randomly arranged
sediments can potentially improve reach-scale predictions based on reduced-
order models.

###### Acknowledgements.

## Acknowledgements

This work was initiated as part of SKK’s internship at Pacific Northwest
National Laboratory. Simulations were performed at the Texas Advanced
Computing Center’s (TACC) Frontera system. Computing resource from Pacific
Northwest National Laboratory’s EMSL (Environmental Molecular Sciences
Laboratory) is also acknowledged.

## Funding

SKK acknowledges support from Pacific Northwest National Laboratory (PNNL) as
part of an internship program. SVA and SKK gratefully acknowledge funding from
US Depart of Energy, Office of Basic Energy Sciences (Geosciences) under award
number DE-SC0021626 as well as US National Science Foundation award #205324.
The computing resources used were made available under NSF’s Leadership
Resources Allocation (LRAC) award. XH and TDS acknowledge funding from the DOE
Office of Biological and Environmental Research, Subsurface Biogeochemical
Research program, through the PNNL Subsurface Science Scientific Focus Area
project (http://sbrsfa.pnnl.gov/).

## Declaration of Interests

The authors report no conflict of interest.

## APPENDIX A. Grid refinement, integral scales, and domain sizes

The solver has been thoroughly verified and validated for a range of cases
(Apte & Patankar, 2008) and has also been used for large-scale, parallel
simulations of oscillatory flow over a layer of sediment particles (Ghodke &
Apte, 2016, 2018a) and flow through porous media (Finn & Apte, 2013; He et
al., 2019).

$Re_{D}$ | 50 | 150 | 350
---|---|---|---
| Drag Coefficient, $C_{D}$
Isotropic, $D_{p}/\Delta y=100$ | 1.54 | 0.9 | 0.65
$\Delta x=\Delta z=3\Delta y$ | 1.58 | 0.91 | 0.66
$\Delta x=\Delta z=4\Delta y$ | 1.56 | 0.9 | 0.66
Table 7: Grid refinement study for flow over an isolated sphere with non-
isotropic, rectilinear grids similar to those used in the present study.

In turbulent flow over a sediment bed, there is a need to use non-isotropic
and high-aspect ratio grids to minimize the total control volumes and yet
provide sufficient resolution needed to capture all scales of turbulence.
Specifically, for DNS of open channel flows, the resolution near the sediment
bed and in the bed-normal direction should be such that $\Delta y^{+}<1$,
where $\Delta y^{+}$ represents resolution in wall units. The code was used to
predict flow over an isolated sphere at different Reynolds numbers using
isotropic and non-isotropic, rectilinear grids. The drag force was compared
with published data (Apte & Patankar, 2008) and is given in Table 7. It is
observed that the high-aspect ratio grids are capable of predicting the drag
forces accurately for $Re_{D}$ up to 350, where $Re_{D}=UD_{p}/\nu$ is
Reynolds number based on the sphere diameter, $D_{p}$, and the uniform
undisturbed upstream velocity, $U$. Also, the effectiveness of the non-
isotropic grids in capturing vortex shedding was verified using the Strouhal
number for vortex shedding at $Re_{D}=350$. The obtained value for Strouhal
number was 0.131, which compared reasonably well with the range of values
between 0.135–0.14 predicted on finer, isotropic grids in literature (Mittal
et al., 2008; Mittal, 1999; Bagchi et al., 2001).

In the present permeable bed cases, the flow velocity near the bed is much
smaller than the free-stream velocity and hence the relevant Reynolds number
for the spherical roughness elements is the roughness Reynolds number
$D^{+}=D_{p}u_{\tau}/\nu$, where $D_{P}$ the particle diameter is a measure of
the roughness height for monodispersed particle beds. $D^{+}$ values for the
PBL, PBM and PBH cases are 77, 156, and 270, respectively, which fall well
within the range of Reynolds numbers of the above grid refinement study. The
grid resolution used in the present DNS is thus sufficient to capture the
inertial flow features around the particle, including flow separation and
wakes.

Eulerian two-point auto-correlations are used to compute the integral length
scales in streamwise and spanwise directions, and are compared with the domain
sizes in those directions. With increase in Reynolds number, the vortical
structures are broken down by the roughness elements, and the integral length
scales in the streamwise and spanwise directions are expected to decrease. The
time averaged Eulerian two-point auto-correlations were computed as

$\displaystyle\rho^{E}_{ij}\left(|\bm{s}|\right)$
$\displaystyle=\dfrac{\overline{\langle u^{\prime}_{i}(\bm{x},t)\
u^{\prime}_{j}(\bm{x+s},t)\rangle}}{\overline{\langle
u^{\prime}_{i}(\bm{x},t)\ u^{\prime}_{j}(\bm{x},t)\rangle}},$

here $\rho^{E}_{ij}$ is the Eulerian auto-correlation, $\bm{s}$ represents the
set of all possible vector displacements for which the auto-correlation is
calculated. The correlations were first computed at 100,000 randomly picked
locations $(\bm{x})$ in the fluid domain at one instant of time and then
spatially averaged to obtain the overall representation. This procedure is
then repeated over several flow through times to get a temporal average of the
spatially averaged values. The Eulerian integral length scales,
$\rho^{E}_{ij}$, are then calculated by integrating the correlations over the
respective abscissas and the results are presented in table 8. The length
scales for the PBL case are comparable to values obtained by Krogstad &
Antonia (1994); Shen et al. (2020) at similar Reynolds number. The domain
sizes for the PBL and PBM cases are approximately 11–13 times the integral
length scale in the streamwise and 20–32 times that in the spanwise direction.
The domain size for PBH is same as PBM. The domain lengths $L_{x}\times L_{y}$
are $12.56\delta\times 6.28\delta$ for PBL and $6.28\delta\times 3.14\delta$
for PBM are sufficient for the periodicity assumption to obtain mean and
turbulence statistics without any domain confinement effects. In comparison,
DNS of turbulent flow over rough, impermeable walls by Ma et al. (2021) used
$4\delta\times 2.4\delta$ for similar Reynolds numbers.

Case | $L_{11}/\delta$ | $L_{33}/\delta$ | $L_{11}/D_{p}$ | $L_{33}/D_{p}$ | $L_{x}/L_{11}$ | $L_{z}/L_{33}$ |
---|---|---|---|---|---|---|---
PBL | 0.953 | 0.2 | 3.33 | 0.7 | 13.18 | 31.45 |
PBM | 0.558 | 0.151 | 1.95 | 0.53 | 11.25 | 20.79 |
Table 8: Eulerian length scales in the streamwise ($\rho^{E}_{11}$) and
spanwise ($\rho^{E}_{33}$) directions normalized by the free-stream height,
$\delta$, and the particle diameter, $D_{P}$. Also shown are the streamwise
($L_{x}$) and spanwise ($L_{z}$) domain lengths normalized by $L_{11}$ and
$L_{33}$, respectively.

## APPENDIX B. Variable volume averaging

Pokrajac & De Lemos (2015) used a variable volume averaging methodology, to
spatially average time-averaged quantities, wherein thinner volumes were used
to average the flow in the free-stream and thicker volumes were used to
average inside the porous region with a smooth transition in volume height
between the two regions. Following a similar approach, in this study thin
volumes are used for averaging in the free-stream region near the crest of the
bed where steep gradients in flow quantities are present. The averaging volume
is gradually coarsened away from this region and deeper inside the bed thicker
averaging volumes are used. For averaging purposes the domain in the vertical
direction is divided into four regions: (Region 1) uniform averaging-volume
height deep inside the bed, (Region 2) a transitioning averaging-volume height
between the uniform region below and the top layer of the bed, (Region 3) a
refined uniform averaging-volume height in the top layer and the crest region,
(Region 4) a transitioning averaging-volume height in the free-stream region.
The variable volume averaging approach across the various segments is given as

$\displaystyle l_{y}$
$\displaystyle=l_{1}\eta_{1}\qquad\qquad\qquad\qquad\qquad\quad\textrm{Region
1}\qquad-1.14\leq y/\delta\leq-0.57$ $\displaystyle l_{y}$
$\displaystyle=\frac{l_{2}\tanh{\left(\gamma_{1}\eta_{2}\right)}}{\tanh{\left(\gamma_{1}\right)}}\qquad\qquad\qquad\
\ \textrm{Region 2}\qquad-0.57<y/\delta\leq-0.28$ $\displaystyle l_{y}$
$\displaystyle=l_{3}\eta_{3}\qquad\qquad\qquad\qquad\qquad\quad\textrm{Region
3}\qquad-0.28<y/\delta\leq 0.031$ $\displaystyle l_{y}$
$\displaystyle=\frac{l_{4}\left(1-\tanh{\left(\gamma_{2}-\gamma_{2}\eta_{3}\right)}\right)}{\tanh{\left(\gamma_{2}\right)}}\qquad\textrm{Region
4}\qquad\quad\ 0.031<y/\delta\leq 1,$

where $l_{1},l_{2},l_{3}\text{ and }l_{4}$ are the vertical heights of each
region, $\gamma_{1}$ and $\gamma_{2}$ control the rate of transitioning of the
averaging volume height in the bed-normal direction,
$\eta_{1}=\eta/w_{1},\eta_{2}=(\eta-
w_{1})/w_{2},\eta_{3}=(\eta-(w_{1}+w_{2}))/w_{3}$, and $\eta=j/gnj$. Here,
values $w$ are weights based on the ratio of number of assigned volumes for
averaging in each region and the total number of volumes, $j$ is the index of
the averaging volume, $gnj$ is the total number of averaging volumes used in
the bed normal direction. Typical $\gamma_{1}$ and $\gamma_{2}$ values are
between $1.5-3$ and $0.7-1.3$, respectively. Similar variable volume averaging
has been carried out in previous studies by Karra et al. (2022b, a)

## Appendix C: Validation Study

Pore-resolved direct numerical simulation of turbulent flow over a sediment
bed (case VV) was first validated with experimental data of Voermans et al.
(2017) as well as DNS data of Shen et al. (2020). Permeable bed case with
porosity of $0.41$, $Re_{K}=2.56$ and $Re_{\tau}\sim 180$ matches with case
L12 in Voermans et al. (2017). In the present work, the numerical algorithm
developed by Dye et al. (2013) is used to generate a random distribution of
monodispersed spheres for a given porosity. It uses the collective
rearrangement algorithm introduced by Williams & Philipse (2003), coupled with
a mechanism for controlling the overall system porosity, providing a periodic
arrangement in the streamwise and spanwise directions. Although the average
porosity in the sediment bed is matched, the actual random configuration of
the sediment particles is likely different compared to both the experimental
(Voermans et al., 2017) and published DNS data (Shen et al., 2020).

Figure 13: Schematic showing positions of the sediment crest ($y=0$), the
zero-displacement plane ($y=-d$), and particle diameter ($D_{p}$).

Voermans et al. (2017) defined the origin of the sediment bed to be the
inflection point in the porosity profile, that is, where
$\partial^{2}_{yy}\phi=0$. Therefore, in order to facilitate comparison of the
current DNS results with the experimental work the origin for case VV is taken
to be the inflection point ($\partial^{2}_{yy}\phi=0$ ) of its porosity
profile. Shown also in the schematic in figure 13 is the zero-displacement
plane, $y=-d$, whose physical meaning and values are given in Appendix D. The
time-space averaged mean velocity profile normalized by channel free-stream
velocity, $U_{\delta}$, is shown in figure 14a. Excellent agreement is seen
between the DNS data, experimental measurements and DNS data from Shen et al.
(2020). Figures 14b, 14c, and 14d show a comparison of turbulence intensities,
namely streamwise, bed-normal and shear stresses. Again very good agreement
between DNS and experiment is observed. The slight deviation in Reynolds
stress between the DNS and experimental in the outer channel flow region can
be attributed to the high measurement uncertainty (between $6-30\%$ Voermans
et al. (2017)) in sampling this variable in the experiment.

Figure 14: Comparison of (a) mean streamwise, velocity and (b) streamwise, (c)
wall-normal, and (d) shear components of spatially averaged Reynolds stress
tensor. Experimental data by Voermans et al. (2017) ( ), DNS by Shen et al.
(2020) ( ), present DNS ( ).

Figure 15: Comparison of (a) streamwise, (b) wall-normal, and (c) shear
components of form induced stress tensor. Experimental data by Voermans et al.
(2017) ( ), emulating experimental sampling ( ), Shen et al. (2020) ( ),
present DNS ( ).

The normalized form-induced or dispersive stresses are shown in figures 15a,
15b, and 15c. The differences between the present DNS and the experimental
results can be explained based on the sampling procedures used. Firstly, as
mentioned in the previous section, spatial averaging is carried out over an
entire $x-z$ volume at a given $y$ location for DNS results. However, for the
experimental data, spatial averaging was performed over three different
spanwise locations over six different measurements. To quantify the
differences in the sampling procedures between the experiments and DNS, the
experimental sampling process is replicated in the DNS data whereby spatial
averaging is carried out at a few finite uncorrelated spanwise locations and
repeated over different streamwise locations. A family of curves, shown by
grey squares, indicates the associated uncertainty with the sampling locations
of the experimental data. The averaged experimental and DNS data are within
this scatter for all streamwise locations. Secondly, it is has been reported
in literature (Nikora et al., 2002; Fang et al., 2018) that the spanwise
averaging is highly sensitive to the geometry at the sediment-water interface.
For the present DNS, only the mean porosity of the randomly distributed
arrangement of monodispersed spherical particles is matched with the
experimental geometry. However, the exact sediment-grain distribution in the
experiments is unknown and is likely different compared to that used in the
DNS. This difference, especially near the top of the bed can also contribute
to differences in the form-induced or dispersive stresses.

In spite of the potential differences in the sediment bed distribution between
DNS and experimental work, the present results reproduce the mean flow and
turbulence stresses observed in the experiment. The form-induced stresses fall
within the uncertainty associated with sampling locations in the experiments.
In addition, turbulence statistics from the current work are compared with DNS
predictions from Shen et al. (2020). Good agreement between the two sets of
DNS results is observed. The consistency in predicted results with the
published experimental and numerical studies persuasively validates the
numerical approach used in this work.

## D. The log-law and zero-displacement thickness

In turbulent flows over rough walls and permeable beds, the log-law has the
following form

$\displaystyle\dfrac{U(y)}{u_{\tau}}$
$\displaystyle=\dfrac{1}{\kappa}\log\left(\dfrac{y+d}{y_{0}}\right),$

where $\kappa$ is the von-Kármán constant, $d$ is distance between the zero-
displacement plane and the top of the sediment crest (see figure 13), and
$y_{0}$ is the equivalent roughness height which is related to the measure of
the size of the roughness elements.

Case | $\kappa$ | $d/{\delta}$ | $d/D_{p}$ | $d^{+}$ | $y_{0}/\delta$ | $y_{0}/D_{p}$ | $y_{0}^{+}$ | $\delta_{b}/D_{p}$ | $\delta_{p}/D_{p}$
---|---|---|---|---|---|---|---|---|---
PBL | 0.325 | 0.1686 | 0.59 | 45 | 0.0244 | 0.085 | 6.57 | 0.76 | 0.76
PBM | 0.31 | 0.1743 | 0.61 | 95 | 0.0298 | 0.10 | 16.22 | 0.94 | 0.94
PBH | 0.2875 | 0.1829 | 0.64 | 172 | 0.0364 | 0.127 | 34.35 | 1.31 | 1.11
Table 9: The von-Kármán constant ($\kappa$), zero-displacement thickness
($d$), and equivalent roughness height ($y_{0}$) normalized by $\delta$,
$D_{p}$ and $\nu/u_{\tau}$. The last two columns shown Brinkman layer
thickness $\delta_{b}$, and shear penetration depth $\delta_{p}$.
$(\hphantom{0})^{+}$ denotes wall units. Results are shown for PBL, PBM, and
PBH cases.

Although several techniques have been used to determine these parameters in
literature (Raupach et al., 1991), the procedure described by Breugem et al.
(2006) is followed here. First, the extent of the logarithmic layer is
determined by plotting $(y+d)^{+}{\partial_{y^{+}}U^{+}}$ against $y^{+}$ for
several values of $d$. From the equation of log-law, it is easy to see that
the value of $(y+d)^{+}{\partial_{y^{+}}U^{+}}$ is a constant equal to
$1/\kappa$ in the logarithmic layer. Therefore, the value of $d$ is the one
that gives a horizontal profile in the logarithmic layer. The values of $d$,
$\kappa$, and $y_{0}$ determined from a least squares fit of log-law equation
to the velocity profile in the logarithmic layer are given in table 9. The
von-Kármán constant ($\kappa$) for the three permeable bed cases is lower than
the value of 0.4 for flows over smooth walls. This decrease in $\kappa$ has
also been observed in flows over permeable walls by Suga et al. (2010); Manes
et al. (2011); Breugem et al. (2006). Both the zero-displacement thickness and
equivalent roughness height show a dependency on $Re_{K}$ and increase with
increasing Reynolds number. Their values in wall units, $d^{+}$ and
$y_{0}^{+}$, compare reasonably well with the studies by Manes et al. (2011);
Suga et al. (2010); Breugem et al. (2006).

From Nikuradse (1933) experiments of flows over impermeable fully rough walls,
the ratio of $y_{0}/D_{p}$ was found be approximately $1/30$. The values of
$y_{0}/D_{p}$ for the PBL, PBM, and PBH cases given in table 9 are
approximately 2-4 times larger than the value observed in fully rough walls.
Importantly, they show a correlation with $Re_{K}$, due to the influence of
the bed permeability. Hinze (1975) (p. 637) report that for fully rough
impermeable walls the ratio of $d/D_{p}$ is approximately $0.3$. This ratio
for the three permeable bed cases is roughly 2 times larger as shown in table
9 and also shows a dependency on $Re_{K}$. Therefore, the permeable bed cases
in the current study show the characteristics of a fully rough wall regime
(based on their roughness Reynolds number ($D^{+}>70$) as shown in table 1)
influenced by permeability. Breugem et al. (2006) in their study for cases
with $1<Re_{K}<10$ found the values of $y_{0}/D_{p}$ and $d/D_{p}$ to be
orders of magnitude greater than $1/30$ and $0.3$ respectively as their
$D^{+}$ values were $<7$ which meant that the effects of surface roughness
were negligible.

## References

  * Anderson et al. (2008) Anderson, Mary P & others 2008 Groundwater: Selection, Introduction and Commentary.. IAHS Press.
  * Apte & Patankar (2008) Apte, SV & Patankar, NA 2008 A formulation for fully resolved simulation (FRS) of particle–turbulence interactions in two-phase flows. International Journal of Numerical Analysis and Modeling 5, Suppl, 1–16.
  * Apte et al. (2009) Apte, S. V., Martin, M. & Patankar, N. A. 2009 A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows. Journal of Computational Physics 228 (8), 2712–2738.
  * Bagchi et al. (2001) Bagchi, Prosenjit, Ha, MY & Balachandar, S 2001 Direct numerical simulation of flow and heat transfer from a sphere in a uniform cross-flow. J. Fluids Eng. 123 (2), 347–358.
  * Bencala et al. (1983) Bencala, KE, Rathbun, RE, Jackman, A P, Kennedy, VC, Zellweger, GW & Avanzino, RJ 1983 Rhodamine wt dye losses in a mountain stream environment. Water Resources Bulletin 19 (6), 943–950.
  * Breugem et al. (2006) Breugem, WP, Boersma, BJ & Uittenbogaard, RE 2006 The influence of wall permeability on turbulent channel flow. Journal of Fluid Mechanics 562, 35.
  * Breugem & Boersma (2005) Breugem, Wim-Paul & Boersma, Bendiks-Jan 2005 Direct numerical simulations of turbulent flow over a permeable wall using a direct and a continuum approach. Physics of fluids 17 (2), 025103\.
  * Briggs et al. (2009) Briggs, Martin A, Gooseff, Michael N, Arp, Christopher D & Baker, Michelle A 2009 A method for estimating surface transient storage parameters for streams with concurrent hyporheic storage. Water Resources Research 45 (4).
  * Chen et al. (2018) Chen, Xiaobing, Cardenas, M Bayani & Chen, Li 2018 Hyporheic exchange driven by three-dimensional sandy bed forms: Sensitivity to and prediction from bed form geometry. Water Resources Research 54 (6), 4131–4149.
  * Chen et al. (2020) Chen, Yunxiang, Bao, Jie, Fang, Yilin, Scheibe, Timothy D, Richmond, Marshall C, Perkins, William A, Ren, Huiying, Song, Xuehang, Duan, Zhuoran, Hou, Zhangshuan & others 2020 Modeling the effects of topography heterogeneity and discharge variations on riverbed hydrodynamics in a 30-kilometer-long reach over a nine-year period using openfoam .
  * D’angelo et al. (1993) D’angelo, DJ, Webster, JR, Gregory, SV & Meyer, JL 1993 Transient storage in appalachian and cascade mountain streams as related to hydraulic characteristics. Journal of the North American Benthological Society pp. 223–235.
  * Dye et al. (2013) Dye, Amanda L, McClure, James E, Miller, Cass T & Gray, William G 2013 Description of non-darcy flows in porous medium systems. Physical Review E 87 (3), 033012.
  * Fang et al. (2018) Fang, Hongwei, Xu, Han, He, Guojian & Dey, Subhasish 2018 Influence of permeable beds on hydraulically macro-rough flow. Journal of Fluid Mechanics 847, 552–590.
  * Finn & Apte (2013) Finn, Justin & Apte, Sourabh V 2013 Relative performance of body fitted and fictitious domain simulations of flow through fixed packed beds of spheres. International journal of multiphase flow 56, 54–71.
  * Finn (2013) Finn, Justin Richard 2013 A numerical study of inertial flow features in moderate Reynolds number flow through packed beds of spheres. PhD thesis, Oregon State University.
  * Ghisalberti (2009) Ghisalberti, Marco 2009 Obstructed shear flows: similarities across systems and scales. Journal of Fluid Mechanics 641, 51–61.
  * Ghodke & Apte (2016) Ghodke, Chaitanya D & Apte, Sourabh V 2016 Dns study of particle-bed–turbulence interactions in an oscillatory wall-bounded flow. Journal of Fluid Mechanics 792, 232–251.
  * Ghodke & Apte (2018a) Ghodke, Chaitanya D & Apte, Sourabh V 2018a Roughness effects on the second-order turbulence statistics in oscillatory flows. Computers & Fluids 162, 160–170.
  * Ghodke & Apte (2018b) Ghodke, Chaitanya D & Apte, Sourabh V 2018b Spatio-temporal analysis of hydrodynamic forces on the particle bed in an oscillatory flow environment. Journal of Fluid Mechanics 841, 167.
  * Grant et al. (2018) Grant, Stanley B., Gomez-Velez, Jesus D. & Ghisalberti, Marco 2018 Modeling the effects of turbulence on hyporheic exchange and local-to-global nutrient processing in streams. Water Resources Research 54 (9), 5883–5889.
  * Harvey et al. (1996) Harvey, Judson W, Wagner, Brian J & Bencala, Kenneth E 1996 Evaluating the reliability of the stream tracer approach to characterize stream-subsurface water exchange. Water Resources Research 32 (8), 2441–2451.
  * He et al. (2018) He, Xiaoliang, Apte, Sourabh, Schneider, Kai & Kadoch, Benjamin 2018 Angular multiscale statistics of turbulence in a porous bed. Physical Review Fluids 3 (8), 084501.
  * He et al. (2019) He, Xiaoliang, Apte, Sourabh V, Finn, Justin R & Wood, Brian D 2019 Characteristics of turbulence in a face-centred cubic porous unit cell. Journal of Fluid Mechanics 873, 608–645.
  * Hester et al. (2017) Hester, Erich T, Cardenas, M Bayani, Haggerty, Roy & Apte, Sourabh V 2017 The importance and challenge of hyporheic mixing. Water Resources Research 53 (5), 3565–3575.
  * Hinze (1975) Hinze, JO 1975 Turbulence. mcgraw-hill publishing co. New York .
  * Jackson et al. (2013a) Jackson, TR, Haggerty, R & Apte, SV 2013a A fluid-mechanics based classification scheme for surface transient storage in riverine environments: quantitatively separating surface from hyporheic transient storage. Hydrol. Earth Syst. Sci 17, 2747–2779.
  * Jackson et al. (2015) Jackson, Tracie R, Apte, Sourabh V, Haggerty, Roy & Budwig, Ralph 2015 Flow structure and mean residence times of lateral cavities in open channel flows: influence of bed roughness and shape. Environmental Fluid Mechanics 15 (5), 1069–1100.
  * Jackson et al. (2013b) Jackson, Tracie R, Haggerty, Roy, Apte, Sourabh V & O’Connor, Ben L 2013b A mean residence time relationship for lateral cavities in gravel-bed rivers and streams: Incorporating streambed roughness and cavity shape. Water resources research 49 (6), 3642–3650.
  * Jeon et al. (1999) Jeon, Sejeong, Choi, Haecheon, Yoo, Jung Yul & Moin, Parviz 1999 Space–time characteristics of the wall shear-stress fluctuations in a low-reynolds-number channel flow. Physics of fluids 11 (10), 3084–3094.
  * Jiménez (2004) Jiménez, Javier 2004 Turbulent flows over rough walls. Annu. Rev. Fluid Mech. 36, 173–196.
  * Karra et al. (2022a) Karra, Shashank, Apte, Sourabh V, He, Xiaoliang & Scheibe, Timothy D 2022a Particle resolved dns study of turbulence effects on hyporheic mixing in randomly packed sediment beds. In 12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12), Osaka, Japan, July 19–22, 2022.
  * Karra et al. (2022b) Karra, Shashank K, Apte, Sourabh V, He, Xiaoliang & Scheibe, Timothy D 2022b Pore-resolved simulations of turbulent boundary layer flow over permeable and impermeable sediment beds. arXiv preprint arXiv:2204.13875 .
  * Kim et al. (2020) Kim, Taehoon, Blois, Gianluca, Best, James L & Christensen, Kenneth T 2020 Experimental evidence of amplitude modulation in permeable-wall turbulence. Journal of Fluid Mechanics 887.
  * Krogstad & Antonia (1994) Krogstad, P-Å & Antonia, RA 1994 Structure of turbulent boundary layers on smooth and rough walls. Journal of Fluid Mechanics 277, 1–21.
  * Kuwata & Suga (2016) Kuwata, Y & Suga, K 2016 Lattice boltzmann direct numerical simulation of interface turbulence over porous and rough walls. International Journal of Heat and Fluid Flow 61, 145–157.
  * Lee et al. (1990) Lee, Moon Joo, Kim, John & Moin, Parviz 1990 Structure of turbulence at high shear rate. Journal of Fluid Mechanics 216, 561–583.
  * Leonardi et al. (2018) Leonardi, Alessandro, Pokrajac, D, Roman, F, Zanello, F & Armenio, V 2018 Surface and subsurface contributions to the build-up of forces on bed particles. Journal of Fluid Mechanics 851, 558–572.
  * Ma et al. (2021) Ma, Rong, Alamé, Karim & Mahesh, Krishnan 2021 Direct numerical simulation of turbulent channel flow over random rough surfaces. Journal of Fluid Mechanics 908.
  * Manes et al. (2011) Manes, Costantino, Poggi, Davide & Ridolfi, Luca 2011 Turbulent boundary layers over permeable walls: scaling and near-wall structure. Journal of Fluid Mechanics 687, 141–170.
  * Manes et al. (2009) Manes, Costantino, Pokrajac, Dubravka, McEwan, Ian & Nikora, Vladimir 2009 Turbulence structure of open channel flows over permeable and impermeable beds: A comparative study. Physics of Fluids 21 (12), 125109.
  * Manes et al. (2012) Manes, C, Ridolfi, L & Katul, G 2012 A phenomenological model to describe turbulent friction in permeable-wall flows. Geophysical research letters 39 (14).
  * Mittal (1999) Mittal, Rajat 1999 A fourier–chebyshev spectral collocation method for simulating flow past spheres and spheroids. International journal for numerical methods in fluids 30 (7), 921–937.
  * Mittal et al. (2008) Mittal, Rajat, Dong, Haibo, Bozkurttas, Meliha, Najjar, FM, Vargas, Abel & Von Loebbecke, Alfred 2008 A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. Journal of computational physics 227 (10), 4825–4852.
  * Moser et al. (1999) Moser, Robert D, Kim, John & Mansour, Nagi N 1999 Direct numerical simulation of turbulent channel flow up to re $\tau$= 590\. Physics of fluids 11 (4), 943–945.
  * Nikora et al. (2013) Nikora, Vladimir, Ballio, Francesco, Coleman, Stephen & Pokrajac, Dubravka 2013 Spatially averaged flows over mobile rough beds: Definitions, averaging theorems, and conservation equations. Journal of Hydraulic Engineering .
  * Nikora et al. (2004) Nikora, Vladimir, Koll, Katinka, McEwan, Ian, McLean, Stephen & Dittrich, Andreas 2004 Velocity distribution in the roughness layer of rough-bed flows. Journal of Hydraulic Engineering 130 (10), 1036–1042.
  * Nikora et al. (2007) Nikora, Vladimir, McEwan, Ian, McLean, Stephen, Coleman, Stephen, Pokrajac, Dubravka & Walters, Roy 2007 Double-averaging concept for rough-bed open-channel and overland flows: Theoretical background. Journal of hydraulic Engineering 133 (8), 873–883.
  * Nikora et al. (2002) Nikora, Vladimir Ivanovich, Koll, K, McLean, SR, Ditrich, A & Aberle, J 2002 Zero-plane displacement for rough-bed open-channel flows. In Fluvial Hydraulics River Flow 2002, pp. 83–92.
  * Nikuradse (1933) Nikuradse, Johann 1933 Stromungsgesetze in rauhen rohren. vdi-forschungsheft 361, 1.
  * Örlü & Schlatter (2011) Örlü, Ramis & Schlatter, Philipp 2011 On the fluctuating wall-shear stress in zero pressure-gradient turbulent boundary layer flows. Physics of fluids 23 (2), 021704.
  * Pokrajac & De Lemos (2015) Pokrajac, Dubravka & De Lemos, MJS 2015 Spatial averaging over a variable volume and its application to boundary-layer flows over permeable walls. Journal of Hydraulic Engineering 141 (4), 04014087.
  * Pokrajac et al. (2006) Pokrajac, Dubravka, Finnigan, JJ, Manes, C, McEwan, I & Nikora, V 2006 On the definition of the shear velocity in rough bed open channel flows. In River flow, , vol. 1, pp. 89–98.
  * Raupach et al. (1991) Raupach, MR, Antonia, RA & Rajagopalan, S 1991 Rough-wall turbulent boundary layers. Applied Mechanics Reviews 44 (1), 1.
  * Rousseau & Ancey (2022) Rousseau, G & Ancey, C 2022 An experimental investigation of turbulent free-surface flows over a steep permeable bed. Journal of Fluid Mechanics 941, A51.
  * Shen et al. (2020) Shen, Guangchen, Yuan, Junlin & Phanikumar, Mantha S 2020 Direct numerical simulations of turbulence and hyporheic mixing near sediment–water interfaces. Journal of Fluid Mechanics 892.
  * Suga et al. (2010) Suga, K, Matsumura, Y, Ashitaka, Y, Tominaga, S & Kaneda, M 2010 Effects of wall permeability on turbulence. International Journal of Heat and Fluid Flow 31 (6), 974–984.
  * Tennekes et al. (1972) Tennekes, Hendrik, Lumley, John Leask, Lumley, Jonh L & others 1972 A first course in turbulence. MIT press.
  * Valett et al. (1996) Valett, H Maurice, Morrice, John A, Dahm, Clifford N & Campana, Michael E 1996 Parent lithology, surface–groundwater exchange, and nitrate retention in headwater streams. Limnology and oceanography 41 (2), 333–345.
  * Voermans et al. (2017) Voermans, JJ, Ghisalberti, M & Ivey, GN 2017 The variation of flow and turbulence across the sediment–water interface. Journal of Fluid Mechanics 824, 413–437.
  * Voermans et al. (2018) Voermans, Joey J, Ghisalberti, Marco & Ivey, Gregory N 2018 A model for mass transport across the sediment-water interface. Water Resources Research 54 (4), 2799–2812.
  * Williams & Philipse (2003) Williams, SR & Philipse, AP 2003 Random packings of spheres and spherocylinders simulated by mechanical contraction. Physical Review E 67 (5), 051301.
  * Wilson et al. (2008) Wilson, Alicia M, Huettel, Markus & Klein, Stephen 2008 Grain size and depositional environment as predictors of permeability in coastal marine sands. Estuarine, Coastal and Shelf Science 80 (1), 193–199.
  * Zagni & Smith (1976) Zagni, Anthony FE & Smith, Kenneth VH 1976 Channel flow over permeable beds of graded spheres. Journal of the hydraulics division 102 (2), 207–222.
  * Zippe & Graf (1983) Zippe, Hans J & Graf, Walter H 1983 Turbulent boundary-layer flow over permeable and non-permeable rough surfaces. Journal of Hydraulic research 21 (1), 51–65.