Annual Journal of Hydraulic Engineering, JSCE, Vol.52, 2008 February

KINEMATIC WAVE FLOW MODELS FOR
RIVER BASIN RUNOFF SIMULATION
Michiharu SHIIBA1 , Yasuto TACHIKAWA2 and Yutaka ICHIKAWA3
1

Member of JSCE, Dr. of Eng., Professsor, Dept. of Urban & Environmental Engineering, Kyoto
University (Nishikyo-ku, Kyoto City, 615-8540, Japan)
2
Member of JSCE, Dr. of Eng., Associate Professsor, Dept. of Urban & Environmental Engineering,
Kyoto University (Nishikyo-ku, Kyoto City, 615-8540, Japan)
3
Member of JSCE, Dr. of Eng., Assistant Professsor, Dept. of Urban & Environmental Engineering,
Kyoto University (Nishikyo-ku, Kyoto City, 615-8540, Japan)

Kinematic wave ﬂow models for simulating river basin runoﬀ are reviewd. After the baic equations of kinematic wave ﬂow are derived, some modiﬁed kimematic wave ﬂows which are devisedfor treating surface-subsurface ﬂow on watersheds are introduced. Then, the routing of kinematicwave ﬂows on digital terrain models are inroduced. Finally, lumping methods of kinemtatic wave
ﬂows for watershed runoﬀ are shown.
Key Words : Kinematic Wave Flow, Surface-Subsurface Flow, Lumping

1.

Introduction

runoﬀ.

Rainfall-runoﬀ processes are water ﬂow phenomena that occur on and beneath the groundsurface, and can be modeled using diﬀerential
equations to describe movement of ﬂowing water.
The kinematic wave ﬂow model is one of rainfallrunoﬀ models which have been constructed basedon such a notion.
By applying the kinematic wave ﬂow model forsurface ﬂow on mountaneous slopes, Ishihara andTakasao1) developed the theory of surface runoﬀ
and initiated modern physically-based hydrology.
In their theory, the discharge of surface ﬂow is assumed to be expressed by the power-law funcionof the depth of the ﬂow. However, as was shown
later by themselves2) , there exists saturated overland ﬂow which appears when the subsurface ﬂow
which ﬂows through higly permeable soil layers on
montain slopes comes over the layer. To express
such ﬂow, it is convenient to express both thesaturated overland ﬂow and the subsurface ﬂow
in an integrated fashion. Thus in Section 2, weintroduce more general form of kinematic wave
ﬂow model and its solution procedure. In Section
3, the applicaions of those general form of kinematic ﬂow model to modeling watershed runoﬀsare introduced. Section 4 explains the routing of
the kinematic wave ﬂow on digital terrain modelsof river basins. Finally, Section 5 shows a lumping method of kinematic wave ﬂows for watershed

2.

BASIC EQUATIONS OF THE
KINEMATIC WAVE FLOWThe kinematic wave model is an unsteady open
channel ﬂow wave model which is often usedin engineering hydrology. The kinematic wave
model was originally developed for river routing,but later was applied to catchment rainfall-runoﬀ
modeling and several variations of the kinematicwave model were developed.
The basic form of the kinematic wave equationis∂A ∂Q+= qL (x, t)(1)∂t∂xQ(x, t) = αA(x, t)m(2)
where A(x, t) is the ﬂow cross-sectional area,Q(x, t) is the ﬂow discharge, qL (x, t) is the lateral inﬂow discharge per unit length, and α andm are constants. A(x, t) and Q(x, t) are abbreviated to A and Q, respectively in Eq.(1).
Equation (2) is derived from Manning’s or
Chezy’s laws, which are ﬂow resistance laws ofopen channel uniform ﬂow. Manning’s law is described as:
1v = R2/3 sin θ1/2(3)nin which v is the average ﬂow velocity, n is Manning’s friction coeﬃcient, R is the hydraulic radius, and θ is the bed slope. By assuming that RK-1

Dlaw and Manning’s law:{ah (0 ≤ h < d)q=α(h − d)m + ah (h ≥ d)rx

H

where a (= k sin θ/γ) is the velocity of water
which ﬂows in the soil layer, k√is the hydraulicconductivity, and α and m are sin θ/n and 5/3,respectively. The upper relation given in Eq.(7)is used for the case where only subsurface ﬂow
occurs, and the lower relation in Eq.(7) is used
for the case where saturation excess overland ﬂow
occurs.
Shiiba et al.6) developed a new kinematic wave
model by which both ﬂood and low ﬂow simulation can be conducted by introducing a concept of ﬁeld capacity of hillslope soil layer tothe kinematic wave model for surface-subsurface
ﬂows. This new kinematic wave model enables
a runoﬀ simulation to be connected continuously
from ﬂood period to low ﬂow period, or conversely
from low ﬂow period to ﬂood period. This model
assumes that a part of water which inﬁltrates intoa surface soil layer is mainly captured by the capillary force of soil particles and lateral water ﬂow
occurs when water content exceeds the ﬁeld capacity of the soil.
The kinematic wave equation considering ﬁeld
capacity is as follows:{ahf (0 ≤ hf < d)(8)q=α(hf − d)m + ahf (hf ≥ d)θLFig. 1 Schematic
representation
subsurface ﬂow on a hillslope

ofsurface-can be generally expressed as z1 Az2 (z1 and z2 areconstants), we obtain the expression of Eq.(2).
By applying the kinematic wave approximation
to the sheet ﬂow on a plane, we can get also
the kinematic wave ﬂow model of surface ﬂow onmountain slopes.
∂h ∂q+= r(x, t)
(4)
∂t
∂x
q = αhm(5)in which h is the water depth of the surface ﬂow,
q is the discharge per unit width of the slope, ris the lateral input.
It is possible to use the more general form ofdischarge and depth relationq = f (h, x)(6)(() )1/N
d − hf N hf + hc 1 −dh=(0 ≤ hf < d)
hf + hc (hf ≥ d)as we show in Section 3.

3.

(7)

SURFACE-SUBSURFACE FLOW
MODEL AND CONSIDERATION
OF FIELD CAPACITY(9)

where q is the discharge per unit width, h is thetotal water content, hf is the free water content,a (= k sin θ/γe ) is the water ﬂow velocity in asurface soil layer, k is the hydraulic conductivity,
γe is the eﬀective porosity of a surface soil layer,d (= γe D) is the depth of pore contained in asurface soil layer,
D is the depth of a surface soil
√layer, α = sin θ/n, θ is the slope gradient, n
is the Manning roughness coeﬃcient, m is a constant ( > 1), N is a constant ( > 1), hc = γc D,
and γc is the ﬁeld capacity of a surface soil layer.

Takasao et al.4),5) developed a kinematic wave
model which considers both surface and subsurface ﬂow, and its interaction. This model candescribe a saturation excess overland ﬂow mechanism.
They assumed that hillslope surface was covered by a soil layer with high surface inﬁltration
capacity and permeability, which is illustrated inFig. 1. In this ﬁgure, D is the depth of the soillayer, L is the length of the hillslope, H is thedepth of water ﬂow, and r is the rainfall intensity. Let γ be the porosity of the soil layer andbe spatially constant. Then the depth of the eﬀective 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 et al. derived the following
piecewise relation between h and q from Darcy’s

4.

ROUTINGOFKINEMATICWAVE FLOW ON DIGITAL TERRAIN MODELSIn recent years, digital elevation maps havebecome available at a resolution ﬁne enough toclearly describe the topography of catchments.

K-2

B

the computational burden of water ﬂow routing.
These lumped models have a physical basis andthe parameter values are calculated using topographic variables obtained from the DEM.
Owing to the limit of space, only the procedure of the lumping of the power-law kinematicwave model is introduced. In this section, thekinematic wave model, A = KQP (A: ﬂow crosssectional area, Q: ﬂow discharge, K and P : constants), is used to describe shallow water ﬂow on
a rectangular plane instead of Eq.(5).
The main purpose of this section is to derivea steady-state relation between storage volumeand outﬂow discharge of a rectangular plane byspatially integrating the kinematic wave model.
When a rainfall-runoﬀ system reaches a steadycondition, the discharge at some point, Q, is expressed as r × M , in which r is the rainfall intensity and M is the upslope contributing area atthe point.
Suppose that a catchment is represented by N
slope units. Let us consider the ith slope unit(i = 1, 2, ..., N ). Let Li be the length, Bi thewidth, Mi the area, θi the gradient, and Ki and Pi
the kinematic constants K and P , of the ith slopeunit. Also let Ui be the upslope contributing areaat the upper end of the ith slope unit. If the ith
slope unit is located on the upper boundary of acatchment, Ui = 0.
Let Ai (x) be the ﬂow cross-sectional area andQi (x) the ﬂow discharge at the distance x fromthe upper end of the ith slope unit. Since a slope
unit is a rectangular plane, the upslope contributing area at the distance x from its upper end of
the ith slope unit, Mi (x), is expressed as:

slope unit

flow direction

Fig. 2 Schematic representation of the grid-based
DEM

Using these data, a number of digital elevation
models (DEMs) have been developed to precisely
and eﬃciently represent catchment topography.
Tachikawa et al.7) developed a distributed hydrological model based on a triangular irregularnetwork DEM.
Ichikawa et al.8) applied the kinematic wave
model considering ﬁeld capacity to a grid-basedDEM9) . Fig. 2 shows a schematic representation of the grid-based DEM which they used.
In this ﬁgure two adjacent grid points are connected by a rectangular plane. They call thisrectangular plane “a slope unit”. Catchment topography is represented by a set of slope units.
For each slope unit, its area, length and gradient can be easily calculated. They applied thekinematic wave model considering ﬁeld capacity
to all the slope units comprising the catchment
topographic model and routed water ﬂows fromupstream slope units to downstream slope units.
The outﬂow discharge from the catchment wasused as the lateral inﬂow condition of a river
routing model. From the result of the simulaion,
Ichikawa et al. showed the fact that the kinematic wave model considering ﬁeld capacity hashigh applicability to low ﬂow periods.

Mi (x) =Then we can obtain(Qi (x) = r
and5.

Mi
x + UiLi
Mi
x + Ui
Li(LUMPING KINEMATIC WAVEFLOW EQUATIONS

Ai (x) = Ki r

Pi(10))

Mi
x + Ui
Li(11))Pi(12)Spatially integrating Ai (x), we can obtain thestorage volume of the ith slope unit, si
∫ Li
si =
Ai (x)dx

The distributed runoﬀ models, such as the kinematic wave models, have the advantage that theyare based on the physical processes of ﬂowing water. However, generally speaking, they require
much computation time and resources. In thissection, lumped kinematic wave runoﬀ models areintroduced. They are derived by spatially integrating the kinematic wave models. The resulting lumped kinematic wave runoﬀ models reduce

0=}rPi Ki Li {(Mi + Ui )Pi +1 − (Ui )Pi +1
Pi + 1 Mi(13)

Summing up si , we can obtain the total storage

K-3

volume of the catchment, S =∑N

6.

i=1 si as

CONCLUDING REMARKSN∑

Mainly the researches done by the authors{}
Ki Li(Mi + Ui )Pi +1 − (Ui )Pi +1 rPi about the kinematic wave ﬂow were introduced.
(Pi + 1)Mii=1
Ichikawa12) , written in English, covers the top(14)ics introduced in this paper.
If we assume that the constant in the kinematicwave model, Pi , is spatially uniform and rewrite
REFERENCES
Pi as P , Eq.(14) can be further simpliﬁed as:

S=S = ErP

1) Ishihara, T. and Takasao, T: Fundamental researches on the unit hydrograph method and its
application, Transactions of the Japan Soc. ofCivil Eng., No. 60, 3-3, pp. 1–34, 1959.
2) Ishihara, T., Ishihara, Y, Takasao, T and Rai,C.: A study on the runoﬀ characteristics of YuraRiver Basin, Dis. Prev. Res. Inst., Kyoto Univ.
Annuals, No. 5, pp. 147-173, 1962.
3) Ikebuchi, S., Shiiba, M., Takara, K. andTachikawa, Y.: Hydrology, Asakura Shoten,
Chapter 7, 2006 (in Japanese).
4) Takasao, T., Shiiba, M. and Kitamura, H.: A
distributed runoﬀ model and lumping scale of ariver basin, Proc. 31st Annu. Conf. JSCE, 2-86,163, 1976 (in Japanese)5) Takasao, T. and Shiiba, M.: Incorporation of theeﬀect of concentration of ﬂow into the kinematicwave equations and its applications to runoﬀ system lumping, J. Hydrol., 102, pp. 301–322, 1988.
6) Shiiba, M., Tachikawa, T., Ichikawa, Y., Hori,
T. and Tanaka, K.: Development of slope runnoﬀ models which consider ﬁeld capacity and pipe
ﬂow, Annuals of Disas. Prev. Res. Inst. Kyoto
Univ., 41, B-2, pp. 229–235, 1998 (in Japanese).
7) Tachikawa, Y., M. Shiiba, and Takasao, T.: Development of a basin geomorphic information system using a TIN-DEM data structure, Wat. Resour. Bull., 30, 1, 9–17, 1994.
8) Ichikawa, Y., M. Murakami, Y. Tachikawa andShiiba, M.: Development of a basin runoﬀ simulation system based on a new digital topographic
model, Journal of Hydraulic, coastal and environmental engineering, 691/II-57, 43–52, 2001 (inJapanese).
9) Shiiba, M., Y. Ichikawa, T. Sakakibara andTachikawa, Y.: A new numerical representation
form of basin topography, Journal of Hydraulic,coastal and environmental engineering, 621/II47, 1–9, 1999 (in Japanese).
10) Ichikawa, Y., T. Ogura, Y. Tachikawa, M. Shiiba
and Takara, K.: Lumping general kinematic wave
equation of slope runoﬀ system, Annual Journal of Hydraulic Engineering, JSCE, 44, 145–150,
2000 (in Japanese).
11) Ichikawa, Y. and Shiiba, M.: Lumping of kinematic wave equation considering ﬁeld capacity,
Proc. 3rd International Conference on Water
Resources and Environment Research, Dresden,Germany, Vol. 1, 56–60, 2002.
12) Ichikawa, Y.: Distributed Runoﬀ Modeling, Chap
10, ”Preicipitation and Water Resources - The
Textbook for Twelfth IHP Training Course in2002, pp. 171–201, 2003.

(15)

where}
1 ∑ Ki Li {E=(Mi + Ui )P +1 − (Ui )P +1P +1Mii=1(16)Let O be the outﬂow discharge from the catchment.
∑ Under the steady-state condition, O =r× N
i=1 Mi , rearranged for,NOr = ∑N

i=1 Mi(17)

Substituting Eq.(17) into Eq.(15), we can obtain)P(O= F OP
(18)
S = E ∑N
Mii=1∑P
where F = E/( Ni=1 Mi ) .
Equation (18) gives the steady-state relation
between the storage volume and outﬂow discharge of a catchment. However we use this relation to calculate outﬂow discharge for the unsteady condition, that is, time variant rainfall
case. Thus we have the relationship between storage and outﬂow as follows:
S(t) = F O(t)P(19)The continuity equation of storage volume isgiven by()dS
S(t) 1/P= QL (t) − O(t) = QL (t) −(20)dtF∑where QL (t) = r(t) × N
i=1 Mi . Since Eq.(20)
is an ordinary diﬀerential equation for S(t), we
numerically solve this equation by using RungeKutta method. We call the model derived here as“the lumped kinematic wave model”. The valueof F in Eq.(18) can be calculated from the topographic quantities and kinematic constants. Thismeans that the lumped kinematic wave model isa lumped model but the parameters can be speciﬁed on a physical basis.
The procedure of lumping the surfacesubsurface model is shown in Ichikawa et al.10) ,and the procedure of lumping the kinematicwave model considering ﬁeld capacity is shown
in Ichikawa et al.11) .

(Received January 31, 2008)

K-4