In this study we aim to evaluate potential controls on the relationship between drainage density and erosion rate, using both numerical modeling and analysis of real landscapes with high‐resolution topographic data. We develop a 1‐D analytical model using linear and nonlinear hillslope sediment flux laws, along with detachment‐limited fluvial incision models, to examine the effect of different parameters on the relationship between drainage density and erosion rate. We then use the Channel‐Hillslope Integrated Landscape Development (CHILD) model [ Tucker et al. , 2001 ] to test our analytical predictions, using both steady state and transient scenarios. We constrain channel network density using a recently developed technique for extracting channel networks that takes advantage of high‐resolution (1 m) light detection and ranging (lidar) data sets in order to test our model predictions on real landscapes. We compare drainage density to basin‐averaged erosion rates obtained from detrital cosmogenic radionuclide (CRN) analyses. In sites where CRN‐derived erosion rates are not available, we calculate the mean hilltop curvature ( C HT ) of each basin. Mean C HT has been demonstrated by Hurst et al. [ 2012 ] to vary linearly with erosion rate in high‐relief soil‐mantled landscapes, providing the basins are eroding at a uniform rate.

Several authors have predicted, based on numerical and analytical models, the relationship between D d and metrics describing valley morphology, including slope, erosion rate, and sediment transport process [ Kirkby , 1980 , 1993 ; Willgoose et al. , 1991 ; Tucker and Bras , 1998 ]. Here we build on previous studies that have attempted to examine the relationship between D d and erosion rate directly [ Montgomery and Dietrich , 1989 ; Oguchi , 1997 ; Talling and Sowter , 1999 ; DiBiase et al. , 2012 ]. These studies used varying methods for identifying the channel to hillslope transition, including slope‐area scaling relationships. However, such methods can lack precision due to noise when extracting slope from digital elevation models (DEMs). Recently developed methods of channel head identification allow the extraction of accurate drainage networks from high‐resolution DEMs [ Passalacqua et al. , 2010a , 2010b , 2012 ; Pelletier , 2013 ; Clubb et al. , 2014 ]. In addition, early studies were limited by the coarse resolution of the DEMs available at the time.

How valley spacing, and the associated landscape properties of hillslope length and drainage density (), change with erosion rates has been predicted to be sensitive to parameters in common fluvial incision models. Fluvial incision can be modeled using a detachment‐limited scenario in which the incision rateis a power law function of upstream drainage areaand channel slope[e.g.,]:whereis an erodibility coefficient (T), andandare constant exponents. Theandexponents in the stream power model have been shown to control the relationship between erosion rate and topographic gradient []. This relationship has important implications for how landscapes respond to changing tectonic forcing.] showed analytically that if the parameters are such that the fluvial incision model forecasts a linear relationship between erosion rate and slope, then river profiles will retain features that reflect changes in erosion rates (such as knickpoints). This is assumed in many studies and means that river profiles can be inverted to obtain uplift histories over millions of years, for example [], under the assumption that knickpoints are not structurally controlled. If the relationship is nonlinear, however, channels will be imperfect recorders of channel uplift history as the rate of knickpoint migration becomes a function of incision rate []. The value of these parameters in fluvial incision models also controls whether drainage density is sensitive to changing erosion rates [].also affects the transit time of water through the landscape and thus exerts a primary control on catchment flood response [].

One of the most distinctive features of soil‐mantled upland landscapes is the repeating patterns of ridges and valleys. The spacing of these ridges and valleys is fundamentally controlled by the competition between creep‐like sediment transport processes, which tend to smooth the landscape, and fluvial processes, which incise the landscape [ Tarboton et al. , 1992 ; Tucker and Bras , 1998 ; Perron et al. , 2012 ]. Perron et al. [ 2008 ] elegantly demonstrated that the spacing of valleys reflects the relative efficacy of advective (e.g., fluvial) and diffusive (e.g., hillslope) transport processes, both of which may be influenced by climate. Sweeney et al. [ 2015 ] used laboratory experiments to further demonstrate that the competition between hillslope and valley‐forming sediment transport processes controls the degree of landscape dissection. The erosion rate of the landscape also plays a major role in controlling the spacing of ridges and valleys [ Tucker and Bras , 1998 ; Perron et al. , 2008 ].

Our analytical model outlined above is the simplest case scenario, with a number of assumptions. For example, we assume a linear relationship between A and x t . Alternative model formulations are possible, such as that of a power law relationship, where [ Pelletier et al. , 2016 ]. However, using this alternative power law relationship predicts the same relationship between D d and E as that of our simpler scenario (see supporting information ). Furthermore, our model scenario neglects colluvial infilling of the valley head (equation 4 ). Pelletier et al. [ 2016 ] present a transport‐limited model for predicting drainage density, calibrated for the Walnut Gulch Experimental Watershed, Arizona, where they assume that the valley head occurs where the fluvial erosion rate is greater than the colluvial deposition rate by an amount equal to the net erosion rate E . In the model outlined above, we follow a similar approach to Tucker and Bras [ 1998 ] by assuming that the valley head occurs where the fluvial erosion rate is greater than colluvial erosion. Including colluvial deposition at the valley head in our analytical model may lead to decreased absolute values of drainage density. In order to test our simple 1‐D predictions we therefore carried out 2‐D numerical modeling, described in section 3.1 .

Analytical predictions of the relationship between drainage density and erosion rate for (a) linear hillslope sediment flux and (b) nonlinear hillslope sediment flux. We set parameters in equations 6 and 8 to the following:= 0.0088 myr= 1×10m yr= 0.5, and=1.25. The values of these parameters are the same as for the numerical modeling runs (Table 1 ). The relationship depends on the value ofin the stream power law: we predict a positive relationship for> 1, a negative relationship for< 1, and no relationship betweenandfor= 1.

We can equate the slopes above and below the valley head, given byand, respectively (equations 3 and 4 ), to obtain the mean length of overland flow:where. This reduces to(T) if we assume a rectangular flow strip of unit width (i.e.,= 1 m). The mean length of overland flow is approximately equal to half the reciprocal of]; and therefore, equation 5 can be converted towhere= (1 −)/(). The relationship betweenand erosion ratedepends on the slope exponent, shown in Figure 1 a. If> 1, the contributing area at the valley head will decrease with increasing erosion rate andwill therefore increase. However, if< 1 thenwill decrease with increasing erosion rate []. Performing a power law regression on a plot of mean length of overland flow against erosion rate allows the calculation of theexponent, assuming theratio is known, asis the gradient of the regression.

Theexponent may be constrained by examining the relationship betweenand erosion rate []. We can represent drainage density using the downslope distance from the hilltop to the valley head,, at which the slopes above and below the valley head are equal. The equilibrium slope for linear hillslope diffusion () can be expressed as [e.g.,whererepresents a diffusivity coefficient (L). We assume thatanddo not vary with erosion rate. In order to equate the channel slope,, given by equation 2 , with, we assume that the contributing areaat the valley head is given by a flow strip of length(L) and width(L), where

Choosing correct values of the exponents m and n is important in landscape evolution studies, because these values control the relationship between landscape steepness and erosion rates, as well as the competition between advective (fluvial) and diffusive (hillslope) processes. Although the m / n ratio has been reported for many landscapes [ Stock and Montgomery , 1999 ; Whipple and Tucker , 1999 ; Anthony and Granger , 2007 ; Perron and Royden , 2013 ; Mudd et al. , 2014 ], relatively few studies have reported individual values for the m and n exponents, as the erosion rate and K coefficient must be known. In particular, the slope exponent n is a critical parameter as it largely controls the timescale and magnitude of fluvial response to perturbations [ Whipple and Tucker , 1999 ]. Howard and Kerby [ 1983 ] studied a gradient in erosion rates across badlands in Virginia and suggested values of m = 0.45 and n = 0.7. Whittaker and Boulton [ 2012 ] examined knickpoint retreat rates above active faults in the Mediterranean and found that increasing channel steepening with fault throw rate may be explained only if the n exponent is greater than 1.3. Lague [ 2014 ] reanalyzed data from a range of studies and suggested that m ≈1 and n ≈2 in most cases, with a few exceptions. Data from the Siwaliks in the Himalayas [ Lavé and Avouac , 2001 ], Eastern Tibet [ Ouimet et al. , 2009 ], and the Mendocino triple junction in the Western USA [ Snyder et al. , 2000 ] suggest values of m = 0.55, 0.85, and 2 and n = 1.1, 1.7, and 4, respectively [ Lague , 2014 ]. Whipple et al. [ 2000 ] argued that the n exponent depends on the dominant erosion process and varies between ∼2/3 and ∼5/3. Royden and Perron [ 2013 ] used transformed river long profiles along with previously determined uplift rates for the Rio Torto in the central Apennines to estimate an n value of approximately 0.5. Mudd et al. [ 2014 ] analyzed the gradient ( M χ ) of these transformed profiles to estimate 0.52 < n < 0.82 for the Rio Torto. A recent global study of the stream power law parameters by Harel et al. [ 2016 ] found that in most landscapes the exponent n is greater than 1.

The relationship betweenand erosion rate can be predicted by combining models of river incision with models of hillslope sediment transport []. Here we model fluvial incision using the stream power model, a common detachment‐limited scenario (equation 1 ). Depending on the values chosen for the exponentsand, this model can represent fluvial erosion rate as a function of shear stress, for example, or unit stream power []. There are significant limitations to this detachment‐limited model formulation. It assumes that channel width scales with contributing area, and it does not take into account the effects of sediment flux or the impact of stochasticity and thresholds, all of which can modulate fluvial incision for a given channel geometry []. However,] found that sediment flux and threshold effects can be cast in the general form of equation 1 . It is often used to examine fluvial response to changing climatic and tectonic conditions, for example, by solving for the relationship between channel slope and contributing area, assuming uniform incision []:

Theoretically, the n exponent in the detachment‐limited incision model (equation 1 ) may be calculated using the relationship between the mean length of overland flow (inversely proportional to drainage density) and the erosion rate, if known. We fitted a power law to the relationship between mean length of overland flow and erosion rate for the two field sites with CRN‐derived erosion rates available: Feather River, California, and Boulder Creek, Colorado. We used the gradient of the regression, α , to calculate the n exponent based on equation 6 .

The size of the window radius is determined through identifying scaling breaks in the interquartile range and standard deviation of the curvature [ Lashermes et al. , 2007 ; Roering et al. , 2010 ; Hurst et al. , 2012 ]. This ensures that curvature is sampled over a length scale characteristic of hillslope to valley transitions. Mean C HT was computed for each third‐order basin. The relationship between C HT and D d was then examined across each field site.

Ridgetops were mapped as the intersecting margins of basins from zeroth stream order and upward, following the methodology of] and]. Only hilltops internal to each study basin were considered, in order to ensure thatwas adjusted to the erosion rate within each basin. Curvature was calculated using polynomial surface fitting with a circular window radius of 7 m []. The polynomial surface has the formwhere curvature () and slope () can be determined from the fitted coefficients

Mean C HT may be used to infer the distribution of erosion rates across the landscape [ Roering et al. , 2007 ; Hurst et al. , 2012 ]. Hurst et al. [ 2012 ] demonstrated that mean C HT continues to vary linearly with erosion rate after hillslope gradient has become insensitive to increasing erosion rate. Mean C HT has been demonstrated to respond rapidly to changing channel steepness in soil‐mantled landscapes and can therefore be used as a proxy for erosion rate in areas where CRN‐derived erosion rates are not available [ Hurst et al. , 2013a ]. In order to ensure our landscapes were dominantly soil mantled, we calculated the percentage of ridgetops that were soil mantled in each field site using the surface roughness algorithm described in Milodowski et al. [ 2015a ]. We detect patches of bedrock from our lidar DEMs, using a surface roughness ratio of 0.015 as the threshold for bedrock, following Milodowski et al. [ 2015a ]. The roughness ratio is a measure of the local variability of the vectors normal to the topographic surface, which has been shown to correspond to the outcrop of bedrock [ Milodowski et al. , 2015a ]. The percentage of ridgetops identified as soil mantled are reported in Table 2 and is ≥70% for every field site. In some circumstances, C HT may not reflect the variability of erosion rates across the landscape. In transient landscapes, some basins may contain knickpoints, leading to differing erosion rates within the same basin. Therefore, hilltops connected to the channel above and below the knickpoint may not be adjusted to the same channel incision rate [e.g., Mudd and Furbish , 2007 ; Reinhardt et al. , 2007 ; Anderson et al. , 2012 ; Hurst et al. , 2012 , 2013a ]. In addition, the presence of landslides in some high‐relief basins may lead to decoupling of hilltops from the channel network. To avoid such issues, we visually excluded basins with landslides or knickpoints evident from the lidar DEM.

For each basin of interest, we then extracted the total length of channels via the DrEICH method and divided it by the basin contributing area to calculate the D d (expressed in m/m 2 ). We extracted the drainage density for two different sets of basins: all basins with CRN‐derived erosion rates where these were available (the Feather River and Boulder Creek field sites), and all third‐order basins for every field site to investigate the relationship between D d and mean C HT .

The DrEICH algorithm first identifies concave portions of the landscape using a curvature threshold, which is calculated using the quantile‐quantile method of Passalacqua et al. [ 2010a ]. First of all, the DEM is smoothed using optimal Wiener filtering, which distinguishes the large‐scale signal of the fluvial‐hillslope system from microtopographic noise [ Pelletier , 2013 ]. After smoothing the DEM, the curvature threshold is calculated based on a quantile‐quantile plot of the distribution of curvature in each landscape (for more details on this methodology, see Passalacqua et al. [ 2010a , 2010b , 2012 ]). The DrEICH algorithm identifies the upstream extent of fluvial incision within the valleys based on χ ‐transformed longitudinal profiles [ Clubb et al. , 2014 ]. It assumes that the channel profile will be made up of two segments in χ ‐elevation space: a linear channel segment and a nonlinear hillslope segment. The channel head in each valley is calculated as the transition point between the best fit linear channel segment and nonlinear hillslope segment. The DrEICH algorithm was tested against 167 field‐mapped channel heads from a variety of landscapes by Clubb et al. [ 2014 ] and was found to accurately reproduce the field‐mapped channel networks when compared to other channel extraction methods. Our analytical model described in section 2 relies on equating the channel and hillslope gradient at the channel head. In order to test whether these were comparable, we extracted the gradient 2 m above and below each channel head for every field site and calculated the percentage difference between the two gradients. We then calculated the mean percentage difference across each landscape (Table 2 ). For each field site there was less than 25% mean difference in the gradients above and below the channel heads, suggesting that our assumption of equating the slopes in our analytical model is valid.

We ran the test described by Mudd et al. [ 2014 ] on two catchments at each field site using 10 m DEMs derived from the United States Geological Survey's National Elevation Dataset. We used the National Elevation Dataset instead of the lidar data at each field site to provide a larger area for calculation of the m / n ratio and to reduce computational cost. The Mudd et al. [ 2014 ] algorithms require four user‐defined parameters: the target skip value, the standard deviation of the elevation data ( σ ), the minimum segment length, and the number of target nodes (for a detailed description of each of these parameters, see Mudd et al. [ 2014 ]). The value of these parameters can influence the result of the m / n analysis. We performed a sensitivity analysis by changing each parameter and examining the variation in m / n ratios extracted. In total we ran 54 combinations of the parameters to determine the best fit m / n . We varied the skip parameter between 1 and 3, the number of target nodes between 80 and 100, and the minimum segment length between 10 and 20 nodes. We used a σ value of 3 m for all field sites, as analyses performed by Mudd et al. [ 2014 ] showed that the most reliable m / n ratios were calculated when σ values were ≤ 3 m. We used the mean value of the sensitivity analyses as the best fit m / n for the sites; mean, median, and standard deviation of the analyses are reported in Table 4 .

D d was calculated for each of the study basins using the DrEICH method, a channel head extraction algorithm [ Clubb et al. , 2014 ]. The DrEICH method extracts channel heads based on transforming river profiles into χ ‐elevation space [ Perron and Royden , 2013 ], identifying the upstream transition between fluvial and hillslope processes [ Clubb et al. , 2014 ]. Perron and Royden [ 2013 ] showed that river profiles are linearized when transformed into χ ‐elevation space. The DrEICH algorithm identifies channel heads as the point at which these profiles become nonlinear, representing the transition to hillslope processes [ Clubb et al. , 2014 ]. In order to extract channel networks for each field site via the DrEICH methodology, the m / n value for the landscape must be calculated. This was done using the independent statistical collinearity tests described by Mudd et al. [ 2014 ] which assume channel profiles are made up of a number of different segments depending on heterogeneities and spatial variations in incision rate within the river profile. The collinearity test loops through all potential m / n values and performs a piecewise linear regression on the profile. For each regression the Akaike Information Criterion (AIC) is calculated [ Akaike , 1974 ], which measures how well the data fit the regression while penalizing for overfitting. The best fit m / n is selected at the minimum AIC value.

Two sets of study basins were used in the analysis. The first set of study basins included all catchments for which there were CRN‐derived erosion rates available (the two sites above). The second set of study basins, extracted for each of the five sites, included all third‐order basins. We chose to use third‐order basins to sample across a large number of catchments at different erosion rates in each site. We obtained the mean hilltop curvature for each of these basins to use as a proxy for erosion rate where CRN‐derived erosion rates were not available, as previous work by Hurst et al. [ 2012 ] demonstrated that hilltop curvature scales linearly with erosion rate.

In order to examine the relationship between D d and erosion rate, published CRN‐derived erosion rates were compiled from two sites: Feather River, CA [ Riebe et al. , 2000 ; Hurst et al. , 2012 ] and Boulder Creek, CO [ Dethier et al. , 2014 ] summarized in Table 3 . No CRN data are available for the other sites. The catchment‐averaged erosion rates for the Feather River field site were derived using 10 Be concentrations from fluvial sands, assumed to have minimal storage in the fluvial system, by Riebe et al. [ 2000 ] and Hurst et al. [ 2012 ]. A total of 21 CRN‐derived erosion rates are available for the Feather River. These erosion rates span an order of magnitude, varying from 12.5 ± 1.4 mm/kyr to 253.8 ± 66.6 mm/kyr. CRN‐derived erosion rates for the Boulder Creek field site were also calculated by measuring 10 Be concentrations from quartz in alluvial channel sediments by Dethier et al. [ 2014 ]. Within the Boulder Creek catchment there are 12 basins for which a CRN‐derived erosion rate is available, ranging from 14.97 ± 1.25 mm/kyr to 62.92 ± 5.96 mm/kyr.

As well as testing our predictions on model landscapes, we report D d and hilltop curvature for five field sites with 1 m resolution lidar data: two sites in California, one site in Colorado, one site in New Mexico, and one site in Idaho (Figure 3 ). These sites were chosen based on the following criteria: (i) the availability of 1 m resolution lidar data, (ii) relatively uniform lithology across the site, and (iii) a gradient in erosion rates across the landscape, either measured or inferred based on highly variable slopes and ridgetop curvatures. Table 2 summarizes the mean annual temperature and mean annual precipitation of each site (PRISM Climate Group, http://prism.oregonstate.edu ), the underlying lithology, and the elevation range. These sites are predominantly soil mantled, as shown in Table 2 , although some of these sites have bedrock outcrops, such as the Guadalupe Mountains and San Gabriel Mountains, which become bedrock dominated in the rapidly eroding parts of the landscape [ DiBiase et al. , 2010 ]. Debris flows are also prevalent at the San Gabriel Mountains site [ DiBiase et al. , 2012 ].

A potential limitation of using the CHILD model is that the hillslope sediment transport term does not account for flow width (equation 3 ) [ Howard , 1994 ; Pelletier , 2010 ]. Pelletier [ 2010 ] suggested that, if the grid resolution of the model ( δ ) is greater than that of the valley width ( w ), the diffusive transport term should be scaled by a ratio of δ / w . This is not accounted for in CHILD, suggesting that our model runs may underestimate the colluvial deposition rate and potentially predict higher drainage densities. In order to test the sensitivity of our model scenarios to grid resolution, we ran our steady state scenarios at grid resolutions of 2.5 m, 7.5 m, and 10 m, along with the original 5 m runs. We found that the predicted relationship between D d and E was independent of grid resolution ( supporting information Figures S3 to S5 ).

In addition to steady state runs, we examined a transient scenario for three values of n : n = 1, n = 2, and n = 0.7 with linear hillslope sediment transport. We ran these scenarios to test whether such landscapes conform to the same theory, and whether spatial changes in drainage density resulting from varying erosion rates can be detected. These were performed with a larger model domain (2000 m by 2000 m) in order to examine the variation in drainage density across different basins in the same landscape (Figure 2 c). The model was run at a low uplift rate (40 mm/kyr) for 20 × 10 6 years, and then the uplift rate was increased to 320 mm/kyr for 1 × 10 6 years. This allowed us to compare basins responding to different uplift rates in the same model landscape. Our transient scenarios were analyzed following the same procedures (sections 3.4 and 3.5 ), with the drainage density and hilltop curvature extracted for different basins in the domain.

In order to test whether the theory outlined by equations 1 – 8 is applicable in 2‐D, we analyze a series of model landscape evolution scenarios created using the Channel‐Hillslope Integrated Landscape Development (CHILD) model [ Tucker et al. , 2001 ]. In the model, topography evolves based on a combination of fluvial incision using the stream power law (equation 1 ), as well as linear and nonlinear diffusive hillslope sediment transport (equations 3 and 7 ). Our model domain is 500 m by 500 m with a node spacing of 5 m, comparable to the size of the catchments extracted and the DEM resolution of the real data sets, respectively. Although the real data sets have a DEM resolution of 1 m, we could not run our numerical models at the same resolution due to computational constraints. Our domain has one boundary set to a fixed elevation ( z = 0 m) and three boundaries set to no flux. We detail the model setup and the values of all parameters used in Appendix A and Table 1 . We ran three different series of scenarios with different values of n . Our first scenario set n = 1 and m = 0.5, where erosion is proportional to specific stream power. We then kept all other parameters constant while changing the value of n . We ran further scenarios with n = 0.4, n = 0.7, and n = 2. We chose these values of n in order to run scenarios with different values of α . For each scenario, several runs were performed with uplift rates varying between 10 and 320 mm/kyr. We ran each simulation for 5 × 10 6 years to allow the topography to reach steady state, which we determined as occurring when the volume of rock above sea level in the modeled domain became constant. Figures 2 a and 2 b show examples of the topography generated during these runs with high and low erosion rates. The CHILD model uses a triangulated irregular network (TIN). We converted the output TINs to rasters and performed the same topographic analysis as with the real data sets, described in the following sections. We extracted the channel network and calculated the drainage density for the whole catchment generated; the methodology for channel extraction is described in section 3.4 . We also computed mean C HT for each run following the methodology outlined in section 3.5 . We further wished to examine the effect of different hillslope transport laws on the relationship between D d and C HT for varying values of n . Therefore, all our steady state scenarios have been run twice, using linear and nonlinear hillslope sediment flux laws (equations 3 and 7 , respectively).

The relationship between D d and erosion rate can theoretically be used to calculate the n exponent. The scatterplots of D d against CRN‐derived erosion rate show a positive relationship for the Feather River and Boulder Creek (Figure 8 ). Furthermore, there is also a positive relationship between D d and mean C HT for four out of the five field sites (Figure 10 ). This suggests that n > 1 at each of these sites. We can use the relationship between the mean length of overland flow (inversely proportional to D d ) and erosion rate to calculate n if the m / n ratio is known using equation 6 . However, we find that varying the gradient of the regression within the range of acceptable values results in a wide variation in the calculated n exponent. Therefore, it was not possible to further constrain the value of the n exponent using this technique.

Mean hilltop curvature was calculated for every third‐order basin in each of the five study sites and compared to drainage density. Figure 9 shows an example of the spatial distribution of hilltop curvature and drainage density for the Guadalupe Mountains field site. Scatterplots were created of D d against mean C HT for each of these basins (Figure 10 ), and the data were also binned with a bin width of 0.005 m −1 . A power law relationship was fit through all of the data points for each field site, and the p value and R 2 were reported (see Figure 10 ). A significant positive relationship between mean C HT and D d was observed for four out of the five field sites analyzed in this study, with the exponent in the power law relationship varying between 0.15 and 0.6. There was no significant relationship observed between C HT and D d for the San Gabriel Mountains field site, with a p value of 0.02 and an R 2 of 0.18. Mean C HT may only be used as a proxy for erosion rate if the ridgetops are soil mantled. Therefore, the percentage of bedrock ridgetops as a function of mean C HT was also plotted for each field site (Figure 11 ). We found a positive linear relationship between the percentage of bedrock ridgetops and mean C HT for each field site. The vast majority of basins in each site had a low percentage of bedrock ridgetops (Figure 11 ), although one basin in the Bitterroot National Forest site had an anomalously high percentage (70%).

For our real landscapes, we created scatterplots of D d against CRN‐derived erosion rates for the two field sites with available CRN data: Feather River, CA, and Boulder Creek, CO. A power law regression was fitted to the raw data for each of the field sites (Figure 8 ). A regression was deemed to be significant if the p value was less than 0.01 (99% confidence interval). Figure 8 shows that for both the Feather River and Boulder Creek field sites there is a positive relationship between erosion rate and D d . The regressions for the Feather River and Boulder Creek both have p values <0.01, and R 2 values of 0.76 and 0.82, respectively. The exponents on the power law relationships ( α ) are 0.91 and 1.37, respectively.

For each of the steady state modeling scenarios ( n = 0.4, n = 0.7, n = 1, and n = 2) the relationship between D d and uplift rate was plotted for both linear and nonlinear hillslope sediment transport (Figures 4 and 5 ). Figure 4 shows that in the scenarios with linear hillslope sediment transport, there is a positive relationship between drainage density and uplift rate (and therefore erosion rate as the scenarios are at steady state) for n = 2; a negative relationship for n = 0.7 and n = 0.4; and that drainage density is invariant with uplift rate for n = 1. The negative relationship between D d and U is steeper for n = 0.4 than n = 0.7, as would be expected from equation 6 and Figure 1 . We fit a linear regression to the relationship between hilltop curvature and uplift rate based on the predictions of erosion rate and mean C HT set out by Roering et al. [ 2007 ] and following the methodology of Hurst et al. [ 2012 ], shown in Figure 6 . We find a significant positive relationship between mean C HT and uplift rate for both the linear and nonlinear hillslope sediment transport scenarios. These results in an ideal landscape mirror those from our theory (section 2 ) and justify the use of C HT as an indicator of erosion rate in soil‐mantled landscapes. The same general trends between D d and U are apparent for nonlinear sediment transport (Figure 5 ). Figure 1 b suggests that the relationship between D d and uplift rate should be steeper for n = 2 for nonlinear sediment transport, which is not evident from our modeling results. Our transient simulations (Figure 7 ) show the same trends as our steady state runs, suggesting that the same theory can be applied to transient landscapes.

5 Discussion

The results of our landscape evolution modeling show that the theoretical concepts outlined in section 2 are applicable in a 2‐D domain. We find that the nature of the relationship between D d and E in our model scenarios varies with the value of the slope exponent in the detachment‐limited incision model, n. Our steady state modeling runs (Figures 4 and 5) show that if n > 1, there is a positive relationship between D d and C HT , whereas if n < 1 there is a negative relationship. Our modeling of transient scenarios also supports this, showing that the theory is still applicable in these landscapes. There was a significant trend between drainage density and mean hilltop curvature for our transient model runs with n = 0.7 and n = 2, with R2 values of 0.85 and 0.58, respectively (Figure 7).

Our modeling results also have implications for examining the impact of nonlinear hillslope sediment transport on length scales in landscapes. As relief increases and hillslopes approach threshold gradients, hillslope sediment transport becomes increasingly nonlinear, as sediment flux becomes dominated by mass wasting and landslides [Roering et al., 1999]. Our landscape evolution modeling scenarios where hillslope sediment transport was nonlinear (Figure 5) exhibited the same relationships between drainage density and erosion rate as scenarios in which hillslope sediment transport was linear. However, our analytical solution (Figure 1) predicted a steeper relationship between D d and uplift rate with nonlinear transport for n = 2. This may be due to noise in our modeling results (Figures 4-7). This noise may be caused by the extraction of hilltop curvature from the model domain, as well as the loss of information when transforming a TIN network onto a regular grid. The model has a grid spacing of 5 m due to computational constraints, but this resolution may not be fine enough to perfectly capture the variation in curvature along the ridgetops. Despite the noise, a clear significant trend between drainage density and uplift rate is observed from the steady state model runs using linear and nonlinear sediment flux laws, with R2 values ranging from 0.88 to 0.96. (Figures 4 and 5).

It may be expected that the effect of the nonlinear sediment flux law will increase with erosion rates in higher‐relief landscapes. In our CHILD model runs we tested a maximum erosion rate of 320 mm/kyr in order to compare these results to our real landscapes with CRN‐derived erosion rates in this order of magnitude (Table 3). At higher erosion rates, where landscapes transition from soil mantled to bedrock dominated, C HT cannot be used as an indicator of erosion rate as soil production can no longer keep pace with transport rates [Hurst et al., 2012]. Therefore, the results of our study are applicable to landscapes with soil‐mantled ridgetops where C HT can be used as a proxy for erosion rate across the landscape. Figure 11 shows that in each field site with the exception of Bitterroot National Forest, the majority of basins had below 20% of ridges identified as bedrock. In the Bitterroot site the majority of basins had less than 35% bedrock ridgetops, although with more variability than the other field sites. This may lead to more noise in the relationship between D d and C HT in this site, although a significant positive relationship was still observed (Figure 10). In regions with much higher erosion rates, a positive relationship between D d and erosion rate may not be observed. Previous studies by Oguchi [1997] in the mountainous region of central Japan, and by Talling and Sowter [1999] in the Southern San Joaquin Valley, California, found lower drainage densities corresponding to higher relief. These authors concluded that the dominance of mass‐wasting processes on steep slopes in these regions resulted in a negative relationship between D d and relief. In contrast, Sangireddy et al. [2016] found that across a wide range of humid landscapes D d was positively correlated with relief.

We also tested our predictions on real landscapes, using lidar‐derived DEMs for five field sites in the USA. Our results show a positive relationship between D d and erosion rate, using CRN‐derived erosion rates for two of the field sites, and mean C HT for four out of the five sites. Drainage density has profound implications for the transit time of runoff through catchments, and maximum storm runoff has been demonstrated to be a function of drainage density [Chorley and Morgan, 1962; Gregory and Walling, 1968]. Our results suggest that increasing erosion rates will, therefore, result in more rapid catchment response to storms or precipitation events. Furthermore, an increase in drainage density with erosion rate may increase the rapidity of sediment supply to the fluvial network. This is an important control on downstream fluvial morphology, influencing the transition between braided and meandering channels, for example, as stable meandering channels are more likely to develop with low rates of sediment transport [Church, 2006]. Furthermore, based on our landscape evolution modeling results, this positive relationship between drainage density and erosion rate is consistent with a value of n in the stream power law (equation 1) greater than 1. In landscapes with linear hillslope sediment transport, if n is greater than 1 and other parameters are constant, as slope increases (in response to an increase in uplift, for example) fluvial processes will out‐compete diffusive processes. This would lead to channel incision occurring further upstream and an increase in D d . However, where hillslope sediment transport is nonlinear, this relative efficiency of advective and diffusive processes may also depend on the critical gradient, S c . We set S c constant in our modeling scenarios but acknowledge that the value of S c may vary spatially [Grieve et al., 2016a, 2016b], which could affect the observed relationship between drainage density and erosion rate.

The value of the n exponent also has important implications for how the landscape responds to transient forcing. The slope of river profiles may be used to extract information about the uplift history of the channel [Pritchard et al., 2009; Roberts and White, 2010]. However, complete uplift histories can only be extracted from channel profiles if n = 1, when knickpoint retreat rates should be independent of erosion rate. Royden and Perron [2013] demonstrated that if n > 1, rapid incision signals should propagate upstream more rapidly than slow incision (with the converse true for n < 1). Steep segments in river profiles are predicted to lengthen when n > 1, consuming lower gradient segments and therefore progressively destroying the record of the preceding uplift history. Our results are consistent with n > 1 for four of the field sites analyzed, and n≥1 for all of the field sites. This agrees with Lague [2014], who found that n≈2 in the majority of cases. Our results therefore imply that channels in these landscapes will be imperfect recorders of tectonic forcing, and complete uplift histories cannot be extracted from these river profiles.

D and K has been shown to exert a first‐order control on valley spacing [e.g., Perron et al., 2008 2009 Sweeney et al., 2015 Perron et al. [ 2008 m, n, and relief (ζ). They define the Peclet number, Pe as (12) l is the horizontal length of a drainage basin. Perron et al. [ 2008 Pe through an increase in relief, ζ, if n > 1. This leads to narrower valley spacing and increased drainage density. However, if n = 1, Pe is independent of relief. This theory is consistent of the results of our study, where we find that a positive relationship between drainage density and erosion rate is consistent with n > 1. The competition between the parametersandhas been shown to exert a first‐order control on valley spacing [e.g.,].] showed that valley spacing is also predicted to vary with the parameters, and relief (). They define the Peclet number,aswhereis the horizontal length of a drainage basin.] suggest that, assuming the other parameters are constant, higher erosion rates will increasethrough an increase in relief,, if> 1. This leads to narrower valley spacing and increased drainage density. However, if= 1,is independent of relief. This theory is consistent of the results of our study, where we find that a positive relationship between drainage density and erosion rate is consistent with> 1.

However, a key assumption of our study is that the D and K parameters in equations 1 and 3 are constant. The competition between these parameters has been shown to exert a first‐order control on valley spacing [e.g., Perron et al., 2008, 2009; Sweeney et al., 2015]. However, the values of D and K may vary both spatially and temporally. The hillslope diffusion coefficient, D, is a function of hillslope sediment properties, such as soil thickness, cohesion, and grain size [Furbish et al., 2009]. Both soil thickness and grain size are thought to vary with erosion rate [Heimsath et al., 1997; Attal et al., 2015; Riebe et al., 2015]. D has also been shown to vary with climate through controls on soil transport processes [e.g., Carson and Kirkby, 1972; Yoo et al., 2005], and lithology, which affects material properties and soil particle sizes [Hurst et al., 2013b]. If soil thickness decreases with erosion rate, the models of depth‐dependent sediment transport suggest that D may also vary with erosion rate [e.g., Braun et al., 2001]. The channel erodibility coefficient, K, is a function of many parameters such as lithology, climate, sediment cover, and channel width [Whipple and Tucker, 1999]. K may vary with erosion rate through channel width adjustments, as channels have been demonstrated to narrow in response to steepened reaches from increased uplift rates [e.g., Finnegan et al., 2005; Amos and Burbank, 2007]. If n = 1, then equation 6 simplifies so that D d is dependent on D/K. This suggests that, if n = 1, a positive relationship between D d and E may result from a decreasing D/K ratio with erosion rate (Figure S2, supporting information). This may be caused by a decrease in D; an increase in K; or K increasing faster than D such that the ratio decreases. However, as no field evidence has demonstrated how K or D vary with E, these three scenarios cannot be distinguished. With these limitations, we suggest that our results are consistent with the hypothesis that n > 1 for four out of the five field sites, although acknowledge that there may be other possible explanations for the observed relationship.

A further assumption of our analytical predictions is that of detachment‐limited fluvial incision (see equation 1). Detachment‐limited incision assumes that the erosion rate is related to the shear stress, velocity, or power of overland flow, and that sediment is transported by the channel without being deposited. It has been assumed in many studies modeling evolution of soil‐mantled landscapes [e.g., Howard, 1994, 1997; Perron et al., 2008, 2009]. Other studies, however, suggest that erosion in soil‐mantled landscapes is transport limited, where erosion rate is proportional to the divergence of sediment flux [e.g., Tucker and Bras, 1998; Simpson and Schlunegger, 2003; Istanbulluoglu et al., 2003]. Pelletier [2012] demonstrated through analysis of field measurements, along with numerical modeling, that at small spatial scales, both detachment‐ and transport‐limited conditions may apply depending on the texture of the eroding soil. The assumption of detachment‐limited conditions is a simplifying one that we make in this study to generate simple predictions that are testable against our real landscapes. However, Tucker and Bras [1998] present a purely transport‐limited model of the drainage density in soil‐mantled landscapes and predict similar relationships between D d and E as we find in our detachment‐limited model.

Previous studies have suggested that the underlying lithology has an effect on D d [Oguchi, 1997]. Three of the sites analyzed (Feather River, Boulder Creek, and the San Gabriel Mountains) were situated on granitic lithologies; the Guadalupe Mountains site was primarily composed of limestone; and the Bitterroot National Forest site was composed of schist and gneissoid bedrock. Despite these variations, the relationship between D d and erosion rate was positive for four of these sites (Figures 8 and 10). The San Gabriel Mountains is the only site to show no relationship between drainage density and erosion rate. DiBiase et al. [2012] analyzed the same DEM and found that fluvial drainage density decreased with increasing erosion rate, while colluvial channels become denser, leading to the total drainage density remaining constant across the landscape. This contrasts with our analysis, as we found that fluvial drainage density was invariant with erosion rate. The difference between these results may be due differences in channel extraction: DiBiase et al. [2012] used slope‐area plots to identify fluvial channels, whereas in our analysis we used the DrEICH algorithm, which identifies channels based on transformed river long profiles. Our results are consistent with the n value in this landscape being close to 1, as implied by our numerical modeling results. As shown by the San Gabriel Mountains site, the presence of colluvial channels in steep landscapes formed through debris flow processes [Stock and Dietrich, 2003] may complicate the results of our analysis. These colluvial channels can impact the results of channel extraction algorithms and, therefore, the calculation of drainage density across the landscape, as the DrEICH algorithm is focused on identifying the extent of the fluvial channel network.