A mechanistic model for river incision into bedrock by saltating bed load

1. Introduction

[2] Surprisingly little is known about the physical controls on rates of erosion of bedrock river beds [Whipple et al., 2000a], given the importance of river incision into bedrock in driving landscape evolution [e.g., Howard et al., 1994] and linking climate and tectonics to topography [e.g., Molnar and England, 1990]. Most numerical models of landscape evolution use bedrock incision rules that are related to the physics of bedrock wear only through simple scaling arguments, typically by assuming that incision rate is proportional to some measure of flow intensity such as unit stream power [Seidl and Dietrich, 1992; Anderson, 1994; Tucker and Slingerland, 1994; Willett, 1999] or average boundary shear stress [Howard and Kerby, 1983; Howard, 1994; Moglen and Bras, 1995; Tucker and Slingerland, 1996]. These models are difficult to apply in real landscapes because they lump the potentially complex influence of fundamental variables such as rock strength, channel slope, discharge, sediment supply and grain size, into a set of poorly constrained parameters unrelated to any particular erosional mechanism. Model parameters are calibrated using longitudinal profiles of rivers with a known incision history [e.g., Seidl et al., 1994; Stock and Montgomery, 1999; Whipple et al., 2000b; van der Beek and Bishop, 2003] or by assuming steady state topography [e.g., Sklar and Dietrich, 1998; Snyder et al., 2000], in effect using landscape form to infer process. Here we take the reverse approach, and develop a physically based model for bedrock incision by a single mechanism, abrasion by saltating bed load. Our goals are (1) to explore whether detailed consideration of the mechanics of bedrock wear will lead to a view of river incision fundamentally different from the simple scaling models, (2) to use the model to refine mechanistic hypotheses that can be tested in the field and in laboratory experiments, and (3) to explain river profile form in terms of the dynamics of bedrock incision.

[3] Many erosional mechanisms operate at the interface between a river and its rock bed, including cavitation, dissolution, abrasion by both bed load and suspended load, and the plucking of rock fragments by fluid shear stresses acting without the aid of sediment impacts. We focus on bed load abrasion for several reasons: (1) all rivers carry sediment in some quantity so abrasion should be a ubiquitous process; (2) other mechanisms occur only under limited conditions, for example cavitation requires extreme flow velocities and plucking may primarily occur in weak or highly fractured rock; (3) particle impacts are an efficient pathway for the transfer of energy and momentum from the flow to the bed; (4) by definition, bed load interacts with the underlying bedrock far more frequently than does suspended load; (5) much is known about the motion of saltating bed load from both theoretical and experimental studies; and (6) we have direct experimental evidence that saltating bed load causes measurable bedrock wear [Sklar and Dietrich, 2001]. Furthermore, by focusing on abrasion by saltating bed load, we are forced to consider in detail the role of sediment in influencing incision rates, not only in providing tools for abrasion, but also in controlling the extent of bedrock exposure on the channel bed. Ultimately, by focusing on bed load abrasion, we hope to gain insight into the potential feedbacks between incision rate, sediment supply and grain size that may strongly influence landscape response to changes in climate or rock uplift rate.

[4] Actively incising channels in hilly and mountainous terrain exhibit a wide range in the extent of alluvial cover and size distribution of bed sediment. Figure 1 shows photographs of two bedrock canyon channel beds in the Coast Range mountains of southern Oregon and northern California. The bedrock bed of the Smith River is nearly devoid of sediment deposits, in contrast to the South Fork Eel River, where bedrock is exposed only in isolated patches amid a nearly continuous mantle of coarse alluvium. Yet both of these channels, which are separated by about 150 km, are actively incising through similar marine sedimentary rocks. To what extent does the size and rate of sediment supply to these two channels contribute to or inhibit the lowering of the bedrock bed?

[5] Recent experimental results from bedrock abrasion mills [Sklar and Dietrich, 2001] support the hypothesis of Gilbert [1877] that the quantity of sediment supplied to the river should influence bedrock incision rates in two essential yet opposing ways: by providing tools for abrasion of exposed bedrock and by limiting the extent of exposure of bedrock in the channel bed. Incision is limited at lower supply rates by a shortage of abrasive tools (what we term the “tools effect”) and at higher supply rates by partial burial of the bedrock substrate beneath transient sediment deposits (“the cover effect”). Thus, as observed experimentally, the maximum rate of river incision should occur when the stream receives a moderate supply of sediment, relative to its sediment transport capacity. It follows that the size distribution of sediment grains supplied to the channel should also influence incision rates because only the coarser fraction is capable of forming an alluvial cover and because the finer fraction is carried in suspension and rarely collides with the bedrock bed. As suggested by the bedrock abrasion mill results [Sklar and Dietrich, 2001], the most efficient abrasive tools are sediments of intermediate size, large enough to travel as bed load rather than suspended load but not so large as to be immobile.

[6] Gilbert's hypothesis is also consistent with numerous field observations. Sediment clearly contributes to incision through the scouring of potholes [e.g., Alexander, 1932] and longitudinal grooves [Wohl, 1992a]. Boulders deposited in extreme flows [Wohl, 1992b] or eroded from steep canyon walls [Seidl et al., 1994] have been shown to shut off the incision process over decadal to millennial timescales. In contrast, incision under a thin alluvial cover has been documented in a badlands eroded into weak shale [Howard and Kerby, 1983]. Bedrock channel beds partially mantled with alluvial deposits are common in actively incising terrain [Howard, 1998] and occur even in lowland coastal streams [Ashley et al., 1988]. Montgomery et al. [1996] and Massong and Montgomery [2000] have shown that, for a given drainage area, a wide range in extent of alluviation can occur over a narrow range in channel slope, due in part to the influence of large woody debris. This suggests that the extent of bedrock exposure in incising channels may be highly dynamic.

[7] Sediment supply and grain size must be explicitly represented in bedrock incision models to adequately capture the potentially complex influence of sediment on incision rates. Although to first approximation the long-term flux of sediment to a given channel reach will be simply the product of the landscape denudation rate times the upstream drainage area, the supply rate of sediment in the bed load size class will be considerably less. Some mass eroded from upstream is transported as dissolved load. How much of the solid load moves as bed load, and thus contributes to bed load abrasion and partial bed cover, will depend on the size distribution of sediment entering the channel network, the rate of particle comminution with downstream travel distance, and the fluid shear stress available to transport sediment. Unlike fully alluvial channels in depositional environments, the sediment load in mixed bedrock-alluvial channels is commonly assumed to be less than the sediment transport capacity, because there is a limited supply of sediment available to be scoured from the channel bed. Therefore the transport capacity of adjacent upstream reaches cannot be used as a proxy for local sediment supply. At the short timescale of individual incision events, sediment delivery can be highly episodic, particularly in steep landslide-dominated terrain [e.g., Benda and Dunne, 1997]. Significant spatial and temporal variability in sediment supply should also occur due to heterogeneous lithology, nonuniform rates of landscape lowering, fans and other upstream sediment storage elements, and the discrete locations of sediment input such as tributary junctions and stream bank failures.

[8] Few bedrock incision models have attempted to explicitly account for the influence of sediment supply and grain size on incision rate, although many landscape evolution models [e.g., Howard, 1994; Tucker and Slingerland, 1994] represent the cover effect as a sudden shift from bedrock to alluvial conditions when sediment supply exceeds transport capacity. Beaumont et al. [1992] proposed that bedrock incision could be treated like alluvial incision, with incision rate scaling with the sediment transport capacity in excess of the sediment supply, in effect modeling the coverage effect but not the tools effect. Foley [1980] modified Bitter's [1963a, 1963b] sandblast abrasion model to apply to wear by saltating bed load, accounting for the tools effect but neglecting the coverage effect. A preliminary version of the full model developed here [Sklar et al., 1996; Sklar and Dietrich, 1998] captured both the coverage and tools effects and has inspired efforts to incorporate sediment into the simple stream power incision rule [Whipple and Tucker, 2002]. Slingerland et al. [1997] used studies of slurry pipeline wear [e.g., de Bree et al., 1982] to support the view that the coverage and tools effects should occur in rivers. Here we describe in detail a model for bedrock incision by saltating bed load based directly on the mechanics of bed load sediment transport and the wear of brittle materials by particle impacts.

[9] The paper is divided into four sections. We begin by deriving an expression for the rate of bedrock incision by saltating bed load, using established theory for wear of brittle materials by low-velocity particle impacts, in which the erosion rate depends on saltation hop length and impact velocity. We then use data from a suite of published studies on saltation trajectories to parameterize the model directly in terms of shear stress. Next, we explore the sensitivity of the predicted bedrock abrasion rates to variations in sediment supply, shear stress and grain size, and then use a nondimensional framework to collapse the full model behavior onto a single graph. Finally, we compare the model predictions to laboratory measurements of bedrock wear rates [Sklar and Dietrich, 2001] to calibrate the one adjustable parameter, which characterizes rock resistance to fluvial abrasion.

2. Model Development

2.2. General Expression

[13] Consider an idealized flux of saltating bed load, moving across an approximately planar river bed composed of a mixture of exposed bedrock and alluvial patches, as depicted in Figure 2. Bedrock erosion occurs as the integrated result of numerous grain collisions with exposed bedrock. Many nondestructive grain impacts occur as sediment travels across local sediment deposits where bedrock is covered. Bedrock incision rate (E) can be expressed simply as the product of three terms: the average volume of rock detached per particle impact (V_i), the rate of particle impacts per unit area per unit time (I_r), and the fraction of the river bed made up of exposed bedrock (F_e)

Below we develop expressions for each of these terms using empirical relationships for wear by particle impacts, bed load sediment transport capacity, and the trajectory of saltating bed load motion.

[14] Saltating grains move in a characteristic trajectory (Figure 3). The magnitude of particle velocities and saltation hop length and height are well described by functions of boundary shear stress in excess of the threshold of motion [e.g., Wiberg and Smith, 1985], as is bed load sediment transport capacity [e.g., Gomez and Church, 1989]. Ultimately, we will obtain an expression for erosion rate in terms of sediment supply rate, grain size, excess shear stress, and rock strength.

2.3. Volume Eroded Per Unit Impact (V_i)

[15] Erosion of brittle materials by low-velocity particle impacts occurs through the formation, growth and intersection of a network of cracks [Engle, 1978]. Cracks occur due to tensile stresses arising from elastic deformation around the point of impact. Detachment of an eroded fragment is generally the result of an accumulation of intersecting cracks due to multiple particle impacts. Erosion of brittle materials by repeated impacts is referred to as “deformation wear” [e.g., Bitter, 1963a]. Although the volume eroded by any single impact will depend on the local fracture density, the wear rate averaged over a large number of impacts scales with the flux of kinetic energy transferred by the impacting grains [Engle, 1978].

[16] Here we treat bedrock as an elastic brittle material, and thus do not consider ductile or nonelastic behavior. The net transfer of kinetic energy (Δɛ_N) in a collision between a small particle and an essentially infinite rock mass is

where M_p is the mass of the impacting particle, U_i and U_r are the velocities of impact and rebound respectively, and (1 − U_i²/U_r²) = f(sin α) where α is the angle of impact [Head and Harr, 1970]. Some of the kinetic energy transferred may contribute to deformation wear of the bedrock (Δɛ_D), while the remainder is lost (Δɛ_L) through transformation into other energy forms such as sound and heat (i.e., Δɛ_N = Δɛ_D + Δɛ_L). Thus, bedrock wear depends on the energy transferred in excess of losses, Δɛ_D = Δɛ_N − Δɛ_L. This forms the basis of the classic impact wear model of Bitter [1963a], which we write as

where V_i is the mean volume eroded per impact, ɛ_v is the total energy required to erode a unit volume of rock, and ɛ_t is the threshold energy that must be exceeded for detachment to occur. Erosion depends on the sine of the impact angle because the magnitude of the peak tensile stress varies with the normal component of the impact velocity [Engle, 1978].

[17] Note that the tangential velocity component can contribute to erosion by scratching or gouging, referred to as “cutting wear” [Bitter, 1963b]. Cutting wear is important in ductile materials and when the impacting particles are highly angular, but is not significant when brittle materials are impacted by rounded grains [Head and Harr, 1970], as is generally the case in bedrock rivers. Note also that this formulation implicitly assumes that the mass of impacting particles is conserved, because it does not account for the energy that may contribute to wear of the impacting grain or for the kinetic energy of any fragments that may be produced by the grain impact.

[18] The linear dependence of impact abrasion rate on excess energy is illustrated in Figure 4 by the experimental data of Head and Harr [1970], who eroded brittle plate glass with a stream of air-blown 0.04-mm diameter glass beads. In these experiments, they varied the energy flux by changing the air velocity and determined the net energy transferred by measuring particle impact and rebound velocities with high-speed photography. Similarly, the influence of impact angle on abrasion rates is shown in Figure 5. Here data from Head and Harr [1970] are plotted along with results from experiments by Bitter [1963a] in which round 0.3-mm diameter iron pellets were dropped on plate glass from a height of five meters. Note that Head and Harr's [1970] results are consistent with a negligible erosion threshold, unlike Bitter's [1963a] data, which do indicate an erosion threshold possibly due to the use of ductile rather than brittle abrasive material.

[19] The resistance of rock and other brittle materials to abrasion by impacting particles (parameterized in equation (3) by ɛ_v) depends on the capacity of the material to store energy elastically [Engle, 1978]. Referred to variously as the material toughness, resilience or strain energy, the capacity to store energy elastically (β) is defined quantitatively as the area under the stress-strain curve at the yield stress, so that

where σ_T is the rock tensile yield strength, Y is Young's modulus of elasticity and k_v is a dimensionless coefficient that will depend in part on the material properties of the impacting particle. Because the variation in modulus of elasticity of rocks is limited [Clark, 1966], to first order Y can be treated as constant. Thus bedrock abrasion resistance scales directly with the square of rock tensile strength, as demonstrated in the bedrock abrasion mill experiments of Sklar and Dietrich [2001]. In principle, the erosion threshold parameter ɛ_t should also scale with rock tensile strength, however, the abrasion mill data [Sklar and Dietrich, 2001] suggest that the energy threshold ɛ_t can be neglected in the case of bed load abrasion because even the relatively low-impact energy of fine sand moving as bed load is sufficient to cause measurable bedrock wear. We henceforth set ɛ_t = 0, an assumption supported by comparison of model predictions with the experimental data (section 5). Thus, for erosion by bed load sediment moving by saltation, equation (3) can be rewritten as

where πρ_sD_s³/6 = M_p, ρ_s is the density of sediment, D_s is the diameter of a spherical sediment grain, and w_si is the vertical component of the particle velocity on impact (w_si = U_i sin α) (Figure 3).

2.4. Particle Impact Rate (I_r)

[20] The rate of particle impacts, per unit time and per unit area (I_r), should be linearly proportional to the flux of bed load particles and inversely proportional to the downstream distance between impacts. The particle flux is simply the mass flux per unit width (q_s) divided by the mass per particle, and the distance traveled between impacts is the saltation hop length (L_s), so that particle impact rate can be written as

2.5. Fraction of Bed Exposed (F_e)

[21] When sediment mantles portions of the bedrock bed, some of the bed load kinetic energy acquired from the surrounding flow is expended in collisions with immobile sediment rather than exposed bedrock. Because transient sediment deposits are common in actively incising rivers, we assume a fraction of saltation hops occur as bed load across local patches of alluvium. The kinetic impact energy is thus lost to intergranular friction, without contributing to bedrock fracture growth, except perhaps when the alluvial cover is only one or two grain diameters thick. Additional energy is lost in collisions between mobile grains.

[22] We parameterize the effect of local alluvial cover and grain-to-grain interaction as a function of the sediment transport capacity in excess of the sediment supply. Consider the end-member cases. When the channel has just enough transport capacity to move the coarse sediment supplied from upstream and from local hillslopes, the bed should tend to be fully alluviated. Any increase in the rate of bed load sediment supply (assuming no change in grain size) will cause aggradation until a steeper slope is achieved so that transport capacity again equals the sediment supply rate. In contrast, when the channel receives no coarse sediment supply, as might occur in mechanically weak lithologies that weather rapidly to suspended load size particles, no coarse sediment is available to form patches of alluvial cover, and bedrock should be exposed across the entire bed. Here we assume that the fraction of the bed composed of exposed bedrock (F_e) varies linearly between these end-members, so that

where q_t is the sediment transport capacity for a fully alluviated bed. In subsequent calculations we use the Fernandez-Luque and van Beek [1976] bed load sediment transport relation

where q_t is the sediment mass transport capacity per unit width, ρ_s is the sediment density, R_b is the nondimensional buoyant density of sediment (R_b = ρ_s/ρ_w − 1), ρ_w is the density of water, g is the acceleration due to gravity, τ* = τ_b/[(ρ_s − ρ_w)gD_s] is the nondimensional form of the boundary shear stress (τ_b), and τ*_c is the value of τ* at the threshold of particle motion.

[23] The linear assumption of equation (7) is supported by flume experiments [Sklar and Dietrich, 2002] in which we observe a generally linear increase in bed cover with increasing sediment supply. The experiments also demonstrate the potential for more complex behavior, such as a threshold supply rate below which alluviation does not occur and instabilities associated with changing bed roughness, particularly for the case of smooth, planar bedrock channels. We do not attempt to model such nonlinear behavior, in effect assuming that it will be suppressed by the irregular topography typical of bedrock channel beds.

2.6. Composite Expression

[24] Our simple expression for bedrock abrasion rate, E = V_iI_rF_e, can now be restated by substituting equations (5), (6), and (7) into equation (1) and simplifying, resulting in the composite expression for bedrock incision rate

3. Saltation Trajectories

[25] In this section we develop expressions for the saltation hop length (L_s) and the vertical component of the particle impact velocity (w_si), so that we can express bedrock abrasion rate (equation (9)) directly in terms of shear stress. The most mechanistically rigorous approach would be to build a model for saltation over both bedrock and alluvial beds that solves the equations of motion for an individual grain throughout the saltation trajectory. Theoretical models have been developed for saltation in water over an alluvial bed [Wiberg and Smith, 1985; Sekine and Kikkawa, 1992] and for saltation in air [e.g., Anderson and Haff, 1988]. Such models need to account for the forces due to hydrodynamic lift and drag, gravity, the virtual mass effect (caused by the relative accelerations of the fluid and particle), the spinning motion of saltating grains, and collisions with the bed [Wiberg and Smith, 1985]. This level of complexity is unnecessary for our purposes here. Instead, we use a large set of published data to obtain a set of empirical expressions that together constitute a self-consistent description of how saltation trajectories vary as a function of shear stress and grain size.

3.2. Mean Downstream Saltation Velocity (u_s)

[29] We find the best data collapse when we nondimensionalize u_s by dividing by (R_bgD_s)^1/2. Figure 6 shows the variation in mean downstream saltation velocity with transport stage (τ*/τ*_c). The best fit log-log linear regression line through the data ensemble is

Goodness of fit for this and subsequent regressions is indicated by the standard error of the fit parameters and is shown in Figures 6, 7, and 9. Note that the regression fit predicts a nonzero mean downstream velocity at the threshold of motion. A more physically consistent approach is to treat u_s as a function of nondimensional excess shear stress (τ*/τ*_c − 1) so that velocity goes to zero at the threshold of motion. Regressing the data in this manner (Figure 7a), we obtain

Equation (11) is plotted in Figure 6 as the thick line.

[30] Equation (11) better represents the trend in the data because there is no strong curvature in the residuals, unlike equation (10); note also that the uncertainty in the parameter estimates is lower. Equation (11) is also more consistent with the fact that rolling and sliding occur in the transition between no motion and full saltation. The vast majority of grains move in saltation mode once τ*/τ*_c reaches about 1.2 [Abbott and Francis, 1977; Hu and Hui, 1996a]. Therefore we henceforth treat downstream velocity and the other components of the saltation trajectory as functions of τ*/τ*_c − 1 (Figure 7) to account for the reduced intensity of particle motion at low excess shear stresses, and to be consistent with bed load sediment transport capacity relations (e.g., equation (8)), which predict zero sediment transport at the threshold-of-motion shear stress.

3.3. Saltation Hop Height (H_s)

[31] Figure 7b shows the variation in hop height with excess shear stress, nondimensionalized by particle diameter. Maximum observed hop heights are about five particle diameters, for large excess shear stresses. The log-log linear regression line that best fits the data ensemble is

3.4. Saltation Hop Length (L_s) and the Transition to Suspended Transport

[32] Figure 7c shows the variation in saltation hop length with excess shear stress, nondimensionalized by particle diameter. The log-log linear best fit relationship is

This relationship, however, should not hold for very large excess shear stresses because of the transition from bed load to suspended particle motion. Equations (11) and (12) similarly should not apply when the transport mode shifts from bed load to suspended load. The detailed observations of Abbott and Francis [1977] show that rising turbulent eddies interrupt the downward acceleration of descending saltating grains even at low excess shear stresses. As shear stress increases, such “suspensive” saltation hops become more frequent. Full suspended transport, when the upward diffusion of momentum balances gravity, occurs when the shear velocity (u* = (τ_b/ρ_w)^1/2) exceeds the particle fall velocity in still water (w_f) [Rouse, 1937]. The transition from bed load to suspension is equivalent to a grain taking a hop of infinite length.

[33] To account for the shift in the mode of particle transport, we divide equation (13) by a nonlinear function of u*/w_f so that hop length grows rapidly as u* approaches w_f and becomes infinite at the threshold of suspension

This nonlinear extrapolation is shown as the thick line in Figure 7c. Note that significant deviation from the log-log linear regression fit occurs only above excess shear stresses of about 20. The fact that no experimental data exist for any of the three saltation trajectory components at excess shear stresses above this range is consistent with a shift from saltation to suspended motion.

[34] The transport stage that corresponds to the threshold of suspension (i.e., the value of τ*/τ*_c when u*/w_f = (τ_b/ρ_w)^1/2/w_f = 1) and thus the excess shear stress at which saltation hop length becomes infinite, depends nonlinearly on grain size

because of the dependence of fall velocity on grain size. Dietrich [1982] provided an empirical method for calculating fall velocity as a function of grain size for varying grain shape (Cory shape factor) and angularity (Powers scale). We calculate w_f using this method, assuming values of Cory shape factor (0.8) and Powers scale (3.5) typical for natural gravel grains [Dietrich, 1982]. Figure 8 shows fall velocity and transport stage at the threshold of suspension plotted against grain size. For grain sizes less than about 0.0001 m, the threshold of suspension occurs at the threshold of motion (τ*/τ*_c = 1), thus silt and clay sized grains travel as wash load. The peak in the (τ*/τ*_c)_susp curve is due to the change in flow around the falling grain, from hydraulically smooth flow for fine grains to hydraulically rough flow for coarse grains.

3.5. Nondimensional Saltation Trajectory

[35] The ratio of hop length to hop height is a nondimensional measure of the shape of the saltation trajectory. We can use this ratio to test whether equations (12) and (14) fit not only the individual components of the trajectory but the overall trajectory shape as well. We do this by comparing the trajectory shape predicted by the ratio of equations (12) and (14)

to the subset of data where measurements of hop length and height were taken from the same hop event. This is shown in Figure 9, along with the best fit log-log linear line through the data pairs. We conclude that equations (12) and (14) represent the saltation trajectory shape well. Importantly, note that saltation trajectories grow more elongated with increasing excess shear stress; thus saltating grains impact less frequently per unit distance traveled as excess shear stress increases.

3.6. Vertical Impact Velocity (w_si)

[36] No direct measurements of the vertical component of particle impact velocity (w_si) are available in the published experimental studies. However, the variation in w_si with excess shear stress should scale with the mean particle descent velocity (w_sd). We can calculate w_sd as simply the hop height (H_s) divided by the descent time (t_d)

where (t_d ≅ L_sd/u_s), and L_sd and L_su are the descending and rising portions of the hop length respectively (Figure 3). Hu and Hui [1996b] report measurements of L_su and L_s that indicate a linear relationship [Sklar, 2003], which we approximate here as

Combining equations (11), (12), (14), (17), and (18) we obtain an expression for mean particle descent velocity

Saltating grains accelerate as they fall and probably reach their maximum descent velocity at the end of the hop. From the scaled plots of individual saltation trajectories of Abbott and Francis [1977, Figure 13] and from the modeled vertical velocity reported by Wiberg and Smith [1985, Figure 15], we estimate that vertical impact velocity can be approximated as twice the mean descent velocity

3.7. Calculating Bedrock Incision Rates

[37] Using expressions for saltation hop length (equation (14)), vertical impact velocity (equation (20)), and bed load sediment transport capacity (equation (8)), we can now express the three component terms of the saltation-abrasion bedrock incision model directly in terms of excess shear stress

Multiplying these three terms, rounding the resulting exponent on excess shear stress to one decimal place (i.e., −0.52 ∼ −0.5), and simplifying we obtain a final form of the saltation-abrasion bedrock incision model

for q_s ≤ q_t and u* ≤ w_f (E = 0 for q_s > q_t or for u* > w_f). If we neglect the suspension effect term, or equivalently consider only excess shear stresses well below the threshold of suspension{i.e., 0 < (τ*/τ*_c − 1) < 10}, and substitute equation (8) for q_t and rearrange, we obtain a simpler expression for bedrock abrasion rate

where the tools effect and cover effect are represented by the first and second terms respectively, k₁ = 0.08YR_bg/k_vσ_T, k₂ = 0.014Y(R_bg)^0.5/(τ*_c)^1.5k_vσ_t²ρ_s, valid for q_s ≤ q_t.

[38] In the model results reported below, we use equation (24a) to calculate bedrock incision rates from seven input variables: channel slope (S), channel width (W), total sediment supply rate (Q_s = q_sW), discharge (Q_w), grain size (D_s), a channel roughness parameter (n), and rock tensile strength (σ_T). Average boundary shear stress (τ_b) is calculated using the hydraulic radius (r_h) in place of flow depth (H_w) to account for drag due to the channel banks

where r_h is found using the Manning equation for mean flow velocity u_w = r_h^2/3S^1/2/n and assuming a rectangular channel cross section

All other quantities (ρ_s, ρ_w, g, τ*_c, and Y) are considered constants. Note that the only remaining free parameter is k_v, the dimensionless proportionality constant relating rock resistance (ɛ_v) to elastic strain energy (β = σ_T²/2Y) in equation (4). In section 6 we use experimental measurements of bedrock wear [Sklar and Dietrich, 2001] to obtain an order of magnitude estimate for k_v of 10¹².

4. Model Results

[39] The saltation-abrasion model (-) suggests that rates of bedrock incision by saltating bed load should depend, for a given rock strength, on three principal variables: shear stress, sediment supply, and grain size. In this section, we explore the behavior of the saltation-abrasion bedrock incision model as we vary each of those three principal variables independently. We then use a nondimensional form of the model to show that the full model behavior can be depicted compactly in a single nondimensional graph. Note that in all the calculations reported below we have been careful not to violate the assumptions underlying the development of the model, for example flow depths are always greater than the grain diameter and particle Reynolds numbers remain in the hydraulically rough regime.

[40] To help guide our selection of model input values we use a particular reach of bedrock river as a reference site, a gauged reach of the South Fork Eel River in northern California (Figure 1) [Seidl and Dietrich, 1992; Howard, 1998]. This actively incising, mixed alluvial-bedrock bedded canyon river motivated many of the assumptions underlying development of this model. We use the reference site as a device to ensure a physically reasonable set of values for the variables held constant as we conduct the following model sensitivity analysis, and to provide perspective on the extent of the range of possible input values encompassed by the model parameter space.

[41] Table 3 lists the reference site values for channel geometry, roughness, rock strength, discharge, sediment supply, and grain size used as model inputs in the following sensitivity analysis. Also listed are some of the calculated hydraulic conditions, such as flow depth and mean velocity, and model output for the reference site. The choice of values of discharge and sediment supply to represent the reference site is developed fully by Sklar [2003]. The reference discharge of 39.1 m³/s has an exceedence probability of 0.013 for mean daily flow, based on analysis of a 24 year record of daily mean discharges, USGS gauge 11475500 “S.Fork Eel near Branscom, CA.” A 24 hour duration of the reference bed load sediment supply of 42.6 kg/s would represent 6% of the annual coarse sediment yield to the reach assuming an average landscape lowering rate of 0.9 mm/yr [Merritts and Bull, 1989] and that 22% of the total load is in the bed load size range.

[42] We assume here that the net effect of the full distribution of discharge and sediment supply events in transporting sediment and eroding bedrock, can be accomplished by the reference discharge and sediment supply acting over a limited duration, which in this case is approximately 22 days per year. We assume that flows too low to mobilize sediment (and thus erode bedrock) occur on the remaining 343 days. For reference site conditions, the saltation-abrasion model predicts an instantaneous bedrock incision rate of 20 mm/yr (0.056 mm/day), which is equivalent to an annual average incision rate of 0.9 mm/yr. Note that in each of the figures showing model results in this section (Figures 10, 11, 12, 13, 14, and 15), the curve passing through the reference site values is plotted with a thick solid line; an open circle indicates the reference site.

4.1. Influence of Sediment Supply

[43] Figure 10 shows how incision rate predicted by the saltation-abrasion model varies with sediment supply, holding all other input variables to constant reference site values. Also plotted are the component terms from equation (1), the volume eroded per impact (V_i), the rate of impacts by saltating grains (I_r), and the fraction of bed area exposed to abrasive wear (F_e). We find a parabolic variation in erosion rate, as expected from -. At low supply rates, incision rate grows with increasing sediment supply, reflecting the growth in impact rate. At high supply rates, incision rate decreases with increasing supply rate, due to the reduction in the fraction of bed area exposed. Erosion rate is zero when the sediment supply equals the transport capacity and the bed is fully alluviated, and in the trivial case when there is no sediment supply. The maximum incision rate corresponds to a critical level of sediment supply where the growth in impact rate is balanced by the reduction in fraction of bed area exposed. Because shear stress and grain size are not varied here, the volume eroded per impact remains constant.

[44] If we repeat this calculation for different values of channel slope, and thus different shear stresses, we generate a set of parabolic erosion rate curves, as shown in Figure 11. As shear stress increases from the threshold of motion, the peak erosion rate increases along with the sediment supply rate required to alluviate the bed and shut off incision. As transport stage rises above a value of about 10, however, the magnitude of the peak erosion rate begins to decline. This occurs because the increase in particle impact energy with greater excess shear stress is more than offset by the reduction in frequency of particle impact as the saltation hop length rapidly lengthens.

4.2. Influence of Transport Stage

[45] The dependence of erosion rate on excess shear stress is more clearly illustrated in Figure 12, where we again plot erosion rate along with the component terms from equation (1), now holding sediment supply and other variables constant. Erosion rate is zero below the transport stage required to transport the supplied sediment load. Even though transport stage exceeds the threshold of motion, the bed is assumed fully alluviated. Erosion rate increases rapidly once there is excess sediment transport capacity, reflecting the growth in both the extent of bedrock exposure in the channel bed and the particle impact energy. The erosion rate peaks due to a decline in the rate of growth in V_i and F_e, so that these two terms no longer offset the continuous reduction in impact frequency (plotted separately on a logarithmic scale) with increasing shear stress. Erosion rates go to zero for high-transport stage as the saltation hop length approaches infinity at the threshold of suspension. Note that the peak erosion rate occurs at an excess shear stress well within the range of experimental saltation data used to parameterize the model, and is not an artifact of the extrapolation we use to account for the transition from bed load to suspended transport mode (equation (14)). The peak in erosion rate with increasing shear stress (hereafter called the “transport peak” as distinguished from the “supply peak” discussed above) occurs at a relatively low transport stage in Figure 12 because sediment supply is held constant. The transport peak would occur at a higher excess shear stress if we allowed sediment supply to increase along with shear stress.

[46] Figure 13 illustrates the variation in incision rate with shear stress for various levels of sediment supply, holding the remaining variables constant. Greater rates of sediment supply require greater excess transport capacity to expose bedrock and initiate incision (i.e., for q_t > q_s). Peak erosion rates occur at greater transport stages for larger sediment supply rates. The magnitude of the peak erosion rate rises and then declines with increasing sediment supply rate. This is similar to the pattern of the erosion rate peaks shown in Figure 11, except that now the roles of sediment supply and shear stress are reversed. Erosion rate goes to zero in all cases as transport stage approaches the threshold of suspension, which for D_s = 0.06 m occurs when τ*/τ*_c ≅ 30.

4.3. Influence of Grain Size

[47] Another way to understand the variation in erosion rate with shear stress is to vary transport stage by varying grain size, holding slope and all other variables constant. This is shown for reference site conditions in Figure 14, where we again plot V_i, I_r and F_e, the component terms of equation (1), as well as the transport stage. Here high-transport stage corresponds to small grain size, because D_s is in the denominator of τ*. Erosion rate is zero for grain sizes greater than 0.065 m because transport capacity drops below sediment supply rate. Erosion rate is zero for grain sizes smaller than about 0.003 m because volume eroded per impact (V_i) and particle impact rate (I_r) are zero when grains are in suspension. For grain sizes above 0.003 m, I_r is initially high and V_i initially low because the very small mass of an individual particle results in a large number of particles, for a fixed sediment supply rate, each with very low kinetic energy. As grain size increases, V_i grows more rapidly than I_r declines, such that the product of V_i and I_r, equivalent to the net flux of kinetic energy, grows larger. Erosion rate peaks at the grain size where the rate of reduction in bed exposure (F_e) equals the rate of increase in kinetic energy flux, and declines to zero as the reduction in excess shear stress approaches the threshold of alluviation. Once F_e = 0 (when Q_s = Q_t), erosion rate is zero, even though larger grains would still be in motion and thus have a nonzero impact energy. Note that as the threshold of motion is approached the decline in particle impact rate (I_r) with increasing grain size is reversed and at the threshold of motion the impact rate is infinite. This singularity in equation (23) is physically equivalent to a transition to rolling motion, but does not affect predicted erosion rates because the threshold of alluviation always occurs at a higher excess shear stress than the threshold of motion. Figure 14 is conceptually a mirror image of Figure 12 because the direction of increase in transport stage along the horizontal axis is reversed. Note that the curve shapes are not exact mirror images because linear variation in grain size is equivalent to hyperbolic variation in transport stage.

[48] The role of grain size in controlling sediment transport capacity for a given slope and discharge and thus influencing the extent of bedrock exposure is illustrated in Figure 15, where we plot erosion rate as a function of grain size for various sediment supply rates. Irrespective of sediment supply, erosion rate is zero for grains too small to move as bed load. As in Figure 14, erosion rate rises with increasing grain size to a peak and then declines to zero as the sediment transport capacity approaches the sediment supply rate. A maximum in the peak erosion rate occurs for an intermediate level of sediment supply, where the tools and coverage effects are balanced. Note that for larger sediment supply rates, full bed alluviation occurs at smaller grain sizes, in effect narrowing the range of grain sizes capable of eroding bedrock. Figure 15 is conceptually a mirror image of Figure 13; each illustrates the effect of variable transport stage on bedrock erosion rates.

4.4. Nondimensional Framework

[49] The predicted dependence of bedrock abrasion rate on sediment supply and shear stress can be described compactly by expressing the saltation abrasion model in a nondimensional form. We define a dimensionless sediment supply (q*_s) as the sediment supply (q_s) relative to the transport capacity (q_t)

Substituting equation (27) and the bed load sediment transport relation (equation (8)) into equation (24a) and rearranging, we obtain an expression for a nondimensional bedrock incision rate (E*)

where E* = Eσ_T²/ρ_sY(gD_s)^3/2 and k₃ = 0.46R_b^3/2τ*_c. Because the ratio u*/w_f is a function of transport stage for a given grain size (Figure 8), the nondimensional erosion rate (E*) can be considered a function of just two nondimensional quantities, the relative sediment supply (q*_s) and the transport stage (τ*/τ*_c). An equivalent nondimensional functional relationship can be derived by formal dimensional analysis using the Buckingham pi theorem [Sklar, 2003].

[50] Using this nondimensional framework the saltation-abrasion model is plotted in Figure 16, which shows contour lines of equal dimensionless erosion rate (E*, contour labels × 10⁻¹⁵) as a function of relative sediment supply and transport stage. As discussed above, relative sediment supply is assumed to be equal to the fraction of bed area covered by transient alluvial deposits (i.e., q*_s = 1 − F_e). Using this framework, the saltation-abrasion model collapses to a unique functional surface for all physically reasonable combinations of the seven model inputs: sediment supply (Q_s), coarse sediment grain size (D_s), discharge (Q_w), slope (S), channel width (W), channel roughness (n), and rock tensile strength (σ_T). The surface has single peak value of E* = 7.7 × 10⁻¹⁵, which occurs at q*_s = 0.50 and τ*/τ*_c ≅ 12.

[51] Erosion rate goes to zero at each of the four boundary conditions. Along the transport stage axis, erosion rate is zero at the threshold of motion and at the threshold of suspension. Along the relative supply axis erosion rate is zero when the bed is fully alluviated and when the bed is composed of bare bedrock due to a negligible supply of tools. The peak nondimensional erosion rate occurs when the tradeoffs between the cover and tools effects, and between the energy and frequency of particle impacts, are exactly balanced.

[52] Figure 17 shows how the three curves corresponding to reference site conditions in Figures 11, 13, and 15 are, in effect, vertical slices through the nondimensional erosion function surface. Varying sediment supply while holding shear stress constant (Figure 17a, as in Figure 11) creates a vertical slice parallel to the relative supply axis. Similarly, varying transport stage while holding supply constant, either by increasing slope (Figure 17b, as in Figure 13) or by decreasing grain size (Figure 17c, as in Figure 15), creates a vertical slice through the erosion function surface that bends away from the thresholds of motion and alluviation, and toward the threshold of suspension and maximum bedrock exposure. The variable transport stage slices (Figures 17b and 17c) bend because transport stage occurs on both axes, as sediment transport capacity, in the denominator of the relative supply axis, depends on transport stage in equation (8).

[53] The inverse relationship between bedrock incision rate and excess shear stress at high-transport stages (equation (24a) and 12-16–17), which results from the tendency of saltation trajectories to grow more elongated with increasing excess shear stress (equation (16) and Figure 9), is an unanticipated feature of the saltation-abrasion model. It suggests that there is a fundamental limit to how quickly bedrock can be eroded by the process of bed load abrasion. The critical shear stress at which erosion rate reaches a maximum for variable excess shear stress corresponds to a critical relative sediment supply rate, which for low excess shear stress is approximately Q_s/Q_t = 0.25. This can be seen most clearly by considering the derivative of equation (24b) (the dimensional form of the saltation-abrasion model without the suspension term) with respect to excess shear stress, which for compactness we represent here by γ (i.e., γ = τ*/τ*_c − 1),

Equation (29) equals zero when

The critical relative sediment supply of 0.25 results from the ratio of the exponents on excess shear stress in equation (24b) (i.e., 0.5/2.0). The expression for the derivative of equation (24a), the saltation-abrasion model including the suspension term, with respect to excess shear stress, is considerably more complex. The critical relative supply is no longer a constant, but becomes an increasing function of excess shear stress.

[54] Figure 18 is a contour map of the derivative of the nondimensional erosion rate (E*) predicted by the full saltation-abrasion model (equation (24a)), with respect to transport stage, showing how the critical relative sediment supply varies between 1.0 at the threshold of suspension to 0.25 at the threshold of motion. The point (τ*/τ*_c, Q_s/Q_t) where the curves of erosion rate versus excess shear stress (Figures 13, 15, and 17) cross the curve of dE*/d(τ*/τ*_c) = 0 (Figure 18) corresponds to the maximum possible erosion rate for the given fixed input values. Figure 18 also shows that erosion rate is most sensitive to changes in excess shear stress near the threshold of alluviation (Q_s/Q_t = 1).

4.5. Implications for the Stream Power Model

[55] The nondimensional framework is useful for understanding the differences between the saltation-abrasion model and other models for river incision into bedrock incision. Here we briefly consider the widely used model [e.g., Whipple and Tucker, 1999] in which incision rate is assumed to be proportional to stream power per unit bed area (ω)

to illustrate how the saltation-abrasion model can be used to disaggregate the many factors lumped into the dimensional erosional efficiency coefficient k_sp. Several studies have attempted to calibrate k_sp at a landscape scale [e.g., Stock and Montgomery, 1999], assuming that k_sp primarily represents the role of rock resistance. However, many other factors are implicitly represented by k_sp that influence the efficiency of the erosion process for a given rate of energy expenditure, most notably sediment supply and grain size.

[56] The power per unit bed area expended by the river in eroding rock (ω_e) can be expressed as the product of erosion rate and the impact energy required to erode a unit volume of rock

We can define a dimensionless stream power erosional efficiency (e_f) as the ratio of ω_e to the total unit power expended by the flow (ω)

For the stream power model, assuming wear primarily by sediment impacts,

and for the saltation abrasion model (equation (24a))

Note that for the saltation-abrasion model, unlike the stream power model, stream power erosional efficiency is independent of rock strength; energy is dissipated by fluid friction and grain impacts in the same proportions whether the rocks are strong or weak. Equations (34) and (35) can be used to explicitly evaluate the roles of sediment supply and grain size in controlling the stream power efficiency of the abrasion process, and in influencing the value of k_sp.

[57] Figure 19 shows contours of equal stream power erosional efficiency as a function of relative sediment supply and transport stage, assuming that all erosion takes place by saltating bed load impacts. The efficiency surface peaks as expected at a relative sediment supply of 0.5, where the tools and cover effects are balanced. However, the efficiency peak occurs at τ*/τ*_c = 2.7, a low transport stage relative to the nondimensional erosion rate peak of τ*/τ*_c = 12 (Figure 16). This suggests that moderate magnitude-frequency discharges are the most efficient events, in terms of energy dissipation, for eroding bedrock by bed load abrasion.

[58] Note that the common assumption that K_sp is constant for a given rock type [e.g., Stock and Montgomery, 1999] is equivalent to assuming that efficiency should plot in this parameter space as horizontal plane. Note also that the peak efficiency is less than 1%, indicating that nearly all stream energy is lost to heat by turbulent and viscous fluid friction while only a negligible fraction of total stream power contributes to rock wear. This result is independent of grain size as well as rock strength.

5. Comparison With Experimental Data

[59] In this section we compare erosion rates predicted by the saltation-abrasion model with laboratory measurements of bedrock wear by saltating bed load in rotational bedrock abrasion mills [Sklar and Dietrich, 2001]. The primary goal of this comparison is calibrating the one free parameter in the model (k_v), the dimensionless constant in equation (4) relating the energy required to erode a unit volume of rock (ɛ_v) to the elastic strain energy of the rock at failure (σ_T/2Y). We also use the experimental data to evaluate two key assumptions: (1) that the threshold energy for detachment (ɛ_t) can be neglected and (2) that the abrasion rate of the bedrock bed becomes negligible as the transport mode shifts from bed load to suspension.

[60] The data from the bedrock abrasion mills are not appropriate, however, for fully testing the saltation-abrasion model because the experimental conditions are different from those assumed for the model in two important respects. First, the rotational geometry of flow in the abrasion mills produces a strong secondary flow, in which water descends along the walls and rises in the center of the vortex, that results in significant cross stream and vertical components of the stress field at the bedrock bed. In contrast, the model derivation above only considered the downstream stress component (τ_b). The secondary circulation sweeps sediment to the center of the bedrock disks where wedge-shaped deposits form and expand outward with increasing sediment loading [Sklar and Dietrich, 2001, Figure 4b]. Second, unlike a reach of river channel, the rotational abrasion mill is a closed system, so that partial bed coverage results from a partitioning of sediment additions into mobile and immobile fractions, rather than from a channel adjustment to maintain steady state between sediment input and output. As a result, partial bed coverage cannot be directly related to sediment flux, as the model assumes, but rather depends on the mass of sediment present in the abrasion mill.

[61] The bedrock abrasion mill data most useful for calibrating k_v are from runs with very low sediment loading, such as where a single 70 gram particle was used in testing the variation in erosion rate with rock strength, with bedrock and sediment composed of the same lithology [Sklar and Dietrich, 2001, Figure 3]. In this case, the differences in the pattern of partial bed coverage between the model and abrasion mills do not affect the results because no significant partial coverage occurred. Similarly, the low sediment loading for runs exploring the variation in erosion rate with grain size, for constant sediment mass [Sklar and Dietrich, 2001, Figure 5], can be used to consider the effects of including the threshold of suspension and neglecting the threshold of detachment in the saltation-abrasion model.

5.1. Calibrating the Rock Strength Parameter (k_v)

[62] Use of the abrasion mill data to constrain the value of k_v requires estimates of the rate and energy of grain impacts, which we estimate from observations of sediment motion made through clear-walled abrasion mills and from the predictions of the saltation-abrasion model.

[63] Because of the uncertainty in these estimates we first conducted simple impact wear tests in which we eroded bedrock disks by repeatedly dropping particles from a known height (1.0 m), in order to know the impact energy and number of impacts precisely. The impact energy of a dropped particle ɛ_i = M_pgH_d where M_p is the mass of the impacting particle and H_d is the height of drop. From equation (4) we obtain an expression for k_v for the case of dropped particles

where N_d is the total number of drops (100 in each case), M_e is the total mass of bedrock eroded, and ρ_s is the density of bedrock. Table 4a lists the conditions and results for drop tests on two artificial sandstones; both tests indicate an order of magnitude estimate of k_v ∼ 1 × 10¹².

Table 4a. K_v Calibration From Wear by Dropped Particle Experiments

σt, MPa	ρs, kg/m3	Mp, kg	Me, kg	ɛi, J	Nd	kv
1.12	2300	0.074	0.0101	0.722	100	1.32 × 1012
1.12	2300	0.074	0.0087	0.729	100	1.54 × 1012
1.12	2300	0.072	0.0138	0.706	100	9.38 × 1011
3.44	2300	0.084	0.0012	0.824	100	1.33 × 1012
3.44	2300	0.087	0.0009	0.853	100	1.84 × 1012
3.44	2300	0.087	0.0012	0.853	100	1.38 × 1012

[64] To apply the saltation-abrasion model to the experimental abrasion mill setting, we use observations of sediment motion to estimate q_s and τ_b, and rearrange equation (24a) to solve for k_v for the case of a single abrasive grain

neglecting the cover term (1 − Q_s/Q_t) because there are no other sediment grains for the single abrasive grain to interfere with. Assuming sediment motion can be represented by travel along a circular path at the mean radial distance of 0.06 m, a single grain covers a circumferential distance (L_c) of 0.38 m per revolution. With the propeller speed set to 1000 rpm, we observed 70 gram particles rotating with an average angular velocity of 70 rpm, so that average “downstream” sediment velocity u_s = 0.44 m/s. Total sediment flux Q_s = M_pu_s/L_c = 0.081 kg/s, and, for a total radial width (r_w) of 0.075 m, sediment flux per unit width q_s = Q_s/r_w = 1.08 kg/m s.

[65] A 70 gram abrasive grain has a nominal diameter of 0.038 m. We observed that grains larger than this size would either not move at all, or would only move sporadically. Grains with a diameter of about 0.045 m appeared to be closest to the threshold of motion, vibrating in place and occasionally sliding or rolling a short distance. Assuming a low value of τ*_c = 0.03 to account for the low friction angle for a single grain on a nearly planar bedrock surface [Buffington and Montgomery, 1997], and ρ_s = 2600 kg/m³, the “downstream” component of the average boundary shear stress τ_b = D_s(ρ_s − ρ_w)gτ*_c = 22 Pa. For the 0.038 m diameter grain, transport stage τ*/τ*_c = 1.2 and the ratio u*/w_f = 0.19. At this transport stage, the model prediction of downstream sediment velocity u_s = 0.49 m/s agrees well with the observed average sediment velocity. The model prediction of impact rate I_r = 90 impacts/m² s is also consistent, for the bedrock disc area of 0.03 m², with the observed frequency of 2 to 3 audible grain impacts per second.

[66] To solve for k_v using equation (37), we use the measured wear rates for the set of runs in which bedrock of various strengths was eroded by a single abrasive grain composed of the same rock type [Sklar and Dietrich, 2001, Figure 3], assuming as before a constant elastic modulus Y = 5 × 10⁴ MPa [Clark, 1966]. Calculated values of k_v average ∼3 × 10¹², and range between 1 × 10¹² and 9 × 10¹² (Table 4b), roughly consistent with an order of magnitude estimate of k_v = 1 × 10¹².

Table 4b. K_v Calibration From Single-Grain Abrasion Mill Experiments

Rock Type	σt, MPa	ρs, kg/m3	E, g/hr	E,a m/s	kv
Artificial (20:1)b	0.163	2300	215	8.69 × 10−7	4.07 × 1012
Artificial (6:1)	0.448	2300	71	2.87 × 10−7	1.63 × 1012
Artificial (4:1)	1.12	2300	5.1	2.06 × 10−8	3.63 × 1012
Sandstone	1.583	2450	2.9	1.11 × 10−8	3.37 × 1012
Graywacke	9.1	2500	0.22	8.18 × 10−10	1.38 × 1012
Limestone	9.78	2600	0.21	7.51 × 10−10	1.30 × 1012
Welded tuff	10.9	2600	0.036	1.26 × 10−10	6.23 × 1012
Quartzite	19	2600	0.008	2.86 × 10−11	9.09 × 1012
Andesite	24.4	2600	0.030	1.06 × 10−10	1.47 × 1012

• a Erosion rate units conversion made with bedrock disc area of 0.03 m².
• b Numbers in parenthesis indicate sand to Portland cement ratios for artificial sandstone.

5.2. Testing Assumptions Regarding Erosion Threshold (ɛ_t) and Hop Length (L_s)

[67] We next use the abrasion mill data to test the assumptions of a negligible threshold of erosion (ɛ_t = 0) and infinite saltation hop lengths above the threshold of suspension (equation (14)). Data from runs in which grain size was varied, holding total sediment mass and all other variables constant, are most appropriate, because no partial coverage occurred at the low sediment loading of 70 g. In the variable grain size runs limestone bedrock was not abraded by limestone gravel, but rather by quartzite gravel (0.01 m < D_s < 0.035 m), vein quartz-rich mixed pebbles (0.003 m < D_s < 0.01 m), and quartz sand (D_s < 0.002 m). Because these sediments are stronger than the limestone bedrock, we use k_v = 3 × 10¹² to account for the experimentally observed threefold increase in erosion rates when quartzite gravels were used instead of gravel composed of the same rock type as the bedrock [Sklar and Dietrich, 2001, Figure 3].

[68] Figure 20 shows that the erosion rates predicted by the saltation-abrasion model (equation (24a)) closely follow the observed variation in bedrock wear as grain size decreases from near the threshold of motion (D_s ∼ 0.035 m) to the threshold of suspension (D_s ∼ 0.001 m). Note that in contrast to Figures 14 and 15, we plot erosion rate on a logarithmic scale in Figure 20 so that the low erosion rates for sand-sized grains can be distinguished from zero. Also plotted in Figure 20 are curves for model predictions using nonzero values of the erosion threshold (ɛ_t), and for the case of ɛ_t = 0 but without the infinite hop length suspension term (i.e., equation (24b)). As ɛ_t is increased from zero by successive orders of magnitude the model predictions increasingly deviate from the measured values. A nonzero erosion energy threshold affects the predicted erosion rates most strongly for the smaller grain sizes, because impact energy depends on the fourth power of grain size (equation (21)).

[69] Erosion rates for sand sizes are greatly overpredicted by the version of the model that lacks the suspension term (i.e., equation (24b)). Neither version of the model adequately predicts erosion rates for medium sand (D_s ∼ 0.0005 m), because the model does not represent the dynamics of suspended sediment transport in which sediment is exchanged between a bed load layer and the water column. In the rotational abrasion mills, sand grains are carried down to and across the bed by the secondary circulation, resulting in some bedrock wear. Note that no wear was measured for fine sand (D_s ∼ 0.0002 m) despite the occurrence of some sand transport across the bed.

6. Discussion

[70] The saltation-abrasion model is intended to represent the dynamics of bed load abrasion for the full range of sediment transport and supply conditions. Therefore all channels incising into bedrock predominantly by bed load abrasion should plot somewhere on the nondimensional erosion function surface depicted in Figure 16. For example, the “reference” reach of the South Fork Eel River shown in Figure 1b, which is mostly mantled with coarse alluvium, should plot near the upper left corner of Figure 16, not far above the threshold of motion during incision events, with the bed mostly buried. In contrast, the reach of the Oregon Smith River shown in Figure 1a, which receives little coarse sediment supply, should plot nearer to the lower right corner, with high bedrock exposure and large excess shear stress during incision events. Because the saltation-abrasion model was developed at the timescale of individual discharge events, the dynamic variation in channel conditions and incision rate due to changes in discharge and sediment supply can be represented by shifts in position on the erosion function surface.

[71] To apply the model at much longer, geomorphically significant timescales, we need to address tradeoffs between magnitude and frequency inherent in the distribution of discharges and sediment supply events that an incising channel might experience. The most frequent discharges typically do not exceed the threshold of motion shear stress and thus move no bed load and cause no bedrock wear. Peak erosion rates predicted by the saltation-abrasion model occur at large excess shear stresses, which in most incising channels are likely to occur only during rare extreme discharges. In our selection of “reference site” conditions (Table 3) we assume for simplicity that the integrated effect of the full range of discharge and sediment supply events can be represented by a dominant discharge that moves all the annual coarse sediment supply and a fraction of the total runoff in a small fraction of the total time. In alluvial channels with wide adjacent floodplains, the geomorphically dominant events have moderate magnitude and recurrence intervals of 1–2 years [Wolman and Miller, 1960]. Baker and Kale [1998] have argued that the magnitude of the dominant discharge in bedrock channels should scale with the rock resistance, with more resistant rocks being eroded primarily by infrequent extreme events.

[72] An important implication of our model is that extreme magnitude events may not be the most effective in eroding bedrock, instead, more frequent moderate magnitude events may be responsible for the majority of bedrock wear. This is because of the reduced frequency of particle impact with increasing shear stress, the same effect that imposes an ultimate limit to the rate of erosion by saltating bed load (Figure 16). In addition, the coupling between discharge and sediment supply will strongly influence which discharges are most important in eroding bedrock. Some extreme events, like the 1964 storm of record on Redwood Creek and other N. California Coast Range rivers, result in meters of aggradation and complete burial of underlying bedrock [Lisle, 1982]. Large magnitude events can also strip away alluvial cover, or mobilize coarse alluvial armor, and erode bedrock that may have been buried for many years. Whether big storms are net importers or exporters of sediment may primarily depend on how much time has passed since the last big event, and thus on how much sediment has been produced by hillslope processes [Howard, 1998].

[73] An inherent limitation of the mechanistic approach adopted here is that we have only considered a single erosional process and thus probably underpredict total incision rate, particularly at high excess shear stresses. However, many of the effects of sediment predicted by the saltation-abrasion model influence the efficiency of other incision mechanisms as well. For example, partial bed cover by transient alluvial deposits will prevent incision by all other erosional mechanisms, with the exception of dissolution. Impacts by saltating bed load will certainly influence and in some lithologies may control the rate of plucking by enlarging and extending fractures along planes of weakness [Whipple et al., 2000a]. Increased fracture density due to multiple bed load impacts exposes fresh mineral surfaces to weathering and dissolution, enhancing rates of chemical weathering and dissolution. Bed load abrasion may also promote cavitation by creating small pits and other irregularities in the bedrock surface that favor flow separation and bubble formation at high-flow velocities [Barnes, 1956; Whipple et al., 2000a].

[74] By focusing only on abrasion by bed load and assuming that the shift to suspended transport is equivalent to grains taking saltation hops of infinite length (equation (14)), we do not account for the possibility of bed lowering by suspended load abrasion. Abrasion by suspended load has been argued to dominate wear of boulders and bedrock outcrops protruding high into the flow [Hancock et al., 1998; Whipple et al., 2000a] and the walls of narrow bedrock slots [Whipple et al., 2001]. However, the long-term rate of channel incision will be controlled by the erosion rate of the lowest portions of the channel bed, where bed load will be concentrated. The only proposed theory for bedrock erosion by suspended load [Whipple et al., 2000a, 2000b] assumes, by analogy with vertical aeolian wear profiles [Anderson and Haff, 1988], that sediment concentration, and thus particle impact rate, scales with mean downstream fluid velocity. However, as Whipple et al. [2000a] acknowledge, this holds only when an unlimited supply of suspendable sediment is available on the bed, as occurs in lowland sand bedded channels but not in incising bedrock channels where suspended sediment concentration is controlled by supply from upstream.

[75] Nonetheless, grains small enough to travel in suspension may contribute to lowering of bedrock channel beds. Turbulent bursts and sweeps, and flow separation over bed roughness elements will bring suspended grains into contact with exposed bedrock. Moreover, contrary to our assumption of infinite hop lengths for suspended sediment, some amount of suspendable sediment will remain as bed load in the near-bed boundary layer due to the exchange of sediment between bed load and suspended load that underlies conventional suspended sediment transport theory [e.g., McLean, 1992]. Although neglecting abrasion by fine sediments transported predominantly in suspension may result in underpredicting abrasion rates at very high excess shear stresses, the low erosional efficiency of fine sediments seen in bedrock abrasion mill experiments [Sklar and Dietrich, 2001] suggests that the extent of underprediction is small.

[76] The saltation-abrasion model could be improved by treating channel roughness as a dynamic variable. In all the calculations reported here we kept the roughness parameter constant for simplicity. A more realistic approach would account for the effect of form drag in reducing the shear stress available for sediment transport and bed load abrasion. The influence of grain size on skin drag should also be included, as well as the roughness effect of bed load sediment motion. Few data are available to help guide a mechanistic parameterization of bedrock channel roughness. In channels with approximately planar beds and rectangular cross sections, form drag should be at a minimum when there are few alluvial deposits due to low relative sediment supply; form drag would grow with increasing bed coverage. Other channels with highly irregular bedrock topography may be most rough at low coverage, so that form drag may reach a minimum at high bed coverage.

[77] Our assumption of a planar bedrock bed neglects the possible feedbacks between differential bedrock wear and the process of partial bed alluviation. At low relative sediment supply, planar beds may be unstable because bed load may preferentially move through local low points, increasing the impact rate and thus the erosion rate in the depression at the expense of neighboring high points. This positive feedback would grow until the form drag associated with the bed depressions reduced the shear stress available to transport bed load enough to force sediment deposition to occur. Bed load impacts would then preferentially erode the exposed higher points, limiting the growth in the small-scale relief of the bedrock bed. Temporal fluctuations in relative sediment supply would have a similar effect, shifting the focus of bed load wear between higher and lower portions of the bedrock bed. Topographic high points that protrude well into the flow, above the reach of bed load abrasion, will be eroded primarily by abrasion by suspended sediments [Whipple et al., 2000a]. Irregular bed topography will also influence the saltation impact angle, and thus the fraction of total kinetic energy that contributes to deformation wear. Thus bedrock high points that do not protrude above the zone of bed load transport will erode preferentially on their upstream faces, leading to downstream migration of rock bed forms [e.g., Tinkler and Parrish, 1998].

[78] Another important simplification in the saltation-abrasion model is our assumption that bed cover is linearly proportional to excess transport capacity. In addition to neglecting the threshold of alluviation, we also do not account for the interactions between bed sediments and the water surface at critical flow. In flume experiments with Froude number ∼1, we have observed the sudden formation of antidunes as sediment supply and partial bed coverage reach a threshold value [Sklar and Dietrich, 2002]. The large increase in form drag associated with the growth of antidunes decreases the transport capacity resulting in rapid bed alluviation without further increases in the rate of sediment supply. In a separate set of experiments [Sklar and Dietrich, 2001], we have also observed bedrock wear under a thin but complete alluvial cover, presumably because sufficient impact energy from mobile grains striking stationary grains is transmitted to the underlying bedrock. We expect that a more sophisticated treatment of the controls on partial bed cover and its insulating effects might alter the shape of the nondimensional erosion function surface (Figure 16), but would not alter the fundamental constraints that the cover and tools effects impose on rates of bedrock erosion.

[79] Perhaps the most important simplifying assumption in the saltation-abrasion model is our use of a uniform grain size to represent the wide distribution of grain sizes supplied to incising channels. Grain size has a profound effect on predicted incision rates, through (1) the threshold of motion, (2) by scaling the effect of shear stress (i.e., nondimensional transport stage), and (3) in influencing the height of the erosion function surface (Figure 16). The grain size distribution of sediment supplied to the channel network will strongly influence the extent of alluvial cover and the availability of bed load abrasive tools. As Howard [1980] has shown for coarse alluvial channels, the critical threshold of motion is set by coarsest grain size fraction supplied in concentrations sufficient to dominate channel spanning deposits. The grain size distribution will largely determine the fraction of the total load that moves as bed load and thus affects the magnitude of the sediment supply variable (Q_s).

[80] Little is known about what controls the size distribution of sediments supplied to channels by adjacent hillslopes. The initial grain size distribution is likely to depend on the bedrock lithology and its tectonic history, the climatically influenced weathering rate and style, and hillslope sediment transport processes. Once sediments reach the channel network, the distribution evolves rapidly due to selective transport [Paola et al., 1992], lateral sorting [Paola and Seal, 1995], and particle breakdown in transport [Kodama, 1994]. Downstream fining can be interrupted by resupply of coarse sediments at junctions with steep low-order tributaries [Rice, 1998].

7. Conclusion

[81] We have developed a mechanistic model for bedrock incision by saltating bed load using well established empirical relationships for impact wear of brittle materials and the motion of saltating bed load grains. The model predicts that maximum incision rates occur at moderate sediment supply rates, relative to sediment transport capacity, because of a tradeoff between the availability of abrasive tools and the partial alluviation of the bed. Incision rates also peak at intermediate levels of excess shear stress due to a tradeoff between particle impact energy and impact frequency. The predicted bed load abrasion rate can be expressed in a compact nondimensional form as a function of excess shear stress and sediment supply relative to sediment transport capacity. We nondimensionalize erosion rate by grain size, and by rock resistance, which scales with the square of rock tensile strength. In this parameter space, the saltation-abrasion model defines a unique functional surface bounded by four essential constraints on incision rate imposed by sediment supply, grain size and the dynamics of sediment motion. The rate of bed load abrasion will be zero at shear stresses below the threshold of motion and above the threshold of suspension, and when the supply of coarse sediment is either negligible or exceeds the sediment transport capacity.

[82] The saltation-abrasion model illustrates the benefits of building geomorphic process models from the detailed mechanics of individual erosional mechanisms. Directly coupling bed load transport with bedrock wear by saltating grain impacts leads to several important insights: (1) grain size is a key control on bedrock incision rate because the threshold of motion introduces a minimum shear stress required to erode bedrock, and because the erosional efficiency of a given shear stress depends on the grain size; (2) incision rate is most sensitive to changes in shear stress when sediment supply is nearly equal to sediment transport capacity because small increases in shear stress can produce large changes in the extent of bedrock exposure (Figure 12); and (3) there is an upper limit to the rate at which bed load abrasion can erode bedrock due to the reduced frequency of bed load impacts with increasing shear stress (Figure 16), which suggests that other erosional mechanisms may be dominant under tools-poor or high excess shear stress conditions.

[83] Several priorities for future research emerge from this work, including understanding: (1) the controls on the grain size distribution supplied to river channels; (2) how coarse sediment is transported through bedrock bedded channels; (3) the role of the grain size distribution and the bedrock bed roughness in influencing the formation and growth of alluvial deposits; and (4) how to integrate incision processes other than abrasion by bed load into a comprehensive mechanistic model for fluvial erosion of bedrock.