京 都 大 学 防 災 研 究 所 年 報 第 50 号 B 平 成 19 年 4 月
Annuals of Disas. Prev. Res. Inst., Kyoto Univ., No. 50 B, 2007

Quantification of Parameter Uncertainty in Distributed Rainfall-Runoff Modeling

Giha LEE*, Yasuto TACHIKAWA*, and Kaoru TAKARA*Graduate school of Urban and Environment Engineering, Kyoto University
Synopsis
In general, hydrological models have several (or a lot of) parameters that cannot bedirectly measured, which only are inferred by calibration procedure against a historicalinput-output data record. While the applications of automatic parameter estimation
techniques have received considerable attention over the last decades, such classical methods
have received criticism for their lack of rigor in handling with uncertainty in the parameter
estimates. This work addresses the calibration of the distributed rainfall-runoff model
KsEdgeFC2D, the quantification of parameter uncertainty and its effect on the prediction of
streamflow for Kamishiiba catchment (211km2). In this study, to analyze the propagation ofparameter uncertainty into prediction, we employ the Shuffled Complex Evolution
Metropolis (SCEM-UA) global optimization algorithm. Moreover, we compare SCEM-UA
derived optimal parameter values to those estimated using deterministic SCE-UA methodwith three different objective functions to account for the structural stability of KsEdgFC2D
model and to demonstrate the capability of the SCEM-UA algorithm to efficiently evolve toparameter posterior distribution.
Keywords: Automatic parameter estimation, Parameter uncertainty, SCEM-UA
1. Introduction
Rainfall-runoff process in heterogeneous real worldis commonly simplified and represented by varioushydrological models (Wagener et al., 2004). Thesemodels are conversion and simplification of reality,
thus no matter how spatially sophisticated and accurate
they may be those models only represent aspects ofconceptualization or empiricism of modelers (orhydrologists). Accordingly, their outputs are as reliable
as hypothesis, structure of models, and quantity andquality of input data, and parameter estimates (Gupta
et al., 1999; Muletha and Nicklow, 2005). Modelparameters, in general, are classified into two (Kuczera
and Franks, 2002). First is physical parameter which is
physically measurable property of the watershed andthe other is conceptual parameter that can only beinferred by some matching process, often called
calibration, between the simulated responses and theobserved ones (Sorooshian and Gupta, 1995). Fromabove
definition,
virtuallyevenso-called‘physically-based’ models, be they lumped or
distributed, would be regarded as being conceptual.
The conceptual parameter values are adjusted
between each run of the model, either manually bymodelers or automatically by some computer-based
optimization algorithm and the corresponding
observation until some optimal parameter set has beenfound. Especially, the application of automaticparameter optimization algorithms has received
considerable attention and has improved over the last
decades (e.g., Sorooshian and Dracup, 1980; Duan etal., 1992; Vrugt et al., 2003). Until the early 1990’s,
the available automated optimization techniques
practically were not able to find the global optimalvalues in a prescribed objective function because of thepresence of multiple local optima, parameterinteraction, and discontinuous and non-convex－ 45 －

response surface defined by the selected objectivefunction. These insights in the response surface led tothe development of global optimization algorithm,
namely Shuffled Complex Evolution (SCE-UA)
optimization algorithm (Duan et al., 1992; 1993; 1994).
However, while remarkable progress has been made in
the development and application of automatedcalibration procedures, such methods have not been
free from criticism for their insufficiency handlingwith parameter uncertainty (Vrugt et al., 2003; 2005).
As many hydrologists pointed out in theirliteratures, estimates of hydrologic model parameters
are generally error-prone because used data (e.g.,
rainfall, streamflow) during calibration containmeasurement errors and model structural imperfection
never explicitly represents the system or fits theobserved data (Schaake, 2003). Manifold studies havebeen conducted to quantify or assess parameter
uncertainty and its propagation into subsequentprediction results. For example, multinormal
approximations (Kuczera and Mroczkowski, 1998),
simple uniform random sampling over the feasible
parameter space (Uhlenbrook et al., 1999), andMarkov Cahin Monte Carlo (MCMC) methods(Kuczera and Parent, 1998; Vrugt et al., 2003) aredeveloped to analyze parameter uncertainty. In
comparison with traditional statistical theory based on
first-orderapproximations
andmultinormal
distributions, MCMC methods have become
increasingly popular as one of the general purpose
approximation methods for complex inference, searchand optimization problem. Shuffled Complex
Evolution Metropolis (SCEM-UA) algorithm is aneffective and efficient evolutionary MCMC sampler
which has enhanced search capability. This method,
modification of the SCE-UA algorithm by developedby Duan et al., (1992), operates by merging thestrengths of the Metropolis algorithm, controlled
random search, competitive evolution, and complex
shuffling. The stochastic characteristic of Metropolis
scheme makes it possible to avoid the tendency of the
SCE-UA algorithm to collapse to the global minimum
on response surface. This powerful algorithm canprovide not only optimal parameter set but also its
underlying posterior distribution within a single
optimization run (Vrugt et al., 2003).
In this study, we employ the SCEM-UA method toestimate confidence limits for model parameters.

Especially, the aim of this work is to explore thecapacity of SCEM-UA to identify the posterior
parameter distribution for the complex KsEdgeFC2D
model applied to the Kamishiiba catchment (211 ㎢),
and to evaluate the effect of parameter uncertainty onhydrograph simulation. Moreover, we compare
SCEM-UA derived parameter values to thoseestimated using deterministic SCE-UA method withthree different objective functions, Simple Least
Square estimation criterion (SLS), Heteroscedastic
Maximum Likelihood Estimator (HMLE) and
Modified Index of Agreement (MIA) to account for the
structural stability of KsEdgFC2D model and todemonstrate the capability of the SCEM-UA algorithm
to efficiently evolve to parameter posterior distribution.
This paper is organized as follows. In section 2, we
describe
thedistributedhydrologic
model,
KsEdgeFC2D. Section 3 presents details on the
SCEM-UA algorithm used to infer the posterior model
parameter distributions as compared with SCE-UA. In
section 4 and 5 the result of case study is illustrated.
Finally, we summarize methodology and analyzed
results in section 6.
2. Applied Hydrologic Model
KsEdgeFC2D is a conceptual distributedhydrologic model developed by Ichikawa et al. (2001)
includingdischarge-stage
relationshipwithsaturated-unsaturated flow (Tachikawa et al., 2004).
The model solves the one-dimensional kinematic wave
equation with the discharge-stage equation using theLax-Wendroff finite difference scheme according to
the flow direction map (see Figure 1). All
geomorphologic information is extracted from 250m
based DEM.
The model assumes that a permeable soil layer
covers the hillslope as illustrated in Figure 1(d). The
soil layer consists of a capillary layer in which
unsaturated flow occurs and a non-capillary layer
where saturated flow occurs. According to thismechanism, if the depth of water is higher than the soil
depth, then overland flow occurs. The discharge-stage
relationship is expressed by equation (2) corresponding
to water levels (see Figure 1(d)) defined as:－ 46 －∂h ∂q+= r (t )
∂t ∂x

(1)

(a)

(b)

Radar Rainfall
or Gage Rainfall
Grid Precipitation
Geo-data
Processing
Rainfall-Runoff
Simulation(c)surface flow
subsurface
flow
Hillsloperouting
Channelrouting(d)(Option)Discharge at
Nodes and Edges

Fig. 1 Schematic representation of KsEdgeFC2D (a)Modular structure of KsEdgeFC2D (b) Distributed
grid rainfall data (c) Slope and channel components
extracted from DEM (d) Model structure for thehillslope soil layer and discharge -stage relationship⎧vc d c ( h / d c ) β ,0 ≤ h ≤ dc⎪
q = ⎨vc d c + va ( h − d c ),dc ≤ h ≤ d s⎪v d + v ( h − d ) + α ( h − d ) m ,ds ≤ h
c
c
a
c
s
⎩

(2)

Flow rate of each slope segment are calculated byabove governing equations combined with thecontinuity equation like equation (1). where, v = k i ;v = k i ; k = k / β ; α = i / n ; i is slope gradient, kis saturated hydraulic conductivity of the capillary soil
layer, kis hydraulic conductivity of the
non-capillary soil layer, n is roughness coefficient, dis the depth of the capillary soil layer and d is soil
depth. Detailed explanations of model structure appear
in Tachikawa et al., (2004).
c

a

a

c

a

c

c

a

c

s

3. Methodology
Estimation of the reliability of model output isimportant to hydrologic engineering and water
resources planning since model output reliability canbe used to assess model verification and the selectionof suitable models (Melching et al., 1990). Beven(1989) has reviewed the limitations of the current
generation of distributed, physically-based models ofwatershed hydrology and has suggested that a possible
way forward requires a realistic estimation ofprediction uncertainty. Only recently methods forrealistic assessment of parameter uncertainty havebegan to appear in the literature. These include, as
stated in the introduction part, the classical use of

first-order approximations to parameter uncertainty,
evaluation of likelihood ratios, and parametricbootsrapping or MCMC methods. SCEM-UA is one ofthe popular and efficient tools to evaluate parameter
uncertainty and predictability of hydrological models(Feyen et al., 2006).
In this study, we examine the applicability orcapability of SCEM-UA to assess hydrologic model
parameter uncertainty and to search for best-fit
parameter values through comparing with those
optimized by deterministic automated global searchalgorithm, SCE-UA. Moreover, the influence ofobjective functions (e.g., SLS, HMLE, MIA) on
optimal parameter and model performance isdemonstrated to account for the model structural
uncertainty (Lee et al., 2007). In the following
subsections 3.1, 3.2, and 3.3 we describe the SCE-UA
based on three different objective functions to searchfor global optimal parameter set and the SCEM-UA
based on MCMC sampler using Metropolis Hasting(MH) strategy respectively.
3.1 Shuffled Complex Evolution (SCE-UA)
Algorithm
The SCE-UA is used to identify the best fitted
parameter set, which is a single-objective global
optimization method designed to handle with high-parameter dimensionality encountered in calibrationof a nonlinear hydrologic simulation models. Thisevolutionary approach method has been performed bya number of researchers on a variety of models withoutstanding positive results and has proved to be an
efficient, powerful method for the automatic
optimization. SCE algorithm is basically synthesized
by following three concepts: (1) combination of asimplex procedure with the concepts of controlled
random search approaches; (2) competitive evolution;and (3) complex shuffling. The integration of thesesteps above mentioned makes the SCE method
effective and flexible.
3.2 Applied Objective Functions for SCE_UA
Calibration
The aim of computer-based automatic calibration is
to find the values of the model parameters thatminimize or maximize the numerical value of theobjective functions. In general, the most commonly
utilized objective functions in hydrological modeling－ 47 －

are variations of the Simple Least Squares (SLS)
function defined as:
N

SLS = ∑ (qtobs − qt (θ )) 2

(3)

t =1

where qtobs is observed stream flow value at time t;
qt (θ ) is model simulated stream flow value at time tusing parameter set θ ; N is the number of flow values
available. SLS has a feature that large discharge isemphasized due to squared errors while low flows are
neglected, thus the parameter set fitting around peakdischarge value is likely to obtain.
A newly proposed measure is the Modified Index
of Agreement (MIA) to reduce the influence of thesquared term during high flows as putting a weight onflow values. This objective function is calculated as:⎛ qtobs − qt (θ ) ⎞⎟qtobs
t =1 ⎝⎠2n∑⎜
MIA = 1 −⎛ qt (θ ) − qtmean + qtobs − qtmean ⎞⎜⎟∑⎟qtmeant =1 ⎜
⎝
⎠

(4)
2

n

where qtmean is mean value of observed time series.
Sorooshian and Dracup (1980) proposed a different
objective function to consider entire behavior ofhydrograph, the Heteroscedastic Maximum Likelihood
Estimator (HMLE), which enables to estimate the most
likely weights through the use of the maximum
estimation theory. This new measure can eliminate
some of the subjectivity involved in the selection oftransformation and / or a weighting scheme byhandling heteroscedastic error, so that it yields a morebalanced performance over the entire hydrograph. It is
calculated as:
wε
1 ∑
N

t

min HMLE =
θ ,λ

N

t =1
N

∏w

t

(5)

t

t =1

where ε t = qtobs − qt (θ ) is the model residual attime t; wt is the weight assigned to time t, computed
as wt = f t 2( λ −1) ; f t = qttrue is the expected true flow attime t;λ is the transformation parameter which

stabilizes the variance.
3.3 Shuffled Complex Evolution Metropolis(SCEM-UA) Algorithm
The Markov Chain Monte Carlo sampler, entitledin Shuffled Complex Evolution Metropolis algorithm
is well suited for the practical assessment of parameter
uncertainty in hydrological models. This sampler
incorporate effective characteristics of SCE-UA suchas controlled random search, competitive evolution,and complex shuffling with the strengths of theMetropolis-Hasting algorithm to evolve a populationof sampled points to an estimation of the stationary
posterior distribution of the parameters. Two big
revisions prevent the convergence to a small attraction
region, instead, and facilitate convergence to astationary posterior target distribution of parameters.
The first modification is replacement of the downhill
simplex method by a Metropolis annealing
covariance-based offspring approach, thereby avoiding
a deterministic drift toward a single mode. Second, the
SCEM-UA does not divide the complex into
subcomplexes during the generation of the offspring
and uses a different replacement way, to terminate
occupations having lower posterior density on theparameter space. General steps for SCEM-UA
implementation are outlined below.
3.3.1 SCEM-UA process
1. Initialize the process including the selection of the
population size s and the number of complexes q.

2. Generate s samples {θ 1 , θ 2 , L , θ s } from thepredefined prior distribution and compute theposterior density p (θ1 y ) , p (θ 2 y ) ,L , p (θ s y ){}of each point, where y = known observed data.
3. Sort the points in order of decreasing posterior
density and store them in array D[1: s,1: n + 1] ,where n is the number of parameters.
4. Initialize the starting points of the parallel
sequences S 1 , S 2 ,L , S q , such that S k is
D[k ,1: n + 1] , where k = 1, 2,L , q .
5. Partition the s points of D into q
complexes C 1 , C 2 ,L , C q , each containing m points.
6. Evolve each sequence. Evolve each of the parallel
sequences according to the Sequence Evolution
Metropolis (SEM) algorithm described in 3.3.2.

－ 48 －

7. Unpack all complexes C back into D and rank thepoints in order of decreasing posterior density.
8. Check the convergence criteria, Gelman and Rubin(GR) statistic. If GR are satisfied, stop; otherwise,
return to step 5.
3.3.2 SEM algorithm, Key component of the
SCEM-UA
This algorithm produces new candidate points ineach of the parallel sequences, S k by generating draws
from an adaptive proposal distribution by using theinformation induced in the m samples of C k .
1. Compute the mean μ k and covariance structureΣ k of the parameters of C k . Sort the m point incomplex C k in order of decreasing posterior
density.
2. Compute α k , the ratio of the mean posterior
density of the m points in C k to the mean posterior
density of the last m generated points in S k .
3. If α k is smaller than a predefined likelihood ratio, T,
generate a candidate point, θ ( t +1) , by using a
multinormal distribution centered on the last draw,θ ( t ) , of the sequence S k , and covariance structure
cn2 Σ k , where, cn is a jumprate. Otherwise, go tostep 4.
4. Generate offspring θ ( t +1) , by using multinormal
distribution with mean μ k and covariance structure
cn2 Σ k .
( t +1)y ) , of θ ( t +1) .
5. Compute the posterior density, p (θIf θ ( t +1) is outside the feasible space,( t +1)y ) to zero.
set p (θ( t +1)y ) / p (θ ( t ) y ) and6. Compute the ratio r = p (θrandomly sample a uniform label Z over interval [0,1].
7. If Z is smaller than or equal to r, accept the new
candidate point. Otherwise, reject the candidate
point and then remain at the current position in thesequence, it means θ ( t +1) = θ ( t ) .
8. Add the point θ ( t +1) to the sequence S k .
9. If the candidate point is proper, replace the best
member of C k with θ ( t +1) , and go to 10; otherwise
replace the worst member of C k with θ ( t +1) .
10. Compute AR k , the ratio of the posterior density ofthe best to the posterior density of the worst
member of C k .
11. Repeat the steps 1-10 until predefined number of

iteration before complexes are shuffled.
A detailed description and explanation of the
SCEM-UA method appears in Vrugt et al. (2003).
4. Case Study
The study site is the Kamishiiba catchment,
upstream area of Kamishiiba dam, which lies within
Kyushu region in Japan and covers an area of 211km2.
The topography of this area is hilly with the elevation
varying from 400m to 1700m. Land use type is forest,
thereby enabling to regard the drainage as a typicalJapanese mountainous area. In the model, thecatchment is based on 250 by 250m grid blocks.
Observed discharge data converted from water level of
dam inflow having 10min temporal resolution areavailable for Kamishiiba dam gauging station (Kyushu
Electric Power Co., Inc). Distributed rainfall (1km by
1km, see figure 1(b)) is gauged on Eshiroyama radar
station. Table 1 shows historical short-term events forthis study.
Table 1 Historical storm events
Peakdischarge(m3/s)

objective of
Event

Date ofoccurrence

Rainfallduration(hr)

Calibration

1999 / 09 / 22

96

590

Validation

1997 / 09 / 15

84

1203

Although recent hydrologic models are based on
physics to a certain extent, some processes are only
represented in a lumped conceptual way (Feyen et al.,2006). Consequently, even some parameters ofphysically based model lack physical basis and cannot
be directly observed from field surveying. In theKsEdgeFC2D, there remain five parameters that need
to be estimated by calibration procedure against
observed streamflow data. Parameters to be optimized
and a range of prior distribution are illustrated in Table
2.
The SCEM-UA algorithm contains two algorithmic
parameters as well as the SCE-UA method that need to
be specified manually by the user: 1) the number ofcomplexes and sequences, q; 2) the population size, s,
which also determine the number of points within eachcomplex. Vrugt et al. (2003) recommended the use oflarger population sizes and larger number of parallel
sequences to be able to precisely capture the complex－ 49 －

shape of the covariance structure. Moreover, the SEM,
sub-process of SCEM-UA algorithm also has threealgorithmic parameters to be chosen carefully: 1) the
number of evolution steps defined as L=m/10; 2)predefined likelihood ratio T; 3) predefined jumprate,
cn = 2.4 / n . Essential algorithmic parameters forrunning SCEM-UA are tabulated in Table 3.
Table 2 Calibration parameters of the KsEdgeFC2D
model with upper and lower bounds of the predefinedparameter space
Parameters optimized

Lower

Upper

(KsEdgeFC2D)

bound

bound

n

0.1

0.5

ka

0.01

0.05

d s (m)0

1.0

d c (m)0.0

0.6

β = k / k (beta )

0

20

a

c

functions to identify the structural inadequacy of usedhydrologic model.
The SCEM-UA algorithm was implemented with apopulation size s=200 and q=10 complexes, 20 points
in each complex. An important issue in MCMC
sampling is convergence of the sampler to thestationary posterior distribution. Practically, Scale
Reduction score ( SR ) developed by Gelman and
Rubin (1992) has been used to be a criterion of
convergence. If the SR is less than 1.2, the Markov
chain is considered to be converged into the target
posterior distribution; otherwise, the evaluation steps
are repeated until obtaining a suitable value. Figure 2
illustrates the calculated SR against the number of
MCMC iterations. All parameters become stable after
approximately 5,000 iterations (i.e., SR <1.2).

Table 3 Algorithmic parameter in SCEM-UA method
Parameterq

Description
the number of complexes
q = 10s

the population size
s = 200

Lthe number of evolution steps
L = m /10 = ( s / q ) /10 = (200 /10) /10 = 2

TFig. 2 Evaluation of the Gelman and Rubin Scale
Reduction score for KsEdgeFC2D parameters

likelihood ratio
T=106

cn

jump rate
cn = 2.4 / n = 2.4 / 5 ≈ 1.07

5. Results
This case study illustrates the efficiency andeffectiveness of the Shuffled Complex Evolution
Metropolis algorithm for realistic assessment ofprediction uncertainty on hydrologic responses
provided by the complex hydrological modeling like
distributed modeling. In addition to analysis ofparameter uncertainty, we compare the optimal
parameter values estimated using SCEM-UA with theones derived using the original SCE-UA global
optimization algorithm with different objective

Figure 3 presents the marginal posterior probability
distributions for the KsEdgeFC2D parameters obtainedusing 5,000 sample parameter sets generated afterconvergence of the SCEM-UA algorithm. The
moments of the posterior parameter distributions andthe most likely parameter combination are presented inTable 4. Here, the first 5,000 simulations of eachparallel sequence were discarded (i.e., SR >1.2). The
limit of the X-axis in Figure 3 corresponds to the range
specified for the estimated parameter uncertainty
through the stochastic optimization procedure. The
histograms for the k , d and β approximate a normal
distribution centered around the optimal parametervalues and then the posterior mean value is close to the
optimal parameter value while other histograms suchas n and d are skewed to the upper value of priorac－ 50 －

srange. The encompassed points around the maximum
value of the prior parameter range presented in Figure
4 (see the plot for correlation between n and β ) infers
that pre-specified prior parameter range for calibration
is not suitable. Although the initial guess (i.e., prioruncertainty range) for model parameters is physicallyreasonable or practically appropriate, unknown hugeparameter interaction leads to the numerically
acceptable values, which are too much different fromhydrologist’s concern.
Figure 4 describes scatter plots in two dimensionsof the parameter space of the 5000 parameter sets
sampled from the posterior distribution afterconvergence to stationary posterior distribution. Due tono contribution of the surface roughness coefficient, n
to subsurface flow simulation, this parameter shows nocorrelation with the other calibration parameters whichis confirmed by the summarizing correlation
coefficients presented in Table 4. On the other hand,dominant parameters of subsurface flow such as k , d ,β show a positive or negative correlation (see the
Figure 4 and Table 4).
Moreover, the optimal parameter set estimated by
the SCE-UA method is included in Table 4. Gupta et al.
(1998) introduced a multi objective analysis
framework to investigate deficiencies in the model
structure, which are reflected in a structural inability tosimultaneously reproduce different aspects of thesystem response with a single set of parameters.
Different parameter combinations are required to fit
different response modes such as low and high flow
(Bastidas, 1998). However, in KsEdgeFC2D model
case, we did not detect a big difference of hydrologicalresponses or a distinguished improvement for different
modes of hydrographs according to objective functions.
As shown in Table 4, optimal parameter values of not
only the SCEM-UA but also the SCE-UA algorithm
are approximately similar regardless of three objectivefunctions. It means that the KsEdgeFC2D model has
stable and appropriate structure to simulate
rainfall-runoff processes. Also, this result indicates that
the Pareto solutions estimated by Multi objective
optimization algorithm for KsEdgeFC2D model might
have a very small and narrow range. Hydrographs
reproduced by two algorithms using each optimal
parameter set are presented in Figure 5. Because theresiduals between the observed streamflow and theassimulated ones are assumed to be mutuallyindependent, Gaussian distributed, with constantvariance, all best-fit parameters estimated by
SCEM-UA are much closer to the calibrated
parameters using SLS objective function than othersfrom HMLE, MIA (See Figure 3).
Finally, probabilistic predictions of the hydrographwere obtained from the ensemble simulation of theKsEdgeFC2D model for 5000 parameter combinations
sampled from the posterior parameter distribution.
Figures 6 and 7 illustrate how the results of the
SCEM-UA algorithm can be translated into estimatesof hydrograph prediction uncertainty. In these figures,
the black line indicates the observed streamflow data,
the blue line, the simulated hydrograph using the most
likely parameter set having the highest posterior
probability and the grey shaded region, 90%hydrograph prediction uncertainty associated with theposterior distribution of the parameter estimates. Here,parameter uncertainty boundary estimated by theposterior distribution is narrow and failed to bracket
the observations for the most calibration period;
particularly, the recession part of simulated hydrographis not matched to the observed one. This result may
come from the imperfect model structure and/or input
data. Figure 7 shows that the parameter sets evaluatedfrom calibration event reproduce the observed
discharges reasonably well during validation period.
6. Summary
Identifiability process including a selection ofsuitable model and an estimation of proper parametervalues is difficult due to a range of uncertainties
involved in a modeling process that are also
unavoidably propagated into the model output. Even
with the robust global automatic optimization
algorithms, classical parameter estimation approaches
are not able to treat parameter uncertainty and its
influence on hydrograph simulation. However, the
SCEM-UA, an automatic Bayesian parameter
inference algorithm based on Markov Chain Monte
Carlo methods, has been proven to be very efficient for
estimation of target posterior parameter distribution. In
this study, we demonstrated the capability of the SCE
M-UA algorithm to calibrate the KsEdgeFC2D modelfor the Kamishiiba catchment. Summarized resultsfrom this study are as follows.

－ 51 －

Fig. 3 Marginal posterior probability distribution of the KsEdgeFC2D model parameters
Table 4 Shuffled Complex Evolution Metropolis Posterior Mean (Mean), Standard Deviation (STDV) andCorrelation Coefficients between the generated samples for the KsEdgeFC2D parameters
Parameter

Mean

Correlation coefficients

STDV

Optimal parameter values

n

ka

ds

dc

β

SCEM

SCE_SLS

SCE_HMLE

SCE_MIA

n

0.497

0.052

1

-0.033

-0.13

0.05

-0.06

0.499

0.5

0.5

0.5

ka

0.037

0.004

-

1

-0.61

0.09

-0.69

0.037

0.038

0.03

0.047

ds

0.684

0.045

-

-

1

0.19

0.87

0.691

0.688

0.682

0.683

dc

0.589

0.050

-

-

-

1

0.14

0.599

0.6

0.6

0.6

β

6.12

1.058

-

-

-

1

6.082

6.046

5.882

6.642

－ 52 －

Fig. 4 Scatter plots in two dimensional parameter space of 5000 sampled parameter sets showing the correlation ofdifferent combinations of parameters (Correlation Coefficients for all combinations are illustrated in Table 4)

Fig. 5 Simulated Hydrographs using two (SCEM-UA, SCE-UA with 3 different objective functions) algorithms－ 53 －

Fig. 6 Hydrograph prediction uncertainty associated with the most probable set derived using the SCEM-UA
algorithm for the calibration period (9/22/1999).

Fig. 7 Hydrograph prediction uncertainty associated with the most probable set of calibration event derived using
the SCEM-UA algorithm for the validation period (9/15/1997).

－ 54 －

1) The SCEM-UA algorithm was able to successfully
explore the feasible parameter space and to converge
into the target posterior parameter distributions after
approximate 5000 iterations.
2) When comparing with the optimal parameter set
calibrated by the original SCE-UA algorithm, the most
likely parameters having the highest frequency werevery similar as those values.
3) The subjectively selected three different objectivefunctions provided approximately constant parameter
sets and the simulated hydrographs based on those
optimal parameter sets are quite similar.
4) Posterior density such as k , d and β followed toa normal distribution while other evaluated parameter
were excessively skewed to the upper value of priorparameter range.
5) Hydrographs predictions based on the posterior
parameter
distributions
demonstrated
that
KsEdgeFC2D model was able to reproduce theobserved discharges with reasonable accuracy for the
Kamishiiba catchment.
6) The parameter uncertainty bounds were narrow anddid not cover all observations. It means that
improvement in model structure or input data may
result in more accurate predictions.
The analysis of reasons why equifinality happen inhydrologic modeling, including the comparison ofposterior parameter distributions in differenthydrologic models and its translation to hydrographsimulation, is on going.
asReferences
Bastidas, L.A. (1998) Parameter estimation forhydrometeorological models using multi-criteria
methods, Ph.D Thesis, University of Arizona, Tucson,
Ariz., U.S.A.
Beven, K.J. (1989) Changing ideas in hydrology – Thecase of physically-based models, Journal of
Hydrology, 105, 157-172
Duan, Q., Sorooshian, S. and Gupta, V.K. (1992)
Effective and efficient global optimization forconceptual rainfall-runoff models, Water Resources
Research 28(4), 1015-1031.
Feyen, L., Vrugt, J.A., Nualláin, B.Ó., Knijff, J. and Roo,
A.D. (2006) Parameter optimization and uncertainty
assessment for large-scale streamflow simulation with
the LISFLOOD model, Journal of Hydrology,

doi:10.1016/j.jhydrol.2007.07.004.
Gelman, A. and Rubin, D.B. (1992) Inference fromiterative simulation using multiple sequences,
Statistical Science 7(4), 457-472.
Gupta, H.V., Sorooshian, S. and Yapo, P.O. (1998)
Toward improved calibration of hydrologic models:Multiple and noncommensurable measures of
information, Water Resources Research 34(4),751-763.
Gupta, H.V., Sorroshian, S and Yapo, P.O (1999) Statusof automatic calibration for hydrologic models:comparison with multilevel expert calibration, Journalof Hydrologic Engineering, ASCE, Vol. 4, No. 2,
135-143.
Ichikawa, Y., Murakami, M., Tachikawa, Y., Shiiba, M.
(2001) Development of a basin runoff simulation
system based on a new digital topographic model,
Journal of Hydraulic, Coastal and Environmental
Engineering, JSCE, no. 691/II-57, 43-52.
Kuczera, G. and Mroczkowski, M. (1998) Assessmentof hydrologic parameter uncertainty and the worth ofmulti-response data, Water Resources Research 34(6),1481-1490.
Kuczera, G. and Parent, E. (1998) Monte-Carloassessment of parameter uncertainty in conceptual
catchment models, Journal of Hydrology 211, 69-85.
Kuczera, G. and Franks, S.W. (2002) Testing Hydrologic
Models: Fortification or Falsification?, Mathematical
Models of Large Watershed Hydrology (edited bySingh, V.P. and Frevert, D.K.), Water Resources
Publications.
Melching, C.S., Yen, B.C. and Wenzel, H.G. (1990) A
reliability estimation in modeling watershed runoff
with uncertainties, Water Resources Research 26(10),2275-2286.
Muletha, M.K. and Nicklow, J.W. (2005) Sensitivity anduncertainty analysis coupled with automatic
calibration for a distributed watershed model, Journalof Hydrology 306, 127-145.
Lee, G., Tachikawa, Y. and Takara, K. (2007)
Identification of model structural stability throughcomparison of hydrologic models, Annual Journal of
Hydraulic Engineering, JSCE, vol. 51, 49-54.
Schaake, J.C. (2003) Introduction of calibration ofwatershed models, Calibration of watershed models(edited by Duan, Q., Gupta., H.V., Sorroshian, S.,
Rousseau, A.N. and Turcotte, R.), Water Science andApplication 6, AGU, Washington, 1-7.

－ 55 －

Sorooshian, S. and Dracup, J.A. (1980), Stochastic
parameter estimation procedures for hydrologic
rainfall-runoff
models:Corrolated
andheteroscedastic error cases, Water Resources
Research 16(2), 430-442.
Sorooshian, S. and Gupta, V.K. (1995) Modelcalibration, Computer Models of Watershed
Hydrology (ed. Singh, V.P.), Water Resources
Publications, Colorado, 23-68.
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,
7-12.
Uhlenbrook, S., Seibert, J., Leibundgut, C. and Rodhe,
A. (1999) Prediction uncertainty of conceptual
rainfall-runoff models caused by problems in

identifying model parameters and structure,
Hydrological Sciences Journal 44, 779-797.
Vrugt, J.A., Gupta, H.V., Bouten, W., Sorroshian, S.
(2003) A Shuffled Complex Evolution Metropolis
algorithm for optimization and uncertainty assessmentof hydrologic parameter estimation, Water Resources
Research 39(8), 1201. doi:10.1029/2002WR001642.
Vrugt, J.A., Diks, C.G.H., Gupta, H.V., Bouten, W. and
Verstraten, J.M. (2005) Improved treatment of
uncertainty in hydrologic modeling: combining thestrengths of global optimization and data assimilation,
Water
Resources
Research
41,
W01017.
doi:10.1029/2004WR003059.
Wagener, T., Wheater., H.S. and Gupta, H.V. (2004)
Rainfall-Runoff Modelling in Gauged and Ungauged
Catchments, Imperial college press, London.

分布型降雨流出モデルにおけるパラメータの不確実性の定量化

Giha LEE*・立川康人*・宝馨* 京都大学大学院工学研究科都市環境工学専攻

要 旨
一般的に水文モデルには直接計測することのできないモデルパラメータが含まれ，それらの値は
過去の水文データに適合するように決定される。これまで数十年に渡って，モデルパラメータの自動推定に関する研究が多数なされてきたが，それらの手法はパラメータ推定に伴う不確かを考慮できないという欠点を有する。本研究では，分布型降雨流出モデルKsedgeFC2Dのモデルパラメータ推定について，パラメータ推定の不確かさとその河川流量予測への影響を，上椎葉流域(211㎢)を対象に分析する。パラメータの不確かさが予測値にどのように伝達するかを分析するために，
Shuffled Complex Evolution Metropolis アルゴリズム(SCEM-UA)を採用した。SCEM-UAを用いることにより，効率的にモデルパラメータの事後分布が得られることがわかった。また，SCEM-UAによって求められたパラメータの値と3つの異なる目的関数を設定したSCE-UA法によって得た値とを比較し，SCEM-UAの適用性を検討した。
キーワード: 自動パラメータ推定，パラメータの不確かさ, SCEM-UA－ 56 －