The Influence of Three-Dimensional Topography on Turbulent Flow Structures Over Dunes in Unidirectional Flows

Dunes are the most prevalent bedform present in sand‐bedded rivers and their morphology typically comprises multiple scales of three‐dimensional topography. However, our understanding of flow over dunes is predicated largely on two‐dimensional models, a condition which is rare in nature. Here, we present results of Large Eddy Simulations over a static, three‐dimensional dune field, using a two‐ and three‐ dimensional topographic realisation, to investigate the interaction between bed topography and turbulent flow structures. We show that flow over two‐dimensional bedforms increases the velocity over the stoss slope and reduces the size of the leeside separation zone as compared to 3D topography. Flow over three‐dimensional bedforms generates twice as many vortices as over two‐dimensional bedforms, and these vortices are longer, wider and taller than flow over their two‐dimensional counterparts. Turbulence is dominated by hairpin‐shaped vortices and Kelvin‐Helmholtz instabilities that interact with the bed in the brink point region of the dune crest and down the lee slope, and generate high shear stresses for long durations. These results are used to propose a new conceptual model showing the differences between flow over two‐ and three‐dimensional bedforms. The findings highlight how the size, morphology and stacking of coherent flow structures into larger flow superstructures may be critical in sediment entrainment, and may dictate the relationship between event duration and magnitude that drives sediment impulses at the bed. This will ultimately lead to an increased in the three‐dimensionality of bedform morphology.

When the lee-side angle is less than 10°, permanent flow separation is not observed, and when the lee-side angles is less than 4° there is no flow separation (Best & Kostaschuk, 2002;Lefebvre & Winter, 2016;Lefebvre et al., 2016); 3. flow reattachment 4-6 dune heights downstream of the crest (Engel, 1981); 4. a bounding shear layer between the separated flow and streamwise flow that expands downstream and includes vorticity generated by Kelvin-Helmholtz instabilities (Abad, 2008;Best, 2005a;Frias & Abad, 2013).Above the bounding shear layer, a quasi-inviscid interior region exists that is affected by spatial accelerations and the wakes of the upstream dunes (Nelson & Smith, 1989); and 5. an internal boundary layer that grows from the reattachment point along the stoss slope of the next downstream dune (Best, 2005a).
Although the simplified model summarized above has been developed from an idealized case that rarely, if ever, exists in natural environments, it forms the basis for interpreting flow over alluvial dunes.However, this conceptual model has a significant oversimplification, in that it does not consider three-dimensional morphology.
The importance of three-dimensional dune topography on flow was shown by Maddux et al. (2003aMaddux et al. ( , 2003b) ) who undertook a series of flume experiments examining flow over a fixed, artificial, asymmetric, sinuous-crested three-dimensional dune.Maddux et al. (2003a) reported that the flow showed similar characteristics to flow over two-dimensional dunes, including the presence of flow separation and generation of a shear layer.However, the turbulence generated over the three-dimensional dunes was much weaker than over their two-dimensional counterparts, despite the increase in friction.Venditti (2007) extended this analysis by comparing a flat bed, and a two-dimensional asymmetric dune with four three-dimensional dune morphologies, including both lobe-and saddle-shaped crestlines.Venditti (2007) showed that the three-dimensional morphology could either increase or decrease turbulence as compared to its two-dimensional counterpart.Lobe-shaped crestlines increased turbulence, whilst saddle-shaped dunes reduced turbulence by suppressing both the wake structure and mixing within the flow separation cell.Hence, three-dimensionality of the topography is intrinsically linked to the intensity of turbulence (Lefebvre et al., 2016;Maddux et al., 2003aMaddux et al., , 2003b;;Venditti, 2007). 10.1029/2021JF006121 3 of 24 Numerical models offer an alternative approach for investigating the dynamics of turbulent flow over riverbeds.Large Eddy Simulation (LES) has been applied to study flow over dunes, although principally over static, two-dimensional dunes of varying complexity (e.g., Abad, 2008;Frias & Abad, 2013;Nabi et al., 2012;Stoesser et al., 2008;Yue et al., 2006Yue et al., , 2005)).More recent studies have applied LES to consider flow over a single three-dimensional dune (Nabi et al., 2012;Omidyeganeh & Piomelli, 2013a, 2013b) and have captured hairpin-shaped vortices originating at the crestline and upwelling through the water column.Khosronejad and Sotiropoulos (2014) developed a coupled LES-morphodynamic model to simulate the evolution of these flow structures starting from a flat sand bed, and reported that as bedforms developed and grew, spanwise roller vortices were generated that were aligned with the dune crestlines.As the crestlines became deformed in the spanwise direction, these roller vortices evolved into hairpin vortices, which rose to the free surface.
To date, studies of flow over static three-dimensional dunes, in either physical (e.g., Maddux et al., 2003aMaddux et al., , 2003b;;Venditti, 2007) or LES numerical (Omidyeganeh & Piomelli, 2013a, 2013b) models have constructed the morphology by deforming a standard two-dimensional dune profile in the downstream direction according to a sine wave function.Yet, these morphologies and studies are not representative of natural dune fields, as they typically consider individual bedforms and do not consider superimposed bedforms, morphologies with irregular scales of topography (spacing, heights and lengths) or sinuous/discontinuous crestlines.As such, a fully three-dimensional spatially and temporally dependent picture of the turbulent flow field over a train of three-dimensional dunes has yet to be obtained.This is achieved herein through application of a LES model over a static dune field morphology that was generated and measured in a flume under equilibrium flow conditions.This strategy permits examination of the time-dependent three-dimensional flow field over an entire dune field.Here we consider two different three-dimensional morphologies constructed from the same topographic data.For the first Digital Elevation Model (DEM) we extract a long profile from the flume DEM but represent it in three-dimensional space by making the topography homogenous across-stream, thereby forming a pseudo-transverse 2D dune field.For the second DEM, we use the three-dimensional topography created in the flume.This experimental design thereby allows for the first time a quantification of: (a) how the three-dimensional topography governs the mean flow; (b) how topographic steering of the flow characterizes the magnitude and dynamics of the turbulent flow; (c) the time-dependent forces acting on the bed and their implications for bedform evolution; and (d) an assessment of the limitations in applying two-dimensional models to understand flow over natural bedforms.

Methodology
The methodology detailed below comprises three sections: (a) a summary of the physical experiments used in the generation of the dunes and an analysis of their topographic characteristics; (b) description of the numerical scheme, and (c) data analysis procedures that enable investigation of the influence of dune morphological dimensionality upon the flow structure.

Bedform Topography
The topography used herein is taken from the experimental data of Reesink et al. (2018), where a full description of the experimental set-up is provided.In summary, a uniform sand with a D 50 of 239 μm was placed in a recirculating flume 16 m long (L), 2 m wide (W) and 0.5 m deep (H) and allowed to develop under flowing water.The flow conditions used in the experiment had a mean depth-averaged velocity of 0.49 ms −1 and depth of 0.15 m, yielding a subcritical (Froude number = 0.36) and fully turbulent (Reynolds number 1.27 × 10 5 ) flow.Once dunes were formed and in equilibrium with the flow, the flow was stopped without change to the bed, the flume was drained carefully to avoid any disturbance to the bed, and the three-dimensional topography mapped with a Leica Scan Station 2 terrestrial LiDAR.These data were used to construct a DEM with errors <1 mm, which provided the necessary topographic boundary conditions for the numerical model (Figure 1a).In this DEM, there are no dune crests perpendicular to the edges (the flume walls) nor are they continuous across the width of the domain.In order to generate the two-dimensional topography, a profile (y/W = 0.25) was extracted and converted into three-dimensional domains by making the lateral topography homogenous.This two-dimensional pseudo-transverse dune field thus still comprises a greater range of topography than previous two-dimensional experiments in that it possesses a range of dune sizes and superimposed bedforms (Figure 1b).Herein, the complete three-dimensional topography is referred to as the 3D DEM and the two-dimensional pseudo-transverse dune field referred to as 2D DEM.The geometric characteristics of the DEMs were determined through manual analysis of the: (a) maximum and mean bedform wavelength (λ and λ mean ) as estimated from the distance between two crests; (b) maximum crest elevation and trough depression (η c and η t ) (van der Mark, 2008), and (c) standard deviation of bed elevations (Be σ ) (Barabasi & Stanley, 1995).Results are summarized in Table 1 for bed elevation profiles extracted from y/W 0.25 and 0.5, where W is flume width.
The data in Table 1 only provide two-dimensional metrics of the topographic surfaces and herein the three-dimensionality of the topography is of importance.There are fundamental differences between the two DEMs: for example the surface area of the 3D DEM is greater (3D: 3.32 m 2 ; 2D: 3.19 m 2 ) while the standard deviation of bed elevation is lower (3D: 0.018 m; 2D: 0.021 m).To characterize further the difference between the two topographies, we consider the topographic distribution, and calculate both the effective roughness and calculate the lee side angles.Topography was classified into 1 mm bin sizes (Figure 1c) and, although the distribution of topography is similar, the 2D DEM is positively skewed with its central tendency dominated by a peak at 0.17 m.This implies that 2D DEM is rougher, which is also confirmed in the effective roughness calculations.The effective roughness (Hardy et al., 2010) calculates a roughness value that is a scale-dependent measure of the 84th percentile (R 84 ) of the elevation differences between pixels as a function of horizontal scale.To maximize the sample size, a maximum horizontal distance of W/2 was used, which yielded a constant R 84 over a search distance of 0.03 m (Figure 1d).The R 84 value is less for the 3D DEM (0.04 m) compared to the 2D DEM (0.07 m).These differences in roughness characteristics are also shown when the lee side angles are calculated by applying a simple finite difference scheme on the resolution of the spatial discretization, and not over the whole length of the dune (Figure 1e).The mean and standard deviation of the lee side angle for the 3D DEM is 16° ± 9° and for the 2D DEM is 11° ± 9.2°. Lee side angle has previously been shown to be significant in determining the amount of flow separation, with permanent flow separation existing when slope angles are greater than 30°, and intermittent flow separation occurring between 10° and 20° (Best & Kostaschuk, 2002;Cisneros et al., 2020;Kwoll et al., 2016).

The Numerical Model
The LES scheme used herein solves the full three-dimensional Navier-Stokes equations discretized using a finite-volume method.The topography (Figure 1) is incorporated in the numerical model using a Mass Flux Scaling Algorithm (MFSA; Hardy et al., 2005;Lane et al., 2004), with the density of topographic data being equal to the grid density and the numerical grid being defined with vertices that are co-located to the DEM (see Lane et al., 2004).Therefore, the DEM elevations map directly onto grid cells.A five-term MFSA algorithm approach is applied (Miori et al., 2012) and improves the representation of topography by modifying each of the four effective cell face areas at the cell representing the bed-water interface and scales the amount of flux that can pass through each cell.This method interpolates between grid points, and any topographic variability in the bed at a scale smaller than the discretization scale of the topography is not represented explicitly.Therefore, a roughness height, using a non-equilibrium wall function (Launder & Spalding, 1974), is used to represent the effect of subgrid scale topography, with the sides of the domain being frictionless boundaries.In this application we use a 0.005 m spatial resolution for a domain 3 m long and 1 m wide and 0.25 m high (nx = 600, ny = 200, nz = 50).
The improved representation of topography in the model set-up allows topographically induced secondary circulation to be predicted, which has been shown to dominate over both inertial and secondary flows (Blanckaert et al., 2010) as turbulence-driven secondary circulations are typically up to 2% of the primary flow velocities (Nezu & Rodi, 1985).This implies that the flow structures are conditioned largely by the bed topography (van Balen et al., 2010).The sides of the domain are treated with symmetry conditions.
The MFSA is incorporated into a LES model and coupled to a sub-grid scale model (SGS) in order to allow the prediction of time-dependent flow.The SGS model calculates the exchange of energy between the grid and subgrid scales.The energy transfer in turbulent flows is predominantly from larger to smaller eddies (energy cascade), although this process may not be continuous, especially where eddies impinge upon either a solid surface or occur in shear layers.Herein, we use the Wall-Adapting Local Eddy-viscosity model (Nicoud & Ducros, 1999) structured on a tensor invariant that reproduces the proper scaling at the wall.This approach is well suited for complex geometries because: (a) no explicit filtering is needed and only local information is required to build the eddy viscosity; and (b) it is sensitive to both the strain and rotation rate of small turbulent structures (Nicoud & Ducros, 1999).A second order central differencing scheme is used to spatially discretize the governing equations with implicit time-stepping, and the pressure and momentum equations are coupled by applying SIMPLEST, a modification of the SIMPLE algorithm (Patankar & Spalding, 1972).Note.The maximum crest elevation (η c ), maximum trough depression (η t ) and standard deviation of bed elevation (Be σ ; Barabasi & Stanley, 1995) were calculated from the mean profile elevation (blue line and red line, Figure 1a), while the maximum bedform wavelength (λ) and average bedform wavelength (λ mean ) have been measured.All metrics are in meters.

Table 1 Geometric Characteristics of the Digital Elevation Model Used in the Experiments
Initially, a logarithmic flow profile from the experimental data is used for the inlet conditions.After two turnover times recirculating boundary conditions in the streamwise direction are used for pressure and the three components of velocity whereby values from the outlet are mapped onto the inlet.Due to the irregular bed topography, the bed profile was aperiodic.Across the inlet/outlet faces, 91.5% of the area is identical (either fluid or bed) and was mapped directly.The remainder consists of areas of fluid at the outlet, which map onto the bed at the inlet and are thus set to zero (7.5%) or the opposite case (2%) where inlet values are linearly interpolated from the nearest fluid cell at the outlet, assuming zero velocity at the bed.A mass balance check is then applied to the inlet, and a global correction made to ensure mass flow rate is constant.After a period of 10 turnover times to generate the initial turbulent flow, the LES simulations were run at a frequency of 10 Hz for a total of 204.8 s.

Data Analysis
Primary analysis of the flow field is undertaken through Reynolds-averaged statistics; the time-averaged means of each component of velocity; the mean turbulent kinetic energy (TKE) to quantify the energy extracted from the flow by eddies; the Reynolds Shear stress and the boundary layer correlation coefficient (BLCC) to provide a measure of the structural coherence within the flow (Nelson et al., 1993).Next, time series of two-dimensional instantaneous flow maps down the midline (y/W = 0.5) of the dune are shown through maps of normalised instantaneous downstream velocity (u-) and vorticity magnitude (ω) where instantaneous velocity is normalised by the mean flow velocity and the vorticity by the square of the mean flow velocity.Quadrant analysis is then used to classify the turbulent events and understand the generation of Reynolds shear stresses (Wallace, 2016).By applying the definition of Lu and Willmarth (1973), four quadrants are defined by the covariation of the downstream (u') and vertical velocity fluctuations (w') around a zero mean: (a) Q1 events or outward-interactions (+u', +w'); (b) Q2 events or ejections (bursts) (-u', +w'); (c) Q3 events or inward-interactions (-u', -w') and; (d) Q4 events or inrushes (sweeps) (+u', -w').Here, each velocity pair was studied using a range of "hole" sizes (Bennett & Best, 1995;Lu & Willmarth, 1973) of varying multiples of the rms of the velocity components (u'w').A sensitivity analysis was undertaken to quantify the percentage time and the contribution to the total Reynolds stress of events in each quadrant over the whole simulation for increasing hole sizes.Parcels of flow, defined herein as a coherent volume of fluid that contains flow with the same attributes, of Q2 and Q4 events are identified.An Eulerian vortex method is used to detect vortices within an instantaneous snapshot of the flow, by analyzing spatial patterns in the velocity gradient field and its invariants (Green et al., 2007).Herein, we apply the Q criterion, the second invariant of divergence of the velocity vector, which implies that a vortex is present if the magnitude of the vorticity is greater than the strain (Hunt et al., 1988).Each quadrant parcel and vortex structure identified is counted and its geometric characteristics are quantified, by identifying whether the structure is present in a numerical cell, with the smallest size that can be identified being equal to the cell size.The influence that the flow has on the bed is analyzed initially through examining the primary (τ ik = -ρu'w', where ρ is the density of water [1,000 kg m −3 ]) and secondary (τ jk = -ρv'w') Reynolds stresses acting on the bed, and this analysis is extended to consider the contribution through each quadrant.

Reynolds-Averaged Flow
The Reynolds-averaged flow (Figure 2) shows the influence of topography on the flow.Two profiles are reported: the midline (y/W = 0.5) of the 2D DEM and 0.25 y/W of the 3D DEM (Figure 1a, blue line) as they represent the same topographic profile.In the third column (Figures 2c,2f,2i,2l,and 2o), the difference between flows over the 3D-and 2D DEM are reported in relation to the value over the 2D DEM (Difference = (3D-2D)/2D).
The u-components of velocity (Figures 2a and 2b) show regions of faster accelerating flow on stoss slopes and slower than average flow on the lee slopes.For the 2D DEM, the downstream velocity is greater over the stoss slope, but smaller over the lee side and trough, than for the 3D DEM.The difference map between flows shows that the largest differences are in the near-bed region (Figure 2c, label i) where flow over the 2D DEM is greater in magnitude (>125%) than over the 3D DEM.This band of flow is not consistent across the whole profile, but on stoss slopes the u-component is of greater spatial extent and magnitude over the two-dimensional topography (Figure 2c).For the w-component (Figures 2d-2f), positive w-velocities are observed on the stoss slopes, with negative w-velocities downslope of the brink point in the lee side, implying potential flow separation.
Reattachment points can be estimated through both the u-and w-components of flow, where the u-component is equal to zero at the bed and the w-component is negative.This analysis reveals no consistent pattern, and possible reattachment points are present on most lee slopes (Figures 2a and 2b, labeled r p ).However, the difference map of the w-component (Figure 2f) shows the greatest change on stoss slopes (≈200%) with the flow over the 2D DEM being greater in magnitude, although it does not extend as far into the overlying boundary layer (see vertical component of flow, Figure 2f, label ii).This observation is consistent for all stoss slopes.The greatest visual difference between DEMs for all three components of velocity is in the spanwise v-component (Figures 2g-2i).For flow over the 3D DEM, there is a clearly identifiable lateral movement, implying steering of the flow, although there is no discernible pattern as regards a specific cross-stream direction on either the stoss or lee slopes.However, for the v-velocity over the 2D (Figure 2h) there is minimal lateral flow, especially in the near-bed region, and flow here is at least an order of magnitude less than the v-component for the same profile along the 3D DEM (Figure 2g).These differences are clearly apparent in Figure 2i that shows distinct regions of difference, denoting the considerable changes in the lateral steering of flow between a 3D-and 2D DEM.High turbulent kinetic energy (Figures 2j-2l) and Reynolds shear stresses (τ ik ; Figures 2m-o) are detected in the near-bed region.Larger values of TKE are present on the stoss slopes of bedforms and linked to flow in the leeside of the dunes, but appear to form higher magnitude, coherent parcels of high TKE flow, when the stoss slope is sheltered behind a larger lee slope in a topographic hollow (Figure 2j, label iii).For both DEMs, higher magnitude TKE can be observed at 0.2 and 0.5 x/l, although high magnitude TKE is less over the stoss slopes located between 0.65 and 0.8 x/l (Figure 2j, label iii).The difference between the 2D DEM and 3D DEM (Figure 2l) shows that regions of flow in the near-bed region are of much lower magnitude and do not protrude as high into the flow over the 2D DEM, and that flow over the 3D DEM has higher TKE throughout much of the domain.Furthermore, over the 2D DEM, the TKE is greater in the region of skimming flow over some of the smaller dunes (Figure 2l, label iv) although TKE is greater in the leeside region associated with leeside flow for the larger dunes.A similar pattern is evident for the Reynolds shear stress (Figures 2m-2o), with high values on stoss slopes that follow regions of lower magnitude on the preceding lee side (Figure 2m, 0.15 x/l, label v).These Reynolds stresses extend into the flow and can extend over several small topographic protrusions downstream (Figure 2m, 0.2 x/l -0.6 x/l).For flow over the 3D DEM, the regions of high-or low-Reynolds stresses extend from the bed through the boundary layer to the free surface, but this is not the case for the flow over the 2D DEM (Figure 2n).Over the 2D DEM, the regions of high Reynolds stresses are smaller and do not extend to the free surface (Figure 2o) but, similar to the trends in TKE, generate a band of higher magnitude Reynolds stresses close to the bed (Figure 2n).
For the 3D DEM, high boundary layer correlation coefficients (Figures 2p and 2q; R uw > 0.6) start at the dune crest, located in the brink region of the lee slope.This flow moves away from the bed and into the boundary layer, decreasing to moderate values (0.6 > R uw > 0.5, [Nelson et al., 1993]), and by a height of >0.4 z/h this flow appears to coalesce to form large regions of flow with moderate values of BLCC.Downstream of these areas of flow moving away from the bed, R uw attains values of c. 0.3.On the stoss slopes, low BLCCs (R uw < 0.1) are observed and these regions are smaller in their extent than the regions of high BLCC values, and do not move away from the bed.By a height of c. 0.7 z/h, consistently low values of BLCC (R uw < 0.1) are found along the whole profile.A similar pattern is observed for flow over the 2D DEM (Figure 2q), although the regions of high and moderate BLCC are larger and the low-magnitude near-bed regions are not clearly detected along the stoss slopes.

Time Dependent Flow
A series of instantaneous normalised downstream (u-) velocity and vorticity magnitude (ω) plots, each separated by 1 s and starting at 0.7 t/T, are shown in Figure 3 (for animations see Supporting Information S1) for the same topographic profile illustrated in Figure 2.For flow over the 3D DEM, the instantaneous normalised u-component of flow (Figure 3a) shows low-magnitude flow in the topographic lows in the bed, with parcels of low u-magnitude starting in the bed region and advecting towards the free surface (Figure 3a, arrows), and separated by regions of higher magnitude flow over the crestal zone.The size of low-magnitude flow parcels in the topographic hollows varies with time (see region 0.15 x/l in Figure 3a), and parcel size grows until it protrudes into the skimming flow (dashed line, labeled i) and then moves away from the bed, reducing the size of the flow parcel in the bed region (labeled ii).This pattern is not easily detectable in flow over the 2D DEM where flow up the stoss slope and over the crest is faster than the 3D DEM (see Figure 2b, up to 150% of the 3D DEM) and extends to the free surface (see also Figure 5 below).
The normalised vorticity magnitude (Figures 3c and 3d) shows an intense band of vorticity in the near-bed region (≈<0.15z/h) for flow over both the 3D and 2D DEMs.This high vorticity is constant along the whole profile, although the area of high vorticity increases where there are topographic hollows.The normalised vorticity magnitude detects skimming flow at a height determined by the relative protrusion of the largest bedform.Furthermore, as with the u-component, for the 3D DEM (Figure 3c), parcels of high vorticity flow generated near the bed (see region at 0.5 x/l) advect through the entire boundary layer.This is not observed for the 2D DEM where the high-magnitude vorticity in the topographic hollows is supressed closer to the bed (see region at 0.6 x/l, Figure 2a) and the movement of low vorticity fluid away from the bed and through the boundary layer is not as intense.

The Effect of Topography on Three-Dimensional Flow
In order to understand the nature of flow steering over three-dimensional dunes and the mechanism of turbulence generation, flow parcels and vortices are identified using streamlines, quadrant analysis and the Q-criterion.Quadrant analysis is applied to the whole domain for every model discretization for every time step, with the hole size defined by the rms values of the u-and w-components.Initially, we consider the total time and total Reynolds stress contribution from each quadrant with increasing hole sizes.Next, we use quadrants to identify fluid parcels with similar flow characteristics and quantify their frequency and geometric characteristics.We then apply the Q criterion to enable the identification of turbulent structures within the flow.

Streamlines
In order to illustrate the effect of bed topography on the flow, a series of streamlines were calculated from the Reynolds-averaged velocity components starting at the inlet between 0.25 y/w and 0.75 y/w in 0.05 y/w increments and at heights of 0.3 and 0.4 z/h (Figures 4a and 4b).The effect of topographic steering on flow is evident.
For the 3D DEM (Figure 4a), flow is forced, at both heights, into confined flow paths, and this forcing is absent over the 2D DEM (Figure 4b).For the streamlines at 0.3 z/h, convergence occurs at 0.15 x/l whilst for streamlines at 0.4 z/h this occurs at 0.7 x/l (Figure 4a).This flow convergence occurs in the lee side of the dunes where the crestlines are lobate and it is noticeable that, once the flow has converged into these preferential streamlines, there is no evidence that the flow diverges downstream (Figure 4a).

Temporal Occurrence and Reynolds Stress of Quadrants Events
The percentage contribution of time that flow is within each quadrant event has been calculated, to compare the nature of turbulent events between the two DEMs and identify the most appropriate quadrant hole sizes for further analysis (Table 2, Figure 5).When applying a hole size of zero, Q4 events occur most frequently, followed by Q2 and then with Q1 ≅ Q3.When any hole size is applied (Table 2), the sequence changes with more Q2 events detected than Q4 events and with minimal contributions from Q1 and Q3 events.For a hole size of 3 or above (Figure 5a), Q2 events are the only detectable quadrant.In terms of total absolute Reynolds stress (Figure 5b) for a threshold of zero, a sequence of Q2, Q4, and Q1≅Q3 is observed in a ratio of 3:2:1:1.This sequence is consistent for increasing thresholds, although the ratio changes at greater hole sizes.By a hole size of 3, Q1 and Q3 events do not contribute significantly to the total Reynolds stress, although there is still a large contribution of Reynolds stress in quadrants Q2 and Q4.However, Q2 always exceeds Q4.This analysis thus demonstrates both the importance of Q2 and Q4 events, in terms of their temporal occurrence and the nature of the turbulent field, as well as the magnitude of Reynolds stress contributed by Q2 and Q4 events.However, this quadrant analysis does not reveal any significant difference between flow over the 2D -and 3D DEMs.

Quadrant Analysis in Three-Dimensional Space
In order to better interpret the flow structures, quadrant parcels were calculated for every time realisation.In light of the above results, only Q2 and Q4 events are considered herein as they dominate both the proportion of total time and contribute most of the Reynolds shear stress.Furthermore, to aid detection and visualisation (Figures 4c  and 4d) only flow structures exceeding a hole size of 2 are considered.
For both topographies there are more Q2 than Q4 parcels identified.Furthermore, the mean quadrant parcel length (l q ), width (w q ) and height (h q ) are larger for Q2 than Q4 events (Table 3).However, there is no significant difference in the number of Q2 and Q4 parcels identified between the topographies (3D/2D, Q2 = 1.03,Q4 = 1.01).There is also a change in the mean dimensions of these structures between the 2D-and 3D DEMs.
The mean length and height of Q2 parcels over 3D DEM are smaller, while the mean width is greater, than over the 2D DEM.For Q4 structures, all mean dimensions are greater over the 3D DEM.This implies that; i) over the 3D DEM the mean Q2 and Q4 structures are more elongated (anisotropic) than over the 2D DEM; and ii) the Q2 structures are flatter than Q4 structures over the 3D DEM, but this pattern reverses over the 2D DEM.To visualize these structures, one time realisation (0.7 t/T) over both topographies is presented (Figures 4c-4f) and reveals similar features: (a) there are more Q2 flow structures than Q4 (Q2/Q4 ≈ 1.2); (b) the volume of Q2 structures is larger for both topographies at the 95th and 50th percentiles, and (c) the spatial location of these structures can be identified, and for the 3D DEM Q2 events can be observed throughout the flow depth.These events start close to crestlines and move upwards to the free surface (Figure 4c, marked with a circle), while other structures that have become detached from the bed continue to move downstream.
Q4 events are also found throughout the flow depth, but the larger parcels of flow are found in the leeside of dunes in the topographic hollows (Figure 4e, marked with a square).Although this approach of identifying parcels of flow through quadrant detection allows an understanding of the ejection (Q2) or inrush (Q4) events in the flow, it does not allow identification of turbulent flow structures, just volumes of flow with similar attributes.

Q Criterion
Applying the Q criterion enables identification of turbulent structures in the flow, and herein an identical process is applied to that used in identifying the quadrant structures (Figures 4g and 4h).For the whole simulation, there were ∼100% more identified Q-criterion vortices over the 3D DEM as compared to the 2D DEM.In terms of the size of these structures, when the 95th and 50th percentile of each dimension are calculated, the vortices over the 2D DEM were shorter (2D/3D: 95th percentile 0.8; 50th percentile 0.7), narrower (2D/3D: 95th percentile 0.4; 50th percentile 0.5), and thinner in height at the 95th percentile (0.7).
In order to visualize the structures throughout the flow depth, Q criterion isosurfaces (threshold = 3.5) have been cropped so only the middle of the flow profile is shown (0.25 < y/w < 0.75; Figure 6; for animations see Data Availability Statement) with the isosurfaces colored by flow depth.For both the 3D and 2D DEMs, vortices are observed (marked with circles) that have geometric characteristics similar to hairpin vortices.For both DEMs, these vortices originate close to the bed, where a high number of structures are present that rise up into the flow.As these structures rise through the boundary layer, they grow in size, but with the number of structures decreasing and thus indicating the coalescence of vortices.These structures show the geometric characteristics of a hairpin shape, but because of the number of vortices and processes of coalescence/interaction, the geometry of the vortices becomes distorted and results in a series of entangled hairpin vortices.In order to present this pattern more clearly, a randomly identified structure is traced through five images, each 1 s apart (Figure 6), over the 3D DEM, with the same structure being circled as it moves downstream.Initially, the structure possesses a classical hairpin shape (T1), which grows upwards, widens and its legs thicken (T2) until it reaches the free surface (T3).At this point, the structure loses its classical shape as the area of the vortex in contact with the free surface increases (T4), before the structure breaks up and is no longer in contact with the free surface (T5).This is similar to the schematic model of dune-generated macroturbulence interacting with the free surface proposed by Best (2005b).

Reynolds Stress and Quadrant Interactions With the Bed
The primary and secondary absolute Reynolds stress (|τ ik |, |τ jk |) (Figure 7) demonstrates the influence of topographically induced flow on the bed shear stress.For |τ ik | on the 3D DEM (Figure 7a), a high |τ ik | is observed starting at the crest and lee slope of the dunes.This pattern is consistent for both the larger and smaller crest heights (Figure 7a, regions i and ii).High |τ ik | is observed on the stoss slopes (Figure 7a, iii), although this is not observed consistently on all stoss slopes (Figure 7a, iv) and the absolute magnitude of the stress is approximately half that identified on the lee slopes.It is also apparent that a large proportion of the bed experienced minimal |τ ik | throughout the whole simulation.Similar observations can be made over the 2D DEM (Figure 7b) with high |τ ik | identified in both lee slopes and on portions of the stoss slopes.Furthermore, over the 2D DEM, less of the bed area experiences high |τ ik |.Similar spatial patterns are observed when the secondary Reynolds stress (|τ jk |) is analyzed (Figures 7c and 7d).Overall, the area of bed experiencing high Reynolds stress decreases over the 2D DEM, and this is especially marked for the crest and lee side regions.High secondary Reynolds stresses are observed on transverse slopes (Figure 7c, v) demonstrating how the flow is steered over the 3D DEM.Over the 2D DEM, the area of bed that experiences high Reynolds stress is reduced with only three of the crests showing high values.There is also minimal lateral flow over the 2D DEM (Figure 2h), thus reducing the possibility for high secondary Reynolds stresses.
Our numerical data set can also be used to explore the Reynolds stresses present on the leeside slopes for the two bed topographies (Figure 8; see Figure 1e for distribution of leeside angle).Figure 8 reveals the greater magnitude of maximum |τ ik | on the 3D DEM, as well as again illustrating the abundance of leeside angles less than 15° for both beds.For example, on leeside slopes of c. 20°, maximum |τ ik | on the 3D DEM is up to 3 times greater than on the 2D DEM, with this difference being less apparent for leeside slopes <c.7° and greater for slopes >20°.This finding supports past research that has suggested the importance of leeside angle upon Reynolds stresses and flow resistance (Best, 2005a;Best & Kostaschuk, 2002;Kwoll et al., 2016;Naqshband et al., 2014;van der Mark, 2008).
In order to investigate where turbulence is generated in the near-bed region and consider how this would affect an erodible mobile bed, the quadrants acting on the bed were investigated (Figure 9).To achieve this, we calculate the percentage of time that each individual quadrant, at a hole size of zero, is active in the cell above the bed (Figure 9, values displayed for percentages >15%).None of the quadrants, over either topography, exceed a value of 60%, but, as with previous analysis, Q2 and Q4 events show the highest percentage of time interacting with the bed.The spatial distribution of quadrants is more easily identifiable on the 3D DEM (Figure 9, left-hand column).Maximum temporal interactions with the bed for Q1 and Q3 events are located on the stoss slopes.On the 3D DEM, the stoss slopes tend to be transverse as the crestlines are not parallel to the domain sides.This feature is seen consistently, regardless of whether it is located on a large dune or a smaller bedform, and either superimposed on a larger dune or located in front of a crest in a topographic hollow (Figure 9, label i).However, the opposite trend is observed for Q2 and Q4 events, where the highest magnitude events are located on the crest in the brink point region and down the lee slope (Figure 9, label ii).Although this observation is only for one 3D DEM, the majority of high-magnitude temporal interactions with the bed appear to occur where the crestline curvature is lobe-shaped.Similar patterns are present for the 2D DEM (Figure 9, right-hand column) with maximum temporal interactions for Q1 and Q3 located on the stoss slopes and for Q2 and Q4 events located on the lee slopes starting in the brink point region.However, over this 2D DEM, the high-magnitude temporal interactions for Q2 and Q4 events extend across the whole domain, but are not continuous and show a lineated pattern (see Figure 7).
In order to identify which quadrants contribute the greatest total Reynolds stress at the bed and over which topography, we combine the primary shear stress (|τ ik |) analysis with quadrant analysis (Figure 10).For Q1 events areas of high magnitude occur on transverse slopes in topographic hollows (Figure 10, label i) that extend up the stoss slope of the next bedform downstream.Areas dominated by high magnitude Q2 events are identified starting at the crest and down the lee slope (Figure 10, label ii).Q3 events contribute least to the total τ ik, with the only contributions associated with smaller-scale topographic features in the topographic hollows, similar to those observed for Q1.Lastly, Q4 events cover the greatest total area of the bed, show the highest magnitude values, and are identified on the lee slopes, both orientated downstream and transverse to flow (Figure 10, label iii).When the 3D 2D N 1.11 1.09 l q 1.28 1.7 w q 1.19 1.18 h q 1.25 1.44

Table 3
Ratios (Q2/Q4) of the Number of Identified Quadrant Structures (N) and the Ratio of the Geometric Characteristics of Quadrant Parcels of Fluid for Their Length (l q ), Width (w q ) and Height (h q ), for the 2D and 3D DEMs 2D DEM is considered (Figure 9, right-hand column), Q1 and Q3 events contribute on the stoss slope, whereas Q2 and Q4 events contribute on the lee slope.Furthermore, over the 2D DEM, the spatial area of high magnitude |τ ik | contributions is reduced compared to the 3D DEM.

Discussion
This study has examined, through the application of LES, the characteristics of mean and turbulent flow generated over a three-dimensional dune field produced under a steady uniform flow, rather than considering a single three-dimensional dune in isolation.We compare this analysis with a hypothetical two-dimensional 10.1029/2021JF006121 15 of 24 pseudo-transverse dune field in order to enable an improved understanding of how bed morphology determines the mean and turbulent flow and how this affects bed shear stresses.

Flow Characteristics
The Reynolds averaged statistics (Figure 2) and the instantaneous velocity and vorticity time series (Figure 3) show many of the general features reported in most generalised schematic models of dunes (e.g., Best, 2005a).The boundary layer correlation coefficient (R uw ; Figures 2p and 2q) allows identification of different flow characteristics within the flow field.Venditti and Bennett (2000) reported a range of R uw values; for flow separation   it is compressed nearer the bed than for the 3D DEM (Figures 2c, 2l, and 2i).The velocity gradient produces a shear layer that is detectable in the TKE (Figures 2j-2l) and vorticity (Figure 3) plots.The vorticity forms an intense band in the near-bed region and is constant along the whole domain, although increases in magnitude and height above the bed of vorticity are found in the topographic hollows and have previously been hypothesized to comprise Kelvin-Helmholtz instabilities (Abad, 2008;Best, 2005a;Bennett & Best, 1995;Frias & Abad, 2013;Kadota & Nezu, 1999;Nelson et al., 1993;Nezu & Nakagawa, 1993;Müller & Gyr, 1986).Above the shear layer, the flow rises away from bed (0.6 > R uw > 0.5) and moves into the outer boundary layer as a wake flow, with these regions of greater vorticity coalescing >0.4 z/h.The wake zone contains fluid ejected from the shear layer from dunes upstream (Nelson & Smith, 1989), which can be identified in the TKE plots (Figures 2j-2l) where turbulence decreases from the near-bed values, and where the regions of larger Reynolds stresses that are advecting from the bed merge into a wake zone.Above the wake region, R uw = 0.2-0.35,identifies the freestream flow, a quasi-inviscid region affected by the remnants of the upstream flow composed of mild spatial accelerations and wakes from upstream dunes (Nelson & Smith, 1989).On lee slopes, low R uw values (R uw < 0.1) are detected, which implies reattachment of turbulent eddies impacting on the bed (Nelson et al., 1993;Müller & Gyr, 1986).
Other regions of low R uw are located on the stoss slopes and identify the presence of uncorrelated wake turbulence interacting with the boundary layer (McLean et al., 1994).High values of TKE are detected in the near-bed region on the stoss slopes, and form high magnitude, coherent parcels of high TKE flow, especially when the stoss slope is sheltered behind a larger lee slope, within a topographic hollow.Venditti and Bennett (2000) made similar observations that showed TKE reaches a maximum at the reattachment point, and that greater TKE is present within the separation cell.Furthermore, high vorticity is generated in similar topographic locations to the high Reynolds stresses (Figures 2 and 3; Bennett & Best, 1995;Nelson et al., 1993Nelson et al., , 1995;;McLean et al., 1999;van Mierlo & de Ruiter, 1988) and advects towards the bed.Such flow would cause the entrainment and suspension of sediment (Best, 2005a;Cellino & Graf, 2000;Jackson, 1976;Nelson et al., 1995) and therefore influences the dune morphology (Bennett & Best, 1996).
When the velocity components are considered, distinct differences between the three-and two-dimensional topography are evident.For the u-and w-components, flow on the stoss slopes is consistently faster over the 2D DEM, but the w-component, although greater in magnitude, does not extend as high into the flow as the u-component.The greatest difference between the 3D and 2D DEMs is identified in the spanwise (v-) component.Minimal lateral flow occurs over the 2D DEM as flow cannot be steered through topographic lows.Over the 3D DEM (Figure 1a), the crestlines are not at a constant elevation, but instead possess scour holes and topographic lows at the end of crests, similar to the horns of barchan dunes (Bristow et al., 2020), thus allowing flow paths around the lower topography.This is shown (Figure 4a) by both streamlines over the 3D DEM in the near-bed region (<0.25 z/h) and also flow higher in the boundary layer (streamlines at 0.3 & 0.4 z/h) that is drawn down in the leeside, where the crestlines are lobate, and into confined flow paths.This concentration of streamlines is absent over the 2D DEM and the increase in velocity over the 2D DEM thus reduces the size of the leeside separation zone as compared to the 3D DEM (Engel, 1981), with R uw thus reflecting a more intense bounding shear layer between the flows.For the 2D DEM, the reduced lateral exchange and increase in magnitude of the w-velocity, which does not extend as far into the boundary layer, has implications for the use of a two-dimensional understanding of flow in the modeling of suspended sediment.Furthermore, the Froude number for the 2D DEM is slightly lower (≈5%) than for the 3D topography.As such, the morphology may be out of equilibrium with the flow, a common occurrence in natural rivers as flow changes more rapidly than morphology (Martin & Jerolmack, 2013), and the flow structures detected may thus cause the morphology to evolve.

Turbulent Flow Structures
Quadrant analysis with different thresholds hole sizes was conducted for every grid cell and every time iteration in order to calculate the total amount of time and Reynolds stresses that could be attributed to each quadrant.The limitation of the approach is one of resolution, whereby: (a) spatially, it was assumed each grid cell was an independent unit and therefore does not account for the fact that flow structures were larger than the model discretization, and (b) temporally, a quadrant is detected for each model time step rather than the duration of the event.
However, in terms of the total amount of time, Q2 and Q4 events dominate over Q1 and Q3 for both the 3D-and 2D DEMs.As soon as any threshold value is applied, Q2 events are seen to dominate over Q4 events: when a hole size of 2 or greater is applied, Q1 and Q3 events are seen to become insignificant.This is in agreement with Nelson et al. (1995) who report that Q2 events are extremely common and Q1 and Q3 events are relatively rare.The present findings contradict Nelson et al. (1995) in the larger contribution of Q4 events detected herein.However, in terms of continuity, if fluid advects away from the bed then it must be replaced, and previous work suggests that Q2 events are replaced by inrushes of high-velocity fluid from above (e.g., Best, 2005aBest, , 2005b;;Grass, 1971;Grass et al., 1991;Nezu & Nakagawa, 1993), but herein we detect the replacement by lower-velocity fluid.In terms of the quadrant contributions to the Reynolds stress for a threshold of zero, the present results indicate a sequence of Q2, Q4 and then Q1 ≅ Q3, in a ratio of 3:2:1:1.
Quadrant analysis, applying a hole size of 2, is used to detect parcels of flow associated with the topography.For both 2D -and 3D DEMs, there are more Q2 high-magnitude parcels than Q4, and the Q2 parcels are larger.However, over the 3D DEM, a greater number of Q2 and Q4 flow parcels are identified, and they are also more anisotropic in their morphology.The high-magnitude Q2 parcels are smaller over the 3D DEM than the 2D-DEM, but the Q4 parcels are larger.Spatially, Q2 flow parcels originate close to dune crests and move upwards towards the free surface, whilst other flow parcels that have become detached from the bed continue to advect downstream moving towards the free surface.These parcels of large-scale vorticity detected in Q2 events are generated by Kelvin-Helmholtz instabilities (Bennett & Best, 1995), and develop both along the shear layer and also at flow reattachment (Müller & Gyr, 1986).Q4 events are also found throughout the flow depth, but the larger parcels of Q4 flow are found in the leeside of dunes in the topographic hollows.It is suggested that these Q4 events are associated with the interaction of hairpin-shaped vortices with the wall, which have been hypothesized as the main mechanism for destabilizing the bed (Best & Kostaschuk, 2002;Coleman & Melville, 1994;Grass, 1970Grass, , 1971)).
Maps of instantaneous flow and vorticity (Figure 3) illustrate the dynamic interaction of flow with the bed, with parcels of low-momentum flow advecting towards the free surface and separated by regions of higher momentum flow.The Q criterion detects hairpin-shaped vortices (Figures 5e, 5f, and 6) that penetrate the entire boundary layer, originating from the crestlines and rising up towards the surface, similar to that modeled by Omidyeganeh andPiomelli (2013a, 2013b).The multiple scales of topography, over both 2D-and 3D DEMs, generate a range of sizes of vortices and increase turbulence, similar to the findings of Fernandez et al. (2006) who showed that the presence of superimposed bedforms increases levels of turbulent kinetic energy and Reynolds stresses.Here, we show the presence of twice as many hairpin vortices over the 3D DEM than over the 2D DEM, but that these vortices are shorter and narrower than their counterparts generated over a 2D DEM.This finding agrees with Khosronejad and Sotiropoulos (2014) who showed that as the bed develops and becomes more three-dimensional, then the hairpin vortices become more frequent.Furthermore, Khosronejad and Sotiropoulos (2014) suggested that roller vortices evolve into a hairpin-shape and appear in the flow as soon as incipient crestlines emerge on the bed.Because of the number of vortices generated over the fully formed topography herein, it is difficult to confirm this hypothesis of transverse roller-hairpin vortex growth.Moreover, similar to Khosronejad and Sotiropoulos (2014) we show (Figure 6) that vortices grow upwards, coalesce, and their legs are stretched in the streamwise direction by velocity gradients in the local flow and therefore intensify.Because of the multiple crests in the dune field, the hairpin vortices become entangled and merge to form "packets" of hairpin vortices (e.g., Christensen & Adrian, 2001) that may extend to the free surface (Figure 6).These surface upwellings are frequently observed in turbulent open channel flow over dune beds, with the aperiodic and highly intermittent appearance of free surface "boils" (Best, 2005a(Best, , 2005b;;Bradley et al., 2013;Jackson, 1976;Kostaschuk & Church, 1993;Müller & Gyr, 1986).Figure 6 shows that once these flow structures have reached the surface, they break up, decay and move away from the free surface, as described in past conceptual models of dune-related macroturbulence (Best, 2005b).The morphology of hairpin vortices has previously been shown to be dependent on the flow Reynolds number (Head & Bandyopadhyay, 1981), with hairpin vortices becoming stretched and narrower at higher Reynolds numbers: herein, we observe a similar trend when the bed topography becomes three-dimensional.
From this analysis, a new conceptual model of flow over three-dimensional dunes is presented (Figure 11).For both the 3D-and 2D-DEMs, the flow characteristics summarized in Best (2005a) can be observed.Flow separation occurs at the dune crest, with the generation of hairpin vortices that also are detected in the brink region of the lee slope.There is a strong a band of vorticity in the near-bed region forming over the wake zone and region of reattachment along the lee slope.However, our study reveals substantial differences from this idealized model.Compared to flow over the 3D bed, over the 2D DEM, the flow (u-& w-) on the stoss slopes is faster, but parcels of near-bed flow do not extend as high into the flow, which determines the near-bed vorticity, and there is a minimal spanwise (v-) component of flow.For the 3D DEM, the crestlines are not at a constant elevation, and possess scour holes and topographic lows at the end of their crests.This morphology steers the flow through the 3D DEM, while over the 2D bed a shear/skimming flow develops and suppresses turbulence closer to the bed.The topographic steering of flow over the 3D DEM results in greater turbulence and our results show the presence of double the number of hairpin vortices over the 3D DEM.However, these vortices are shorter and narrower than those generated over a 2D DEM, and quadrant analysis reveals the presence of a greater number of Q2 and Q4 flow parcels associated over the 3D DEM but that possess a greater anisotropy in their shape than over the 2D DEM.

Flow Interactions With the Bed and Implications for Sediment Transport
The total Reynolds shear stresses (|τ ik | and |τ jk |) show areas of potential erosion on a mobile bed (Figure 7).For the 3D DEM, a consistent observation, on all crests and lee slopes, is the presence of high total |τ ik| .Although high τ ik is observed on the stoss slopes, it is not consistently observed on all stoss slopes, and the absolute magnitude of the Reynolds shear stress is approximately half that identified on the lee slopes.The absolute magnitude of the secondary Reynolds shear stress is generally lower, especially in the crest and lee side regions but is high on transverse slopes, showing how flow is steered over the 3D DEM (Venditti et al., 2005).In addition a large proportion of the bed experiences minimal absolute magnitude of |τ ik| or | τ jk| throughout the whole simulation.
Quadrant analysis for events detected at the nearest cell to the bed shows distinct regions where events interact with the bed.Q1 and Q3 events are dominant on the transverse stoss slopes, with Q2 and Q4 events dominating the crest/brink point and lee slopes.The percentage of time that quadrant events dominate is greater for Q2 and Q4 events and these events reach a maximum where crestline curvature is at an apex.Such information concerning the flow may aid an understanding of how sediment is transported and how the bed morphology evolves.Nelson et al. (1995) report experiments examining two-dimensional dunes and show that Q2 events contribute positively to the mean bed shear stress and collectively move the majority of the sediment, primarily because they are extremely common.Nelson et al. (1995) also report that Q1 events, although relatively rare, individually move as much sediment as Q4 events of comparable magnitude and duration.Finally, Q4 events, associated with the interaction of Q2 events and hairpin vortices with the flow surface, have been hypothesized as a mechanism for destabilizing the bed (Best & Kostaschuk, 2002;Coleman & Melville, 1994;Grass 1970Grass , 1971).As such, Q1 and Q4 events, produced downstream of points of flow separation, may interact with the sediment on stoss slopes, causing sediment transport up slope to the brink point followed by sediment avalanching down the lee slope.Significant peaks in the spatial pattern of bedload transport have been observed downstream of points of flow separation, thus lending credence to the idea (Best, 2005a;Nelson et al., 1993;McLean et al., 1994;Raudkivi, 1966) that flow separation plays a central role in the development of bedforms.The high percentage of Q2 events associated with the crestline entrain the finer sediment into suspension, and which may dominate sediment transport in sand-bedded rivers (Best, 2005b;Kostaschuk & Church, 1993;Rood & Hickin, 1989;Shugar et al., 2010;Venditti & Bennett, 2000).The frequency and duration of Q2 events (Figure 10) may also highlight how the size, form, and stacking of coherent flow structures into larger flow superstructures may be critical in sediment entrainment and may dictate the relationship between event duration and magnitude that drive sediment impulses at the bed (e.g., Diplas et al., 2008).

Conclusions
The present study has examined the characteristics of turbulent flow generated over a static train of three-dimensional dunes produced under steady flow conditions and compared this to a hypothetical pseudo-transverse two-dimensional topography.We developed a conceptual model (Figure 11) highlighting the differences in flow structure between 2D and 3D topographies that has shown that: 1.The velocity patterns differ between flow over 2D and 3D dunes.Over 2D dunes the magnitude of downstream and vertical flow is greater (u->125%, w-≈ 200%) over the stoss slopes, although the w-component does not extend as far into the boundary layer.Furthermore, there is minimal lateral flow over a 2D bed, which generates a skimming flow and thus reduces the size of the leeside separation zone.2.Over the 3D dunes, 100% more vortices are detected using a Q-criterion analysis, and they are longer, wider and taller than those detected over the 2D dunes.3. Quadrant analysis shows that Q2 and Q4 events dominate the flow and contribute the greatest proportion of Reynolds stresses.However, there is no significant difference between the number of quadrant parcels identified between the 2D and 3D beds, although Q2 parcels (Kelvin-Helmholtz instabilities) are larger in length, width and height than Q4 parcels (hairpin-shaped vortices).Quadrant parcels become more anisotropic in shape over the 3D bed compared to over 2D dunes.4. The coherent flow structures directly influence the bed.Maximum temporal interactions with the bed for Q2 and Q4 events are located on the dune crest in the brink point region and down the lee slope, while Q1 and Q3 events are detected on the stoss slopes.High primary Reynolds stresses (|τ ik |) are observed at the crest and lee slope of the dunes, while high secondary Reynolds stress (|τ jk |) are observed on transverse slopes over the 3D topography.The primary Reynolds stresses are generated by both Q2 and Q4 events and cover the greatest spatial area of the bed.Q1 events influence transverse slopes over the 3D dunes.However, large proportions of the bed experience minimal Reynolds stresses.The high shear stresses that occur for the longest duration are located along crestlines, and highlight how the size, form, and stacking of coherent flow structures into larger flow superstructures are likely critical for sediment transport.These larger-scale flow structures may also determine how the duration and magnitude of turbulent events control impulses of sediment movement at the bed, and that will lead to an increase in the three-dimensionality of bedform morphology.

Figure 1 .
Figure 1.The Digital Elevation Models (DEMs) used as boundary conditions in the numerical models.(a) The three-dimensional morphology collected using terrestrial LiDAR; (b) The 2D DEM from the longitudinal profile shown as red line in (a).The topography is characterized through; (c) the topographic height of each cell in 1 mm bins, (d) the R 84 of the elevation differences between pixels; (e) the distribution of the lee side angles for 3D DEM, and (f) the distribution of the lee side angles for the 2D DEM (f).In c and d, the solid and dashed lines represent three-and two-dimensional topography respectively.

Figure 2 .
Figure 2. The Reynolds-averaged statistics for flow over dunes for (a and b) mean downstream (u-) flow; (d and e) mean vertical (w-) flow; (g and h) mean crossstream (v-) flow; (j and k) turbulent kinetic energy; (m and n) Reynolds shear stress; (p and q) boundary layer correlation coefficient; for the profile down quartile 1 (y/W = 0.25) of the three-dimensional domain (a, d, g, j, m and p) and (b, e, h, k, n and q) midline (y/W = 0.5) of the two-dimensional domain as the profile topography is identical.Column3 (c, f, i, l, and o) represents the percentage difference between the two cases ((3D-2D)/2D); red and blue thus depict areas in which the 3D dunes have increased or decreased values respectively, compared to the 2D dunes, of the given variable.Note a vertical exaggeration of x2 on the ordinate axis.

Figure 3 .
Figure 3.A time series of normalised downstream velocity (u-) and normalised vorticity magnitude (ω) with each image separate by 1 s starting at t/T = 0.7.Panel (a) downstream velocity over the three-dimensional dunes (y/W = 0.25); (b) downstream velocity over the two-dimensional dunes (y/W = 0.5); (c) vorticity magnitude over the three-dimensional dunes; (d) vorticity magnitude over the two-dimensional dunes (y/W = 0.5).The downstream velocity normalised by the mean U while vorticity magnitude normalised by squared U mean.Animations can be found in the Supporting Information S1.

Figure 4 .
Figure 4.The effect of topography on flow.Streamlines (a and b) initiated at inlet between 0.25 y/W and 0.75 y/w in 0.05 y/W increments, at 0.3 z/h (blue) and 0.4 (red) z/h.Three-dimensional quadrant analysis(c-f) where a hole size of two is applied for the threshold for Quadrants 2 (c, d; red) and 4 (e, f; blue) for t/T 0.75.The Q criterion (g and h) colored according to depth for t/T 0.75, where the Q criterion has been masked so only the region between 0.25 < y/W < 0.75 is shown.

Figure 5 .
Figure 5. Quadrant composition for the whole simulation depending on the hole size applied, for (a) time in each quadrant and (b) the total Reynolds stress contribution for the 3D domain (solid lines) and 2D (dashed lines) Digital Elevation Models.

Figure 6 .
Figure 6.A series of Q criterion plots each separated by 1 s in a sub-region of the numerical domain (0.16 < x/l < 0.5, 0.25 < y/w < 0.75).An individual horseshoe vortex is circled, showing its evolution over the three-dimensional topography.

Figure 7 .
Figure 7.The total Reynolds stresses on the bed calculated, for (a and b) Primary Reynolds stress (∑|-ρu'w'|) and (c and d) secondary Reynolds stress (∑|-ρv'w'|) for the (a and c) 3D and (b and d) 2D topography.

Figure 8 .
Figure 8.The relation between lee side angle and total Reynolds stress (|τ ik |) shown through a two-dimensional probability distributions for (a) the 3D DEM and (b) 2D DEM.

Figure 9 .
Figure 9.The total percentage of time that quadrants interact with bed for the 3D DEM (Column 1) and 2D DEM (Column 2).

Figure 11 .
Figure11.A conceptual model of the difference between turbulent flow over a two and three-dimensional dune field.Over both topographies, flow separation (intermittent or permanent) occurs at the dune crest with the generation of hairpin vortices (yellow).Over the lower lee slope (gray region) a band of vorticity is present in the near-bed region (brown).A wake zone associated with flow separation or leeside flow expansion, is formed where turbulent eddies impact on the bed (black solid arrows).Over 2D dunes, flow acceleration over the stoss side, and downstream flow, is greater than over the 3D dunes (as shown by the dark blue shading).Steering of flow by the 3D topography generates stronger lateral-component of flow and more vortices that are more anisotropic in shape over the 3D dunes than over the 2D dunes.