Free surface flow over square bars at intermediate relative submergence

ABSTRACT Results from large-eddy simulations and complementary flume experiments of turbulent open channel flows over bed-mounted square bars at intermediate submergence are presented. Scenarios with two bar spacings, corresponding to transitional and k-type roughness, and three flow rates, are investigated. Good agreement is observed between the simulations and the experiments in terms of mean free surface elevations and mean streamwise velocities. Contours of simulated time-averaged streamwise, streamfunction and turbulent kinetic energy are presented and these reveal the effect of the roughness geometry on the water surface response. The analysis of the vertical distribution of the streamwise velocity shows that in the lowest submergence cases no logarithmic layer is present, whereas in the higher submergence cases some evidence of such a layer is observed. For several of the flows moderate to significant water surface deformations are observed, including weak and/or undular hydraulic jumps which affect significantly to the overall streamwise momentum balance. Reynolds shear stress, form-induced stress and form drag are analysed with reference to the momentum balance to assess their contributions to the total hydraulic resistance of these flows. The results show that form-induced stresses are dominant at the water surface and can contribute significantly to the overall drag, but the total resistance in all cases is dominated by form drag due to the presence of the bars.


Introduction
Although important advances have been made in recent decades, a comprehensive understanding of the turbulence structure and the flow resistance in shallow rough-bed flows remains elusive. In civil and environmental engineering applications the situation is further complicated by the fact that most shallow flows have steep slopes and free surface deformations can be significant. Bathurst (1985) summarized that most standard resistance equations are applicable only to rivers with gentle slopes, and those that have been developed for steep stretches are either empirical or valid only for deep flows. In many shallow flows of practical interest the mean flow depth, H, is of the same order of magnitude as the roughness height, k, and the well-known concept of the logarithmic boundary layer (LBL) may no longer be applicable (Jiminez, 2004;Raupach, Antonia, & Rajagopalan, 1991). The LBL is nonetheless widely used to study such flows, often leading to erroneous and misleading interpretations of the flow physics.
Several early attempts at modifying the LBL for application to gravel bed rivers with steep slopes and low submergences resulted in a number of empirical or semi-empirical formulations for the well known Darcy-Weisbach friction factor, f, which is defined for flow in an open channel as follows: where U b is the bulk flow velocity, g is acceleration due to gravity and S is the bed slope. A brief review of these approaches was given by Bathurst (1985), who considered a semi-empirical equation proposed by Hey (1979) to be the most complete. However, based on a thorough analysis of field and flume data from his own measurement campaigns and those performed by others, Bathurst concluded that the Hey (1979) equation was subject to large errors in f for shallow flows. Bathurst therefore proposed an alternative equation that was more applicable to low submergence flows but was nonetheless characterized by a possible error of between ± 25% and ± 35%. Bathurst (1985) was left to conclude that the complicated nature of the flow resistance processes, coupled with the lack of available data, currently prevent the development of a satisfactory practical method for predicting flow resistance in steep mountain rivers. More recently, Pagliara & Chiavaccini (2006) and Pagliara, Das, and Carnacina (2008) investigated flow resistance in steep chutes with large-scale roughness that comprised protruding boulders in some cases. Bed slope, relative submergence and boulder arrangement were systematically altered and a logarithmic expression for f was obtained and was shown to fit the available experimental data well. Furthermore, it was observed that f increased with increasing slope for a given relative submergence and roughness arrangement and the influence of larger boulders on the overall resistance was shown to decrease with increasing submergence. In order to achieve a rigorous understanding of low submergence flows over rough beds, Nikora, Goring, McEwan, and Griffiths (2001) and Nikora et al. (2007a,b) applied the doubleaveraging methodology to the governing Navier Stokes equations. The follows the concept of the Reynolds decomposition, e.g. θ =θ + θ , for an instantaneous variable θ, where the overline denotes the temporally-averaged value and the prime is used to denote a fluctuation due to turbulence. The double-averaging decomposition, e.g.θ = θ +θ , is used for temporallyaveraged variables where the angular brackets denote the spatially-averaged value and the tilde denotes the spatial fluctuation of that variable. The spatial averaging procedure is generally carried out across a plane parallel to the bed. Nikora et al. (2001) and Nikora et al. (2007a,b) identified four distinct flow regimes based on relative submergence. Flow type I (high submergence) comprises the outer, logarithmic, form-induced and interfacial sub-layers; flow type II (intermediate submergence) is characterized by the absence of a logarithmic layer because H /k is too small to sustain it; in flow type III (low submergence) the roughness layer ( = interfacial + form-induced) extends to the free surface; and in flow type IV (partiallyinundated roughness) only the interfacial sub-layer is present. Nikora (2009) extended the use of the double-averaging methodology to develop a theoretical expression that explicitly showed that the friction factor can be accounted for by six additive components, of which only three are present in two-dimensional uniform spatially-averaged flow without secondary currents: viscous stress, turbulent stress and forminduced stress. The contribution of the form-induced stress, ũw , has recently attracted substantial attention as it has been shown to play an important role in the overall streamwise momentum balance in the near-bed region (Aberle, Koll, & Dittrich, 2008;Dey & Das, 2012;Ferreira et al., 2010;Gimenez-Curto and Corniero, 1996;Gimenez-Curto and Corniero Lera, 2003;Manes, Pokrajac, Coceal, & McEwan, 2008;Manes, Pokrajac, & McEwan, 2007;Manes, Pokrajac, McEwan, & Nikora, 2009;Mignot, Barthelemy, & Hurther, 2009a,b;Nikora et al., 2007b;Sarkar & Dey, 2010). Indeed, Gimenez-Curto and Corniero (1996) and Gimenez-Curto and Corniero Lera (2003) suggested that at very low submergences the form-induced stress could become the dominant component in the streamwise momentum balance. Interestingly Manes et al. (2007) found that the relative contribution of form-induced stress increased as submergence decreased in their experiments on open channel flow over closely-packed spheres.
Alongside the question of relative submergence, the geometrical characteristics of the roughness itself clearly play an important role. Perry, Schofield, and Joubert (1969) were the first to identify two types of roughness, d-type and k-type, which give rise to fundamentally different roughness functions. Their pipe flow experiments showed that for k-type roughness the roughness function depended on the size of the roughness elements, whereas in d-type roughness it depended on the pipe diameter. For a flow of given bulk Reynolds number, R, over two-dimensional roughness elements, the transition between d-and k-type roughness depends solely on the streamwise spacing of the elements. Many researchers have chosen to focus on square bar roughness to investigate the effect of spacing on turbulence structure and mean flow characteristics, both experimentally (Coleman, Nikora, McLean, & Schlicke, 2007;Djenidi, Antonia, Amielh, & Anselmet, 2008;Djenidi, Elavarasan, & Antonia, 1999;Krogstad, Andersson, Bakken, & Ashrafian, 2005;Okamoto, Seo, Nakaso, & Kawai, 1993;Roussinova & Balachandar, 2011) and numerically (Cui, Patel, & Lin, 2003;Ikeda & Durbin, 2007;Stoesser & Nikora, 2008;Stoesser & Rodi, 2004). Simpson (1973), Tani (1987), Jiminez (2004) and Coleman et al. (2007) all proposed that the transition from d-to k-type roughness occurs at around l/k = 5, where l is the crest-to-crest bar spacing and k is the roughness height. Leonardi, Orlandi, Smalley, Djenidi, and Antonia (2003) carried out a series of direct numerical simulations (DNS) of flow over 2D bars aligned transverse to the main flow direction, varying l/k between 1.33 and 20, and found that bars are considered isolated (flow completely reattaches to bed before next roughness structure) when l/k is greater than 8. Stoesser and Nikora (2008) also confirmed this finding in a large-eddy simulation (LES) study, revealing how wall shear stress increases from d-to k-type roughness. Roussinova and Balachandar (2011) investigated the effect of altering the submergence for two different k-type spacings, l/k = 9 and 18, and found that for the larger spacing case the effects of roughness are felt only in the region 0 ≤ z/k ≤ 3 while for the closer spacing case they are felt throughout most of the flow in the outer layer. It is generally accepted for high submergence flows that d-type roughness is characterized by stable separated vortices occupying the entire cavity between roughness elements while k-type roughness entails a mean recirculation bubble in the wake of the roughness elements, with reattachment to the bed occurring between successive elements (Stoesser & Rodi, 2004).
For isolated (k-type) roughness at small submergences the local Froude number may become large enough to produce a hydraulic jump or standing wave at the free surface. These features are characterized by extremely vigorous turbulence production, formation of large-scale turbulent structures, air spray, air entrainment and energy dissipation (Chanson, 2009;Chanson & Brattberg, 2000). Although such features are ubiquitous in civil and environmental engineering scenarios, their contribution to the overall streamwise momentum balance in open channel flows has not been studied with reference to hydraulic resistance.
The aim of the present paper is to shed light on the hydrodynamics of shallow flows that are strongly influenced by water surface deformation, and to quantify the effect of roughness spacing and relative submergence on hydraulic resistance in such flows. Six flow cases, with varying relative submergence and k-type and transitional (between d-and k-type) roughness types, are investigated using LES and complimentary laboratory flume experiments. The double-averaging methodology is applied to the simulated flow to assess qualitatively the contributions to the overall momentum balance of the various additive components such as form-induced and turbulent stresses. The remainder of the paper is organized as follows: Section 2 outlines the experimental set-up, Section 3 presents the numerical methods and computational details, and Section 4 presents and discusses the results. Finally some conclusions are drawn in Section 5.

Experimental set-up
Experiments were carried out in a 10 m long, 30 cm wide glasswalled recirculating flume in the Hyder Hydraulics Laboratory at Cardiff University. A series of plastic square bars 30 cm wide with cross-section 12 mm × 12 mm were installed along the length of the flume, perpendicular to the direction of mean flow (Fig. 1). The roughness height, k, was therefore 12 mm. Two different bar spacings were investigated, l = 62.5 mm and l = 125 mm, corresponding to the normalized spacings l/k = 5.2 and l/k = 10.4 respectively. According to Coleman et al. (2007), the l/k = 5.2 case should be classified as transitional roughness as it is very close to the boundary between d-and k-type roughnesses, while the l/k = 10.4 case constitutes k-type roughness. Bed slope was fixed at 1:50 for all experiments and three flow rates, Q, were tested for each bar spacing (1.7, 2.5, 4.0 l s −1 ), giving a total of six experimental cases. Each case had a different relative submergence, H /k, where H is the double-averaged depth, defined as the distance between the spatial and temporal mean water surface position and the spatial mean bed elevation, z mb . The doubleaveraged bulk velocity, U b (= Q/H ), ranged from 0.20 m s −1 to 0.36 m s −1 . Note that the angular brackets and overbar have not been used for H and U b in the interest of simplicity. Table 1 provides a summary of the flow conditions. In all cases the aspect ratio, B/H, was sufficiently large to ensure either that the flow was two-dimensional or that any secondary flows were weak and that their influence on the mean flow could be neglected.
Measurements of instantaneous velocity and free surface position were taken in a section of the flume where the flow was considered to be uniform and fully-developed. The flow was also considered to be spatially periodic with wavelength l in the streamwise direction, that is to say the temporal mean values of all flow variables in successive cavities between bars were considered to be the same. The bulk Reynolds number (= U b H /ν) was in the range 5700 ≤ R ≤ 13, 000 and the friction Reynolds number R τ (= u * H /ν) where u * is the global friction velocity based on the bed shear stress, was in the range 1.8 × 10 3 ≤ R τ ≤ 3.7 × 10 3 . The global Froude number of the flows, F = (U b / √ gH ), was in the range 0.38 ≤ F ≤ 0.59; local values based on local depths and velocities can be much higher.
A Nixon Flowmeters (Cheltenham, UK) propeller meter was used to measure streamwise velocities along the channel centreline, at two streamwise positions for the l/k = 5.2 case and four streamwise positions for the l/k = 10.4 case. At each streamwise location the velocity was measured at between four and 11 discrete depths, to reveal vertical velocity profiles. At each depth measurements were taken during 120 s at a sampling frequency of 1 Hz; 120 samples of instantaneous velocity were therefore available, from which the temporal mean was calculated.
The water surface elevation was measured using high speed photography with illuminated seeded particles. A Baumer TXG14F CCD camera (Pune, India) was used in conjunction with a Polytec BUS-11 Wotan Flash stroboscope (Waldbronn, Germany) and a halogen lamp to capture the dynamic free surface. Seeding particles that had the same approximate density   Figure 2a presents an example of one such image that was taken for the l/k = 10.4, H /k = 2.9 flow case. Each recording therefore comprised 1500 individual frames, which were then digitally processed using pixel recognition software to reveal the x -z coordinates of the free surface at the channel sidewall at each time instant. These data were then averaged to give a temporal mean of the water surface position. It should be noted that the free surface elevation recorded by this technique was somewhat different from the elevation observed in the centre of the channel due to sidewall effects and the fact that light source illuminated more than just one very thin sheet as would be the case in "proper" particle image velocimetry. Hence the free surface profiles are somewhat smeared and as can be seen from Fig. 2a the water surface appears as a thick or double line rather than a very sharp interface. Thus, in addition to the data that was derived from the photographs, a point gauge was used to measure water surface elevations at several locations along the channel centreline. Measurements were taken over a length spanning two or more cavities, at streamwise intervals of between 2.5 mm and 10 mm. Point gauge measurements were taken for all flow cases except C6, which was characterized by an extremely dynamic free surface with significant spanwise wandering of the hydraulic jump, and it was therefore impossible to measure the water surface with the point gauge.

Numerical framework and computational considerations
The in-house HYDRO3D LES code is used to solve the filtered Navier-Stokes equations for an unsteady, incompressible, viscous flow (Bomminayuni & Stoesser, 2011;Fraga, Stoesser, Lai, & Socolofsky, 2016;Kara, Stoesser, & Sturm, 2012;Kim, Kim, & Stoesser, 2013;Ouro, Harrold, Stoesser, & Bromley, 2017;Stoesser, 2010;Stoesser, McSherry, & Fraga, 2015). LES is an eddy-resolving technique in which the energetic portion of the flow is simulated directly and only the small-scale turbulence is modelled (Stoesser, 2014), and is therefore capable of explicitly predicting turbulence and unsteadiness in flows of engineering importance (Koken & Constantinescu, 2009;McCoy, Constantinescu, & Weber, 2007). The effects of the small-scale turbulence on the large eddies are calculated using the wall-adapting local eddy-viscosity (WALE) subgrid scale model introduced by Nicoud and Ducros (1999). The diffusive terms are approximated by fourth-order central differences while convective fluxes in the momentum and level-set equations are approximated using a fifth-order weighted, essentially non-oscillatory (WENO) scheme. Rodi, Constantinescu, and Stoesser (2013) state that the accuracy and credibility of a code can be justified by the use of high-order spatial discretization schemes together with sufficiently fine grids. A fractional-step method is used with a Runge-Kutta predictor and the multigrid method is used as a corrector of the final step when solving the Poisson pressure-correction equation. The immersed boundary method (Kara, Stoesser, & McSherry, 2015), which maps Eulerian velocities onto Lagrangian point-based representations of non-fluid bodies in the flow, is used to define the geometries of the roughness elements. The position of the free surface is tracked using the level set method (Kang & Sotiropoulos, 2012;Kara, Stoesser, Sturm, & Mulahasan, 2015;Sussman, Smereka, & Osher, 1994), which treats the air-water interface as a sharp interface across which the density and viscosity transitions smoothly. Air entrainment is not taken into account: any air bubbles that appear in the water phase as a result of plunging waves at hydraulic jumps are removed by the diffusivity of the method away from the water surface.
In total 10 LES of flows over transverse square bars have been performed. Computational details are provided in Table 1. Cases C1 to C6 correspond to the six experimental flow cases that are described in the previous section and they were performed on a reasonably coarse grid. Cases F1 to F3 correspond to three selected flow cases and were performed on a fine grid to investigate the sensitivity of the results to grid resolution. Figure 3 presents the computational domains that were used for these first nine simulations; for the l/k = 5.2 cases the domain spanned two cavities while in the l/k = 10.4 cases it spanned only one cavity. For both l/k values the length of the domain, L x , was 10.4k while the width, L y , and height, L z , were both 5k. In the coarse grid the number of uniformly-spaced computational points in the streamwise direction, N x , was 128; in the spanwise direction the number of points, N y , was 64, and the number of points in the vertical direction, N z , was 120. The total number of computational points was therefore approximately 9.8 × 10 5 . The fine grid was discretized with 256 × 128 × 240 (= 7.8 × 10 6 ) computational points. Figure 3 includes isosurfaces showing the instantaneous free surface positions at two arbitrary moments in time, for two example simulations (C2 and C5). The figure shows that the domain extended higher than the free surface: the volume above the surface was occupied by the air phase, and the volume below was occupied by the water phase. A free-slip boundary condition was applied to the top of the domain while a no-slip condition was stipulated on the channel bed. The bars were represented by immersed boundaries, which also achieve an effective no-slip boundary condition on their surfaces (refer to Cevheri, McSherry, & Stoesser, 2016 for a rigorous validation).
Periodic boundary conditions were applied in the streamwise direction, and the flow was driven by the component of gravitational acceleration parallel to the channel bed, based on the bed slope that was applied in the flume experiments (1:50). The global shear velocity, u * (= √ gSH ), was therefore exactly the same as in the experiments. As discussed in Section 2, the aspect ratio of the experimental cases was large enough to ensure that the influence of the vertical side walls on the mean flow would be negligible, and that secondary currents would be negligibly weak. Periodic boundary conditions were therefore applied on the lateral faces of the computational domain to produce a two-dimensional mean flow with no secondary currents.
The domain dimensions, in terms of the mean height of the water surface, H, for each of the cases, are provided in Table 1. A tenth simulation, D1, was performed using a coarse grid in a domain that was twice as long and twice as wide as the original domain, so that the effect of domain size on the computation of turbulent structures could be assessed.
In all cases the simulations were started with a planar rigid lid applied at the mean free surface position that was recorded in the experiments. A free-slip boundary condition was stipulated at the rigid lid and the simulation was run for 100,000 time steps, which corresponded to between 10 and 15 flow through periods, T f (= L x /U b ), to allow the flow to develop fully. The simulation was then restarted without the rigid lid but with the level set algorithm now activated to track the free surface. Averaging of the flow quantities was begun after 4-6 more flow through periods, and continued for between 40 and 60 further flow throughs to ensure that the turbulence statistics were well converged. Further averaging was performed in the homogeneous spanwise direction to obtain a smooth distribution of turbulence statistics. The variables averaged in this way are denoted with an overbar above the corresponding symbols.
The multi-phase flows under investigation herein requires the introduction of the volume fraction, φ, defined as the temporal mean of the fraction of the volume in a plane parallel to the mean flow direction that is occupied by water, i.e. φ = V w /V o , where V w is the temporal mean of the volume of air in the plane and V o is the overall volume of the plane. In the roughness layer the volume fraction is identical to the roughness geometry function defined by Nikora et al. (2001), and takes a value in the range 0.0 < φ < 1.0 (V w < V o ), depending on the roughness 830 R. McSherry et al. Journal of Hydraulic Research Vol. 56, No. 6 (2018) Figure 4 Volume fraction versus height above bed for the six LES cases distribution. In the region immediately above the roughness tops the water occupies the entire volume of the plane (V w = V o ) and therefore φ = 1.0. Due to the spatial heterogeneity of the mean water surface in some of the flow cases (refer to Figs 2 and 3b), a portion of the volume of a plane parallel to the flow direction may be occupied by air and the volume of water in the plane may therefore be less than the volume of the plane. In this region the plane volume may be expressed as where V a is the temporal mean of the volume of air in the plane. It follows that φ < 1.0 when V a = 0.0. Figure 4 presents vertical distributions of the volume fraction for the LES cases C1 to C6. The plots show that, since the roughness elements have a square cross-section, φ is constant in the interfacial sub-layer and undergoes a step change as z increases beyond the roughness tops at z/k = 1. As z increases towards the water surface, two of the cases (Fig. 4b and c) display another abrupt step change as φ drops from 1.0 to 0.0. This is because the mean water surface in these two cases is flat, so the volume of a plane parallel to the mean flow direction is occupied entirely by water on one side of the surface and entirely by air on the other side. In contrast, in the other four cases φ transitions gradually from 1.0 to 0.0 due to the spatial heterogeneity of the mean water surface: the fraction of air in the plane increases with increasing z, until finally no water volume is present.

Results and discussion
Figure 5 plots the variation of the friction factor, presented as (8/f ) 1/2 and calculated from Eq. (1), with relative submergence for the six cases that have been studied. Flume data from Bathurst (1985) are included for comparison. Note that the Bathurst data correspond to gravel bed roughness and were originally presented in terms of H /D 84 , where D 84 corresponds to the 84th percentile used to represent the coarse gravel fraction (i.e. 84% of the gravel elements are smaller than D 84 ). The data have been re-plotted in terms of the equivalent grain roughness, k s , and the relationship k s = 3.5D 84 (Dietrich & Whiting, 1989) has been used to equate the two measures of roughness height. With regards to the data from the present simulations, the equivalent grain roughness is assumed to be equal to the bar height (k s = k); this relationship is likely to be subject to some error but has been used to permit a comparison with the Bathurst flume data. The data lie within the spread of Bathurst's data and agree very well with the general trend of decreasing resistance with increasing relative submergence. Similar to flow over square bars in ducts, for a given submergence, the small bar spacing produces less resistance than the large bar spacing (Leonardi et al., 2003). Figure 2 enables a qualitative comparison of an instantaneous free surface profile recorded by the high speed photography with one computed in the corresponding LES. Although an exact match cannot be expected due to the unsteadiness of the flow and the aforementioned "sidewall smearing" that affected the experimental image (refer to Section 2), the agreement is nonetheless very encouraging. The figure shows the presence of a standing wave, both in the experiment and in the simulation, and that its height and streamwise position are fairly accurately predicted by the LES.
A more quantitative evaluation of the match between LES and experiment can be achieved with the help of Fig. 6, which presents water surface profiles for the six flow cases. Plotted are the spanwise mean of the temporal mean water surface predicted by the LES and the temporal means that were measured using the digitally-processed photographs and data points from the point gauge. Coarse and fine grid LES results are plotted, where available. The agreement between experiments and LES is particularly good for the small spacing (l/k = 5.2) cases, and the accuracy of the free surface computation seems not to be affected by grid resolution (Fig. 6b). The mean free surface in all three small spacing cases is relatively flat, although some undulation is noticeable in the lowest submergence case (Fig. 6a) and this is predicted very well by the LES. In fact two unsteady undular jumps are present in this flow case, slightly upstream of the bars. These jumps establish and disappear (quasi-) periodically. In the large spacing cases well-defined weak hydraulic jumps are visible between the bars, where the local Froude number approaches unity and the flow becomes critical.  For the low submergence, large spacing case (Fig. 6d) the elevation predicted by the coarse grid simulation, C4, is noticeably higher than the measured elevations. This may be the result of numerical diffusion and the imperfect approximation of density and viscosity gradients across the boundary that can occur when the level set method is used to simulate complex free surface topologies on coarse grids. Indeed the fine grid simulation, in which the numerical diffusion is lower and the computation of gradients is more precise, produces a pretty good agreement with the experimental data. The agreement in the medium submergence, large spacing case (Fig. 6e) appears to be very convincing for both coarse and fine grids. It is noteworthy 832 R. McSherry et al. Journal of Hydraulic Research Vol. 56, No. 6 (2018) Figure 7 Measured and computed vertical profiles of temporal and spanwise mean streamwise velocity,ū, for l/k = 5.2 that the coarse grid simulation, C5, appears not to suffer from the same problem as C4, even though the complexity of the free surface profiles is similar for both cases.
As stated in Section 2, point gauge measurements were not possible in the high submergence, large spacing case (Fig. 6f) due to the extremely dynamic nature of the surface and very significant three-dimensionality of the water surface. This dynamic, three-dimensional behaviour is also thought to be the reason for the relatively poor agreement between the LES and the photography data in this case. The periodic boundary conditions that have been applied to the lateral faces of the domain ensure that the mean flow is, by definition, two-dimensional. The spanwise effects that were observed in the experiments due to the presence of the side walls cannot therefore be simulated using this approach. It must therefore be stated that, although the wide channel assumption appears to be valid for the other five flow cases, it may not be valid for C6. Nevertheless some general observations and comparisons can be made from examination of the results. It should also be noted that the spanwise behaviour at the water surface was very dynamic in nature, and may have been the result of an onset of three-dimensional instability at the jump. As no measurements of spanwise and vertical velocity components were possible it is impossible to state with any certainty if significant secondary currents were present, nor indeed whether or not they were related to the three-dimensionality of the water surface.
Further validation of the LES is provided by vertical profiles of the spanwise and temporal mean of the streamwise velocity from the LES in Figs 7 and 8, alongside the experimental data, which were recorded at the channel centreline using the Nixon propeller meter. The measurement locations are illustrated in Fig. 1: two profiles were measured in the l/k = 5.2 cases and four were measured in the l/k = 10.4 cases. In general, the LES data from all six cases show good agreement with the experiments. It should be noted that the experimental data points at z/k = 1 have been given a filled symbol to highlight the fact that they are possibly less reliable than the others: these points are located in the shear layer at the top of the roughness elements, and it is known that the Nixon probe is not well suited to regions of high shear and rotation. In addition some of the discrepancy between experimental and numerical results is likely to be due to spatial resolution: in the experiments this is defined by the propeller diameter (comparable to the bar height) while in the LES the spatial resolution is defined by the grid and is therefore much finer than the resolution in the experiments. In the three l/k = 5.2 cases significant negative spanwise-and temporallyaveraged velocity is observed close to the bed in the middle of the cavity, while there is relatively little streamwise variation in the profiles above the roughness top, indicating a skimming type flow with stable or quasi-stable recirculations in the cavities. In each of the larger spacing cases a region of negative mean velocity is observed near the bed at l/k = 0.25 but not at l/k = 0.5, suggesting that reattachment to the bed occurs in the upstream half of the cavity. The fine grid simulations, F1, F2 and F3, have been included in the plots and the profiles suggest that the mean velocity predictions are not particularly sensitive to the grid resolution.
A more detailed understanding of the flow field can be achieved by considering the contours of spanwise and temporal mean streamwise velocity,ū, that are presented in Fig. 9. For the small spacing, low submergence case (Fig. 9a) the contours reveal that the mean water surface is relatively flat. Non-zero contours of mean velocity are observed above the mean water surface, due to the periodic establishment and disappearance of undular hydraulic jumps just upstream of the bars. When such jumps are present the local water surface rises and a non-zero velocity is observed in the near surface region, contributing to the non-zero mean. Importantly, Fig. 9a also shows significant streamwise variation of velocity throughout the depth, with local acceleration and deceleration occurring between and above the bars, respectively. It is worthwhile to note that computation of these phenomena is permitted by the accurate free-surface tracking algorithm employed by the numerical method; the physics would be incorrectly predicted if a standard free-slip rigid slid assumption was used to represent the water surface. In contrast to the shallow small-spacing case, the two deeper small-spacing cases (Figs 9b and c) are characterized by relatively little streamwise variation of velocity above the bars. The plots for the largespacing cases, Figs 9d-f, reveal regions of locally-accelerated flow that give rise to local Froude numbers close to 1.0, which in turn lead to dramatic deformations of the water surface and the formation of the hydraulic jumps. The hydraulic jumps are characterized by a small region of negative streamwise velocity produced by the plunging and breaking motion of the standing waves. The plots also confirm the presence of a mean recirculation bubble in the wake of the upstream bar, as well as a small corner vortex at the leading face of the downstream bar. Figure 10 presents contours of the spanwise and temporal mean of the streamfunction, ψ, for the six cases. Note that positive streamfunction is denoted by solid contour lines while dashed lines denote negative streamfunction. The plots confirm the observations noted above regarding the velocity profiles and contours, and also offer some important additional insights. Firstly, the plots clearly show that the close spacing cases are characterized by large recirculations in the cavities, while the mean flow does reattach to the bed in the large spacing cases. The contours also confirm that the l/k = 5.2 cases should indeed be classified as transitionally rough, rather than d-type, as the streamlines in the shear layer do not quite connect the top of the elements to create a pseudo-smooth wall, as in Stoesser and Rodi (2004) for example. In fact the ψ = 0 contour impinges on the leading faces of the roughness elements, indicating that some form drag will be experienced at the top of the bars and flow in this region will undergo some degree of unsteadiness. Therefore, although the velocity profiles in Fig. 7 suggest relatively little streamwise variation in the streamwise velocity, the mean flow field above the roughness crests is not quite homogeneous, as would be the case in a true skimming or pseudo-smooth wall flow. Instead, it can be considered that the bars present a pseudo-rough wall at the crest height, resulting in the streamwise undulations in the streamfunction contours that are visible in the plots. The lowest submergence case (Fig. 10a) in particular reveals a somewhat wavy streamfunction field and suggests the presence of an undular hydraulic jump at the free surface, as discussed above with reference to the contours of mean streamwise velocity. In the l/k = 10.4 cases the mean flow reattaches to the bed between the bars and can be unequivocally categorized k-type roughness. The length of the mean recirculation bubble decreases as submergence increases: this will be discussed in more detail below. Figure 11 illustrates how the wall-normal distance of the first computational point for streamwise velocity, z +l (= zu l * /ν), varies along the length of the computational domain. Here u l * is the local friction velocity, calculated using the local shear stress, τ l (= μ∂ū/∂z), and therefore does not account for the form drag acting on the bars. The global friction velocity, u * (= √ gSH ), does account for the form drag and is expected to be higher than u * l . Note that, because the LES are performed on a staggered grid in which the velocities are assigned to the centre of the cell faces, the normalized distance between the wall and the closest u grid point to it is half the cell height, i.e. z + /2. Based on the cell dimensions given in Table 1, this means that the distance from the wall to the closest u grid point is in the range 18 ≥ z + ≥ 23 for the coarse grids and 9 ≥ z + ≥ 11 for the fine grids. The plots show, as expected, that the values of z +l are much lower than the corresponding z + values. The plots also provide confirmation that at least one computational point is well inside the viscous sub-layer along the length of the domain, indicating that wall-resolved LES have been achieved in every case. Peaks in z +l coincide with peaks in wall shear stress, and for both bar spacings the magnitude of the peak increases with increasing relative submergence. In the small spacing cases one significant peak in z +l is observed on the channel bed, the result of the large quasi-steady vortex in the cavity which produces highly negative velocities in the near-wall region. In the large spacing cases a significant peak is observed in the region that is subject to negative near-wall velocities due to the presence of the recirculation in the wake. It is noteworthy that the peak values in the recirculation region are higher than those in the region 4.3 ≤ x/k ≤ 9.0, which corresponds to the portion of the bed in which the flow has reattached to the wall. The minima that are visible approximately half way between the bars correspond to the position at which wall shear stress switches from negative to positive, and give an alternative measure of the location of the reattachment point that was discussed above with reference to the streamfunction. The variation of z +l suggests that in the lowest submergence case the flow reattaches at x = 5.3k and in the highest case it reattaches at x = 4.4k, a decrease of approximately 17%. Figure 12 presents contours of the spanwise and temporal mean of the turbulent kinetic energy, tke, normalized on the square of the global friction velocity u * , for the six flow cases. The mean free surface position is overlaid for reference. A highly turbulent shear layer is present above the roughness crests in the small spacing cases, and in the two higher submergence flows for that spacing (Figs 12b and c) there appears to be relatively little turbulent activity near the water surface. The low submergence, small spacing case (Fig. 12a) displays a different character, however, showing high levels of turbulent energy near the surface and a noticeably weaker turbulent shear layer near the roughness crests compared to the higher submergence cases. The turbulent energy observed at the free surface is produced by the unsteady undular jumps that are indicated in the corresponding streamfunction plot (Fig. 10a) and involve the generation of large scale turbulent structures that are convected downstream with the mean flow. The complexity of the flow in the large spacing cases is evident from the tke plots. Each case displays a very localized peak in the turbulent kinetic energy at the location of the hydraulic jumps that occur between the bars.
All cases under investigation were expected to belong to the class of intermediate submergence (type II) flows as defined by Nikora et al. (2007a). The boundary layer was therefore not expected to comprise the logarithmic layer that occurs in smooth bed flows at high Reynolds numbers and in rough bed flows at 836 R. McSherry et al. Journal of Hydraulic Research Vol. 56, No. 6 (2018) Figure 13 Variation of double-averaged streamwise velocity with wall-normal distance. Dashed lines represent fitted log law lines high submergence. Figure 13 presents profiles of the doubleaveraged streamwise velocity, u + (= ū /u * ), where u * is the global shear velocity, the overbar denotes temporal averaging and the triangular brackets denote streamwise and spanwise spatial averaging. In all cases spatial averaging is performed over the entire computational area. The origin of the vertical axis and the quantification of the zero-plane displacement, d, are defined following the Clauser method (Wei, Schmidt, & McMurty, 2005). Cases with the same flow rate (Table 1) and therefore similar relative submergence are plotted together to enable meaningful comparisons. The plots show that the streamwise velocity is generally higher in the small spacing cases than in the large spacing cases, over the entire flow depth.
To analyse the vertical distribution of the double-averaged streamwise velocity in greater detail it is useful to consider the universal formulation for the logarithmic layer in turbulent smooth bed flows: where κ is the von Kármán constant, typically taken to be between 0.40 and 0.42 and A is a constant typically taken to be 5.2 (Bailey, Vallikivi, Hultmark, & Smits, 2014). For rough bed flows the double-averaged velocity distribution may be used, and Eq. (2) becomes: where and u + is the downward shift of the profile due to bed roughness, relative to the corresponding smooth bed case. Figure 13 includes fitted log-law lines using Eq. (3), with κ = 0.41, A = 5.2 and u + determined using a best-fit procedure. The plots show that some of the profiles do comprise relatively straight sections that are fairly well aligned to the loglaw lines, and it could therefore be argued that a logarithmic layer is present to some extent in these cases. In general this straight section occupies a larger portion of the water depth in the small spacing cases (Fig. 13a) than in the large spacing cases (Fig. 13b), suggesting that the logarithmic region, if it exists at all, is more established in the small spacing cases.
(3) and rearranging after making substitutions for u + and z + , it is possible to show that: In a logarithmic layer κ, and therefore (d ū /dz)(z/u * ), would be constant, i.e. would not vary with distance from the bed. Figure 14 plots (d ū /dz)(z/u * ) as a function of z/k and reveals that for none of the six cases could it be definitively stated that this quantity remains constant over a significant portion of the flow depth. It is therefore not possible to conclude that any of the flows comprise a logarithmic layer. Nonetheless, the profiles for the highest submergence cases do comprise sections (for z/k 1.4) in which (d ū /dz)(z/u * ) varies relatively little and may therefore be considered quasi-constant. It is interesting to note that these quasi-constant values are different for the two spacings: approximately 4 for l/k = 5.2 and approximately 2.4 for l/k = 10.4. These values correspond to κ = 0.25 and κ = 0.41 respectively (Eq. (4)) and suggest that for low submergence flows of this type the constant of proportionality in the logarithmic layer, if indeed it is present at all, may differ markedly from the values that have previously been observed for larger submergence flows. Figure 15 presents plots of the double-averaged streamwise velocity distribution in the near-bed region. Figure 15a presents the small spacing cases while Fig. 15b does the same for the large spacing cases. The plots reveal that, when the normalized velocity u + is used, the profiles for each spacing are very similar, particularly in the interfacial sub-layer (0 ≤ z ≤ k). As such, it is possible to assign a functional form to the distributions in this region. As can be seen in the figure, the distribution appears to follow an exponential form in the small spacing case whereas the distribution may be described as linear in the large spacing cases. This is consistent with the findings of Coleman et al. (2007). Above the roughness tops the flow profiles begin to deviate from the functional forms that define them in the interfacial region, and for neither spacing could the velocity profile Vertical profiles of normalized spatially-averaged Reynolds shear stress, u w , are shown in Fig. 16. The plots reveal markedly different trends for the small and large spacing cases. For the two small spacing, high submergence cases (Figs 16b and c), u w increases approximately linearly with distance from the free surface, reaching a peak at the height of the roughness crests. This is consistent with other studies of turbulence in rough-bed flows (Bomminayuni & Stoesser, 2011;Manes et al., 2007;Singh, Sandham, & Williams, 2007;Stoesser & Nikora, 2008). Below the roughness crests u w decreases approximately linearly from the peak until it reaches zero at the channel bed, an observation also made by Stoesser and Nikora (2008). In the corresponding large spacing cases the variation of the Reynolds stress with height is clearly not 838 R. McSherry et al. Journal of Hydraulic Research Vol. 56, No. 6 (2018) Figure 17 Vertical profiles of spatially-averaged form-induced primary shear stress linear, either above the crests or below, and displays significant curvature throughout the depth.
An interesting feature of the Reynolds stress variation in all three large spacing cases, and also in the small spacing, low submergence case, is the portion of significantly positive u w that occurs close to the water surface. In all cases where positive near-surface u w is observed, the depth at which the u w line changes sign corresponds to the depth at which maximum double-averaged streamwise velocity is observed. These positive stresses are the result of turbulence production at the unsteady surface, which has already been discussed with reference to the tke contours (Fig. 12). This finding has implications for the streamwise momentum balance and overall hydraulic resistance in the channel. Nikora (2009) showed that negative spatially-averaged Reynolds stress adds to the friction, while positive stress acts to reduce it. These results therefore suggest that in open channel flows with unsteady, turbulent free surfaces, turbulent stresses near the surface produce negative streamwise drag. Figure 17 presents vertical profiles of the normalized spatially-averaged form-induced shear stress, ũw , which arises from possible correlations in spatial fluctuations in the mean flow field due to the spatially heterogeneous nature of the rough bed. As with the Reynolds shear stress, a negative spatiallyaveraged form-induced stress will contribute positively to the overall hydraulic resistance in the channel. The plots show that in all cases ũw increases nonlinearly with distance from the bed until it reaches a peak at the height of the roughness crests. This is consistent with previous studies by Manes et al. (2007) and Stoesser and Nikora (2008), for example, who also showed that ũw then decreased fairly rapidly above the roughness crests and remained at zero for the upper part of the water column. Fig. 17b and c confirm that this general trend is also observed in the present study for the two small spacing, high submergence cases. This suggests that these two cases should be classified as type I (large submergence) or type II (intermediate submergence) flows, since there clearly exists a portion of the depth in which form-induced stresses are not significant.
In the other four cases, however, significant non-zero forminduced stresses are observed throughout the whole water depth, indicating that the form-induced sub-layer extends all the way to the free surface. Furthermore, the most striking features of Fig. 17 are the very large negative stresses in the near-surface region for the four cases in which the surface was observed to be turbulent and unsteady. This result indicates that the standing wave-type responses that were observed at the free surface will make a significant contribution to the overall resistance in the flow.
It is noteworthy that the magnitudes of these negative nearsurface peaks in form-induced stress are always larger than the corresponding positive Reynolds stress peaks that are observed in the same region (Fig. 16): the net contribution of the free surface response to the overall drag is therefore always positive, as confirmed by Fig. 18, which presents profiles of the combined (form-induced + Reynolds) stress. The stress variation above the roughness crests is shown to be approximately linear, even in the large spacing cases, except in the near-surface region. The plots suggest that extrapolation of the slope from the height of the roughness crests to the x-axis would result in a stress of approximately −1, which is consistent with other studies of flows over rough beds, for example in Manes et al. (2007). The reason that the stress does not reach this value at the bed is the dominance of the form drag in the interfacial sub-layer (0 ≤ z ≤ k), which produces an approximately linear reduction in the contribution of Reynolds and form-induced stresses until they disappear completely at the bed.
Also included in Figs 16-18 are results from the fine grid and large domain simulations, where available. The plots show that the agreement is largely good and that, at least in terms of the identification of trends and qualitative analysis, the coarse grid simulations on the smaller domains produce adequate results. One notable exception is the small spacing, low submergence case, which was shown to produce inaccurate predictions for the mean surface elevation (Fig. 6d). As would be expected given the fact that the surface elevations were over-predicted in the coarse grid simulation, the near-surface stress peaks for this case are noticeably lower in the fine grid simulation, but the characteristics is consistent (Figs 16a and 17a). Figure 19 presents the variation with distance from the bed of all terms that appear in the double-averaged streamwise momentum balance for steady two-dimensional uniform open channel flows: Reynolds stress, form-induced stress, skin friction and pressure drag. The total line corresponds to the sum of all of these terms, while the gravity line is also included for reference. For convenience only two cases, one for each bar spacing (C2 and C5), have been included. It is noteworthy that for both spacings the total momentum flux follows the gravity line rather well over the majority of the flow depth above the roughness tops. Also evident in both cases is the significant contribution to the momentum balance of the form drag, which dominates in the region below the roughness tops and reduces to zero above the roughness. The character of the form drag distribution differs slightly between the cases: for small bar spacing the distribution follows a slightly concave shape whereas for large spacing it is characterized by an almost linear slope. This reflects the fact that in the small-spacing case there is significant recirculation inside the cavities, possibly causing a nonlinear distribution of the form drag. Figure 19b, the momentum balance for the large spacing case, shows that form-induced stress at the undular jump dominates in the near-surface region and results in an overall momentum flux that exceeds the gravity flux. In both cases steps in the Reynolds stress and form-induced stress profiles are evident at the location of the roughness tops: this is because in the momentum balance these terms are multiplied by the roughness geometry function φ which, as shown in Fig. 4, is characterized by a step change as z increases beyond the roughness tops. In both cases the sums of the four selected contributors to the momentum balance do not fall perfectly on the gravity line because various other minor contributions are not included. These are subgrid-scale stresses, flow non-uniformity in the streamwise direction, normal stresses and contributions from secondary currents (due to not perfectly converged statistics of the LES), however their individual magnitudes are all very small.
The total momentum flux is plotted in Fig. 20, using the roughness height k to normalize the vertical coordinate. This normalization permits meaningful comparison of the data with two cases from the experimental study of Coleman et al. (2007), which were characterized by much higher relative submergence (H /k = 11) than was the case for the present study. It should also be noted that the bar spacings from the Coleman et al. (2007)  different trend in the interfacial region: a smaller bar spacing (l/k = 8) produced a total flux that never surpassed the gravity flux, whereas in the large spacing case it surpassed the theoretical gravity flux by a considerable margin. In Coleman et al. (2007) no explanation for this result was proposed, and it is assumed to be due to a combination of limitations of experimental methods, measurement accuracy and unaccounted-for non-uniform flow effects. As mentioned above in the present study the total flux for all of the LES cases is slightly lower than the theoretical gravity flux over most of the interfacial sub-layer, echoing to some extent the behaviour of the l/k = 8 of Coleman et al. (2007). The Coleman et al. (2007) data are also characterized by rather pronounced concave-upward curves above the roughness tops, which is due to the fact that the flows were accelerating. A further difference is noted in Fig. 20b near the surface, where neither experimental case exhibits any significant flux. This is not surprising given the high submergence of the experimental tests, and it is assumed that the water surface was flat and unaffected by the bed topology.

Conclusions
Results from large-eddy simulations and complementary flume experiments of turbulent open channel flows over bed-mounted square bars at intermediate submergence have been presented. In total six flow cases were investigated, comprising two roughness spacings that correspond to transitional and k-type roughness and three flow rates. The bed slope was held constant for all cases, and relative submergence therefore increased with flow rate. In the experiments the water surface was observed to be very complex and turbulent for the large spacing cases, and comprised a single hydraulic jump between the bars. The streamwise position of the jump varied between the cases, with the distance of the jump from the previous upstream bar increasing with flow rate. The free surface was observed to be less complex in the small spacing cases, particularly for the two higher flow rates, in which the flow resembled a classic skimming flow. The Darcy-Weisbach friction factor was calculated for all six cases from a simple momentum balance, and it has been shown that for a given flow rate the larger bar spacing produces higher resistance.
The predictions of the LES have been shown to be in reasonably good agreement with the experiments in terms of mean free surface position and mean streamwise velocity. The position of the hydraulic jumps have been well represented. Contours of spanwise and temporal mean streamlines revealed that the small spacing cases are characterized by classic cavity flow with a quasi-rough wall presented at the height of the roughness crests. In the large spacing cases the length of the mean recirculation bubble was observed to decrease with relative submergence. Analysis of the vertical distribution of streamwise velocity has revealed that in the lowest submergence cases no logarithmic layer is present, whereas in the higher submergence cases some evidence of a logarithmic layer is observed.
Analysis of spatially-averaged Reynolds shear and forminduced stress profiles revealed that, in cases in which significant turbulent activity occurs at the free surface, the Reynolds stress may be positive near the surface and can therefore reduce the overall drag experienced by the flow. Secondly, vertical profiles of form-induced stress revealed significant negative peaks in these complex near-surface regions, resulting from spatial fluctuations in the mean flow due to the presence of hydraulic jumps. The magnitudes of these negative form-induced stresses are larger than those of the positive Reynolds stresses, and therefore act to produce a net contribution to the overall drag at the surface.

Funding
This work was supported by the UK Engineering and Physical Sciences Research Council (EPSRC; grant number EP/k041088/1).

Notation f
= Darcy-Weisbach friction factor g = acceleration due to gravity k = roughness height tke = turbulent kinetic energy u * = global friction velocity u = double-averaged streamwise velocity ũw = dispersive shear stress u w = turbulent shear stress z = bed elevation z mb = spatial mean bed elevation A = constant B = channel width D 84 = the 84 th percentile used to represent the coarse gravel fraction H = mean water depth L x = channel length L y = channel width L z = channel height N x = # of grid points in streamwise direction N y = # of grid points in spanwise direction N z = # of grid points in the vertical direction Q = flow rate U b = bulk flow velocity S = bed slope T f = flow through periods V a = temporal mean of the volume of water in a plane V o = overall volume of a plane V w = temporal mean of the volume of air in a plane F = bulk Froude number R = bulk Reynolds number R τ = bulk friction Reynolds number φ = volume fraction κ = von Kármán constant l = crest-to-crest bar spacing ψ = streamfunction θ = an instantaneous variable τ l = local shear stress ORCID T. Stoesser http://orcid.org/0000-0001-8874-9793