A numerical investigation into the importance of bed permeability on determining flow structures over river dunes

Although permeable sediments dominate the majority of natural environments past work concerning bed form dynamics has considered the bed to be impermeable, and has generally neglected flow between the hyporheic zone and boundary layer. Herein, we present results detailing numerically modeled flow which allow the effects of bed permeability on bed form dynamics to be assessed. Simulation of an isolated impermeable bed form over a permeable bed shows that flow is forced into the bed upstream of the dune and returns to the boundary layer at the leeside, in the form of returning jets that generate horseshoe‐shaped vortices. The returning flow significantly influences the leeside flow, modifying the separation zone, lifting the shear layer adjoining the separation zone away from the bed. Simulation of a permeable dune on a permeable bed reveals even greater modifications as the flow through the dune negates the formation of any flow separation in the leeside. With two dunes placed in series the flow over the downstream dune is influenced by the developing boundary layer on the leeside of the upstream dune. For the permeable bed case, the upwelling flow lifts the separated flow from the bed, modifies the shear layer through the coalescence with vortices generated, and causes the shear layer to undulate rather than be parallel to the bed. These results demonstrate the significant effect that bed permeability has on the flow over bed forms that may be critical in affecting the flux of water and nutrients.


Introduction
Dunes are widespread in alluvial channels and generated from sediment that range in size from fine sands to gravels [Dinehart, 1992;Seminara, 1995;Best, 1996;Carling, 1999;Kleinhans, 2001Kleinhans, , 2002Carling et al., 2005;Best, 2005;Bradley et al., 2013]. Their presence significantly influences the mean and turbulent flow, and therefore sediment dynamics [Best, 2005;van der Mark et al., 2008;Naqshband et al., 2014], as well as flow resistance and energy losses within open channels [Garc ıa, 2008;Lefebvre et al., 2014]. However, much of our current understanding of the flow dynamics of these bed forms stems from the simplest scenariothat of an impermeable two-dimensional dune that is asymmetric in cross section, with an angle-of-repose leeside, and formed under equilibrium conditions in uniform, unidirectional flow. These flow conditions are described by Best [2005] who proposed a schematic model from a synopsis of previous work [e.g., McLean andSmith, 1979, 1986;Nelson and Smith, 1989;McLean, 1990;Nelson et al., 1993;McLean et al., 1994McLean et al., , 1999aMcLean et al., , 1999bBennett and Best, 1995;Maddux et al., 2003a;Bridge, 2003;Kleinhans, 2004] that consisted of five distinct flow zones: (i) accelerating flow over the stoss side; (ii) flow separation or deceleration [Nelson et al., 1993;McLean et al., 1994;Best and Kostaschuck, 2002] from the crest in the leeside; (iii) flow reattachment at 4-6 dune heights downstream; (iv) a shear layer between the separated flow and streamwise flow, which expands downstream, and (v) an internal boundary layer that grows from the reattachment point and along the stoss slope of the next downstream dune. These flow dynamics have been inferred to strongly influence the morphodynamics of the river bed. For example, turbulence formed by leeside flow separation, by Kelvin-Helmholtz instabilities produced along the bounding shear layer [Best, 2005;Abad, 2008;Frias and Abad, 2013], generates large instantaneous Reynolds stresses [van Mierlo and de Ruiter, 1988;Nelson et al., 1993Nelson et al., , 1995Bennett and Best, 1995;McLean et al., 1999a]. These high shear stresses are significant in determining the suspension of sediment [Jackson, 1976;Nelson et al., 1995;Best, 2005aBest, , 2005b, local bed load transport rates, the formation of dunes from ripples [Bennett and Best, 1995], and therefore influence the overall dune morphology [Bennett and Best, 1995]. However, this conceptual model assumes that the river bed, including the bed forms, is impermeable. Yet, porous bed forms migrating over permeable beds (e.g., gravel/sand dunes) are abundant in fluvial environments, and as such develop hyporheic flow-a flux of water between the boundary layer flow and the river bed [Harvey et al., 1996;Winter et al., 1998;Storey et al., 2003]. Such hyporheic flows are driven by hydrostatic and hydrodynamic forces [Boano et al., 2014]. Hydrostatic forces are induced by pressure determined by the depth of the overlying water, and typically influences hyporheic flows over spatial scales of hundreds of meters. In contrast, hydrodynamic forcing tends to be most important in driving shallow localized hyporheic flows, where such flows are formed by local velocity and pressure gradients that are generated by flow around, over, or through bed forms and that result in momentum transfer into the bed [Boano et al., 2014]. The flow typically enters the bed in a region of high velocity upstream of the bed form and then reemerges downstream of the crest within the recirculation zone [Thibodeaux and Boyle, 1987;Savant et al., 1987;Blois et al., 2014]. These hydrodynamic hyporheic flows scale positively with mean boundary layer velocity and mean bed permeability Cardenas et al., 2004;Sawyer and Cardenas, 2009], and negatively with bed form wavelength [Boano et al., 2014]. The wavelength of the bed form also directly influences the depth and average residence times of such hydrodynamically induced flows. The depth of hyporheic flow typically scales on the order of 80% of the bed form wavelength [Elliott and Brooks, 1997;Cardenas and Wilson, 2007a], although both the depth of hyporheic flow [Stonedahl et al., 2013;K€ aser et al., 2013] and the decay rate of the water flux into the bed [Boano et al., 2014] scale positively with the depth to the free surface. Finally, the average residence time of fluid within the bed scales positively with the bed form wavelength and negatively with streambed permeability [Elliott and Brooks, 1997;Stonedahl et al., 2010;Gualtieri, 2012Gualtieri, , 2014.
Hydrodynamically induced flow can also be generated by boundary layer turbulence in open channel flow [e.g., Blois et al., 2012Blois et al., , 2014, especially in gravels where pore water velocities near the bed interface can reach magnitudes in the order of 0.2 m s 21 [Nagaoka and Ohgaki, 1990]. These flow velocities are generated by coherent flow structures in the turbulent boundary layer, as well as flow separation induced by bedtopography, and may be expected to generate significant pressure gradients across the bed and into the pores beneath [Detert et al., 2004]. These flows result in intense turbulent fluctuations in the hyporheic flow within the pore spaces in the Brinkman layer (defined as the transition layer between the fully turbulent flow and the deeper groundwater Darcian flow, Goharzadeh et al. [2006]). As such, the role of bed forms in both inducing and driving hyporheic flow has been recognized previously [Stonedahl et al., 2012[Stonedahl et al., , 2013 and investigated through the use of numerical models [Boano et al., 2007;Cardenas andWilson, 2007a, 2007b;Packman et al., 2004;Hester et al., 2013].
In these previous simulations neither turbulence in the Brinkman layer nor the advective (hydrostatic) pumping is correctly predicted from the pressure field generated in the boundary layer. However, it is clear that this is too simple a representation of flow in many coarse porous beds, as recent work has shown the recirculation zone and downstream reattachment point reported in the standard model of flow over dunes is replaced by strong vertical upwelling from the porous bed [Blois et al., 2014]. An understanding of these flows is essential to understand how the boundary layer hydraulics are altered, and how these flows may influence the eco-hydraulics of the river bed through: (i) solute exchange that are dependent on hyporheic exchange and residence time in the subsurface [Zarnetske et al., 2011;Bardini et al., 2012;Marzdari et al., 2012;Arnon et al., 2013;Gomez-Velez et al., 2015], and (ii) the structure of the bed, since flows across the boundary can cause dilation of framework gravels [Allan and Frostick, 1999], thus acting as a mechanism that can introduce fine sediment into the river bed, and thus potentially degrading its quality as a habitat for species such as salmonids [Greig et al., 2007].
The present paper develops and applies a novel Computational Fluid Dynamics (CFD) model that is able to simulate flow both in the boundary layer flow and Brinkman layer. This provides a new methodological approach to predict simultaneously surface and subsurface flow in one numerical scheme. As the same discretization of the Navier-Stokes equations is applied throughout the solution domain, the pressure term, as well as the velocity and turbulence terms, is maintained without requiring any coupling. This allows the pressure-driven advective terms, discussed above, to be correctly predicted. To assess and validate this methodology, the physical experiments of Blois et al. [2014] are utilised. The analysis is then extended to undertake a numerical sensitivity experiment to assess how bed permeability influences flow over bed forms (Table 1). This enables a reassessment of the influence of bed and bed form permeability on the Water Resources Research 10.1002/2016WR019662 generation of flow structures and their pathways within the bed and boundary layer. The coupled approach presented herein thus allows both hydrostatically and hydrodynamically induced hyporheic flow to be modeled simultaneously, providing insight into the nature of turbulence-driven hyporheic flow.

Numerical Model Description
The hydraulics of open-channel flow are governed by the Navier-Stokes equations, a set of conservation laws for mass and momentum. Here we apply a Reynolds-averaged Navier-Stokes (RANS) model to predict the time-averaged flow field. This set of equations written in Cartesian tensor form is: for i 5 1, 2, 3 (x, y, z), where u is the time-averaged velocity, u 0 is the turbulent velocity fluctuation (relative to the time-averaged value), p is pressure, q is the fluid density, m m is the fluid molecular kinematic viscosity, g i is the gravitational acceleration constant (g 1 5 g 2 5 0 and g 3 5 g), and u 0 i u 0 jj is the Reynolds stress (a product of time-averaging the Navier-Stokes equations).
To close the governing equations, we apply a two-equation renormalization group theory (RNG) j-e turbulence model [Yakhot et al., 1992]. This methodology has previously been shown to outperform the standard j-e turbulence model in areas of flow separation and reattachment [Yakhot et al., 1992;Lien and Leschziner, 1994a;Dargahi, 2004].
The numerical model is a three-dimensional finite volume model. The pressure and momentum equations are coupled using SIMPLEST, a variation on the SIMPLE algorithm of Pantankar and Spalding [1972]. Finally, the convective terms were computed using the second-order upstream monotonic interpolation for scalar transport (UMIST) differencing scheme [Lien and Leschziner, 1994b].

Experimental Methods
Five numerical experiments (Table 1) were conducted based on the experimental data of Blois et al. [2014]. Two initial experimental configurations were studied: (i) an impermeable dune over a smooth impermeable bed, and (ii) an impermeable dune over a rough permeable bed, as these represent the conditions reported by Blois et al. [2014, see below] and provided boundary conditions and validation data for the numerical model. In the third experiment, we modeled a permeable dune over a permeable rough bed, as also reported by Blois et al. [2014], but herein we modeled the permeable dune as comprising regular spheres of 0.01 m diameter. Blois et al. [2014] used natural heterogeneous gravels of the same size but did not record the packing density or the morphology. Finally, the first two experiments were repeated but with the presence of two bed forms.
In all the present numerical experiments, we employed the same channel dimensions as Blois et al. [2014]; 4.8 m in length (L), 0.35 m in width (W) with a height (H) of 0.60 m. In experiments 2, 4, and 5 that simulated a permeable bed (Table 1), we modeled a permeable bed covering the entire flume and that comprised uniform spheres of diameter (D) 0.038 m packed in rectangular fashion, identical to the physical experiments of Blois et al. [2014]. The permeable section of the domain in the numerical model comprised six layers of uniform spheres, producing a depth of permeable section (h bed ) of 0.228 m giving a bed porosity of 50%. A flow depth (h w ) of 0.19 m above the bed was maintained, and hence the ratio of flow depth to permeable bed thickness was 0.8. For experiments 1, 2, 4, and 5, an impermeable dune was modeled that was 0.41 m in wavelength (k), 0.056 m in amplitude and with a leeside angle a lee 5 278. The dune spanned the entire channel width. This geometry was replicated as closely as possible with 0.01 m spheres for the permeable dune in experiment 3. All topography-the bed form and porous bed-was represented in the numerical model using a mass-flux scaling algorithm [Hardy et al., 2005[Hardy et al., , 2007 to represent the dune and bed morphology in a structured grid solution, without the need for boundary-fitting grids [Lane et al., 2004].

Numerical Model Boundary Conditions
A streamwise freestream velocity (U o ) of 0.16 m s 21 was applied in all the numerical simulations, with the inlet flow velocity (U o ) and flow depth (h w ) yielding a freestream flow Reynolds number of 3.0 3 10 4 and a Froude number of 0.12. The outlet defined at the downstream end of the domain used a fixed-pressure boundary condition, although mass was allowed to enter and leave the domain. The walls of the channel were specified as nonslip conditions. At the free surface a rigid-lid approximation was applied, through the introduction of an additional surface pressure term, but requiring a mass correction of the adjacent cells [Bradbrook et al., 2000].
The resolution of the computational mesh (Table 1) was identified after conducting a grid sensitivity test [e.g., Hardy et al., 2003] to ensure that the model was both mesh independent and able to capture the pertinent flow features. The convergence criterion was set such that mass and momentum flux residuals were reduced to 0.1% of the inlet flux. The mesh dimensions used for the different simulations are given in Table 1.

Analysis of Flow
Results are presented and discussed in their dimensionless form, where downstream (U xt ) and vertical (W xt ) velocities are nondimensionalized with respect to the inlet velocity (U o ), turbulent kinetic energy (TKE xt ) is nondimensionlized with respect to the square of U o , and spatial units (x/d in downstream or y/d in vertical) are nondimensionalized with respect to the diameter of the spheres (D 5 0.038 m) used to represent the porous bed. For the single dune case, profiles of streamwise velocity and turbulence intensity were extracted at the dune crest and then every 2 x/d downstream, while for the experiments with two dunes, profiles were extracted at both dune crests and one 2 x/d downstream from the crests. In the vertical axis, the top of the spheres comprising the channel bed is always set as zero height, and thus in experiments 2, 4, and 5 where the Brinkman layer flow is investigated, the flow depth also possesses a negative height.
Once the Reynolds-averaged flow fields and velocity profiles have been reported, the analysis is extended to consider flow pathways, the generation of turbulent structures and the prediction of near-bed shear stresses. Flow pathways are calculated through an array of streamlines, coloured by W xt , and demonstrate potential flow paths through the domain. Turbulent structures are visualised applying the Q criterion [Hunt et al., 1988] that identifies a vortex if the magnitude of the vorticity tensor is greater than that of the rate of strain tensor and there exists a localized pressure minimum. Finally, near-bed shear stress (s) is predicted (s 5 0.19 q TKE, Biron et al. [2004]) to assess how the modification in local flow conditions may influence potential morphodynamic evolution of the bed.
We first report results for the two cases of an asymmetric dune placed on an impermeable and then a permeable bed, and validate the model using the data of Blois et al. [2014] by comparing experimental and modeled at-a-point velocity. In order to minimize any error introduced through the geolocation and matching of exact points between the two set ups, we have taken the five nearest points within a 0.05 m search radius of the measured points and averaged the modeled data weighted by the inverse of the squared distance (the interpolation method of Shepard [1968]). Second, the spatial patterns in the modeled velocity fields are compared with key hydraulic observations (Table 2). Once the two validation cases have been presented, we then extend the analysis to consider cases of a permeable dune upon a permeable bed and two dunes.

Validation Against at-a-Point Experimental Data
Validation of the scheme is undertaken when the experimental and modeled at-a-point downstream (u2) velocity data are compared (Figures 1a and 1b) for: (i) an impermeable dune on an impermeable bed, and (ii) an impermeable dune on a permeable bed. In both cases, there is a good overall agreement with a coefficient of determination of 0.89 (correlation coefficient 5 0.94) for case 1 and a coefficient of determination of 0.55 (correlation coefficient 5 0.74) for case 2. For both cases, the level of agreement is statistically significant (n 1200, r > 0.058). This is extended to consider the prediction of flow within the bed (Figure 1b, red dots). Again good overall agreement is demonstrated when experimental and modeled at-a-point downstream (u2) velocity data and turbulent kinetic energy are compared (correlation coefficients of 0.6 and 0.65) considering possible geolocation errors and the strong spatial gradients in shear flow within the pore volume. Furthermore, this shows a similar level of agreement for this numerical scheme to that which has accurately predicted only boundary layer flow in a series of previous geomorphological applications [e.g., Bradbrook et al., 1998;Ferguson et al., 2003;Lane et al., 2004;Hardy et al., 2013]. This thus demonstrates confidence in the scheme to explore boundary layer-hyporheic flow interactions.

Flow Fields Over a Single Impermeable Dune on an Impermeable Bed
The U xt interpolated experimental data and modeled data (Figures 1c and 1e) both show the known flow characteristics associated with flow over a single two-dimensional dune and can be interpreted using the schematic model of Best [2005]. In both cases, accelerating flow over the stoss side of the dune is observed, with flow separation or deceleration from the crest in the dune leeside and flow reattachment at 4-6 x/d downstream. In the modeled simulation, the recirculation zone behind the dune is stronger with contours more closely spaced, although this may be a function of the interpolation of the experimental data. Furthermore, the recirculation zone extends further downstream with flow reattachment at 9 x/d downstream of the dune, which is equivalent to a scale of 6 bed form heights, consistent with previous observations [Best, 2005;Venditti et al., 2013]. Finally, in both measured and modeled cases, a shear layer exists between the separated flow and streamwise flow, and which expands downstream. This is further identified in W xt (Figure 1e) with positive velocity on the stoss side and negative flow above the flow separation in the leeside. Finally, a zone of maximum TKE xt is contained within the expanding shear layer in the downstream direction ( Figure 1i) and is formed due to the high vertical gradient of U xt generated across the recirculation zone. The similarities in both the magnitude and spatial distribution of the flow structures confirm that the numerical scheme is accurately predicting the flow in agreement with the schematic model of Best [2005] and the experimental validation data given in Blois et al. [2014].

Flow Fields Over a Single Impermeable Dune on a Permeable Bed
A comparison between the U xt interpolated experimental data and modeled data (Figures 1d and 1f) provides a second validation test of the numerical scheme, and a comparison of the main flow features is presented in Table 2. In both the experimental and modeled data, the flow downstream of the crest is modified with a decrease in the size of the recirculation zone that is located in the leeside dune (x/d < 0). Furthermore, the recirculation zone is no longer homogenous and possesses more than one velocity minima in the velocity contours (for modeled case at x/d 5 0 and 3), which are potentially formed by reemerging flow out of the bed. This flow modifies the shape of the shear layer, which is no longer parallel to the bed as in the impermeable bed case (Figures 1e and 1g) but shows an undulating nature (Figure 1f). This disruption to the classical recirculation leeside separation cell and shear layer can be explained through analysis of the modeled W xt (Figure 1h). In the leeside of the dune at 3 and 6 x/d, jets of flow are detected leaving the hyporheic zone and entering the overlying boundary layer. This jet causes the formation of two counter-clockwise rotating vortices immediately behind the dune, as previously observed by Blois et al. [2014].
To enable a comparison between the permeable and impermeable bed conditions, a series of streamwise velocity profiles were extracted commencing at the crest and moving 4 x/d downstream, and then two further profiles spaced at 2 x/d ( Figure 2). Generally, U xt is greater over the impermeable bed than the permeable bed, with a maximum at a height >8 y/d. The only discrepancy from this trend is in the leeside of the dune at y/d 6, where U xt is more strongly negative due to the strength of the recirculation zone. In fact, by 4 x/d, U xt is always positive for the permeable bed case, demonstrating that no recirculation is predicted  Table 2 for comparison with the other numerical experiments. The flow is from right to left.

Water Resources Research
10.1002/2016WR019662 that far downstream. When W xt is analyzed, the greatest magnitude W xt , although negative, is again present for the impermeable case, thus showing the intensity of the recirculation zone. However, in the velocity profiles located in the leeside of the dune (profiles b-d, Figures 2b-2c), the proportion of the profile that is characterized by negative flow increases; for example, in Figure 2c, the primary negative flow is between 6 and 8 y/d for the permeable bed but 6 and 10 y/d for the impermeable bed. This would suggest, as demonstrated in Figure 1, that the intensity of the recirculation zone is reduced by the presence of hyporheic flow. When the TKE xt profiles are examined (Figure 2, row iv), the influence of hyporheic flow on the generation of turbulence in the boundary layer can be observed. For the first profile in the leeside of the dune ( Figure  2, IV: row iv, b), the turbulence generated from the recirculation zone over the impermeable bed exceeds that of the permeable bed, but further downstream in profiles c and d, the TKE xt for the impermeable bed exceeds the permeable case by a factor of 4 (Figures 1i and 1j). This demonstrates the considerable influence of turbulence generated by hyporheic flow in the leeside of the bed form, and how it generates more turbulence in the near-bed region than due to flow separation and reattachment.

10.1002/2016WR019662
In order to further investigate the influence of an impermeable bed, we consider streamlines, the Q criterion and bed shear stress (Figure 3). For the streamlines (Figures 3a and 3b), the tracer has been colored by the vertical (w2) component of velocity. For the impermeable case, the results are as expected following the schematic model of Best [2005], in which the streamlines pass over the crest of the dune and then either follow the shear layer or enter the recirculation zone. However, when the permeable bed is considered, the influence of both hydrostatically and hydrodynamically induced flows are observed. Just upstream of the bed form (210 x/d), flow is forced downward due to the local velocity gradients and then reemerges downstream of the crest in the recirculation zone [Thibodeaux and Boyle, 1987;Savant et al., 1987;Blois et al., 2014]. These flow paths extend to 2 to 3 y/d beneath the bed form, which is shallower than previous observations which suggest that the depth of flow penetration typically scales on the order of 80% of the bed form wavelength [Elliott and Brooks, 1997;Cardenas and Wilson, 2007a] (8 y/d), although this depth of flow penetration is also related to the free surface flow depth [Stonedahl et al., 2013;K€ aser et al., 2013]. These shallow hyporheic flow paths reemerge downstream in the recirculation zone but then pass through, and above, the streamlines formed by flow passing over the dune crest, potentially causing an instability that generates an undulating shear layer. A second set of streamlines are also identified much lower in the bed, at a depth of 4-5 y/d, and are potentially formed by hydrostatically induced forces and that upwell to the bed at 10, 17, and 20 x/d. Once the flow reaches the bed, it follows a similar pathway to the other upwelling flow in that it moves through the shear flow and up into a secondary shear layer located higher in the flow. These upwelling flows also influence the near-bed turbulence. The Q criterion (Figures 3c and 3d) detects areas of vortex generation and evolution. For the impermeable bed case (Figure 3c), these regions are located on the stoss side of the dune and at the crest where flow separation occurs. However, a different pattern is observed for the permeable bed case (Figure 3d), where the upwelling flow appears to generate turbulent structures with classical horseshoe geometry and that are associated with the upwelling flow predicted by the streamlines (Figure 3b). Furthermore, the results also display a possible coalescence of flow structures, with three bands of flow structures present from 0 to 8 x/d that then appear to merge into one large structure downstream. These flow structures have significant impact on the near-bed shear stress (Figures  3e and 3f). For the permeable case, patches of greater shear stress are located between these turbulent structures that possess a horseshoe geometry (Figure 3f). These patches have linear downstream striations (that become dominant >12 x/d) within them, and that are located above the contact points in the permeable bed where more upwelling flow is able to pass back into the boundary layer.

Permeable Dune Over a Porous Bed
The simulation of a permeable dune on a permeable bed (Figures 4 and 5) generally shows that all components of the flow are lower in magnitude than for the impermeable dune ( Table 2). The reason for this  Table 2 for comparison with the other numerical experiments. Flow is from left to right.

10.1002/2016WR019662
decrease in the magnitude of the peak flow velocity is attributed to the observation that although the permeable dune confines the flow and increases the velocity, the porous nature of the dune allows the passage of some flow through it, which subsequently reduces the peak velocity in the boundary layer above the dune. Furthermore, the permeable nature of the dune inhibits the formation of any recirculation zone in the leeside of the dune. As the flow is able to pass through the dune, U xt in the leeside of the dune is greater than zero, unlike the case of an impermeable dune placed on a porous bed (Figure 4a). This phenomenon has previously been observed in numerical simulation of laminar flow over a set of five permeable bed forms [Gualtieri, 2012]. The absence of flow separation and the presence of the partial passage of flow through the dune also prohibits the formation of a pronounced shear layer as observed in the previous cases. In the downstream direction, and due to the gradient of streamwise momentum in the vertical direction, upwelling occurs over an irregular spacing (Figure 4b). The permeable nature of the dune has previously been suggested to enhance turbulence [Huq and Britter, 1995], although high TKE xt in the dune area is not observed (Figure 4c).
The streamlines (Figure 5a) show a different pattern to those observed in the previous cases, with more flow paths being detected both through the bed and through the bed form. Some flow is forced downwards upstream of the bed form (x/d 5 11), although this is not of the same W xt magnitude and thus does not penetrate so deeply into the bed as that previously observed. Streamlines are also detected lower in the flow (between 22 and 24 y/d), and this flow travels further into the bed before it upwells back toward This change in the flow dynamics between the permeable and impermeable bed forms influences the generation of flow structures ( Figure 5b) and bed shear stresses (Figure 5c). The horseshoe-shaped structures in the leeside of the impermeable dune generated through reemerging flow are smaller in magnitude and spatial scale as compared with the permeable bed form. There are two zones (11 and between 16 and 24 x/ d; Figure 5b) where the return flow is associated with turbulent structures, although these are smaller in size than those detected over an impermeable dune on a permeable bed (Figure 3d). Furthermore, some flow structures can be seen upstream and along the stoss side of the dune (214 to 210 x/d). The smallerscale turbulent structures also result in lower magnitude shear stresses (up to 16 x/d downstream) as compared to an impermeable dune over a permeable bed, and with only a slight increase in shear stress at the end of the dune where the hyporheic flow upwells into the boundary layer flow.

Flow Associated With Two Asymmetric Dunes
In the fourth and fifth numerical experiments (Table 1), the flow fields around two bed forms in sequence were considered, over an impermeable and permeable bed, respectively (Figures 6-8). When flow over two dunes upon an impermeable bed is analyzed (Figures 6a, 6c, and 6e), flow over the first dune has an identical length flow separation zone to the isolated dune (Figure 1), with the point of reattachment being at c. Six bed form heights. However, the influence of the bed form affects the boundary layer flow on the next dune downstream, as the boundary layer develops on the stoss side of the second bed form. In particular, the shear layer, most easily identified in the TKE xt (Figure 6c), skims over the crest of the second dune, and this increases the length of flow reattachment downstream.
When two dunes are considered on a permeable bed, the flow on the stoss side of the first dune is again identical to the isolated dune (Figures 1d and 6b, 6d). However, flow separation in the leeside of the dune does not form two separate circulation cells as observed for the singular dune ( Figure 1d). However, similar to the single dune over a permeable bed (Figure 1d), the shear layer is not parallel to the bed, due to the hyporheic flow upwelling into the boundary layer in the dune leeside. Jetting of the flow can be observed at the end of the bed form (Figure 6d, x/d 5 0) and more significantly at x/d . 5 8. This nonuniform interaction of the shear layer with the second dune causes a complex flow pattern, with the recirculation zone in the lee of the downstream dune being located just below the crest, and being bounded by the dune leeface and returning water. The magnitude of flow in the lee of the second dune is greater in both intensity

10.1002/2016WR019662
and spatial extent than the upstream dune, with regions of high W xt (Figure 6d) leaving the bed at the end of the dune (x/d 5 26) and generating a region of TKE xt (Figure 6f).
To enable comparison between the permeable and impermeable bed conditions, a series of velocity profiles were extracted for locations at the crests and in the lee of both dunes (Figure 7). When the u-component is considered (Figure 7, row 1), flow over the impermeable bed is faster, as previously observed in the single bed form experiments (Figure 2). In fact, there is very little difference in the first three (Figures 7a-7c) profiles to those observed over the single dune experiments. Over the crest of the first dune, the difference is minimal (<0.05 U xt ), with this difference increasing only very slightly over the crest of the second dune. The U xt profile in the lee of the second dune (profile d) appears different, but this is largely due to a geolocation difference in that the dune length is not an exact integer multiple of x/d, and as such there is a pore space at the end of the dune. The profile here is thus altered due to the jetting of flow out of the bed, with the flow in the leeside being 0.25 U xt greater. This, as previously observed, is because the classical recirculation zone associated with dunes is not detected in the leeside (Figure 6d).
At the crest of the dune, W xt (Figure 7) again shows that the negative flow velocity is greatest for the permeable bed case. The location of this negative flow, especially in the lee of the first dune, is also higher in the flow showing this consistent pattern of the recirculation zone being lifted away from the bed due to the forcing of the returning flow (Figure 7, row 2). A similar pattern of flow is detected over the crest of the second dune but a more complex pattern is observed in the lee of the second dune. The increased upwelling flow further splits the recirculation cell into several smaller ones. The magnitude of this upwelling/jetting (profile d where y/d < 6) is high (W xt up to 0.4), explaining the disruption of the classical leeside flow cell. The magnitude of W xt decreases with depth within the bed (at depths <3 y/d W xt < 0.05) implying that the shallower driven flow has a greater influence on the flow dynamics. In all profiles, TKE xt is greater in the permeable bed (Figure 7, row 3).
Finally, the flow fields of the two dunes are analysed through streamlines, Q criterion and bed shear stress plots (Figure 8). The streamlines over the two impermeable dunes show very similar patterns to those observed over the single impermeable dune, and identify an area of recirculation spanning a distance 6 x/d in the lee of the dune. The recirculation cell behind the second dune appears not to touch the bed and the streamlines are closer together (Figure 8a, between 26 and 32 x/d). However, when the flow of the permeable bed is considered, a far more complex flow pattern is observed (Figure 8b). Although flow is forced down into the bed in front of both dunes, it is the distance downstream where the upwelling flow returns to the boundary layer that influences the overall flow field. Flow in front of the first dune is forced down to a depth of 3-4 y/d. Some of the upwelling flow in this case returns to the boundary layer flow between the two dunes (at 4-10 x/d) and influences the shear layer, as was predicted in the case of the single dune, and forms an undulating shear layer (Figure 3b). However, some of this flow continues under the second dune, and is added to the flow that is being hydrodyamically forced down into the bed in front of the second dune (at 16 x/d). The flow forced down in front of the second dune appears to be forced into the bed at a much steeper gradient (i.e., it goes deeper into the bed over a shorter distance) both into bed and when it upwells in the leeside of the second dune. The return flow appears at the end of the dune but also at 26 x/d where it has a high vertical velocity and further disrupts the flow in the lee of the dune (Figure 8b). The flow appears to coalesce with the shear layer and increases the vertical thickness of the shear layer. This can be confirmed through the Q criteria ( Figure 8c) where horseshoe vortices are again detected. Flow behind the first dune on a permeable bed shows a similar pattern to that observed behind a single dune (Figure 3d), although a far more complex pattern exists behind the second dune. The jetting flow in the lee of the dune forms a vortex structure and this continues downstream, with the number of structures identified by the Q criterion increasing and the spacing between these structures decreasing. Finally, this flow field has an influence on the bed shear stress, which increases in magnitude downstream of the dune, with values approximately 5 times higher in the case of a permeable bed (compare Figures 8e and 8f).

Discussion
The paper presents and applies a new numerical model that is capable of predicting flow both in the boundary layer and within the hyporheic zone. The method enables a numerical representation of both the hydrostatic and hydrodynamic forcing's that drive flow exchanges between the two flow regions. Validation between numerical and experimental results showed good quantitative agreement (experiment 1: r 5 0.94; experiment 2: r 5 0.74) and produced realistic flow field predictions, thus giving confidence in the applicability of the model.
In this application, we chose to investigate flow over dune bed forms constructed of coarse (gravel sized) particles. These are morphologically similar to a riffle and pool system which has previously been studied [e.g., Harvey and Bencala, 1993] including ones in gravel-bed rivers [Tonina and Buffington, 2009]. In this application similar dynamic hyporheic exchange flows (HEF) [Wondzell and Gooseff, 2013] are predicted such as curved flow paths tending away from the bed under the bed form. However, due to the process representation of the model more complex flow dynamics can be predicted. Therefore, a new conceptual model of flow over dunes has been proposed ( Figure 9) that summarizes flow downstream of coarsegrained dunes on permeable beds and demonstrates that this is fundamentally different to that over impermeable beds. For the control case, that of an impermeable dune on an impermeable bed (Figure 9a), the five distinct flow zones reported in the schematic model of Best [2005] are unchanged. However, this flow field is significantly altered when the bed is permeable (Figure 9b). In this case, the flow is hydrodynamically forced into the bed on the stoss side of the dune by local velocity gradients that are generated by flow around, or over, the bed form and that result in momentum transfer into the bed [Boano et al., 2014]. These flow paths are 2-3 y/d beneath the bed form and are qualitatively consistent with earlier studies of mixing zones beneath dunes [Sawyer and Cardenas, 2009;Jin et al., 2010;Hester et al., 2013]. The flow reemerges downstream of the crest within the recirculation zone, as observed in previous studies [Thibodeaux and Boyle, 1987;Savant et al., 1987;Blois et al., 2014], and this upwelling causes a reduction in the size and intensity of the recirculation zone. Furthermore, the recirculation zone is no longer homogenous, and jets of flow from the reemerging flow from two counter-clockwise rotating vortices immediately behind the dune, replicating well observations from physical experiments [Blois et al., 2014]. These jets of upwelling flow occur at several points downstream of the dune. Analysis of the numerical data using the Q criterion suggests that these jets generate coherent flow structures with a horseshoe-shaped geometry that coalesce with other flow structures downstream of the dune, and that these weaken/dissipate the shear layer so that it adopts Water Resources Research 10.1002/2016WR019662 an undular boundary. The flow structures detected by the Q criterion also have a significant impact on the bed shear stresses, with high patches of shear stress being located between these turbulent structures, perhaps due to downwelling fluid that, for continuity, has to replace the upwelling horseshoe-shaped vortices. These upwellings are likely driven by pressure gradients formed around the dune. Separation of flow in the leeside of a dune generates a low-pressure region [Maddux et al., 2003b], with the sudden expansion of flow causing the hydrostatic pressure to increase in the direction of flow [Motamedi et al., 2014]. Over an impermeable bed, the generation of this adverse pressure gradient causes flow to separate forming vortices along the shear layer [e.g., Buckles et al., 1984] that increases the turbulence intensity and Reynolds stresses [Raudkivi, 1966;Engelund and Fredsoe, 1982;McLean and Smith, 1986;Mendoza and Shen, 1990;Nelson et al., 1993;McLean et al., 1994;Bennett and Best, 1995;Best, 2005]. However, over a permeable bed, the low pressure in the lee of the dune generates significant pressure gradients both across the bed and into the pores beneath [Detert et al., 2004;Klar et al., 2004]. These potentially act to draw flow back from the Brinkman layer through upwelling into the boundary layer above. This causes the flow separation cell in the lee of the dune to split into two, as previously observed by Blois et al. [2014], with the separation cell closest to the leeside of the dune being displaced up the lee slope (Figure 9b).
When two bed forms are considered in sequence over both impermeable ( Figure 9c) and permeable ( Figure  9d) beds, the trends observed above are enhanced. For the impermeable bed, the boundary layer is still developing by the time it interacts with the second dune, although the wake turbulence decays and diffuses up into the boundary layer [Maddux et al., 2003a] that has the influence of slightly increasing the size of the recirculation zone. However, in general there is little modification to the overall flow field. The same is not true for the permeable case (Figure 9d), especially behind the second bed form where upwelling jets of flow: (i) lift the separated flow away from the bed; (ii) modify the shear layer through coalescence with vortices generated by flow reemerging from the bed, which acts to increase the vertical extent of the shear layer and; (iii) causes the shear layer to undulate rather than be parallel to the bed.
When a permeable dune overlying a permeable bed is considered, the flow field is even more substantially modified. An impermeable dune may be viewed as causing a blockage to the flow that creates a pressure difference between the stoss and lee slopes, thus generating form drag [Maddux et al., 2003]. However, in the case of a highly permeable bed form, flow is retarded but not blocked by the dune, and thus some of the flow passes through the bed form. This flow through the dune further modifies the leeside flow structure, which for the highly permeable dune considered herein does not possess flow separation in its leeside.
These findings improve our understanding the dynamics of gravel bed rivers as any natural gravel river consists of some nature of bed form (e.g., transverse clast dams, longitudinal clast ridges, transverse ribs, unit bars, and clusters). Furthermore, real rivers are more heterogeneous than the simple examples reported. The size, spacing, and sequence of bed forms along the reach also affect hyporheic exchange patterns [Anderson, 2005;Gooseff et al., 2006;Boano et al., 2014] including the location and magnitude of upwelling and downwelling. Furthermore, flow in gravel bed rivers is unsteady, and typically shallow, with the ratio of mean flow depth to typical roughness height seldom exceeding 10-20 in floods and less than 5 during normal flow conditions. During base flow conditions bed-forms cause form (pressure) drag and thereby drive hyporheic exchange flow. However, Storey et al. [2003] examined hyporheic exchange in a pool-riffle sequence for both flood and normal flows and demonstrated that at higher flows the effect of the riffle diminished substantially reducing the hyporheic exchange. This current work needs to be expanded to investigate; (i) the effect of size, spacing, and sequence of bed forms; (ii) the influence of particle size; and (iii) how the depth and velocity of boundary layer flow effects the magnitude of the processes driving hyporheic exchange.
The model can also be applied to understand the sediment dynamics of a gravel bed river where there are sequences of bed forms. A coarser matrix can form in the topographic lows between bed forms through the action of winnowing. The winnowing would be enhanced through the increased shear stresses predicted through the generation of turbulent vortices in the lee of the dune (Figure 8f). However, this coarser matrix, especially if located near to the stoss of the next downstream dune, would be a preferential location for flow to be driven into the bed. This process information, therefore, may be used to understand how surface textures change. For example, Singh et al. [2012] hypothesize that finer material entrained from a dune crest may be deposited in the trough and infiltrate into the subsurface sediment pores, causing the subsurface material to become finer. This may also affect the bed surface texture through colmation, clogging of the top layer of the channel sediments [Brunke and Gonser, 1997], which would lead to a reduction of pore volume, consolidation of the sediment matrix, and decreased permeability of the bed. This would hinder both hyphoeric exchange processes and also modify the bed roughness and potential depositional environment across the riverbed [Buffington and Montgomery, 1999], which would lead to a spatial variation in the hydraulic conductivity [Genereux et al., 2008]. Previous flume experiments have demonstrated patterns of fine-sediment deposition, with fine sediment preferentially deposited in downwelling zones [Packman and MacKay, 2003;Rehg et al., 2005], which has been suggested as an explanation in the differences in hydraulic conductivity between upzones and downzones [Scordo and Moore, 2009]. Furthermore, the action of the returning hyporheic flow into the dune leeside in the form of the strong jets of flow, that create horseshoeshaped vortices revealed herein, may act to both flush fines from the bed or delay their deposition within the bed, this potentially creating more openwork gravel textures in the immediate dune leeside. It is also worth considering how these returning horseshoe-shaped vortices could affect the boundary layer flow. Recent experimental work looking at coherent flow structures over gravel beds [Hardy et al., 2016] identified large individual packets of fluid, that contain several smaller scales of fluid motion, were initiated at the bed through shear that generate a bursting mechanism. No consideration in these experiments was taken of flow passing through the bed and the possible influence these structures have on the overall nature of the boundary layer flow.
Although there is a need to test a range of scenarios considering both bed form and hydraulic conditions there is also a need to scale up these findings. It is not computationally efficient to apply this type of model to the whole river profile, however, modeling individual morphological features such as a meander bend would be feasible. This would enable the prediction of flow pathways [e.g., Cardenas and Zlotnik, 2003], and the speed of flow where water flowing along deeper pathways would travel at much lower velocities than shallower flow and thus predict a range of travel times [Stonedahl et al., 2010;Cardenas, 2008;Zarnetske et al., 2011;Bardini et al., 2012;Marzdari et al., 2012;Arnon et al., 2013]. This form of process understanding has both ecological [Bardini et al., 2012;Hester et al., 2013] and morphodynamic implications [Harrison and Clayton, 1970]. For example, detailed quantification of the flow field around bed forms resting on a permeable bed could inform development of more complete models of solute and particle transport [Packman et al., 2000;K€ aser et al., 2013;Hester et al., 2013], nutrient and carbon cycling [Jin et al., 2010[Jin et al., , 2011Bardini et al., 2012;Gomez-Velez et al., 2015], and finally help understand the spawning habitat of some fish [Baxter and Hauer, 2000]. Furthermore, numerical studies such as that detailed herein can allow more complete investigations of the influence of bed form morphology in controlling a significant proportion of such mass, momentum, and particulate exchange between the flow and bed [Thibodeaux and Boyle, 1987;Lu and Chiew, 2007].

Conclusions
Through the development and application of a numerical model which allows hyporheic flow to be predicted, the effects of bed permeability on bed form dynamics have been assessed. Our key findings are: 1. Flow over of an isolated impermeable bed form over a permeable bed is forced into the bed upstream of the dune and returns to the boundary layer flow at the leeside, in the form of returning jets that generate horseshoe-shaped vortices. The reemerging flow significantly influences the leeside flow, modifying the separation zone, lifting the shear layer adjoining the separation zone away from the bed. The returning flow in the leeside appear as jets of flow that generate horseshoe vortices 2. Flow over a permeable dune on a permeable bed generally shows that all the velocity components of the flow are lower in magnitude than for the impermeable dune. The decrease in the magnitude of the peak flow velocity is attributed to the observation that although the permeable dune confines the flow and increases the velocity, the porous nature of the dune allows the passage of some flow through it, which subsequently reduces the peak velocity in the boundary layer above the dune. This prevents the formation of any recirculation zone in the leeside of the dune and the absence of flow separation also prohibits the formation of a pronounced shear layer. 3. When there are multiple bed forms the flow over the downstream dune is influenced by the developing boundary layer on the leeside of the upstream dune. For the permeable bed case, the upwelling flow lifts the separated flow from the bed, modifies the shear layer through the coalescence with vortices generated by the returning flow, and causes the shear layer to undulate rather than be parallel to the bed.
These results demonstrate the significant effect that bed permeability has on the flow over bed forms that may be critical in affecting the flux of water and nutrients between the boundary flow and surface. Data presented in this paper can be obtained by contacting RJH. The work was carried out under NERC grant NE/ K003194/1 awarded to RJH and GHSS. We are grateful to the Associate Editor Christophe Ancey and three anonymous referees for providing helpful comments that have led to significant improvements in this manuscript.