2.1 Theoretical background

The derivation of a general relationship for the friction factor for fixed- and mobile-rough-bed flows starts with the double-averaged (in time and in volume of a thin slab parallel to the bed) momentum equation proposed in Nikora et al. (Reference Nikora, Ballio, Coleman and Pokrajac2013). Here, for simplicity (but without losing conceptual generality) we consider a simpler case of fixed-rough-bed OCF for which the double-averaged momentum equation for the streamwise velocity component can be written as

(2.1) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\langle \overline{u}_{1}\rangle }{\unicode[STIX]{x2202}t}+\langle \overline{u}_{j}\rangle \frac{\unicode[STIX]{x2202}\langle \overline{u}_{1}\rangle }{\unicode[STIX]{x2202}x_{j}} & = & \displaystyle gS_{b}-\frac{1}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}\langle \overline{p}\rangle }{\unicode[STIX]{x2202}x_{1}}-\frac{1}{\unicode[STIX]{x1D719}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}\langle \overline{u_{1}^{\prime }u_{j}^{\prime }}\rangle }{\unicode[STIX]{x2202}x_{j}}-\frac{1}{\unicode[STIX]{x1D719}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}\langle \widetilde{\overline{u}}_{1}\widetilde{\overline{u}}_{j}\rangle }{\unicode[STIX]{x2202}x_{j}}+\frac{1}{\unicode[STIX]{x1D719}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\unicode[STIX]{x1D719}\left\langle \overline{\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}u_{1}}{\unicode[STIX]{x2202}x_{j}}}\right\rangle \nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{V_{f}}\iint _{S_{int}}\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}\overline{u}_{1}}{\unicode[STIX]{x2202}x_{j}}m_{j}\text{d}S+\frac{1}{\unicode[STIX]{x1D70C}}\frac{1}{V_{f}}\iint _{S_{int}}\overline{p}m_{1}\text{d}S.\end{eqnarray}$$

In (2.1), the Reynolds decomposition $\unicode[STIX]{x1D703}=\overline{\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D703}^{\prime }$ for instantaneous variables and the decomposition $\overline{\unicode[STIX]{x1D703}}=\langle \overline{\unicode[STIX]{x1D703}}\rangle +\widetilde{\overline{\unicode[STIX]{x1D703}}}$ for time-averaged variables are used where $\unicode[STIX]{x1D703}$ is a flow variable (i.e. velocity or pressure) defined in the fluid but not in the space occupied by the roughness elements; angle brackets denote spatial (volume) averaging, overbar denotes time averaging and prime and tilde denote temporally and spatially fluctuating components, respectively; $u_{i}$ is the $i$ th component of the velocity vector; $p$ is pressure; $m_{i}$ is a unit vector normal to the bed surface and directed into the fluid; $S_{int}$ is the extent of water–bed interface bounded by the averaging domain; $\unicode[STIX]{x1D719}=V_{f}/V_{0}$ is the roughness geometry function or ‘spatial porosity’ (Nikora et al. Reference Nikora, Goring, McEwan and Griffiths2001), where $V_{f}$ is volume occupied by fluid within the averaging domain of the total volume $V_{0}$ ; and $x_{1}=x$ (streamwise), $x_{2}=y$ (spanwise), $x_{3}=z$ (bed normal) are orthogonal coordinates. In wall-bounded flows such as rivers, there are strong gradients in flow properties in the wall-normal direction, especially near the bed, and therefore the volume-averaging domain should be a thin slab parallel to the mean bed. The dimensions of the averaging domain in the plane parallel to the mean bed define a scale of consideration and should typically exceed relevant correlation length scales of bed roughness. Thus, the double-averaged quantities in (2.1) depend, in general, on all three spatial coordinates ( $x$ , $z$ and $y$ ) and time $t$ , considering length and time scales well exceeding characteristic scales of bed roughness and turbulence, respectively (Nikora et al. Reference Nikora, McEwan, McLean, Coleman, Pokrajac and Walters2007). For simplicity, however, we will omit the explicit indication ( $x$ , $y$ , $z$ , $t$ ) of this dependency, presenting double-averaged quantities $\langle \overline{u}_{i}\rangle (x,y,z,t)$ , $\langle \overline{p}\rangle (x,y,z,t)$ , $\langle \overline{u_{i}^{\prime }u_{j}^{\prime }}\rangle (x,y,z,t)$ and $\langle \widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{j}\rangle (x,y,z,t)$ as $\langle \overline{u}_{i}\rangle$ , $\langle \overline{p}\rangle$ , $\langle \overline{u_{i}^{\prime }u_{j}^{\prime }}\rangle$ and $\langle \widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{j}\rangle$ . Compared to the conventional ensemble-(time-) averaged momentum equation, (2.1) contains two groups of additional terms: (i) dispersive stresses $\langle \widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{j}\rangle$ due to potential correlations of spatial variations in time-averaged velocity fields (also known as form-induced stresses); and (ii) viscous drag per unit fluid volume $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D708}}=(1/V_{f})\iint _{S_{int}}(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}\unicode[STIX]{x2202}\overline{u}_{1}/\unicode[STIX]{x2202}x_{j})m_{j}\text{d}S$ and form (pressure) drag per unit fluid volume $\unicode[STIX]{x1D6F7}_{p}=-(1/V_{f})\iint _{S_{int}}\overline{p}m_{1}\text{d}S$ , which together define the total drag per unit volume $F_{D}=\unicode[STIX]{x1D719}(\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D708}}+\unicode[STIX]{x1D6F7}_{p})$ (e.g. Nikora et al. Reference Nikora, McEwan, McLean, Coleman, Pokrajac and Walters2007, Reference Nikora, Ballio, Coleman and Pokrajac2013). The development of the double-averaged hydrodynamic equations for rough-bed flows was initiated by atmospheric scientists for describing turbulent flows within and above terrestrial canopies such as forests or bushes (Wilson & Shaw Reference Wilson and Shaw1977; Raupach & Shaw Reference Raupach and Shaw1982; Finnigan Reference Finnigan2000) and later this approach was adopted in studies of rough-bed open-channel flows.

In a first step of derivation, (2.1) is appropriately re-arranged and integrated from the roughness troughs $z_{t}$ to the water surface $z_{ws}$ , to obtain an equation that explicitly involves the local bed shear stress $\unicode[STIX]{x1D70F}_{0}$ , defined here as

(2.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}_{0}=\int _{z_{t}}^{z_{c}}\frac{1}{V_{0}}\left[\iint _{S_{int}}\left(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}\overline{u}_{1}}{\unicode[STIX]{x2202}x_{j}}\right)m_{j}\,\text{d}S-\iint _{S_{int}}\overline{p}m_{1}\,\text{d}S\right]\text{d}z=\int _{z_{t}}^{z_{c}}F_{D}\,\text{d}z, & & \displaystyle\end{eqnarray}$$

where $z_{c}$ is the elevation of the roughness crests. Then, a triple integration $\int _{z_{t}}^{z_{ws}}\text{d}z\int _{z_{t}}^{z}\text{d}z\int _{z_{t}}^{z}\text{d}z$ is applied to the difference between (2.1) and the obtained equation, similar to Fukagata et al.’s (Reference Fukagata, Iwamoto and Kasagi2002) procedure for smooth-bed flows. The first two integrations are required to obtain an explicit equation for the local double-averaged velocity $\langle \overline{u}_{1}\rangle$ while the third integration over the flow depth gives a depth-averaged velocity $\widehat{U}=H^{-1}\int _{z_{t}}^{z_{ws}}\unicode[STIX]{x1D719}\langle \overline{u}_{1}\rangle \,\text{d}z$ , where $H$ is the mean flow depth, defined as $H=\int _{z_{t}}^{z_{ws}}\unicode[STIX]{x1D719}\,\text{d}z=z_{ws}-z_{c}+\int _{z_{t}}^{z_{c}}\unicode[STIX]{x1D719}\,\text{d}z$ . As follows from this definition, the depth-averaged velocity $\widehat{U}$ and ‘local’ friction factor $f=8\unicode[STIX]{x1D70F}_{0}/(\unicode[STIX]{x1D70C}\widehat{U}^{2})$ may depend, in general, on the streamwise and transverse coordinates and time (at time scales well exceeding turbulent scales and at length scales well exceeding roughness scales). For steady, uniform, and two-dimensional flow (in the double-averaged sense) these quantities do not depend on spatial coordinates and time.

Subsequent re-arrangement of the integrated equation leads to the following decomposition of the friction factor

(2.3)

where the Reynolds number is defined as $Re=\widehat{U}H/\unicode[STIX]{x1D708}$ ; $Q=\widehat{U}H$ is a specific flow rate (i.e. per unit width); $N=3(L_{\unicode[STIX]{x1D70F}}/H)^{2}+(\unicode[STIX]{x1D70F}_{02D}/\unicode[STIX]{x1D70F}_{0})(L_{\unicode[STIX]{x1D719}}/H)^{3}-(H_{m}/H)^{3}$ is a parameter characterising flow–rough-bed interactions; $\unicode[STIX]{x1D70F}_{02D}=\unicode[STIX]{x1D70C}gS_{b}H$ is the bed shear stress that would act if the flow was steady, two-dimensional (2-D) and uniform in terms of the double-averaged quantities (Nikora et al. Reference Nikora, McEwan, McLean, Coleman, Pokrajac and Walters2007); $H_{m}$ is the maximum flow depth, i.e.  $H_{m}=z_{ws}-z_{t}$ ; $L_{\unicode[STIX]{x1D70F}}$ is a ‘drag’ length scale, $L_{\unicode[STIX]{x1D70F}}=(\int _{z_{t}}^{z_{c}}(z_{ws}-z)^{2}F_{D}\,\text{d}z/\int _{z_{t}}^{z_{c}}F_{D}\,\text{d}z)^{0.5}$ with $\int _{z_{t}}^{z_{c}}F_{D}\,\text{d}z=\unicode[STIX]{x1D70F}_{0}$ ; $L_{\unicode[STIX]{x1D719}}$ is a ‘depth-roughness’ length scale, $L_{\unicode[STIX]{x1D719}}=(3\int _{z_{t}}^{z_{ws}}(z_{ws}-z)^{2}(1-\unicode[STIX]{x1D719})\,\text{d}z)^{1/3}$ ; and $F_{in}^{\prime \prime }$ , $F_{s12}^{\prime \prime }$ and $F_{drp}^{\prime \prime }$ represent deviations of the inertial force ( $F_{in}^{\prime \prime }$ ), fluid forces in the $x$ – $y$ plane ( $F_{s12}^{\prime \prime }$ ) and pressure force ( $F_{drp}^{\prime \prime }$ ) per unit volume from their depth-averaged counterparts, i.e.

(2.4) $$\begin{eqnarray}\displaystyle F_{in}^{\prime \prime } & = & \displaystyle \unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}\left(\frac{\unicode[STIX]{x2202}\langle \overline{u}_{1}\rangle }{\unicode[STIX]{x2202}t}+\langle \overline{u}_{j}\rangle \frac{\unicode[STIX]{x2202}\langle \overline{u}_{1}\rangle }{\unicode[STIX]{x2202}x_{j}}\right)-\frac{1}{H}\int _{z_{t}}^{z_{ws}}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}\left(\frac{\unicode[STIX]{x2202}\langle \overline{u}_{1}\rangle }{\unicode[STIX]{x2202}t}+\langle \overline{u}_{j}\rangle \frac{\unicode[STIX]{x2202}\langle \overline{u}_{1}\rangle }{\unicode[STIX]{x2202}x_{j}}\right)\text{d}z\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle F_{s12}^{\prime \prime } & = & \displaystyle \unicode[STIX]{x1D70C}\left(-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}\langle \overline{u_{1}^{\prime }u_{j=1,2}^{\prime }}\rangle }{\unicode[STIX]{x2202}x_{j=1,2}}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}\langle \widetilde{\overline{u}}_{1}\widetilde{\overline{u}}_{j=1,2}\rangle }{\unicode[STIX]{x2202}x_{j=1,2}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j=1,2}}\unicode[STIX]{x1D719}\left\langle \unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}\overline{u}_{1}}{\unicode[STIX]{x2202}x_{j=1,2}}\right\rangle \right)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{H}\int _{z_{t}}^{z_{ws}}\unicode[STIX]{x1D70C}\left(-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}\langle \overline{u_{1}^{\prime }u_{j=1,2}^{\prime }}\rangle }{\unicode[STIX]{x2202}x_{j=1,2}}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}\langle \widetilde{\overline{u}}_{1}\widetilde{\overline{u}}_{j=1,2}\rangle }{\unicode[STIX]{x2202}x_{j=1,2}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j=1,2}}\unicode[STIX]{x1D719}\left\langle \unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}\overline{u}_{1}}{\unicode[STIX]{x2202}x_{j=1,2}}\right\rangle \right)\text{d}z\qquad\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle F_{drp}^{\prime \prime } & = & \displaystyle -\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}\langle \overline{p}\rangle }{\unicode[STIX]{x2202}x_{1}}-\frac{1}{H}\left(-\int _{z_{t}}^{z_{ws}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}\langle \overline{p}\rangle }{\unicode[STIX]{x2202}x_{1}}\,\text{d}z\right).\end{eqnarray}$$

It should be noted that similar to decomposition (2.3), the effects of turbulent and dispersive stresses on the friction at the wall have been accounted for in a recent work of Jelly & Busse (Reference Jelly and Busse2018) in an expression for the Hama roughness function, a parameter in the logarithmic velocity distribution for rough-wall flows (Hama Reference Hama1954). Unlike their study, our approach is more general as it does not require any assumptions regarding the applicability of the velocity log and defect laws and associated conditions imposed by these laws.

2.2 Friction factor decomposition and flow classification

Equation (2.3) explicitly shows that the Darcy–Weisbach friction factor $f$ can be split into at least five additive components, due to: (i) viscous stress; (ii) turbulent stress; (iii) dispersive stress; (iv) unsteadiness and spatial heterogeneity of double-averaged velocities and pressure; and (v) spatial heterogeneity of fluid stresses in the $x$ – $y$ plane. The proposed decomposition (2.3) can be viewed as a generalisation of the previous momentum-based decompositions (Fukagata et al. Reference Fukagata, Iwamoto and Kasagi2002; Peet & Sagaut Reference Peet and Sagaut2009; Bannier et al. Reference Bannier, Garnier and Sagaut2015), which are essentially focused on hydraulically smooth-wall flows, i.e. without accounting for important dynamic effects such as dispersive stresses and pressure drag that typically emerge in flows at a high roughness Reynolds number. It is therefore worth highlighting some key advantages of decomposition (2.3). First, equation (2.3) for the friction factor $f$ can be directly used for assessing relative contributions of particular mechanisms to bed resistance, which would be more difficult when using alternative resistance formulations such as Manning’s $n$ or Chezy’s $C$ due to their nonlinear relations to the Darcy–Weisbach $f$ . The second important feature of (2.3) is that it explicitly accounts for the scale of consideration, through the ‘filtering’ effect of the spatial averaging procedure, i.e. the ‘resolution’ scale is defined by the size of the averaging domain in (2.1). Third, all terms in (2.3) depend on the roughness geometry (via parameters $N$ and $\unicode[STIX]{x1D719}$ ) and highlight the significance of momentum transfer mechanisms in the near-bed flow region due to the role of weighting factors $(z_{ws}-z)$ and $(z_{ws}-z)^{2}$ . Fourth, equation (2.3) is fully consistent with the definitions of the conventional regimes of flat-bed laminar flow; hydraulically fully smooth-bed turbulent flow; and hydraulically fully rough-bed turbulent flow and transitions between them. Indeed, for steady, 2-D and uniform flat-bed laminar flow, equation (2.3) reduces to $f=24/Re$ , which is a classical relationship for flat-bed laminar OCF (e.g. Graf & Altinakar Reference Graf and Altinakar1998). For steady, 2-D uniform fully smooth-bed turbulent OCF (2.3) gives

(2.7) $$\begin{eqnarray}\displaystyle f=\frac{24}{Re}+\frac{24}{Q^{2}}\int _{z_{t}}^{z_{ws}}(z_{ws}-z)\langle -\overline{u_{1}^{\prime }u_{3}^{\prime }}\rangle \,\text{d}z, & & \displaystyle\end{eqnarray}$$

reflecting the emergence of the turbulent momentum flux in this flow type, especially pronounced in the near-bed region. For this flow type, the roughness height $\unicode[STIX]{x1D6E5}$ is much smaller than the thickness of the viscous sublayer $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}\sim \unicode[STIX]{x1D708}/u_{\ast }$ , i.e. the surface roughness in such flows is fully submerged within the viscous sublayer ( $\unicode[STIX]{x1D6E5}\ll \unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ ). For steady, 2-D uniform fully rough-bed turbulent flow ( $\unicode[STIX]{x1D6E5}\gg \unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ ) the resulting resistance equation incorporates, in addition to viscous and turbulent momentum fluxes, potential roughness effects in the form of the dispersive stress $\langle -\widetilde{\overline{u}}_{1}\widetilde{\overline{u}}_{3}\rangle$ , i.e.

(2.8) $$\begin{eqnarray}\displaystyle f=\frac{48}{Re}\frac{1}{N}+\frac{48}{Q^{2}N}\int _{z_{t}}^{z_{ws}}(z_{ws}-z)\unicode[STIX]{x1D719}\langle -\overline{u_{1}^{\prime }u_{3}^{\prime }}\rangle \,\text{d}z+\frac{48}{Q^{2}N}\int _{z_{t}}^{z_{ws}}(z_{ws}-z)\unicode[STIX]{x1D719}\langle -\widetilde{\overline{u}}_{1}\widetilde{\overline{u}}_{3}\rangle \,\text{d}z. & & \displaystyle\end{eqnarray}$$

Note that the dispersive stresses in fully rough-bed turbulent flow ( $\unicode[STIX]{x1D6E5}\gg \unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ ) contribute to the friction factor mainly in the domain outside the viscous sublayer, i.e. between and above the roughness elements. It is worth noting that the potential effect of dispersive stresses in fully rough-bed turbulent flows ( $\unicode[STIX]{x1D6E5}\gg \unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ ) follows from (2.3) as a prediction, i.e. their effect has not been involved previously in conventional considerations of such flows. Fifth, equation (2.3) also suggests the possibility of two additional flow regimes: rough-bed laminar flow and marginally rough-bed turbulent flow ( $\unicode[STIX]{x1D6E5}\approx \unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ ). For a steady, 2-D uniform rough-bed laminar flow, the turbulent stress contribution in (2.8) vanishes hence leading to

(2.9) $$\begin{eqnarray}\displaystyle f=\frac{48}{Re}\frac{1}{N}+\frac{48}{Q^{2}N}\int _{z_{t}}^{z_{ws}}(z_{ws}-z)\unicode[STIX]{x1D719}\langle -\widetilde{\overline{u}}_{1}\widetilde{\overline{u}}_{3}\rangle \,\text{d}z & & \displaystyle\end{eqnarray}$$

and thus exposing potential effects of dispersive stresses that are absent in conventional flat-bed laminar flows. For a steady, 2-D uniform marginally rough-bed turbulent flow ( $\unicode[STIX]{x1D6E5}\approx \unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ ), the expression for the friction factor is the same as (2.8) but with dispersive stresses now expected to act within the viscous sublayer, being non-existent above it. The marginally rough-bed turbulent flow ( $\unicode[STIX]{x1D6E5}\approx \unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ ) may be viewed as a distinct form of a conventional transitional regime between fully smooth-bed flow and fully rough-bed flow. It may occur if the dispersive stresses effect in (2.8) becomes significant. This flow type may well correspond to the riblet-bed flow, where the roughness height is known to be comparable to the thickness of the viscous sublayer, i.e. $\unicode[STIX]{x1D6E5}\approx \unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ .

Equation (2.3) can be further expanded to cover OCFs with mobile boundaries (both bed and water surfaces). Examples include river flows with surface waves and/or moving bed particles or waving aquatic plants on the bed. For the case of OCF with mobile boundaries, equation (2.3) will include additional terms describing potential correlations between time-averaged velocities and so called ‘time porosity’ (Nikora et al. Reference Nikora, Ballio, Coleman and Pokrajac2013). This extended equation for the friction factor further suggests the possibility of two additional flow types: mobile-rough-bed laminar flows and mobile-marginally rough-bed turbulent flows. However, the topic of mobile-boundary flows is outside the scope of this paper, which focuses on exploring the significance of viscous, turbulent, and dispersive contributions to the total bed friction for the case of fixed-rough-bed flows.

To explore decomposition (2.3), we have undertaken an extensive study of steady, uniform OCFs over self-affine rough beds that resemble a wide range of natural and man-made surfaces (see § 3.1). The key equation underpinning these tests is a reduced form of (2.3) that includes only the first three terms on the right-hand side, where the parameter $N$ is expressed as $N=3(L_{\unicode[STIX]{x1D70F}}/H)^{2}+(L_{\unicode[STIX]{x1D719}}/H)^{3}-(H_{m}/H)^{3}$ bearing in mind that $\unicode[STIX]{x1D70F}_{02D}/\unicode[STIX]{x1D70F}_{0}=1$ for the studied flows, which are essentially 2-D flows in terms of the double-averaged flow quantities. In addition, the third term (dispersive stress contribution) in (2.3) is subdivided into two terms, as explained in § 2.3 below.

Figure 1. Roughness-scale and SC-scale averaging domains in a rough-bed open-channel flow. Circles indicate secondary currents cells.

2.3 Partitioning of the total dispersive stress

There may be situations when there is more than one source of mean flow heterogeneity. For example, the roughness effects on the flow may be superimposed with the effects of the secondary currents (SCs) which are typically scaled with the flow depth. Under secondary currents in straight channels we consider helical motions visible in time- and streamwise-averaged flow fields (e.g. Nezu & Nakagawa Reference Nezu and Nakagawa1993). In cross-sectional planes they are typically expressed with counter-rotating motions preserved along the flow. In this case, the total dispersive stresses are formed due to combined effects of bed roughness and SCs. Close to the bed the main contribution to the dispersive stresses is expected from the roughness effects which typically quickly die out away from the bed, where the SCs may become the dominant contributor (e.g. figures 6–8 to be discussed in § 4.1). Thus, in principle the total dispersive stress in uniform double-averaged flow over (statistically) homogeneous in streamwise direction rough bed can be subdivided into two terms: (i) a term responsible for the roughness-scale effects; and (ii) a term responsible for the secondary current effects. This subdivision is motivated by the observations that OCFs over rough beds (e.g. Nikora & Roy Reference Nikora, Roy, Church, Biron and Roy2012) are often filled with helical secondary currents which may serve as an additional source of dispersive stresses typically occurring at a larger scale compared to the roughness-generated dispersive stresses. The subdivision makes use of two different averaging domains: (i) a thin elongated domain with a streamwise dimension that well exceeds the roughness correlation length, while the spanwise dimension is much smaller than the flow depth $H$ (or in other words much smaller than the transverse size of the SC cells; subsequently called the ‘roughness-scale’ domain); and (ii) a thin spanwise-extended domain that is several flow depths wide and several roughness correlation lengths long (subsequently called the ‘SC’ domain). A sketch in figure 1 illustrates both domains. Use of the first domain allows roughness effects to be isolated while with the second domain both roughness and SC contributions to the total dispersive stress are accounted for. If we use the double-averaged equations based on the SC domain, then by additionally employing a ‘roughness-scale’ averaging we can distinguish contributions to the total dispersive stresses due to the bed roughness and secondary currents (see appendix A for details). This approach yields the following decomposition of the total dispersive stress

(2.10)

where $\widetilde{\overline{u}}_{i}=\overline{u}_{i}-\langle \overline{u}_{i}\rangle =\widetilde{\langle \overline{u}_{i}\rangle }_{r}+\widetilde{\overline{u}}_{ir}$ is a dispersive fluctuation based on the SC averaging domain; $\langle \overline{u}_{i}\rangle$ is the double-averaged velocity based on the SC averaging domain, defined by angular brackets $\langle \cdots \rangle$ with no subscript; $\widetilde{\overline{u}}_{ir}=\overline{u}_{i}-\langle \overline{u}_{i}\rangle _{r}$ is a dispersive component based on the ‘roughness-scale’ averaging domain; $\langle \overline{u}_{i}\rangle _{r}$ is the double-averaged velocity based on the ‘roughness-scale’ averaging domain, defined by angular brackets $\langle \cdots \rangle$ with subscript $r$ ; $\widetilde{\langle \overline{u}_{i}\rangle }_{r}=\langle \overline{u}_{i}\rangle _{r}-\langle \overline{u}_{i}\rangle$ is a dispersive component due to the secondary currents alone. The final equation underpinning our estimates in §§ 4.2–4.4 can then be obtained from (2.3) and (2.10) as

(2.11)

It should be noted that the dispersive stress partitioning, proposed in this section, stems from the definition of the secondary currents in OCFs as helical motions visible in time- and streamwise-averaged flow fields (e.g. Nezu & Nakagawa Reference Nezu and Nakagawa1993). However, this conventional definition is not unique, as in principle there may be ‘meandering’ time-averaged SC cells, among other possible scenarios. After streamwise averaging of such meandering cells they may appear invisible in time- and streamwise-averaged flow fields, invalidating the proposed partitioning. Furthermore, the streamwise wavelengths of SC meandering may be comparable to the roughness scales, making the separation of roughness-induced and SC-induced contributions problematic, if possible at all. This later case is suggested by a recent paper of Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2018), who reported a comprehensive study of flow structure in pipes covering a wide range of three-dimensional roughnesses of the staggered type. The emergence of near-wall ‘meandering’ SCs, with scales comparable to the scales of roughness elements, is clearly seen in their figure 14 (Chan et al. Reference Chan, MacDonald, Chung, Hutchins and Ooi2018). That said, the SCs observed in our work, as will be shown in § 4.1, fit the conventional definition very well making the partitioning proposed here appropriate.

Table 1. Flow and roughness characteristics for LES data set. $\unicode[STIX]{x1D6FD}$ is spectral slope, $H$ is mean flow depth, $U$ is cross-sectionally averaged mean velocity, $u_{\ast }$ is shear velocity computed using streamwise pressure gradient as explained in § 3.3, $H/\unicode[STIX]{x1D6E5}$ is relative submergence, $\unicode[STIX]{x1D6E5}=4\unicode[STIX]{x1D70E}_{z}$ , $B/H$ is aspect ratio, $B$ is flow width, $\unicode[STIX]{x1D6E5}^{+}=u_{\ast }\unicode[STIX]{x1D6E5}/\unicode[STIX]{x1D708}$ is roughness Reynolds number, $Re_{b}=UH/\unicode[STIX]{x1D708}$ is bulk Reynolds number, $Fr=U/(gH)^{0.5}$ is Froude number, $\unicode[STIX]{x1D708}$ is fluid kinematic viscosity and $g$ is acceleration due to gravity.

3.1 Roughness characteristics

For our experiments and simulations we have selected single-valued self-affine fractal surfaces as they seem to represent a fairly general model of natural and man-made rough surfaces. Indeed, the studies of gravel-bed and sand-dune-bed rivers (e.g. Hino Reference Hino1968; Nikora, Sukhodolov & Rowinski Reference Nikora, Sukhodolov and Rowinski1997; Nikora, Goring & Biggs Reference Nikora, Goring and Biggs1998; Butler, Lane & Chandler Reference Butler, Lane and Chandler2001), the ocean floor (e.g. Bell Reference Bell1975), glaciated bedrocks and sub-glacial surfaces (e.g. Hubbard, Siegert & McCarroll Reference Hubbard, Siegert and McCarroll2000; Mankoff et al. Reference Mankoff, Gulley, Tulaczyk, Covington, Liu, Chen, Benn and Głowacki2017), the surfaces of other planets (e.g. Turcotte Reference Turcotte1997; Nikora & Goring Reference Nikora and Goring2004) and a variety of machined surfaces (e.g. Majumdar & Tien Reference Majumdar and Tien1990) have demonstrated that all these surfaces can be described with one-dimensional power spectra (or second-order structure functions) of a similar shape. Typically, the observed spectra include three scale ranges with distinctly different behaviour: (i) a ‘saturation’ (white spectrum) range at low wavenumbers ( $k_{x}<k_{xc1}$ ); (ii) a scaling region at intermediate wavenumbers ( $k_{xc1}<k_{x}<k_{xc2}$ ) where the spectrum $S(k_{x})$ decays as a power function with a spectral exponent $\unicode[STIX]{x1D6FD}$ , i.e.  $S(k_{x})\propto k_{x}^{-\unicode[STIX]{x1D6FD}}$ ; and (iii) a ‘smooth’ region at high wavenumbers ( $k_{x}>k_{xc2}$ ) where the spectrum rapidly declines to zero (and thus within this region the surface is ‘smooth’, i.e. differentiable). Note that this schematisation assumes that the surface is a continuous single-valued field, which is a fair assumption even for granular surfaces if grain spacings are much smaller than the prevailing scales of surface fluctuations. The cutoff wavenumbers $k_{xc1}$ and $k_{xc2}$ define the upper and lower boundaries of the scaling region, respectively, with $(k_{xc1})^{-1}$ providing a measure of the roughness correlation length scale (which also depends on $\unicode[STIX]{x1D6FD}$ ). The surface with such features may be classified as representative of self-affine fractal surfaces that exhibit, at least within a certain scale range, a power dependence of the amplitude $A$ of a disturbance on its wavelength $L$ , i.e.  $A\propto L^{\unicode[STIX]{x1D6FC}}$ (e.g. Turcotte Reference Turcotte1997). The exponent $\unicode[STIX]{x1D6FC}$ , known as the Hurst (or Hölder) exponent, relates to the spectral exponent $\unicode[STIX]{x1D6FD}$ as $\unicode[STIX]{x1D6FD}=2\unicode[STIX]{x1D6FC}+1$ . The exponent $\unicode[STIX]{x1D6FC}$ can take a value between 0 and 1, with $\unicode[STIX]{x1D6FC}=1$ corresponding to the special case of self-similarity, thus yielding a range for $\unicode[STIX]{x1D6FD}$ from 1 to 3 for a surface to be classified as self-affine fractal (e.g. Turcotte Reference Turcotte1997). The described three-range spectral model is adopted in our work as most representative of naturally formed and man-made surfaces. It may be viewed as an extension of the approaches used in the studies of Anderson & Meneveau (Reference Anderson and Meneveau2011) and Barros, Schultz & Flack (Reference Barros, Schultz and Flack2015, Reference Barros, Schultz and Flack2018), where spectra of rough surfaces had a single range corresponding to our range (ii). For our study, we have chosen isotropic self-affine Gaussian surfaces with one-dimensional spectra as described above. Three self-affine fractal patterns, referred to hereafter as R1, R2 and R3, were generated using the inverse discrete Fourier transform of the designed spectra. Full details of the applied spectral synthesis procedure and the algorithms used are provided in Stewart et al. (Reference Stewart, Cameron, Nikora, Zampiron and Marusic2019).

For all three surface types, the wavenumbers defining the location and extent of the scaling range were the same, i.e.  $k_{xc1}=0.02~\text{mm}^{-1}$ (i.e.  $l_{x}=1/k_{xc1}=50~\text{mm}$ ), $k_{xc2}=0.20~\text{mm}^{-1}$ (i.e.  $l_{sm}=1/k_{xc2}=5~\text{mm}$ ), as well as the surface standard deviation $\unicode[STIX]{x1D70E}_{z}=1.5~\text{mm}$ . The exponent $\unicode[STIX]{x1D6FD}$ in the scaling range (ii) is selected to be 1 for R1 (resembling ‘ $-1$ ’ scaling in near-bed velocity spectra), $5/3$ for R2 (resembling the inertial subrange exponent) and 3 for R3 (resembling scaling in sand-dune spectra). Although the selection of the exponents was rather intuitive, they do cover the whole range of possibilities. Figure 2 shows the resulting roughness patterns demonstrating the strong influence of the exponent $\unicode[STIX]{x1D6FD}$ while keeping scaling ranges and surface standard deviation the same for all three cases. The patterns are periodic in both $x$ and $y$ directions and thus can continuously cover the bed of our open-channel facility as well as to serve as a ‘numerical bed’ for LES runs (see §§ 3.2 and 3.3). The roughness height $\unicode[STIX]{x1D6E5}$ is defined as four standard deviations of the bed elevations to represent a 95.4 % range of their fluctuations, i.e.  $\unicode[STIX]{x1D6E5}=4\unicode[STIX]{x1D70E}_{z}=6~\text{mm}$ . We do not employ directly the equivalent roughness height $k_{s}$ to characterise bed roughness in our study as its use requires the presence of the genuine logarithmic layer and validity of the logarithmic velocity distribution. Although our data do exhibit quasi-logarithmic velocity distribution above the roughness tops (see § 4.1), the relative submergence $H/\unicode[STIX]{x1D6E5}=5.7{-}20$ in our experiments is well below a suggested threshold value of 40 for the logarithmic layer to emerge (Jiménez Reference Jiménez2004). That said, it should be noted that the correlations obtained by Flack & Schultz (Reference Flack and Schultz2010) and Barros et al. (Reference Barros, Schultz and Flack2018) suggest that for the Gaussian rough surfaces the following approximate relation is valid: $k_{s}\approx (3.5{-}4.5)\unicode[STIX]{x1D70E}_{z}$ . Its equivalent presentation $k_{s}\approx (0.88{-}1.13)\unicode[STIX]{x1D6E5}$ provides additional support for choosing $\unicode[STIX]{x1D6E5}=4\unicode[STIX]{x1D70E}_{z}$ for characterising roughness heights in our study. Besides, Stewart et al. (Reference Stewart, Cameron, Nikora, Zampiron and Marusic2019) also showed that the spectral slope $\unicode[STIX]{x1D6FD}$ appears to be closely related to the roughness ‘effective slope’ ES, which is equal to the double solidity $\unicode[STIX]{x1D706}$ (i.e. the ratio of the projected frontal area of the roughness elements to the wall-parallel projected area; Napoli, Armenio & De Marchis Reference Napoli, Armenio and De Marchis2008; Flack & Schultz Reference Flack and Schultz2010). Estimates of solidity $\unicode[STIX]{x1D706}$ corresponding to $\unicode[STIX]{x1D6FD}=$ 1, 5/3 and 3 were found to be $\unicode[STIX]{x1D706}=0.327$ , 0.219 and 0.102, respectively. The obtained data suggest an approximate relation $\unicode[STIX]{x1D706}=0.5ES=0.58\exp (-0.58\unicode[STIX]{x1D6FD})$ which, however, is not general and is limited to the roughness spectral shapes explored in our experiments.

Figure 2. Self-affine rough tiles used in the experiments and LES (a); (b) one-dimensional streamwise transects; the origin of the $z_{b}$ axis corresponds to the mean bed level. Tiles manufactured for laboratory experiments are shown.

For the laboratory experiments, the digital patterns have been used to manufacture the roughness tiles employing a three-step mould and cast procedure. First, a master plate was created for each of R1, R2 and R3 from acetal copolymer using a high-precision three-axis CNC milling machine. Second, moulds of each master plate were produced using a two-part addition cure room-temperature-vulcanizing (RTV) silicone rubber. Third, replicas of the master plates were then cast in the silicone moulds from epoxy resin. Owing to the large bed surface area of our laboratory facility it was necessary to manufacture at least 150 of these casts (‘tiles’) for each roughness pattern. As part of the final adjustments, the dimensions of the plates were reduced from the initial design value of 392 mm  $\times$  392 mm to a final value of 388 mm $\times$ 388 mm ( $\pm$ 0.1 mm). The prepared plates were aligned on the flume bed in a 50 $\times$ 3 array creating a continuous rough bed of repetitive patterns. The magnitude of potential discontinuities at the joints between the plates was appreciably less than $l_{sm}$ . A detailed assessment and analysis of the manufactured plates, as well as extensive measurements of the friction factor for a wide range of hydraulic conditions, are reported in Stewart et al. (Reference Stewart, Cameron, Nikora, Zampiron and Marusic2019). Specifically, Stewart et al. (Reference Stewart, Cameron, Nikora, Zampiron and Marusic2019) found that the exponent $\unicode[STIX]{x1D6FD}$ should be considered as a significant factor of the overall hydraulic resistance, in addition to the conventionally used length scales such as characteristic roughness heights and wavelengths. The results demonstrated that with all other parameters being the same, a decrease in $\unicode[STIX]{x1D6FD}$ leads to a visible increase in the friction factor. This effect was observed to be as great as 49 % between the R1 and R3 beds and it was not dependent on the channel aspect ratio. The present paper is, to a certain degree, a follow up of the work reported in Stewart et al. (Reference Stewart, Cameron, Nikora, Zampiron and Marusic2019), with a particular focus on the mechanisms contributing to the overall friction factor.

3.2 Laboratory experiments

Experiments were conducted in the Aberdeen Open Channel Facility (AOCF). The AOCF consists of an 18 m long and 1.18 m wide recirculating open-channel flume, an instrumental carriage that traverses the length of the flume with an integrated particle image velocimetry (PIV) system, along with data capture and analysis hardware and software systems. Flow is fed to the flume entrance tank by a pair of centrifugal pumps which are speed controlled in a closed loop system, with feedback provided by magnetic encoders attached to the pump drive shafts. The entrance tank has been optimised with honeycomb panels and guide vanes to ensure that the flow enters the channel uniformly distributed and free from large-scale turbulence and surface waves. The tail water level in the channel is controlled by adjusting a vertical slat weir at the channel exit to ensure flow uniformity. The flow depth distribution along the channel is calculated and adjusted before each experiment from profiles of the bed surface (using a laser triangulation sensor) and water surface (using a confocal displacement sensor). The relative standard errors of measured flow depth $H$ , bed slope $S_{b}$ and depth-averaged flow velocity $\widehat{U}$ (from PIV, see below) did not exceed 0.3 %, 1.4 % and 0.3 %, respectively, resulting in the maximum relative error in the friction factor of 1.5 %. This estimate gives the upper bound and corresponds to the overlap (with LES) validation cases that had smallest bed slope. The relative errors for all high Reynolds number cases are 5–8 times smaller.

Two configurations of the AOCF PIV system were used in this study: a streamwise orientated ‘mini’ mode for near-bed measurements, including the region below roughness tops and a transverse orientated ‘macro’ mode for measuring the flow region above the roughness tops. The ‘mini’ mode (figure 3 a) is a two-camera stereoscopic system with a measurement domain of $50~\text{mm}\times 25~\text{mm}$ ( $x\times z$ ). Three streamlined glass bottom boats sit at the water surface (one for each camera, and one for the laser sheet) to allow optical access to the flow. The camera boats contain water filled prisms to limit the distortion and internal reflection at the air–water interface. Flow disturbance from the boats is limited to a thin boundary layer ( ${<}1.5$  mm thick) which does not interact with the measured region of the flow. We use Dalsa 4M180 cameras equipped with Nikon 60 mm lenses with aperture set to f/16, and capture data directly to a 24 TB solid state disk array. The light sheet was 0.7 mm thick and generated by an Oxford Lasers 532 nm Nd:YAG laser with pulse energy up to 100 mJ. The seeding particles for the ‘mini’ mode were neutrally buoyant 10–20 micron titania coated hollow glass microspheres, produced by Microsphere Technology Ltd.

Figure 3. ‘Mini’-mode PIV configuration (a); ‘macro’-mode PIV configuration (b); footprint of the PIV measurement planes, with lines across the flow showing ‘macro’ mode and lines along the flow presenting ‘mini’-mode measurement planes, black-border squares show the roughness tiles (c); and transfer function $T$ of the ‘mini’ and ‘macro’ PIV measurement modes, with $\unicode[STIX]{x1D706}$ , $\unicode[STIX]{x1D706}_{x}$ , $\unicode[STIX]{x1D706}_{y}$ and $\unicode[STIX]{x1D706}_{z}$ denoting a generic wavelength and directional wavelengths in $x$ , $y$ and $z$ directions, respectively (d).

The ‘macro’ mode (figure 3 b) is a four-camera stereoscopic PIV system with a measurement domain of $225~\text{mm}\times 120~\text{mm}$ ( $y\times z$ ). For this mode we have designed and constructed semi-immersible camera lenses that penetrate approximately 10 mm through the water surface. The lenses have a 7 mm right angle prism mirror attached to the tip, which orientates the viewing direction of the downward facing cameras in the upstream direction. The cameras are placed 500 mm downstream of the transverse orientated light sheet so that the small wakes from the semi-immersed lenses cannot disturb the upstream measurement domain. The lenses have a fixed aperture of f/16. Larger silver coated seeding particles with diameters of 20–32 microns were used for the ‘macro’ mode. The light sheet was 2.8 mm thick and entered the flow through the glass side wall of the channel.

Both ‘mini’ and ‘macro’-mode PIV systems are attached to a robotic 3-axis platform that can programmatically position the respective measurement domains anywhere in the channel without manual adjustments or recalibration. Axis positions are recorded with sub-micron noise and resolution by linear-scale optical encoders attached to each axis. Calibrations of the stereoscopic PIV systems are performed by capturing images of a backlit 50 micron diameter pinhole from different camera positions. The experimental scenarios covered by PIV measurements are shown in table 2.

A series of 128 ‘mini’-mode and 32 ‘macro’-mode measurements were carried out for each flow scenario and roughness configuration to provide sufficient data for double averaging. The measurements were made 12 m (100–300 flow depths) from the flume entrance ensuring fully developed flow conditions. The ‘mini’-mode measurement planes were configured in an 8  $\times$  16 ( $x\times y$ ) grid with a step size of $40~\text{mm}\times 15~\text{mm}$ . Each measurement plane was recorded for 100 s at a sampling rate of 10 Hz. The combined footprint of the 128 planes is shown in figure 3(c). The 32 ‘macro’ mode measurement planes were configured on a grid of 8  $\times$  4 ( $x\times y$ ) with a grid step of $40~\text{mm}\times 165~\text{mm}$ (figure 3 c). Each of the ‘macro’ measurement planes were sampled for 240 s at a sampling rate of 10 Hz.

The PIV images were analysed using an iterative deformation method (IDM) PIV algorithm (Cameron Reference Cameron2011). For our analysis we have selected 64  $\times$  64 pixel Blackman weighted interrogation regions with 52 pixel overlap between neighbouring interrogation regions, and a stretched sin-cardinal (sinc) kernel to low-pass filter the velocity field between each of the 8 iterations of the algorithm. The transfer functions that reflect the low-pass filtering of the velocity field from the combined effects of the interrogation region size and the light sheet thickness is illustrated in figure 3(d). The mini-mode data have been used to obtain the double-averaged quantities using the roughness-scale averaging domain ( $0.7\times 0.7\times 388$ mm). The macro-mode data were employed to compute the SC-scale double-averaged quantities within the domain ( $2.5~\text{mm}\times 388~\text{mm}\times 388~\text{mm}$ ) at the channel centre (central tile in figure 3 c), including the depth-averaged velocity $\widehat{U}=H^{-1}\int _{z_{t}}^{z_{ws}}\unicode[STIX]{x1D719}\langle \overline{u}_{1}\rangle \,\text{d}z$ and ‘local’ friction factor $f=8\unicode[STIX]{x1D70F}_{0}/(\unicode[STIX]{x1D70C}\widehat{U}^{2})$ defined in §§ 1 and 2.1. The averaging times in the data analysis significantly exceeded turbulence length scales, typically being well above $70H/u_{\ast }$ . Similar domain sizes and averaging times have been used in LES data analysis, which is described below.

3.3 Large-eddy simulations

To carry out the LESs for selected scenarios (table 1), the in-house code Hydro3D was employed. The code has been previously validated and applied to study several flows of similar complexity (Bomminayuni & Stoesser Reference Bomminayuni and Stoesser2011; Kara, Stoesser & Sturm Reference Kara, Stoesser and Sturm2012; Bai, Fang & Stoesser Reference Bai, Fang and Stoesser2013; Kim, Stoesser & Kim Reference Kim, Stoesser and Kim2013; Kara et al. Reference Kara, Kara, Stoesser and Sturm2015b ,Reference Kara, Stoesser, Sturm and Mulahasan c ; Fraga et al. Reference Fraga, Stoesser, Lai and Socolofsky2016; McSherry et al. Reference McSherry, Chua, Stoesser and Mulahasan2018). Hydro3D is based on a finite-difference method employing a staggered Cartesian grid. The code solves the filtered Navier–Stokes equations for incompressible unsteady flow, i.e.

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{i}}=0 & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}v_{i}v_{j}}{\unicode[STIX]{x2202}x_{j}}=-\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}{\dot{p}}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}(2\unicode[STIX]{x1D708}\unicode[STIX]{x1D61A}_{ij})}{\unicode[STIX]{x2202}x_{j}}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{j}}, & \displaystyle\end{eqnarray}$$

where $v_{i}$ is the resolved velocity in the $i$ th direction ( $i=1$ , 2 and 3 correspond to $x$ -, $y$ - and $z$ -directions, respectively) and $x_{i}$ is the space coordinate vector; ${\dot{p}}$ is resolved pressure; and $\unicode[STIX]{x1D61A}_{ij}$ denotes the filtered strain-rate tensor $\unicode[STIX]{x1D61A}_{ij}=1/2(\unicode[STIX]{x2202}v_{i}/\unicode[STIX]{x2202}x_{j}+\unicode[STIX]{x2202}v_{j}/\unicode[STIX]{x2202}x_{i})$ . The subgrid-scale (SGS) stress $\unicode[STIX]{x1D70F}_{ij}$ is defined as $\unicode[STIX]{x1D70F}_{ij}=-2\unicode[STIX]{x1D708}_{t}\unicode[STIX]{x1D61A}_{ij}$ where $\unicode[STIX]{x1D708}_{t}$ is the eddy viscosity. In this study, the wall-adapting local eddy viscosity (WALE) proposed by Nicoud & Ducros (Reference Nicoud and Ducros1999) is used to calculate $\unicode[STIX]{x1D708}_{t}$ and to obtain the subgrid stress model.

In Hydro3D, the convection and diffusion terms in the Navier–Stokes equations are approximated by fourth-order accurate central differences. An explicit 3-step Runge–Kutta scheme is used to integrate the equations in time, providing second-order accuracy. A fractional step method is employed, i.e. within the time step the convection and diffusion terms are solved explicitly first in a predictor step that is then followed by a corrector step during which the pressure and divergence-free velocity fields are obtained through a Poisson equation. The latter is solved iteratively through a multi-grid procedure. Hydro3D features a local mesh refinement (LMR) method allowing a finer mesh in the vicinity of the rough beds. As a result of using the fine mesh near the bed, the SGS eddy viscosity is of the same order of magnitude as the molecular viscosity even for the high-Reynolds-number validation (overlap) cases 5.a, 10.a and 15.a (table 1). Details of the LMR method and its validation are reported in Cevheri, McSherry & Stoesser (Reference Cevheri, McSherry and Stoesser2016).

Twenty-six LES runs were carried out in total. Details of the computational set-ups are provided in table 3, where the first column presents the run number, with a, b, c and d identifying different flow submergences $H/\unicode[STIX]{x1D6E5}$ . The second column shows the relative dimensions of the computational domain ( $L_{x}$ , $L_{y}$ , $L_{z}$ ) in terms of the maximum flow depth $H_{m}$ for each simulation. In general, the simulation domain covered two roughness tiles in the streamwise direction and one roughness tile in the spanwise direction; for the deepest flows this results in a domain size of the order $2\unicode[STIX]{x03C0}H_{m}\times \unicode[STIX]{x03C0}H_{m}\times H_{m}$ , typically accepted as adequate for resolving the largest structures in the flow. The third column of table 3 provides the total number of grid points in the streamwise, $n_{x}$ , spanwise, $n_{y}$ , and wall-normal, $n_{z}$ , directions, respectively. The domain is discretised with a locally refined mesh, a section of which is depicted in figure 4(a). The mesh consists of two parts: a coarser mesh (half of the resolution of the fine mesh, which makes it a 2:1 LMR) is employed in the upper region (80 % of the domain) whereas a finer grid is used in the near-wall region (20 % of the domain), to (i) allow the accurate representation of the roughness geometry and (ii) sufficiently resolve the flow features around individual roughness elements and hence their contributions to the drag. In the fourth column of table 3, the grid spacings in terms of wall units $\unicode[STIX]{x0394}x^{+}$ , $\unicode[STIX]{x0394}y^{+}$ , $\unicode[STIX]{x0394}z^{+}$ are listed for the finest grid resolution. The flow in the simulations is driven by a streamwise pressure gradient that is adjusted at every time step so that the constant bulk velocity is maintained. The LES-computed double- and depth-averaged pressure gradient $\text{d}\widehat{P}/\text{d}x$ is then used to calculate the domain-averaged wall shear stress and shear velocity, i.e. $\unicode[STIX]{x1D70F}_{0}=(\text{d}\widehat{P}/\text{d}x)H$ and $u_{\ast }=(\unicode[STIX]{x1D70F}_{0}/\unicode[STIX]{x1D70C})^{0.5}=[(1/\unicode[STIX]{x1D70C})(\text{d}\widehat{P}/\text{d}x)H]^{0.5}$ .

Figure 4. Part of the computational grid for $\unicode[STIX]{x1D6FD}=3$ (a); LES-predicted friction factors as a function of the bulk Reynolds number (symbols) and known analytical curves (lines) for laminar flows and smooth-bed turbulent flow (b); comparison of the friction factor estimates using experimental and LES data (see § 3.3 for details) (c,d); and comparison of velocity patterns from LES and PIV measurements; time- and streamwise-averaged data are shown (e,f). The coordinate $y=0$ corresponds to the right side wall of the flume.

In all LES runs, periodic boundary conditions are applied in the streamwise and spanwise directions, mimicking fully developed flow in both these directions. The water surface is treated as a frictionless rigid lid at which the free-slip condition is applied. The use of the rigid lid condition is deemed a reasonable treatment given the relatively low Froude number (table 1) and the experimentally observed fairly flat water surface. The no-slip condition is applied at the rough surfaces, employing a refined version (Kara, Stoesser & McSherry Reference Kara, Stoesser and McSherry2015a ; Ouro et al. Reference Ouro, Harrold, Stoesser and Bromley2017; Ouro & Stoesser Reference Ouro and Stoesser2017) of the direct-forcing immersed-boundary method (IBM) first proposed by Uhlmann (Reference Uhlmann2005). The maximum grid spacing near the (rough) wall is $\unicode[STIX]{x0394}z^{+}/2$ being less than 5 wall units. Thus, the grid is fine enough for resolving, at least approximately, the viscous sublayer. The LES were run for at least 50 eddy turnover times (calculated as the ratio of water depth to the shear velocity) to collect the required statistics. In addition, the LES data were averaged spatially and convergence of the quantities was carefully monitored. All simulations are considered to be sufficiently converged.

In order to verify the code’s ability to predict the friction factor accurately, three series of runs for turbulent and laminar OCFs were carried out, i.e. (i) laminar flows; (ii) smooth-bed turbulent flow; and (iii) rough-bed turbulent flows. Figure 4(b) shows the LES-computed friction factors for laminar flows (open circles) as a function of the Reynolds number for five runs, together with the analytical solution $f=24/Re_{b}$ (solid line, $Re_{b}=UH/\unicode[STIX]{x1D708}$ ; Graf & Altinakar Reference Graf and Altinakar1998). As can be seen in this figure, the simulations reproduce the friction factor change with increasing Reynolds number fairly accurately. A simulation of turbulent smooth-bed flow ( $Re_{b}=10\,400$ ) shows that the friction factor values, evaluated using $\text{d}\widehat{P}/\text{d}x$ (open circle) and from (2.7) (black triangle), are essentially the same. They both match very well the Blasius empirical equation $f=0.224Re_{b}^{-1/4}$ , which was initially proposed for smooth-wall turbulent pipe flows and later widely adopted for wide open-channel flows (e.g. Yen Reference Yen2002). The third validation data set includes LES runs 5a, 5b, 10a, 10b, 15a and 15b that resemble experimental runs 24, 20, 29, 25, 37 and 30, respectively (tables 1 and 2). Figure 4(c,d) compares the estimates of the friction factors computed from $f_{H}=8(gHS_{b})/\widehat{U}^{2}$ (figure 4 c) and from (2.8), $f$ , (figure 4 d). Figure 4(c,d) demonstrates that the LES-predicted friction factor estimates match the experimentally obtained values reasonably well, with the data points being within a 10 % envelope except for the two shallow flows. The LES data for deeper flows (smaller friction factor) tend to match the experimental data better. A possible reason for this is that the rigid lid assumption is more appropriate for deeper flows, where water surface deformations are sufficiently small, while their effect may not be negligible for shallow flows. Also, the spanwise flow boundaries for LES and in the experiments are different (‘infinitely wide channel’ in the LES and smooth side walls in the experiments). This difference may explain why the experimental values of $f$ are always slightly larger than the LES values. Figure 4(e,f) presents contours of the time- and streamwise-averaged velocity $\langle \overline{u}_{1}\rangle _{r}$ with velocity vectors in the transverse plane for LES (e) and PIV (f) data. The agreement between simulations and experiments is favourable, with both revealing cells of the secondary currents located at approximately the same spanwise positions. The discrepancies in the patterns are mainly due to the spanwise boundaries and their effect on the development of the secondary currents in terms of strength, position and size. Overall, from figure 4 it can be concluded that the LES are able to predict flows over rough beds with sufficient accuracy, allowing extension of the experimental scenarios towards smaller Reynolds numbers using LES.