Hostname: page-component-7479d7b7d-m9pkr Total loading time: 0 Render date: 2024-07-11T22:34:04.806Z Has data issue: false hasContentIssue false

Particle-resolved simulation of antidunes in free-surface flows

Published online by Cambridge University Press:  24 April 2023

Christoph Schwarzmeier*
Affiliation:
Chair for System Simulation, Friedrich-Alexander-Universität Erlangen-Nürnberg, Cauerstraße 11, 91058 Erlangen, Germany
Christoph Rettinger
Affiliation:
Chair for System Simulation, Friedrich-Alexander-Universität Erlangen-Nürnberg, Cauerstraße 11, 91058 Erlangen, Germany
Samuel Kemmler
Affiliation:
Chair for System Simulation, Friedrich-Alexander-Universität Erlangen-Nürnberg, Cauerstraße 11, 91058 Erlangen, Germany
Jonas Plewinski
Affiliation:
Chair for System Simulation, Friedrich-Alexander-Universität Erlangen-Nürnberg, Cauerstraße 11, 91058 Erlangen, Germany
Francisco Núñez-González
Affiliation:
Department of Civil and Environmental Engineering, Universitat Politécnica de Catalunya – BarcelonaTech (UPC), Jordi Girona 1-3, 08034 Barcelona, Spain
Harald Köstler
Affiliation:
Chair for System Simulation, Friedrich-Alexander-Universität Erlangen-Nürnberg, Cauerstraße 11, 91058 Erlangen, Germany
Ulrich Rüde
Affiliation:
Chair for System Simulation, Friedrich-Alexander-Universität Erlangen-Nürnberg, Cauerstraße 11, 91058 Erlangen, Germany CERFACS, 42 Avenue Gaspard Coriolis, 31057 Toulouse CEDEX 1, France
Bernhard Vowinckel
Affiliation:
Leichtweiß-Institute for Hydraulic Engineering and Water Resources, Technische Universität Braunschweig, 38106 Braunschweig, Germany
*
Email address for correspondence: christoph.schwarzmeier@fau.de

Abstract

The interaction of supercritical turbulent flows with granular sediment beds is challenging to study both experimentally and numerically. This challenging task has hampered advances in understanding antidunes, the most characteristic bedform of supercritical flows. This article presents the first numerical attempt to simulate upstream-migrating antidunes with geometrically resolved particles and a liquid–gas interface. Our simulations provide data at a resolution higher than laboratory experiments, and they can therefore provide new insights into the mechanisms of antidune migration and contribute to a deeper understanding of the underlying physics. To manage the simulations’ computational costs and physical complexity, we employ the cumulant lattice Boltzmann method in conjunction with a discrete element method for particle interactions, as well as a volume-of-fluid scheme to track the deformable free surface of the fluid. By reproducing two flow configurations of previous experiments (Pascal et al., Earth Surf. Process. Landf., vol. 46, issue 9, 2021, pp. 1750–1765), we demonstrate that our approach is robust and accurately predicts the antidunes’ amplitude, wavelength and celerity. Furthermore, the simulated wall shear stress, a key parameter governing sediment transport, is in excellent agreement with the experimental measurements. The highly resolved data of fluid and particle motion from our simulation approach open new perspectives for detailed studies of morphodynamics in shallow supercritical flows.

Type
JFM Rapids
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2023. Published by Cambridge University Press.

1. Introduction

Antidunes are short-wave periodic disturbances that develop on the surface of loose granular beds in response to the interaction with supercritical and near-critical shallow turbulent flows (Kennedy Reference Kennedy1963), where criticality is defined by the Froude number. These types of bedforms arise in fluvial, coastal and submarine environments and are closely tied to the resulting flow resistance, turbulence and sediment transport, which are quantities that are of great interest in geology as well as in hydraulic and environmental engineering. Among the bedforms at supercritical flow conditions, antidunes are the shortest in wavelength (of the order of the flow thickness), and they are also the only bedforms that can migrate either upstream or downstream, or even remain stationary. This work deals with upstream-migrating antidunes (UMAs). The longitudinal profile of a UMA is symmetrical, exhibiting a sinusoidal shape that is mirrored by the undulating free water surface. Consequently, UMAs are coupled to the fluid surface and are, thus, dependent on the presence of a free surface or an interface where stationary waves can develop. The movement of UMAs in the direction opposite to the main flow is counterintuitive, and it is believed to result from the rhythmic aggradation of particles on the upstream flank of the bedform and erosion on the lee side. However, little is known about the migration mechanism in connection to turbulence, bed morphology and sediment transport.

Theoretical studies have made significant advances to understanding the characteristics of antidunes, by defining antidune instability regions (e.g. Kennedy Reference Kennedy1963; Colombini & Stocchino Reference Colombini and Stocchino2012; Bohorquez et al. Reference Bohorquez, Cañada-Pereira, Jimenez-Ruiz and del Moral-Erencia2019) and identifying hydraulic constraints for antidune migration regimes (Núñez-González & Martín-Vide Reference Núñez-González and Martín-Vide2011). Unfortunately, due to the challenging supercritical flow conditions often associated with low submergences, experimental datasets to validate those theories are extremely sparse. The limited number of systematic experimental datasets with antidune configurations is in part related to the inherent technical challenges in reproducing rapid flows over an erodible bed in laboratory flumes (Slootman et al. Reference Slootman, Ventra, Cartigny, Normandeau and Hubbard2021), as well as to the difficulties in performing non-intrusive measurements in flows that are generally unstable and shallow. Numerical simulations of supercritical flows above an erodible bed, as a methodological alternative for the study of antidunes, are challenging because of the strongly varying water surface of rapid flows. Existing simulations of supercritical bedforms have mainly considered the Saint-Venant shallow-water equations (e.g. Parker & Izumi Reference Parker and Izumi2000; Balmforth & Vakil Reference Balmforth and Vakil2012) and the Reynolds-averaged Navier–Stokes approach to model cyclic steps (Vellinga et al. Reference Vellinga, Cartigny, Eggenhuisen and Hansen2018) or downstream-migrating antidunes (Olsen Reference Olsen2017, Reference Olsen2022), but, to the best of the authors’ knowledge, simulations of UMAs have not been reported so far.

Recently, Kidanemariam & Uhlmann (Reference Kidanemariam and Uhlmann2014) have successfully simulated downstream-migrating ripples in subcritical flows, with particle-resolved direct numerical simulations in which the fluid–particle interaction is resolved by considering the motion of the flow and of each individual particle. Owing to the subcritical flow conditions in their study, Kidanemariam & Uhlmann (Reference Kidanemariam and Uhlmann2014) could simplify the flow geometry by keeping the fluid depth constant. In this work we propose to use particle-resolved simulations in conjunction with a deformable fluid surface, to simulate the formation and propagation of UMAs in supercritical flows. We aim to numerically reproduce an experimental campaign recently reported by Pascal, Ancey & Bohorquez (Reference Pascal, Ancey and Bohorquez2021), who managed to measure the propagation of UMAs with high spatial and temporal resolution. For this, we use the lattice Boltzmann method (LBM) to simulate the fluid hydrodynamics (Rettinger & Rüde Reference Rettinger and Rüde2022) and extend it with a volume-of-fluid scheme (Schwarzmeier et al. Reference Schwarzmeier, Holzer, Mitchell, Lehmann, Häusl and Rüde2023) to track the strongly deformable free fluid surface in supercritical flows. In turn, we couple the fluid hydrodynamics with the particle dynamics via the momentum-exchange method to simulate the antidune formation with high fidelity. The parameter choices of Pascal et al. (Reference Pascal, Ancey and Bohorquez2021), with coarse sediment grains and low particle relative submergence, allow for a direct overlap of experimental conditions with simulations. In this manner, the present work closes the gap between river morphodynamics’ experiments and numerical simulations, and demonstrates the capabilities offered by large-scale simulations of a mobile sediment bed comprising thousands of particles in unidirectional supercritical turbulent flows, to study coupled bedform and free-surface interactions in great detail.

2. Numerical methods

All numerical simulations in this work were performed using the open-source multi-physics software framework waLBerla (https://www.walberla.net/; Bauer et al. Reference Bauer2021). In what follows, we briefly summarize the simulation's key features, namely the fluid–gas, the particle–particle and the fluid–particle interactions. Figure 1(a) presents a sketch of our numerical scheme.

Figure 1. (a) Schematic representation of the simulations, coupling the liquid (l) with a particle (p) and gas (g) phase, and (b) the simulation set-up in its initial condition. Panel (b) is a zoom into the computational domain covering only 16 % of its streamwise extent.

2.1. Free-surface lattice Boltzmann method

The LBM (Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017) discretizes the Boltzmann equation from kinetic theory and describes the evolution of particle distribution functions on a uniform computational grid to simulate a fluid flow. Macroscopically, the LBM approximates the Navier–Stokes equations with second-order accuracy in space $\boldsymbol {x}$ and time $t$. The cell-local fluid density $\rho _{{l}}(\boldsymbol {x},t)$ and velocity $\boldsymbol {u}_{{l}}(\boldsymbol {x},t)$ are given by the zeroth- and first-order moments of the cell's particle distribution functions. In the remainder of this article, all quantities denoted with the superscript LU are specified in a normalized LBM unit system. We set the cell size $\Delta x^{{LU}} = 1$, time step size $\Delta t^{{LU}} = 1$ and the reference density of the fluid $\rho ^{{LU}}_{{l},{ref}}=1$. We use the D3Q27 cumulant collision model proposed by Geier et al. (Reference Geier, Schönherr, Pasquali and Krafczyk2015) and set all relaxation rates to unity, except for $\omega ^{{LU}}$, which defines the kinematic viscosity of the fluid, $\nu _{{l}}^{{LU}} = (1/3) (1/\omega ^{{LU}} - \Delta t^{{LU}}/2)$. We model no-slip boundary conditions at solid obstacles using the bounce-back rule (Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017).

The waLBerla framework includes an implementation of the free-surface lattice Boltzmann method (Schwarzmeier et al. Reference Schwarzmeier, Holzer, Mitchell, Lehmann, Häusl and Rüde2023, FSLBM,). The FSLBM is based on a volume-of-fluid scheme and simulates the interface between two immiscible fluids. The model assumes that the liquid fluid (higher density) governs the flow dynamics of the system so that the fluid flow in the gaseous fluid (lower density) is neglected. A free-surface boundary condition captures the motion of the liquid–gas interface. The boundary condition incorporates the gas pressure $p_{{g}}$ and the Laplace pressure $p_{{l},{L}}(\boldsymbol {x},t)=2\sigma _{{l}}\kappa _{{l}}(\boldsymbol {x},t)$, where $\sigma _{{l}}$ is the surface tension of the liquid and $\kappa _{{l}}(\boldsymbol {x},t)$ is the local curvature of the free surface.

2.2. Particle–particle and fluid–particle coupling

The temporal evolutions of the translational and rotational velocities, $\boldsymbol {u}_{{p},i}$ and $\boldsymbol {\omega }_{{p},i}$, of an individual rigid particle $i$ are described by the Newton–Euler equations for rigid-body motion. The total force and torque acting on a particle comprise contributions from inter-particle collisions, hydrodynamic interactions and an external force attributed to the gravitational acceleration $\boldsymbol {g}$. The inter-particle collisions are modelled via the discrete element method (DEM) specifically designed for a four-way coupled lattice Boltzmann method of particle-laden flows. The model is fully parametrized by the dry coefficient of restitution $e_{{dry}} = 0.97$, the collision time $T_{{c}}^{{LU}} = 4 d_{50}^{{LU}} \Delta t^{{LU}} /\Delta x^{{LU}}$, Poisson's ratio $\nu _{{p}} = 0.22$ and the friction coefficient $\mu _{{p}}=0.5$ to simulate silica beads. Details on the particle–particle and fluid–particle coupling can be found in Rettinger & Rüde (Reference Rettinger and Rüde2022).

For fluid–particle interactions, we employ a geometrically resolved coupling between the fluid and the particulate phase. We impose no-slip boundary conditions for the fluid along the particle surface and compute the hydrodynamic forces and torques acting on the particles directly. This geometrically resolved coupling is achieved by the LBM-specific momentum-exchange method (MEM). The MEM requires an explicit mapping of the particles onto the underlying computational grid (Aidun, Lu & Ding Reference Aidun, Lu and Ding1998; Wen et al. Reference Wen, Zhang, Tu, Wang and Fang2014). Following Rettinger & Rüde (Reference Rettinger and Rüde2022), we additionally augment the DEM with lubrication corrections in the normal and tangential directions to model these strong subgrid-scale hydrodynamic forces for situations where gaps between particles become too small to be properly resolved by the LBM mesh.

3. Simulation set-up

3.1. Physical scenario and simulation parametrization

The simulations carried out for this article aim to reproduce the laboratory experiments of Pascal et al. (Reference Pascal, Ancey and Bohorquez2021), who performed their experiments with water and a sediment bed of natural gravel. The gravel particles had a density of $\rho _{{p}} = 2550~{\rm kg}~{\rm m}^{3}$ and characteristic sieving diameters of $d_{16}=2.5$ mm, $d_{50}=2.9$ mm and $d_{84}=3.3$ mm, where the subscripts denote the respective percentile of the particle size distribution. The mean flow velocity $U_{x,{l}} = q_{{w}} / h_{0}$ is computed from the measured flow discharge in the surface flow layer $q_{{w}}$ and reference flow depth $h_{0}$ (figure 1b). Three dimensionless numbers specify the flow conditions: the Reynolds number ${Re} = U_{x,{l}} h_{0} / \nu _{{l}}$, the Froude number ${Fr} = U_{x,{l}}/\sqrt {g h_{0}}$ and the Weber number ${We} = \rho _{{l}} U_{x,{l}}^{2} h_{0} / \sigma _{{l}}$, where $\nu _{{l}} = 1 \times 10^{-6}~{\rm m}^{2}~{\rm s}^{-1}$ and $\sigma _{{l}} = 7.2 \times 10^{-2}~{\rm kg}~{\rm s}^{-2}$ are the kinematic viscosity and surface tension of water, respectively. The values of these are listed in table 1, where we refer to the simulation set-ups with the same labels, E1 and E4, as chosen by Pascal et al. (Reference Pascal, Ancey and Bohorquez2021).

Table 1. Summary of the main simulation conditions. The labels are chosen as in the reference experiments from Pascal et al. (Reference Pascal, Ancey and Bohorquez2021).

The simulation domain was discretized with a mesh size $\Delta x = 2.5 \times 10^{-4}$ m on a uniform computational grid of size $L^{{LU}}_x \times L^{{LU}}_y \times L^{{LU}}_z = 3200 \times 60 \times 160$ cells, corresponding to $0.8~\text {m} \times 0.015~\text {m} \times 0.04~\text {m}$ in the streamwise, spanwise and vertical directions, respectively. We have used periodic boundary conditions in the $x$- and $y$-directions. The boundaries in the $z$-direction and on the particle surface were modelled with a no-slip condition. The sediment bed consisted of $6528$ spherical particles with a mean diameter of $d^{{LU}}_{50} = 11.6$, where the actual diameters were sampled from a uniform distribution of the interval $[0.9,1.1] d_{50}$. We have generated the initially flat particle bed of height $h_{{b}} = 4.14 d_{50}$ using a rigid-body dynamics simulation for a sedimentation process in dry conditions. To this end, particles were dropped onto a single layer of fixed particles that were positioned on a horizontal hexagonal grid with small random perturbations in their $z$-positions. We fix all particles with a vertical centre coordinate $z_{{p},i}< h_{{b,f}} = 1.5 d_{50}$ to generate a rough bottom. We have carried out sensitivity studies to verify that the chosen number of particles and $h_{{b,f}}$ were sufficient to exclude finite-size effects of the domain extent.

The liquid was initialized with a height of $h_{{b}}+h_{0}$, where $h_{0}$ represents the initial height of the liquid on top of the sediment bed, and with a parabolic velocity profile on top of the bed. The hydrostatic density was initialized, such that the liquid pressure was equal to the constant atmospheric gas pressure at the free surface. The flow was driven in the $x$-direction via an externally imposed body force $F_x$, which was controlled during the simulation to maintain the desired bulk flow rate. The mean fluid velocity was set to $U^{{LU}}_{x,{l}} = 0.02$, leading to the time-step sizes $\Delta t$ specified in table 1.

3.2. Evaluation quantities and grid resolution study

The temporal and spatial variability of the bed topography is one of the key quantities to characterize the evolution of antidunes. We follow the approach of Kidanemariam & Uhlmann (Reference Kidanemariam and Uhlmann2014) and compute the $y$-averaged solid volume fraction $\phi _{{b}}(x,z,t)$ for each computational cell in the $x,z$-plane. The bed height $h_{{b}}(x,t)$ is then given as the linearly interpolated $z$-value, for which $\phi _{{b}}(x,z,t) = 0.3$. The bedload transport rate $q_{{b}}(t)$ is computed via

(3.1)\begin{equation} q_{{b}}(t)=\frac{\sum_{i \in {P}} u_{x,{p},i}(t) V_{{p},i}}{L_{x}L_{y}}, \end{equation}

where $P$ is the set of all particles and $u_{x,{p},i}(t)$ is the streamwise particle velocity of particle $i$.

We carried out a grid sensitivity study to verify that our results are independent of numerical parameters. To this end, we focused on the more turbulent flow conditions of E4 (see table 1). In addition to the reference case from § 3.1 with $\Delta x_{50,{medium}} = 2.5 \times 10^{-4}$ m as a medium resolution, we have tested two additional cases with a two-times coarser ($\Delta x_{50,{coarse}} = 5 \times 10^{-4}$ m) and finer resolution ($\Delta x_{50,{fine}} = 1.25 \times 10^{-4}$ m). These computational grid resolutions lead to $d^{{LU}}_{50,{coarse}} = 5.8$, $d^{{LU}}_{50,{medium}} = 11.6$ and $d^{{LU}}_{50,{fine}} = 23.2$. For the grid sensitivity study, we shortened the simulation domain by a factor of $10$ in the $x$-direction to keep the computational cost manageable. We have scaled the time-step size $\Delta t$ proportionally to $\Delta x$ in our study with reference $\Delta t_{50,{medium}} = 1.08 \times 10^{-5}$ s, as listed in table 1. Figure 2 compares the bedload transport rates (3.1) for the different grid resolutions. From this study, we conclude that the ‘medium’ grid is sufficiently independent of the grid resolution and can be used for further analysis.

Figure 2. Simulated bedload transport rate $q_{{b}}(t)$ with different computational grid resolutions using a shortened simulation domain. The time averages were computed only in the time ranges conforming to the lengths of the plotted lines to ensure that the fluid flow was adequately developed.

4. Results

Our particle-resolved simulations provide the position of each particle as well as the location of the water surface in space and time, and full fluid information (figure 3). Bed elevation perturbations arose spontaneously from the initial flat bed right from the start of the simulations. From then on, a regular longitudinal pattern of quasi-periodic bedforms, with gentle slopes up- and downstream of the crests, prevailed throughout the simulation time. The bed patterns were approximately in phase with respect to the free surface and they slowly migrated upstream (see animations in the supplementary movies available at https://doi.org/10.1017/jfm.2023.262). Therefore, these bedforms can be unequivocally classified as UMAs that compare well to the observations of the experiments of Pascal et al. (Reference Pascal, Ancey and Bohorquez2021). A dataset with such a high resolution as the one generated here numerically has so far not been available from physical experiments on movable beds under supercritical flows.

Figure 3. Visualization of the simulated velocity in streamwise direction $u_{x,{l}}(\boldsymbol {x},t)$. The undulations of the sediment bed and free surface of the liquid are in phase, conforming to the definition of antidunes. The figure is a zoom into the computational domain covering only $25$ % of its streamwise extent.

4.1. Bed elevation perturbations

Space–time evolution diagrams of the bed surface, for the simulated and the experimental bedforms, are plotted in figure 4. Such diagrams result from laterally averaged bed surface profiles at every 1000 time steps of the simulation, and from footage recorded from the transparent channel walls in the experiments of Pascal et al. (Reference Pascal, Ancey and Bohorquez2021). From the original experimental data, we have selected the intervals corresponding to the lower left parts of figure 3(a,d) in that reference, to reproduce a part of their bed elevation plots to compare with the scales of our numerical data. The alternating blue and yellow (dark and light) diagonal strips denote the troughs and crests of the bedforms, respectively. The experimental and numerical diagrams with the evolution of the bottom elevation in the spatial–temporal domain are in good qualitative agreement. A similar upstream migration trend of the bedforms is clearly visible due to the negative slope of the strips. In some regions the strips bend, indicating acceleration/deceleration, or even stationarity of the bedforms. Particularly for run E1, these regions are related to local perturbations that migrate downstream. Note that similar perturbations can be identified in both experimental and simulation plots. Overall, the simulated bedforms appear stable and compare well to their experimental counterparts.

Figure 4. Sediment bed elevation $h_{{b}}(x,t)$ for E1 (a,b) and E4 (c,d): (a,c) experimental data from Pascal et al. (Reference Pascal, Ancey and Bohorquez2021), and (b,d) simulation data from the numerical simulations. We consider the system fully developed after $t=30$ s, as illustrated by the dashed black line.

As pointed out by Pascal et al. (Reference Pascal, Ancey and Bohorquez2021), the definition of a dominant bedform migration speed seems inappropriate given the non-uniformity of the bedform celerities observed. The same reasoning is valid for the amplitude and the wavelength of the antidunes. In our simulations, the amplitude of the simulated bed undulations ranged approximately from $1$ to $3$ times the particle median grain size $d_{50}$, and the wavelength ranged from approximately $10$ to $15$ times the water depth $h_0$. These values are in good quantitative agreement with the experimental data, where the dune amplitudes ranged from one particle median grain size to the mean flow depth ($\approx 3d_{{50}}$), and typical wavelengths varied from $0.05$ to $0.15$ m ($\approx 6h_{0}$ to $15h_{0}$).

To obtain a more precise comparison of the numerical and experimental data in terms of their bedform size and fluctuations, we have followed Pascal et al. (Reference Pascal, Ancey and Bohorquez2021) and computed the power spectral density (PSD) from the square of the two-dimensional discrete Fourier transform of $h_{{b}}(x,t)$ with respect to the streamwise position $x$ and time $t$. We have normalized the PSD by the total number of samples available for $x$ and $t$. The results are shown in figure 5, along with the experimental PSD. Although there is not a perfect match between the spectra for the experimental and numerical bedforms, the range of wavelengths and periods agree reasonably well. For the numerical data, the ranges of the periods and wavelengths in the spectra are slightly narrower than in the experiments, which can be attributed to the much shorter data series from the simulations ($45$ s simulation time versus more than $4600$ s in E1 and $2800$ s in E4 in the experiments).

Figure 5. PSD of wavelength $\lambda$ and period $T$ for (a) E1 and (b) E4 from Pascal et al. (Reference Pascal, Ancey and Bohorquez2021). The orange line is the outermost hull of the simulation results’ PSD.

We compared experimental and simulated results further by obtaining the spectra in the celerity–wavelength domain (figure 6). The celerity $c = \lambda /T$, that is, the movement speed of the dunes, was computed as the ratio between wavelength $\lambda$ and period $T$ from the spectra in figure 5. Similarly to the PSD in the $\lambda$$T$ plane, a fairly good agreement was found between the simulations and the experiments. More precisely, the celerities in the simulations range from $2$ to $12~{\rm mm}~{\rm s}^{-1}$ in E1 (experiment: $\approx 2\unicode{x2013}15~{\rm mm}~{\rm s}^{-1}$) and from $2$ to $18~{\rm mm}~{\rm s}^{-1}$ in E4 (experiment: $\approx 5\unicode{x2013}30~{\rm mm}~{\rm s}^{-1}$). It is important to note that the simulations reproduced the trends observed experimentally, that, for a given wavelength, celerity increases with bedload transport rate; while, for a constant bedload transport rate, celerity increases with increasing antidune wavelength. This latter trend is opposite to that commonly observed for morphodynamics in subcritical flows (e.g. Coleman & Nikora Reference Coleman and Nikora2011).

Figure 6. PSD contour of celerity $c$ and wavelength $\lambda$ for (a) E1 and (b) E4 from Pascal et al. (Reference Pascal, Ancey and Bohorquez2021).

4.2. Hydraulic variables and sediment transport

To further compare the experimental and simulated flow and sediment transport conditions in a quantitative sense, the normalized average water depth, the bottom slope, the Shields number and the Einstein bedload number are defined and summarized in table 2. It must be noted that the bottom profile in the simulations was intentionally horizontal, and thus the slope $\psi$ considered in table 2 is the tangent to the acting force imposed in $x$$z$ direction. Assuming a steady uniform flow where the fluid gravitational force is balanced by boundary friction, the average bed shear stress was computed with the depth–slope product $\rho _1 g R_{{b}}\tan \psi$, where $R_{b}$ is the hydraulic radius.

Table 2. Summary of experimental (exp.) and simulated (sim.) variables, where $\psi$ is the mean bed slope, $\varTheta =R_{{b}}\tan \psi /[d_{{50}}(\rho _{{p}}/\rho _1 - 1)]$ is the Shields number and $q_b^\ast =q_{{b}}/[(\rho _{{p}}/\rho _1 - 1)gd_{{50}}^{{3}}]^{{1/2}}$ is the Einstein bedload number.

Comparing the numbers in table 2, the numerical and experimental Shields numbers are very similar. We note that in our simulations the horizontal boundary conditions were periodic, that is, the domain was not bound by lateral walls as in the physical experiments. Therefore, for the simulations, the boundary shear stress in the Shields number was computed based on the average water depth ($R_{{b}}=\bar {h}$) from $t=30$ s to $t=75$ s; whereas, for the experiments, the boundary shear stress considered the hydraulic radius as a correction for the effect of the lateral wall roughness. The simulated sediment transport rates were lower than the experimental mean values by a factor of approximately $2$. Nevertheless, the increase in sediment transport rate from case E1 to E4 was similar in the simulations and the experiments (roughly an increment by a factor of $3$).

Figure 7 shows the evolution of the simulated sediment transport rates with time for E1, together with a representative time window of similar duration in the experiments. Compared to the simulation data, the experimental transport rates show larger fluctuations. These variations are likely due to sampling, which captured all particles ejected from the downstream end of the channel in the experiment, as opposed to the averaged sediment transport rates over the entire simulation domain. It is, therefore, remarkable that the simulated and experimental transport rates match well in the last $15$ s of the time series. Some reasons for the underestimation of the average sediment transport rates in the simulations can be related to the shape of the particles (natural gravel in the experiments, spheres in the simulations (Deal et al. Reference Deal, Venditti, Benavides, Bradley, Zhang, Kamrin and Perron2023, cf.)), to experimental uncertainties as, for instance, additional moment added by the sediment feeding system, and to the strong nonlinear dependence of sediment transport rates on boundary conditions and the wall shear stresses. These notwithstanding, the comparison between numerical and experimental hydraulic variables and sediment transport rates must be rated as satisfactory, and, given the good correspondence between the simulations and the experiments, we conclude that our simulation approach is adequate to investigate UMAs with high fidelity.

Figure 7. Bedload transport rate $q_{{b}}(t)$ as measured in the experiment E1 from Pascal et al. (Reference Pascal, Ancey and Bohorquez2021) and in the numerical simulations. Time-averaged values were computed for a time range of $2804$ s and $45$ s for the experiment and simulation, respectively.

5. Conclusion

We have performed particle-resolved numerical simulations of supercritical turbulent flows (${Fr}>1$) over an erodible granular bed of spherical particles. Our goal was to numerically reproduce the experimental work of Pascal et al. (Reference Pascal, Ancey and Bohorquez2021), in which upstream-migrating antidunes developed on a bed of natural gravel sheared by a shallow-water flow (water depth $\approx 3$ times the sediment grain size). Supercritical shallow flows over erodible beds are extremely unstable due to the strong feedback between the dynamics of the bottom and the wavy free surface. Our results demonstrate that one can accurately reproduce the bedform dynamics at different supercritical flow conditions by using a massively parallel simulation approach relying on fully coupled fluid–particle–gas interactions. In our simulations, no suspended load was present and the bed instability was naturally excited by the simulated flow. Hence, the range of amplitudes, wavelengths and celerities of the self-generated upstream-migrating bedforms were results of the simulations and agree well with the experimental patterns reported by Pascal et al. (Reference Pascal, Ancey and Bohorquez2021). Particularly, the averaged bed Shields stress from the simulations, as a key quantity to predict sediment transport rates, matched very well with the experimental values.

The vast amount of data generated by simulations in a non-intrusive manner such as those presented in this work can be of great use to supplement experimental measurements under challenging supercritical flow conditions. The detailed data of particle and fluid motion available through our simulation approach opens a multitude of possibilities. The present study, therefore, encourages further simulation campaigns to understand antidune mechanics and the physical controls on bedform initiation and morphodynamics in supercritical flows.

Supplementary material

Animations of the simulations are available online as supplementary movies available at https://doi.org/10.1017/jfm.2023.262. The source code of the software and the simulation results are available at https://doi.org/10.5281/zenodo.7794080.

Acknowledgements

Francisco Núñez-González is María Zambrano Fellow in the Department of Civil and Environmental Engineering, Universitat Politécnica de Catalunya – BarcelonaTech. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for providing computing time on the GCS Supercomputer SuperMUC at the Leibniz Supercomputing Centre (www.lrz.de).

Funding

The authors thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for funding the projects 408062554, 434946896, 433735254 and 428445330. This work was supported by the SCALABLE project. This project has received funding from the European High-Performance Computing Joint Undertaking (J.U.) under Grant Agreement No. 956000. The J.U. receives support from the European Union's Horizon 2020 research and innovation programme and France, Germany and the Czech Republic.

Declaration of interests

The authors report no conflict of interest.

Author contributions

C.S.: conceptualization, methodology, software, formal analysis, validation, investigation, data curation, writing the original draft, review and editing, visualization, supervision, project administration. C.R.: conceptualization, methodology, software, formal analysis, validation, investigation, data curation, writing the original draft, visualization, supervision, project administration. S.K.: software, validation, investigation, writing the original draft. J.P.: software, validation, writing the original draft. F.N.-G.: conceptualization, writing the original draft, review and editing. H.K.: resources, funding acquisition. U.R.: resources, funding acquisition, writing the original draft. B.V.: conceptualization, writing the original draft, review and editing.

References

Aidun, C.K., Lu, Y. & Ding, E.-J. 1998 Direct analysis of particulate suspensions with inertia using the discrete Boltzmann equation. J. Fluid Mech. 373, 287311.CrossRefGoogle Scholar
Balmforth, N.J. & Vakil, A. 2012 Cyclic steps and roll waves in shallow water flow over an erodible bed. J. Fluid Mech. 695, 3562.CrossRefGoogle Scholar
Bauer, M., et al. 2021 waLBerla: a block-structured high-performance framework for multiphysics simulations. Comput. Maths Applics. 81, 478501.CrossRefGoogle Scholar
Bohorquez, P., Cañada-Pereira, P., Jimenez-Ruiz, P.J. & del Moral-Erencia, J.D. 2019 The fascination of a shallow-water theory for the formation of megaflood-scale dunes and antidunes. Earth-Sci. Rev. 193, 91108.CrossRefGoogle Scholar
Coleman, S.E. & Nikora, V.I. 2011 Fluvial dunes: initiation, characterization, flow structure. Earth Surf. Process. Landf. 36, 3957.CrossRefGoogle Scholar
Colombini, M. & Stocchino, A. 2012 Three-dimensional river bed forms. J. Fluid Mech. 695, 6380.CrossRefGoogle Scholar
Deal, E., Venditti, J.G., Benavides, S.J., Bradley, R., Zhang, Q., Kamrin, K. & Perron, J.T. 2023 Grain shape effects in bed load sediment transport. Nature 613 (7943), 298302.CrossRefGoogle ScholarPubMed
Geier, M., Schönherr, M., Pasquali, A. & Krafczyk, M. 2015 The cumulant lattice Boltzmann equation in three dimensions: theory and validation. Comput. Maths Applics. 70 (4), 507547.CrossRefGoogle Scholar
Kennedy, J.F. 1963 The mechanics of dunes and antidunes in erodible-bed channels. J. Fluid Mech. 16 (4), 521544.CrossRefGoogle Scholar
Kidanemariam, A.G. & Uhlmann, M. 2014 Direct numerical simulation of pattern formation in subaqueous sediment. J. Fluid Mech. 750, R2.CrossRefGoogle Scholar
Krüger, T., Kusumaatmaja, H., Kuzmin, A., Shardt, O., Silva, G. & Viggen, E.M. 2017 The Lattice Boltzmann Method: Principles and Practice. Springer.CrossRefGoogle Scholar
Núñez-González, F. & Martín-Vide, J.P. 2011 Analysis of antidune migration direction. J. Geophys. Res. 116 (F2), F02004.CrossRefGoogle Scholar
Olsen, N.R.B. 2017 Numerical modelling of downstream migrating antidunes. Earth Surf. Process. Landf. 42 (14), 23932401.CrossRefGoogle Scholar
Olsen, N.R.B. 2022 Explaining the formation of sedimentary structures under antidunes using a 2D width-averaged numerical model. Norwegian J. Geol. 102, 202210.Google Scholar
Parker, G. & Izumi, N. 2000 Purely erosional cyclic and solitary steps created by flow over a cohesive bed. J. Fluid Mech. 419, 203238.CrossRefGoogle Scholar
Pascal, I., Ancey, C. & Bohorquez, P. 2021 The variability of antidune morphodynamics on steep slopes. Earth Surf. Process. Landf. 46 (9), 17501765.CrossRefGoogle Scholar
Rettinger, C. & Rüde, U. 2022 An efficient four-way coupled lattice Boltzmann – discrete element method for fully resolved simulations of particle-laden flows. J. Comput. Phys. 453, 110942.CrossRefGoogle Scholar
Schwarzmeier, C., Holzer, M., Mitchell, T., Lehmann, M., Häusl, F. & Rüde, U. 2023 Comparison of free-surface and conservative Allen–Cahn phase-field lattice Boltzmann method. J. Comput. Phys. 473, 111753.CrossRefGoogle Scholar
Slootman, A., Ventra, D., Cartigny, M., Normandeau, A. & Hubbard, S. 2021 Supercritical-flow processes and depositional products: introduction to thematic issue. Sedimentology 68 (4), 12891296.CrossRefGoogle Scholar
Vellinga, A.J., Cartigny, M.J.B., Eggenhuisen, J.T. & Hansen, E.W.M. 2018 Morphodynamics and depositional signature of low-aggradation cyclic steps: new insights from a depth-resolved numerical model. Sedimentology 65 (2), 540560.CrossRefGoogle Scholar
Wen, B., Zhang, C., Tu, Y., Wang, C. & Fang, H. 2014 Galilean invariant fluid–solid interfacial dynamics in lattice Boltzmann simulations. J. Comput. Phys. 266, 161170.CrossRefGoogle Scholar
Figure 0

Figure 1. (a) Schematic representation of the simulations, coupling the liquid (l) with a particle (p) and gas (g) phase, and (b) the simulation set-up in its initial condition. Panel (b) is a zoom into the computational domain covering only 16 % of its streamwise extent.

Figure 1

Table 1. Summary of the main simulation conditions. The labels are chosen as in the reference experiments from Pascal et al. (2021).

Figure 2

Figure 2. Simulated bedload transport rate $q_{{b}}(t)$ with different computational grid resolutions using a shortened simulation domain. The time averages were computed only in the time ranges conforming to the lengths of the plotted lines to ensure that the fluid flow was adequately developed.

Figure 3

Figure 3. Visualization of the simulated velocity in streamwise direction $u_{x,{l}}(\boldsymbol {x},t)$. The undulations of the sediment bed and free surface of the liquid are in phase, conforming to the definition of antidunes. The figure is a zoom into the computational domain covering only $25$ % of its streamwise extent.

Figure 4

Figure 4. Sediment bed elevation $h_{{b}}(x,t)$ for E1 (a,b) and E4 (c,d): (a,c) experimental data from Pascal et al. (2021), and (b,d) simulation data from the numerical simulations. We consider the system fully developed after $t=30$ s, as illustrated by the dashed black line.

Figure 5

Figure 5. PSD of wavelength $\lambda$ and period $T$ for (a) E1 and (b) E4 from Pascal et al. (2021). The orange line is the outermost hull of the simulation results’ PSD.

Figure 6

Figure 6. PSD contour of celerity $c$ and wavelength $\lambda$ for (a) E1 and (b) E4 from Pascal et al. (2021).

Figure 7

Table 2. Summary of experimental (exp.) and simulated (sim.) variables, where $\psi$ is the mean bed slope, $\varTheta =R_{{b}}\tan \psi /[d_{{50}}(\rho _{{p}}/\rho _1 - 1)]$ is the Shields number and $q_b^\ast =q_{{b}}/[(\rho _{{p}}/\rho _1 - 1)gd_{{50}}^{{3}}]^{{1/2}}$ is the Einstein bedload number.

Figure 8

Figure 7. Bedload transport rate $q_{{b}}(t)$ as measured in the experiment E1 from Pascal et al. (2021) and in the numerical simulations. Time-averaged values were computed for a time range of $2804$ s and $45$ s for the experiment and simulation, respectively.

Schwarzmeie et al. Supplementary Movie 1

Animation E1 initial phase - Simulation of the initial phase of the test case E1. The animations show 30 s of simulated time, starting from an initially flat bed and free surface. The animations play 2 times slower to aid observation. The animation at the top shows a side view, where particles are colored by their diameter, and the flow field is colored by its velocity in streamwise direction. The animation in the center shows the sediment bed from the top view, where particles are colored by their elevation. Note that the periodic simulation domain is repeated three times in spanwise direction to provide a better visual impression of the bedforms. The animation at the bottom shows the free-surface elevation from the top view. To keep computational costs for post-processing manageable, we only stored a two-dimensional slice of the flow field. The free-surface elevation animation is extruded from this slice to conform with the simulated domain width.

Download Schwarzmeie et al. Supplementary Movie 1(Video)
Video 44.2 MB

Schwarzmeie et al. Supplementary Movie 2

Animation E1 developed system - Simulation of the developed phase of the test case E1. The animations show 45 s of simulated time, starting from a developed system. The animations play 2 times slower to aid observation. The animation at the top shows a side view, where particles are colored by their diameter, and the flow field is colored by its velocity in streamwise direction. The animation in the center shows the sediment bed from the top view, where particles are colored by their elevation. Note that the periodic simulation domain is repeated three times in spanwise direction to provide a better visual impression of the bedforms. The animation at the bottom shows the free-surface elevation from the top view. To keep computational costs for post-processing manageable, we only stored a two-dimensional slice of the flow field. The free-surface elevation animation is extruded from this slice to conform with the simulated domain width.

Download Schwarzmeie et al. Supplementary Movie 2(Video)
Video 45.7 MB

Schwarzmeie et al. Supplementary Movie 3

Animation E4 initial phase - Simulation of the initial phase of the test case E4. The animations show 30 s of simulated time, starting from an initially flat bed and free surface. The animations play 2 times slower to aid observation. The animation at the top shows a side view, where particles are colored by their diameter, and the flow field is colored by its velocity in streamwise direction. The animation in the center shows the sediment bed from the top view, where particles are colored by their elevation. Note that the periodic simulation domain is repeated three times in spanwise direction to provide a better visual impression of the bedforms. The animation at the bottom shows the free-surface elevation from the top view. To keep computational costs for post-processing manageable, we only stored a two-dimensional slice of the flow field. The free-surface elevation animation is extruded from this slice to conform with the simulated domain width.

Download Schwarzmeie et al. Supplementary Movie 3(Video)
Video 50.4 MB

Schwarzmeie et al. Supplementary Movie 4

Animation E4 developed system - Simulation of the developed phase of the test case E4. The animations show 45 s of simulated time, starting from a developed system. The animations play 2 times slower to aid observation. The animation at the top shows a side view, where particles are colored by their diameter, and the flow field is colored by its velocity in streamwise direction. The animation in the center shows the sediment bed from the top view, where particles are colored by their elevation. Note that the periodic simulation domain is repeated three times in spanwise direction to provide a better visual impression of the bedforms. The animation at the bottom shows the free-surface elevation from the top view. To keep computational costs for post-processing manageable, we only stored a two-dimensional slice of the flow field. The free-surface elevation animation is extruded from this slice to conform with the simulated domain width.

Download Schwarzmeie et al. Supplementary Movie 4(Video)
Video 49.8 MB