Hillslopes in humid regions are typically convex to concave in profile and have a relatively thick, continuous regolith cover. Conversely, hillslopes in arid regions are typically cliff-dominated and have a relatively thin, discontinuous regolith cover. The difference between these two end-member slope forms is classically attributed to climate, but climate, tectonics, and lithology all play a role. In this paper, we describe a mathematical model for hillslope gradient and regolith thickness using basic climatic and tectonic input data for a given rock type. The model first solves for the regolith thickness on a planar slope segment in topographic steady state using the soil production function and a prescribed uplift/incision rate. The climatic and lithologic controls on soil production rates are quantified using an empirical energy-based model for the physical weathering of bedrock. The slope gradient is then computed by balancing uplift/incision rates with sediment fluxes calculated using a nonlinear depth- and slope-dependent sediment transport model. The model quantifies the ways in which, as aridity and uplift/incision rates increase, regolith thicknesses decrease and hillslope gradients increase nonlinearly until a threshold condition is reached, beyond which bare, cliff-dominated slopes form. The model can also be used to estimate long-term erosion rates and soil residence times using basic input data for climate and regolith thickness. Model predictions for erosion rates closely match cosmogenically derived erosion rates in granitic landscapes. This approach provides a better quantitative understanding of the climatic and tectonic controls on slope form, and it provides a simple, widely applicable method for estimating long-term erosion rates and the thickness of regolith cover on hillslopes.


It is a “textbook” observation in geomorphology that hillslopes in arid regions tend to be cliff-dominated and have a thin, discontinuous regolith cover, while hillslopes in more humid regions tend to be convex to concave in profile with a relatively thick, continuous regolith cover (Ritter et al., 2002) (Fig. 1). These differences in slope form are also associated with distinctly different modes of evolution: cliff-dominated, bare-regolith slopes tend to evolve by lateral slope retreat, thus maintaining steep slopes over time, even in the absence of uplift, while convex to concave slopes with continuous regolith cover tend to evolve by slope replacement, evolving to gentler slopes over time (Fig. 1). While the formation of these two end-member slope types is classically attributed to climate, uplift/incision rates and lithology must also play a role. In arid climates where weathering rates are very low, for example, regolith production can, nevertheless, keep pace with uplift rates if those rates are sufficiently low. Rock type also plays a role in controlling regolith cover and hillslope gradient, because rock types differ in their degree of weatherability. The geologic conditions required for the formation of bare, cliff-dominated slopes must, therefore, be a function of climate, tectonics, and lithology. In addition, the precise mechanisms of climate control on slope form (i.e., what are the respective roles of temperature and precipitation?) are not well understood.

Recently, three major advances have been made that bear on the link between slope form and climate and tectonics. First, Heimsath et al. (1997, 1999) used in situ cosmogenic isotopes to constrain the soil production function to be an exponential function of regolith thickness. Heimsath et al.'s work directly links the rate of physical weathering of bedrock with regolith thickness. Second, Roering and coworkers quantified the nonlinear dependence of sediment flux on slope gradient using experimental, field, and modeling studies (e.g., Roering et al., 1999; 2001; Roering, 2004). Third, several recent studies have documented the dependence of sediment flux on regolith thickness or depth (e.g., Gabet, 2000; Heimsath et al., 2005). By combining the nonlinear slope-dependent transport model of Roering et al. with the depth-dependent model of Gabet, Heimsath et al., and others, we get the following equation for unit sediment flux qs:  
where κ is a transport coefficient, hn is the regolith thickness normal to the surface, S is the slope gradient, and Sc is the tangent of the angle of stability. Equation 1 is not appropriate for very thick soils (i.e., more than several meters) because the sediment flux cannot continue to increase indefinitely as soil thickness increases (unless the hillslope is near the angle of stability and capable of transporting the entire soil profile by mass movement). More complex, depth-dependent models are available that capture this depth-saturation effect (e.g., Roering, 2008). In a comparison of several end-member sediment transport models proposed in the literature, Roering (2008) found that a depth- and slope-dependent transport model similar to that of Equation 1 was the most accurate model for quantifying sediment transport on the steep upland slopes of the Oregon Coast Range.

The exponential soil production function implies that bedrock lowering is a maximum for bare bedrock slopes and decreases exponentially with increasing regolith thickness. This exponential relationship has been inferred from in situ cosmogenic isotope measurements in upland soil profiles (Heimsath et al., 1997, 1999). Conceptually, this relationship reflects the fact that regolith cover acts as a buffer for the underlying bedrock, protecting it from the diurnal temperature changes and infiltrating runoff that drive subsurface weathering. The exponential soil production function may not capture the full complexity of regolith production, however. As regolith thickness decreases below a critical value in arid and semiarid climates, the landscape may be unable to store enough water to promote weathering or support plant life. Therefore, in some environments, weathering rates may actually increase as regolith thickness increases (for relatively thin regolith cover), leading to a “humped” production function (Fig. 1D) (Ahnert, 1977; Cox, 1980; Dietrich et al., 1995; Anderson and Humphrey, 1989; Furbish and Fagherazzi, 2001; Anderson, 2002; Furbish, 2003; Minasny and McBratney, 1999, 2006). Recent cosmogenic radionuclide data from granitic landscapes in Australia provide at least preliminary support for a humped production function (Heimsath et al., 2006). Also, depending on whether tors are included or separated from analyses of soil production data sets, earlier cosmogenic analyses may also be consistent with a humped production function. Tors are small bedrock hills characterized by thin regolith cover and relatively low weathering rates. Because tors are associated with locally more resistant bedrock, Heimsath et al. (1997, 1999) excluded tors from their analysis of the soil production function. Strudley et al. (2006), however, argued that tor formation, although triggered by locally resistant bedrock, also involves a feedback process between slope gradient, regolith cover, and erosion that they associated with a humped production function. In this paper, we consider both the exponential and humped production functions in our analyses, comparing and contrasting the model predictions for these two possible end-member forms of the soil production function.


Hillslopes in upland landscapes are composed of a system of two interacting surfaces: the topographic surface, with elevations given by z(x,y), and the underlying bedrock surface or weathering front, given by b(x,y). The difference between these two surfaces is the soil or regolith thickness, h(x,y). The topographic and weathering-front surfaces are strongly coupled because the shape of the topography controls erosion and deposition, which, in turn, changes the values of h(x,y) (e.g., Furbish and Fagherazzi, 2001). The values of h(x,y), in turn, control weathering rates. In the following analysis, we will assume topographic steady state at the hillslope scale, i.e., that bedrock uplift and soil erosion rates are equal after the effect of the relative densities of rock and soil are taken into account. The appropriateness of this assumption and the robustness of the model predictions with respect to this assumption will be addressed in the discussion section. If topographic steady state is assumed, the rate of bedrock lowering must equal the rate of bedrock uplift:  
assuming the exponential form of the soil production function, and  
assuming the humped form of the soil production function plotted in Figure 1, where P0 is the maximum bedrock lowering rate, h0 is a characteristic regolith depth, and U is the uplift/incision rate. Equations 2 and 3 assume that regolith production is controlled by the regolith thickness measured vertically. It is more reasonable to assume that regolith production is controlled by the regolith thickness measured normal to the surface, since it is that latter thickness that buffers the bedrock from the elements (Heimsath et al., 2001). This “normal-directed” soil production function has the effect of introducing a cos θ term (where θ is the hillslope angle) into the analysis. We will consider that effect later in this section. The humped production function represented by Equation 3 is a linear function of h for small h that rises to a maximum value of approximately one-half the maximum rate of the exponential form and then matches the exponential form for larger values of h. This specific form of the humped production function is not unique, and other forms that increase from a low value at h = 0 to an exponential form at large values of h are possible. We refer to U as the “uplift/incision” rate because it is the rate of channel incision into bedrock at the base of the slope, not the regional uplift rate, that controls hillslope response to uplift. Local incision rates may be higher or lower than regional uplift rates, depending on how the channel network responds to regional uplift.
Equations 2 and 3 can be solved for the regolith thickness h along a planar slope segment in topographic steady state given values for P0, U, and h0. For the exponential production function, regolith thickness is obtained by solving Equation 2 for h to obtain  
If P0/U < 1, the model predicts that a bare bedrock slope will form. For the humped production function, solving Equation 3 for h gives  
With the humped production function, two solutions for h are possible. The solution that applies depends on whether uplift rates are increasing or decreasing over time. We will return to this point later in this section. Given a solution for regolith thickness (h), the soil residence time can be computed as  
where E is the rate of soil erosion, and ρb and ρs are the densities of the bedrock and soil, respectively.
Climate and rock type enter into the model through their control on the maximum rate of bedrock lowering, P0 (or P0/2 in the humped model). Recently, Rasmussen et al. (2005) and Rasmussen and Tabor (2007) proposed a new variable, called “effective energy and mass transfer” (EEMT), for quantifying the role of climate in pedogenic processes. Unlike previous attempts to empirically relate rates of pedogenesis to climatic indices such as mean annual temperature and precipitation, the EEMT variable is explicitly based on a conceptual model of energy input and flow through the soil profile, including thermal (temperature) and material (water and biomass) forms of energy. These authors demonstrated the effectiveness of the EEMT framework by documenting significant correlations between EEMT and soil order, regolith thickness, and a variety of pedogenic indices using data from study sites over a broad range of climates and rock types. Plotting these data versus mean annual temperature (MAT) and mean annual precipitation (MAP), in contrast, yielded no significant correlations. As such, the EEMT variable appears to provide an effective quantitative measure of the controlling influence of climate and lithology on rates of regolith production. EEMT may be computed using an empirical function of MAT and MAP (Rasmussen and Tabor, 2007):  
Biotic controls on EEMT are implicitly included in Equation 7 because they are empirical functions of MAP and MAT. Rasmussen and Tabor (2007) found that regolith thickness on stable (low gradient) slopes increased exponentially with EEMT, suggesting that rates of regolith production also increase exponentially with EEMT. This hypothesis suggests a relationship of the form  
where a (units of m/k.y.) and b (units of m2/[kJ yr]) are empirical coefficients. In order to constrain the parameters a and b for granitic landscapes, we used published data from Riebe et al. (2004) (Figs. 2A and 2B). Physical erosion rates measured by Riebe (E) were used to estimate the rate of peak regolith production using the equation  
assuming the exponential production function, where hn is the soil thickness normal to the surface, and  
assuming the humped production function, with h0 = 0.5 m. The resulting estimates for P0 are plotted in Figures 2A and 2B as a function of EEMT on logarithmic scales, assuming ρbs = 1.35. The data points plotted in Figures 2A and 2B include all of the data published in Riebe et al. (2004) for which regolith thicknesses were measured. The continuous curve plotted in Figures 2A and 2B corresponds to Equation 8 with bestfit coefficients a = 0.037 m/k.y. and b = 0.00003 m2/(kJ yr). The two data points plotted with open circles represent outliers in this and subsequent analyses of the data set from Riebe et al. (2004). Potential limitations of the model with respect to these outlying points will be discussed at the end of this section. The coefficients a and b vary with rock type, and the variation in these parameters represents the control lithology has on slope and soil thickness in the model. In this paper, we focus on model results for granitic landscapes only, because it is in these landscapes where regolith production rates are currently most well constrained. Rasmussen and Tabor (2007), however, provided preliminary values of a and b for basalt and andesite as well. Future studies that apply the model to those rock types, therefore, may use the preliminary values provided by Rasmussen and Tabor (2007).

Figures 3A and 3B plot model predictions for regolith thickness h given by Equations 3 and 4 as a function of EEMT and uplift/incision rate for both the exponential (Figs. 3A and 3B) and humped production functions (Figs. 3C and 3D), assuming a = 0.037 m/k.y., b = 0.00003 m2/(kJ yr), and = 0.5 m. For the exponential production function, regolith thickness h0 increases linearly with EEMT and decreases exponentially with uplift/incision rate until the threshold condition P0/U < 1 is met, beyond which, bare slopes are formed. For example, assuming an EEMT of 50,000 kJ/(m2 yr) (i.e., a moist, temperate climate) and an uplift/incision rate of 0.01 m/k.y., Figure 3A predicts a regolith thickness of 1.4 m. For the humped production function, the threshold condition for bare slopes to form requires lower values for U and/or higher values for P0 compared to the exponential case, i.e., P0/U < 2, meaning that bare slopes will form under a wider range of climatic and uplift/incision values compared to the exponential case. Also, in the humped production case, there is a “left” and a “right” branch of solutions corresponding to relatively thin and thick regolith, respectively. The “thin” solutions are plotted in Figure 3B, while the “thick” solutions are the same as those of the exponential production function plotted in Figure 3A for P0/U < 2 (otherwise h = 0). The solution that is applicable depends on whether the value of P0/U is increasing or decreasing with time. If a hillslope is initially in a state of tectonic quiescence (i.e., thick soils have had time to form) and the hillslope experiences steadily increasing uplift rates, regolith thickness will decrease from its initially high value along the right branch of the humped production function until uplift outpaces production and a bare slope forms. Alternatively, if a hillslope begins as a bare slope (e.g., formed during a period of high uplift rates) and uplift rates steadily decrease, the regolith thickness will increase along the left branch until equilibrium is achieved with the new uplift rate. Regolith thicknesses are much lower along the left branch compared to the right branch for similar uplift rates. As such, if a hillslope is governed by a humped production function, regolith thicknesses will be lower for uplift rates decreasing below the threshold condition for bare, cliff-dominated slopes to form compared to regolith thicknesses in conditions of increasing uplift rates. This behavior suggests a type of “hysteresis,” in which thin soils persist on a landscape subject to decreasing uplift rates, while thick soils persist on a landscape subject to increasing uplift rates. In this example, we described the hysteresis effect in terms of increasing or decreasing uplift rates, but climate change also exerts an influence through its effect on P0. Hence, climate change could thicken or thin the regolith cover hysteretically in a similar manner.

The slope gradient of a hillslope can be estimated using the solutions for regolith thickness plotted in Figures 3A and 3B together with the relationship  
Equation 11 applies the sediment continuity equation to the base of a slope of horizontal length L with a basal slope gradient S and basal regolith thickness normal to the surface hn, undergoing constant erosion at a rate E. Equation 11 is general and does not place any limitations on the variations in S and hn upslope. However, when Equation 11 is used together with input values for hn obtained by solving Equations 9 and 10 (which assume a planar slope segment in topographic steady state), then the predictions of Equation 11 are exactly applicable only to planar slopes and thus provide only approximate slope gradients for more general hillslope shapes. Equation 11 is a quadratic equation with solution  

Equation 12 is plotted in Figures 3C and 3D for the exponential and humped production functions, assuming the same parameters as in Figures 3A and 3B and κ = 0.1 m/k.y., L = 50 m, and Sc = 1. Figures 3C and 3D illustrate that slope gradients are inversely related to regolith thickness. For all climatic and tectonic conditions for which bare slopes are predicted (i.e., h = 0), the model predicts slope gradients that are equal to the threshold of stability (i.e., S/Sc = 1). However, although the maximum value for S/Sc is 1 in the model, actual maximum slope gradients may be much larger because topographic steady state is impossible to achieve above the threshold condition for bare slopes. As such, the slope gradient will continually increase in such cases as long as the uplift/incision rate outpaces soil production. Figure 3C illustrates that the model prediction that slope gradient increases with increasing uplift/incision rate and decreasing EEMT values until bare, cliff-dominated slopes form. If the humped production function is assumed, there are two solutions for slope gradient corresponding to the two solutions of regolith thickness plotted in Figures 3A and 3B. Figure 3D shows that, in cases where thin regolith forms on the left branch of the humped production function, the slope will compensate with very steep slopes over a wide range of climates and uplift rates. This result underlines the fact that bare, cliff-dominated slopes can persist under a wide range of conditions once they form in a given landscape (if the humped production function applies). This model behavior may help explain the persistence of cliff-dominated slopes even in areas not undergoing active tectonic uplift (e.g., Monument Valley and other cliff-dominated landscapes of the Colorado Plateau that are far enough from the Colorado River to have been unaffected by late Cenozoic incision). Qualitatively, the ways in which the slope gradients predicted by the model vary with climate and uplift rate do not depend on the specific values of κ, L, and Sc. Quantitatively, predicted slope values decrease as the ratio κSc/L increases. The threshold for cliff-dominated slopes to form, however, given by P0/U < 1 and P0/U < 2 for the exponential and humped production functions, respectively, does not depend on κ, L, or Sc.

Equations 2 and 3 relate regolith production to the regolith thickness in the vertical direction. As noted at the beginning of this section, it is more appropriate to relate these variables in a coordinate system normal to the surface (Heimsath et al., 2001). If we adopt this approach, then Equations 2 and 3 become  
assuming the exponential form of the soil production function, and  
assuming the humped form. If slope gradients are known, the analytic solutions to Equations 13 and 14 are  
Equations 11 and either 13 or 14 form a set of two coupled equations for two unknowns that cannot be solved analytically. The simplest numerical method for solving these equations is to use the analytic solutions for Equations 12 and either 15 or 16 to iteratively converge on the solution. In this method, Equation 15 or 16 is used to estimate h assuming S = 0. This estimate for h is then used in Equation 12 to calculate a value for S that is then used in Equation 15 or 16 to refine the value of h. This process can be repeated and converges to the correct solution in all cases. Figure 4 illustrates the numerical solutions for h and S obtained in this way for the exponential case. For S/Sc < 0.4, the analytic solutions predicted by Equations 3, 4, and 12 for the original, vertically directed production model differ from the numerical solutions for the revised, normally directed production model by less than 1%. For steeper slopes, however, the revised model predicts significantly thicker regolith for the same values of U and EEMT compared to the original model. Slope gradients predicted by the revised model are not significantly different from those of the original model as a function of uplift/incision rate and EEMT.

Thus far, we have assumed uplift rate to be a prescribed input to the model. The model framework can also be inverted to estimate erosion rates by solving Equations 13 and 14 using site-specific measurements for h and S and empirical data for P0 and h0. In contrast to the previous results in this paper, this approach does not require topographic steady state and hence can be expected to be more generally applicable. Figures 5A and 5B plot the erosion rate as a function of regolith thickness normal to the surface (hn) and EEMT for the exponential and humped production functions, respectively. The corresponding soil residence times, calculated using Equation 6, are plotted in Figures 5C and 5D. The surfaces plotted in Figures 5C and 5D dip down and away from the reader, so that the thickest soils have the largest residence times. As a test of the model predictions for erosion rate, Figures 2C and 2D illustrate the predicted erosion rates versus the measured rates using the data of Riebe et al. (2004), assuming exponential and humped production functions, respectively. The straight line in this figure represents an exact match between the predicted and measured values. Figure 2C shows that the model does an excellent job of predicting erosion rates using both the exponential and humped production functions, except for two sites (FR6 and FR8) on the low-relief, high-elevation plateau of the Sierra Nevada (shown as open circles), where predicted rates are several times higher than actual rates. Persistent snow cover is one reason why the predicted erosion rates are substantially higher than the measured erosion rates at sites FR6 and FR8. Persistent snow cover affects both the accuracy of cosmogenically derived erosion rates, and it serves to buffer the subsurface bedrock from diurnal temperature changes and, hence, may decrease erosion rates relative to model predictions. The empirical function for EEMT (i.e., Eq. 7) does not account for seasonality and timing of precipitation inputs because it is based on mean annual temperature and precipitation only. As such, higher MAP values drive higher weathering rates regardless of precipitation type, i.e., rain or snow. In areas of persistent snow cover, such as the low-relief, high-elevation zones of the Sierra Nevada, the empirical function in Equation 7 may be of limited accuracy because higher MAP values in such cases may equate to lower weathering rates resulting from the buffering effect provided by the persistent snow layer. The original formulation of EEMT (Rasmussen et al., 2005) is based on monthly data and mitigates this problem to an extent because precipitation that falls during below-freezing temperature conditions does not contribute to weathering.


The framework of this paper relies upon the assumption of topographic steady state at the hillslope scale. Given that topographic steady state is difficult to prove and may also be rarely achieved in nature, it is appropriate to ask how robust the model predictions are with respect to this assumption. Topographic steady state assumes that uplift and erosion are in precise balance, implying that elevation values do not change over time. If, however, uplift and erosion rates are in approximate balance (e.g., values of E and ρbU/ρs are within 20% of each other), then the model predictions can be expected to apply with comparable accuracy. Therefore, topographic steady state need not apply precisely in order for the model to provide useful predictions. For example, even in a landscape with no active tectonic uplift, isostatic adjustment generates rock uplift rates that are equal to ~80% of erosion rates. Therefore, an approximate balance between uplift and erosion rates can be expected to hold even in landscapes that are not subject to active tectonic uplift.

In the arid and semiarid regions of the southwestern United States, cliff formation is especially prevalent in the relatively undeformed sedimentary rocks of the Colorado Plateau. In order to understand that prevalence within the context of this paper, it is important to note that the erosion of bedded or banded rocks on cliff-dominated slopes is a special case because the erosion of one stratum may depend on the erodibility of strata exposed above it, because weaker strata cannot erode faster than more slowly eroding, resistant cap rock units above them. For this reason, cliff formation is more common in layered strata because the effective resistance to erosion of the entire exposed sequence is controlled by the most resistant units within that sequence. In the model of this paper, the role of rock type is quantified using lithologically dependent coefficients linking regolith production rate to climate. As more quantitative data become available on the values of these parameters for different rock types and as the model becomes more broadly applied, it will be important to remember that in bedded and banded rocks, the erodibility of one rock type can exert a rate-limiting effect on the erosion of rocks exposed below it.

The method used in this paper provides a quantitative tool for estimating the regolith cover and slope gradient based on readily available indices for climate (i.e., MAT and MAP) and tectonics (i.e., uplift/incision rate). Roering et al. (2007) took a broadly similar approach, utilizing the nonlinear slope–dependent transport model to derive a generally applicable analytic relationship between erosion rate and hillslope relief. This approach by Roering et al. provides a very useful theoretical basis for quantitatively linking slope gradient and erosion rates. Here, we extend the results of Roering et al. by including both transport-limited (regolith covered) and weathering-limited (bare) slopes (and the transition between them) and the effects of climate and lithology on erosion rates and slope form.


In this paper, we presented a new method for estimating regolith thickness, slope gradient, and erosion rates as a function of basic tectonic and climatic input data. The model provides a more complete understanding of the controlling factors that influence the formation of cliff-dominated versus convex to concave slopes. Model results highlight the dynamic complexity inherent in landscapes governed by a humped soil production function. Application of the model framework is currently limited by the relative paucity of data on regolith thickness and long-term erosion rates. The model framework would therefore greatly benefit from additional calibration and test data, particularly in humid climates, areas of persistent snow cover, and in study sites with nongranitic lithologies.