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
\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.
Open access
