京 都 大 学 防 災 研 究 所 年 報 第47号 C 平成16年3月
Annuals of Disas. Prev. Res. Inst., Kyoto Univ., No. 47 C, 2004

Distributed Runoff Model Linking Surface with Groundwater Processes

Amin NAWAHDA*, Toshiharu KOJIRI, and Shuichi IKEBUCHI* Graduate School of Engineering, Kyoto UniversitySynopsis
A method with multi-layer and mesh-typed runoff model using Hydro-BEAM(Hydrological River Basin Environment Assessment Model) is proposed to analyze theintegrated hydrological processes. The spatiotemporal simulation is calculated with thekinematic wave model for surface runoff, Richard’s equation for unsaturated subsurface
flow and the unconfined flow for groundwater. The initial loss of rainfall due tointerception by depression storage reprocess is considered here. Moreover the basin
division and land use dynamics are introduced to encounter reservoir operation andland utilization with human activities. The proposed model is calibrated for different
initial conditions and parameters, and applied into the Yasu River to verify the dynamic
linkage between surface and groundwater.
Keywords: Distributed runoff model, Parameterization, Saturated and unsaturated flow,
Spatio-temporal distribution.

1. Introduction

equations are formulated as follow;

The hydrological cycle involves complicatedinteractions between atmosphere layer, surface andground layers. The integrated hydrological modelingis achieved by dynamic linking and simultaneous
calculation of the water cycle related processes.
Traditionally the lumped models simulate theformation of the surface runoff. These models arebased on the assumptions of uniformity and linearityof watersheds. In reality these assumptions encountera high percentage of errors. The distributed rainfallrunoff (Hydro-BEAM) was introduced in order toovercome the shortcomings of lumped models Kojiri(2000). At each mesh the heat balance method andthen the kinematic wave model are used for
calculating the surface runoff. The linear storage
model is used for groundwater flow modelling. The
ground of the whole basin is divided into four
uniform layers. The mass balance and the momentum∂A ∂Q+= q ( t, x )∂t ∂x
A = f ( x,Q ) = αQβ

(1)
(2)

where A is the discharge area, Q is the flow rate, q is
the spatiotemporal effective rainfall, and α and β arethe parameters of Manning equation. The generalmomentum equation from the first to forth layer isformulated with Darcy’s law. In this study thewatershed is treated extensively with more details.
These details encounter distributed criteria for eachmesh such as; geological cross sections, soil maps,
dynamics of land uses, sinks and sources, and otherhydroclimatic processes. The flow chart of thecalculation procedure for the integrated HydroBEAM is shown in Fig. 1.

parameter repre-sent the space of variability, and an
optimized value must lye within the boundaries ofthis space.

Large time step (hr)Total mesh number
Small time step (min)Water source
Interception

Soil 1

Evapotranspiration

Utilization

Soil 2

City

Mountain

Water body

Small time step (sec)Soil moisture

Agriculture
Soil nInfiltration

Groundwater recharge
River discharge

RunoffGroundwater level

Fig. 1 The main flowchart of integrated HydroBEAM

2. Integrated hydrological modeling
2.1 Spatiotemporal variability of hydrologicalparameters
The surface and subsurface hydrological parameters vary in space and/or time; this variability
causes a significant change in the shape of the riverhydrograph. Therefore the watershed is divided intocatchments, each catchement is girded into squaregrids, and each grid is connected to a drainage network represented by a channel segment. The variability of the hydrological parameters for each grid is
integrated, and then the grid is divided into elementsof uniform parameters. The number of elements isdecided by the surface properties and the availabilityof hydrological data. In kinematic wave model thesoil roughness affects both the speed of runoff andinfiltration rate, infiltration affects both the depth of
the surface runoff and the groundwater recharge.
Therefore for each type of soil and land uses there
are a preliminary set of roughness and infiltration
coefficients. These coefficients are evaluated byusing the scaling methods and the scaling factors arederived by two steps; first is soil sampling fromdifferent locations in the watershed and second is to
use the probability densities and auto correlation
structures for the calibration. The rational method for
calculating the representative infiltration and roughness coefficients for each mesh is shown in Fig.2. It
shows the soil map based parameterization n*, andthe land use based parameterization n, where n is asurface hydrological parameter, Lj is the area
percentage of the sub parameter nj, and A is the areaof the ith mesh. The two values for a hydrologicaln*i =∑n L*jPaddy
jAi,

ni =

∑n L
jj

AiFig. 2 Integrated hydrological parameters based on
distributed soils and land uses

2.2 Subsurface flow process
Richards equation Eq. (3) for unsaturated flow in
one dimension Richards (1931) is used to simulate
the distribution of the soil moisture and to estimate
the spatiotemporal groundwater recharge.

C∂Ψ ∂Ψ  ∂Ψ  =+ 1 K (Ψ) ∂t∂z  ∂z

(3)

where, C is the water capacity, ψ is the water
potential head, K is the hydraulic conductivity, z is
the thickness of the unsaturated layer, and t is time.
2.3 Groundwater flow process
Quasi three-dimensional equation for thetransient unconfined groundwater flow with Dupuitassump-tions as illustrated in Eq. (4) is used;1∂2h2∂2h2 ∂hK+K= S − R ( x, y, t ) xy22 2∂x∂y ∂t(4)where h is the groundwater level, K is the hydraulicconductivity, S is the specific yield, R is thedistributed groundwater recharge, x and y are thedisplacement coordinates, and t is time.
2.4 Interception process
A fraction of rainfall is not contributing with therunoff formation, caused by interception. Similar tothe interception by vegetation as it was derived byHorton (1919), the interception caused by residentialbuildings can be derived as follows;
I = C + VEt

(5)

.

where I is the interception loss, C is the interception
storage depth, V is the ratio of building’s surfaces
area to its projected area on the ground, E is theevaporation rate, and t is the storm duration.

Therefore Hydro-BEAM allows for automatic
adjustment of the distribution of land uses for eachmesh. On the other hand the hydrological parameters
are also subjected to the same dynamic behavior.

2.5 Evapotranspiration process (ET)Evapotranspiration is a description for both evaporation and transpiration by plants, the temporal
behaviours of weather conditions and thegroundwater level are affecting the rate of ET. Thereare many methods to estimate ET according to theamount of available metrological data. The heat andmass balance methods are used to evaluate ET foreach mesh in the watershed as a function of the net
solar radiation Kojiri (2000). The snowfall andsnowmelt processes are simulated by using theenergy balance at the earth’s surface. An empiricalformula for the threshold temperature Tc is evaluatedfor each mesh in order to predict the formation ofsnow Park et al (2003);

3. Case study: The Yasu River Basin

Tc = 11.01 − 1.5e a(6)where ea is the vapor pressure. If T ≤ Tc; thensnowfall. The snowmelt is simulated with an updated
energy balance equation that includes the snowmelt
heat, temperature near the snow surface, snow coverstored heat, rainfall heat, and the latent heat for thesnow layer.
2.6 Reservoir operation and change of basinsituation
The basin division approach can be defined as
dividing the watershed that has a reservoir into upperand lower parts. Compared with the lower catchment
the upper catchments of the watershed is subjected todry and wet precipitations that differ in time andamount, also the stratigraphy of the lower part differsfrom the upper part. The wave propagation from theupper catchments down to the lower catchments will
be subjected to the well known shock wave, becauseof the flow of water from high levels to low levels in
the lower catchments Shin (1979). On the other hand,
the reservoir operation can be considered as a strong
down stream control for the developed kinematicwave from the upper catchments, so the kinematicwave approximation will be no longer valid if thereservoir operation is not included Ishihara (1959).
The seasonality of the reservoir operation, and theextremes of the seasonal water demands and safety
operation are set in the hydrological model, and thenlinear storage model is used for routing the reservoir.
The changes of land uses within the watershed
will occur dynamically in space and/or in time.

The Yasu River is located in Shiga prefecture,
Japan. The area of the whole basin is 445 km2 . The
hydrological data for the Yasu River basin wereprepared. Spatial data includes, geological formation,soil map, elevation, land uses, slope, flow direction,and existing sewers. Temporal data includes; meteorological data, land utilization, groundwater level,
operation of existing sinks, water consumption, andseasonal measurements of the hydrological parameters. There are two reservoirs; the Yasu and theAoto reservoirs. The reservoirs are used to supply
domestic water and to satisfy industrial water
demand. The geological cross sections from the Yasu
River basin are analysed; there are four main types ofrocks: granite; chert; sandstone; and limestone. The
soil and the geological maps are used in order to
distribute the soil types and geological formations foreach mesh. There are five main soil types; gray low
land, gray, peat, not matured soil, and dry brown
forest soils. According to the boring logs andgeological cross sections of the study area, there aremainly unconfined and confined layers respectively,
the unconfined layer has an average thickness of 11m. The hydraulic parameters of the unconfined layer
are determined by the rational method for the surface
sub-layers. The groundwater divide in south, north,and east of the Yasu Basin forms the basinboundaries. These boundaries are assumed to be noflow boundaries, the western boundary will be used
for calibrating the groundwater model. Threedimensional simulations for the unconfinedgroundwater levels were carried out. A one-year
simulation was done by using these data and by
applying physical parameters obtained from field.
These data was set as initial conditions for thegroundwater flow model.
4. Simulation results and discussion
Integrated Hydro-BEAM is applied to hourly datafrom the upper and lower catchments. Initialsaturated depth profiles were evaluated according tothe distri-buted groundwater levels. Precipitation andthe corresponding runoff in the Yasu River basin
outlet are shown in Fig.3. The groundwater level foreach mesh is simulated by the dynamic linking
between the surface water models and groundwatermodel in the integrated Hydro-BEAM. The
distributed ground-water levels for the lower
catchments are shown in Fig.4.
The soil moisture profile for each mesh issimulated by considering the spatiotemporalparameterization which is caused by the land use
dynamics. At the start of the rain seasons the soil isassumed to be dry, and according to the soil types
and depths of groundwater the potential head at theboundary of the unsaturated layer is evaluated. AtYasu River basin the groundwater potential headvalue ranges from -40 to -100 m, the soil layer isdivided into (100) divisions between the surface andthe groundwater level. Then the soil moisture and theRainfall

distribution of the potential head is simulated every
10 sec, during this small time step the groundwaterlevel is assumed to be steady and is not affected by
the soil moisture. At every large time step thegroundwater level is updated in order to set the lower
boundaries for the unsaturated flow model. Then thecalculated soil moistures will be accumulated forevery large time step and used for estimating the
distributed recharge to groundwater model. Also themodel simulates the distributed discharge fromgroundwater storage, because the accumulated soil
moisture can be negative in the case of dry weather
conditions and extensive water uptake from thesubsurface layers.

Observed

Simulated
05080

10060

150

40

Rainfall (mm/day)

1002002002501316191121 151 181 211 241 271 301 331 361
Time (day)Fig. 3 Simulated and observed daily river discharges at Minakuchi station (1997)

60

1.0
2.0

40

3.0
20

4.00

5.0
11001

2001

3001

4001

5001

6001

7001

Time (hr)Fig. 4

Simulated and observed (thin line) ground-water levels in lower Yasu basin

Groundwater level (m)0.0Rainfall (mm/hr) .

Discharge (m3/s)

120

5. Conclusions
A general framework for the integrated hydrological modeling has been presented. The framework
includes three main linked models: distributed
rainfall-runoff model, unsaturated flow model, andunconfined groundwater model. The importance oflinking the three models has been illustrated andemphasized. Manual calibration has been conducted
for the case study of the Yasu River basin. The
simulated groundwater levels and the riverdischarges show a good agreement with the observed
values, but in certain catchments there was a poor
agreement with the corresponding observed data.
The manual calibration procedure for the model
parameters doesn’t succeed in finding global optimalvalues; the spatiotemporal fluctuations of the basincharacteristics might harden the attempts for finding
global optimal parameters for every catchment within
the basin. The integrated hydrological modelingincluding groundwater and surface water might

require long time periods for automatic calibrationand more research is needed in this subject.
References
Ishihara T. (1959): Fundamental researches on theunit hydrograph method and its application, Japan
Society of Civil Engineers, Vol. 60, No. 3-3, pp. 19.
Horton R. E. (1919): Rainfall Interception, Monthly
Weather Review, Vol. 47, No. 1, pp. 603-623.
Kojiri T. (2000): GIS-based environment model forwater quantity and quality with river basin
simulation, Annual DPRI, Kyoto University, Vol. 44,No. 1, pp. 150-160.
Park J., Kojiri T., Tomosugi K. (2003): Developmentof GIS based distributed runoff model forenvironmental assessment, Japan Soc. Hydrol. &Water Resources, Vol. 16, No. 1, pp. 541-555.
Richards L. A. (1931): Capillary conduction ofliquids through porous media. Physics 1, pp. 318333.
Shin H. (1979): Modeling of rivers, 1st edition, Wiley.

要 旨
統 合 的 な 水 文 過 程 を 解 析 す る た め に 、 Hydro‑BEAM(Hydrological River Basin Environment
Assessment Model)をベースとした多層メッシュ型流出モデルを提案する。表面流出にはキネマティ
ックウェーブモデルを、不飽和地中流にはリチャーズ式を、地下水流には非圧地下水モデルを用いて時空間的な解析がなされる。ここでは窪地貯留や樹冠遮断に伴う降雨の初期損失が考慮されている。さらに、貯水池操作や人間活動に伴う土地利用変化に対応するために、流域分割や土地利用のダイナミクスが導入されている。ここで提案するモデルは異なる初期条件やパラメータを用いてキャリブレーションされ、地表水と地下水の直接的な結合（dynamic linkage）を検証するために、野洲川流域に適用された。
キーワード: 分布型流出モデル，パラメタリゼーション，飽和不飽和地中流，時空間分布