Bedrock incision by bedload: insights from direct numerical simulations

. Bedload sediment transport is one of the main processes that contribute to bedrock incision in a river and is therefore one of the key control parameters in the evolution of mountainous landscapes. In recent years, many studies have addressed this issue through experimental setups, direct measurements in the ﬁeld, or various analytical models. In this article, we present a new direct numerical approach: using the classical methods of discrete-element simulations applied to granular materials, we explicitly compute the trajectories of a number of pebbles entrained by a turbulent water stream over a rough solid surface. This method allows us to extract quantitatively the amount of energy that successive impacts of pebbles deliver to the bedrock, as a function of both the amount of sediment available and the Shields number. We show that we reproduce qualitatively the behaviour observed experimentally by Sklar and Dietrich (2001) and observe both a “tool effect” and a “cover effect”. Converting the energy delivered to the bedrock into an average long-term incision rate of the river leads to predictions consistent with observations in the ﬁeld. Finally, we reformulate the dependency of this incision rate with Shields number and sediment ﬂux, and predict that the cover term should decay linearly at low sediment supply and exponentially at high sediment supply.


Introduction
The incision of bedrock channels is one of the key processes that govern the formation and evolution of mountain ranges (Anderson, 1994;Howard, 1994;Whipple and Tucker, 1999).Long-term averaged incision rates can take values from 0.02 to 14 mm yr −1 (see, for instance, the review by Lague, 2014).It has also been observed that, under rarely reached conditions, the short-term incision rate can reach even higher values, up to a few metres per day (Hartshorn et al., 2002;Lamb and Fonstad, 2010;Cook et al., 2014).In order to model the long-term evolution of the morphology of bedrock rivers, and more generally of mountainous landscapes, it is often necessary to adopt a simple macroscopic law to take into account the process of bedrock incision.One of the most commonly used approaches, the stream-power incision model, assumes that the incision rate within a river channel varies as a power law of both its local slope and its drainage area (which is equivalent to introducing a depen-dence in the water discharge) (Seidl et al., 1994;Whipple and Tucker, 1999).The suitability of this model to adequately reproduce several features of bedrock channels has recently been reviewed extensively by Lague (2014).One of its main restrictions is that it does not take into account more detailed parameters such as the dynamics of the alluvial cover in the channel.
In a bedrock mountain river, various processes contribute to incision: chemical dissolution, cavitation, abrasion (or wear) by both bedload and suspended load, plucking, and macroabrasion (Whipple et al., 2000;Chatanantavet and Parker, 2009).Amongst those, abrasion, plucking, and macroabrasion depend directly on the amount of material that is removed from the bedrock by rolling, sliding, or impacting particles transported by the flow, which itself depends mainly on the amount of energy that is transmitted to the bedrock by moving particles (Foley, 1980).This transfer of energy from the impacting particle to the bedrock has been Published by Copernicus Publications on behalf of the European Geosciences Union.recently directly measured experimentally by Turowski and Bloem (2015), as a function of the thickness of the sediment layer covering the bedrock.As can be expected, the fraction of the incipient kinetic energy that is effectively transmitted to the bedrock decreases when the sediment thickness increases.This confirms that, as proposed early on by Gilbert (1877) or Shepherd (1972), the amount of sediment available in the river channel should influence the downcutting rate in two opposite ways: incision should be first enhanced by an increase in the number of impacts of abrasive tools on the bedrock ("tool effect"), but if the supply rate becomes too high, the bedrock should become partially or completely protected from these impacts ("cover effect").
The first direct measurement of the effect of sediment transport on abrasion was performed by Sklar and Dietrich (2001): by measuring the mass loss of a rock disk eroded by a bedload layer of saltating grains in a rotating flow, they confirmed that a maximum abrasion rate is observed for a critical amount of sediment above the bedrock.Following this experimental work, Sklar and Dietrich (2004) developed a mechanistic approach in order to derive the saltationabrasion model.In this model, the incision rate I is written as the product of three terms: the volume of rock eroded by each impact, V i ; the number of impacts by unit time and surface, n i ; and the probability that a saltating grain impacts an exposed area of the bedrock, F : (1) Sklar and Dietrich (2004) derive the frequency of impacts from a mechanistic description of saltation trajectories, all pebbles being assumed to have the same dynamics.They consider that F is the fraction of bedrock not shielded by immobile particles, and varies linearly with the amount of available sediment Q s when below the transport capacity of the stream, and vanishes when the transport capacity Q t is reached: A further model was developed by Turowski et al. (2007), where each saltating grain has a given probability to impact an exposed or covered region of the bedrock, which leads to an exponential expression for F : with ϕ a constant.In both models, the bedload layer is described as made of two distinct populations: static particles that cover and protect the bedrock, and moving particles that all have the same trajectories.Therefore, this analytical approach does not take into account the fact that, when the amount of sediment increases, moving pebbles will interact with each other and with static ones, which should modify their trajectories.An impact is likely to be qualitatively different depending on whether it hits a covered or exposed region of the bedrock: a saltating particle impacting the raw bedrock can bounce and continue its saltating trajectory.Conversely, when it impacts an area already covered with immobile pebbles, it is likely that more energy will be dissipated in the collision, and the impacting particle might not bounce back.The two populations (static and saltating particles) should therefore be permanently interacting, with static particles being ejected by an impact and becoming mobile, while mobile particles can get trapped in asperities of the static cover (Charru et al., 2004).This implies that Eq. ( 1) is somewhat ill-defined, since the number of impacts n i should also be a function of F .
In this article, we propose a new numerical approach of abrasion, based on the discrete-element method.This method allows us to model the individual trajectories of all particles within the bedload layer and therefore to obtain a physically based value of the amount of energy transmitted to the bedrock by impacts.The article in organized as follows.In Sect. 2 we present our numerical setup and the physical laws implemented in our simulations.In Sect. 3 we expose the numerical results regarding the sediment transport rate, the energy delivered to the bedrock and the influence of bedrock roughness.Finally, in Sect. 4 we discuss the implications of our results on the influence of both the Shields number and the sediment supply on the incision rate.We propose a new definition for the cover function F , compare our results to available experimental and analytical models and estimate the long-term incision rate predicted by our simulations.

Numerical setup
We use the discrete-element method to simulate the individual dynamics of pebbles entrained by a turbulent water flow over a fixed bedrock.The same method would allow for modelling of non-spherical particles by considering composite particles made of two or more "glued" spheres, but for the sake of simplicity and to limit the number of control parameters, we restrict our study to the dynamics of spherical particles.
The computational domain is a parallelepipedic box of length L = 2 m, width W = 1 m, and height H = 2 m (see Fig. 1 and Table 1 for the list of all physical parameters used in the simulation).Periodic boundary conditions are used in both horizontal directions x and y: any pebble coming out of the box on one side is reinjected with the same velocity on the other side.The bedrock is modelled as a horizontal surface located at z = 0, over which we simulate a natural roughness of the bedrock by glueing N b spheres of radius R = 5 cm, centred at a height z = z r .These protruding spheres have the same mechanical properties as the pebbles entrained by the flow but are fixed and considered part of the bedrock.In all the presented results except in Sect.3.4, these spheres protrude by a height h b = R+z r = 4 cm and their surface density is given by Considering that we only model spherical particles, the presence of this roughness on the bedrock is necessary to allow for the existence of patches of immobile pebbles (particles roll towards asperities and can get trapped) and also to enhance vertical motion, that is, saltation of mobile pebbles (without any roughness and within a purely horizontal flow, particles tend to simply roll along the smooth surface).The aim of the simulation is to compute the amount of energy that is transmitted to the bedrock when it is hit by saltating pebbles, and to evaluate the erosion of the bedrock induced by these impacts.However, let us note that we simulate the bedload dynamics over a timescale of the order of 1 min, whereas significant abrasion of the bedrock only happens over at least the duration of a flood, that is, a few days, and in many cases over years.Therefore, we can admit that the bedrock (both the horizontal surface and the fixed spheres) remains unchanged and immobile within the timescale of the simulation (in particular, its altitude remains z = 0 throughout the whole duration of the simulation).This can be seen as an advantage compared to experimental studies such as conducted by Johnson and Whipple (2010), where incision affects preferential areas of the bedrock, which rapidly leads to the formation of a narrow inner channel, thus making both shear stress and alluvial highly non-uniform within the whole channel.
The bedload consists of N pebbles that are modelled as spheres of radius R = 5 cm and density ρ s = 2500 kg m −3 .Even if repeated impacts might lead to a slow comminution of the pebbles, their size is considered constant over the duration of the simulation.The number of pebbles that can be disposed in a single layer over the bedrock is of the order of W L/(π R 2 ).Therefore, we quantify the sediment supply by defining the dimensionless surface density σ as the surface of bedrock covered by pebbles, divided by the total available surface: A horizontal turbulent water flow in the x direction puts the pebbles into motion (see Sect. 2.3).Their trajectories are then driven by their immersed weight (W ), fluid friction (drag force F and torque M), and contact forces exerted by other pebbles (N, T ).The evolution with time of the position r and rotational velocity of a pebble is given by Newton's equations of motion: with m = 4π 3 ρ s R 3 the mass and J = 8π 15 ρ s R 5 the angular momentum of a pebble.The immersed weight of a pebble is

Contacts between pebbles
When two pebbles are in contact, the exact deformation of each solid particle is not explicitly computed but spheres are instead allowed to overlap slightly (see, for instance, Pöschel and Schwager, 2005).We assume that two pebbles i and j are in collision if the distance between their centres is lower than the sum of their radii, that is, if δ = 2R − |r i − r j | > 0 (see Fig. 2).When two pebbles make contact, they experience an inelastic rebound that can be modelled by the sum of an elastic and a viscous force (Cundall and Strack, 1979).The elastic force is linear in the overlap δ.The viscous dissipation is proportional to the temporal variation in this overlap: the normal force experienced by a pebble i in contact with a pebble j is then where k is the elastic constant, is the effective viscosity, and is the unit normal vector of the collision.
The value of the elastic constant is related to the material's Young's modulus and the pebble size: we adopt the value k = 2 × 10 8 N m −1 , which corresponds to an elastic modulus Y ∼ k/R = 4 GPa.This value is quite low for rocks, but increasing the elastic modulus would imply reducing the numerical time step too much.Let us note, however, that the pebbles that we model are nevertheless very rigid: the deformation of a pebble under its own weight is only 60 nm.The effective viscosity is a numerical parameter that is responsible for the inelasticity of the collision but does not have a direct physical equivalent.When a pebble impacts another one in water, energy is dissipated in plastic deformations or micro-fractures within the rock (which are responsible for wear), as well as in the viscous interstitial flow.Within our model, the effective inelasticity of the collision can be quantified by the coefficient of restitution e, which compares the velocity of the pebble before and after a collision: e = 1 corresponds to an elastic collision and e = 0 to total dissipation.If the force is given by Eq. ( 8), e is expressed as the typical duration of a collision.We choose the value = 2 × 10 4 kg s −1 for the effective viscosity, which leads to collisions of duration T coll = 10 −4 s and a coefficient of restitution e = 0.3, which means that a pebble loses 1−e 2 = 90 % of its incident kinetic energy during an impact.This value of e lies in the lower range of the experimental observations of Schmeeckle et al. (2001) for natural sediment at high Stokes number, and below the prediction of Davis et al. (1986) for elastic spheres (e = 0.65).However, let us note that, in our simulation, the energy loss is due to not only the viscous dissipation in the film of water that appears between the two particles in elastic collision but also the dissipation induced by the mechanical damage that both particles experience upon impact.In the following (see Eq. 19), we will assume that all the energy lost during an impact with the bedrock contributes to its abrasion.
The tangential force T ij generated at a contact between two pebbles is described by the regularized Coulomb's law of solid friction, as in Cundall and Strack (1979).This force opposes the tangential motion and is expressed as where v s is the sliding velocity at the contact, which is a function of the two pebbles' translational and angular velocities, µ = 0.6 is the local friction coefficient, and G = 5 kg s −1 is the slope of the regularization of the Coulomb's law.This regularization prevents the indetermination of the friction force when the two particles in contact have a zero sliding velocity.
Finally, each collision between a pebble and the horizontal surface of altitude z = 0 is treated as a collision with a pebble of infinite size, and same mechanical properties.

Turbulent water flow
A stationary turbulent flow over a rough bedrock follows the average velocity profile: where x is the direction of the flow and κ = 0.41 is the von Kármán constant.z 0 is the bedrock roughness and depends on pebble size: a bedrock made of pebbles of radius R has a roughness z 0 = R/15 (Nikuradse, 1933;Valance, 2005).U * is the shear velocity, whose expression is given by the relationship between turbulent shear stress and velocity gradient: where ρ w = 1000 kg m −3 is the density of water.The ability of the stream to put pebbles into motion is described by the Shields number.This dimensionless quantity is proportional to the ratio between the drag force on a pebble and its immersed weight: Pebbles are put into motion by the flow if exceeds a threshold value c : measurements of this threshold, both in the field and in experiments, give values in the range 0.01 < c < 0.2 (Buffington and Montgomery, 1997;Lamb et al., 2008).In abrasion experiments conducted in a flume setup by Attal and Lavé (2009) the maximum fluid velocity is 4 m s −1 for a water height H = 60 cm and pebbles of size 10 to 80 mm.As reported in data reviewed by Rickenmann and Recking (2011), flow velocity in mountain streams varies typically between 0.3 and 4 m s −1 , for a water height between 0.1 and 3 m.In order to be consistent with these observations, in our simulations we adopted mean water velocities up to 5.0 m s −1 , which corresponds to varying from 0 to 0.11.
If the bedrock is covered with a layer of mobile pebbles, as in our simulations, the turbulent velocity profile is modified.Recently, Duran et al. (2012) developed a quasi-2-D mechanistic approach that takes into account the retroaction of the pebbles on the water flow by assuming the conservation of total horizontal momentum in horizontal slices.They showed that the fluid velocity vanishes where the local solid fraction is high enough (that is, at a depth of one or two grain diameters within the bedload layer), and tends towards a logarithmic profile in the "clear water" region.The "intermediate" region where the velocity goes from zero to the logarithmic profile is very thin (usually of the order of one grain diameter).Therefore, we simplified the treatment of this retroaction by modelling only two different regions in the flow: at each time step, we compute the average solid fraction φ in horizontal slices.If φ < φ b = 0.5, the velocity profile is logarithmic and not affected by the presence of pebbles.If φ φ b , the velocity of water vanishes: V f = 0 (see the resulting velocity profile in Fig. 5).This approximation would be too simplistic if we were to study the exact flux of grains at the surface of a thick layer, but remains relevant enough in our simulations, where in most cases the bedload layer remains relatively thin.A more detailed description in our geometry would require taking into account not only the average retroaction of grains on the fluid flow but also the horizontal variations in pebble density.

Interactions between pebbles and the flow
The turbulent flow exerts on each mobile pebble a drag force given by where U = V (z)− dr dt is the relative velocity between the local flow and the pebble.The drag coefficient C D of a sphere can be expressed semi-empirically as a function of the particle Reynolds number (Clift et al., 1978).In the present study, we use the following approximation: where Re p = 2ρ w U R η w , with η w = 10 −3 Pa s the dynamic viscosity of water, is the particle Reynolds number.This approximation is equivalent to the Stokes formula for the drag force at low Re p .When a sphere is rotating in a viscous fluid such as water, its angular velocity induces a diffusion of momentum in a boundary layer.This results in a viscous torque applied to the pebble (Liu and Prosperetti, 2010), which opposes its rotation:

Computational methods
Pebbles are initially disposed on a regular lattice at a height z = 8 cm and released with no initial velocity at t = 0.At the same time, the fluid is set into motion and pebbles start to move, driven by both gravity and the drag force.We use the classical numerical methods of molecular dynamics to compute the positions (r) and rotational velocities ( ) of the pebbles as a function of time: at each time step, all forces acting on each pebble are computed, and Newton's equations of motion (both translational and rotational) are integrated simultaneously for all pebbles by the Verlet method, of fourth order (Cundall and Strack, 1979;Pöschel and Schwager, 2005).The time step used in the simulation is t = 10 −6 s = T coll /100, which ensures that the trajectories during a collision are computed with sufficient accuracy.The "instantaneous" sediment flux q s (t) is computed over temporal windows of duration δt = 100 ms: www.earth-surf-dynam.net/4/327/2016/Earth Surf.Dynam., 4, 327-342, 2016 Within the alluvial cover, some pebbles are almost immobile, either because they got trapped by the bedrock roughness or because they belong to bottom layers of the cover and are therefore not entrained by the water flow.These pebbles constitute a static cover that contributes to protecting the bedrock from rapid impacts by saltating pebbles.We quantify the cover fraction in the following way: the bedrock surface is divided into square cells of side 2R.At each time step, we compute the velocity distribution of pebbles.If a pebble centred in a given cell has a velocity lower than 1/10th of the maximum velocity, this cell is considered "covered" by an immobile pebble.If n cells are covered at a time t, the time-averaged static cover fraction is then defined as With this definition, C = 1 when one layer of immobile particles is completely shielding the bedrock.
As are collisions between mobile particles, each impact of a pebble on the bedrock is inelastic: the impacting pebble loses a fraction of its incipient kinetic energy during the collision (due to the dissipative term in Eq. 8).This energy loss can result in the erosion of a small volume of bedrock (see Sect. 4.4).We can evaluate the energy lost during an impact by computing the work of the repulsive force during the collision, that is If we consider all impacts on bedrock occurring over a duration T , the total energy delivered to the bedrock by unit time and surface can then be expressed as 3 Results

Sediment transport
Let us first investigate the structure and dynamics of the bedload layer.Figure 3 shows the evolution of the flux of sediment q s with time when the Shields number is beyond the threshold of motion.After the pebbles are released in the flow, the bedload flux increases regularly during a transient phase.The duration of this transient phase depends on the Shields number but is always of the order of a few seconds.The bedload flux then reaches a relatively steady value: from this moment we consider that the system is in the permanent regime.Let us note that we still observe relatively large fluctuations of the bedload flux, consistent with what has been observed experimentally (see, for instance, Böhm et al., 2004, andAncey et al., 2006), and which motivated stochastic approaches to bedload transport on top of a sediment layer ( Roseberry et al., 2012;Ancey and Heyman, 2014;Fan et al., 2014).However, we will not investigate these fluctuations in further detail and in the rest of the discussion all results are computed in the permanent regime and only concern the average values of the sediment and energy fluxes.
Transport of sediment only occurs if the fluid drag force on a pebble is large enough to overcome solid friction, that is, if the Shields number exceeds a critical value c .In Fig. 4 we plot the flux of sediment, averaged over time (in the permanent regime) Q s = q s (t) as a function of .We observe that the threshold of motion corresponds to a critical value of the Shields number c = 0.012.Below this threshold, the average sediment flux vanishes after a short transient.The exact value of the threshold is somewhat difficult to assess, since some transport can occur with intermittency even below c .Within the margin of error of this definition, the threshold c depends only very slightly on the sediment supply σ and we shall assume in the following that c is a constant.The value found in our simulations is relatively low compared to many experimental observations (see, for example, the data compiled in the Fig. 1 of Lamb et al., 2008).However, let us note that our particles are perfectly spherical and therefore easy to put into motion, plus most observations of the transport threshold concern the incipient motion of particles over a thick sediment layer, where isolated moving particles are more likely to get trapped and stop their motion at low Shields number.
In Fig. 5, we plot the flow velocity, the mean velocity of pebbles and the solid volume fraction as a function of height and for a relatively large sediment supply (σ = 1.6).The local volume fraction is computed in horizontal slices of height R/3.5.Its profile presents two local maxima and vanishes for z ∼ 2R, which shows that the bedload is structured (in this example) into two rather compact layers.As evidenced by the water velocity profile, the flow only penetrates the upper layer of pebbles (which in this case is incomplete   and Müller, 1948;Métivier and Meunier, 2003) is plotted for comparison.By analogy with most sediment transport laws, which give the flux of sediment at saturation, we fit the evolution of Q s with ( − c ) by a power law, whose prefactor depends on the amount of sediment available, that is, on the sediment density σ : It has to be noted, though, that we only explore a limited range of Shields number, considering that the sediment supply σ is the main control parameter in this study, which does not guarantee a high precision in the determination of the index n.The index of the best fit roughly increases with the sediment supply, and we have n(σ ) < 1 for σ < 1 and n(σ ) > 1 otherwise.If the sediment supply is low, the sediment flux increases slowly with : pebbles can be transported at a higher speed when the flow accelerates, but the amount of available pebbles remains below the transport capacity.If the sediment supply is high enough, a rapid flow is able to put more pebbles into motion, which leads to a rapid variation in Q s with .
Let us now focus on the effect of the sediment supply on the bedload flux.In Fig. 7a, we plot the flux of sediment as a function of the sediment supply for a few values of the Shields number.For a given stream velocity, the amount of available pebbles and therefore the sediment flux both increase with σ .However, if the total number of pebbles is too high, the transport capacity of the flow is reached and the flux of sediment saturates.We remark that, in most cases, a local minimum in the bedload flux is reached around σ = 1, which can be seen both as a geometrical effect (mobile pebbles can then form a compact layer, consolidated by the bedrock roughness, and become hard to dislodge) and as an artificial effect of our fluid model: if the volume fraction in the only bottom layer is low, the fluid velocity is larger than zero; it vanishes as soon as the first layer is dense enough.If the sediment supply slightly increases, some pebbles will pop up above the first layer and be easily entrained by the flow.Similar but less pronounced local minima can also be observed around σ = 2 and σ = 3. Figure 7b confirms that the function f (σ ) = Q s ( − c ) n(σ ) , though not monotonous, is indeed independent of .The variation in Q s with σ shows that the sediment flux first increases linearly before saturating beyond a critical sediment supply σ 0 .For the sake of simplicity and in order to evaluate this critical value, we dismiss the local minima reached at σ = 1, 2, and 3 and fit the curve Q s (σ ) by a simple exponential function: Given this model, the surface density 3σ 0 corresponds to the maximum quantity of sediment that can be transported by a flow of given Shields number.The variation in σ 0 with is plotted in Fig. 7c: it can be well fitted by the affine function This finding is once again consistent with most transport models (see, for instance, the review by Lajeunesse et al., 2010) in the limit of high sediment supply: if the amount of sediment is large enough, the number of particles put into motion increases linearly with , whereas their velocity is proportional to the shear velocity, that is, to 0.5 .Finally in Fig. 8 we plot the maximum sediment transport rate Q t ( ), reached for large values of the sediment supply σ σ 0 .The best fit by a power law suggests Q t ∼ ( − c ) 1.2 .However, let us note that we do not have many data points, and that modelling spherical particles is likely to enhance transport close to the threshold.As shown in the same fig- ure, our data could also be consistent with the power law Q t ∼ ( − c ) 1.5 .

Static cover
In Fig. 9, we plot the evolution of the static cover fraction C as a function of σ for different values of the Shields number.As expected, below the threshold of motion, the static cover fraction first increases linearly with the sediment supply.Because of the roughness of the bedrock, it departs from the function C = σ beyond σ 0.5: this is due to the fact that the distribution of pebbles on the surface is not strictly homogeneous.A fraction of the bedrock is then covered by two layers of immobile pebbles, while other areas are bare.If exceeds the threshold of motion, the immobile cover fraction is very low for σ < 1 (and strictly zero when the Shields number is high enough, all particles being entrained

Incision process
The flux of energy that is delivered to the bedrock by the impacts is given by the work of the dissipative normal force during each collision between a mobile pebble and the bedrock (whether it is the flat surface or one of the glued spheres).In Fig. 10, we plot the variation in this flux of energy E with for different values of the sediment supply.As can be expected, E becomes larger when the velocity of the stream increases: if pebbles move at a higher speed in the bedload layer, their incoming velocity at the impact on the bedrock increases, as does the energy dissipated during the impact.As illustrated in Fig. 10a, the variation in E with can be fitted by a power law whose index varies with the sediment supply: As shown in Fig. 10b, the exponent m (σ ) increases roughly linearly with σ , with m(0) 1, and more generally with 1 ≤ m(σ ) ≤ 4.
Let us now investigate more precisely the effect of the sediment supply on the energy transfer, which is plotted in Fig. 11 for a few values of the Shields number.The shape of the curve is similar in all cases: E first increases with σ until it reaches a maximum value, and then decays to zero for large sediment supplies.The estimate of the incision rate corresponding to a given flux of energy delivered to the bedrock is obtained though the procedure described in Sect.4.4.The inset in Fig. 11 shows the variation in the maximum flux of energy with the Shields number.This variation can be well fitted by an affine relationship.This demonstrates that the process of incision only happens if the Shields number exceeds an incision threshold i 0.025 > c .Therefore, sediment transport can occur on a bedrock while not contributing to river incision if c < < i : in this regime, pebbles are rolling along the bedrock without impacting, and therefore do not contribute significantly to its erosion.

Influence of bedrock roughness
The roughness of the bedrock can be modified by varying two parameters: the surface density χ of fixed spheres on the bedrock and the protruding height of those pebbles over the horizontal bedrock (roughness height h b = R+z r ). Figure 12 illustrates the aspect of the bedrock, viewed from above, for two different values of χ .
In order to study the variation in energy transfer with the bedrock roughness, we plot the flux of energy delivered to the bedrock with respect to the sediment supply (σ ) for a single value of − c = 0.071, and for different roughness configurations.Figure 13a shows the same plot for increasing www.earth-surf-dynam.net/4/327/2016/Earth Surf.Dynam., 4, 327-342, 2016 Figure 11.Flux of energy delivered to the bedrock as a function of the dimensionless sediment supply.The energy transfer increases when the sediment supply increases in the range 0 < σ < 0.5 ("tool effect").It then decays for larger sediment supplies ("cover effect") and vanishes when the sediment supply reaches σ ≈ 3. The maximum flux is reached for a critical value σ m of the sediment supply, that depends only slightly on the Shields number.Plain lines are the best empirical fits by Eq. ( 27).The equivalent incision rate (rightside axis) is computed using Eq. ( 34) (see Sect. 4.4).The inset plots the maximum value of the energy flux as a function of − c .The dotted line is the best affine fine, which reveals the existence of an incision threshold i > c .
values of χ , and a roughness height h b = 4 cm.The aspect of the erosion curve is the same whatever the roughness density.However, a higher density of protruding spheres systematically leads to a decay of the energy received by the bedrock at high sediment supplies: the cover effect is more efficient if the bed possesses a dense roughness.In Fig. 13b, we plot the same evolution for a given roughness density χ = 0.36 but for spheres protruding more or less within the flow.At high sediment supply, the energy received by the bedrock appears to be higher if the bedrock is smooth.

Sediment transport rate
Although the model that we adopt for the interaction between the water flow and the pebbles is rather simple, the dynamics of the bedload layer appears to be consistent with experimental observations: as shown by Fig. 6, the order of magnitude of the sediment transport rate is comparable to the one predicted by the Meyer-Peter-Müller law.When the transport rate is fitted by a power law, the exponent that we obtain depends on the dimensionless sediment supply σ and is in general smaller than 1.5.At low Shields number, the sediment flux is higher for low supplies: indeed, it is easier for pebbles to roll along a flat bedrock than on the rough surface of a sediment layer.When is increased, the bedload flux increases faster for large sediment supply.For low sup-plies, once all pebbles have been put into motion, only their velocity can increase with .If more pebbles are available, an increase in leads to both more pebbles moving and an increase in their velocity, which implies n(σ ) > 1.The main discrepancy between our results and experimental observations is the value of the motion threshold, which is much lower ( c 0.012) in our simulations than in most reported measurements (Lamb et al., 2008).This is likely due to the fact that we model pebbles as spheres, which can easily roll along the bedrock.Let us note that Duran et al. (2012), who also used spheres but modelled more sophisticated fluidgrain interactions, obtained a critical Shields number of 0.12.However, like in most experimental studies, they were interested in the threshold of motion for grains rolling along the surface of a thick sediment layer, rather than on a relatively smooth surface, as is the case in our simulations.It would be necessary to consider more irregular pebbles (for instance by considering cohesive clusters of several spherical particles) in order to perform more realistic simulations.Once rescaled by the value of the threshold ( c = 0.012), the sediment flux that we obtain is, however, quantitatively consistent with common measurements (Meyer-Peter and Müller, 1948).

The role of sediment supply, Shields number, and bed roughness on incision
As shown in Fig. 11, the results of our simulations are qualitatively consistent with experimental observations by Sklar and Dietrich (2001): for a given Shields number , energy dissipated in the bedrock first increases with the sediment supply and reaches a maximum for σ m 0.5 before decaying, and vanishing at high sediment supply.This can be understood as the result of a competition between a "tool effect" and a "cover effect": as long as the bedrock is still exposed to impacts, the more pebbles are put into motion, the more energy they provide to the bedrock.When the sediment supply increases, immobile (or slowly rolling) pebbles start to accumulate on the bedrock, thus partially protecting it from direct impacts by saltating pebbles.The total energy transferred to the bedrock vanishes entirely beyond σ 3: at this point the bedrock is completely protected by the bedload layer.In Fig. 14, we compare more quantitatively our numerical results to the experiments by Sklar and Dietrich (2001) and the predictions of the linear (Sklar and Dietrich, 2004) and exponential (Turowski et al., 2007) cover models.Since these results were originally expressed in terms of the mass m s of gravel of diameter d, within a cylindrical container of diameter D, we derive the corresponding dimensionless sediment supply through the expression The comparison shows that our simulations are able to predict the right tendency for the flux of energy delivered to the Earth Surf.Dynam., 4, 327-342, 2016 www.earth-surf-dynam.net/4/327/2016/bedrock, that is, for the incision rate.However, let us note that because of the cylindrical geometry of the experiment by Sklar and Dietrich (2001), the fluid shear stress and the cover fraction are not uniform in their setup: immobile particles tend to accumulate in the centre of the disk, while saltating grains can still impact the bedrock around this protected area.This discrepancy might explain why the value of σ corresponding to maximum incision differs slightly between our simulations and the experiments.The quantitative description of the cover effect will be discussed further in Sect.4.3.Mudstone (Sklar & Dietrich, 2001) Limestone (Sklar & Dietrich, 2001) Andesite (Sklar & Dietrich, 2001) Linear cover model (Sklar & Dietrich, 2004) Exponential cover model (Turowski et al., 2007) Figure 14.We plot on the same graph our numerical prediction for the flux of energy delivered to the bedrock (full circles), erosion rates measured experimentally by Sklar and Dietrich (2001) (empty symbols), and the best fit for these experimental values by the linear cover model of Sklar and Dietrich (2004) and the exponential cover model of Turowski et al. (2007) (plain line).The sediment mass used in the experiments is converted into a dimensionless sediment supply using Eq. ( 26).The scale for E has been chosen so that its maximum coincides with the experimental maximum erosion rate.
We also quantified the influence of the Shields number on the abrasion process and showed a power-law dependency of E on the excess shear stress − c .In particular, our fits (see Fig. 10b) show that the exponent m is always greater than the index n of the transport law.This implies that if Q s and σ are kept constant, the incision rate increases with the Shields number.This result contradicts the prediction of the saltation-abrasion model, where the incision rate scales like ( − c ) −0.5 (Sklar and Dietrich, 2004), and findings by Chatanantavet and Parker (2009) and Johnson and Whipple (2010), where there is no explicit dependency of E in (for given values of Q s and σ ).However, it is consistent with the prediction of the shear stress and stream-power incision models (Whipple and Tucker, 1999).The variation in E with is nonlinear, with an exponent m > 1, which means that for a given sediment load, the effect of large hydraulic events in the long term will be amplified.However, let us remark that such events can also result in a sudden increase in the sedi- ment load, which could result in total inhibition of erosion, as evidenced by Lague (2010).
Our results also show that abrasion only occurs beyond a given threshold which is higher than the threshold of motion of pebbles, which can be explained by the fact that rolling or sliding pebbles do not contribute significantly to the erosion of bedrock.This is inconsistent, however, with observations by Sklar and Dietrich (2001), who report abrasion occurring as soon as the flow is able to put sediment into motion.Let us note, however, that we vary the flow velocity for a given pebble size, whereas Sklar and Dietrich (2001) vary the sediment size for a given velocity.The fact that there is only a small difference between the two values c and i could explain that the discrepancy was not observed experimentally.If the existence of an incision threshold i > c were to be confirmed, it would mean that this value should be taken into account in models of river incision instead of the critical Shields number c .This is also consistent with the possible existence of an energy threshold necessary to effectively erode a small volume of material at impact, as observed by Bitter (1963).
Our simulations indicate that the roughness of the bedrock does not affect the general evolution of the energy delivered to the bed with respect to the cover fraction.However, incision appears to be enhanced if the surface density of asperities is low and if they are not too high.These two effects can be related to the geometrical explanation of the cover effect: if the roughness is denser or higher, mobile pebbles are more likely to get trapped and immobilized along the bedrock, therefore protecting it from further impacts by rapid pebbles.This enhanced cover effect will disappear if the roughness density χ is too large: indeed, if the bedrock was entirely covered with glued spheres, it would become equivalent to a smooth bedrock.Let us note that recently Huda and Small (2014) modified the saltation-abrasion model in order to take into account bedrock roughness, and found the opposite result: the incision rate is considerably increased (by more than 1 order of magnitude) by the presence of long-scale bed topography.However, the roughness that we implement in our model does not modify the local slope, and has a length scale comparable to the pebble size.
Finally, let us remark that the influence of the coefficient of restitution on the results of our simulations should be of importance and will be the object of further investigation.Increasing the coefficient of restitution would certainly facilitate the saltating motion of pebbles, whereas they only roll along the bedrock at low e.Increasing e could therefore decrease the incision threshold by narrowing the rolling/sliding regime.In addition, a high coefficient of restitution means that a lower fraction of the kinetic energy of the projectile is delivered to the bedrock.However, this implies that the impactor rebounds with a higher kinetic energy, and is then more likely to impact again the bedrock at high speed.It is therefore not trivial to assess in which way the total energy delivered to the bedrock (that is, the number of impacts multiplied by the energy given at each impact) will evolve with the value of e.

Cover effect
By analogy with both the linear (Sklar and Dietrich, 2004) and the exponential (Turowski et al., 2007) cover models, we first isolate in E the influence of the sediment supply, and fit the flux of energy by the empirical function E ( , σ ) = ( ) σ p e −σ/γ . ( The values for p and γ , corresponding to the best fits plotted in Fig. 11, are reported in Table 2.The parameter γ appears to be roughly independent of the Shields number, and we consistently find p > 2, whereas p = 1 in the exponential cover model.The fact that p depends on underlines that we cannot express the incision rate as a simple product of a function of and a function of σ .Furthermore, in our simulations, forcing a fit with p = 1 always leads to underestimate the maximum incision rate and overestimate it at high sediment supply.This is likely due to the fact that, in the saltation-abrasion model, the number of impacts is proportional to the sediment supply, whereas our simulations show that the sediment flux does not vary linearly in σ .Therefore, the incision rate increases faster than linearly at low sediment supply. Following the approach for incision rate by Sklar and Dietrich (2004) and Turowski et al. (2007) (see Eq. 1) and for the impulse rate by Turowski and Rickenmann (2009), let us express the flux of energy delivered to the bedrock, defined by Eq. ( 20), as the product of the energy provided by each impact (E i ), the number of impacts per unit time and surface (n i ), and a cover function (F ): Earth Surf.Dynam., 4, 327-342, 2016 www.earth-surf-dynam.net/4/327/2016/In the saltation-abrasion model, n i is proportional to the sediment supply while the cover function F is interpreted as the probability that an impact hits the bedrock.In the original model by Sklar and Dietrich (2004), F is simply the fraction of exposed bedrock and decays linearly with the relative sediment flux Q s /Q t .In the stochastic model of Turowski et al. (2007), F is expressed as an exponential function of Q s /Q t , leading to a prediction closer to experimental observations (see Fig. 14).As noted in the Introduction, this equation implies that the number of impacts and the cover function F are independent parameters.However, if we take into account the fact that mobile pebbles are not isolated from one another, nor from the static cover, they should be interdependent.Let us now reformulate Eq. ( 28) in order to take into account only the Shields number and the sediment flux, which are easier to compute, or to measure experimentally, than the frequency of impacts and their energy.Indeed, if a saltating particle hits a mobile particle shielding the bedrock, a small fraction of its incipient energy could still be transmitted to the bedrock (Turowski and Bloem, 2015), though this event would not be counted as an eroding impact with Eq. (1).From the simulations we can then compute a cover function F (σ, ) that does not require a geometrical or stochastic description of the alluvial cover of the bedrock.The energy of an impact E i is expected to scale like the typical kinetic energy of moving pebbles, which is itself proportional to U * 2 (Fan et al., 2014).Using Eq. ( 13), and assuming that the impact energy vanishes for < c , we can therefore write, dimensionally, As proposed by Foley (1980), the number of impacts per unit time and surface is expected to be proportional to the sediment transport rate Q s , which can be written dimensionally as Therefore, the flux of energy should read with K a dimensionless constant, which depends neither on nor on σ .This allows us to redefine the cover function as In Fig. 15 we plot F in a log-linear scale.Let us first observe that K is indeed independent of , since all curves converge to F = 1 for σ = 0, which implies that K = 1 and validates our dimensional analysis.On the one hand, as emphasized in the inset of Fig. 15, the cover function F decays slower than exponentially, and rather linearly, for σ 0.75.On the other hand, F can be well fitted by an exponential decay for σ 1.This behaviour is similar to the case that Hodge and Hoey (2012) refer to as "sigmoidal" in their cellular automaton model (though as a function of Q s /Q t and not σ , see their Fig. 7, and although they only consider static cover), and to some observations by Chatanantavet and Parker (2008), where the cover function appears to decay below 1 only for Q s /Q t 0.25 − 0.75 (see their Fig.13).Our result therefore confirms that the "exponential cover model" overestimates the cover effect at low sediment supplies but fits correctly at high enough sediment supplies.In the exponential regime, following the stochastic approach by Turowski et al. (2007), we can fit the cover function by with 3σ 0 an estimate of the normalized sediment mass transport capacity (see Fig. 11) and ϕ the "cover factor" of the probabilistic approach (Turowski et al., 2007).Our results show systematically ϕ > 1 (see Table 3), which implies that it is more probable for a pebble to impact on an uncovered zone of the bedrock than a covered one.This is consistent with observations by Chatanantavet and Parker (2008) (in flume experiments) and Turowski and Rickenmann (2009) (in the field).In Fig. 15, we also plot the function 1 − C(σ ), where C is the static cover fraction computed in Sect.3.2.If the evolution of both functions with σ is similar, the cover function F is systematically smaller than 1 − C: this implies that the cover effect is due to not only immobile but also mobile (rolling or saltating) pebbles.The latter can either directly shield the underlying bedrock (which was referred to as "dynamic cover effect" by Turowski et al., 2007) or hit other saltating pebbles and slow them down.

Estimation of the incision rate
We can estimate the rate of incision induced by the impacts on the bedrock, based on the flux of energy delivered.Fol-G.Aubert et al.: Numerical study of bedrock incision by the bedload layer lowing Engle (1978) and Sklar and Dietrich (2004), we express the incision rate as where v is the energy required to incise a unit volume of rock.The value that is used in Sklar and Dietrich (2004) is derived from the mechanical properties of rocks: σ T = 7 × 10 6 Pa is the rock tensile strength, Y = 5.0 × 10 10 Pa is the rock elastic modulus, and k v = 10 6 is a dimensionless rock resistance parameter.The incision rate obtained is plotted as a function of the dimensionless sediment supply and for different values of the Shields number in Fig. 11.
For instance, the incision rate corresponding to the value E = 40 W m −2 (obtained at − c = 0.049 and σ = 0.6) is I = 2.6 m yr −1 = 7 mm day −1 .Let us note that this value corresponds to an instantaneous incision rate and not to the average incision rate over a year: it may be reached for a high water discharge and for a particular value of the sediment supply.In a river, these conditions may be verified only during a few days per year, while the instantaneous incision rate would be very low the rest of the time, when the discharge is small, and the sediment supply is either very low or very high.The order of magnitude of the incision rate that we predict is comparable to values measured in rapidly eroding rivers in Taiwan during storm events (Hartshorn et al., 2002).As we have predicted the value of the instantaneous incision rate for a wide range of both Shields number and sediment supply, it would be possible to compute the longterm average incision rate for a given stream, provided that the probability distribution functions of both water discharge and sediment supply are known.

Conclusions
We have presented the results of a new model for incision of a river bedrock based on the direct simulation of physically based trajectories of pebbles in a stream.In this model we solved the equations of motion for a large number of pebbles entrained by a turbulent water flow, with a simplified retroaction of the presence of the pebbles on the flow.This allowed us to explicitly compute the trajectories of pebbles transported by the flow, and therefore to quantify the energy dissipated during collisions between the bedload and the bedrock, which is directly responsible for the incision of the bedrock.We found that the sediment transport rate can be fitted by a power law of the Shields number, similar to most classical transport laws at saturation.However, we also evidenced the influence of the sediment supply: the exponent of the transport law increases with the quantity of available pebbles.For a given Shields number, we showed that the bedload flux increases with the sediment supply until it reaches its saturated value.This allowed us to compute the sediment mass that the flow is able to transport.However, extracting a unique general expression for the flux of sediment as a function of both the Shields number and the sediment supply remains non-trivial.
The amount of energy that impacts of saltating pebbles deliver to the bedrock can be directly computed from the simulation data.This flux of energy, which is expected to be proportional to the incision rate, shows the same qualitative variations with sediment supply as observed in experiments by Sklar and Dietrich (2001): it first increases with the amount of sediment available (as the number of impacts increases) before decaying when there is too much sediment and the bedrock becomes shielded.We also showed that the energy delivered to the bedrock increases as a power law of the Shields number, and is zero below a given incision threshold, higher than the motion threshold, which was not observed in experiments.Finally, by extracting a cover function from our data, we showed that the classical linear and exponential models for the cover effect lead to underestimation of the incision rate, at high and low sediment supplies respectively.The shape of our cover function resembles experimental observations by Chatanantavet and Parker (2008) and some numerical results by Hodge and Hoey (2012).If defined as in Eq. ( 33), the cover function appears instead to decay linearly at low sediment supply and exponentially at high sediment supply.This underlines the fact that the amount of sediment available contributes not only to shield the bedrock but also to change the dynamics of saltating particles.Finally, we evaluated the rate of incision predicted by our simulations as a function of the hydraulic conditions (the Shields number) and the amount of sediment available (dimensionless sediment supply).Its order of magnitude appears to be consistent with long-term observations made in mountain streams.
Though our results are qualitatively consistent with experimental observations (Sklar and Dietrich, 2001) and another type of model (Sklar and Dietrich, 2004;Turowski et al., 2007), the quantitative aspect is probably affected by our numerical method: indeed, the fact that we model the flow by a horizontally averaged and purely horizontal velocity profile is likely, on the one hand, to have a (negative) impact on the possibility for pebbles to gather into immobile patches, and therefore on the efficiency of the cover effect.On the other hand, taking into account turbulent velocity fluctuations, and in particular local bursts of vertical velocity, could enhance saltating motion.A better explicit model of the dynamics of the alluvial cover would therefore require accounting for a spatially non-uniform velocity field, and ideally the exact velocity field around each mobile particle, which would be much more time-consuming numerically.Such stochastic effects are probably better accounted for in models such as the cellular automaton by Hodge and Hoey (2012).The assumption of spherical particles is also likely to have an effect on our predictions: indeed, angular or irregularly shaped sediments would probably be less easily put into motion by the flow (for instance, a flat pebble would be harder to dislodge from the bedrock; Rust, 1972;Komar and Li, 1986).In con-trast, once entrained by the flow they could slide along the bedrock (instead of simply rolling) and therefore contribute to its wear via solid friction.
Finally, let us note that in the prediction of the long-term evolution of a river bed (see, for example, Lague, 2010), incision of the bedrock is not the only relevant parameter.Our numerical approach would also be relevant for the study of the incision of lateral walls (if we add sidewalls to the computational domain) and the comminution of mobile particles (since our simulations also give us access to the energy lost by mobiles pebbles, not only in their impacts with the bedrock but also in all their contacts with one another).

Figure 1 .
Figure 1.Snapshot of the numerical simulation.The colour scale codes the horizontal velocity of each pebble.The grey plane and (immobile) grey spheres constitute the rough bedrock.
with g = −g e z the gravitational acceleration (g = 9.8 m s −2 ).

G
. Aubert et al.: Numerical study of bedrock incision by the bedload layer

Figure 2 .
Figure 2. Two pebbles in contact, located respectively at r i and r j , with an overlap δ = 2R − |r i − r j |. v ij and ij are, respectively, the translational and angular relative velocities of pebble j with respect to pebble i.Normal N ij and tangential T ij forces apply at the contact.

Figure 3 .Figure 3 .
Figure3.Flux of sediment q(t) as a function of time, for a Shields number Θ = 0.083 and a sediment supply σ = 0.6.A transient state is observed for a few seconds before the steady state is reached.

Figure 4 .Figure 5 .
Figure 4. Average flux of sediment Q s = q s as a function of the Shields number and for different values of the sediment supply σ .Q s becomes significantly larger than 0 when the Shields number exceeds c ≈ 0.012.

Figure 6 .
Figure 6.Flux of sediment as a function of the relative excess shear stress ( − c )/ c for different values of the sediment supply, in a log-log scale.The critical Shields number is c = 0.012.Plain lines are best fits by a power law (see Eq. 22).The bold line represents the Meyer-Peter-Müller transport law.Inset: exponent of the best fit by a power law, and its standard deviation, as a function of the sediment supply σ .The dotted line represents n = 1.5.

Figure 7 Figure 8 .Figure 7 .
Figure 7. a) Flux of sediment as a function of the sediment supply for different values of the Shields number.The flux of sediment increases linearly for low sediment supplies, but tends to a saturated value when transport capacity is reached.Plain lines are best fits by the equation Qs (Θ, σ) = Qt (Θ)×(1−e −σ/σ 0 ).b) The function f (σ) (as defined by equation (22)) is shown to be roughly independent of the excess shear stress Θ ′ = Θ − Θc. c) Evolution of the critical sediment supply σ0 with Θ − Θc.The dashed line is the best linear fit.

Figure 8 .
Figure8.Saturated value Q t of the sediment transport rate as a function of ( − c )/ c .The plain line is the best fit by a power law; the dashed line is the fit by a power law of index n = 1.5.

Figure 9 .
Figure 9. Static cover fraction as a function of the sediment supply, for different values of the Shields number.The dashed line represents the function C = σ .

Figure 10 .
Figure 10.(a) Flux of energy received by the bedrock as a function of the relative excess shear stress, in a log-log scale.Plain lines are best fits of E with power laws (see Eq. 25).(b) Exponent m(σ ) of the power law as a function of the dimensionless sediment supply σ .Data points are fitted by a linear function.

Figure 12 .
Figure 12.Two different cases of bedrock roughness: the positions of the glued spheres are plotted as seen from above.Left: χ = 0.08.Right: χ = 0.52.The flow is in the x direction, from left to right.

Figure 13 .
Figure 13.Flux of energy delivered to the bedrock as a function of the sediment supply for − c = 0.071: (a) for increasing values of the roughness density χ and (b) for different values of the roughness height.

FFigure 15 .
Figure 15.Empty symbols: cover function F , as defined by Eq. (32), plotted as a function of σ , in a log-linear scale.Plain lines are best fits by exponential tails for σ > 1. Full symbols plot 1 − C, with C the static cover fraction as computed in Sect.3.2.Inset: Zoom on the cover function for low values of σ , in a linear scale.

Table 1 .
List of the physical parameters used in the model.
L length of the box (m) W width of the box (m) H height of the box (m) N number of mobile pebbles R radius of pebbles (m) ρ s pebble density (kg m −3 ) g gravitational acceleration (m s −2 ) k elastic constant of collisions (kg s −2 ) effective viscosity of collisions (Pa s)

Table 2 .
Coefficients used in the empirical fit of the flux of energy as a function of sediment supply (see Eq. 27).

Table 3 .
Cover factor ϕ extracted from the exponential tail of the cover function F (see Eq. 33).