Journal of Japan Society of Civil Engineers, Ser.
Annual Journalof HydraulicEngineering,
2013, February
B1 (HydraulicEngineering),
Vol. JSCE,
69, No.Vol.57,
4, I_163-I_168,
2013.

ESTIMATING THE 2011 LARGEST FLOODDISCHARGE AT THE KUMANO RIVER
USING A 2D DYNAMIC WAVE MODEL
AND PARTICLE FILTERS
Yeonsu KIM1, Yasuto TACHIKAWA2, Sunmin KIM3,
Michiharu SHIIBA4, Seong Jin NOH1, and Takashi MENO1,
1Student Member of JSCE, Dept. of, Civil and Earth Resources Eng., Kyoto Univ. (Kyoto 615-8540, Japan)
2Member of JSCE, Assoc. Prof., Dept. of Civil and Earth Resources Eng., Kyoto Univ. (Kyoto 615-8540, Japan)
3Member of JSCE, Lecturer, Dept. of Civil and Earth Resources Eng., Kyoto Univ. (Kyoto 615-8540, Japan)
4Member of JSCE, Prof., Dept. of Civil and Earth Resources Eng., Kyoto Univ. (Kyoto 615-8540, Japan)

The largest-ever flood event was recorded at the Kumano River in the Kii peninsula of Japan from
August 31 to September 4, 2011. Water level data exceeded the observation range in many places except
in the Ouga water stage station at downstream. This study aimed to estimate the peak discharge of theevent by using the estimation method, which is a 2D dynamic wave model combined with particle filters,
and by considering the water level observed at the Hitari and Ouga stations and a river discharge
estimated by a hydrological model. The estimation method was applied to the three historic flood events
for quantifying the uncertainties of the Manning roughness coefficient (Manning’s n) of the subject river
channel. Based on the optimized range of Manning’s n, the filtering method was applied to estimate thelargest peak discharge of the 2011 flood using the discharge estimated by a hydrological model. The
possible range of the largest peak discharge was successfully evaluated through the comparison of theobserved flood marks. Finally, a rating curve established by the estimation results at the Ouga station is
examined.
Key Words: 2D dynamic wave model, particle filters, Kumano River, 2011 flood, rating curve

1. INTRODUCTION
From August 31 to September 4, 2011, the largest
ever flood event occurred in the Kii peninsula, Japandue to the 12th typhoon in 2011. Although the peakdischarge was able to be estimated using a rating curveor a hydrological model1), its accuracy is controversial
because of the uncertainties in the rating curve at the
Ouga station, and prediction uncertainties of ahydrologic model for the largest flood we have everexperienced. Historic flood events utilized in thehydrologic model calibration1) are all located in theextrapolation ranges in the rating curve.
Prediction uncertainties of discharge come fromthe various characteristics of the open channel flowsuch as channel geomorphology, channel roughness,initial and boundary conditions and so on. Consideringsuch factors properly, it makes sense that a hydraulicmodel and a hydrologic model should be fully utilized
to estimate river discharge combined with filtering
techniques and all available observed information.
To estimate river discharge using the hydraulicmodel, one essential prerequisite is the boundary
condition such as inflow from the upper reach. The
upper boundary inflow is possible to estimate using ahydrologic model; however, the uncertainty includedin the estimated boundary inflow should be treatedproperly with filtering techniques. To consider theuncertainties of the boundary conditions, Tachikawa etal.2) and Kim et al.3) introduced particle filters (PFs)
into 1-dimensional and 2-dimensional hydraulicmodels, respectively, and successfully provided
reliable river stages and discharges within theuncertain upstream inflow data. In the estimation
methods to deal with the uncertainties of the discharge,
available information such as water level was used as
reference data to evaluate the uncertainties and update
the state variables. Kim et al.3) adopted a 2D hydraulicmodel to consider the channel geomorphology more
precisely, which incorporated reference water level
data more effectively than a 1D model.
Simultaneous estimation of two variables such as
channel roughness and upper boundary inflow
discharge is a delicate problem because their

I_163

Annual Journal of Hydraulic Engineering, JSCE, Vol.57, 2013, February

analysis based on the previous study7) is applied firstto quantify a proper range of channel roughness of thestudy river using the events occurring in August, 2003,
August, 2004 and September, 2004. Based on the
optimized model parameters, a distributed hydrologic
model1) and a 2D hydraulic model are fully utilizedwith limited water level information and flood marks
to estimate the discharge of the largest flood in 2011.

2. APPLICATION PROCESSAs shown in Fig. 1, the whole process is composedof two parts: a 2D hydraulic model and particle filters.
The connector “A” and “B” indicate reading boundary
conditions and reading reference data, respectively.
The basic frame shown in Fig. 1 is the method
proposed by Kim et al. 7), which is used to estimate thespatial distribution of the channel roughness and theupper boundary inflow discharge at the same time.
Kim et al.7) applied the method to estimate channelroughness using discharge converted from a rating
curve adding noise as the upper boundary condition.
At the upper boundary of the channel reach in thissubject area, the Hitari station, water level data isobserved, although a rating curve is not established.
Therefore, the upper boundary inflow discharge wasestimated by a distributed hydrological model1) as analternative.
The initial range of Manning’s n is sequentially
updated as the process progresses. The analysis
process consists of three applications for the different
events in August 2003, August 2004, and September
2004, respectively. At the first application for theAugust 2003 flood, the initial range of Manning’s n isset by certain values for the main channel, flood plain,
and inundation area. Then, from the second application
for the August 2004 flood, the initial range ofManning’s n is set by the 90% interval estimated at the
final step of the previous application.
Fig. 1 also shows the process to estimate thedischarge of 2011 flood. The estimation method is thesame as that applied to the three historic floods.
However, the resampling step stops after the upstream
water level used as the reference data does not exist
due to exceed the observation range. In case of no data
as a reference data of the estimation process, theprocess proceeds without resampling. The distributionof the inflow boundary discharge and Manning’s nfollow the distribution of the final step result when
there is resampling.
For the application of the method to the largest
flood in 2011, an upper boundary inflow is used as anaveraged hydrograph estimated by a distributedhydrological model with different parameter sets
optimized for different events1). A noise term is added
to the inflow, and it is used for the upper boundary
condition. Manning’s n is used for the finallyestimated value from the third event occurring in

Fig. 1 Flow chart of the estimation process.

combination can be almost unlimited for a given
observed water level. Thus Hsu et al.4), Ding et al.5),
and Kim et al.6) only focused on one of two factors:
either channel roughness or boundary inflow.
However, the boundary inflow and the channelroughness should be estimated simultaneously
considering that the channel roughness is affected
continuously by the evolving channel geomorphology
and the vegetation distribution of the flood plain.
In our previous study7), it was shown that boundary
inflow and channel roughness were able to be
estimated simultaneously in a reasonable range using a
2-D hydraulic model combined with particle filters.
Based on several historic flood events, the uncertain
ranges of channel roughness can be identified and themost reasonable boundary inflow can be estimated. Byutilizing this information, it is possible to estimate theproper river discharge of the following event, eventhough it has the largest ever flood peak.
In September 2011, the largest flood happened at
the Kumano River basin resulting in an exceeded
observation range at various water stage stations withthe exception of the Ouga station. The rating curve at
the Ouga station could have large uncertainties
because the observed water level in the 2011 flood far
exceeds the applicable range of the existing rating
curve. Therefore, this study aimed to estimate thedischarge of the flood event in September, 2011 by
fully utilizing a 2D hydraulic model, a distributedhydrologic model and filtering techniques with all
available information. A sequential uncertainties

2

I_164

Annual Journal of Hydraulic Engineering, JSCE, Vol.57, 2013, FebruaryZ  by N uN vN gh s   (uvh)  (v2 h)y  x
y
t
x
y

(3)

where h is water depth; Z s is water stage; t istime; u and v are depth averaged velocity in x andy direction; M and N denote discharge flux in x andy direction;  is density of water; and  bx and  by
are the bed shear stress in x and y direction.

(3) Particle Filters(PFs)PFs perform the sequential Monte Carlo (SMC)
estimation based on particle representations ofprobability densities within Bayesian theorem9). The
main purpose of PFs is to track a target variable as it
evolves over time, typically with a non-Gaussian andmulti-modal
probability
densityfunction10).
Particularly, PFs have special advantages in dealing
with non-linear systems. Among the various PFs,sequential importance resampling (SIR) is introducedin this study9).
In general, PFs require a number of particles for
estimation and prediction, which sometimes causes
equifinality problems. Kim et al.7) introduced the
reduction factor and the variance tuition factor into theperturbation step to reduce the effect of the equifinality
problems which are induced by the interaction of thevariables(e.g., discharge and Manning’s n) whenutilizing the water level data as reference data. The
perturbation step consists of two steps: one is at theinitial stage and another is at the updating stage. At theinitial stage, the perturbation equations for eachvariable are as below:

Fig. 2 Study area and calculation domain.

September 2004, which is used for the initial range ofthe channel roughness to estimate the discharge of thelargest event in 2011.

3. SIMULTANUOUS ESTIMATION OF
DISCHARGE AND CHANNEL ROUGHNESS(1) Study area
The study area is the Kumano River located in theKii peninsula, Japan. The total reach length is about
11km from the Hitari water stage station to the Ouga
water stage station, which are located at upstream anddownstream boundaries, respectively. The calculation
domain is composed of 212 longitudinal points and 88
cross sectional points for a 2D hydraulic model (Fig.
2). The channel roughness is classified into three parts:
main channel, flood plain, and inundation area. These
parts are classified by the aerial photo. As upperboundary inflow conditions, the upstream discharge
from the tributary (the Akagi River) and the main
channel (the Kumano River) are considered. The
lateral inflow to the channel reach is disregarded
because the largest volume per second is smaller than
600m3/s even in the largest flood in 2011 as estimatedby the distributed hydrologic model.

H ti , down  H tobs , down   t ,Qi , int

Qhydro , intn

 t ,
i , maint(4)

 t ~ U(-0.1,0.1 ) t ~ U(0.7 Q t,

obs , int

,1.3Q

obs , in
t

)

(5)
(6)

 t ~ U( ,  )where U indicates uniform distribution; Hi , downt

andare the downstream water level of eachparticle and the observed water level at downstream,
respectively; Qti ,in and Qthydro ,in are the upstreamH

obs , downtdischarge of particle i and the discharge obtained
from a hydrological model, respectively; n ti ,main is theManning’s n of particle i at time t for the main
channel; and  and  are determined by 0.02 to
0.06, 0.02 to 0.08, and 0.04 to 0.1 to the main channel,
flood plain, and inundation area, respectively.
At the updating stage the perturbation equations
for each variable are as follows:(2) 2D dynamic wave model
The 2D dynamic wave model is composed of thecontinuity equation (Eq. (1)) and the momentum
equation (Eq. (2) and (3)). The equations are
numerically solved with the Finite Volume Method(FVM) on a structure grid8). The simple first order
upwind scheme is utilized for the convection term. The
Adams-Bashforth method is introduced for time
integration8).

(7), downH ti,down H tobs
  t 1 ,  t 1 ~ N(0,0.1 2 )
1
1

,in
,in 2 (8)Qti,in1  Qti ,in  a t (Qthydro ,in  Qthydro)   t 1 ,  t 1 ~ N(0, (0.1Qthydro
) )
1
1

nti1  nti ,main   t 1 ,  t 1 ~ N(0, s 2Varth M N(1)
 0
t x yZ M uM vM gh s  bx  (u2 h)  (uvh) (2)txyx  xyin case Varn, main 0.002 , n
2

i
t 1

n

i , maintn , main(9))

  t 1 ,  t 1 ~ N(0,0.002 )2(10)where N indicate normal distribution; at is therecursive correction factor, which is determined by
3

I_165

(a) EV_2003(1)(b) EV_2004A(2)(c) EV_2004B(3)Fig. 3 The comparison between the estimated water level and the observed water level at upstream.

(a) EV_2003(1)(b) EV_2004A(2)(c) EV_2004B(3)Fig. 4 The estimated Manning’s n at the main channel.

(a) EV_2003(1)(b) EV_2004A(2)(c) EV_2004B(3)Fig. 5 The estimated Manning’s n at the flood plain.
,inat  (Qtin1 / Qthydro Qtin / Qthydro ,in ) / 2 ;1sis the tuition

factor for the variance reduction proposed by
Moradkhani et al.11); and Vart n ,main is the variancecalculated at the current step for the weight averaged
values.
In case that Vart n ,main is smaller than the 0.0022,
Eq. (10) is applied to maintain the minimum tracking
ability and to consider the change of Manning’s n as
the channel evolves. In addition, the same equation isapplied to the flood plain and inundation area.
Manning’s n for the specific area is excluded from thefiltering target when the number of wetted grids issmaller than 10% of the number of total grids of eachclassified zone.
At the resampling step, the likelihood of eachparticle against the observed water stage is calculatedusing the Gaussian function:w  p(h
i
t

obs
t

|h )
i
t

1

 2

 ( hti  htobs ) 2

e

2 2

(11)

where wti indicates the likelihood of the i th particle

at time t ;  is the standard deviation associatedwith the observed stage( htobs ) and determined by 0.2m;

and hti is the simulated water stage in each particle.
Then, the normalized weight ( W i ) is calculated byt

the likelihood as follows:
wi
t
Wi 
t
N
i
 wt
i 1

(12)

(4) Estimation results
For the purpose of identifying the channelroughness of the river channel, the method was applied
sequentially to the flood events in August 2003,
August 2004, and September 2004. Through thesequential applications with 120 particles, Manning’s
n is continuously quantified. Figs. 3, 4 and 5 show
estimated upstream water levels and Manning’s n forthe main channel and flood plain, respectively. Figures(a), (b), and (c) in each graph indicate the eventsoccurring in August 2003, August 2004, andSeptember 2004, respectively. In each figure, the cyan
and the yellow color area indicate the 70% and 90%
intervals of the particles. The black line means weightaveraged values, and red points indicate the value ofeach particle. The black points in Fig. 3 are theobserved water level used for reference data in the

I_166

Annual Journal of Hydraulic EnginAeering, JSCE, Vol.57, 2013, February

estimation process to calculate the weight of eachparticle. The estimated water level follows theobserved water level as the process progresses asshown in Fig. 3.
In terms of Manning’s n, the initial range at main
channel, flood plain and inundation area is set at 0.02
to 0.06, 0.02 to 0.08 and 0.04 to 0.1, respectively.
From the second application, the initial range ofManning’s n is set by the 90% interval of the final stepin the previous application (Fig. 4 and 5). As shown in
Fig. 4 and 5, the variation of Manning’s n is
continuously reduced, but the minimum errors aremaintained to keep the tracking ability. In addition, theestimated Manning’s n for each section shows good
agreement with the values reported by Chow12) and
DEFRA/EA13).

Fig. 6 Discharge estimated from the estimation process at the
Ouga station.

4. ESTIMATIING 2011 FLOOD DISCHARGE
The method shown in Fig. 1 is applied to estimate
the discharge of the event occurring in September,2011. The observed water level at upstream exists until
5am on September 3, which is roughly 24 hours before
flood peak. Then the water level at upstream exceeds
the measurable range of the station. Therefore theestimation process is proceeding by the method
proposed by Kim et al.7) before 5am on September 3.
The initial range of Manning’s n is set by the values,
which are 0.021~0.033, 0.040~0.047, and 0.057~0.08
at the main channel, flood plain, and inundation area,
respectively, from the previous analysis process. After
6am on September 3, there is no reference data, so theweight of each particle is fixed by the values at the lasttime step of the update estimation, and is maintained
for the subsequent simulation.
In terms of an upstream boundary condition, an
averaged hydrograph for the main channel and thetributary (the Akagi River) from the hydrologicalmodel with different optimized model parameters1) areconsidered. The difference between each hydrograph issmaller than 15% of the given discharge at each time
step.
The estimated discharge at the Ouga station isplotted in Fig. 6. The updating procedure is
implemented in the section where the percentile
interval is located. For the verification of the peakdischarge in 2011 flood, we utilized the observed flood
marks provided by Wakayama prefecture. The flood
marks are compared with the highest water level ineach particle of the estimation process in Fig. 7. As
shown in Fig. 7, almost all of the flood marks arewithin the lower parts of the highest water level ofparticles. This tendency is induced by the effect of thehydrological model outputs because there is alimitation to predicting discharge uncertainties fromthe time without updating. The level of flood marks atthe Hitari station is about 41.4m, and the particles
including the highest water level within 0.5m from41.4m are tracked. From the corresponding particles,

Fig. 7 Comparison between the flood marks and the highest
water level estimated from the estimation method.

Fig. 8 Comparison of an inundation area, which is the upper
part of the study area(as shown in Fig. 2 by a blue line) and ismarked by a red broken line to compare a simulated inundation
area with an observed inundation area(bottom right).

the peak discharge at the Ouga station is estimated in
the range from 22,500m3/s to 25,500m3/s when
considering errors by ±0.5m from the flood marks(Fig. 6 and 7). Then the inundation areas estimated by
the 2D hydraulic model to have the highest water
surface profile similar with flood marks in Fig. 7 arecompared with the inundation map provided by
Wakayama prefecture. The red broken line indicatesthe inundation area, and it shows good agreement withthe simulated results (Fig. 8).

5

I_167

Annual Journal of Hydraulic Engineering, JSCE, Vol.57, 2013, February

range from 22,500m3/s to 25,500m3/s. In addition, therating curve established from the estimation process iscompared to the currently updated rating curve.
In further research, the estimation process would
be combined with a hydrological model to improve thepredictability of the water stage and the uncertainties
of the hydrologic parameters.
ACKNOWLEDGEMENT
This works is based on data from Wakayama prefecture and
J-POWER CO. Ltd. We gratefully acknowledge this support.

REFERENCES
1) Meno, T., Tachikawa, Y., Kim, S., Shiiba, M., Yorozu, K., :Analysus of the Kumano river flood by the Typhoon 12, 2011,
Proc. of the 67th Annual Conference of the JSCE, Ⅱ-160, 2012.
2) Tachikawa, Y., Sudo, J., Shiiba, M., Yorouzu, K., and Kim, S.:
Development of a real-time rivers stage forecasting method
using particle filter, Annual Journal of Hydraulic Engineering,
JSCE, Vol.55, pp. S_511-S_516, 2011.
3) Kim, Y., Tachikawa, Y., Kim, S., Shiiba, M., Yorozu, K. andNoh, S.J.: Short term prediction of water level and dischargeusing a 2D dynamic wave model with particle filters, Journal ofJapan Society of Civil Engineering, JSCE, 68(4), pp. I_25-I_30,
2012.
4) Hsu, M.H., Fu, J.C., and Liu, W.C.: Dynamic routing model withreal time roughness updating for flood forecasting, Journal of
Hydraulic Engineering, 132, pp.605-619, 2006.
5) Ding, Y., Jia, Y., Wang, S.S.Y.: Identification of Manning's
roughness coefficients in shallow water flows, Journal of
Hydraulic Engineering, ASCE, 130( 6), pp. 501-510, 2004.
6) Kim, W., Kim, Y.S. and Woo, H.Y.: Estimation of channelroughness coefficients in the Han River using unsteady flow
model, Journal of Korea Water Resources Association, KWRA,
28(6), pp. 133-146, 1995.
7) Kim, Y., Tachikawa, Y. Shiiba, M., Kim, S., Noh, S.J.: (In press)
Simultaeneous estimation of inflow and channel roughness
using 2D hydraulic model and particle filters. Journal of floodrisk management, accepted in 2012.
8) Nagata, T.: Hydraulics formulae: Hydraulics worked examples
with CD-ROM, Hydraulic committee, JSCE, pp.16-19, 2002.
9) Ristic, B., Arulampalam, S., and Gordon, N.: Beyond the Kalman
Filter, Artech House, pp. 35-65, 2004.
10) Ioannis Rekleitis.: A Particle Filter Tutorial for Mobile Robot
Localization. Technical Report TR-CIM-04-02, Centre forIntelligent Machines, McGill University, Montreal, Quebec,
Canada, 2004.
11) Moradkhani, H., Hsu, K. L. and Sorooshian, S.: Uncertainty
assessment of hydrologic model states and parameters,
Sequential data assimilation using the particle filter, Water
Resour.Res., Vol.41, W05012, doi:10.1029/2004WR00364,
2005.
12) Chow, V. T.: Open-channel hydraulics, New York, McGrawHill Book Co. 1959.
13) DEFRA/EA.: Reducing uncertainty in river flood conveyance,
Roughness Review, R&D Technical Report to DEFRA /Environment Agency, Project W5A- 057, HR Wallingford Ltd.,
United Kingdom, pp. 144–147, 2003.

Fig. 9 Comparison between an established rating curve and theexisting rating curve at the Ouga station.

In addition, the discharge estimated from theestimation process and the observed water level atdownstream are utilized to examine the rating curve at
Ouga station. The rating curve is treated as thefollowing power function:(13)
Q  c(a  h)
where c and a are the constants, and  is anexponent. They define the unique relation between
discharge and water level of the open channel.
The rating curves ((1), (2) and (3) in Fig. 9) areestablished from the estimated results of the three
events. Then the rating curve considering the estimatedresults of 2011 flood and the three events, the old
rating curve and the current rating curve developed by
MLIT are compared. Through the comparison, it isconfirmed that the rating curve established from theestimated results is similar to the currently used rating
curve. It shows that sequential applications of theestimation process are able to establish a rating curve
comparable with the current rating curve. Using theestablished rating curve and the observed water level,
the peak discharge at the Ouga station is estimated byabout 25,300m3/s.

5. CONCLUSION
An estimation method using a 2D dynamic wave
model and particles filters is introduced to estimate thedischarge of the event occurring in September, 2011.
Before estimating the discharge of the event, theuncertainties of the channel roughness for the main
channel, flood plain, and inundation area arequantified by the sequential application of the method.
With the quantified channel roughness, inflow to theupper boundary is also estimated from distributedhydrological model outputs. The discharge of thelargest event is estimated by considering theuncertainties of roughness coefficients and upperboundary inflows. For the verification of the peakdischarge, the highest water level of each particle iscompared with flood marks. By the comparison to theflood marks, the peak discharge is estimated in the(Received September 30, 2012)6

I_168