1. Introduction
According to global numerical predictions, the frequency of extreme flooding will double across 40 % of the globe by 2050 (Arnell & Gosling Reference Arnell and Gosling2016). With this work, we hope to pave the way for the development of future management technologies that can effectively predict and control these complex environmental events. This is possible only through a thorough understanding of the physical mechanisms governing riverine processes, which is the core objective of our research. Riverine environments are complex systems comprising various phases (Silva-Arya et al. Reference Silva-Arya, Elansary, Mohapatra and Chaudhry2008) whose interactions are interconnected and intricate. Understanding river dynamics requires a comprehensive grasp of sediment-transport mechanisms, which influence the morphodynamics of the system, meaning the spatial and temporal evolution of the riverbed. In recent years, for example, small-scale sediment formations have drawn increasing attention as a close tie with large-scale formations was identified (Best Reference Best1996, Reference Best2005; Venditti et al. Reference Venditti, Church and Bennett2005b ). The most accessible evidence of sediment morphodynamics is the formation of ripples on beaches. These are just one example of the patterns that can arise from the interaction of a fluid and a loose bed (a bed composed of solid granular material). However, a wide variety of large-scale morphological patterns are observed in nature: mega-dunes, ripples, dunes, anti-dunes, bars, chevrons. These bedforms differ in many aspects: size, dimensionality, orientation (Venditti et al. Reference Venditti, Church and Bennett2005a ; Andreotti et al. Reference Andreotti, Claudin, Devauchelle, Durán and Fourrière2011; Coleman & Nikora Reference Coleman and Nikora2011), formation process (Williams & Kemp Reference Williams and Kemp1971; Best Reference Best1992; Blondeaux et al. Reference Blondeaux, Vittori and Mazzuoli2016) and migration patterns (Mohrig & Smith Reference Mohrig and Smith1996; Charru et al. Reference Charru, Andreotti and Claudin2013). Almost all bedforms observed in nature are three-dimensional, even if quasi-two-dimensional bedforms have also been detected in specific flow conditions (Venditti et al. Reference Venditti, Church and Bennett2005a ; Coleman & Nikora Reference Coleman and Nikora2011). The length scales of bedforms vary greatly: from large barrier islands and mega-dunes to dunes and ripples and down to small sand waves. Multiple scales are often present in the same environment. In fact, superimposition of bedforms has been identified as an important aspect of dune evolution and interaction (Galeazzi et al. Reference Galeazzi, Almeida, Mazoca, Best, Freitas, Ianniruberto, Cisneros and Tamura2018). With the term ‘superimposition’ we refer to the simultaneous presence of bedforms characterised by a significant range of scales, e.g. small morphological structures (ripples and dunes) that are hosted by larger-scale bedforms (larger dunes or bars). Reesink et al. (Reference Reesink, Parsons, Ashworth, Best, Hardy, Murphy, McLelland and Unsworth2018) investigated the response and adaptation of alluvial dunes to changing flow conditions typical of flooding events. While the characteristics of large-scale bedforms depend on macroscopic flow parameters, such as flow depth and discharge, the appearance of superimposed bedforms in response to increases in flow depth and velocity was also observed. The authors linked the appearance and disappearance of superimposed bedforms to changing flow rate and depth. They found that an increase in flow discharge or flow depth induces an amplification of the scale range of the bedforms. However, the associated changes in the local flow field were not measured, thus leaving unanswered questions on the nature of the interaction of small-scale bedforms and the flow itself.
In this work, we focus our attention on riverine environments with loose granular beds. Laboratory and field experiments have been extensively performed to examine aspects such as bed-load transport, suspension, settling and bedform evolution. Bedforms can be generated from bed instabilities or pre-existing bed defects (Venditti et al. Reference Venditti, Church and Bennett2005a ). Both analytical (Andreotti et al. Reference Andreotti, Claudin, Devauchelle, Durán and Fourrière2011; Colombini & Stocchino Reference Colombini and Stocchino2011) and experimental (Coleman & Melville Reference Coleman and Melville1996) approaches have been adopted to investigate their formation process. Once formed, bedforms evolve over time and space, growing, migrating and interacting through merging, sediment exchange and linking (Kocurek et al. Reference Kocurek, Ewing and Mohrig2010). These interactions generate complex flow fields (Best Reference Best2005b ).
These smaller structures typically appear over the host dune and, given their smaller size, migrate at a faster speed than the host dune itself. The formation and evolution of large-scale bedform depend on small-scale sand waves, for the mechanism of dune formation (Coleman & Nikora Reference Coleman and Nikora2011) involving near-wall turbulence and localised imperfections and for the mechanism of dune migration (Venditti et al. Reference Venditti, Church and Bennett2005a ). The faster-moving superimposed bedforms transport the sediment particles across the host dune and deposit them on its lee side (downwind-facing portion of the dune). Unsteady flow conditions play a role in the intensity of the turbulent eddies and, therefore, in the host dune composition as well. Furthermore, the presence of superimposed bedforms effectively changes the slope of the stoss side (the upwind-facing portion of the dune) of the host dune, changing the characteristics of the crest (Reesink & Bridge Reference Reesink and Bridge2011). Although these processes might seem secondary and not relevant for the global equilibrium of river systems, some works have shown that the evolution of bedforms is directly linked to the total sediment transport and that the interaction of dunes can cause spikes of sediment flux (Terwisscha van Scheltinga et al. Reference Terwisscha van Scheltinga, Coco, Kleinhans and Friedrich2022). However, Ashley (Reference Ashley1990) questions the significance of bedform superposition, proposing that superimposed bedforms may not fundamentally differ from large-scale bedforms.
Despite the progress made, achieving a complete understanding of the mechanisms that govern bed morphodynamics has been and remains a challenging task from an experimental perspective. This is due to several factors, such as the difficulty of detecting the early stages of pattern formation, the limited amount of flow data that can be collected or the limited extent of the observation time window (Scherer et al. Reference Scherer, Kidanemariam and Uhlmann2020). Additionally, field measurements of the sediment flux during morphogenetic events are difficult to obtain (Garcia Reference Garcia2008).
Numerical simulations are a useful complementary tool to overcome these difficulties and can be used to analyse the interactions between sediment particles and the carrying flow, especially in the near-bed region (Rodi Reference Rodi2017). The great advantage of numerical simulations is the ability to predict quantities such as pressure and vorticity, which are difficult to measure in the laboratory and the field. Lastly, numerical measurements are non-intrusive, as opposed to experimental measurements which often require instruments to be submerged in the fluid. Particle image velocimetry is an experimental technique that provides more detailed information about the flow (two- or three-dimensional measurements) and has no feedback on the flow; however, its implementation is challenging, especially in sediment-laden flows.
In the past, mostly models based on the solution of the Reynolds-averaged Navier–Stokes (RANS) equations have been applied to the study of sediment transport over sand bedforms, highlighting the influence of the bedform geometrical properties on the mean flow characteristics (Parsons et al. Reference Parsons, Walker and Wiggs2004; Lefebvre Reference Lefebvre2019; Yamaguchi et al. Reference Yamaguchi, Giri, Shimizu and Nelson2019). Parsons et al. (Reference Parsons, Walker and Wiggs2004) focused on idealised two-dimensional transverse dunes, while Lefebvre (Reference Lefebvre2019) and Yamaguchi et al. (Reference Yamaguchi, Giri, Shimizu and Nelson2019) studied natural three-dimensional bedforms. More recently, Chiodi et al. (Reference Chiodi, Claudin and Andreotti2014) and Ahadi et al. (Reference Ahadi, Bergstrom and Mazurek2018) performed RANS of sediment transport in open-channel flow using a two-fluid model to relate the bed-load sediment transport to the mean flow quantities. Salimi-Tarazouj et al. (Reference Salimi-Tarazouj, Hsu, Traykovski and Chauchat2024) applied this methodology to study oscillatory flows over fixed and mobile ripples, including both suspended and bed-load sediment transport. Due to the moderate computational requirements, the RANS approach allows one to simulate and study large areas with relatively realistic geometries (Lefebvre Reference Lefebvre2019; Salimi-Tarazouj et al. Reference Salimi-Tarazouj, Hsu, Traykovski and Chauchat2024). However, it fails to explain some of the hydrodynamic and morphodynamic mechanisms due to the parametrisation of both large- and small-scale eddies (Piomelli Reference Piomelli2008; Khosronejad et al. Reference Khosronejad, Kang, Borazjani and Sotiropoulos2011; Coleman et al. Reference Coleman, Garbaruk and Spalart2015; Wu et al. Reference Wu, Soligo, Marchioli, Soldati and Piomelli2017).
An alternative is represented by direct numerical simulations (DNS) of the Navier–Stokes equations, in which all scales of motion are captured by the grid. The increased spatial and temporal accuracy allows the explicit calculation of the distributions of shear stress, sediment concentration and sediment flux that is essential to understand the micro-mechanics of sediment transport, namely the physical processes occurring at the scale of small fluid motions (Kidanemariam & Uhlmann Reference Kidanemariam and Uhlmann2014, Reference Kidanemariam and Uhlmann2017; Akiki & Balachandar Reference Akiki and Balachandar2020; Mazzuoli et al. Reference Mazzuoli, Blondeaux, Vittori, Uhlmann, Simeonov and Calantoni2020; Vittori et al. Reference Vittori, Blondeaux, Mazzuoli and Simeonov2020).
A DNS provides detailed insights but is computationally expensive, often requiring reduced sediment particles or lower Reynolds numbers. Large-eddy simulations (LES) can achieve higher Reynolds numbers but face similar limitations. Direct computation of flow around particles may not be feasible with DNS or LES, depending on the smallest resolved flow scale and sediment grain size. In such cases, a point-particle Lagrangian approach is used, modelling particles as mass points influenced by flow forces (Marchioli et al. Reference Marchioli, Armenio, Salvetti and Soldati2006; Vowinckel Reference Vowinckel2021). A LES with Lagrangian models also requires additional modelling for sub-filter-scale advection velocity (Marchioli Reference Marchioli2017).
The limitations of certain approaches hinder the study of flow configurations with multiple bedform scales. An Eulerian approach, modelling sediment as a continuum, offers lower computational demands than a Lagrangian one. Large-eddy simulations, particularly wall-modelled LES, are commonly used for sediment-transport studies due to their efficiency at high Reynolds numbers (Zedler & Street Reference Zedler and Street2001, Reference Zedler and Street2006; Khosronejad & Sotiropoulos Reference Khosronejad and Sotiropoulos2014, Reference Khosronejad and Sotiropoulos2017). Chou & Fringer (Reference Chou and Fringer2008) employed LES to examine sediment-transport time scales in turbulent flow and later studied sand-ripple formation (Chou & Fringer Reference Chou and Fringer2010), focusing on fixed beds. Khosronejad & Sotiropoulos (Reference Khosronejad and Sotiropoulos2014) developed a hydromorphodynamic model for sand-wave formation in open-channel flow, later exploring Barchan dune formation (Khosronejad & Sotiropoulos Reference Khosronejad and Sotiropoulos2017). Their work, combining fluid and sediment dynamics models, identified key mechanisms in dune formation, though findings are more applicable to desert than riverine environments.
The objective of this work is to investigate superimposed bedforms and advance the understanding of their effects on large-scale morphodynamics and hydrodynamics. Starting from the experimental observations of Reesink et al. (Reference Reesink, Parsons, Ashworth, Best, Hardy, Murphy, McLelland and Unsworth2018), we aim to reproduce numerically the flow over a river bed where the appearance of superimposed bedforms was observed. We simulate a laboratory-scale flow over a realistic river-bed geometry with a coupled hydraulic and morphodynamic system. We analyse the formation process of the bedforms, and the flow dynamics and turbulence characteristics in the presence of such morphological structures. The bed geometry, which is complex and develops over time, is modelled using an immersed-boundary method (IBM). The sediment phase is treated via an Eulerian approach. We analyse the flow features and how they are affected by the presence of superimposed bedforms. From the standpoint of morphodynamics, we analyse the characteristics of the bedforms, the dynamics and interaction with the host dune and their formation mechanism.
In summary, we aim to answer the following research questions: (i) What are the interaction dynamics between superimposed bedforms and the host dune in a coupled hydromorphodynamic system? (ii) How do superimposed bedforms influence the flow dynamics and turbulence characteristics of a laboratory-scale realistic river flow?
In the following, § 2 describes the numerical set-up (governing equations, boundary conditions and simulation parameters). In § 3 we present the predictive results in a realistic application. Lastly, in § 4, concluding remarks and recommendations for future work are also made. In Appendix A, we show the validation procedure of the different methods implemented in this study, namely the IBM, the wall model and the morphological model.
2. Problem formulation
In this section, we summarise the main features of the models presently considered. Firstly, the numerical models for the solution of the fluid phase are presented: these include the spatial and temporal discretisation of the Navier–Stokes equations, the IBM and the wall model. Secondly, we introduce the morphological model implemented for the evolution of the sediment phase.
2.1. Fluid transport
2.1.1. Numerical methods
We solve the Navier–Stokes equations (2.1) using LES. A spatial filter, represented by the mesh itself in our case, is applied to the flow field and all scales smaller than the filter width are modelled. The filtered Navier–Stokes equations read as
\begin{align} &\qquad\qquad\qquad\qquad\qquad\qquad\boldsymbol \nabla \cdot {\tilde {\mathbf u}} =0 \;, \end{align}
\begin{align} &\frac {\partial {\widetilde {\mathbf{u}}} }{\partial t} +\boldsymbol \nabla \cdot \left ( {\widetilde {\mathbf{u}}} {\widetilde {\mathbf{u}}} \right ) = -\frac {1}{\rho }\boldsymbol \nabla {\widetilde {p}} +\boldsymbol \nabla \cdot \left (2\nu {\widetilde {\mathbf{S}}} \right ) -\boldsymbol \nabla \cdot {\widetilde {\boldsymbol{\tau }}} +\boldsymbol \Pi +{\mathbf F}^{\textit{IBM}}+{\mathbf F}^{\textit{WM}} \; . \end{align}
Here we indicate vectors in bold font and scalars in italics. A tilde is used to indicate filtered quantities. The instantaneous filtered velocity is indicated by
${\tilde {\mathbf u}}$
, with
$(\tilde {u},\tilde {v},\tilde {w})$
being the velocity components along the horizontal, vertical and spanwise coordinate directions
$\mathbf{x}=(x,y,z)$
, respectively. The dynamic pressure is
$\tilde {p}$
,
$\rho$
is the fluid density,
$\nu$
is the fluid kinematic viscosity,
$t$
is the time and
$\mathbf{S}$
is the strain-rate tensor:
$\mathbf{S} = (\boldsymbol \nabla \mathbf{u} +\boldsymbol \nabla {\mathbf u}^T )/2$
. Force
$\boldsymbol \Pi$
is the body force used to drive the flow, and adjusted at each time step to maintain a constant flow rate. Terms
${\mathbf F}^{\textit{IBM}}$
and
${\mathbf F}^{\textit{WM}}$
are two forcing terms used to implement the IBM and the wall model, which are discussed in §§ 2.1.2 and 2.1.3, respectively.
The sub-filter-scale stress tensor,
$\boldsymbol{ {\widetilde {\tau }} }= {\widetilde {\mathbf{u}\mathbf{u}}} - {\widetilde {\mathbf{u}}} {\widetilde {\mathbf{u}}}$
, is modelled using an eddy-viscosity approximation,
$\boldsymbol{ {\widetilde {\tau }} }=-2 {\nu _{{}_{\textrm{SFS}}} } {\widetilde {\mathbf S}}$
. To reduce the already significant cost of the calculations, we employed the Smagorinsky model (Smagorinsky Reference Smagorinsky1963) to parametrise the eddy viscosity. This model has been shown to perform adequately far from solid walls and in quasi-equilibrium conditions (Balaras et al. Reference Balaras, Benocci and Piomelli1995). The validation of the fluid-dynamical model, discussed in Appendix A, shows that this choice is sufficiently accurate for the present cases.
Equations (2.1a ) and (2.1b ) are solved on a structured, non-staggered grid by a finite-volume technique (Silva Lopes et al. Reference Silva Lopes, Piomelli and Palma2006). A fractional step method (Kim & Moin Reference Kim and Moin1985) is used, and the solution is advanced in time using the second-order-accurate Adams–Bashforth scheme. The spatial derivatives are discretised employing a combination of a standard second-order central-difference scheme and a quadratic upwind interpolation for convective kinematics (QUICK) scheme (Leonard Reference Leonard1987). The code gradually switches from the upwind scheme near the wall to the central-difference scheme in the outer flow via a hyperbolic-tangent function. More details about this procedure and the reasons behind this choice are given in Appendix A.1. Periodic boundary conditions are applied in the streamwise and spanwise directions and the top boundary is modelled as a fixed free surface:
\begin{equation} \frac {\partial u}{\partial y} = \frac {\partial w}{\partial y} = 0 \; \hbox{ and } \; v=0 \; . \end{equation}
Free-surface fluctuations are neglected, since (i) the Froude number is low (
$Fr=0.37$
), indicating subcritical free-surface deformations; (ii) the primary focus of this paper is to study the influence of superimposed bedforms on river flow, making the behaviour of the free surface of secondary importance; and (iii) Reesink et al. (Reference Reesink, Parsons, Ashworth, Best, Hardy, Murphy, McLelland and Unsworth2018) noted that free-surface deformations did not correlate with the channel bedform morphology. In the remainder of this paper, the tilde is omitted. All spatial averages are intrinsic, i.e. flow quantities are averaged solely within the fluid volume,
$\Omega _f$
, at times when the grid cell is occupied by fluid (Nikora et al. Reference Nikora, Ballio, Coleman and Pokrajac2013). The intrinsic volume average is defined as
\begin{equation} \langle f \rangle = \frac {1}{\Omega _f} \int _{\Omega _f} f(\mathbf{x},t) {d}\Omega = \frac {\int _{\Omega _0}f(\mathbf{x},t)\phi (\mathbf{x},t){d}\Omega } {\int _{{d}\Omega _0}\phi (\mathbf{x},t){d}\Omega } \;, \end{equation}
where
$\phi$
is the volume of fluid,
$\Omega _0$
the total volume and
\begin{equation} \Omega _f =\int _{\Omega _0}\phi (\mathbf{x},t) {d}\Omega \; . \end{equation}
Analogously, the intrinsic spanwise average is defined as
\begin{equation} {\left \langle {f}\right \rangle } _z(x,y,t) =\frac {\int _{W}f(\mathbf{x},t)\phi (\mathbf{x},t)\textrm{d}z} {\int _0^W\phi (\mathbf{x})\textrm{d}z} \;, \end{equation}
where
$W$
is the spanwise length of the domain. The intrinsic time- and spanwise-averaged quantities are computed as
\begin{equation} {\overline {f}} (\mathbf{x})=\frac {\int _{T_0}f(\mathbf{x},t)\phi (\mathbf{x},t)\textrm{d}t} {\int _{T_0}\phi (\mathbf{x},t)\textrm{d}t} \;, \quad {\left \langle {f}\right \rangle } _z(x,y,t) =\frac {\int _0^{W}f(\mathbf{x},t)\phi (\mathbf{x},t)\textrm{d}z} {\int _0^W\phi (\mathbf{x})\textrm{d}z} \;, \end{equation}
where
$T_0$
is the averaging period and
$W$
the width of the domain. The volume-and-time and the spanwise-and-time averages are defined as
\begin{equation} {\left \langle \overline {f} \right \rangle } =\frac {\int _{T_0}\int _{{d}\Omega _0}f(\mathbf{x},t) \phi (\mathbf{x},t){d}\Omega \,\textrm{d}t} {\int _{T_0}\int _{{d}\Omega _0}\phi (\mathbf{x},t){d}\Omega \,\textrm{d}t} \;, \qquad {\left \langle \overline {f} \right \rangle _z} (x,y) =\frac {\int _{T_0}\int _W f(\mathbf{x},t)\phi (\mathbf{x},t)\textrm{d}z\,\textrm{d}t} {\int _{T_0}\int _W\phi (\mathbf{x},t)\textrm{d}z\,\textrm{d}t} \; . \end{equation}
Finally, the average over a cross-section is
\begin{equation} {\left \langle f \right \rangle _{yz}} (x,t)=\frac {\int _{y_1(x)}^{L_y}\int _0^Wf(\mathbf{x},t) \phi (\mathbf{x},t)\textrm{d}y\,\textrm{d}z} {\int _{y_1(x)}^{L_y}\int _0^W\phi (\mathbf{x},t)\textrm{d}y\,\textrm{d}z} \;, \end{equation}
where
$y_1(x)$
is the
$y$
coordinate of the bottom surface of the computational domain, as shown in figure 1(a), and
$L_y$
is the top-surface ycoordinate.
The vector
$\mathbf{u}$
is the instantaneous resolved (filtered) flow velocity, such that
$\overline {\mathbf{u}}$
is the time-averaged velocity,
$\langle \overline {\mathbf{u}} \rangle _z$
is the time- and span-averaged velocity field and
${{\mathbf{u}}^{\prime }} =\mathbf{u}- {\langle \overline {\mathbf{u}} \rangle _z}$
indicates the resolved velocity fluctuations since all flow configurations are statistically homogeneous in the spanwise direction.
(a) Computational grid in the
$xy$
plane. Only every eighth point is shown for clarity. 
, Initial bed shape; 
, zoom area of (b). (b) Close-up view of the bed surface in the area denoted by the magenta rectangle in (a) (
). 
, Fluid nodes; 
, IBM interface nodes; 
, solid nodes. 
denotes the inner-/outer-layer interface, at a distance
$\varDelta _n$
from the IBM surface, and the arrow indicates the velocity vector
${{\mathbf U}^\vartheta _{LES}}$
. (c) Schematic representation of the sediment phase surface and the morphological model.
2.1.2. Immersed boundary method
To treat complex geometries and mobile geometries, we apply an IBM based on the volume-of-fluid approach in the formulation proposed by Scotti (Reference Scotti2006). The flow equations are solved on a fixed grid, and the presence of the solid phase, mobile or fixed, is imposed through the source terms
${\mathbf F}^{\textit{IBM}}$
added to the momentum equation (2.1b
):
\begin{equation} {{\mathbf F}^{\textit{IBM}}}=(1-\phi )\left ( \frac {{{\mathbf u}^{\textit{IBM}}}-\mathbf{u}}{\Delta t} \right ) \;, \end{equation}
where
${{\mathbf u}^{\textit{IBM}}}$
is the velocity of the immersed boundary. The immersed-boundary velocity is set to
\begin{equation} {{\mathbf u}^{\textit{IBM}}}(\mathbf{x})= \begin{cases} {{\mathbf u}^{IS}} \text{ for solid and boundary nodes,}\\ \mathbf{0} \text{ for fluid nodes.} \end{cases} \end{equation}
Solid nodes are identified as the grid cells that are completely below the immersed surface (with
$\phi =0$
); fluid nodes are the set of all grid cells that are completely outside of the immersed surface (with
$\phi =1$
); finally, boundary nodes are the set of all grid cells cut by the immersed surface (with
$0\lt \phi \lt 1$
). The immersed-surface velocity is zero if the geometry is fixed and is equal to the surface velocity if the bed is mobile:
${{\mathbf u}^{IS}}=\textrm{d}{{\mathbf x}^{IS}}/\textrm{d}t$
. Figure 1(a) shows the immersed surface as brown-shaded area. Note that the entire sediment mass can be mobilised if the flow exerts sufficient force. The initial sediment layer, therefore, must be thick enough so that the sediment interface never reaches the computational boundary within the time frame of the simulation. As long as this condition is satisfied, the results are independent of the sediment-layer thickness.
2.1.3. Wall model
Because of the coarseness of the grid, in wall-modelled LES the wall shear stress cannot be obtained directly by differentiating the velocity profile; a model must be provided to link the outer-layer velocity to the wall stress. We use the generalised Moody diagram proposed by Meneveau (Reference Meneveau2020). In this model, it is assumed that the near-wall region is in equilibrium. The flow in the inner layer, under these conditions, is governed by a simplified equation, derived by setting advection and wall-parallel derivatives in the momentum equation to yield
\begin{equation} 0 = -\frac {1}{\rho }\left [\frac {\partial p}{\partial \vartheta }\right ]_{LES} +\frac {\partial }{\partial n} \left [\left (\nu +\nu ^{\textit{WM}}\right ) \frac {\partial \mathbf{u}^{\textit{WM}}}{\partial n}\right ] +\boldsymbol \Pi \; . \end{equation}
Here,
$\hat {\boldsymbol \vartheta }$
and
${\hat {\mathbf n}}$
are the tangent and normal unit vectors to the instantaneous immersed surface,
${{\mathbf S}_{\mathbf b}}$
, and
$\varDelta _n$
is the distance on the inner-/outer-layer interface from the surface. The pressure gradient
$(\partial p/\partial \vartheta )_{LES}$
is taken from the inner-/outer-layer interface and assumed to be constant through the wall layer (consistent with the boundary-layer approximation),
$u^{\textit{WM}}$
is the wall-parallel velocity in the inner layer and
$\nu ^{\textit{WM}}$
a turbulent viscosity that accounts for all the turbulent motions in the inner layer (and which is not the same as
$\nu _{SFS}$
). A mixing length model is used for
$\nu ^{\textit{WM}}$
that includes the effect of surface roughness, and the governing ordinary differential equation (2.11) is solved numerically on an embedded grid (sufficiently fine to resolve the velocity gradient at the wall). The wall stress is then obtained as
\begin{equation} \boldsymbol{{\tau _w^{\textit{WM}}} }=\nu \left (\partial {\mathbf u}^{\textit{WM}}/ \partial n\right )_w . \end{equation}
Boundary conditions are
${\mathbf u}^{\textit{WM}}={\mathbf u}^{\textit{IBM}}$
at the solid/liquid interface and
${\mathbf u}^{\textit{WM}}= {{{\mathbf U}^\vartheta _{LES}}}$
at the inner-/outer-layer interface
$n=\varDelta _n$
(where
${{{\mathbf U}^\vartheta _{LES}}}$
is the outer-layer velocity component parallel to the solid surface, interpolated to
$n=\varDelta _n$
, as shown in figure 1
b).
A fit to the solution is then calculated (analogous to the Moody diagram) from which
$\boldsymbol{\tau }_w^{\textit{WM}}$
can be obtained, knowing
${{\mathbf U}^\vartheta _{LES}}$
,
$(\partial p/\partial \vartheta )_{LES}$
, and the equivalent sand-grain roughness height
$k_s$
, which here is taken to be equal to five times the sand-grain diameter (van Rijn Reference van Rijn1993). Additional details about the model can be found in the original paper (Meneveau Reference Meneveau2020).
The wall stress is then used to generate the two components of a forcing function that is added to the Navier–Stokes equations. The forcing is calculated to ensure that the momentum flux at the solid is equal to the wall stress at the solid/liquid interface:
\begin{equation} {{\mathbf F}^{\textit{WM}}} -\boldsymbol \nabla \cdot \left (\mathbf{u}\mathbf{u}\right )_w +\boldsymbol \nabla \cdot \left [\left (\nu + {\nu _{{}_{\textrm{SFS}}} } \right ) \mathbf{S}\right ]_w = {\tau _w^{\textit{WM}}} . \end{equation}
Note that
${{\mathbf F}^{\textit{WM}}}=0$
at all other fluid points.
2.2. Morphological model
We now discuss the model that describes the time evolution of the channel bed, i.e. the morphodynamics. First, we introduce the conservation equation, and then the model used to represent the sliding of sand that occurs when the surface slope exceeds its repose angle (the angle beyond which spontaneous failure of the slope occurs) (Chien & Wan Reference Chien and Wan1999; García Reference García2008).
2.2.1. Sediment mass conservation
The transport equation for the sediment phase is based on the Exner–Polya equation (Exner Reference Exner1920), revisited by Paola & Voller (Reference Paola and Voller2005). The equation reads as follows:
\begin{equation} \left (1-\gamma \right )\frac {\partial y_b}{\partial t} +\boldsymbol{\nabla }\cdot {{{\mathbf q}_{\mathbf{bl}}}} +\Omega _E-\Omega _D=0 \;, \end{equation}
where
$\gamma$
is the sediment porosity, set to the typical value
$\gamma =0.35$
in this work (Graton & Fraser Reference Graton and Fraser1935; Martin & Aral Reference Martin and Aral1971; Wheatcroft Reference Wheatcroft2002),
${{\mathbf y}_{\mathbf{b}}}$
is the channel-bed elevation,
${{\mathbf q}_{\mathbf{bl}}}$
is the bed-load sediment-flux vector (discussed later in this section),
$\Omega _E$
is the erosion sediment flux and
$\Omega _D$
is the deposition flux. Equation (2.14) represents a sediment mass balance equation for a loose granular bed. Suspended transport entails full lift-off and steady re-entrainment of the sediment particles. This phenomenon takes place only if the particles are small enough to be advected by the turbulent motions. In the present work, we choose a flow–sediment configuration where the particle size and density are not conducive to such a mode of transport, hence bed-load transport is the dominant transport mode. Further discussion of this assumption can be found in Appendix A.3. Therefore, we neglect the suspended-sediment transport, and the balance equation (2.14) assumes the following simpler form:
\begin{equation} \left (1-\gamma \right )\frac {\partial y_b}{\partial t} +\boldsymbol{\nabla }\cdot {{{\mathbf q}_{\boldsymbol{bl}}}} =0 \; . \end{equation}
Equation (2.15) is discretised using a central-difference scheme and advanced in time with the Adams–Bashforth scheme.
Mesh representation of a sample bedform field: flow dynamics mesh slices are shown (black Cartesian mesh; every other point is shown for clarity). The inset shows a zoom on the region
$9\lt x/H\lt 10$
and
$0\lt z/H\lt 0.5$
. The triangular mesh of the bed surface is coloured by the vertical distance from a reference level
$\eta _b=(y_b-y_{ref})/(y_{max}-y_{ref})$
, with
$y_{ref}/H=0.4$
(every point is shown).
The immersed surface is discretised by triangular elements connecting the surface nodes, and one surface node for each grid cell is defined. Figure 2 shows a representation of the mesh used to discretise the bedform field developed in an open-channel flow (see Appendix A.3), overlapped to slices of the Cartesian flow mesh. The position of the node is advanced in time employing equation (2.15), and subsequently the sand-slide algorithm is applied, as described in the following section.
The bed-load transport process is influenced by the balance of forces on settled particles, where transport occurs when eroding forces exceed settling forces. The critical bottom shear stress for initiation of motion is determined by factors such as particle Reynolds number and size (Shields Reference Shields1936; van Rijn Reference van Rijn1984a
; Nielsen Reference Nielsen1992; Niño et al. Reference Niño, Lopez and García2003; Debnath & Chaudhuri Reference Debnath and Chaudhuri2010). We adopt an empirical model for the bed-load transport rate. The vector
${{\mathbf q}_{\mathbf{bl}}}$
is tangent to the channel bed with non-null components in the local streamwise and spanwise directions (van Rijn Reference van Rijn1987; García Reference García2008). It is defined as the product of sediment volumetric concentration
$C_{bl}$
, sediment velocity
${{\mathbf U}^\vartheta _{\mathbf{bl}}}$
and bed-load layer thickness
$\delta _{bl}$
:
\begin{equation} {{{\mathbf q}_{\mathbf{bl}}}} ={\delta _{bl}} {C_{bl}} {{{\mathbf U}^\vartheta _{\mathbf{bl}}}} \; . \end{equation}
We set the bed-load layer thickness to twice the particle diameter,
$d$
, as suggested in the literature (Charru et al. Reference Charru, Andreotti and Claudin2013; Khosronejad et al. Reference Khosronejad, Kang, Borazjani and Sotiropoulos2011). Notice that the immersed surface represents the location of the settled particles, while the edge of the bed-load layer is a fictitious location where we extrapolate the velocity
${{\mathbf U}^\vartheta _{\mathbf{bl}}}$
(see figure 1
c). The velocity vector of the sediment particles within the bed-load layer,
${{\mathbf U}^\vartheta _{\mathbf{bl}}}$
, is approximated as the resolved velocity parallel to the bed surface, interpolated to the edge of the bed-load layer. The concentration in the bed-load layer is given by the empirical equation suggested by van Rijn (Reference van Rijn1987):
\begin{equation} {C_{bl}} =0.18C_0\frac {T}{d_*} \;, \end{equation}
where
$C_0=0.65$
is the maximum sediment concentration and
$d_*$
and
$T$
are defined as
\begin{align} d_*=d\left [\frac {(s-1)g}{\nu ^2}\right ]^{1/3} \;, \quad T=\frac {\theta -\theta _{cr}}{\theta _{cr}} \;, \quad \theta =\frac {|\tau _w|}{(s-1)\rho gd}=\frac {u_\tau ^2}{(s-1)gd} \;, \end{align}
where
$s=\rho _s/\rho$
is the ratio between particle density,
$\rho _s$
, and fluid density,
$\rho$
, and
$g$
is the acceleration of gravity. The Shields parameter,
$\theta$
, represents the bottom shear stress normalised by the specific gravity and grain diameter of the sediments. In this formulation, the threshold shear stress is
$\theta _{cr}$
, and the critical Shields parameter is given by
\begin{equation} \theta _{cr}=\dfrac {\tau _{w,cr}}{(s-1)\rho g d} \,, \end{equation}
where
$\tau _{w,cr}$
is defined according to the empirical formula suggested by van Rijn (Reference van Rijn1993), and corrected to take into account the influence of gravity on erosion from sloping beds with the following correction factor:
\begin{equation} f_{\theta }=\frac {\sin (\alpha +\beta )}{\sin (\alpha )} \,, \end{equation}
where
$\beta$
represents the bed slope angle.
The system formed by equations (2.16)–(2.19) provides the bed-load transport-rate vector, using the wall shear stress and the extrapolated bed-load-layer velocity as inputs. Then, by computing its gradient, we can determine the bed elevation displacement through (2.15). Finally, a further step is required to ensure that the computed bed elevation does not exceed the threshold for landslides, as explained in the following section.
2.2.2. Sand-slide algorithm
The sand-slide algorithm used in this work is similar to that applied by Khosronejad et al. (Reference Khosronejad, Kang, Borazjani and Sotiropoulos2011). Whenever the angle between two adjacent immersed-surface nodes is larger than the angle of repose,
$\alpha$
, the elevation of the node and its neighbours is adjusted so that the local volume of sediment is conserved and all slopes are lower than the repose angle. The procedure is iterative, since the avalanche can cause neighbouring nodes to exceed the angle of repose, and the iteration stops when all steep slopes have been eliminated. The algorithm takes into account only the action of gravity, and neglects any dependence on the flow conditions (Beakawi Al-Hashemi & Baghabra Al-Amoudi Reference Beakawi Al-Hashemi and Baghabra Al-Amoudi2018). Experimental studies show that the difference between the static and dynamic angle of repose is less than 10 % in our range of sediment characteristics (Cheng & Zhao Reference Cheng and Zhao2017); therefore we neglect this difference. To avoid an accumulation of error, each iteration starts from a different corner of the domain. The system of equations defining the algorithm is
\begin{align} & \frac {c_j\left (y^*_{b,P} +\Delta y_P\right )-\left (c_jy^*_{b,j} +\Delta y_j\right )}{\Delta l_{Pj}} =c_j\tan (\alpha ) \quad \text{ for } j=E,W,B,F,\notag\\ & A_P\Delta y_P+\sum _{j=E,W,B,F} A_j\Delta y_j=0 \;, \end{align}
with
\begin{equation} c_j=\begin{cases} 1\quad \text{ if } {\left | {\left (y^*_{b,P} -y^*_{b,j}\right )/\Delta l_{Pj}} \right |}\gt \tan (\alpha ), \\ 0\quad \text{ otherwise,} \end{cases} \end{equation}
where
$y^*_{b,P}$
and
$y^*_{b,j}$
are the bed elevation values predicted by the morphological model at point P and its neighbour
$j$
, respectively. We use four neighbouring points to the east (E), west (W), back (B) and front (F) of each point (P). Here
$\Delta y_P$
and
$\Delta y_j$
represent the changes in elevation to be applied to ensure that all angles are lower than
$\alpha$
. Distance
$\Delta l_{Pj}$
is the distance between point
$P$
and its neighbour
$j$
. Therefore, equation (2.21) is a system of five equations in five unknowns
$(\Delta y_P,\Delta y_E,\Delta y_W,\Delta y_B,\Delta y_F)$
to be solved for each node of the immersed surface.
2.3. Summary of the numerical model
In this section, we have presented all the techniques required for the solution of the flow dynamics and morphodynamics of the selected flow configurations. To summarise, the coupled fluid–sediment system is advanced in time as follows:
-
(i) Given a bed geometry, the volume-of-fluid fraction,
$\phi$
, and the IBM forcing term,
${{\mathbf F}^{\textit{IBM}}}$
, are calculated from (2.9). -
(ii) The wall shear stress,
$\boldsymbol{\tau _w^{\textit{WM}}}$
, is calculated from (2.12). -
(iii) The wall-model forcing term,
${{\mathbf F}^{\textit{WM}}}$
, is calculated from (2.13). -
(iv) The flow field is updated using the information from the two previous time steps (
$\mathbf{u}$
,
$p$
,
$ {\nu _{{}_{\textrm{SFS}}} }$
,
${{\mathbf F}^{\textit{IBM}}}$
,
${{\mathbf F}^{\textit{WM}}}$
) using (2.1). -
(v) The bed-load sediment transport,
${{\mathbf q}_{\mathbf{bl}}}$
, is calculated using
$\boldsymbol{\tau _w^{\textit{WM}}}$
and (2.16). -
(vi) The bed elevation,
${{\mathbf y}_{\mathbf{b}}}$
, is updated using (2.15). -
(vii) The slope angle is adjusted using (2.21).
When the bed geometry is fixed the first step is performed only once at the beginning of the simulation, and the last three steps are not performed. The full flow field and the bedform field from this procedure are collected for analysis, the results of which are presented in § 3.
2.4. Simulation set-up
In § 3, we show the results of wall-modelled LES over a realistic bed geometry extracted from the laboratory measurements of Reesink et al. (Reference Reesink, Parsons, Ashworth, Best, Hardy, Murphy, McLelland and Unsworth2018), who considered an open-channel flow where the bottom surface is composed of loose granular material and the bed is initially flat. They analysed the response of the channel bed to flow variations, such as changes in flow velocity and water depth. The flow conditions were varied suddenly and the bed morphology was recorded for a prolonged time (
$1.5$
hours). Reesink et al. (Reference Reesink, Parsons, Ashworth, Best, Hardy, Murphy, McLelland and Unsworth2018) tested 23 flow configurations with the bulk Reynolds number varying between
$90\,000$
and
$135\,000$
and the Froude number varying between
$0.3$
and
$0.5$
. They could vary the flow conditions by increasing or decreasing both the discharge and the water depth. The study revealed different morphological changes in response to variations in flow depth and velocity: (i) an increase in bedform superimposition following an increase in depth and velocity, (ii) the flattening of dunes in response to decreased flow depth and (iii) an increase in scour, the removal of sediment from the bed, in response to increased flow velocity. We selected this study because of the emergence of superimposed bedforms and the desire to understand the physics behind this multiscale phenomenon. Superimposition occurred in nine (out of 23) cases. We chose one configuration at a bulk Reynolds number
$Re_b=\widehat {u_b}\widehat {H}/\nu =1.25\times 10^{5}$
and Froude number
$Fr=u_b/\sqrt {gH}=0.37$
, where
$\widehat {H}$
represents the bulk channel depth and
$\widehat {u_b}$
the bulk velocity (defined thoroughly in the following).
Some differences exist between the laboratory experiment and the numerical simulation presented here. First, we extract a section of the measured centreline bed profile and extrude it in the spanwise direction, thus generating a bed geometry that is initially two-dimensional, shown as a shaded brown area in figure 1(a). This choice is dictated by the fact that no data are available for the spanwise shape of the bedform. However, since the flow field is highly three-dimensional and turbulent, the bed starts to develop three-dimensional features immediately.
Second, we apply periodic boundary conditions in the streamwise and spanwise directions for both the fluid and sediment phases, while the laboratory flume is bounded laterally by solid walls, and the inlet conditions are imposed. Following the studies by Ojha & Mazumder (Reference Ojha and Mazumder2008), Omidyeganeh & Piomelli (Reference Omidyeganeh and Piomelli2011, Reference Omidyeganeh and Piomelli2013) and Hardy et al. (Reference Hardy, Best, Marjoribanks, Parsons and Ashworth2021), we chose the domain length and width to be large enough to guarantee that the fluid structures can form and develop naturally. Within the time frame of the present calculation the main dunes move only a small fraction of the domain’s total length, which justifies the use of a single dune, rather than a sequence.
The simulation is initialised by integrating the Navier–Stokes equations over the fixed initial geometry until a statistically steady state is reached. Then, at time
$t=0$
, we allow the bed to deform and observe its subsequent development. Table 1 shows the sediment and flow parameters used. This choice neglects the previous history of the flow effects on the bedform, but cannot be overcome, unless the experimental measurements provided the three-dimensional bedform as a function of time, at closely spaced time instants. The simulation is run at a constant flow rate to resemble the experimental set-up.
Flow, roughness and sediment parameters. Here
$h$
is the dune crest height,
$H(x)$
is the local channel depth and the superscript ‘crest’ indicates quantities evaluated at the crest location.
The cross-sectional average of channel depth, the volume-averaged channel depth and the cross-sectional average of the velocity (i.e. the bulk velocity) are defined as follows:
\begin{align} &H(x,t)=\frac {1}{W}\int _0^W\int _{y_1(x)}^{L_y}\phi (\mathbf{x},t)\textrm{d}y\,\textrm{d}z \;, \end{align}
\begin{align} &\qquad\qquad\widehat {H}=\frac {\int _{V_0}\phi (\mathbf{x},t)\textrm{d}V}{WL} \;, \end{align}
\begin{align} &\qquad\qquad u_b(x,t)= {\left \langle u(\mathbf{x},t) \right \rangle _{yz}} \;, \end{align}
\begin{align} &\qquad\quad\widehat {u}_b= {\left \langle {u(\mathbf{x},t)}\right \rangle } =\hbox{constant} . \end{align}
In this flow configuration
$u_b(x,t)$
and
$H(x,t)$
vary along the streamwise direction and in time (if the bed is mobile); however, the product
$u_bH$
is constant at all cross-sections and all times and equal to
$\widehat {u}_b\widehat {H}$
. Hereinafter,
$\widehat {u}_b$
will be used as reference velocity for normalisation. If the geometry is fixed, the cross-sectional averages of the velocity and the channel depth are only a function of the streamwise direction,
$x$
. As a consequence, the bulk Reynolds number,
$Re_b$
, shown in table 1, is a constant global parameter, while the Froude number,
$Fr$
, varies such that its maximum is found over the dune crest. The flow is fully turbulent and subcritical (
$Fr\lt 1$
). We use the friction velocity over the dune crest to evaluate a posteriori the equivalent sand-grain roughness, which is
$k_s^+=53$
, in wall units, hence in the transitional regime (Jiménez Reference Jiménez2004; Schultz & Flack Reference Schultz and Flack2007).
The grid has a region with fine wall-normal spacing near the surface, then is stretched to the top surface. The number of points in the fine-resolution region is chosen so that the morphology variation in time is well captured by the grid. The domain has length
$L/h=13$
and width
$W/h=8$
and height in the range
$1.7\lt L_y/h\lt 2.8$
, where
$h$
is the height of the larger crest of the initial bed geometry at the initial time.
Since no flow measurements are available for the geometry considered in this study, we performed a careful validation study, which is described in Appendix A. We used several reference cases to validate the various parts of the numerical model. Since the Reynolds number of the actual dune is higher than in any of the validation cases, the near-wall flow is closer to equilibrium than in the test cases, and we expect the error to be no greater than that observed in Appendix A.3.
2.5. Grid-convergence study
We performed a grid-convergence study with a fixed bed to determine the resolution needed to capture accurately the flow dynamics. Table 2 summarises the numerical parameters of the grids used. To decouple the wall-modelling and resolution errors (Kawai & Larsson Reference Kawai and Larsson2012), the grid-resolution study was performed with a fixed wall-model interface location,
$\varDelta _n=0.05h$
.
Grid parameters in wall units. Wall units are computed using the value of the time-averaged friction velocity at the crest (
$x=0$
). Both the minimum value (near the immersed surface) and the maximum value (near the top surface) of the spanwise grid spacing,
$\Delta y^+$
, are reported.
Figure 3 shows the streamwise distribution of the time- and span-averaged friction coefficient:
\begin{equation} {\left \langle \overline {C_f} \right \rangle _z} =\frac { {\left \langle \overline {\tau _w} \right \rangle _z} }{1/2\rho \widehat {u}_b^2} \; . \end{equation}
Note that the
$C_f$
predicted by the wall model includes the contribution of the form drag of the roughness elements, through the rough-wall modification of the generalised Moody diagram.
Figure 3 shows the presence of two recirculation regions behind the two crests at
$x/h=0$
and
$x/h=5$
. After separation, the flow reattaches to the bottom boundary at
$x/h\simeq 3$
over the smaller dune, and at
$x/h\simeq 7$
over the larger dune, with recirculation bubbles of average length
${\overline {l_1}} /h=2.002$
and
${\overline {l_2}} /h=1.012$
, respectively. Once the flow reattaches, the friction coefficient increases over the stoss side of both dunes (
$3\lt x/h\lt 5$
and
$7\lt x/h\lt 13$
). The three grids give similar results. The integral difference between the skin friction coefficients on the medium and fine grids is 4 % of the maximum skin friction coefficient on the fine grid.
Time- and spanwise-averaged friction coefficient,
${\langle \overline {C_f} \rangle _z}$
. 
, Coarse grid; 
, medium grid; 
, fine grid.
Vertical profiles of the time- and span-averaged (a) streamwise velocity
${\langle \overline {u} \rangle _z} /\widehat {u_b}$
, (b) Reynolds shear stress
${\langle \overline {{{u}^{\prime }} {{v}^{\prime }} } \rangle _z} /\widehat {u_b}^2$
and (c) turbulent kinetic energy
$2{\langle \overline {\mathcal{K}} \rangle _z} /\widehat {u_b}^2$
. The profiles are shown at seven streamwise locations:
$x/h=0.1$
, 2, 4.5, 6.5, 8, 10 and 12. 
, Coarse grid; 
, medium grid; 
, fine grid; 
, reference level for each section.
Figure 4 shows the vertical profiles of the time- and span-averaged streamwise velocity,
${\langle \overline {u} \rangle _z} /\widehat {u_b}$
(figure 4
a), shear Reynolds stress,
${\langle \overline {{{u}^{\prime }}{{v}^{\prime }} } \rangle _z} /\widehat {u_b}^2$
(figure 4
b) and the turbulent kinetic energy,
$2{\langle \overline {{\mathcal{K}} } \rangle _z} /\widehat {u_b}^2= {\langle \overline {u'_iu'_i} \rangle _z} /\widehat {u_b}^2$
(figure 4
c). Only the fluid domain is shown in this and all the following figures.
The streamwise velocity in figure 4(a) shows very good convergence between the three grids, with negligible difference over most of the domain. The second-order statistics in figures 4(b) and 4(c) also show good convergence towards the solution on the fine grid over most of the fluid domain. The largest difference is observed in the separated region off of the larger crest. In the remainder of the paper, the medium grid resolution will be used for the simulations of dunes with a mobile bed.
Figure 1(a) shows the grid (medium) used for the mobile-bed calculations; the parameters of this grid are summarised in table 2. The grid in the
$xy$
plane is extruded in
$z$
, the direction in which the spacing is constant. As mentioned, the bottom boundary of the computational grid is positioned deep enough within the sediment volume to prevent mobilisation of the entire sediment mass.
3. Results
In this section, we present the results of the simulation of a laboratory-scale, realistic river geometry. The numerical methodology illustrated in the previous sections is employed to perform high-Reynolds-number simulations of the flow over a bed composed of loose sediment. With the aid of the experiments by Reesink et al. (Reference Reesink, Parsons, Ashworth, Best, Hardy, Murphy, McLelland and Unsworth2018), we study a configuration where the bed presents both large morphological structures and superimposed bedforms.
3.1. Bed geometry
Bed surface evolution over time: (a)
$t/T_s\times 10^{-6}=0$
, (b)
$t/T_s\times 10^{-6}=0.3$
, (c)
$t/T_s\times 10^{-6}=0.55$
, (d)
$t/T_s\times 10^{-6}=0.8$
.
Figure 5 shows the bed surface at four different time instants. The initial bed surface (figure 5
a) is uniform in the spanwise direction. Figures 5(b)–5(d) show the bed surface at times
$t/T_s\times 10^{-6}=0.3,\ 0.55,\ 0.8$
, respectively,
$T_s=d/w_s$
being the characteristic sediment time scale. The formation of superimposed bedforms takes place over the entire bed surface: the morphology of the stoss side of the first dune (
$8\lt x/h\lt 13$
) is characterised by trains of transverse superimposed bedforms, while that of the second dune (
$1\lt x/h\lt 5$
) shows local crests and troughs that are strongly three-dimensional and tend to align with the flow. The bedforms at the beginning of the first slope (
$6\lt x/h\lt 8$
) show similar characteristics. On the other hand, the bedforms that form over the stoss side of the second crest are transverse to the flow and extend over the entire width of the channel. They display a dune-like triangular section with a lee side steeper than the stoss side. A comparison of the location of the crest lines in figures 5(c) and 5(d) highlights the increasing three-dimensionality of the bedforms. The bedform on the stoss side of the first dune shows closed-loop crests that travel up the stoss side and develop into transverse crests, as also observed by Venditti et al. (Reference Venditti, Church and Bennett2005a
). In the following, we refer to the large-scale bedforms, namely the two large crests and troughs present in the initial bed geometry, as the ‘main’ bedforms, to distinguish them from the superimposed small-scale bedforms.
Figure 6 shows the spanwise-averaged bed elevation profiles at the same four times shown in figure 5, and highlights that the main features of the bedform do not vary considerably: we notice a slightly sharper crest angle of the first main crest (
$x/h\simeq 0$
), a slightly eroded second main crest (
$x/h\simeq 5$
) and mild deposition in the trough downstream of the first main crest (
$1\lt x/h\lt 5$
). As the superimposed bedforms migrate downstream, they reach the main crest and form a rugged, irregular brink (
$0\lt x/h\lt 1$
). Subsequently, the sediment slides down and deposits into the trough. By comparing the lee side of the first main dune at the four time instants, we notice that it is migrating forward, and so does the second main dune.
To determine the celerity of the crests, we calculate the auto-correlation coefficient:
\begin{equation} R (\xi, \tau ) = \frac { \left \langle \langle {y_b}\rangle _z(x,t)\langle {y_b}\rangle _z(x+\xi, t+\tau )\right \rangle _z} {\left \langle \langle {y_b}\rangle _z(x,t)^2\right \rangle _z} \; . \end{equation}
At each instant, we use an iterative algorithm to determine the spatial shift,
$\xi ^*(t)$
, that maximises the autocorrelation:
$R(\xi ^*,t)=R(\xi, t)\big |_{\max }$
. Shifting the bedform by
$\xi ^*(t)$
gives the best possible overlap of the geometries as they migrate. The celerity can finally be calculated as the derivative of the forward shift of the dune crests over time,
$c_s={d}\xi ^*(t)/\textrm{d}t$
. Figure 6(b) shows the shifted spanwise-averaged profiles
$\langle {\widetilde {y_b}}\rangle _z=\langle {y_b}\rangle _z(x+\xi ^*(t),t)$
.
Flow over dunes with mobile bed. (a) Spanwise-averaged profiles of the channel bed,
$\langle {y_b}\rangle _z$
at 
,
$t/T_s=0$
; 
,
$t/T_s=0.3\times 10^6$
; 
,
$t/T_s=0.55\times 10^6$
; and 
,
$t/T_s=0.8\times 10^6$
. (b) Spanwise-averaged and shifted profiles of the channel bed; 
,
$\widetilde {y}_b^{ref}(x)$
is the averaged profile during the quasi-steady-state period.
The evolution of the channel bed undergoes an initial transient phase during which the superimposed bedforms start to form and grow in size and a second phase in which the bed morphology is statistically invariant. The celerity of the main crests settles, after a relatively long transient (
$0.3\times 10^6$
$T_s$
), on a value of approximately
$c/w_s\approx 0.003$
, which does not vary considerably until the end of the simulated time window.
Contours of the spanwise-averaged profiles of the channel bed,
$\langle {y_b}\rangle _z$
, over time. The dashed line identifies the location of a dune crest over time, and its slope indicates its celerity.
Time evolution of the total drag coefficient,
$C_D$
. 
, Instantaneous drag coefficient; 
, average value of the drag coefficient in the quasi-steady-state interval,
$T_{qss}/T_s\times 10^{-6}=[0.3{-}0.8]$
.
Figure 7 shows contours of the spanwise-averaged profiles of the channel bed,
$\langle {y_b}\rangle _z$
, over time. The first main crest, the second main crest and the superimposed bedforms have different celerities, as shown in the figure. The celerity shows an inverse relationship with dune size, as expected. It is interesting to notice that the velocity of the superimposed dunes reduces during the simulation as they grow in size, to reach a quasi-steady state.
Figure 8 shows the time evolution of the total drag coefficient of the surface, which is the sum of the friction drag coefficient,
$C_{D,f}$
, and the pressure drag coefficient,
$C_{D,p}$
:
\begin{equation} C_D=C_{D,f}+C_{D,p} =\frac {1}{S_b}\int _{S_b}C_f\hat {\boldsymbol \theta }_b \cdot \hat {\mathbf{x}}\textrm{d}A +\frac {1}{S_b}\int _{S_b}-C_p{\hat {\mathbf n}_{\mathbf{b}}} \cdot \hat {\mathbf{x}}\textrm{d}A \; . \end{equation}
The instantaneous force coefficients,
$C_f$
and
$C_p$
, are defined as
\begin{align} C_f=\dfrac {\tau _w}{\rho \widehat {u_b}^2/2} \;, \qquad C_p=\dfrac {p_w}{\rho \widehat {u_b}^2/2} \;, \end{align}
where
$p_w$
is the pressure at the immersed-boundary location. The time evolution of the drag coefficient highlights the presence of the initial transient phase and of the steady-state phase for
$t/T_s=T_{qss}/T_s=[0.3{-}0.8]\times 10^6$
.
During this time interval, the large features of the dunes translate almost as a solid body, while smaller bedforms grow and decay on a shorter time scale. On the basis of these observations, we can assume that a dynamic equilibrium between the flow and the main bedform has been achieved. Therefore, we consider the system to be in a quasi-steady state in the time interval
$T_{qss}$
, during which we take an ensemble average of the shifted bedforms to define the reference profile:
\begin{align} \widetilde {y}_b^{ref}(x) &=\frac {1}{T_{qss}}\int _{T_{qss}} \langle {y_b}\rangle _z(x+\xi ^*(t),t)\textrm{d}t \;, \end{align}
\begin{align} \widetilde {x}_b^{ref}(x) &=\frac {1}{T_{qss}}\int _{T_{qss}} \langle {x_b}\rangle _z(x+\xi ^*(t),t)\textrm{d}t \;, \end{align}
where the integral is extended over the quasi-steady-state period,
$T_{qss}$
, only. The reference profile given by (3.4) is also shown in figure 6.
Knowing
$\xi ^*$
also allows us to perform time averages of the flow properties, by shifting each instantaneous field by the appropriate distance, so that all the bedforms have maximum overlap. The shift previously computed is then applied backwards to the entire field and the reference surface, before calculating any instantaneous quantity. Note that we use the reference profile calculated via (3.4) for the quasi-steady period, while for the transient part of the bed evolution we use the initial profile as a reference. These geometrical definitions will aid and streamline our discussion of the results and their interpretation in the context of a changing morphology.
With these definitions in place, we can now explore how the migration of the main bedforms is influenced by the superimposed bedforms, which transport packets of sediment across the dune and down into the trough. The sediment is then trapped in regions where transport is lower (
$\theta$
is lower) causing the dune to move. Both main crests migrate downstream but at different migration rates (figure 7). This indicates that the superimposed bedforms still travel downstream and deposit sediment inside the trough.
(a–d) Instantaneous contours of the elevation variation normalised by the sand grain size,
$(y_b-\widetilde {y}_b^{ref})/d$
, at
$t/T_s\times 10^{-6}=0$
, 0.3, 0.55, 0.8. 
,
$(y_b-\widetilde {y}_b^{ref})/d=20$
; 
, first main crest; 
, first main trough; 
, second main crest; 
, second main trough.
To elucidate this point we show top-view contours of
$(y_b-\widetilde {y}_b^{ref})/d$
in figure 9 at the same times as figure 5. The locations of the main crests and troughs (identified as the local maxima and minima of the spanwise-averaged profiles) are indicated by the red and blue lines, and the black contours identify the crests of the superimposed bedforms. The first main crest moves forward to
$x/h\approx 0.5$
as soon as the first superimposed bedform reaches the main crest and later travels forward more slowly. The downstream motion of the transverse superimposed bedforms over the stoss side of the first main dune is readily visible. On the other hand, it is less clear how the second main dune is migrating forward.
(a–d) Instantaneous contours of the elevation variation normalised by the sand grain size,
$(y_b-\widetilde {y}_b^{ref})/d$
, at
$t/T_s\times 10^{-6}=0.12$
, 0.17, 0.21, 0.25. 
,
$(y_b-\widetilde {y}_b^{ref})/d=20$
; 
, second main crest; 
, fixed rectangular box for reference.
To explain this mechanism, in figure 10 we show a zoom of the region
$4\lt x/h\lt 6$
. The magenta boxes in the figure highlight highly irregular and short-lived crest lines near the second main crest (figure 10
a) and travelling downstream for a short distance (barely visible in figure 10
b,c) before cascading down into the trough (figure 10
d). The fairly regular streamwise contour lines in the regions
$5\lt x/h\lt 6$
represent the sediment cascading down the lee side of the two main crests.
To further characterise the time development of the superimposed bedforms, we identify the crests and troughs by locating the local maxima and minima of the superimposed bedforms with respect to the reference profile (3.4). Figure 11 shows the evolution in time of height and wavelength of the superimposed bedforms. After the transient, the superimposed bedform height settles on a statistically steady value of
$\varDelta /d\simeq 31.5$
and the wavelength on
$\lambda /d\simeq 981$
. These measurements are in reasonable agreement with the profiles measured by Reesink et al. (Reference Reesink, Parsons, Ashworth, Best, Hardy, Murphy, McLelland and Unsworth2018), where
$\varDelta /d\approx 45$
and
$\lambda /d=500{-}1200$
.
Spanwise-averaged bedform wavelength averaged over all bedforms,
$\lambda /d$
, on the left-hand axis (
); span-averaged bedform height,
$\varDelta /d$
, on the right-hand axis (
).
3.2. Local skin-friction and pressure coefficient
The wall pressure (not shown) decreases over the stoss side of each superimposed bedform, reaches a local minimum near the crest line and then remains nearly constant over the local trough (as expected in separated flows).
(a–d) Instantaneous contours of the friction coefficient,
$C_f(x,z,t)$
, at
$t/T_s\times 10^{-6}=0$
, 0.3, 0.55, 0.8. 
,
$(y_b-\widetilde {y}_b^{ref})/d=20$
; 
, first main crest; 
, first main trough; 
, second main crest; 
, second main trough.
Figure 12 shows contours of the instantaneous friction coefficient,
$C_f$
. Notice that the friction due to sand-grain roughness is included in the friction coefficient through the wall model, by including a roughness displacement (
$y_0\simeq 0.17d$
) into the wall-model equations. On the other hand, the form drag due to the separation behind the superimposed bedforms is explicitly obtained from the immersed-boundary force. Coefficient
$C_f$
increases over the stoss side of the superimposed bedforms and then suddenly drops.
Figure 12(a) shows the friction coefficient at the initial time of the simulation, when the bed morphology is uniform in the spanwise direction, but the velocity field, although also statistically homogeneous, is three-dimensional. The streaky structure of
$C_f$
over the stoss side is similar to that in boundary layers with a zero pressure gradient or mild favourable pressure gradient. The recirculation bubbles behind the two main crests are highly three-dimensional and time dependent.
As time progresses, the skin-friction coefficient indicates the presence of recirculation bubbles downstream of the crest lines of the superimposed bedforms. Also, on average,
$C_f$
on the two stoss sides (
$4\lt x/h\lt 5$
and
$8\lt x/h\lt 13$
) is smaller in the presence of the superimposed bedforms. This is particularly evident on the stoss side of the first main crest (
$8\lt x/h\lt 13$
). The first stoss side alternates streaks of low and high
$C_f$
, corresponding respectively to the troughs and crests of the superimposed bedforms.
The total drag coefficient over the mobile bed, however, is larger than that over a fixed bed. Figure 13 shows the time evolution of the drag for the fixed- and mobile-bed configurations. Both friction and pressure components of drag are shown. When the bed is fixed, most of the drag is due to pressure, and the friction component is only 11 % of the total drag. In the mobile-bed simulation the fraction of drag due to friction is even lower (
$\simeq 3\,\%$
) due to the flow reversal, which produces a negative drag (thrust) on the wall. Furthermore, the pressure-drag component is roughly two times larger (at a quasi-steady state) than in the fixed-bed case.
Time evolution of the friction- and pressure-drag coefficients. Fixed bed: 
,
$\,C^{fix}_{D,p}$
; 
,
$C^{fix}_{D,f}$
. Mobile bed: 
,
$C^{mob}_{D,p}$
; 
,
$C^{mob}_{D,f}$
; 
,
$\Delta C_{D,p}^{SB}$
; 
,
$\Delta C_{D,p}^{main}$
.
To determine the cause of this increase of form drag, we define two additional pressure coefficients that separate the added form drag due to the larger recirculation bubble on the main geometry,
$\Delta C_{D,p}^{main}$
, from the form drag due to the recirculation behind the superimposed bedforms,
$\Delta C_{D,p}^{SB}$
, such that
$C_{D,p}^{mob}(t)= {\overline {C_{D,p}^{fix}}} +\Delta C_{D,p}(t)^{main}+\Delta C_{D,p}^{SB}(t)$
, with
\begin{align} &\Delta C_{D,p}^{main}(t) =\frac {1}{S_b}\int _0^W\int _{x_C}^{x_R}\left ( -C_p^{mob}+ {\overline {C_p^{fix}}} \right ){\hat {\mathbf n}_{\mathbf{b}}\cdot \hat {x}}\,\textrm{d}x\,\textrm{d}z \;, \end{align}
\begin{align} &\Delta C_{D,p}^{SB}(t) =\frac {1}{S_b}\int _0^W\int _{x_R}^{L+x_C}\left ( -C_p^{mob}+ {\overline {C_p^{fix}}} \right ){\hat {\mathbf n}_{\mathbf{b}}\cdot \hat {x}}\,\textrm{d}x\,\textrm{d}z \;, \end{align}
where the superscripts ‘mob’ and ‘fix’ refer to the mobile-bed and fixed-bed cases, respectively, while
$x_R$
is the reattachment location and
$x_C$
is the location of the first main crest in the mobile-bed case.
The time evolution of these quantities is also shown in figure 13. In the initial phase of the bed evolution, the additional form drag is entirely due to the superimposed bedforms, since
$\Delta C_{D,p}^{main}$
oscillates around zero. In the subsequent phases of the bedform evolution, the superimposed bedforms generate larger local recirculation regions that contribute significantly to the total drag. The time- and spanwise-averaged pressure coefficient due to the superimposed bedforms in the quasi-steady-state interval is
${\overline {\Delta C_{D,p}^{SB}}} =0.032$
, contributing 42 % of the total form drag. It is also interesting to notice that the form drag due to the main recirculation bubbles increases roughly at the time the quasi-steady state is reached (
$t/T_s\times 10^{-6}\gt 0.3$
), showing a correlation with the presence of superimposed bedforms. This suggests that the presence of recirculation bubbles behind the superimposed bedforms promotes larger separation on the main geometry as well, a mechanism that we will further elucidate momentarily. The contribution of
${\overline {\Delta C_{D,p}^{main}}} =0.019$
in the quasi-steady-state interval is 25 % of the total form drag.
3.3. Dynamics of the main recirculation bubbles
To support the previous findings, we analysed the variation in time of the two recirculation bubbles (downstream of the two main crests). The lengths of the two recirculation bubbles,
$\ell _1$
and
$\ell _2$
, are identified by the locus of points where the spanwise-averaged skin-friction coefficient,
$\langle {C_f}\rangle _z(x,t)$
, is less than or equal to zero. The average and root-mean-square of the length variation in time are summarised in table 3. As anticipated in § 3.2, the recirculation bubbles on the mobile bed are larger, by 20 % and 10 %, respectively.
Recirculation-bubble characteristics. Average and relative root mean square of the two main recirculation bubbles (labelled as 1 and 2, respectively) for the simulations with fixed and mobile beds.
Contours of the instantaneous streamwise velocity in the middle plane near the two main crests:
$u/\widehat {u}_b$
at (a)
$t/T_s\times 10^{-6}=0$
, (b)
$t/T_s\times 10^{-6}=0.3$
, (c)
$t/T_s\times 10^{-6}=0.55$
and (d)
$t/T_s\times 10^{-6}=0.8$
. The black dashed lines 
denote the maximum and minimum bed level during the calculation.
All separated flows exhibit some level of unsteadiness. The reattachment point, in particular, always oscillates. Unless separation occurs at a sharp corner, the separation point may oscillate as well. In this case, an additional cause of unsteadiness, although with a longer time scale, is present, namely the fact that the geometry of the crest is itself unsteady. Figure 14 shows instantaneous contours of the streamwise velocity in the middle plane for the fixed bed and for the mobile bed at three widely spaced times. The panels focus on the crests of the two main dunes. By comparing figures 14(b) and 14(c), we can see that the separation point of the first main dune shifts upstream. This shift occurs because the main recirculation region merges with the local recirculation behind the superimposed bedform, which reaches the crest. The interconnection between the train of superimposed bedforms reaching the main crest and the variation in the size of the recirculation bubbles is linked to the increase in form drag generated by the main geometry,
$\Delta C_{D,p}^{main}$
. Figure 14 also highlights the motion of the superimposed bedforms, shown in more detail in supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.334. Interestingly, superimposition seems to occur at multiple scales. In fact, smaller bedforms can be seen travelling on top of the superimposed bedforms.
The streamwise velocity over the crest is larger over the mobile bed than over the fixed bed due to the flow constriction caused by the superimposed bedforms, which can significantly decrease the cross-sectional area (compare for instance, figures 14 a and 14 d). Another reason for the increased average velocity is the creation of a low-momentum zone near the wall, due to the flow separation behind the superimposed bedforms on the stoss side. Similar phenomena can be observed on both main crests, although on the second crest, the flow is also affected by the recirculation bubble downstream of the first dune. As a consequence of the increased mean velocity, the flow outside the low-momentum region has more inertia, and the angle formed by the separated shear layer with respect to the horizontal is larger. This can be also observed from the mean flow: figures 15 and 16 show contours of the time- and spanwise-averaged streamwise velocity and spanwise vorticity:
\begin{equation} {\left \langle \overline {\omega _z} \right \rangle _z} =\frac {\partial {\left \langle \overline {u} \right \rangle _z} }{\partial y}-\frac {\partial {\left \langle \overline {v} \right \rangle _z} }{\partial x} \; . \end{equation}
As mentioned in § 3.1, the time averaging is performed after shifting the geometry to account for the motion of the large-scale dune. Figure 16 shows an increase in spanwise vorticity near the bed over the stoss side of the first main dune, determined by the flow separation over the crests of the superimposed bedforms, confirming the results of Hardy et al. (Reference Hardy, Best, Marjoribanks, Parsons and Ashworth2021). As observed by Zgheib et al. (Reference Zgheib, Fedele, Hoyal, Perillo and Balachandar2018), regions of high vorticity remain confined within the wall region, extending only to a distance approximately equal to the height of the superimposed bedforms. Figures 15 and 16 show the mean flow modifications due to the bed mobility: the different inclination of the separated shear layer, the larger region of reversed flow and the formation of the low-momentum region on the stoss side of the two main dunes.
Contours of time- and spanwise-averaged streamwise velocity
${\langle \overline {u} \rangle _z} /\widehat {u_b}$
. (a) Fixed bed. (b) Mobile bed. The dashed lines 
denote the maximum and minimum bed level during the calculation.
Contours of the time- and spanwise-averaged spanwise vorticity,
${\langle \overline {\omega _z} \rangle _z} {\langle {H}\rangle } /\widehat {u_b}$
. (a) Fixed bed. (b) Mobile bed. The dashed lines 
denote the maximum and minimum bed level measured during the simulation.
Contours of the instantaneous turbulent kinetic energy in the middle plane,
$2 {\mathcal{K}} /\widehat {u_b}^2$
, at (a)
$t/T_s\times 10^{-6}=0$
, (b)
$t/T_s\times 10^{-6}=0.3$
, (c)
$t/T_s\times 10^{-6}=0.55$
and (d)
$t/T_s\times 10^{-6}=0.8$
. The black dashed lines 
denote the maximum and minimum bed level measured during the simulation.
The instantaneous contours of the normalised turbulent kinetic energy,
$2{\mathcal{K}} /\widehat {u_b}^2=u_iu_i/\widehat {u_b}^2$
, shown in figure 17, highlight the role that the shear layers play in the generation of turbulent fluctuations. In the separated shear layer that departs from the crest the turbulent production is large, and the Reynolds stresses are amplified due to the shear-layer instability. In addition to the main shear layers, however, smaller ones are separating from the crests of the superimposed dunes, contributing to an increase in the turbulence levels near the wall when the bed is mobile, which is reflected in the mean values (figure 18).
Contours of the time- and spanwise-averaged turbulent kinetic energy
$2{\langle \overline { {\mathcal{K}} } \rangle _z} /\widehat {u_b}^2$
. (a) Fixed bed. (b) Mobile bed. The black dashed lines 
denote the maximum and minimum bed level measured during the simulation.
3.4. Roughness effects
As previously pointed out, the superimposed bedforms create separated shear layers that amplify the turbulent fluctuations, cause form drag and generate a region of low momentum near the wall. These are typical phenomena observed in flows over rough surfaces, and the present flow is characterised by two scales of roughness: the sand-grain size and the superimposed bedforms.
Acceleration parameter,
$K$
. (a) Bed time- and spanwise-averaged geometry: 
, fixed-bed; 
, mobile-bed reference profiles; 
, locations of sections I, II, III and IV. (b) Acceleration parameter along the fixed-bed dune. (c) Acceleration parameter along the mobile-bed dune. The zones of 
APG and 
FPG are highlighted.
The flow is characterised by a succession of favourable pressure gradient (FPG) and adverse pressure gradient (APG) regions. Figure 19 shows the acceleration parameter:
\begin{equation} K=\frac {\nu }{u_\infty ^2}\frac {\textrm{d}u_\infty }{\textrm{d}x} \;, \end{equation}
where
$u_\infty = {\langle \overline {u} \rangle _z} (x,y=L_y)$
is the time- and spanwise-averaged streamwise velocity at the top boundary. A positive
$K$
indicates an FPG and a negative
$K$
indicates an APG. Although
$K$
is not large enough to cause relaminarisation of the flow (Spalart Reference Spalart1986), significant APG and FPG effects can be noticed.
Profiles of the wall-parallel velocity at (a–d) the locations shown in figure 19. 
, Fixed bed; 
, mobile bed; 
, logarithmic law of the wall, with slope
$0.12H/d_n$
.
Figure 20 shows the mean-velocity profiles at four streamwise locations as a function of the wall-normal distance,
$d_n$
. We observe a slight decrease of the velocity gradient in the outer layer in FPG sections (Piomelli & Yuan Reference Piomelli and Yuan2013), and an increase of the gradient in APG regions (Knopp et al. Reference Knopp, Reuther, Novara, Schanz, Schülein, Schröder and Kähler2021). As a consequence of the superimposed bedforms, which act as large roughness elements, the momentum deficit in the near-wall region is significantly larger when the bed is mobile. The slope
$0.12H/d_n$
, shown in figure 20, is the same that characterises the logarithmic region in a channel flow with zero pressure gradient (Bradshaw & Huang Reference Bradshaw and Huang1995; Pope Reference Pope2000), and serves as a reference.
In the APG region (sections I and IV in figure 20) the velocity profile shows mild features of accelerating flows, such as a decreased slope of the velocity profile. Conversely, in sections II and III, located in the APG region, the slope of the velocity profile away from the wall increases, as is typical of APG flows. The slope change is more significant when the bed is mobile (probably due to the near-wall momentum deficit).
Note that, as mentioned above, the sand-grain-size roughness is included through the generalised Moody diagram (Meneveau Reference Meneveau2020). Therefore, the differences in velocity distribution (e.g. near-wall momentum deficit, wall-normal velocity gradient, acceleration parameter etc.) between the fixed and mobile cases are entirely attributable to the three-dimensional morphology.
3.5. Sediment transport
The non-dimensional bed-load transport rate is calculated as
\begin{equation} \varPhi =|{{{\mathbf q}_{\mathbf{bl}}}} |/\sqrt {(s-1)gd^3} \; . \end{equation}
Figure 21 shows that the bed-load transport (and its gradient) are largest on the crests of the superimposed bedforms, highlighting that they control most of the sediment transport. In the separated flow region (
$2\lt x/h\lt 6$
), we observe small areas where the wall shear stress is sufficiently high to erode the bed and transport minor amounts of sediment upstream, as illustrated in figure 22, which focuses on the recirculation zones.
(a–d) Contours of the instantaneous non-dimensional bed-load sediment transport rate,
$\varPhi$
, at
$t/T_s\times 10^{-6}=0$
, 0.3, 0.55, 0.8. 
,
$(y_b-\widetilde {y}_b^{ref})/d=20$
; 
, first main crest; 
, first main trough; 
, second main crest; 
, second main trough.
(a) Contours of the instantaneous magnitude of the non-dimensional bed-load sediment transport vector,
$\varPhi$
, at
$t/T_s\times 10^{-6}=0.2$
. (b) Contours of the instantaneous streamwise component of the bed-load-layer velocity,
${{{\mathbf U}^\vartheta _{\mathbf{bl}}}} /\widehat {u_b}$
, at
$t/T_s\times 10^{-6}=0.8$
. 
,
$(y_b-\widetilde {y}_b^{ref})/d=20$
; 
, first main trough; 
, second main crest; 
, upstream transport; 
, downstream transport.
The non-dimensional bed-load transport rate is shown in figure 22(a), and the streamwise component of the velocity at the edge of the bed-load layer in figure 22(b), at
$t/T_s\times 10^{-6}=0.2$
. The green box highlights regions where sediment is eroded by the reverse flow and transported upstream, as indicated by the velocity contours, towards the lee side of the first main dune. This process, in addition to landslides from the crests, contributes to sediment deposition in the trough of the first main dune. The recirculation region is highly three-dimensional; while most sediment in the trough is transported upstream by the recirculating flow, on the stoss side of the smaller main crest (
$4\lt x/h\lt 5$
), sediment can be transported either upstream or downstream.
Similarly to the bed-load transport rate, we compute the non-dimensional suspended-sediment transport rate:
\begin{equation} \varPhi _s=\frac {|{{\mathbf q}_{\boldsymbol{s}}}|} {\sqrt {(s-1)gd}}=0.00033d^{0.3}T^{1.5} \;, \end{equation}
where
$q_s$
is the suspended-sediment erosion rate (van Rijn Reference van Rijn1984a
). The erosion rate of suspended sediment is two orders of magnitude smaller than the bed-load transport, thus validating our choice to neglect it in the present simulation.
4. Conclusions
We have studied the flow over a realistic dune geometry at a laboratory scale, to investigate the response of a loose sediment bed to a turbulent flow at a Reynolds number of
$121\,900$
, intending to analyse a flow configuration that resulted in the formation of superimposed bedforms. Our results show that the morphology of the superimposed bedforms varies, depending on the local state of the flow: the bedforms exhibit a very three-dimensional shape in the recirculation regions, and tend to be oriented in the streamwise direction; on the other hand, they tend to be a transverse over the stoss side of the larger crest and they may become large enough to cause local detachment of the flow and recirculation. As the bedforms migrate downstream over time, they become more three-dimensional. The shape of the main crest line and the configuration of the brink also change as the superimposed bedforms reach the top of the stoss side from upstream.
The bedform reaches a quasi-steady state (statistically) after a long transient: the size and speed of the superimposed bedforms cease to change significantly, but fluctuate around an average height of
$\varDelta /d=31.5$
and length
$\lambda /d=981$
, and a speed of
$c/w_s\approx 0.003$
. Although the mean profile of the dune does not vary considerably during this quasi-steady period, we observe slight changes in crest angle, erosion of the second main crest and mild deposition in the trough downstream of the first main crest. Nonetheless, these changes are sufficient to induce strong variations of the flow features, indicating that the presence of superimposed bedforms plays an important role in determining the global dynamics of the system.
Our results also show that the presence of superimposed bedforms leads to an increase in the size and temporal oscillation of the recirculation bubbles. The recirculation behind the superimposed bedforms causes pressure drops that contribute significantly to the form drag. The shear layer near the wall is weakened in the presence of superimposed bedforms, due to the near-wall momentum deficit over the stoss side of the preceding main dune and to the irregular shape of the brink. We also find that the shear layer deviates towards the top surface when the recirculation bubble is very large due to the sharper crest and irregular brink. Overall, the main effect of the superimposed bedforms on the stoss side of the first main dune is to reduce the momentum near the wall and decrease the velocity gradient at the wall, where form drag supplies most of the resistance to the flow.
Using the present model, which combines wall-modelled LES, an IBM-based bedform description and a transport equation for the sediment motion, we were able to simulate successfully a realistic hydromorphodynamic system and investigate the main effects that superimposed bedforms have on the global flow dynamics. Although the investigation focuses on a single specific morphology, we suggest that several observations regarding the influence of superimposed bedforms on flow dynamics may be applicable across a range of Reynolds numbers, provided that similar bedforms develop and generate local flow separation. The findings discussed in this paper contribute to our understanding of how such features affect large-scale flow and sediment dynamics in natural and engineered channels, potentially supporting a more precise quantification of their impact on overall flow dynamics. The knowledge and quantification of such impact are crucial to improve the accuracy of riverine system simulations on a large scale. Our study addresses the entire spectrum of spatial scales that bedforms may exhibit and discusses the inaccuracies that may be expected if these scales are underestimated or neglected, thus paving the way for improved predictions of sediment transport and bedform morphodynamics.
In addition, the study contributes to the assessment of the scaling behaviour that these phenomena exhibit with the flow and sediment parameters. Reesink et al. (Reference Reesink, Parsons, Ashworth, Best, Hardy, Murphy, McLelland and Unsworth2018) revealed the relation between bedforms and flow parameters such as flow depth and flow regime. Reesink et al. (Reference Reesink, Parsons, Ashworth, Best, Hardy, Murphy, McLelland and Unsworth2018) noted that superimposed bedforms appear in response to an increase in water depth (
$Fr$
increase); on the other hand, increasing the flow regime (
$Re$
increase) flattens the bedform field. This aligns with our observations, as the superimposed bedforms form and migrate through the mechanism of cascade and deposition in the troughs, which would be inhibited by a stronger shear layer over the superimposed bedform crests. Estimating the limits of this transition would require further simulations.
The results discussed in this paper also highlight an intriguing aspect of superimposition, specifically its occurrence at multiple scales. The bedform elevation contours (see figure 9) reveal at least three scales of bedforms: the main dunes, the superimposed bedforms extensively described in this study and very small-scale bedforms that travel over the superimposed bedforms at a higher migration rate. This observation raises the question of how small these structures can become and what the smallest scale is at which sediment transport can occur through this mechanism. The scale of these structures approaches the grid resolution, which means the current simulation could not capture the dynamics of anything smaller. Investigating these small-scale structures further would be desirable but would require more advanced simulations to fully understand their behaviour and impact.
The numerical toolbox assembled in our study appears suitable for a wide range of wall-bounded flows, Reynolds numbers and flow configurations. A natural expansion of the present work would be to simulate other flow configurations studied by Reesink et al. (Reference Reesink, Parsons, Ashworth, Best, Hardy, Murphy, McLelland and Unsworth2018) to understand further the response of superimposed bedforms to the flow and vice versa. For instance, it would be interesting to study the phenomenon of flattening of the dunes when the flow depth is reduced, as understanding the cause of the disappearance of superimposed bedforms could bring to light different features of the flow that foster bedform formation.
Another possible route to extend the study would be to consider realistic Reynolds numbers, rather than the laboratory-scale ones used here. If we were to simulate a natural river flow, the Reynolds number would be considerably larger; however, the computational cost of the simulation would still be acceptable, since the grid size required for wall-modelled LES does not scale in viscous units. For instance, a realistic river configuration with an average water depth and width of
$\mathcal{O}(10\,\textrm {m})$
and
$\mathcal{O}(100\,\textrm {m})$
, respectively, with an average annual flow discharge of
$\mathcal{O}(1 \times 10^3 \textrm {m}^3\,\textrm {s}^{-1})$
, results in a bulk Reynolds number of
$\mathcal{O}(10^7)$
, two orders of magnitude larger than in the present laboratory-scale simulation. The time and length scales influence the numerical requirements of the simulation, in terms of both grid resolution and time advancement iterations. Since only the large scales of the flow need to be captured in wall-modelled LES, whose size is a weak function of the Reynolds number, the number of points required would increase approximately by a factor of 16. The current simulation required approximately 150 000 CPU hours to cover a physical time of 4 hours, while the realistic river flow would require 2.5 million CPU hours, a number that is large, yet manageable.
Supplementary movie
The supplementary movie is available at https://doi.org/10.1017/jfm.2025.334.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Validation and assessment of the model
In this appendix, we validate the numerical methods implemented, namely the IBM for the treatment of the complex and movable wall boundary, the wall model for the treatment of the wall boundary condition and the morphological model for the simulation of the evolving channel bed. To test the wall model and the IBM, we first compare our results with reference data for the flow in an open channel at a high Reynolds number and over periodic hills. For the morphological model we reproduce numerically a laboratory experiment (Venditti et al. Reference Venditti, Church and Bennett2005a ) of an open-channel flow over a loose, initially flat, granular bed.
A.1. Channel flow
First, we validate the implementation of the wall model coupled with the IBM by performing simulations of a plane channel. The domain is a cuboid of length
$L=4H$
, width
$W=2H$
and height
$1.1H$
, and it is discretised with a Cartesian grid. The immersed surface is a flat horizontal surface, located at
$y_b=0.1H$
, such that the effective flow depth is
$H$
. The flow is forced with a constant mean-pressure gradient,
$\Pi$
, resulting in a Reynolds number
$Re_\tau \simeq 5200$
, for which high-quality DNS data are available for comparison (Lee & Moser Reference Lee and Moser2015). We use periodic boundary conditions in the streamwise and spanwise directions, the top of the domain is a free-slip surface and at
$y=0.1H$
(the solid surface modelled by the IBM) the wall model discussed earlier is used. The computational domain is discretised with a uniform grid in the streamwise and spanwise directions with 128 points, respectively. The spacing in wall units of the grid in these two directions is
$\Delta x^+=156$
and
$\Delta z^+=78$
, respectively. In the vertical direction, the grid has a region of uniform fine resolution near the immersed boundary and is then stretched upwards through a hyperbolic tangent distribution, with a total of 128 points. The resulting grid spacing near the immersed surface is therefore
$\Delta y^+_{min}=20$
and
$\Delta y^+_{max}=56$
near the top surface.
The main issue we faced was the presence of numerical oscillations near the wall and the free surface, as shown in figure 23. At the free surface, the grid is coarse and neither the molecular nor the sub-filter-scale dissipation is sufficient to dampen the saw-tooth oscillations introduced by the solution of (2.1b ) on a staggered grid (Ferziger & Peric Reference Ferziger and Peric2002). Below the inner-/outer-layer interface, the oscillations were due to the application of the IBM forcing. The refinement of the grid reduces their intensity but does not remove them completely. To remove them, we blended the central-difference scheme with a second-order-accurate upwind method, the QUICK method:
\begin{equation} \frac {{\textrm{d}}\bullet }{\textrm{d}x} =f\left (\frac {d\bullet }{\textrm{d}x}\right )_{QUICK} +(1-f)\left (\frac {{\textrm{d}}\bullet }{\textrm{d}x}\right )_{CD}. \end{equation}
Figure 23 shows that the blending function
$f$
restricts the upwind method to be active only below the interface and near the top boundary, and is effective in removing the unphysical oscillations.
Figure 24 shows the wall-normal profiles of the time- and spanwise-averaged streamwise velocity,
${\langle \overline {u} \rangle }$
, shear Reynolds stress,
${\langle \overline {{{u}^{\prime }} {{v}^{\prime }} } \rangle }$
, and normal Reynolds stresses,
${\langle \overline {{{u}^{\prime }} {{u}^{\prime }} } \rangle }$
and
$ {\langle \overline { {{v}^{\prime }} {{v}^{\prime }} } \rangle }$
, in wall units. The vertical black dashed line (
) in figure 24 indicates the location of the wall-model interface, i.e. the location at which the
${{\mathbf U}^\vartheta _{LES}}$
is extracted from the LES solution and fed to the wall model. The use of a fine grid coupled with the blended method gives reasonably accurate results. In wall-modelled LES, the coarseness of the grid causes low Reynolds stresses because only part of the fluctuation spectrum is resolved. Since we simulate a channel flow, the fluctuations of the vertical velocity component,
$v'$
, are driven to zero at the surface while the streamwise and spanwise ones increase. Note that the region below the interface is effectively a sponge layer in which the turbulent eddies required in the LES region develop. Also, in this region, the integral scale of turbulence and the grid spacing are of the same order: refining the grid reduces both at the same rate and thus cannot improve the results.
Channel flow: time- and spanwise-averaged vertical velocity profile,
${\langle \overline {v} \rangle } /u_b$
, and QUICK blending function,
$f$
. 
, QUICK blending function (right-hand axis); 
, CD grid D; 
, blended grid D; 
, wall-model interface location.
Channel flow,
$Re_\tau =5200$
: mean flow statistics. 
, Lee & Moser (Reference Lee and Moser2015) (DNS); 
, CD grid D; 
, blended grid D; 
, wall-model interface location. (a) Mean velocity,
${\langle \overline {u} \rangle } /u_\tau$
. (b) Reynolds shear stress,
${\langle \overline { {{u}^{\prime }} {{v}^{\prime }} } \rangle } /u_\tau ^2$
. (c) Normal Reynolds stresses,
${\langle \overline {{{u}^{\prime }}{{u}^{\prime }} } \rangle } /u_\tau ^2$
and
$-{\langle \overline {{{v}^{\prime }} {{v}^{\prime }} } \rangle } /u_\tau ^2$
.
A.2. Flow over periodic hills
Periodic hills: flow streamlines (white) and mean streamwise velocity, normalised by the bulk velocity (colour contours).
Next, we performed wall-modelled LES of flow over periodic hills, a benchmark case for heavily separated flows. The Reynolds number was
$Re_h=37\,000$
, and the data are compared with experimental results by Rapp & Manhart (Reference Rapp and Manhart2011). Figure 25 shows the general features of the flow, which separates from the curved surface generating a recirculation bubble. The height of the channel constriction,
$h$
, is one-third of the total channel depth, the distance between consecutive crests is
$L=9h$
and the spanwise extension of the channel is
$W=4.5h$
. Rapp & Manhart (Reference Rapp and Manhart2011) showed that these dimensions are sufficient to allow periodic boundary conditions to be applied. At the top of the domain, we use a free-slip condition.
At the bottom boundary, the approximated wall-model boundary condition is applied as discussed in § 2.1. The wall-model interface was initially set at a fixed distance from the boundary, namely
$\varDelta _n=0.05H^{crest}$
, where
$H^{crest}$
is the water depth at the crest.
We performed a grid-convergence study with three meshes with increasing resolution, with 64, 96 and 128 cells in each direction, respectively. While the grid is uniform in
$x$
and
$z$
, the wall-normal grid spacing varies due to channel constriction. The finest grid (F) has spacing equal to
$\Delta x/h=0.057{-}0.10$
,
$\Delta y/h=0.0015{-}0.038$
and
$\Delta z/h=0.035$
, in the streamwise, wall-normal and spanwise directions, respectively. The standard wall model was used. We then ran two additional cases with the finest grid in which the wall model was modified first by moving the inner-/outer-layer interface close to the wall (F-NW), and then by adding the non-equilibrium correction (F-NW-NEQ), suggested by Meneveau (Reference Meneveau2020).
Flow over periodic hills at
$Re_h=37\,000$
: friction coefficient normalised by the density and bulk velocity,
${\langle \overline {C_f} \rangle _z} (x/h)=2\tau _w/(\rho \widehat {u_b}^2)$
. 
, Experimental data by Rapp & Manhart (Reference Rapp and Manhart2011); 
, F; 
, F-NW; 
, F-NW-NEQ.
Figure 26 shows the time- and span-averaged friction coefficient,
${\langle \overline {C_f} \rangle _z} (x/h)$
, and compares it with the experiments of Rapp & Manhart (Reference Rapp and Manhart2011). We observe significant error for simulation F, in the recirculation region and, to some extent, in the upward slope downstream. To reduce the errors in the recirculation region, we moved the inner-/outer-layer interface to the first fluid point. This modification is justified by the fact that the rationale for having the interface away from the wall is to create a buffer layer, between the wall and the interface, in which eddies can be generated (Kawai & Larsson Reference Kawai and Larsson2012). In a separated flow, however, the existence of an inflectional velocity profile and the inherent unsteadiness of the separation and reattachment lines provide a vigorous stirring of the flow that accelerates the generation of eddies. Moving the interface closer to the wall improves the results significantly (simulation F-NW in figure 26). The addition of the non-equilibrium model correction (Meneveau Reference Meneveau2020) was also beneficial, and the wall shear stress predicted by the F-NW-NEQ simulation agrees much better with the data. Similar conclusions can be drawn from an examination of the mean velocity and Reynolds stresses (not shown). These modifications were used in the production runs described in the paper.
A.3. Morphological model
The final test we carried out had the purpose of testing the accuracy of the morphological model. To this end, we simulated one of the cases studied experimentally by Venditti et al. (Reference Venditti, Church and Bennett2005a
). Those authors considered an open-channel flow where the bottom surface was composed of loose granular material and the bed was initially flat, and analysed five flow regimes, spanning some cases in which sediment transport is minimal, and others in which the sediment begins to be transported immediately, causing bedform deformation. We focus on the experiment at the highest Reynolds number, where the bedform develops most rapidly (FLOW A in Venditti et al. Reference Venditti, Church and Bennett2005a
). Table 4 shows the sediment and flow parameters used to match the experiment. The flow is fully turbulent and subcritical (
$Fr\lt 1$
). The equivalent sand-grain roughness associated with the sediment dimensions was
$k_s=30y_0$
(van Rijn Reference van Rijn1993), and lies at the intersection between the transitional and fully rough regimes. This configuration falls into the gravitational-settling regime (Finn & Li Reference Finn and Li2016), in which the particle motion is mainly driven by settling.
Experimental parameters for the validation of the morphological model (Venditti et al. Reference Venditti, Church and Bennett2005a ).
The simulations match all the flow and sediment parameters. We use the methodology discussed in § 2.1 with the combination of the IBM and the wall model and the blended central difference–QUICK scheme for the convective terms. The domain is a cuboid of length
$L=32H$
, width
$W=4H$
and height
$1.5H$
, and the initial immersed surface in this case is a horizontal plane, located at
$y_b=0.5H$
, such that the effective flow depth is
$H$
. The streamwise and spanwise directions are periodic and the top surface is treated as a fixed free-slip surface. The grid is Cartesian with uniform spacing in the streamwise and spanwise directions. The wall-normal distribution of grid points is similar to that shown in figure 1(b). There is a band of cells with uniform spacing to cover the area across the flow–sediment interface; above this layer the grid is stretched to the top surface. Table 5 shows the numerical parameters. Three grids were used (coarse, medium and fine) to verify the grid convergence of the results.
Open-channel flow over loose granular bed: simulation grid parameters.
Open-channel flow at
$Re_\tau =4590$
, fixed bed. Time- and spanwise-averaged profiles of (a) streamwise velocity,
${\langle \overline {u} \rangle } /u_\tau$
, and (b) Reynolds stress profiles in wall units,
${\langle \overline {{{u}^{\prime }} {{v}^{\prime }} } \rangle }$
and
$ {\langle \overline {{{u}^{\prime }} {{u}^{\prime }} } \rangle }$
. 
, Venditti et al. (Reference Venditti, Church and Bennett2005a
); 
, coarse grid; 
, medium grid; 
, fine grid; 
, wall-model interface location.
Initially, the simulation is run with a fixed bed to allow the flow to develop fully and reach a statistically steady state. Then, statistics for the fixed-bed case were collected for
$150T_f$
, where
$T_f=\delta /u_\tau$
is the large-eddy turnover time. Experimental measurements collected before the bed started to deform are available, and are compared in figure 27 with the numerical results. Since the medium and fine grids are in good agreement with each other and with the experimental data, the medium grid is used in the following.
The simulation starts from an instantaneous flow field obtained from the fixed-bed case. We use the bed-load transport model by van Rijn (Reference van Rijn1984b
), which was shown (D’Alessandro et al. Reference D’Alessandro, Hantsis, Marchioli and Piomelli2021) to yield accurate predictions in this transport regime
$ ({\langle \overline {\theta } \rangle } \simeq 0.1\gt \theta _{cr} )$
. Starting from the flat bed, Venditti et al. (Reference Venditti, Church and Bennett2005a
) reported the emergence of several, progressively larger and more complex features on the bed:
-
(i) Longitudinal striations with spacing of the order of the hairpin streaks, which form and disappear quickly.
-
(ii) A cross-hatch pattern (
$t/T_s=6720$
). -
(iii) Chevron-shaped features develop at the nodes of the cross-hatch (
$t/T_s=8400$
). -
(iv) Incipient crest lines.
-
(v) Quasi-two-dimensional crest lines for
$15\,120\lt t/T_s\lt 50\,400$
. -
(vi) Quasi-equilibrium two-dimensional bedforms for
$t/T_s=2\times 10^6$
.
Contours of the bedform-elevation displacement
$\widehat {y_b}$
. Five time instants from flat bed to quasi-two-dimensional transverse dunes are shown at (a)
$t/T_s=0$
, (b) 12 000, (c) 15 000, (d) 35 000 and (e) 50 000. Only a portion of the domain is shown for clarity.
The structures corresponding to points (i)–(iii) are too small to be captured by a wall-modelled LES grid and, therefore, we focus our attention on the structures corresponding to points (iv)–(vi), namely on the later stages of development.
Figure 28 shows contours of the bed elevation displacement,
$\widehat {y_b}$
, normalised by the grain diameter,
$d$
, at five time instants:
$t/T_s=0$
, 12 000, 15 000, 35 000 and 50 000, corresponding to the phase of development of quasi-two-dimensional crest lines. The bed-elevation displacement from the initial flat bed is defined as
$\widehat {y_b}(x,z,t)=y_b(x,z,t)- {\langle {y_b(t=0)}\rangle }$
. Supplementary movie 2 shows the full evolution of the bedforms.
Figure 28(a) shows the initial flat bed with
$\widehat {y_b}=0$
, while figure 28(b–e) shows the elevation variation during the formation and growth of the two-dimensional crest lines. Compared with the experimental findings, our bedform field exhibits a greater spanwise variation (i.e. sinuosity). Similar results were observed in the numerical simulations by Khosronejad & Sotiropoulos (Reference Khosronejad and Sotiropoulos2014) which also reproduced the configuration of the laboratory experiment conducted by Venditti et al. (Reference Venditti, Church and Bennett2005a
). They suggested that this issue might be related to the different inflow boundary conditions and the length of the domain compared with the experiment. In this work, instead of inlet conditions, we adopt periodic boundary conditions in the streamwise direction, which recirculates all the initial three-dimensional bed structures.
To compare quantitatively the bedform evolution to the experiment, we determined the wavelength and amplitude of the bedforms, by smoothing the bed profiles and identifying the local maxima. Figure 29 shows the wavelength,
$\lambda (z,t)$
, and the height,
$\varDelta (z,t)$
, of the bedforms. Each symbol (circle and star) represents the average wavelength and height of a profile at a specified spanwise location and time instant. The grey area represents the 90 % confidence range of the experimental measurements in figure 29. The numerical results are slightly below the average height and wavelength but fall well within the uncertainty of the experiments. Determining the height and wavelength of the ripples is challenging both experimentally and numerically, leading to significant variability in the experimental (and numerical) data. Given these factors, we find the agreement satisfactory.
Channel flow over loose granular bed: comparison of (a) bedform wavelength,
$\lambda$
, and (b) bedform height,
$\varDelta$
, with experimental data. 
, Venditti et al. (Reference Venditti, Church and Bennett2005a
); 
,
$90\,\%$
confidence range; 
, best-fit law by Venditti et al. (Reference Venditti, Church and Bennett2005a
); 
, coarse grid; 
, medium grid. Computed data at multiple spanwise locations are shown for each time instant.
As anticipated in § 2.2, we neglect the suspended-sediment transport in the present simulations. Khosronejad & Sotiropoulos (Reference Khosronejad and Sotiropoulos2014) performed a numerical simulation with the same flow and sediment configuration and a numerical model similar to the one used here, including suspended-sediment transport and stratification effects. They observe that fine grids are required to achieve bed destabilisation unless stratification is included, and that suspended sediment accounts only for about 17 % of the total sediment transport rate.
In the present work, the sand bed is destabilised even with the coarse grid, probably because of the less dissipative nature of the numerical scheme compared with that of Khosronejad & Sotiropoulos (Reference Khosronejad and Sotiropoulos2014). Furthermore, the agreement of the bedform height and wavelength indicates that the role of suspended-sediment transport, in the present simulations, is not crucial. The results shown in this appendix prove that the methodology applied in this work can effectively model strongly separated flows, and can accurately predict the morphodynamics of a mobile bed.
Open access
