The influence of Holocene vegetation changes on topography and erosion rates: A case study at Walnut Gulch Experimental Watershed, Arizona

. Quantifying how landscapes have responded and will respond to vegetation changes is an essential goal of geomorphology. The Walnut Gulch Experimental Watershed offers a unique opportunity to quantify the impact of vegetation 10 changes on landscape evolution over geologic time scales. The Walnut Gulch Experimental Watershed (WGEW) is dominated by grasslands at high elevations and shrublands at low elevations. Paleovegetation data suggest that portions of WGEW higher than approximately 1430 m a.s.l. have been grasslands and/or woodlands throughout the late Quaternary, while elevations lower than 1430 m a.s.l. changed from a grassland/woodland to a shrubland c. 2-4 ka. Elevations below 1430 m a.s.l. have decadal time-scale erosion rates approximately ten times higher, drainage densities approximately three 15 times higher, and hillslope-scale relief approximately three times lower than elevations above 1430 m. We leverage the abundant geomorphic data collected at WGEW over the past several decades to calibrate a mathematical model that predicts the equilibrium drainage density in shrublands and grasslands/woodlands at WGEW. We use this model to test the hypothesis that the difference in drainage density between the shrublands and grassland/woodlands at WGEW is partly the result of a late Holocene vegetation change in the lower elevations of WGEW, using the upper elevations as a control. Model 20 predictions for the increase in drainage density associated with the shift from grasslands/woodlands to shrublands are consistent with measured values. Using modern erosion rates and the magnitude of relief reduction associated with the transition from grasslands/woodlands to shrublands, we estimate the timing of the grassland-to-shrubland


Introduction 2 Problem statement
Understanding how climate change controls landscape evolution is a central problem in geomorphology.Climate changes are multifaceted, with changes in temperature (mean and variability), precipitation (mean and variability), and vegetation cover (type and density) often occurring simultaneously.
The multifaceted nature of climatic changes can make it difficult to identify which aspects of climate change are most important in driving landscape modification in specific cases.However, given the accelerated climatic changes expected to occur in the coming decades, understanding how landforms are likely to respond to specific climatic drivers, acting alone or in concert, is critically important to society (e.g., Pelletier et al., 2015).
Published by Copernicus Publications on behalf of the European Geosciences Union.
In the southwestern USA the existence of an extensive, regionally correlative fan and valley floor depositional unit (the Q3a unit of Bull, 1991) suggests that the Pleistoceneto-Holocene transition was a major perturbation in what are now semi-arid shrubland-dominated landscapes but which were predominantly pinyon-juniper woodlands at the Last Glacial Maximum (LGM).Where rates of aggradation in modern shrubland drainage basins have been measured, the Pleistocene-to-Holocene transition was associated with more than an order-of-magnitude increase above either LGM or mid-to-late Holocene rates (e.g., Fig. 3.11 of Weldon, 1986).The retreat of grasslands and woodlands to higher elevations and their replacement by shrublands likely exposed elevations of the southwestern USA between approximately 800 m and at least 1400 m a.s.l. to significant increases in percent bare ground, thus modifying the rainfall-runoff partitioning of hillslopes and their resistance to fluvial/slope-wash erosion.In especially arid areas of the southwestern USA such as the central Mojave Desert, the range of elevations affected by late Quaternary conversions of grasslands/woodlands to shrublands extends to elevations as high as 1800 m a.s.l.(Pelletier, 2014).
Given the strong correlation between percent bare ground and drainage density in the southwestern USA (Melton, 1957), it has been hypothesized that modern shrublands sufficiently high in elevation to have been grasslands or woodlands at LGM underwent large increases in drainage density during the Pleistocene-to-Holocene transition.Such an expansion of the fluvial drainage network could have mobilized hillslope deposits stored as colluvium during the last glacial epoch, mobilizing large volumes of sediment into the fluvial system during the transition to the present interglacial (Bull, 1991;Pelletier, 2014).This hypothesis is consistent with measured ages of the onset of aggradation in valley floor and alluvial fan depositional zones in the central Mojave Desert, in which aggradation occurs earliest (ca.15 ka) in depositional zones fed by source regions with relatively low elevations in the 800-1800 m a.s.l.range and progressively later in areas fed by eroding regions at higher elevations (Pelletier, 2014).Alternative explanations for the punctuated nature of late Quaternary aggradation in the southwestern USA invoke changes in the frequency and/or magnitude of floods (Antinao and McDonald, 2013b) and argue that hillslope vegetation changes played a limited role (Antinao and McDonald, 2013a).To date, tests of the paleovegetation change hypothesis in the southwestern USA have been limited to studies of the timing of deposition on valley-floor channels and alluvial fans.Erosion of the source drainage basins themselves has been relatively understudied.
The Walnut Gulch Experimental Watershed (WGEW) provides an excellent opportunity to test the paleovegetation change hypothesis in a drainage basin that has been extensively monitored for decades.The western, low-elevation portion of WGEW is currently a shrubland while the eastern portion is predominantly a grassland (a small area of woodland occupies the highest elevations).Paleovegetation data, however, indicate that all of WGEW was a grassland or woodland until just a few thousand years ago.Scientists working at WGEW have measured rates of sediment export from watersheds using sediment samplers (e.g., Nichols et al., 2008) and rates of sediment redistribution within watersheds using anthropogenic radionuclides (e.g., 137 Cs) and rare-earth element tracers (e.g., Nearing et al., 2005;Polyakov et al., 2009).These data, in addition to the results of analyses of airborne lidar data presented here, make it possible to calibrate every parameter of a mathematical model that predicts the equilibrium drainage density of the landscape under different dominant vegetation types.In this paper we use this mathematical model to test the hypothesis that grassland/woodland-to-shrubland vegetation changes in the lower elevations of WGEW drove large increases in drainage density and erosion rates and a decrease in hillslope-scale relief.

Geology and soils
The WGEW is part of the USA Department of Agriculture (USDA) Agricultural Research Service's (ARS) Southwest Watershed Research Center (SWRC).WGEW is located at the boundary between the Chihuahuan and Sonoran deserts and elevations of between approximately 1300 and 1600 m a.s.l.The approximately 150 km 2 watershed has a planar geometry at large spatial scales, dipping to the WSW with a slope of approximately 1.5 %, i.e., approximately 230 m of elevation change over a distance of 15 km.
The bedrock geology of WGEW includes sedimentary, plutonic, and volcanic rocks of Precambrian to late Cenozoic age (Fig. 1).Due to the complex nature of the rock types exposed in the southern portion of the watershed, we focus this study on the northern portion of the watershed, which is dominated by the Gleeson Road Conglomerate (GRC).
The GRC was derived primarily from erosion of the Dragoon Mountains to the east of WGEW and is estimated to be Plio-Pleistocene in age by Osterkamp (2008), who noted that "The upper part of the Gleeson Road Conglomerate is probably equivalent both stratigraphically and in age to the Plio-Pleistocene upper basin fill of Brown et al. (1966).To the west and northwest, along the axis of the San Pedro River Valley, the upper part of the Gleeson Road Conglomerate grades into fine-grained fluviatile and lacustrine beds of the Plio-Pleistocene St. David Formation (Gray, 1965;Melton, 1965)."The GRC dips gently (5-8 • ) to the north and northwest (Gilluly, 1956, plate 5), a fact which may contribute to the generally longer and more gently sloping south-facing hillslopes relative to north-facing hillslopes in the watershed.The tilted strata of the GRC were beveled to a gently sloping topographic surface and incised into during Quaternary time.Incision of the GRC was driven by uplift/tilting of the fan and/or by incision of the San Pedro River, triggering headward erosion of Walnut Gulch and its tributaries.
The northern portion of WGEW is composed of the whetstone pediment of Bryan (1926) and is divided into two parts: the Dissected Whetstone Pediment (DWP) at elevations below approximately 1430 m a.s.l. and the Upper Whetstone Pediment (UWP) at higher elevations (Fig. 1).The DWP is distinguished from the UWP by its lower relief and less well-developed soils.Differences between the DWP and UWP have been attributed primarily to headward extension of tributaries resulting from late Quaternary incision of the San Pedro River (Cooley, 1968) as well as renewed river and tributary incision following livestock grazing beginning ca.1880 (Renard et al., 1993;Osterkamp, 2008).However, the boundary between the DWP and UWP coincides with the transition from the modern shrubland to the grassland (Fig. 2).As  (Breckenfeld, 1995).The boundary between these two soil groups coincides with the boundary between the DWP and UWP and the transition from the modern shrubland to the grassland.

Climate and vegetation
Mean annual temperature at Tombstone (located in the western portion of WGEW at an elevation of 1384 m a.s.l.) is 17.6 • C and mean annual precipitation (MAP) is approximately 300 mm.There is both a winter and summer rainy season, but approximately 70 % of rainfall occurs in the summer months and the greater intensity of summer storms means that runoff occurs almost exclusively during the summer season of July through September.Mean annual precipitation is approximately 10 % higher and rainfall erosivity (which is based on rainfall intensity at 30 min duration) is approximately 20 % higher at the highest elevations of the watershed relative to the lowest elevations (Nearing et al., 2015).
Modern vegetation cover in the USA-Mexico borderland region closely follows elevation, with desert scrub occurring below approximately 1500 m, grasslands from 1400 to 1700 m, lower encinal (encina is Spanish for "oak") from 1700 to 2600 m, upper encinal from 1900 to 2600 m, and for-est from 2200 to 2600 m (Wagner, 1977).For each zone, the highest elevations are reached on dry, south-, and west-facing slopes and the lowest elevations on north-facing slopes and valley bottoms.
Paleovegetation records from the borderland region suggest that the western portions of WGEW have transitioned from a grasslands/woodland to a shrubland over the past few thousand years, while the eastern half of the watershed has been a grassland/woodland for at least the past 40 000 years and likely since the penultimate interglacial period.Low stalagmite 18 O values at the LGM in the Cave of the Bells paleoclimate record indicate that conditions were much wetter and cooler during LGM (Wagner et al., 2010), in agreement with paleovegetation studies (Betancourt et al., 1990(Betancourt et al., , 2001;;Anderson, 1993;Arundel, 2002;Holmgren et al., 2005;Holmgren, 2005).Packrat midden records indicate the presence of grasslands and/or pinyon-juniper woodlands at LGM in what is currently Chihuahuan desert scrub at elevations of 1200-1400 m a.s.l.(Betancourt et al., 2001).Holmgren (2005) documented the presence of the primary grass species at WGEW (Bouteloua eriopoda) (currently abundant only at elevations above approximately 1430 m a.s.l.) at an elevation of 1287 m ca.4750 14 C years BP.Elements of the modern shrubland, such as Larrea tridentate, appeared as late as 2190 14 C years BP at 1287 m a.s.l.based on macrofossils, but they may have been present as early as 4095 14 C years BP based on pollen.As such, the latest Quaternary transition from grasslands/woodlands to shrublands in WGEW occurred gradually and was completed only within the past few thousand years.

Intensive monitoring sites
WGEW is home to two intensive monitoring sites, one in the shrubland of the DWP and the other in the grassland of the UWP.These sites provide the data necessary, in conjunction with the topographic analyses presented here, to calibrate a mathematical model that predicts the equilibrium drainage density as a function of vegetation cover and to test the hypothesis that differences in landscape morphology and erosion rates between the northwestern and northeastern portions of WGEW are partly the result of a transition from grassland/woodland to shrubland within the past few thousand years in the northwestern portion of the drainage basin that did not occur in the northeastern portion.
Watersheds 102, 103, 104, 105, and 106 are located in shrublands at an elevation of approximately 1370 m a.s.l. in what is referred to as "Lucky Hills," which has been the site of a variety of intensive scientific studies since the 1960s.Cover during the rainy season at Lucky Hills is approximately 25 % bare soil, 25 % canopy, and 50 % erosion pavement (rocks) (Nearing et al., 2007).Dominant vegetation includes Creosote (Larrea tridentata, shrub) and Whitethorn (Acacia constricta, shrub), with lesser populations of Desert Zinnia (Zinnia acerosa, shrub), Tarbush (Flourensia cernua, shrub), and sparse Black Grama (Bouteloua eriopoda, grass) (King et al., 2008).The matrix material of surface layer is composed of 60 % sand, 25 % silt, and 15 % clay (Nearing et al., 2007).Sediment from the watershed is monitored with a supercritical flume with an automatic traversing slot sampler (Simanton et al., 1993).The shrubland sites have a range of slope aspects and exhibit sediment yields approximately 30 times higher than the grassland sites across all slope aspects.
Watershed 112 is located in the Kendall grassland site at WGEW, approximately 10 km east of Lucky Hills and at an elevation of approximately 1525 m.Kendall drains to the west and has roughly equal areas of north-and southfacing hillslopes.The site is predominantly covered by grass and forbs with some shrubs and succulents with a combined canopy cover of approximately 35 %.Ground cover during the rainy season has been measured at 28 % rock, 42 % litter, and 14 % plant basal cover (Nearing et al., 2007).Compared to the 25 % bare soil at Lucky Hills, the bare soil exposed at Kendall is negligible (i.e., a few percent or less).Historically, the dominant desert grassland bunchgrasses at the site have been black grama (Bouteloua eriopoda), sideoats grama (B.curtipendula), three-awn (Aristida sp.), cane beardgrass (Bothriochloa barbinodis) (King et al., 2008), and, more recently, Lehmann lovegrass (Eragrostis lehmanniana) (Moran et al., 2008).Nearing et al. (2005) used spatially distributed 137 Cs measurements to quantify fluvial/slope-wash erosion and deposition rates within and from watersheds 103 (Lucky Hills) and 112 (Kendall).Nearing et al. (2005) found that in Lucky Hills the fraction of the drainage basin experiencing erosion was much higher (85 %) than in Kendall (53 %), where erosion and deposition rates were lower and approximately balanced such that the rate of net sediment exported from the Kendall watershed was more than a factor of 10 lower than from the Lucky Hills watershed.These observations are consistent with sediment yields measured from 1995 to 2005 that are more than 10 times higher in drainage basins of similar size at Lucky Hills (231 t km −2 yr −1 in watershed 102, area of 1.46 Ha) than at Kendall (7 t km −2 yr −1 in watershed 112, area of 1.86 Ha) (Nearing et al., 2007).Assuming a soil bulk density of 1500 kg m −3 , these sediment yields correspond to mean erosion rates of approximately 1.5 × 10 −4 m yr −1 in the shrublands of Lucky Hills and 5 × 10 −6 m yr −1 in the grasslands of Kendall (Table 1).Hillslopes are slightly steeper in Kendall, so if anything we would expect erosion rates to be higher at Kendall than at Lucky Hills if slope gradient were the dominant factor in controlling erosion rates.Nearing et al. (2005) interpreted the differences in erosion rates between Lucky Hills and Kendall to be primarily a function of vegetation cover; i.e., "hydrologic response differences as a function of vegetation differences are probably largely responsible for the differences in hillslope erosion rates between the two watersheds.When flows are more concentrated and vegetative cover is less, as on the Lucky Hills site, flow shear stresses and stream power will tend to be greater, resulting in a greater hydrologic potential for erosion.Also important is probably the higher litter cover and plant basal area cover on the grassland site that would have a direct protective effect against erosion."This interpretation is consistent with the conclusions of Abrahams et al. (1995) and Parsons et al. (1996), who emphasized the role of vegetation cover in controlling fluvial/slope-wash erosion rates in their plot-scale studies at WGEW.In this paper we explore the implications of these erosion rate differences for landscape morphology and topographic evolution of WGEW over geologic timescales.

Topographic analysis
In this section we describe the methods used to quantify the similarities and differences in landscape morphology between the modern grassland and shrubland sites, with an eye toward providing the data necessary to calibrate the mathematical model as described in Sect.3.2.In all of the topographic analyses described in this section we used a 1 m pixel −1 digital elevation model (DEM) for WGEW derived from airborne laser swath mapping (Heilman et al., 2008).Prior to the analyses, the DEM was smoothed with an optimal Wiener filter (OWF), following the approach of Pelletier ( 2013), to remove small-scale variability related to errors related to distinguishing ground from vegetation points and other imperfections of the lidar-derived DEM.The smoothing did not significantly alter slope gradients but did significantly reduce anomalous curvature values related to DEM imperfections.The drainage density in the shrubland areas appears to be much higher than in the grassland areas (Fig. 3).To quantify this difference we used the method developed by Pelletier (2013) to identify the network of valley bottoms (i.e., where water flow is localized and fluvial processes are responsible for the majority of erosion) in the vicinity of Lucky Hills and Kendall.In this method, the DEM is filtered using the OWF, the contour curvature is computed at every pixel, and valley heads are identified as the areas closest to the divides where the contour curvature exceeds a userdefined threshold value.In Pelletier (2013) and in this paper a threshold contour curvature of 0.1 m −1 was used for valley head identification (i.e., the method of Pelletier (2013) was used without modification or parameter tuning).Once valley heads are identified, the multiple-flow-direction routing algorithm of Freeman et al. (1991) is used to route a unit of runoff from each valley head to identify the valley bottoms downstream.We mapped the longest distance upslope from every hillslope pixel along steepest-descent flow lines.We then averaged this flow distance value for all pixels adjacent to valley heads.This mean value can be compared to the prediction of a mathematical model that computes the mean distance along flow lines from topographic divides to valley heads as a function of colluvial and fluvial transport coefficients.Hillslope length measured in this way is inversely related to drainage density (Horton, 1945).Its mean value contains information equivalent to drainage density, but it has the advantage that it is a mappable quantity (Tucker et al., 2001).The drainage network analysis was performed on representative 1.6 × 1.6 km areas in the vicinity of Lucky Hills and Kendall to quantify the difference in drainage density between shrubland and grassland areas within WGEW.
Relief was mapped as the difference between the highest elevation upstream from each pixel along flow lines.The results of the drainage network identification procedure described above were used to limit this analysis to the hillslope pixels only.We then computed the mean hillslope relief within 10 m elevation bins from 1320 to 1550 m a.s.l.The resulting graph quantifies how the mean hillslope relief varies with elevation across the shrubland-to-grassland transition.
We also computed the mean value of the along-channel slope gradient and curvature (i.e., the Laplacian) as functions of contributing area.Differences in mean curvature as a function of contributing area provide a quantitative signature of how late Holocene vegetation changes have modified the landscape morphology in the vicinity of the hillslope-tovalley-bottom transition.

Mathematical modeling
In this section we describe the mathematical model used to quantify erosion over geologic timescales and its dependence on landscape morphology at WGEW.The mathematical model is used to predict the equilibrium drainage density, quantified as the mean distance along flow lines from divides to valley bottoms, in both shrublands and grasslands, in order to test the hypothesis that the difference in drainage density observed between the shrublands and grasslands at WGEW can be attributed, in part, to late Holocene vegetation changes in the shrubland portion of the watershed.
Erosion at WGEW over geologic timescales can be approximated as the sum of erosion due to colluvial and fluvial/slope-wash processes.Sediment transport by colluvial processes increases approximately linearly with slope based on short-term monitoring studies (Kirkby, 1967) and cosmogenic radionuclide analyses (McKean et al., 1993) when slopes are uniformly soil-mantled and topographic gradients are modest.Linear slope-dependent transport, combined with conservation of mass, leads to a diffusion equation for topography (Culling, 1960): where E c is the erosion rate by colluvial processes (defined to be positive when the landscape is eroding), D is diffusivity in m 2 yr −1 , and z is elevation in m.Equation (1) assumes that colluvial sediment flux is equal to the product of a coefficient (i.e., the diffusivity D) and the local slope gradient.
This assumption is reasonable for WGEW, where hillslopes are uniformly soil mantled and the mean hillslope gradient is 7 %.We assume that fluvial/slope-wash processes in WGEW can be approximated as transport limited.We use the term fluvial/slope-wash to refer to all sediment transport by flowing water, wherever it occurs along the continuum from hillslopes (i.e., as sheet and rill flow) to channels (i.e., fluvial erosion).A transport-limited condition applies to landscapes in which most of the sediment is transported as bedmaterial load and sediment is readily deposited when the shear stress by flowing water declines with increasing distance along flow pathways.In the alternative detachmentlimited end-member model of fluvial/slope-wash erosion, the shear stress required to detach sediment is much larger than the shear stress required to transport it; hence sediment redeposition is rare or nonexistent once detachment/entrainment occurs.Pelletier (2012) addressed the geomorphic conditions under which transport-limited vs. detachment-limited conditions are likely to occur, taking into account data for the relative proportion of sediment transported as bed-material load vs. wash load, among other factors, using WGEW as a case study.He concluded that among these two end-member models, WGEW is most accurately considered to be transport limited.
The assumption of transport-limited conditions implies that fluvial/slope-wash erosion, E f , is equal to the divergence of the fluvial/slope-wash volumetric unit sediment flux, q s (the volumetric sediment flux per unit width of water flow): where ε 0 is the grain packing density (assumed here to be 0.55, e.g., equivalent to a bulk density of 1500 kg m −3 for a grain density of 2700 kg m −3 ).The divergence of the fluvial/slope-wash sediment flux can, within the valley network, be approximated by the derivative of the volumetric unit sediment flux in the along-valley direction, q s , with respect to the distance downstream along flow lines, x.
Adopting this approximation and summing the colluvial and fluvial/slope-wash components, the total erosion rate at any point on the landscape is thus given by In this paper we use Eq.(3) to predict the equilibrium drainage density, quantified as the mean distance from divides to valley bottoms, associated with grasslands and shrublands at WGEW.Specifically, we substitute an empirical power-law relationship for the mean distance, x, along flow lines from divides to valley bottoms vs. contributing area, A, into Eq.( 3) and solve for the value of x where the fluvial erosion rate exceeds the colluvial deposition rate by an amount equal to E.
Next, we further parameterize the three terms in Eq. (3) in terms of measured proxies for erosion rate (e.g., decadalscale sediment fluxes) and topographic parameters.The erosion rate, E, is constrained using the ratio of the volumetric sediment flux, Q s , reported by Nearing et al. (2007) and contributing area, A: The fluvial erosion term can be written in terms of Q s using where w is the valley-bottom width.We approximate the mean rate of colluvial deposition at valley heads by where S h is the mean slope gradient of toe slopes as they intersect the valley bottom at valley heads (Fig. 5).Equation ( 6) derives from the mass balance of a square segment of the valley bottom in the vicinity of the valley head (Fig. 5).According to the assumption that soil-mantled hillslopes evolve diffusively, valley bottoms in the vicinity of the valley head receive a unit sediment flux equal to DS h by colluvial transport processes from each of three adjacent hillslope segments, i.e., the two hillslopes in the cross-sectional direction and the one entering the valley head from upslope.In the cross-sectional profile, the difference in sediment flux across the profile is equal to DS h and that difference occurs over a distance of w.
As such, the colluvial deposition rate, computed in a manner that is independent of DEM or pixel resolution, includes the valley width in the denominator (Pelletier, 2010a).Similarly, the divergence in the along-valley direction is approximately DS h /w.The factor of 2 does not appear in the along-valley derivative because colluvial sediment flux enters the valley bottom segment only from the upslope direction.Colluvial sediment flux leaving the valley head is assumed to be negligible because the slopes of the valley bottom immediately downslope from the valley head (approx.5 %) tends to be much smaller than that of the hillslopes entering it from upslope (approx.15-20 %) (data presented in Sect.3.1).Next, we introduce three power-law relationships that relate volumetric sediment flux to contributing area, channel width to contributing area, and contributing area to distance along flow lines from topographic divides.Sediment flux is a power-law function of contributing area at WGEW: Equation ( 7) will be calibrated in Sect.3.2 based on data from Nearing et al. (2007).The coefficient k Q s in Eq. ( 7) is a sediment transport efficiency parameter that is a function of runoff, sediment texture, vegetation cover, and potentially other factors.It takes on different values in the shrubland and the grassland (i.e., k Q s s and k Q s g ), reflecting the fact that sediment yield is a function of vegetation cover.Its value in shrublands, along with that of p, is obtained by a leastsquares regression of Eq. ( 7) to Q s values for the five shrubland drainage basins of Lucky Hills reported by Nearing et al. (2007).Its value for the grasslands is constrained by assuming that the same value of p applies to both shrublands and grasslands and using the single Q s data point available for grasslands (i.e., for watershed 112) to constrain k Q s using Eq. ( 7).Bedload sediment flux is a function of both contributing area and along-channel slope.However, along-channel slopes (Table 1) vary only weakly with contributing area at the monitoring sites; i.e., along-channel slope varies by only a factor of 2 as contributing area varies by more than a factor of 20.As such, including slope does not significantly modify or improve the regression of measured data against its controlling factors (see Sect. 3.1).Hence, we did not include it explicitly in Eq. ( 7).Miller (1995) measured 222 cross-sectional channel profiles in the field at WGEW and used those data to calibrate a power-law relationship between channel width and contributing area: where k w = 0.023 m 0.32 and l = 0.34.Contributing area has a power-law relationship with the mean distance, x, from divides along flow lines: In Table 2 we report the values of k Q s s and k Q s g , p, k w , l, k A , and c for WGEW.Substituting Eqs. ( 4)-( 9) into Eq.( 3) yields Equation ( 10) is applied using vegetation-specific values for k Q and S h to solve for the mean distance from divides to valley heads in shrublands and grasslands, i.e., x s and x g , respectively.That is, x s is obtained by solving and x g is obtained by solving Equation ( 3) is generally applicable within the fluvial network.Once the colluvial deposition rate is approximated using Eq. ( 6) (which makes use of S h , the mean gradient of the hillslope toe slopes at valley heads), subsequent equations, including Eqs. ( 10)-( 12), apply only to valley heads.Equations (10)-( 12) are a mathematical representation of the conceptual model, first proposed by Tarboton et al. (1992), that erosion at valley heads is a competition between transportlimited fluvial erosion and colluvial deposition.That is, fluvial erosion rates must exceed colluvial deposition rates by an amount equal to the net erosion rate on the landscape in order to maintain an equilibrium drainage density.
To estimate the time required for low-order channels to grow headward in response to hillslope vegetation changes, we modeled the channels as a diffusive system (e.g., Begin, 1988).In diffusive systems, the timescale, t r , required for a response over a length scale, L, can be estimated using where D is a diffusivity.The value of D in the diffusive model for alluvial channels is given by the ratio of the unit volumetric sediment flux to the along-channel slope:

Topographic analyses
Figure 3 illustrates the results of the drainage network identification for 1.6 km × 1.6 km examples of the landscape in the vicinity of the Lucky Hills and Kendall sites.The mean distance along flow lines from divides to valley bottoms is 18 m in the shrubland area and 50 m in grassland area, using the Pelletier (2013) algorithm.Figures 4a and b illustrate that hillslopes in the shrubland area are more finely dissected with rills and gullies than in the grassland area.
As part of this analysis we also measured the mean slope gradient of hillslopes immediately adjacent to valley heads, i.e., S h in Eq. ( 10).We obtained S hs = 0.17 ± 0.03 m m −1 in shrublands and S hs = 0.19 ± 0.04 m m −1 in grasslands (uncertainty is the standard deviation).In first-order valleys, along-channels slopes are approximately 0.05 ± 0.03 m m −1 .Mean hillslope relief increases substantially across the shrubland-to-grassland transition (Fig. 6).Between elevations of approximately 1320 and 1430 m a.s.l., mean relief is uniformly low (approx.0.5-1 m).Above elevations of approximately 1450 m, hillslope relief increases abruptly and continues to increase with increasing elevation.
Contributing area follows a piecewise power-law function of distance along flow lines from topographic divides (Fig. 7  the valley network.Above contributing areas of approximately 50 m 2 (or, equivalently, distances from the divide equal to approximately 15 m), contributing area increases as the 2.5 power of distance from the divide for both grassland and shrubland areas, i.e., k A = 0.3 m −0.5 and c = 2.5 in Eq. ( 9).Below contributing areas of 50 m 2 , k A = 2 m −0.75 and c = 1.25.We used k A = 0.3 m −0.5 and c = 2.5 when solving Eqs. ( 11) and ( 12), since these values are most applicable to points within the valley network (i.e., at valley heads and points downstream).This is a self-consistent approach because the solutions to Eqs. ( 11) and ( 12) are larger than 15 m (Sect.3.2).Plots of mean topographic curvature (i.e., the Laplacian of z) vs. contributing area (Fig. 8) indicate that mean curvatures are nearly identical at small and large contributing areas but differ substantially within the range of contributing areas from ∼ 10 to 300 m 2 .Plots of mean topographic curvature as a function of contributing area (distance from divide also shown along x axis using data in Fig. 7).Error bars represent the standard error of the mean for each bin.Topographic curvatures are similar in shrublands and grasslands at small and large contributing areas, with a minimum of 0.03 m −1 near divides.Within a range of contributing areas from ∼ 10 to ∼ 300 m 2 the data show significant differences in mean curvature between shrublands and grasslands.

Mathematical modeling
In this section we use the results of Sect.3.1.,together with analyses of the sediment yield reported by Nearing et al. (2007), to constrain the terms in Eqs. ( 11) and ( 12) in order to solve for the mean distance from divides to valley heads in shrublands and grasslands.In order to constrain the absolute values of k Q s s and p, we performed a least-squares minimization of Eq. ( 11) to the decadal-scale sediment yields reported by Nearing et al. (2007) for watersheds 102-106 (Lucky Hills) (Table 1).This regression yields k Q s s = 2 × 10 −6 ± m 1.56 yr −1 (with a range of values including 1 standard error from 2 × 10 −7 to 2 × 10 −5 ), p = 1.44 ± 0.2, and R 2 = 0.93 (Fig. 9).Assuming that the value of p derived from the shrubland watersheds also applies to the grassland watershed, the value of k Q s for the grassland is estimated to be k Q s g = 6 × 10 −8 m 1.56 yr −1 .For D we adopt the value of 1 × 10 −3 m 2 yr −1 commonly inferred from scarp degradation studies in the southwestern USA (e.g., Hanks, 2000; Table 2 cites D values for the basin and range of the western USA of between 6.4 × 10 −4 and 2 × 10 −3 m 2 yr −1 based on eight published studies).The full list of model parameters and their values is provided in Table 2.
We used Eqs.( 11) and ( 12) to predict the mean distance along flow lines from divides to valley bottoms in shrublands and grasslands.The predicted values are x s = 16 m and x g = 66 m (Table 3).In the topographic analysis shown in Fig. 3 7), from which the values of k Q s s and p were constrained.
lands.As such, the model predicts mean distances from divides to valley bottoms within 10 % of measured values.
Figure 10 plots the magnitude of the three terms for a range of possible mean distances from divide to valley bottom.The fluvial erosion rate is ∼ 10 −2 m yr −1 in shrublands and ∼ 10 −3 m yr −1 in grasslands, increasing with distance downstream, reflecting the nonlinear relationship between sediment flux and drainage area (Fig. 9).The colluvial deposition rate is ∼ 10 −3 m yr −1 for both shrublands and grasslands and decreases modestly with increasing distance from the divide as a result of the increase in channel width with increasing contributing area (Eq.8). Figure 10 demonstrates that the fluvial erosion rate must be quite large (approximately 2 orders of magnitude larger than the net erosion rate in this case) in order to counteract the effects of colluvial deposition and thus maintain a valley head.
We used the mean curvature vs. contributing area data plotted in Fig. 8 to reconstruct an average longitudinal profile from divides to valley bottoms in shrublands and grasslands in order to infer the approximate relief reduction associated with the late Holocene shift from grasslands/woodlands to shrublands in WGEW. Figure 11 illustrates the results of this integration.Integrating the curvature vs. contributing area data in Fig. 8 twice leads to two constants of integration: one is constrained by the requirement that the slope along flow pathways at divides is zero and the other by using a constant reference elevation at a contributing area of approximately 300 m 2 .We chose A ≈ 300 m 2 as the location to enforce the reference elevation because this is the contributing area below which the mean curvature begins to deviate significantly between the grassland and shrubland areas.The difference in elevation between the two profiles plotted in Fig. 11 provides an estimate of the minimum erosion or  relief reduction associated with the shift from grassland to shrubland in the lower elevations of WGEW.We consider the results of Fig. 11 to be a minimum estimate because systematic differences in mean slope between the grassland and shrubland across a wide range of scales are not reflected (or not fully reflected) in curvature or any reconstruction of the longitudinal profile based on integrating the curvature.The results in Fig. 11 suggest that divides have lowered by a minimum of approximately 0.3 m as a result of Holocene vegetation changes.Given hillslope-scale erosion rates of approximately ∼ 10 −4 m yr −1 in shrublands and the much smaller erosion rate of ∼ 10 −5 m yr −1 in grasslands, it would take approximately 3 kyr following a transition from grasslands/woodlands to shrublands to erode the landscape by an amount 0.3 m greater in shrublands compared to grasslands/woodlands.This estimate is comparable to the age of the grassland/woodland-to-shrubland transition in the region at an elevation of 1287 inferred from paleovegetation studies in the region, i.e., 2-4 14 C kyr BP (Holmgren, 2005).
In order to test the hypothesis that 2-4 kyr is sufficient time for drainage density to have fully responded to the recent grassland/woodland-to-shrubland transition, we used Eqs.( 13) and ( 14) to estimate the response time of loworder channels.We used L = 50 m to represent the approximate distance that the valley heads migrated upslope in response to late Holocene vegetation changes.We used data associated with the outlets of watersheds 105 and 106 as representative of the low-order channels in the shrublanddominated portions of WGEW.Using Eq. ( 14), the D value for these channels is approximately 3 × 10 3 m 2 kyr −1 based on the ∼ 0.1 m 3 yr −1 sediment flux of the 105 and 106 watersheds (Fig. 10), a width of approximately 1 m, and slope of approximately 0.03.Substituting these values for L and D into Eq.( 13) yields an estimate for t r equal to 0.83 kyr, or ∼ 1 kyr in keeping with the approximate nature of this calculation.This timescale is significantly shorter than the age of the transition (2-4 kyr), suggesting that sufficient time has occurred for drainage density to have fully responded to the changes in vegetation cover in the lower elevations of WGEW.

Uncertainty in parameter values and their impact on the model results
Equations ( 11) and ( 12) predict mean distances from divides to valley bottoms that are broadly similar to measured values.Several factors limit the accuracy of the comparison between measured and predicted values.First, we relied upon a regional value for the diffusivity D because we do not have a reliable means of calibrating this value locally.Second, decadal-scale erosion rates computed from sediment samplers may differ from long-term erosion rates.Nearing et al. (2007) found that for the six watersheds considered here, half of the sediment yield measured from 1995 to 2005 was derived from 6 to 10 events and that the largest events contributed between 9 and 11 % of the total yield.Thus, while transport is highly episodic in WGEW, the most effective flood events have a sufficient number of recurrences to provide an estimate of the yield that does not depend sensitively on the timescale.This conclusion is consistent with the fact that erosion rates inferred from sediment sampling from 1995 to 2005 are similar to erosion rates measured over the postbomb period using 137 Cs (Nearing et al., 2005).That said, the sediment yields reported by Nearing et al. (2007) may not include the most extreme drought conditions or other disturbances that could cause long-term sediment yields to be larger than those reported in Table 1.
In particular, Kirchner et al. (2001) demonstrated the potential pitfall of applying sediment yields or erosion rates measured at one timescale to another timescale by demonstrating that rates over different timescales can differ by orders of magnitude.We agree that this is a significant potential concern.However, we believe that our decadal-scale sediment yields are an appropriate estimate for millennial-scale sediment yields based on three lines of argument.First, as noted above, published analyses have shown that fluvial sediment transport in WGEW, while highly episodic, is not dominated by just a few extreme events.The effective discharge (i.e., the discharge above which half of the total load is transported) occurs many times within a 30-year record.In a study of sediment transport in the 1995-2005 period, for example, Nearing et al. (2007) addressed this issue as follows: "For six of the seven watersheds, between 6 and 10 events produced 50 % of the total sediment yields over the 11-year period."That is, the effective discharge has a recurrence interval of between approximately 1 and 2 years.Second, sediment yields calculated from 1995 to 2005 by Nearing et al. (2007) closely match yields measured over approximately 50 years using 137 Cs (Nearing et al., 2005).This 50-year timescale includes significant droughts at WGEW; i.e., time periods when vegetation cover would have been anomalously low and erosion rates high.Third, recent work suggests that the order-of-magnitude increase in sediment yields measured by Kirchner et al. (2001) between interannual and millennial timescales may partly reflect the influence of high-severity (i.e., stand-replacing) forest fires in their study area in Idaho.Orem and Pelletier (2016) measured wildfire-affected and non-wildfire-affected erosion rates measured over interannual timescales together with erosion rates measured over millennial and million-year timescales in the forested Valles Caldera, New Mexico.Orem and Pelletier (2016) observed a similar increase in erosion rates to that of Kirchner et al. (2001).They were also able to demonstrate that the increase they measured was due to the episodic effects of highseverity wildfire.In contrast to forested areas, wildfires are of limited size and severity in shrublands.Similarly, grassland fires in Arizona typically result in modest (if any) increases in runoff and erosion rates (e.g., Stone et al., 2003).As such, there is a basis for concluding that the large increase in erosion rates measured by Kirchner et al. (2001) between interannual and millennial timescales is unlikely to occur in WGEW.

Further discussion of the hypothesis of a vegetation-change-driven increase in drainage density in the shrublands of WGEW
We also propose that the difference in mean curvatures between shrublands and grasslands within the range of contributing areas from ∼ 10 to ∼ 300 m 2 partly reflects recent expansion of the drainage network in the shrublands of WGEW.This hypothesis is consistent with the fact that the deviation of curvature values between the two sites begins at a contributing area comparable to the support area (i.e., the contributing area required to form a valley head) in the grasslands and complementary studies (e.g., Yetemen et al., 2010) that have quantified the influence of vegetation cover on curvature-area relationships.As such, we propose that the shrubland areas previously had support areas comparable to the grassland areas and that drainage network expansion has influenced the drainage network and the morphology of the adjacent hillslopes down to spatial scales corresponding to contributing areas of ∼ 10 m 2 .The fact that curvature values are very similar between shrubland and grassland below spatial scales ∼ 10 m 2 is consistent with the hypothesis that hill-slopes in the lower elevations of WGEW have not yet fully adjusted to the increase in drainage density associated with the grassland-to-shrubland transition.However, as the analysis of Eqs. ( 13) and ( 14) reveal, the drainage density itself (which is the primary focus of this paper) has likely had sufficient time to fully adjust to the vegetation changes.Figure 5 demonstrates that mean hillslope relief increases substantially across the shrubland-to-grassland transition.We propose that some of this difference in hillslope relief is a consequence of the difference in fluvial/slope-wash erosion rates, i.e., that higher erosion rates in the lower-elevation shrublands have caused relief reduction in the past few thousand years relative to the higher-elevation grasslands, and that this difference in fluvial/slope-wash erosion rates is the result of a geologically recent increase in drainage density in the shrublands of WGEW.However, it is likely that a portion of the difference in mean hillslope relief across the study site also reflects variable uplift rates, i.e., the fact that the uplift of any piedmont or foothill region tends to increase towards the mountain range.Flexural-isostatic response to erosion (which has been proposed to be an important component of late Cenozoic landscape evolution in southern Arizona; Menges and Pearthree, 1989;Pelletier, 2010b) of the Dragoon Mountains can be expected to have caused eastward tilting of WGEW, i.e., higher uplift rates in the higher elevations of WGEW compared to the lower elevations.Tilting would not explain the abrupt increase in relief at elevations just above 1430 m a.s.l., however, since no faulting occurs in the vicinity of this contour.Therefore, it is likely that some of the difference in hillslope-scale relief across the shrublandto-grassland transition at WGEW is a result of the difference in erosion rates between the shrublands and grasslands.While we can be certain of grassland-to-shrubland shift only during the present interglacial period, the Quaternary period has seen many interglacial periods broadly similar in climate to the current period; hence it is likely that the lower elevations of WGEW have seen grassland-to-shrubland conversions more than once over the past approximately 2 Myr.Each of these episodes could have contributed to relief reduction of the lower elevations of the study site relative to the higher elevations.
Previous studies at WGEW have emphasized the role of base-level lowering and vegetation changes within the past 130 years on the differences in erosion rate between Lucky Hills and Kendall (Nearing et al., 2007).However, recent paleovegetation studies have provided a new perspective.Specifically, Holmgren's (2005) documentation of shrubland species in the region several thousand years before present at an elevation less than 100 m lower than Lucky Hills suggests that the lower elevations of WGEW likely shifted from a grassland/woodland to a shrubland prior to the 1880s.While base-level lowering has steepened hillslopes and channels close to the main-stem channel of Walnut Gulch, hillslopescale relief is clearly larger at Kendall than at Lucky Hills (Fig. 6), indicating that base-level lowering may be a domi-nant factor only for those areas within close proximity to the main channel.The magnitude of the differences in topography (i.e., drainage density and the magnitude of erosion than can be inferred from the change) is difficult to fit into a period as short as 130 years.Given erosion rates measured over the past 60 at Lucky Hills, approximately 2 cm of erosion can be expected to have occurred over the past 130 years.Figure 11 suggests that erosion associated with a recent increase in drainage density is likely at least 10 times this value and thus more consistent with a vegetation change that occurred several thousand years before present.

Implications for our understanding of the controls on drainage density
The model of this paper contributes to our broader understanding of the controls on drainage density and it provides a mathematical model for predicting drainage density that may be useful in other study sites.
Previous studies have demonstrated that drainage density is controlled by relief (e.g., Montgomery and Dietrich, 1992;Tarboton, 1992;Tucker and Bras, 1998), climate (Melton, 1957;Abrahams and Ponczynski, 1984), parent material (e.g., Ray and Fischer, 1960;Day, 1980), and time (e.g., Ruhe, 1952;Dohrenwend et al., 1987).While many studies have demonstrated the importance of individual factors on drainage density, we lack a comprehensive model for drainage density that integrates all of these factors.Equation (10) represents one possible candidate for such a model in soil-mantled, transport-limited landscapes.Time is not included in the model because it is based on an equilibrium mass-balance framework.Nevertheless, Eq. ( 10) predicts the drainage density to which a transient landscape will likely approach over time following a perturbation.Relief enters the model via the erosion rate, E (quantified for the case of WGEW using multiple measurements of Q s /A), and the mean toe slope gradient near valley heads, S h .Climate and vegetation cover enters the model through the parameters D and k Q s .The value of D increases with increasing soil moisture and temperature cycling around 0 • C (which together drive soil creep) and likely also with increasing vegetation cover.Evidence for vegetation control of D values comes from Hanks (2000), who compiled data on D values estimated from morphological analysis of dated shorelines from the Negev in Israel, the semi-arid southwestern USA, and sub-humid to humid areas in California and Michigan.The available data suggest that D values increase from dry to wet climates and/or from areas of low to high vegetation cover: D values are ∼ 0.1-0.3m 2 kyr −1 in arid areas, ∼ 1 m 2 kyr −1 in semi-arid areas, and ∼ 10 m 2 kyr −1 in sub-humid and humid areas.Values of k Q s increase with rainfall and decrease with vegetation cover.In addition, Eq. ( 10) explicitly includes channel width and its scaling with contributing area, factors that, to our knowledge, have not been included in previous mathematical models of drainage density.
Drainage density is most commonly found to be an inverse function of mean annual precipitation or effective precipitation.This finding is consistent with the conceptual model of this paper that vegetation cover is the predominant climate-related variable that influences drainage density and that vegetation cover and drainage density are inversely related.Melton (1957), for example, documented an inverse correlation between drainage density and the precipitationeffectiveness (P-E) index at over 80 sites in the southwestern USA, including arid (low-elevation) and humid (highelevation) climates.A similar negative correlation between drainage density and MAP was found by Abrahams and Ponczynski (1984).Naively, one might expect more precipitation to result in greater channel incision and hence less contributing area between divides and valley heads (and hence a higher drainage density), all else being equal.Melton (1957), however, proposed that greater aridity results in a lower vegetation density and, hence, a reduction in the cohesive strength protecting soils on hillslopes, thus leading to a higher drainage density.One might also expect a lower vegetation density to increase the runoff-to-rainfall ratio and hence also increase drainage density, but runoff intensity varied by only a factor of two across Melton's sites while drainage density varied by nearly two orders of magnitude.Istanbulluoglu and Bras (2007) provided theoretical support for Melton's interpretation, illustrating that a lower vegetation densities can lead to higher drainage densities through the cohesive or anchoring effect of plant roots.

Implications for our understanding of erosion-climate linkages
There has been an ongoing debate in the geomorphic literature regarding the importance of climate (not limited to but often defined as MAP) on erosion rates.Given the significant correlation between MAP and erosion rates in many studies within individual mountain ranges (e.g., Reiners et al., 2003;Bookhagen and Strecker, 2012), it is perhaps surprising how little correlation exists between MAP and erosion rates in global compilation/synthesis studies (e.g., von Blanckenburg, 2005;Portenga and Bierman, 2011).Even studies that emphasize the climatic control on erosion rates note that R 2 values between erosion rates and MAP are quite small (e.g., Yanites and Kesler, 2015).
Recent work on the role of vegetation, and its changes through time, can provide a basis for understanding the relatively low correlation between erosion rates and MAP in unglaciated areas outside the dominant influence of periglacial processes and the complex relationship between erosion rates and climate in such areas more generally.For example, Torres Acosta et al. (2015) recently documented a negative correlation between erosion rates and both vegetation cover and MAP in Kenya.They proposed that the primary effect of more humid conditions is to increase vegetation cover on hillslopes, thereby reducing erosion rates on otherwise similar slope gradients.This concept is consistent with the classic Langbein and Schumm (1958) curve.Langbein and Schumm (1958) proposed that sediment yields are maximized in semi-arid climates (all else being equal) because such climates generate sufficient rainfall to detach and transport soil in overland/rill flow but insufficient vegetation cover to protect/anchor the soil.As MAP increases in this conceptual model, more precipitation is available to drive erosion, but this effect is more than offset by a decrease in the susceptibility of soil to erode due to the increased anchoring effect associated with greater plant cover/biomass.The results of this paper demonstrate further complexity in the erosion-climate relationship, i.e., that the change in climate (and hence of vegetation cover) can be as important or more important that its mean state.It is important to emphasize that the effect of vegetation on the rate of erosion by colluvial processes may be entirely different than its effect on fluvial/slope-wash processes.All else being equal, increased vegetation cover is likely to increase erosion rates by colluvial processes, since more plants can be expected to drive higher rates of bioturbation (e.g., Osterkamp et al., 2012).As such, it is crucial to consider colluvial and fluvial/slope-wash processes separately when considering the effects of vegetation on hillslope erosion.

Conclusions
In this study we leveraged all relevant data from a uniquely well-studied semi-arid watershed to test the hypothesis that late Holocene vegetation changes can modulate drainage density, hillslope-scale relief, and watershed-scale erosion rates.We documented that areas below 1430 m a.s.l.have decadal-scale erosion rates approximately 10 times higher, drainage densities approximately 3 times higher, and hillslope-scale relief approximately 3 times lower than elevations above 1430 m.We calibrated all the terms of a mathematical landscape evolution model and used the model to predict the equilibrium drainage density associated with shrublands and grasslands.Model predictions for the increase in drainage density associated with the shift from grasslands/woodlands to shrublands are broadly consistent with measured values.Using modern erosion rates and the magnitude of relief reduction associated with the transition from grasslands/woodlands to shrublands, we also estimated the timing of the grassland-to-shrubland transition in the lower elevations of WGEW to be approximately 3 ka, i.e., broadly consistent with constraints from paleovegetation studies.Our work provides a mathematical model for predicting equilibrium drainage density in transport-limited fluvial environments that may be applicable in other study sites.
www.earth-surf-dynam.net/4/471/2016/Earth Surf.Dynam., 4, 471-488, 2016 The DEM used in this paper can be obtained upon request from the corresponding author.All other data used in the paper are in the published literature.

Figure 1 .
Figure 1.Maps of the bedrock geology and geomorphology of Walnut Gulch Experimental Watershed (WGEW).(a) Bedrock geology from Osterkamp (2008): rectangle identifies the portion of WGEW that is the focus of this study.(b) Geomorphic map from Osterkamp (2008): the boundary between the Dissected Whetstone Pediment and Upper Whetstone Pediment marks a key transition in landscape morphology, soil type, and vegetation cover (Fig. 2).

Figure 2 .
Figure 2. Relationships among landscape morphology and vegetation cover in the study area.(a) Shaded relief image of the topography, illustrating the significant increase in hillslope-scale relief from the western to the eastern portion of the study area.(b) Grayscale map of topographic curvature (i.e., Laplacian), demonstrating generally lower absolute hillslope curvatures (i.e., more gray) in the western (shrubland) portion of the study area relative to the eastern (grassland) portion (more black).(c) Vegetation map, after Skirvin et al. (2008), identifying the areas that are primarily shrubland, grassland, and transitional between the two.

Figure 4 .Figure 5 .
Figure 4. Detailed shaded relief images (locations shown in Fig. 3) illustrating the presence of parallel hillslope rills and gullies in the shrubland areas (shown in a).Grassland areas (shown in b) have fewer such features.

Figure 6 .Figure 7 .
Figure 6.Plot of mean hillslope relief as a function of elevation, illustrating the marked increase in relief above elevations of approximately 1430 m a.s.l. in the study area.Each data point represents the mean hillslope relief in 10 m wide elevation bins.
Figure8.Plots of mean topographic curvature as a function of contributing area (distance from divide also shown along x axis using data in Fig.7).Error bars represent the standard error of the mean for each bin.Topographic curvatures are similar in shrublands and grasslands at small and large contributing areas, with a minimum of 0.03 m −1 near divides.Within a range of contributing areas from

Figure 9 .
Figure9.Plot of the decadal timescale volumetric sediment fluxes in shrublands measured at the five watersheds of Lucky Hills as a function of contributing area.The straight line is the result of a least-squares regression to the logarithms of both sides of Eq. (7), from which the values of k Q s s and p were constrained.

Figure 10 .
Figure10.Plots of the total erosion rate, E, the fluvial/slope-wash erosion rate, E f , and the colluvial deposition rate, −E c , as a function of distance along flow lines from divides, x.The values of x s and x g (where E = E f + E c ) are also shown.

Figure 11 .
Figure 11.Plots of the mean longitudinal profile of hillslopes and valley bottoms in shrubland and grassland areas, constructed by integrating the mean curvature data in Fig. 8.

Table 1 .
Sediment yield data and watershed characteristics.The erosion rate calculation assumes a soil bulk density of 1500 kg m −3 .

Table 2 .
List of model parameters and values.

Table 3 .
Measured (from DEM analysis) and predicted values (from Eqs.11 and 12) for the mean distance from divides to valley bottoms in shrublands (x s ) and grasslands (x g ).