CE74.55 Modeling of Water Resources Systems
February 24th and 25th, 2011

DISTRIBUTED RAINFALL-RUNOFF MODELING
Yasuto TACHIKAWA
Dept. of Civil and Earth Resources Eng., Kyoto University, Kyoto 615-8540, Japantachikawa@hywr.kuciv.kyoto-u.ac.jp

Modeling the actual hydrological response of a basin is a difficult task as the actual basin processes
are complicated and basin properties change temporally as well as spatially. However, this difficulty can beaddressed to some extent by using a distributed hydrological modeling approach, which is capable ofincorporating the basin heterogeneity to model the basin hydrological responses. The paper introduces abasic idea of distributed hydrological modeling, which is organized as follows:
z
z
z
z
z
z

1

Section 1 explains a basic concept of distributed hydrological modeling.
Section 2 presents a typical water flow model embedded in a distributed hydrological model. As aflow model, the kinematic wave flow on a rectangular plane is explained.
Section 3 explains numerical solution methods of the kinematic wave flow. A typical explicit method
and an implicit method are introduced.
Section 4 represents extensions of the kinematic wave model for surface and subsurface flow
modeling.
Section 5 shows applications of a distributed hydrologic model to basins with different catchment andclimatic characteristics.
Section 6 shows applications of a distributed hydrologic model for future discharge projection underclimate change.

BASIC CONCEPT OF DISTRIBUTED HYDROLOGIC MODELING

1.1 Open-book type catchment modeling

A simple distributed hydrologic representation is a rainfall-runoff model based on an open-book
catchment model. Figure 1 illustrates a schematic representation of an open-book catchment model. As its
name implies, an open-book catchment model consists of two rectangular planes and a stream. Rain waterflows on the planes and drains into the stream. The stream flow drains out of the catchment outlet. The
flows on the rectangular planes and the stream are routed using a flow model such as the kinematic wave
model. This is one of the simplest approximations of catchment hydrology using a distributed hydrologic
model. A more complex structure of a catchment can be modeled as a cascade of the open-book catchment
model as shown in Figure 2. Spatial distribution information of topography, soil characteristics, land coverand rainfall intensity are used for each sub-catchment.

Figure 1: Schematic drawing of an open-book catchment
model.
1

Figure 2: Cascade of the open-book catchment models.

1.2 Catchment modeling using Digital Elevation Models

Digital elevation models are now available at any catchments with a high spatial resolution enough todescribe local catchment topography. For example, HydroSHED (Hydrological data and maps based on
SHuttle Elevation Derivatives at multiple Scales, http://hydrosheds.cr.usgs.gov/index.php) provides
hydrographic information for regional and global-scale, which includes digital elevations, drainage
directions, and flow accumulation with the spatial resolutions of 3 arc-second (about 100m), 15 arc-second(about 500m) and 30 arc-second (about 1km). Figure 3 is an example of flow direction information derived
from a grid-based DEM using the eight-direction method. After deriving the flow direction information, it
is easy to define a rectangular plane formed by two adjacent grid points as shown in Figure 4. Catchment
topography is represented by a set of slope units. For each slope unit, its area, length and gradient used for aflow model are easily calculated. A flow model such as the kinematic wave model is applied to all slope
units and runoff is routed according to the flow direction information.

Figure 3: Flow direction information using a DEM at the Maruyama Basin in Japan.

2

Figure 4: Schematic drawing of a catchment modelling using a DEM, which is a set of slope units
connected according to the flow direction.

2

KINEMATIC WAVE MODEL

The kinematic wave model was originally developed for river routing. Then, it was applied tocatchment rainfall-runoff modeling.
2.1 Basic theory

The basic form of the kinematic wave flow equation is:∂A ∂Q+= q L ( x, t )
∂t ∂x

(1)

Q = αAm(2)

where A(x, t) is the flow cross-sectional area, Q(x, t) is the flow discharge, qL(x, t) is the lateral inflow perunit length and α and m are constants. Eq. (1) is the continuity equation and is derived from the principleof mass conservation within a control volume. Eq. (2) is derived from Manning's or Chezy's laws which areflow resistance laws of open channel uniform flow. Manning's law is described as:v=

1 1/ 2 2 / 3
I R
n

(3)

in which v is the average flow velocity, n is the Manning's roughness coefficient, R is the hydraulic radius,
and I is the friction slope. By assuming that R is expressed as z1Az2 (z1 and z2 are constants), we can easily
obtain the expression of Eq. (2).

3

Figure 5: Shallow water flow on a rectangular plane.

2.2 Kinematic wave model for shallow surface flow on a rectangular plane

A catchment is often modeled as a rectangular plane or a set of rectangular planes as shown in Figure
2 and 4. The rainfall-runoff on a rectangle slope is modeled as a shallow water flow as shown in Figure 5.
The hydraulic radius of shallow water flow can be approximated by the water depth h. Then, Eq. (3) iswritten as:v=

1 1/ 2 2 / 3
I h
n

(4)

Then, discharge per unit width q is expressed as:

q = vh =1 1/ 2 5 / 3
I h = αh mn(5)in which α = I / n , m = 5/3. The value of I is usually approximated as the slope gradient of the rectangularplane, and the value of n is calibrated from observation or determined empirically. Eq. (5) is the momentum
equation of the kinematic wave model on shallow water flow on a rectangular plane. The continuity
equation becomes:∂h ∂q+= r = qL / B
∂t ∂x

(6)

by dividing both sides of Eq. (1) with the width of a rectangular plane B, where r is rainfall intensity. Eqs.
(5) and (6) are the basic equations of the kinematic wave model for shallow surface flow on a rectangularplane.

4

3

NUMERICAL SOLUTIONS OF KINEMATIC WAVE MODEL

The combination of Eqs. (5) and (6) derives a nonlinear first-order partial differential equation which
describes the change of q and h in time and space. This equation can be solved by using numerical methods.
Several methods are available, for example, the finite difference method, the finite element method, and thefinite volume method. Here, the finite difference methods are introduced to solve the equation.
Finite difference methods are classified into two types, explicit and implicit methods. Explicit methods
evaluate all spatial derivatives in a partial differential equation at the time level where the solution is known.
Implicit methods evaluate all spatial derivatives at the time level where the solution is unknown. Generally,
implicit methods are more complex than explicit methods; however the former is better in numericalstability to the latter.

3.1 Explicit method (The Lax-Wendroff scheme)

The Lax-Wendroff scheme is a typical explicit method to solve the kinematic wave flow numerically.
Figure 6 shows the process to obtain the unknown variables defined at the grid points on discretized x-t
plane. The subscripts i and j denote incremented time and space levels, respectively. This scheme calculateshji+1 from hj-1i, hji , and hj+1i as stated below.
Expanding h( x, t + Δt ) using the Taylor's series up to the second order, we can obtain:(7)
where Δt is the time step size. Eq. (6) can be rewritten as:(8)Using Eq. (8), ∂ h / ∂t can also be rewritten as:
2

2

(9)
The right hand side of Eqs. (8) and (9) do not have derivatives of q and h with respect to time. We
replace all derivatives in these expressions by difference operators as:

(10)

5

Unknown value
Known value

Figure 6: Schematic drawing of the Lax-Wendroff scheme on discretization of x-t plane.

(11)
where Δx is the distance between two adjacent grid points, Δt is the time step size and f '= dq / dh . Let
the right hand side of Eq. (10) be H1 and that of Eq. (11) be H2, then we obtain:(12)This equation means hji+1 (and also qji+1) can be calculated from hj-1i, hji, and hj+1i as shown in Figure 6.
For j=n, the following backward difference approximation:(13)
is used. When the upper boundary condition h0i, or q0i (i = 0, 1, …), and the initial conditions hji or qji, (j =0, 1, …, n) are given, hji+1 can be calculated sequentially.
3.2 Implicit method (The four point implicit scheme)

The celerity of the kinematic wave is determined as:(14)
By using Eq. (14), the continuity equation (6) can be rewritten as:(15)
6

Unknown value
Known value

Figure 7: Schematic drawing of the four point implicit scheme on discretization of x-t plane.

Using the discretized x-t plane shown in Figure 7, each term of Eq. (15) can be written as:

(16)

(17)
in which Δx is the distance between two adjacent grid points, Δt is the time step size, θ is the timeweighting parameter, and λ is the space weighting parameter. The default value of θ is 0.6 and λ is
0.5.
i +1
In Eq. (16), the unknown variable is q j +1 . Collecting unknown terms to the left hand side of theequation, we obtain:

(18)
i +1

i +1

Eq. (18) is a nonlinear equation of q j +1 because c* is dependent on the value of q j +1 . To solve theequation, an iterative method is used as follows:
i +1
i
i +1
step 1: Let the initial estimate of qˆ j +1 be 0.5( q j +1 + q j ).
i +1

i +1

step 2: Obtain hˆ j +1 corresponding to qˆ j +1 .
step 3: Calculate c* using Eq. (17).
i +1step 4: Calculate an estimate of q j +1 using Eq. (18).
step 5: If q j +1 − qˆ j +1 < ε ( ε is a tolerance level), q j +1 is obtained. Otherwise, set q j +1 as a new
i +1

i +1

i +1

i +1

estimated value qˆ j +1 and go back to the step 2.
7i +1

4 KINEMATIC WAVE MODEL FOR SUBSURFACE FLOW
Several variations of the kinematic wave model were developed to improve the prediction accuracy ofrainfall-runoff from catchments. These models realize various flow mechanisms, characteristics of soil andland use and these spatial distributions within a catchment.
4.1 Kinematic wave model for surface-subsurface saturated flow (Takasao & Shiiba, 1988)

The kinematic wave model described in 2.2 is applied to infiltration excess overland flow (Hortonian
overland flow). To extend the discharge-depth relationship to represent saturation excess overland flow,
Takasao and Shiiba (1988) developed a new discharge-depth relationship which includes both surface andsubsurface flow and its interaction. They assumed that hillslope surface was covered by a soil layer withhigh surface infiltration capacity and permeability illustrated in Figure 8. In the figure, D is the depth of thesoil layer, L is the length of the hillslope, H is the depth of water flow and r is the rainfall intensity. Let γ
be the porosity of the soil layer and be spatially constant. Then the depth of the effective pore of the soil
layer, d, becomes γD . Also, the actual water depth (water volume in height), h, becomes γH (0 < H < D)or H-D+ γD (H > D). Takasao and Shiiba (1988) derived the following relation between h and q from
Darcy's law and Manning's law:

⎧ah,
q=⎨
m
⎩α (h − d ) + ah,

0≤h≤d

(19)

h≥d

where a = k sinθ / γ is the velocity of water which flows in the soil layer, k is the hydraulic conductivity,and α = sinθ / n and m=5/3, respectively. The upper relation given in Eq. (19) is used for the case where
only subsurface flow occurs and the lower relation in Eq. (19) is used for the case where saturation excess
overland flow occurs.

Figure 8: Schematic drawing of surface-subsurface flow on a hillslope.

4.2 Kinematic wave model for surface-subsurface saturated-unsaturated flow (Tachikawa et al., 2005)

The kinematic wave models described in 2.2 and 4.1 are basically used to simulate flood runoff which
occurs in heavy rainfall events where the gravity force is dominant as a driving force of water flow.
However these discharge-stage relations cannot be used in dry conditions because the capillary force of soil
particles cannot be negligible in comparison with the gravity force in such conditions. To realize long term
runoff simulations from dry to humid seasons continuously, Tachikawa et al. (2005) further extended thedischarge-stage relation to realize both high flow and low flow simulation introducing an unsaturated flowmechanism in capillary pore of hillslope soil layer.
8

q(x)

vm d m + va ( h − d m ) + α ( h − d a ) m

q = g(h)vm d m + va ( h − d m )daSurface flowdm

Non-capillary pore
Saturated flow
Unsaturated flowDqavm d m ( h / d m ) β

F(q(x))

qm

dm

Soil

da

h(x)

Figure 9: Schematic drawing of the discharge-depth relationship for saturated-unsaturated flow.

The discharge depth relation realizes three lateral flow mechanisms including 1) subsurface flow
through capillary pore; 2) subsurface flow through non-capillary pore; and 3) surface flow on the soil layer(see Figure 9). At a slope segment, when the water depth is lower than the equivalent water depth forunsaturated flow (0 < h < dm), flow is simulated by Darcy law with an unsaturated hydraulic conductivity
km. If the water depth exceeds the equivalent depth for unsaturated flow, the exceeded water flows as
saturated subsurface flow that is simulated by Darcy law with saturated hydraulic conductivity ka. Once thewater depth is greater than the effective soil layer d a = γD , the water flows as surface flow, which issimulated by the Manning’s equation.
These three flow processes are represented by the following single set of stage-discharge relationship
Eq. (20) and the continuity equation (6) (Tachikawa et al., 2005; Sayama and McDonnel, 2009):⎧vm d m (h / d m ) β ,0 ≤ h ≤ dm⎪
q = ⎨vm d m + va (h − d m ),dm ≤ h ≤ da⎪m⎩vm d m + va (h − d m ) + α (h − d a ) , d a≤ h(20)

ν m = k m i , ν a = k a i , km = ka / β , α = i / n⎧βvm (h / d m ) β −1 ,∂q ⎪= ⎨va ,
c* =
∂h ⎪
m −1
⎩va + mα (h − d a ) ,0 ≤ h ≤ dm
dm ≤ h ≤ da(21)

d a≤ hwhere q is discharge per unit width, h is water depth, i is the slope gradient, km is the saturated hydraulicconductivity of the capillary soil layer, ka is the hydraulic conductivity of the non-capillary soil layer, dm isthe depth of the capillary soil layer, da represents the depths of the capillary and non-capillary soil layers,
β is the exponent constant of unsaturated flow, vm and va are the flow velocities of unsaturated andsaturated subsurface flows, respectively, and n is the Manning’s roughness coefficient which varies as afunction of land use type. Model parameters in the stage-discharge relationship are dm, km, da, ka, n, andβ (= ka/ km). The kinematic wave celerity for each layer is formulated in Eq. (21).

95 APPLICATION OF DISTRIBUTED MODEL
A distributed hydrologic model using the kinematic wave flow with the discharge-depth relationship(20) is applied to basins with different catchment and climatic characteristics.

Illinois Basin, US

Mae Chaem Basin,
ThailandUpper Kotmale Basin, Sri Lanka

Maruyama Basin, Japan
Figure 10: Study basins and slope distributions.

5.1 Simulation results at Maruyama basin in JapanThe Maruyama River basin is situated in the central part of Japan and the area of the basin is about
909 km2. The catchment land use consist of 85% forest, 7 % rice, 2% buildings, 1% inland water and 5%
others. The average annual rainfall of the basin is about 2,000 mm. Hourly flow data at the Fuichiba station
for three flood events in 2004 and one event in 2003 are used in this study. Hourly rainfall data at ten
stations within the Maruyama basin is also used. To represent the basin topography, the SRTM 90m digital
10

elevation model (USGS 2008, http://seamless.usgs.gov/) was used. Considering the computational time
taken for automatic calibration of the distributed hydrological model, the original digital elevation model(DEM) was resampled to 250m grid size. The derived slope variation of the basin is shown in Figure 10.

Calibration - Event 1

Validation - Event 2

Validation - Event 3

Validation - Event 4

Figure 11: Observed and model predicted hydrographs for the Maruyama basin.

Figure 11 shows the observed and predicted flood hydrographs for the Maruyama basin for model
calibration and validation. The distributed model is calibrated at the basin outlet for flood Event 1 using the
SCE algorithm; then the calibrated model parameter set is used to predict the other flood events forvalidation. The hydrological model with the optimized parameter set for Event 1 predicted the peak flow
and the time to peak of all flood events reasonably. To investigate further on these findings, the model wascalibrated and validated for other three events independently. The model structure captures the governinghydrological processes in the steep and wet Maruyama basin.

5.2 Simulation results at the Upper Kotmale basin in Sri Lanka

The Upper Kotmale basin of 304 km2 above the Thalawakelle station is the upper most basin of the
Mahaweli River, one of the most important water resources in Sri Lanka. The terrain of the area isgenerally mountainous with moderate slope (Figure 10). The basin is situated in the wet zone of thecountry and the average annual rainfall varies from 2,200 mm to 2,600 mm. The catchment land use
consists of tea 44%, forest 36%, home garden 7%, grass 5%, cash crops 5% and scrub. The main soil type11

of the basin is silt loam. The Upper Kotmale basin was set up as an experimental basin in Sri Lanka for
PUB studies and six tipping bucket rain gauges were installed in the basin. The data collection has been
carried out since 2003 and we used hourly rainfall data for the study. Hourly flow data for two flood eventsfrom April to September, 2003 at the basin outlet Talawakelle is also used. The SRTM DEM of 90m wasresampled to 250m grid size to represent the basin topography in the hydrological model.

Calibration-Event 1

Validation-Event 2

Calibration-Event 2

Validation-Event 1

Figure 12: Observed and model predicted hydrographs of the Upper Kotmale basin.

Figure 12 shows the observed and predicted hydrographs for model calibration and validation. Twoparameter sets are obtained for the basin by calibrating the model for each flood event independently using
the SCE algorithm. Then the parameter set obtained for flood Event-1 was used to predict the flood Event-2
and vice-versa. The observed and predicted hydrographs for both events reasonably agreed with each otherin the model calibration stage as well as in the model validation stage. The Nash coefficient value for both
events in the model calibration is over 0.96 and the percentage error in the peak flow prediction is less than
4 %. Moreover, the percentage error in the flow volume is less than 10 % and also showed a high Nashcoefficient value for both events in the model validation stage.
The model for the Upper Kotmale basin showed similar performances that we have experienced inJapanese basins. The identified values of the model parameters are different from Japanese basins, which
would reflect the differences of land use and soil types. Wet and moderate slope basin or wet and steep
slope basin tends to produce saturated subsurface flow and overland flow quickly after rainfall. The applied
model structure fits to reproduce such quick response hydrological behaviors of these catchments. Onlychanging the model parameter values, the hydrological model reproduces hydrographs well in both
calibration and validation. This suggests the model structure is transferable to these catchments.
12

6 APPLICATION OF DISTRIBUTED MODEL FOR RIVER FLOW PROJECTION
UNDER CLIMATE CHANGE
Future river discharge projection was conducted using a distributed hydrologic model with a futureclimate projection data generated by the general circulation model (MRI-AGCM3.1S) developed by theMeteorological Research Institute in Japan. MRI-AGCM3.1S realizes 1920×960 of grid cells of about 20
km spatial resolutions (Kitoh et al., 2009). The products of MRI-AGCM3.1S consists of various
atmospheric and hydrologic variables of the present climate experiment (1979-2003), the near futureclimate experiment (2015-2039), and the future climate experiment (2075-2099), which were simulated
under the SRES A1B scenario.
The river discharge for the Chao Phraya River basins is predicted by feeding the future climate
projection data into the 1km-spatial resolution distributed hydrologic model. The inputted data to the
distributed hydrologic model as lateral flow is daily surface runoff generation and daily sub-surface runoff
generation data, which are simulated by the land-surface process model embedded in the MRI-AGCM3.1S.
Runoff simulation data of hourly maximum and daily mean are stored for each day with about 1km spatial
resolution. The simulated discharge data were analyzed to discuss the change of the flood risk and water
resources.
6.1 Change of floods

Annual maximum hourly discharge data was compiled and the statistical characteristics wereanalyzed. Figure 13(a) shows the change ratio of the mean of the annual maximum hourly discharge for thepresent climate experiment and the near future climate experiment. Generally, the annual maximum
discharge of the main stream of the Chao Phraya River does not change, however the one of the tributaries
changes location to location. Figure 13(b) shows the change ratio of the 10-year return period riverdischarge in the future climate experiment with respect to the future climate experiment using the Gumbel
distribution. The spatial pattern is similar to the change of the mean of the annual maximum riverdischarge. It is important to recognize that the change of the discharge would appear especially at thetributaries.

(a)

(b)

Figure 13: (a) Change of the ratio of the mean of the annual maximum hourly discharge in the futureclimate experiment with respect to the present climate experiment. (b) Change of the ratio of the 10-years
return period river discharge in the future climate experiment with respect to the present climate
experiment.
13

Figure 14: (a) Change of the ratio of the mean October discharge in the future climate experiment withrespect to the present climate experiment. (b) Change of the monthly discharge at the lower reach of the
Pasak River represented the circle in (a).

6.2 Change of water resources

Figure 14(a) shows the change ratio of the mean October discharge in the future climate experiment withrespect to the present climate experiment. Generally, the monthly discharge on October decreases in the
middle and lower part of the Chao Phraya River basin. Figure 14(b) shows the projected monthly discharge
at the lower reach of the Pasak River. The figure shows the clear pattern of the decrease of the October
discharge.

7 SUMMARY
Distributed hydrological modeling based on the kinematic wave flow was introduced. A distributedhydrologic model is capable of incorporating the basin heterogeneity to model the basin hydrologicalresponses. The paper was organized as follows:
z
z
z
z
z
z

Section 1 explains the basic concept of distributed hydrological modeling.
Section 2 presents a water flow model embedded in a distributed hydrological model. As a typicalflow model, the basic theory of the kinematic wave flow on a rectangular plane is explained.
Section 3 explains numerical solution methods of the kinematic wave flow. A typical explicit method
and an implicit method are introduced.
Section 4 represents extensions of the kinematic wave model for surface and subsurface flow
modeling.
Section 5 shows applications of a distributed hydrological model to basins with different catchment
and climatic characteristics.
Section 6 shows applications of a distributed hydrologic model for future discharge projection underclimate change.
14

REFERENCES
Kitoh, A., Ose, T., Kurihara, K., Kusunoki, S., Sugi, M. and KAKUSHIN Team-3 Modeling Group (2009):
Projection of changes in future weather extremes using super-high-resolution global and regional
atmospheric models in the KAKUSHIN Program: Results of preliminary experiments, HydrologicalResearch Letters, Vol. 3, pp. 49-53.
Sayama, T. and McDonnel, J. J. (2009) A new time-space accounting scheme to predict stream water
residence time and hydrograph source components at the watershed scale, Water Resources Research,
45, pp. 1-14.
Tachikawa Y., Nagatani, G. and Takara K. (2004) Development of stage-discharge relationship equation
incorporating saturated-unsaturated flow mechanism, Annual Journal of Hydraulic Engineering, JSCE,48, pp. 7-12 (in Japanese).
Takasao T and Shiiba M. (1988) Incorporation of the effect of concentration of flow into the kinematicwave equations and its application to runoff system lumping, Journal of Hydrology, 102, pp. 301-322.

REFERENCES BOOKS
Beven, K. J. (2000) Rainfall-runoff Modeling The Primer, Wiley.
Bras, R. L. (1989) Hydrology: An Introduction to Hydrologic Science, Addison-Wesley.
Brutsaert, W. (2005 Hydrology: An Introduction, Cambridge University Press.
Chow, V. T., D. R. Maidment and L. W. Mays (1988) Applied Hydrology, McGraw-Hill.
Eagleson, P. S. (1970) Dynamic Hydrology, McGraw-Hill.
Hornberger, G. M., Raffensperger, J. P., Wiberg, P. L. and Eshleman, K. N. (1998) Elements of PhysicalHydrology, The Johns Hopkins University Press.
Maidment, D. R. (ed.) (1993) Handbook of Hydrology, McGraw-Hill.
Vieux, B. E. (2004) Distributed Hydrogic Modeling Using GIS, second edition, Kluwer Academic
Publishers.

15