How long is a hillslope?

Hillslope length is a fundamental attribute of landscapes, intrinsically linked to drainage density, landslide hazard, biogeochemical cycling and hillslope sediment transport. Existing methods to estimate catchment average hillslope lengths include inversion of drainage density or identification of a break in slope–area scaling, where the hillslope domain transitions into the fluvial domain. Here we implement a technique which models flow from point sources on hilltops across pixels in a digital elevation model (DEM), based on flow directions calculated using pixel aspect, until reaching the channel network, defined using recently developed channel extraction algorithms. Through comparisons between these measurement techniques, we show that estimating hillslope length from plots of topographic slope versus drainage area, or by inverting measures of drainage density, systematically underestimates hillslope length. In addition, hillslope lengths estimated by slope–area scaling breaks show large variations between catchments of similar morphology and area. We then use hillslope length–relief structure of landscapes to explore nature of sediment flux operating on a landscape. Distinct topographic forms are predicted for end‐member sediment flux laws which constrain sediment transport on hillslopes as being linearly or nonlinearly dependent on hillslope gradient. Because our method extracts hillslope profiles originating from every ridgetop pixel in a DEM, we show that the resulting population of hillslope length–relief measurements can be used to differentiate between linear and nonlinear sediment transport laws in soil mantled landscapes. We find that across a broad range of sites across the continental United States, topography is consistent with a sediment flux law in which transport is nonlinearly proportional to topographic gradient. © 2016 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd.

A fundamental measure of hillslopes is their down-slope length. Hillslope length exerts a primary control on the rate of sediment transport through a catchment, which modulates flows of pollutants, nutrients and water, and in addition influences rates of soil erosion at both the hillslope and catchment scale (Mathier et al., 1989;Dunne et al., 1991;Gabriels, 1999;Liu et al., 2000;Kinnell, 2009;Thompson et al., 2010). Consequently, in a catchment which has short hillslope lengths, the signal of a storm event in the catchment's headwaters will be transmitted more rapidly to the outlet of the catchment than would be the case in an otherwise similar catchment with longer hillslopes (Dunne et al., 1991).
Work to quantify the rate of sediment accumulation on hillslopes by estimating colluvial hollow filling rates requires accurate constraint on the length of the hillslope upslope of the accumulation zone (Reneau et al., 1989;Reneau and Dietrich, 1991). Hillslope length also acts to constrain debris flow hazard in soil mantled landscapes: a scaling relationship between the area of a landslide scar and the volume of its deposit has been demonstrated (Guzzetti et al., 2009;Klar et al., 2011) and the length of a hillslope exerts a fundamental constraint on the maximum source area of landslides (Hurst et al., 2013c).
Hillslope length is governed by the amount of dissection of a landscape by a channel network, a process governed by a range of factors including climate, erosion rate and hillslope morphology (Chorley, 1957;Abrahams and Ponczynski, 1984;Montgomery and Dietrich, 1989;Rinaldo et al., 1995;Oguchi, 1997;Moglen et al., 1998;Tucker and Bras, 1998;DiBiase et al., 2012) and additionally by the roughness of the topographic surface (Dunne et al., 1991;Thompson et al., 2010). On highly diffuse hillslopes, with little vegetation cover to increase surface roughness, the length of overland flow will be reduced, when contrasted with a hillslope characterized by a large amount of microtopography or vegetation (Dunne et al., 1991;Thompson et al., 2010).
The importance of quantifying hillslope length has been recognized for many decades, but doing so has been a challenge: Horton (1932) commented, 'The exact determination of the average distance which rain must travel overground to reach the stream-channels is impracticable.' Owing to advances in both topographic data resolution and computational power Horton's assertion is no longer true.
A number of methods have been proposed to measure hillslope length. It can be estimated through interpretation of slope-area plots  and analysis of drainage density data (Tucker et al., 2001;DiBiase et al., 2012) and by modelling flow paths along hillslopes (Tucker et al., 2001;DiBiase et al., 2012;Hurst et al., 2012). An estimate of hillslope length averaged over a landscape can also be approximated through spectral analysis of characteristic length scales found in high resolution topography (Perron et al., 2008), however we do not consider this method, as it yields a single length scale for a landscape, which cannot be readily evaluated against the hillslope and basin scale measurements of other techniques.
Since sediment flux and overland flow are in part governed by hillslope length (Dunne et al., 1991), we reason that characterizing hillslope length as a flow path at the hillslope, rather than catchment scale, should be the most reliable description. We therefore set out to compare the length of flow paths against other previously reported measures of hillslope length and consider the distribution of hillslope length measurements across four field sites in the continental United States.
The rate of hillslope sediment flux is a crucial component in understanding landscape evolution, as it is integral in determining both the geometry of hillslopes (Carson and Kirkby, 1972;Dietrich et al., 2003) and their response to climatic and tectonic forcing (Armstrong, 1980(Armstrong, , 1987Rosenbloom and Anderson, 1994;Arrowsmith et al., 1996;Fernandes and Dietrich, 1997;Roering et al., 2001Roering et al., , 2004Mudd and Furbish, 2005). Here we demonstrate that our method of measuring hillslope length via flow paths can be used to test two common sediment flux laws, by comparing the hillslope lengthrelief structure of landscapes with theoretical predictions made by each flux law.
Spatially continuous soil thickness data is not available for our field sites and particle motion based models make a wide range of topographic predictions Tucker and Bradley, 2010;Furbish and Roering, 2013), as they have no analytical solution, meaning that neither of these groups of flux laws can be falsified using topographic data. However, it is possible to falsify the nonlinear or linear models of sediment flux as these two laws predict distinct relief structures which can be measured using high resolution topographic information. Previous tests of the appropriate sediment flux law for a given landscape typically combine topographic analysis and numerical modelling (Roering et al., 2001;Roering, 2008), topographic analysis and field measurements (Roering et al., 1999;Heimsath et al., 2005), field measurements and numerical modelling (Braun et al., 2001;Herman and Braun, 2006;), or a combination of all three of these techniques (Pelletier et al., 2011). Here we demonstrate a test of flux laws using topography alone.
One approach to examining the topographic outcome of different sediment flux laws is to calculate the predicted relief structure of a landscape. Roering et al. (2007) solved a statement of mass conservation along a one dimensional hillslope for hillslopes obeying a nonlinear flux law of the form (Andrews and Bucknam, 1987;Roering et al., 1999Roering et al., , 2001Roering et al., , 2007: where q s (in dimensions Length 2 Time À1 , dimensions henceforth denoted in square brackets as [L]ength, [M]ass and [T] ime) is a volumetric sediment flux per unit contour length, K [L 2 T À1 ] is a sediment transport coefficient, S [dimensionless in L/L] is hillslope gradient and S c [dimensionless] is a critical hillslope gradient; as S approaches S c , q s asymptotically progresses towards infinity (Figure 1(a)). Relief, R [L], defined as the elevation difference between the hilltop and the channel can be determined from Equation 1 in the case of steady-state hillslopes, where erosion balances uplift and soil thickness does not change over time . Here, we employ the steady state definition of Mudd and Furbish (2004), whereby a hillslope which retains a steady form with regard to baselevel (the channel at the toe of the hillslope) is deemed to be in steady state (also known as steady topographic form). Roering et al. (2007) demonstrated that on hillslopes with steady topographic form, R is a function of the erosion rate E [L T À1 ], the sediment transport coefficient K and the hillslope length L H [L] where sediment transport can be described by Equation 1. Furthermore, Roering et al. (2007) found by 1040 S. W. D. GRIEVE ET AL.
normalizing the equation describing hillslope form for these hillslope properties, all steady state hillslopes obeying Equation 1 should fall on a single curve when nondimensional relief R * = R/(S c L H )) is plotted against nondimensional erosion, E * = (À2C HT L H )/S c , where C HT [L À1 ] is the curvature of the hilltop. Hurst et al. (2012) was able to show that in a landscape with a wide range of erosion rates, along the Feather River in northern California, ridgeline-averaged E * vs R * data of hilltops fell on the curve predicted by Roering et al. (2007), thus providing a test of the nonlinear sediment flux law. One drawback of this approach is that it requires a distribution of erosion rates in a landscape; Roering et al. (2007) tested their theory in landscapes with more spatially homogeneous erosion rates (the Oregon Coast Ranges and Gabilan Mesa, California) and these two landscapes fell upon single points in E * vs R * space. Other efforts at constraining flux laws have examined manually selected hillslope profiles (Rosenbloom and Anderson, 1994;Arrowsmith et al., 1996;Furbish and Roering, 2013) or a combination of manually selected hillslope profiles and high resolution topography (Roering et al., 1999). Such profiles are both time consuming to identify and difficult to reproduce without transect coordinates.
Here we use the dimensional relationship between hillslope length, L H , and relief R, combined with our algorithmic extraction of ridge to valley profiles to examine if topography is more consistent with a particular sediment flux law. One advantage of using dimensional instead of nondimensional length vs relief data is that hillslopes should have a wide range of hillslope lengths so relief and hillslope length may be tested on a large dataset in steadily eroding landscapes. Unlike previous tests of sediment flux laws, this procedure can be performed without field measurements or numerical modelling.
A simple statement of mass conservation on a one-dimensional hillslope, assuming negligible chemical weathering processes, is: where z [L] is the elevation of the surface, U [L T À1 ] is tectonic uplift and ρ s and ρ r [M L À3 ] are soil and rock densities, respectively. At steady state, U = E, ∂z/∂t = 0, Equation 2 reduces to (ρ r /ρ s )E = ∂q s /∂x. To close this equation, we must define a sediment flux law, q s . One end member sediment flux law is a linear sediment flux law (Figure 1(a)) (Gilbert, 1909;Culling, 1960;Ahnert, 1987): which was first used primarily for its simplicity but has been supported by field studies (McKean et al., 1993;Small et al., 1999) and has been demonstrated to describe sediment transport by rainsplash Dunne et al., 2010). We also consider the nonlinear flux law of Equation 1. We do not consider depth dependent sediment flux because we do not have soil thickness data for all our study sites; thus we cannot truly evaluate all possible flux laws. Rather our goal is to assess whether one of two popular flux laws are consistent with the data and show that, at a minimum, flux laws are falsifiable. At steady state, the solution of Equation (2) for linear sediment flux is: where the hilltop is located at x = 0 and the base of the hillslope is at x = L H . The steady state solution of Equation (2) for nonlinear sediment flux is : The relief is found by subtracting the elevation at the hilltop from the elevation at the hillslope base. In the linear case, relief is: and in the nonlinear case steady state relief is: Equations (7) and (8) describe contrasting relationships between hillslope length and relief. If sediment transport can be described by a linear flux law, hillslope relief should go as the square of the hillslope length given constant densities, transport coefficients and erosion rates. On the other hand, relief will be approximately linear as a function of hillslope length if sediment flux can be described by a nonlinear sediment flux law of the form in Equation (1) (Figure 1(b)).
In order to test if theoretical predictions of relief as a function of hillslope length are consistent with real topography, we must first ensure that our topographic measurements are consistent with the assumptions contained within the equations of sediment flux and mass conservation. Equations (7) and (8) are formulated for one-dimensional, soil mantled hillslopes, therefore we must restrict our analysis to hillslope profiles that can be approximated as one-dimensional, and where we can be confident of a persistent soil mantle.
Quasi-one-dimensional hillslope traces are identified by comparing the measured hillslope length to the Euclidean distance between the start and end point of the flow path. If these two values are approximately equal, within a tolerance of 0.02%, we conclude that a trace is planar, and can be described using Equations (7) and (8). In addition, any traces which were initiated on hilltops with a positive hilltop curvature or a slope value of greater than 1.2 are excluded from this analysis, ensuring that hillslopes studied are soil mantled. This filtering can introduce a bias towards shorter hillslopes, as the longer traces have more opportunity to diverge from planarity, and this limitation should be considered when evaluating hillslope length-relief relationships.
Identifying sections of a landscape which approximate onedimensional hillslope profiles is a key challenge of attempting to test one-dimensional models with topographic data, because significant proportions of upland, soil mantled landscapes are characterized by planform curvature, forcing topographic profiles to deviate from a one-dimensional profile (Roering et al., 1999). Our method generates a large amount of traces for a landscape, and can filter them automatically to find all of the quasi-one-dimensional traces in an area. However, such traces can be rare, limiting the validity of tests of flux laws using this method, particularly where only a small area of a landscape is covered by the filtered topographic data.

Topographic processing
In order to implement and test techniques for extracting hillslope length from high resolution topography we must ensure consistency in how topographic processing is performed. The topographic derivatives of slope, aspect and curvature are required for all three of the hillslope length extraction techniques. We calculate these derivatives following techniques outlined by Lashermes et al. (2007); Roering et al. (2010) and Hurst et al. (2012) to fit a polynomial surface to elevation values within a circular window, passed across every cell in the DEM. The radius of the circular window is determined by identifying breaks in values of the interquartile range and standard deviation of curvature as a function of window size. This ensures that the slope values are representative of slope at the hillslope scale rather than a function of the combination of microtopographic variations, such as the roughening of hillslopes from tree throw, and measurement noise produced during LiDAR data acquisition (Roering et al., 2010;Hurst et al., 2012).
The foundation of the hilltop flow routing and drainage density inversion methods is being able to define the location of channels. Without an accurate constraint of the location of the channel network neither method will be accurate. Each of these techniques therefore require precise mapping of the location of the channel head, either through field exploration (Julian et al., 2012;Jefferson and McGee, 2013), a process based method (Clubb et al., 2014), geomorphometric analysis of a DEM (Passalacqua et al., 2010;Orlandini et al., 2011;Pelletier, 2013) or contour map (Oguchi, 1997). Here, we extract channels from the DEM following the DrEICH method outlined by Clubb et al. (2014), using the parameters in Supplementary Table 2. In addition, large floodplains were identified in each landscape as patches of low positive local relief and low slope. Local relief was calculated as the range in elevations in a circular moving window, of the same radius as the window used to fit the polynomial surface, passed over the landscape. These patches were then manually combined to create contiguous patches of floodplain. This floodplain mask was then combined with the channel network to include the floodplain extent, ensuring that only the hillslope portion of each landscape is sampled.
Basins are extracted by identifying junctions in the channel network where the Strahler stream order increases and then identifying all of the upslope pixels for that junction. Drainage area is calculated using the D-infinity algorithm (Tarboton, 1997).
Extracting hillslope length using slope-area analysis Relationships between local slope and drainage area have been used to identify process transitions in landscapes (Montgomery, 2001;Stock and Dietrich, 2003;Roering et al., 2007;Tarolli and Dalla Fontana, 2009;Booth et al., 2013;Tarolli, 2014;Tseng et al., 2015). In areas where hillslope processes dominate, slope increases with drainage area and in areas where fluvial processes dominate slope decreases with increasing drainage area (Montgomery and Foufoula-Georgiou, 1993). Therefore, the transition between the hillslope and fluvial domain can be identified as the first inflection point in slope-area space (Figure 2) (Montgomery and Foufoula-Georgiou, 1993;Hancock and Evans, 2006;Tarolli and Dalla Fontana, 2009). As drainage area increases beyond this transition point, further scaling breaks have been demonstrated to correspond with large-scale landsliding (Tarolli and Dalla Fontana, 2009;Booth et al., 2013;Tarolli, 2014;Tseng et al., 2015).
Slope [dimensionless in L/L] and drainage area [L 2 ] are sampled on a per pixel basis for each catchment in the landscape and the resulting values are placed in logarithmically spaced bins with a width of 0.1 in base 10 logarithmic space, with area measured in units of m 2 . Binned values are then plotted in loglog space and the algorithm searches for the binned point with the maximum slope value. This maximum slope bin is identified as the inflection point and the corresponding drainage area is recorded. We also record the drainage area of the maximum slope of a cubic spline fitted to the binned data following Roering et al. (2007). We have found these two methods pro- Figure 2. Schematic diagram of a slope-area plot showing the predicted relationship between slope and drainage area in a catchment in a soil mantled landscape. Slope increases with increasing drainage area on hillslopes, while in channels, slope decreases as drainage area increases. The inflection point between these two domains, denoted by the dashed line, is identified as the characteristic drainage area, which is transformed into a representative hillslope length for the particular catchment. The second inflection point corresponds to the signature of large-scale landsliding (Tseng et al., 2015). duce indistinguishable results and thus results presented in figures use maximum binned slope values as this is less computationally expensive. The resulting characteristic drainage area, often used to identify the threshold area for channel initiation (Montgomery and Dietrich, 1992;Montgomery and Foufoula-Georgiou, 1993) is then converted into a characteristic hillslope length by dividing the area by the unit contour width , which can be approximated as the DEM resolution (Montgomery and Dietrich, 1989;Moore et al., 1991;Costa-Cabral and Burges, 1994;Mitasova et al., 1996;Tarboton, 1997).
Inversion of drainage density to extract hillslope length Drainage density (D D [L]) is a fundamental landscape parameter which has been shown to vary with climate (Chorley, 1957;Abrahams and Ponczynski, 1984;Rinaldo et al., 1995;Moglen et al., 1998), relief (Schumm, 1956;Oguchi, 1997) and dominant sediment transport process (Montgomery and Dietrich, 1989;Willgoose et al., 1991;Tucker and Bras, 1998). It is defined as the total length of the channel network (L T [L]) divided by the catchment area (A [L 2 ]) (Horton, 1932(Horton, , 1945: The drainage density of a catchment provides a measure of the level of dissection of a landscape, the inverse of which will reflect the catchment average hillslope length. This parameter is described by Horton (1932Horton ( , 1945 and Tucker et al. (2001) as the mean distance water must travel from a random point in a catchment to reach a channel and is considered by Schumm (1956) to be equivalent to the catchment average length of overland flow.
This measurement can be quantified as (Horton, 1932;Schumm, 1956;Tucker et al., 2001): where L H [L] is hillslope length; thus drainage density can be transformed into a measurement of the average flow path length of hillslopes in the catchment.
Using hilltop flow routing to measure hillslope length Building on the work of Hurst et al. (2012), we have, in addition, developed a measure of hillslope length in which hillslope length is defined as the typical travel distance from divide to channel . The simplest method employed to model flow paths is the D8 algorithm, which distributes flow at 45°azimuth angles into the eight neighboring grid cells (Mark, 1984;O'Callaghan and Mark, 1984;Tarboton et al., 1991). This technique operates well when constraining channelized flow (Shelef and Hilley, 2013) and is a commonly used method for extracting a channel network from a DEM (Braun and Willett, 2013). However, when using this algorithm in the hillslope domain, the direct flow paths do not reflect the dispersive nature of overland flow (Tarboton, 1997;Pelletier, 2010;Shelef and Hilley, 2013) and as such the resulting measured hillslope lengths are biased towards short hillslopes. Hurst et al. (2012) modified this approach and modelled flow from hilltop to channel as a point source crossing each DEM cell.
Hilltops are extracted by identifying the edges of basins which share a stream order following techniques outlined by Hurst et al. (2012). From each of these selected hilltop pixels a trace can be run to measure its flow path to the channel.
We then modify the algorithm of Lea (1992) to model the flow of a point source of water from the ridge line down to the channel. Flow is modelled within each DEM cell, flowing from an inlet point on each cell edge to an outlet point on another cell edge, with the flow path within each cell following that cell's aspect angle. Flow is then routed across each cell from the inlet point (or cell center for hilltop pixels where tracing begins) across to its outlet point. This in-pixel routing is repeated recursively, with the flow outlet becoming the subsequent new inlet as the routing proceeds, until the trace has reached a channel or floodplain pixel. If two cells flow into each other the algorithm calculates the secondary aspect of these two pixels and flow is pushed along the cell edge which corresponds to this secondary aspect until it enters another cell whereby flow can continue.
When using aspect angle calculated from a polynomial fit, areas of high topographic noise can generate large sinks in the topography, where flow cannot be routed to the channel. These large areas commonly correspond to low gradient areas of a landscape, such as terraces, valley fills or areas filled by DEM pre-processing to remove pits. It is not possible to allow flow to pass across these zones without extensive smoothing of the landscape or the trace paths themselves, which will destroy any topographic signal which could otherwise be measured. Therefore these traces, which account for <5% of the total number of traces in the nosiest landscapes and <0.5% in smoother landscapes are excluded from all final analysis.
The initial starting point of the trace is given by Lea (1992) as the cell center for computational simplicity. In high resolution topographic data this placement has little influence on the final results as hillslope length is effectively sampled at intervals along each ridgeline equal to the data resolution.
When our routing algorithm is used to measure hillslope length we are implicitly choosing to quantify sediment transport over millenial timescales, because the calculation of aspect averages the land surface over several cycles of geomorphic roughening (i.e. tree-throw pits are considered transient and do not affect millennial scale sediment transport directions). Our objective for this study is to compare hillslope lengths with predictions of geomorphic flux laws operating over millenial time scales (i.e. the time required to change the relief of the hillslope, sensu Mudd and Furbish, 2007) rather than event-based routing of water, so we route flow over topography which is effectively smoothed using an aspect driven routing technique.
The hilltop flow routing algorithm records a small number of basins with short average hillslope lengths, which is consistent with the typical implementation of the algorithm which will measure hillslope lengths near the basin outlet, where the hillslope length approaches the effective channel width. In some basins the only hilltops that are sampled are these interfluves, due to the filtering of hilltops outlined in this section, to ensure that hillslopes sampled are soil mantled. These values are valid measurements; however their bias in smaller basins requires careful analysis of basin average results and motivates the use of median and median absolute deviation values to minimize the influence of these exceptional measurements on the dataset.
Our flow tracing method not only can be used to measure hillslope length, but it can also be used to measure local relief at the hillslope scale, by calculating the difference in elevation between the start and end point of each trace. We use data collected for hillslope length measurements to explore the relationships between hillslope length and relief at the hillslope scale.

Study Sites
We selected four locations in the continental United States to evaluate methods which measure hillslope length in a range of tectonic, climatic and lithologic settings (Figure 3). Each site is covered by high resolution LiDAR data and has been the subject of previous work constraining the hillslope diffusivity in addition to the tectonic and erosional setting. Each dataset was collected by the National Center for Airborne Laser Mapping (NCALM) with point densities ranging from 5.56 to 9.84 points per m 2 . These data were collected with horizontal accuracies ranging from 0.06 m to 0.13 m and vertical accuracies ranging from 0.05 m to 0.35 m (for full accuracy information see Supplementary Table 1). These point clouds were each gridded to 1 m resolution DEMs following standard workflows generated by Kim et al. (2006). This grid resolution was selected as it was the highest resolution that the point clouds could be gridded to, ensuring no variation in results driven by grid resolution changes. Lower resolution data could not be utilized for this study, as the DrEICH algorithm used to extract the channel networks cannot be used on lower resolution data (Clubb et al., 2014). For each study site detailed below, we compared L H data generated using each of the methods detailed earlier.

Coweeta, North Carolina
The Coweeta Hydrologic Laboratory is located in the Southern Appalachian Mountains in North Carolina, USA (Figure 3(a)). The landscape is underlain by Coweeta and Otto Group metasediments (Price et al., 2005) and is characterized by ridge hollow topography (Hack and Goodlett, 1960) with soils approaching 2 m thickness in hollows (Hales et al., 2012). Northern Hardwood forests dominate at high elevations alongside fire intolerant species such as mountain laurel (Kalmia latifolia) and Rhododendron maxima, spread due to anthropogenic fire suppression techniques in the last 80 years (Elliott and Swank, 2008;Hales et al., 2009). Precipitation in Coweeta is dominated by small convective storms with occasional high intensity storms driven by hurricanes (Swift et al., 1988;Hales et al., 2009Hales et al., , 2012, such storms drive slope failures and debris flows within the catchment (Hursh, 1941;Clark, 1987;Wieczorek et al., 2000;Witt, 2005;Hales et al., 2009).
The Southern Appalachians are tectonically quiescent (Oyarzun et al., 1997;Faill, 1998), however, the preservation of relief leads to debate regarding whether the landscape is in steady state (Baldwin et al., 2003;Matmon et al., 2003; , 2011, 2013). Catchment average erosion rates for the southern Appalachians have been measured from cosmogenic radionuclides at between 0.023 and 0.031 mm yr À1 (Matmon et al., 2003) and a pilot study exploring hillslope erosion rates inside the Coweeta watershed calculated the hillslope erosion rate to be between 0.051 and 0.111 mm yr À1 (Hales et al., 2012). These data indicate that in steeply incised terrain, such as the Coweeta Hydrologic Laboratory, the rates of hillslope erosion are governed by the rate of soil production (Hales et al., 2012).

Oregon coast range, Oregon
The Oregon Coast Range in Oregon, USA (Figure 3(b)) is a humid, forested landscape with regularly spaced and heavily incised valleys (Roering et al., 1999). Soils in the area range in thickness from 0.4 m on hilltops to 2 m in unchanneled valleys , identified as common debris flow source areas (Dietrich and Dunne, 1978;Montgomery et al., 2000;Heimsath et al., 2001). The study area is predominantly underlain by Eocene sedimentary rocks of the Tyee Formation (Baldwin, 1956;Snavely et al., 1964). Uplift rates in the OCR have been measured at between 0.1 and 0.3 mm yr À1 from marine terrace data (Kelsey et al., 1996) which corresponds with erosion rates measured using cosmogenic radionuclides of 0.07 to 0.15 mm yr À1 (Beschta, 1978;Reneau and Dietrich, 1991;Bierman et al., 2001;Heimsath et al., 2001). These data have been taken to suggest that the Oregon Coast Range is approximately in steady state (Reneau and Dietrich, 1991;Roering et al., 1999Roering et al., , 2007Roering et al., , 2010Montgomery, 2001;Roering, 2008).

Gabilan mesa, California
Gabilan Mesa is located within the Central Coast Ranges in California, USA (Figure 3(c)) and is characterized by large linear canyons running NE to SW fed by smaller, parallel tributaries typically orientated perpendicular to the canyons (Dohrenwend, 1978(Dohrenwend, , 1979Perron et al., 2008). The area has a semi-arid Mediterranean climate with most rainfall occurring between October and March (Dohrenwend, 1978) and an oak savannah ecosystem . Soil depths are reported by Roering et al. (2007) to be less than 1 m thick on hilltops and thicker in unchanneled valleys. These soils overlie tectonically undeformed marine and continental sedimentary rocks of the Pancho Rico and Paso Robles Formations (Galehouse, 1967;Durham, 1974;Dohrenwend, 1978). Longterm erosion rates in the Gabilan Mesa have been estimated using surface exposure dating and the scale of valley incision into the Mesa at between 0.14 and 0.74 mm yr À1 . The regularity of the valley spacing indicates that the landscape is in approximate topographic steady state (Perron et al., 2009), an observation that is supported by the uniformity of the hilltop curvature across the landscape .

Northern sierra Nevada, California
The northern Sierra Nevada of California, USA (Figure 3(d)) consists of a combination of pine and oak forest, with chaparral vegetation at lower elevations  underlain by thick soils in a semi-arid climate (Yoo et al., 2011;Hurst et al., 2012). The topography is low relief with unevenly spaced valleys incising into a relict surface and exhibiting smooth drainage divides (Gabet et al., 2015). The local geology is primarily ophiolitic, volcanic and sedimentary units of the Fiddle Creek Complex, in addition to intrusive granitoid bodies (Day and Bickford, 2004). A wide range of erosion rates have been measured in the northern Sierra Nevada, with low erosion rates found on the relict surfaces of 0.02 mm yr À1 and higher erosion rates of 0.25 mm yr À1 in canyons (Riebe et al., 2000;Hurst et al., 2012). However, these rates have been shown to vary through time due to accelerating tectonic uplift rates (Stock et al., 2004). This results in a complex soil mantled landscape with spatially varying erosion rates driven by climatic and tectonic variations (Riebe et al., 2000;Stock et al., 2004;Hurst et al., 2012). Hurst et al. (2012) and Milodowski et al. (2015) demonstrated that although the northern Sierra Nevada exhibits a wide range of erosion rates, due to the long response times of the channels in this landscape, the hillslopes have adjusted to local channel lowering and retain a steady topographic form. This was further demonstrated by Hurst et al. (2012) as this landscape plots on the steady state curve in E * vs R * space, whereas Hurst et al. (2013a) showed that transient landscapes should plot above or below the curve.

Basin average measurements of hillslope length
To ensure techniques are compared against equivalent datasets, the hilltop flow routing method results are averaged over the same second-order basins used for results from slope-area and drainage density analyses. Both drainage density inversion and slope-area plots provide information about the landscape and sediment transport over millennial timescales, as both methods generalize topography through surface fitting and spatial averaging. Aspect angle is utilized in the hilltop flow routing algorithm and the surface fitting algorithm from which aspect is computed provides an effective smoothing of the topography without altering the elevation values needed for subsequent analysis. Therefore we can gain a sensible comparison between the three measures, without the influence of small-scale topographic noise on the results.
In Coweeta the basin median hillslope length from hilltop flow routing is 123 m with a median absolute deviation (MAD) of 50 m (Figure 4(a)), the median absolute deviation is a measure of statistical dispersion that is robust, i.e. performs well for probability distributions that are non-normal (Hampel et al., 2011) such as our measured distributions of hillslope lengths. The large outliers seen in this dataset correspond to the large basins present in Coweeta. The slope-area data exhibit underestimation of hillslope length with a median of 28 m with a MAD of 14 m (Figure 4(b)). Few of the longer hillslopes seen in the hilltop flow routing data are reflected in the slope-area data with a strong clustering around the median value. Hillslope length values from drainage density measurements in the southern Appalachians (Figure 4(c)) reflect the smaller area covered by the Coweeta study site with a sparser distribution than seen in the other landscapes. The median lies at 58 m with a MAD of 33 m.
In the Oregon Coast Range the hilltop flow routing algorithm produces basin average hillslope lengths with a median of 89 m and a MAD of 17 m (Figure 5(a)). Results derived from slopearea measurements for this location show a skewed distribution where the majority of measurements are below 10 m, with a median value of 6 m and a MAD of 1 m (Figure 5(b)). The Oregon Coast Range drainage density derived hillslope length 1045 HOW LONG IS A HILLSLOPE? data is normally distributed about a median of 39 m with a MAD of 10 m (Figure 5(c)).
For Gabilan Mesa the median hillslope length from hilltop flow routing is 103 m with a MAD of 26 m (Figure 6(a)). The slope-area data for Gabilan Mesa show fewer measurements overestimating hillslope length than in the other landscapes with a median value of 112 m with a MAD of 56 m (Figure 6(b)). Similar to the slope-area data, the drainage density derived measurements do not form a continuous distribution, with a lower median hillslope length of 42 m and a MAD of 14 m (Figure 6(c)) In the northern Sierra Nevada the hilltop flow routing data show a long tail extending to 400 m with a median of 104 m, and a MAD of 35 m (Figure 7(a)). Slope-area measurements for this location exhibit a skewed distribution of hillslope lengths, extending towards 2.75 x 10 5 m with a median value of 177 m and a MAD of 173 m (Figure 7(b)). The northern Sierra Nevada drainage density data has a similar median to the Oregon Coast Range at 33 m, but with a larger MAD of 13 m (Figure 7(c)), driven by the presence of a number of outlying basins exceeding a median hillslope length of 200 m.

Individual hillslope measurements
The Oregon Coast Range data has a median value of 98 m (Figure 8(b)) and the least variance of the four datasets, with a (MAD) of 35 m. The Coweeta and Sierra Nevada data have similar right skewed distributions but their medians differ by 100 m from 228 to 127 m, respectively (Figure 8(a), (d)). The Coweeta data shows a much larger MAD of 140 m compared with the Sierra Nevada data at 71 m. The Gabilan Mesa data is less right skewed than the Coweeta and the Sierra Nevada data with a median of 140 m and a MAD of 60 m (Figure 8(c)). The dataset shows several smaller peaks which correspond to the different length scales that can be identified in the Gabilan Mesa data.

Discussion
Comparison of hillslope length measurement techniques Figure 9 presents the median and MAD hillslope length value for each technique across the four test landscapes. The slopearea derived hillslope length values for the Oregon Coast Range are distinctly different to the expected range of measurements for this locality, with most basins being underestimated by at least an order of magnitude. However the underestimation does not appear to be systematic as some of the basins have values outside the range of the expected hillslope lengths, with one basin showing a predicted median hillslope length exceeding 2 × 10 4 m. In both the northern Sierra Nevada and Coweeta slopearea derived datasets the large range within the data is driven by the large variation in basin areas found in the landscape with the largest measurements skewing the median data. As with the Oregon Coast Range the modal values are at the lowest end of the distribution, showing a similar underestimation in hillslope lengths to the Oregon Coast Range data. In Gabilan Mesa the slope-area measurements appear to correlate more strongly with the hillslope length values measured from the landscape using hilltop flow routing, suggesting that the regularity of valley spacing and smoothness of the landscape in this location may facilitate more accurate slope-area measurements.
These data show the difficulty in measuring hillslope length using slope-area plots, and although in Gabilan Mesa the estimated values show good agreement with measured hillslope lengths, this pattern is not repeated in the other three landscapes. This failure of slope-area plots to predict hillslope length is not systematic, with both over-and under-estimates observed. This suggests that in steep, complex terrain with a range of processes acting upon it the slope-area technique fails to accurately capture the hillslope to valley transition point.
A fundamental problem in using slope-area plots to infer hillslope length is that the value identified by the kink in the curve is an area, not a length. This value is converted into a length by dividing the area by the unit contour width, which assumes a uniform shape of the upstream contributing area. As can be observed in nature, zero order basins are rarely this uniform, resulting in the derived length value possibly having little bearing on the actual geometry of the landscape.
In order to generate a slope-area plot following techniques described by Hancock and Evans (2006); Roering et al. (2007) and Tarolli and Dalla Fontana (2009), the raw slope and area data must be binned and smoothed in order to highlight the trend. This introduces free parameters and reduces the repeatability of the technique across different landscapes where different parameters may be needed to extract the inflection point (Figure 2). Although the technique can be shown to work in selected cases, it needs user supervision and modification for each individual basin it is applied to and therefore caution is required if it is to be used at the landscape scale. In this analysis the identification of an inflection point is performed automatically and in each location there are several basins where the inflection point identified far exceeds the possible range of hillslope lengths which can be observed in the landscape.
The technique is also limited by a minimum basin size, as a sufficient number of data points are required in order to extract a signal from the data; such a limitation does not exist in the hilltop flow routing technique, which can extract hillslope lengths for individual hillslopes.
Measurements of hillslope length from drainage density systematically underestimate the value measured through hilltop flow routing in all four landscapes. Each location also has a very low variance when compared with the hilltop flow routing measurements. This highlights the inability of this method to capture the true variability in hillslope lengths which can be observed in a landscape, with the majority of basin measurements falling within a very small range of values, even in locations such as Coweeta which has a range of basin areas that scale over an order of magnitude.
The drainage density technique only requires generation of a channel network, which is routinely generated for most forms . Comparison between the three hillslope length measurement techniques, hilltop flow routing (black), slope-area analysis (blue) and drainage density inversion (red). Error bars indicate ±1 MAD around the median value for each landscape. Black dashed lines separate measurements from different landscapes. This figure is available in colour online at wileyonlinelibrary.com/journal/espl of topographic analysis, although the value will be sensitive to the choice of method for identifying channel heads. However, the drainage destiny approach can only provide a single basin average value whereas measurement techniques using flow paths demonstrate that hillslope length can vary substantially within a basin and can therefore provide a range of values for a single basin, allowing more complex patterns of hillslope length variations to be studied across smaller spatial areas. Another benefit of using drainage density to quantify hillslope length is that it is computationally efficient and unambiguous in its implementation compared with other techniques that have many free parameters.
We observe, however, that using drainage density systematically underestimates hillslope length values when compared with flow path derived measurement techniques. Such systematic underestimation is driven by the technique's inability to capture small-scale local variations in hillslope length. This is a necessary outcome of using a basin average method and one which cannot be overcome; the urge to use smaller and smaller basins to increase the measurement resolution consequently will decrease the accuracy of the measurements as fewer streams will be included in each basin's drainage density value, making the resulting value less robust.

Spatial patterns of hillslope length from hilltop flow routing
The large range observed in both the individual hillslope and basin average hilltop flow routing data for Coweeta and the northern Sierra Nevada are driven by variations in hillslope and basin morphology. These landscapes appear transient, with incision into plateaus creating populations of steep, short hillslopes and long hillslopes running across lower gradient plateaus. The large range of basin areas, varying over an order of magnitude, produces the potential for longer hillslopes, particularly in Coweeta. When contrasting Coweeta and northern Sierra Nevada it is clear that the Coweeta dataset is not as well constrained, due to the smaller size of the area covered by the dataset, limiting the number of basins that can be studied.
The more uniform valley spacing of the Gabilan Mesa and Oregon Coast Range Perron et al., 2009) results in a less skewed distribution of hillslope length measurements when contrasted with the other two landscapes. This is particularly apparent in Gabilan Mesa where the more rounded nature of the landscape, with less sharp transitions between hillslope and channel regimes, minimize variation in hillslope length. In the Oregon Coast Range the majority of basin average hillslope length data are within 1 MAD of the median, emphasizing the uniformity in valley spacing and basin area in this landscape. However, with these uniform landscapes a wide range of individual values are observed, highlighting the utility of the hilltop flow routing algorithm to capture topographic information at the hillslope scale.
Comparing relief with predictions from linear and nonlinear flux laws The data from Coweeta is sparse (Figure 10(a)) with only 82 traces identified as planar across the landscape, and all of these traces occurring in lower relief sections of the landscape, yet still follows a broadly linear trend. The reduction in planar traces in this landscape may be connected to the general morphology of the Coweeta catchment, which has less regularly spaced valleys than other similar sites such as the Oregon Coast Range, presenting more opportunity for longer traces to diverge from planarity. This is exacerbated by the high incidence of landsliding reported in the study area (Hales et al., 2009) which leads to increased topographic roughening.
The Oregon Coast Range data are tightly grouped along a linear trend (Figure 10(b)), showing a broad range of L H and R values which correspond to observations of the range of valley spacing and relief made in this location (Roering et al., 1999). These data highlight that the planar trace filtering method is not biasing the sampling towards certain spatial locations or topographic patterns and that the trends of L H and relief reflect the study area as a whole.
The Gabilan Mesa is a smaller geographic area than the other three sites and correspondingly has fewer data points than the Oregon and Sierra Nevada data. However, as with the Oregon Figure 10. Scatter plots of the relationship between hillslope length and relief for each landscape. Each data point corresponds to a single hillslope trace. Red lines reflect Equation 8 fitted to the data using parameters in Supplementary Table 4 The critical slope for this location is computed using a range of erosion rates, producing a maximum and minimum critical gradient for this landscape which differ by 0.003. As these two best fit lines are visually indistinguishable, only the line generated using the upper erosion rate bound is displayed. This figure is available in colour online at wileyonlinelibrary.com/journal/espl 1049 HOW LONG IS A HILLSLOPE?
Coast Range data, a broad range of values have been sampled, highlighting the ability of this technique to operate across a range of landscape morphologies, and again the data shows a linear trend. However, at longer L H , R values stop increasing, reaching an apparent limit at approximately 100 m (Figure 10(c)) suggesting that in the highly diffusive landscape relief greater than this threshold cannot be sustained.
In the northern Sierra Nevada (Figure 10(d)) the data follows a similar pattern to the Oregon Coast Range, although there is more scatter and a larger range of values. This increased scatter is indicative of the wider range of topography in the northern Sierra Nevada, as this data captures hillslope traces from the rapidly incising canyons and from the more slowly eroding plateaus. Being able to identify a linear trend in such variable data, with no need for special processing highlights the robustness and fundamental character of the L H -R relationship and its ability to predict nonlinear sediment flux.
In both the Oregon Coast Range and the northern Sierra Nevada the data at low R and low L H values show much less variance than the remainder of the data for each site. These data points correspond to hilltops located close to the outlet of catchments, which produce very short, planar traces with little variation in the measurements. Such a lack of variation has no impact on the overall patterns in the data and as such there is no requirement to exclude these points from further analysis.
This independent evidence for nonlinear flux is not dependent on collecting erosion rate data or an erosion gradient in a landscape, as is required when analyzing E * R * data to constrain sediment flux laws Hurst et al., 2012). This technique provides clear evidence for the sediment flux law operating on a range of landscapes and highlights that hillslopes are able to record information about the sediment transport process which can be quantified through topographic analysis.
Through the use of published parameters K, ρ r , ρ s and E (see Supplementary Table 4 for parameter values) Equation (8) can be fitted to each dataset to test how well this model explains real topography. The best fit of Equation (8) to the Coweeta dataset uses S c = 0.57 (R 2 = 0.82). This predicted critical gradient is low and suggests that processes other than the creep-like sediment transport described by Equation (1) may be operating in Coweeta that limit relief of hillslope traces. We can say with confidence, however, due to this strong relief limit, that a linear sediment flux law such as that described by Equation (3) is not appropriate in this landscape.
In the Oregon Coast Range the predicted S c is 0.79 with R 2 = 0.96, this gradient falls below the value reported by Roering et al. (2007) of 1.2 ± 0.1, however, it is within published ranges of critical slope values for similar landscapes (Hurst et al., 2012). This predicted critical gradient is lower than a large proportion of the hillslopes in the Oregon Coast Range, and this discrepancy may indicate that Equation (8) does not fully capture the range of sediment transport processes occurring in the Oregon Coast Range.
A similar critical gradient of 0.7245 ± 0.0015 (R 2 = 0.94) is predicted for the northern Sierra Nevada, with the small error bounds generated by performing the fit using the erosion rate from the rapidly eroding canyon and from the more slowly eroding plateau, which is similar to the predicted value of 0.8 from Hurst et al. (2012) which was estimated using E * R * data. The estimate of Hurst et al. (2012) requires averaging relief and ridgetop curvature along ridgetop segments, which are then themselves binned and averaged. In addition the E * R * technique requires a range of erosion rates to constrain the value. The predicted critical gradient for the Gabilan Mesa is equal to the lower bounds of the previously published value of 1.2 ± 0.4 (Roering et al., 1999) at 0.8 with R 2 = 0.88. However, due to the highly diffusive nature of Gabilan Mesa the model does not explain the L H -R relationship as well as in other landscapes.
However, the published values of erosion rate and diffusivity are less well constrained for Gabilan Mesa than in the other landscapes . If the diffusivity value is reduced toward the lower bound of the published range, or the erosion rate is increased toward the upper bound, this fit considerably improves, suggesting that these values may have been over-or under-estimated, respectively.
A limitation of estimating the critical gradient for a landscape using a best fit relationship is that approximately 50% of hillslopes in the landscape will have a gradient in excess of this best fit value. This can explain the discrepancy between critical gradients predicted in this study and those of Roering et al. (1999) which were generated by estimating the erosion rate for planar hillslopes from their morphology, and identifying the values of K and S c which minimized the deviation from the landscape average erosion rate.
A landscape can be expected to exhibit a probability density of critical gradients, controlled by local scale variations in soil, vegetation and topographic properties (Hales et al., 2009). The best fit S c value generated using LH-R data is an averaged S c value, representative of an average hillslope in the landscape, whereas the technique of Roering et al. (1999) may better constrain an upper bound for a landscape critical gradient.

Conclusions
We implement a technique to measure the length of individual hillslopes from digital topography, modelling flow as a point source travelling through DEM cells from hilltop to channel. Through the implementation of this technique two previous methods used to generate basin average hillslope lengths are evaluated. Inverting drainage density produces normally distributed hillslope lengths which systematically underestimate the true lengths of hillslopes found in a landscape. Analyzing slope-area plots for their inflection point in order to extract a characteristic hillslope length is found to be very sensitive to catchment morphology and user defined parameters, producing large amounts of uncertainty in their accuracy.
Hilltop flow routing produces robust results across a range of landscape morphologies, coping with erosional gradients, vegetation changes and differing tectonic regimes. The technique allows for exploration of the variability of hillslope length within a landscape for the first time. By using this flow path method the relationship between relief and hillslope length is explored, providing evidence of nonlinear sediment transport processes in four test landscapes. It is also shown that the critical gradient used in the nonlinear flux model can be constrained. In three of the four test landscapes the predicted critical gradient is shown to be lower than previously published, suggesting that this value may represent the average of a population of critical gradients across a landscape, highlighting the heterogeneity of landscape properties which can be found in steady state landscapes.
Acknowledgements-The topographic data used in this paper is freely available from http://www.opentopography.org. All the code used in this analysis is open source and can be downloaded from http:// csdms.colorado.edu/wiki/Model:Hilltop_flow_routing. This work was supported by NERC grant NE/J009970/1 and US Army Research Office contract number W911NF-13-1-0478. This paper is published with the permission of the Executive Director of the British Geological Survey and was supported in part by the Climate and Landscape Change research programme at the BGS. We thank David Milodowski, Fiona Clubb, T.C. Hales and Robert Parker for discussions and advice which has shaped the final form of this manuscript. The Associate Editor, Joshua Roering, and two anonymous reviewers provided extensive comments and insights which have allowed us to considerably improve this manuscript.