Effects of catchment discretization on topographic index distributions

Effects of catchment discretization on topographic index distributions

Journal of Hydrology (2008) 359, 150– 163 available at www.sciencedirect.com journal homepage: www.elsevier.com/locate/jhydrol Effects of catchment...

1MB Sizes 2 Downloads 34 Views

Journal of Hydrology (2008) 359, 150– 163

available at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/jhydrol

Effects of catchment discretization on topographic index distributions Santosh K. Aryal a b

a,*

, Bryson C. Bates

b

CSIRO Land and Water, Floreat, Western Australia 6014, Australia CSIRO Marine and Atmospheric Research, Floreat, Western Australia 6014, Australia

Received 24 January 2008; received in revised form 13 June 2008; accepted 24 June 2008

KEYWORDS Digital elevation models; Hillslopes; Terrain analysis; Topographic index

Digital elevation models (DEMs) are widely used in topographically-based hydrology models. Most researchers dealing with DEMs have worked exclusively with either grid or contour based techniques. Intercomparisons of these techniques are empirical in nature, few in number and limited to a small number of catchments. We examine the relative effects of these techniques on the distribution of the topographic index (TI) for nine idealised hillslope types and eight natural catchments (4.4–83 km2) in southeast Australia. In particular we investigate the underlying causes for the frequently observed increase in the TI with grid size. A simple analytical model for a parallel-planar hillslope indicates that median TI is asymptotic with respect to the number of grid cells and a linear function of grid size. We demonstrate that the latter result holds for hillslopes with more complex geometry. Using a multiple flow direction algorithm, we find that hillslope geometry and flow partitioning affect the spatial pattern of the TI and the rate of change (M) of median TI with grid size. This makes it impractical to recommend general guidelines for the value of the directional parameter in flow partitioning algorithms. For the natural catchments we show that median TI is also a linear function of grid size and that there is some evidence for a linear relationship between M and the inverse of the median topographic slope. Throughout the investigation we find that contour DEMs are less sensitive to changes in resolution than grid DEMs. ª 2008 Elsevier B.V. All rights reserved. Summary

Introduction Grid and contour DEMs * Corresponding author. Tel.: +61 8 93336318; fax: +61 8 93336211. E-mail addresses: [email protected] (S.K. Aryal), Bryson. [email protected] (B.C. Bates).

The prediction of landscape saturation patterns using terrain attributes derived from digital elevation models (DEMs) is still an evolving area of science. The DEMs most commonly

0022-1694/$ - see front matter ª 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2008.06.025

Effects of catchment discretization on topographic index distributions used in hydrological applications are grid based and contour based (Moore and Grayson, 1991; Costa-Cabral and Burges, 1994; Tarboton, 1997; Wilson and Gallant, 2000). More recently, Vivoni et al. (2004, 2005) have investigated the use of triangulated irregular networks (TIN) algorithm to study catchment response. Although the use of grid DEMs is more prevalent due to their simplicity and ease of implementation (Moore et al., 1986,1991; Beven et al., 1988; O’Loughlin, 1990; Iorgulescu and Jordan, 1994; Wolock, 1995; Moore and Thomson, 1996; Freer et al., 1997; Woods et al., 1997; Woods and Sivapalan, 1997; Western et al., 1999; Scanlon et al., 2000), neither approach is usable in catchments with large flat areas as contours cannot describe their topography efficiently. Grid DEMs exploit the fact that digital terrain data are made available in the form of either a three-dimensional elevation map on a regular grid or digital contour maps which are easily transformed into grid DEMs. Several methods are available to estimate flow directions from grid DEMs, including D8 (O’Callaghan and Mark, 1984), Rho8 (Fairfield and Leymaire, 1991), DEMON (Costa-Cabral and Burges, 1994), and D1 (Tarboton, 1997). Single flow direction methods such as D8 have been criticised because they cannot model flow dispersion: they only allow outflow from a grid cell to its nearest neighbour in the direction of steepest descent (Quinn et al., 1991; Costa-Cabral and Burges, 1994; Desmet and Govers, 1996; Moore, 1996; Zhou and Liu, 2002; Pan et al., 2004). Researchers have attempted to remove this limitation by using multiple flow direction algorithms where the outflow from one grid cell passes not only to its neighbour of steepest descent but also to other neighbours (Freeman, 1991; Holmgren, 1994; Quinn et al., 1991, 1995; Wolock and McCabe, 1995; Tarboton, 1997; Gallant and Wilson, 2000). Contour DEMs divide landscapes into irregularly shaped polygons (elements) based on contour lines and their orthogonals. They use a one-dimensional flow principle where flow accumulates downslope along flow tubes or flow trajectories (O’Loughlin, 1986; Moore and Grayson 1991; Moore, 1996; Maunder, 1999; Wilson and Gallant, 2000). Although flow convergence (or divergence) due to topographic control is efficiently handled, the computation of terrain attributes for large catchments can be problematic as the irregular element sizes caused by an irregular contour spacing require frequent manual intervention (Moore, 1996). Such intervention includes the addition of extra contours in between regular contour intervals and localised modification of contour shapes to avoid invalid flow lines (Maunder, 1999). In large catchments the extent of these interventions can become large and unmanageable (Moore, 1996; Mizukoshi and Aniya, 2002). Also, contour DEMs tend to produce large elements in gully areas resulting in high contributing areas spread across wide valley bottoms. Gallant and Wilson (2000, p. 82) describe an ad hoc method for dealing with this situation. Nevertheless contour DEMs are regarded by some as the superior approach for hillslope flow computation (Wise, 2000).

Topographic index The topographic index (TI) is an important terrain attribute as it parameterizes the effect of topography on runoff generation (Beven and Kirkby, 1979; O’Loughlin, 1981; Moore,

151

1996; Woods et al., 1997; Hjerdt et al., 2004). Three common forms of the TI are 8 > < as = tan b ðaÞ TI ¼ lnðas = tan bÞ ðbÞ ð1Þ > : as =K tan b ðcÞ where as = a/b is the specific contributing area in which a is the upslope contributing area and b is the contour or grid cell length subtending a at a given location, tan b is the local slope, ln is the Napierian logarithm, and K is the lateral saturated hydraulic conductivity of the soil. The TI is a physically-based indicator of a landscape’s propensity to develop and maintain steady-state saturation at a given location for a given soil type. Areas within a catchment possessing the same value of the TI are assumed to have the same hydrologic response to rainfall, regardless of their position in the landscape.

Effects of DEM attributes on terrain analysis and hydrological modelling It is widely reported that grid size selection has a marked influence on the spatial pattern and statistical distribution of the TI (Quinn et al. 1991, 1995; Wolock and Price, 1994; Zhang and Montgomery, 1994; Band et al. 1995; Bruneau et al., 1995; Woods and Sivapalan, 1997; Wilson and Gallant, 2000; Wilson et al., 2000; Grayson and Blo ¨schl, 2001). The works of Wolock and Price (1994), Zhang and Montgomery (1994), Band et al. (1995) and Bruneau et al. (1995) indicate that differences between the frequency distributions of TI data obtained using various DEM resolutions are due primarily to differences in location (central tendency), and that any differences in shape are slight. In general, an increase in grid size leads to an increase in as, a decrease in tan b, and an increase in the mean of ln (as/ tan b). Contour DEMs for natural catchments also result in variable mean TI when the number of elements is changed (Vertessy et al. 1993; Band et al., 1995; Chirico et al. 2003). Different topographically-based hydrology models or modelling concepts such as MIKE-SHE, TOPMODEL, THALES and TOPOG use different DEM techniques. There is evidence that DEM selection and resolution can have noticeable effects on model predictions of streamflow hydrographs (Band et al., 1995), overland flow and water table depth (Wolock and Price, 1994) and catchment average soil moisture (Chirico et al., 2003) if the model parameters are not adjusted to compensate for the effects of grid size. Most of the above knowledge has been obtained through empirical investigations. Comparisons of contour and grid DEMs are few in number and typically limited to a few catchments with drainage areas less than 20 km2 (see e.g. Band et al., 1995; Chirico et al., 2003).

Objectives and outline of this paper Clearly there is a need for a systematic evaluation and intercomparison of current DEM techniques. Since landscape saturation patterns vary at hillslope and catchment scales, it is sensible to explore the issues described above using elevation data from a comprehensive set of idealised hillslopes with different geometries, and small- to moderately-sized

152 natural catchments (1–100 km2). Thus the objectives of our investigation are to: (1) identify the inherent reasons why the TI varies with grid size, (2) assess the extent to which this variation changes from hillslope to hillslope and from catchment to catchment, and identify the causative factors for these changes, (3) examine the effect of flow partitioning on the spatial pattern of TI and the variation of the TI with grid size, and (4) compare the performances of grid and contour DEMs. We focus on the median of TI (referred to as median TI in the rest of the paper) since changes in grid size primarily affect the location of its frequency distribution and our data sets contained a few extreme observations (close to infinity) that would unduly influence the mean. The outline of the paper is as follows. We describe our methods and the idealised hillslopes and natural catchments used in our study in the section ‘‘Methods and data’’. The section ‘‘Results for idealised hillslopes’’ begins with the development of an analytical model that describes the variation of median TI with the number of grid cells and grid size for the simplest hillslope considered and a single flow direction algorithm. This is followed by investigations of the extent to which this description applies for the combined effects of more complex hillslope geometries and a multiple flow direction approach, and the effect of flow partitioning on the spatial pattern of TI. The section ‘‘Results for natural catchments’’ describes the results obtained from the extension of the study to natural catchments. We also assess the strengths of empirical relationships between the rate of change (M) of median TI with grid size, and M and median topographic slope. We discuss the relevance of our research to prevailing ideas on flow partitioning and the calibration of topographically-based hydrology models in the section ‘‘Results’’. Finally, our conclusions are presented in the section ‘‘Conclusions’’.

Methods and data DEM techniques Initially, we attempted to use three DEM techniques in our investigation: one contour-based method and two gridbased methods (a single flow direction algorithm and the D1 approach). The contour-based technique described by O’Loughlin (1986) and Vertessy et al. (1993) was used to extract topographic attributes from the idealised hillslopes and the natural catchments. The single flow direction algorithm was used solely to facilitate the development of an analytical model that describes the variation of median TI with the number of grid cells and grid size for the simplest hillslope geometry considered. The D1 method was selected because of its abilities to minimise unrealistic dispersion in the calculation of as and handle ridges, pits and flat areas in natural catchments (Tarboton, 1997; Zhou and Liu, 2002). Moreover, Endreny and Wood (2003) found that D1 and 2D-Lea, an algorithm which provided the building blocks for DEMON, achieved highest accuracy ranking among five flow direction algorithms for spatial congruence between observed and DEM-delineated overland flow networks for two agricultural hillslopes. However, application of D1 to the idealised hillslopes considered herein gave rise to severe

S.K. Aryal, B.C. Bates edge contamination. This phenomenon arises because flow direction cannot be calculated for grid cells at hillslope edges. This problem extended further into the domain when the flow was inward from a hillslope edge. Consequently, we applied a multiple flow direction algorithm to the idealised hillslopes because of its superior performance for simple landscape elements (Pan et al., 2004); and the opportunities it provided to test the generality of the analytical model and to explore the interplay between the TI, grid resolution and flow partitioning. We applied the D1 method to the natural catchments for the reasons given above. The parameterization of TI given by Eq. (1)a was used throughout. We used a multiple flow algorithm based on the formulae proposed by Freeman (1991) and Holmgren (1994) for apportioning outflow from a central grid cell to a lower neighbour. maxð0; Svi Þ ui ¼ Pr v j¼1 maxð0; Sj Þ

ð2Þ

Here, i and j are indices for neighbouring cells with lower elevation, /i is the flow proportion (0–1) assigned to cell i, Si is the slope between the central cell and cell i, r is the number of neighbouring cells that outflow can be apportioned to, and v is a directional parameter such that outflow is fully multidirectional as v ! 0 and fully unidirectional as v ! 1. Quinn et al. (1995) note that use of v > 10 gives an approximation to a single flow direction algorithm. We set v = 1, 3, 5, 7 and 10, which encompass the values of v considered by Freeman (1991), Quinn et al. (1991, 1995), Holmgren (1994), Wolock and McCabe (1995) and Gallant and Wilson (2000). Several definitions of the flow width used to compute as have been proposed (see e.g. Quinn et al., 1991; Costa-Cabral and Burges, 1994; Wolock and McCabe, 1995). Chirico et al. (2005) evaluated different alternatives for defining flow width and concluded that invariant flow width equal to the grid size is generally the best approach in most practical circumstances. We adopted this approach using equal flow widths in cardinal and diagonal directions (see Tarboton, 1997). We used the TauDEM software to implement the D1 approach (http://hydrology.neng.usu.edu/taudem/). TauDEM uses cubic interpolation to derive regular grids with different resolutions. The D1 approach represents flow direction as a single angle taken as the steepest downward slope on the eight triangular facets formed in a 3 · 3 grid cell window centred on the cell of interest. Flow is never proportioned between more than two downslope cells. The flow proportion from the upslope cell to a downslope cell is determined by how close the flow angle is to the angle formed by a line connecting their centres. The DEM routine embedded in the TOPOG model was used to implement a contour-based method (http:// www.per.clw.csiro.au/topog/). In this method the number of elements can be changed by adjusting either the contour interval or the flow tube spacing. We used the latter approach throughout our investigation (see also, e.g., Vertessy et al., 1993). The element network is topologically sound as it is arranged around peaks, saddles, and ridge and valley axes. Flow diffusion between adjacent flow tubes is not permitted. Since element size changes along a flow tube, we

Effects of catchment discretization on topographic index distributions compare the performances of contour- and grid-based techniques on the basis of the number of computational elements used (i.e., the number of grid cells in a grid DEM or the number of elements in a contour DEM).

Idealised hillslopes We investigated the effects of DEM resolution and variations in v on the TI for idealised hillslopes with concave, convex and planar longitudinal profiles, and convergent, divergent and parallel planform shapes (Fig. 1). Apart from the highly idealised hillslope used in the section ‘‘Effect of grid size on median TI for a hillslope with simple geometry’’ for an analytical derivation (Fig. 1(V)), we improved the degree of physical realism of our hillslopes by using the cosine function c1 cos (c2p/2) to define the shape of their transverse cross-sections. Here, c1 is the maximum depression midway across the widest cross-section and 1 6 c2 6 3 gives the half cycle of the cosine curve (Fig. 2). These more complex hillslopes were used throughout the remainder of our investigation. The mean transverse slope (St, defined in Fig. 2) was set to 10% for parallel hillslopes, and reduced from 10% at the widest cross-section to 4% at the narrowest cross-section for convergent and divergent hillslopes. Although r in Eq. (2) is set to 8 in most real-world applications, we set r = 4 because of the simplicity of the hillslopes considered herein (Fig. 3). Given the use of r = 4 and cosine transverse cross-sections, loss of flow due to edge effects occurs only for one and two directions for divergent and convergent hillslopes, respectively. The loss ranges from less than 0.6% (2.5 m grid) to less than 5% (20 m grid) with an average loss of less than 2.3%. Thus the impacts of edge effects on the outcomes of our research are negligible.

Analyses of hillslopes from a number of catchments in Australia revealed that longitudinal profiles given by a ‘‘bow’’ shape, as used in this paper, are much more prevalent than the ‘‘S’’ shaped profiles with midslope points of contraflexure used by O’Loughlin (1981). Table 1 lists the equations used to calculate the longitudinal profiles of concave and convex hillslopes, and expressions for the longitudinal profile function f(X) defined by S ¼ SfðXÞ

Convergent IV

VII

Convex

I

ð3Þ

where S* = dh/dx is the local slope in which x is distance from the toe of the hillslope and h is the surface elevation at x; S = H/L is the average slope in which H is the overall hillslope height and L is the hillslope length; X = x/L; and f(X) is obtained by differentiating the profile equations given in Table 1 with respect to X. The profile factor B determines the degree of concavity (B > 0) or convexity (B < 0) of the longitudinal profile. The profile equations in Table 1 yield a planar longitudinal profile (Y = h/H = X; f(X) = 1) as B ! 0. The non-dimensional longitudinal profiles described by the profile equations are shown in Aryal et al. (2002, 2003). The planform width is assumed to vary linearly with distance upslope; this shape is described by the convergence ratio CR = br/bo where br and bo are the widths of the hillslope ridge and toe, respectively. CR < 1, CR = 1 and CR > 1 give divergent, parallel and convergent planform shapes, respectively (see Aryal et al., 2002). Although these hillslopes may be somewhat different to those of the real world, they are amenable to experimentation. For our investigation we set: L = 1000 m; longitudinal slope S = 0.2; br = bo = 400 m for parallel planform; br = 200 m and bo = 600 m for divergent planform; and br = 600 m and bo = 200 m for convergent planform. This

Parallel

Divergent

153

V

VIII

Planar

II

VI

IX

Concave

III

Figure 1

Schematic of idealised hillslope types (after R. Suzuki, as discussed by Tsukamoto and Ohta, 1988).

154

S.K. Aryal, B.C. Bates

br

S t = 0.04

bo

S t = 0.1

a .1

St

br

04

St

=

=0

0.

bo

b Figure 2 Schematic of transverse cross-sections for two idealised hillslope types: (a) divergent-planar and (b) convergent-planar (br and bo are the widths of the hillslope ridge and toe, respectively, c1 the amplitude of the cosine curve used to model crosssectional shape, and St the mean transverse slope).

Figure 3

Possible flow directions on an idealised hillslope.

resulted in CR = 0.33, 1 and 3 for divergent, parallel and convergent planform shapes, respectively. Longitudinal profiles given by B = ±2 were used, as well as the planar profile (B  0). These values are based on the hillslope characteristics of the Arkins and Bruthen catchments (section ‘‘Study catchments’’). We used the ANUDEM software to handle elevation data generated using equations in Table 1 to create digital elevation data sets for different grid resolutions (Hutchinson, 1989; http://cres.anu.edu.au/outputs/anudem.php). This software uses a thin plate spline interpolation algorithm.

37

Study catchments 1

38

Melbourne 6 7

3 5

8

4

39

Latitude south

2

Tasman Sea 143

145

147

Longitude east Figure 4 Location of study catchments in southeast Australia (for key to numerals see Table 2).

We considered eight natural catchments located in southeast Australia (Fig. 4), for which digitised elevation contours mapped at the scale of 1:25,000 are available. The vertical separation of the contour lines is 10 m. Details of the catchments are given in Table 2. Catchment shapes range roughly from circular to oblong. For each catchment we compared the TI distributions for different DEM resolutions and confirmed that their differences are due primarily to changes in location rather than scale or shape. Consider Fig. 5, for example, which illustrates the changes in the distribution of TI induced by changes in DEM resolution for the Traralgon catchment. Observe that for both grid and contour DEMs the upper tails of the TI distributions (at which surface saturation is likely to occur) as well as their central parts and lower tails shift to the right as spatial resolution is decreased, and that no cumulative distribution curve crosses another. Further

Effects of catchment discretization on topographic index distributions Table 1

Longitudinal profile equations and functions for concave and convex hillslopes

Longitudinal profile Equations Functions

Table 2

155

Concave (B > 0)

Convex (B < 0)

1 1 ðBÞ Y ¼ tan BðX1Þþtan tan1 ðBÞ B fðXÞ ¼ f1þB2 ðX1Þ2 g tan1 ðBÞ

ðBXÞ Y ¼ tan tan1 ðBÞ

1

fðXÞ ¼ f1þB2 X 2 Bg tan1 ðBÞ

Characteristics of study catchments

No.

Name

Drainage area (km2)

Relief (m)

1 2 3 4 5 6 7 8

Wannon River at Jimmy Creek Jimmy Creek at Jimmy Creek Lardner Creek at Gellibrand Cumberland River at Lorne Arkins Creek West Branch at Wyelangta Traralgon Creek at Koornalla Bruthen Creek at Carrajung Lower Jack River at Jack River

47 22 52 40 4.4 83 16 51

830 820 470 630 180 630 490 630

investigation showed that the 0.97 quantile and higher values of TI correspond to computational elements that lie in the streambed. Thus examination of changes in median TI with grid size provides insight across the range of the TI data.

Results for idealised hillslopes Effect of grid size on median TI for a hillslope with simple geometry Consider a parallel-planar hillslope of finite length with St = 0 (Fig. 1(V)). Let Dn denote the grid size after n (=0, 1, 2, . . .) grid reductions, c ( = 2, 3, 5, 6, 7, . . .) the grid reduction parameter, a the area of the entire hillslope, and b the original grid size (equal to the hillslope width). After n grid reductions, the number of grid cells derived from the original cell is given by c2n and Dn = b/cn. Assume that the outflow from a grid cell passes only to its neighbour of steepest descent. For n = 0, D0 = b and as = a/b = b. Halving the original grid size once c = 2; n = 1 yields 22 = 4 grid cells each with area b2/4, D1 = b/2 and as = b, b, b/2, b/2. Thus median as = 3b/4 < b. Halving the grid size once more (n = 2) yields 24 = 16 grid cells with D2 = b/4 and median as = 5b/8 < 3b/4 < b (Fig. 6). Given that tan b = constant for this hillslope, it can be shown by mathematical induction and simple algebraic manipulation that ( ðcn þ1Þb ðaÞ 2cn tan b median TI ¼ bþD ð4Þ n ðbÞ 2 tan b Here, Eq. (4a) is a nonlinear function of the number of grid cells and there is a horizontal asymptote at median as/tan b = b/2tan b. Eq. (4b) is a linear function of Dn with slope M = (2tan b)1 and intercept b/2tan b. M is an important quantity since it is the rate of change of median TI with grid size. Thus, Eqs. (4)(a,b) describe simple relationships in which median TI decreases as DEM resolution increases, and give the theoretical median TI as b/2 tan b for planar parallel hillslopes with St = 0 .

Effects of DEM resolution and v on median TI Fig. 7 illustrates the variation of median TI with number of computational elements for the nine idealised hillslope types described in the section ‘‘Idealised hillslopes’’ using the grid-based method (v = 1, 5 and 10) and the contourbased method. Four features are apparent: (1) median TI is constant for the contour DEMs for all but the coarsest resolutions; (2) for all grid DEMs, there is a systematic variation in median TI with reductions in grid size and asymptotic behaviour is only apparent, if at all, at high resolutions; (3) median TI values for v = 1 are quite different to those for v = 5, 10; and (4) in seven cases out of nine, median TI values obtained from contour DEMs are noticeably different to those obtained from grid DEMs. For the contour-based method, reductions in flow tube spacing should not theoretically affect as since they reduce both the upslope contributing area and the flow width by the same factor. Slight deviations from this theoretical argument can be seen at very coarse resolutions for convergent-concave, parallel-planar and divergent-concave hillslopes. Since the computational element network in the contour DEMs is created by natural flow lines intersecting the contour lines, and given the three-dimensional curvature (plan, profile and cross-section) of the hillslopes considered, slight irregularities in grid shapes and sizes are inevitable. The effects of these irregularities are pronounced at low DEM resolution leading to largely non-systematic fluctuations in the median value of TI. For the grid DEMs, least-squares linear regressions of median TI on grid size for the nine hillslope types and v = 1, 5, 10 (not shown) explained more than 95% of the variance for all but one case (convergent-concave with v = 10 which accounted for 80%. These results are in keeping with Eq. (4b) despite the use of a multiple flow direction algorithm, the complexity of the hillslope geometries considered and the fact that tan b 5 constant. Table 3 lists the estimates of M for the idealised hillslopes and v = 1, 3, 5, 7, 10. In all but two cases (divergent-convex and divergent-concave), M reduces as the

150 60

0.4

Cumulative Frequency

75

0.90 0.92 0.94 0.96 0.98 1.00

0.6

100

0.2

Cumulative Frequency

0.8

a

S.K. Aryal, B.C. Bates 1.0

156

5e+03

5e+04

5e+05

0.0

as/tanβ

1e+02

5e+02

5e+03

5e+04

5e+05

Figure 6 Schematic of successive grid reduction for c = 2 and n = 0, 1 and 2.

4400 3700

0.8

Cumulative Frequency

0.90 0.92 0.94 0.96 0.98 1.00

21400

0.4

0.6

7300 10200

0.2

Cumulative Frequency

b

1.0

as/tanβ

5e+03

5e+04

5e+05

0.0

as/tanβ

1e+02

5e+02

5e+03 as/tanβ

5e+04

5e+05

Figure 5 Cumulative frequency curves for Traralgon TI data (a) numerals indicate grid size (m) for grid DEMs. and (b) numerals indicate number of computational elements for contour DEMs. Insets show magnified view of the upper tail of the TI distribution. Curves above the horizontal dashed lines in the insets refer to streambed elements.

flow gradually becomes more unidirectional. The intercepts of the regression lines also decrease with increasing v in every case. Therefore, it is very difficult to provide guidance as to what value or range of v should be used to construct DEMs for natural catchments. The effect of grid size on M is also a function of planform and profile shape. The M values for concave hillslopes generally exhibit the smallest sensitivities to variations in planform shape and v, unlike those for convex and planar hillslopes where M is more sensitive to variations in planform and in v. For these hillslope types the ones with convergent planform have the highest M. Similarly for low v, convex hillslopes have the largest M and concave hillslopes the lowest. Finally, recall from the section ‘‘Effect of grid size on median TI for a hillslope with simple geometry’’ that M = (2tan b)-1 when a single

flow direction algorithm is applied to a planar parallel hillslope of finite length with zero transverse slope. Thus M = 2.5 for tan b = S = 0.2. For high v we observe that M  2.5 for the planar parallel hillslope with the cosine transverse cross-section described in the section ‘‘Idealised hillslopes’’ and that M > 2.5 for the other profile-planform combinations considered. Therefore, it is unclear whether the linear relationship between M and (tan b)1 does not hold generally or whether different proportionality constants apply for different hillslope geometries. Further investigation of this issue is postponed until the section ‘‘Relationship between M and median topographic slope’’, where the strength of the relationship between M and (tan b)1 is examined in the more important context of natural catchments.

Effect of v on spatial pattern of TI Fig. 8 illustrates the effects of hillslope geometry and v on the spatial distribution of the TI for four profile-planform combinations with cosine transverse cross-sections and v = 1, 5 and 10. For hillslopes with convex longitudinal profiles (bottom two rows) the TI exhibits patterns that are diffused both longitudinally and transversely. The largest values of TI occur at midslope positions. This implies that convex hillslopes are prone to saturation at those locations giving rise to the hanging saturation phenomenon (Aryal et al., 2002). For hillslopes with concave longitudinal profiles, moderate to high values of TI are concentrated in narrow longitudinal bands. The value of v also manipulates the spatial patterns. For example, high values of v tend to enhance the proportion of transverse flow and diminish the proportion of longitudinal flow. This causes the concentration of high TI in the longitudinal centres of the hillslopes, which is particularly noticeable for the hillslopes with concave longitudinal profiles. One of the artefacts of grid DEMs is that the edges of convergent hillslopes with relatively flat transverse cross-sections can experience concentrations of high TI due to smaller diagonal flow, even when a multidirectional flow algorithm is employed (not shown). Under-

Effects of catchment discretization on topographic index distributions v=1 v=5 v = 10 Contour based

Divergent convex

1200

1200

1100 1000 900 800

1900 1800 1700

Parallel convex

1100 Median as/tanβ (m)

Median as/tanβ (m)

1300

1000

Median as/tanβ (m)

1400

157

v=1 v=5 v = 10 Contour based

900 800 700

700 5000

10000

15000

5000 10000 15000 Number of computational elements

950 850 750

v =10

650

Contour based

900 800 700 600

550

500

0

5000

10000

15000

20000

0

4000

Number of computational elements 1300

1000 900 800 700 600 5000

10000

16000

15000

1200 1100 1000 900 800 700 0

5000 10000 15000 Number of computational elements

1700

20000

Convergent concave

1500

1000

20000

20000

v=1 v=5 v =10 Contour based

1400 1300

20000

1100

500 0

12000

Parallel concave

1200

Median as/tanβ (m)

1100

8000

15000

Convergent planar

Number of computational elements v=1 v = 10 v=5 Contour based

Divergent concave

1200 Median as/tanβ (m)

1600 1500

v=5

1000

10000

Number of computational elements v=1

Parallel planar

1100

5000

20000

900 800 v=1 v=5 v = 10 Contour based

700 600 500

Median as/tanβ (m)

Median as/tanβ (m)

1050

1200

Median as/tanβ (m)

v=1 v=5 v =10 Contour based

Divergent planar

1200 1100 1000 0

0

1150

1600 1500 1400 1300

600

20000

Number of computational elements

Median as/tanβ (m)

0

v=1 v=5 v = 10 Contour based

Convergent convex

1300 1100

v=1 v=5 v = 10

900

Contour based

700

400 500

300 0

Number of computational elements

5000 10000 15000 Number of computational elements

20000

0

5000

10000

15000

20000

Number of computational elements

Figure 7 Effect of grid size on median TI for idealised hillslopes for grid-based method for v = 1, 5 and 10 and contour-based method. Results for v = 3 and 7 not shown for clarity.

Table 3

Estimates of M for idealised hillslopes

Profile shape

v

Planform shape Divergent

Parallel

Convergent

Convex

1 3 5 7 10

12.1 6.30 5.30 4.90 4.20

9.60 5.70 4.70 4.60 3.80

24.0 11.1 8.10 6.80 7.20

Planar

1 3 5 7 10

8.30 5.00 4.40 3.90 3.70

7.80 3.50 4.40 2.10 2.30

15.4 7.10 5.20 4.70 4.10

Concave

1 3 5 7 10

8.20 6.70 6.00 5.60 5.80

6.00 5.40 4.90 4.80 4.40

5.90 4.10 3.80 3.70 3.20

standably, for hillslopes with large mean transverse slopes flow is concentrated midslope forming a longitudinal band of high TI for all v.

Results for natural catchments Effect of DEM resolution on median TI Fig. 9 illustrates the variation of median TI with the number of computational elements for contour and grid DEMs for

eight natural catchments. The maximum number of elements used was constrained by the resolution of the digital elevation data (section ‘‘Study catchments’’). For contour DEMs, median TI is constant for moderate to high resolutions despite the topographic complexity of the natural catchments. This behaviour suggests the existence of a threshold resolution above which there is no added value for the estimation of median TI. For the grid DEMs, the asymptote behaviour of median TI is absent or weak across the range of resolutions considered. This is consistent with the results obtained for the nine idealised hillslopes (Fig. 7). Fig. 10 depicts the relationships between median TI and grid size for the study catchments. The regression lines account for 97–98% of the variance of the data, which are again in keeping with Eq. (4b) (cf. section ‘‘Effects of DEM resolution and v on median TI’’). This is despite the use of the D1 approach to analyse the catchments rather than the single flow direction approach used to derive Eq. (4b), the topographic complexity of the natural catchments and tan b 5 constant. Note also that the estimates of the regression parameters vary widely from catchment to catchment.

Relationship between M and median topographic slope In this investigation we used the DEM routine embedded in the TOPOG model to compute local topographic slope so as to avoid the sensitivity of tan b to the resolution of grid-based DEMs. A plot of M against 1/median(tan b) is shown in Fig. 11. The least-squares regression line accounts for 49% of the variance (significant at the 5.2% level). The

158

S.K. Aryal, B.C. Bates

Figure 8 Distribution of TI = as/tan b for convergent-concave (top row), divergent-concave (second row), convergent-convex (third row) and divergent-convex (bottom row) hillslopes with cosine transverse cross-sections; v = 1, 5 and 10 (left to right) and S = 0.2. Flow direction is from top to bottom. Scale intervals are irregular.

data point for Jimmy Creek is distant from the remainder of the data set and thus may draw the fitted line towards it. Application of Cook’s distance test (Cook and Weisberg, 1982) revealed that the data points for the Jimmy and Arkins catchments exert a large influence on the regression parameters. While no grounds could be found for omitting the data point for Arkins, when compared with its counter-

parts the topography of Jimmy Creek is quite peculiar in that it is very steep in the western half and northeast corner of the catchment and relatively flat in the remainder. The TOPOG-generated element network for Jimmy Creek is shown in Fig. 12. The number of computational elements (and hence the number of tan b values) in the steeper portion of the catchment is very high. This will cause the in-

Effects of catchment discretization on topographic index distributions

4000

flow tube

1000

0

3000 2000 1000

grid

3000 2000 1000

10000

15000

20000

0

0

Number of computational elements

2000

4000

6000

8000

10000

12000

0

Number of computational elements

Cumberland

Median as /tanβ

flow tube grid

2000

1000

5000

10000

15000

grid

1000

0

0

20000

flow tube

2000

0

0

10000

Traralgon

grid

1000

8000

3000

3000

flow tube

0

2000 4000 6000 Number of computational elements

Arkins

2000

Median as/tanβ

5000

flow tube

4000

0 0

500

1000

1500

2000

2500

3000

3500

0

Number of computational elements

Number of computational elements Bruthen

5000 10000 15000 Number of computational elements

20000

25000

Jack

3000

flow tube

2000

flow tube

grid

grid

Median as /tanβ

Median as/tanβ

Median as/tanβ

grid

2000

5000

flow tube

grid

Median as /tanβ

Median as/tanβ

Lardner

Jimmy

Wannon 3000

Median as/tanβ

159

2000

1000

1000

0

0 0

1000

2000

3000

4000

5000

6000

7000

0

5000

10000

15000

Number of computational elements

Number of computational elements

Figure 9

Effect of DEM resolution on median TI for study catchments.

5000 Jimmy Creek y = 25.9x + 179

Median as/tanβ (m)

4000

Lardner y = 23.8x - 501

Wannon y = 16.8x + 530

3000

Arkins y = 32.6x - 403 Bruthen y = 21.9x - 142

2000

Jack River y = 16.4x - 242

Traralgon y = 11.9x - 173

1000 Cumberland y = 15.7x - 74.0

0 0

Figure 10

50

100 Grid size (m)

150

200

250

Effect of grid size on median TI for study catchments.

verse median depicted in Fig. 11 to fall at a lower value than what may be considered to be a characteristic value of (tan b)1. Thus the data point for Jimmy Creek might be considered as an outlier. Application of robust regression methods led to regression lines that were very similar to the least squares fit obtained by removing the data point for Jimmy Creek. For example, the root mean square difference between the fits of the least-squares and least trimmed squares regression lines is 0.26. As shown in Fig. 11, the least-squares regression line for the reduced data set accounts for 86% of the variance (significant at the 0.26% level).

Discussion The results of our experiments with idealised hillslopes complement one of the findings of Holmgren (1994). Holmgren recommends 4 6 v 6 6 based on experiments with flow overlaps in three DEMs with three combinations of horizontal and vertical resolution for a 3 · 4 km area in northern Sweden and a fourth DEM in which the grid cells are filled with random numbers for spot heights. We found that use of a similar range of v gives consistent values of median TI, but only for a given grid size. This finding holds over a wide range of grid sizes and hillslope types. Moreover,

160

S.K. Aryal, B.C. Bates 35 Arkins

30 Jimmy

25

Lardner

Bruthen 2

R = 0.49

M

20

Wannon

Jack Cumberland

15

Traralgon

2

R = 0.86

10 5 0

0.0

1.0

2.0

3.0

4.0

5.0

6.0

7.0

8.0

9.0

1/median(tanβ)

Figure 11 Plot of M versus 1/median(tan b) for study catchments. Solid line is the linear regression line for the full data set. Dashed line is the regression line with Jimmy Creek omitted.

N

2000 m

Figure 12

TOPOG-generated element network for Jimmy Creek.

Holmgren (1994) and Quinn et al. (1995) suggest that v should be increased as grid resolution is increased. Our work indicates that the distribution of TI is sensitive to grid size and that the above strategy will lead to a distribution of TI that is quite different to what would have been obtained otherwise (Fig. 7). Our results suggest the alternative strategy of choosing smaller v for increasing grid resolution (or

vice versa) to partially reduce the effect of the selection of v on the determination of median TI. Kim and Lee (2004) have proposed and demonstrated a spatially distributed flow apportioning algorithm in which the degree of flow apportioning to neighbouring grid cells can be modulated according to its position relative to the hillslope ridge and toe. This approach obviates the need to select a value

Effects of catchment discretization on topographic index distributions of v prior to a DEM analysis, and is worthy of further investigation in the context of assessing its impact on the sensitivity of TI to changes in grid size. Our results have interesting implications for the calibration of topographically-based hydrology models based on grid DEMs. This topic has received some attention in recent years. For example, Franchini et al. (1996) found that the hydraulic conductivity parameter within TOPMODEL increases with grid size, and can take values that are physically infeasible. This phenomenon led Ibbitt and Woods (2004) to suggest ways to compensate objectively for the distorting effects of grid size on saturated hydraulic conductivity so that its measured value can be used in TOPMODEL. Quinn et al. (1995) suggest that if only the (catchment outlet) hydrograph is required then the use of a coarse DEM employing any of the algorithms they describe will suffice, and that any inaccuracy incurred will be compensated by the calibrated parameters (see e.g., Wolock and McCabe, 1995). They also recommend the use of fine resolution DEMs if information about internal state variables is required. We interpret these coping strategies as a symptom of an inherent problem with grid DEMs that has been clearly elucidated in this study: the distribution of the TI data is a strong function of grid size. In our view the internal adjustment of model parameter estimates to compensate for grid size effects undermines the physical credibility of any hydrology model: a better strategy is to use a contour DEM with sufficiently fine resolution (section ‘‘Effect of DEM resolution on median TI’’) whenever practicable. If users continue to use gridbased methods then they must accept the limitation that their model calibrations will be scale dependent. A possible way forward for the grid-based approach would be to always use the finest possible resolution consistent with the accuracy of the original data. This would reduce the error due to grid discretization but not completely remove it.

3.

4.

5.

6.

7.

Conclusions We examined the construction of digital elevation models (DEMs) using grid- and contour-based techniques, studying the effects of these techniques and DEM resolution on the distribution of the topographic index (TI) defined by TI = as/tan b for idealised hillslopes and natural catchments. We focused on the median of TI since changes in the distributions of TI due to variations in DEM resolution were well characterised by changes in location rather than changes in scale or shape. The following conclusions can be drawn from our study: 1. A simple analytical model for a parallel-planar hillslope and a single flow direction algorithm indicates that median TI is asymptotic with respect to number of grid cells (Eq. (4a)) and a linear function of grid size (Eq. (4b)). A linear variation in median TI with grid size was also found for nine idealised hillslopes with more complicated geometries and eight natural catchments. 2. For the natural catchments and their corresponding grid DEMs, the asymptotic behaviour of median TI with respect to the number of grid cells was either absent or weak across the range of resolutions considered. This suggests that either the resolution of the original elevation data (10 m) is too coarse for asymptotic behaviour

161

to be obtained or Eq. (4a) may not be applicable to the complex topography of natural catchments. However our results indicate Eq. (4b), which was obtained for a highly idealised hillslope geometry, is applicable to natural catchments across a wide range of resolutions. In multiple flow direction algorithms variations in the directional parameter v affect estimates of TI and the rate of change (M) of median TI with grid size. The value of M is also dependent on hillslope geometry. These findings make it impractical to recommend appropriate values of v for general use. Thus conventional flow partitioning will remain an open subject and a potential source of model error in hydrologic applications. The effect of DEM resolution on the performance of algorithms that accommodate spatially varying v needs to be investigated. Values of TI derived from contour DEMs are far less sensitive to variations in the number of computational elements than those derived using grid DEMs. Moreover, unlike grid-based methods, contour-based techniques provide consistent estimates of TI at moderate to high DEM resolutions. Our results for natural catchments suggest that grid and contour DEMs do not often provide similar estimates of TI even with resolutions approaching that of the original elevation data. The extent to which this phenomenon can be attributed to the resolution of the original elevation data or differences in interpolation algorithms and their parameter settings needs to be investigated with very-high resolution data sets. For the eight natural catchments considered herein, there is some evidence for a relationship between M and median topographic slope that has a similar functional form to that derived from the simple analytical model. A stronger finding could not be reached due to the presence of one apparent outlier and second data point with high leverage in a small sample setting. Given the above findings, interpretations of saturated areas from grid DEMs are likely to be unreliable. Provided the need for manual intervention is not prohibitive (section ‘‘Grid and contour DEMs’’), contour DEMs should be the method of choice for identifying and interpreting hillslope saturation patterns.

These results are encouraging but should be verified against additional experimental data. For the sake of completeness, and with a view to applying the same DEM techniques to both the idealised hillslopes and the natural catchments, these additional studies should also include the use of triangulated irregular network algorithms. Wilson and Gallant (2000, p. 108) list seven key assumptions and limitations of the topographic index approach and the details will not be repeated here. Perhaps the most important consideration is that use of the TI assumes steady-state conditions and thus our investigation will have limited application in studies of the transient behaviour of catchments.

Acknowledgements We thank David Tarboton and Warwick Dawes for their assistance with the TauDEM and TOPOG softwares, and John

162 Gallant, Rodger Grayson and Ross Woods for their helpful comments and suggestions. The constructive comments of Stephen J. Burges and three anonymous reviewers led to significant improvements to the manuscript.

References Aryal, S.K., O’Loughlin, E.M., Mein, R.G., 2002. A similarity approach to predict landscape saturation in catchments. Water Resources Research 38 (10), 1208, doi:10.1029/2001 WR000864. Aryal, S.K., Mein, R.G., O’Loughlin, E.M., 2003. The concept of effective length in hillslopes: assessing the influence of climate and topography on the contributing areas of catchments. Hydrological Processes 17 (1), 131–151. Band, L.E., Vertessy, R.A., Lammers, R., 1995. The effect of different terrain representations and resolution on simulated watershed processes. Zeitschrift fuer Geomorphologie, Supplementbaende 101, 187–199. Beven, K.J., Kirkby, M.J., 1979. A physically-based variable contributing area model of basin hydrology. Hydrological Sciences Bulletin 24, 43–69. Beven, K.J., Wood, E.F., Sivapalan, M., 1988. On hydrological heterogeneity – catchment morphology and catchment response. Journal of Hydrology 100, 353–375. Bruneau, P., Gascuel-Odoux, C., Robin, P., Merot, Ph., Beven, K., 1995. Sensitivity to space and time resolution of a hydrological model using digital elevation data. Hydrological Processes 9 (1), 69–81. Chirico, G.B., Blo ¨schl, G., Grayson, R., Western, A., Woods, R., 2003. Contour-based and grid-based modelling: a comparison. Geophysical Research Abstract. European Geophysical Society 5, 13370. Chirico, G.B., Western, A., Grayson, R., Blo ¨schl, G., 2005. On the definition of the flow width for calculating specific catchment area patterns from gridded elevation data. Hydrological Processes 19 (13), 2539–2556. Cook, D., Weisberg, S., 1982. Residuals and Influence in Regression. Chapman and Hall, New York. Costa-Cabral, M.C., Burges, S.J., 1994. Digital terrain model networks (DEMON): A model of flow over hillslopes for computation of contributing and dispersal area. Water Resources Research 30 (6), 1681–1692. Desmet, P.J.J., Govers, G., 1996. Comparison of routing algorithms for digital elevation models and their implications for predicting ephemeral gullies. International Journal of Geographical Information Systems 10, 311–333. Endreny, T.A., Wood, E.F., 2003. Maximizing spatial congruence of observed and DEM-delineated overland flow networks. International Journal of Geographical Information Science 17 (7), 699– 713. Fairfield, J., Leymaire, P., 1991. Drainage networks from grid digital elevation models. Water Resources Research 27 (5), 709– 717. Franchini, M., Wendling, J., Obled, C., Todini, E., 1996. Physical interpretation and sensitivity analysis of the TOPMODEL. Journal of Hydrology 175, 293–338. Freeman, T.G., 1991. Calculating catchment area with divergent flow based on a regular grid. Computers & Geosciences 17 (3), 413–422. Freer, J., McDonnell, J., Beven, K.J., Brammer, D., Burns, D., Hooper, R.P., et al, 1997. Topographic controls on subsurface storm flow at the hillslope scale for two hydrologically distinct small catchments. Hydrological Processes 11 (9), 1347–1352. Gallant, J.C., Wilson, J.P., 2000. Primary topographic attributes. In: Wilson, J.P., Gallant, J.C. (Eds.), Terrain Analysis: Principles and Applications. John Wiley and Sons Inc., New York, pp. 51– 85.

S.K. Aryal, B.C. Bates Grayson, R., Blo ¨schl, G., 2001. Spatial modelling of catchment dynamics. In: Spatial Patterns in Catchment Hydrology: Observations and Modelling. Cambridge University Press, Cambridge, UK. Hjerdt, K.N., McDonnell, J.J., Seibert, J., Rodhe, A., 2004. A new topographic index to quantify downslope controls on local drainage. Water Resources Research 40 (5), W05602, doi:10.1029/2004WR003130. Holmgren, P., 1994. Multiple flow direction algorithms for runoff modeling in grid based elevation models: an empirical evaluation. Hydrological Processes 8 (4), 327–334. Hutchinson, M.F., 1989. A new method for gridding elevation and stream line data with automatic removal of spurious pits. Journal of Hydrology 106, 211–232. Ibbitt, R., Woods, R., 2004. Re-scaling the topographic index to improve the representation of physical processes in catchment models. Journal of Hydrology 293, 205–218. Iorgulescu, I., Jordan, J.P., 1994. Validation of TOPMODEL on a small Swiss catchment. Journal of Hydrology 159, 225–273. Kim, S., Lee, S., 2004. A digital elevation analysis: a spatially distributed flow apportioning algorithms. Hydrological Processes 18, 1777–1794. Maunder, C.J., 1999. An automated method of constructing contour based digital elevation models. Water Resources Research 35 (12), 3931–3940. Mizukoshi, H., Aniya, M., 2002. Use of contour-based DEMs for deriving and mapping topographic attributes. Photogrammetric Engineering and Remote Sensing 68 (1), 83–93. Moore, I.D., 1996. Hydrologic modelling and GIS. In: Goodchild, M.F., Steyaert, L.T., Parks, B.O., Johnson, C., Maidment, D., Crane, M., et al. (Eds.), GIS and Environmental Modelling: Progress and Research Issues. GIS World Books, Fort Collins, CO, pp. 143–148. Moore, I.D., Grayson, R.B., 1991. Terrain-based catchment partitioning and runoff prediction using vector elevation data. Water Resources Research 27 (6), 1177–1191. Moore, I.D., Grayson, R.B., Ladson, A.R., 1991. Digital terrain modeling: a review of hydrological, geomorphological and biological applications. Hydrological Processes 5 (1), 3–30. Moore, I.D., Mackay, S.M., Wallbrink, P.J., Burch, G.J., O’Loughlin, E.M., 1986. Hydrologic characteristics and modeling of a small forested catchment in southeastern New South Wales: prelogging condition. Journal of Hydrology 83, 307–335. Moore, R.D., Thomson, J.C., 1996. Are water table variations in a shallow forest soil consistent with the TOPMODEL concept? Water Resources Research 32 (3), 663–669. O’Callaghan, J.F., Mark, D.M., 1984. The extraction of drainage networks from digital elevation data. Computer Vision, Graphics and Image Processing 28 (3), 323–344. O’Loughlin, E.M., 1981. Saturation regions in catchments and their relations to soil and topographic properties. Journal of Hydrology 53, 229–246. O’Loughlin, E.M., 1986. Prediction of surface saturation zones in natural catchments by topographic analysis. Water Resources Research 22 (5), 794–804. O’Loughlin, E.M., 1990. Modelling soil moisture status in complex terrain. Agricultural and Forest Meteorology 50, 23–38. Pan, F., Peters-Lidard, C.D., Sale, M.J., King, A.W., 2004. A comparison of geographical informations systems-based algorithms for computing the TOPMODEL topographic index. Water Resources Research 40 (6), W06303. doi:10.1029/2004WR00306. Quinn, P., Beven, K., Chevallier, P., Planchon, O., 1991. The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models. Hydrological Processes 5 (1), 59–79. Quinn, P.F., Beven, K.J., Lamb, R., 1995. The ln(a/tanb) index: how to calculate it and how to use it within the TOPMODEL framework. Hydrological Processes 9 (2), 161–182.

Effects of catchment discretization on topographic index distributions Scanlon, T.M., Raffensperger, J.P., Hornberger, G.M., Clapp, R.B., 2000. Shallow subsurface stormflow in a forested headwater catchment: Observation and modeling using a modified TOPMODEL. Water Resources Research 36 (9), 2575–2586. Tarboton, D.G., 1997. A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research 33 (2), 309–319. Tsukamoto, Y., Ohta, T., 1988. Runoff process on a steep forested slope. Journal of Hydrology 102, 165–178. Vertessy, R.A., Hatton, T.J., O’Shaughnessy, P.J., Jayasuriya, M.D.A., 1993. Predicting water yield from a Mountain Ash forest catchment using a terrain analysis based catchment model. Journal of Hydrology 150, 665–700. Vivoni, E.R., Ivanov, V.Y., Bras, R.L., Entekhabi, D., 2004. Generation of triangulated irregular network based on hydrologic similarity. Journal of Hydrologic Engineering 9 (4), 288–302. Vivoni, E.R., Ivanov, V.Y., Bras, R.L., Entekhabi, D., 2005. On the effects of triangulated terrain resolution on distributed hydrologic model response. Hydrological Processes 19, 2101–2122. Western, A.W., Grayson, R.B., Blo ¨schl, G., Willgoose, G.R., McMahon, T.A., 1999. Observed spatial organization of soil moisture and its relation to terrain indices. Water Resources Research 35 (3), 797–810. Wilson, J.P., Gallant, J.C., 2000. Secondary topographic attributes. In: Wilson, J.P., Gallant, J.C. (Eds.), Terrain Analysis: Principles and Applications. John Wiley and Sons Inc., New York, pp. 87–131. Wilson, J.P., Repetto, P.L., Snyder, R.D., 2000. Effect of data source, grid resolution, and flow-routing method on computed topographic attributes. In: Wilson, J.P., Gallant, J.C. (Eds.),

163

Terrain Analysis: Principles and Applications. John Wiley and Sons Inc., New York, pp. 133–161. Wise, S., 2000. Assessing the quality for hydrological applications of digital elevation models derived from contours. Hydrological Processes 14 (11–12), 1909–1929. Wolock, D.M., 1995. Effects of subbasin size on topographic characteristics and simulated flow paths in Sleepers River watershed, Vermont. Water Resources Research 31 (8), 1989– 1997. Wolock, D.M., Price, C.V., 1994. Effects of digital elevation model map scale and data resolution on a topography-based watershed model. Water Resources Research 30 (11), 3041–3052. Wolock, D.M., McCabe, G.J., 1995. Comparison of single and multiple flow direction algorithms for computing topographic parameters in TOPMODEL. Water Resources Research 31 (5), 1315–1324. Woods, R.A., Sivapalan, M., 1997. A connection between topographically driven runoff generation and channel network structure. Water Resources Research 33 (12), 2939– 2950. Woods, R.A., Sivapalan, M., Robinson, J.S., 1997. Modeling the spatial variability of subsurface runoff using topographic index. Water Resources Research 33 (5), 1061–1073. Zhang, W., Montgomery, D.R., 1994. Digital elevation model grid size, landscape representation, and hydrologic simulations. Water Resources Research 30 (4), 1019–1028. Zhou, Q., Liu, X., 2002. Error assessment of grid-based flow routing algorithms used in hydrological model. International Journal of Geographical Information Science 16 (8), 819–842.