Sand waves in environmental flows: Insights gained by coupling large-eddy simulation with morphodynamics

Sand waves in fluvial, coastal, and Aeolian environments arise as the result of sediment transport processes driven by the coupled interaction between turbulent environmental flow and a sand-covered surface (see Fig. 1). Depending on flow governing parameters (Reynolds number, Froude number, mean flow characteristics, etc.), bathymetric features (e.g., meandering river geometric features), and other factors (e.g., sand properties and sand supply), this interaction can give rise to the emergence of a broad range of organized patterns of varying shapes and sizes. For instance, ripples, transverse dunes, alternate bars, or anti-dunes^1,2 are common bed forms in subaqueous (riverine and coastal) environments where sand supply is plentiful and continuous. Other spectacular formations such as crescent-shaped barchan dunes or star-shaped dunes sculpt desert landscapes when sand supply is limited over the bedrock and their emergence depends largely on the mean wind direction.³ Barchan dunes typically occur in regions where the wind flow is unidirectional whereas star dunes have been associated with 360° large-scale variability in the wind direction.⁴ The scale (amplitude) of these earth features range from cm-scale ripples to tens of meters large mega-dunes in large rivers, while aeolian mega-dunes in deserts can reach heights of up to hundreds of meters. The equilibrium size of these dunes has been shown to be limited by the vertical (bed-normal) extent of the flow domain, i.e., the water depth in rivers and the thickness of the atmospheric boundary layer in Aeolian settings.^5,6 Dunes are complex and dynamic earth features that grow and migrate continuously at speeds (or sand wave celerity) that are very slow (from many hours to years for mega-dunes in rivers and deserts) compared to the integral time scale of the turbulent flow that generates them.^7,8 This continuous migration in conjunction with their large size (relative to characteristic flow length scales) has a profound effect on the turbulent flow and drives landscape evolution. Dunes in rivers, for instance, change streambed roughness, give rise to slowly evolving coherent structures with distinct signatures at the water surface (e.g., surface boils), modify the spectra of turbulence, and provide the primary mechanism for transporting large amounts of sediment during floods, thus, impacting the morphology, streambank stability, and ecology of waterways.^9,10

Most of what we know today about the genesis and long-term evolution of sand dunes has been derived from small-scale laboratory experiments and field observations,^11–17 theoretical studies based on linearized stability theory,^1,18 or simplified computational models based on cellular automata or other macroscopic approaches.^19,20 Understanding gained from experiments is inherently limited by the complexity and scale of the underlying phenomena. Stability analysis studies are valid during early stages of bed deformation when the amplitude of sand waves is very small and as such they are not able to yield insights at later stages of the process. Simple computational models, on the other hand, can produce realistic, albeit overly simplified, emergent landscapes, but they are inherently limited by their lack of coupling the turbulent flow with the motion of the sand bed and the emerging morphodynamics. Consequently, and in spite of considerable progress, many fundamental questions about dunes and their effect on turbulent environmental flows remain unanswered. Among them are questions like: What are the key mechanisms via which the near-bed turbulence destabilizes the initially flat sand bed? How do organized large-scale dunes emerge following this initial turbulence-sand interaction? Why do different shapes of dunes emerge in different environments? How does dune-induced coherent structures modify the structure of turbulence in the ambient flow?

Computational fluid dynamics models coupling eddy-resolving approaches (detached eddy simulation (DES) and large-eddy simulation (LES)) for the flow with bed deformation provide the only feasible approach for tackling such questions.^10,21,22 Previous work, however, has either focused on carrying out LES over fixed dunes^9,23,24 or employed computational techniques that are only capable of resolving small amplitude dunes^21,22 (e.g., the so-called Arbitrary Lagrangian-Eulerian, of ALE, approach). The development of computational models that are able to resolve arbitrarily large and arbitrarily complex dunes, have the resolution required for tackling fundamental questions such as those posed above, and at the same time are also practical enough to tackle morphodynamics at large scales has been a major focus of our work at the St. Anthony Falls Laboratory during the last ten years. The most recent results of our work were presented in the invited talk given by the lead author of this paper during the 67th APS DFD meeting in San Francisco, CA. In this vision paper, we present a much broader summary of our work in this field and outline our vision for the future.

The paper is organized as follows. In Sec. II, we present the macroscopic mathematical equations of bed morphodynamics and discuss approaches for parameterizing the bedload flux, which is a key quantity in these equations. In Sec. III, we present the computational framework we have developed for solving these equations coupled with LES or unsteady Reynolds-averaged Navier-Stokes (URANS) models for the turbulent flow. In Sec. IV, we present new insights we have gained from highly resolved LES into the mechanisms that cause the initial bed instability and lead to the emergence of organized large dunes. In Sec. IV, we also discuss the effects of equilibrium dunes on the turbulent flow both for transverse dunes and barchan dunes. Finally in Sec. V, we summarize our key contributions, present our vision for developing data-driven, site-specific hydro-morphodynamic models of river flooding, and outline key modeling challenges that need to be addressed to realize this vision.

The non-linear interaction of sand particles with a turbulent flow is governed by the coupled interactions of sand grains with each other and the instantaneous turbulent flow.^9,12,25–27 Thus, a first principal mathematical description of the underlying phenomena must employ the Lagrangian equations of motion for each sand grain accounting for the instantaneous hydrodynamic forces imparted by the turbulent flow on each particle and the coupled interactions of particles with other particles, the flow, and the bed. Therefore, such an approach should resolve the flow patterns around each individual sand grain and account for particle collisions and the effect of sand grains on the carrier flow. Various levels of approximations resolving particle-scale phenomena have been proposed in the literature. Escauriaza and Sotiropoulos²⁸ developed a one-way coupled approach, accounting only for the effects of the flow on sand particles, to simulate the interaction of the turbulent horseshoe vortex at a cylindrical pier with sand grains.

A much more sophisticated particle-resolving model was reported by Uhlmann²⁹ who carried out the first two-way coupled direct numerical simulation of sand wave initiation. They employed an immersed boundary (IB) method to resolve flow patterns around individual particles accounting for particle-particle and flow-particle interactions to simulate the early stages of destabilization of an initially flat sand bed in an open channel. Their simulations demonstrated the ability of such a DNS approach to simulate the formation of small-scale ripples but also underscored that such an approach is computationally intractable if the goal is to simulate long-term evolution of bed forms at high Reynolds numbers.^28,29

Simulation of sand wave initiation at realistic scales requires adopting macroscopic approaches, which are based on the continuity equation for bedload sediment. Such equation requires formulating expressions for the bedload flux vector, the rate of entrainment of bedload sediment into suspension, and the rate of deposition of suspended sediment. In what follows, we present the mathematical equations governing the sediment motion and bed deformation and highlight the key computational challenges we have addressed to enable us coupling these equations with the equations governing the turbulent flow.

Macroscopic morphodynamic models employ the Exner-Polya, sediment material mass balance equation,³⁰ to simulate the transient mobile-bed evolution. In its most general form, the Exner-Polya equation for the evolution of bed elevation, z_b, reads as follows:

where

γ is the sediment material porosity, Θ_b is the bedload concentration across the bedload layer, $δ b$ is the bedload layer thickness, Φ is the bedload flux vector, Ω_E is the rate of entrainment of bed material into suspended load, Ω_D is the rate at which suspended load is deposited on the mobile bed (see Fig. 2), u_b is the sediment particle velocity vector within the bedload layer, F_f represents the forces exerted by turbulent flow on the sediment particles within the bedload layer (see Sec. II C), τ₀ is the bed shear stress, ρ_f and ρ_s are the fluid and sediment material density, respectively, and D is the median size of bed material.

A common assumption in most existing computational models is that the variation of bedload concentration in the vertical direction (within the thin bedload layer) is considered to be negligible and Θ_b is assumed to be locally constant across the thickness of the bedload layer. Thus by replacing Θ_b with its vertically average value within the bedload layer, $Θ ¯ b$⁠, Eq. (3) can be rewritten as follows:^21,22

Depending on the turbulence level and sediment material characteristics, one might take the suspended sediment load into account or neglect it in Eq. (1). For most cases, however, sediment material suspension plays a crucial role in bed evolution and thus one needs to compute the entrainment, Ω_E, and deposition, Ω_D, fluxes (for details on how to compute these fluxes see Refs. 31–33). We note, however, that the most critical part of the Eq. (1) for solving mobile bed evolution is to obtain the bedload flux vector, which, among other factors, is a function of sediment particle velocity vector, u_b, and the mean bedload concentration, $Θ ¯ b$⁠. Depending on the method employed to compute bedload flux vector, various coupled hydro-morphodynamics modeling approaches arise, which we will discuss in more details in Sec. II C.

By knowing the mean bedload concentration and sediment particle velocity, one can compute the bedload flux vector (Eq. (6)). The mean bedload concentration, which is in terms of mean flow quantities and bed material characteristics, can be computed using a number of well-known experiment-based equations (see Ref. 34 for more details). Almost all of the existing numerical studies (see, e.g., Refs. 21 and 35–37) have employed some form of such mean-flow based equations in which the mean bedload concentration is obtain as a function of mean-flow characteristics and the Galileo number (G_a), which is a measure of the ratio between gravitational and viscous forces acting on a sediment particles and can be written in the form of a Reynolds number,^34,38

where

where u_g is the gravitational (settling) velocity scale of the sediment particles, T is the effective bed shear stress, and τ_cr is the critical bed shear stress of the sediment material.

Among the most commonly used mean-flow based equation for the bedload concentrations are van Rijn³⁹ and Fredsoe^11,12,40 pickup functions. Major disparity among computational models corresponds to the method they employ to obtain the bedload sediment particle velocity. Various strategies, with different level of sophistication, have been adopted to solve for the sediment particle velocity as we discuss below.

In a Lagrangian approach, the trajectories of individual sediment grains are calculated by solving the particle momentum equations by considering the forces acting on the grain and the interactions with other grains and surrounding structures. The mathematical equations that govern the transport of sediment particles read as follows:²⁸ 

where x and u_b are the position and the velocity vectors of the sediment particle, m is the particle mass, $f → G$ is the total gravitational force, $f → F$ is the flow-particle interactions force that (i.e., drag, lift, added mass, and fluid stress induced forces), $f → P$ is the resultant force due to interactions of the particle with solid walls and other particles, g is the gravitational acceleration vector; C_D, C_L, and C_M are the empirically obtained coefficients for drag, lift, and added mass, respectively, u_r is the particle velocity relative to the local flow velocity u_f (i.e., u_r = u_b − u_f), ω is the vorticity vector (=∇ × u_f), $D D t$ is the Lagrangian time derivative; p is the pressure; and μ is the fluid dynamic viscosity. For more details regarding the formulas and empirical correlations, the reader is referred to Ref. 28. We note that for direct numerical simulation approaches²⁹ that resolve the flow details around sediment particles and their impact on the flow, the lift, drag and added mass coefficients do not need to be modeled as the corresponding forces can be directly computed.

Eqs. (12)–(15) incorporate most of the underlying physics for sediment particle transport. Thus, coupled solution of turbulence and these equations can obtain a physical description of flow and bed deformation. Such an undertaking, however, is not feasible for practical purposes as it would require computing the trajectories and accounting for the interactions of billions of individual sediment particles among themselves and with the flow over time sufficiently long for bed evolution under live-bed conditions to reach equilibrium (see Sec. II C 2). A significantly simpler case to model, on the other hand, would be clear-water scour in which the flow could be treated as dilute mixture (low concentration of particles in the flow) and particle-particle interactions and the impact of particles on the flow can be neglected. Even for such a case, however, using Lagrangian equations to model the initiation and evolution of scour toward equilibrium is not feasible since the computational cost would still be excessive due to the large number of particles and long integration times. A Lagrangian approach could be very effective, however, for gaining insights into the fundamental mechanisms of initiation of motion and sediment transport near hydraulic structures. Escauriaza and Sotiropoulos²⁸ applied a model similar to that described by the Eqs. (12)–(15) to simulate the evolution of 10⁵ inertial particles initially placed on the bed upstream of cylindrical pier mounted in a rectangular open channel driven by flow field in the junction region. As shown in Fig. 3, Escauriaza and Sotiropoulos²⁸ simulated the initial stages of clear-water scour and showed that in accordance to experimental observations the transport of sediment grains is highly intermittent and exhibits essentially all the characteristics of bedload sediment transport observed in experiments, including random bursting events, sliding, saltation, and particle clustering along the bed. They also showed that the bedload flux is characterized by scale-invariance and multifractality (see Ref. 28 for details).

In summary, while Lagrangian models can, at least in principle, incorporate and accurately simulate most of the physics of sediment particle transport they are not practical engineering tools. This is because simulations at Reynolds number of practical interest across the long time scales that characterize sand wave development and evolution are out of reach of present day computing power.^28,29 Such models, however, can provide insight into fundamental processes^28,29 and also be used to develop macroscopic modeling approaches that incorporate the effects of turbulence on bedload sediment transport. The latter was illustrated in the work of Escauriaza and Sotiropoulos²² who employed the Lagrangian momentum equations for sediment grains (Eq. (13)) to propose an Eulerian transport equation for the local particle velocity in the bedload layer. This equation incorporated in its right hand side all the instantaneous forces and was used to calculate the local sediment velocity in the flatbed layer, which is required to calculate the instantaneous bedload flux (see Eq. (6)). Escauriaza and Sotiropoulos²² applied this approach to calculate the early stages of sand wave development in the scour hole developing around a cylindrical bridge pier (see Fig. 4) and demonstrated good agreement with experimental observations and scaling laws previously observed in nature. However, this approach requires very fine mesh resolution in the near-wall layer since the source terms in the particle velocity transport equation involve instantaneous forces that need to be accurately computed at the particle scale. For that Escauriaza and Sotiropoulos²² employed this method in conjunction with an arbitrary Lagrangian Eulerian (ALE) approach to allow for fine near-wall resolution at all times. ALE methods, however, are restricted to simulation of small amplitude sand waves^21,22 and cannot be readily extended to simulate initiation and long term evolution of arbitrarily large and arbitrarily complex dunes. In what follows, we present a macroscopic approach that does not suffer from such limitations and can be used to simulate sand waves in a broad range of complex aquatic environments.

The underlying assumption in this approach is that the sediment particle velocity in the bedload layer depends on the flow field at the edge of the bedload layer. Therefore, the particle momentum equations (either in Lagrangian or Eulerian forms) do not need to be solved and the sediment particle velocity, u_b, is obtained as a function of near-bed flow velocity, shear velocity (u_∗), and the thickness of the bedload layer,^{35,36,41–47}

The flow velocity in such models is reconstructed at the edge of the bedload layer (generally via a wall model) to calculate the sediment particle velocity. Once the sediment particle velocity is calculated, it can be used to obtain the bedload flux vector and solve the sediment mass balance Eq. (1). Most coupled hydro-morphodynamics models use this approach within an ALE framework.^{21,22,36,44,48,49} As mentioned above, however, the amplitude of the sand waves that can be captured with ALE-based models is limited because of the ALE’s intrinsic shortcoming in handling large deformations. In a series of numerical studies, Sotiropoulos and co-workers^{10,17,45–47,50,51} used the flow velocity interpolation method to obtain the sediment particle velocity in conjunction with the curvilinear immersed boundary (CURVIB) method (instead of ALE) to simulate micro- (6 cm high) to macro-scale (1.5 m high) sand waves. In Sec. III, we present an overview of the computational methodology we have developed to enable simulations of arbitrarily complex sand waves and discuss the new physics that this simulation approach has yielded for multi-scale sand waves.

We have developed a novel fully coupled hydro-morphodynamics numerical method that solves the spatially filtered (LES) and/or unsteady Reynolds-averaged continuity and Navier-Stokes equations with the Boussinesq assumption to simulate incompressible, stratified, turbulent flow of dilute water-sediment mixture over evolving mobile sediment beds of arbitrarily complex morphology. The governing equations for the flow field and suspended sediment concentration, ψ, read in curvilinear coordinates $ξ i$ as follows (i, j, k, l = 1, 2, 3):

where J is the Jacobian of the geometric transformation given by $J= ∂ ( ξ 1 , ξ 2 , ξ 3 ) / ∂ ( x 1 , x 2 , x 3 )$⁠, x₃ is the vertical coordinate, $ξ l i =∂ ξ i /∂ x l$ are the transformation metrics, u_i is the ith Cartesian velocity component, $U i = ( ξ m i / J ) u m$ is the contravariant volume flux, $G j k = ξ l j ξ l k$ are the components of the contravariant metric tensor, τ_ij is the subgrid stress tensor for the LES model (dynamic Smagorinsky model^52,53) or the Reynolds stress tensor for the URANS model,^54,g is the gravitational acceleration, δ is the Kronecker delta, $ρ ¯$ is the density of sediment-water mixture (=ρ_f(1 − ψ) + ρ_sψ), W^j is the contravariant volume flux of suspended sediment (⁠$= ( ξ 3 j / J ) w s$⁠), σ^* is the Schmidt number (=0.75),⁴⁹ μ_t is the eddy viscosity, and w_s is the settling velocity of sediment particles.

The governing equations are discretized on a hybrid staggered/non-staggered grid using second-order accurate finite-difference formulas and integrated in time with an efficient fractional-step method (see Refs. 10 and 55 for details). The sharp-interface CURVIB method is employed to simulate the dynamic evolution of the sediment bed and take into account the arbitrary geometric complexity of real-life waterways with embedded hydraulic structures. To simulate the coupling between the flow and morphodynamics, we employ a partitioned fluid-structure interaction approach. We partition the problem into the fluid and bedload sediment domains and solve the governing equations for the stratified flow and the bed morphodynamics separately in each domain accounting for the interaction of the two domains by applying boundary conditions at the sediment/water interface. The dynamically evolving shape of the bed is determined by solving the Exner-Polya equation (see below) coupled with the flow and suspended sediment equations. A mass-conservative sand-slide algorithm is incorporated to account for sudden avalanche events wherever the bed surface slope exceeds the sediment material angle of repose. The numerical model is parallelized using the message passing interface (MPI) communication standard and can be run on parallel high-performance supercomputers. Details of stratified turbulent flow field equations, numerical methods and discretization algorithms can be found in a series of papers form our group.^{10,47,50,55–58} A summary of the key elements of our approach is given below.

Geometrically complex computational domains, like those arising in real-life waterways with migrating sand waves, are handled using the sharp-interface CURVIB method.⁵⁷ The computational domains we consider herein are open channels (straight and/or curved). We employ a background grid that outlines the actual channel geometry and is sufficiently deep to contain the sediment/water interface at all times. The background grid is discretized with a boundary-fitted Cartesian mesh, or curvilinear mesh for curved channels, while the sediment/water interface and channel side walls are discretized with an unstructured triangular grid and treated as sharp interface immersed boundaries in accordance with the CURVIB formulation. The governing equations are solved at the background grid nodes in the fluid phase (fluid nodes) with boundary conditions specified at fluid nodes in the immediate vicinity of the sediment/water interfaces (denoted as IB nodes). The nodes inside the rigid objects are blanked out from the computations. Boundary conditions for the velocity components, turbulence quantities, and sediment concentration are reconstructed at all IB nodes using the wall functions presented above implemented in a manner consistent with the CURVIB method.⁵⁰ A crucial task in immersed boundary method is the classification of the background grid nodes into fluid, IB, and external nodes whenever the shape and/or location of the IB are modified as the evolving sand waves migrate. For this purpose, we employ the efficient ray-tracing algorithm developed by Borazjani et al.⁵⁹ to re-classify the grid nodes every time step as the bed shape deforms.

In the context of the CURVIB method, we discretize the sediment/water interface with an unstructured triangular grid (Fig. 5) and the Exner-Polya sediment mass balance equation (Eq. (1)) is solved on this unstructured mesh using a finite volume method. The bedload layer thickness has been shown to be a function of the flow regime and vary between D to 10D (see Ref. 34 for details). For the flow regimes in this work, we assume that the bedload layer thickness is equal to $δ b =2D$⁠.³⁴ Using the divergence theorem, Eq. (1) can be discretized at a given bed-surface triangular element (Fig. 5) as follows:¹⁰ 

where Δz_b is the bed elevation change over time step Δt and, as schematically illustrated in Fig. 5, A_h is the horizontal projected area of the triangular cell face, the summation is carried out over the three edges of the triangular element (i_e = 1, 2, 3), S is the vector along a cell edge with a length equal to the edge length (⁠$S$⁠), and δ_b is the vertical vector along the z-axis with a length equal to the bedload layer thickness (⁠$δ b$⁠).

We employ a partitioned fluid-structure interaction (FSI) approach⁵⁹ to simulate the coupling between the stratified flow and bed-morphodynamics. We account for the interaction of the two domains by applying boundary conditions at the interface of the sediment and water. Namely, to solve the flow and suspended sediment concentration equations, boundary conditions are specified on the bed in terms of the bed surface location, rate of bed elevation change, and sediment concentration on the bed. On the other hand, to solve the bed elevation equation (Eq. (20)), the near bed stratified flow velocity components, bed shear stress, and sediment concentration from the stratified flow domain are required to calculate: (i) the bedload sediment flux (Φ) and (ii) the net suspended sediment fluxes (Ω_D and Ω_E) at the sediment/water interface (i.e., at the leading edge of bedload layer). For the geophysical simulations we present in this work, the loose-coupling FSI approach is found to be robust enough. For a complete detail of implemented FSI method and how we prescribe the boundary conditions in the two domains reader is referred to Refs. 10 and 50.

The coupled hydrodynamic and bed morphodynamic model might lead to non-physical surface-slopes at the sediment/water interface, i.e., regions on the bed where the calculated surface slope is greater than that of angle of repose of sediment material, φ. In our approach, we examine and correct the calculated bed surface slopes at the end of each hydro-morphodynamics time step. To do so, we use a mass-conservative sand-slide model, which finds the non-physical local slopes and redistributes the sediment mass among neighboring bed cells so that no local bed surface slope exceeds the angle of repose.

As shown in Fig. 6, provided that the angle between center point “P” and any one of the three adjacent cell centers “i” exceeds φ, sediment material is redistributed among neighboring cells to form a new slope equal to the angle of repose (Fig. 6(b)). This is implemented by applying corrections δz_P and δz_i to the bed elevations obtained from the mass balance Eq. (20) at point “P” (z_P) and its ith neighbor (z_i), respectively, for which the bed slope between the two cells whose centroids are δl_Pi apart exceeds the material angle of repose using the following equation:⁵⁰ 

The elevation corrections are computed by balancing the sediment mass among cell “P” and its three adjacent cells (whose horizontal projected areas are equal to A_hP and A_hi, respectively) as follows:⁵⁰ 

Once the Exner-Polya equation (Eq. (20)) is solved, the entire bed is swept to identify elements whose maximum surface slope is larger than the angle of repose. The sand-slide algorithm is then applied only to those identified elements. Since a single sweep is not sufficient to identify and correct all element surfaces whose surface slope exceed the angle of repose, we apply the algorithm until all surface slopes are smaller than the material angle of repose. The iterative process is declared converged when the bed cell steepest surface angle is less than or equal to 0.99 of the angle of repose of the bed material. To eliminate bias in the final solution, the iterative algorithm is implemented through successive sweeps of alternating direction. Usually 2–4 iterations are enough for the sand-slide algorithm to converge.

In this section, we review new findings into the initiation and long term evolution of sand waves derived via numerical simulations using the computational approach we described in Sec. III. We employ LES coupled with the morphodynamics model described above to carry out high fidelity simulation of sand waves in a straight laboratory flume under two conditions: (i) continuous sand supply, via sediment recirculation, which leads to the formation of transverse dunes¹⁰ and (ii) finite sand supply on bedrock, which leads to the formation of barchan dunes.⁶⁰ The specific laboratory flume we have used in our simulations corresponds to that studied experimentally by Venditti and Church.¹³ The dimensions of the flume, sediment properties, and all relevant numerical details can be found in Ref. 10. For quantitative and qualitative validation studies of sand waves simulations using laboratory and field data across scales, the reader is referred to Refs. 10, 17, and 47. In Refs. 17 and 47, we systematically investigated the ability of unsteady RANS and LES-based coupled models to capture the initiation and development of sand waves in different streams and rivers with disparity of scales. Unlike LES, unsteady RANS models are not able to capture the formation of sand waves in small- to intermediate-scale streams.^17,47 However, unsteady RANS models are shown to be able to reproduce migrating dunes in large rivers. We argued that in large rivers, the perturbations imposed by bed movement and the underlying flow structures are presumably sufficiently large for unsteady RANS model to excite the bed instability and give rise to bed forms that are in good agreement with bed forms observed in nature (see Refs. 17 and 47 for details).

Beginning from a flat sand bed, Venditti and Church¹³ and Venditti et al.¹⁵ reported that during the first 150 s of their laboratory experiments, the sand bed was deformed in the following sequence of events (see Fig. 7): (i) spontaneously occurring cross-hatch patterns, which consist of oblique linear striations with amplitude of about 2D and wavelength of about 4–7 cm and an angle of 35°–45° between them, over the entire bed surface; (ii) chevron-shaped features at the nodes of the cross-hatch; (iii) and incipient sand wave crestlines arising as the result of chevron migration. Our numerical simulations¹⁰ captured all the experimentally observed features during the initial stage of sand waves formation. In Fig. 8, we plot contours of the sediment/water interface elevation along the full length of the channel at several instants in time. As shown in this figure, the simulated oblique striations are about 2D in amplitude and about 4 cm in wavelength with angles ranging from 40° to 60° (Fig. 8(a)). At t = 30 s, the cross-hatch patterns have grown in length and amplitude sufficiently to start forming more organized linear striations with small nodes of cross-hatch, which grow in size as the striations migrate downstream to form well-defined chevron-shaped features on the bed (Fig. 8(b)). Subsequently, these chevrons grow and form crestlines that are about 1.5 cm high and 10 cm long (Fig. 8(c)). These processes and computed values for wave amplitude and wavelength are well within the range of the experimental observations of Venditti and Church¹³ (see Ref. 10 for more details).

By conducting two different numerical experiments, we provided compelling evidence to support the role of near-bed turbulence versus interfacial Kelvin-Helmholtz (KH) type waves (i.e., the stratified near-bed sediment/water layer induced by the suspension of sediment material in the flow) as the key mechanism for the spontaneous destabilization of the bed and the onset of cross-hatch patterns.¹⁰ Specifically, we performed two different simulations: (i) one that takes into account the density stratification (i.e., with the Boussinesq term in Eq. (18) activated) in the coupled simulations and (ii) by neglecting the density stratification effects (i.e., with the Boussinesq term in Eq. (18) set equal to zero) but retaining all other aspects of the coupled hydro-morphodynamic model. These simulations showed that during the initial stages of sand wave initiation the two approaches produce essentially identical results, clearly showing that density stratification effects do not play a significant role in initiating the bed instability. At later times, however, the results from the two simulations started growing apart with the stratified flow simulation yielding results in good agreement with the measurements.

These results also support the findings of Venditti and Church¹³ and suggest that near-bed flow structures over the initially flatbed play a key role in destabilizing the bed. To further investigate this hypothesis, we employed our simulated data sets to calculate the fluctuating bed shear stress and also quantify the extend of the initially flatbed occupied by turbulent sweep events. We note that the turbulent sweep events are characterized by fourth-quadrant (Q4) velocity fluctuations with u′ > 0 and w′ < 0, where u′ and w′ are the streamwise and vertical components of flow velocity fluctuations, associated with the interaction of near-wall coherent structures in the turbulent boundary layer with the mobile bed. From the simulated velocity fields, we can readily calculate the velocity fluctuations (see Ref. 10 for details) and detect sweeps by defining the following sweep detection function (SDF), F_s, as follows:

The SDF and fluctuating bed shear stress along with mobile bed features at t = 2 s (when the bed is first destabilized in our simulations) are shown in Fig. 9. It is evident from this figure that both the SDF and the fluctuating bed shear stress exhibit cross-hatch patterns similar in length, spacing, and overall shape to those appearing on the bed elevation contours during the early stage of instability. The length of most of these cross-hatches in both the bed shear stress fluctuations and SDF fields is of the order of ∼6 cm. The angles the lines in the SDF and fluctuating bed shear stress contours form with the axis of channel is about 40°, which is again very similar to that of cross-hatch features on the bed (see Fig. 9). These results make a conclusive case that near-bed coherent structures are responsible for destabilizing the bed and for giving rise to the cross-hatch and chevron shapes during the very early stage of bed evolution. Our numerical experiments have further shown, however, that stratification effects do become important but only later in the process when the initial interaction of sediment and near-bed turbulence has created a stratified layer of suspended sediment above the bed (for a detailed analysis see Ref. 10).

Our simulation results present the first complete description of the initial deformation of the bed but also clarify the fundamental mechanism that leads to organized, large-scale bed forms to ultimately emerge. Turbulent sweeps transport high momentum fluid toward the bed and lead to a spatially heterogeneous fluctuating shear stress field. Concentrated regions of high shear stress fluctuations cause local micro-scour events and associated areas where the eroded material will be deposited nearby. This process is very rapid and occurs spontaneously across the entire bed within the first few seconds of the process as sweep events in the turbulent boundary layer interact with the mobile bed. As shown in Fig. 8, the so resulting bed elevation contours are highly heterogeneous due to the heterogeneity of the bed shear stress fluctuations and the so-resulting spatial variability in the intensity of the local micro-scour and deposition events. As a result the emerging chevron nodes on the bed do not all have the same amount of mass. Smaller nodes tend to travel slightly faster than their somewhat bigger neighbors and ultimately catch up with them merging to form bigger sediment nodes. Fig. 8 clearly illustrates this process and demonstrates how this process biases the bed evolution toward large, organized bed forms to ultimately emerge. Therefore, these results clearly show that the seeds of the emergent structure in sand wave fields are random turbulent flow motions near the initially flatbed.

The specific type of sand waves that will ultimately emerge depends on a variety of factors, including the availability of sand and the characteristics of the turbulent flow that is transporting the sediment. In what follows, we discuss the characteristics of (i) transverse dunes, which emerge when the sediment is continuously recirculated in the flow domain and (ii) barchan dunes, which arise when the supply of sediment is limited over a bedrock. For both dune types, the initiation process we have discussed above is identical and differences become pronounced at later stages when the availability of sediment material starts affecting the emerging dune fields.

To achieve the simulation of transverse dunes, we numerically recirculate the sediment in the flow domain to replicate the sediment recirculation in the laboratory experiment of Venditti and Church.¹³ The details of how this is accomplished can be found in Ref. 10. Following the initial deformation of the bed, organized transverse dunes of complex morphology emerge that span the entire width of the simulated channel. As discussed in detail in Ref. 10, the so emerging dunes migrate at continuously decreasing speeds as they grow in size toward the quasi-equilibrium state. When this state is reached the dune migration speed asymptotes to approximately 5% of the bulk flow velocity and the resulting dune amplitude is established at approximately 35% of the flow depth. Khosronejad and Sotiropoulos¹⁰ presented detailed quantitative comparisons between the measured and calculated temporal evolution of dune amplitude, wavelength, and celerity and reported excellent agreement.

In Fig. 10, we plot a representative snapshot of the simulated bed morphology at t = 700 s during the quasi-equilibrium state of transverse sand waves. The captured bed morphology in this figure shows that the simulation has captured a number of the complex topological features observed in the laboratory, including, among others: the highly distorted and sinuous geometry of the sand waves in the spanwise direction, the complex bifurcations of crest-lines, and the formation of distinct distorted O-shaped features connecting adjacent crest-lines. The latter feature is further highlighted in Fig. 11, where we plot an instant in time when such complex patterns are realized in the simulation and the experiment. The similarity between the two patterns is striking, which further underscores the ability of our simulations to capture sand wave topologies realized experimentally (for further validations see Ref. 10).

In order to show the linkage between flow and the sand waves geometry, we show in Fig. 12 snapshots of an instantaneous iso-surface of the Q-criterion, which are colored with the distance from the bed. As shown in this figure, our simulations capture the emergence of large-scale horseshoe vortices, which appear as soon as the chevrons-like bed features emerge on the mobile bed (see Fig. 12(a)). The appearance of such bed features causes local separation in the lee side of sand waves and the emergence of distorted Kelvin-Helmholtz rollers in the spanwise direction (due to the complex sand wave morphology), which quickly evolve into the horseshoe (or hairpin) structures observed in Fig. 12. The top part of the horseshoes grows upwards while their bottom part is stretched in the streamwise direction by velocity gradients in the local flow and intensifies. As sand waves grow in amplitude, so do the strength and size of horseshoes, which ultimately become strong enough not to be dissipated in the water column and are able to rise all the way to the free surface, causing the distinct “boil” events, which in turn become more frequent and pronounced when sand waves become fully 3D in shape during the quasi-equilibrium stage of bed evolution (for more details see Ref. 10). We note that the role of horseshoes as the key coherent structure responsible for surface boils has long been hypothesized by Best and others^9,61,62 based on field observations. Omidyeganeh and Piomelli^63,64 carried out LES for turbulent flow over frozen dunes and reproduced horseshoes rising to the free surface. Our work, however, presents the first coupled simulation of the emergence and long term evolution of horseshoes starting from the flat sand bed and continuing up until large equilibrium scale dunes emerge and affect the turbulent flow by shedding energetic horseshoes that rise to and interact with the water surface.

A major feature of equilibrium sand dunes is that, as we already mentioned above, they migrate very slowly with an average velocity of about 5% of the mean flow velocity. Such slowly moving sand waves modulate the turbulent flow in the channel at frequencies much lower than the integral time scale of the flow, which is proportional the channel depth and the bulk velocity. Singh et al.⁶⁵ illustrated experimentally that this distinct feature of turbulent flow over mobile beds manifests itself in the flow velocity spectra in the form of a clearly defined spectral gap, which delineates the low-frequency fluctuations associated with the migrating dunes from those of higher-frequency fluctuations associated with flow turbulence. Singh et al.⁶⁵ reported velocity measurements over a gravel bed at only one elevation and hypothesized that the width of spectral gap should be a function of the distance of the velocity sensor from the mobile bed.

In our simulations, we were able to resolve for the first time this important feature of the turbulent flow for a sand bed and demonstrated the existence of a clear spectral gap in the computed velocity spectra (See Fig. 13). Fig. 13 shows computed spectra at three distances above the sand bed and it can be clearly seen that all three spectra exhibit scaling regimes separated by a region of flat spectral gap. The low-frequency f⁻¹ scaling in Fig. 13 is due to migrating sand waves, while the f^−5/3 regions are associated with the turbulent flow in the channel. Fig. 13 clearly shows that the width of the spectral gap increases with distance from the mobile bed. We argue that this is because the right end of the spectral gap (i.e., high-frequency side) is associated with the integral time scale of turbulence (i.e., the scale of the largest turbulent eddies in the flow field). The left (low frequency) end of the spectral gap, on the other hand, corresponds to the velocity fluctuations caused by the smallest migrating sand waves that can be resolved by the velocity sensor. As discussed in Ref. 10, turbulence structures with vertical size smaller than the distance Δ_z between the velocity sensor and the surface, however, cannot produce a signature in the velocity spectra at height Δ_z. Therefore, the smaller the Δ_z the better the velocity sensor is able to resolve velocity fluctuations by the smallest migrating sand waves causing the left end of the spectral gap to move to the right (toward higher frequencies) as the bed is approached. This is clearly evident in the calculated spectra in Fig. 13, which show that, while the right end of the spectral gap is fixed at approximately 1 Hz (the frequency of the largest eddies in the open channel flow), its left end is displaced further to the right into higher frequencies as the distance from the bed is decreased. To our knowledge, our simulations are the first to resolve numerically this important feature of turbulent flows with slowly migrating sand waves.

Barchans are a subcategory of sand waves that are characterized by their distinct crescentic shape, which consists of a windward side, two horns, and a slip face on the lee side⁶⁶ (see Fig. 1). Barchan fields observed in nature consist of many individual barchans spaced at distances that have been shown to scale with either the thickness of the atmospheric boundary layer, for Aeolian barchans, or the flow depth for subaqueous barchans.^6,67

Barchans arise when the flow is unidirectional and the sand supply is limited over the bedrock. In Aeolian environments, they can grow to be hundreds of meters high and migrate at very slow speeds, ranging from centimeters per year for mega-barchans to tens of meters per year for small barchans.^68,69 To date Aeolian barchans in deserts have been studied either by limited satellite observations,^70–73 theory,^7,67,72,74 or simplified computational model that do not account for the impact of turbulence and yield simplified barchan shapes.^75–77 In subaqueous environments, however, scaling arguments show that barchans evolve much faster and scale with the water depth^4,6,7,72 which make the numerical simulation of their genesis and long term evolution tractable from the computational standpoint. For that, we set up a numerical experiment using the same open channel configuration we employed above for the simulation of transverse dunes but with finite sediment supply. More specifically, we prescribed in the simulation a 2 cm thick sand bed and allowed this finite amount of sediment to exit the computational domain without recirculating through the inlet as we did in the simulation of transverse dunes. The sediment/water interface is discretized with an unstructured triangular mesh consisting of over 3.5 × 10⁵ cells and the channel, within which the sediment/water interface is immersed, is discretized with a Cartesian grid, which consists of over 100 × 10⁶ uniformly spaced grid nodes. The minimum distance for the first node off the bed in the background grid system is approximately equal to 16 wall units and the computational time step is set to about one millisecond.⁶⁰ 

Our numerical simulations⁶⁰ reproduced for the first time the entire process of genesis, from a flat sand bed, to long term evolution of equilibrium barchan fields. As discussed above, during the early stages of bed instability the response of the bed is not sensitive to how much sediment is available for transport. Therefore, the same mechanisms described above for transverse dunes are also at work initially in the barchan dune simulation. As the bed organizes to bigger and slower moving sediment masses, however, the limited sediment supply ultimately causes crescent shape features to develop on the bed instead of dunes that span the entire channel width. Field observations^6,72 have shown that when large barchan dunes emerge, transverse sand waves on the windward side of the dunes arise that transport sediment toward the horns. These waves destabilize the horns via a process known as calving, causing new small-scale barchans to be shed from the horns of large parent dunes. These mechanisms, which so far has been observed only via satellite images, have been proposed as the key mechanisms for sustaining equilibrium field scale dunes and preventing all the available sand to ultimate form a single super-large barchan.⁵ Our numerical simulations have been able to capture for the first time the entire process of barchan field evolution. Fig. 14 shows a snapshot of the simulated barchans at equilibrium. The onset of transverse sand waves on the windward side of large dunes, the genesis of small dunes through calving from the horns of the parent dune, and collision and merging of small, faster traveling dunes into larger dunes are all phenomena that are evident in Fig. 14 as well as in video animations of the simulated surface evolution we have created. Most importantly, the dimensions of the simulated dunes scaled with the so-called saturation length scale⁷² are shown to be in excellent agreement with field observations and laboratory experiments. For more details, the reader is referred to Ref. 60. Finally, we note that the barchan dune simulations we present in this work (Fig. 14) are carried out under subaqueous conditions but the present morphodynamics model can be extended to simulate the aeolian barchans shown in Fig. 1 (see Refs. 78–80 for more details).

The computational tools we have developed have enabled us to couple high-fidelity LES with morphodynamic models to simulate sand waves of arbitrarily large size and topological complexity in turbulent flows. The numerical simulations reproduced experimental observations with remarkable accuracy, uncovered new physical phenomena, and elucidated numerous questions of fundamental importance. For example, our results clarified the role of near-bed turbulent sweeps as the mechanism for the experimentally documented spontaneous initial deformation of the sand bed, showed that the seeds of the emergent structure of large dunes reside in the spatial heterogeneity of bed shear stress fluctuations on the flatbed, and elucidated the coherent dynamics and spectral characteristics of turbulent flow over slowly migrating sand waves. Furthermore, we have been able to reproduce for the first time via a coupled simulation the formation of barchan dunes that exhibit all the complexities observed in nature in an environment with unidirectional flow over bedrock with limited supply of sand.

The morphodynamic models we have successfully coupled with LES have been for the most part derived from classical experimental correlations linking near-bed transport processes with mean flow quantities (i.e., bed shear stress). Clearly, therefore, computational models have outpaced experimental work in this field. Further improvements of computational approaches need to be driven by fundamental laboratory scale experiments that will systematically link the instantaneous characteristics of the turbulent flow with near-bed sediment transport, entrainment, and deposition processes.⁸¹ 

Computational models can also provide new insights that will lead to improved macroscopic models of transport. Lagrangian approaches resolving the interactions of sand grains with the turbulent flow, the bed and each other can be a powerful tool in this regard. While such models are still too expensive and impractical for simulating long-term dune evolution at practical Reynold numbers, they can be used to enable simulation-based discovery of new and improved macroscopic correlations linking near-bed turbulence with near-bed transport.^22,29

Most work in this field so far has focused on simulation of single-size sand while in nature sediment material contains multiple classes of sand grains and/or also involves cohesive materials. Incorporating such properties into LES computational models is far from trivial and requires further insights from highly resolved laboratory experiments and fundamental Lagrangian calculations coupling particles with turbulence. This is a major fundamental challenge that we anticipate will be at the center of experimental and computational work in this field.

The computational tools we have developed have also enabled the simulation of mega-dunes in large rivers with complex hydraulic structures.^51,82 At such scale, we have used unsteady Reynolds-averaged models coupled with a technique for decoupling the time steps of the flow from that of the sediment transport to circumvent the large disparity of time scales between the flow and dune migration in nature (note that large dunes can take weeks or even months to migrate in real-life waterways).^17,47 Fig. 15 shows an example from our recent simulations of mega-dunes in a large meandering river with an embedded stream-restoration rock structure. Most of the required computational infrastructure for carrying out such simulations at unprecedented resolution and on a site-specific basis is now in place. Further refining the computational tools to incorporate adaptive mesh refinement techniques, coupling them with site-specific data from multi-beam sonars, LiDAR systems, and satellites, and enabling computer codes to scale to take advantage of emerging extreme massively parallel computing platforms are some of the major challenges for the future. Addressing these challenges will lead to a powerful simulation-based engineering science framework that can be used to virtually create extreme flooding events, understand their potential impact to communities, infrastructure, and economic development, and develop effective flood protection and mitigation strategies to address the impacts of global environmental change.

This work was partially supported by the National Science Foundation Grants EAR-0738726 and IIP-1318201. Computational resources were provided by the University of Minnesota Supercomputing Institute.