The role of velocity , pressure , and bed stress fluctuations in bed load transport over bed forms : numerical simulation downstream of a backward-facing step

Bed load transport over ripples and dunes in rivers exhibits strong spatial and temporal variability due to the complex turbulence field caused by flow separation at bedform crests. A turbulence-resolving flow model downstream of a backward-facing step, coupled with a model integrating the equations of motion of individual sand grains, is used to investigate the physical interaction between bed load motion and turbulence downstream of separated flow. Large bed load transport events are found to correspond to low-frequency positive pressure fluctuations. Episodic penetration of fluid into the bed increases the bed stress and moves grains. Fluid penetration events are larger in magnitude near the point of reattachment than farther downstream. Models of bed load transport over ripples and dunes must incorporate the effects of these penetration events of high stress and sediment flux.


Introduction
The details of turbulent flow over dunes and ripples in rivers and oceans have been described by field and laboratory experiments (see Best, 2005, for an extensive review), as well as high-resolution, turbulence-resolving numerical simulations (Shimizu et al., 1999(Shimizu et al., , 2001;;Nelson et al., 2006;Zedler and Street, 2001;Omidyeganeh and Piomelli, 2011;Grigoriadis et al., 2009;Stoesser et al., 2008;Chang and Constantinescu, 2013).However attempts to couple turbulence to the transport of sediment over bedforms have usually relied on empirical formulas, wherein the sediment flux is either a direct function of boundary shear stress or indirectly through entrainment rate and deposition rate formulas (Niemann et al., 2011;Nguyen and Wells, 2009;Giri and Shimizu, 2006;Chou and Fringer, 2010;Paarlberg et al., 2009;Kraft et al., 2011) (although see Nabi et al., 2013, andPenko et al., 2013).Unlike suspended sediment fields, experiments detailing the spatiotemporal pattern of bed load transport over ripples and dunes have not been reported.Grass and Ayoub (1982) hypothesized that the mean bed load sediment flux could be calculated as the integral of the probability density of bed stress due to turbulence times the sediment flux as a function of stress.They further hypothesized that the bed stress distribution could be determined from the distribution of near-bed downstream velocity.The experiments of Nelson et al. (1995) simultaneously measuring sediment flux and near-bed fluid velocity over a flat bed and downstream of a backward-facing step showed that the relationship between near-bed fluid Reynolds stress and bed load transport was not simple, and the spatially varying distribution of velocity fluctuations relative to the shear velocity must be considered in formulating transport relationships over ripples and dunes.They also found that there was not a simple monotonic relationship between instantaneous, downstream, near-bed velocity and sediment flux.Specifically, longer duration positive fluctuations of near-bed velocity were found to transport more sediment per unit time than shorter duration events.As such, the hypothesis of Grass and Ayoub (1982) needs significant modification to be useful downstream of separated flows.
Published by Copernicus Publications on behalf of the European Geosciences Union.
Experimental measurement of turbulence is often limited to a time series of fluid velocity components at a single point or, more rarely, several points.In such instances, the detection of spatially and temporally evolving turbulence structures is difficult.Quadrant analysis has been used to detect certain types of turbulence structures (Lu and Willmarth, 1973;Bogard and Tiederman, 1986).Quadrant analysis involves joint examination of the fluctuating components of fluid velocity in the downstream, x, and bedperpendicular, z, directions.u and w are the downstream and bed-perpendicular fluctuating components of fluid velocity.With u and w measurements plotted on a twodimensional graph, the first quadrant (Q1) is a point with u > 0 and w > 0; this is also known as an outward interaction.Quadrant 2 events (Q2, u < 0 and w > 0) are termed bursts or ejections, quadrant 3 events (Q3, u < 0 and w < 0) are called inward interactions, and quadrant 4 events (Q4, u > 0 and w < 0) are known as sweeps.Q2 and Q4 events transport downstream momentum toward the bed, and are thus positive contributions to the Reynolds stress component, −ρu w , whereas Q1 and Q3 events are negative contributions.
Schmeeckle ( 2014) used a coupled turbulence-resolving numerical model of flow and a particle model of sediment motion to simulate the interaction between turbulence and sediment movement over a flat bed.Vortical structures embedded within broader sweep structures were found to bring fluid into and out of the bed, and were sites of sediment entrainment and transport.In this article I extend the model of Schmeeckle (2014) and Furbish and Schmeeckle (2013) to the case of bed load transport downstream of a backwardfacing step, largely matching the experiments of Nelson et al. (1995).Quadrant analysis is extended to include sediment flux, bed stress, and fluid pressure.Flow over a backwardfacing step, like that over bed forms, causes flow separation, but does not have the complicating effect of flow acceleration by an upstream sloping bed.

Methodology
The fluid is modeled by the large eddy simulation (LES) technique in which the spatially filtered Navier-Stokes equations are integrated using the finite-volume method.The equations of motion of each sediment grain are integrated over time using the distinct element method (DEM).The sand grains are assumed to be spheres, and forces between particles are calculated when grain boundaries overlap.The LES and DEM models are coupled in momentum.The flow field is interpolated to the particle centers and used to derive fluid forces on the particles.In turn, each fluid force acting on the particles is given as a resistance term to the fluid momentum equations at the fluid cell containing the center of the particle.Only drag, pressure gradient, and buoyancy forces are included as fluid-particle forces.The bed of par-ticles is about three to four grain diameters thick above the lower fluid boundary, and non-moving particles of the particle bed rapidly damp the fluid velocity.The details of the numerical model reported here are the same as reported in Schmeeckle (2014).
The flow magnitude, step height (x step ), particle diameter (D) and density (ρ s ), and flow depth of the numerical simulations are specified to nearly match the experiments of Nelson et al. (1995).The computational domain extends 0.2 m upstream of the 0.04 m backward-facing step and 1.2 m (30 step heights) downstream of the step, and 0.1 m (2.5) step heights across stream.The fluid domain in the vertical dimension extends 0.16 m above the step, and downstream of the step the vertical domain is 0.2025 m.At rest, the topmost particles are roughly 0.0025 m above the lower wall.Thus, the flow height downstream of the step is about 0.2 m from the bed of particles to the top of the numerical domain.The grid used in this study is a structured mesh of 4 655 000 hexagonal cells.The grid is evenly spaced in the downstream and cross-stream directions.The downstream and cross-stream grid lengths are 0.002 and 0.00143 m, respectively.The vertical grid spacing is nonuniform, with smaller grid cells containing the particles and near-bed flow.The vertical grid dimension is also significantly reduced in a zone containing the separation bubble shear layer.The vertical dimension of cells containing particles at the bottom of the numerical domain is 0.00025 m.The downstream and cross-stream grid dimensions are slightly larger than the diameter of the particles, but the particle diameters are about 3.6 times larger than the vertical grid dimension.As such there is rarely more than one particle in a grid cell.
If fluid and particles are coupled in mass and momentum, the fluid solver becomes unstable when the particles are of the same size as the fluid grid cells.Each particle is treated as a source of resistance in the cell where the particle(s) center is located.It is possible to smooth the effect of a particle over a broader number of cells to achieve a stable algorithm, but this smooths the sharp interface between the bed and the overlying relatively sediment-free fluid.Smoothing of the bed interface was deemed not appropriate, because most of the fluid momentum is damped within a grain diameter below the bed interface.Here, fluid and particles are coupled in momentum but not in mass.The flow around each particle is not directly modeled; only the damping of flow by the integrated force of each particle is modeled.Thus, the flow separation and turbulence generated by flow around particles and in the interstices of particles is not modeled.Further, the calculated pressure field within the bed may be different than reality, because the fluid continuity equation does not account for the reduced volume of fluid within the bed.It is difficult to predict how these assumptions degrade the fidelity of the simulations.However, the integrated momentum effect of each particle on each fluid cell (and vice versa) is calculated in the model and should lead to relatively accurate fluid simulations at the scale of the grid.
The domain of the DEM model begins at the step in the x direction but otherwise coincides with the fluid domain boundaries.The particle domain is periodic in the downstream and cross-stream directions.The diameters of the 415 000 particles in this simulation are randomly drawn from a normal distribution with a median of 0.9 mm and a standard deviation of 0.1 mm.The diameters are varied to avoid close packing arrangement of the bed during the simulation.The particle parameters for the DEM model are the same as in Schmeeckle (2014) except that the Young's modulus is increased to 5 × 10 6 Pa.
Boundary conditions for fluid velocity are no shear at the upper boundary, periodic conditions in the cross-stream direction, and zero gradient at the outlet.The no-slip condition is applied at the lower boundary, but the fluid velocity becomes negligible before reaching this boundary because of the presence of a bed of particles above it.The inlet boundary condition is specified as the velocity 0.15 m downstream of the inlet.This is similar to a periodic boundary condition wherein the inlet and outlet are the same velocity, but the recycled velocity is taken before the backward step, thus ensuring fully developed boundary layer turbulence upstream of the backward-facing step.

Results
Prior to recording simulation results, the flow reached dynamic equilibrium after about 30 s of simulated time.Results reported here are for 20 s of simulated time.This length of time provided adequate statistics, but it was not so long that the bed of particles developed bedforms and areas of the bed without sufficient numbers of particles.However, bed elevation changes of one to two particle diameters were apparent in response to the passage of individual turbulent events.The position, velocity, and fluid force of each particle of known diameter is recorded at 40 Hz.Similarly, the fluid velocity and pressure are saved simultaneously at 40 Hz along two near-bed horizontal slices at z = 1 mm and z = 5 mm and along a slice perpendicular to the cross-stream at the center of the numerical domain, y = 0.05 m.The lower boundary of the fluid and particle domain is at z = −0.0025m and the topmost particles of the bed at rest are at approximately z = 0.A local depth-integrated downstream sediment flux, q sx , is calculated by summing the product of each particle volume and velocity that is found in a 0.01 × 0.01 m horizontal area of the bed and then dividing by the local bed area (i.e., division by 0.0001 m 2 ).Downstream bed shear stress, τ bx is calculated by summing the downstream component of fluid force acting on all particles with centers contained in the same 0.01 × 0.01 m horizontal grid, and then dividing by the local bed area.The saved fluid velocity and pressure at z = 1 mm are extracted at points in the center of the grid of local bed areas used to calculate boundary stress and sediment flux.In this manner, the data examined in this article It should be noted that the bed stress, as defined here, is the sum particle force per bed area.It is not the near-bed Reynolds stress component, −ρu w .Averaged over sufficient time, these two quantities are equivalent by a balance of forces.However, at a particular instant in time, there are fluid accelerations, and τ bx = −ρu w , except by coincidence.
Simultaneous visualization of the fluid pressure fluctuation, p = p − p; downstream fluid velocity, u; and particle velocity magnitude, |U|, reveals a positive covariance of particle motion with near-bed, downstream fluid velocity and fluid pressure (Fig. 1).Low-frequency, cross-channel variations in pressure are apparent in Fig. 1, which were also noted in the direct numerical simulations of Le et al. (1997).Figure 1 shows that positive pressure fluctuations are associated with both large near-bed downstream fluid and particle velocity magnitude.There are small areas of the bed with large negative vertical velocity (red areas of Fig. 1d).Fluid that penetrates the bed leads to neighboring areas where fluid exits the bed (blue areas of Fig. 1d).These large fluctuations in vertical velocity are associated with significant sediment motion (Fig. 1c).
Figure 2 shows the temporal statistics (mean, 10th, and 90th percentile) vs. downstream distance for u, w, p, τ bx , and q sx .The position of the point of reattachment is plotted at a vertical dotted line in the five plots of Fig. 2. Figure 2a www.earth-surf-dynam.net/3/105/2015/Earth Surf.Dynam., 3, 105-112, 2015 shows that the mean near-bed velocity (at z = 1 mm) increases rapidly near the point of reattachment, and it increases, albeit much less rapidly, all the way to the downstream outlet.Interestingly, the difference between the 90th and 10th percentile of velocity is larger near the outlet than in the reattachment zone.However, the difference between the 90th and 10th percentile of bed stress (Fig. 2d) is smaller near the outlet.The 10th percentile transport rates are essentially zero (or slightly negative) downstream of flow reattachment (Fig. 2e).The mean sediment transport rate (Fig. 2e) increases rapidly downstream of reattachment and does not show a peak in transport at 20 step heights as do the results of Nelson et al. (1995).However, they suggested that the peak could be the result of sampling error.The fluid pressure (Fig. 2c) rises rapidly from the recirculation region through the zone of reattachment (from about x/h step = 3 to 7).The largest magnitude of this upstreamdirected pressure gradient is about 400 Pa m −1 , which leads to a stress of about −0.5 Pa at x/h step = 5.This "pressure gradient stress" is about one-third to one-quarter of the negative bed stress in the recirculation and reattachment zone.However, this stress is distributed throughout the bed of particles, and the resulting pressure force on individual grains is more than an order of magnitude smaller than is required to entrain the topmost grains that are able to move.
While Fig. 1 qualitatively shows the spatial covariance of some of the fluid and particle variables, Fig. 3 shows some of the significant temporal correlation pairs of variables u, |w |, τ bx , and p.The absolute magnitude of the vertical velocity fluctuation, |w | is used rather than w because transport was found to peak when the fluctuations of vertical velocity were high.It is perhaps unsurprising that u is positively correlated with τ bx and q sx , but Fig. 3 also shows the positive correlation with fluid pressure, p.

Permeable splat events
Given that the force on bed grains results primarily from fluid drag, it is somewhat paradoxical that the temporal variance in the bed stress is much larger near reattachment, despite the smaller variance in downstream velocity at reattachment, relative to farther downstream.This apparent paradox is due to the fluctuations in vertical velocity being much larger near the point of reattachment than farther downstream (Fig. 2b).A large negative vertical velocity brings high downstream fluid velocity into the bed, thus creating peak bed stresses.Consider the plots in Fig. 2 between about x/h step = 8 and x/h step = 12. Figure 2a shows that the 90th percentile of the near bed velocity continues to increase downstream.However, Fig. 2d shows that the 90th percentile in boundary shear stress is as high or higher than farther downstream.Figure 2e shows that the transport, similarly, is as large as farther downstream, despite having a lower near-bed downstream velocity.Figure 2b shows that the magnitude of the 10th percentile of vertical velocity is larger in this zone than farther downstream.These large negative vertical velocities bring highmomentum fluid into the bed, increasing the force on grains and causing transport.When a localized volume of fluid approaches and impinges on an impermeable boundary, the boundary-normal velocity must stagnate, and the fluid gets redirected to move parallel to the wall.Perot and Moin (1995) refer to these wall impingements as "splat events", and Stoesser et al. (2008) note the occurrence of splats near flow reattachment in their simulations of turbulence over dunes.In the simulations reported here, the bed is a permeable boundary, and splats can penetrate the bed.To satisfy fluid continuity, infiltration of the bed by a splat must be accompanied by exfiltration of the bed surrounding the splat.Permeable splat events are apparent near reattachment in Fig. 1d (areas of intense red and blue) and the dynamics of the splat events are apparent in the Supplement animation of Fig. 1.Schmeeckle (2014) remarked that significant entrainment of bed load grains occurs on the boundaries between areas of bed infiltration and exfiltration.
The very large negative stresses in the recirculation region which peak at about x/h step = 4 in Fig. 2d are due to a negative mean vertical velocity, w, at the particle bed (Fig. 2b) and the large negative vertical velocity fluctuations, w 10 .There is a mean penetration of fluid into the bed, and there are permeable splat events.The downstream fluid velocity is also negative in the bed of particles due to the adverse pressure gradient.However, once again, the drag forces produced on the grains in this region are more broadly distributed through the bed, in contrast to the bed well downstream of flow reattachment, where the boundary shear stress is concentrated on only the topmost particles.This set of conditions also explains why the mean transport rate and nearbed downstream velocity are negligible even though Fig. 2d shows that the mean boundary shear stress is negative at the point of reattachment.

Quadrant analysis
Recall that the simulation data were collected for u, w, q sx , p, and τ bx simultaneously at a horizontal grid of points.In Fig. 4 all of the data were aggregated from all points downstream of x/h step = 12.5.Figure 4a shows the frequency of u -w paired bins, and the predominance of burst and sweep events is apparent.In Fig. 4b, q sx is summed for each uw bin.The bins are then normalized by dividing all bins by the bin with the maximum sum of q sx .Figure 4b shows that most of the transport (about 80 %) takes place during sweeps and outward interactions.This result is consistent with Nelson et al. (1995).In Fig. 4c, τ bx is summed for each u -w bin, and each bin is normalized by the largest magnitude bin.Percentages for each quadrant in Fig. 4c are given by the sum τ bx and divided by the total deviation, |τ bx |.Sweeps are associated with high bed stress, and bursts are associated with low bed stress.The pressure deviation is summed in each u -w bin in Fig. 4d and normalized by the magnitude of the bin with the largest magnitude.Percentages for each quadrant are given by the sum p and divided by the total deviation, |p |.Sweeps and outward interactions are associated with high-pressure events and bursts and inward interactions are associated with low-pressure events.
The spatial correlation between sweeps and outward interactions and between bursts and inward interactions is apparent in Fig. 5. Areas of the bed occupied predominately by sweeps and outward interactions are also areas with high fluid pressure, large particle forces, and large sediment fluxes.Conversely, bursts and inward interactions are associated low pressure, small particle forces, and small sediment fluxes.Sweeps and inward interactions occur together when a broad volume of fluid moves toward the bed, bringing with it high downstream velocity.Such a situation is apparent in Figs. 1 and 5 at x/h step ≈ 17.When a broad sweep impinges on the permeable bed, there is infiltration and exfiltration at spatial scales smaller than the broader sweep structure.Areas of exfiltration are apparent as outward interactions.Downstream-elongated structures of high-and low-speed fluid begin to emerge downstream of flow reattachment (Fig. 5a) (as also noted by Le et al., 1997).These emerging streaks also produce streaks of high particle forces (Fig. 4b) and particle motion (Fig. 4c).

Conclusions
Temporally averaged bed stress is not sufficient to specify the rate of bed load transport downstream of separated flow (compare Fig. 2d and e).Most of the transport takes place at high-stress events that are associated with both high downstream velocity and high-magnitude vertical velocity events (Fig. 4b).The temporal distribution of bed stress is broader near flow reattachment than farther downstream (Fig. 4d), even though the temporal distribution of near-bed downstream velocity is less broad near flow reattachment than downstream."Near-bed" and "in the bed" fluid velocities are different.In this study near-bed was specified at z = 1 mm, which is about one sand grain diameter above the top of the bed.Negative vertical velocity events (splats) bring high downstream momentum fluid into the bed, and those bed penetration events are stronger near flow reattachment (Fig. 2b).Consequently, the 90th percentile of stress and the mean sediment flux reaches a peak in a relatively short distance downstream of reattachment.This provides a probable explanation of the findings of Nelson et al. (1995) that instantaneous, near-bed downstream velocity was not sufficient to specify the instantaneous sediment flux; the actual force on bed particles is also dependent on the penetration of turbulence structures into the bed.The upstream inclination of the stoss of bed forms, relative to the flat bed considered here, is expected to increase the intensity of fluid penetration events near flow reattachment.

Figure 1 .
Figure 1.Visualization of the downstream fluid velocity, u; vertical velocity, w; particle velocity magnitude, |U|; and fluid pressure fluctuation, p ; at an instant in time.(a) Downstream velocity on a vertical slice at the middle of the cross-stream domain.(b) Downstream velocity on a horizontal slice at z = 1 mm.(c) Particle velocity magnitude.(d) Vertical velocity on a horizontal slice at z = 1 mm.(e) Fluid pressure fluctuation on a horizontal slice at z = 1 mm.(f) Fluid pressure at the middle of the cross-stream domain.An animation associated with this figure can be found in the Supplement.

Figure 2 .
Figure 2. Temporal mean and the 10th and 90th percentile flow and transport parameters are plotted against distance downstream relative to step height, x/h step .(a) Downstream fluid velocity at z = 1 mm, (b) vertical fluid velocity at z = 1 mm, (c) fluid pressure at z = 1 mm, (d) bed shear stress, and (e) depth-integrated downstream sediment flux.N95 is the measured sediment flux ofNelson et al. (1995).A smoothed line of a moving average of q * of all points within 0.025 m upstream and downstream is also shown.

Figure 3 .
Figure 3. Correlation coefficients vs. time lag for various pairs of flow and transport variables as indicated in the legend in (c).The lag is of the second variable relative to the first variable in the legend shown in (c).

Figure 4 .Figure 5 .
Figure 4. Flow and transport data binned by downstream and vertical velocity fluctuation pairs, u -w .Data are aggregated for all points downstream of x/h step = 12.5.(a) Bin counts normalized by the largest bin.The total percentage of counts for each quadrant are shown.(b) The sum of downstream sediment transport in each bin, q sx , normalized by the bin with the largest transport sum.The percentage of transport for each quadrant is shown.(c) The sum of the downstream bed stress fluctuation for each bin, normalized by the magnitude of the bin with the largest magnitude.Percentages shown in each quadrant are for the sum of the stress fluctuation, τ bx , divided by the total absolute deviation, |τ bx |.(d) The sum of the fluid pressure fluctuation for each bin, p , normalized by the magnitude of the bin with the largest magnitude.Percentages shown in each graph are for the sum of the pressure fluctuation divided by the total absolute pressure deviation, |p |.