Annual Journal of Hydraulic Engineering, JSCE, Vol.53, 2009, February
Annual Journal of Hydraulic Engineering, JSCE, Vol.53, 2009, February

EXTENDING A STORAGE-DISCHARGERELATIONSHIP FOR SUBSURFACE FLOW MODELINGIN DRY MILD-SLOPE BASINS
Hunukumbura P.B.1, Yasuto TACHIKAWA2 , Yutaka ICHIKAWA3
and Michiharu SHIIBA4
1Student Member, Graduate student, Dept. of Urban and Environmental Eng., Kyoto University(Kyoto 615-8540, Japan) E-mail: hunu@kyoto.mbox.media.kyoto-u.ac.jp
2Member of JSCE, Dr. Eng., Associate Professor, Dept. of Urban and Environmental Eng., Kyoto University(Kyoto 615-8540, Japan)
2Member of JSCE, Dr. Eng., Assistant Professor, Dept. of Urban and Environmental Eng., Kyoto University(Kyoto 615-8540, Japan)
3 Member of JSCE, Dr. Eng., Professor, Dept. of Urban and Environmental Eng., Kyoto University(Kyoto 615-8540, Japan)

In this study, a new storage-discharge equation for unsaturated zone considering the soil moisturesuction is proposed for dry mild-slope basins. For wet steep basins a distributed hydrologic modelwith kinematic wave approximations including surface-subsurface runoff, OHDIS-kwmss has been wellstudied. However, when we apply this model for a dry mild-slope basin, the Illinois basin (2,400 km2) in
U.S, we found that the model predictions were not good. Especially, we observed that the model was notcapable to reproduce the hydrograph well even for the calibrated event itself by parameter tuning . Soil
moisture suction effect is important in dry mild-slope basins and the present equation does not include theeffect of soil moisture suction. Therefore, we derive a new equation by including the soil moisture suction
effect. Finally, the new equation is tested for different hillslopes and it is found that the new equationshows different flow characteristics from the existing equation.
Key Words: subsurface flow, distributed hydrological model, storage-discharge relationship, mild- slope
basins, OHDIS-kwmss, prediction in ungauged basins

1. INTRODUCTION
Understanding the governing hydrologicalprocesses in different regions of the world isnecessary to realize predictions in ungauged basins.
In this regards, physically realistic models can playa major role. However, in practice simplicity of theprocess is also needed because even at hillslope
scales complicated models such as three
dimensional Richards equation are not feasible1).
Furthermore, measuring physical properties of soils,
vegetation etc. in finer spatial and temporal
resolutions is not realistic. Therefore, lumping
processes in some extent is necessary for basin scale
applications. Kinematic approximation is one of thesimplified approaches, which is used to model thehydrological process in both surface and subsurface
flow2),3),4),5).
OHDIS-kwmss is a distributed hydrologicalmodel with kinematic wave approximations
including surface-subsurface runoff developed at
Kyoto University3),6). In general, Japanese basins arewet steep and this model has been successfully
applied for many basins in Japan. Hunukumbura etal.7) applied the model to the Maruyama River basin(909 km2) in Japan and found that the model
calibrated for small flood event can be used to
predict large flood events as well with goodprediction accuracy. Further they concluded that thecalibrated model using the runoff data at the basin
outlet is capable of predicting the inside locations
with reasonable accuracy. However, this model is
not tested in other places in the world. Applicationsof hydrological models and testing their capabilities
in different basins in the world are vital tounderstand the governing hydrological processes ineach region.
In this study, we investigate the prediction
capability of OHDIS-kwmss model for the Illinois
River basin in U.S, which is one of the experimental
basins of the Distributed Model Inter-comparison
Project, DMIP 8),9). In DMIP phase 1, seven
distributed hydrological models (SWAT, NOAH
Land Surface model, HRCDHM, HL-RMS, VIC-

- 25 -

3L, TOPNET and WATFLOOD) were applied to
the Illinois basin at Tahlequah. Results show thatthe overall Nash coefficient for each model varies
from 0.27 to 0.85 for the calibration period. Eventabsolute error for all the models for the same period
is greater than 20%. The best Nash-coefficient value
for the calibration period (0.85) was obtained fromthe HL-RMS model with Sacramento Soil Moisture
Accounting Model (SAC-SMA), and hillslopechannel routing employs the kinematic wave. Thismodel uses larger spatial resolution (16 km2) grids10),11),12).
Our experience in OHDIS-kwmss in wet steep
basins in Japan shows higher Nash-coefficient value(greater than 0.95) for model calibration period. As
the model was not tested yet in dry mild-slope
basins, main objectives in this study are to check
whether the present model structure is good enoughto capture the hydrological process in dry mildslope basin and if it is not sufficient, which
modification is necessary to make good prediction
in dry mild-slope basins.
The model is calibrated for the basin by using
SCE-UA algorithm13),14) and after applying thecalibrated model parameters for other events, it isobserved that the model predictions are not good in
the Illinois River basin. Even for the calibrated
event, it was not possible to reproduce the observed
outflow hydrograph well by parameter tuning. It
means that the present model equations representing
the surface and subsurface flow with kinematicwave approximation, is not sufficient to capture thegoverning hydrological processes in the basin. The
Illinois River basin is mild-slope dry basin, andhence the effect of soil moisture suction in theunsaturated zone would be significant.
In this study a new storage-discharge equation
for unsaturated zone considering the soil moisturesuction is proposed. The new equation is based on
the Darcy’s equation, and soil moisture suction headterm is replaced using Campbell’s simplified
model15) representing the soil moisture andhydraulic conductivity characteristics. The equationis possible to include the effect of soil moisturesuction by adding only one additional parameter tothe existing model. Finally, the new storagedischarge equation is tested for single hill slopes
with different parameters and the slope flow
characteristics of the new equation is compared withthe existing equation.
Section 2 describes the existing model used inthis study and its applications to the Illinois basin.
Section 3 presents the derivation of the new
equation and numerical solution. Application andcomparison of the existing and new equation are
given in section 4.

2. EXISTING MODEL APPLICATIONS(1) The existing storage-discharge equation
OHDIS-kwmss is based on the one dimensional
kinematic wave theory and developed by Ichikawa
et al.6). The basin topography in the hydrologic
model is represented according to the methodology
described in Shiiba et al.16). In the distributedhydrologic model, it is considered that the basin
consists of number of rectangular slope elements
which drains to the deepest gradient of its
surrounding.
It is assumed that the each slope element of thebasin is covered with a permeable soil layer and astorage-discharge relationship for unsaturated,saturated and overland flow defined by the threshold
dm and da3) is adopted. The stage discharge
relationship used in the OHDIS-kwmss model are
given in Eq.(1). In this equation, the soil moisturesuction head in the unsaturated zone is not included.
Therefore, it is supposed that the equation isbasically applicable for hill-slopes with steep slope
angle and wet condition:⎧k m d m i ( h / d m ) β
0 ≤ h ≤ dm⎪
q = ⎨k m d m i + k a i (h − d m )dm < h ≤ da⎪k d i + k i (h − d ) + α (h − d )m d < h
a
m
a
a
⎩ m m

(1)

where, q is the discharge per unit width; h is thestage; i is the slope; km and ka are saturated
hydraulic conductivities in the capillary pore andnon-capillary pore respectively.
(2) Application of OHDIS-kwmss model for the
Illinois River basin
To check the capability of the present storagedischarge equation to capture the hydrologicalprocesses in different type of basins, we haveapplied the model to the Illinois River basin in US(Fig. 1). It is one of the Distributed Model Intercomparison Project experimental basins8),9). The
Illinois basin at Tahlequah is about 2,400 km2 insize and situated in both Oklahoma and Arkansas
state, U.S. The average annual flow rate at the basin
outlet and the average annual rainfall of the basin
are about 29 m3/s and 1142 mm, respectively. In thebasin, rocky soils and outcrops tend to remain either
forest or pasture12). Hourly precipitation data
derived from gauged-adjusted radar and hourly
stream flow data at Tahlequah are obtained fromDMIP website9). Altogether 15 flood eventsoccurred between 1996 to 2001 were selected forthis study. The summery of all selected events are
given in Table 1. The 30 m digital elevation modelis also downloaded from the same website and it
was resampled to 300 m grid size. The OHDISkwmss model was then set up to the basin.

- 26 -

Table 1 Detail of the selected flood events at Tahlequah.

DMIP ID
Event-6
Event-7
Event-8
Event-9
Event-10
Event-11
Event-12
Event-13
Event-14
Event-15
Event-16
Event-17
Event-18
Event-19
Event-20

Rainfall /(mm)Rainfall /(mm)Fig. 1 The Illinois River basin.

Event 7

Rainfall /(mm)Rainfall /(mm)Discharge/(m3/s)

Discharge/(m3/s)

Event 6

Event 9Rainfall /(mm)Rainfall /(mm)Discharge/(m3/s)

Discharge/(m3/s)

Event 8

Event 13

Discharge/(m3/s)

Discharge/(m3/s)

Event 10

Fig. 2 Model calibration and validation at the Illinois Riverbasin.
120
100

Area/(%)

80

Maruyama
Illinois

60
40
20
0
0

10

20

30

40

50

Slope (degree)

Fig. 3 Slope-area comparison for the Maruyama and the Illinois
River basins.

Total RF
(mm)
58
120
92
110
87
75
135
67
105
37
15
95
108
39
56

Total Flow
(mm)
18.2
35.1
32.9
63.2
38.9
5.0
81.6
48.5
17.0
28.4
17.4
35.8
48.5
5.8
14.4

Initial Flow
(mm)
0.026
0.005
0.049
0.059
0.020
0.011
0.041
0.054
0.007
0.043
0.048
0.046
0.051
0.011
0.013

(3) Results and discussion for Illinois Riverbasin
The model was calibrated for the flood Event 6
in Table 1 using SCE-UA algorithm and thecalibrated model parameters were then use topredict the other 14 events. Prediction results ofsome events including the calibration event are
shown in Fig. 2.
According to our past experience in Japan,
calibrated parameters for any event can be used to
predict the other events with good prediction
accuracy7). In contrast, the model predictions are notgood for all events in Illinois basin (Fig.2).
Therefore, it is clear that the calibrated parameters
for Event 6 are not suitable to predict the otherevents. For further clarification on this, we havecalibrated the model using Event 9 data and used thecalibrated parameters to reproduce other events. We
observed similar results. It means that, the existingstorage discharge relationship is not well capture thegoverning hydrological process in the Illinois basin.
Furthermore, it is found that the objective function
value (Nash coefficient) for model calibration is0.68 and it is much lower value than that of usuallyexperienced in Japanese basins. Moreover, Event 9
and 13 shows little better prediction than others.
From Table 1, it can be noted that the initial flowrate of the Event 9 and 13 are higher than that ofother events. In other words, initial condition of thebasin for these two events is wetter than the otherevents. Therefore, it is concluded that thehydrological responses of this basin are highly
nonlinear for dry conditions.
Considering all these facts and figures, the model
is not applicable to reproduce the observedhydrographs even for calibrated events by changing
the model parameters. Thus we concluded that thepresent storage-discharge equation is not sufficientfor flood prediction of dry and mild-slope basinssuch as the Illinois basin. Figure 3 shows the slope-- 27 -

area distributions of the Illinois basin and the
Maruyama basin in Japan. According to this, theaverage slope of the Illinois basin is about 5 deg
while that of the Maruyama basin is about 20 deg.
In the derivation of the stage discharge equation
for the unsaturated zone, kinematic wave
approximations were used in the present model andthe soil moisture suction head term was neglectedassuming the steep slope. Illinois basin is a dry
mild-slope basin. Therefore the effect of soil
moisture suction head is important. We believe thatthe extension for the present model equations
considering the effect of soil moisture suction headis necessary for modeling highly nonlinear dry mildslope basins.

3. DERIVATION OF NEW STORAGEDISCHARGE
EQUATIONFORUNSATURATED ZONE
(1) Derivationofnewstorage-discharge
equation
The subsurface flow is conceptualized as two
types of flows, saturated flow through non-capillarypores and unsaturated flow through capillary pores.
First, unsaturated flow happens and after the water
storage in the soil becomes greater than a thresholddm, the saturated flow starts. Figure 4 shows thevertical section of the hill-slope, and Figure 5
shows the water storage h profile along the slope.
The total pore volume, capillary pore volume andwater storage in the soil at each cross section of theslope are represented as height da, dm and h.
For the unsaturated flow (0 ≤ h < dm) , byapplying the Darcy’s equation along the slope,average moisture flux velocity vm to the down slope
direction can be written as;v m = −k∂H∂( z + ϕ )⎛ ∂ϕ ⎞= −k= k⎜i −⎟∂x∂x
⎝ ∂x ⎠1andk ⎛α ⎞=⎜ ⎟k m ⎜⎝ φ ⎟⎠Fig. 4 Cross section of the soil layer.

xFig. 5 Soil water storage variation along the slope.

moisture suction curve. The parameter β has
relationship of β =3+2b.
According to the setting shown in Fig. 4,
α h
=
φ dm

Therefore;⎛ h ⎞⎟⎟⎝ dm ⎠

ϕ = H a ⎜⎜

−b

⎛ h ⎞
and k = k m ⎜⎜ ⎟⎟⎝ dm ⎠

β(5)Using Eq.3 and Eq.5 we can get;( − b −1)⎛∂h ⎞⎟1 ⎛ h ⎞⎜⎜ ⎟⎟vm = k ⎜ i + bH a⎜∂x ⎟dm ⎝ dm ⎠⎝⎠(6)

Therefore, the moisture flux rate per unit width, q
can be written as a function of storage and its
gradient;
β

⎛ h ⎞
⎛ h ⎞
q = k m d mi⎜⎜ ⎟⎟ + k mbH a ⎜⎜ ⎟⎟⎝ dm ⎠⎝ dm ⎠( β − b −1)

∂h
∂x

(7)

The difference between this new equation and theexisting equation for the unsaturated zone is that the
new equation has an additional term with the storage
gradient, the second term of Eq.7. Interestingly, it ispossible to incorporate the effect of soil moisturesuction by using only one additional parameter Ha
to the existing equation. Furthermore, by setting theparameter Ha to zero in the Eq.7, we can obtain
exactly the same equation as in the current OHDISkwmss model. Continuity equation of the water flowin the slope element can be written as:β∂h ∂ q+= r (t )
∂t
∂x

(4)

where, km, φ , α and Ha are saturated hydraulicconductivity, porosity, moisture content and air
entry suction pressure respectively (Ha < 0). The
parameter b is a constant and defines the soil

h

q

(3)

where k, i, H, z and ϕ , hydraulic conductivity, slope,
total head, elevation head and the moisture suction
head respectively.
In the case of steep, the difference of the soil
moisture suction head in the flow direction can beassumed very small as compared to the difference ofthe elevation head. Therefore the suction head term
can be neglected. However, for mild slope this could
significant and should be included. To incorporate
the soil moisture suction term, Campbell’s model(Eq.4) for soil moisture and hydraulic conductivity
is used:α ⎛ Ha ⎞ b⎟=⎜
φ ⎜⎝ ϕ ⎟⎠

Soil

(8)

where r(t) is the rainfall intensity to the slope
element.
Finally substituting the Eq.7 to Eq.8 we canobtain;- 28 -

2

∂h
∂h
∂ 2h
⎛ ∂h ⎞+ A + B⎜ ⎟ + C 2 = r (t )∂t∂x∂x
⎝ ∂x ⎠

where;

60 mm/hr(9)( β −1)⎛ h ⎞k bH (β − b − 1) ⎛ h ⎞⎜⎜ ⎟⎟,B = m a
A = k miβ ⎜⎜ ⎟⎟dm⎝ dm ⎠⎝ dm ⎠( β −b −1)⎛ h ⎞C = k mbH a ⎜⎜ ⎟⎟⎝ dm ⎠( β −b − 2 )Fig. 6 Test case 1 – rainfall applied to the entire hillslope.

(2) Numerical solution
This non-linear equation (Eq.9) needs to besolved for h to calculate q. In this regard, A, B and C
in the Eq.9 are calculated using known time level
while the derivatives are calculated implicitly17).
Following boundary conditions are assumed to solve
the above equation.
For x = 0; h = 0 for all t andFor t = 0; h = 0 for all x

60 mm/hrVmFig. 7 Test case 2 – rainfall applied at the upper half.
Table 2 Parameter values used in the models.

Using the five point implicit scheme, partial
derivatives can be written as;

Parametervaluehn −1, j − hn −1, j −1
hn , j − hn , j −1∂h= (1 − θ )
+θ
∂x
Δx
Δx

Km (cm/hr)
β (‑)

0.05
4

dm ( mm)600

hn, j −1 − hn −1, j −1
hn, j − hn −1, j∂h= (1 − λ )
+λ
∂t
Δt
Δt

da ( mm)

650

Ha− 2hn −1, j −1 + hn −1, j − 2h∂ 2h= (1 − θ ) n −1, j
+
∂x 2
Δx 2
hn, j − 2hn , j −1 + hn , j − 2

θ

Δx 2

Substituting above derivatives into Eq.9, h isrepresented as:
⎧
⎫
⎪
⎪
⎪
⎪
⎪r (e) − A⎡θ hn , j − hn , j −1 + (1 − θ ) hn −1, j − hn −1, j −1 ⎤ ⎪
⎢
⎥
⎪
⎪
Δx
⎢⎣
⎦⎥ ⎪
⎪
2
⎪
⎪
Δt ⎪ ⎛⎜ θ hn, j − hn , j −1 + (1 − θ ) hn −1, j − hn −1, j −1 ⎞⎟⎪hn , j =⎨− B ⎜
⎬
⎟
Δx
λ ⎪ ⎝
⎪(10)
⎠
⎪
⎪
⎤
⎪ ⎡θ hn , j − 2hn , j −1 + hn , j − 2 +
⎪
⎢
⎥
2
⎪
⎪
Δx
⎢
⎥
⎪− C
⎪
n
⎢
⎥
(1 − θ ) hn −1, j − 2hn −1, j −1 + hn −1, j − 2
⎪
⎪
⎥
⎪⎩ ⎢⎣
⎪⎭
Δx 2
⎦
(1 − λ ) 1−
hn , j −1 − hn −1, j −1 + hn −1, j

(

(

)

)

(

(

)

(

)

)

(

λ

(

)

)

where, n and j are current time step and current
spatial step, respectively. The parameters λ and θ areweights factors. We use 1 and 0.5 for λ and θ
respectively. Eq.10 is solved for hn,j iteratively andthen find the corresponding the value of q.

4. APPLICATION AND COMPARISION
OF THE NEW EQUATION WITH
EXISTING EQUATION
One simple hillslope is used to study thecharacteristics of the new equation and to compare itwith the existing equation. The same parametervalues were select for both models except the Ha

(mm)

-100

Slope Length (m)100parameter which is the additional parameter used in
the new equation. The parameter values are given in
Table 2. It is assumed that the hill-slope is initiallyfully dried. Constant rainfall intensity (60mm/hr for
10 hours) is applied for the hillslope and two cases
are tested; 1) rainfall is given to the entire hillslope,(this represent slope elements near to the river oredges with more water content in the downstream)(Fig. 6); and 2) to the upper half of the slope (thisrepresent when the rain occurred on hill mountains
while the downstream is dry)(Fig. 7) tounderstand the characteristics of flow for spatiallydistributed rainfall. Outflow hydrographs are
simulated for both equations for different slope
angles. With the parameters which we have chosen,
the soil remains unsaturated throughout thesimulation period.
Figures 8 and 9 show the simu latedhydrographs of the new and existing equations fordifferent hillslope angles for Case 1 and Case 2,
respectively. When apply the rainfall for entire slope(Case 1), it is found that difference of the simulatedhydrographs from new and existing equation is not
significant. It means if the soil moisture distribution
pattern looks like Fig.6, the differences ofdischarges by the two equations are small. In Case2, for higher slope angles, the difference between
two equations is also not significant. Even though
we consider the soil moisture suction in the new
equation, we cannot see difference because theelevation head difference for steep slopes is
considerably larger than the difference of soil
moisture suction. In contrast, we can see a clear- 29 -

it is found that the effect of lateral soil moisturesuction difference has significant effect on the
discharge.

Runoff height (mm/hr)

Slope 30￮
Slope 15￮

Slope 5￮Further research is to apply the distributedhydrologic model with the new storage-discharge
equation to the Illinois River basin and confirm theapplicability of the new equation for dry and mildslope basins.

Slope 3￮

REFERENCES

Runoff height (mm/hr)

Fig. 8 Test case 1 – Simulated hydrographs fordifferent slope angles.

Slope 30￮

Slope 15￮
Slope 5￮
Slope 3￮Fig. 9 Test case 2 – Simulated hydrographs fordifferent slope angles.

difference of the simulated hydrographs for mild
slopes in Case 2. For mild slopes, when the downslope is dry, the effect of soil moisture suction
becomessignificant
andthekinematics
approximation is inapplicable. In general, when thedown slope soil is dry, moisture flux velocity tends
to have higher value due to the effect of soil
moisture suction. The new equation captures thisphenomenon and starts the flow at the end of theslope earlier than that of the present equation.
The higher velocity of unsaturated flow in the
new equation leads to identify smaller hydraulicconductivity. This means if heavy rainfall isprovided, the soil is easily saturated and large flood
will happen. We believe the high nonlinear
characteristics of the new equation will provide
better simulation results at the dry mild-slopes
basins.

5. CONCLUSIONSzzz

Even though the existing storage discharge
equation is capable of giving good predictions
for wet and steep basins, flood predictions indry mild-sloped basins were not good.
New storage discharge relationship for
unsaturated zone considering the soil moisturesuction was developed.
Simulated flow by using both equations forhillslopes with different slope angles wereobtained and compared. For steep slopes, bothmodels give nearly the same results. However
for mild slopes with dry soil in the down slope,1) Troch, P. A., Loon, E.E.V and Hilberts, A.G.J: Hillslopestorage Boussnessq model for subsurface flow and variable
sources area along complex hillslopes:2Intercomparison
with a three-dimentional Rechards equation model, Water
Resours. Res., Vol.39, No.11, 2003.
2) Beven,K. : On subsurface stromflow: predictions withsimple kinematic theory for saturated and unsaturated
flows, Water Resours. Res., Vol.18, No.6, 1982.
3) Tachikawa, Y., Nagatani, G. and Takara, K: Development ofstage-discharge relationship equation incorporating
saturated–unsaturated flow mechanism, Ann. J. HydraulicEngng. JSCE 48, 7-12, 2004.
4) Troch, P. A., Loon, E.V and Hilberts, A: Analytical solution
to hillslope-storage kinematic wave equation for subsurface
flow, Advance in Water Resources, 25, 637-649, 2002.
5) Ichikawa, Y., Murakami, M., Tachikawa, T., and Shiiba, M:
Development of a basin runoff simulation system based ona new digital topographic model, J. Hydraulic, Coastal andEnviron. Engng. JSCE 691 (II-57)43-52, 2001.
6) Shiiba, M., Tachikawa, T., Ichikawa, Y., Hori, T, and Tanaka
K: Development of a slope runoff models which consider
field capacity and pipe flows, Annuals of Disas. Prev. Res.
Inst., Kyoto Univ., No. 41B-2, 1998.
7) Hunukumbura, P.B., Tachikawa, Y. and Takara, K:
Improvement of internal behavior in distributedhydrological model, Annual Journal of HydraulicEngineering, JSCE, Vol.52, 2008.
8) NOAA: Distributed Model Inter-comparison Project
http://www.nws.noaa.gov/oh/hrl/dmip/update_log.html,
January, 2008.
9) NOAA: Distributed Model Inter-comparison Project
http://www.nws.noaa.gov/oh/hrl/dmip/2/index.html,
February, 2008.
10) Reed, S, Koren, V, Smith, M, Zahang, Z, Moreda, F and
Seo, D : Overall distributed model intercomparison project,
J. Hydrol., Vol.298, 27-60, 2004.
11) Koren, V, Reed, S, Smith, M, Zahang, Z, and Seo, D :Hydrology laboratory research modeling system (HLRMS)of the US national weather service, J. Hydrol.,Vol.291, 297-318, 2004.
12) Vieux, E.B, Cui, Z and Gaur, A: Evaluation of a distributedhydrological model for flood forecasting, J. Hydrol.,Vol.298, 155-177, 2004.
13) Duan, Q., Sorooshian, S. and Gupta, V.K: Effective andefficient global optimization for conceptual rainfall-runoff
models, Water Resours. Res., Vol.28, No.4, 1015-1031,
1992.
14) Duan, Q., Sorooshian, S. and Gupta, V.K: Optimal use of
the SCE-UA global optimization method for calibratingwatershed models, J. Hydrol., Vol.158, 265-284, 1994.
15) Maidment, D.R: Handbook of Hydrology, McGRAW-HILL,
INC, 5.6-5.7, 1993.
16) Shiiba, M., Ichikawa, Y., Sakakibara, T. and Tachikawa, Y:
A new numerical representation form of basin topography,
J. Hydraulic, Coastal and Environ. Engng. JSCE 621 (II-47)
1-9,1999.
17) Singh, V.P: Kinematic wave modeling in water resources,
John Wiley & Sons, 691-965, 1996.
(Received September 30, 2008)- 30 -