WANtaroHP (f90: Run-off Analysis)


解析に用いる方程式

貯留関数法は以下の式で構成されている．

基礎方程式

\begin{align} &\cfrac{dS}{dt}=r_e(t-T_L)-q_c(t) & \text{連続の式}\\ &S(t)=k\cdot \{q_c(t)\}^p & \text{貯留量と流出量の関係式}\\ &r_e(t)=f\cdot R(t) & \text{有効雨量算定式} \end{align}
$S(t)$ : 貯留量(mm/hr)     $k$ : 定数
$r_e(t)$ : 有効雨量(mm/hr)     $p$ : 定数
$T_L$ : 遅延時間(hr)     $f$ : 定数(流出率)
$q_c(t)$ : 直接流出量(mm/hr)     $R(t)$ : 観測雨量(mm/hr)

直接流出量算定式

\begin{align} &q_c(t)=q(t)-q_B(t) & \text{直接流出量算定式}\\ &q(t)=\cfrac{3.6\cdot Q(t)}{A} & \text{総流出量算定式}\\ &q_B(t)=q_1+(i-n_1)\cfrac{q_2-q_1}{n_2-n_1} & \text{基底流出量算定式} \end{align}
$q_c(t)$ : 直接流出量(mm/hr)     $Q(t)$ : 観測流量(m$^3$/s)
$q_B(t)$ : 基底流出量(mm/hr)     $A$ : 流域面積(km$^2$)
$q(t)$ : 総流出量(mm/hr)
$n_1$ : 直接流出開始点     $q_1$ : 直接流出開始点流出量
$n_2$ : 直接流出終了点     $q_2$ : 直接流出終了点流出量
$i$ : $n_1$から$n_2$までの数値

流出率算定式

\begin{equation} f=\cfrac{\sum q_c}{\sum R_T} \end{equation}
$f$ : 流出率     $q_c$ : 直接流出量(mm/hr)     $R_T$ : 観測雨量(mm/hr)

ここに，観測雨量とは，時間シフトを考慮した直接流量と同時間内における雨量である．



未知数の求め方と流量再現

連続の式を以下のように差分近似する．

\begin{equation} \cfrac{dS}{dt}=r_e(t-T_L)-q_c(t) \quad\Rightarrow\quad S(t+\Delta t)=S(t)+\Delta t\cdot r_e(t+\Delta t -T_L)-\cfrac{\Delta t}{2}\left\{q_c(t+\Delta t)+q_c(t)\right\} \end{equation}

上式のうち，$q_c$と$r_e$は全時刻で既知であるため，時刻$0$での$S$の初期値$S(0)=0$を設定することにより，各時刻での貯留量$S$を求めることができる．

遅延時間$T_L$については，その値を$T_L=0, 1, 2, \cdots$と変化させ，適切な値を定める．

貯留量$S$と直接流出量$q_c$には以下の関係がある．

\begin{equation} S=k\cdot \{q_c\}^p \quad\Rightarrow\quad \ln S=\ln k+p\cdot \ln q_c \end{equation}

$T_L$を変化させながら両対数軸上で貯留量$S$と直接流出量$q_c$の関係を調べ，流量増加時と流量現象時で挙動の一致度が最も高い時の$T_L$を遅延時間として定める．

$k$および$p$の値は，定めた$T_L$に対し，両対数軸上の直線回帰により定めることができる．

$k$，$p$，$T_L$が定まれば，以下の微分方程式をRunge-Kutta法などにより，数値積分することにより，$q_c$の時刻歴が定まる．

\begin{equation} \cfrac{dq_c}{dt}=\cfrac{1}{k\cdot p}\cdot(r_e-q_c)\cdot q_c{}^{1-p} \end{equation}

なお，上式における$q_c$は，$q_c$-$S$関係を算出した時に用いた$q_c$とは無関係である． この段階では，$r_e$，$p$，$k$が既知量であり，$q_c$は未知量である． このため，直接流出量の分離点$n_1$に相当する時刻に小さな$q_c$の値を入力し，以降逐次数値積分により次回の値を予測していくようにする．

$q_c$の時刻歴が定まれば，以下の式により流量の再現ができる．

\begin{equation} Q=\cfrac{1}{3.6}\cdot A \cdot (q_c+q_B) \end{equation}

ここで，上式適用にあたり，単位には注意を要する． 流量$Q$の単位は m$^3$/s であり，直接流出量$q_c$および基底流出量$q_B$の単位は mm/hr，流域面積$A$の単位は km$^2$ である．



計算事例

以下の6つのケースについて試計算を行い，流量曲線の再現を試みた．

  • Case 1 : 鵡川 (09/Aug/1992)
    (Source: 若手水文学研究会：現場のための水文学(1) -- 流出解析　その1 -- in Japanese character)
  • Case 2 : 教科書の例題
    (Source: 椎葉充晴・立川康人・市川温，例題で学ぶ水文学，森北出版株式会社，20120年5月21日in Japanese character)
  • Case 3 : 忠類川 (05/Aug/1975)
    (Source: 北海道開発局土木試験所河川研究室：実用的な洪水流出計算法，昭和62年3月in Japanese character)
  • Case 4 : 落部川 (29/Jun/1960)
    (Source: 実用的な洪水流出計算法 (洪水データ集) in Japanese character)
  • Case 5 : 落部川 (03/Apr/1961)
    (Source: 実用的な洪水流出計算法 (洪水データ集) in Japanese character)
  • Case 6 : 古舞川 (28/Jun/1966)
    (Source: 実用的な洪水流出計算法 (洪水データ集) in Japanese character)

以降に，貯流量-直接流出量の関係図および流量再現結果図を示す．



Case 1: 鵡川 (1992.08.09)

png png
TL=0TL=3.6hr (adopted)
qc-S relationship
png
nd n1 n2 Delay time f p k
72 3 72 3.6 0.736 0.763 14.364
Hydrograph for Case 1


Case 2: 教科書例題

png png
TL=0TL=3.4hr (adopted)
qc-S relationship
png
nd n1 n2 Delay time f p k
24 3 23 3.4 0.476 1.000 2.868
Hydrograph for Case 2


Case 3: 忠類川 (1975.08.05)

png png
TL=0TL=0.9hr (adopted)
qc-S relationship
png
nd n1 n2 Delay time f p k
36 4 36 0.9 0.231 0.743 6.150
Hydrograph for Case 3


Case 4: 落部川 (1960.06.29)

png png
TL=0TL=0.5hr (adopted)
qc-S relationship
png
nd n1 n2 Delay time f p k
29 2 28 0.5 1.405 1.000 3.218
Hydrograph for Case 4


Case 5: 落部川 (1961.04.03)

png png
TL=0TL=1.0hr (adopted)
qc-S relationship
png
nd n1 n2 Delay time f p k
30 2 24 1.0 1.205 0.783 6.329
Hydrograph for Case 5


Case 6: 古舞川 (1966.06.28)

png png
TL=0TL=3.9hr (adopted)
qc-S relationship
png
nd n1 n2 Delay time f p k
72 8 72 3.9 0.410 0.656 9.598
Hydrograph for Case 6


参考文献

  • 角屋睦，永井昭博：流出解析法 (その10) -- 4.貯留法-貯留関数法による洪水流出解析 --，農業土木学会誌第48巻第10号，pp43-50，1980年10月
    (http://www.scj.go.jp/ja/member/iinkai/bunya/doboku/takamizu/pdf/haifusiryou04-1.pdf)
  • 若手水文学研究会：現場のための水文学(1) -- 流出解析　その1 --
    (http://kankyou.ceri.go.jp/genba/suimon.html)
  • 椎葉充晴・立川康人・市川温，例題で学ぶ水文学，森北出版株式会社，2010年5月21日
    (ISBN 978-4-627-49631-6 published 05.2010)
  • 北海道開発局土木試験所河川研究室：実用的な洪水流出計算法，昭和62年3月
    (http://river.ceri.go.jp/contents/tool/kouzui.html)
  • 実用的な洪水流出計算法 (洪水データ集)
    (http://river.ceri.go.jp/contents/tool/kouzui.html)


プログラムリスト

FilenameDescription
a_f90.txt実行用シェルスクリプト
f90_sfm.f90貯留関数法パラメータ予測
a_gmt_ctl_sfm.txt洪水波形再現用GMTコマンド(制御部)
a_gmt_plt_sfm.txt洪水波形再現用GMTコマンド(作図部)
a_gmt_ctl_qcss.txt直接流出量-貯留量関係作図GMTコマンド(制御部)
a_gmt_plt_qcss.txt直接流出量-貯留量関係作図GMTコマンド(作図部)
dat_sfm1.csv入力データ事例 (1)
dat_sfm2.csv入力データ事例 (2)
dat_sfm3.csv入力データ事例 (3)
dat_sfm4.csv入力データ事例 (4)
dat_sfm5.csv入力データ事例 (5)
dat_sfm6.csv入力データ事例 (6)

プログラミング

Fortran プログラムの実行コマンドは以下のとおり． 

./f90_sfm  fnameR  fnameW  n1  m  md
fnameR入力ファイル名
fnameW出力ファイル名
n1 基底流量分離開始点
m 基底流量分離方法 (0: 水平分離[固定値]，1: 傾斜直線分離[最適値検索])
md 時間増分の分割数 (1 以上の整数)
  • 直接流出開始点 $n_1$ はコマンドライン引数として入力する．
  • 直接流出終了点 $n_2$ の扱いに関しては，コマンドライン引数$m$が入力される．
    • $m=0$の場合，基底流出量 $q_B$ は一定値であり，その値は直接流出開始点での値に等しい．
    • $m=1$の場合，基底流出量 $q_B$ は直接流出開始点 $n_1$ での値から直接流出終了点 $n_2$までの間で変化する．
    • 直接流出終了点 $n_2$ は自動的に検索・設定される．
    • $n_2$ 検索・設定の条件は，ハイドログラフにおける観測流量と再現流量の差分が最小となることである．
    • 実際の計算では，$n_2$ は大きな値から小さな値に変化させられる．
    • 直接流出量 $q_c$ は0以上の値であるべきなので，計算された$q_c$が0未満の場合は，その値は0に修正される．
  • 遅延時間の算定法は以下のとおり．
    • 直接流出量$q_c$と貯留量$S$の関係はループ状になることが知られているけれども，パラメータ$p$および$k$は，両対数グラフ上の回帰分析により計算することができる．
    • 有効雨量の時刻をシフトさせることにより，遅延時間は，観測データと回帰推定値の残差が最小になった時の時間として得ることができる．
    • $p$および$k$を計算するための回帰分析では，$q_c$-$S$関係のループ形状を考慮するため特殊な方法を用いている．そのデータ選別方法を下に示す．
      • $q_c$軸上の領域を小さな領域に分割する．
      • 各々の領域において，$S$の最大値と最小値を抽出する．
      • 各々の領域から抽出された全データを用いて回帰分析を行う．
    • $p$の値の上限値は1.0に規定している．もし$p$の値が計算上1.0を超えた場合，$p$を1.0に固定し，$k$は，$p=1.0$の条件での回帰分析により求める．
  • 数値積分のための時間間隔は変化させることができる．その値$md$は，コマンドライン引数として入力される． この時，流量データはスプライン補間により内挿され，雨量データは1時間の中で同じ値とされる．なぜなら雨量の単位は'mm/h'だからである．
  • $T_L$, $k$ および $p$の値が得られた後，ハイドログラフ再現のための数値計算には，Runge-Kutta法が用いられている．


与えられたパラメータに対するハイドログラフ

プログラム

FilenameDescription
a_f90.txt実行用シェルスクリプト
f90_sfmq.f90洪水波形計算
a_gmt_ctl_sfmQ.txtハイドログラフ作成用GMTコマンド(制御部)
a_gmt_plt_sfmQ.txtハイドログラフ作成用GMTコマンド(作図部)
dat_sfm1Q.csv入力データ事例 (1)
dat_sfm2Q.csv入力データ事例 (2)
dat_sfm3Q.csv入力データ事例 (3)
dat_sfm4Q.csv入力データ事例 (4)
dat_sfm5Q.csv入力データ事例 (5)
dat_sfm6Q.csv入力データ事例 (6)

作図事例

png png png
fig_SFM1Q.pngfig_SFM2Q.pngfig_SFM3Q.png
png png png
fig_SFM4Q.pngfig_SFM5Q.pngfig_SFM6Q.png


inserted by FC2 system