Role of turbulent kinetic energy modulation by particle–fluid interaction in sediment pick-up

Abstract

Reliable prediction of the erosion rate of sediment beds is important for many applications in coastal and river engineering. Theoretical understanding of empirically derived scaling relations is still lacking. This applies in particular for the scaling anomaly between low and high Shields number conditions. In this work, the erosion process is studied from the perspective of the phase-averaged turbulent kinetic energy (TKE) equations. The multi-phase TKE equations are written in a form that allows for a direct comparison with the TKE equation that appears for a stratified single-phase flow under the Boussinesq approximation. This reveals that next to buoyancy destruction, several other TKE modulation mechanisms become important at high Shields numbers and concentrations. Two scaling laws are derived for both moderate and high Shields numbers, and are tested against a wide range of experimental data.

1. Introduction

Reliable prediction of the erosion of non-cohesive sediment is important for assessment of the safety level of levees, prediction of coastal erosion, and optimization of dredging methods and equipment. Also, this information is important for the understanding of river systems and fluvial landforms. Although there exist several theoretical attempts to address the erosion problem (Cao 1997; Zhong, Wang & Ding 2011), understanding is still far from satisfactory, and the proposed scaling relations are mainly empirical.

The vertical particle flux from a sediment bed is defined as $F(y)=\rho _p\,\langle \alpha v_p\rangle$, where $\rho _p$ is the density of the solid particles, $\alpha$ is the instantaneous solid phase fraction, $v_p$ is the instantaneous vertical velocity of the particle phase obtained by volume averaging, and the brackets $\langle \cdot \rangle$ denote the Reynolds average, i.e. an ensemble average over multiple realizations of the flow (Pope 2000). It is helpful to decompose $F(y)$ into a pick-up flux $E_p(y)\geqslant 0$ and a settling flux $S(y)\geqslant 0$:

(1.1)\begin{equation} F=E_p-S. \end{equation}

Figure 1 illustrates the pick-up process of particles from an erodible bed. It is assumed that a horizontal flow profile of the fluid, $\langle u_f\rangle$, has been established before it reaches the erodible bed. As a result of particle entrainment, a solid phase fraction profile $\left \langle \alpha \right \rangle$ develops, while the interface gradually moves downwards with the erosion velocity $v_e$. It is also possible to consider a co-moving frame of reference that translates vertically with $v_e$. By mass conservation, $v_e$ can be related to the pick-up flux and settling flux at a given reference height $y_{r}$ according to

(1.2)\begin{equation} v_e=\frac{E_p(y_{r})-S(y_{r})}{\left(1-n_i-\left\langle \alpha(y_{r}) \right\rangle\right)\rho_p}, \end{equation}

where $n_i$ is the in situ porosity of the bed, and $\left \langle \alpha (y_{r}) \right \rangle$ is the mean solid phase fraction at the reference height. The reference height can be considered as the seabed boundary of a large-scale sediment transport model that typically incorporates turbulent suspension and particle settling, but strategically ignores near-bed transport phenomena. It is assumed that the pick-up flux is governed by the stress exerted on individual grains at the bed $\tau$, a critical stress required for the initiation of motion of particles $\tau _{cr}$, the particle size $d_p$, the viscosity of the fluid $\nu _f$, gravity $g$, the density of the solid particles $\rho _p$, and the density of the carrier fluid $\rho _f$, such that

(1.3)\begin{equation} E_p(y_{r})=f_1( \tau,\tau_{cr},d_p, \nu_f, g, \rho_p, \rho_f). \end{equation}

It is also assumed that the settling flux can be related to the terminal settling velocity of a single grain in stagnant water $v_p^\infty$, and $\left \langle \alpha (y_{r}) \right \rangle$ and $\rho _p$:

(1.4)\begin{equation} S(y_{r})=f_2(v_p^\infty, \left\langle \alpha(y_{r}) \right\rangle,\rho_p). \end{equation}

A very frequently used pick-up function was proposed by Van Rijn (1984):

(1.5)\begin{equation} \phi=\frac{E_p(y_{r})}{\rho_p\left(\varDelta g d_{50}\right)^{0.5}}=0.00033 D_*^{0.3}f_D\left(\frac{\theta-\theta_{cr}}{\theta_{cr}}\right)^{1.5}, \end{equation}

where $\varDelta =\rho _p/\rho _f-1$, $D_*$ is a dimensionless particle size defined as $D_*=d_{50}\sqrt [3]{{\varDelta g}/{\nu _f^2}}$, i.e. the square root of the ratio between the submerged weight of the particle and viscous stresses (Galileo number), $\theta =\tau /(\rho _p-\rho _f)gd_{50}$ is the Shields parameter (ratio of shear stress and submerged weight of particles), $\theta _{cr}$ is the critical Shields parameter related to the minimum shear stress that is required for the initiation of particle motion, and $f_D$ is a correction factor, with $f_D=1.0$ in the original formulation for $\theta _{cr}\leqslant \theta <1.0$. The normalization of the pick-up flux with velocity scale $\sqrt {\varDelta g d_{50}}$ originates from Einstein (1950). The critical Shields parameter $\theta _{cr}$ can be determined from the empirical Shields diagram (Shields 1936). Buffington & Montgomery (1997) presented an overview of different techniques to determine $\theta _{cr}$. There exists quite some scatter in the measurement data that can be attributed to the preparation of the bed, particle size distributions and the degree of exposure of the particles at the top of the bed to the fluid flow (Miedema 2012). The criteria for initiation of motion have also been defined in terms of other parameters. For example, Kudrolli, Scheff & Allen (2016) considered a critical shear rate condition, and Maldonado & de Almeida (2019) studied the impulse acting on the particles at the bed. In the regime $\theta \gg \theta _{cr}$, (1.5) implies that the pick-up flux becomes proportional to $\theta ^{1.5}$ or in terms of the friction velocity, $u_*=\sqrt {\tau /\rho _f}$, $E_p\propto u_*^3$, i.e. linearly proportional to the power of the eroding flow. The fact that the critical Shields number $\theta _{cr}$ is present in the denominator of the term in parentheses on the right-hand side of (1.5) suggests that the critical Shields parameter remains important in the limit of high Shields numbers, $\theta \gg \theta _{cr}$. Other pick-up relations were proposed that do not reflect this behaviour in the limit of high Shields values (Fernandez Luque & Van Beek 1976; Nakagawa & Tsujimoto 1980). The pick-up model of Einstein (1950) does not contain a critical shear stress, but contains a probability function $P(\theta )$ that describes the fraction of time that the lift force on a particle at the bed exerted by the turbulent flow is sufficient to overcome the weight of the particle. The normalized pick-up flux according to Einstein (1950) is

(1.6)\begin{equation} \frac{E_p}{\rho_p\sqrt{\varDelta d_{50}g}}=\beta\,P(\theta), \end{equation}

where $\beta$ is a non-dimensional parameter. The probability function $P(\theta )$ decays rapidly in the limit $\theta \rightarrow 0$ and converges quickly to unity for $\theta \gtrsim 0.15$. Yalin (1977) followed a similar approach but considered the friction velocity $u_*$ as the relevant velocity scale for sediment pick-up instead of the velocity scale $\sqrt {\varDelta d_{50}g}$ as considered by Einstein (1950):

(1.7)\begin{equation} \frac{E_p}{\rho_p u_*}=\beta_2\,P_2(\theta), \end{equation}

where $P_2$ is a modified probability function. The value $\theta =0.15$ can be considered as a reasonable separator between the regime for which the probability function is smaller than 1 and the regime for which the probability function is equal to 1, since $P(\theta =0.15)\approx 0.9$ in the model of Einstein (1950) and $P_2(\theta =0.15)\approx 0.99$ in the model of Yalin (1977). Based on force measurements on a single grain near a smooth bed, Schmeeckle, Nelson & Shreve (2007) showed the resemblance of the measured force with the Saffman lift force, which supports the analysis of Einstein (1950) and Yalin (1977). Also, Garcia & Parker (1993) proposed a scaling law that does not contain a critical Shields condition; their pick-up relation contains the terminal settling velocity of the particles, the friction velocity and the particle Reynolds number. They observed a very steep scaling of the pick-up flux with the friction velocity $E_p\propto u_*^5$. In addition, Cheng & Emadzadeh (2016) proposed an alternative scaling relation without a critical Shields parameter.

Figure 1. Definition sketch of pick-up flux $E_p$, settling flux $S$, and total particle flux $F$.

Winterwerp et al. (1992) executed experiments with hyper-concentrated flows over an erodible bed ($\theta >1.0$). They observed a significant deviation from the scaling relation (1.5). Based on their data, they proposed a scaling law $E_p\propto \theta ^{0.5}$. The deviations were attributed to the high near-bed concentrations that are absent in the experiments of Van Rijn (1984). Bisschop (2018) conducted experiments at extremely high Shields parameters, up to $\theta =329$, and high near-bed concentrations. Recently, van Rijn, Bisschop & van Rhee (2019) performed a re-analysis of these data and proposed a correction of scaling (1.5) by setting $f_D={1}/{\theta }$ for $\theta >1$. According to van Rijn et al. (2019), this factor takes three processes into account that are relevant for the high Shields range: damping of turbulence in the near-bed region where the solid concentrations are extremely high, the increase of the mixture kinematic viscosity, and the increase of apparent shear resistance as a result of dilatancy of the sand bed near the interface. Van Rhee (2010) demonstrated that this dilatancy effect can become important for a large ratio between the erosion velocity and the hydraulic conductivity of the bed, $v_e/k \gg 1$, in combination with densely packed beds (low porosity $n_i$).

To elucidate the key difference in the observations, it is helpful to reconsider the studies mentioned above in more detail. First, the erosion experiments of Van Rijn (1984), Okayasu, Fujii & Isobe (2010), Cheng & Emadzadeh (2016) and Bisschop (2018) cover a wide range of flow conditions, sediment sizes and measurement techniques (table 1). Van Rijn (1984) and Cheng & Emadzadeh (2016) used sediment lifts to keep the bed aligned with the floor of the flumes. The equilibrium speed of the lifts can be interpreted as the erosion velocity $v_e$, defined in (1.2). Okayasu et al. (2010) constructed a flume with two sections. In the first section, particles were fixed to the flume floor, and in the second section, a sand bed could freely erode. The position of the bed was tracked near the starting point of the erodible section with a laser pointer. The bed shear stress was derived from velocity measurements in the first section and the application of standard boundary layer theory for rough walls. Bisschop (2018) used a closed-circuit flume where the flow velocity can be increased up to $10\ {\rm m}\ {\rm s}^{-1}$. The position of the bed was measured at approximately 3 m from the starting point of the erodible section using a set of conductivity probes mounted in the side wall of the flume and high-speed camera images. The conductivity probes were calibrated to obtain the absolute concentration profiles above the bed. The shear stress was estimated by measuring the pressure difference over the erodible section (over a length $\approx$6 m), which was corrected for acceleration effects and the wall shear stress contributions from the side and top walls of the flume.

Table 1. Parameter ranges in the erosion experiments ($d_{50}$, $v_p^\infty$, $u_*$), corresponding dimensionless numbers ($\theta$, $D_*$), entry flow at erosion measurement section being either clear-water flow (CWF) or hyper-concentrated flow (HCF), and typical concentrations in the near-bed region $\alpha _r$.

Figure 2 shows the normalized pick-up flux versus the Shields parameter of all the experimental studies mentioned above. It is seen that the scaling at relatively low Shields values differs markedly from the high near-bed concentration regimes addressed by Winterwerp et al. (1992) and Bisschop (2018). Hence it is useful to distinguish clear-water flow (CWF) and hyper-concentrated flow (HCF), with a maximum near-bed concentration <0.01 and much larger than 0.01, respectively. The normalization of the pick-up flux as suggested by Einstein (1950), $\phi =E_p/\rho _p\sqrt {gd_p\varDelta }$, does not collapse all the data onto a single master curve $f(\theta )$. For Shields parameters lower than unity, a power-law trend with exponent 1.5 seems acceptable even though not perfect. In the HCF case, a power law with exponent 0.5 better describes the experimental trend. The fact that (i) the data do not collapse perfectly on master curves, and (ii) the power law exhibits an abrupt modification between CWF and HCF cases, indicates that some physical processes are not captured by the scaling law $\phi =f(\theta )$, and that further analysis is needed to explain and understand the pick-up flux. This work addresses the observed scaling behaviour of the pick-up flux from the perspective of the turbulent kinetic energy (TKE) budgets of the particle phase and the fluid phase. The main focus is to explain the striking difference in the scaling of the pick-flux in the CWF and HCF regimes by expected changes in the relative contributions of various turbulence modulation terms in the TKE equations. To analyse the CWF case, we model the fluid–particle system as a stratified single-phase flow under the Boussinesq approximation. The HCF case is analysed by using the phase-averaged TKE equations for fluid–particle systems. New scaling laws are retrieved from this analysis, which are compared with existing experimental data.

Figure 2. Scaling anomaly in the normalized pick-up flux and Shields number between CWF and HCF cases. CWF experiments: Van Rijn (1984), Okayasu et al. (2010) and Cheng & Emadzadeh (2016). HCF experiments: Winterwerp et al. (1992) and Bisschop (2018).

2. Analysis of CWF experiments at moderate Shields numbers

For the experimental cases with CWF over erodible beds at moderate Shields numbers $0.15<\theta <1$, the near-bed concentration is relatively low. The Stokes number in the near-bed region can be defined as $St=t_p/t_f$, with a fluid time scale $t_f=\kappa y_{r}/u_*$, that is, the ratio between the characteristic length scale of eddies in the near-bed region and the friction velocity, with $\kappa$ representing the von Kármán constant, and $t_p$ the particle drag time scale; see Appendix A. A basic estimate for $t_p$ can be obtained by considering the terminal settling velocity of a single particle in stagnant water $v_p^\infty$ (Greimann & Holly 2001):

(2.1)\begin{equation} t_p=\frac{v_p^\infty \rho_p}{|g|\,(\rho_p-\rho_f)}. \end{equation}

Using this estimate, we found that for the CWF experiments, the Stokes number is typically smaller than 1 for $y_{r}>2d_p$. This implies that the relative velocity between the phases is negligible, and it suffices to consider a single-phase formulation. Because the concentrations are low, the Boussinesq approximation holds.This means that the effect of density differences can be modelled by source terms that appear on the right-hand sides of the momentum and TKE equations, while the density in the inertia terms is considered to be constant. For more details, see Appendix A.

To study the effect of self-stratification, three local characteristic numbers need to be examined: the flux Richardson number, buoyancy Froude number and buoyancy Reynolds number. First, the local flux Richardson number (Turner 1979) is defined as the ratio between the buoyancy destruction and the shear production of TKE:

(2.2)\begin{equation} Ri_f(y)={-}\frac{|g|\,\langle{v^\prime \rho^\prime}\rangle }{\rho_0\, \langle{u^\prime v^\prime}\rangle\,\partial {\left\langle u \right\rangle}/\partial y}={-}\frac{\varDelta\,|g|\,\langle{v^\prime \alpha^\prime}\rangle }{\langle{u^\prime v^\prime}\rangle\,\partial \left\langle u \right\rangle/\partial y}, \end{equation}

where $\left \langle u \right \rangle$ is the mean horizontal velocity, $u'$ and $v'$ are the velocity fluctuations in the horizontal and vertical directions, respectively, $\rho '$ is the density fluctuation of the mixture, $\alpha '$ is the fluctuation of the particle phase fraction, and $\rho _0\approx \rho _f$ in the dilute limit.

The gradient of the mean velocity can be modelled by the Monin–Obukhov similarity theory (Monin & Obukhov 1954)

(2.3)\begin{equation} \frac{\partial \langle u\rangle}{\partial y}=\frac{u_*}{\kappa y}\,\psi_m(\zeta), \end{equation}

where $\zeta =y/L$ is the non-dimensional height above the bed, $L$ is the Obukhov length (Obukhov 1946)

(2.4)\begin{equation} L=\frac{u_*^3}{\varDelta\,|g| \left\langle v'\alpha' \right\rangle\kappa}, \end{equation}

and $\psi _m$ is the stability function for the momentum flux. The Obukhov length is a profile-dependent parameter that characterizes the strength of the stratification. It represents the theoretical height where $Ri_f=1$ based on a neutrally stratified $\psi _m=1$ velocity profile. Substitution of (2.3) in (2.2) gives

(2.5)\begin{equation} Ri_f(\zeta)=\frac{\varDelta\,|g|\,\langle v'\alpha'\rangle\,\kappa y}{u_*^3\,\psi_m({\zeta})}=\frac{\zeta}{\psi_m(\zeta)}. \end{equation}

The second control parameter is the buoyancy Froude number (Lesieur 1990)

(2.6)\begin{equation} \displaystyle Fr=\left(\frac{l_B}{l_e}\right)^{2/3}, \end{equation}

where $l_e=u_*^3/\varepsilon$ is a turbulent flow scale associated with the most energetic eddies in the boundary layer at height $y$, and $l_B=(\varepsilon /N^3)^{1/2}$ is the Osmidov length scale (Osmidov 1975), with $\varepsilon =u_*^3\,\psi _m(\zeta )/\kappa y$ the turbulent dissipation, and $N$ the Brunt–Vaïsala frequency (Turner 1979; Lesieur 1990) given by

(2.7)\begin{equation} \displaystyle N^2={-} \frac{|g|}{{\rho_0}}\,\frac{{\rm d} \langle\rho\rangle}{{{\rm d} y}}. \end{equation}

Parameter $N$ is also known as the buoyancy frequency, which is the frequency of an internal wave in a stably stratified fluid. The Osmidov scale $l_B$ is a local scale for the stratification, and it represents the maximum eddy length scale that can develop in the presence of internal waves such that the turnover time of the eddies at that scale equals $1/N$.

The density profile $\left \langle \rho \right \rangle =\left \langle \alpha \right \rangle \rho _p + (1-\left \langle \alpha \right \rangle )\, \rho _f$ is a function of the local phase fraction such that

(2.8)\begin{equation} \displaystyle \frac{{\rm d} \langle\rho\rangle}{{\rm d} y} = \rho_f \varDelta\,\frac{{\rm d} \langle\alpha\rangle}{{{\rm d} y}}. \end{equation}

The gradient of the phase fraction can also be described by the Monin–Obukhov similarity theory:

(2.9)\begin{equation} \frac{\partial \langle\alpha\rangle}{\partial y}=\frac{\langle{v'\alpha'}\rangle}{u_*}\,\frac{\psi_\alpha(\zeta)}{\kappa y}, \end{equation}

where $\psi _\alpha$ is the stability function for $\alpha$. Since for the clear-water case $\rho _0\approx \rho _f$, it follows that

(2.10)\begin{equation} N^2=\frac{|g|\,\varDelta\,\langle{v'\alpha'}\rangle\,\psi_\alpha(\zeta)}{\kappa y u_*}. \end{equation}

This implies that the local buoyancy Froude number and local flux Richardson number can be related by

(2.11)\begin{equation} Fr=\frac{1}{\sqrt{Ri_f}}\,\sqrt{\frac{\psi_m}{\psi_\alpha}}, \end{equation}

where it is expected that $\psi _m\approx \psi _\alpha$. This signifies that self-stratification in the near bed region can be interpreted in terms of both buoyancy destruction of TKE and stability constraints imposed by internal waves associated with the vertical density gradient above the bed.

The third control parameter of stratified flows is the buoyancy Reynolds number (Venayagamoorthy & Koseff 2016) defined as

(2.12)\begin{equation} Re_b=\frac{\varepsilon}{\nu_f N^2}. \end{equation}

The buoyancy Reynolds number can be related to the ratio of the Osmidov length scale $l_B$ and the Kolmogorov length scale $l_\eta =(\nu _f^3/\varepsilon )^{1/4}$ since $l_B/l_\eta =Re_b^{3/4}$. If the separation between these scales is too small, then the turbulent cascade cannot develop, which in turn drives the flow towards a laminar state and reduces particle entrainment. Figure 3 illustrates this argument. Using (2.12), $\varepsilon =u_*^3\,\psi _m(\zeta )/\kappa y$ and (2.10), it follows that

(2.13)\begin{equation} Re_b=\frac{Re_{\tau}(y)}{Ri_f(y)}\,\frac{\kappa}{\psi_{\alpha}(\zeta) }, \end{equation}

where $Re_{\tau }=u_* y/\nu _f$ is the friction Reynolds number. In the regime of minimum scale separation between $l_B$ and $l_\eta$, $Re_b$ reaches a critical value, and according to (2.13), this will impose a constraint on the local flux Richardson number $Ri_f^c\propto Re_\tau \kappa /\psi _\alpha$. By virtue of (2.5), the definition of the pick-up flux, $\phi =\langle {v'\alpha '}\rangle /\sqrt {\varDelta g d_p}$ and $\psi _\alpha \approx \psi _m$, for small $\zeta$, it then follows that

(2.14)\begin{equation} \phi\propto\lim_{\zeta \to 0}\frac{\psi_m(\zeta)}{\psi_\alpha(\zeta)}\,Re_p\,\theta^{1.5}=Re_p\,\theta^{1.5}, \end{equation}

where $Re_p=d_pu_*/\nu _f$ is the particle Reynolds number based on the friction velocity.

Figure 3. Constraint on the ratio between the Osmidov length scale $l_B$ and the Kolmogorov scale $l_\eta$. Increasing the pick-up flux results in a reduction of $l_B$ by buoyancy until a critical value $l_B/l_\eta$ (i.e. $Re_b$) is reached. This argument corresponds to scaling (2.14).

In dilute suspensions, the inter-particle distance allows enough space for TKE dissipation events to evolve, as in CWF. However, it is possible that a finite-size effect sets in when $l_B\approx d_p$: the energy carrying eddies might no longer be able to induce a net lift or drag force on the particles at the interface, at least not in a statistical-averaged sense; see figure 4. From this constraint on the Osmidov length scale and (2.13), it follows that

(2.15)\begin{equation} Ri_f^c\propto \left(\frac{\kappa y_r}{d_p} \right)^{4/3}\frac{1}{\psi_\alpha\psi_m^{1/3}}, \end{equation}

and for the non-dimensional pick-up flux at small $\zeta$,

(2.16)\begin{equation} \phi\propto \frac{\psi_m^{2/3}}{\psi_\alpha}\left(\frac{\kappa y_r}{d_p} \right)^{1/3}\theta^{1.5}\approx \left(\frac{\kappa y_r}{d_p} \right)^{1/3}\theta^{1.5}. \end{equation}

This leaves a weak dependence on the choice of the reference height.

Figure 4. Constraint on the ratio between the Osmidov length scale $l_B$ and the particle size $d_p$. Increasing the pick-up flux results in a reduction of $l_B$ by buoyancy until a critical value $l_B/d_p$ is reached. This line of reasoning results in scaling (2.16).

In conclusion, in the CWF experiments at moderate Shields numbers $0.15<\theta <1.0$, by a self-stratification argument, two scaling laws can be expected, (2.14) and (2.16) for sub-Kolmogorov and finite-size particles, respectively. We will test the scaling relations, as well as the empirical relation given by (1.5), in § 4. Note that they all have the same $\theta ^{1.5}$ proportionality, but differ in the pre-factor in front of this term, notably the dependence on the particle diameter $d_p$.

3. Analysis of HCF experiments at high Shields numbers

For HCF over erodible beds at high Shields numbers, it is required to address more generic equations for the TKE of dispersed multi-phase flows.

These equations appeared in various forms in literature. In this work, we use the framework of Fox (2014) and adopt the same notations here for convenience. The two-fluid model of Fox (2014) is restricted to drag forces only, but we think that this is not a limitation for this study. Although the Saffman lift force is important for the pick-up of individual particles from the bed (see § 1), we assume that it is sufficient to restrict our near-bed flow analysis to the drag force. Jha & Bombardelli (2010) demonstrated that adding the lift force in a two-fluid model for particle-laden channel flows modifies the flow dynamics only marginally.

To study the departures from stratified single-phase flow under the Boussinesq approximation, it is helpful to combine the balances for the TKE of the fluid and particle phases. In the equations for the fluid and particle TKE as derived by Fox (2014), the buoyancy term, as it appears under the Boussinesq approximation, is difficult to recognize. In Appendix A, it is shown that after multiple manipulations of the equations of Fox (2014), the combined TKE balance of the two phases can be written as follows:

(3.1)\begin{align} &\underbrace{\frac{\partial\rho_f\left\langle 1-\alpha \right\rangle k_f+\rho_p\left\langle \alpha \right\rangle k_p}{\partial t}}_{\textit{time change}} +\boldsymbol{\nabla}\boldsymbol{\cdot}\biggl[\underbrace{\rho_f\left\langle 1-\alpha \right\rangle\left\langle \boldsymbol{u}_f \right\rangle_fk_f+ \rho_f\left\langle 1-\alpha \right\rangle\frac{1}{2}\left\langle \boldsymbol{u}_f'''\boldsymbol{u}_f'''\boldsymbol{\cdot}\boldsymbol{u}_f''' \right\rangle_f}_{\textit{transport by mean fluid velocity and fluid velocity fluctuations}} \nonumber\\ &\quad +\underbrace{ {\left\langle 1-\alpha \right\rangle}\left\langle p_f'\boldsymbol{u}_f''' \right\rangle_f-\left\langle \boldsymbol{\sigma}_f\boldsymbol{\cdot}\boldsymbol{u}_f''' \right\rangle}_{ \textit{transport by fluid stresses}} \nonumber\\ &\quad +\underbrace{\rho_p\left\langle \alpha \right\rangle\left\langle \boldsymbol{u}_p \right\rangle_p k_p+\rho_p\left\langle \alpha \right\rangle\frac{1}{2}\left\langle \boldsymbol{u}_p''\boldsymbol{u}_p''\boldsymbol{\cdot} \boldsymbol{u}_p'' \right\rangle_p}_{\textit{transport by mean particle velocity and particle velocity fluctuations}} \nonumber\\ &\quad +\underbrace{\left\langle \alpha \right\rangle\left\langle p_f'\boldsymbol{u}_p'' \right\rangle_p+\rho_p\left\langle \alpha \right\rangle\left\langle \boldsymbol{\mathsf{P}}_p \boldsymbol{\cdot} \boldsymbol{u}_p'' \right\rangle_p\biggr]}_{\textit{transport by fluid pressure fluctuations and particle stresses}} \nonumber\\ &\quad =\underbrace{-\rho_f\left\langle 1-\alpha \right\rangle\left\langle \boldsymbol{u}_f'''\boldsymbol{u}_f''' \right\rangle_f: \boldsymbol{\nabla}\left\langle \boldsymbol{u}_f \right\rangle_f-\rho_p\left\langle \alpha \right\rangle\left\langle \boldsymbol{u}_p''\boldsymbol{u}_p'' \right\rangle_p:\boldsymbol{\nabla}\left\langle \boldsymbol{u}_p \right\rangle_p}_{\textit{shear production in both phases, conversion MKE-TKE}} \nonumber\\ &\qquad +\underbrace{\left\langle p_f'\boldsymbol{\nabla}\boldsymbol{\cdot}\alpha'\left(\left\langle \boldsymbol{u}_f \right\rangle_f- \left\langle \boldsymbol{u}_p \right\rangle_p\right) \right\rangle}_{\textit{compression effect, conversion MKE-TKE}} \nonumber\\ &\qquad +\underbrace{\frac{\rho_p\left\langle \alpha \right\rangle}{ t_p}\,\frac{\left\langle \alpha'\boldsymbol{u}_f' \right\rangle}{\left\langle \alpha \right\rangle\left\langle 1- \alpha \right\rangle}\boldsymbol{\cdot}\left(\left\langle \boldsymbol{u}_p \right\rangle_p-\left\langle \boldsymbol{u}_f \right\rangle_f-\frac{\left\langle \alpha'\boldsymbol{u}_f' \right\rangle}{\left\langle \alpha \right\rangle\left\langle 1-\alpha \right\rangle}\right)}_{\textit{drag conversion MKE-TKE, contains buoyancy via momentum equations, see (A38)}} \nonumber\\ &\qquad \underbrace{-\left\langle \boldsymbol{\nabla} \boldsymbol{u}_f''':\boldsymbol{\sigma}_f \right\rangle}_{\textit{viscous dissipation, conversion TKE-heat}}+\underbrace{\rho_p\left\langle \alpha \right\rangle\left\langle \boldsymbol{\mathsf{P}}_p: \boldsymbol{\nabla} \boldsymbol{u}_p'' \right\rangle_p}_{\textit{conversion TKE-$\varTheta_p$}}\nonumber\\ &\qquad - \underbrace{\frac{\rho_p\left\langle \alpha \right\rangle}{ t_p}\left\langle \left|\boldsymbol{u}_f''-\boldsymbol{u}_p''\right|^2 \right\rangle_p}_{ \textit{drag dissipation, conversion TKE-heat}}, \end{align}

where $k_f$ and $k_p$ are the TKE, in kinematic units, of the fluid and particle phases, respectively, $\boldsymbol {u}_f$ and $\boldsymbol {u}_p$ are the velocities of each phase, $\left \langle \cdot \right \rangle$ represents standard Reynolds averaging, $\left \langle \cdot \right \rangle _f$ and $\left \langle \cdot \right \rangle _p$ are conditional averages over each phase, $A'$, $A''$ and $A'''$ denote fluctuations with respect to the Reynolds average, particle phase average and fluid phase average, respectively, $\boldsymbol{\mathsf{P}}_p$ is the particle stress, $\boldsymbol {\sigma }_f$ is the viscous stress, $p_f$ is the fluid pressure, and $t_p$ is a constant time scale associated with drag.

The left hand-side of (3.1) contains the time change of TKE and transport of TKE by the mean flows, turbulent fluctuations, pressure fluctuations, collisional stress and viscous stress. Similar transport terms appear in the TKE equation for stratified single-phase flow.

The first two terms on the right-hand side of (3.1) represent the shear production of the fluid and particle phases. Similar to single-phase flow, they describe the transfer from kinetic energy of the mean flow (or mean kinetic energy, MKE) to TKE. The next two terms describe the transfer from MKE to TKE by pressure coupling and drag forces. Using the momentum balances of the fluid and particle phases, the drag contribution can be related to the buoyancy term as it appears in the TKE equation under the Boussinesq approximation; see (A5) in Appendix A. The last three terms describe the conversion from TKE into heat and granular temperature. After the conversion into granular temperature, energy can be further converted into heat by drag and inelastic collisions; see Fox (2014). Figure 5 gives an overview of the energy flows between the different variables. In the context of this study, the relation between the conversion of MKE to TKE by the mean shear, the conversion of TKE to MKE via drag, and the conversion of MKE to potential energy (PE) is important. For increasing particle concentrations, granular dissipation $\epsilon _p$, and drag dissipation become more relevant, next to the viscous dissipation of TKE $\epsilon _f$.

Figure 5. Energy flows between PE, MKE, TKE, granular temperature $\varTheta _p$, temperature of the particles $T_p$, pseudo-TKE $\varTheta _f$ (turbulence in the wake of the particles) and temperature of the fluid $T_f$, where $e$ is the coefficient of restitution between the particles, $\epsilon _f$ is the viscous dissipation of fluid TKE, $\epsilon _p$ is the granular dissipation of particle TKE, and $\epsilon _D$ is the small-scale viscous dissipation associated with pseudo-turbulence (Fox 2014).

Now the TKE budget is written in a convenient form to study the departures from the Boussinesq approximation, there is an opening to retrieve possible scaling laws for the pick-up flux in the HCF regime at high Shields numbers. In the previous section, it is anticipated that the pick-up flux at moderate Shields numbers and low near-bed concentrations is controlled by the balance between the shear production, viscous dissipation and buoyancy destruction terms. For increasing near-bed concentrations and Shields numbers, the magnitude of each term in (3.1) will likely change. For simplicity, we assume that although the relative contribution of each transport term will vary in the HCF regime, the key deviations that control the pick-up flux stem from differences in the source terms on the right-hand side of (3.1).

First, we need to estimate the combined shear production of both phases. The assumption is that in the flow direction, the macroscopic particle/fluid slip velocity is small at the length scale ($\kappa y$) of the mean flow and of the turbulent fluctuations in the energy-containing range. In then follows that $\left \langle u_p \right \rangle _p\approx \left \langle u_f \right \rangle _f$ and $\left \langle u_p''v_p'' \right \rangle _p\approx \left \langle u_f''v''_f \right \rangle _f$. The combined shear production term from both phases then becomes

(3.2)\begin{align} &-\rho_f\left\langle 1-\alpha \right\rangle\left\langle u_f'''v_f''' \right\rangle_f\frac{\partial \left\langle u_f \right\rangle_f}{\partial y}-\rho_p\left\langle \alpha \right\rangle\left\langle u_p''v_p'' \right\rangle_p\frac{\partial \left\langle u_p \right\rangle_p}{\partial y}\nonumber\\ &\quad \approx \rho_m\left\langle u_f'''v_f''' \right\rangle_f\frac{\partial \left\langle u_f \right\rangle_f}{\partial y}. \end{align}

Next, we estimate the turbulent shear stress by

(3.3)\begin{equation} -\rho_m\left\langle u_f'''v_f''' \right\rangle_f={-}\rho_m\left\langle u_f'v_f' \right\rangle_f\approx \rho_m u_*^2, \end{equation}

which is expected to hold in the logarithmic (or constant-stress) layer as long as the Reynolds number is sufficiently high. Furthermore, assuming that the fluid velocity profile obeys the standard logarithmic scaling as in single-phase flow, the velocity gradient becomes equal to

(3.4)\begin{equation} \frac{\partial \left\langle u_f \right\rangle_f}{\partial y}\approx \frac{u_*}{\kappa y}. \end{equation}

The validity of the logarithmic scaling for the HCF regime has been confirmed by velocity measurements by Winterwerp et al. (1990) and Bisschop (2018). Our final estimate of the shear production term becomes

(3.5)\begin{equation} -\rho_f\left\langle 1-\alpha \right\rangle\left\langle u_f'''v_f''' \right\rangle_f\frac{\partial \left\langle u_f \right\rangle_f}{\partial y}-\rho_p\left\langle \alpha \right\rangle\left\langle u_p''v_p'' \right\rangle_p\frac{\partial \left\langle u_p \right\rangle_p}{\partial y}\approx \rho_m\frac{u_*^3}{\kappa y}, \end{equation}

with $\rho _m=\rho _f+(\rho _p-\rho _f)\left \langle \alpha \right \rangle$; for sediment-water flows, it typically holds that $\rho _p=2.6\rho _f$ so $\rho _m=\rho _f(1+1.6\left \langle \alpha \right \rangle )$. The shear production is balanced by the other terms: viscous dissipation, transfer to granular temperature, conversion of MKE to TKE, modified buoyancy destruction and drag dissipation. We expect that in the HCF regime, the transfer to granular temperature and drag dissipation will start to consume the largest portion of the TKE production, while changes in the MKE–TKE conversion have secondary effects.

From this conjecture, let us try to find some estimates for the pick-up flux in the HCF regime. The transfer to granular temperature term can be decomposed in two parts:

(3.6)\begin{equation} {\rho_p\left\langle \alpha \right\rangle\left\langle \boldsymbol{\mathsf{P}}_p:\boldsymbol{\nabla} \boldsymbol{u}_p'' \right\rangle_p}=\left\langle p_p\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}_p'' \right\rangle-\left\langle \boldsymbol{\sigma}_p:\boldsymbol{\nabla}\boldsymbol{u}_p'' \right\rangle, \end{equation}

where $p_p$ is the particle pressure, and $\boldsymbol{\sigma} _p$ is the particle shear stress. The pressure dilatation term in (3.6) can be modelled as

(3.7)\begin{equation} \left\langle p_p\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_p'' \right\rangle=\left\langle [1+4\alpha\, g_0(\alpha)]\varTheta\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_p'' \right\rangle_p, \end{equation}

where $g_0(\alpha )$ is the radial distribution function, and $\varTheta$ is the granular temperature (Fox 2014). In the dilute case, this term is usually neglected, with the argument that $\varTheta$ varies over the integral length scales and thus does not correlate with gradients of velocity fluctuations that live on a smaller scale. It is expected that a diverging flow field $\boldsymbol {u}_p''$ is associated with a lower concentration ($\alpha '<0$). As a result, the dilation term is likely negative at high particle concentrations. The second term on the right-hand side of (3.6) can be recognized as the dissipation rate of particle TKE by the particle shear stress. Combining the effect of pressure dilatation and shear provides the following definition of the dissipation rate of particle TKE:

(3.8)\begin{equation} \left\langle \boldsymbol{\sigma}_p:\boldsymbol{\nabla} \boldsymbol{u}_p'' \right\rangle-\left\langle p_p\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_p'' \right\rangle=\left\langle \alpha \right\rangle\rho_p{\epsilon}_p, \end{equation}

where it should be noted that this dissipation term is not a direct conversion into heat, but the particle TKE first converts into granular temperature and then converts into heat by particle collisions and particle–fluid drag on the level of individual granular fluctuations. On the other hand, the drag dissipation term on the right-hand side of (3.1) is a direct conversion into heat by particle–fluid drag; see figure 5.

To estimate the relative contribution of the drag dissipation and $\epsilon _p$, it is helpful to consider the viscous time scale $t_\eta$, the particle-drag time scale $t_p$ and the turnover time of the most energetic eddies in the near-bed region $t_e=\kappa y/u_*$. The viscous time scale can be estimated as

(3.9)\begin{equation} t_\eta=\left(\frac{\mu_{m}}{\rho_m \epsilon_m}\right)^{1/2}, \end{equation}

where the subscript $m$ indicates variables defined for the particle–fluid mixture. Typically, the effective shear viscosity of the mixture can be modelled as $\mu _m=\mu _f(1-{\left \langle \alpha \right \rangle }/{\left \langle \alpha \right \rangle _m})^{-2}$ (Maron & Pierce 1956; Guazzelli & Pouliquen 2018), with $\left \langle \alpha \right \rangle _m\approx 0.6$ the maximum flowable packing fraction. Using $\epsilon _m=u_*^3/\kappa y$ and $\rho _m=\rho _f(1+\varDelta \left \langle \alpha \right \rangle )$, it follows that

(3.10)\begin{equation} t_{\eta}=\left(\frac{\left(1-\dfrac{\left\langle \alpha \right\rangle}{\left\langle \alpha \right\rangle_m}\right)^{{-}2}}{1+\varDelta\left\langle \alpha \right\rangle}\right)^{1/2}\frac{\kappa y}{u_*}\,\frac{1}{\sqrt{\kappa \,Re_\tau}}. \end{equation}

Now if $t_{\eta }< t_p< t_e$, then it can be anticipated that the turbulent cascade breaks down at the time scale $t_p$, and drag dissipation is the most important dissipation mechanism. If $t_p< t_\eta < t_{e}$, then the turbulent cascade is halted by viscous dissipation, and at high particle concentrations it is expected that $\epsilon _p$ is at least of the same order as $\epsilon _f$ or higher (higher for higher concentration). By using (3.10), the constraint $t_\eta < t_p$ can be written as

(3.11)\begin{equation} \left(\frac{\left(1-\dfrac{\left\langle \alpha \right\rangle}{\left\langle \alpha \right\rangle_m}\right)^{{-}2}}{1+ \varDelta\left\langle \alpha \right\rangle}\right)^{1/2}\frac{1}{\sqrt{\kappa \,Re_\tau}}< St, \end{equation}

which sets a lower bound on the Stokes number, and likewise the constraint $t_p< t_\eta$ sets an upper bound on the Stokes number. This implies that the Stokes number controls the relative contribution of the drag dissipation (high Stokes numbers) and $\epsilon _p$ (low Stokes numbers) in the HCF regime.

Using estimates obtained for statistically steady-state, homogeneous and isotropic turbulence (Fox 2014), we derive in Appendix B the scaling for the dissipation of particle TKE

(3.12)\begin{equation} \epsilon_p=St\left(\frac{2}{St+2}\right)^2\frac{k_f}{t_p}, \end{equation}

and for the dissipation by drag forces,

(3.13)\begin{equation} \frac{1}{t_p}\left\langle \left|\boldsymbol{u}_f''-\boldsymbol{u}_p''\right|^2 \right\rangle_p=2\,\frac{k_f}{t_p}\left(\frac{St}{2+St}\right)^2. \end{equation}

The ratio $R$ of (3.12) and (3.13) is

(3.14)\begin{equation} R=\frac{2}{St}, \end{equation}

which indicates that at high Stokes numbers, drag dissipation is dominant, and at low Stokes numbers, turbulent dissipation of particle TKE overcomes the drag dissipation. This is in line with estimate (3.11) obtained from the analysis of the relevant time scales.

For the flow in the near-bed region, the Stokes number can be estimated as $St={u_*t_p}/{\kappa y}$. Estimate (3.14) then suggests that close to the bed, drag dissipation is dominant, while the importance of particle TKE dissipation increases when moving from the bed. At small $y$ and thus high Stokes numbers, $St\gg 1$, the dominant balance is between drag dissipation and shear production, yielding

(3.15)\begin{equation} \frac{\rho_p\left\langle \alpha \right\rangle}{t_p}\left\langle \left|\boldsymbol{u}_f''-\boldsymbol{u}_p''\right|^2 \right\rangle_p=2\rho_p\left\langle \alpha \right\rangle\frac{k_f}{t_p}\left(\frac{St}{2+St}\right)^2=\rho_m\,\frac{u_*^3}{\kappa y}. \end{equation}

Since $k_f\propto u_*^2$, and using $St={u_*t_p}/{\kappa y}$, it follows that for $St\gg 1$,

(3.16)\begin{equation} \left\langle \alpha \right\rangle\propto \frac{(2+St)^2}{St}\approx St, \end{equation}

which indicates that the near-bed concentration increases with the Stokes number. At higher regions in the profile, such that $St\ll 1$, the key balance is between dissipation of particle TKE and shear production,

(3.17)\begin{equation} \left\langle \alpha \right\rangle\rho_p\,St\left(\frac{2}{St+2}\right)^2\frac{k_f}{t_p}=\rho_m\,\frac{u_*^3}{\kappa y}. \end{equation}

This indicates that for $St\ll 1$,

(3.18)\begin{equation} \left\langle \alpha \right\rangle\propto (St+2)^2\propto St^0. \end{equation}

From the limiting cases (3.16) and (3.18), the following scaling can be anticipated for the near-bed region:

(3.19)\begin{equation} \left\langle \alpha \right\rangle\propto St^{\xi}, \end{equation}

with $0<\xi <1$, where the scaling exponent should be determined from concentration measurements in the near-bed region.

The following scaling of the pick-up flux is then expected:

(3.20)\begin{equation} E_p\propto \rho_p \left\langle \alpha \right\rangle u_*\propto{\rho_p}\,{St}^{\xi}\,u_*. \end{equation}

In the traditional non-dimensional form, this becomes

(3.21)\begin{equation} \phi\propto St^{\xi}\,\theta^{0.5}. \end{equation}

4. Results and discussion

In this section, we reconsider the experimental data introduced in § 1. Figure 6 shows the Obukhov length versus the Shields parameter. A minimum in the Obukhov length is observed for the CWF experiments conducted in the range $0.15<\theta <1$. This indicates that stratification effects, as discussed in § 2, could be present in this range, but buoyancy destruction is likely negligible for the other cases.

Figure 6. Obukhov length over channel height versus the Shields parameter. The dashed vertical line indicates the region ($\theta < 0.15$) where the pick-up probability of a single grain strongly influences the pick-up rate.

In the CWF experiments with a lower Shields number, $\theta <0.15$, particle contact load is the dominant transport form; see figure 7. For this reason, it is very unlikely that the pick-up flux for those cases is controlled by TKE budget considerations. It is more reasonable to expect that the pick-up flux in this regime is controlled by the probability that an individual particle is lifted from the bed by the turbulent flow that developed above the fixed bed section upstream of the erodible bed section in the CWF experiments.

Figure 7. Ratio of the terminal settling velocity and friction velocity plotted against the Shields parameter to identify the expected transport mode. Symbols are as defined in figure 2. The horizontal lines indicate the threshold values for the different transports (Dey 2014).

In the experiments of Bisschop (2018), where $\theta \gg 1$, the near-bed concentration was measured by conductivity probes with vertical separation 1 cm. The first measurement point at 1 cm height above the bed shows that $\left \langle \alpha \right \rangle >0.1$ for all cases; see figure 8. In this case, the Boussinesq assumption is no longer valid. As detailed in § 3, mechanisms other than stratification are probably controlling the scaling of the pick-up flux. This is consistent with the observed high values of the Obukhov length in figure 6. It is expected that the near-bed concentration depends only weakly on the Stokes number, which could be modelled by (3.19). Figure 8 demonstrates that this type of scaling is reasonably accurate, using a simplified Stokes number $St={t_p u_*}/{\kappa y}$, with a constant $t_p$ modelled by (2.1).

Figure 8. Near-bed concentration versus Stokes number (see (2.1)), measured at (a) $1$ cm and (b) $3$ cm height above the bed; data of Bisschop (2018). Symbols are as defined in figure 2. The dashed line represents a least squares regression for both heights separately, $\left \langle \alpha \right \rangle =0.41\,St^{0.20}$. The dash-dotted line is a least squares regression for the combined measurements at $1$ and $3$ cm height above the bed.

Now we will examine the scaling of the pick-up flux for both the CWF and HCF experiments. Regarding the CWF experiments, figure 9 shows the two scaling laws for the dimensionless pick-up flux (2.14) and (2.16), as derived from buoyancy arguments in § 2. It is observed that (2.14) gives a reasonably accurate fit to the experimental data. This indicates that stratification effects play an important role in the pick-up flux in the CWF regime at moderate Shields numbers. The assumption that the Osmidov length scale is halted by the particle size, which leads to (2.16), could be valid in some cases. The limited quality of the scaling in figure 9(b) suggests, however, that this is not generally true.

Figure 9. Measured pick-up flux versus computed pick-up flux using (a) (2.14) and (b) (2.16). The Shields parameter falls in the range $0.15<\theta <1.36$ for the selected data points. The outer dashed lines indicate a 30 % error margin.

For the HCF experiments, the scaling of near-bed concentration as observed in figure 8 can be used to predict the pick-up flux; see figure 10. Although there is some scatter in the data, the overall trend is quite well captured. It is also satisfying that the pick-up flux in the Winterwerp et al. (1992) data follows the same trend as obtained from the concentration and pick-up measurements of Bisschop (2018). These results support the scaling hypothesis that the pick-up flux at high Shields numbers, with associated high near-bed concentrations, is controlled by additional TKE modulation terms, instead of the standard buoyancy destruction term, as discussed in § 3.

Figure 10. Scaling of the pick-up flux at a high Shields number $\theta >8.7$ of Bisschop (2018) and data points obtained from Winterwerp et al. (1992). Pick-up flux (a) evaluated at a fixed reference height $y_r=3$ cm above the bed ($y_r/h=0.15$) as in Bisschop (2018), and (b) evaluated at a fixed non-dimensional reference height $St=1$ or $y_r=t_pu_*/\kappa$ and $\left \langle \alpha \right \rangle _r=0.41$. The dashed line indicates perfect correlation.

This does not mean that other explanations are not possible, such as the dilatancy hardening effect; see § 1. The pore pressure just below the bed and prior to pick-up is difficult to measure and therefore not available. However, using the estimate for the critical Shields parameter $\theta _{cr}$ corrected for the dilatancy effect by the theory of Van Rhee (2010), it is verified that $\theta _{cr}\ll \theta$. This suggest that dilatancy hardening could be present, but it is unlikely that it affected the pick-up rate in the experiments considered here. Another explanation proposed by van Rijn et al. (2019), as discussed briefly in § 1, is that the pick-up flux at high Shields numbers is reduced by a higher effective viscosity of the mixture of water and particles. In the context of this study, we need to substitute the effective mixture viscosity in (2.14). In summary, changes in the relative importance of terms in the TKE budget under high Shields conditions are a plausible explanation for the scaling anomaly observed in figure 2.

5. Conclusion

In this paper, the physical origin of several empirical scaling laws for the sediment pick-up flux has been investigated by considering the turbulent kinetic energy (TKE) balances for a stratified single-phase flow under the Boussinesq approximation and a full two-phase flow formulation. The scaling relations derived from the aforementioned balances are tested against existing erosion experiments.

At low Shields number, $\theta <0.15$, the probability of pick-up of a single particle by the overlying turbulent flow strongly depends on the Shields number. Buoyancy destruction of TKE is presumably not the dominant process in this regime.

In the range $0.15 < \theta \lesssim 1.5$, stratification effects likely dominate the sediment pick-up process. The non-dimensional pick-up flux $\phi$ can be scaled with the particle Reynolds number based on the shear velocity and the Shields number: $\phi \propto Re_p\,\theta ^{1.5}$.

For high Shields number flow, $\theta \gg 1$, a high near-bed concentration develops. Based on the observation of the Obukhov length in this regime, the buoyancy destruction of TKE is not the dominant process. It is anticipated that dissipation of TKE by collisions and TKE dissipation by drag become the dominant sinks and balance the shear production of TKE. From this hypothesis, a non-dimensional scaling relation is retrieved for the non-dimensional pick-up flux, $\phi \propto St^{\xi }\,\theta ^{0.5}$. From the observed sediment concentrations in the high Shields regime, a value $\xi =0.2$ is found. It is observed that the non-dimensional pick-up rate scales in accordance with this value, $\phi \propto St^{0.2}\,\theta ^{0.5}$. This scaling is confirmed by two independent data sets.

In order to further advance our understanding of erosion processes, new high-resolution measurements of near-bed sediment transport phenomena as well as interface-resolving DNS simulations of this problem are needed. These data could be used to formulate closures for both the instantaneous and Reynolds-averaged two-phase equations. The long-term goal would be to develop an accurate pick-up flux formula for large-scale morphodynamical models for riverine, coastal and offshore applications.