Stochastic bedload transport in mountain streams

. Describing bedload transport as a stochastic process is an idea that emerged in the 1930s with the pioneering work of Einstein. For a long time, the stochastic approach attracted marginal attention, but the situation has radically changed over the last decade with the recent advances in the theory of bedload transport. In parallel, the implementation of bedload monitoring techniques at high temporal resolution has produced a wealth of interesting results showing, among other things, that classic empirical bedload transport equations do not capture neither the mean behavior of sediment transport rates q s nor its order of magnitude, especially at low sediment transport rates (a case that is most frequent in mountain streams). We have developed a stochastic model, which takes in-spiration from population dynamics and provides a stochastic partial di ﬀ erential equation for the number of moving particles. Taking the ensemble average leads to a fairly simple advection di ﬀ usion equation for particle activity (i.e., the number of moving particles per unit streambed area). The model has a number of unique features. For instance, it yields the probability distribution of the bed-load transport rate and predicts bedform formation for a wide range of Froude numbers.


Introduction
In recent years, sediment transport has been monitored in a number of mountain streams using geophones, which has made it possible to measure bedload transport rate continuously at a high sampling rate [1][2][3][4].Measurements have shown that bedload transport rates q s exhibit considerable temporal variability even at constant water discharge q w , with the fluctuation scale often exceeding one order of magnitude.Surprisingly this has not breached the trust placed in empirical equations that express bedload transport rate as a one-to-one scalar function of water discharge q s = f (q w ).In his historical perspective of bedload transport theory, Gomez already noted the discrepancy between the great variability in transport rates and the simple bedload transport equations commonly used [5].The paradox could be merely apparent and easily resolved.Indeed, many physical systems are affected by fluctuations (noise), but when averaged over appropriate space or time scales, these fluctuations cancel out, and it is permitted to consider only the mean values as the driving dynamical parameters of the system.It does not matter whether or not this has been intentional, but the fact remains that the issue of rate fluctuations has gone unnoticed for years.
Yet, reality is hard-headed.The issue of rate fluctuations reemerged when one tried to measure beadload transport rates using samplers or baskets.In these techniques, the key point in the protocol is the time during which sediment is collected in the basket: how long should one wait for the measured transport rate to be fairly constant?There is no easy solution to this issue as the measurements showed considerable dependence on sampling time [6].A similar problem arises when trying to compare field data (e.g., volumes accumulated over timescales ranging from second to year) and empirical bedload equations [7]: integrating bedload transport equations over time does not yield the right value of transported volume.In both cases, as sediment transport results from many nonlinear fluctuating processes, time integration provides random outcomes that can deviate significantly from the mean trend.
In this paper, we review the approach we have developed to cope with fluctuations in bedload transport, summarize the main achievements and outline the perspectives.

Governing equations
There are various approaches to stochastic sediment transport modelling.A pragmatic approach is to start from the commonly accepted governing equations (Saint-Venant and Exner) expressing mass and momentum conservation at the bulk scale, and then see how they can be extended to include fluctuation-driven contributions [8,9].Another one, probably more demanding and risky, could also be more rewarding: if one is able to model particle transport at the microscale, then one should be able to infer their time-or volume-averaged behaviour and investigate its interplay with the water stream.This idea has attracted growing interest over the years, with the development of Lagrangian and Eulerian models [10][11][12].This is the approach we have taken and we will summarize here.

Mass conservation for bedload transport
The initial idea is to consider a fixed control volume of length Δx, and focus on mass conservation in a statistical manner.We consider two populations of particles: moving and resting particles.At time t, there are N particles in the control volume.This number can be increased when particles enter the volume or are entrained from the bed; N decreases when particles leave the volume or are deposited.To describe these time variations, a convenient framework is the theory of birth-death Markov processes, widely used in population-dynamics models or chemical kinetics [13]: entrainment or arrival is equivalent to birth or immigration, respectively, whereas deposition or departure corresponds to death or emigration.
If there are N moving particles in the control volume, the probability of deposition within the time increment δt is σNδt, with σ the deposition rate.For entrainment, we assume that there are two processes referred to as individual and collective entrainment resulting in a probability of entrainment P = (λ � + µN)δt, where λ � and µ denote the individual and collective entrainment rates, respectively.Collective entrainment acts as a feedback loop: µ is a key parameter, which controls the development and strength of wide fluctuations.A caveat is in order: here, collective entrainment implies that the probability of entrainment depends not only on the flow conditions (through λ � ), but also on the number of moving particles (through µ) as these can impact the bed and impart momentum to the bed particles, favouring their entrainment.It does not mean that there are massive departures of particles (avalanches) within short time spans.For subsequent use, we also introduce a volumetric particle entrainment rate per unit length λ = λ � � p /Δx and the differential rate κ = σ − µ between deposition and collective entrainment, with � p the particle volume per unit width.
Under these assumptions and neglecting particle motion (emigration and immigration) for the moment, we end up with the "forward master equation," that describes the evolution of the number of moving particles resulting from entrainment and deposition [11]: A stumbling block in this approach is that the governing equation for N involves discrete probabilities.To generalize the model and derive a continuum formulation, we wish to replace the discrete variable N with a continuous variable a.To that end, we use a trick that involves using a kind of Fourier transform.Fourier transforms are reversible operations that map the time and frequency domains in spectral theory of signals.On many occasions when working with times series, it is easier to work in the frequency domain than the original time domain.Similarly here, we can use the Poisson transform that maps the discrete and continuous probability domains [14] where integration is made over a certain domain C and f (a, t) is a positive real-valued function.Using this transform, we have shown that the master equation ( 1) can be transformed into a second-order nonlinear parabolic diffusion equation, which has the same structure as that of a Fokker-Planck equation [11] ∂ with f (a) the transform of P in the a-space.The continuous random variable a is called the Poisson rate for the control volume of size Δx.We can now derive a continuum variant of this equation by introducing the Poisson density b = lim Δx→0 a/Δx.We now also include particle motion by assuming that particles move at a velocity of mean ūp and variance D u .We end up with a generalized form for (3) in the form of a stochastic partial differential equation where The decisive advantage of this formulation is that whereas we cannot easily calculate the inverse Poisson transform, there exists a simple relation between the moments of a and N, or between b and particle activity (the volume of moving particles per unit bed length) γ(x, t) = N� p /Δx: �a� = �a� and �b� = �γ�.By taking the ensemble average of (4), we obtain the governing equation for �γ� that can be seen as the mass conservation for bedload: This is a linear advection diffusion equation with a source term.Albeit of very common structure, this equation yields many interesting insights into the physics of sediment transport.We can make a link between the Exner equation and this equation.Indeed, (5) can be cast in the following form with Q = ūp �γ� − ∂ x (D u �γ�), E = λ + µ�γ�, and D = σ�γ�.Interestingly, if we borrow the definition of the sediment flux rate from Furbish et al. [12] and refer to Q as the macroscopic sediment transport rate, then ( 6) is the generalized Exner equation: where y b (x, t) denotes the bed elevation, x the downstream position, t time, ζ b the bed porosity, and where we identified qs as Q.Note that equation ( 7) does not usually include the time variation in the particle activity ∂ t �γ� as this term is vanishingly small.

Mass and momentum balance equations for the water stream
For the water stream, the governing equations comprise the conservation of mass and momentum balance equations (Saint-Venant) [15]: in which h(x, t) = y s − y b denotes the flow depth, y s (x, t) the free surface position, v the depthaveraged velocity, � the water density, τ b is the bottom shear stress, g is the gravitational acceleration and ν the eddy viscosity.The bed slope is defined as tan θ = −∂ x y b .

Closure equations
We need to specify all the parameters involved in the governing equations.The bedload transport model includes 5 parameters: the entrainment and deposition rates (λ, µ, σ), the mean particle velocity ūp , the particle diffusivity D u .Preliminary tests have made it possible to relate these parameters to the flow variables (h, v) [10,16]: λ and σ vary almost linearly with v, whereas µ is nearly constant for the flow range explored.For the bottom shear stress, we use the Darcy-Weisbach equation and a simple mixing length model for ν [17].

Numerical solutions
Our framework enables two kinds of computation: (i) if we use the stochastic partial differential equation ( 4) in combination with the Saint-Venant-Exner equations, then we have to solve a stochastic system of equations; (ii) we can take the ensemble-averaged formulation, in which (4) is replaced with (6), and thus have to solve a deterministic set of hyperbolic differential equations.In both cases, we used the finite-volume library SharpClaw [15,17].The stochastic equation (4) was solved using an Euler-Maruyama scheme [17].Note that integrating (4) required special attention when finite-volume techniques are used (pseudo-diffusive terms arise due to the correlation between neighbouring cells) [18].

Model predictions
We illustrate our approach's potential by focusing on two applications.[19] and Birnbaum-Saunders distribution [20].We also report the probability density function P ṅ( ṅ) for ζ = 5 (ζ being a free dimensionless parameter describing velocity fluctuation strength, see [11]): when the number of moving particles follows the negative binomial distribution (solid red line) or the Poisson distribution (dashed red line).Except for the Poisson distribution (whose variance equals the mean), the coefficient of variation is √ 2 and all of the distributions have the same mean (�ṅ� = 1 bead s −1 ).Drawn from [11].

Fluctuations in the bedload transport rate
When knowing N, we can define the instantaneous particle flux as the number of moving particles per unit time in the control volume ṅ = � N(t) i=1 U p,i /Δx, where both the number of moving particles N and particle velocities U p,i are random variables.The bedload transport rate q s is related to ṅ: q s = � p ṅ.For a stationary process, the probability distribution of N is given by a negative binomial distribution-or a Poisson distribution if µ = 0-, while that of U p is given by a truncated normal distribution [11].After a little bit of work, we can derive a closed-form expression for ṅ [11].For the sake of conciseness, we will merely illustrate its properties in Fig. 1, and compare it with the Hamamori and Birnbaum-Saunders distributions [19,20].High-resolution data confirm (i) the significant proportion of zero values of the particle flux and (ii) the highly fluctuating nature of time series, two features that are consistently described by our model [10,[21][22][23].

Bedform development
A longstanding issue in morphodynamics is related to bedforms, their formation and migration.It can be shown that the bed is unconditionally stable for Froude numbers Fr < 2 when the Saint-Venant equations ( 8)-( 9) are coupled with the standard Exner equation (1 − ζ b )∂ t y b = −∂ x qs , with qs given by a scalar relation qs = f (q w ) [24].By contrast, the SVE equations ( 8)-( 9) together with ( 7) and ( 6) are unstable and catch the most unstable wavelength: as shown by Fig. 2 (supercritical shallow water flow over a sloping erodible gravel bed), initially planar beds become unstable, and very quickly, antidunes develop and migrate upstream.Linear stability analysis of these equations reveals that the flow is absolutely unstable.

Concluding remarks
A number of theoretical or numerical problems have been studied in our recent papers to determine how well the model performs at predicting bedload transport [11,15,17,18]: (i) the derivation of the stochastic evolution equation for the number of moving particles over fixed plane beds, which leads to exact analytical solutions of the particle activity and velocity fluctuations [11]; (ii) the nonlocal pseudo-diffusive behaviour of bedload transport due to particle entrainment and deposition [18]; ; (iii) the nonlinear simulation of the stochastic SVE equations, which is considered a major challenge, and stability analysis of the stochastic SVE showing the absolute character of bed instability [15]; (iv) the nonlinear simulation of the ensemble-averaged Saint-Venant-Exner equations (SVE), which successfully captures the anti-dune regime observed experimentally [17]; (v) numerical simulation of aggradation in good agreement with available data [17].
What we have learned from these studies is that the stochastic approach holds promise for predicting bed morphology and transport rates under various flow conditions provided that sediment transport is not too intense.A remarkable achievement has been the proof that the SVE equations are absolutely unstable and develop bedforms as soon as the Shields number exceeds a threshold.In other models, this result can be achieved only by significantly affecting the structure of the governing equations or employing ad hoc assumptions.
Yet, more work is required to explore the full potential of our stochastic approach to bedload transport.Among other things, this involves studying the coupling between bedload transport and water stream by determining how the model parameters (λ, µ, σ, ūp , D u ) are related to the flow variables (h, v).Experiments have been carried out [10,16,25,26], but more work is needed.Experiments have also revealed a number of critical points in the current framework.First, a key assumption of our developments based on Markov processes is that at most, a single event occurs over a short time increment δt and particle movements are statistically independent.Yet there is growing evidence that even at low sediment transport rates, correlated motion and bursts, i.e., several particles entrained and moving together, can occur [23].This means that the statistical framework has to be extended to cope with multiple events and correlated motion.
Second, we have mainly focused on one-dimensional flows, but most flows of practical importance are two-dimensional.In two-dimensional flows, e.g., flows around alternate bars, sediment transport does not occur necessarily in the same direction as the water stream [23].Vortices can also cause the particles to move against the stream [27].This makes the coupling between flow and particles much more complicated.
Third, a number of processes such as hyporheic flows and grain sorting have been omitted.We are currently studying them in isolation.While grain sorting (segregation) can be modelled using continuum models [28], its inclusion into a stochastic approach is not straightforward.
Fourth, our framework enables computation of individual realizations (by solving the stochastic partial differential equation ( 4) in combination with the Saint-Venant-Exner equations) and ensemble-averaged properties (by solving the advection diffusion equation ( 6) together with the Saint-Venant-Exner equations).If we want to obtain more statistical information (e.g., the variance or the probability distribution of a key variable), then we have to run the stochastic model and compute a large number of realizations, which is time consuming.Comparing the model outcomes with field data (that can be seen as realizations) is not as direct as with deterministic models.This will probably require thinking of specific techniques as those used in uncertainty propagation [29].
Last, other approaches have been developed [12,30].The Lagrangian description proposed by Furbish and coworkers leads to similar results for the definition of the bedload transport rate, but some of their experiments show features that have been not in our experiments (e.g., the probability distribution of particle velocity).This also calls for further analysis [16].

E3SFigure 1 .
Figure1.Comparison of the probability density function P(ṅ) in a log-linear plot: Hamamori's equation (with �ṅ� = 1 bead s −1 )[19] and Birnbaum-Saunders distribution[20].We also report the probability density function P ṅ( ṅ) for ζ = 5 (ζ being a free dimensionless parameter describing velocity fluctuation strength, see[11]): when the number of moving particles follows the negative binomial distribution (solid red line) or the Poisson distribution (dashed red line).Except for the Poisson distribution (whose variance equals the mean), the coefficient of variation is √ 2 and all of the distributions have the same mean (�ṅ� = 1 bead s −1 ).Drawn from[11].

E3SFigure 2 .
Figure 2. (a) Snapshots of the bed and free-surface elevations together at t = 0 and at t = 200 s: dashed lines show the uniform base flow at t = 0; solid lines show the anti-dune train in the numerical simulation at t = 200 s.The gray and black lines correspond to the free surface and bed elevation, respectively.(b) Evolution of the maximum perturbation in the bed elevation.(c) Convection contribution to the bedload transport rate in the plane {x, t} scaled by the uniform background flow's steady-state transport rate qss .Drawn from [15].