Large-eddy simulation of the flow over a realistic riverine geometry

Abstract

In this paper, we discuss the transport of sediment and the formation of bedforms in turbulent river flows, under flow conditions typical of flooding events. Through the implementation of an immersed boundary method, a wall model and a morphological model, we were able to simulate complex and mobile geometries under high Reynolds numbers at an affordable computational cost. In particular, we examined the evolution of bedforms on a loose sediment bed under turbulent flow conditions, using input parameters obtained from laboratory measurements. Over time, the bedforms become more three-dimensional and irregular in shape, leading to changes in the shear layer, crest angle and separation patterns. The bedforms continue to evolve until a quasi-steady equilibrium is reached. Our simulations highlight the crucial role played by the small-scale bedforms, which significantly affect the flow dynamics: an increase in the total drag is observed, related to the form drag generated by the local recirculation and the increased size of the large-scale recirculation bubble. Furthermore, a stronger turbulent activity ensues from the shear layers forming on the crests of the small-scale bedforms. Finally, a wider shedding angle of the shear layer is caused by the irregular crest line.

1. Introduction

According to global numerical predictions, the frequency of extreme flooding will double across 40 % of the globe by 2050 (Arnell & Gosling 2016). 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. 2008) 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 1996, 2005; Venditti et al. 2005b ). 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. 2005a ; Andreotti et al. 2011; Coleman & Nikora 2011), formation process (Williams & Kemp 1971; Best 1992; Blondeaux et al. 2016) and migration patterns (Mohrig & Smith 1996; Charru et al. 2013). 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. 2005a ; Coleman & Nikora 2011). 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. 2018). 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. (2018) 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. 2005a ). Both analytical (Andreotti et al. 2011; Colombini & Stocchino 2011) and experimental (Coleman & Melville 1996) 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. 2010). These interactions generate complex flow fields (Best 2005b ).

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 2011) involving near-wall turbulence and localised imperfections and for the mechanism of dune migration (Venditti et al. 2005a ). 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 2011). 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. 2022). However, Ashley (1990) 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. 2020). Additionally, field measurements of the sediment flux during morphogenetic events are difficult to obtain (Garcia 2008).

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 2017). 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. 2004; Lefebvre 2019; Yamaguchi et al. 2019). Parsons et al. (2004) focused on idealised two-dimensional transverse dunes, while Lefebvre (2019) and Yamaguchi et al. (2019) studied natural three-dimensional bedforms. More recently, Chiodi et al. (2014) and Ahadi et al. (2018) 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. (2024) 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 2019; Salimi-Tarazouj et al. 2024). However, it fails to explain some of the hydrodynamic and morphodynamic mechanisms due to the parametrisation of both large- and small-scale eddies (Piomelli 2008; Khosronejad et al. 2011; Coleman et al. 2015; Wu et al. 2017).

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 2014, 2017; Akiki & Balachandar 2020; Mazzuoli et al. 2020; Vittori et al. 2020).

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. 2006; Vowinckel 2021). A LES with Lagrangian models also requires additional modelling for sub-filter-scale advection velocity (Marchioli 2017).

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 2001, 2006; Khosronejad & Sotiropoulos 2014, 2017). Chou & Fringer (2008) employed LES to examine sediment-transport time scales in turbulent flow and later studied sand-ripple formation (Chou & Fringer 2010), focusing on fixed beds. Khosronejad & Sotiropoulos (2014) developed a hydromorphodynamic model for sand-wave formation in open-channel flow, later exploring Barchan dune formation (Khosronejad & Sotiropoulos 2017). 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. (2018), 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

(2.1a) \begin{align} &\qquad\qquad\qquad\qquad\qquad\qquad\boldsymbol \nabla \cdot {\tilde {\mathbf u}} =0 \;, \end{align}

(2.1b) \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 1963) 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. 1995). 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. 2006). A fractional step method (Kim & Moin 1985) 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 1987). 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:

(2.2) \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. (2018) 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. 2013). The intrinsic volume average is defined as

(2.3) \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

(2.4) \begin{equation} \Omega _f =\int _{\Omega _0}\phi (\mathbf{x},t) {d}\Omega \; . \end{equation}

Analogously, the intrinsic spanwise average is defined as

(2.5) \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

(2.6) \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

(2.7) \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

(2.8) \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.

Figure 1. (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 (2006). 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 ):

(2.9) \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

(2.10) \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 (2020). 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

(2.11) \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

(2.12) \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 1993). Additional details about the model can be found in the original paper (Meneveau 2020).

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:

(2.13) \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 1999; García 2008).

2.2.1. Sediment mass conservation

The transport equation for the sediment phase is based on the Exner–Polya equation (Exner 1920), revisited by Paola & Voller (2005). The equation reads as follows:

(2.14) \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 1935; Martin & Aral 1971; Wheatcroft 2002), ${{\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:

(2.15) \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.

Figure 2. 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 1936; van Rijn 1984a ; Nielsen 1992; Niño et al. 2003; Debnath & Chaudhuri 2010). 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 1987; García 2008). 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}$ :

(2.16) \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. 2013; Khosronejad et al. 2011). 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 (1987):

(2.17) \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

(2.18) \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

(2.19) \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 (1993), and corrected to take into account the influence of gravity on erosion from sloping beds with the following correction factor:

(2.20) \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. (2011). 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 2018). 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 2017); 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

(2.21) \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

(2.22) \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:

1. (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).

2. (ii) The wall shear stress, $\boldsymbol{\tau _w^{\textit{WM}}}$ , is calculated from (2.12).

3. (iii) The wall-model forcing term, ${{\mathbf F}^{\textit{WM}}}$ , is calculated from (2.13).

4. (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).

5. (v) The bed-load sediment transport, ${{\mathbf q}_{\mathbf{bl}}}$ , is calculated using $\boldsymbol{\tau _w^{\textit{WM}}}$ and (2.16).

6. (vi) The bed elevation, ${{\mathbf y}_{\mathbf{b}}}$ , is updated using (2.15).

7. (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. (2018), 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. (2018) 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 (2008), Omidyeganeh & Piomelli (2011, 2013) and Hardy et al. (2021), 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.

Table 1. 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:

(2.23) \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}

(2.24) \begin{align} &\qquad\qquad\widehat {H}=\frac {\int _{V_0}\phi (\mathbf{x},t)\textrm{d}V}{WL} \;, \end{align}

(2.25) \begin{align} &\qquad\qquad u_b(x,t)= {\left \langle u(\mathbf{x},t) \right \rangle _{yz}} \;, \end{align}

(2.26) \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 2004; Schultz & Flack 2007).

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 2012), the grid-resolution study was performed with a fixed wall-model interface location, $\varDelta _n=0.05h$ .

Table 2. 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:

(2.27) \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.

Figure 3. Time- and spanwise-averaged friction coefficient, ${\langle \overline {C_f} \rangle _z}$ . , Coarse grid; , medium grid; , fine grid.

Figure 4. 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. (2018), we study a configuration where the bed presents both large morphological structures and superimposed bedforms.

3.1. Bed geometry

Figure 5. 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. (2005a ). 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:

(3.1) \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)$ .

Figure 6. 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.

Figure 7. 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.

Figure 8. 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}$ :

(3.2) \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

(3.3) \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:

(3.4a) \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}

(3.4b) \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.

Figure 9. (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.

Figure 10. (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. (2018), where $\varDelta /d\approx 45$ and $\lambda /d=500{-}1200$ .

Figure 11. 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).

Figure 12. (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.

Figure 13. 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

(3.5a) \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}

(3.5b) \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.

Table 3. 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.

Figure 14. 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:

(3.6) \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. (2021). As observed by Zgheib et al. (2018), 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.

Figure 15. 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.

Figure 16. 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.

Figure 17. 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).

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.

Figure 19. 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:

(3.7) \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 1986), significant APG and FPG effects can be noticed.

Figure 20. 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 2013), and an increase of the gradient in APG regions (Knopp et al. 2021). 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 1995; Pope 2000), 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 2020). 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

(3.8) \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.

Figure 21. (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.

Figure 22. (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:

(3.9) \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 1984a ). 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. (2018) revealed the relation between bedforms and flow parameters such as flow depth and flow regime. Reesink et al. (2018) 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. (2018) 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.