Spatial statistical properties and scale transform analyses on the topographic index derived from DEMs in China

ARTICLE IN PRESS Computers & Geosciences 35 (2009) 592–602

Contents lists available at ScienceDirect

Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo

Spatial statistical properties and scale transform analyses on the topographic index derived from DEMs in China Bin Yong a,, Wan-Chang Zhang b, Guo-Yue Niu c, Li-Liang Ren a, Cheng-Zhi Qin d a

State Key Laboratory of Hydrology, Water Resources and Hydraulic Engineering, Hohai University, No. 1 Xikang Road, Nanjing 210098, China Key Laboratory of Regional Climate-Environment Research for Temperate East Asia, Institute of Atmospheric Physics, CAS, Beijing 100029, China John A. and Katherine G. Jackson School of Geosciences, Department of Geological Sciences, University of Texas, Austin, TX, USA d State Key Laboratory of Resources and Environmental Information System, Institute of Geographical Sciences and Natural Resources Research, CAS, Beijing 100101, China b c

a r t i c l e i n f o

abstract

Article history: Received 19 September 2007 Received in revised form 27 March 2008 Accepted 31 March 2008

The topographic index (TI), frequently used in approximately characterizing the spatial distribution of variable source areas within a watershed, has been widely applied in topography-based land-surface process schemes coupled in regional or global climatic models. The TI concept, however, was originally developed for studying hill-slope scale hydrological processes and was most commonly used in simulations from small- to medium-sized watersheds. It is still questionable whether the TI computed from coarseresolution digital elevation models (DEMs) for large-scale hydrology and climate studies can effectively reflect the spatial distribution of soil moisture, surface saturation, and runoff generation processes in most areas. In this study, we first proposed an improved multiple flow direction algorithm (IMFD) for accurately computing the TI distribution. We then evaluated the IMFD algorithm quantitatively by using four types of artificial mathematical surfaces. Subsequently, we conducted statistical analyses on the TI distributions computed with IMFD from 90  90 m2 and 1000  1000 m2 resolution DEM blocks sampled from across the whole of China. We found there are linear relationships between the mean TI values computed from the two different resolution DEMs in three sampled blocks of different sizes, i.e., 0.11  0.11, 0.51  0.51 and 11 11. Systematic analyses further suggested that the forms of these linear relationships are evidently affected by the algorithm used for the TI computation, while the size, location, and number of the selected TI samples have minor effects on them. Finally, we investigated the influence of DEM resolution on the spatial statistical properties of TI. From the viewpoint of terrain discretization and smoothing effects, we also discussed the mechanism and the reasons causing the similarity on TI at different spatial resolutions. & 2008 Elsevier Ltd. All rights reserved.

Keywords: Topographic index TOPMODEL DEM RCM GCM Downscaling Multiple flow direction algorithm

1. Introduction Originally defined in TOPMODEL (Beven and Kirkby, 1979) to approximate the likely spatial distribution of variable source areas (Quinn et al., 1995) and predict local variations in water table depths (Brasington and Richards,

 Corresponding author. Tel.: +86 25 83787361; fax: +86 25 83787364.

E-mail address: [email protected] (B. Yong). 0098-3004/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2008.03.006

1998) within a watershed, the topographic index (TI, also called wetness index), has been widely used to study the effects of topography on hydrologic processes at a basin scale (Beven et al., 1984; Wolock, 1993; Ambroise et al., 1996; Beven, 1997). The TOPMODEL concept, owing to its advantages of a simple structure and physically based nature as well as its potential to couple with groundwater variations in time and space, has been frequently implemented to topography-based land-surface process schemes used in regional climate models (RCMs) or global

ARTICLE IN PRESS B. Yong et al. / Computers & Geosciences 35 (2009) 592–602

climate models (GCMs) (Famiglietti and Wood, 1994; Stieglitz et al., 1997; Ducharne et al., 2000; Koster et al., 2000; Chen and Kumar, 2001; Niu and Yang, 2003; Niu et al., 2005). For such large-scale models, directly using high-resolution digital elevation models (DEMs) to compute the spatial distribution of TI is very difficult due to the enormous volume of the data, especially for global applications (Kumar et al., 2000). In addition, it is hard to acquire high-resolution DEM (e.g., 30 m) data for all the land areas of global continents. Therefore, the coarseresolution DEMs, such as GTOPO30 (30 arcsec, about 1000 m resolution) provided by the US Geological Survey (USGS), are often used for the retrieval of the TI distribution at regional or global scales for land-surface hydrological process simulations (Chen and Kumar, 2001; Niu and Yang, 2003; Niu et al., 2005; Nourani and Mano, 2007). It is questionable whether the TI computed from coarse-resolution DEMs can be reasonably used for characterizing the likely spatial distribution of variable source areas and determining runoff generation processes in regional or global climate simulations. To assess the impacts of DEM resolution on the TI computation, a number of studies have been carried out in the past two decades and some valuable results have been achieved (Hutchinson and Dowling, 1991; Jenson, 1991; Quinn et al., 1991; Iorgulescu and Jordan, 1994; Wolock and Price, 1994; Zhang and Montgomery, 1994; Bruneau et al., 1995; Mendicino and Sole, 1997; Saulnier et al., 1997; Brasington and Richards, 1998). However, there is a lack of researches on the possibility of using coarse-resolution DEMs (e.g., GOTOP30, about 1000 m resolution) to obtain the TI distributions with accuracies comparable to those derived from high-resolution DEMs (e.g., STRM DEM3, about 90-m resolution, available from USGS). To answer the above question, we focus on: (1) developing an improved multiple flow direction algorithm (IMFD) for more accurate computation of TI, which will be evaluated with an objective assessment approach, (2) investigating the statistical properties of TI computed from different resolution DEMs, (3) establishing a simple but operational downscaling scheme for the TI applications with different spatial resolutions in China, and (4) examining the impacts of DEM resolution on the retrieved spatial characteristics of the TI distribution and exploring the mechanism and the reasons causing the similarity on TI at different spatial resolutions. This study will be beneficial not only for application of TOPMODEL concept to large-scale hydrological models and landsurface schemes for climate studies, but also for understanding the fractal properties of the TI fields at different DEM resolutions.

2. The IMFD algorithm used for the TI computation The TI is defined as ln(a/tan b), where a is the accumulative upslope area per unit contour length (or specific catchment area) and b is the local slope angle. Current GIS-based computation methods of TI mainly include algorithms of single flow direction (SFD), bi-flow direction (BFD), and multiple flow direction (MFD)

593

(Pan et al., 2004). The MFD algorithm, originally proposed by Quinn et al. (1991) (referred to as ‘MFD-Quinn’ in this paper), is recognized as the most accurate (Wolock and McCabe, 1995; Giuseppe and Aurelia, 1997; Pan et al., 2004). On the basis of MFD-Quinn, we develop an IMFD for the accurate computation of the TI distribution. 2.1. The multiple flow direction algorithm (MFD-Quinn) Different from SFD, the MFD algorithm assumes that flow from a current grid cell could drain into more than one downslope neighboring grid cells. The fraction of the area draining through each current cell to each downslope direction is proportional to the slope gradient of each downhill flow path, so that steeper gradients will attract more water from the upslope contributing area (Quinn et al., 1991). In MFD-Quinn, the calculation of accumulative upslope area (i.e., a) can be concisely summarized as following four important equations (i.e., Eqs. (1)–(4)). According to Quinn et al. (1991), a was defined as

a ¼ A=

n X

Lj

(1)

j¼1

where A is the total area draining into each cell, Lj is the effective contour length of the cell boundary between current cell and its jth downslope neighboring cell. The value of Lj is set to a half of grid size for cells in cardinal directions as illustrated as L1 and L3 (Fig. 1), and 0.354 for cells in diagonal directions, such as L2 and L4 (Fig. 1). The local slope angle in the downslope direction, tan b, is computed as tan b ¼

n n X X ðtan bj Lj Þ= Lj j¼1

(2)

j¼1

where tan bj is the slope gradient between the current cell and its jth downslope neighboring cell. Combining Eqs. (1) and (2), we can obtain the ln(a/tan b) expression 0 1 n X @ ðtan bj Lj ÞA (3) lnða= tan bÞ ¼ ln A= j¼1

K1 10 L4 8

13

12

11

K2

K3

10

L1

L3

L2

7

9

8

Fig. 1. An example of a simulated flow partition on DEM using multiple flow direction algorithm.

ARTICLE IN PRESS 594

B. Yong et al. / Computers & Geosciences 35 (2009) 592–602

The total upslope area of the current cell is distributed to its downslope neighboring cells as 0 1 n X DAi ¼ Aðtan bi Li Þ=@ ðtan bj Lj ÞA (4) j¼1

where DAi is the drainage area passed from the current cell onto its downslope neighboring cell i. 2.2. The improved multiple flow direction algorithm (IMFD) In MFD-Quinn, a is the effective cumulative area draining from the upslope neighboring cells to the current cell, which reflects the tendency of upslope water to accumulate at the current cell in the catchment. While tan b is the local slope angle of the current cell, which reflects the tendency for gravitational forces to move upslope water from the current cell to downslope neighboring cells. Eqs. (1) and (2) suggest that the effective contour lengths of the grid-cell boundaries between current cell and its downslope neighboring cells (e.g., L1, L2, L3, and L4 in Fig. 1) are used to determine both a and tan b. However, in IMFD, we consider that the contour lengths used to determine a should be the effective contour lengths of the grid-cell boundaries between upslope cells and current cell (e.g., K1, K2, and K3 in Fig. 1) rather than that boundaries between current cell and downslope cells. Thus, Eq. (1) should be amended as following:

a ¼ A=

m X

Ki

(5)

i¼1

where Ki is the effective contour length of the grid-cell boundary between the current cell and its ith upslope neighboring cell. Using Eq. (5) to substitute (1), Eq. (3) then becomes 0 1 n X lnða= tan bÞ ¼ ln@ A= ðtan bj Lj ÞA j¼1

0 1 n m X X þ ln@ Lj = K iA j¼1

complex terrain conditions of real watershed include both convergence and divergence because the exponent cannot be altered in response to local terrain conditions. In this study, we adopt a simple but adaptive scheme proposed by Qin et al. (2007) to dynamically determine the downslope flow-partition exponent and accurately compute the TI distribution under various complex terrain conditions. According to Qin et al.’s (2007) approach, Eq. (7) can be substituted to the following expression for computing DAi: 0 1 n X f ðeÞ f ðeÞ 2 (8) DAi ¼ ðA þ Dx Þððtan bi Þ Li Þ=@ ðtan bj Þ Lj A j¼1

Under MFD-Quinn, which sets f(e) ¼ 1, the partition of downslope flow cannot adapt to varying local terrain conditions. Here, e is the tangent value of the maximum downslope angle; f(e) is a linear function of e, rather than a fixed constant value. More details about how to construct the function f(e) are in Qin et al. (2007). Because the ‘real-world’ DEMs often contain depressions or flat areas where the TI cannot be directly calculated, a preprocessor procedure for raw DEM data is necessary for various TI computation algorithms. IMFD uses the simple but operative method of Jenson and Domingue (1988) to deal with depressions in DEM. With regard to the computation of slopes in flat areas, Wolock and McCabe (1995) suggested that the slope gradient of current flat cell equals to (0.5  vertical resolution)/ (horizontal resolution). This method is typical and simple enough, but not too reasonable in some very low relief areas. In IMFD, we use the tracking flow direction (TFD) method recommended by Pan et al. (2004) to process the flat cells in DEM. It has been tested that TFD is more accurate than Wolock and McCabe’s method under most terrain conditions because TFD can produce smaller slope value of flat area. IMFD should compute the spatial distribution of TI more accurately than MFD-Quinn, because the above improved equations and approaches can compensate for some serious deficiencies in MFD-Quinn.

(6)

i¼1

Additionally, the calculated total drainage area passed onto the downslope cells does not include the area of current cell in MFD-Quinn. This shortcoming will bring about some errors in the iterative programming for the computation of TI. Thus, it is necessary to modify Eq. (4) as 0 1 n X 2 @ DAi ¼ ðA þ Dx Þðtan bi Li Þ= ðtan bj Lj ÞA (7) j¼1

where Dx is the grid size of DEM, and the definition of other symbols is the same as that in Eq. (4). Although some computing deficiencies of the effective cumulative area have been remedied through the above amendments, another important problem still remains unresolved, i.e., the fixed-exponent strategy to model the flow-partition proportion used in MFD-Quinn. In practice, the fixed-exponent strategy is unsuitable when the

2.3. Quantitative evaluation on IMFD To further investigate the accuracy and detect the effectiveness of IMFD, we should quantitatively evaluate and compare the errors generated by MFD-Quinn and IMFD. The most commonly used method of evaluating error is to apply various algorithms for measuring on ‘real-world’ DEMs. However, such an evaluation method frequently takes ‘common sense’ as the standard and cannot provide the ‘true’ value of the topographic parameters, so that it is uncertain and subjective to judge ‘reasonable’ or ‘unreasonable’ for proposed algorithms (Zhou and Liu, 2002). Furthermore, we should not neglect the fact that no ‘real-world’ DEM is perfect because of the errors of DEM itself such as sample errors, interpolation errors, and representation errors. Thus, the errors produced by the algorithm could be masked by the DEM errors (Zhou and Liu, 2002). Finally, the results derived from some selected study areas may be unsuitable

ARTICLE IN PRESS B. Yong et al. / Computers & Geosciences 35 (2009) 592–602

elsewhere. So using a ‘real-world’ DEM to compare algorithms is data dependent and difficult to quantify (Qin et al., 2007). While using a selection of mathematical models with controlled parameters, digital terrain analysis (DTA) algorithms can be objectively compared and evaluated, independently from data and human analyst’s bias (Zhou and Liu, 2004). For all these reasons, we adopt some standard and artificial mathematical surfaces to compare different algorithms of TI quantitatively. In this paper, we use the evaluation method developed by Zhou

595

and Liu (2002) to assess IMFD. The kernel of Zhou and Liu’s (2002) approach is to construct four types of artificial surfaces (i.e., ellipsoid, inverse ellipsoid, saddle and plane), which can be represented by mathematical models, thus the theoretical ‘true’ value of a specific catchment area (i.e., a in this paper) on any point of each surface can be inferred from mathematical formulations. Comparing the theoretical values and the computed values, the errors at each grid cell caused by different algorithms will be estimated easily. This evaluation

Fig. 2. Contour maps of artificial surfaces, theoretical a distributions and theoretical ln(a/tan b) distributions of four typical mathematical surfaces: 2 2 2 ðx400Þ2 ðx400Þ2 Z2 þ ðy300Þ þ 20 þ ðy300Þ þ ðZ300Þ ¼ 1 (z4300);0oxo800;0oyo600, 2 ¼ 1 (z40);0oxo800;0oyo600, (b) inverse ellipsoid: 4002 3002 4002 3002 3002 ðx400Þ2 ðy400Þ2 z þ ¼ 0oxo800; 0oyo800 and (d) plane: z ¼ 0.5x+0.25y+1000ox500;0oxyox500; slope ¼ 27.91. 2 2 0:002 2 1

(a) ellipsoid: (c) saddle:

ARTICLE IN PRESS 596

B. Yong et al. / Computers & Geosciences 35 (2009) 592–602

method can simulate various terrain conditions including convergent, divergent, ridge, planar, and so on. Furthermore, its four types of artificial DEM surfaces could be regarded as a good surrogate for real topography. The approach is described in more detail in Zhou and Liu (2002). Fig. 2 gives an example of the four mathematical models. Based on the artificial DEM surfaces and the theoretical a and ln(a/tan b) distributions (see Fig. 2), we can compute the root-mean-square errors (RMSEs) of a and ln(a/tan b) for IMFD and MFD-Quinn, respectively. The RMSE of a (RMSEa) can be computed as following: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pn 2 i¼1 ðaTi  aCi Þ RMSEa ¼ n

(9)

where subscript Ti and Ci stand for theoretical and computed a at grid i, and n is the number of grid cells for evaluation. The RMSE of the calculated TI (RMSETI) is sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pn 2 i¼1 ½lnða= tan bÞTi  lnða= tan bÞCi  RMSETI ¼ n

(10)

Table 1 lists RMSEa and RMSETI computed for the four mathematical surfaces using MFD-Quinn and IMFD. The RMSEs of a show that IMFD yields lower errors than MFD-Quinn under all tested artificial terrain conditions. While for the RMSEs of ln(a/tan b), except for the inverse ellipsoid surface, IMFD produces greater accuracy for other three kinds of terrain conditions than MFD-Quinn does. We conclude that IMFD performs reasonably well and produces better results in most cases. In summary, the quantitative evaluation suggests that the specific catchment area values and the TI values computed by IMFD are closer to the theoretical ‘true’ values than those computed by MFD-Quinn. Compared with MFD-Quinn, IMFD can more accurately express the spatial distribution of the depth to the water table and more reasonably reflect the hydrological similarity within a watershed. In the following analyses and discussions, the IMFD algorithm will be used to compute the spatial distribution of TI for all the selected ‘real-world’ DEMs rather than SFD as adopted by previous researches (Kumar et al., 2000; Wolock and McCabe, 2000).

3. Data description and sampling strategies Table 1 RMSEs for a and ln(a/tan b) computed from four mathematical surfaces by using MFD-Quinn and IMFD Surfaces

Algorithms RMSETI

RMSEa

Ellipsoid Inverse ellipsoid Saddle Plane

60°E

MFD-Quinn

IMFD

MFD-Quinn

IMFD

6.602 500.513 109.779 412.448

5.892 394.534 102.817 32.564

0.063 0.081 0.573 0.176

0.048 0.096 0.459 0.158

90°E 100°E

110°E

70°E

80°E

To perform statistical analyses on the TI distributions computed from different resolution DEMs, two datasets of DEM covering the whole of China were considered: STRM DEM3 (90 m resolution) and GTOPO30 (1000 m resolution) available for public use released by USGS. These two datasets were firstly geo-referenced and projected from a geographical coordinate system to the Albers Equal Area coordinate system. A total of 100 DEM blocks covering a 11 11 latitude–longitude area (grid-cell size of GCMs) were randomly selected from both STRM DEM3 and GTOPO30 datasets (see Fig. 3). The 100 blocks from 90 and 1000 m resolution DEMs representative of the general

120°E

130°E

140°E

40°N

1°

1°

DEM (m) 1421 669

30°N 0.5°

0.5° DEM (m) 1232 732

20°N 0.1°

0.1° DEM (m) 1206 801

Fig. 3. Geographic locations of DEM blocks in three different sizes (11 11, 0.51  0.51, 0.11  0.11) sampled over two different resolution DEMs (90 and 1000 m) all over China as well as sampling strategies used in this study.

ARTICLE IN PRESS B. Yong et al. / Computers & Geosciences 35 (2009) 592–602

parameter used for calculating surface runoff and characterizing the water table regimes of groundwater is the mean value of TI computed in each fundamental grid cell of the models (Chen and Kumar, 2001; Niu and Yang, 2003; Niu et al., 2005). The most representative methods for parameterization of the TI regimes in land-surface hydrological process simulations are a three-parameter gamma distribution function (Stieglitz et al., 1997; Ducharne et al., 2000; Chen and Kumar, 2001; Niu and Yang, 2003) and a two-parameter exponential function (Niu et al., 2005). In these two methods, the mean TI value in each grid-cell plays an important role in connecting the actual spatial distribution of TI with the mathematical expression of TI in regional- and global-scale simulations. Considering the importance of the mean TI value for hydrological processes represented in RCMs and GCMs, we focus on establishing the transformation functions between different resolution DEMs for various typical grid-cell sizes of RCMs and GCMs. For each of the 100 11 11 DEM blocks covering the whole of China, the mean TI values were computed with the IMFD algorithm from 90 and 1000 m resolution DEM datasets. The scatter diagram in Fig. 4(a) suggests a linear relationship between the mean TI values obtained from 90 m resolution DEMs (l90) and those derived from the 1000-m resolution DEMs (l1000) in the case of 11 11 block size. We also find very similar correlation relationships between the mean TI

terrain complexity were subjectively selected for the TI computations. To examine the relationships between the TIs computed from the DEM blocks in different sizes and locations, we clipped other 100 0.51  0.51 DEM blocks (grid-cell size of RCMs) from the center of each 11 11 DEM block as shown in Fig. 3. Following the same procedure, we take the same amount 0.11  0.11 DEM blocks from the center part of each 0.51  0.51 blocks (see the right part of Fig. 3). 4. Spatial scaling analyses on TI We systematically analyzed the statistic relationships between the mean TI values derived from the six DEM datasets (two DEM datasets with spatial resolution of 90 and 1000 m, respectively, at three different sizes: 11 11, 0.51  0.51, 0.11  0.11) computed with the IMFD algorithm. In addition, we carefully investigated the transformation functions of TI between different resolutions for different block sizes of DEM. 4.1. Spatial statistical properties and scale transformation functions In topography-based land-surface process schemes implemented in RCMs or GCMs, the most important

11

y = 0.519x + 2.579

y = 0.518x + 2.588

R2 = 0.762

mean value of TI (90-m)

mean value of TI (90-m)

11

10

9

8

R2 = 0.754

10

9

8

7

7 9

11

10

11 12 13 mean value of TI (1000-m)

14

9

15

R2 = 0.727

10

10

11 12 13 mean value of TI (1000-m)

14

15

11

y = 0.586x + 1.887 mean value of TI (90-m)

mean value of TI (90-m)

597

9 8

10 9 IMFD(1°×1°) IMFD(0.5°×0.5°) IMFD(0.1°×0.1°) SFD(Wolock,2000) SFD(Kumar,2000)

8 7 6

7 9

10

11 12 13 14 mean value of TI (1000-m)

15

9

10

11 12 13 14 mean value of TI (1000-m)

15

Fig. 4. Scatter plots showing fitted linear relationships between the mean TI values computed from two different resolution DEMs (90 and 1000 m) over China. (a) 11 11 grid-cell size, (b) 0.51  0.51 grid-cell size, (c) 0.11  0.11 grid-cell size, and (d) the fitted lines obtained from different algorithms.

ARTICLE IN PRESS 598

B. Yong et al. / Computers & Geosciences 35 (2009) 592–602

values computed from the two different resolution DEMs for the cases of 0.51  0.51 and 0.11  0.11 block sizes as shown in Fig. 4(b) and (c), respectively. The simple linear relationships between the two DEMs of different resolutions for different block sizes differ very slightly. Statistic analyses on the transformation functions between the mean TI values derived from the different resolution DEMs in different locations and sampling sizes indicate that the block size and location has minor effects on the linear transformation relationships between different resolution DEMs. If we regard the sampled DEM blocks as sub-basins, our conclusion is consistent with the finding of Wolock (1995) that the mean TI values are almost constant for sub-basins in the case that their sizes are greater than 5 km2. Comparing to Wolock and McCabe (2000) and Kumar et al. (2000), we find that the algorithm used for the TI computation distinctly affects the linear downscaling equations between different resolution DEMs (see Fig. 4(d)). The slopes of fitted lines computed with the IMFD algorithm are lower than those with the SFD algorithm, while the intercepts tend to be larger. Statistical analyses were performed on the mean TI values computed from the two resolution DEMs (90 and 1000 m) over the 300 selected DEM blocks of three different sizes (11 11, 0.51  0.51, and 0.11  0.11) evenly distributed in the whole of China, a linear downscaling equation, l90 ¼ 0.543l1000+2.332, was obtained. Since we have demonstrated that the size and location of the DEM blocks only slightly affect the transformation functions,

Table 2 Simple linear downscaling scheme on TI computed with IMFD all over China Grid size

Linear downscaling equation

R2

11 11 0.51  0.51 0.11  0.11 Other sizes

l90 ¼ 0.519l1000+2.579 l90 ¼ 0.518l1000+2.588 l90 ¼ 0.586l1000+1.887 l90 ¼ 0.543l1000+2.332

0.762 0.754 0.727 0.738

60°E

70°E

80°E

therefore, this equation may be assumed to represent the downscaling algorithm for a wide range of sizes of simulated grid-cell conditions. The simple linear downscaling equations on the TI regimes obtained in this study for China are summarized in Table 2. Previous researches by Wolock et al. (1990) suggested that an increase in the mean TI value generally leads to an increase in the ratio of overland flow to the total stream flow with a decrease in subsurface flow. Incorporating the terrain information at a higher resolution (90 m), these simple linear downscaling equations on the TI regimes (Table 2) will facilitate coupling TOPMODEL to GCMs and RCMs to improve runoff production and the partitioning of runoff between surface and subsurface runoff components. Thus, this study will be beneficial for producing the same accurate simulations by using coarse-resolution DEMs as using finer-resolution DEMs for large-scale hydrological models and providing a way for more accurate terrain analyses for land-surface process modeling in China.

4.2. The validity of the scale transformation relationships in local areas For detecting the possible effects of the number, location and sampling densities of the DEM blocks on the transformation relationships obtained, we clipped 50 0.51  0.51 DEM blocks covering 301401N–10011101E from the central part of China (Fig. 5) to compare with those computed from 100 0.51  0.51 blocks all over the whole of China. Using the same method and procedure as those used for the DEM blocks selected from across the whole of China, we obtain a downscaling equation, l90 ¼ 0.518l1000+2.488, with R2 being about 0.916 for the central part of China (Fig. 6(a)). Comparing with those computed from 100 0.51  0.51 DEM blocks all over China, we find a minor difference (see Fig. 6(b)). From Fig. 6(b), it is evident that the number, location, and sampling density of selected blocks do not have much

90°E 100°E 110°E 120°E 130°E 140°E

40°N Loess Plateau 30°N Weihe Plain 20°N

Fig. 5. Geographic locations of 50 0.51  0.51 DEM blocks chosen from central China.

ARTICLE IN PRESS

11

y = 0.518x + 2.488 R2 = 0.916

mean value of TI (90-m)

mean value of TI (90-m)

B. Yong et al. / Computers & Geosciences 35 (2009) 592–602

10 9 8 7 9

10 11 12 13 14 mean value of TI (1000-m)

15

599

10 9 8 50 blocks( 0.5°×0.5°) 100 blocks( 0.5°×0.5°) 100 blocks( 1°×1°) 100 blocks( 0.1°×0.1°)

7 6

9

13 10 11 12 14 mean value of TI (1000-m)

15

Fig. 6. (a) Scatter diagram showing a linear regression line for mean TI values of 50 0.51  0.51 DEM blocks clipped from central China and (b) fitted lines obtained from different number or size DEM blocks.

effect on the transformation relationships for TI. In other words, the simple linear downscaling scheme proposed in this paper is applicable to the whole of China. 4.3. Analyses of similarity on TI at different resolutions Although a number of previous studies (e.g., Quinn et al., 1991; Wolock and Price, 1994; Zhang and Montgomery, 1994) have shown that the spatial distribution of TI is strongly dependent on the DEM resolution, they did not further explore the origins and nature of such scaledependence. Recently, Kumar et al. (2000) suggested that the simple linear relationship between the mean TI values derived from different resolution DEMs might be a result of similarity in topographic characteristics at different resolutions. Gallant and Hutchinson (1996) and Wolock and McCabe (2000) indicated that the effects of DEM resolution on topographic characteristics can be separated into two components: terrain discretization and terrain smoothing, which work together to contribute to such topographic similarity. The terrain discretization effect is induced by dividing the same terrain into different numbers of grid cells. While terrain smoothing is caused by using a coarse-resolution DEM, which results in a loss of detailed topographic information. By comparing the slope, specific catchment area, and TI values derived from two different resolution (100 and 1000 m) DEMs for 50 locations in the coterminous USA, Wolock and McCabe (2000) concluded that terrain discretizaiton, and not terrain smoothing, is the primary mechanism by which DEM resolution affects the computation of topographic characteristics. 4.3.1. Terrain discretization and terrain smoothing effects For analyzing the terrain discretization and smoothing effects of DEM resolution on TI, three 0.51  0.51 DEM blocks located at Qilian Mountain, Loess Plateau, and Weihe Plain, representing three typical terrain types of cliffy, undulating, and flat respectively, were selected from the 90 and 1000 m resolution DEM datasets from the central part of China (Fig. 5). To discriminate the exact contribution of terrain discretization and terrain smoothing on topographic similarity, we bi-linearly interpolated the 1000 m resolution DEMs to a new 90 m resolution

DEMs to isolate these two influential factors. Thus, the 1000 m and new 90 m DEMs have the same amount of terrain information but different terrain discretization, while the new 90 m DEMs have the same degree of terrain discretization as the original 90 m DEMs but less terrain information. From Fig. 7, it can be seen clearly that the terrain characteristics emerged on the spatial distributions of TI in three different terrain types derived from original 90 m resolution DEMs is the most detailed. Although the new 90 m DEMs seem consistent with the original 90 m DEMs in spatial resolution, their topographic features are quite different because the terrain smoothing effect causes terrain information losses in the new 90 m DEMs. With respect to the spatial distribution of TI computed from the same resolution DEM, the cliffy Qilian Mountain with high peaks and low valleys exhibits much terrain information than those relatively flat areas, such as Loess Plateau and Weihe Plain. The discrepancy between the 1000 and new 90 m TI regime is induced only by the discretization effect because of the same terrain information they have. Additionally, it is very clear that the terrain discretization effect of DEM resolution is more remarkable on relatively flat terrain (see the charts of 1000 m minus new 90 m in Fig. 7). For example, the TI differences of the flat Weihe Plain are relatively large. Accordingly, the TI differences between new 90 m and original 90 m DEMs caused by terrain smoothing are more apparent on the relatively steep terrain of Qilian Mountain (see the last row of Fig. 7). From the above discussions, we conclude that both terrain discretization and terrain smoothing effects are closely related to the local terrain conditions. For relatively cliffy terrain areas, terrain smoothing is mainly responsible for the losses of terrain information, so the smoothing effect becomes more dominant. While for relatively flat terrain, in contrast, terrain discretization becomes a major factor in governing the effect of DEM resolution on the TI distribution, because the effect of discretization on relatively flat terrain areas is stronger than that of terrain information losses. Taking the deviation of the mean TI values between each 1000 m DEM block and new 90 m DEM block as a measure of discretization degree for characterizing the effect of terrain discretization on the original corresponding grid, we find that the deviation is about 1.114 (Table 3)

ARTICLE IN PRESS 600

B. Yong et al. / Computers & Geosciences 35 (2009) 592–602

Weihe Plain (flat)

TI value

TI value

new90-m minus 90-m

1,000-m minus new90-m

new 90-m

90-m

Loess Plateau (middle)

1,000-m

Qilian Mountain (steep)

3- 5 5- 7 7- 9 9 - 11 11 - 13 13 - 15 15 - 17 17 - 19

3- 5 5- 7 7- 9 9 - 11 11 - 13 13 - 15 15 - 17 17 - 19

TI value 3- 5 5- 7 7- 9 9 - 11 11 - 13 13 - 15 15 - 17 17 - 19

TI value below -2 -2 - -1 -1 - 0 0 - 1 1 - 2 2 - 3 3 - 4 above 4

TI value below -2 -2 - -1 -1 - 0 0 - 1 1 - 2 2 - 3 3 - 4 above 4

Fig. 7. Comparison of TI spatial distributions computed from three DEMs in three different topographic conditions. (a) 0.51  0.51 DEM block located at Qilian Mountain with cliffy terrain, (b) 0.51  0.51 DEM block sampled from Loess Plateau with undulating terrain, (c) 0.51  0.51 DEM block taken from relatively flat Weihe Plain. (90 m showing TI distribution computed from the original 90 m DEM of STRM DEM3 datasets, 1000 m computed from 1000 m DEM of GOTOP30 datasets, and new 90 m from 90 m DEM bi-linearly interpolated from 1000 m DEM.)

for all 50 0.51  0.51 mesh grids chosen in Fig. 5. Similarly, the smoothing effect can be defined as the deviation of the mean TI values between new 90 m and original 90 m, and it was calculated as 1.437 (Table 3). The sum of the discretization degree and the smoothing degree can be expressed as the deviation of the mean TI values between the 1000 m DEM blocks and the original 90 m DEM blocks

as 2.551 given in Table 3. However, using the same method as that proposed by Wolock and McCabe (2000), we get a different result that the discretization component only reaches about 43.63% of the total effects. In other words, the effects of DEM resolution on TI are almost evenly shared with discretization and smoothing in central China, which is differed from 77.3% shares for

ARTICLE IN PRESS B. Yong et al. / Computers & Geosciences 35 (2009) 592–602

90 m DEMs can indicate the terrain smoothing effect. The scatter plots as shown in Fig. 8(a) suggest that a rather good linear relationship exists between the mean TI values computed from the 1000 m DEMs and those obtained from the new 90 m DEMs, which implies that the TI variation at different resolutions caused solely by the terrain discretization effect still follows the linear relationship. Similarly, a linear relationship was also found in the terrain smoothing effect (Fig. 8(b)). From this study, we confirm that terrain discretization and terrain smoothing follow the linear regularity fairly well so that we can conclude that the general effects of DEM resolution on the TI distribution also keep the linear relationship. This is most likely the primary reason to explain the similarity on topographic characteristics at different DEM resolutions.

terrain discretization derived by Wolock and McCabe (2000) in the coterminous USA. The most important reason for causing such obvious differences is attributed to the discrepancy of the terrain type and the sampling location of DEM blocks selected. In this study, 37 out of 50 blocks were chosen in the mountainous areas with undulating terrain conditions, and only 13 blocks were sampled from the flat plain areas. In general, the area selected for the present study is clearly more undulated than the whole of coterminous USA. Therefore, such relatively cliffy landforms will certainly result in much topographic details being lost with the decrease of DEM spatial resolutions. Consequently, the proportion of terrain smoothing effect is predominant in the case of the landforms in central China.

4.3.2. The reason of similarity on TI at different DEM resolutions As the analyses and discussions presented in the above sections, two mechanisms, namely terrain discretization and terrain smoothing, govern the impact of DEM resolution on the spatial distribution of TI. For further investigating the cause of similarity in topographic characteristics at different resolutions, we conducted correlation analyses to examine the relationships between the mean TI values of the 1000, 90, and new 90 m DEMs within the 50 0.51  0.51 blocks selected in central China. The correlation between the mean TI values of the 1000 m and new 90 m DEMs can reflect the effect of terrain discretization on the TI variations at such two resolutions, while the relationship between the new 90 m and original

5. Summary and conclusions Based on the MFD-Quinn algorithm (Quinn et al., 1991), we proposed an improved scheme (i.e., IMFD) for more accurate computation of the TI distribution. The improvements mainly include: (1) amending the equations of the classical MFD algorithm for reasonably computing the specific catchment area a, (2) adopting a varying exponent strategy for accurately modeling flow partition, and (3) using a flow direction tracking method (i.e., TFD) to meticulously deal with flat grid cells. We quantitatively evaluated IMFD for four types of mathematical surfaces against their theoretical values of the specific catchment area and the TI. Assessment of the results indicates that the spatial distribution of TI computed by IMFD is more accurate than the result of MFD-Quinn. We systematically analyzed the spatial statistical properties of TI computed with IMFD over the whole of China and then determined a simple linear relationship between the mean TI values derived from different resolution DEMs. A crucial finding of this study is that, on the basis of the IMFD algorithm and with the coarseresolution DEM data usually employed in GCMs or RCMs simulations, more accurate mean values of TI can still be obtained by using a simple linear downscaling scheme. We found that the algorithm used for the TI computation

Median

Maximum

90 m 1000 m New 90 m 1000 m minus new 90 m New 90 m minus 90 m 1000 m minus 90 m

4.373 7.218 4.724

8.291 10.842 9.728 1.114 1.437 2.551

16.771 18.523 16.914

mean value of TI (new90-m)

Minimum

12

y = 0.547x + 3.702 R2 = 0.895

11 10 9 8 9

10 11 12 13 14 mean value of TI (1000-m)

15

mean value of TI (new90-m)

Table 3 Minimum, median, and maximum values of mean TI computed from three different DEMs (90, 1000, new 90 m) for 50 mesh grids chosen from central China and median differential values of TI between every two different DEMs DEM

601

12

y = 1.002x + 1.526 R2 = 0.879

11 10 9 8

7

8 9 mean value of TI (90-m)

10

Fig. 8. Scatter plots showing a relationship between mean TI values computed from different DEMs. (a) 1000 m and new 90 m DEMs, (b) original 90 m and new 90 m DEMs (50 0.51  0.51 DEM blocks sampled from central China).

ARTICLE IN PRESS 602

B. Yong et al. / Computers & Geosciences 35 (2009) 592–602

importantly affects the linear downscaling equations, while the size, number, and location of selected DEM blocks have a minor influence. The downscaling scheme proposed in this study has been shown to be applicable to the whole of China and demonstrated to have an extensive validity. The simple linear downscaling scheme of TI between the two resolutions (90 and 1000 m) will facilitate applications of TOPMODEL concept to hydrology and climate studies at regional and global scales. Another important finding of this work is that the effects of DEM resolution on topographic characteristics can be attributed to terrain discretization and terrain smoothing. In addition, the terrain discretization and smoothing effects are closely related to the local terrain conditions. The statistical analyses suggest that both terrain discretization and terrain smoothing follow a linear law at different DEM resolutions, and such linear properties can primarily explain the similarity of topographic characteristics at different resolutions. Acknowledgements This study is financially supported by the National Key Basic Research Program of China (Grant No. 2006CB400502). Also the paper is a part of results supported from the 111 Project under Grant B08048, Ministry of Education and State Administration of Foreign Experts Affairs, P.R. China. The authors would like to thank two anonymous reviewers for their comments on an earlier version of this paper. References Ambroise, B., Freer, J., Beven, K.J., 1996. Towards a generalization of the TOPMODEL concepts: topographic indices of hydrological similarity. Water Resource Research 32, 2135–2145. Beven, K.J., 1997. TOPMODEL: a critique. Hydrological Processes 11, 1069–1085. Beven, K.J., Kirkby, M.J., 1979. A physically based, variable contributing area model of basin hydrology. Hydrological Science Bulletin 24 (1), 43–69. Beven, K.J., Kirkby, M.J., Schofield, N., Tagg, A.F., 1984. Testing a physically-based flood forecasting model (TOPMODEL) for three UK catchment. Journal of Hydrology 69, 119–143. Brasington, J., Richards, K., 1998. Interactions between model predictions, parameters and DTM scales for TOPMODEL. Computers and Geosciences 24 (4), 299–314. Bruneau, P., Gascuel-Odoux, C., Robin, P., Merot, P., Beven, K.J., 1995. Sensitivity to space and time resolution of a hydrological model using digital elevation data. Hydrological Processes 9, 69–82. Chen, J., Kumar, P., 2001. Topographic influence on the seasonal and interannual variation of water and energy balance of basins in North America. Journal of Climate 14, 1989–2014. Ducharne, A., Koster, R.D., Suarez, M.J., Stieglitz, M., Kumar, P., 2000. A catchment-based approach to modeling land surface processes in a general circulation model: 2. Parameter estimation and model demonstration. Journal of Geophysical Research 105 (D20), 24823–24838. Famiglietti, J.S., Wood, E.F., 1994. Multiscale modeling of spatially variable water and energy balance processes. Water Resource Research 30 (11), 3061–3078. Gallant, J.C., Hutchinson, M.F., 1996. Towards an understanding of landscape scale and structure. In: Proceedings The Third International Conference/Workshop on Integrating GIS and Environmental Modeling, Santa Fe, New Mexico, USA, CD ROM, [online], available from /http://www.ncgia.ucsb.edu/conf/SANTA_FE_CD-ROM/sf_papers/ gallant_john/paper.htmlS (accessed 2 July 2008). Giuseppe, M., Aurelia, S., 1997. Information content theory for the estimation of the topographic index distribution used in TOPMODEL. Hydrological Processes 11, 1099–1114.

Hutchinson, M.F., Dowling, T.I., 1991. A continental hydrological assessment of a new grid-based digital elevation model of Australia. Hydrological Processes 5, 45–58. Iorgulescu, I., Jordan, J.P., 1994. Validation of TOPMODEL on a small Swiss catchment. Journal of Hydrology 159, 255–273. Jenson, S.K., 1991. Applications of hydrological information automatically extracted from digital elevation models. Hydrological Processes 5, 31–44. Jenson, S.K., Domingue, J.O., 1988. Extracting topographic structure from digital elevation data for geographic information system analysis. Photogrammetric Engineering and Remote Sensing 54, 1593–1600. Koster, R.D., Suarez, M.J., Ducharne, A., Stieglitz, M., Kumar, P., 2000. A catchment-based approach to modeling land surface processes in a general circulation model: 1. Model structure. Journal of Geophysical Research 105 (D20), 24809–24822. Kumar, P., Verdin, K.L., Greenlee, S.K., 2000. Basin level statistical properties of topographic index for North America. Advances in Water Resources 23, 571–578. Mendicino, G., Sole, A., 1997. The information content theory for the estimation of the topographic index distribution used in TOPMODEL. Hydrological Processes 11, 1099–1114. Niu, G.-Y., Yang, Z.-L., 2003. The versatile integrator of surface and atmosphere processes. Part 2: Evaluation of three topography-based runoff schemes. Global Planet Change 38, 191–208. Niu, G.-Y., Yang, Z.-L., Dickinson, R.E., Gulden, L.E., 2005. A simple TOPMODEL-based runoff parameterization (SIMTOP) for use in global climate models. Journal of Geophysical Research 110, D21106. Nourani, V., Mano, A., 2007. Semi-distributed flood runoff model at the subcontinental scale for southwestern Iran. Hydrological Processes 21, 3173–3180. Pan, F., Peters-Lindard, C.D., Sale, M.J., King, A.W., 2004. A comparison of geographical information systems-based algorithms for computing the TOPMODEL topographic index. Water Resource Research 40, W06303. Qin, C., Zhu, A.-X., Pei, T., Li, B., Zhou, C., Yang, L., 2007. An adaptive approach to selecting a flow-partition exponent for a multipleflow-direction algorithm. International Journal of Geographical Information Science 21, 443–458. Quinn, P., Beven, K.J., Planchon, O., 1991. The prediction of hillslope flow paths for distributed hydrological modeling using digital terrain models. Hydrological Processes 5, 59–79. Quinn, P., Beven, K.J., Lamb, R., 1995. The ln(a/tan b) index: how to calculate it and how to use it in the TOPMODEL framework. Hydrological Processes 9, 161–185. Saulnier, G.-M., Obled, Ch., Beven, K.J., 1997. Analytical compensation between DTM grid resolution and effective values of saturated hydraulic conductivity within the TOPMODEL framework. Hydrological Processes 11, 1331–1346. Stieglitz, M., Rind, D., Famiglietti, J., Rosenzweig, C., 1997. An efficient approach to modeling the topographic control of surface hydrology for regional and global modeling. Journal of Climate 10, 118–137. Wolock, D.M., 1993. Simulating the variable-source-area concept of stream flow generation with the watershed model TOPMODEL. US Geological Survey Water-Resources Investigation Report 93-4124, Lawrence, KS, USA, 33pp. Wolock, D.M., 1995. Effects of subbasin size on topographic characteristics and simulated flow paths in Sleepers River watershed, Vermont. Water Resource 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 Resource Research 30, 1665–1680. Wolock, D.M., McCabe, G.J., 1995. Comparison of single and multiple flow direction algorithms for computing topographic parameters in TOPMODEL. Water Resource Research 31 (5), 1315–1324. Wolock, D.M., McCabe, G.J., 2000. Differences in topographic characteristics computed from 100- and 1000-m resolution digital elevation model data. Hydrological Processes 14, 987–1002. Wolock, D.M., Hornberger, G.M., Musgrove, T.J., 1990. Topographic effects on flow path and surface water chemistry of the Liyn Brianne catchments in Wales. Journal of Hydrology 115, 243–259. Zhang, W., Montgomery, D.R., 1994. Digital elevation model grid size, landscape representation, and hydrologic simulations. Water Resource Research 30, 1019–1028. Zhou, Q., Liu, X., 2002. Error assessment of grid-based flow routing algorithms used in hydrological models. International Journal of Geographical Information Science 16, 819–842. Zhou, Q., Liu, X., 2004. Analysis of errors of derived slope and aspect related to DEM data properties. Computers and Geosciences 30, 369–378.