Statistical analysis of drainage density from digital terrain data

Geomorphology 36 Ž2001. 187–202 www.elsevier.nlrlocatergeomorph

Statistical analysis of drainage density from digital terrain data Gregory E. Tucker ) , Filippo Catani 1, Andrea Rinaldo 2 , Rafael L. Bras Ralph M. Parsons Laboratory, Department of CiÕil and EnÕironmental Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Received 10 August 1999; received in revised form 18 September 2000; accepted 20 September 2000

Abstract Drainage density Ž Dd ., defined as the total length of channels per unit area, is a fundamental property of natural terrain that reflects local climate, relief, geology, and other factors. Accurate measurement of Dd is important for numerous geomorphic and hydrologic applications, yet it is a surprisingly difficult quantity to measure, particularly over large areas. Here, we develop a consistent and efficient method for generating maps of Dd using digital terrain data. The method relies on Ži. measuring hillslope flow path distance at every unchanneled site within a basin, and Žii. analyzing this field as a random space function. As a consequence, we measure not only its mean Žwhich is half the inverse of the traditional definition of drainage density. but also its variance, higher moments, and spatial correlation structure. This yields a theoretically sound tool for estimating spatial variability of drainage density. Averaging length-to-channel over an appropriate spatial scale also makes it possible to derive continuous maps of Dd and its spatial variations. We show that the autocorrelation length scale provides a natural and objective choice for spatial averaging. This mapping technique is applied to a region of highly variable Dd in the northern Apennines, Italy. We show that the method is capable of revealing large-scale patterns of variation in Dd that are correlated with lithology and relief. The method provides a new and more general way to quantitatively define and measure Dd , to test geomorphic models, and to incorporate Dd variations into regional-scale hydrologic models. q 2001 Elsevier Science B.V. All rights reserved. Keywords: drainage patterns; DEM; Apennines; morphometry

1. Introduction The degree to which a landscape is dissected by stream channels has long been recognized as a fun-

damental property of natural terrain. Horton Ž1932, 1945. defined drainage density as the total length of stream channels divided by the area they occupy Dd s L T rA

Ž 1.

)

Corresponding author. Present address: School of Geography and the Environment, Oxford University, Mansfield Road, Oxford, OX1 3TB, UK. Tel.: q44-1865-271-913; fax: q44-1865271-929. E-mail address: [email protected] ŽG.E. Tucker.. 1 Dipartimento di Scienze della Terra, Universita’ di Firenze, Firenze, Italy. 2 Dipartimento di Ingegneria Idraulica, Marittima e Geotecnica, Universita’ di Padova, Padova, Italy.

where L T is the total length of channels within area A. Drainage density, which has dimensions of inverse length, is physically related to the mean distance one has to walk from a random location before encountering a channel; that is, L ; 1rŽ2 Dd . ŽHorton, 1932.. Drainage density is also closely related to several other measures of landscape dissection, in-

0169-555Xr01r$ - see front matter q 2001 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 9 - 5 5 5 X Ž 0 0 . 0 0 0 5 6 - 8

188

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

cluding valley density ŽMontgomery and Dietrich, 1994; Rinaldo et al., 1995; Tucker and Bras, 1998. and channel-head source area ŽMontgomery and Dietrich, 1989.. Drainage density is known to vary with climate and vegetation ŽChorley, 1957; Chorley and Morgan, 1962; Melton, 1958; Gregory and Gardiner, 1975; Moglen et al., 1998., soil and rock properties ŽSmith, 1958; Carlston, 1963; Cotton, 1964; Wilson, 1971; Madduma Bandara, 1974; Kelson and Wells, 1989., relief ŽSchumm, 1956; Montgomery and Dietrich, 1988; Oguchi, 1997., and landscape evolution processes. The theoretical basis for these observed variations has been explored within the context of dynamic process models ŽKirkby, 1987; Montgomery and Dietrich, 1989; Willgoose et al., 1991b; Tarboton et al., 1992; Howard, 1997; Tucker and Bras, 1998., though testing of these models has been hampered by the difficulties involved in measuring Dd and characterizing its spatial variations. The question of whether Dd varies with spatial scale Že.g., basin size. ŽGardiner et al., 1977. remains unresolved, also largely for lack of an appropriate method of analyzing spatial variability in Dd . The close association between Dd and climatically driven erosional processes implies that Dd can act as a signature of process and landscape history ŽKirkby, 1987; Willgoose et al., 1991b; Rinaldo et al., 1995; Howard, 1997; Tucker and Slingerland, 1997; Tucker and Bras, 1998.. Furthermore, the association between Dd and rates of erosion or tectonic uplift means that Dd can also potentially be used to infer tectonic and geomorphic history. However, decoding these various signatures is not straightforward, in part because drainage density may not appropriately reflect the large spatial variability of hillslope lengths. Unlike continuous terrain properties such as elevation or slope gradient, Dd Žas conventionally defined. is a spatially averaged quantity. Yet when seen as a measure of hillslope flow path distance, it is a property that can and does vary in space; and these spatial variations are of great interest in reading signatures of various geomorphic processes, testing theoretical models of drainage density, and providing input to hydrologic models. In this paper, we introduce a consistent and objective method for analyzing and mapping drainage density using digital terrain data. The method starts

by defining a local and easily measurable property, the downslope distance to the nearest channel or valley from a given point. This quantity is then treated as a random space function. When hillslope flow path distance is averaged in space over a length scale equal to its autocorrelation scale, the result is a digital map of drainage density at a meaningful landscape scale. These methods are applied to the upper Reno drainage basin in the northern Apennines of Italy, for which detailed field data are available. We show that the resulting maps reveal systematic, basin-wide variations in Dd that are associated with properties such as relief and lithology. By combining such maps with maps of lithology, soils, topography, rainfall, and other properties, it becomes possible to test quantitative theories of drainage density and its controls.

2. Methods In conventional terrain analysis, three scalars are ordinarily computed for each site, say i, of a digital elevation model ŽDEM.: the local slope and the local curvature of the topographic field of elevations, and the total upstream contributing area. Such quantities are instrumental in defining two other scalar properties that provide the basis for our proposed mapping of drainage density, and therefore it seems appropriate to review some technicalities pertinent to their calculation. The local slope, <= z < i , of the topographic field z Žx. is easily computed for each site of a given landscape. In the discrete form naturally employed by DEMs, at a given site i, the slope, <= z < i , can be computed from the drop in elevation between the ith site and the lowest among its eight nearest neighbors. This maximum drop, identified from the steepest descent direction, is divided by either the grid size, D x Žif the orientation occurs along the grid coordinate directions., or by '2 D x Žif the orientation is diagonal.. The curvature at the ith site, =z i2 , is also an important geomorphic indicator. It is computed by suitably accurate numerical averaging schemes Žsee, e.g., Rodriguez-Iturbe and Rinaldo, 1997; Roering et al., 1999., typically involving all of the eight near-

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

est-neighbors in a rectangular lattice. The resulting accuracy is O Ž D x ) 2 ., where D x ) is the dimensionless size of the DEM Ži.e., divided by the size of the map.. Curvatures critically distinguish landscape sites. If =z i2 - 0, the site is topographically convex, or topographically divergent, and is conventionally defined as a hillslope. If =z i2 G 0, the site is topographically concave, or convergent, and typically corresponds to a valley, hollow, or basal hillslope sculpted by advective processes ŽKirkby, 1971; Smith and Bretherton, 1972; Dunne, 1980.. It is customary ŽDietrich et al., 1992, 1993; Howard, 1994. to consider topographic concavity of a site as a necessary but not sufficient condition for channelization to occur, and indeed unchanneled valleys have important geomorphic implications, in particular with reference to signatures of climate changes ŽDietrich and Dunne, 1993; Rinaldo et al., 1995; Tucker, 1996; Tucker and Slingerland, 1997.. Total contributing area, A i , at the ith site is also a significant local property. It is commonly computed by the recursive equation A i s Ý wi , j A j q a

Ž 2.

j

where A j is the contributing area of site j, a is the elementary area of a pixel, and wi, j is the element of the connectivity matrix, W, defined as wi, j s 0 if j is unconnected to i, and by 0 - wi, j F 1 if j is connected to Ži.e., drains to. i. Notice that, by flow continuity, S i wi, j s 1 for all j. Furthermore, wi, j s 1 if a single-flow direction issues from the ith node Žfor a review, see e.g., Rodriguez-Iturbe and Rinaldo, 1997.. Where the connectivity matrix W has only Boolean Žzero or unit. elements and continuity allows only for one element of a row to differ from zero, the matrix represents a tree Ži.e., a unique pattern links any site to the seed.. Where W is made up of real numbers, in the range 0–1 and constrained by flow continuity to add to the unit value for every row, we have two possible cases. If the maximum eigenvalue of W is zero, we deal with a network; and Eq. Ž2. is still valid. In any other case, circuits arise and total contributing area becomes meaningless ŽRodriguez-Iturbe and Rinaldo, 1997.. It is quite important to compute total contributing area properly

189

when dealing with divergent topography, and this involves suitable algorithms for handling multipleflow directions in areas of topographic divergence Ži.e., where =z 2 - 0. Že.g. Quinn et al., 1991; Costai Cabral and Burges, 1994; Tarboton, 1997.. In extracting a channel network from a DEM, the first task one faces is classifying the channelization status of each pixel. To do so, one needs to have assessed the issue of where channels begin ŽMontgomery and Dietrich, 1988, 1992.. Several criteria have been proposed, employing threshold support areas Že.g., Band, 1986; Tarboton et al., 1991., threshold slope-area values Že.g., Montgomery and Dietrich, 1988, 1992; Montgomery and Foufoula-Georgiou, 1993; Moore et al., 1988a,b; Dietrich et al., 1992, 1993; Prosser and Dietrich, 1995., or critical curvatures ŽHoward, 1994.. In our study area, we found through comparison with channel locations mapped in the field that in order to capture the strong heterogeneity of surface properties multiplethreshold criteria work best. This approach involves identifying channelized sites through a combination of drainage area, slope-area product, and curvature criteria ŽCatani et al., in preparation.. Regardless of the criteria chosen for channel identification, the first step in deriving drainage density from a DEM is to apply an automated procedure that separates channeled from unchanneled sites. Once this classification has been completed, it becomes possible to define a local Ži.e., at-a-point. measure of drainage density based on the distance of each unchanneled site from the network of channeled sites. Suppose the ith site is unchanneled. We define two characteristic lengths for the ith unchanneled pixel. The local hillslope length to the channel network, L i . The length L i is defined by following the steepest descent path downslope until a channeled pixel is first reached ŽFig. 1.. Such definition, however, hinges on the proper identification of the channel network and may be affected by uncertainty either in the mapped streams ŽMontgomery and Foufoula-Georgiou, 1993. or in the criteria used to estimate channel network extent from a DEM Že.g., in the values of the geomorphic thresholds A t , Dt or A i <=z < i .. Such uncertainty will generally be greatest in low relief areas or when using relatively coarse or noisy DEMs. v

(

190

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

steepest descent direction. before encountering a channel ŽHorton, 1945.. Since LŽx. and L Žx. can be reasonably considered random space fields, however correlated, we can focus not only on average distances but also on higher moments and spatial correlation properties. In particular, from the probability distribution functions of the lengths, say pŽ L. and pŽ L . Žthus considered as random space functions., one can compute the means, ² L: and ² L :, variances, sL2 and sL2 , and higher-order moments. The qth moment of the distribution is ² Lq : s

`

q

q j

H0 L p Ž L . d L A Ý L

Ž 3.

j

Fig. 1. Schematic illustration of hillslope flow path distance, L.

The local hillslope length to the valley network, Li . Analogous to L i , this length is defined by following the steepest descent path downslope until a concave pixel is first reached Ži.e., the first site encountered for which =z i2 G 0, or more generally =z i2 G C0 , where C0 is a suitably chosen minimum curvature.. The uncertainty on the length computation depends only on the accuracy of the estimation of the sign of the curvature, and this proves generally robust ŽRodriguez-Iturbe and Rinaldo, 1997.. On the other hand, because curvature is the second derivative of elevation, it can be quite sensitive to noise or other errors in DEM data and thus requires either high-quality data or appropriate smoothing procedures. Needless to say, although L is a measure of hydrologic significance Že.g., Dietrich and Dunne, 1993; Rodriguez-Iturbe and Rinaldo, 1997., it does not necessarily portray the actual extent of unchanneled flow paths due to possible noteworthy differences between valley and drainage densities Že.g., Rinaldo et al., 1995; Tucker and Slingerland, 1997.. Using the above procedures, two scalar fields are identified: L i Ža discrete sample of the continuous field LŽx.. and Li Žwith continuous counterpart L Žx... We argue that drainage density may be mapped using these scalar fields. In fact, drainage density is conventionally defined as the inverse of twice the average length one has to walk Žin the v

where the summation spans all unchanneled sites. Notice that analogous expressions are valid for the field L . Interestingly, one may color-code the distances L, L and then map them, thereby in essence mapping AlocalB drainage density. In fact, focusing on hillslope lengths to define local drainage density allows a pointwise determination of hydrologic paths —the larger the length to the network, the smaller the local drainage density, and conversely. It is also possible to define a hillslope length covariance function, which represents the degree of average correlation in L or L between points separated by a given distance l. Under the assumption of isotropy and second-order stationarity, the hillslope covariance function is defined by CLŽ l . A ² L Ž x . L Ž x q l . :ls< l < .

Ž 4.

The assumption of second-order stationarity deserves some comment. Height–height Ži.e., topographic. correlation functions have been theoretically proven stationary for particular surface evolution models and for many fractal geometries Že.g., Barabasi and Stanley, 1995. and have been empirically shown to be reasonably stationary over basin Žor larger. scales ŽChan and Rothman, 1999.. The stationarity of length covariances, which are derived from topographic properties, defy theoretical assessment; but our preliminary testing supports the validity of the stationarity assumption. Interestingly, CLŽ l. is a measure of the heterogeneity of drainage density within a given area; in particular, it is related to the variance by CLŽ0. A sL2 Žand similarly for C L .. It is of particular interest that, following random function theory, from

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

the covariance function, we can define AmicroscalesB and AmacroscalesB of the processes, which describe objectively the resulting heterogeneity ŽBras and Rodriguez-Iturbe, 1985.. The microscale, l) , defines the distance one may walk in any direction without significantly affecting the distance from the nearest channel Žmathematically, l) is the intercept of an osculating parabola centered in CLŽ0... Conversely, the macroscale, L, is the distance one has to walk, on average, away from any site to obtain a completely uncorrelated measure of the distance of any site to the first occurrence of a channel. The macroscale is defined by a suitable integral of the covariance functions, i.e. `

LA

H0 C Ž l. d l . L

Ž 5.

The two above scales, l) and L, define objectively the degree of heterogeneity in channel network den-

191

sity and the spatial variability in the characteristic scale of landscape dissection. The concept of hillslope-to-channel Žor valley. length generalizes the traditional definition of drainage density, which is based on mean values. Treating hillslope lengths as random functions indeed embeds the traditional definition Žaccording to which Dd ; 1rŽ2² L:.. but has the advantage of incorporating additional information about spatial variability that can be characterized in a statistical framework and directly compared with other spatial fields such as slope, aspect, vegetation, soils, and relief. Finally, we note that in our results we will be focusing on properties of both hillslope-to-channel length and hillslope-to-valley length. Similarities and departures between these two measures will be discussed, in particular in view of the robust identification characteristics of L .

Fig. 2. Map of study area showing bedrock lithology and drainage patterns. Lithologies are classified on the basis of internal structure, from well stratified to highly deformed ŽAchaoticB ..

192

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

3. Study area The study area is located in the northern sector of the Apennines Mountains ŽItaly., close to the principal drainage divide which separates Tuscany from Emilia-Romagna ŽFig. 2.. The region covers the upper half of the Reno River basin and is characterized by strong spatial variations in Dd that are correlated with bedrock geology, relief, and processes ŽFigs. 2 and 3.. In this region, the Apennines chain is predominantly made up of arenaceous and calcareous turbidite sequences and highly deformed argillaceous units of sedimentary and tectonic origin belonging to the Ligurian Formations ŽMerla, 1951; Abbate and Sagri, 1970; Kliegfield, 1979; Treves, 1984.. The area covered by the 20-m DEM used in the analysis, approximately 50 = 40 km wide, is roughly equally subdivided between regular well-stratified turbidite sequences Žbelonging to the Sestola-Vidiciatico Unit, the Porreta-Suviana sandstones, and the Mont

Cervarola sandstones. and the highly deformed, melange-like shaly rocks of the Basal Complexes ´ ŽBettelli and Panini, 1987, 1991.. The Northern Apennines have been subject to intense and nearly continuous uplift since the Pliocene and continue to be tectonically active ŽBartolini et al., 1982; Consiglio Nazionale delle Ricerche and Progetto Finalizzato Geodinamica, 1992.. During the middle Pleistocene, the southwestern sector Žto which the study area belongs. was progressively uplifted above sea level and subjected to intense erosion, particularly within the uplands along the Apenninian Divide ŽBoccaletti et al., 1985; Bertotti et al., 1997; Sorgi et al., 1998.. The landscape morphology displays a strong correlation with the underlying geology ŽFigs. 2 and 3.. Where the well-stratified units prevail, the relief tends to be characterized by steep hillsides and narrow deep valleys, especially in the southern area close to the main drainage divide ŽFigs. 2 and 3.. In the central portion of the study area where the

Fig. 3. Shaded relief image of digital elevation model of study area, with lighting from the northwest. Inset shows contour map. ŽTerrace-like features are artifacts of DEM generation..

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

structurally complex units are widespread, valleys are wider and hillsides are smoothly and gently sloped. Hillslope transport is dominated by mass failures, principally in the form of large, roto-translational landslides ŽCasagli and Catani, 1998.. Field mapping indicates that active or recently active landslides cover as much as 40% of the central region of the study area. The footslopes and the narrow irregular alluvial flats are shaped by the toes of such landslides, which often divert river courses. The type and frequency of slope movements are influenced by the relief, the underlying lithologies, and by the rather severe meteorological conditions: as the principal drainage divide splitting the Adriatic and Tyrrhenian sides of Italy, this part of the Apennine chain receives a mean annual precipitation of about 2800 mm. The area covered by the DEM used in this analysis is characterized by an unusually complicated drainage density pattern: the southern, high-relief portion of the Reno basin shows high drainage density; while much of the remaining area, especially where the highly deformed units prevail, exhibits lower values of Dd Žsee Fig. 2.. These differences are clearly connected with the bedrock geology, relief, and the different geomorphic processes cited above. The regional variability in factors associated with Dd makes this region well suited to analysis of spatial variations in drainage density and demonstrates the need for a more flexible parameterization of this important measure. 3.1. Digital terrain data A digital elevation model of the study region was generated through interpolation of a 1:25,000-scale contour map Ž50-m contour interval.. DEM grid spacing is 20 m. In addition to the DEM, a digitized 1:25,000-scale map of stream networks was used as the basis for channel locations. This stream network map was then mapped onto a 20-m grid using a standard vector-to-raster conversion procedure. Drainage directions and drainage areas were computed using a single-direction flow algorithm, based on the steepest downslope neighbor of each pixel Že.g., Jenson and Domingue, 1988; Tarboton et al., 1991.. The hillslope-to-channel Ž L. and hillslope-to-

193

valley Ž L . lengths were then computed by tracing the downslope path from each unchanneled or convex pixel to the first occurrence of a channelized or concave pixel, respectively, and recording the total path distance at each site.

4. Results The digital map of distances to the nearest channel for the upper Reno basin is shown in Fig. 4. Not surprisingly, hillslope-to-channel length is large near ridgelines and small along the bases of hillslopes and within valleys. Thus, the map highlights the ridgevalley network. In addition to these fine-scale patterns, large-scale variations in drainage density are apparent, with the darker gray regions Žlarger hillslope-to-channel distances. indicating areas of lower drainage density Žcompare with Fig. 2.. Fig. 5 indicates that L has an approximately exponential frequency distribution Žnote the semi-log axes. with a mean of 300 m and a standard deviation, sL , of 230 m. Just as hillslope-to-channel length, L, can be thought of as a local measure of drainage density, its mean value directly reflects the regional average drainage density. The exponential frequency distribution ŽFig. 5. implies that there are many short pathways from sites near channels and few long, contorted pathways Že.g., Fig. 1.. At the scale of an individual hillslope, L is strongly autocorrelated, as might be expected ŽFig. 4.. At larger scales, this correlation breaks down. The correlation scale is reflected in the covariance function ŽFig. 6., which indicates that L becomes essentially independent Žuncorrelated. beyond length scales of about 1000 m. Notice that this length scale, L, is closely related to regional average drainage density: the higher the drainage density, the shorter the scale over which L covaries. Thus, L provides a natural length over which L can be spatially averaged in order to map drainage density at a regional scale. Spatial averaging of L using a circular moving window of radius L s 1000 m results in a map of drainage texture at a scale between that of the purely local field, L, and the regional average drainage density ŽFig. 7.. The resulting map of LŽ L. highlights systematic variations in drainage texture that

194

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

Fig. 4. Map of hillslope-to-channel distances, L, computed from DEM and digitized streams. Note that L tends to be low near ridgelines and high near channels, thus highlighting drainage patterns. Large-scale variations in L, associated with variations in relief and bedrock geology, are evident.

are correlated with lithology and relief Žcompare with Figs. 2 and 3.. Just as L reflects the density of channels in a given area, L reflects the density of concave ele-

ments which, in most landscapes, are primarily associated with valleys. Thus, L can be thought of as a local measure of the degree of landscape dissection into concave valleys. Typically, L will be smaller

Fig. 5. Frequency distribution of L within the study area. The approximately straight trend on a semi-log plot suggests an exponential distribution.

Fig. 6. Covariance functions of L and L , where r is the length scale in meters Žsee text..

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

195

Fig. 7. Map of spatially averaged L generated by applying a circular moving average of radius L s 1000 m to the map shown in Fig. 4.

Fig. 8. Map of topographic curvature for the southern portion of the study area, derived from the DEM.

196

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

Fig. 9. Ža. Digital map of hillslope length-to-valley, L , for the southern portion of the study area. Valley locations estimated using curvature threshold Ct G 0.01. Žb. Map of spatially averaged L for study area, obtained using circular moving window of 1 km radius.

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

than L for three reasons: Ži. valleys are wider than the channels they contain, Žii. unchanneled valleys often occur upslope of channel heads Že.g., Dietrich and Dunne, 1993., and Žiii. the lower portions of hillslopes are often concave-upward and unchanneled Že.g., Dunne, 1980.. Depending on the quality of the DEM, operational problems can arise when calculating L . Small errors in the data can create sites of apparent convexity or concavity that do not exist in the field. In addition, interpolation procedures that are sometimes used in DEM generation can create concave Žand fictitious. artifacts such as AterracesB associated with contour lines Žsuch features appear in parts of the Reno dataset; see Fig. 3.. Finally, roto-translational landsliding on hillslopes can sometimes create subtle concavities that do not represent true valley forms. For these reasons, care must be taken in selecting a minimum threshold for concavity. A map of topographic curvature for the study area is shown in Fig. 8. Hillslope length to valley, L , was calculated from the curvature field by tracing flow paths from each convex or planar point to the nearest occurrence of a site with a curvature greater than a threshold value of Ž= z 2 . min s 0.01 ŽFig. 9a.. The mean value of L , 87 m, is smaller than that of L, as expected. The correlation scale is also shorter ŽFig. 6.. When L is averaged using a moving window of 1-km radius, the resulting map ŽFig. 9b. highlights both the promise and the challenge in mapping valley density based on L . Within the high-relief block in the southern portion of the DEM, the curvature criterion successfully delineates valleys and ridgelines very clearly ŽFigs. 8 and 9a.. In the lower-relief areas to the north, DEM processing artifacts in the form of contour-line terraces partly obscure the valley-ridge topography and apparently lead to underestimates of true hillslope-to-valley distance. Spatial averaging ŽFig. 9b. yields a map that reveals systematic differences in valley density between the northern and southern portions of the study region, though the map fails to distinguish between the high-Dd southern massif and the low-Dd central region. Operational problems notwithstanding, however, hillslope length-to-valley maps like Fig. 9 have the important advantage of being based solely on topography and do not require independent measurement of channel networks. Moreover, the difference between L and

197

L is likely to contain information about properties such as characteristic valley width and the magnitude of past variations in climate, as manifest in the extent of unchanneled valleys Že.g., Rinaldo et al., 1995; Tucker and Slingerland, 1997.. Drainage density is but one of several landscape attributes whose large-scale variations we might wish to map. Properties such as gradient, elevation, and soil moisture vary both on a fine scale, in relation to ridge-valley topography, and on a larger scale, reflecting lithology, soils, structurertectonics, and spatial variations in hydroclimatology Že.g., orographic effects.. One can conceive of applications, for example, in which it would be desirable to map average slope gradient on a scale larger than that of individual hillslopes but small enough to highlight systematic gradient variations related to lithology and structure. Fig. 10 shows a plot of hillslope gradient Žcomputed using a 3 = 3 cell weighted average. averaged using a moving window of radius L Žin this case, 1000 m.. The gradient map highlights the massif in the south of the basin and the locally steep topography associated with the NE–SW trending outcrop pattern ŽFig. 2.. Thus, the drainage correlation macroscale L provides a natural length scale for averaging out the fine-scale terrain roughness associated with ridge-valley topography in order to better resolve regional-scale trends in various landscape attributes, including but not limited to drainage density. One way to test our method is to apply it to terrain in which drainage density and its spatial variations are known independently to a high degree of accuracy. Needless to say, this can be difficult to achieve in the field. As an alternative, however, we can exploit recent advances in process-based modeling Že.g., Willgoose et al., 1991a,b; Howard, 1994; Tucker and Slingerland, 1997. to generate simulated landscapes in which valley-head source area Žas a surrogate for drainage density; Montgomery and Dietrich, 1989. is exactly known ŽTucker and Bras, 1998.. Fig. 11 shows a simulated terrain that was generated using three simple rules: Ži. the rate of vertical fluvial erosion at any point is equal to K'A S, where A is drainage area, S is gradient in the downslope direction, and K is an erosion coefficient that reflects climate, hydrology, and lithology ŽWhipple and Tucker, 1999; Tucker and Bras, 2000.;

198

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

Fig. 10. Image showing average surface gradient obtained by applying a circular moving average Žradiuss 1000 m. to the local gradient field Žlocal gradients calculated on a 3 = 3 cell basis..

Žii. landsliding imposes an upper limit to gradient, Scr , that cannot be exceeded ŽTucker and Slingerland, 1994.; and Žiii. the landscape is subject to a steady rate of uplift Žor equivalently, base level lowering. until equilibrium between uplift and erosion is reached. Under this set of conditions, Howard Ž1997. and Tucker and Bras Ž1998. have shown that valley-head source area, A vh , is determined by

A vh s

Fig. 11. Hypothetical terrain with spatially varying drainage density, simulated using the GOLEM model ŽTucker and Slingerland, 1996, 1997.. Diagonal band of low drainage density reflects a reduced value of the critical slope threshold Scr Žfrom 508 to 308.. Channels are drawn in at the landslide-fluvial transition Žsee Tucker and Bras, 1998 for details.. Side boundaries are connected to form a periodic boundary condition. The spatial scale is arbitrary; nominally, it covers a 2560=2560 m region Ž128=128 20-m pixels..

U

ž / KScr

2

.

Ž 6.

In Fig. 11, spatial variations in drainage density are simulated by varying the rock strength parameter Scr : a diagonal band of AweakB rocks Ž Scr s 308. is inset into a region of AstrongB rocks Ž Scr s 508.. Fig. 12a shows length-to-channel Ž L. for the simulated terrain. Averaging L over a moving window of width L s 20 pixels Žequal to the covariance macroscale. results in a map of drainage density that successfully recaptures the main features of the imposed pattern of lithologies ŽFig. 12b..

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

199

Fig. 12. Ža. Map of hillslope length-to-channel derived from Fig. 11. Žb. Contour map of average drainage density, LŽ L. obtained by applying a moving average of width L s 20 pixels to the map in Ža.. Units are in pixels. Compare with Fig. 11.

200

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

5. Discussion and conclusions The hillslope length method provides a simple and straightforward way to analyze drainage density statistically and to map variations in landscape texture across a given region using DEM data. The two proposed metrics, hillslope length-to-channel and length-to-valley, are complementary. Of the two, hillslope length-to-channel is closest to the traditional definition of drainage density, but as we have suggested, our treatment of L as a random space function with spatial correlation introduces a new and more general definition. Computing L requires both elevation data and either accurate stream network data or a reliable procedure for extracting the channelized portion of the DEM. This drawback, however, is compensated for by the increasing quality, resolution, and availability of DEM data; and given reliable channel network data, the hillslope length-to-channel statistics are quite robust and straightforward to compute. The second statistic, hillslope length-to-valley, measures valley density rather than channel network density ŽDietrich and Dunne, 1993; Montgomery and Dietrich, 1994.. An important difference is that valley density is a morphological rather than a hydrological property per se. Valley density reflects the long-term integrated effects of advective Žvalley-forming. and diffusive Žvalley-filling. processes over periods of perhaps tens to hundreds of thousands of years rather than the present-day extent of active channel features, which can fluctuate in response to short-term climate change or landscape disturbance. Thus, differences between the length-to-channel and length-to-valley statistics may contain information about paleoclimate and landscape history. Hillslope length-to-valley relies on the proper identification of the sign of topographic curvature, which can be sensitive to DEM errors or noise. With sufficiently high quality data, however, curvature can be reliably calculated from DEMs Že.g., Roering et al., 1999; Fagherazzi et al., 1999; Rinaldo et al., 1999., suggesting that operational problems with the current generation of DEM data will disappear as radar- and laser-based terrain data become increasingly available. Distributed maps of drainage density can be used, among other things, to test theoretical predictions

of landscape texture. When Horton’s Ž1932, 1945. well-known papers inaugurated the field of quantitative morphometry in the 1950s and 1960s, drainage density was one of the key variables that received attention. At that time, however, practitioners encountered two significant challenges: the difficulty in making measurements by hand from maps over large areas; and the difficulty in separating key controlling factors such as climate, lithology, and vegetation. The automated methods introduced here go a long way toward solving the first problem. They may also mitigate the second problem by making it possible to collect statistically significant data sets within regions of uniform climate, relief, or lithology. Moreover, by treating drainage density as a random field, it becomes possible to analyze variability in Dd within a quantitative framework. Recently, a number of theoretical models of drainage density have been developed. These theories link the mechanics of different hillslope and channel processes to the expected variations in drainage density Žor surrogates of drainage density such as channel-head source area. as a function of variables such as slope, relief, hydrology, and climate. The proposed automated procedures for measuring channel network and valley density from DEM data provide a natural way to test these theories over large areas, where the key controlling factors can also be mapped. By linking small-scale process mechanics with large-scale morphologic outcomes, such tests promise to place new constraints on models of landscape evolution.

Acknowledgements This project was supported by the Consiglio Nazionale delle Ricerche ŽItalian National Research Council. through the CNR-MIT Cooperative Agreement for the Study of Climatic Change and Hydrogeologic Risks. AR is funded by MURST 40% Morfodinamica Fluviale e Costiera, and by Universita’ di Padova. GT and RB gratefully acknowledge additional support from the US Army Research Office ŽDAAH 04-95-1-0181. and the US Army Corps of Engineers ŽDACA 88-95-C-0017.. The views expressed in this paper are not an official position of the Department of Defense.

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

References Abbate, E., Sagri, M., 1970. The eugeosynclinal sequences. Sediment. Geol. 4, 251–340. Band, L.E., 1986. Topographic partition of watersheds with digital elevation models. Water Resour. Res. 22, 15–24. Barabasi, A.L., Stanley, H.E., 1995. Fractal Concepts in Surface Growth. Cambridge Univ. Press, NY. Bartolini, C., Bernini, M., Carloni, G.C., Costantini, A., Federici, P.B., Gaspari, G., Lazzarotto, A., Marchetti, G., Mazzanti, R., Papani, G., Pranzini, G., Rau, A., Sandrelli, F., Vercesi, P.L., Castaldini, D., Francavilla, F., 1982. Carta neotettonica dell’Appennino settentrionale. Note Illus., Boll. Soc. Geol. Ital. 101, 523–549. Bertotti, G., Capozzi, R., Picotti, V., 1997. Extension controls Quaternary tectonics, geomorphology and sedimentation of the N-Apennines foothills and adjacent Po plain ŽItaly.. Tectonophysics 282 Ž1–4., 291–301. Bettelli, G., Panini, F., 1987. I melanges dell’Appennino settentri´ onale dal T. Tresinaro al T. Sillaro. Mem. Soc. Geol. Ital. 39, 187–214. Bettelli, G., Panini, F., 1991. Liguridi, melanges e tettoniti nel ´ Complesso Caotico lungo la linea del Sillaro. Mem. Descr. Carta Geol. Ital. 46, 387–415. Boccaletti, M., Coli, M., Eva, C., Ferrari, G., Giglia, G., Lazzarotto, A., Merlanti, F., Nicolich, R., Papani, G., Postpischl, D., 1985. Considerations on the seismotectonics of the Northern Apennines. Tectonophysics 117, 7–38. Bras, R.L., Rodriguez-Iturbe, I., 1985. Random Functions and Hydrology. Addison-Wesley, Reading, MA. Carlston, C.W., 1963. Drainage Density and Stream Flow: U.S.G.S. Professional Paper No. 422-C, pp. 1–8. Casagli, N., Catani, F., 1998. Fractal dimension and mechanical properties of melanges in the Northern Apennine. In: De Frede ´ ŽEd.., Proc. 4th Cong. of the Int. Assoc. Math. Geol. Neaples, vol. 2, pp. 561–571. Chan, K., Rothman, D.H., 1999. Coupled length scales in eroding landscapes. preprint. Chorley, R.J., 1957. Climate and morphometry. J. Geol. 65, 628–638. Chorley, R.J., Morgan, M.A., 1962. Comparison of morphometric features, Unaka Mountains, Tennessee and North Carolina, Dartmoor, England. Geol. Soc. Am. Bull. 73, 17–34. Consiglio Nazionale delle Ricerche and Progetto Finalizzato Geodinamica, 1992. Structural model of Italy and gravity map. Quad. Ric. Sci. 114. Costa-Cabral, M., Burges, S.J., 1994. Digital elevation model networks: a model of flow over hillslopes for computation of contributing and dispersal areas. Water Resour. Res. 30, 1681–1692. Cotton, C.A., 1964. The control of drainage density. N. Z. J. Geol. Geophys. 7, 348–352. Dietrich, W.E., Dunne, T., 1993. The channel head. In: Beven, K., Kirkby, M.J. ŽEds.., Channel Network Hydrology. Wiley, New York, pp. 175–219. Dietrich, W.E., Wilson, C.J., Montgomery, D.R., McKean, J.,

201

Bauer, R., 1992. Erosion thresholds and land surface morphology. Geology 20, 675–679. Dietrich, W.E., Wilson, C.J., Montgomery, D.R., McKean, J., 1993. Analysis of erosion thresholds, channel networks, and landscape morphology using a digital terrain model. J. Geol. 101, 259–278. Dunne, T., 1980. Formation and controls of channel networks. Progress Phys. Geogr. 4, 211–239. Fagherazzi, S., Bortoluzzi, A., Dietrich, W.E., Adami, A., Lanzoni, S., Marani, M., Rinaldo, A., 1999. Tidal Networks 1. Automatic network extraction and preliminary scaling features from digital terrain maps. Water Resour. Res. 35, 3891–3904. Gardiner, V., Gregory, K.J., Walling, D.E., 1977. Further notes on the drainage density-basin area relationship. Area 9, 117–121. Gregory, K.J., Gardiner, V., 1975. Drainage density and climate. Z. Geomorphol. N.F. 19, 287–298. Horton, R.E., 1932. Drainage basin characteristics. Am. Geophys. Union, Trans. 13, 348–352. Horton, R.E., 1945. Erosional development of streams and their drainage basins; hydrophysical approach to quantitative morphology. Geol. Soc. Am. Bull. 56, 275–370. Howard, A.D., 1994. A detachment-limited model of drainage basin evolution. Water Resour. Res. 30, 2261–2285. Howard, A.D., 1997. Badland morphology and evolution: interpretation using a simulation model. Earth Surf. Processes Landforms 22, 211–227. Jenson, S.K., Domingue, J.O., 1988. Extracting topographic structure from digital elevation data for geographic information system analysis. Photogramm. Eng. Remote Sens. 54 Ž11., 1593–1600. Kelson, K.I., Wells, S.G., 1989. Geologic influences on fluvial hydrology and bedload transport in small mountainous watersheds, northern New Mexico, USA. Earth Surf. Processes Landforms 14, 671–690. Kirkby, M.J., 1971. Hillslope process–response models based on the continuity equation. Slopes, Form and Process, Institute of British Geographers Special Publication, vol. 3, pp. 15–30. Kirkby, M.J., 1987. Modelling some influences of soil erosion, landslides and valley gradient on drainage density and hollow development. Catena Suppl. 10, 1–14. Kliegfield, R., 1979. The Northern Apennines as a collisional orogen. Am. J. Sci. 279, 676–691. Madduma Bandara, C.M., 1974. Drainage density and effective precipitation. J. Hydrol. 21, 187–190. Melton, M.A., 1958. Correlation structure of morphometric properties of drainage systems and their controlling agents. J. Geol. 66, 442–460. Merla, G., 1951. Geologia dell’Appennino settentrionale. Boll. Soc. Geol. Ital. 70 Ž1., 95–382. Moglen, G.E., Eltahir, E.A., Bras, R.L., 1998. On the sensitivity of drainage density to climate change. Water Resour. Res. 34, 855–862. Montgomery, D.R., Dietrich, W.E., 1988. Where do channels begin? Nature 336, 232–234. Montgomery, D.R., Dietrich, W.E., 1989. Source areas, drainage density and channel initiation. Water Resour. Res. 25, 1907– 1918.

202

G.E. Tucker et al.r Geomorphology 36 (2001) 187–202

Montgomery, D.R., Dietrich, W.E., 1992. Channel initiation and the problem of landscape scale. Science 255, 826–830. Montgomery, D.R., Dietrich, W.E., 1994. A physically based model for the topographic control on shallow landsliding. Water Resour. Res. 30, 1153–1171. Montgomery, D.R., Foufoula-Georgiou, E., 1993. Channel network source representation using digital elevation models. Water Resour. Res. 29, 3925–3934. Moore, I.D., Burch, G.J., MacKenzie, D.H., 1988a. Topographic effects on the distribution of surface soil water and the location of ephemeral gullies. Am. Soc. Agric. Eng. Trans. 31, 1098–1107. Moore, I.D., O’Loughlin, E.M., Burch, G.J., 1988b. A contourbased topographic model for hydrological and ecological applications. Earth Surf. Processes Landforms 13, 305–320. Oguchi, T., 1997. Drainage density and relative relief in humid steep mountains with frequent slope failure. Earth Surf. Processes Landforms 22, 107–120. Quinn, P., Beven, K.J., Chevallier, P., Planchon, O., 1991. The prediction of hillslope flow paths for distributed hydrological modeling using digital terrain models. Hydrol. Processes 5, 59–80. Prosser, I.P., Dietrich, W.E., 1995. Field experiments on erosion by overland and their implication for a digital terrain model of channel initiation. Water Resour. Res. 31, 2867–2876. Rinaldo, A., Dietrich, W.E., Rigon, R., Vogel, G.K., RodriguezIturbe, I., 1995. Geomorphological signatures of varying climate. Nature 374, 632–636. Rinaldo, A., Fagherazzi, S., Lanzoni, S., Marani, M., Dietrich, W.E., 1999. Tidal networks 3: studies in empirical geomorphic relationships. Water Resour. Res. 35, 3919–3929. Rodriguez-Iturbe, I., Rinaldo, A., 1997. Fractal River Basins: Chance and Self-Organization. Cambridge Univ. Press, New York. Roering, J.J., Kirchner, J.W., Dietrich, W.E., 1999. Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology. Water Resour. Res. 35, 853–870. Schumm, S.A., 1956. Evolution of drainage systems and slopes in badlands at Perth Amboy, New Jersey. Geol. Soc. Am. Bull. 67, 597–646. Smith, K.G., 1958. Erosional processes and landforms in Badlands National Monument, South Dakota. Geol. Soc. Am. Bull. 69, 975–1008. Smith, T.R., Bretherton, F.P., 1972. Stability and the conservation of mass in drainage basin evolution. Water Resour. Res. 8, 1506–1529.

Sorgi, C., Deffontaines, B., Hippolyte, J.C., Cadet, J.P., 1998. An integrated analysis of transverse structures in the northern Apennines, Italy. Geomorphology 25, 193–206. Tarboton, D.G., 1997. A new method for the determination of flow directions and contributing areas in grid digital elevation models. Water Resour. Res. 33, 309–319. Tarboton, D.R., Bras, R.L., Rodriguez-Iturbe, I., 1991. On the extraction of channel networks from digital elevation data. Hydrol. Processes 5, 81–100. Tarboton, D.R., Bras, R.L., Rodriguez-Iturbe, I., 1992. A physical basis for drainage density. Geomorphology 5, 59–76. Treves, B., 1984. Orogenic belts as accretionary prism: the example of the Northern Apennine. Ofioliti 9 Ž3., 577–618. Tucker, G.E., 1996. Modeling the large-scale interaction of climate, tectonics and topography. Tech. Rep. 96-003, Earth System Science Center, Pennsylvania State Univ., University Park, 223 pp. Tucker, G.E., Bras, R.L., 1998. Hillslope processes, drainage density and landscape morphology. Water Resour. Res. 34 Ž10., 2751–2764. Tucker, G.E., Bras, R.L., 2000. A stochastic approach to modeling the role of rainfall variability in drainage basin evolution. Water Resour. Res. 36, 1953–1964. Tucker, G.E., Slingerland, R.L., 1994. Erosional dynamics, flexural isostasy, and long-lived escarpments: a numerical modeling study. J. Geophys. Res. 99, 12229–12243. Tucker, G.E., Slingerland, R.L., 1996. Predicting sediment flux from fold and thrust belts. Basin Res. 8, 329–349. Tucker, G.E., Slingerland, R.L., 1997. Drainage basin responses to climate change. Water Resour. Res. 33, 2031–2047. Whipple, K.X., Tucker, G.E., 1999. Dynamics of the stream power river incision model: implications for height limits of mountain ranges, landscape response timescales and research needs. J. Geophys. Res. 104, 17661–17674. Willgoose, G.R., Bras, R.L., Rodriguez-Iturbe, I., 1991a. A physically based coupled network growth and hillslope evolution model: 1. Theory. Water Resour. Res. 27, 1671–1684. Willgoose, G.R., Bras, R.L., Rodriguez-Iturbe, I., 1991b. A physically based coupled network growth and hillslope evolution model: 2. Nondimensionalization and applications. Water Resour. Res. 27, 1671–1684. Wilson, L., 1971. Drainage density, length ratios, and lithology in a glaciated area of southern Connecticut. Geol. Soc. Am. Bull. 82, 2955–2956.