Macro-roughness model of bedrock–alluvial river morphodynamics

The 1-D saltation–abrasion model of channel bedrock incision of Sklar and Dietrich (2004), in which the erosion rate is buffered by the surface area fraction of bedrock covered by alluvium, was a major advance over models that treat river erosion as a function of bed slope and drainage area. Their model is, however, limited because it calculates bed cover in terms of bedload sediment supply rather than local bedload transport. It implicitly assumes that as sediment supply from upstream changes, the transport rate adjusts instantaneously everywhere downstream to match. This assumption is not valid in general, and thus can give rise to unphysical consequences. Here we present a unified morphodynamic formulation of both channel incision and alluviation that specifically tracks the spatiotemporal variation in both bedload transport and alluvial thickness. It does so by relating the bedrock cover fraction to the ratio of alluvium thickness to bedrock macro-roughness, rather than to the ratio of bedload supply rate to capacity bedload transport. The new formulation (MRSAA) predicts waves of alluviation and rarification, in addition to bedrock erosion. Embedded in it are three physical processes: alluvial diffusion, fast downstream advection of alluvial disturbances, and slow upstream migration of incisional disturbances. Solutions of this formulation over a fixed bed are used to demonstrate the stripping of an initial alluvial cover, the emplacement of alluvial cover over an initially bare bed and the advection–diffusion of a sediment pulse over an alluvial bed. A solution for alluvial–incisional interaction in a channel with a basement undergoing net rock uplift shows how an impulsive increase in sediment supply can quickly and completely bury the bedrock under thick alluvium, thus blocking bedrock erosion. As the river responds to rock uplift or base level fall, the transition point separating an alluvial reach upstream from an alluvial–bedrock reach downstream migrates upstream in the form of a “hidden knickpoint”. A tectonically more complex case of rock uplift subject to a localized zone of subsidence (graben) yields a steady-state solution that is not attainable with the original saltation–abrasion model. A solution for the case of bedrock–alluvial coevolution upstream of an alluviated river mouth illustrates how the bedrock surface can be progressively buried not far below the alluvium. Because the model tracks the spatiotemporal variation in both bedload transport and alluvial thickness, it is applicable to the study of the incisional response of a river subject to temporally varying sediment supply. It thus has the potential to capture the response of an alluvial–bedrock river to massive impulsive sediment inputs associated with landslides or debris flows. Published by Copernicus Publications on behalf of the European Geosciences Union. 114 L. Zhang et al.: Macro-roughness model of bedrock–alluvial river morphodynamics


Introduction
The pace of river-dominated landscape evolution is set by the rate of downcutting into bedrock across the channel network. The coupled process of river incision and hillslope response is both self-promoting and self-limiting (Gilbert, 1877). Although there are multiple processes that can lead to erosion into bedrock, we here focus on incision driven by abrasion of a bedrock surface as moving particles collide with it. Low rates of incision entail some sediment supply from upstream hillslopes, which provides a modicum of abrasive material in river flows that further facilitates bedrock channel erosion. Faster downcutting leads to higher rates of hillslope sediment supply, boosting the concentration of erosion "tools" and bedrock wear rates, but also leading to greater cover of the bedrock bed with sediment (Sklar and Dietrich, 2001, 2006Turowski et al., 2007;Lamb et al., 2008;Turowski, 2009). Too much sediment supply leads to choking of the channels by alluvial cover and the retardation of further channel erosion (e.g., Stark et al., 2009). This competition between incision and sedimentation leads long-term eroding channels to typically take a mixed bedrock-alluvial form in which the pattern and depth of sediment cover fluctuate over time in apposition to the pattern of bedrock wear.
Theoretical approaches to treating the erosion of bedrock rivers have shifted over recent decades (see Turowski, 2012, for a recent review). The pioneering work of Howard and Kerby (1983) focused on bedrock channels with little sediment cover; it led to the detachment-limited model of Howard et al. (1994), in which channel erosion is treated as a power function of river slope and characteristic discharge, and the "stream-power-law" approach, in which the powerlaw scaling of channel slope with upstream area underpins the way in which landscapes are thought to evolve (Whipple and Tucker, 1999;Whipple, 2004;Howard, 1971, foreshadows this approach). At the other extreme, sediment flux came into play in the transport-limited treatment of mass removal from channels of, for example, Smith and Bretherton (1972), in which no bedrock is present in the channel and where the divergence of sediment flux determines the rate of lowering. Whipple and Tucker (2002) blended these approaches, and imagined a transition from detachment limitation upstream to transport-limited behavior downstream. They also discussed, in the context of the stream-powerlaw approach, the idea emerging at that time (Sklar and Dietrich, 1998) of a "parabolic" form of the rate of bedrock wear as a function of sediment flux normalized by transport capacity. Laboratory experiments conducted by Sklar and Dietrich (2001) corroborated this idea, and they led to the first true sediment flux-dependent model of channel erosion of Dietrich (2004, 2006). This saltation-abrasion model was subsequently extended by Lamb et al. (2008) and Chatanantavet and Parker (2009). It was explored experimentally by Chatanantavet and Parker (2008) and Chatanantavet et al. (2013); evaluated in a field context by Johnson et al. (2009), Chatanantavet and Parker (2009), Hobley et al. (2011) and Turowski et al. (2013); adapted to treat alluvial intermittency by Lague (2010); and given a stochastic treatment by Turowski et al. (2007), Turowski (2009) and Lague (2010), the latter of whom introduced several new elements. Howard (1998) presents an alternative formulation for incision that relates bedrock wear to the thickness of alluvial cover rather than sediment supply, in a form that can be thought to be a predecessor of the present work.
At the heart of their saltation-abrasion model lies the idea of a cover factor p corresponding to the areal fraction of the bedrock bed that is covered by alluvium (Sklar and Dietrich, 2004). This bedrock bed is imagined as a flat surface on which sediment intermittently accumulates and degrades during bedload transport over it. The fraction of sediment cover is assumed to be a linear function of bedload transport relative to capacity. Bedrock wear takes place when bedload clasts strike the exposed bedrock. In the simplest form of the saltation-abrasion model, the subsequent rate of bedrock wear is treated as a linear function of the impact flux and inferred to be proportional to the bedload flux, which leads to the parabolic shape of the cover-limited abrasion curve.
The saltation-abrasion model is considerably more sophisticated and flexible Dietrich, 2004, 2006) than this sketch explanation can encompass. It does, however, have three major restrictions. First, it is formulated in terms of sediment supply rather than local sediment transport. The model is thus unable to capture the interaction between processes that drive evolution of an alluvial bed and those that drive the evolution of an incising of bedrock-alluvial bed. Second, for related reasons, it cannot account for bedrock topography significant enough to affect the pattern of sediment storage and rock exposure. Such a topography is illustrated in Fig. 1 for the Shimanto River, Japan. Third, it is designed for quasi-steady conditions, and thus cannot account for the effects of cyclic variation in sediment supply on channel development downstream of the point of sediment supply.
Here we address all three of these issues in a model that allows both alluvial and incisional processes to interact and coevolve. We do this by relating the cover factor geometrically to a measure of the vertical scale of elevation fluctuations of the bedrock topography, here called macroroughness, rather than to the ratio of sediment supply rate to capacity sediment transport rate. Our model encompasses downstream-advecting alluvial behavior (e.g., waves of alluvium), diffusive alluvial behavior and upstream-advecting incisional behavior (e.g., knickpoint migration). In order to distinguish between the model of Dietrich (2004, 2006) and the present model, we refer to the former as the CSA (Capacity-based Saltation-Abrasion) model, and the latter as the MRSAA (Macro-Roughness-based Saltation-Abrasion-Alluviation) model. We point out here that the first and third issues indicated above have also been addressed by Lague (2010), although in a substantially different way than a b Figure 1. Views of the Shimanto River, a mixed alluvial-bedrock river in Shikoku, Japan. (a) Upstream view at low stage. (b) Macroscopic roughness of the bed and alluvial patches. Channel width is about 100 m. presented here. The notation used in this paper is defined in Table A1.

Capacity-based Saltation-Abrasion (CSA) geomorphic incision law and its implications for channel evolution: upstream-migrating waves of incision
2.1 CSA geomorphic incision law Dietrich (2004, 2006) present the following model, referred to here as the Capacity-based Saltation-Abrasion (CSA) model, for bedrock incision in mixed bedrockalluvial rivers transporting gravel. Defining E as the vertical rate of erosion into bedrock, q a as the volume gravel transport rate per unit width (specified in their model solely in terms of a supply, or feed rate q af ) and q ac as the capacity volume gravel transport per unit width such that q a < q ac , where β is an abrasion coefficient with dimension L −1 . By introducing a cover factor parameter p, this equation can be rewritten as This cover factor p is defined (Sklar and Dietrich, 2006) as the areal fraction of bedrock surface covered with alluvium and is given by where this fraction is calculated by averaging over a window larger than a characteristic macroscale of bedrock elevation variation. We refer to this formulation for cover factor p as "capacity based" because Eq.
(2) dictates that p is determined in terms of the ratio of sediment supply to its capacity value in the CSA model. In the above formulation, it is assumed that the gravel transport rate q a over a bedrock surface can be estimated by simply multiplying the capacity rate q ac by the areal cover fraction p. While this is the simplest first-order assumption, it should be recognized that the roughness of the bedrock itself can change the flow resistance, leading to a relationship that is more complex than Eq. (1b) (Inoue et al., 2014;Johnson, 2014).
Before introducing the relation of Sklar and Dietrich (2006) for abrasion coefficient β, it is of value to provide an interpretation for this parameter not originally given by Dietrich (2004, 2006), but which plays a useful role in the analysis below. The abrasion coefficient has a physical interpretation in terms of Sternberg's law (Sternberg, 1875) for downstream diminution of grain size (Parker, 1991(Parker, , 2008Chatanantavet et al., 2010). The analysis leading to this interpretation is given in Appendix A; salient results are summarized here. Consider a clast of material that is of identical rock type to the bedrock being abraded. Sternberg's law is where D is gravel clast size, D u is the upstream value of D, x is downstream distance and α d is a diminution coefficient. If all diminution results from abrasion, α d is related to β by In the case of constant β, and therefore constant α d , the distance L half for such a clast to halve in size is given by This interpretation of abrasion coefficient β in terms of diminution coefficient α d allows for comparison of the experimental results of Sklar and Dietrich (2001) with values of α d previously obtained from abrasion mills (Parker, 2008: see Fig. 3-41 therein;Kodama, 1994). The relations of Dietrich (2004, 2006) to compute β and q ac can be cast in the following form: In the above relations, D corresponds to the characteristic size of the gravel clasts that are effective in abrading the bedrock; ρ g is the material density of the grains; R is their submerged specific gravity (∼ 1.65 for quartz); g is gravitational acceleration, τ * is the dimensionless Shields number of the flow; τ * c is the threshold Shields number for the onset of significant bedload transport; α a and n a denote, respectively, a relation-specific dimensionless coefficient and an exponent; v f is the fall (settling) velocity corresponding to grain size D, Y is the bedrock modulus of elasticity; σ t is the rock tensile strength; and k is a dimensionless coefficient of the order of 10 −6 . (In the above two relations and the text, several misprints in Dietrich (2004, 2006) have been corrected on the advice of the authors.) Equation (5c) corresponds to the bedload transport relation of Fernandez Luque and van Beek (1976) when α a = 5.7 and n a = 1.5; Dietrich (2004, 2006) used this relation with the assumed value τ * c = 0.03. The relations above define a 0-D formulation. It must be augmented with other parameters and relations, including channel width, relations for hydraulics, quantification of flow discharge or flow duration curve, etc., to allow application at the river reach scale.
It is useful to cast Eq. (5a) in the form where β ref is a reference value of β, either computed from known values of the parameters Y , k, σ t , R f , etc., or estimated indirectly.

Embedding of CSA into a model of bedrock surface evolution
A relation for the evolution of bedrock surface elevation η b is obtained by substituting the CSA geomorphic law for incision of Eq. (1b) into a simplified 1-D mass conservation equation for bedrock material subjected to piston-style rock uplift or base level fall (Sklar and Dietrich, 2006): Here t denotes time, υ denotes the relative vertical velocity between the rock underlying the channel (which is assumed to undergo no deformation) and the point at which base level is maintained, and I denotes a flood intermittency factor to account for the fact that only relatively rare flow events are likely to drive incision (Chatanantavet and Parker, 2009). Also, I is assumed to be a prescribed constant; a more generalized formulation for flow hydrograph is given in Sklar and Dietrich (2006) and DiBiase and Whipple (2011). In interpreting Eq. (6a), it should be noted that υ denotes a rock uplift rate (in the sense of England and Molnar, 1990) for the case of constant base level, or equivalently a rate of base level fall for rock undergoing neither uplift nor subsidence. Below we use the term "rock uplift" as shorthand for the relative vertical velocity between the rock and the point of base level maintenance. Substituting Eq. (1b) into Eq. (6a) yields (Sklar and Dietrich, 2006)

Character of the CSA model: upstream waves of incision
The MRSAA model (introduced below) has several new features as compared to CSA. These are best illustrated by first characterizing the mathematical nature of CSA in the context of Eq. (6). Let denote the streamwise bedrock-surface slope. Reducing Eq. (6b) with Eq. (7) the CSA model of Eq. (1) reveals itself as a nonlinear kinematic wave equation with a source term: Here c b denotes the wave speed associated with bedrock incision. The form of Eq. (8a) dictates that disturbances in bedrock elevation always move upstream. We will see later that these disturbances can take the form of upstream-migrating knickpoints (e.g., Chatanantavet and Parker, 2009). Any solution of Eqs. (8a) and (8b) subject to the cover relation of Eq. (2) requires specification of a flow model. In mountain streams, backwater effects are likely to be negligible (e.g., Parker, 2004). The normal (steady, uniform) flow assumption allows for simplification. Let Q denote water discharge during (morphodynamically active) flood flow taking place with intermittency I , and let H denote flood depth and g denote acceleration due to gravity. Momentum and mass balance take the forms where τ is boundary shear stress at flood flow, U is the corresponding mean flow speed, B is channel width and ρ is water density. The dimensionless Shields number τ * and dimensionless Chézy resistance coefficient Cz are defined as As shown in Parker (2004) and Chatanantavet and Parker (2009), reducing Eqs. (7), (9) and (10) yields the following relations for H and τ * : A comparison of Eqs.
(2), (5c) and (11b) indicates that even for constant values of other parameters, the functional forms for q ac and thus p are such that c b is in general a nonlinear function of S b = −∂η b /∂x.

Limitations of the CSA model
The CSA model Dietrich, 2004, 2006) was a major advance in the analysis of bedrock incision due to abrasion because it (a) accounts for the effect of alluvial cover and tool availability on the incision rate through the term p(1 − p) in Eq. (1b) and (b) provides a physical basis for incision due to abrasion as gravel clasts collide with the bedrock surface. The CSA model been used, modified, adapted and extended by a number of researchers (Crosby et al., 2007;Lamb et al., 2008;Chatanantavet and Parker, 2009;Turowski, 2009;Lague, 2010). The model does, however, have a significant limitation in that it specifically does not include either alluvial morphodynamics or the morphodynamics of transitions between bedrock and alluvial zones. Here we study this limitation, and how to overcome it, in terms of the highly simplified configuration of a reach (HSR, highly simplified reach) with constant width; fixed, non-erodible banks; constant water discharge; and sediment input only from the upstream end. For simplicity, we also neglect abrasion of the gravel itself, so that grain size D is a specified constant. (This condition, while introduced arbitrarily here, can be physically interpreted in terms of clasts that are much more resistant to abrasion than the bedrock.) The means to relax these constraints is available (e.g., Chatanantavet et al., 2010;DiBiase and Whipple, 2011), and indeed many of them have been implemented in the SSTRIM model of Lague (2010). Such a relaxation, however, obscures the first-order physics underlying the rich patterns of interaction between completely and partially alluviated conditions illustrated herein.
In the CSA model, the bedload transport rate q a is specified as a "supply". That is, the bedload transport rate is constrained so that it cannot change in the downstream direction, and is always equal to the bedload feed rate (supply) q af at the upstream end. When the feed rate q af increases, q a must increase simultaneously everywhere. That is, a change in bedload supply is felt instantaneously throughout the entire reach, regardless of its length.
We illustrate this behavior in Fig. 2. The reach has length L. The gravel feed rate at x = 0 follows a cyclic "sedimentograph" (in analogy to a hydrograph) with period T = T h + T l , in which the sediment feed rate has a constant high feed rate q afh for time T h , and a subsequent constant low feed rate q afl for time T l . According to the CSA model, at x = L corresponding to the downstream end of the reach, the temporal variation in bedload transport rate must precisely reflect the feed rate. That is, the model was not designed to route sediment in the downstream direction.
In a more realistic model, the effect of a change in bedload feed rate q af would gradually diffuse and propagate downstream, so that the bedload transport rate at the downstream end of the reach would show more gradual temporal variation. This effect is illustrated in Fig. 2. This same diffusion and propagation can be expected in the cover fraction p, which in general should vary in both x and t. The change in cover fraction in turn should affect the incision rate as quantified in Eq. (1a). To capture this effect, however, Eq. (1b) must be coupled with an alluvial formulation that routes sediment downstream over the bedrock.
A second limitation concerns alluviation of the bedrock surface. Consider a wave of sediment moving over this surface, as shown in Fig. 3. We characterize the vertical scale of the geometric roughness of the bedrock surface (as seen in Fig. 1) in terms of a vertical macro-roughness L mr . For simplicity, we assume that the bottom of the bedrock relief has a specific elevation η b , an assumption that will be relaxed later in favor of a probabilistic formulation. The alluvial thickness above this basal elevation η a represents an average value of local bed elevation over an appropriately defined window. It can be seen in Fig. 3 that the surface undergoes both partial (η a < L mr ) and then complete (η a ≥ L mr ) alluviation, only to be excavated later as the wave passes through.
Bed elevation η is given as Figure 3 shows that, in the case of complete alluviation, the elevation of the bed η can be arbitrarily higher than the elevation η b of the bedrock, the difference between the two  Here η b denotes the elevation of the bottom of the bedrock relief (black line), L mr denotes the bedrock macroroughness thickness (dashed black line), η a denotes the thickness of the alluvial cover (which may be less than or greater than L mr ) and η = η b + η a denotes the elevation of the top of the alluvium. corresponding to the thickness η a . The CSA model does not describe the variation in bed elevation η when the bed undergoes transitions between partial and complete alluviation; it simply infers that incision is shut down by the complete alluvial cover.
The goal of this paper is the development and implementation of a model that overcomes these limitations by (a) capturing the spatiotemporal coevolution of the sediment transport rate, alluvial cover thickness and bedrock incision rate, and (b) explicitly enabling spatiotemporally evolving transitions between bedrock-alluvial morphodynamics and purely alluvial morphodynamics. The form of the model presented here is simplified in terms of the HSR outlined above, including a constant-width channel and a single sediment source upstream.

Formulation for alluvial sediment conservation and cover factor
The geomorphic incision law of the MRSAA model is identical to that of CSA, i.e., Eq. (1b). The essential differences are contained in (a) a formulation for the cover factor p that differs from Eq. (2) and (b) the inclusion of alluvial morphodynamics in a way that tracks the spatiotemporal evolution of the bedload transport rate, and allows for smooth spatiotemporal transitions between the bedrock-alluvial state and the purely alluvial state. The specific case we consider here is one for which (a) the bedrock surface is rough in a hydraulic sense (as opposed to a hydraulically smooth or transitional surface; see Schlichting, 1979), and (b) the characteristic vertical scale of bedrock elevation fluctuation about a mean value based on an appropriately defined window, here denoted as the macro-roughness L mr of the bedrock, is large compared to the characteristic size of the clasts constituting the alluvium. We use the term "macro-roughness" so as to clearly distinguish it from hydraulic roughness, which is specifically defined in terms the logarithmic velocity profile. Inoue et al. (2014) introduced the terms "clast-rough" and "clast-smooth", the former referring to a bedrock surface roughness that is large compared to the characteristic size of the alluvium, and the latter referring to a bedrock surface macro-roughness that is small compared to the size of the alluvium. Here we specifically consider the clast-rough case.
We formulate the problem by considering a conservation equation for the alluvium, in standard Exner form, appropriately adapted to include below-capacity transport over a non-erodible surface. The first model of this kind is due Figure 4. Illustration of the statistical structure or local hypsometry of the bedrock surface topography (dark-grey line). Here z denotes an elevation above some arbitrary datum (dark-red line) deep in the bedrock and p(z ) is the probability (green line, marker) that a point at elevation z is in water or alluvium rather than bedrock, i.e., above the local bedrock surface. The effective "bottom" of the bedrock relief is located at elevation z 0 = η b , where p takes an appropriately selected low value p 1 (e.g., 0.05); the effective "top" of the bedrock relief is located at elevation z 1 , where p takes an appropriately selected high value p 1 (e.g., 0.95); and the macro-roughness L mr is given as z 1 − z 0 . The coordinate z = z − z 0 is referenced to the effective "bottom" of the bedrock relief.
to Struiksma (1999), and further progress has been made by Parker et al. (2009Parker et al. ( , 2013, Izumi and Yokokawa (2011), Izumi et al. (2012), Tanaka and Izumi (2013) and . These models are expressed in continuous form; Lague (2010) presents a discrete version based on a series of reaches of finite length that allows for generalization to a continuous form. None of the above models is specifically designed to handle the clast-rough case, in particular that shown in Fig. 1, where the elevation of the bedrock roughness has a random element. Here we handle the clast-rough case by first characterizing the statistical nature of the bedrock surface alone. As noted in Fig. 4, z denotes elevation above an arbitrary datum deep in the bedrock, and p(z ) denotes the probability that a point located at elevation z is located in alluvium or water rather than bedrock. Conversely, 1 − p denotes the probability that a point at elevation z is in bedrock (rather than water or alluvium above). As seen on the right-hand side of the figure, p(z ) → 0 as z → −∞ (pure bedrock) and p(z ) → 1 as z → +∞. This statistical structure function (a hypsometric curve for local bedrock topography) which we use here to characterize bedrock elevation fluctuations is analogous to that used in Parker et al. (2000) for alluvial beds. It should be noted that "−∞" is shorthand for "far below the bedrock surface" and "+∞" is shorthand for "far above the bedrock surface". It should also be emphasized that the tilde in the parameter p indicates it is not a cover factor, but rather a statistical parameter referring to the bedrock relief itself.
In such a statistical formulation, bedrock relief has neither a precise "bottom" nor a precise "top". Rather, the "bottom" and "top" of the bedrock topography, as well as the macroroughness L mr , are here defined in a statistical sense. This  3) and (4), z is the elevation above the effective "bottom" of the bedrock relief, L mr is macro-roughness height, η b is elevation at base of the bedrock and η a is the thickness of alluvium (yellow fill). The diagram shows alluvium-water interfaces (blue line) that are at spatially constant elevations. This is for illustrative purposes only; the interfaces should instead be those that result from averaging over an appropriate spatial window, i.e., alluvial fill levels are expected to vary from one pocket in the bedrock relief to another.
can be done using moments or exceedance probabilities; here we use the latter. Let p 0 denote some low reference value of p (e.g., p 0 = 0.05, or deep into the bedrock relief) and p 1 denote a corresponding high reference value of p (e.g., p 1 = 1 − p 0 = 0.95, or near the upper portion of the bedrock relief), and z 0 and z 1 denote the corresponding bed elevations. An effective "base" of the bedrock relief can be set at z 0 , a macro-roughness height L mr defined as and an effective "top" of the bedrock specified as z 0 + L mr . The clast-rough condition considered here satisfies the constraint that L mr /D 1. The problem can now be rephrased in terms of a vertical coordinate z with its origin located at the effective bottom of the bedrock topography: As noted in Fig. 4, the statistical variable p(z) with shifted coordinate z now has the following properties: Here we define the thickness of the alluvial cover η a as the elevation difference between the locally averaged top of the alluvium and the elevation z = 0. The cover fraction p associated with any alluvial thickness η a relative to the macroroughness height L mr is then given as

Exner equation of alluvial sediment conservation over a bedrock surface
The alluvial sediment is taken to have constant porosity λ.
As illustrated in Fig. 5, the volume of alluvial sediment per unit area between elevations z and z + z is (1 − λ) p(z) z, and the corresponding volume bedload transport rate per unit width q a is estimated as p q ac , where again, according to Eq. (16), p = p(η a ).
For the case of sediment of constant density, the Exner equation for mass balance of alluvial sediment can be expressed as where the factor I accounts for the fact that morphodynamics are active only during floods. Reducing Eq. (17) using Leibniz's rule, The above formulation for the conservation of alluvial sediment over a bedrock surface differs in one essential way from the earlier forms due to Struiksma (1999), Parker et al. (2009Parker et al. ( , 2013, Izumi and Yokokawa (2011), Izumi et al. (2012), Tanaka and Izumi (2013) and Zhang et al. (2013). Specifically, in Eq. (18), the cover fraction p is present on the lefthand side of the equation as well as the right-hand side. It is shown below that this feature dictates a strong nonlinearity in the speed of propagation of alluvial waves over a bedrock surface, such that wave speed increases with decreasing wave amplitude. This feature is specifically captured by means of Leibniz's rule, as implemented between Eqs. (17) and (18). The combination of Eqs. (6b) and (18) delineates a formulation encompassing both mixed bedrock-alluvial rivers and alluvial rivers.

Closure model for cover relation
In the present formulation, the cover fraction p is free to vary in both x and t, i.e., p = p(x, t). In order to complete the problem, however, it is necessary to specify a closure model for p. We characterize the local variation in bedrock elevation in terms of the macro-roughness, i.e., the vertical length scale L mr of Fig. 4. Here we seek a formulation that averages over a window capturing a statistically relevant sample of this local variation. In general, we assume a cover relation that characterizes to what extent the alluvial cover "drowns" the bedrock roughness elements. More specifically, we assume the form The precise details of the relation can be expected to vary from case to case, but the overall characteristics that we hypothesize are illustrated in Fig. 6a-c. It is seen therein that p takes the residual value p = p 0 (e.g., 0.05) at χ = 0, increases monotonically to p = p 1 (e.g., 0.95) at χ = 1, and then takes the asymptotic value p → 1 as χ becomes sufficiently large. The first of these conditions corresponds to a bedrock surface that is bare of alluvium except in deep pockets of the macro-roughness elements, the second to a bedrock surface that is nearly completely alluviated but with some parts of the macro-roughness elements exposed, and the third to a bedrock surface that is deeply alluviated. Note that the cover relation of Fig. 6 and Eq. (19) is based on the macro-roughness height scale L mr rather than the transport capacity q ac of Eq. (2). This is the motivation for referring to the new model presented here as the Macro-Roughness-based Saltation-Abrasion-Alluviation (MRSAA) model.
In applying the MRSAA model to general cases, it is useful to delineate the simplest functional form for the closure relation for cover fraction that satisfies the constraints of Eq. (19) and Fig. 6. This relation is the piecewise-linear form or rephrasing to emphasize the dependence of p on the thickness of alluvial cover η a , The above relation is illustrated in Fig. 7a with the sample evaluations Equations (20) and (21) (2008), but with a specific focus on various forms of macro-roughness that mimic those in the field. The form for the derivative of Eq. (20) with respect to χ , which is given below, will prove useful in succeeding analysis.
Our model is specifically meant to apply to the case for which the characteristic size of the roughness elements is large compared to the size of clasts transported as bedload. With this in mind, it should be noted that in Fig. 6, no bed elevation variations are shown over parts of the bed that are covered with alluvium. This is done only for simplicity, and reflects  Figure 6. Illustration of the MRSAA model relation between areal fraction of alluvial cover of bedrock p(χ) (green curves) and χ = η a /L mr , where η a is the thickness of alluvium (yellow fill) and L mr is the macro-roughness height, for (a) low cover, (b) intermediate cover (c) complete alluviation above the top of the bedrock. The diagrams show alluvium-water interfaces at spatially constant elevations (blue lines). This is for illustrative purposes only: the interfaces should instead be those that result from averaging over an appropriate spatial window. the condition that in the clast-rough case considered here, grain size is small compared to macro-roughness height. Figure 6 also contains another simplification, in that all pockets are assumed to be filled to the same level by alluvium. While this condition is not likely to be true at the local scale, it is a reasonable first approximation when averaging over an appropriately defined window.
The formulation presented here has an obvious limitation. Since it is a 1-D expression of sediment conservation over a bedrock surface, it cannot capture 2-D variation, which will result in a more complex pattern than that shown in Fig. 6, and in particular will provide more connectivity between adjacent pockets. This two-dimensionality is known to have an effect on the pattern of incision, as illustrated by Johnson and Whipple (2007). The extension of the formulation to the 2-D case represents a future goal; some relevant comments can be found in the "Discussion" section (Sect. 8).
Sections 3.4 and 3.5, immediately below focus, on the mathematical interpretation of the MRSAA problem in terms of diffusion and wave characteristics. The reader whose primary interest is in applications may jump directly to Sect. 3.6, with the two exceptions of Eqs. (27) and (28) in Sect. 3.5. Equation (27) is a version of Eq. (6b) in which incisional morphodynamics are recast into a kinematic wave equation, revealing upstream-migrating waves of incision with wave speed c b . Equation (28) where The left-hand side of Eq. (23a) thus takes a kinematic wave form, such that c a is the wave speed of downstream-directed alluviation.
It is important to realize that alluvial wave speed c a is a nonlinear function of alluvial thickness η a . Using the example functions of Eqs. (20) and (22), the speed c ai of an alluvial wave of infinitesimal height η a = 0 is given from Eq. (23b) as The ratio c a /c ai of the wave speed c a at height η a to the corresponding value for an infinitesimal wave is then found to be In Fig. 7b, the ratio c a /c ai is plotted against χ = η a /L mr in accordance with the specifications of Eqs. (20), (21), and (22). It is seen therein that alluvial wave speed takes its maximum value for the limit η a → 0, and decreases to a vanishing value as the bed becomes completely alluviated (η a /L mr = 1.056). That is, waves of alluvium run fastest over a nearly bare bed, and wave-like behavior ceases to exist under conditions of complete alluviation. This latter result is in accordance with Lisle et al. (2001) and Cui et al. (2003a, b). It is of interest to inquire as to how the model would behave if the clast-rough condition, i.e., L mr /D 1, were not satisfied here. The limiting case of clast-smooth conditions would correspond to a probability distribution p(z ) in Fig. 4 that obeys a step function; p(z ) would be vanishing up to a smooth, horizontal bedrock surface, and would take the value unity above it. Consequently, the cover fraction p would converge precisely to zero as η a → 0, thus resulting, according to Eq. (23b), in an infinite speed of propagation of an alluvial wave of infinitesimal height. This is not entirely unrealistic: the physical realization would consist of clasts rolling rapidly over a smooth bed with no alluviation (Inoue et al., 2014). The presence of such a singularity, would, however, preclude the modeling of the migration of a pulse of alluvium of finite extent over the otherwise bare bed schematized in Fig. 3. In the present clast-rough formulation, the presence of deep pockets within the bedrock relief where alluvium can be stored without transport ensures that the wave speed of alluvium never displays a singularity.
The form of Eq. (23a) can be further clarified by rewriting it as where In the above relation, κ a has the physical meaning of a kinematic diffusivity. In general, q ac , q ac /S and thus κ a are nonlinear functions of S. The alluvial problem thus takes the form of a nonlinear advective-diffusive problem with a source term arising from a bedrock term.

Full
where c b denotes the speed of upstream-migrating incisional waves. Equation (18) can be cast in conjunction with Eqs. (19) and (20) as and where c a denotes the speed of downstream-migrating alluvial waves, and κ a is the kinematic diffusivity of alluvium. In this way, upstream-migrating incisional waves are combined with downstream-migrating alluvial waves and alluvial diffusion.
In MRSAA, then, the spatiotemporal variation in the cover fraction p(x, t) is specifically tied to the corresponding variation in η a through Eq. (19), e.g., the specific example of Eq. (29) above. This variation then affects incision through Eq. (27). Consider the simplified case of a wave of alluvium of finite extent illustrated in Fig. 3. There is no incision upstream of the wave because p = 0 and there is no sediment in motion over the bed. At the peak of the wave, η a > L mr , so p = 1; the bed is entirely covered with sediment, and again there is no incision. Incision can only occur on the rising and falling parts of the wave, where bedrock is partially exposed and sediment is in motion over it, i.e., 0 < p < 1. It can thus be expected that the spatiotemporal variation in cover thickness η a will affect the evolution of the long profile of an incising river that undergoes transitions between alluvial and mixed bedrock-alluvial states.

Amendment of the flow component of the MRSAA model
The flow model, and in particular Eqs. (9a) and (11), must be modified to include the alluvial formulation, so that bedrock slope S b is replaced with slope S of the top of the bed, where Thus Eqs. (9a) and (11a, 11b) are amended to and The purely alluvial case, i.e., p = 1, df/dχ = 0 and η b = const < η a in Eq. (28), results in the purely diffusional relation in which the diffusivity κ a is a function of S a = −∂η a /∂x.

How the governing equations connect to each other
In the numerical analysis below, the actual equations used to solve for morphodynamic evolution are not those of Sects. 3.4 and 3.5, but rather the primitive forms presented earlier. The unknowns to be solved are η b , η a , η, p, q ac , τ * , S, S b and S a . These nine parameters are connected to each other via nine equations, i.e., Eqs. (5c), (6b), (12), (18), (20), (30a, b, c) and (32b).

Equivalence of the MRSAA and CSA models at steady state
In the restricted case of the highly simplified reach (HSR) configuration constrained by (a) temporally constant, belowcapacity sediment feed (supply) rate q af , (b) bedload transport rate q a everywhere equal to the feed rate q af , and (c) a steady-state balance between incision and rock uplift, η a , p and S b become constant and S a vanishes, so that Eq. (28) is satisfied exactly. Equation (18) integrates to give so that η a can then be back-calculated from Eq. (20). In this case, then, the MRSAA model reduces to Eqs. (27) and (34), i.e., the CSA model.

The below-capacity steady-state case common to the CSA and MRSAA models
The steady-state form of Eq. (6) under below-capacity conditions (p < 1) can be expressed with the aid of Eq.
(2) in the form where p ss , β ss and q acss denote steady-state values of p, β and q ac , respectively. Equations (35a)-(35c) describe a balance between the incision rate and relative vertical rock velocity (e.g., constant rock uplift rate at constant base level or constant rock elevation with constant rate of base level fall). CSA and MRSAA yield the same solution for this case, which must be characterized before showing how the models differ. Equation (35a) has an interesting character. When the value of the dimensionless number ϕ exceeds unity, p falls below zero and no steady-state solution exists. Equation (35b) reveals that ϕ can be interpreted as a dimensionless rock uplift rate. Thus when the rock uplift rate is sufficiently large for ϕ to exceed unity, incision cannot keep pace with rock uplift. The model thus implicitly predicts the formation of a hanging valley. This issue was earlier discussed in Crosby et al. (2007).
In solving for this steady state, and in subsequent calculations, we use the bedload transport relation of Wong and Parker (2006a), a development and correction of the semiempirical relation of Meyer-Peter and Müller (1948), rather than the similar formulation of Fernandez Luque and van Beek (1976); in the case of the former, α a = 4, n a = 1.5 and τ * c = 0.0495. We consider two cases: one for which β ss = β is a specified constant, and one for which only a reference value β ref is specified, and β ss is computed from Eq. (5d).
In the case of a specified constant abrasion coefficient β, specification of υ, I and q af allow computation of ϕ, p ss and q acss from Eqs. (35a)-(35c). Further specification of R (here chosen to be 1.65, the standard value for quartz) and D allows the steady-state Shields number τ * ss to be computed from Eq. (5c). Steady-state bedrock slope S bss can then be computed from Eq. (11b) upon specification of flood discharge Q, Chézy resistance coefficient Cz and channel width B. In the case of β ss calculated according to Eq. (5d) using a specified reference value β ref , the problem can again be solved with Eqs. (35), (5c) and (11b), but the solution is implicit.
We performed calculations for conditions loosely based on (a) field estimates for a reach of the bedrock Shimanto River near Tokawa, Japan (Fig. 1), for which bed slope S is about 0.002 and channel width is about 100 m and (b) estimates using relations in Parker et al. (2007) for alluvial gravel-bed rivers with similar slopes, and reasonable choices for otherwise poorly constrained parameters. The input parameters, Cz = 10, Q = 300 m 3 s −1 , B = 100 m, are loosely justified in terms of bankfull characteristics of alluvial gravel-bed rivers of the same slope (Parker et al., 2007;  Wilkerson and Parker, 2011) as shown in Fig. 8a and b. The value D = 20 mm represents a reasonable characteristic size of the substrate (and thus the bedload) for gravel-bed rivers; a typical size for surface pavement is 2 to 3 times this (e.g., Parker et al., 1982). Flood intermittency I is estimated at 0.05, i.e., 18 days per year, and thus a reasonable estimate for a river subject to frequent heavy storm rainfall. Alluvial porosity is λ = 0.35. Two sediment feed rates were considered. The high feed rate was set at 3.5 × 10 5 t yr −1 , which corresponds to the following steady-state parameters at capacity conditions: Shields number τ * = 0.12, depth H = 1.5 m, steady-state alluvial bed slope S ass = 0.0026 and Froude number Fr = 0.51, where The low feed rate was set at 3.5 × 10 4 t yr −1 , corresponding to the following parameters at capacity conditions: Shields number τ * = 0.064, depth H = 2.1 m, steady-state alluvial bed slope S ass = 0.0010 and Froude number Fr = 0.32.
The value β ss = 0.05 km −1 was used for the case of a constant steady-state abrasion coefficient. This corresponds to a value of α d of 0.017 km −1 , which falls in the middle of the range measured by Kodama (1994) for chert, quartz and andesite (see Fig. 3-41 of Parker, 2008). For the case of a variable abrasion coefficient, Eq. (1a) was used with β ref set to 0.05 km −1 and τ * ref set to 0.12, i.e., the value for the high feed rate. This value of τ * ref is about 2.5 times the threshold value of Wong and Parker (2006a).
For the high feed, predicted relations for a steady-state abrasion coefficient β ss versus rock uplift rate υ are shown in Fig. 9a; the corresponding predictions for S bss versus υ are shown in Fig. 9b; the corresponding predictions for p ss and ϕ are shown in Fig. 9c. Both the cases of constant and variable β ss are shown. There are five notable aspects of these figures: (a) in Fig. 9a, the predictions for variable β ss are very similar to the case of constant, specified β ss , and indeed are nearly identical for υ ≤ 3.3 mm yr −1 (corresponding to ϕ ≤ 0.05 in Fig. 9c). (b) In Fig. 9b and c, the predictions for S bss , p ss and ϕ for variable β ss are again nearly identical to those for constant β ss , and again essentially independent of υ for υ ≤ 3.3 mm yr −1 . (c) In Fig. 9c, p ss is only slightly below unity (i.e., p ss ≥ 0.95), and ϕ ≤ 0.05 for υ ≤ 3.3 mm yr −1 ). (d) For υ > 3.3 mm yr −1 , the predictions for S bss and p ss become dependent on υ, such that S bss increases, and p ss decreases, with increasing υ. The values for constant β ss diverge from those for variable β ss , but are nevertheless close to each other up to some limiting value. (e) This limiting value corresponds to ϕ = 1 and thus p ss = 0 from Eq. (35a); larger values of ϕ lead to hanging valley formation. Here ϕ = 1 for the very high values υ = 65 mm yr −1 for constant β ss and υ = 30 mm for variable β ss .
These results require interpretation. It can be seen from Eqs. (35a-35c) that when υ/(I β ss q af ) = ϕ 1, p becomes nearly equal to unity (very little exposed bedrock), in which case q af is constrained to be only slightly smaller than q ac . From Eqs. (5c) and (11), then, S bss is only slightly above the steady-state alluvial bed slope S ass . Note that the steadystate bedrock slope decouples from rock uplift rate under these conditions: the predictions for υ = 0.2 mm yr −1 are nearly identical to this for υ = 3.3 mm yr −1 . This behavior is a specific consequence of the condition ϕ 1 corresponding to a low ratio of uplift rate to reference incision rate E ref = I β ss q af . They imply a wide range of conditions for which (a) very little bedrock is exposed, and (b) bedrock slope is independent of uplift rate.
The results for the low feed rate are very similar. The values for variable steady-state abrasion coefficient β ss differ from the constant value β ss in Fig. 10a, but this is because the constant value β ss = 0.05 was set based on the high feed rate. The results in Fig. 10b and c are qualitatively the same for Fig. 9b and c; the uplift rate below which ϕ < 0.05 is υ < 0.33 mm yr −1 for the case of constant β ss ,  Figure 9. Variation at steady state (black curves) of (a) abrasion coefficient β ss , (b) bedrock slope S bss and (c) cover fraction p ss and parameter ϕ on rock uplift or base lowering rate υ, for a high bedload feed rate of 3.5 × 10 5 t yr −1 .  Figure 10. Variation at steady state (black curves) of (a) abrasion coefficient β ss , (b) bedrock slope S bss and (c) cover fraction p ss and parameter ϕ with rock uplift or base lowering rate υ, for a low bedload feed rate of 3.5 × 10 4 t yr −1 . The cases of constant, specified β ss and S bss varying according to Eq. (34) are shown as blue curves. The vertical dashed lines denote the incipient conditions for the formation of a hanging valley. The predictions are the same for the CSA and MRSAA models. and υ < 0.73 mm yr −1 for the case of variable β ss . The critical value of υ beyond which a hanging valley forms is υ ≥ 6.8 mm yr −1 for constant β ss and υ ≥ 7.1 mm yr −1 for variable β ss .
The lack of dependence of steady-state bedrock slope S bss on rock uplift rate υ below a threshold value for the steadystate solutions of the CSA model (and thus the MRSAA model as well) is in stark contrast to earlier work for which the incision rate E is assumed to have the following dependence on slope S b and drainage area A ("slope-area" formulation, Howard and Kerby, 1983): where A denotes drainage area, n and m are specified exponents, and K is a constant assumed to decrease with increasing rock hardness. In order to compare the steady-state predictions of the slope-area relation in Eq. (37) for constant υ with CSA, drainage area A must be taken to be a constant value A o so as to correspond to the HSR configuration used here. The steady-state slope S bss corresponding to a balance between incision and rock uplift is found from Eq. (37) to be The issue as to the values of m and n has been considered by many researchers, including Whipple and Tucker (2002) and Lague (2014).
In their Table 1, Whipple and Tucker (2002) quote a range of values of n, but their most quoted value is n = 2. We compare the results for CSA for S bss with the predictions from Eq. (38) with n = 2 by normalizing against a reference value S bref that corresponds to a reference rock uplift rate υ ref of 0.2 mm yr −1 . Equation (38) yields In Fig. 11, Eq. (39) is compared against the CSA predictions of Figs. 9b and 10b (high and low feed rate, respectively) for both constant and variable β ss . In order to keep the plot within a realistic range, only values of υ between 0.2 and 10 mm yr −1 (the upper limit corresponding to Dadson et al., 2003) have been used in the CSA results. The remarkable insensitivity of the CSA predictions for steady-state slope S bss on rock uplift rate is readily apparent from the figure.
One more difference between the CSA and slope-area formulations is worth noting. If the slope-area relation is installed into Eq. (6) in place of CSA, it is readily shown that bedrock slope gradually relaxes to zero in the absence of rock uplift. CSA does not obey the same behavior under the constraint of constant sediment feed rate: Figs. 9b and 10b indicate that bedrock slope converges to a constant, nonzero value as rock uplift declines to zero. This is not necessarily a shortcoming of CSA; the sediment feed rate can be expected to decline as relief declines.   Figure 11. Normalized steady-state bedrock slope S bss /S bref versus normalized rock uplift rate υ as predicted by the CSA model for a low feed rate (orange-brown curves) and a high feed rate (black curve and dashed red curve), and for constant and variable abrasion coefficient. The results are the same for the MRSAA model. Also shown is the prediction of a model for which the incision rate is specified in terms of bedrock slope and upstream drainage area (green curve). Note that the predictions for steady-state bedrock slope of the CSA model are insensitive to the rock uplift rate over a wide range.

Boundary conditions and parameters for numerical solutions of the MRSAA model
Having conducted a fairly thorough analysis of the steady state common to the CSA and MRSAA models, it is now appropriate to move on to examples of behavior that can be captured by the MRSAA model, but are not captured by models that assume a relation for cover based on the ratio of sediment supply to capacity transport rate, i.e., Eq. (2). Before doing so, however, it is necessary to delineate the boundary conditions and other assumptions used in the MRSAA model. Let L denote the length of the reach. Equation (27a) indicates the formulation for bedrock incision is first order in x and so requires only one boundary condition. The example considered here is that of a downstream bedrock elevation, i.e., base level, set to zero: According to Eq. (18), or alternatively Eq. (28a), the alluvial formulation is second order in x and thus requires two boundary conditions. The following boundary condition applies at the upstream end of the reach, where q af (t) denotes a feed rate that may vary in time, At the downstream end, a free boundary condition is applied for η a /L mr < 1, and a fixed boundary condition is applied for η a /L mr ≥ 1 as follows: Here, Eq. (42a) specifies a free boundary in the case of partial alluviation, thus allowing below-capacity sediment waves to exit the reach. Equation (42b), on the other hand, fixes the maximum downstream elevation at η = η a = L mr . In order to illustrate the essential features of the new formulation of the MRSAA model for the morphodynamics of mixed bedrock-alluvial rivers, it is useful to consider the most simplified case that illustrates its expanded capabilities compared to the CSA model. Here we implement the HSR simplification. In addition, based on the results of the previous section, we approximate β ss as a prescribed constant. Finally, we assume that the clasts of the abrading bedload are sufficiently hard compared to the bedrock so that grain size D can be approximated as a constant. These constraints are easily relaxed.
In the numerical solution of the differential Eqs. (6b) and (18), spatial derivatives have been computed using an upwinding scheme for short timescales (so as to capture downstream-migrating alluvial waves) and a downwinding scheme for long timescales (so as to capture upstreammigrating incisional waves). Time derivatives have been computed using the Euler step method.

Sediment waves over a fixed bed: stripping and emplacement of alluvial layer and advection-diffusion of a sediment pulse
Three numerical solutions of the MRSAA model are studied here: (a) stripping of an alluvial cover to bare bed, (b) emplacement of an alluvial cover over a bare bed and (c) advection-diffusion of an alluvial pulse over a bare bed. Reach length L is 20 km. As the time for alluvial response is short compared to incisional response, β ss and υ are set equal to zero for these calculations. In addition, flood intermittency I is set to unity so as to illustrate the migration from the feed point to the end of the reach under the condition of continuous flow. The macro-roughness L mr is set to 1 m based on observation of the Shimanto River near Tokawa, Japan. The values for Cz, Q, B, D and λ are the same as in Sect. 5, i.e., Cz = 10, Q = 300 m 3 s −1 , B = 100 m, D = 20 mm and λ = 0.35. Bedrock slope S b , which is constant due to the absence of abrasion, is set to 0.004. The above numbers combined with Eqs (5c) (using the constants of the formulation of Wong and Parker, 2006a), (32a) and (32b) yield the following values: depth H = 1.32 m, Froude number Fr = 0.63, Shields number τ * = 0.016 and capacity bedload transport rate q ac = 0.0017 m 2 s −1 .
None of these three cases can be treated using models that assume a relation for cover based on the ratio of sediment supply to capacity transport rate, i.e., Eq. (2). They thus illustrate capabilities unique to MRSAA.

Alluvial stripping
The case of stripping of an initial alluvial layer to bare bedrock is considered here. In this simulation, the bedload feed rate q af = 0 and the initial thickness of alluvial cover η a is set to 0.8 m, i.e., 80 % of the macro-roughness length L mr . To drive stripping of the alluvial layer, the feed rate is set equal to zero. Figure 12a shows how the alluvial cover is progressively stripped off from upstream to downstream as a wave of alluvial rarification migrates downstream. The alluvial layer is completely removed (except for residual sediment in deep pockets, as specified by Eq. 21) after a little more than 0.12 years.
Of interest in Fig. 12a is the fact that the wave of stripping maintains constant form in spite of the diffusive term in Eq. (28a) which should cause the wave to spread. The reason the wave does not spread is the nonlinearity of the wave speed c a in Eq. (28b): since p enters into the denominator on the right-hand side of the equation, wave speed is seen to increase as p decreases, and thus η a decreases. As a result, the lower portion of the wave tends to migrate faster than the higher portion, sharpening the wave and opposing diffusion.

Emplacement of an alluvial layer over an initially bare bed
In this simulation, the initial thickness of alluvium η a is set to zero and the sediment feed rate is set to 0.0013 m 2 s −1 , i.e., 80 % of the capacity value. The result of the calculation is shown in Fig. 12b. Here nonlinear advection and diffusion act in concert to cause the wave of alluviation to spread. The steady-state thickness of alluvium is 0.83 m; by 0.1 years it has been emplaced only down to about 5 km from the source. This steady-state condition, and only this condition, corresponds to a convergence of results from MRSAA and CSA.

Propagation of a pulse of alluvium over an initially bare bed
In this example the initial bed is bare of sediment. The sediment feed rate is set equal to 0.0012 m 2 s −1 , i.e., 70 % of the capacity value for 0.05 years from the start of the run, and then dropped to zero for the rest of the run. Figure 12c shows the propagation of a damped alluvial pulse through the reach, with complete evacuation of the pulse in a little more than 0.15 years. Nonlinear advection acts against diffusion to suppress the spreading of the upstream side of the pulse, but advection acts together with diffusion to drive spreading of the downstream side of the pulse.

Comparison of evolution to uplift-driven steady state for the CSA and MRSAA models
Here we consider three cases of channel profile evolution to steady state that include both rock uplift and incision. In the first case, the initial bedrock slope is set to a value below the steady-state value, and the sediment feed rate is set to a value that is well above the steady-state value for the initial bedrock slope, causing early-stage massive alluviation. The configuration for the second case is a simplified version of a graben with a horst upstream and a horst downstream. The configuration for the third case is such that there is an alluviated river mouth downstream and a bedrock-alluvial transition upstream. In all cases, MRSAA predicts evolution that cannot be predicted by models that assume a relation for cover based on the ratio of sediment supply to capacity transport rate, i.e., Eq. (2).

Evolution of bedrock profile with early-stage massive alluviation
Here we set Q, B, Cz, D and λ to the same values as Sect. 6. The reach length L is 20 km, the flood intermittency I is set to 0.05, macro-roughness L mr is set to 1 m, initial alluvial thickness η a | t=0 = 0.5 m, downstream bed elevation η b | x=L = 0 and the abrasion coefficient β ss is 0.05 km −1 . The initial bed slope is 0.004. The feed rate is set to twice the capacity rate for this slope, i.e., q af = 0.0033 m 2 s −1 . The uplift rate is set to the very large value of 5 mm yr −1 . It should be noted, however, that as shown in Fig. 10b, the steady-state bedrock slope for this feed rate is independent of the uplift rate for υ ≤ 5 m yr −1 . This is because the steady-state value of ϕ is 0.019, i.e., ϕ 1.
The results for the CSA model are shown in Fig. 13a. The bed slope evolves from the initial value of 0.004 to a final steady-state value of 0.0068. Evolution is achieved solely by means of an upstream-migrating knickpoint. Only the first 4000 years of evolution are shown in the figure. Figure 13b shows the results of the first 400 years of the calculation with MRSAA. By 100 years, the bed is completely alluviated, and by 400 years, the thickness of the alluvial layer at the upstream end of the reach is 52 m. This massive alluviation is, unsurprisingly, not predicted by CSA, which was designed to treat incision only. Figure 13c shows the results of the first 4000 years of evolution. The upstreammigrating knickpoint takes the same form as CSA, but it is nearly completely hidden by the alluvial layer. The knickpoint gradually migrates upstream, driving the completely alluviated layer out of the domain, but this process is not complete by 4000 years. A comparison of Fig. 13a and c show that a knickpoint that is exposed in CSA is hidden in MRSAA. Models that assume a relation for cover based on the ratio of sediment supply to capacity transport rate cannot predict the presence of a hidden knickpoint.

Evolution of horst-graben configuration
In this example, Cz, Q, B, D, λ, I , β ss , L mr , L, η a | t=0 and η b | x=L are set to the values used in Sect. 7.1. The sediment feed rate q af = 0.00083 m 2 s −1 , and the initial bedrock slope S b is set to the steady-state value for a rock uplift rate of 1 mm yr −1 , i.e., 0.0027. The model is then run for a rock uplift rate of 1 mm yr −1 for the domains 0 ≤ x ≤ 8 km and 12 km ≤ x ≤ 20 km and a rock subsidence rate of 1 mm yr −1 for the domain 8 km < x < 12 km. This configuration corresponds to a simplified 1-D configuration of a graben bounded by two horsts, one upstream and one downstream.
This case cannot be implemented in models that assume a relation for cover based on the ratio of sediment supply to capacity transport rate, i.e., Eq. (2). This is because such models are not designed to handle the case of alluvial fill of accommodation space created by subsidence. The results for MRSAA are shown in Fig. 14 domains evolve to a steady state in terms of both bedrock elevation and alluvial cover. The bedrock elevation of the subsiding domain never reaches steady state because it is completely alluviated. The profile at the top of the alluvium in this domain has indeed reached steady state by 15 000 years, with a bed slope that deviates only modestly from the steadystate bedrock slope driven by a uniform υ = 1 mm yr −1 .  Figure 14. Evolution predicted by the MRSAA model for localized subsidence at a narrow graben superimposed on broader uplift. Note the bedrock-alluvial and alluvial-bedrock transitions at the margins of the graben. By 15 kyr, the bed top has reached steady state, even though the bedrock surface in the graben continues to subside. The regional rock uplift rate and graben subsidence rate are assumed constant for simplicity. Here η denotes elevation at the bed top, η b denotes elevation of the bottom of bedrock relief, and x denotes streamwise distance.

Evolution of river profile with alluviated zone at river mouth
In this example Cz, Q, B, D, λ, β ss , L mr , L and η a | t=0 are again set to the values chosen in Sect. 7.1. The bedload feed rate is 0.00083 m 2 s −1 ; the steady-state bedrock slope S b associated with this feed rate is 0.0026 for υ < 5 mm yr −1 (Fig. 10b). The initial bedrock slope is set, however, to the higher value of 0.004. The rock uplift rate υ for this case is set to zero, for which the steady-state slope is again 0.0026. The result of CSA for this case, with base level η b | x=L pinned at zero elevation, is shown in Fig. 15. As in the case of Sect. 7.1, the bedrock slope evolves from the initial value of 0.004 to the steady-state value 0.0026 by means of an upstream-migrating knickpoint. Only 4000 years of evolution are shown in the figure, by which time the knickpoint is 4.8 km from the feed point.
MRSAA is implemented with somewhat different initial and downstream boundary conditions in order to model the case of a bed that remains alluviated at the downstream end. This condition thus corresponds to an alluviated river mouth. The initial bedrock slope is again 0.004, and the downstream bedrock elevation η b | x=L is again 0 m. The downstream alluvial elevation η a | x=L , however, is held at 10 m, so that the downstream end is completely alluviated. The initial slope S for the top of the bed is 0.0021, a value chosen so that the bed elevation equals the bedrock elevation at the upstream end. Results of the MRSAA simulation are shown in Fig. 16. Figure 16a-c show the early-stage evolution, i.e., at t = 0, 10 and 100 years. Over this period, a bedrock-alluvial transition (from mixed bedrock-alluvial to purely alluvial) migrates downstream from the feed point to x = 13.6 km, i.e., 6.4 km upstream of the terminus. Bedrock incision is negligible over this period. Figure 16d-f show the bedrock and top bed profiles for 1000, 2000 and 4000 years. Over this period, the bedrock-alluvial transition migrates upstream. As it does so, the bedrock slope downstream of x = 13.6 km remains alluviated and does not change. The bedrock slope between the transition and x = 13.6 km evolves to the steady-state value of the case in Sect. 7.2, and the top bed slope downstream of the transition evolves to the same slope as the steadystate bedrock slope (because with υ = 0, ϕ vanishes). The figures show that the upstream-migrating bedrock knickpoint is located at the bedrock-alluvial transition. By 4000 years, the transition has migrated out of the domain and the bed is completely alluviated. The thickness of the alluvial cover upstream of x = 13.6 km is, however, only 1.05 m, i.e., only slightly larger than the macro-roughness height of 1 m. This means that although the reach is everywhere alluvial at 4000 years, the bedrock is only barely covered.

Discussion
The MRSAA model is a direct descendant of the model of Sklar and Dietrich (2004) in terms of the formulation for bedrock incision, and the model of Struiksma (1999) in terms of the formulation of the conservation alluvium over a partly covered bedrock surface. In terms of its capabilities, however, it shares much in common with the previous work of Lague (2010), and in particular with his SSTRIM model. These include (a) the melding of incision and allu- viation into a single model, (b) the inclusion of a cover relation that is based on geometric bed structure, and (c) the ability to track simultaneously the spatiotemporal variation in both incision rate and alluvial cover. Priority should accrue to Lague (2010) in regard to these features. The present model has the following advantages: (a) the Exner equation of sediment conservation is specifically based on a formulation of the statistics of partial and complete cover over a rough bedrock surface; (b) the formulation yields a specific relation for alluvial wave velocity as a function of cover, ranging to a maximum value for minimum cover to 0 for complete alluviation; and (c) it allows for explicit description of the nonlinear advective-diffusive physics of the problem in terms of an alluvial diffusivity and two wave celerities, one directed upstream and associated with bedrock incision, and one directed downstream and associated with alluviation. The form of the MRSAA model presented here has been simplified as much as possible, i.e., to treat a HSR (highly simplified reach) with constant grain size D. This has been done to allow for a precise and complete characterization of the behavior of the governing equations. It can relatively easily be extended to: (a) abrasion of the clasts that abrade the bed, so abrasional downstream fining is captured (Parker, 1991); (b) size mixtures of sediment (Wilcock and Crowe, 2003); (c) multiple sediment sources (Lague, 2010;Yanites et al., 2010); (d) channels with width variation downstream (Lague, 2010); (e) discharge varying according to a flow duration curve (Sklar and Dietrich, 2006;Lague, 2010), or fully unsteady flow ; and (f) cyclically varying hydrographs (Wong and Parker, 2006b) or "sedimentographs", the latter corresponding to events for which the sediment supply rate first increases, and then decreases cyclically . In addition, the model can and should be extended to include the stochasticity emphasized by Lague (2010). Sections 6 and 7 illustrate features captured by MRSAA but not by models that assume a cover relation based on the ratio of sediment supply to capacity transport rate, i.e., Eq. (2).
The MRSAA model presented here is applied to several 1-D cases with spatiotemporal variation. The model can easily be generalized to 2-D simply by expressing Eq. (18) in 2-D form. In any such implementation, however, the effect of 2-D connectivity between deep holes should be considered in the relation to the cover factor.
The MRSAA model in the form presented here has a weakness in that the flow resistance coefficient Cz is a prescribed constant. The recent models of bedrock incision of Inoue et al. (2014) and Johnson (2014) provide a much more detailed description of flow resistance. In addition to characterizing macro-roughness, their models use two micro-roughnesses, one characterizing the hydraulic roughness of the alluvium and the other characterizing the hydraulic roughness of the bedrock surface. Their models can thus discriminate between (a) "clast-smooth" beds, for which bedrock roughness is lower than clast roughness and (b) "clast-rough" beds, for which bedrock roughness is greater than clast roughness. This characterization allows for two innovative features: (a) both bed resistance and fractional cover become dependent on the ratio of bedrock micro-roughness to alluvial micro-roughness and (b) incision can result from throughput sediment passing over a purely bedrock surface with no alluvial deposit. The models of Inoue et al. (2014) and Johnson (2014) use modified forms of the capacity-based form for cover of Eq. (2) in order to capture these phenomena. Their models are thus unable to capture the coevolution of bedrock-alluvial and purely alluvial processes of MRSAA. Amalgamation of their models and the one presented here, however, appears to be feasible and is an attractive future goal.
Because MRSAA tracks the spatiotemporal variation in both bedload transport and alluvial thickness, it is applicable to the study of the incisional response of a river subject to temporally varying sediment supply. It thus has the potential to capture the response of an alluvial-bedrock river to massive impulsive sediment inputs associated with landslides or debris flows. A preliminary example of such an extension is given in Zhang et al. (2013). When extended to multiple sediment sources, it can encompass both the short-and long-term responses of a bedrock-alluvial river to intermittent massive sediment supply due to landslides and debris flows. As such, it has the potential to be integrated into a framework for managing sediment disturbance in mountain rivers systems such as those affected by the 2008 Wenchuan earthquake in Sichuan, China. Over 200 landslide dams formed during that event (Xu et al., 2009;Fu et al., 2011). A similar potential application is the case of drastic sediment supply to, and evacuation from, rivers in Taiwan due to typhoon-induced or earthquake-induced landsliding (e.g., Yanites et al., 2010).

Conclusions
We present a 1-D model of alluvial transport and bedrock erosion in a river channel whose bed may be purely alluvial, or mixed bedrock-alluvial, or may transition freely between the two morphologies. Our model, which we call the Macro-Roughness-based Saltation-Abrasion-Alluviation (MRSAA) model, specifically tracks not only large-scale bedrock morphodynamics but also the morphodynamics of the alluvium over it. The key results are as follows: 1. The transport of alluvium over a bedrock surface cannot in general be described simply by a supply rate that instantaneously affects the entire river reach downstream as it is varied in time. Here we track the alluvium in terms of a spatiotemporally varying alluvial thickness.
2. The area fraction of cover p enters into both incisional and alluvial evolution. The alluvial part allows for the downstream propagation and diffusion of sediment waves, so that at any given time the alluvial bed can be above or below the top of the bedrock. The model thus allows for spatiotemporal transitions between complete cover, under which no incision occurs; partial cover, for which incision may occur; and no cover, for which no incision occurs.
3. The MRSAA model captures three processes: downstream alluvial advection at a fast timescale, alluvial diffusion, and upstream incisional advection at a slow timescale. Only the third of these processes is captured by models that assume a relation for cover based on the ratio of sediment supply to capacity transport rate rather than a measure of the thickness of alluvial cover itself. The CSA model can be thought of as a 0-D model that applies locally. The MRSAA model lends itself more directly to application to long 1-D reaches because it embeds the elements necessary to route sediment down the reach.
4. The MRSAA model reduces to the CSA model under the conditions of steady-state incision in balance with rock uplift and below-capacity cover. The steady-state bedrock slope predicted by both models is insensitive to the rock uplift rate over a wide range of conditions. This insensitivity is in marked contrast to the commonly used incision model in which the incision rate is a power function of bedrock slope and drainage area upstream. The two models can differ substantially under transient conditions, particularly under those that include migrating transitions between the bedrock-alluvial and purely alluvial state.
5. In the MRSAA model, inclusion of alluvial advection and diffusion lead to the following phenomena: (a) a wave-like stripping of antecedent alluvium over a bedrock surface in response to cessation of sediment supply, (b) advection-diffusional emplacement of a sediment cover over initially bare bedrock and (c) the propagation and deformation of a sediment pulse over a bedrock surface.
6. In the case of transient imbalance between rock uplift and incision with a massive increase in sediment feed, MRSAA captures an upstream-migrating transition between a purely alluvial reach upstream and a bedrock-alluvial reach downstream (here abbreviated as a alluvial-bedrock transition). The bedrock profile shows an upstream-migrating knickpoint, but this knickpoint is hidden under alluvium. Models that assume a relation for cover based on the ratio of sediment supply to capacity transport rate, i.e., Eq.
(2), capture only the knickpoint, which is completely exposed, and thus miss the thick alluvial cover predicted by MRSAA.
7. MRSAA captures the mixed incisional-alluvial evolution for the case of a simplified 1-D subsiding graben bounded by two uplifting horsts. It captures alluvial filling of the graben, and thus converges to a steady-state top-bed profile with a bedrock-alluvial transition at the upstream end of the graben and an alluvial-bedrock transition at the downstream end.
8. In the case studied here of an uplifting bedrock profile with an alluviated bed at the downstream end modeling a river mouth, MRSAA predicts an upstream-migrating bedrock-alluvial transition at which the bedrock undergoes a sharp transition from a higher to a lower slope. MRSAA further predicts a bedrock long profile under the alluvium that has the same slope as the top bed. It also predicts that the cover is thin, so that the purely alluvial reach is only barely so. The steady state for this case is purely alluvial.
9. The new MRSAA model provides an entry point for the study of how bedrock-alluvial rivers respond to occasional large, impulsive supplies of sediment from landslides and debris flows. It thus can provide a tool for forecasting river-sedimentation disasters associated with such events. An example application would be treatment of the aftereffects of the 2008 Wenchuan earthquake, which triggered massive alluviation and the formation of over 200 landslide dams.

Appendix A: Interpretation of the abrasion coefficient β
Consider a clast or grain of size D and volume V g ∼ D 3 causing abrasion over an exposed bedrock surface. The bedload transport rate q a is given as where E g denotes the volume rate per unit time per unit area at which clasts are ejected from the bed into saltation, and L g denotes the saltation length. Each clast ejected into saltation collides with the bed a distance L g later; therefore, the number of clasts that collide with the bedrock (rather than with other bed particles) per unit time per unit area is found using Eq. (A1), The volume lost from the striking clast per strike is defined as β * g V g , and the volume lost from the stricken bedrock per strike is similarly defined as β * V g . The parameters β * g and β * could be expected to be approximately equal if the striking grain is the same rock type as the bedrock.
The rate at which a grain strikes the bed per unit distance moved is 1/L g ; hence the rate at which grain volume decreases downstream is given by Using V g ∼ D 3 , Eq. (A3) reduces to Equation (A4a) is the differential form of Sternberg's law; α d is a diminution coefficient with units L −1 . The exponential form of Eq.
(3) corresponds to the case of spatially constant α d .
The incision rate of the bedrock E is the number of grains that collide with bedrock per unit area per unit time (Eq.A2) multiplied by the volume lost per strike β * V g , which gives the relation The above relation is identical to Eq. (1b).