Articles | Volume 10, issue 4
Earth Surf. Dynam., 10, 817–832, 2022
Earth Surf. Dynam., 10, 817–832, 2022
Research article
03 Aug 2022
Research article | 03 Aug 2022

Probabilistic description of bedload fluxes from the aggregate dynamics of individual grains

Probabilistic description of bedload fluxes from the aggregate dynamics of individual grains
J. Kevin Pierce1, Marwan A. Hassan1, and Rui M. L. Ferreira2 J. Kevin Pierce et al.
  • 1Department of Geography, The University of British Columbia, Vancouver, Canada
  • 2Department of Civil Engineering, Instituto Superior Técnico, University of Lisbon, Lisbon, Portugal

Correspondence: J. Kevin Pierce (


We formulate the bedload sediment flux probability distribution from the Lagrangian dynamics of individual grains. Individual particles obey Langevin equations wherein the stochastic forces driving particle motions are switched on and off by particle entrainment and deposition. The flux is calculated as the rate of many such particles crossing a control surface within a specified observation time. Flux distributions inherit observation time dependence from the on–off motions of particles. At the longest observation times, distributions converge to sharp peaks around classically expected values, but at short times, fluctuations are erratic. We relate this scale dependence of bedload transport rates to the movement characteristics of individual sediment grains. This work provides a statistical mechanics description for the fluctuations and observation-scale dependence of sediment transport rates.

1 Introduction

Bedload transport refers to conditions when grains bounce and skid along the riverbed (Church2006). Numerous applications require predictions of bedload transport rates: ecological restoration, river engineering, and landscape evolution modeling provide examples. Predicting bedload fluxes is notoriously difficult, in part because they display wide fluctuations (Bunte and Abt2005; Recking et al.2012). Variability in transport rates is strongest in weak transport conditions, wherein instantaneous fluxes can deviate from mean values by several orders of magnitude (Ancey et al.2008; Benavides et al.2022). Transport fluctuations lend scale dependence to measured transport signals, whereby the statistical moments of transport depend on the averaging timescales (Singh et al.2009; Saletti et al.2015; Ancey2020). Because weak bedload transport conditions are common in gravel-bed rivers (Andrews1994; Oldmeadow and Church2006; Church and Hassan2002), predictions of bedload transport should include estimates of variability and the temporal scales over which averaged values converge.

Particle-based, stochastic approaches have been developed from which mean values, probability distributions, and the dependence of measured values on averaging scales can all be obtained (e.g., Ancey and Pascal2020; Turowski2010). These approaches provide a more complete description of bedload transport than deterministic models (e.g., Kalinske1947; Bagnold1966). However, existing stochastic models remain largely kinematic because they do not model the Newtonian forces generating particle transport (Furbish et al.2021). In this paper, we develop a new statistical mechanics formulation of the bedload flux which is based on a Newtonian model for the dynamics of individual grains. This effort provides new understanding of how transport rate fluctuations and scale dependence originate from the grain-scale dynamics.

The original stochastic description of bedload displacement is due to Einstein, who calculated particle trajectories as a random sequence of rests interrupted by instantaneous steps (Einstein1937), in what amounts to a pioneering application of the continuous-time random walk (e.g., Montroll1964; Weiss1994). Many subsequent works have since generalized Einstein's original model, either by introducing different distributions for resting times and step distances than Einstein originally considered (Hubbell and Sayre1964; Schumer et al.2009; Zhang et al.2012) or by incorporating resting of particles within the bed subsurface (Wu et al.2020; Pierce and Hassan2020b, a). A conceptually different modification of Einstein's model replaces his instantaneous steps with intervals of particle motion at a constant velocity (Lisle et al.1998; Lajeunesse et al.2017). This approach can be summarized with the stochastic equation

(1) x ˙ ( t ) = V σ ( t ) ,

where x(t) is the sediment position, V is the constant velocity of moving grains, and σ(t) is dichotomous noise which flips randomly between σ=0 (rest) and σ=1 (motion) (e.g., Bena2006). Owing to the constant velocity assumption, this model for sediment displacement applies only at long timescales when the details of individual particle movements become irrelevant (Weiss1994). Despite this timescale limitation, Eq. (1) establishes a starting point to describe bedload displacements with additional resolution compared to the random walk approach. To calculate displacements at short timescales, we should replace V with a time-dependent velocity derived from the forces generating particle motions.

The velocities of moving grains fluctuate due to hydrodynamic forcing and particle–bed collisions (Heyman et al.2016; Fathel et al.2015). Experiments indicate that downstream velocity distributions of particles can be exponential (Fathel et al.2015; Lajeunesse et al.2010), Gaussian (Heyman et al.2016; Ancey and Heyman2014), or gamma-like (Liu et al.2019; Houssais et al.2015). This range of observations has been attributed to differences in the balance of long and short movements within the particle velocity statistics (Wu et al.2020; Wu et al.2021) or to differences between experiments in the amount of momentum dissipated by particle–bed collisions (Pierce2021). Several studies have successfully modeled the downstream velocity fluctuations of moving particles by making an analogy to Brownian motion (Fan et al.2014; Ancey and Heyman2014). These studies represent exponential and Gaussian velocities with the Langevin equation

(2) v ˙ ( t ) = F ( v ) + ξ ( t ) ,

where v(t) is the downstream sediment velocity, F(v) is a deterministic force per unit mass, and ξ(t) is a random force per unit mass. ξ(t) has been modeled as Gaussian white noise for simplicity, although this is not an essential requirement. The choice F(v)=-μv/|v|+Δ produces exponential velocities, while F(v)=γ(V-v) produces Gaussian velocities. In these expressions, μ is a Coulomb friction-like coefficient, Δ is a constant force per unit mass, γ is a relaxation rate per unit mass, and V a steady-state velocity. Equation (2) provides a reasonable description of bedload velocities over short timescales, but it cannot describe particle displacements over longer timescales when entrainment and deposition also moderate the particle transport.

Motion–rest alternation and the fluctuating movement velocities of individual particles both lend variability to the sediment flux (Böhm et al.2004; Roseberry et al.2012). Sediment fluxes have been defined with both surface and volume definitions (Ballio et al.2014, 2018). The ensemble-averaged flux has been formulated using a surface definition by Furbish et al. (2012, 2017):

(3) q ( x , t ) = 0 d x 0 d t R ( x , t ) E ( x - x , t - t ) .

This “nonlocal formulation” generalizes earlier approaches based on aggregating particles from upstream locations (Nakagawa and Tsujimoto1976; Parker et al.2000). In this equation, E(x-x,t-t) is the average entrainment rate a distance x upstream and a time t in the past, while R(x,t) is the probability that a just-entrained particle displaces at least a distance x in time t. For the steady case E(x,t)=E wherein particles step an average distance  between entrainment and deposition, the nonlocal flux becomes q〉=E, in accord with Einstein (1950). The nonlocal formulation has not yet been extended to characterize the sediment flux as a stochastic process.

The sediment flux has been described as a stochastic process using both renewal theory and population modeling approaches. These methods rely on additional Eulerian characteristics of particle transport, which apply to an ensemble of particles occupying a volume or crossing a surface. Renewal models introduce inter-arrival time distributions ψ(t), which characterize the time intervals between successive arrivals of particles to a surface (Turowski2010; Ancey and Pascal2020). The flux is calculated as the rate of particle crossings over an observation time T:

(4) q ( T ) = N ( T ) T .

Here, 𝒩(T) is the number of particle crossings by T – related to ψ(t). For exponential inter-arrival times ψ(t)=λe-λt, the flux becomes Poissonian with rate constant λ, so the mean flux is q〉=λ. For other ψ(t) values, the flux depends on the observation time T in a type of scale dependence. Instead of counting particle arrivals, population modeling approaches introduce Eulerian entrainment and deposition rates to count the number of moving particles in a volume (Ancey et al.2006, 2008). These approaches compute the flux by summing the velocities of all moving particles in a volume (Heyman2014; Ancey and Pascal2020). Both of these stochastic approaches provide excellent correspondence with experimental data. However, we still have little understanding of how to relate their Eulerian input parameters to the underlying grain-scale mechanics (Heyman et al.2016). This limitation exposes a need for further research into how the trajectories of individual grains ultimately produce a stochastic sediment flux and control its dependence on observation scale.

In this paper, our objective is to formulate the stochastic sediment flux directly from the Lagrangian dynamics of individual grains rather than by introducing additional Eulerian quantities, such as volumetric entrainment and deposition rates or inter-arrival time distributions. To achieve this, in Sect. 2 we extend the motion–rest alternation model for particle displacement in Eq. (1) to include Newtonian velocities from Eq. (2). This produces a mechanistic–stochastic model of particle displacement which is valid across a wider range of timescales than earlier descriptions. We then construct the sediment flux in Sect. 3 by accumulating the displacements of individual particles. The resulting formulation shares elements from both the nonlocal and renewal approaches summarized in Eqs. (3) and (4). In Sect. 4 we solve the formulation to derive the probability distributions of particle position and the sediment flux, and in Sects. 5 and 6 we discuss the new features and limitations of our approach, summarize its relationship to earlier work, and suggest several directions for further development.

2 Stochastic description of bedload transport

The starting point for our analysis is an idealized one-dimensional domain populated with sediment particles on the surface of a sedimentary bed. Particles are set in motion by the flow and move downstream until they deposit, and the cycle repeats. The downstream coordinate is x so that x˙=dx/dt describes a velocity in the downstream direction. The flow is considered weak enough that interactions among moving grains are uncommon, although interactions between moving particles and the bed occur regularly. These conditions are characteristic of “rarefied” bedload transport conditions (e.g., Kumaran2006; Furbish et al.2017). Particles are considered to have similar enough shapes and sizes so as to have nearly identical mobility characteristics. These conditions allow all particles to be described as independent from one another but governed by the same underlying dynamical equations. Any of these conditions could be relaxed in exchange for additional mathematical difficulty.

2.1 Mechanistic formulation of intermittent transport

From these assumptions, we propose an equation of motion for the individual sediment grain including two features. First, particles should alternate between motion and rest, similar to the earlier motion–rest models summarized by Eq. (1). Second, the velocities of moving particles should evolve according to the Newtonian equation (Eq. 2). These dynamics can be represented as

(5) x ˙ ( t ) = v ( t ) σ ( t ) , v ˙ ( t ) = [ F ( v ) + ξ ( t ) ] σ ( t ) .

In these equations, F(v) is a deterministic forcing term whose structure can be chosen to produce the desired velocity distribution for moving particles (exponential, Gaussian, or others), ξ(t) is Gaussian white noise with a correlation function

(6) ξ ( t ) ξ ( t + τ ) = 2 Γ δ ( τ )

that represents the velocity fluctuations of moving particles, and Γ is a velocity diffusivity [L2 T−3] which controls the intensity of velocity fluctuations. Figure 1 displays particle velocities and displacements derived from these equations for the choice of Gaussian movement velocities.

Figure 1Panel (a) shows the velocity from Eq. (5) for a particular realization of the noises ξ(t) and σ(t), while panel (b) shows the position derived in Eq. (5) alongside other possible trajectories. Keys in panel (a) indicate the average movement time 1/kD, rest time 1/kE, and motion–rest alternation process of σ(t), while keys in panel (b) show the average movement velocity V and evolution of the displacement for different values of σ(t). Motion–rest alternation with velocity fluctuations produces tilted stairstep trajectories with unsteady slopes in the xt plane.


These equations represent Newtonian dynamics that are randomly paused as σ(t) alternates. The quantity σ(t) is dichotomous Markov noise (e.g., Horsthemke and Lefever1984; Bena2006) which produces alternation between motion (x˙ generally nonzero) and rest (x˙ zero). This noise takes on values σ=1 (motion) and σ=0 (rest). The entrainment rate (σ=0σ=1) is labeled kE, and the deposition rate (σ=1σ=0) is labeled kD. Times spent in motion and rest are respectively distributed as P(t)=kDexp(-kDt) and P(t)=kEexp(-kEt), so the mean movement time is kD-1, and the mean resting time is kE-1. The notation k=kE+kD is a shorthand used throughout the paper. 1/k represents the correlation time of the dichotomous noise.

The transport process described by Eq. (5) is intermittent because the particle velocity x˙ and acceleration v˙ are randomly switched on and off by entrainment and deposition. Notably, the presence of σ(t) in these two equations decouples x˙ and v. In general, x˙v, in contrast to conventional mechanics. This decoupling means the relevant velocity scale for particle displacement is the “virtual velocity” x˙. This velocity is “virtual” because it contains the particle velocity during motions in addition to the intervening rests (Hassan et al.1991; Bradley and Tucker2013). The second velocity scale v represents the velocities of particles if they are moving. During motions, v evolves as a stochastic process. This evolution is paused during particle rests, when deposition shuts off the driving forces in Eq. (5).

The formulation of Eq. (5) is conceptually similar to the work of Fan et al. (2016), which simulates particle transport by manually switching Eq. (2) on and off to represent entrainment and deposition. We have replaced this manual switching with dichotomous noise to obtain an analytical representation of particle transport by motion–rest alternation. This description is similar to several intermittent transport models in the stochastic physics literature, as it involves both (1) interrupted diffusion (e.g., Laskin1989; Łuczka et al.1993; Balakrishnan et al.2001) and (2) dichotomous noise entering at second order in time (e.g., Masoliver1992, 1993). To our knowledge, Eq. (5) is the first stochastic transport model to include both of these components simultaneously.

2.2 Phase space description of motion–rest alternation

The time evolution of Eq. (5) for a particular realization of ξ(t) and σ(t) maps a trajectory through the phase space spanned by x and v. The conditional probability density W(x,v,t|x0,v0) represents the likelihood that a phase trajectory reaches (x,v) at time t provided it passed through (x0,v0) at t=0. This density characterizes the stochastic evolution of the particle position and velocity.

A master equation for the phase space density can be formulated by noting that the combined process (x,v,σ) is Markovian (Horsthemke and Lefever1984). In Appendix A, we demonstrate that Eq. (5) implies the master equation

(7) t t + k W ( x , v , t | 0 ) = t + k E L ^ K W ( x , v , t | 0 ) ,

using the abbreviation W(x,v,t|x0,v0)=W(x,v,t|0). In this equation, L^K is the Kramers operator, familiar from the description of Brownian particles in an external force field (e.g., Risken1984; Kubo et al.1978):

(8) L ^ K = - v x + v - F ( v ) + Γ v .

The Kramers operator evolves the particle position and velocity. The term with x represents drift, while terms containing v represent accelerations, with advective contributions from F(v) and diffusive ones from Γ (Duderstadt and Martin1979).

To understand the structure of Eq. (7), it is useful to recall the classical Kramers equation for particles driven by a fluctuating force without intermittency: tW(x,v,t|0)=L^KW(x,v,t|0), with the same operator L^K as above (Kramers1940). A comparison suggests that the additional terms in Eq. (7) weave intermittency into the distribution function. In fact, there are two limiting behaviors contained in Eq. (7) at opposite extremes of t (Bena2006). At short times, higher orders of t dominate, giving the classical Kramers equation tWL^KW after an integration over time. At long times, the lower-order terms in t dominate instead, giving tW(kE/k)L^KW. This is still a Kramers-type equation, although the evolution of x and v is slowed by the “intermittency factor”, kE/k. This factor represents the average fraction of time particles have spent in motion. Thus, the higher-order time derivatives in Eq. (7) encode a transition between conventional and intermittent regimes of motion. This structure reflects the increasing influence of entrainment and deposition on the evolution of x and v with time.

2.3 The displacement statistics of bedload grains

To calculate the sediment transport rate we will use the probability distribution of position for bedload particles, defined as

(9) P ( x , t | 0 ) = - d v W ( x , v , t | 0 ) .

Unfortunately, even for the classical Kramers equation without intermittent motion, it is extremely difficult to obtain a governing equation for P(x,t|0) without first solving the phase space master equation for W(x,v,t|0) (e.g., Brinkman1956; Olivares-Robles and García-Colin1996). For the case of F(v) associated with exponential velocities (e.g., Fan et al.2014), the Kramers equation has only been solved numerically (Menzel and Goldenfeld2011).

Fortunately, bedload experiments often show Gaussian velocities for moving particles (e.g., Martin et al.2012; Ancey and Heyman2014; Heyman et al.2016), corresponding to the forcing term

(10) F ( v ) = γ ( V - v )

in Eq. (7). With this force, the classical (non-intermittent) Kramers equation analogous to Eq. (7) can be solved exactly for the position distribution P(x,t|0) (e.g., Wang and Uhlenbeck1945; Chandrasekhar1943).

When motions are intermittent as in Eq. (7), an analytical solution appears to not yet be possible. However, we can still solve Eq. (7) with the force (Eq. 10) approximately by using the same “overdamped” approximation often applied to the classical Kramers equation (e.g., Risken1984; Gardiner1983). This approximation is possible whenever a term like Eq. (10), linear in velocity, is present in a governing Langevin equation. In our context, the overdamped approximation applies whenever moving particles attain their steady-state velocities relatively quickly after entrainment. Campagnol et al. (2015) suggest that the “acceleration phase” following entrainment will only affect a small portion of the trajectory between motion and deposition when the flow is sufficiently strong. In addition, earlier Langevin models have successfully described the velocity distributions of bedload particles by neglecting the transient acceleration phase following entrainment (Fan et al.2014; Ancey and Heyman2014). For these reasons, we expect that the overdamped approximation is a reasonable first approach for solving Eq. (7).

We can actually obtain the overdamped approximation for the phase space equation (Eq. 7) using the same method originally introduced by Kramers (1940).

We integrate Eq. (7) along the straight line x+vγ-1=const. from v- to v→∞. Because γ−1 is small for fast relaxation, we can take the line integral along dv only (Coffey et al.2004), providing the overdamped master equation:

(11) t 2 + k t + V x t + k E V x - D x 2 t - k E D x 2 P ( x , t ) = 0 .

The spatial diffusivity D is defined as D=γ-2Γ with units [L2 T−1]. Hereafter we suppress the explicit dependence of the position distribution on its initial conditions [P(x,t|x0,v0,t0)=P(x,t)].

Equation (11) interleaves two different diffusion processes: one associated with motion–rest alternation and another with velocity fluctuations during motions. Terms involving V represent advection, while those involving D represent diffusion from velocity fluctuations. Terms involving kE and k represent diffusion due to motion–rest alternation. Mixed orders of t encode the aforementioned transition between conventional and intermittent transport.

3 Mechanistic formulation of the sediment flux

To phrase the probability distribution of the sediment flux in terms of these particle dynamics, we apply a method very similar to the one developed by Banerjee et al. (2020) to describe currents of active particles in condensed matter physics. The basic idea, as depicted in Fig. 2, is to initially distribute N particles in all states of motion along the domain -Lx0. Later, the number of particles N and the size L of the domain will be extended to infinity such that their ratio ρ=N/L remains constant. This limit produces a configuration similar to the one considered in the nonlocal formation of Eq. (3).

Figure 2Panel (a) indicates the configuration for the flux. Here, time increases from the bottom to the top. Particles begin their transport with positions -Lx0 to the left of S at observation time T=0, and the flux is calculated in panel (b) as the rate N(T)/T of particles crossing x=0 over the observation time T. Particle crossing events are indicated in (b) by color-coded lines. The probability distribution of q(T) is determined from all possible realizations of the trajectories and initial positions as N and L tend to infinity, while the density of particles ρ=N/L to the left of S is held constant.


From this initial configuration, the flux is calculated as the average rate of particles crossing to the right of the control surface at x=0 after the sampling time T, analogous to the renewal formulation of Eq. (4):

(12) q ( T ) = 1 T i = 1 N I i ( T ) .

In this equation, Ii(T) is an indicator function which equals 1 if the ith particle has crossed the control surface (x=0) by T and 0 otherwise. Particles which have not crossed the surface (or which have crossed and then crossed back) do not contribute to the flux.

The probability density F(q|T) of the flux, conditional on the sampling time T, is then an average involving Eq. (12) across all possible initial configurations of particles and their trajectories:

(13) F ( q | T ) = δ q - 1 T i = 1 N I i ( T ) .

Taking the Laplace transform F̃(s|T)=0dqe-sqF(q|T) (forming the characteristic function) produces

(14) F ̃ ( s | T ) = i = 1 N 1 - 1 - e - s / T I i ( T ) .

This formula relies on the independence of averages for each particle (so that the average of a product is the product of averages) and the observation that eαI=1-(1-eα)I if I=0,1.

The average over initial conditions and possible trajectories of the indicator for the ith particle involved in this characteristic function is

(15) I i ( T ) = 1 L - L 0 d x 0 d x P ( x - x , T ) .

P(x,T) is the probability distribution of position at time T, either derived exactly from Eqs. (7) and (9) or from the overdamped approximation Eq. (11). The integral over x evaluates the probability that the position of a particle at T exceeds x=0, while the integral over x averages across a uniform distribution of possible initial positions. These Ii(T)〉 terms are the components of the flux that depend on the particle dynamics.

Inserting Eq. (15) into Eq. (14) and taking the limits L→∞ and N→∞, as the density of particles ρ=N/L is held constant, provides

(16) F ̃ ( s | T ) = exp - 1 - e - s / T Λ ( T ) .

Λ(T) is the central parameter of the sediment flux probability distribution:

(17) Λ ( T ) = ρ 0 d x 0 d x P ( x + x , T ) .

The quantity Λ(T)/T is a rate function. This ratio describes the rate of particle arrivals to the control surface at x=0, and it explicitly depends on the observation time T.

Equation (16) is the characteristic function of a Poisson distribution (Cox and Miller1965). Expanding in e-s/T and inverting the Laplace transforms provides the probability distribution of the flux conditional on the sampling time T:

(18) F ( q | T ) = n = 0 Λ ( T ) n n ! e - Λ ( T ) δ q - n T .

This equation implies that the mean flux is q(T)=0qF(q|T)dq=Λ(T)/T. Similarly, the variance is σq2(T)=Λ(T)/T2. For the case when Λ(T) is proportional to the observation time (Λ∝T), these formulas become identical to the renewal approach with exponential inter-arrival times.

Equation (18) formulates the flux probability distribution directly from the particle dynamics set out in Eq. (5). This equation is a scale-dependent Poisson distribution. The Poisson form originates primarily from the assumption that particles undergo independent dynamics. The distribution is scale-dependent through the displacement statistics of individual particles ingrained in Eq. (17).

4 Results

4.1 Displacement by intermittent transport

The overdamped master equation (Eq. 11) describes the displacement statistics of particles alternating between motion and rest with Gaussian movement velocities. This equation is founded on the approximation that just-entrained particles attain their steady-state (but fluctuating) velocities rapidly.

The overdamped master equation (Eq. 11) is solved in Appendix B with transform calculus, obtaining

(19) P ( x , t ) = A ^ 0 t d t e - k E ( t - t ) - k D t × I 0 2 k E k D t ( t - t ) G ( x - V t , t ) .

Here 0 is a modified Bessel function, A^=-Dx2+Vφx+k+t is a differential operator, and G(x,t)=exp[-x2/4Dt]/4πDt is a Gaussian propagator. Within the integral, the Bessel term represents the proportion of time u/t the particle has spent in motion, while the Gaussian term describes the distribution of displacements achieved in time u. This distribution is compared with numerical simulations of the exact distribution from Eq. (5) in Fig. 3a. The general decreasing trend of mean transport with observation time is qualitatively consistent with laboratory observations (Singh et al.2009; Saletti et al.2015) and the renewal approach summarized earlier (Turowski2010; Ancey and Pascal2020).

Figure 3Panel (a) indicates the overdamped probability distribution of particle position (Eq. 19) as it evolves through time, while panel (b) shows the resulting particle diffusion from Eq. (20). All results are scaled by the mean hop length =V/kD and the timescale τ=1/k of the motion–rest alternation. Curves represent the analytical results, while the points are results of Monte Carlo simulations of the exact Eq. (5) produced by evaluating the cumulative transition probabilities on a small time step (e.g., Barik et al.2006). In panel (a) from the initial mixture of motion and rest states, particles advect downstream as they diffuse apart from one another due to motion–rest alternation and velocity variations. The initial distribution persists as a delta function spike for the t=10τ curve. In panel (b) at timescales t1/k the diffusion is approximately ballistic since particles have not had time to exchange between motion and rest. For t1/k particles undergo normal diffusion as particles become well-mixed among motion and rest states. Small discrepancies are visible at short times between the simulations and analytical approximations due to our neglect of velocity fluctuations in deriving Eq. (20).


The moments of particle position from Eq. (5) are challenging to obtain. An approximation for the mean can be obtained by considering that velocity fluctuations during motions approximately cancel out, since these are symmetrical around V. Therefore, we set Γ=0 in Eq. (5), multiply by x, and integrate to find the mean position x(t)=kEVt/k, which is Vt scaled by the expected fraction of time spent in motion. Similarly, we can approximate the variance by reasoning that motion–rest alternation, and not velocity fluctuations during motions, is the primary source of particle diffusion. Again setting Γ=0 we find the variance of position (σx2=x2-x2):

(20) σ x ( t ) 2 = 2 k E k D V 2 k 3 t + 1 k e - k t - 1 k .

This describes a two-range particle diffusion process, whereby the rate of particle spreading depends on how long the dynamics have been ongoing (Taylor1920). Figure 3b compares this variance to the numerical simulations. The variance exhibits a crossover between ballistic and normal scaling at τ=1/k. The approximate variance calculation is good apart from small undershooting at short times. Here, small contributions to the variance from velocity fluctuations in the motion state become visible, consistent with the expectation that particle velocity fluctuations during motions will enhance the particle diffusion at short timescales.

4.2 The flux rate function for overdamped transport

The formalism in Sect. 3 provides the central parameter Λ(T) of the sediment flux distribution. Using the overdamped probability distribution of position in Eq. (19), we evaluate Eq. (17) in Appendix C, providing

(21) Λ ( T ) = ρ 0 T d T I 0 2 k E k D T ( T - T ) e - k E ( T - T ) - k D T × D π T T + k T - k D 2 k e - V 2 T / 4 D + V 2 T + k T - k D k erfc - V 2 T 4 D .

In this equation, the notation T means that the partial time derivative acts from the left of all terms in which it is involved, as in f(T)Tg(T)=T[f(T)g(T)], and erfc(⋅) denotes the complementary error function.

This result indicates a nuanced observation-scale dependence in the sediment flux. We can better understand Eq. (21) by investigating extreme cases of the observation time. As shown in Appendix D, Eq. (21) takes on simple forms at extreme values of T:

(22) Λ ( T ) = ρ k E k D T π , T k D P e - 1 , ρ k E V T k , T k D P e - 1 .

Here, the quantity Pe=V2/(2DkD) is a Péclet number which measures the relative importance of advection and diffusion during particle motions. After entrainment, particles typically move for a duration 1/kD before deposition, over which they will displace LA=V/kD by advection and LD=±2D/kD by diffusion. Pe is composed as a squared ratio of these length scales: Pe=(LA/LD)2.

Figure 4Panel (a) plots the mean sediment flux for different values of the Péclet number Pe=V2/(2kDD), characterizing the relative significance of particle velocity fluctuations during motions. Plotted lines show the analytical result (Eq. 21), while the points are Monte Carlo simulations. Panel (b) displays the evolution of the flux distribution (Eq. 18) across observation times. The flux is normalized by the prediction q0=E of the Einstein model. In panel (a), the convergence time of the mean flux is controlled by attributes of individual particle motions. In all cases, the mean flux converges for T1/(kDPe), as expected by Eq. (22). Discrepancies between numerical simulations and analytical calculations emerge for timescales Tγ-1, indicating that the overdamped approximation cannot account for this “acceleration phase” period. A plateau emerges at intermediate times when velocity fluctuations are strong (Pe≪1). This plateau is likely associated with particles whose displacements lag behind due to diffusion “catching up” to other crossing grains. In panel (b), the Einstein limit F(q|T)=δ(q-E) is approached as the observation time T grows. Stronger velocity fluctuations (smaller Pe) slow the convergence. In all plots, Pe is modified by changing D, while all other parameters are held constant.


The limiting form of Λ(T) implies that for T(kDPe)-1 the mean flux converges to the eventual value q(T)=q0:

(23) q 0 = ρ k E V / k .

This can be understood as the result q0=E of the nonlocal formulation (Eq. 3) in the case of steady transport conditions. Thus, at T→∞, our formulation becomes equivalent to that of Einstein (1950). Here, E=ρkE is an averaged areal entrainment rate and =V/kV/kD is the mean step length of particles. kkD holds since the mean duration of a single motion (1/kD) is typically much smaller than the duration of a single rest (1/kE).

Figure 4a shows the rate constant decaying toward its asymptotic value in Eq. (22) for different values of Pe, with all other parameters fixed. Numerical simulations of the exact equations (Eq. 5) are superimposed. The overdamped approximation pursued in this paper provides a valid characterization of the sediment flux for T1/γ, but for T1/γ, the approximate result overshoots. Thus, we expect that the particle acceleration phase immediately after entrainment (Campagnol et al.2015) slows the observation-scale dependence of the flux at short timescales.

Figure 4b demonstrates the adjustment of the flux distribution (Eq. 12) across observation times. At the shortest times, the flux approaches a uniform-like distribution due to the (apparent) divergence of Λ(T). At very long observation times, the flux adopts the deterministic (Einstein) form

(24) F ( q | T ) δ ( q - E ) .

This limit can be seen by taking large T in Eq. (18) and using the correspondence between Poisson and Gaussian distributions for large Poisson rates (e.g., Cox and Miller1965). The Einstein (1950) result for the sediment transport rate becomes exact only when kDTPe-1. Otherwise, the flux has intrinsic fluctuations related to the unpredictability of particle arrivals in a finite observation window, as characterized by Eq. (12).

5 Discussion

In this paper, we have formulated a mechanistic description of the bedload sediment flux using a detailed stochastic model of individual particle displacements. The resulting sediment flux distribution shows Poissonian fluctuations that depend on observation scale. Our displacement model applies over a wider range of timescales than earlier formulations because it includes both Newtonian velocities and motion–rest alternation. In appropriate simplified limits, the displacement model Eq. (7) reduces to many earlier descriptions of grain-scale transport (e.g., Einstein1937; Lisle et al.1998; Lajeunesse et al.2017; Ancey and Heyman2014; Fan et al.2014). In this sense, the grain-scale transport component of this paper generalizes and unifies all of these earlier works.

We solved the displacement model analytically to obtain the displacement probability distribution. This derivation relied on the “overdamped” approximation that particles accelerate rapidly following entrainment. This approximation is only possible when the velocities of moving grains are Gaussian. We then formulated the stochastic sediment flux using the resulting particle displacement statistics. The obtained flux distribution mimics earlier renewal theory descriptions of the bedload flux (Turowski2010; Ancey and Pascal2020), except its input parameters relate directly to the Lagrangian transport characteristics of individual grains. The scale dependence of the flux is mainly controlled by the Péclet number and deposition rate of grains, indicating that the scale dependence originates from the velocity fluctuations and motion–rest alternations of individual grains. The flux is enhanced at short observation times because it is dominated by uniquely fast-moving particles. As time grows, a plateau may emerge if the diffusion is sufficiently strong. This plateau can likely be attributed to slow-moving particles “catching up” with grains which were faster to cross the control surface, slowing the flux decay. When the time spent in motion goes to zero (idealized steps) or velocity fluctuations vanish (Pe→∞), the flux loses its scale dependence. In other conditions, as the observation time becomes large, the flux approaches the deterministic result q=E of Einstein (1950) and later nonlocal formulations (Furbish et al.2012, 2017).

5.1 Newtonian description of bedload displacement

Our description of individual particle displacements in Sect. 2 provides an analytically solvable alternative to computational physics models of grain-scale transport (e.g., Schmeeckle2014; Ji et al.2014; Clark et al.2017). Analytical progress was possible largely because the forces on moving particles were modeled as linear in the particle velocity (Eq. 10) with Gaussian white noise fluctuations (Ancey and Heyman2014). This forcing structure is a crude approximation for the actual hydrodynamic and granular forces acting on bedload particles. In reality, flow forces are nonlinear in the flow velocity and history-dependent (Michaelides1997; Schmeeckle et al.2007), while collision forces are velocity-dependent and episodic (Brach1989; Schmeeckle and Nelson2003; Pierce2021). Fluctuations in turbulent flow forces are colored noise, not white (Schmeeckle et al.2007; Dwivedi et al.2010; Celik et al.2014). Although we cannot hope to include completely realistic forces in an analytically solvable model, it should still be possible to introduce some level of additional complexity into the forcing structure of Eq. (5). A first step is to solve Eq. (7) with F(v) tailored to produce exponential velocities for moving particles (e.g., Fan et al.2014). The result would lend additional insight into how the scale dependence of the flux depends on the particle velocity statistics.

A simplified element of our approach is the representation of entrainment and deposition by instantaneous alternation between motion and rest states (e.g., Einstein1937). In actuality, motion and rest are not perfectly defined because “resting” particles undergo slow downstream creep (Houssais et al.2015; Allen and Kudrolli2018) and can shuttle in their pockets without any net translation (Diplas et al.2008; Celik et al.2010). Our dichotomous representation of entrainment and deposition provides no direct information on the timescales of particle motion and rest, and relating the motion and rest timescales of grains to the underlying fluid and granular mechanics remains an important open problem. It is challenging to imagine an analytically solvable model of particle displacement which does not discriminate between multiple states of transport, but improvements can nonetheless be made to the dichotomous representation of entrainment and deposition we employed. One option is to render slower-moving particles more likely to deposit, which can be viewed as making the dichotomous noise “state-dependent” (e.g., Laio et al.2008; Bartlett and Porporato2018). This would preferentially cull slow-moving particles and positively skew the particle velocity distribution (Williams and Furbish2021), probably affecting the sediment flux.

5.2 Mechanistic interpretation of transport fluctuations

The sediment flux probability distribution derived in Sect. 3 represents a Poisson distribution with a scale-dependent rate. Poisson distributions have relatively thin tails which reflect narrow sediment transport fluctuations. In reality, sediment flux distributions are only Poissonian at high transport rates, whereas in other conditions they have wide tails representing the possibility of large fluctuations (Ancey et al.2008; Saletti et al.2015; Turowski2010), which appear as bursts (e.g., Goh and Barabási2008) in sediment flux time series (Dhont and Ancey2018; Singh et al.2009). This observation highlights a need to improve the mechanistic description we developed here to produce wider transport rate fluctuations.

A vast set of processes generates transport rate fluctuations in real channels. At the shortest timescales, fluctuations arise from the intermittent arrivals of individual grains, as we have described here. Over longer timescales, activity waves (Heyman et al.2014), cluster entrainment (Papanicolaou et al.2018; Strom et al.2004), bedform migration (Guala et al.2014; Hamamori1962), grain-size sorting (Cudden and Hoey2003; Iseya and Ikeda1987), flow variations (Mao2012; Wong and Parker2006), and sediment supply perturbations (Lisle et al.1993; Madej et al.2009) all contribute to sediment transport variability. It should be possible to include particle–particle interactions in our description to capture the subset of these processes which originate at the grain scale. Activity waves, clusters, and bedforms might result from including interactions between particles in the entrainment and deposition rates, such as collisions (Lee and Jerolmack2018), the stabilization of bed particles by neighbors (Church et al.1998), or coordinated deposition based on the locations of sedimentary deposits (McDowell and Hassan2020). We might formulate the resulting joint distribution of particle positions and velocities by analogy to reaction–diffusion problems (e.g., Pechenik and Levine1999; Cardy2008) or other interacting particle systems available in the physics literature (e.g., Escaff et al.2018; Hernández-García and López2004).

5.3 Observation-scale dependence of the flux distribution

Phrasing the transport rate in terms of the Lagrangian dynamics of individual grains produces a flux distribution that adjusts with the observation time. According to Eq. (22), this adjustment is largely controlled by the deposition rate and Péclet number Pe. The deposition rate controls the duration of motions, while the Péclet number characterizes the significance of velocity fluctuations. Typical values for the Péclet number can be estimated from experimental publications as the ratio of mean and root mean square streamwise velocities of grains, giving Pe=0.8 (Fathel et al.2015), Pe=1.4 (Heyman et al.2016), and Pe=3.7 (Martin et al.2012). For this range of Péclet numbers, Fig. 4 leads us to expect a monotonic decay of the mean flux with no obvious plateau at moderate observation times. Such a decay has been observed in several studies. Both Bunte and Abt (2005) and Singh et al. (2009) noted a power-law decay to a constant flux with observation time in field and laboratory data, respectively. However, these studies also observed that at high transport rates, the flux increased with observation time rather than decreased. Possibly, at higher transport rates, collective rather than grain-scale processes may principally control the scaling of the mean flux with observation time.

Both Singh et al. (2009) and Saletti et al. (2015) investigated higher statistical moments of the bedload flux beyond the mean in laboratory data. They found that higher moments also shift with the observation time T, but Singh et al. (2009) identified statistical multiscaling, wherein the flux distribution changes shape with T, while Saletti et al. (2015) identified monoscaling, wherein the distribution does not change shape with T. Our analysis predicts monoscaling because the flux distribution (Eq. 18) is always Poissonian, even though all of its moments scale together with T in a nontrivial way (see Fig. 4). Probably, the Poissonian and monoscaling characteristics of the flux distribution both originate from the assumed independence of individual particle motions, not from the distributions of velocities, movement times, or resting times. In particular, Eq. (18) indicates that changing the force F(v) to produce exponential velocities will only modify Λ(T), not the shape or scaling type (mono, multi) of the flux distribution. Possibly, wider-tailed, multiscaling sediment flux distributions will derive from generalizations of our approach to include particle–particle interactions. Turowski (2010) demonstrated that renewal theories with certain non-exponential inter-arrival times produce multiscaling, and it is likely that non-exponential distributions will originate from particle–particle interactions (Ancey et al.2008). To some extent, a flux distribution which is wide-tailed at short observation times must be multiscaling, since it should approach the thin-tailed Einstein distribution (Eq. 24) as the observation time becomes large, changing shape as it adjusts.

6 Conclusions

We have formulated the bedload flux probability distribution from the statistical mechanics of individual grains in transport. This formulation produces Poissonian flux distributions having scale-dependent rates, meaning transport rate fluctuations are relatively narrow, and transport characteristics shift with the timescales over which they are observed. In laboratory experiments, sediment transport fluctuations are typically wider than Poissonian. Notably, we can assert that the Poisson flux distribution derived in this paper originates exclusively from the independence of individual grains: the Poisson form is completely indifferent to the forces driving particles downstream so long as these forces do not introduce correlations between particles. In the future, it will be necessary to refine the statistical mechanics formulation presented here to produce wider transport fluctuations. We expect that introducing any component in Eq. (5) which couples one particle to another will achieve flux distributions wider than Poissonian. The severe challenge will be evaluating the average in Eq. (13) when grains are not independent.

Appendix A: Derivation of the phase space master equation

Because the joint process α=(x,v,σ) is Markovian, its phase space distribution function for a particular realization of the white noise ξ(t) obeys the Chapman–Kolmogorov equation (Cox and Miller1965; Van Kampen2007):

(A1) W ξ α , t + δ t | α 0 = d α W ξ ( α , t + δ t | α , t ) W ξ α , t | α 0 .

Here, dα=σ=0,1-dx-dv. The Chapman–Kolmogorov equation relates the phase space distribution function at t+δt to its value at t through the transition amplitudes Wξ(x,v,σ,t+δt|x,v,σ,t). The distribution Wξ is a functional of the white noise ξ(t) (Hänggi1985; Łuczka2005).

In the limit of vanishing δt, the transition amplitudes in Eq. (A1) can be directly evaluated from the dynamical equations (Eq. 5) using a method analogous to Horsthemke and Lefever (1984). The transition rates involve δ-function terms in x and v. These terms are expanded in δt to first order, and the integrals in Eq. (A1) are conducted over the δ functions. This produces the pair of equations

(A2) t M ξ = k E R ξ - k D M ξ + - x v + v [ - F ( v ) + ξ ( t ) ] M ξ , t R ξ = k D M ξ - k E R ξ .

In these equations, the shorthand expressions are Mξ=Wξ(x,v,1,t|x0,v0,σ0,t0) and Rξ=Wξ(x,v,0,t|x0,v0,σ0,t0). We now average Eq. (A2) over realizations of the white noise and compute the correlator ξ(t)Mξ=ΓvM using the Furutsu–Novikov theorem (e.g., Van Kampen2007; Gitterman2005). Incorporating this correlator into Eq. (A2) obtains

(A3) t M = k E R - k D M + L ^ K M , t R = k D M - k E R ,

where L^K=-xv+v{-F(v)+Γv} is the Kramers operator, M=〈Mξ, and R=〈Rξ. Summing the above equations, noting W(x,v,t|0)=M+R, and eliminating M from the sum finally produces Eq. (7).

Appendix B: Solution for the position probability distribution

The position probability distribution can be obtained from Eq. (11) provided we have a pair of initial conditions on P. We consider particles to have a probability kD/k=φ to start from rest, so they have a probability 1-φ=kE/k to start from motion. Particles are initially located at x=0, and particles that start from motion are considered to have a random initial velocity selected from the steady-state distribution

(B1) f ( v ) = γ 2 π Γ e - γ v 2 / 2 Γ .

With these assumptions, the initial state can be written M(x,v,0)=(1-φ)δ(x)f(v) and R(x,v,0)=φδ(x)δ(v). Summing these two equations and integrating out the velocity provides P(x,0)=δ(x). Plugging these two equations into Eq. (A3), summing the result, then integrating out the velocity provides tP(x,0)=-kEVkδ(x). This produces the required pair of initial conditions. A similar calculation is available in Weiss (2002).

Now we take Fourier transforms over space and Laplace transforms over time of the overdamped master equation (Eq. 11), obtaining

(B2) P ̃ ( g , s ) = s + k + D g 2 - i g V φ s ( s + k ) + D g 2 - i g V s + k E .

The Fourier transform can be inverted by partial fraction decomposition and contour integration to obtain

(B3) P ̃ ( x , s ) = - φ D x 2 + V φ x + s + k V Q s + k E exp V x 2 D - V | x | 2 D Q ,


(B4) Q = 1 + 4 D V 2 s ( s + k ) s + k E .

The Laplace transform can then be inverted with the shift property L-1{f̃(s+k)}=e-ktf(t), the derivative property (Arfken1985)

(B5) L - 1 s f ̃ = δ ( t ) + t f ( t ) ,

the Bessel function identity (Prudnikov et al.1992, p. 5, formula

(B6) L - 1 1 s g ̃ ( s - a / s ) = 0 t I 0 2 a t ( t - t ) g ( t ) d t ,

and known Laplace transform pairs (Prudnikov et al.1992), eventually giving Eq. (19).

Appendix C: Calculation of the scale-dependent rate function

The Laplace transform of Eq. (17) over T provides

(C1) Λ ̃ ( s ) = ρ 0 d x i 0 d x P ̃ x + x i , s .

Noting that x+xi is always positive, inserting Eq. (19), and integrating twice gives

(C2) Λ ̃ ( s ) = - ρ φ D V Q s + k E + 2 ρ D φ V Q ( 1 - Q ) s + k E + 4 ρ D 2 ( s + k ) V 3 Q ( 1 - Q ) 2 s + k E .

Taking the inverse transform, applying Eq. (B5), and using the shift property develops

(C3) Λ ( T ) = ρ e - k E T L - 1 - φ D V Q * s + 2 D φ V Q * 1 - Q * s + 4 D 2 T + k V 3 Q * 1 - Q * 2 s .

Here, the notation T means the derivative acts from the left on all terms multiplying it, and

(C4) Q * = 1 + 4 D k D - k E V 2 + 4 D V 2 s - k E k D s .

Laplace inverting Eq. (C3) with Eq. (B6) and tabulated Laplace transforms (e.g., Arfken1985; Prudnikov et al.1992) eventually provides Eq. (21).

Appendix D: Asymptotic limits of the flux rate function

The behavior of Eq. (21) at extreme values of T can be obtained with Tauberian theorems by inverting the Laplace-transformed rate function (Eq. C2) at the opposite extreme of s (Weiss1994). At short times, expanding Eq. (C2) as s→∞ gives

(D1) Λ ̃ ( s ) = ρ k E V 2 k s 2 + ρ k E k D 4 s 3 ,

which inverts to

(D2) Λ ( t ) ρ k E V T 2 k + ρ k E k D T π ,

giving the small T behavior. This has two scaling limits within it. Provided that T2D/V2, the scaling goes as Λ(T)T-1/2. But if T2D/V2, it goes as Λ(t)∼T. For large times, taking s→0 gives

(D3) Λ ̃ ( s ) = ρ k E V k s 2 ,

and this inverts to Λ(T)=ρkEVT/k. These limits are summarized in Eq. (22).

Code availability

Python scripts used for the Monte Carlo simulation of Eq. (5) and to generate all figures have been made publically available at (Pierce et al.2022). These scripts contain comments detailing the stochastic simulation methods.

Data availability

All of the data presented in the paper is freely available at (Pierce et al.2022).

Author contributions

All authors (JKP, MAH, and RMLF) contributed equally to ideation and paper preparation. JKP performed all calculations and constructed all figures.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We would like to thank the two anonymous reviewers for their constructive reviews of the paper.

Review statement

This paper was edited by Jens Turowski and reviewed by two anonymous referees.


Allen, B. and Kudrolli, A.: Granular bed consolidation, creep, and armoring under subcritical fluid flow, Phys. Rev. Fluids, 7, 1–13,, 2018. a

Ancey, C.: Bedload transport: A walk between randomness and determinism. Part 1. The state of the art, J. Hydraul. Res., 58, 1–17,, 2020. a

Ancey, C. and Heyman, J.: A microstructural approach to bed load transport: Mean behaviour and fluctuations of particle transport rates, J. Fluid Mech., 744, 129–168,, 2014. a, b, c, d, e, f

Ancey, C. and Pascal, I.: Estimating mean bedload transport rates and their uncertainty, J. Geophys. Res.-Earth, 125, 1–29,, 2020. a, b, c, d, e

Ancey, C., Böhm, T., Jodeau, M., and Frey, P.: Statistical description of sediment transport experiments, Phys. Rev. E, 74, 1–14,, 2006. a

Ancey, C., Davison, A. C., Böhm, T., Jodeau, M., and Frey, P.: Entrainment and motion of coarse particles in a shallow water stream down a steep slope, J. Fluid Mech., 595, 83–114,, 2008. a, b, c, d

Andrews, E. D.: Marginal bedload transport in a gravel bed stream, Sagehen Creek, California, Water Resour. Res., 30, 2241–2250,, 1994. a

Arfken, G.: Mathematical Methods for Physicists, Academic Press, Inc., San Diego,, 1985. a, b

Bagnold, R. A.: An approach to the sediment transport problem from general physics, Tech. Rep. 4, US Geological Survey, Washington, DC,, 1966. a

Balakrishnan, V., Van Den Broeck, C., and Bena, I.: Stochastically Perturbed Flows: Delayed and Interrupted Evolution, Stoch. Dynam., 01, 537–551,, 2001. a

Ballio, F., Nikora, V., and Coleman, S. E.: On the definition of solid discharge in hydro-environment research and applications, J. Hydraul. Res., 52, 173–184,, 2014. a

Ballio, F., Pokrajac, D., Radice, A., and Hosseini Sadabadi, S. A.: Lagrangian and Eulerian description of bed load transport, J. Geophys. Res.-Earth, 123, 384–408,, 2018. a

Banerjee, T., Majumdar, S. N., Rosso, A., and Schehr, G.: Current fluctuations in noninteracting run-and-tumble particles in one dimension, Phys. Rev. E, 101, 1–16,, 2020. a

Barik, D., Ghosh, P. K., and Ray, D. S.: Langevin dynamics with dichotomous noise; Direct simulation and applications, J. Stat. Mech.: Theory Exp., 2006, P03010,, 2006. a

Bartlett, M.  S. and Porporato, A.: State-dependent jump processes: Itô–Stratonovich interpretations, potential, and transient solutions, Phys. Rev. E, 98, 1–16,, 2018. a

Bena, I.: Dichotomous Markov noise: Exact results for out-of-equilibrium systems, Int. J. Modern Phys. B, 20, 2825–2888,, 2006. a, b, c

Benavides, S. J., Deal, E., Rushlow, M., Venditti, J. G., Zhang, Q., Kamrin, K., and Perron, J. T.: The Impact of Intermittency on Bed Load Sediment Transport, Geophys. Res. Lett., 49, e2021GL096088,, 2022. a

Böhm, T., Ancey, C., Frey, P., Reboud, J. L., and Ducottet, C.: Fluctuations of the solid discharge of gravity-driven particle flows in a turbulent stream, Phys. Rev. E, 69, 061307,, 2004. a

Brach, R. M.: Rigid body collisions, Am. Soc. Mech. Eng., 56, 133–138,, 1989. a

Bradley, D. N. and Tucker, G. E.: The storage time, age, and erosion hazard of laterally accreted sediment on the floodplain of a simulated meandering river, J. Geophys. Res.-Earth, 118, 1308–1319,, 2013. a

Brinkman, H. C.: Brownian motion in a field of force and the diffusion theory of chemical reactions, Physica, 22, 149–155,, 1956. a

Bunte, K. and Abt, S. R.: Effect of sampling time on measured gravel bed load transport rates in a coarse-bedded stream, Water Resour. Res., 41, 1–12,, 2005. a, b

Campagnol, J., Radice, A., Ballio, F., and Nikora, V.: Particle motion and diffusion at weak bed load: Accounting for unsteadiness effects of entrainment and disentrainment, J. Hydraul. Res., 53, 633–648,, 2015. a, b

Cardy, J.: Reaction-diffusion processes, in: Non-Equilibrium Statistical Physics and Turbulence, edited by: Nazarenko, S. and Zaboronski, O. V., Cambridge University Press, Cambridge, 108–161,, 2008. a

Celik, A. O., Diplas, P., Dancey, C. L., and Valyrakis, M.: Impulse and particle dislodgement under turbulent flow conditions, Phys. Fluids, 22, 1–13,, 2010. a

Celik, A. O., Diplas, P., and Dancey, C. L.: Instantaneous pressure measurements on a spherical grain under threshold flow conditions, J. Fluid Mech., 741, 60–97,, 2014. a

Chandrasekhar, S.: Stochastic problems in physics and astronomy, Rev. Mod. Phys., 15, 1–89, doi10.1103/RevModPhys.15.1, 1943. a

Church, M.: Bed material transport and the morphology of alluvial river channels, Annu. Rev. Earth Planet. Sci., 34, 325–354,, 2006. a

Church, M. and Hassan, M. A.: Mobility of bed material in Harris Creek, Water Resour. Res., 38, 1237,, 2002. a

Church, M., Hassan, M. A., and Wolcott, J. F.: Stabilizing self-organized structures in gravel-bed stream channels: Field and experimental observations, Water Resour. Res., 34, 3169–3179,, 1998. a

Clark, A. H., Shattuck, M. D., Ouellette, N. T., and O'Hern, C. S.: Role of grain dynamics in determining the onset of sediment transport, Phys. Rev. Fluids, 2, 034305,, 2017. a

Coffey, W., Kalmykov, Y., and Waldron, J.: The Langevin equation, 1st Edn., World Scientific, London,, 2004. a

Cox, D. R. and Miller, H. D.: The Theory of Stochastic Processes, Wiley, New York,, 1965. a, b, c

Cudden, J. R. and Hoey, T. B.: The causes of bedload pulses in a gravel channel: The implications of bedload grain-size distributions, Earth Surf. Proc. Land., 28, 1411–1428,, 2003. a

Dhont, B. and Ancey, C.: Are bedload transport pulses in gravel bed rivers created by bar migration or sediment waves?, Geophys. Res. Lett., 45, 5501–5508,, 2018. a

Diplas, P., Dancey, C. L., Celik, A. O., Valyrakis, M., Greer, K., and Akar, T.: The role of impulse on the initiation of particle movement under turbulent flow conditions, Science, 322, 717–720,, 2008. a

Duderstadt, J. J. and Martin, W. R.: Transport Theory, 1st Edn., Wiley Interscience, New York, ISBN 0-471-04492-X, 1979. a

Dwivedi, A., Melville, B., and Shamseldin, A. Y.: Hydrodynamic forces generated on a spherical sediment particle during entrainment, J. Hydraul. Eng., 136, 756–769,, 2010. a

Einstein, H. A.: Bed load transport as a probability problem, PhD thesis, ETH Zurich, Zurich, 1937. a, b, c

Einstein, H. A.: The bedload function for sediment transportation in open channel flows, Tech. rep., United States Department of Agriculture, Washington, DC, ISBN TB1026, (last access: 24 July 2022), 1950. a, b, c, d

Escaff, D., Toral, R., Van Den Broeck, C., and Lindenberg, K.: A continuous-time persistent random walk model for flocking, Chaos, 28, 1–11,, 2018. a

Fan, N., Zhong, D., Wu, B., Foufoula-Georgiou, E., and Guala, M.: A mechanistic-stochastic formulation of bed load particle motions: From individual particle forces to the Fokker-Planck equation under low transport rates, J. Geophys. Res.-Earth, 119, 464–482,, 2014. a, b, c, d, e

Fan, N., Singh, A., Guala, M., Foufoula-Georgiou, E., and Wu, B.: Exploring a semimechanistic episodic Langevin model for bed load transport: Emergence of normal and anomalous advection and diffusion regimes, Water Resour. Res., 52, 2789–2801,, 2016. a

Fathel, S. L., Furbish, D. J., and Schmeeckle, M. W.: Experimental evidence of statistical ensemble behavior in bed load sediment transport, J. Geophys. Res.-Earth, 120, 2298–2317,, 2015. a, b, c

Furbish, D. J., Ball, A. E., and Schmeeckle, M. W.: A probabilistic description of the bed load sediment flux: 4. Fickian diffusion at low transport rates, J. Geophys. Res.-Earth, 117, 1–13,, 2012. a, b

Furbish, D. J., Fathel, S. L., and Schmeeckle, M. W.: Particle motions and bedload theory: The entrainment forms of the flux and the Exner equation, in: Gravel-Bed Rivers: Process and Disasters, chap. 4, 1st Edn., edited by: Tsutsumi, D. and Laronne, J. B., John Wiley & Sons Ltd., 97–120,, 2017. a, b, c

Furbish, D. J., Roering, J. J., Doane, T. H., Roth, D. L., Williams, S. G., and Abbott, A. M.: Rarefied particle motions on hillslopes – Part 1: Theory, Earth Surf. Dynam., 9, 539–576,, 2021. a

Gardiner, C. W.: Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, Springer-Verlag, Berlin,, 1983. a

Gitterman, M.: The noisy oscillator: The first 100 years, from Einstein until now, World Scientific,, 2005. a

Goh, K. I. and Barabási, A. L.: Burstiness and memory in complex systems, Europhys. Lett., 81, 48002,, 2008. a

Guala, M., Singh, A., Badheartbull, N., and Foufoula-Georgiou, E.: Spectral description of migrating bed forms and sediment transport, J. Geophys. Res.-Earth, 119, 123–137,, 2014. a

Hamamori, A.: A theoretical investigation on the fluctuations of bed load transport, Technical Report No. R4, Tech. rep., Delft Hydraulics Laboratory, Delft, (last access: 24 July 2022), 1962. a

Hänggi, P.: The functional derivative and its use in the description of noisy dynamical systems, in: Stochastic processes applied to physics, edited by: Pesquera, L. and Rodriguez, M. A., World Scientific Publishing Company, Santander, Spain, 69 pp., ISBN 9971978202, 1985. a

Hassan, M. A., Church, M., and Schick, A. P.: Distance of movement of coarse particles in gravel bed streams, Water Resour. Res., 27, 503–511,, 1991. a

Hernández-García, E. and López, C.: Clustering, advection, and patterns in a model of population dynamics with neighborhood-dependent rates, Phys. Rev. E, 70, 016216,, 2004. a

Heyman, J.: A study of the spatio-temporal behaviour of bed load transport rate fluctuations, PhD thesis, École Polytechnique Fédérale de Lausanne, Lausanne,, 2014. a

Heyman, J., Ma, H. B., Mettra, F., and Ancey, C.: Spatial correlations in bed load transport: Evidence, importance, and modeling, J. Geophys. Res.-Earth, 119, 1751–1767,, 2014. a

Heyman, J., Bohorquez, P., and Ancey, C.: Entrainment, motion, and deposition of coarse particles transported by water over a sloping mobile bed, J. Geophys. Res.-Earth, 121, 1931–1952,, 2016. a, b, c, d, e

Horsthemke, W. and Lefever, R.: Noise-Induced Transitions: Theory and Applications in Physics, Chemistry and Biology, in: 1st Edn., Springer-Verlag, Berlin,, 1984. a, b, c

Houssais, M., Ortiz, C. P., Durian, D. J., and Jerolmack, D. J.: Onset of sediment transport is a continuous transition driven by fluid shear and granular creep, Nat. Commun., 6, 1–8,, 2015. a, b

Hubbell, D. and Sayre, W.: Sand Transport Studies with Radioactive Tracers, J. Hydraul. Div., 90, 39–68,, 1964. a

Iseya, F. and Ikeda, H.: Pulsations in bedload transport rates induced by longitudinal sediment sorting: A flume study using sand and gravel mixtures, Geograf. Ann. A, 69, 15–27,, 1987. a

Ji, C., Munjiza, A., Avital, E., Xu, D., and Williams, J.: Saltation of particles in turbulent channel flow, Phys. Rev. E, 89, 1–14,, 2014. a

Kalinske, A. A.: Movement of sediment as bed load in rivers, EOS Trans. Am. Geophys. Union, 28, 615–620,, 1947. a

Kramers, H. A.: Brownian motion in a field of force and the diffusion model of chemical reactions, Physica, 7, 284–304,, 1940. a, b

Kubo, R., Toda, M., and Hashitsume, N.: Statistical Physics II: Nonequilibrium Statistical Physics, 1st Edn., Springer-Verlag, Berlin,, 1978. a

Kumaran, V.: Granular flow of rough particles in the high-Knudsen-number limit, J. Fluid Mech. 561, 43–72,, 2006. a

Laio, F., Ridolfi, L., and D'Odorico, P.: Noise-induced transitions in state-dependent dichotomous processes, Phys. Rev. E, 78, 1–7,, 2008. a

Lajeunesse, E., Malverti, L., and Charru, F.: Bed load transport in turbulent flow at the grain scale: Experiments and modeling, J. Geophys. Res.-Earth, 115, F04001,, 2010. a

Lajeunesse, E., Devauchelle, O., and James, F.: Advection and dispersion of bed load tracers, Earth Surf. Dynam., 6, 389–399,, 2017. a, b

Laskin, N.: Non-Gaussian diffusion, J. Phys. A, 22, 1565–1576,, 1989. a

Lee, D. B. and Jerolmack, D.: Determining the scales of collective entrainment in collision-driven bed load, Earth Surf. Dynam. 6, 1089–1099,, 2018. a

Lisle, I. G., Rose, C. W., Hogarth, W. L., Hairsine, P. B., Sander, G. C., and Parlange, J. Y.: Stochastic sediment transport in soil erosion, J. Hydrol., 204, 217–230,, 1998. a, b

Lisle, T. E., Iseya, F., and Ikeda, H.: Response of a channel with alternate bars to a decrease in supply of mixed‐size bed load: A flume experiment, Water Resour. Res., 29, 3623–3629,, 1993. a

Liu, M. X., Pelosi, A., and Guala, M.: A statistical description of particle motion and rest regimes in open-channel flows under low bedload transport, J. Geophys. Res.-Earth, 124, 2666–2688,, 2019. a

Łuczka, J.: Non-Markovian stochastic processes: Colored noise, Chaos, 15, 026107,, 2005. a

Łuczka, J., Niemiec, M., and Piotrowski, E.: Linear systems with randomly interrupted Gaussian white noise, J. Phys. A, 26, 4849–4861,, 1993. a

Madej, M. A., Sutherland, D. G., Lisle, T. E., and Pryor, B.: Channel responses to varying sediment input: A flume experiment modeled after Redwood Creek, California, Geomorphology, 103, 507–519,, 2009. a

Mao, L.: The effect of hydrographs on bed load transport and bed sediment spatial arrangement, J. Geophys. Res.-Earth, 117, 1–16,, 2012. a

Martin, R. L., Jerolmack, D. J., and Schumer, R.: The physical basis for anomalous diffusion in bed load transport, J. Geophys. Res.-Earth, 117, 1–18,, 2012. a, b

Masoliver, J.: Second-order processes driven by dichotomous noise, Phys. Rev. A, 45, 706–713,, 1992. a

Masoliver, J.: Second-order dichotomous processes: Damped free motion, critical behavior, and anomalous superdiffusion, Phys. Rev. E, 48, 121–135,, 1993. a

McDowell, C. and Hassan, M. A.: The influence of channel morphology on bedload path lengths: Insights from a survival process model, Earth Surf. Proc. Land., 45, 2982–2997,, 2020. a

Menzel, A. M. and Goldenfeld, N.: Effect of Coulombic friction on spatial displacement statistics, Phys. Rev. E, 84, 1–9,, 2011. a

Michaelides, E. E.: Review – The transient equation of motion for particles, bubbles, and droplets, J. Fluids Eng. T. ASME, 119, 233–247,, 1997. a

Montroll, E. W.: Random walks on lattices, J. Math. Phys., 6, 193–220,, 1964. a

Nakagawa, H. and Tsujimoto, T.: On probabilistic characteristics of motion of individual sediment particles on stream beds, in: Hydraulic Problems Solved by Stochastic Methods: Second International IAHR Symposium on Stochastic Hydraulics, Water Resources Publications, Lund, Sweden, 293–320, ISBN 0918334225, 1976. a

Oldmeadow, D. F. and Church, M.: A field experiment on streambed stabilization by gravel structures, Geomorphology, 78, 335–350,, 2006. a

Olivares-Robles, M. A. and García-Colin, L. S.: On different derivations of telegrapher's type kinetic equations, J. Non-Equilib. Thermodynam., 21, 361–379,, 1996. a

Papanicolaou, A. N., Tsakiris, A. G., Wyssmann, M. A., and Kramer, C. M.: Boulder array effects on bedload pulses and depositional patches, J. Geophys. Res.-Earth, 123, 2925–2953,, 2018. a

Parker, G., Paola, C., and Leclair, S.: Probabilistic Exner sediment continuity equation for mixtures with no active layer, J. Hydraul. Eng., 126, 818–826,, 2000. a

Pechenik, L. and Levine, H.: Interfacial velocity corrections due to multiplicative noise, Phys. Rev. E, 59, 3893–3900,, 1999. a

Pierce, J. K.: The stochastic movements of individual streambed grains, PhD thesis, The University of British Columbia,, 2021. a, b

Pierce, J. K. and Hassan, M. A.: Joint stochastic bedload transport and bed elevation model: Variance regulation and power law rests, J. Geophys. Res.-Earth, 125, 1–15,, 2020a. a

Pierce, J. K. and Hassan, M. A.: Back to Einstein: Burial-induced three-range diffusion in fluvial sediment transport, Geophys. Res. Lett., 47, e2020GL087440,, 2020b. a

Pierce, J. K., Hassan, M. A., and Ferreira, R. M. L.: Code for ESurf submission “Stochastic description of intermittent transport and aggregate derivation of the bedload flux”, Zenodo [code],, 2022. a, b

Prudnikov, A. P., Brychkov, I. A., Marichev, O. I., and Gould, G. G.: Integrals and Series Vol. 5: Inverse Laplace Transforms, Gordon and Breach Science Publishers, Paris, ISBN 2881246826, 1992. a, b, c

Recking, A., Liébault, F., Peteuil, C., and Jolimet, T.: Testing bedload transport equations with consideration of time scales, Earth Surf. Proc. Land., 37, 774–789,, 2012. a

Risken, H.: The Fokker-Planck Equation: Methods of Solution and Applications, 2nd edn., Springer-Verlag, Ulm,, 1984. a, b

Roseberry, J. C., Schmeeckle, M. W., and Furbish, D. J.: A probabilistic description of the bed load sediment flux: 2. Particle activity and motions, J. Geophys. Res.-Earth, 117, F03032,, 2012. a

Saletti, M., Molnar, P., Zimmermann, A., Hassan, M. A., and Church, M.: Temporal variability and memory in sediment transport in an experimental step-pool channel, Water Resour. Res., 51, 9325–3337,, 2015. a, b, c, d, e

Schmeeckle, M. W.: Numerical simulation of turbulence and sediment transport of medium sand, J. Geophys. Res.-Earth, 119, 1240–1262,, 2014. a

Schmeeckle, M. W. and Nelson, J. M.: Direct numerical simulation of bedload transport using a local, dynamic boundary condition, Sedimentology, 50, 279–301,, 2003. a

Schmeeckle, M. W., Nelson, J. M., and Shreve, R. L.: Forces on stationary particles in near-bed turbulent flows, J. Geophys. Res.-Earth, 112, 1–21,, 2007. a, b

Schumer, R., Meerschaert, M. M., and Baeumer, B.: Fractional advection-dispersion equations for modeling transport at the Earth surface, J. Geophys. Res.-Earth, 114, 1–15,, 2009. a

Singh, A., Fienberg, K., Jerolmack, D. J., Marr, J., and Foufoula-Georgiou, E.: Experimental evidence for statistical scaling and intermittency in sediment transport rates, J. Geophys. Res.-Earth, 114, 1–16,, 2009. a, b, c, d, e, f

Strom, K., Papanicolaou, A. N., Evangelopoulos, N., and Odeh, M.: Microforms in gravel-bed rivers: formation, disintegration, and effects on bedload transport, J. Hydraul. Eng., 130, 554–567,, 2004. a

Taylor, G. I.: Diffusion by continuous movements, Proc. Lond. Math. Soc., s2–20, 196–212,, 1920. a

Turowski, J. M.: Probability distributions of bed load transport rates: A new derivation and comparison with field data, Water Resour. Res., 46, W08501,, 2010. a, b, c, d, e, f

Van Kampen, N. G.: Stochastic Processes in Physics and Chemistry, in: 3rd Edn., Elsevier B. V.,, 2007. a, b

Wang, M. C. and Uhlenbeck, G.: On the theory of Brownian motion II, Rev. Modern Phys., 17, 323–342,, 1945. a

Weiss, G. H.: Aspects and Applications of the Random Walk, North Holland, Amsterdam,, 1994. a, b, c

Weiss, G. H.: Some applications of persistent random walks and the telegrapher's equation, Physica A, 311, 381–410,, 2002. a

Williams, S. G. W. and Furbish, D. J.: Particle energy partitioning and transverse diffusion during rarefied travel on an experimental hillslope, Earth Surf. Dynam., 9, 701–721,, 2021. a

Wong, M. and Parker, G.: Reanalysis and correction of bed-load relation of Meyer-Peter and Müller using their own database, J. Hydraul. Eng., 132, 1159–1168,, 2006.  a

Wu, Z., Furbish, D., and Foufoula-Georgiou, E.: Generalization of hop distance-time scaling and particle velocity distributions via a two-regime formalism of bedload particle motions, Water Resour. Res., 56, e2019WR025116,, 2020. a, b

Wu, Z., Singh, A., Foufoula-Georgiou, E., Guala, M., Fu, X., and Wang, G.: A velocity-variation-based formulation for bedload particle hops in rivers, J. Fluid Mech., 912, A33,, 2021. a

Zhang, Y., Meerschaert, M. M., and Packman, A. I.: Linking fluvial bed sediment transport across scales, Geophys. Res. Lett., 39, 1–6,, 2012. a

Short summary
We describe the flow of sediment in river channels by replacing the complicated details of the turbulent water with probability arguments. Our major conclusions are that (1) sediment transport can be phrased in terms of the movements of individual sediment grains, (2) transport rates in river channels are inherently uncertain, and (3) sediment transport in rivers is directly analogous to a number of phenomena which we understand relatively well, such as molecules moving in air.