Proc. of Monitoring, Prediction and Mitigation of Water-Related Disasters (MPMD2005), pp. 91-96, 2005

Scale Issues in Hydrological Geomorphology and Developing Scale Invariancein Surface Flow Hydrology
Nawa Raj Pradhan1, Yasuto Tachikawa 2 and Kaoru Takara2
1

Graduate School of Civil Engineering, Kyoto University, Japan. e-mail: pradhan@rdp.dpri.kyoto-u.ac.jp
Disaster Prevention Research Institute, Kyoto University, Uji, Kyoto 611-0011, Japan

2

ABSTRACT: It is found that basin hydrological response in relations with the drainage basin dominating
geomorphological parameters is directly influenced by the scale of DEM resolution. A Scale Invariant model forthe topographic index distribution (Pradhan et al., 2004a) has fulfilled a part of this gap. A scale independent
relationship in flood routing models in a distributed hydrological model is yet to be developed. To overcome thisproblem, scale laws that govern the relation in digital elevation data resolution on upslope contributing area has
been analyzed and a mathematical formulation has been derived that successfully downscaled the upslope
contributing area from coarse resolution DEM to target fine resolution DEM. The method to downscale theupslope contributing area is used to obtain the similar distribution of depth, cross-section and kinematic wave
celerity from different DEM resolutions in Kamishiiba catchment (210 km2) and to develop a scale invariant
model in the surface flow hydrology.
Key Words: Scale invariant, hydrological geomorphology, topographic index, contributing area, flow routing
1

INTRODUCTION

Hydrological geomorphology is literally the interface between hydrology, the science of water andgeomorphology, the study of landforms and their causative processes. Despite the enormous capacity of today’s(and tomorrow’s) information technologies, the complexity of the Earth’s surface is such that the mostvoluminous descriptions are still only coarse generalizations of what is actually present. The need for continued
and sustained research on scale issues in hydrological geomorphology is therefore self-evident.
Much of the spatial variability can be ignored at “small” spatial scales on the order of 0.1- 1.0m. Indeed, thescientific understanding of individual hydrologic processes at laboratory scales, such as flow through saturated
and unsaturated columns of porous media, is fairly well advanced. In particular, one wants to know how thelaboratory-scale equations can be spatially integrated so as to describe the hydrologic cycle over a hillside. As
the spatial scale under consideration increases to that of a single hillside, spatial variability becomes important,and new elements begin to influence the hydrologic mass balance, such as the topography of the hillside. With
the development of the scale invariant model for the topographic index distribution, Pradhan et al. (2004a)
showed how the laboratory-scale equations can be spatially integrated to provide a consistent hydrologic massbalance in a topography driven model, TOPMODEL.
Beyond a single hillside, a river basin can be viewed as a channel-network-hills system. The hydrologic
cycle for larger sub-basins involves the spatially integrated behaviors of several hills along a channel network.
An understanding of the spatial variability among hillsides and their interactions through a channel network isnecessary for this integration. Thus, at this point scale invariance in surface wave models finds an importantcomponent of the hydrologic cycle in river-basin hydrology.
The basic guide line set by this research to obtain physically based hydrological relationships independent ofregions and scales is to develop an effective translation method of the scale dependence relations of thedominating hydrological and geomorphologic processes linked to typical properties of the catchment into
effective hydrological model. Thus, this research is focused on the development of the scale invariance incatchment hydrology to develop a model consistent with observations. The model can be a potential tool topredict ungauged basins in an effective way.
2

SCALE INVARIANT MODEL FOR TOPOGRAPHIC INDEX (Pradhan et al., 2004a)

It is particularly surface water hydrology that interacts with geomorphology although recently there has been an
increasing convergence between research in geomorphology and in groundwater hydrology. In addition tocontents page

relations between drainage basin characteristics and basin hydrological response, geomorphologists have madeparticular contributions in the investigation of runoff producing areas and the dynamic ways in which such areas
contribute to the generation of stream hydrographs, including headwater drainage systems and the modeling of
their role in runoff production, TOPMODEL (Beven and Kirkby, 1979).
TOPMODEL allows for spatial heterogeneity by making calculations on the basis of the distributionfunction of an index of hydrological similarity, soil topographic index, given by Equation (1):

=

SI

⎧
⎨
⎩ To

ln

⎫
⎬
β ⎭

a
tan

(1)

where SI is the soil topographic index, a is the area draining through a point, tanβ is the local slope angle atthat point and To is the local down-slope transmissivity at soil saturation. Points with the same value of the index
will be predicted as having the same hydrological responses. Specifying a spatial distribution for To being much
more problematic in most applications, it has been assumed to be spatially homogeneous. In the case thesimilarity index defined by Equation (1) reduces to Equation (2).
⎧
⎫
a
(2)
TI
= ln ⎨⎬
tan
β⎩⎭where TI is the topographic index. Topographic index is scale dependent (Zhang and Montgomery, 1994) whichleads identified parameter values to be dependent on DEM resolution. When the resolution of DEM change for acatchment, the spatial distribution of topographic index also change. In this scenario keeping the same effective
parameter value of To for changed topography index cannot fulfill the hydrological similarity concept of
TOPMODEL across different DEM resolutions in the catchment.
Figure 1(a) shows the blunder in simulated result (with Nash efficiency of -45%) of 1000m DEM resolution
TOPMODEL, with effective parameter values identified at 50m-grid resolution DEM (refer Table 1). For thispurpose notions of scale transformation and scale invariance are needed. We developed a method to downscale
topographic index which is given by equation (3):⎛Ci
TI scaled = ln ⎜⎜W R θ
⎝ i f scaled

⎞
⎟
⎟
⎠

(3)

67

73

79

85

91

97

115

61

109

55

103

97

49

91

43

85

37

79

31

73

25

67

19

61

13

55

49

7

43

1

37

400

31

/s )
3

600

(m

(m

800

D i s c h a rg e

3

1000

D i s c h a rg e

/s )

1200

25

0.04

1400

19

0.03

13

0.02

7

0
0.01

1

115

109

103

97

91

85

79

73

67

61

55

49

43

37

31

25

19

13

7

P re c ip i ta ti o n ( m /h r)

P re c ip i ta ti o n ( m /h r)

1where TIscaled is the scaled topographic index, Ci is the upslope contributing area of the coarse resolution DEM,
Wi is the unit contour length of coarse resolution DEM, i is a location in catchment, Rf is a resolution factor andθscaled is the downscaled steepest slope of the target resolution DEM by fractal method. Details of the derivationof Equation (3) is given in Pradhan et al. (2004a). Equation (3) is coupled with TOPMODEL to develop the
Scale Invariant TOPMODEL (Pradhan et al., 2004b). The Scale Invariant TOPMODEL is applied to theKamishiiba catchment (210 km2).
0
0.01
0.02
0.03
0.04
800

700
600
500
400
300
200
100
0

200
0
1

7

13

19 25 31 37 43 49 55

61 67 73 79 85

91 97 103 109 115

103 109 115

Ti me (hour )Obs er ved dis char ge

Time ( hour )Calculated by us ing TOP MODEL with 50m gr id r es olution DEM and ef f ective par ameter scalibr ated at this r es olution DEM

Ob se r v e d d isc h a r g e
Ca lc u la t e d b y u sin g TOP MODEL wit h 5 0 m g r id r e so lu t io n DEM a n d e f f e c t iv e
p a r a m e t e r s c a lib r a t e d a t t h is r e so lu t io n DEM
Ca lc u la t e d b y u sin g TOP MODEL wit h 10 0 0 m g r id r e so lu t io n DEM a n dp a r a m e t e r s d e r iv e d f r o m 5 0 m g r id r e so lu t io n DEM

Calculated by downs caling topogr aphic index f r om 1000m gr id r es olution DEM to 50 m gr id
r es olution DEM and par ameter s der ived f r om 50m gr id r es olution DEM(a)(b)
Figure 1. Comparison of simulation results at Kamishiiba catchment (210 km2) for different resolution DEM byusing the same effective parameter value identified at 50m DEM resolution; (a) without scale invariant modelfor topographic index distribution (b) with Scale Invariant model for topographic index distribution.

92

Table 1. Effective parameter values identified at 50m DEM resolution in Kamishiiba catchment (210 km2).
Lateral transmissivity of soil atsaturation
condition,To[m2/hr]9.8

decay factor of lateraltransmissivity with respect tosaturation deficit, m [m]0.07

Maximumroot
storage, Rzmax [m]

zone

0.001

Manning’s roughness
coefficient n
0.037

Figure 1(b) shows that the simulated runoff from Scale Invariant TOPMODEL applied at 1000m-grid
resolution DEM, with the same effective parameter values derived from 50m-grid resolution DEM in Table 1,
has matched with the simulated runoff of 50m DEM resolution TOPMODEL and also with the observed runoff
with Nash efficiency of 90%. The method to downscale topographic index distribution has solved two problems.
Firstly, it has given consistency to the physically reliable effective parameter value of To with observations at thescale of interest; it can’t be denied the fact that physically realistic value of To also finds equally important ininfiltration excess dominating hydrological processes. Secondly, it has successfully reduced the uncertainty
associated with a scale in prediction.
Although the relative importance of the components in the rainfall-runoff transformation process depends
both on its working scale and on the geographical, climatic and environmental conditions of the site underconsideration, the relative importance of routing phenomenon in surface flow hydrology cannot be ignored for acomplete process model that offers a detailed description of the rainfall-runoff transformation. In Figure 1 (b),the same effective parametric value of Manning’s roughness coefficient n identified at 50m DEM resolution (seeTable 1) is used to 1000m DEM resolution without analyzing the scale effects in routing. Thus, the hydrologicalresponse of scale invariant TOPMODEL alone in Figure 1(b) still suffers discrepancy from that of 50m DEM
resolution hydrological response and the observed discharge by overestimating the peak flows.
3

SCALE EFFECTS IN THE INTERFACE BETWEEN SURFACE FLOW HYDROLOGY AND
GEOMORPHOLOGY: ANALYSIS, SOLUTION METHOD, RESULTS AND DISCUSSIONS

Flow routing in channels has been a subject of much discussion for over half a century and more especially since
the advent of digital computers. Flow routing is a technique for determining the propagation of flow from one
point in the channel to another. Flow routing in open channels entails wave dispersion, wave attenuation or
amplification and wave retardation or acceleration. These wave characteristics constitute the hydraulics of flowrouting or propagation and are greatly affected by the geometric characteristics of channels. The flow variables
whose propagation characteristics are of interest are discharge, velocity, depth, cross-section, volume andduration. In catchment hill slope channel routing these flow variables is a function of upslope contributing area.
In the DEM based distributed hydrological model, the smaller contributing area less than a grid area of a coarse
DEM resolution used is completely lost. In Figure 2(a), the smaller contributing area less than a km2 that appears
over 97% in 50m DEM resolution is seen completely lost when 1000m DEM resolution is used.
3.1 Development of Scale Invariant model for the upslope contributing area
The density of the small contributing area is higher in a catchment. In Figure 3(a), it is observed that this small
1.2

1.05
1

0.8

Fraction of catchment area

Frac tion of c atc hment area1
At 1000m DEM r esolut ion

0.6
0.4

At 50m DEM r esolut ion

0.2

0.95

At 50m DEM r esolut ion

0.9
Downsc aled f r om 1000m DEM
r esolut ion t o 50m DEM r esolut ion

0.85
0.8

0
1

17

33

49

65

81

97

113

129

145

161

177

193

1

209

17

33

49

65

81

97

113

129

145 161

177 193 209

Upslope c ont r ibut ing ar ea ( k m 2 )Upslope c ont r ibut ing ar ea ( k m 2 )(b)(a)
Figure 2. Comparison of upslope contributing area distribution function from different DEM resolutions in
Kamishiiba catchment (210 km2) (a) without downscaling method for upslope contributing area (b) withdownscaling method for upslope contributing area.

93

contributing area less than a grid area of the coarse resolution DEM used is completely lost. In fact the smallest
contributing area derived from a DEM resolution is a single grid of the DEM at that resolution. Thus area
smaller than this grid resolution is completely lost as the larger sampling dimensions of the grids act as filter. But
as we use finer resolution DEM, the smaller contributing area - that is the area of finer grid resolution isachieved. From this point of view, we introduced number of sub grids Ns to derive scaled upslope contributing
area as shown by Equation 4.

⎛ Ci ⎞⎟
C i scaled = ⎜⎜ N I ⎟
⎝ s f ⎠

(4)

where Ciscaled is the scaled upslope contributing area at a point i, If is a influence factor. Ns is the total number ofsubgrids within a coarse resolution grid. i is a location in a catchment. The area of a coarse resolution grid shown
itself is the smallest contributing area for that coarse DEM resolution. When this area of coarse resolution DEM
is divided by the number of sub grids (sub grid as target resolution DEM) that together adds up to make thecoarse resolution grid, area of a sub grid as smallest contributing area for the target DEM resolution is obtained.
In a catchment as the upslope contributing area gets bigger and bigger, the contributing area values given bycoarse and fine resolution DEM at the points downstream becomes closer and closer; thus the influence of Ns on
Ci must gradually decrease in Equation (4). For this reason we introduced influence factor If in Equation (4) and
If is described as;⎧ (1 − N i ) H ⎫⎨I f = e⎩

⎬
⎭

No

(5)

where, Ni is the number of the coarse resolution grids contained in the contributing area at a location i in thecatchment, N0 is the number of the coarse resolution grids contained in the contributing area at the outlet of thecatchment. H in Equation (5) is introduced as harmony factor. Considering the influence of Ns on Ci in Equation(4) is almost negligible at the outlet of the catchment, the value of H can be obtained from Equation (6) as;

N se −H = 1(6)Finally, we developed a Scale Invariant model for the upslope contributing area from Equations (4) and (5) as;

C i scaled

⎞
⎛
⎟
⎜
Ci
⎟
= ⎜
⎧ (1 − N i ) H ⎫⎜⎨⎬ ⎟⎜ N e ⎩ No ⎭ ⎟
⎠
⎝ s

(7)

Using Equation (7), the upslope contributing area is downscaled from 1000m DEM resolution to 50m DEM
resolution. In contrast to Figure 2(a), Figure 2(b) shows the similar distribution of upslope contributing area from
50m DEM resolution and downscaled from 1000m DEM resolution to 50m DEM resolution.
3.2 Development of Scale Invariance in surface flow hydrology
A wave is a variation in flow, such as a change in flow rate or water surface elevation, and the wave celerity isthe velocity with which this variation travels along the channel. The kinematic wave celerity, ck, can be definedin terms of flow depth by Equation (8).
1
5 ⎛⎜ S i 2 ⎞⎟ 2 3
(8)
ck =
yi
3 ⎜⎝ n ⎟⎠where Si is the slope and n is the Manning’s roughness coefficient. yi is the depth of flow and is expressed as;
3

⎞ 5
⎛ nQ
i
⎟
⎜
yi =
⎜ S 12 B ⎟
i ⎠
⎝ i

(9)

where Qi and Bi are the flow rate and channel width respectively at a point i.
In a distributed system routing the flow is calculated as a function of space and time through the system.
Like lateral transmissivity of soil To in Equation (1), the Manning’s roughness coefficient n in Equation (8) and(9) is also an effective parameter. Figure 3(a) shows much difference in distribution of depth yi from 50m DEM
resolution and from 1000m DEM resolution when keeping the same value of Manning’s roughness coefficient

94

1.2

1
0.8

At 1000m DEM resolution

0.6

Downscaled from 1000m DEM to
50m resolution
At 50m DEM resolution

0.4
0.2
0

Fraction of catchment area

Fraction of catchment area

1.2

1

From 1000m DEM resolution

0.8
0.6

From 50m DEM resolution

0.4

Downscaled from 1000m DEM
resolution to 50m DEM resolution

0.2
0

0.1 0.4 0.7

1

1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4 3.7

4

4.3 4.6 4.9 5.2 5.5 5.8 6.1 6.4 6.7

7

0.3

0.6

0.9

Depth (m)

1.2

1.5

1.8

2.1

2.4

2.7

3

3.3

3 .6

3.9

4 .2

4.5

4.8

5.1

5.4

5.7

6

6.3

6.6

6.9

Kinematic wave celerity (m/s)(a)(b)
Figure 3. Comparison of (a) flow depth distribution function from different DEM resolutions in Kamishiiba
catchment (210 km2) without and with it’s downscaling method (b) Kinematic wave celerity distribution fromdifferent DEM resolutions in Kamishiiba catchment (210 km2) without and with it’s downscaling method.

for the 1000m DEM resolution, identified at 50m DEM resolution (0.037 in Table 1). This scale problem in
depth has serious impact on the distribution of kinematic wave celerity as shown in Figure 3(b).
The root of this problem originates from the scale problem on upslope contributing area as discussed earlier
in section 3.1 and in figure 2(a), and also from the underestimation of slope in coarse resolution DEM. Upslope
contributing area is a key variable because of its intrinsic capability to describe the nested aggregation structure
embedded in the fluvial landforms and its important physical implications (e.g., Rodriguez-Iturbe and Rinaldo,1997). Flow rate Qi being a function of upslope contributing area (Strahler, 1964), from Equation (7) weformulate to obtain the flow rate given by target resolution DEM at point i as;
Q i t arg

et

⎛
⎜
= ⎜
⎜
⎜ N
⎝

Q i

s

e

⎧ (1 − N i ) H
⎨
N o
⎩

⎫
⎬
⎭

⎞
⎟
⎟
⎟
⎟
⎠

(10)

The downscaled channel width at a location i, Biscaled is given as;b(11)
B i scaled= aC iscaled
where a and b are the coefficients whose values are independent of the scale of DEM.
Substituting Qi ,Bi and Si in Equation (9) by Qitarget and Biscaled from Equations (10) and (11) and by θscaled (refer
Pradhan et al., 2004a) respectively we develop the method to downscale the flow depth, y iscaled as;y iscaled⎛nQ it arg et= ⎜1⎜θ⎝ scaled 2 B iscaled

3

⎞ 5
⎟
⎟
⎠

(12)

Finally, substituting yi and Si in Equation (8) by y iscaled and θscaled we developed the method to downscale thewave celerity distribution. In Figure 3 (a) (all the simulation results in Figures 3(a) and 3(b) are made at time
step 43 hours of the rainfall event), it is shown that the distribution of downscaled flow depth from 1000 m DEM
resolution to 50m DEM resolution and that from 50m DEM resolution (Manning’s roughness coefficient n used
is identified at 50m DEM resolution) has matched. Thus, by using Equation (12) we successfully reduced theover estimation of depth given by 1000m DEM resolution (also shown in Figure 3(a)). Figure 3(b) shows theperfect match of the distribution of downscaled wave celerity ck from 1000 m DEM resolution to 50m DEM
resolution and that from 50m DEM resolution (n value used is identified at 50m DEM resolution; refer Table 1).
Several variations of the kinematic wave routing method have been proposed. These routing methods can beeasily coupled with Equation (11) and Equation (12) to develop a scale invariant routing method. Cunge (1969)
proposed that Muskingum method can be considered an approximate solution of a modified diffusion equation.
Coupling this routing method with the method to downscale depth and kinematic wave celerity in Equation (11)and (12) we further analyzed the simulation results with the Scale Invariant TOPMODEL as shown in Figure 4.
In Figure 4, the simulated result downscaled from 1000m DEM resolution (with downscaling routing method
and Scale Invariant TOPMODEL, using the same parameters in Table 1 identified at 50m DEM resolution) has
perfectly matched with the simulation result of 50m DEM resolution with further increment in Nash efficiency
from 90% (in Figure 1(b)) to 93%.
4

CONCLUSION

There is a long tradition in geomorphology of seeking generalizable rules for landscape evolution such that reallandscapes, and particularly their scale-dependent attributes, can be modeled. However, basin hydrologicalresponse in relations with the geomorphological parameters are influenced by DEM resolution. We developed

95

1 16

1 11

1 06

1 01

96

91

86

81

76

71

66

61

56

51

46

41

36

31

26

21

16

11

6

1

Precipitation (m/hr)

0
0.01
0.02
0.03
0.04

3

Discharge (m /s)

800
700
600
500
400
300
200
100
0
1

6
11 16 21 26 31 36 41 46 51 56 61 66 71 76 81 86 91 96 101 106 111 116
Time (hour)
Calculated by using T OPMODEL with 50m grid resolution DEM and effective parameters calibrated at this resolution DEM
Calculated by downscaling topographic index from 1000m grid resolution DEM to 50 m grid resolution DEM and parameters
derived from 50m grid resolution DEM
Calculated by downscaling kinematic wave celerity distribution and topographic index distribution from 1000m grid resolution
DEM to 50 m grid resolution DEM and parameters derived from 50m grid resolution DEM

Figure 4. Comparison of simulation results at Kamishiiba catchment (210 km2) at 50m DEM resolution, at 1000m DEM
resolution with downscaled topographic index distribution and at 1000m DEM resolution with downscaled topographic index
distribution and kinematic wave celerity. For all the simulations, same parameters identified at 50m DEM resolution is used.

the concept of resolution factor and a fractal method for scaled steepest slope to account for the effect of scale on
topographic index distribution. The method to downscale the topographic index is then coupled with
TOPMODEL to develop the Scale Invariant TOPMODEL. In this research we analyzed the scale laws thatgovern the relation in digital elevation data resolution on upslope contributing area and developed amathematical formulation to downscale the upslope contributing area. The method to downscale the contributing
area is successfully applied to downscale the flow variables and developed a scale invariant model in surface
flow hydrology. We propose the scale invariant models in basin hydrology as a method to enhance theconsistency of a complete process model in rainfall-runoff transformation. This research aims to offer applied
opportunities to identify physically based hydrological relationships independent of scale and region (Pradhan etal., 2005).
5 ACKNOWLEDGEMENTS
A part of this study is supported by “System Modeling Approaches for Assessment of Interaction Between Social
Changes and Water Cycle (PI: Prof. Kaoru Takara, Disaster Prevention Research Institute, Kyoto University)”,which is conducted under the CREST (Core Research for Evolutional Science and Technology) Program by theJapan Science and Technology Agency.
REFERENCES:
Beven, K. J. and Kirkby, M. J. (1979). A physically based, variable contributing area model of basin hydrology,
Hydrol. Sci Bull., 24, pp. 43-69.
Cunge, J.A. (1969). On the subject of a flood propagation computation method (Muskingum Method). Journel of
Hydraulic Research, 7, pp. 205-230.
Pradhan, N. R., Tachikawa, Y., and Takara, K. (2004a). A scale invariance model for spatial downscaling of
topographic index in TOPMODEL, Annual Journal of Hydraulic Engineering, JSCE, vol. 48, pp. 109-114.
Pradhan, N. R., Tachikawa, Y., and Takara, K. (2004b). Scale Invariant TOPMODEL for the solution of DEM
resolution effects and parameter inconsistency in hydrological modeling, Proc. of the 2nd Asia Pacific
Association of Hydrology and Water Resources, Singapore, Vol. 2, pp. 641-648.
Pradhan, N. R., Tachikawa, Y., and Takara, K. (2005). Development of a transferable hydrologic modeling in
TOPMODEL framework across scale and region, Annual Journal of Hydraulic Engineering, JSCE,(submitted).
Rodriguez-Iturbe, I., and Rinaldo, A. (1997). Fractal River Basins: Chance and Self-organization. Cambridge
University Press.
Strahler, A.N. (1964). Geology-Part II, In Handbook of Applied Hydrology. Chow VT (ed.), McGraw-Hill, New
York, pp. 39-76.
Zhang, W. and Montgomery, D. R. (1994). Digital elevation model grid size, landscape representation, andhydrologic simulations, Wat. Resour. Res., 30(4), pp. 1019-1028.

96

contents page