Bedrock incision by bedload: insights from direct numerical simulations

Abstract

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 field, 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 field. Finally, we reformulate the dependency of this incision rate with Shields number and sediment flux, and predict that the cover term should decay linearly at low sediment supply and exponentially at high sediment supply.

Earth Surf. Dynam., 4, 327342, 2016 www.earth-surf-dynam.net/4/327/2016/ doi:10.5194/esurf-4-327-2016 Author(s) 2016. CC Attribution 3.0 License.

Guilhem Aubert, Vincent J. Langlois, and Pascal Allemand

Laboratoire de Gologie de Lyon, Universit Claude Bernard Lyon 1/ENS de Lyon/CNRS UMR5276, Villeurbanne, France

Correspondence to: Vincent J. Langlois ([email protected])

Received: 21 December 2015 Published in Earth Surf. Dynam. Discuss.: 19 January 2016 Revised: 21 March 2016 Accepted: 5 April 2016 Published: 20 April 2016

Abstract. 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.

1 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 andTucker, 1999). Long-term averaged incision rates can take values from 0.02 to 14 mm yr1 (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 ow, 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.

Bedrock incision by bedload: insights from direct numerical simulations

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-dened, since the number of impacts ni 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 inuence of bedrock roughness. Finally, in Sect. 4 we discuss the implications of our results on the inuence of both the Shields number and the sediment supply on the incision rate. We propose a new denition 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.

2 Description of the numerical model

2.1 Numerical setup

We use the discrete-element method to simulate the individual dynamics of pebbles entrained by a turbulent water ow over a xed 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 = 2m, width W = 1m, and height H = 2m (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 Nb spheres of radius R = 5cm,

centred at a height z = zr. These protruding spheres have the

same mechanical properties as the pebbles entrained by the ow but are xed and considered part of the bedrock. In all the presented results except in Sect. 3.4, these spheres protrude by a height hb = R+zr = 4cm and their surface density

Earth Surf. Dynam., 4, 327342, 2016 www.earth-surf-dynam.net/4/327/2016/

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

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 conrms that, as proposed early on by Gilbert (1877) or Shepherd (1972), the amount of sediment available in the river channel should inuence the downcutting rate in two opposite ways: incision should be rst 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 rst 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 ow, they conrmed 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 saltation abrasion model. In this model, the incision rate I is written as the product of three terms: the volume of rock eroded by each impact, Vi; the number of impacts by unit time and surface, ni; and the probability that a saltating grain impacts an exposed area of the bedrock, F :

I = Vi [notdef] ni [notdef] 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 Qs when below the transport capacity of the stream, and vanishes when the transport capacity Qt is reached:

F =

1 Qs/Qt for Qs < Qt

0 otherwise . (2)

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 :

F = exp(' Qs/Qt), (3) 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

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

Flow direction

z

y

x

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.

is given by

[notdef] = Nb [notdef]

R2

W L = 0.36. (4)

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 ow, 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 signicant abrasion of the bedrock only happens over at least the duration of a ood, that is, a few days, and in many cases over years. Therefore, we can admit that the bedrock (both the horizontal surface and the xed spheres) remains unchanged and immobile within the timescale of the simulation (in particular, its altitude remains z = 0 through

out 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 = 5cm and density s = 2500kg m3.

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 WL/(R2). Therefore, we quantify the sediment supply by dening the dimensionless surface density as the surface of bedrock covered by pebbles, divided by the total available surface:

= N [notdef]

A horizontal turbulent water ow in the x direction puts the pebbles into motion (see Sect. 2.3). Their trajectories are then driven by their immersed weight (W), uid 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 [Omega1] of a pebble is given by Newtons equations of motion:

8 >

>

>

<

>

>

>

:

m d2rdt2 = W + F + N + T

J d[Omega1]

dt = R [notdef] T + M,

(6)

4

8

15 sR5 the angular

momentum of a pebble. The immersed weight of a pebble is

W = (s w)

with m =

3 sR3 the mass and J =

4R3

3 g, (7)

with g = g ez the gravitational acceleration

(g = 9.8ms2).

www.earth-surf-dynam.net/4/327/2016/ Earth Surf. Dynam., 4, 327342, 2016

R2

W L . (5)

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 m3)

g gravitational acceleration (m s2)

k elastic constant of collisions (kg s2)

[Gamma1] effective viscosity of collisions (Pa s)

e coefcient of restitution of collisions overlap between two pebbles (m)

[notdef] local friction coefcient normalized sediment supply0 movable sediment supply at a given Shields number m sediment supply for maximal incision rate solid volume fraction shear stress (Pa)

U shear velocity (m s1)

von Krmn constantw density of water (kg m3)

w dynamic viscosity of water (Pa s)

[Theta1] Shields number[Theta1]c critical Shields number for transport of sediment [Theta1][prime] = [Theta1] [Theta1]c dimensionless excess shear stress

[Theta1]i critical Shields number for incision of the bedrock Nb number of spheres glued on the bedrockhb height of the bedrock roughness (m)

[notdef] surface density of bedrock roughnessTcoll typical duration of a contact between 2 pebbles (s)

[Delta1]t time step in the simulations (s)

q(t) instantaneous sediment ux (kg m1 s1)

Qs average sediment ux (kg m1 s1) [Phi1]E ux of energy (W m2)

I incision rate (m s1)

kv dimensionless rock resistance coefcient Y Youngs modulus of rock (Pa)

[epsilon1]v energy required to erode a unit volume of rock (J m3) T tensile strength of rock (Pa)

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

2.2 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, Pschel 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 [notdef]ri rj [notdef] > 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

Nij =

k [Gamma1]d dt

[parenrightbigg]

nij , (8)

where k is the elastic constant, [Gamma1] is the effective viscosity,

and nij =

rj ri

|rj ri[notdef]

is the unit normal vector of the collision.

The value of the elastic constant is related to the materials Youngs modulus and the pebble size: we adopt the value k = 2 [notdef] 108 N m1, which corresponds to an elastic modu

lus Y k/R = 4GPa. 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 60nm. 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 ow. Within our model, the effective inelasticity of the collision can be quantied by the coefcient of restitution e, which compares the velocity of the pebble before and after a collision: e = 1 cor

responds to an elastic collision and e = 0 to total dissipation.

If the force is given by Eq. (8), e is expressed as

e = exp

Tcoll[Gamma1] 2m

[parenrightbigg]

Figure 2. Two pebbles in contact, located respectively at ri and rj , with an overlap = 2R [notdef]ri rj [notdef]. vij and [Omega1]ij are, respectively,

the translational and angular relative velocities of pebble j with respect to pebble i. Normal Nij and tangential T ij forces apply at the contact.

simulation, the energy loss is due to not only the viscous dissipation in the lm 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 Coulombs law of solid friction, as in Cundall and Strack (1979). This force opposes the tangential motion and is expressed as

T ij = min(G[notdef]vs[notdef];[notdef][notdef]Nij [notdef])

vs

|vs[notdef]

pm/k

p1 [Gamma1]2/(4mk)

,(9)

the typical duration of a collision. We choose the value [Gamma1] = 2 [notdef] 104 kg s1 for the effective viscosity, which leads to

collisions of duration Tcoll = 104 s and a coefcient of resti

tution e = 0.3, which means that a pebble loses 1e2 = 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

, (10)

where vs is the sliding velocity at the contact, which is a function of the two pebbles translational and angular velocities, [notdef] = 0.6 is the local friction coefcient, and G = 5kg s1 is

the slope of the regularization of the Coulombs 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 innite size, and same mechanical properties.

2.3 Turbulent water ow

A stationary turbulent ow over a rough bedrock follows the average velocity prole:

V (z) = Vf (z)ex with Vf (z) =

U

ln

with Tcoll =

[parenleftbigg]

z z0

[parenrightbigg]

(11)

Earth Surf. Dynam., 4, 327342, 2016 www.earth-surf-dynam.net/4/327/2016/

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

where x is the direction of the ow and = 0.41 is the von

Krmn constant. z0 is the bedrock roughness and depends on pebble size: a bedrock made of pebbles of radius R has a roughness z0 = 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:

U = [radicalbigg]

w with = w 2z2[parenleftbigg]

2, (12)

where w = 1000kg m3 is the density of water. The abil

ity 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:

[Theta1] =

[notdef] (2R)2

(s w)g [notdef] (2R)3 =

@Vf @z

R2

2 CD [notdef]U[notdef]U, (14)

where U = V (z)

drdt is the relative velocity between the lo

cal ow and the pebble. The drag coefcient CD 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:

CD =

24Rep + 0.4, (15)

where Rep =

2w UR

w , with w = 103 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 Rep. When a sphere is rotating in a viscous uid 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:

M = 8 wR3 [Omega1]. (16)

2.5 Computational methods

Pebbles are initially disposed on a regular lattice at a height z = 8cm and released with no initial velocity at t = 0. At the

same time, the uid 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 ([Omega1]) of the pebbles as a function of time: at each time step, all forces acting on each pebble are computed, and Newtons equations of motion (both translational and rotational) are integrated simultaneously for all pebbles by the Verlet method, of fourth order (Cundall and Strack, 1979; Pschel and Schwager, 2005).

The time step used in the simulation is [Delta1]t = 106 s =

Tcoll/100, which ensures that the trajectories during a collision are computed with sufcient accuracy. The instantaneous sediment ux qs(t) is computed over temporal windows of duration t = 100ms:

qs(t) =

1 t [notdef]

w(U )2

2(s w) [notdef] gR

. (13)

Pebbles are put into motion by the ow if [Theta1] exceeds a threshold value [Theta1]c: measurements of this threshold, both in the eld and in experiments, give values in the range 0.01 < [Theta1]c < 0.2 (Bufngton and Montgomery, 1997; Lamb et al., 2008). In abrasion experiments conducted in a ume setup by Attal and Lav (2009) the maximum uid velocity is 4m s1 for a water height H = 60cm and pebbles of

size 10 to 80mm. As reported in data reviewed by Ricken-mann and Recking (2011), ow velocity in mountain streams varies typically between 0.3 and 4m s1, for a water height between 0.1 and 3m. In order to be consistent with these observations, in our simulations we adopted mean water velocities up to 5.0m s1, which corresponds to [Theta1] varying from 0 to 0.11.

If the bedrock is covered with a layer of mobile pebbles, as in our simulations, the turbulent velocity prole is modied. Recently, Duran et al. (2012) developed a quasi-2-D mechanistic approach that takes into account the retroaction of the pebbles on the water ow by assuming the conservation of total horizontal momentum in horizontal slices. They showed that the uid 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 prole in the clear water region. The intermediate region where the velocity goes from zero to the logarithmic prole is very thin (usually of the order of one grain diameter). Therefore, we simplied the treatment of this retroaction by modelling only two different regions in the ow: at each time step, we compute the average solid fraction in horizontal slices. If < b = 0.5, the velocity prole is logarith

mic and not affected by the presence of pebbles. If [greaterorequalslant]b, the velocity of water vanishes: Vf = 0 (see the resulting velocity

prole in Fig. 5). This approximation would be too simplistic if we were to study the exact ux 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 uid ow but also the horizontal variations in pebble density.

2.4 Interactions between pebbles and the ow

The turbulent ow exerts on each mobile pebble a drag force given by

F = w

N

1 WL [notdef]

t+ t

[integraldisplay]

t

dt[prime]. (17)

www.earth-surf-dynam.net/4/327/2016/ Earth Surf. Dynam., 4, 327342, 2016

Xi=1 mvxi(t[prime])

[parenrightBigg]

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

Within the alluvial cover, some pebbles are almost immo-bile, 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 ow. 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 dened as

C = [angbracketleft]n(t)[angbracketright] [notdef]

70

60

11)

s = 0.6 Q = 0.083

Flux of sediment ( kg m s

50

40

30

20

Transient state

Steady state

10

0 0 2 4 6 8 10 Time (s)

Figure 3. Flux of sediment qs as a function of time, for a Shields number [Theta1] = 0.083 and a sediment supply = 0.6. A transient state

is observed for a few seconds before the steady state is reached.

(Roseberry et al., 2012; Ancey and Heyman, 2014; Fan et al., 2014). However, we will not investigate these uctuations 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 uxes.

Transport of sediment only occurs if the uid drag force on a pebble is large enough to overcome solid friction, that is, if the Shields number exceeds a critical value [Theta1]c. In Fig. 4 we plot the ux of sediment, averaged over time (in the permanent regime) Qs = [angbracketleft]qs(t)[angbracketright] as a function of [Theta1]. We observe that

the threshold of motion corresponds to a critical value of the Shields number [Theta1]c = 0.012. Below this threshold, the aver

age sediment ux vanishes after a short transient. The exact value of the threshold is somewhat difcult to assess, since some transport can occur with intermittency even below [Theta1]c.

Within the margin of error of this denition, the threshold [Theta1]c depends only very slightly on the sediment supply and we shall assume in the following that [Theta1]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 ow 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 lo

cal volume fraction is computed in horizontal slices of height R/3.5. Its prole 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 prole, the ow only penetrates the upper layer of pebbles (which in this case is incomplete). Most of the pebbles lying in the bottom layer are therefore completely

Earth Surf. Dynam., 4, 327342, 2016 www.earth-surf-dynam.net/4/327/2016/

4R2WL. (18)

With this denition, C = 1 when one layer of immobile par

ticles 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

[Delta1]E =

[integraldisplay]

collision

Figure 3. 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.

24

d dt (t[prime])ndt[prime]. (19)

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

[Phi1]E =

1 WL [notdef]

N t[prime]

[parenrightbig]

[notdef]

1T [notdef]

Ximpacts[Delta1]E. (20)

3 Results

3.1 Sediment transport

Let us rst investigate the structure and dynamics of the bed-load layer. Figure 3 shows the evolution of the ux of sediment qs with time when the Shields number is beyond the threshold of motion. After the pebbles are released in the ow, the bedload ux 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 ux 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 uctuations of the bedload ux, consistent with what has been observed experimentally (see, for instance, Bhm et al., 2004, and Ancey et al., 2006), and which motivated stochastic approaches to bedload transport on top of a sediment layer

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

0 0 0.01 0.02 0.03 0.04Shields number ()

60

Average sediment flux (kg ms)

-1-1

= 0.4 = 0.6 = 0.8 = 1.2 = 2.0 = 3.2

= 0.48 = 0.96 = 1.4 = 2.4 = 4.0

100

Flux of sediment ( kg m s)

11

40

1.8

c = 0.012 0.002

10

1.6

20

1.4

1.2

1

n()

0.8

0.6

0 1 2 3 4 5 6

1

1 10 (c)/c

Figure 4. Average ux of sediment Qs = [angbracketleft]qs[angbracketright] as a function of the

Shields number [Theta1] and for different values of the sediment supply . Qs becomes signicantly larger than 0 when the Shields number exceeds [Theta1]c 0.012.

0 0 0.5 1 1.5 2 2.5Velocity (m s )

1

Figure 6. Flux of sediment as a function of the relative excess shear stress ([Theta1] [Theta1]c)/[Theta1]c for different values of the sediment supply, in

a loglog scale. The critical Shields number is [Theta1]c = 0.012. Plain

lines are best ts by a power law (see Eq. 22). The bold line represents the Meyer-PeterMller transport law. Inset: exponent of the best t by a power law, and its standard deviation, as a function of the sediment supply . The dotted line represents n = 1.5.

whose prefactor depends on the amount of sediment available, that is, on the sediment density :

Qs (,[Theta1]) = f ()([Theta1] [Theta1]c)n(). (22)

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 t 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 ux increases slowly with [Theta1]: pebbles can be transported at a higher speed when the ow accelerates, but the amount of available pebbles remains below the transport capacity. If the sediment supply is high enough, a rapid ow is able to put more pebbles into motion, which leads to a rapid variation in Qs with [Theta1].

Let us now focus on the effect of the sediment supply on the bedload ux. In Fig. 7a, we plot the ux 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 ux both increase with . However, if the total number of pebbles is too high, the transport capacity of the ow is reached and the ux of sediment saturates. We remark that, in most cases, a local minimum in the bedload ux is reached around = 1, which can be seen both as a geometrical effect (mo

bile pebbles can then form a compact layer, consolidated by the bedrock roughness, and become hard to dislodge) and as an articial effect of our uid model: if the volume fraction in the only bottom layer is low, the uid velocity is larger than zero; it vanishes as soon as the rst layer is dense enough.

www.earth-surf-dynam.net/4/327/2016/ Earth Surf. Dynam., 4, 327342, 2016

Solid volume fraction

= 1.6 ; = 0.061

3

Flow velocityGrains volume fraction Grains mean velocity

z / (2R)

0 0.5 1 1.5 2 2.5

2

1

Figure 5. Solid volume fraction and velocity of the uid and the pebbles as a function of height for a Shields number [Theta1] = 0.061 and

a sediment supply = 1.6. Two distinct layers of pebbles can be

observed. The velocity of the ow and of pebbles vanishes in the bottom layer.

immobile or only slightly entrained by the upper mobile particles. The average velocity of pebbles exposed to the ow increases with their vertical position in the bedload, but remains lower than the velocity of water.

Let us now investigate the transport law by varying the Shields number [Theta1]. In Fig. 6, we plot the variation in the average ux of sediment Qs with the reduced Shields number ([Theta1][Theta1]c)/[Theta1]c, and for different values of the sediment supply.

The classically used Meyer-PeterMller law (Meyer-Peter and Mller, 1948; Mtivier and Meunier, 2003)

Qsat = 8

r ww g (2R)3 [notdef] ([Theta1] 0.047)3/2 (21)

is plotted for comparison. By analogy with most sediment transport laws, which give the ux of sediment at saturation, we t the evolution of Qs with ([Theta1] [Theta1]c) by a power law,

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

3

(b)

(c)

f(s) x 10 4

(a)

2

Q = 0.003 Q = 0.015 Q = 0.031 Q = 0.049 Q = 0.071

200

QQ = 0.015 QQ = 0.031 QQ = 0.049 QQ = 0.071

1

11

Flux of sediment ( kg m s)

150

0 0 1 2 3 4 5 6Sediment supply (s)

2

100

1.5

s0 = 20 (QQc)

50

s 0

1

0.5

0 0 1 2 3 4 5 6Sediment supply (s)

0 0 0.02 0.04 0.06 0.08 QQc

Figure 7. (a) Flux of sediment as a function of the sediment supply for different values of the Shields number. The ux of sediment increases linearly for low sediment supplies, but tends to a saturated value when transport capacity is reached. Solid lines are best ts by the equation Qs ([Theta1],) = Qt ([Theta1]) [notdef] (1 e

/0). (b) The function f () (as dened by Eq. 22) is shown to be roughly independent of the excess shear stress [Theta1][prime] = [Theta1] [Theta1]c. (c) Evolution of the critical sediment supply 0 with [Theta1] [Theta1]c. The dashed line is the best linear t.

If the sediment supply slightly increases, some pebbles will pop up above the rst layer and be easily entrained by the ow. Similar but less pronounced local minima can also be observed around = 2 and = 3.

Figure 7b conrms that the function f () =

Qs([Theta1] [Theta1]c)n()

, though not monotonous, is indeed inde

pendent of [Theta1]. The variation in Qs with

shows that the sediment ux rst 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 t the curve Qs() by a simple exponential function:

Qs ([Theta1],) = Qt ([Theta1]) 1 e/0

[parenrightbig]

Figure 7. a) Flux of sediment as a function of the sediment supply for different values of the Shields number.

The ux of sediment increases linearly for low sediment supplies, but tends to a saturated value when transport

capacity is reached. Plain lines are best ts by the equation Qs ( ,) = Qt ( )(1e/

0 ). b) The function

f () (as dened 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 t.

11

Transport capacity Q ( kg m s )

t

100

Slope 1.2

Slope 1.5

10

-1 )

-1 s

100

transport capacity Q t(kg m

1

slope 1.2

slope 1.5

10

.

(23)

Figure 8. Saturated

value Qt of the sediment transport rate as a function of ([Theta1] [Theta1]c)

/[Theta1]c. The plain line is the best t by a power law; the dashed line

is the t by a power law of index n = 1.5.

ure, our data could also be consistent with the power law Qt ([Theta1] [Theta1]c)1.5.

3.2 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 rst increases linearly with the sediment supply. Because of the roughness of the bedrock, it departs from the function C = beyond [similarequal] 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 [Theta1] 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

Earth Surf. Dynam., 4, 327342, 2016 www.earth-surf-dynam.net/4/327/2016/

Given this model, the surface density 3

0 corresponds to the maximum quantity of sediment that can

be transported by a ow of given Shields number. The variation in 0 with [Theta1] is plotted in Fig. 7c: it can be well tted by the afne function

0 = 20([Theta1] [Theta1]c). (24) This nding 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 [Theta1], whereas their velocity is proportional to the shear velocity, that is, to [Theta1]0.5. Finally in Fig. 8 we plot the maximum sediment transport rate Qt([Theta1]), reached for large values of the sediment supply 0. The

best t by a power law suggests Qt ([Theta1] [Theta1]c)1.2. How

ever, 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 g-

0.5 1 5 (c)/c

1

0.5 1 5 (QQc)/Qc

Figure 8. Saturated value of the sediment transport rate Qt as a function of (

c)/ c. The plain line is the

best t by a power law; the dashed line is the t by a power law of index n = 1.5.

28

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

0 0 1 2 3 4 5 6 Sediment supply

1

0.75

Cover fraction C

0.5

= 0.007 < c c = 0.015 c = 0.049

0.25

Figure 10. (a) Flux of energy received by the bedrock as a function of the relative excess shear stress, in a loglog scale. Plain lines are best ts of [Phi1]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 tted by a linear function.

The inset in Fig. 11 shows the variation in the maximum ux of energy with the Shields number. This variation can be well tted by an afne relationship. This demonstrates that the process of incision only happens if the Shields number exceeds an incision threshold [Theta1]i [similarequal] 0.025 > [Theta1]c. Therefore,

sediment transport can occur on a bedrock while not contributing to river incision if [Theta1]c < [Theta1] < [Theta1]i: in this regime, pebbles are rolling along the bedrock without impacting, and therefore do not contribute signicantly to its erosion.

3.4 Inuence of bedrock roughness

The roughness of the bedrock can be modied by varying two parameters: the surface density [notdef] of xed spheres on the bedrock and the protruding height of those pebbles over the horizontal bedrock (roughness height hb = R+zr). Figure 12

illustrates the aspect of the bedrock, viewed from above, for two different values of [notdef].

In order to study the variation in energy transfer with the bedrock roughness, we plot the ux of energy delivered to the bedrock with respect to the sediment supply () for a single value of [Theta1] [Theta1]c = 0.071, and for different roughness

congurations. Figure 13a shows the same plot for increasing

www.earth-surf-dynam.net/4/327/2016/ Earth Surf. Dynam., 4, 327342, 2016

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 = .

by the ow), where mobile pebbles have the freedom to roll along the bedrock. If the sediment supply increases, an incomplete static layer develops over the bedrock. In all cases the bedrock becomes entirely covered by a static layer if the sediment supply exceeds [similarequal] 3. Local maxima in the static

cover fraction are, once again, due to the intrinsic discontinuity of our uid model: when the bottom volume fraction exceeds b = 0.5, the uid velocity suddenly decreases within

the bottom layer.

3.3 Incision process

The ux 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 at surface or one of the glued spheres). In Fig. 10, we plot the variation in this ux of energy [Phi1]E with [Theta1] for different values of the sediment supply. As can be expected, [Phi1]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 [Phi1]E with [Theta1] can be tted by a power law whose index varies with the sediment supply:

[Phi1]E([Theta1],) = [Phi1]1 ()([Theta1] [Theta1]c)m(). (25) As shown in Fig. 10b, the exponent m() increases roughly linearly with , with m(0) [similarequal] 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: [Phi1]E rst 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 ux of energy delivered to the bedrock is obtained though the procedure described in Sect. 4.4.

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

0 0 1 2 3 4 Dimensionless sediment supply ()

70

60

max(W m)

2

60

50

i = c + 0.013

40

10

Incision rate (mm day)

1

50

30

20

2)

10

40

E ( W m

30

0 0 0.02 0.04 0.06 0.08

c

5

= 0.015 = 0.031 = 0.049 = 0.071

20

10

0

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 maxi

mum ux 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 ts by Eq. (27). The equivalent incision rate (right-side axis) is computed using Eq. (34) (see Sect. 4.4). The inset plots the maximum value of the energy ux as a function of [Theta1][Theta1]c. The

dotted line is the best afne ne, which reveals the existence of an incision threshold [Theta1]i > [Theta1]c.

values of [notdef], and a roughness height hb = 4cm. 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 efcient if the bed possesses a dense roughness. In Fig. 13b, we plot the same evolution for a given roughness density [notdef] = 0.36 but

for spheres protruding more or less within the ow. At high sediment supply, the energy received by the bedrock appears to be higher if the bedrock is smooth.

4 Discussion

4.1 Sediment transport rate

Although the model that we adopt for the interaction between the water ow 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-PeterMller law. When the transport rate is tted 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 ux is higher for low supplies: indeed, it is easier for pebbles to roll along a at bedrock than on the rough surface of a sediment layer. When [Theta1] is increased, the bedload ux increases faster for large sediment supply. For low sup-

plies, once all pebbles have been put into motion, only their velocity can increase with [Theta1]. If more pebbles are available, an increase in [Theta1] 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 ([Theta1]c [similarequal] 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 uid grain 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 ([Theta1]c = 0.012), the sed

iment ux that we obtain is, however, quantitatively consistent with common measurements (Meyer-Peter and Mller, 1948).

4.2 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 [Theta1], energy dissipated in the bedrock rst increases with the sediment supply and reaches a maximum for m [similarequal] 0.5 before decay

ing, 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 [similarequal] 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 ms of gravel of diameter d, within a cylindrical container of diameter D, we derive the corresponding dimensionless sediment supply through the expression

= N [notdef]

d2/4

D2/4 =

6mss d D2 . (26)

The comparison shows that our simulations are able to predict the right tendency for the ux of energy delivered to the

Earth Surf. Dynam., 4, 327342, 2016 www.earth-surf-dynam.net/4/327/2016/

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

Figure 12. Two different cases of bedrock roughness: the positions of the glued spheres are plotted as seen from above. Left: [notdef] = 0.08.

Right: [notdef] = 0.52. The ow is in the x direction, from left to right.

Normalized erosion rate (g hMPa)

1 2

= 0.049 (this study)

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)

30

40

30

2 )

20

E(W m

20

10

10

0

Figure 14. We plot on the same graph our numerical prediction for the ux of energy delivered to the bedrock (full circles), erosion rates measured experimentally by Sklar and Dietrich (2001) (empty symbols), and the best t 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 [Phi1]E has been chosen so that its maximum coincides with the experimental maximum erosion rate.

We also quantied the inuence of the Shields number on the abrasion process and showed a power-law dependency of [Phi1]E on the excess shear stress [Theta1] [Theta1]c. In particular, our ts

(see Fig. 10b) show that the exponent m is always greater than the index n of the transport law. This implies that if Qs and are kept constant, the incision rate increases with the Shields number. This result contradicts the prediction of the saltationabrasion model, where the incision rate scales like ([Theta1][Theta1]c)0.5 (Sklar and Dietrich, 2004), and ndings by

Chatanantavet and Parker (2009) and Johnson and Whipple (2010), where there is no explicit dependency of [Phi1]E in [Theta1] (for given values of Qs and ). However, it is consistent with the prediction of the shear stress and stream-power incision models (Whipple and Tucker, 1999). The variation in [Phi1]E with [Theta1] 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 amplied. However, let us remark that such events can also result in a sudden increase in the sedi-

www.earth-surf-dynam.net/4/327/2016/ Earth Surf. Dynam., 4, 327342, 2016

Figure 13. Flux of energy delivered to the bedrock as a function of the sediment supply for [Theta1][Theta1]c = 0.071: (a) for increasing values of

the roughness density [notdef] and (b) for different values of the roughness height.

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 uid 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.

0 0 1 2 3 Dimensionless sediment supply ()

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

0 1 2 3 4 5 Dimensionless sediment supply ()

(2014) modied the saltationabrasion 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 inuence of the coefcient of restitution on the results of our simulations should be of importance and will be the object of further investigation.Increasing the coefcient 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 coefcient 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.

4.3 Cover effect

By analogy with both the linear (Sklar and Dietrich, 2004) and the exponential (Turowski et al., 2007) cover models, we rst isolate in [Phi1]E the inuence of the sediment supply, and t the ux of energy by the empirical function

[Phi1]E ([Theta1],) = [Psi1] ([Theta1])pe/ . (27)

The values for p and , corresponding to the best ts plotted in Fig. 11, are reported in Table 2. The parameter appears to be roughly independent of the Shields number, and we consistently nd p > 2, whereas p = 1 in the exponen

tial cover model. The fact that p depends on [Theta1] underlines that we cannot express the incision rate as a simple product of a function of [Theta1] and a function of . Furthermore, in our simulations, forcing a t with p = 1 always leads to underes

timate the maximum incision rate and overestimate it at high sediment supply. This is likely due to the fact that, in the saltationabrasion model, the number of impacts is proportional to the sediment supply, whereas our simulations show that the sediment ux 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 ux of energy delivered to the bedrock, dened by Eq. (20), as the product of the energy provided by each impact (Ei), the number of impacts per unit time and surface (ni), and a cover function (F ):

[Phi1]E = Ei [notdef] ni [notdef] F. (28)

Earth Surf. Dynam., 4, 327342, 2016 www.earth-surf-dynam.net/4/327/2016/

100

1

0.8

F

0.6

0.4

10-1

0.2

0 0.25 0.5 0.75 1

F

1-C

10-2

= 0.015 = 0.031 = 0.049 = 0.071

10-3

Figure 15. Empty symbols: cover function F , as dened by Eq. (32), plotted as a function of , in a log-linear scale. Plain lines are best ts 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.

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 signicantly to the erosion of bedrock. This is inconsistent, however, with observations by Sklar and Dietrich (2001), who report abrasion occurring as soon as the ow is able to put sediment into motion. Let us note, however, that we vary the ow 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 [Theta1]c and [Theta1]i could explain that the discrepancy was not observed experimentally. If the existence of an incision threshold [Theta1]i > [Theta1]c were to be conrmed, it would mean that this value should be taken into account in models of river incision instead of the critical Shields number [Theta1]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 asperi-ties 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 [notdef] 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

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

Table 2. Coefcients used in the empirical t of the ux of energy as a function of sediment supply (see Eq. 27).

[Theta1] [Theta1]c p 0.015 2.11 0.210.031 2.77 0.180.049 2.70 0.190.071 2.00 0.24

In the saltationabrasion model, ni 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 ux Qs/Qt. In the stochastic model of Turowski et al.

(2007), F is expressed as an exponential function of Qs/Qt, 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 ux, 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 (,[Theta1]) that does not require a geometrical or stochastic description of the alluvial cover of the bedrock. The energy of an impact Ei 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 [Theta1] < [Theta1]c, we can therefore write, dimensionally,

Ei (sR3) [notdef] U 2 (sR3) [notdef] g R ([Theta1] [Theta1]c). (29)

As proposed by Foley (1980), the number of impacts per unit time and surface is expected to be proportional to the sediment transport rate Qs, which can be written dimensionally as

ni

, (33)

with 30 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 ume experiments) and Turowski and Rickenmann (2009) (in the eld). 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.

4.4 Estimation of the incision rate

We can estimate the rate of incision induced by the impacts on the bedrock, based on the ux of energy delivered. Fol-

www.earth-surf-dynam.net/4/327/2016/ Earth Surf. Dynam., 4, 327342, 2016

QssR4 . (30)

Therefore, the ux of energy should read

[Phi1]E = K [notdef] g [notdef] Qs(,[Theta1]) [notdef] ([Theta1] [Theta1]c) [notdef] F (,[Theta1]), (31) with K a dimensionless constant, which depends neither on [Theta1] nor on . This allows us to redene the cover function as

F (,[Theta1]) =

[Phi1]Eg Qs ([Theta1] [Theta1]c)

. (32)

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

[Theta1] [Theta1]c '0.015 2.86 [notdef] 0.03

0.031 4.65 [notdef] 0.06

0.049 6.32 [notdef] 0.18

0.071 8.43 [notdef] 0.26

In Fig. 15 we plot F in a log-linear scale. Let us rst observe that K is indeed independent of [Theta1], 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 [lessorsimilar] 0.75. On the other hand, F can be well tted by an exponential decay for [greaterorsimilar]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 Qs/Qt 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 Qs/Qt [greaterorsimilar] 0.25 0.75 (see their Fig. 13). Our result therefore

conrms that the exponential cover model overestimates the cover effect at low sediment supplies but ts correctly at high enough sediment supplies. In the exponential regime, following the stochastic approach by Turowski et al. (2007), we can t the cover function by

F = A exp

'

30

340 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

I =

2Y , (34)

where [epsilon1]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 [notdef] 106 Pa is the rock tensile strength, Y = 5.0 [notdef] 1010 Pa

is the rock elastic modulus, and kv = 106 is a dimension

less 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 [Phi1]E = 40W m2 (obtained at [Theta1] [Theta1]c = 0.049 and = 0.6)

is I = 2.6myr1 = 7mm day1. 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 veried 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 long-term average incision rate for a given stream, provided that the probability distribution functions of both water discharge and sediment supply are known.

5 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 ow, with a simplied retroaction of the presence of the pebbles on the ow. This allowed us to explicitly compute the trajectories of pebbles transported by the ow, 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 tted by a power law of the Shields number, similar to most classical transport laws at saturation. However, we also evidenced the inuence 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 bed-load ux increases with the sediment supply until it reaches its saturated value. This allowed us to compute the sediment mass that the ow is able to transport. However, extracting a

unique general expression for the ux 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 ux 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 rst 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 dened 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 ow by a horizontally averaged and purely horizontal velocity prole 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 efciency of the cover effect. On the other hand, taking into account turbulent velocity uctuations, 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 eld, and ideally the exact velocity eld 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 ow (for instance, a at pebble would be harder to dislodge from the bedrock; Rust, 1972; Komar and Li, 1986). In con-

Earth Surf. Dynam., 4, 327342, 2016 www.earth-surf-dynam.net/4/327/2016/

[Phi1]E

T 2

[epsilon1]v with [epsilon1]v = kv

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

trast, once entrained by the ow 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).

Acknowledgements. Data supporting the results of this article can be obtained by contacting Vincent Langlois ([email protected]).

The authors would like to acknowledge the insightful input of four reviewers, amongst whom P. Chatanantavet, on an earlier version of this article.

Edited by: E. Lajeunesse

References

Ancey, C. and Heyman, J.: A microstructural approach to bed load transport: mean behaviour and uctuations of particle transport rates, J. Fluid Mech., 744, 129168, 2014.

Ancey, C., Bhm, T., Jodeau, M., and Frey, P.: Statistical description of sediment transport experiments, Phys. Rev. E, 74, 011302, doi:http://dx.doi.org/10.1103/PhysRevE.74.011302

Web End =10.1103/PhysRevE.74.011302 http://dx.doi.org/10.1103/PhysRevE.74.011302

Web End = , 2006.

Anderson, R. S.: Evolution of the Santa Cruz Mountains, California, through tectonic growth and geomorphic decay, J. Geophys.Res.-Sol. Ea., 99, 2016120179, 1994.

Attal, M. and Lav, J.: Pebble abrasion during uvial transport: Experimental results and implications for the evolution of the sediment load along rivers, J. Geophys. Res., 114, F04023, doi:http://dx.doi.org/10.1029/2009JF001328

Web End =10.1029/2009JF001328 http://dx.doi.org/10.1029/2009JF001328

Web End = , 2009.

Bitter, J.: A study of erosion phenomena part I, Wear, 6, 521, 1963. Bhm, T., Ancey, C., Frey, P., Reboud, J.-L., and Ducottet,C.: Fluctuations of the solid discharge of gravity-driven particle ows in a turbulent stream, Phys. Rev. E, 69, 061307, doi:http://dx.doi.org/10.1103/PhysRevE.69.061307

Web End =10.1103/PhysRevE.69.061307 http://dx.doi.org/10.1103/PhysRevE.69.061307

Web End = , 2004.

Bufngton, J. M. and Montgomery, D. R.: A systematic analysis of eight decades of incipient motion studies, with special reference to gravel-bedded rivers, Water Resour. Res., 33, 19932029, doi:http://dx.doi.org/10.1029/96WR03190

Web End =10.1029/96WR03190 http://dx.doi.org/10.1029/96WR03190

Web End = , 1997.

Charru, F., Mouilleron, H., and Eiff, O.: Erosion and deposition of particles on a bed sheared by a viscous ow, J. Fluid Mech., 519, 5580, 2004.

Chatanantavet, P. and Parker, G.: Experimental study of bedrock channel alluviation under varied sediment supply and hydraulic conditions, Water Resour. Res., 44, doi:http://dx.doi.org/10.1029/2007WR006581

Web End =10.1029/2007WR006581 http://dx.doi.org/10.1029/2007WR006581

Web End = , 2008.

Chatanantavet, P. and Parker, G.: Physically based modeling of bedrock incision by abrasion, plucking, and macroabrasion, J.Geophys. Res., 114, F04018, doi:http://dx.doi.org/10.1029/2008JF001044

Web End =10.1029/2008JF001044 http://dx.doi.org/10.1029/2008JF001044

Web End = , 2009.

Clift, R., Grace, J., and Weber, M.: Bubbles, Drops, and Particles,

Academic Press, New York, 1978.

Cook, K. L., Turowski, J. M., and Hovius, N.: River gorge eradication by downstream sweep erosion, Nat. Geosci., 7, 682686, 2014.

Cundall, P. A. and Strack, O. D. L.: A discrete numerical model for granular assemblies, Gotechnique, 29, 4765, 1979.

Davis, R. H., Serayssol, J.-M., and Hinch, E.: The elastohydrodynamic collision of two spheres, J. Fluid Mech., 163, 479497, 1986.

Duran, O., Andreotti, B., and Claudin, P.: Numerical simulation of turbulent sediment transport, from bed load to saltation, Phys.Fluids, 24, 103306, doi:http://dx.doi.org/10.1063/1.4757662

Web End =10.1063/1.4757662 http://dx.doi.org/10.1063/1.4757662

Web End = , 2012.

Engle, P. A.: Impact Wear of Materials, Elsevier Sci., New York,

1978.

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.

Foley, M. G.: Bed-rock incision by streams, Geol. Soc. Am. Bull.,

91, 21892213, 1980.

Gilbert, G.: Report on the Geology of the Henry Mountains, US Geogr. Geol. Surv., Rocky Mountains Region, Washington, 1877.Hartshorn, K., Hovius, N., Dade, W. B., and Slingerland, R. L.:

Climate-driven bedrock incision in an active mountain belt, Science, 297, 20362038, 2002.

Hodge, R. A. and Hoey, T. B.: Upscaling from grain-scale processes to alluviation in bedrock channels using a cellular automaton model, J. Geophys. Res.-Earth, 117, doi:http://dx.doi.org/10.1029/2011JF002145

Web End =10.1029/2011JF002145 http://dx.doi.org/10.1029/2011JF002145

Web End = , 2012.

Howard, A. D.: A detachment-limited model of drainage basin evolution, Water Resour. Res., 30, 22612285, 1994.

Huda, S. A. and Small, E. E.: Modeling the Effects of Bed Topography on Fluvial Bedrock Erosion by Saltating Bed Load, J. Geophys. Res.-Earth, 119, 12221239, doi:http://dx.doi.org/10.1002/2013JF002872

Web End =10.1002/2013JF002872 http://dx.doi.org/10.1002/2013JF002872

Web End = , 2014.

Johnson, J. P. and Whipple, K. X.: Evaluating the controls of shear stress, sediment supply, alluvial cover, and channel morphology on experimental bedrock incision rate, J. Geophys. Res.-Earth, 115, doi:http://dx.doi.org/10.1029/2009JF001335

Web End =10.1029/2009JF001335 http://dx.doi.org/10.1029/2009JF001335

Web End = , 2010.

Komar, P. D. and Li, Z.: Pivoting analyses of the selective entrainment of sediments by shape and size with application to gravel threshold, Sedimentology, 33, 425436, 1986.

Lague, D.: Reduction of long-term bedrock incision efciency by short-term alluvial cover intermittency, J. Geophys. Res., 115, F02011, doi:http://dx.doi.org/10.1029/2008JF001210

Web End =10.1029/2008JF001210 http://dx.doi.org/10.1029/2008JF001210

Web End = , 2010.

Lague, D.: The stream power river incision model: evidence, theory and beyond, Earth Surf. Proc. Land., 39, 3861, 2014.Lajeunesse, E., Malverti, L., and Charru, F.: Bed load transport in turbulent ow at the grain scale: Experiments and modeling, J.Geophys. Res., 115, F04001, doi:http://dx.doi.org/10.1029/2009JF001628

Web End =10.1029/2009JF001628 http://dx.doi.org/10.1029/2009JF001628

Web End = , 2010.Lamb, M. P. and Fonstad, M. A.: Rapid formation of a modern bedrock canyon by a single ood event, Nat. Geosci., 3, 477 481, 2010.

Lamb, M. P., Dietrich, W. E., and Venditti, J. G.: Is the critical Shields stress for incipient sediment motion dependent on channel-bed slope?, J. Geophys. Res.-Earth, 113, doi:http://dx.doi.org/10.1029/2007JF000831

Web End =10.1029/2007JF000831 http://dx.doi.org/10.1029/2007JF000831

Web End = , 2008.

www.earth-surf-dynam.net/4/327/2016/ Earth Surf. Dynam., 4, 327342, 2016

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

Liu, Q. and Prosperetti, A.: Wall effects on a rotating sphere, J. FluidMech., 657, 121, doi:http://dx.doi.org/10.1017/S002211201000128X

Web End =10.1017/S002211201000128X http://dx.doi.org/10.1017/S002211201000128X

Web End = , 2010.

Mtivier, F. and Meunier, P.: Input and output mass ux correlations in an experimental braided stream. Implications on the dynamics of bed load transport, J. Hydrol., 271, 2238, 2003.Meyer-Peter, E. and Mller, R.: Formula for the bedload transport, in: Proceedings of the 3rd Meeting of the International Association of Hydraulic Research, Stockholm, 3964, 1948.Nikuradse, J.: Gesetzmigkeiten der turbulenten Strmung in glatten Rohren (Nachtrag), Forsch. Ingenieurwes., 4, 4444, 1933.Pschel, T. and Schwager, T.: Computational Granular Dynamics,

Springer, Berlin, 2005.

Rickenmann, D. and Recking, A.: Evaluation of ow resistance in gravel-bed rivers through a large eld data set, Water Resour.Res., 47, doi:http://dx.doi.org/10.1029/2010WR009793

Web End =10.1029/2010WR009793 http://dx.doi.org/10.1029/2010WR009793

Web End = , 2011.

Roseberry, J. C., Schmeeckle, M. W., and Furbish, D. J.: A probabilistic description of the bed load sediment ux: 2.Particle activity and motions, J. Geophys. Res.-Earth, 117, doi:http://dx.doi.org/10.1029/2012JF002353

Web End =10.1029/2012JF002353 http://dx.doi.org/10.1029/2012JF002353

Web End = , 2012.

Rust, B. R.: Pebble orientation in uvial sediments, J. Sediment.Res., 42, 384388, 1972.

Schmeeckle, M. W., Nelson, J. M., Pitlick, J., and Bennett, J. P.: Interparticle collision of natural sediment grains in water, Water Resour. Res., 37, 23772391, 2001.

Seidl, M. A., Dietrich, W. E., and Kirchner, J. W.: Longitudinal prole development into bedrock: an analysis of Hawaiian channels,J. Geol., 457474, 1994.

Shepherd, R.: Incised river meanders: Evolution in simulated bedrock, Science, 178, 409411, 1972.

Sklar, L. S. and Dietrich, W. E.: Sediment and rock strength controls on river incision into bedrock, Geology, 29, 10871090, 2001. Sklar, L. S. and Dietrich, W. E.: A mechanistic model for river incision into bedrock by saltating bed load, Water Resour. Res., 40, doi:http://dx.doi.org/10.1029/2003WR002496

Web End =10.1029/2003WR002496 http://dx.doi.org/10.1029/2003WR002496

Web End = , 2004.

Turowski, J. M. and Bloem, J.-P.: The inuence of sediment thickness on energy delivery to the bed by bedload impacts, Geodin. Acta, 28, 110, doi:http://dx.doi.org/10.1080/09853111.2015.1047195

Web End =10.1080/09853111.2015.1047195 http://dx.doi.org/10.1080/09853111.2015.1047195

Web End = , 2015.

Turowski, J. M. and Rickenmann, D.: Tools and cover effects in bedload transport observations in the Pitzbach, Austria, Earth Surf. Proc. Land., 34, 2637, 2009.

Turowski, J. M., Lague, D., and Hovius, N.: Cover effect in bedrock abrasion: A new derivation and its implications for the modeling of bedrock channel morphology, J. Geophys. Res., 112, F04006, doi:http://dx.doi.org/10.1029/2006JF000697

Web End =10.1029/2006JF000697 http://dx.doi.org/10.1029/2006JF000697

Web End = , 2007.

Valance, A.: Formation of ripples over a sand bed submitted to a turbulent shear ow, Eur. Phys. J. B, 45, 433442, doi:http://dx.doi.org/10.1140/epjb/e2005-00201-9

Web End =10.1140/epjb/e2005-00201-9 http://dx.doi.org/10.1140/epjb/e2005-00201-9

Web End = , 2005.

Whipple, K. X. and Tucker, G. E.: Dynamics of the stream-power river incision model: Implications for height limits of mountain ranges, landscape response timescales, and research needs, J. Geophys. Res., 104, 17661, doi:http://dx.doi.org/10.1029/1999JB900120

Web End =10.1029/1999JB900120 http://dx.doi.org/10.1029/1999JB900120

Web End = , 1999.

Whipple, K. X., Hancock, G. S., and Anderson, R. S.: River incision into bedrock: Mechanics and relative efcacy of plucking, abrasion, and cavitation, Geol. Soc. Am. Bull., 112, 490503, 2000.

Earth Surf. Dynam., 4, 327342, 2016 www.earth-surf-dynam.net/4/327/2016/