2.4.1. Data compilation for turbulent bedload transport of spherical grains

In order to isolate channel geometry effects on bedload flux, we compile existing experimental datasets for spherical grains, assuming that grain-shape effects are negligible for them. However, surprisingly, even though turbulent bedload transport has been studied for more than a century (Ancey Reference Ancey2020), flux measurements for spherical grains with a nearly uniform size distribution seem to be exceedingly scarce. (By contrast, such kinds of measurements are the rule rather than exception when studying bedload transport driven by laminar flow; e.g. Ouriemi et al. Reference Ouriemi, Aussillous, Medale, Peysson and Guazzelli2007; Aussillous et al. Reference Aussillous, Chauchat, Pailha, Médale and Guazzelli2013; Houssais et al. Reference Houssais, Ortiz, Durian and Jerolmack2015; Cúñez et al. Reference Cúñez, Franklin, Houssais, Arratia and Jerolmack2022.) After an extensive literature search, we have managed to find only a handful of datasets (Frey et al. Reference Frey, Dufresne, Böhm, Jodeau and Ancey2006; Frey Reference Frey2014; Ni & Capart Reference Ni and Capart2018; Armanini & Cavedon Reference Armanini and Cavedon2019; Carrillo et al. Reference Carrillo, Petrie, Timbe, Pacheco, Astudillo, Padilla and Cisneros2022; Deal et al. Reference Deal, Venditti, Benavides, Bradley, Zhang, Kamrin and Perron2023). Four of these were subsequently disregarded because of either partially crystallised beds and very large scatter (Frey et al. Reference Frey, Dufresne, Böhm, Jodeau and Ancey2006; Frey Reference Frey2014), or the presence of substantial suspended load in addition to bedload (Armanini & Cavedon Reference Armanini and Cavedon2019), or self-inconsistent data (Carrillo et al. Reference Carrillo, Petrie, Timbe, Pacheco, Astudillo, Padilla and Cisneros2022) (almost no bedload flux change with important variation of flow strength: increasing bed slope at constant flow depth). The only two remaining experimental datasets for spherical grains are those by Ni & Capart (Reference Ni and Capart2018) and Deal et al. (Reference Deal, Venditti, Benavides, Bradley, Zhang, Kamrin and Perron2023) discussed in the Introduction. We therefore supplement the data compilation with additional datasets from numerical simulations that couple the DEM for the particle phase with a CFD method for the fluid phase that resolves the sub-grain scale: two datasets based on DNS (Kidanemariam & Uhlmann Reference Kidanemariam and Uhlmann2017; Jain et al. Reference Jain, Tschisgale and Fröhlich2021) for infinitely wide channels ( $b/h=\infty$ , i.e. periodic boundaries), and two datasets based on grain-resolved large-eddy simulations (LES) (Zhang et al. Reference Zhang, Deal, Perron, Venditti, Benavides, Rushlow and Kamrin2022): one for $b/h=\infty$ , and one for a narrow-deep-channel configuration ( $b/h\approx 0.1$ ). Only such high-end methods, because they do not coarse-grain the particle phase and do not rely on empiricism for fluid-particle interactions, can be considered as similarly reliable as controlled experiments. The above spherical-grain datasets are denoted as D23EXPs, NC18EXP, Z22LESn, Z22LESw, J21DNS and KU17DNS, and summarised in table 1 together with other datasets used later. Note that to Z22LESn and Z25LESn (used later), we are unable to apply our sidewall correction method, since the bulk flow velocity $U$ , or data from which it can be inferred, was not reported. In addition, wall-modelled LES such as these are known to underestimate the friction factors of channels with smooth walls (Nikitin et al. Reference Nikitin, Nicoud, Wasistho, Squires and Spalart2000; Yang et al. Reference Yang, Park and Moin2017), rendering (2.10) unusable. We therefore use the empirical sidewall correction that Zhang et al. (Reference Zhang, Deal, Perron, Venditti, Benavides, Rushlow and Kamrin2025) applied to Z22LESn and Z25LESn, $\tau _b=3\tau _R$ , corresponding to $f_w/f_b\approx 0.29$ , which is justified because it leads to an overlap with the sidewall-free simulations Z22LESw and Z25LESw, respectively, for otherwise almost the same conditions (table 1).

Table 1.

Summary of bedload flux data shown in this study. The top, middle and bottom entries correspond to experiments, grain-resolved CFD-DEM simulations, and grain-unresolved CFD-DEM simulations, respectively. The data of the datasets BLRANS1 and BLRANS2 with $s=2.65$ and $\alpha =0$ are from Pähtz & Durán (Reference Pähtz, Liu, Xia, Hu, He and Tholen2020). The datasets Z25LESn and Z25LESw do not include data from simulations based on the ‘artificial- shrinkage method’ by Zhang et al. (Reference Zhang, Deal, Perron, Venditti, Benavides, Rushlow and Kamrin2025), since this method has been falsified (Chen et al. Reference Chen, Durán and Pähtz2026). The asterisk $\ast$ indicates that the friction factor ratio $f_w/f_b$ was not calculated but determined empirically. Note that the listed parameter values of $C_D$ , ${\textit{Ga}}$ , $\varTheta$ and $h^\ast$ are based on $d=c_g$ .

Table 1. Long description

Table with three main sections: experiments, grain-resolved CFD-DEM simulations, and grain-unresolved CFD-DEM simulations. Each section lists datasets with corresponding parameters. The top section includes datasets BLRANS1 and BLRANS2. The middle and bottom sections list datasets Z25LESn and Z25LESw. The table also notes the exclusion of data based on the artificial-shrinkage method by Zhang et al. (2025) due to its falsification by Chen et al. (2026). An asterisk indicates empirically determined friction factor ratios. The table provides parameter values for various datasets, with specific calculations and assumptions detailed for certain datasets.

For the spherical-grain datasets discussed in the previous paragraph, we determine the bed surface elevation using the method described in § 2.2, provided that the required data are available from the studies. This is the case for the dataset NC18EXP (figure 2), but not for D23EXPs, Z22LESn, Z22LESw, J21DNS and KU17DNS. For the latter, we use the values reported in the respective studies. This approximation is justified for D23EXPs, Z22LESn, J21DNS and KU17DNS, since $h^\ast =h/d\gg 1$ (table 1). In the case of Z22LESw, the authors of the respective study (Zhang et al. Reference Zhang, Deal, Perron, Venditti, Benavides, Rushlow and Kamrin2022) stated that they employed the method by Pähtz & Durán (Reference Pähtz and Durán2018a ) to calculate the bed surface elevation, which is the same as our method described in § 2.2. For more details and the calculation of the bed shear stress uncertainty, see Appendix B.

Figure 2.

Determination of bed surface elevation for experimental dataset NC18EXP (Ni & Capart Reference Ni and Capart2018). The vertical profiles of $\dot \gamma \varSigma \phi$ exhibit several local maxima due to layering. We choose the top-most local maximum as the bed surface elevation, since the vertical profile of the fluid shear stress $\sigma ^f_{zx}(z)$ exhibits a focal point there, $\sigma ^f_{zx}(0)\approx \tau _t$ (§ 2.4.3), consistent with the theoretical expectation (Pähtz & Durán Reference Pähtz and Durán2018b ).

2.4.2. Test of channel geometry correction methods against data compilation

In this subsubsection, we test our existing channel geometry correction methods against the spherical-grain data compilation discussed above. Once the bed shear stress $\tau _b$ is determined using a given channel geometry correction method, we are able to fully characterise the hydrodynamic bedload transport conditions by dimensionless similarity parameters. We choose the particle–fluid density ratio $s\equiv \rho _{\!p}/\rho _{\!f}$ , Galileo number ${\textit{Ga}}\equiv d\sqrt {s\tilde g_zd}/\nu$ , Shields number $\varTheta \equiv \tau _b/(\rho _{\!p}\tilde g_zd)$ , dimensionless flow depth $h^\ast \equiv h/d$ , and bed slope angle $\alpha$ (table 1). All spherical-grain datasets (and also those for non-spherical grains considered later) are well within the range of what is typically categorised as bedload transport ( $s\lesssim 10$ ; Pähtz & Durán Reference Pähtz, Liu, Xia, Hu, He and Tholen2020), and most of them also satisfy $s^{1/2}Ga\gtrsim 70$ , here termed ‘rough’ transport conditions. In rough turbulent bedload, particle inertia dominates viscous fluid–particle interactions, and the non-dimensionalised bedload flux $Q_\ast \equiv Q/(\rho _{\!p}d\sqrt {s\tilde g_zd})$ is therefore expected to become essentially independent of $s$ and ${\textit{Ga}}$ for sufficiently small $\alpha$ (Ancey Reference Ancey2020; Pähtz & Durán Reference Pähtz, Liu, Xia, Hu, He and Tholen2020; Zhang et al. Reference Zhang, Deal, Perron, Venditti, Benavides, Rushlow and Kamrin2022). We use this expectation as a tool to evaluate the different sidewall corrections outlined above (figure 3). It can be seen that the data for rough turbulent bedload transport of spherical grains indeed collapse, within measurement uncertainty, on a universal $Q_\ast (\varTheta )$ behaviour for our method (figure 3 d), but do not collapse at all when employing a classical method or left uncorrected (figures 3 a–c), or when using our sidewall correction but a different bed surface definition (NC18EXP $\ast$ in figure 3 d). Note that the spherical-grain dataset Z22LESn is not included in figure 3(c), since the data for the bulk flow velocity $U$ , required to calculate the Einstein–Johnson sidewall correction, are unavailable. Furthermore, Z22LESn in figure 3(d) is based not on our sidewall correction (since it also requires $U$ as an input), but on the correction that Zhang et al. (Reference Zhang, Deal, Perron, Venditti, Benavides, Rushlow and Kamrin2025) applied to this dataset, $\tau _b=3\tau _R$ .

Figure 3.

Channel geometry corrections. Non-dimensionalised bedload flux $Q_\ast$ versus Shields numbers for spherical grains: (a) uncorrected for sidewall effects ( $\varTheta _0$ ); (b) hydraulic-radius correction (Guo Reference Guo2015) ( $\varTheta _R$ ); (c) Einstein–Johnson correction (Einstein Reference Einstein1942; Johnson Reference Johnson1942; Vanoni & Brooks Reference Vanoni and Brooks1957; Guo Reference Guo2017) ( $\varTheta _E$ ); (d) our universal correction ( $\varTheta$ ). Symbols correspond to data from experiments and grain-resolved simulations (table 1). An error bar indicates the standard error and/or uncertainty range. In the absence of error bars, uncertainties are smaller than the symbol size. Except for KU17DNS ( $s^{1/2}\,Ga\approx 44.86$ , smooth), all data are in the rough transport regime ( $s^{1/2}\,Ga\gt 70$ ). The values of $h$ used to determine $\tau _b$ in Z22LESw and NC18EXP use a granular-physics-based bed surface definition (Pähtz & Durán Reference Pähtz and Durán2018b ), (2.7), and in NC18EXP $\ast$ are based on the original definition reported by Ni & Capart (Reference Ni and Capart2018). For all other datasets, reported values of $h$ are used as it has a negligible effect on the Shields numbers due to $h\gg d$ . Case Z22LESn in (d) is based not on our sidewall correction, but on the correction that Zhang et al. (Reference Zhang, Deal, Perron, Venditti, Benavides, Rushlow and Kamrin2025) applied to this dataset, $\varTheta =3\varTheta _R$ (see text).

Figure 3. Long description

The image contains four scatter plots comparing non-dimensionalized bedload flux versus Shields numbers for spherical grains under different correction methods. Panel A shows uncorrected data for sidewall effects. Panel B presents data corrected using the hydraulic-radius correction method. Panel C displays data corrected with the Einstein-Johnson correction. Panel D illustrates data corrected using a universal correction method. Each panel includes data points from experiments and grain-resolved simulations, represented by different symbols and colors. Error bars indicate standard error and uncertainty ranges, with some data points lacking visible error bars due to their small size. The data spans both rough and smooth transport regimes, with specific datasets identified by unique symbols. The graphs aim to clarify the impact of channel geometry corrections on bedload transport measurements.

2.4.3. Further independent tests of our channel geometry correction method

Figure 4.

Channel geometry corrections. (a) Plot of $\varTheta$ versus $2.41\varTheta _R$ (using $d=c_g$ for non-spherical-grain data, as explained in § 3.1.1), based on $\tau _b=2.41\tau _R$ , the empirical correction in Deal et al. (Reference Deal, Venditti, Benavides, Bradley, Zhang, Kamrin and Perron2023). (b) Log-profile shear velocity $u_\tau \approx \sqrt {\tau _b/\rho _{\!f}}$ predicted by our sidewall correction versus measured one for experimental data from particle-free flows in open channels consisting of artificially created hydraulically rough beds ( $r=d_0$ or $r=k$ ) and hydraulically smooth sidewalls (Song Reference Song1994; Song & Chiew Reference Song and Chiew2001; Auel et al. Reference Auel, Albayrak and Boes2014) (data as summarised in the tables of Guo Reference Guo2015). For non-uniform flows, the bed slope $\tan \alpha$ is corrected using (1b) of Guo (Reference Guo2015).

Figure 5.

Independent evidence supporting our bed surface elevation definition. Vertical profiles (solid lines) of the non-dimensionalised fluid shear stress $\varTheta ^f_{zx}\equiv \sigma ^f_{zx}/(\rho _{\!p}\tilde g_zd)$ for the dataset NC18EXP (Ni & Capart Reference Ni and Capart2018). The data exhibit a focal point at the bed surface elevation $z=0$ defined by (2.7): $\varTheta ^f_{zx}(0)\approx 0.04$ , a property required for Bagnoldian-type bedload models (Pähtz & Durán Reference Pähtz and Durán2018b ). The fluid shear stress $\sigma ^f_{zx}$ is calculated from local quantities at the channel centre as described in Chauchat (Reference Chauchat2018), assuming that the mixture viscosity obeys the closure measured by Boyer et al. (Reference Boyer, Guazzelli and Pouliquen2011) for viscous suspensions. The dashed lines correspond to the bed surface elevations determined by Ni & Capart (Reference Ni and Capart2018).

Figure 5. Long description

A line graph displays vertical profiles of the non-dimensionalized fluid shear stress for the dataset NC18EXP. The x-axis represents the non-dimensionalized fluid shear stress, labeled as Theta f subscript xz, ranging from 0 to 0.06. The y-axis represents the normalized height, labeled as z divided by d, ranging from -2.5 to 2.5. The graph includes multiple colored lines, each representing different values of Theta, ranging from 0.119 to 0.207. The solid lines indicate the fluid shear stress calculated from local quantities at the channel center, while the dashed lines correspond to the bed surface elevations determined by Ni & Capart 2018. The data exhibit a focal point at the bed surface elevation defined by the equation z divided by d equals negative 1.5, a property required for Bagnoldian-type bedload models. The fluid shear stress is calculated assuming the mixture viscosity obeys the closure measured by Boyer et al. 2011 for viscous suspensions.

Figures 4(a) and 4(b) present further independent support for our sidewall correction. First, our correction predicts almost the same Shields numbers as those that Deal et al. (Reference Deal, Venditti, Benavides, Bradley, Zhang, Kamrin and Perron2023) empirically determined as $\tau _b=2.41\tau _R$ via fitting (2.1) to their flow velocity data, for both spherical and non-spherical grains (figure 4 a). Second, it is also consistent with experimental data from particle-free flows in open channels ( $b/h=2.94{-}7.79$ ) consisting of hydraulically rough beds and hydraulically smooth sidewalls (Song Reference Song1994; Song & Chiew Reference Song and Chiew2001; Auel et al. Reference Auel, Albayrak and Boes2014) (figure 4 b).

Furthermore, the strong quantitative difference between the spherical-grain data by Ni & Capart (Reference Ni and Capart2018) and Deal et al. (Reference Deal, Venditti, Benavides, Bradley, Zhang, Kamrin and Perron2023) discussed in the Introduction (figure 1) is found to be predominantly the consequence of the latter authors’ bed surface definition (NC18EXP $\ast$ in figure 3 d). In fact, figure 5 shows that the bed surface elevations resulting from their definition are much below ours. It can also be seen that the vertical profiles of the non-dimensionalised local fluid shear stress $\varTheta ^f_{zx}$ for various flow strengths exhibit a focal point at the bed surface elevation ( $z=0$ ) determined through our method (2.7), with a focal value $\varTheta ^f_{zx}(0)\approx 0.04$ close to the transport threshold Shields number $\varTheta _t$ . As discussed in § 2.2, the occurrence of such a focal point is precisely one of the properties that we expect the bed surface to have.

3.1.1. Bedload flux

Using a numerical model coupling the RANS equations with the DEM, Pähtz & Durán (Reference Pähtz, Liu, Xia, Hu, He and Tholen2020) simulated non-suspended sediment transport of spherical grains, and derived a general transport rate expression holding for a large range of conditions, including turbulent bedload with $s^{1/2}\,Ga\gtrsim 80$ and windblown sand:

(3.1) \begin{equation} Q_\ast =2\kappa ^{-1}\sqrt {\varTheta _t}\,M_\ast (1+c_MM_\ast ), \end{equation}

with $M_\ast \equiv M/(\rho _{\!p}d)$ . Based on the definition $M\equiv \rho _{\!p}h\langle \phi \rangle _z$ of the transport load in § 2.1, $M_\ast$ can be interpreted as the number $N$ of transported grains above a bed surface area of size $\sim d^2$ . That is, $M_\ast$ is a measure for the particle activity and therefore proportional to the excess Shields number $\varTheta -\varTheta _t$ (Ancey Reference Ancey2020); the proportionality factor is derived shortly. Hence for bedload transport, where $1+c_MM_\ast \approx 2\sqrt {c_MM_\ast }$ (arithmetic mean $\approx$ geometric mean) is a good approximation as typically $c_MM_\ast \sim O(1)$ , (3.1) predicts a bedload flux scaling similar to the classical formula by Meyer-Peter & Müller (Reference Meyer-Peter and Müller1948).

In (3.1), $2\kappa ^{-1}\sqrt {\varTheta _t}$ represents the dimensionless average particle velocity at $\varTheta _t$ , deduced from assuming that $\varTheta _t$ corresponds to the smallest Shields number at which a quasi-continuous hopping-rebound motion (saltation) of a test grain along a stationary bed in a grain-motion-undisturbed logarithmic fluid velocity profile can be sustained, a so-called ‘rebound threshold’ (Pähtz et al. Reference Pähtz and Durán2020). Furthermore, the term $c_MM_\ast$ , with $c_M=1.7$ a constant parameter determined from the simulations, encodes the rate of collisions between transported grains ( $\sim M_\ast ^2$ ) relative to that of grain-bed collisions ( $\sim M_\ast ^1$ ). In fact, (3.1) expresses the balance between global energy production due to drag (left-hand side) and energy dissipation due to grain contacts (right-hand side). In order to apply it to non-spherical grains, it is important to generalise $d$ in a manner that retains $M_\ast \sim N$ . Since grains transported at elevations close to the bed surface, where most transport takes place, tend to align their largest projected area $A^{\textit{max}}_{\!p}$ parallel to the bed due to torque (Jain et al. Reference Jain, Tschisgale and Fröhlich2020; Zhang et al. Reference Zhang, Deal, Perron, Venditti, Benavides, Rushlow and Kamrin2025), implying $M\sim N\rho _{\!p}V_{\!p}/A^{\textit{max}}_{\!p}$ , this requirement begets $d\equiv ({3}/{2})V_{\!p}/A^{\textit{max}}_{\!p}$ . From approximating grains as ellipsoids, we estimate $V_{\!p}=\pi a_gb_gc_g/6$ and $A^{\textit{max}}_{\!p}=\pi a_gb_g/4$ , leading to $d=c_g$ (already used in figure 4 a).

In steady, homogeneous sediment transport, the average vertical force acting on transported grains is dominated by gravity, buoyancy and potentially fluid lift, since drag forces during the grains’ upward motion tend to cancel those during their downward motion (Pähtz & Durán Reference Pähtz and Durán2018a ). When fluid lift is neglected, as in the simulations by Pähtz & Durán (Reference Pähtz and Durán2018b , Reference Pähtz, Liu, Xia, Hu, He and Tholen2020), this implies that the particle pressure at the bed surface $-\sigma ^p_{zz}(0)$ is approximately equal to the submerged particle weight $M\tilde g_z$ . The bed friction coefficient $\mu _b\equiv -\sigma ^p_{zx}(0)/\sigma ^p_{zz}(0)$ then links $M\tilde g_z$ to the surface particle shear stress $\sigma ^p_{zx}(0)$ . When further separating $\sigma ^p_{zx}(0)$ into gravity and drag contributions, and assuming that particle drag on the fluid reduces the surface fluid shear stress $\sigma ^f_{zx}(0)$ to a value that is just sufficient to sustain transport, one obtains (after some rearrangements) Bagnold’s famous expression (Bagnold Reference Bagnold1956; Pähtz & Durán Reference Pähtz and Durán2018b )

(3.2) \begin{equation} M_\ast =\frac {\varTheta -\varTheta _t}{\mu _b-\tan \alpha }, \end{equation}

which is a key ingredient of numerous bedload transport and windblown sand models in the literature, including the one by Pähtz & Durán (Reference Pähtz, Liu, Xia, Hu, He and Tholen2020). Historically, it has been assumed that $\mu _b$ is equal to $\mu _s$ , the tangent of the static angle of repose of the bulk bed material (Bagnold Reference Bagnold1956). However, RANS-DEM simulations of turbulent bedload transport of spherical grains suggest that although constant, $\mu _b$ is actually substantially larger than $\mu _s$ ( $\mu _b\approx 0.7$ versus $\mu _s\approx 0.4$ ) (Pähtz & Durán Reference Pähtz and Durán2018b ). This observation can be understood as a bed surface strengthening effect (Clark et al. Reference Clark, Shattuck, Ouellette and O’Hern2017; Allen & Kudrolli Reference Allen and Kudrolli2018; Cúñez et al. Reference Cúñez, Franklin, Houssais, Arratia and Jerolmack2022). Bedload transport increases the bed surface’s ability to resist shear stress, since temporarily mobilised grains can settle again in more stable bed surface pockets (Clark et al. Reference Clark, Shattuck, Ouellette and O’Hern2017). Whenever this happens, the bed surface becomes stronger relative to the bulk, i.e. $r_b\equiv \mu _b/\mu _s$ increases. This process continues until a state is reached at which the pockets that emerge at the places where grains are mobilised are equivalent to the most stable pockets that transported grains can reach (Clark et al. Reference Clark, Shattuck, Ouellette and O’Hern2017). Then $r_b$ has acquired its equilibrium value, and will no longer increase. In viscous bedload transport, where all pockets are reachable as a transported grain’s kinetic energy is typically much smaller than the potential energy wells of the bed surface, $r_b\approx 3.4$ (Pähtz et al. Reference Poggiali, Mastrogiuseppe, Hayes, Seu, Birch, Lorenz, Grima and Hofgartner2021). However, in this study on turbulent bedload transport, where these two energy scales are of comparable size, $r_b$ attains a smaller value as grains are too fast to probe the entire phase space of bed surface pockets (Clark et al. Reference Clark, Shattuck, Ouellette and O’Hern2017). We use the value $r_b=1.8$ suggested by the simulations by Pähtz & Durán (Reference Pähtz and Durán2018b ), where $\mu _b\approx 0.7$ and $\mu _s\approx 0.4$ , to calculate $\mu _b=r_b\mu _s$ from the measured $\mu _s$ .

In addition to using $\mu _b=r_b\mu _s$ , we also generalise (3.2) to account for lift forces, defined as non-buoyancy fluid forces that act in the bed-normal direction. Neglecting particle–sidewall interactions in the particle momentum balance, a derivation analogous to that in Pähtz & Durán (Reference Pähtz, Liu, Xia, Hu, He and Tholen2020) yields (Appendix C)

(3.3) \begin{equation} M_\ast =(\varTheta -\varTheta _t)/\mu ^\dagger ,\quad \text{with} \quad \mu ^\dagger \equiv r_b\mu _s\left [\frac {1+\left ({s\tan \alpha }/{(s-1)}+\chi ^\ast \right )\overline {f_L}/\overline {f_D}}{1+r_b\mu _s\overline {f_L}/\overline {f_D}}\right ]-\tan \alpha , \end{equation}

where $\chi ^\ast \equiv \chi /(\rho _{\!p}\tilde g_z)$ , $f_D$ is the drag force per unit mass, $f_L$ is the lift force per unit mass divided by the fluid volume fraction $1-\phi$ , and the overbar, $\overline {{\cdot }}\equiv \langle \phi\, {\cdot }\rangle _z/\langle \phi \rangle _z$ , denotes the $\phi$ -weighted height average, i.e. the average over all transported grains. The expression for $\mu ^\dagger$ takes into account that Reynolds stress gradients do not contribute to the buoyancy force (Maurin et al. Reference Maurin, Chauchat and Frey2018; Zhu et al. Reference Zhu, Chen, Tholen, He and Pähtz2026), a recent finding that implies that there is no buoyancy force counteracting the gravitational and pressure gradient forces in the streamwise direction.

Equation (3.3) contains the single additional calibratable parameter $\overline {f_L}/\overline {f_D}$ , and reduces to (3.2) if fluid lift is neglected. In real-world scenarios, one may expect $\overline {f_L}/\overline {f_D}$ to be of order unity, as suggested by measurements of drag and lift forces acting on surface roughness elements, including grains in bed surface pockets (Chepil Reference Chepil1958, Reference Chepil1961), and by DNS-DEM of sheared mobile beds (Ji et al. Reference Ji, Munjiza, Avital, Ma and Williams2013). We use $\overline {f_L}/\overline {f_D}=1.0$ because it optimises the agreement between the final model and the entire data compilation.

3.1.2. Transport threshold

We reiterate that the prefactor $2\sqrt {\varTheta _t}/\kappa$ in (3.1) resulted from assuming that $\varTheta _t$ is the smallest Shields number that permits a quasi-continuous hopping-rebound motion (saltation) of a test grain along a stationary bed in a grain-motion-undisturbed logarithmic fluid velocity profile (Pähtz & Durán Reference Pähtz, Liu, Xia, Hu, He and Tholen2020). In such a scenario, lift forces acting on the moving test grain are typically much smaller than drag forces due to the formers’ rapid decrease with distance from the bed (Chepil Reference Chepil1961; Li et al. Reference Li, Balachandar, Lee and Bai2019). Then it can be shown that for conditions with $s^{1/4}\,Ga\gtrsim 200$ , $\varTheta _t$ in (3.1) obeys (Pähtz et al. Reference Poggiali, Mastrogiuseppe, Hayes, Seu, Birch, Lorenz, Grima and Hofgartner2021)

(3.4) \begin{equation} \varTheta _t=c_t\left (r_b\mu _s-\frac {s}{s-1}\tan \alpha -\chi ^\ast \right )/C_D, \end{equation}

where $c_t$ is a dimensionless parameter, and $C_D$ is an effective drag coefficient parametrising streamwise drag force acting on grains. However, in spite of the similarity of (3.4) with scaling laws classically associated with the incipient motion of a grain resting in a bed surface pocket, we emphasise again that $\varTheta _t$ is not an incipient motion threshold, but a so-called rebound threshold, as defined and discussed by Pähtz et al. (Reference Pähtz and Durán2020, Reference Poggiali, Mastrogiuseppe, Hayes, Seu, Birch, Lorenz, Grima and Hofgartner2021). In particular, incipient motion thresholds strongly depend on turbulent fluctuations and therefore exhibit an additional dependence on $h^\ast$ that is responsible for their often observed increase with bed slope $\tan \alpha$ (Lamb et al. Reference Lamb, Dietrich and Venditti2008). By contrast, the rebound threshold $\varTheta _t$ is insensitive to turbulent fluctuations (Pähtz et al. Reference Pähtz and Durán2020), as was also directly shown in figure S1 of Pähtz & Durán (Reference Pähtz and Durán2018a ) based on numerical simulation data by Maurin et al. (Reference Maurin, Chauchat, Chareyre and Frey2015), and therefore always decreases with $\tan \alpha$ .

The drag coefficient $C_D$ in (3.4) is smaller than that for settling grains, $C_{D\textit{Settle}}\equiv 4(s-1)gd/(3v_s^2)$ (Deal et al. Reference Deal, Venditti, Benavides, Bradley, Zhang, Kamrin and Perron2023), where $v_s$ is the terminal grain settling velocity, due to the smaller projected grain area exposed to the flow. Assuming again that transported grains tend to align their largest projected area $A^{\textit{max}}_{\!p}=\pi a_gb_g/4$ parallel to the bed, the effective projected area $A^{\textit{eff}}_{\!p}$ in the streamwise direction is estimated as the geometric mean of the remaining two projected areas $\pi a_gc_g/4$ and $\pi b_gc_g/4$ : $A^{\textit{eff}}_{\!p}=\pi \sqrt {a_gb_g}\,c_g/4$ , implying $A^{\textit{eff}}_{\!p}/A^{\textit{max}}_{\!p}=c_g/\sqrt {a_gb_g}\equiv S_{\!f}$ , which is the Corey shape factor. Hence we obtain $C_D=S_{\!f}C_{D\textit{Settle}}$ or

(3.5) \begin{equation} C_D=4S_{\!f}(s-1)gd/(3v_s^2). \end{equation}

In the case of spherical grains ( $S_{\!f}=1$ ), typical values for the quantities in (3.4) are $\mu _s\approx 0.46$ , $C_D\approx 0.43$ and $\varTheta _t\approx 0.05$ (Deal et al. Reference Deal, Venditti, Benavides, Bradley, Zhang, Kamrin and Perron2023), implying $c_t\approx 0.03$ . We use $c_t=0.032$ because it optimises the agreement between the final model and the entire data compilation. On the other hand, when $s^{1/4}\,Ga\lesssim 200$ , viscous-sublayer effects render the scaling of $\varTheta _t$ much more complicated (Pähtz et al. Reference Poggiali, Mastrogiuseppe, Hayes, Seu, Birch, Lorenz, Grima and Hofgartner2021). Therefore, we use (3.4) to determine $\varTheta _t$ , except for the few numerical conditions with $s^{1/4}\,Ga\lt 200$ (BLRANSx in table 1, where ‘x’ stands for an arbitrary character), for which we determine $\varTheta _t$ directly from the simulations, as described by Pähtz & Durán (Reference Pähtz, Liu, Xia, Hu, He and Tholen2020). Note that the datasets J21DNS and KU17DNS in table 1, which also satisfy $s^{1/4}\,Ga\lt 200$ , are excluded from this discussion as they violate the precondition $s^{1/2}\,Ga\gtrsim 80$ required for (3.1).

3.2.1. Experiments and grain-resolved CFD-DEM simulations

To evaluate our model, we require bedload flux datasets for shape-controlled grains, meaning that they must provide the values of $v_s$ and $S_{\!f}$ , required to calculate $C_D$ in (3.5), and the value of $\mu _s$ , which affects $Q_\ast$ both via the calculation of $M_\ast$ in (3.3) and via the calculation of $\varTheta _t$ in (3.4). The datasets by Deal et al. (Reference Deal, Venditti, Benavides, Bradley, Zhang, Kamrin and Perron2023), denoted as D23EXPx in table 1, and by Zhang et al. (Reference Zhang, Deal, Perron, Venditti, Benavides, Rushlow and Kamrin2022, Reference Zhang, Deal, Perron, Venditti, Benavides, Rushlow and Kamrin2025), denoted as Z22LESx and Z25LESx in table 1, satisfy this requirement. However, we excluded from Z25LESx simulation data based on the ‘artificial-shrinkage method’ by Zhang et al. (Reference Zhang, Deal, Perron, Venditti, Benavides, Rushlow and Kamrin2025), since this method has been falsified (Chen et al. Reference Chen, Durán and Pähtz2026). Furthermore, in order to use the spherical-grain dataset NC18EXP by Ni & Capart (Reference Ni and Capart2018), who did not report $\mu _s$ , we acquired the same spheres from the same company, and measured $\mu _s$ using the funnel method (Beakawi Al-Hashemi & Baghabra Al-Amoudi Reference Beakawi Al-Hashemi and Baghabra Al-Amoudi2018) in the manner described by Deal et al. (Reference Deal, Venditti, Benavides, Bradley, Zhang, Kamrin and Perron2023): by slowly pouring them onto an elevated disk bounded by a $2d_o$ high rim (Appendix D). Moreover, we scoured the literature and found two more usable experimental datasets for non-spherical grains: the bedload flux measurements for cylinders (R22EXPc) and lenses (R22EXPl) by Rebai et al. (Reference Rebai, Berzi, Ballio and Matousek2022). For both grain shapes, $v_s$ (Matoušek & Zrostlík Reference Matoušek and Zrostlík2020) and the grain dimensions (required for $S_{\!f}$ ) have been reported, while $\mu _s$ was measured using the funnel method in a box filled with water, which was sufficiently large to rule out potential box–sidewall friction effects (private correspondence with the authors). We confirmed with DEM simulations that different procedures of conducting the funnel method have only a small effect on $\mu _s$ ( $\lt 5\,\%$ ) provided that a sufficiently large number of grains is used and that there is some means of preventing grains from rolling away (such as a rim or rough bed; see Appendix D). By contrast, most other experimental datasets reported in the literature, even well-controlled ones (Lajeunesse et al. Reference Lajeunesse, Malverti and Charru2010), do not report $S_{\!f}$ and/or  $\mu _s$ .

For those datasets mentioned in the previous paragraph that were not already discussed in § 2.4.1 on the spherical-grain data compilation, we use the values of $h$ reported in the respective studies, as justified by $h^\ast =h/d\gg 1$ (table 1). For more details and the calculation of the bed shear stress uncertainty, see Appendix B.

3.2.2. Grain-unresolved CFD-DEM simulations

Since our bedload flux model generalises (3.2), the transport load expression by Bagnold (Reference Bagnold1956), to (3.3), involving the additional parameter $\overline {f_L}/\overline {f_D}=1.0$ , we test it also for conditions to which (3.2) applies (i.e. $\overline {f_L}/\overline {f_D}=0$ ). For this purpose, we add data from grain-unresolved CFD-DEM simulations to the compilation, both from existing studies (Maurin et al. Reference Maurin, Chauchat and Frey2018; Monthiller Reference Monthiller2019; Pähtz & Durán Reference Pähtz, Liu, Xia, Hu, He and Tholen2020) and from our own new simulations, all of which were conducted using either of the numerical models described by Durán et al. (Reference Durán, Andreotti and Claudin2012), Maurin et al. (Reference Maurin, Chauchat, Chareyre and Frey2015) and Xie et al. (Reference Xie, Hu, Pähtz, He and Cheng2022) (table 1; for details, see Appendix E). They differ from grain-resolved simulations in the fact that they do not resolve the flow around single grains, and they therefore incorporate empirical relations for the buoyancy, drag and lift forces acting on grains, the latter of which are typically neglected (also here, i.e. $\overline {f_L}/\overline {f_D}=0$ ). They also allow us to perform thought experiments in which $C_D$ is treated as just another control parameter, since we can artificially modify $C_D$ through changing the drag force parameters (Appendix E); see LESCD and BLRANSCD in table 1. Another purpose for adding these data is to confront our bedload flux model with transport conditions across a much larger range of the control parameters $\varTheta$ , $\alpha$ and $h^\ast$ than that spanned by the experiments and grain-resolved simulations alone.

3.2.3. Streamwise fluid momentum balance in BLRANSx

The numerical model of Durán et al. (Reference Durán, Andreotti and Claudin2012) behind the datasets BLRANSx considers a flow in an inner turbulent boundary layer of infinite height (no sidewalls). To be consistent with (2.5), we slightly modified the streamwise fluid momentum balance employed by Durán et al. (Reference Durán, Andreotti and Claudin2012) to

(3.6) \begin{equation} {\mathrm{d}}_z\sigma ^f_{zx}=\rho _{\!f}\phi g\sin \alpha +\rho _{\!p}\phi f^{f\rightarrow p}_x,\quad \text{with} \quad \sigma ^f_{zx}(\infty )=\tau _b, \end{equation}

before conducting simulations, where $\rho _{\!p}\phi \boldsymbol{f}^{f\rightarrow p}$ is the fluid–particle interaction force per unit volume of the mixture. In fact, integrating (3.6) from zero to infinity, and combining the result with (C5) and $\sigma _{zx}(0)=\sigma ^f_{zx}(0)+\sigma ^p_{zx}(0)$ , leads precisely to (2.5). Note that for the numerical model of Durán et al. (Reference Durán, Andreotti and Claudin2012), the Shields number $\varTheta \equiv \tau _b/(\rho _{\!p}\tilde g_zd)$ based on $\tau _b$ in (3.6) is a prescribed fixed parameter.

3.2.4. Determination of $\varTheta$ for M18RANS, M19RANSx and LESx

Since the raw data required to determine the bed surface elevation, and thus $h$ and $\tau _b$ , are not available for the datasets M18RANS (Maurin et al. Reference Maurin, Chauchat and Frey2018) and M19RANSx (Monthiller Reference Monthiller2019) in table 1, based on the model of Maurin et al. (Reference Maurin, Chauchat, Chareyre and Frey2015), we infer $\varTheta$ from their reported transport-driving Shields number $\varTheta _\ast$ , as explained in Appendix F.