こ う えいフ ォ ー ラ ム 第 27 号 / 2019 . 3

多地点水位情報と一次元不定流解析を用いたデータ同化モデルの開発
DEVELOPMENT OF DATA ASSIMILATION MODEL USING ONE-DIMENSIONAL UNSTEADYFLOW MODEL AND WATER LEVEL DATA AT MULTIPLE-POINTS
西口MTI-HandBook

MTI*研究領域における計算機科学亮太
・壇 鉄也 ∗ ・泉 典洋 **Ryota NISHIGUCHI, Tetsuya DAN and Norihiro IZUMIFor real-time prediction of river water level, data assimilation model using one-dimensional unsteady flow model
was studied. Adjoint sensitivity analysis was applied to a numerical method
データ同化の考え方とその方法of optimization of data assimilation. Optimization variables included upstream discharge and dis[ 講演The: 中野慎也
統計数理研究所] an actual flood event of Tama River. The
charge in tributary.
developed (modelwas applied) toassimilation result using water level data at 6 points had good performance. The prediction result
which set assimilated flow distribution as initial condition were also good.

細川 Analysis,
敬祐
Keywords ：電気通信大学
Adjoint Sensitivity
Water Level Prediction, Data Assimilation

1.

はじめに

1

はじめに

近年, 地球惑星科学の研究においてデータ同化という手法が広く用いられるようになってきた.
2017 年 7 月 九州北部豪雨水害、と毎海道・東北豪雨水害、
この文章は, MTI 研究領域において近年行われ年全国各地で大規模な水害が発生している。記憶に新しいているデータ同化の試みを紹介するのではなく,データ同化の原理・方法論に関してより基礎的な年) 7 月に発生した西日本豪雨によところでは、昨年 (2018
知識を概観することを目的とする. 具体的には,り多くの人命が損なわれた。このように、計画で想定して近年、2015 年 9 月 関東・東北豪雨水害、2016 年 8 月 北いる規模を超えるような出水が頻発するなか、河道改修や• データ同化とはそもそも何なのか?
ダム・調節池の整備などによる対策には限界があり、官民• どのような問題を解くのか?

をあげた減災対策の充実が急務となっている。

• 実際にどのようにして問題を解くのか?
これを受けて、国土交通省は「住民目線のソフト対策」を

図 1: データ同化によるアプローチの模式図

図- 1 データ同化の概念 (出典:中野 (2009)pp.X-14) )ミュレーションに実測データを取り入れる手法のことである.
丸はある点の観測データの時間変化を表す。また、曲線群

数値シミュレーションを行うためには, 初期条は計算条件を変えて実施した数値シミュレーション結果を掲げ、水害ハザードマップなどのリスク情報の周知、事前, 境界条件, パラメータなどを与える必要があ件といった基礎的な部分を解説する. 個々の適用事表す。曲線群の内、ピンク色のそれが観測データを最も良. しかし, それらの値は往々にして正確に決めるの行動計画・訓練、臨時情報の即時提供に重点的に取組む例に関しては, この文章で概説した基礎知識を元く再現している。このように、様々な計算条件の内、観測られない場合が多い. したがって, 初期条件・境界に個別に調べていただきたい.
こととしている。また、新たな施策として、洪水時の観測データに馴染む解析値を探索するのがデータ同化である。
条件などをどのように与えるかによって, 数値シに特化した「危機管理型水位計」を向こう 3 年間で 1 級河図では極く単純な例で示しているが、観測点が１点であミュレーションは様々なシナリオを導きうる. そ川の指定区間および 2 級河川に 5,800 箇所、1 級河川の大臣の様子を図 1 に模式的に示す. 横軸は時間を示
2 データ同化とはなんなのか? る必要はなく、観測データが１種類である必要もないこと
管理区間に 3,000 箇所設置することを発表している。
し, 同じ数値シミュレーションを用いても, 初期は言うまでもない。要はそれぞれの観測値に対応する数値
2.1 データ同化の概念
さらに、水防活動や住民等の避難喚起にとって重要な洪条件・境界条件の与え方によって様々なシナリオ解が得られれば良い。観測実績を説明する数値解を求める（たくさんの灰色のライン）が結果として現れる
データ同化は, data assimilation の訳語である.
水予警報につながる水位予測の精度向上が求められている。
ことができれば、それを時間的に外挿した将来予測は高精ことが示されている. この図には, 赤い丸で観測Assimilationという言葉が,“ある民族が移民としこれについては、かねてより様々な検討がなされていると度であることが期待できる。これがデータ同化を行う目的
データが示されているが, データ同化はこれらのて他の民族に同化する” ことを表現する際に用いころであるが、未だ解決すべき課題が多い。
である。
データを数値シミュレーションに埋め込み, 馴染られることからも分かるように, データ同化は「数そこで、著者らは上記のように今後情報量が増える水位, 実際の観測データをうませていくことによって値シミュレーションに実測データを埋め込み, 馴 この際、無限に採り得る計算条件の組み合わせから最適観測データを活用して、水位予測の精度の向上を図ることまく説明する, より尤もらしい推定（ピンクのラ染ませること」を意味する. 簡単に言えば数値シな計算条件を効率的に探索するためには、数学的な手法を1),2).
イン）を探しだすことを目指している。
を目的に検討を進めてきている
平成 21 年度 MTI 研究会 サイエンスセッション 用いて行う必要がある。このための有力な方法の一つとしc⃝
Mesosphere
Thermosphere
Ionosphere(MTI)
ところで、気象の分野では予測の前処理として行う「デーて、変分法を応用した「随伴変数法」があり、気象庁の現Research Group, Japanタ同化」を洗練することによって、予報精度を大きく向上
業の予測では、この方法で気象観測値と計算値のデータ同させている 3) 。
X–1化を行っている。ちなみに、気象の数理モデルは空間三次
4)
図-1 はデータ同化の概念 を示すものである。図中の赤
元の気圧・気温・風速など変量の時間変化を計算するもの

なので、気象場のデータ同化手法は「四次元変分法」と呼* コンサルタント国内事業本部 インフラマネジメントセンター** 北海道大学教授

環境フィールド工学専攻

ばれている。

15

多地点水位情報と一次元不定流解析を用いたデータ同化モデルの開発

これに倣って、本研究では観測水位と計算水位のデータ同化手法として随伴変数法を採り上げた。
我々以外の河川分野における随伴変数法の研究実績とし
5)

6)(= H(A))、ie：摩擦勾配、qL：単位幅横流入量である。運動
方程式において横流入に伴う運動量の変化は無視した。摩
擦勾配はマニング則により計算した。なお、マニングの粗

7)ては、Sanders et al. 、Ding et al. 、Elhanafy et al. 、 度係数は時間的に変化させずに一定値として取り扱う。
吉田ら 8),9),10) 、渡邊ら 11),12) などの研究があるが、何れも式 (2) の解 U = (A, Q)（以下、「状態量」という）は初モデル・パラメータの推定を目的としており、予測問題に焦期条件と境界条件 D（以下,「制御変数」という）によって点を当てたものは見当たらない。また、これらのいくつか変化するが、ここでは常流の流れ場を想定して、上流端流では多次元の計算モデルを対象に検討を行っている。これ

量 QB 、下流端の河積 AB を境界条件とする。これに初期に対して、我々は現実的な計算機資源と予測更新時間間隔の流量 Q0 、水位 A0 および単位幅横流入量 qL を加えた 5で実行可能な方法であることを念頭において、一次元モデ変数を制御変数として、()
D = QB (t), AB (t), Q0 (x), A0 (x), qL (x, t)

ルを選び、洪水流の予測問題に適用した。以下、本稿では洪水流に対する随伴変数法によるデータ同化に必要な数式の導出を行い、実河川における適用例を示したうえで、今後の展開について記す。

(3)

である。

()
この問題は、等式制約条件 R D, U (D) = 0 のもとで、評価関数 F (D) を最小化する制御変数 D を求める最適化問題である。この種の問題はラグランジュの未定乗数法に解析手法

2.

計算負荷を抑えて予測を迅速に行うため、流れ場のモデルは一次元不定流モデルとした。
前記の通り、データ同化手法は随伴変数法を用いる。同法の適用に当たって、観測流量は観測水位の換算値であることから、これを用いることなく観測水位のみでデータ同化を行うこととした。
以下、随伴変数法による一次元不定流モデルのデータ同化について記す。

(1)

より、次の J を最小化する無制約問題に置き換えられる。
∫ T∫ X()
J=L D, λ, U (D) dxdt
(4)00L = f + λ · R,

λ = (A∗ , Q∗ )(5)

λ はラグランジュの未定乗数あるいは「随伴変数」という。さらに、J は関数を実数又は複素数に変換する関数（「汎関数」という）であるから、「変分法」によって最適解を求めることができる。すなわち、次式で定義される J の第一変分 δJ をゼロで等値することで最適解を求めることがで

随伴変数法

ここでの評価関数 F は、式 (1) の通り、全地点の観測水位より算出した河積 Aobs と計算結果による河積 A の差の二乗和で与える。これを最小化することでデータ同化が成される。ここで、水位ではなく河積を用いる理由は、各観測きる。

δJ =∫ T∫ X(00∂L∂L
δD +δU∂D∂U)

dxdt ≡ 0(6)L が D に関する非線形関数なので、J の最小化は反復地点の河道形状 (川幅) による重みを付加するためである。 計算によるが、D を更新する際に感度係数 δJ/δD が必要∫ T∫ Xになる。そこで、)2)w(∫ T∫ X(A − Aobs
F =(1)
f dxdt , f =∂L2
δUdxdt = 0(7)0 0∂U0 0
ここに、w は観測データがある時間・空間座標で 1、そ�
れ以外で 0 をとる。なお、ここでは観測点間で観測精度に とおいて λ を求め、これを用いて次式で感度係数を計算する。
)]∂ [ (� �δJ
差がないことを仮定した。
�=L D, λ, U (D)
(8)
δD
∂D
上式中の計算水位は次式で与えられる一次元不定流モデ
例えば、D の更新を最急降下法で行う場合、ルの支配方程式を満足する必要がある。

∂A ∂Q+−qL∂t∂xR==0( 2) ∂Q∂ Q∂H++ gA+ gAie∂t∂x A∂x(2)ここに、t：時間座標（0 ≤ t ≤ T ）、x：空間座標（0 ≤x ≤ X ）、Q：流量、A：河積、g : 重力加速度、H ：水位
16

� (k) − α� (k+1)= DD(δJδD)(k)(9)である。k は反復回数、α はステップ・サイズとよばれる正の数である。
次節で記すように、式 (7) は「随伴方程式」とよばれる偏微分方程式に変換される。随伴方程式の求解（「逆解析」こ う えいフ ォ ー ラ ム 第 27 号 / 2019 . 3
制御変数を仮定)()(Q2 ∂Q∗∂zb∂ie∂A∗+ g h̄ − 2+ Q∗ g −− ie − A∂tA∂x∂x∂A)(= w A − Aobs

制御変数を更新

不定流計算
（順解析）

感度係数を計算NO

評価関数が最小？

∂A∗
Q ∂Q∗∂ie∂Q∗++2− Q∗ gA=0∂t∂xA ∂x∂Q

随伴変数方程式の求解（逆解析）

YESt = T の境界では、次式 (12) の条件が設定可能である。

データ同化

データ同化完了A∗ (x, T ) = Q∗ (x, T ) = 0

図- 2 随伴変数法の流れう）と同程度であり、D の全ての要素に対する感度係数が一度に計算されるので、要素毎に摂動を与える直接法に対して、大幅に演算時間が短縮できる。このように、感度係数を効率的に計算するためのアルゴリズムが「随伴変数法」である。
図-2 は随伴変数法の計算の流れを示したものである。
随伴方程式・感度係数

式 (7) は、部分積分法を適用して展開すると、以下のように変形される。

)
∂L∂L
δA +
δQ dxdt =∂A∂Q0 0∫ T∫ X{ () }
w A − Aobs δA dxdt

∫ T∫ X(0+0∫ X[0A∗ δA∫ T∫ X[{]T0dx +∫ T[02(Q∗ g h̄ −伴方程式の求解を「逆解析」と称する所以である。
上下流端の境界条件は、ディリクレ条件、ノイマン条件のいずれも与えることができない。そのため、本稿では、無反射境界条件を設定した。
式 (11) により求めた随伴変数 A∗ および Q∗ を用いて、J の変分は下式で計算される。
)]∫ T[ (∫ X()Q2A∗ δA0 dx +δJ = −Q∗ g h̄ − 2 δAB dt
AB00])∫ X∫ T[(()Q∗ QB−A∗ + 2
δQB dtQ∗ δQ0 dx −A00
∫ T∫ X()−A∗ δqL dxdt
(13)0

したがって、それぞれの制御変数の感度係数は以下で与えられる。
2QA2)

δA]X

δJ= −A∗ (x, 0)δA0 (x)][Q2 (X, t)
δJ= Q∗ (X, t) g h̄(X, t) −δAB (t)A2B (t)δJ= −Q∗ (x, 0)δQ0 (x)δJQ∗ (0, t) QB (t)= −A∗ (0, t) − 2
δQB (t)A(0, t)δJ= −A∗ (x, t)
δqL (x, t)

dt0

ここに、zb は河床高、h̄ は平均水深である。
式 (10) のうち、時空間の境界を除く、内部の状態量 Uに関する変分の項を 0 とおくことにより、式 (11) が導出される。これが、一次元不定流モデルに対する随伴方程式である。

随伴方程式は時間的に退行して解くことになる。これが随0Q ∂Q∗∂Q∗∂A∗− 2+ g h̄∂tA∂x∂x0 0)} ](∂ie∂zb+ Q∗ g −− ie − A
δA dxdt∂x∂A) ]X∫ X[∫ T [(]TQ∗ Q+A∗ + 2Q∗ δQ dx +δQ dt
A0000∫ T∫ X {(∂Q∗Q ∂Q∗∂A∗−++2∂x∂tA ∂x0 0) }∂ie− Q∗ gA
δQ dxdt = 0(10)∂Q−(12)

これは、随伴方程式 (11) の初期条件にあたる。そのため、という）に要する計算規模は式 (2) の求解（「順解析」とい(2)(11)

(3)

(14)

数値解法

式 (2) および式 (11) の離散化は、定義の煩雑さを避けるため、全ての変数を計算格子中央に配置する集中格子を採用した。
式 (2) の離散化は、自然河道断面および集中格子の計算に適した Ying et al.13) によるスキームを採用した。
一方、式 (11) の離散化は、式 (2) との整合性を考慮して、
式 (2) の離散式を 状態量 U で微分することにより導出し

17

多地点水位情報と一次元不定流解析を用いたデータ同化モデルの開発

'HQHQFKRIXXQGHUZLHUNS
:DWHU/HYHO>$3P@awaAkikakawAsa'HQHQFKRIXRYHUZLHUNS
:DWHU/HYHO>$3P@Range:0.0k-60.0k

7DPDJDZDNS
:DWHU/HYHO>$3P@

Map tiles by Geospatial Information Authority of Japan図- 3 解析範囲位置図
た。式 (11) の境界条件は、任意の値を設定できないため、


上下流端ともに無反射境界条件により設定した。

,VKLKDUDNS
:DWHU/HYHO>$3P@また、感度係数による D の更新には降下法を用いる。降下法には、共役勾配法を使用した。
実河川での適用3.

+LQREDVKLNS


:DWHU/HYHO>$3P@

多摩川の下流 60 km を計算区間とした (図-3)。この区間には約 10km 間隔で 7 地点と、他河川と比較して高密度・多点で水位観測がなされている。検討対象の出水は、過去10 年間で最も規模が大きい 2007 年 9 月の台風 9 号による







それを選んだ。

&KRIXEDVKLNS

データ同化の条件

最適化の対象は、上流端流量の時間分布および初期の流量・流積の空間分布に加えて、流域面積が大きい浅川およ2EVHUYHG

$VVLPLODWHG河口の水位は実績で与えた。よって、データ同化には河

び秋川の横流入を考慮した。

(1)

:DWHU/HYHO>$3P@

口を除く 6 地点の水位データを用いることとした。また、 図- 4 毎正時の水位の同化結果 (同化時刻以前 T 時間の結同化時間 (T ) は、本川上流端・河口間の洪水到達時間に相果を図化)当する 4 時間とした。観測値は 10 分間隔で取得されるので、1 地点当たり 25 時点（= 4 × 6 + 1）のデータで評価関数を計算した。

河道形状は 2007 年度の横断測量成果で与えた。本川には複数の可動堰が存在するが、起算時以降は何れも全開状態にある。
不定流計算の時間刻み、距離刻みはそれぞれ 5 秒、約 200(2)

データ同化の結果毎正時のデータ同化の結果を取り出して観測水位と比較したものが図-4 である。観測水位と計算水位の差は、最大でも 30cm 程度であり、極めて良好な結果といえる。
同様に、図-5 は石原地点・日野橋地点・調布橋地点の観測流量と計算流量を比較したものである。

m で与えた。マニングの粗度係数は、河川計画での採用値

石原地点の実測値は最高水位発生時の 3 時間前から他のを参考に、32.4km より下流では 0.025、これより上流では観測所では見られない振動が含まれている。これに追随し0.030 とした。てデータ同化後の流量が振動している。観測値の質はデー18

こ う えいフ ォ ー ラ ム 第 27 号 / 2019 . 3

,VKLKDUDNS

:DWHU/HYHO>73P@'LVFKDUJH>

m3/s

@+LQREDVKLNS)ORRGSODLQ(OHYDWLRQ

$VVLPLODWHG5LVLQJ$VVLPLODWHG)DOOLQJ

2EVHUYHG'LVFKDUJH>

m3/s

@'LVFKDUJH>P V@図- 6 石原地点の同化流量と水位・流量関数の比較 1)'HQHQFKRIXXQGHUZLHUNS


:DWHU/HYHO>$3P@

&KRIXEDVKLNS

2EVHUYHG@m3/s3

'LVFKDUJH>$VVLPLODWHG


図- 5 毎正時の流量の同化結果 (同化時刻以前 T 時間の結果を図化)

'HQHQFKRIXRYHUZLHUNS

:DWHU/HYHO>$3P@

タ同化の結果に直結するので、観測データの検定を前処理現地の状況にそぐわないものとなる。

粗度係数を 0.025 とする必要があり、礫河床となっている3UHGLFWLRQ7LPHとも観測値に符合させようとすると、全区間で マニングの3UHGLFW

だちに違算とは判じ難い。ちなみに、同地点で水位・流量2EVHUYHG観測水位から換算した値であることを考え合わせると、たが極めて高く、観測流量が実際は水位・流量関数を用いてる。ただし、図-4 で見た通り同地点の計算水位の再現精度調布橋地点では観測値と数 100 m3 /s の差異が生じてい:DWHU/HYHO>$3P@として行うことが必要と考えられる。

7DPDJDZDNS


図- 7 同化時刻から 3 時間先の水位の予測結果

一方、図-7 は石原地点のデータ同化後の計算流量のヒステリシスと水位・流量関数を比較したものである。流量のおよび横流入量を流出予測の計算結果から与える必要があ

再現精度が調布橋地点より高いこの地点では高水敷高以上るが、ここでは便宜的に、これらを予測開始時点のデータの水位のレンジで両者の乖離が小さい。こういった比較を同化の値で一定として与えた。また、河口の水位は実績値行うことで、水位・流量関数の適用性を評価できる。これで与えた。予測時間は 3 時間とした。

はデータ同化を行うことの派生的な効用といえる。

図-7 は、田園調布地点（堰上・堰下）および玉川地点の 2時間毎の予測水位を取り出して図化したものである。デー

(3)

予測計算

水位予測の模擬実験を行った。すなわち、前節で同化した計算水位および計算流量を初期条件・境界条件とする不定流計算を行った。本来であれば、将来時点の本川上流端

タ同化時の精度には劣るものの、ほぼ良好な結果が得られている。
主な精度の悪化要因は支川の扱いによるものと考えられる。浅川および秋川でも水位観測が実施されているので、こ19

多地点水位情報と一次元不定流解析を用いたデータ同化モデルの開発

れを取り込んだ河道網でデータ同化を行うことにより、さらに精度の向上が期待できる。
R
一連の計算に要する時間は、シングルスレッド (Intel⃝Core

TMi7-3770 プロセッサ 動作周波数 3.40 GHz) で 10 分程度であった。これは観測データの入信間隔とほぼ同等なので、今後、演算の高速化が必要である。

参考文献

1) 西口亮太・壇鉄也: 随伴変数法による水位縦断分布のリアルタイム予測に関する研究, 河川技術論文集, 第 23 巻, pp.275-280, 2017.
2) 西口亮太・壇鉄也・泉典洋: 随伴変数法による水位縦断分布のリアルタイム予測に関する研究 (2), 河川技術論文集, 第 24 巻, pp.227-232,2018.

3) 気象庁 報道発表資料: 全球数値予報モデルの改善について～高度な初期値解析手法「４次元変分法」の導入～, http://www.jma.go.jp/jma/press/0502/16a/4jigen.pdf, 2005.

以上、本稿で得られた成果を列挙すると以下のようである。

4) 中 野 慎 也 : デ ー タ 同 化 の 考 え 方 と そ の 方 法 ，平 成 21 年
度 MTI 研 究 会 サ イ エ ン ス セッション MTI 研 究 領 域に お け る 計 算 科 学 口 頭 講 演 資 料 ，pp.X-1 - X-9 http://mti.nict.go.jp/MTI_symposium/mti-handbook/top/h21/
MTI_Handbook_Nakano_Assimilation_Ver_1_2.pdf.

• 気象予測のデータ同化手法として採用され、予測精度

5) Sanders, Brett F and Katopodes, Nikolaos D : Adjoint sensitivity analysis for shallow-water wave control, Journal of Engineering Mechanics, 126-9, pp.909-919, 2000.

適用するため、一次元不定流の支配方程式に対する随6) Ding, Yan and Wang, Sam SY : Identiﬁcation of Manning ’sroughness coeﬃcients in channel network using adjoint analysis,
International Journal of Computational Fluid Dynamics, vol.191, pp.3-13, 2005.

おわりに4.の向上に貢献している随伴変数法を河川の水位予測に伴方程式および感度係数を導いた。

• 上条の結果を用いて、多摩川の大臣管理区間で至近 10年間の最大出水を対象としたデータ同化を試みた。この結果、わずか 6 点の観測水位から高精度な同化結果が得られた。

• この同化結果と流量換算式を比較して、同換算式の適用性の評価に用いることができることを示した。

• データ同化の結果を初期条件・境界条件として水位予測の模擬実験を行った。この結果、今後改善の余地はあるものの、2 ないし 3 時間先までの高精度な水位予測の可能性が示唆された。

• 小さい計算負荷でデータ同化と水位予測の演算を行う実用的な水位予測の枠組みがほぼ確立した。

上述のように、今後は支川の扱いなど河道の数理モデルの改良を加えるとともに、本稿では省略した実時間流出予測の計算を加えて予測時間の延長を図る所存である。
また、最適化に当たって、今回与えた評価関数に変更を加えることで、随伴変数法によるダム放流量の最適化などに応用が可能である。このテーマにも取り組みたい。
謝辞 : 本研究で用いた河川測量成果および水位データは国
土交通省関東地方整備局京浜河川事務所より提供を受けた。
ここに記して謝意を表す。

207) Elhanafy, Hossam, Copeland, Graham JM and Gejadze, Igor Yu: Statistical Modelling of Uncertainties in Flood Wave Propagation Using Adjoint Sensitivity Analysis, Flood Risk Assessment
II: IMA Conference Proceedings: University of Plymouth 4-5 ,pp141, 20078) 吉田圭介・石川忠晴: 変分法と浅水流モデルを併用した河床粗度の逆推定法に関する研究, 水工学論文集 54 号, pp.991-996, 2010.
9) 吉田圭介・石川忠晴: Adjoint 法による流量ハイドログラフ推定法に関する研究, 土木学会論文集 B1(水工学) Vol.68, No.4, pp.I 1264I 1266, 2012.
10) Yoshida,Keisuke and Ishikawa,Tadaharu: Flood hydrographestimation using an adjoint shallow-water model, Journal ofHydro-environment Research, vol.9-3, pp.429–440, 2015.
11) 渡邊明英・見上哲章・小島崇・松延和彦・鈴田裕三・富澤慎二郎: 平
面二次元流解析とアジョイント法に基づいた点観測の水位情報に対する縦断水面形時間変化の同化手法の検討, 河川技術論文集, 第 23 巻,
pp.197-202, 2017.
12) 渡邊明英・見上哲章・小島崇・松延和彦・鈴田裕三: 不確実な入力条件に対する河川縦断水面形の同化解析推定量とその分布, 土木学会論文集 B1(水工学) Vol.74, No.4, pp.I 727-I 732, 2018.
13) Ying, Xinya, Khan, Abdul A and Wang, Sam SY : Upwind
conservative scheme for the Saint Venant equations, Journal of
hydraulic engineering, vol.130-10, pp.977–987, 2004.