Annual Journal of Hydraulic Engineering, JSCE, Vol.54, 2010, February

Annual Journal of Hydraulic Engineering, JSCE, Vol.54, 2010, February

VALIDITY ASSESSMENT OF INTEGRATED KINEMATIC WAVE
EQUATIONS FOR HILLSLOPE RAINFALL-RUNOFF MODELING
Hyunuk AN 1 , Yutaka ICHIKAWA 2 , Kazuaki YOROZU 3 ,Yasuto TACHIKAWA 4 and Michiharu SHIIBA 5
1 Student member of JSCE, Graduate student, Graduate School of Eng., Kyoto University(Kyoto-Daigaku-Katsura, Nishikyo-Ku, Kyoto 615-8540, Japan)

2 Member of JSCE, Dr. Eng., Associate Professor, Interdisciplinary Graduate School of Medicine and Eng.,

University of Yamanashi (4-3-11, Takeda, Kohu 400-8511, Japan)

3 Member of JSCE, Dr. Eng., Asssitant Professor, Faculty of Eng., Kyoto University(Kyoto-Daigaku-Katsura, Nishikyo-Ku, Kyoto 615-8540, Japan)

4 Member of JSCE, Dr. Eng., Associate Professor, Faculty of Eng., Kyoto University
5 Member of JSCE, Dr. Eng., Professor, Faculty of Eng., Kyoto University

The integrated kinematic wave model has been successfully applied to Japanese basin as an element
model of distributed hydrological models in many researches. However, a recent research found that themodel performance could be worse for dry and mildly sloped basins. To assess the model validity forvarious hydrological conditions, we conducted numerical experiments using the integrated kinematic wave
model and a physically-based model combining 2D Richards equation for saturated-unsaturated ﬂow andkinematic wave model for overland ﬂow (2D model). The numerical experiments suggested that in case ofhigher rainfall intensity, thinner soil depth, steeper slope angle and wet initial condition, the simulation
results of the integrated kinematic wave model showed good agreement with those of the 2D model.

Key Words : integrated kinematic wave model, hillslope runoﬀ, distributed hydrological modeling, numerical experiment

1.

INTRODUCTION

is still unclear which hydrological condition eﬀects theapplicability of the integrated kinematic wave model,

In hydrological modeling, the kinematic wave models have been broadly used for river ﬂood simulations.
Takasao & Shiiba1) proposed an integrated kinematicwave model which can handle two types of ﬂows, matrix ﬂow and overland ﬂow. Their model can be applied to hillslope runoﬀ simulation. Tachikawa et al.2)
extended this kinematic wave model and proposed a
new integrated kinematic wave model which can deal
with three diﬀerent ﬂows, unsaturated, saturated, andoverland ﬂows. These integrated kinematic wave models have been successfully applied to Japanese basins
as an element model of distributed hydrological models in many researches. Hunukumbra3) extended thismodel application and applied the models to several
places in the world. He found a tendency that themodel performances were good in wet and steep basinssuch as Japanese basins, while, in dry and mildly
sloped basins, the integrated kinematic wave models could not give well agreed results with observationdata.

because only four diﬀerent places were used in the
study of Hunukumbra3) . To assess a correlation between the model validity and the each hydrologicalconditions, it is required to analyze the model performance in various places with many events. However,
it is diﬃcult to conduct a research using observation
data because of limited available data and uncontrollable nature conditions. In this study, we used numerical experiment to assess the validity of the integrated
kinematic wave model.
In this study, the results simulated by a physicallybased 2D model was considered as a surrogate of observations. The 2D model consists of a 2D subsurface ﬂow model and a 1D overland ﬂow model. Even
though 2D Richards equation cannot perfectly simulate a runoﬀ process of real hillslope in all aspect
because of many uncertainty of soil properties, structures, unknowness of hydrological processes and etc.,
the same model concept has been used in several pre-Hunukumbra’s ﬁnding is obviously important andvious researches for simulating hillslope runoﬀ and it
successfully reproduces the runoﬀ discharge and thecritical in distributed hydrological modelings, but itother characteristics of hillslope. For example, Hopp

- 505 -z(m)Rainfall

No flowD
No flow

Seepage
w

0
0

x (m)
100

Fig.1 Concept of integrated kinematic wave model 2) .

Fig.2 2D model, the subsurface ﬂow domain.

& McDonnell4) performed numerical experiments controlling storm size, slope angle, soil depth and bedrock

ﬁve parameters (n, ka , ds , dc , β). If dc , the water
depth corresponding to maximum unsaturated ﬂow,

permeability to study the eﬀect for hillslope runoﬀ using a model which is based on the same concept withis zero, this equation results in the model proposed by

the 2D model of this study. Keim et al.5) also performed virtual experiments to investigate the processof evaporation and canopy interception using a model
based on the same concept. Hence, the 2D model isconsidered to be useful to evaluate the eﬀect of severalhydrological conditions on runoﬀ discharge.

(2)

Under several conditions with controlling slope angle, soil depth, total rainfall and initial condition, thesimulations were conducted using the 2D model. The
integrated kinematic wave model was calibrated to reproduce the results of the 2D model.

2.

MODEL DESCRIPTION(1)

Integrated kinematic wave model

Takasao & Shiiba1) . Combining Eq. 1 and the following continuity equation describes slope runoﬀ system.
∂h ∂q+= r(t)(2)∂t∂xwhere t is time, x is horizontal dimension, r is rainfall
intensity.
2D modelThe physically-based 2D model consists of a 2D subsurface ﬂow model and a 1D overland ﬂow model.
Subsurface ﬂow is described by 2D Richards equation
as follows



 
∂θ
∂
∂ψ
∂ψ
∂
=
+1
(3)
K
+
K
∂t
∂x
∂x
∂z
∂z
where ψ is the pressure head, θ is the volumetric moisture content, K is hydraulic conductivity, x denotes

proposed an integrated kine-

the horizontal dimension, and z denotes the verticaldimension, assumed to be positive upwards. In thismatic wave model which can handle matrix ﬂow andoverland ﬂow in one system based on relation between

study, a coordinate transformation technique was used
to express non-orthogonal slope as Fig. 2 but a de-water depth and discharge. Tachikawa et al.2) modiﬁed this relation and proposed an extended model

tailed numerical scheme would not be shown in thispaper due to page limitation.

describing three types of ﬂows. The extended model
equation is as follows⎧⎪⎪
v (h/dc )β ,0 ≤ h ≤ dc⎪
⎨ c
q = vc dc + va (h − dc ),dc < h ≤ ds (1)⎪⎪⎪⎩vc dc + va (h − dc ) + α(h − ds )m , ds < h

The upper (x = 100 m) and bottom (z − x sin w =0 m) boundaries were set to be no ﬂow. The

where h is water depth, q is ﬂow discharge, vc (= kc i)ary) as long as ψ is negative, otherwise, it becomes
Dirichlet boundary (ψ = 0). If the inﬁltration rate

Takasao & Shiiba

1)

is velocity of unsaturated ﬂow, va (= ka i) is velocity of√saturated ﬂow, α(= i/n) and β are parameters, andthe shallow rectangular cross-section was assumed. To
keep continuity of ﬂow velocity, Kc = Ka /β is assumed, where Kc and Ka are hydraulic conductivitiesof the unsaturated and saturated ﬂows, respectively.
Fig. 1 shows the concept of this model. There arelower boundary treated as a seepage face. Seepage
face length was controlled automatically according to
Neuman6) . At the ground surface, water can enter thesoil domain at the rainfall intensity (Neumann bound-excesses the rainfall intensity while the surface node is
Dirichlet boundary, it turns Neumann boundary condition again. This control is conducted in iterative
procedure.
When the surface condition is saturated, the inﬁltration rate was calculated by the subsurface ﬂow- 506 -

rainfall intensity (mm/hr)5

Table.1 Simulation condition

4
3

slope angle

5, 20, 35 degree

soil depth
total rainfall

0.25, 0.5, 1.0, 2.0 m
10, 20, 40, 70, 100 mm

initial condition

wet, dry

2
1
0
0

2

4

6

8

10

time (hr)Fig.3 Rainfall intensity, total rainfall is 10 mm.

(2)Model parameter

The equation of Van Genuchten9) for the soil watermodel. The amount of diﬀerence between the rainfall
intensity and the calculated inﬁltration rate is the lateral inﬂow for the overland ﬂow from soil domain andprecipitation. Otherwise, when the surface condition
is Neumann boundary, the inﬁltration rate is same as
rainfall intensity and the lateral inﬂow for the overland ﬂow would be zero. The overland ﬂow equation
is written as follows∂h ∂q+= ql (t)∂t∂xq = αhm

retention curve and that of Mualem10) for the unsaturated hydraulic conductivity function were used forthe subsurface ﬂow model. The soil water retention
curve is given by1−1/n1θ − θrSe ==(6)
θs − θr
1 + (α|ψ|)n
where Se is the eﬀective saturation; θr and θs are theresidual and saturated water contents, respectively;(4)and α and n are Van Genuchten parameters, whose
values depend upon the soil properties. Based on the(5)model of Mualem10) , the unsaturated hydraulic conductivity function is given by

where ql is lateral ﬂow rate into surface ﬂow from sub-K = Ks Se1/2 1 − (1 − Sen/(n−1) )1−1/n2(7)surface ﬂow and precipitation which is calculated by
the subsurface ﬂow model. The 2D subsurface ﬂow

where Ks is the saturated hydraulic conductivity. The

model and the 1D overland ﬂow model are combined
explicitly.

parameter values were set as θs = 0.475, θr = 0.28,
Ks = 2.5 m/s, α = 4 m−1 , n = 2 by referring to Hopp
& McDonnell 4) .

3.

NUMERICAL EXPERIMENT(1)

Simulation condition

As the model of Tachikawa et al.2) includes that ofTakasao & Shiiba1) , the former one was used as the integrated kinematic wave model in this study. The ﬁve

Under sev-

parameters (n, ka , ds , dc , β) of the integrated kinematic wave model were calibrated for each simulations

eral conditions, the simulations were conducted using the 2D model and the integrated kinematic waveof the 2D model by using the SCE algorithm7) . The
Nash Sutcliﬀe coeﬃcient8) was used as the objectivemodel was calibrated to reproduce the results of the
2D model as better as possible.

function of parameter optimization.

We considered the results simulated by the 2Dmodel as a surrogate of observations.

A simple slope was considered as Fig. 2. The slope

4.

RESULT AND DISCUSSION

length was 100m. The slope angle and soil depth werevaried in simulation as Table. 1. Rainfall intensity was(1)

Runoﬀ discharge

given as Fig. 3 when total rainfall is 10 mm. When
total rainfall is 20, 40, 70, and 100 mm, the rainfall in-

Fig. 4 describes four examples of runoﬀ discharge
simulated by the 2D model and the integrated kine-

tensity was multiplied by 2, 4, 7 and 10, respectively.
Two types of initial conditions were considered. A

matic wave model. The left column of the ﬁgure shows
the results of the 2D model and the right column

wet initial condition was simulated as soil moisture
contents and the corresponding pressure head ﬁeld af-

shows those of the integrated kinematic wave model.
When soil depth is 0.25 m and 1 m, the runoﬀ dis-

ter three days drainage from a saturated soil domain
without rainfall. A dry initial condition was obtained

charges simulated by the two models are well matched.
On the other hand, the discharge patterns are compar-

after a week drainage.

atively not agreed each other with 2m soil depth. The- 507 -

5

5 degree, wet
20 degree, wet
35 degree, wet
5 degree, dry
20 degree, dry
35 degree, dry

43discharge rate (mm/hr)

discharge rate (mm/hr)

5100mm, D = 2 m
2D model

2
1
0
1

2
Time (days)

3

100mm, D = 2 m
Kinematic wave model

2
1

4

0

20
15

100mm, D = 0.25 m
2D model

10

1

30

5 degree, wet
20 degree, wet
35 degree, wet
5 degree, dry
20 degree, dry
35 degree, dry

25discharge rate (mm/hr)

discharge rate (mm/hr)

30

5
0

2
Time (days)

3

4

5 degree, wet
20 degree, wet
35 degree, wet
5 degree, dry
20 degree, dry
35 degree, dry

25
20
15

100mm, D = 0.25 m
Kinematic wave model

10
5
0

0

1

1.8
1.6

2
Time (days)

3

4

0

1.4
1.2
1
0.8
0.6

40mm, D = 1 m
2D model

0.4
0.2
0
0

1

1.8

2
Time (days)

3

4

3

4

1.2
1
0.8

40mm, D = 1 m
Kinematic wave model

0.6
0.4
0.2
0
1

2

3

4

Time (days)
1.8
discharge rate (mm/hr)

0.8
0.6

2
Time (days)

5 degree, wet
20 degree, wet
35 degree, wet
5 degree, dry
20 degree, dry
35 degree, dry

1.6
1.4

0

5 degree, wet
20 degree, wet
35 degree, wet
5 degree, dry
20 degree, dry
35 degree, dry

1.6
1.4
1.2
1

1

1.8

5 degree, wet
20 degree, wet
35 degree, wet
5 degree, dry
20 degree, dry
35 degree, dry

discharge rate (mm/hr)

discharge rate (mm/hr)

300discharge rate (mm/hr)

5 degree, wet
20 degree, wet
35 degree, wet
5 degree, dry
20 degree, dry
35 degree, dry

420mm, D = 2 m
2D model

0.4
0.2
0

1.6
1.4
1.2
1

5 degree, wet
20 degree, wet
35 degree, wet
5 degree, dry
20 degree, dry
35 degree, dry

0.8
0.6

20mm, D = 2 m
Kinematic wave model

0.4
0.2
0

0

1

2

3

4

0

Time (days)

1

2

3

4

Time (days)

Fig.4 Hydrographs simulated by 2D model and integrated kinematic wave model.

peak time of the integrated kinematic wave model appeared earlier than that of the 2D model with 2m soil
depth. This diﬀerence is supposed to be arised fromvertical inﬁltration which cannot be treated in the integrated kinematic wave model. However, the overall
runoﬀ pattern was well reproduced entirely. The slopeof 5 degree results in earlier peak ﬂow than the otherssometime. This is caused by the diﬀerence of initialcondition which of 5 degree slope is wetter that theothers. As previously mentioned in 3.(1), the initial conditions were made by simulating three days
drainage from a saturated soil domain without rainfall. Hence, the initial conditions of 5 degree slope isgenerally wetter than the other two conditions.

(2)

Nash Sutcliﬀe coeﬃcientFig. 5 shows the Nash Sutcliﬀe coeﬃcients with
diﬀerent soil depths. If the integrated kinematic wave
model perfectly reproduces a runoﬀ discharge pattern
simulated by the 2D model, the corresponding Nash
Sutcliﬀe coeﬃcient becomes 1. Obviously, Nash Sutcliﬀe coeﬃcients become worse as soil depth becomes
thicker. It could be arised from the eﬀect of vertical inﬁltration. As soil depth is thicker, vertical inﬁltration
gives larger eﬀect to runoﬀ process. Because the integrated kinematic wave model does not have a mechanism describing vertical inﬁltration in soil domain, it
fails to reproduce the runoﬀ discharge calculated by
the 2D model which can express eﬀect of vertical inﬁltration in thick soil condition. Furthermore, Nash
Sutcliﬀe coeﬃcients tend to be higher with wet initial- 508 -

0.8

0.8

0.6

0.6

0.4

0.4
NSE

(a) 1

NSE

(a) 1

0.2

0.2

0

0

-0.2

-0.2

10 mm
20 mm
-0.4 40 mm
70 mm
100 mm
-0.6
0
0.5

-0.4
-0.6

1

1.5

10 mm
20 mm
40 mm
70 mm
100 mm0

2

Soil depth (m)

10

40

0.8

0.8

0.6

0.6

0.4
NSE

0.4
NSE

40

(b) 1

(b) 1

0.2

0.2
0

0

-0.2

-0.2

10 mm
20 mm
40
mm
-0.4
70 mm
100 mm
-0.6
0
0.5

10
20
30
Slope angle (degree)

-0.4
-0.6

10 mm
20 mm
40 mm
70 mm
100 mm
0

1

1.5

2

20

30

Slope angle (degree)

Soil depth (m)Fig.5 Nash Sutcliﬀe coeﬃcient with diﬀerent soil depths:(a) wet, (b) dry initial condition.

Fig.6 Nash Sutcliﬀe coeﬃcient with diﬀerent slope angles: (a) wet, (b) dry initial condition.

posed to be same as that of the initial condition eﬀect.
conditions and vice versa with dry initial condition.
This reason also might be two dimensional ﬂow eﬀect.
As the initial condition is wetter, the soil domain iseasily and quickly saturated and the vertical inﬁltration gives relatively small eﬀect to runoﬀ process.

In high rainfall intensity, vertical inﬁltration gives less
eﬀect to overall runoﬀ process.
(3)

Discharge components

Fig. 6 shows the Nash Sutcliﬀe coeﬃcients withAs previously mentioned, if dc = 0 in the model
proposed by Tachikawa et al.2) , the model results in

diﬀerent slope angels. It seems like that there is noapparent correlation between slope angle and validitythat proposed by Takasao & Shiiba1) .
Fig. 7 shows separate components of discharges ofof the integrated kinematic wave model. However, ex-three models with 100 mm total rainfall, D = 0.25

cept considerably wrong results (Nash Sutcliﬀe coeﬃcient less than 0.7), reproducibility tends to be betterm, 5 degree sloep angle and wet initial condition. It
can be found that the 2D model and the model pro-with larger slope angle. This might be arised fromthat ﬂow condition of the lower part does not aﬀect

posed by Tachikawa2) give diﬀerent discharge components even though the overall discharges are wellthat of upper part in the kinematic wave model because it assumes supercritical ﬂow.

agreed each other. The two integrated kinematic wave
models gave diﬀerent results and the hydrograph sim-

These two characteristics shown in Fig. 5 and 6
about validity of the integrated kinematic wave modelulated by Takasao & Shiiba’s model1) is closer to the
2D model’s result. However, as Takasao & Shiiba’s

are well agreed with the research of Hunukumbra3) .
He found that the performance of the integrated kine-model1) can be considered as a variant of Tachikawa’s
model2) , the result simulated by Takasao & Shiiba’s

matic wave model was good in wet and steep basinswhile it becomes worse in dry and mildly slopedmodel1) can be considered as the result of Tachikawa’s
model2) . This is a typical parameter identiﬁcation

basins.
Another point is that the validity is better with highproblem. It means that another information is required to reproduce discharge components using therainfall intensity in Fig. 5 and 6. The reason is sup-integrated kinematic wave model, e.g. speciﬁc param-- 509 -

discharge rate (mm/hr)

35

corresponds with the study of Hunukumbura3) . In-

total flow
matrix flow
overland flow

30

cluding vertical inﬁltration mechanism, the integrated
kinematic wave model may give better result for dry

25

100mm, D =0.25m
5 degree, wet
2D model20
15

and mildly sloped basins.
Another important ﬁnding is that if the integrated

105kinematic wave model perfectly reproduces the over-

0
0

0.5

1
Time (days)

discharge rate (mm/hr)

35

1.5

2

all discharge from hillslope, it does not guarantee thegood agreement with respect to discharge components.

total flow
matrix flow
overland flow

30

To obtain reliable result of separated discharge components using the integrated kinematic wave model,

25

100mm, D = 0.25 m
5 degree, wet
Kinematic wave model

20
15

the other information which limits the range of model
parameter should be needed in terms of hydrological

proposed by Tachikawa

10
5characteristics of the study basins.

0
0

0.5

1
Time (days)

discharge rate (mm/hr)

35

1.5

2

REFERENCES

total flow
matrix flow
overland flow

30
25100mm, D = 0.25 m
5 degree, wet
Kinematic wave model

20
15

proposed by Takasao & Shiiba

10
5
0
0

0.5

1
Time (days)

1.5

2

Fig.7 Matrix and overland ﬂows simulated by the 2Dmodel and the integrated kinematic wave model
proposed by Tachikawa et al.2) and Takasao &Shiiba1) (100 mm, D = 0.25 m, 5 degree, wet).

eter range or discharge rate of speciﬁc component.

5.

CONCLUSION

To assess the validity of the integrated kinematicwave model, we conducted numerical experiments using the physically-based 2D model and the integrated
kinematic wave model. The results simulated by the
2D model were assumed to be a surrogate of observations and the kinematic wave model was calibrated
to reproduce the result of the 2D model. It was supposed that rainfall intensity, soil depth, slope angle,
and initial condition could eﬀect to the validity of theintegrated kinematic wave model and the numericalexperiments were conducted with changing these factors. According to the results of the numerical experiments, in condition of higher rainfall intensity, thinner

1) Takasao, T. and Shiiba, M.: Incorporation of the effect of concentration of ﬂow into the kinematic wave
equations and its applications to runoﬀ system lumping, Journal of Hydrology, Vol.102, pp.301–322, 1988.
2) Tachikawa, T., Nagatani, G. and Takara, K.: Developement of stage-discharge relationship equation
incorporating saturated-unsaturated ﬂow mechanism,
Annual Journal of Hydraulic Engineering, JSCE,
Vol.48, pp.7–12, 2004. (in Japanese)
3) Hunukumbura, J.M.P.B.: Distributed hydrologicalmodel transferability across basins with diﬀerent
physio-climatic characteristics, Phd thesis, Kyoto
University, 2009.
4) Hopp, L. and McDonnell, J.J.: Connectivity of hillslope scale: identifying between interactions storm
size, bedrock permiability, slope angle and soil depth,
Journal of Hydrology, Vol.376, pp.378–391, 2009.
5) Keim, R.F.: Tromp-van Meerveld, H.J. and McDonnell, J.J., A virtual experiment on the eﬀects of evaporation and intensity smoothing by canopy interception on subsurface stormﬂow generation, Journal of
Hydrology, Vol.327, pp.352–364, 2006.
6) Neuman, S.P.: Saturated–unsaturated seepage by ﬁnite elements, Journal of Hydraulics Division, ASCE
HY12, pp.2233–2250, 1973.
7) Duan, Q., Sorooshian, S. and Gupta, V.K.: Eﬀective and eﬃcient global optimization for conceptual
rainfall-runoﬀ models, Water Resources Research,
Vol.28(4), pp.1015–1031, 1992.
8) Nash, J.E. and J.V., Sutcliﬀe: River ﬂows forecasting through conceptual models, part 1 a discussion ofprinciples, Journal of Hydrology, Vol.27(3), pp.282–
290, 1970.
9) Van Genuchten, M.: A closed-form equation forpredicting the hydraulic conductivity of unsaturated
soils, Soil Science Society of America, Vol.44, pp.892–
898, 1980.
10) Mualem, Y.: A new model for predicting hydraulicconductivity of unsaturated porous media, Water Resources Research, Vol.12, pp.513–522, 1976.

soil depth, steeper slope angle, and wet initial condition, the integrated kinematic wave model could give
good agreement with the 2D model. This analysis- 510 -

(Received September 30, 2009)