1. Introduction
Steep streams are waterways whose slope exceeds a few per cent (Church Reference Church2010; Comiti & Mao Reference Comiti and Mao2012; Buffington & Montgomery Reference Buffington and Montgomery2013). Their bed usually contains coarse particles (cobbles and boulders whose size exceeds 10 cm). Water velocity is high (normally in excess of 1 
${\rm m}\,{\rm s}^{-1}$), while the flow is shallow relative to the bed roughness (typically, the flow depth is of the same order as the bed's coarsest particles) and reaches a supercritical regime most of the time. Turbulence is more strongly affected by bed roughness and varied bedforms (e.g. cascades, steps and pools, riffles) in steep streams than in low-gradient streams: a steep stream's velocity profile rarely exhibits a logarithmic layer (Manes, Pokrajac & McEwan Reference Manes, Pokrajac and McEwan2007), flow resistance is usually much higher (Ferguson Reference Ferguson2013), turbulent kinetic energy is more easily transported from the surface to the subsurface flow (Manes et al. Reference Manes, Pokrajac, McEwan and Nikora2009), and turbulent fluctuations in the near-bed layer are less pronounced (Lamb, Brun & Fuller Reference Lamb, Brun and Fuller2017a; Cooper et al. Reference Cooper, Ockleford, Rice and Powell2018). Steep-stream bed topography and high bed permeability promote momentum and mass exchanges with the hyporheic flow through the bed (Buffington & Tonina Reference Buffington and Tonina2009; Tonina & Buffington Reference Tonina and Buffington2009; Boano et al. Reference Boano, Harvey, Marion, Packman, Revelli, Ridolfi and Wörman2014).
Scale invariance is a fundamental concept that has been used widely in hydraulics (Heller Reference Heller2017). Among other things, it implies that a river's overall dynamics can be determined from certain macroscopic variables (such as bed slope and flow depth) using power-law equations. The Manning–Strickler law, used for modelling flow resistance, provides a typical example of scale invariance; this empirical equation holds for a wide range of flow conditions and stream configurations, and it has been justified on dimensional grounds and connected to Kolmogorov's theory of turbulence (Gioia & Bombardelli Reference Gioia and Bombardelli2002; Bonetti et al. Reference Bonetti, Manoli, Manes, Porporato and Katul2017; Katul, Li & Manes Reference Katul, Li and Manes2019). When applied to steep streams, it underestimates flow resistance by one order of magnitude (Ferguson Reference Ferguson2013; Powell Reference Powell2014). Various equations have thus been developed specifically to deal with flow resistance in steep streams by proposing new power-law equations adjusted on steep-stream data (Rickenmann Reference Rickenmann2016), but because of the wide range of temporal and spatial scales associated with these streams, flow resistance data are scattered, and no unique scaling law emerges from them. There is thus growing evidence that steep streams exhibit no scale invariance (Ferguson Reference Ferguson2021), and thus, faced with this failure of scale invariance, we need to dig down into the issue of spatial scales by looking into their flow dynamics at a mesoscopic scale (that is, at the bed roughness scale).
Steep streams are a particular case of flows of Newtonian fluid over a porous interface. Models describing these flows can be split into three classes depending on the scale of observation (Chandesris & Jamet Reference Chandesris and Jamet2006). Microscopic models resolve the flow on a length scale finer than the mean bed particle size. They usually involve direct numerical simulations of the Navier–Stokes equations or large-eddy simulations (Breugem, Boersma & Uittenbogaard Reference Breugem, Boersma and Uittenbogaard2006; Fang et al. Reference Fang, Han, He and Dey2018; Chen et al. Reference Chen, DiBiase, McCarroll and Liu2019; Kuwata & Kawaguchi Reference Kuwata and Kawaguchi2019; Shen, Yuan & Phanikumar Reference Shen, Yuan and Phanikumar2020). At the other end of the observation scale, macroscopic models assume that the free flow and porous bed are two distinct continua separated by a clear interface. Flow dynamics in each layer can be described using mass and momentum balance equations supplemented by relevant jump conditions at the interface (Goyeau et al. Reference Goyeau, Lhuillier, Gobin and Velarde2003; Jamet & Chandesris Reference Jamet and Chandesris2009). This is a two-part problem: determining how turbulence is affected by the rough, porous interface, and relating the jump conditions to bulk flow conditions. An example of the macroscopic approach applied to steep streams was given by Lamb, Brun & Fuller (Reference Lamb, Brun and Fuller2017b), who developed a model in which the surface flow is described using Boussinesq's eddy-viscosity equation, and the subsurface flow is modelled using the Darcy–Forchheimer equation, including a Brinkman correction. These equations are closed using ad hoc algebraic equations.
In the third approach, referred to as mesoscopic, the system is regarded as a two-phase medium. The governing equations for the fluid phase are obtained by averaging the Navier–Stokes equations in time and space over a control volume aligned with the flow direction (Nikora et al. Reference Nikora, Goring, McEwan and Griffiths2001, Reference Nikora, McEwan, McLean, Coleman, Pokrajac and Walters2007a; Dey & Das Reference Dey and Das2012). This technique, called double-averaging, combines temporal averaging (used in turbulence to derive the Reynolds-averaged Navier–Stokes equation) and spatial averaging (used to derive the governing equations of flows in porous media) (Whitaker Reference Whitaker1999). The presence of the solid phase is reflected by the porosity function 
$\epsilon$ (defined as the ratio of the fluid volume to the total volume in the control volume). The surface and subsurface flows are separated by a transition layer within which 
$\epsilon$ varies from a mean bulk porosity to unity. Compared to the macroscopic approach, mesoscopic models require no jump conditions at the interface between surface and subsurface flows. In addition to the Reynolds stress tensor (accounting for turbulence), mesoscopic models involve the dispersive stress tensor (representing the momentum transfer induced by the spatial heterogeneities in the velocity field created by bed obstacles). Both tensors need to be closed. Although increasingly complex procedures for closing the Reynolds stress tensor have been proposed, and several recent numerical studies have focused on dispersive stresses (Fang et al. Reference Fang, Han, He and Dey2018; Jelly & Busse Reference Jelly and Busse2018; Kuwata & Kawaguchi Reference Kuwata and Kawaguchi2019; Shen et al. Reference Shen, Yuan and Phanikumar2020), we know of few attempts to date to close the dispersive stress tensor (Kuwata & Suga Reference Kuwata and Suga2013, Reference Kuwata and Suga2015).
The double-averaging technique has been applied to study, for instance, turbulence in gravel-bed rivers (Nikora et al. Reference Nikora, Koll, McEwan, McLean and Dittrich2004, Reference Nikora, McLean, Coleman, Pokrajac, McEwan, Campbell, Aberle, Clunie and Koll2007b; Manes et al. Reference Manes, Pokrajac and McEwan2007; Franca, Ferreira & Lemmin Reference Franca, Ferreira and Lemmin2008; Mignot, Barthélemy & Hurther Reference Mignot, Barthélemy and Hurther2009; Dey, Sarkar & Ballio Reference Dey, Sarkar and Ballio2011; Cameron, Nikora & Stewart Reference Cameron, Nikora and Stewart2017; Papadopoulos et al. Reference Papadopoulos, Nikora, Vowinckel, Cameron, Jain, Stewart, Gibbins and Fröhlich2020), mass and momentum exchanges across the bed–stream interface (Giménez-Curto & Corniero Reference Giménez-Curto and Corniero2002; Manes et al. Reference Manes, Pokrajac, McEwan and Nikora2009; Voermans, Ghisalberti & Ivey Reference Voermans, Ghisalberti and Ivey2017, Reference Voermans, Ghisalberti and Ivey2018), the Darcy–Weisbach friction factor for open-channel flows (Nikora et al. Reference Nikora, Stoesser, Cameron, Stewart, Papadopoulos, Ouro, McSherry, Zampiron, Marusic and Falconer2019), and turbulence modification in sediment-laden flows (Revil-Baudard et al. Reference Revil-Baudard, Chauchat, Hurther and Eiff2016). We applied this technique to study experimentally the dynamics of steep streams and interpret the flow data that we acquired in a tilted flume. We revisited the problem initially addressed by Lamb et al. (Reference Lamb, Brun and Fuller2017b) at the macroscopic scale. However, working on the mesoscopic scale required that we be able to capture the velocity field within the entire flow domain. To that end, we used a technique called refractive index matching, which involved using fluid and particles with the same refractive index (Wiederseiner et al. Reference Wiederseiner, Andreini, Épely-Chauvin and Ancey2011; Dijksman et al. Reference Dijksman, Rietz, Lörincz and van Hecke2012). This experimental technique has been used increasingly in recent years to reconstruct the three-dimensional velocity fields far from the sidewalls of free-surface flows (Voermans et al. Reference Voermans, Ghisalberti and Ivey2017; Ni & Capart Reference Ni and Capart2018; Kim et al. Reference Kim, Blois, Best and Christensen2020; Shih & Wu Reference Shih and Wu2020; Trewhela & Ancey Reference Trewhela and Ancey2021).
This paper focuses on the steady uniform flow in a turbulent supercritical regime over a sloping permeable bed. In § 2, we define the physical system that was studied experimentally, we recall how the momentum balance equation for a fluid phase can be obtained using the double-averaging method, and we introduce the closure equations. We emphasise algebraic closure equations based on the mixing-length concept. Section 3 is devoted to the experimental protocol. Although the experimental set-up was not in complete similarity with flows observed in real-world mountain streams, it retained their main features: a steep slope, a permeable bed and a shallow flow in a supercritical regime. Section 4 presents the different experimental profiles used (velocity and turbulent and dispersive shear stresses) and compares them with the theoretical model's predictions. Section 5 discusses the study's limitations and proposes possible extensions to the model. An electronic supplement (available at https://doi.org/10.1017/jfm.2022.310) accompanies the paper and provides additional proofs for § 2.2.
2. Theoretical developments
2.1. Geometry studied
We consider a Newtonian fluid's steady uniform flow in the turbulent supercritical regime over a random packing of stationary spherical particles of diameter 
$d_p$. The fluid's dynamic viscosity is denoted by 
$\mu$, and its density is 
$\varrho$. The particle layer is of uniform depth and rests on an impervious, solid planar boundary inclined at angle 
$\theta$ (
$i=\tan \theta$ is the bed slope). The inclination angle is shallow, making it possible to assume that 
$i \sim \sin \theta$. The bed surface is initially flat at the macroscopic scale, but because of the random particle arrangement, the bed–stream interface is bumpy at the particle scale. Figure 1 shows a sketch of the problem studied. We use a Cartesian frame where the 
$x$-axis is aligned with the impervious wall, the 
$z$-axis is normal to that bed, and the 
$y$-axis points in the cross-stream direction. The flow's free surface is located at 
$z=z_{surf}$. We define a porosity function 
$\epsilon$ as the fraction of the pore volume over a control volume 
$V$ (specified in the next subsection). We assume that 
$\epsilon$ depends solely on 
$z$ because the bed is flat and of uniform thickness. We also assume that because the bed is stationary, the function 
$\epsilon (z)$ is prescribed.
Diagram of a gravity-driven turbulent flow over a rough permeable bed with shallow relative submergence. Here, 
$U_x(z)$ is the double-averaged velocity, 
$\epsilon (z)$ is the porosity, 
$z_{rc}$ is the roughness crest elevation, and 
$z_t$ is the elevation at which the bulk porosity 
$\epsilon _b$ is reached. The roughness layer is bounded by 
$z_{rc}$ and 
$z_t$, while the subsurface layer is below 
$z_t$. The surface layer is above 
$z_{rc}$, and in this region 
$\epsilon =1$.
We split the flow domain into three layers depending on 
$\epsilon$:
(i) The surface layer is the flow region above the roughness crest
$z_{rc}$, where $\epsilon =1$.(ii) The subsurface layer is assumed to be a homogeneous porous medium with a constant porosity
$\epsilon _b$.(iii) The roughness layer is a transition region, in which the porosity increases from
$\epsilon = \epsilon _b$ at $z=z_t$ to $\epsilon = 1$ at $z=z_{rc}$.
Nikora et al. (Reference Nikora, Goring, McEwan and Griffiths2001) proposed a slightly different partitioning, with the roughness layer split into interfacial and form-induced sublayers (between 
$z_t$ and 
$z_{rc}$, and above 
$z_{rc}$, respectively). Within the form-induced sublayer, the dispersive stresses gradually drop to zero. Because the measured dispersive shear stress in our experiments was vanishingly small for 
$z>z_{rc}$, our study did not consider sublayers. For the free stream, flow depth can be defined as 
$h=z_{surf}-z_{rc}$ as a first approximation, but subsequently, we refine this definition and instead use 
$h_f = z_{surf}-z_{\epsilon =0.8}$, where 
$z_{\epsilon =0.8}$ is the elevation for which 
$\epsilon =0.8$. Here, we are concerned with shallow flows – flows with shallow relative submergence, 
$h/d_p=O(1)$. The flow depth is, however, sufficiently large relative to 
$d_p$ for the free surface to remain flat and parallel to the bed.
Different flow regimes in the subsurface layer can be created depending on the flow velocity and pore size (Wood, He & Apte Reference Wood, He and Apte2020). In this paper, we focus on beds sufficiently coarse to be permeable and subject to momentum and mass exchanges with the free stream, but also sufficiently dense to quickly dampen fluid turbulence. The flow dynamics across the sediment–water interface can be described using the permeability Reynolds number 
${\textit {Re}}_K=\sqrt {K} u_* / \nu$, where 
$u_*$ denotes the shear velocity, 
$\nu =\mu /\varrho$ is the kinematic viscosity, and 
$K$ is the permeability (Voermans et al. Reference Voermans, Ghisalberti and Ivey2017). Here, we consider flows for which 
${\textit {Re}}_K=O(1)$.
We are interested in determining the mean velocity's streamwise component 
$U_x(z)$. This function is derived by averaging the local velocity field 
$\boldsymbol {u}=(u, v, w)$ using an appropriate averaging procedure (see next subsection). Fluid pressure 
$p$ is assumed to be hydrostatic, on average, throughout the flow domain.
2.2. Governing equations
A turbulent flow over and through a porous medium exhibits time and/or space fluctuations of the velocity and pressure fields. To derive the governing equations, we use the double-averaging technique, which can be seen as an extension of time-averaging for the Reynolds-averaged Navier–Stokes equations (Nikora et al. Reference Nikora, McEwan, McLean, Coleman, Pokrajac and Walters2007a). We first take the time average of the Navier–Stokes equations and use the Reynolds decomposition 
$\psi =\bar \psi +\psi '$, where 
$\psi$ denotes a flow variable (velocity or pressure), 
$\bar \psi$ is its temporal mean, and 
$\psi '$ is its fluctuation.
We then take the volume average of these Reynolds-averaged Navier–Stokes equations. Volume averaging is achieved by considering a mesoscopic control volume in the form of a thin parallelepiped: its length 
$L_*$ and width 
$W_*$ are much larger than the particle diameter 
$d_p$, whereas its depth 
$H_*$ is shallow relative to 
$d_p$. At any point 
$\boldsymbol {x}=(x, y, z)$, the integration volume is 
$V=(x-L_*/2, x+L_*/2)\times (y-W_*/2, y+W_*/2)\times (z-H_*/2, z+H_*/2)$. This volume 
$V$ comprises fixed solid (
$V_s$) and fluid (
$V_f$) volumes. We refer to 
$\boldsymbol {n}$ as the vector pointing from the solid particle to the fluid phase. We split the surface bounding the fluid volume 
$V_f$ into 
${\mathcal {A}}_s$ (the surface of the interface with the solid particles) and 
${\mathcal {A}}_f$ (the fluid surface of the volume boundaries, 
${\mathcal {A}}_f= {\mathcal {S}}_f$). See the electronic supplement for further information. For any arbitrary function 
$\psi$ related to the fluid phase, Gray (Reference Gray1975) defined the intrinsic phase average
\begin{equation} \langle \psi\rangle=\frac{1}{V_f}\int_{V_f}\psi \,\mathrm{d} V. \end{equation} We define the spatial fluctuation 
$\tilde \psi$ with respect to the double (time–space) average (Gray Reference Gray1975): 
$\bar {\psi }=\langle \bar {\psi }\rangle +\tilde \psi$. Using the decomposition based on the double-averaged variables, we end up with the following expression for the velocity field:
\begin{equation} \boldsymbol{u}(x,y,z,t) = \begin{bmatrix} u_x \\ u_y \\ u_z \\ \end{bmatrix} = \begin{bmatrix} \bar{u}_x +u_x^{\prime} \\ \bar{u}_y + u_y^{\prime} \\ \bar{u}_z+u_z^{\prime} \\ \end{bmatrix} = \begin{bmatrix} \langle\bar{u}_x \rangle + \tilde{u}_x +u_x^{\prime} \\ \langle\bar{u}_y \rangle + \tilde{u}_y + u_y^{\prime} \\ \langle\bar{u}_z \rangle + \tilde{u}_z+u_z^{\prime} \\ \end{bmatrix}. \end{equation}Based on developments by Whitaker (Reference Whitaker1996) and Nikora et al. (Reference Nikora, McEwan, McLean, Coleman, Pokrajac and Walters2007a), we can express the double-averaged Navier–Stokes equations in the form
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}(\epsilon\langle\boldsymbol{u}\rangle)=0, \end{gather}
\begin{gather} \frac{\partial \langle \boldsymbol{\bar u}\rangle}{\partial t}+ \epsilon^{{-}1}\,\boldsymbol{\nabla}\boldsymbol{\cdot}(\epsilon\langle\boldsymbol{\bar u}\rangle \langle\boldsymbol{\bar u}\rangle)= \boldsymbol{g}- \frac{1}{\varrho }\,\boldsymbol{\nabla} \langle \bar p\rangle+ \frac{\nu}{\epsilon} ({\nabla}^{2}(\epsilon \langle\boldsymbol{\bar u} \rangle) - \boldsymbol{\nabla} \langle\boldsymbol{\bar u}\rangle \boldsymbol{\cdot}\boldsymbol{\nabla}\epsilon )\nonumber\\ \hspace{4.5pc}+ \frac{1}{\epsilon\varrho}\,\boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{\tau}_t+ \boldsymbol{\tau}_d)+ \frac{\boldsymbol{\bar f}}{\varrho}, \end{gather}
where 
$\boldsymbol {g}$ denotes gravitational acceleration, 
$\boldsymbol {\tau }_t=-\epsilon \varrho \langle \overline {\boldsymbol { u}' \boldsymbol { u}'}\rangle$ is the Reynolds stress tensor, 
$\boldsymbol {\tau }_d=-\epsilon \varrho \langle \boldsymbol {\tilde u} \boldsymbol {\tilde u} \rangle$ is called the dispersive stress tensor, and 
$\boldsymbol {\bar f}$ is the drag force density resulting from fluctuating pressure and velocity gradients on the particle surface:
\begin{equation} \boldsymbol{\bar f}= \frac{1}{V_f}\int_{\mathcal{A}_s}(\mu\boldsymbol{\nabla} \boldsymbol{\tilde u}-\tilde p \boldsymbol{\mathsf{I}}) \boldsymbol{\cdot} \boldsymbol{n}\,\mathrm{d} S, \end{equation}
where 
$\boldsymbol{\mathsf{I}}$ is the identity tensor. For steady, uniform, one-dimensional open-channel flow over a porous bed, the 
$x$-component of the momentum equation reduces to
\begin{equation} 0 =\epsilon \varrho g i + \frac{\mathrm{d} \tau_{d}}{\mathrm{d} z} + \frac{\mathrm{d} \tau_{t}}{\mathrm{d} z} + \frac{\mathrm{d} \tau_v}{\mathrm{d} z}+\epsilon \bar f - \varrho \mu\,\frac{\mathrm{d} \epsilon}{ \mathrm{d} z}\,\frac{\mathrm{d} U_x}{\mathrm{d} z}, \end{equation}
where 
$U_x(z)=\langle \overline {u}_x \rangle$ is the streamwise component of the double-averaged velocity, 
$\bar f$ is the streamwise component of the drag force density 
$\boldsymbol {\bar f}$, 
$\tau _{t}=-\varrho \epsilon \langle \overline {u_x' u_z'} \rangle$ and 
$\tau _{d}=- \varrho \epsilon \langle \tilde {u}_x \tilde {u}_z \rangle$ are the turbulent and dispersive shear stresses, respectively, and 
${\tau _v}= \mu \,\mathrm {d}( \epsilon U_x)/ \mathrm {d} z$ is the viscous shear stress. The last term on the right-hand side of (2.6) is known as the second Brinkman correction (Ochoa-Tapia & Whitaker Reference Ochoa-Tapia and Whitaker1995). This term becomes zero in the surface and subsurface layers. The mesoscopic formulation's advantage over the macroscopic one is that the governing equation (2.4) holds everywhere. The individual contributions in (2.4) may vary significantly from one layer to another.
Equation (2.6) is subject to the usual boundary conditions for free surface flows down a sloping bed:
\begin{equation} U_x=0\ \text{at }z=z_b\quad \text{and}\quad \frac{\mathrm{d} U_x}{\mathrm{d} z}=0\ \text{at }z=z_{surf}. \end{equation}The first boundary condition in (2.7a,b) reflects the no-slip condition at the bottom wall, while the second assumes that the ambient air exerts no shear stress on the free surface.
2.3. Closure equation for the drag forces
The permeable bed is assumed to behave like a homogeneous porous medium for which the drag force density 
$\boldsymbol {\bar f}$ is approximated by a Darcy equation with a Forchheimer correction (Whitaker Reference Whitaker1996). Here we follow Breugem et al. (Reference Breugem, Boersma and Uittenbogaard2006) and use Ergun's (Reference Ergun1952) equation to close the drag force density 
$\boldsymbol {\bar f}$. Its streamwise component 
$\bar f$ can be expressed as
\begin{equation} \bar f={-}\frac{\mu
\epsilon}{ K(U_x)
}\,U_{x}={-}\overbrace{\underbrace{\frac{A_E
(1-\epsilon)^{2}}{ \epsilon d^{2}_p}\,\mu
U_{x}}_{\text{Kozeny\text{--}Carman}} -\underbrace{\frac{B_E
(1-\epsilon) }{ d_p }\, \varrho
U_{x}^{2}}_{\mathrm{Forchheimer~term}}}^{\mathrm{Ergun}},
\end{equation}
where 
$K(U_x)$ is a velocity-dependent effective permeability, and 
$A_E=180$ and 
$B_E=1.75$ are Ergun's empirical constants. In Ergun's nonlinear model, the force exerted by the fluid depends on the mean grain size 
$d_p$ and porosity 
$\epsilon$, and it is a quadratic function of the mean flow velocity 
$U_x$. When the quadratic term in Ergun's equation is negligible, we end up with Darcy's law, and in that case, the permeability 
$K(U_x)$ becomes a constant whose value is set by the Kozeny–Carman relation (
$K_{KC}=d_p^{2} \epsilon _b^{3}/(A_E(1-\epsilon _b)^{2}$), with 
$A_E=180$). When a steady state is reached in the subsurface layer, the stress gradients in the momentum balance equation (2.6) are small, and the drag force 
$f$ balances the driving force 
$\epsilon \varrho g i$. Solving (2.8) for 
$U_x$ provides the steady-state velocity 
$U_{x,SSL}$, whose leading-order estimate is
\begin{equation} U_{x,SSL}=\frac{ g i K_{KC}}{\nu}, \end{equation}
where 
$K_{KC}$ is the Kozeny–Carman permeability given by the first contribution on the right-hand side of (2.8). Note that this expression (2.8) is assumed to be valid everywhere in the fluid (and the same applies to the closure equations presented below).
2.4. Turbulent stress closure
To model the turbulent stress, we opt for Prandtl's mixing-length equation,
\begin{equation} \tau_t={-}\varrho \epsilon \langle \overline{u_x^{\prime} u_z^{\prime}} \rangle =\varrho \epsilon {\ell_{t}}^{2} \left( \frac{\mathrm{d} U_x}{\mathrm{d} z} \right)^{2}, \end{equation}
where, following Li & Sawamoto (Reference Li and Sawamoto1995), we assume that the mixing length 
$\ell _t$ can be parametrised using an integral approach. In their original formulation, Li & Sawamoto (Reference Li and Sawamoto1995) defined the mixing length as
\begin{equation} \ell_{t,LS}=\kappa\,Z_{LS}\quad \text{with } Z_{LS}=\int_{-\infty}^{z}\frac{\epsilon-\epsilon_b}{1-\epsilon_b}\,\mathrm{d} z, \end{equation}
where 
$\kappa$ denotes the von Kármán constant and 
$Z_{LS}$ is an integral depth, which avoids having to fix the position of the bed–stream interface. This equation has been used by Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013) and Maurin et al. (Reference Maurin, Chauchat, Chareyre and Frey2015), among others.
In open-channel flows, turbulence damping causes the mixing length to decrease in the buffer layer along a solid boundary (Nezu & Nakagawa Reference Nezu and Nakagawa1993; Pope Reference Pope2000; Dey Reference Dey2014). For smooth solid boundaries, Van Driest's equation is commonly used to model this damping on an empirical basis (Pope Reference Pope2000). According to Durán, Andreotti & Claudin (Reference Durán, Andreotti and Claudin2012), damping is also expected when the bed is made up of coarse particles. Following these authors, we combine the mixing length proposed by Li & Sawamoto (Reference Li and Sawamoto1995) with the one developed by Van Driest (Reference Van Driest1956):
\begin{equation} \ell_{t,\mu} = \kappa\,Z_{LS} \left(1- \exp \left(-\frac{\sqrt{Z_{LS}\,U_x(z)/ \nu}}{A^{*}}\right) \right), \end{equation}
where 
$A^{*} =26$ is an empirical constant calibrated by Van Driest. As will be shown in the simulation section, § 5, its application to the present context was mostly driven by its improvements to the model's performance.
2.5. Dispersive stress closure
Although there is a large body of work on dispersive stresses 
$\boldsymbol {\tau }_d$ based on experimental measurements (Mignot et al. Reference Mignot, Barthélemy and Hurther2009; Detert, Nikora & Jirka Reference Detert, Nikora and Jirka2010; Dey & Das Reference Dey and Das2012; Voermans et al. Reference Voermans, Ghisalberti and Ivey2017; Rousseau & Ancey Reference Rousseau and Ancey2020) or direct Navier–Stokes simulations (Fang et al. Reference Fang, Han, He and Dey2018; Kuwata & Kawaguchi Reference Kuwata and Kawaguchi2019; Shen et al. Reference Shen, Yuan and Phanikumar2020), we know of very few attempts to provide a scaling law or a closure equation for the dispersive stress tensor 
$\boldsymbol {\tau }_d$. Based on earlier investigations into dispersive stresses (Breugem et al. Reference Breugem, Boersma and Uittenbogaard2006; Rousseau Reference Rousseau2019), we suggest an empirical relationship for the dispersive shear stress closure:
\begin{equation} \tau_{d} ={-}\varrho \epsilon\langle\tilde{u}_x \tilde{u}_z\rangle =\varrho \epsilon \nu_d\,\frac{\mathrm{d} U_{x}}{\mathrm{d} z} =\varrho \epsilon \lambda d_p\,\frac{1-\epsilon(z)}{1-\epsilon_{b}}\,U_{x}\,\frac{\mathrm{d} U_{x}}{\mathrm{d} z}, \end{equation}which involves an effective dispersive viscosity
\begin{equation} \nu_d=\lambda d_p\,\frac{1-\epsilon(z)}{1-\epsilon_{b}}\,U_{x}, \end{equation}
where 
$\lambda$ is a proportionality constant that needs to be calibrated. We can also recast the dispersive shear stress in the form
\begin{equation} \tau_{d} = \varrho\epsilon\ell_d^{2}\,\frac{U_x}{ d_p}\,\frac{\mathrm{d} U_x}{\mathrm{d} z},\quad \text{where }\ell_d= d_p\sqrt{\lambda}\sqrt{ \frac{1-\epsilon(z)}{1-\epsilon_{b}} }, \end{equation}
where 
$\ell _d$ plays the role of the mixing length and is subsequently called the dispersive length.
3. Experimental protocol
3.1. Set-up
Experiments were performed in a 6 cm-wide tilting flume, as depicted in figure 2. A constant head tank provided a steady fluid discharge into the flume. Equal proportions of borosilicate beads of two diameters (7 and 9 mm for A runs, and 13 and 15 mm for B runs) were packed randomly onto the flume bottom, forming the coarse bed. Thus median diameters were 
$d_{p,A}=8$ mm and 
$d_{p,B}=14$ mm for runs A and B, respectively. A bimodal size distribution was used because otherwise, beads formed parallel layers causing undesirable biases in the averaged porosity and velocity profiles, as observed, for instance, by Ni & Capart (Reference Ni and Capart2015). Before each run, the upper layer was flattened out to form a uniform bed layer of height 
$h_s=5$ cm.
Sketch of the experimental set-up. See Movie 1 and Movie 2 in the online supplementary material available at https://doi.org/10.1017/jfm.2022.310 for an overview of the set-up and a video sample taken by the PIV-RIMS imaging system, respectively.
We used the refractive index matched (RIM) technique to visualise what happened within the stream and bed far from the sidewalls (Wiederseiner et al. Reference Wiederseiner, Andreini, Épely-Chauvin and Ancey2011). This technique involved matching the fluid's refractive index 
$n_f$ with that of the glass beads, making it possible to determine bead positions and probe interstitial flow velocities. The fluid was prepared by mixing volume concentrations of 60 % benzyl alcohol and 40 % ethanol (BAE). The advantage of combining borosilicate, ethanol and benzyl alcohol was that mixtures had bulk densities close to those found in rivers, and this allowed us to run experiments on sloping beds with no sediment transport. Because of the fixed volume of fluid in the reservoir, a steady state was maintained for approximately 40 s. Table 1 recaps the fluid and sediment characteristics.
Fluid and sediment characteristics measured at a controlled temperature 20 
$^{\circ }$C: slope 
$i$, fluid refractive index 
$n_f$, mean particle diameters 
$d_{p,A}$ and 
$d_{p,B}$ for runs A and B, respectively, bed thickness 
$h_s$, particle density 
$\varrho _s$, fluid density 
$\varrho$, and kinematic viscosity 
$\nu$.
3.2. Velocimetry and scanning
We employed a method combining particle image velocimetry (PIV) and refractive index matched scanning (RIMS). This method was presented in great detail in an earlier publication, to which we refer interested readers (Rousseau & Ancey Reference Rousseau and Ancey2020). We outline the method briefly below.
The PIV-RIMS method coupled continuous scanning and PIV; it was thus able to provide more information on turbulence than previous measurements based on fixed laser sheets. We followed a methodology inspired by Dijksman et al. (Reference Dijksman, Rietz, Lörincz and van Hecke2012), van der Vaart et al. (Reference van der Vaart, Gajjar, Epely-Chauvin, Andreini, Gray and Ancey2015) and Ni & Capart (Reference Ni and Capart2015): the porous bed and stream were scanned by moving a laser sheet across the stream's transverse axis. The laser sheet was displaced from 
$y=1$ cm to 
$y=4$ cm (the origin 
$y=0$ of the cross-stream axis is placed at the front sidewall). As shown by figure 3, the fluid volumes (
$y<1$ cm) adjacent to the sidewalls were influenced by the sidewall boundary layers and were thus not representative of the flow conditions far from the sidewalls – see figure 11 and related text in Rousseau & Ancey (Reference Rousseau and Ancey2020). For this reason, they were not scanned.
Three-dimensional view of the streamwise velocity component for run A1. (a) Side view of a flow slice with an indication of the velocity field at the transverse position 
$Y=25$ mm (where 
$Y$ is the distance from the sidewall). (b) Front view of the flow slice. This view enables us to appreciate how the sidewall influences the flow. Different bead colours (purple and pink) were used to distinguish between the two different diameters (
$d_p=7$ mm and 9 mm). The origins of 
$Z$ and 
$X$ are arbitrary in this plot. See Movie 3 in the online supplementary material, illustrating the velocimetry procedure applied to a raw video sample.
For PIV, we seeded the flow with micrometric tracers (made of hollow borosilicate glass spheres with diameters in the range of 8–12 
$\mathrm {\mu }$m) and lit them using the mobile laser sheet. The camera was placed 30 cm from the sidewall, filming an area 
$\Delta x\,\Delta z= 73.8 \times 34.5$ mm
$^{2}$ and operated at a rate of 420 frames per second. To estimate velocities from these images, we used a method called ‘feature tracking’ (Miozzi, Jacob & Olivieri Reference Miozzi, Jacob and Olivieri2008), which is included in the opyf Python package (available from the public GitHub repository https://github.com/groussea/opyflow).
First- and second-order turbulence statistics were computed by adjusting the scanning travel speed so that it was consistent with the constraints imposed by the laser sheet's thickness and space–time velocity fluctuations (see Rousseau & Ancey (Reference Rousseau and Ancey2020) for further explanations and validation of the PIV-RIMS method). We were able to produce the following volume-averaged quantities: the mean streamwise velocity component 
$U_x$, and the turbulent (
$\tau _t$) and dispersive (
$\tau _d$) shear stresses within a mesoscopic sampling volume whose dimensions were 
$\Delta x\,\Delta y\,\Delta z =73.8 \times 30\times 34.5$ mm
$^{3}$ (the mesoscopic scale 
$L_*=\Delta x$ over which spatial averaging was done was approximately ten bead diameters). The porosity (
$\epsilon$) profile was obtained by averaging over a thinner volume (1 px thick volume), so virtually 
$\epsilon$ was a plane-averaged quantity. The profiles were obtained by first averaging in time the velocity field at any point (voxel) within the control volume. We measured the time fluctuations and the spatial disturbance for each voxel. Finally, we deduced the resulting Reynolds and dispersive stresses, and averaged these quantities over the mesoscopic sampling volume 
$\Delta x\,\Delta y\,\Delta z$ parallel to the bed.
4. Results
Table 2 summarises the key features of the nine runs presented in this paper. Each run was made using the same flow discharge per unit width (
$q_f = 3 \times 10^{-3}$ m
$^{2}$ s
$^{-1}$). The bed slope 
$i$ varied from 0.5 % to 8 %. For the sake of comparison, we also report the key features of runs L12, L14 and L15 conducted by Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017), which we will compare with ours below. The main difference was that we investigated bed slopes as steep as 8 %, which led to faster surface and subsurface velocities than those observed by Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) in horizontal beds. In our experiments, the permeability Reynolds number 
${\textit {Re}}_K$ ranged from 2 to 9, whereas Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) investigated a 0.3–6 range. Figure 4 shows the typical velocity and shear stress profiles determined for runs A2 and A3, on which we comment in the following subsections.
Flow characteristics for each of the nine experimental runs (A and B runs) conducted and, for the sake of comparison, for three runs made by Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) (L12, L14 and L15). The variables are as follows: 
$d_p$ is the bimodal size distribution's mean particle diameter, 
$i$ is the bed slope, 
$q_f$ is the water discharge per unit width, 
$\epsilon _b$ is the mean bulk porosity, 
$h_f = z_{surf}-z_{\epsilon =0.8}$ is the flow depth, 
$u_*$ is the shear velocity, 
$u_{*,V}$ is the shear velocity following the Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) definition, as explained in § 4.2, 
$U_{b}$ is the surface layer's mean velocity, 
$U_{surf}$ is the free surface velocity, 
$U_{SSL}$ is the subsurface layer's mean velocity, 
$u_p=\sqrt {g d_p i}$ is a velocity scale, 
${\textit {Re}}_{b}= U_b h_f/ \nu$ is the bulk Reynolds number, 
$h_+= u_* h_f/\nu$ is the friction Reynolds number, 
${\textit {Re}}_{p}=U_{SSL} d_p/\nu$ is the particle Reynolds number, 
$K_{KC}$ is the permeability predicted by the Kozeny–Carman equation (2.8), 
${\textit {Re}}_K=\sqrt {K_{KC}}\,u_*/\nu$ is the permeability Reynolds number, 
${\textit {Re}}_*=u_{*} d_p/\nu$ is the roughness Reynolds number, and 
$Fr=U_{surf}/\sqrt {g h_{f}}$ is the Froude number.
Flow profiles for runs A2 (
$i=1\,\%$) and A3 (
$i=2\,\%$). We increased the bed slope from 1 % to 2 % while keeping a constant flow discharge. (a) Streamwise velocity profiles; the inset shows the subsurface velocities in a log–linear plot. The horizontal error bars represent uncertainties at the free surface level. (b) Porosity profile. (c,d) Shear stresses compared with the driving force density 
$G$ defined by (4.2). The total shear stress 
$\tau _{Tot}=\tau _v+\tau _t+\tau _d$ includes turbulent stress 
$\tau _t$, dispersive stress 
$\tau _d$, and viscous stress 
$\tau _v$. The bed-normal distance 
$z$ is scaled by the mean particle diameter 
$d_p$ and defined relative to the elevation 
$z_{\epsilon =0.8}$: 
$z^{\prime }=z-z_{\epsilon =0.8}$.
4.1. Porosity profiles and bed-normal origin definition
To compare the various runs’ velocity profiles, we needed to define a reference level that could be used to locate the bed–stream interface. We considered different possibilities. Some authors have suggested defining this reference level by fitting the differential form of the logarithmic velocity profile to the experimental profile (Cameron et al. Reference Cameron, Nikora and Stewart2017; Suga et al. Reference Suga, Okazaki, Ho and Kuwata2018; Shen et al. Reference Shen, Yuan and Phanikumar2020). We did not choose this option because logarithmic velocity profiles are not expected in shallow submergence conditions (Jiménez Reference Jiménez2004). One alternative was the definition given by Pokrajac et al. (Reference Pokrajac, Finnigan, Manes, McEwan and Nikora2006), who suggested that the roughness crest – defined as the elevation for which porosity is 0.99 (
$z_{rc} = z_{\epsilon =0.99}$) – could be considered the reference level. However, we found that this definition created significant scatter between the velocity profiles when the bed was rearranged. Indeed, 
$z_{\epsilon =0.99}$ was strongly influenced by individual grains that were slightly higher than the average bed level. We deduced that a better option was to set the reference level at a bed height where the scatter between porosity profiles was minimal. We found that the optimal value was 
$z_{ \epsilon =0.8}$ – see figure 13 in Rousseau & Ancey (Reference Rousseau and Ancey2020), which compared porosity profiles. In our experiments, 
$z_{ \epsilon =0.8}$ was located at approximately 
$0.3 d_p$ below the roughness crest 
$z_{rc}$. Using a bimodal size distribution for the bed packing, we observed no clear inflexion point at 
$0.3 d_p$ below the roughness crest, in contrast to what was observed for beds made up of equally sized beads (Voermans et al. Reference Voermans, Ghisalberti and Ivey2017; Shen et al. Reference Shen, Yuan and Phanikumar2020). In the following subsections, we will use a relative elevation:
\begin{equation} z^{\prime}=z-z_{ \epsilon=0.8}. \end{equation}4.2. Main flow characteristics
Figure 4 shows the bed slope's influence on velocity and shear-stress profiles for runs A2 and A3. As the typical flow depth did not vary much between runs (
$h_f\sim 10$ mm), changing the bed slope 
$i$ was the most effective way of varying the bottom shear stress (
$\tau _b\sim \varrho g h i$) from one run to another. The slope was increased from 
$1\,\%$ to 
$2\,\%$, and all other control parameters were held identical. Figures 4(a–d) show the velocity, porosity, and turbulent shear stress 
$\tau _t=-\varrho \epsilon \langle \overline {u_x^{\prime } u_z^{\prime } } \rangle$ for runs A2 and A3, respectively. If the flow is steady, one-dimensional and uniform, then the total shear stress counterbalances the driving gravitational force density:
\begin{equation} G(z)=\int_z^{z_{surf}} \epsilon (z) \varrho g i \,\mathrm{d} z = \tau_t +\tau_v+\tau_d. \end{equation} We computed 
$G$ by integrating 
$\epsilon \varrho g i$ between 
$z$ and 
$z_{surf}$, and plotted it in figures 4(c,d). These figures show that there was a good match between total shear stress 
$\tau _{Tot}=\tau _t +\tau _v+\tau _d$ and the driving stress 
$G$ in the surface layer, confirming that the flow was close to the steady state (the slight deviation in figure 4(d) was within the acceptable range of uncertainty). We also found that the turbulent shear stress 
$\tau _t$ followed closely the 
$G$ variations between 
$z_{surf}$ and 
$z^{\prime } \sim 0.3 d_p$, that is, at approximately the roughness crest 
$z^{\prime }_{rc}$. Within the roughness layer, the turbulent shear stress decreased significantly. We also observed that the dispersive stress and viscous stress (
$\tau _d =- \epsilon \varrho \langle \tilde {u}_x \tilde {u}_z \rangle$ and 
$\tau _v =\epsilon \mu U'_x(z)$, respectively) were one order of magnitude lower than the turbulent shear stress in the surface layer. The turbulent and viscous stresses reached their maximum near the roughness crest, which was accompanied by a change in convexity (reflected by an inflexion point) in the velocity profile. The dispersive shear stress reached its maximum at a lower elevation 
$z^{\prime } = 0$.
These observations in the surface and roughness layers were consistent with the Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) measurements of flows over rough horizontal beds when the permeability Reynolds numbers were close to ours. Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) estimated that the sediment–water interface was located at 
$z^{\prime } = 0$ (as explained in § 4.1). Although flow velocities were much lower in the Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) experiments than in ours, their velocity profiles, when scaled by the shear velocity 
$u_*$, were close to ours, as shown by figure 5(b). In our experiments on a sloping bed, the shear velocity was estimated as 
$u_*=\sqrt {g h_f i}$, whereas Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) inferred the shear velocity from the total stress maximum: 
$u_{*,V}=\sqrt { \tau _{Tot,max} / \varrho }$, where 
$\tau _{Tot,max}$ was the maximum of the 
$\tau _{Tot}$ profile shown in figure 4. Previous studies investigating flows adjacent to rough beds have provided varied definitions of shear velocity, and they are known to affect data interpretation, especially under shallow submergence conditions (Pokrajac et al. Reference Pokrajac, Finnigan, Manes, McEwan and Nikora2006).
Comparison of the streamwise velocity profiles from our A1, A2 and B2 runs with the L12, L14 and L15 runs obtained by Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) (
$d_p=25$ mm). (a) Streamwise velocity profiles in physical units. (b) Scaled streamwise velocities (velocities were scaled by the shear velocity 
$u_*$).
4.3. Subsurface flows
Figures 6(a,b) show the scaled streamwise velocity profile 
$U_x/u_{p}$ in a log–linear plot, where 
$u_p=\sqrt {g d_p i}$. Subsurface velocities were extracted directly from the profiles by defining the subsurface region between elevations 
$z^{\prime }\approx - 1.7 d_p$ and 
$z^{\prime } \approx - 0.7 d_p$ for 
$d_p=8$ mm, and between 
$z^{\prime }= - 0.4 d_p$ and 
$z^{\prime }= - d_p$ for 
$d_p=14$ mm (these positions are shown in figure 6). Subsurface velocity values were extracted from the A and B runs. We also considered the other runs to assess how sensitive the results were to bed arrangement (to that end, we repeated experiments by creating a new bed layer before each run). We refer the reader to figure 13 in Rousseau & Ancey (Reference Rousseau and Ancey2020): overall, we found that in the subsurface layer, particle rearrangement caused a relative variability of approximately 25 % in the velocity profiles.
Mean streamwise velocity profiles for our nine runs in a log–log plot. (a) Profiles for runs A1, A2, A3 and A4. (b) Profiles for runs B1, B2, B3, B4 and B5. The vertical line locates the bed-normal coordinate range over which the mean subsurface velocity 
$U_{SSL}$ was computed. The dotted horizontal line denotes the roughness crest 
$z'_{rc}$ in the two sets of experiments, while the solid horizontal line shows the reference elevation 
$z'_{\epsilon =0.8}=0$.
We measured the subsurface velocities and compared them with the velocities estimated by the Ergun and Kozeny–Carman equations (2.9). Figures 7(a,b) show how the mean subsurface velocity 
$U_{x,SSL}$ varies with bed slope 
$i$ for A runs and B runs, respectively. As a first approximation, we assumed the drag force density to be 
$f=\epsilon \varrho g i$, and using this assumption we could see the data in figure 7 as the dependence of the drag force density 
$f$ on subsurface velocity 
$U_{x, SSL}$. For 
${\textit {Re}}_{p}\leq 20$, bed slope and subsurface velocity seemed to be linked linearly, as predicted by the linear Kozeny–Carman relationship. Because of velocity fluctuations and uncertainties, it was difficult to be more assertive about the linear dependence of 
$f$ on 
$U_{x,SSL}$. At higher particle Reynolds numbers (for 
${\textit {Re}}_{p}>20$), nonlinear dependence was observed clearly, and it was consistent with Ergun's equation (2.8). This equation is known to perform well at describing flows in homogeneously packed beds (with 
$\epsilon \approx 0.4$).
(a) Measured subsurface velocities as a function of bed slope for 17 A-type runs. (b) Measured subsurface velocities as a function of bed slope for 11 B-type runs . We tested a variety of bed arrangements and slopes (
$i=0.5$ %–8 %). The data were compared with the linear Kozeny–Carman equation (dashed line) and the Ergun equation (dashed-dotted line), computed with the median grain diameter and the measured averaged bulk porosity.
4.4. Turbulence intensities
Figures 8(a–c) show the scaled turbulence intensities in the streamwise (
$\sigma _{u_x}=\langle \overline {u^{\prime 2}_x} \rangle ^{1/2}$) and bed-normal (
$\sigma _{u_z}=\langle \overline {u^{\prime 2}_z} \rangle ^{1/2}$) velocity components, and the scaled turbulent shear stresses for runs A1, A4, B1, B4, L12 and L15. Turbulence intensities in the 
$z$- and 
$x$-directions increased when going from the free surface to the sediment–water interface, and they reached a maximum between the roughness crest and 
$z^{\prime }=0$. We compared our measurements with the empirical equations provided by Nezu & Nakagawa (Reference Nezu and Nakagawa1993) that give the variation in 
$\sigma _{u_x}/u_*$ and 
$\sigma _{u_z}/u_*$ with the distance from smooth boundaries:
\begin{gather} \sigma_{u_x}/ u_* = 2.3 \exp (- (z-z_b)/ h_f ), \end{gather}
\begin{gather} \sigma_{u_z} / u_* = 1.27 \exp (- (z-z_b)/ h_f ), \end{gather}
where 
$z_b$ is the position of the smooth wall. These curves are plotted in figures 8(a,b), with 
$z-z_b$ given by 
$z^{\prime }$ (although this choice cannot be considered universal for rough beds).
Turbulence statistics from our A1, A4, B1 and B4 runs, and from the L12 and L15 runs obtained by Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017). (a) Double-averaged streamwise velocity fluctuations scaled by the shear velocity 
$u_*$. The dashed lines are the empirical curves for turbulence intensities given by Nezu & Nakagawa (Reference Nezu and Nakagawa1993) in (4.3) and (4.4). (b) Double-averaged normal velocity fluctuations scaled by 
$u_*$. (c) Scaled turbulent shear stress profile. The dotted lines represent the driving force density 
$G$ scaled by 
$\varrho \, u_*^{2}$. (d) Scaled streamwise fluctuating velocity 
$\langle \tilde {u}^{2}_x \rangle ^{1/2}$. (e) Scaled normal fluctuating velocity 
$\langle \tilde {u}^{2}_z \rangle ^{1/2}$. ( f) Scaled dispersive stress. The green dashed and dotted lines are the L12 and L15 run profiles measured by Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) and scaled by 
$u_{*,V}$.
Overall, the turbulence variations predicted by (4.3) and (4.4) captured the observed trends and provided the right order of magnitude for all the runs. Ghisalberti (Reference Ghisalberti2009) suggested that the following scaling held at the roughness crest for obstructed shear flows: 
$\sigma _{u_x}(z=z_{rc}) / u_* \sim 1.8$ and 
$\sigma _{u_z}(z=z_{rc}) / u_* \sim 1.1$. Here, we found that the measured streamwise turbulence intensities 
$\sigma _{u_x}(z=z_{rc}) / u_*$ ranged from 1.5 to 1.7, and thus were slightly lower than 1.8. The values of the bed-normal intensities 
$(\sigma _{u_z}(z=z_{rc})/ u_*$ were approximately 0.7. Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) also found that 
$\sigma _{u_z}(z=z_{\epsilon =0.8})/ u_*$ ranged from 0.5 to 0.7, which was consistent with our observations at 
$z_{\epsilon =0.8}$.
Figure 8(c) shows that turbulent stress maxima ranged from 0.5 to 1, in agreement with earlier studies (Nikora et al. Reference Nikora, Goring, McEwan and Griffiths2001; Mignot et al. Reference Mignot, Barthélemy and Hurther2009; Voermans et al. Reference Voermans, Ghisalberti and Ivey2017). We observed that the scaled turbulent shear stress 
$\langle \overline {u_x^{\prime } u_z^{\prime }} \rangle / u_*^{2}$ decreased with increasing bed slope 
$i$ and Reynolds number 
${\textit {Re}}_K$. This behaviour contrasted with that observed by Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017), who reported an increase in turbulent shear stress with 
${\textit {Re}}_K$. As the scaled shear stress depended on how 
$u_*$ was defined, this decrease might have been related to the shear velocity's definition . We observed that scaling the turbulent shear stress by 
$u_{*,V}$, instead of by 
$u_*$, modified the peak value of the turbulent shear stress, but regardless of the scale chosen, the peak value decreased as 
${\textit {Re}}_K$ increased (see Rousseau (Reference Rousseau2019) for a discussion about the definition of shear velocity). Another feature of shallow submergence flows was noticeable: for runs A4 and B4, the peak values of the turbulent shear stress were localised below 
$z_{rc}$, confirming previous observations that the permeable bed was more subject to turbulent shear stress when 
${\textit {Re}}_K$ was increased (Voermans et al. Reference Voermans, Ghisalberti and Ivey2017).
4.5. Dispersive stresses
Figures 8(d–f) show the scaled streamwise and bed-normal fluctuating components (
$\langle \tilde {u}^{2}_x \rangle ^{1/2}/u_*$ and 
$\langle \tilde {u}^{2}_z \rangle ^{1/2} /u_*$), and the scaled dispersive shear stress 
$\langle \tilde {u}_x \tilde {u}_z \rangle /u_*^{2}$, respectively, for our A and B runs, and the L12 and L15 runs from Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017). The scaled profiles had maxima at elevations below the roughness crest. This was consistent with observations by Nikora et al. (Reference Nikora, Goring, McEwan and Griffiths2001). We found that 
$\langle \tilde {u}^{2}_x \rangle ^{1/2}$ and 
$\langle \tilde {u}^{2}_z \rangle ^{1/2}$ reached their maxima at two different positions. The distance between these positions was approximately 
$0.5 d_p$. The Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) experiments led to similar results, but when comparing their results with ours, we noted differences: for L12, the peaks of 
$\langle \tilde {u}^{2}_x \rangle ^{1/2}$ and the dispersive stress 
$\langle \tilde {u}_x \tilde {u}_z \rangle /u_*^{2}$ occurred at 
$z_{rc}$. For runs A and B, the dispersive stresses were vanishingly small at 
$z_{rc}$, whereas the peak values of 
$\langle \tilde {u}_x \tilde {u}_z \rangle$ were located at 
$z^{\prime } = 0$ (that is, where the porosity was equal to 0.8). These fluctuating velocity intensities and dispersive stresses were similar to those shown in the numerical simulation performed by Fang et al. (Reference Fang, Han, He and Dey2018), who found that the peaks of 
$\langle \tilde {u}^{2}_x \rangle ^{1/2}$ and 
$\langle \tilde {u}^{2}_z \rangle ^{1/2}$ were positioned below the roughness crest.
4.6. Dispersive stress closure
For 
${\textit {Re}}_K>2$, the dispersive shear stress was concentrated within a thin layer below the roughness crest. The closure equation (2.13) was used and compared with the dispersive stress data. This equation involved the velocity 
$U_{x}$ and porosity 
$\epsilon$ profiles, as well as a 
$\lambda$ parameter that needed to be calibrated. When using 
$\lambda =0.015$, we found a good match between computed and measured dispersive shear stresses, as shown in figure 9 for all but one run: for run A4, the closure equation (2.13) underestimated the maximum dispersive stress by a factor of 2, whereas for other runs, the relative error did not exceed 25 %.
Comparison between the measured and computed dispersive shear stresses. The solid lines show the experimental data, while the dotted lines show the dispersive shear stress computed using the closure equation (2.13). (a) Runs A1 to A4 (
$d_p=8$ mm). (b) Runs B1 to B4 (
$d_p=14$ mm).
Voermans et al. (Reference Voermans, Ghisalberti and Ivey2018) suggested that the effective dispersive viscosity was related to the effective transverse diffusive coefficient 
$D_T$ through the dispersive Schmidt number, which is the ratio of the dispersive diffusion to the effective dispersive viscosity (
$Sc_d= \nu _d/D_T$). In the literature devoted to dispersion in porous media (Sahimi Reference Sahimi2011), the transverse dispersion is often found to vary linearly with 
$U_x$:
\begin{equation} D_T=\alpha_T\,\frac{d_p}{2}\,U_x. \end{equation}
Within our computational framework, (2.14) defines the effective dispersive viscosity, which takes the following form in the subsurface layer (where 
$\epsilon =\epsilon _b$):
\begin{equation} \nu_{d}= \lambda d_p U_{x}. \end{equation}
Equations (4.5) and (4.6) are structurally similar. Furthermore, we found that the fitted value 
$\lambda =0.015$ was consistent with the 
$\alpha _T$ values found to be in the 0.01–0.05 range in the related literature (Sahimi Reference Sahimi2011, p. 360). This finding supports the Voermans et al. (Reference Voermans, Ghisalberti and Ivey2018) statement that the dispersive Schmidt number is close to unity in the subsurface layer.
4.7. Mixing length
The mixing length can be estimated experimentally by inverting Prandtl's equation:
\begin{equation} \ell_{t}= \frac{\sqrt{-\langle \overline{u_x^{\prime} u_z^{\prime}} \rangle}}{\mathrm{d} U_x/\mathrm{d} z}. \end{equation} Figure 10 shows that the scaled mixing lengths 
$\ell _t/h_f$ for the surface layer (
$z>z_{rc}$) fell onto the same curve for runs A and B. The empirical distribution was similar to the one determined by Nezu & Rodi (Reference Nezu and Rodi1986). For 
$z>z^{\prime }=0.7 h_f$, the mixing length decreased when approaching the free surface. For runs A1, A2, A3, B1 and B2, the mixing length tended toward zero at the free surface – a behaviour that was consistent with the velocity defect effect (Coles Reference Coles1956; Nezu & Rodi Reference Nezu and Rodi1986). In the vicinity of the roughness crest 
$z_{rc}$, the mixing length increased nonlinearly, which evoked the damping effect observed by Van Driest (Reference Van Driest1956) in turbulent shear flows near a smooth solid boundary. All the profiles exhibited a local minimum near the roughness crest. This local minimum, as well as the overall behaviour of the mixing-length profiles, was consistent with the observations of Ghisalberti & Nepf (Reference Ghisalberti and Nepf2004).
Scaled mixing length evaluated using (4.7). (a) Runs A1 to A4 (
$d_p=8$ mm). (b) Runs B1 to B5 (
$d_p=14$ mm). We have also plotted the mixing-length profile predicted by Coles (Reference Coles1956), with 
$\varPi =0$ and 
$\varPi =0.1$, together with (4.8) proposed by Nezu & Rodi (Reference Nezu and Rodi1986), which includes both turbulence damping and Coles’ law of the wake.
Figure 10 shows two empirical mixing-length curves: the equations of Coles (Reference Coles1956) and Van Driest (Reference Van Driest1956) account for a velocity defect (also called Coles’ law of the wake) and turbulence damping, respectively. Combining these two effects, Nezu & Rodi (Reference Nezu and Rodi1986) suggested expressing the mixing length 
$\ell _t$ as
\begin{equation} \ell_{t}=h_f\, \kappa \sqrt{1-\frac{z^{\prime}}{h_f} } \left[ \frac{h_f}{z^{\prime}} + {\rm \pi}\, \varPi\, \sin \left({\rm \pi}\, \frac{z^{\prime}}{h_f} \right) \right] \varGamma, \end{equation}where
\begin{equation} \varGamma=1- \exp \left( \frac{-u_* z^{\prime} }{\nu \, A^{*}} \right) \end{equation}
is van Driest's damping function, and 
$\varPi$ is the Coles parameter expressing the strength of the wake function – see Nezu & Rodi (Reference Nezu and Rodi1986) or Pope (Reference Pope2000, pp. 305–308) for further information. Here, the distance from the wall is given by 
$z^{\prime }$. Accounting for the turbulence damping effect and the wake effect in the 
$\ell _{t}$ definition provided a better agreement than the original expression of the mixing length based on a constant von Kármán coefficient.
For runs B1 to B5, the agreement between empirical curves and experimental data was less pronounced. We observed that the empirical curves were higher, thereby suggesting that the reference 
$z_{\epsilon =0.8}$ was not the optimal value for B runs. The failure to find the same reference level for the sediment–water interface across all the runs led us to think that models relying on a vertical origin were not suitable for shallow submergence flows. This belief led us to develop closure equations based on an integral length presented in § 2.4 and used, among others, by Li & Sawamoto (Reference Li and Sawamoto1995), Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013) and Maurin et al. (Reference Maurin, Chauchat, Chareyre and Frey2015).
Figure 11 shows the scaled mixing-length profile. By scaling the mixing length 
$\ell _t$ by 
$d_p \sqrt {(1-\epsilon )/(1-\epsilon _b)}$, we were able to collapse all experimental 
$\ell _t$ values onto the vertical straight line 
$\lambda '= \ell _{t}/(d_p \sqrt {(1-\epsilon )/(1-\epsilon _b)}) \sim 0.1$ in the subsurface layer (
$z'<0$). This showed that the mixing length was very close to the dispersive length defined in (2.15); because 
$\lambda '\sim \sqrt {\lambda }$, we have 
$\ell _t\sim \ell _d$ in the subsurface layer, whereas 
$\ell _t>\ell _d$ in the roughness layer. This match-up shows that the dispersive and turbulent stresses had the same characteristic length, which was a small fraction of the particle diameter and thus related to the typical pore size.
Mixing length profiles evaluated from (4.7) in the roughness layer and scaled by the dispersive length 
$\ell _d$. (a) Runs A1 to A4 (
$d_p=8$ mm). (b) Runs B1 to B5 (
$d_p=14$ mm).
5. Simulations and discussion
5.1. Numerical scenarios
We solved the double-averaged momentum equation (2.6) supplemented by the closure equations presented in §§ 2.3–2.5. Drag forces in the porous bed were modelled using Ergun's equation (2.8). For Prandtl's mixing-length equation (2.10), we used each of the three algebraic expressions for 
$\ell _t$: 
$\ell _{t,LS}$ given by (2.11), 
$\ell _{t,\mu }$ given by (2.12), and 
$\ell _{t,D}$ given by (5.1). We considered three computational scenarios, as follows.
(i) The L&S scenario was based on the Li & Sawamoto (Reference Li and Sawamoto1995) mixing-length closure equation for
$\ell _{t,LS}$, (2.11).(ii) The damp. + L&S scenario used the
$\ell _{t,\mu }$ mixing length closure given by (2.12), which includes a damping correction.(iii) The disp. + damp. + L&S scenario combined turbulent and dispersive stresses. Since the mixing and dispersive lengths were similar in the subsurface layer, and because the dispersive length dropped to zero in the surface layer, we used a generalised mixing length:
(5.1)where\begin{equation} \ell_{t,D} = \sqrt{\ell_{t,\mu}^{2} + \ell_d^{2} }, \end{equation}$\ell _d$ is the dispersive length given by (2.15), and $\ell _{t,\mu }$ is the mixing length given by (2.12), including a damping correction.
None of these scenarios included the velocity defect effect. Indeed, although we found (see § 4.7) that it explained the decrease in the mixing length 
$\ell _t$ near the free surface, consistent with earlier measurements (Nezu & Rodi Reference Nezu and Rodi1986; Pope Reference Pope2000), its influence on the velocity profile was negligible because as 
$\ell _t$ decreased to zero, the shear stress and shear rate also dropped to zero.
The system of equations was solved using a finite-difference scheme. We considered the hydraulic conditions reported in table 2. One parameter was adjusted from experiments: 
$\lambda =0.015$. These values held for slopes ranging from 1 % to 8 %. Runs A2 and B5 were compared to the simulated profiles in figures 12 and 13, respectively. We also applied the model to runs L12, L14 and L15 conducted by Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017); as these authors did not provide a porosity profile, we synthesised it based on their figure 6.
Simulated and measured double-averaged profiles for run A2. Vertical profiles of (a) mixing length, (b) turbulent 
$\tau _t$, (c) dispersive 
$\tau _d$, (d) drag force density 
$f$, and (e) streamwise velocity. Panel ( f) shows the different contributions to the momentum balance equation (2.6). Three closure scenarios are compared: (L&S) for the Li & Sawamoto (Reference Li and Sawamoto1995) mixing length, (damp. + L&S) accounting for turbulence damping, and (disp. + damp. + L&S) including dispersion.
Simulated and measured double-averaged profiles for run B5. Vertical profiles of (a) mixing length, (b) turbulent 
$\tau _t$, (c) dispersive 
$\tau _d$, (d) drag force density 
$f$, and (e) streamwise velocity. Panel ( f) shows the different contributions to the momentum balance equation (2.6). Three closure scenarios are compared: (L&S) for the Li & Sawamoto (Reference Li and Sawamoto1995) mixing length given, (damp. + L&S) accounting for turbulence damping, and (disp.+ damp. + L&S) including dispersion.
5.2. Observations on experimental and simulated profiles
We tested the three scenarios and compared their numerical predictions with experimental data in figures 12 and 13. The L&S scenario overestimated systematically the mixing length 
$\ell _t$ in the roughness layer. As a consequence, the streamwise velocity was underestimated.
In the damp. + L&S scenario, which accounted for turbulence damping near the bed surface, the agreement was better. However, the increase in the mixing length within the roughness layer was not predicted. When including dispersive mechanisms – the disp. + damp. + L&S scenario – we obtained a better agreement between numerical solutions and experimental data for the turbulent mixing length. We also observed a substantial improvement in the velocity profile's shape for the roughness layer. In the surface layer, none of the three scenarios reproduced the decrease in the mixing length near the free surface. This suggested that including the velocity defect law could improve the predictive capacity. Although (5.1) led to a better agreement between computed and experimental profiles, it could not capture all the patterns exhibited by the mixing-length profiles.
This comparison highlighted how turbulence damping and dispersion could produce complementary effects. Indeed, the damping effect, which increased with decreasing Reynolds number, tended to reduce turbulent vertical exchanges. Dispersion was an additional contributor to vertical exchange in momentum. Then damping effects might become negligible as the Reynolds number increases, but an opposite trend was expected for the dispersive shear stress.
Figures 12( f) and 13( f) show the scaled contributions to the momentum balance equation (2.6) predicted by the disp.+ damp. + L&S scenario. As expected, this scenario predicts the prevalence of the drag force density in the subsurface layer and the predominance of the turbulent shear stress gradient in the surface layer. There was a sharp decay in the drag force density and the turbulent shear-stress gradient in the roughness layer. Furthermore, for a shallow region of thickness 
$\sim 0.2d_p$ at the bed interface, the dominant balance was between the drag force and the turbulent and dispersive shear stresses. The second Brinkman correction term 
$-\varrho \nu \,(\mathrm {d}\epsilon /\mathrm {d} z)\, (\mathrm {d} U _x/\mathrm {d} z)$ was vanishingly small across the entire layer, an observation that supports a posteriori the working assumption made by Nikora et al. (Reference Nikora, Koll, McEwan, McLean and Dittrich2004) and Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) that this term plays a negligible role in the momentum balance equation in the roughness layer of flows in a transitional regime (
${\textit {Re}}_K=O(1$)).
5.3. Flows over rough horizontal beds
To assess how robust the computational framework presented in § 2 was, we applied the model to runs conducted by Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) over horizontal beds. Figure 14 compares the model predictions with their experimental data for run L14 (see the Zenodo repository indicated in the Acknowledgements for the other runs). Since the bed was horizontal in the Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) experiments, the driving force 
$f_d$ was inferred from the information that they provided: 
$f_d=\tau _{max}/(\delta -z_{U})$, where 
$\delta$ was the estimated boundary layer thickness and 
$z_{U} \sim 0.4 d_p$ was related to the bed-normal coordinate at which 
$\tau _{Tot}=\tau _{max}=\varrho u_{*,V}^{2}$ (see Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) for more details). We used a synthetic porosity profile based on figure 16 in Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017), where the roughness layer thickness typically scaled with the bed grain size. We observed good overall agreement between the model output and experimental data. Including both the damping and dispersive mechanisms increased the agreement between the theory and the experiment. This might come as no surprise since, like us, Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017) worked at intermediate Reynolds numbers. This comparison provided further evidence that both damping and dispersive shear stress were important when predicting flows over rough beds in a transitional regime (
${\textit {Re}}_K=O(1$)).
Simulated and measured double-averaged profiles for run L12 conducted by Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017). Vertical profiles of (a) mixing length, (b) turbulent 
$\tau _t$, (c) dispersive 
$\tau _d$, (d) drag force density 
$f$, and (e) streamwise velocity. Panel ( f) shows the different contributions to the momentum balance equation (2.6). Three closure scenarios are compared: (L&S) for the Li & Sawamoto (Reference Li and Sawamoto1995) mixing-length model, (damp. + L&S) including turbulence damping using an empirical damping-correcting factor, and (disp.+ damp. + L&S) including turbulence damping and dispersive shear stress.
6. Conclusion
This paper examined the flow dynamics of shallow supercritical turbulent flows over sloping permeable beds made up of coarse grains – a situation typical of steep streams in mountainous terrain but which may also be in line with flows over canopies. To gain new physical insights into this problem, we used innovative imaging techniques based on refractive index matched scanning. This enabled us to reconstruct the velocity field throughout the bed and stream, far from the sidewalls. Experiments were conducted at the same flow discharge rate. We varied particle sizes (from 8 mm to 14 mm) – but kept the grain-size distribution narrow – and bed inclination (from 0.5 % to 8 %).
To move a step further in our understanding of steep streams, we used the double-averaged Navier–Stokes equations (2.4) supplemented by two closure equations (2.10) and (2.13) for the turbulent and dispersive shear stresses. In both cases, the closure was based on mixing lengths. For the turbulent shear stress, we used the variant of Prandtl's equation proposed by Li & Sawamoto (Reference Li and Sawamoto1995), and we also accounted for turbulence damping near the bed–flow interface by taking inspiration from the Van Driest (Reference Van Driest1956) equation. Coles’ law of the wake was considered a second-order effect, and was thus not included in the model. We knew of no algebraic closure for the dispersive shear stress, and we thus proposed an empirical equation that was structurally close to the definition of turbulent shear stress and consistent with earlier investigations. For the drag force exerted by the fluid on the bed particles, we used Ergun's equation (2.8). In total, our model involved a single free parameter 
$\lambda$, which was fitted to our data (
$\lambda =0.015$).
Overall, the model was able to capture the flow features in our experiments, namely the mean flow velocity in the bed, the velocity profile across the stream, and the turbulent and dispersive shear stresses. Furthermore, our model interpreted the non-logarithmic velocity profile as the consequence of two additive processes: on the one hand, there was noticeable turbulence damping near the bed interface (which we were able to account for by using a damping correcting function in the mixing-length equation (2.11)), and on the other hand, the dynamics of the roughness layer was controlled by the drag force and the dispersive and turbulent shear stresses.
From this perspective, the reason why the Manning–Strickler equation fails to predict flow resistance in steep streams is related directly to the existence of the roughness layer, which is associated closely with the production of dispersive stresses and drag forces.
Supplementary material and movies
Supplementary material and movies are available at https://doi.org/10.1017/jfm.2022.310.
Acknowledgements
We would like to thank P. Frey and I. Pascal for their feedback and all the discussions we had during Rousseau's thesis. We are grateful to J.J. Voermans, who graciously agreed to share his data published in Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017). We are grateful to the anonymous reviewers for the numerous astute remarks that helped to improve the paper.
Funding
This research was supported by EPFL. All the data and Python scripts used for plotting figures are archived at Zenodo (https://doi.org/10.5281/zenodo.5801325).
Declaration of interests
The authors report no conflict of interest.
Open access
