雑記

日々の記録になります．

【流出解析】タンクモデルのプログラム作成

はじめに

先週は貯留関数法を用いて河川流域からの流出量を計算するプログラムを作成した（【流出解析】貯留関数法のプログラムを作ってみる - 雑記）．流出量を計算するモデルは貯留関数法以外にもあり，タイトル中のタンクモデルもその一つである．貯留関数法に続いて，本稿ではタンクモデルのプログラムを作成した．

タンクモデルについて

タンクモデルは，流域を，側面にいくつかの流出孔を持つタンクに置き換えて考える流出計算手法である．降雨はタンクモデルの最上段のタンクに入り，2段目以下のタンクは1段上のタンクの孔から水を受ける．各段のタンク内の水の一部は側面の孔から外部へ流出し，これらの和が河川流域からの流出量となる．

f:id:okome29taberu:20200502201131p:plain:w300


タンクモデルでは一般に直列3～4段のものが利用される．最上段が洪水の表面流出，2段目が中間流出，3，4段目が地下水流出に対応すると考えられている．この理由から，洪水解析用タンクモデルは2～3段，低水解析用タンクモデルは3～4段で構成されることが多い．

基礎式

上図で考えると，まず，1段目のタンクに降った雨 r mm/hはタンクに一時的に貯留され，e_1 mm/hで蒸発散し，貯留量 S_1 mmと流出孔の高さ h_1 および h_2 mmに応じた流量 q_{11} mm/h，q_{12} mm/h，g_1 mm/hが，タンクからそれぞれ外部もしくは下部タンクへ流出する．このことから，Δt 時間における1段目のタンクの貯留量変化 ΔS_1は，


ΔS_{1}=rΔt-e_{1} Δt-q_{11} Δt-q_{12} Δt-g_{1} Δt

と表せられるので，以下の式が得られる．


\frac{dS_1}{dt}=r-e_1-Q_1-g_1

ただし，Q_1=q_{11}+q_{12} であり，これが1段目のタンクの基礎式である．一方，2段目以降のタンクでは，水の供給源は降雨の代わりに上部タンクからの流入（浸透）g_iであるので，上式と同様に考えると以下の基礎式が得られる．


\frac{dS_i}{dt}=g_{i-1}-e_{i}-Q_{i}-g_{i}


Q_{i}=q_{i1}+q_{i2}+⋯+q_{ij}

ここで，S_i は i 段目のタンクの貯留量， e_{i} は i 段目のタンクからの蒸発散量， Q_{i} は i 段目のタンクからの流出量， q_{ij} は i 段目のタンクにおける j 番目の孔からの流出量である． 一方で，流出量 q_{ij} と貯留量 S_i および孔の高さ h_{ij} の間には以下の式の関係があると仮定する．

{\displaystyle q_{ij}= \begin{eqnarray} \left\{ \begin{array}{l} 0\hspace{38mm}(S_{i}≤h_{ij}) \\ a_{ij}(S_{i}-h_{ij})\hspace{12mm}(S_i>h_{ij}) \end{array} \right. \end{eqnarray} }


ここで，a_{ij} は i 段目のタンクにおける j 番目の孔のパラメータである．この式から，i 段目のタンクにおいて，貯留量 S_{i} が孔の高さ h_{ij} よりも低い場合，j 個目の孔からは流出しないということが分かる．同様に下部タンクへの流出量（浸透量） g_i は，


g_{i}=b_{i} S_{i}

であらわされる．ここで，b_{i} は i 段目のタンクにおけるパラメータである．

パラメータ

本稿では，下図の菅原ら（1986）に記載されている河川用パラメータを用いる．

f:id:okome29taberu:20200502191032p:plain:w150


流出孔の高さ h については，以下の通りである．

流出孔の位置 h [mm]
1段目タンクの上側孔 25~60
1段目タンクの下側孔 5~15
2段目タンクの孔 0~30
3段目タンクの孔 0~60
4段目タンクの孔 0
計算方法

流出量 Q_{i} および q_{ij}，浸透量 g_{i} の計算式は前述したので，貯留量 S_{i} の計算式だけ書き記す．S_{i} についての基礎式において時間方向に前進差分をとると，


\frac{S_{i}^{t+1}-S_{i}^{t}}{Δt}=g_{i-1}^{t}-e_{i}^{t}-Q_{i}^{t}-g_{i}^{t}

となる．ここで，t はタイムステップである． S_{i}^{t+1} ついてまとめると，


S_{i}^{t+1}=S_{i}^{t}+(g_{i-1}^{t}-e_{i}^{t}-Q_{i}^{t}-g_{i}^{t})Δt

となる．ただし，1断面目のタンクにおいては g^{t}_{i-1} の代わりに降水量 r^{t} を用いる． 実際の計算の際は，1段目のタンクから順に各孔からの流出量を計算した後に，各タンクの貯留量を計算するという手順となる．

プログラムと計算結果

本稿では，【流出解析】貯留関数法のプログラムを作ってみる - 雑記の際と同様に，降雨量データとして，財）北海道河川防災研究センター・研究所（2006）に記載されている湧別川・丸瀬布地点における2001年9月10日～15日に発生した洪水のデータを利用した．ただしモデルでの蒸発散量は0とした（用いた降雨量データは有効雨量であるから，蒸発散がすでに実際の降雨量から引かれていると仮定した）．また各流出孔のパラメータは菅原ら（1986）に従い，流出孔の高さはその中央値とした．

プログラムは以下の通りである．一応C++で書いた．

#include <iostream>
#include <iomanip>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

using namespace std;

// Parameter setting
// 各タンクの流出孔の係数
double a11=0.20;    // 1段目タンクの上側孔
double a12=0.20;    // 1段目タンクの下側孔
double b1=0.20;     // 1段目タンクの下向き孔
double a21=0.050;   // 2段目タンクの孔
double b2=0.050;    // 2段目タンクの下向き孔
double a31=0.010;   // 3段目タンクの孔
double b3=0.010;    // 3段目タンクの下向き孔
double a41=0.001;   // 4段目タンクの孔

// 各タンクにおける流出孔の高さ(mm)
double h11=42.5;    // 1段目タンクの上側孔
double h12=7.5;     // 1段目タンクの下側孔
double h21=15.;     // 2段目タンクの孔
double h31=2.5;     // 3段目タンクの孔
double h41=0.;      // 4段目タンクの孔

// その他のパラメータ
double dt=1.;       // time step [h]

// 4段タンクモデル
void TankModel4steps(double& r, double& q, double& s1, double& s2, double& s3, double& s4){

  double q11=0.;                     // 1段目タンクの上側孔
  double q12=0.;                     // 2段目タンクの下側孔
  double g1=0.;                      // 1段目タンクの下向き孔
  double q21=0.;                     // 2段目タンクの孔
  double g2=0.;                      // 2段目タンクの下向き孔
  double q31=0.;                     // 3段目タンクの孔
  double g3=0.;                      // 3段目タンクの下向き孔
  double q41=0.;                     // 4段目タンクの孔

  // 各流出孔からの流出量の計算
  if(s1>h11) q11=a11*(s1-h11);    // 1段目タンクの上側孔
  if(s1>h12) q12=a12*(s1-h12);    // 1段目タンクの下側孔
  if(s2>h21) q21=a21*(s2-h21);    // 2段目タンクの孔
  if(s3>h31) q31=a31*(s3-h31);    // 3段目タンクの孔
  if(s4>h41) q41=a41*(s4-h41);    // 4段目タンクの孔
  g1=b1*s1;                       // 1段目タンクの下向き孔
  g2=b2*s2;                       // 2段目タンクの下向き孔
  g3=b3*s3;                       // 3段目タンクの下向き孔

  // 合計流出量の計算
  q=q11+q12+q21+q31+q41;

  // 各タンクの貯留量の計算
  s1+=dt*(r-q11-q12-g1);          // 1段目タンク
  s2+=dt*(g1-q21-g2);             // 2段目タンク
  s3+=dt*(g2-q31-g3);             // 3段目タンク
  s4+=dt*(g3-q41);                // 4段目タンク

}

int main (){

  char fileName[256];     // File name
  string ss;                  // Dammy line
  double s1,s2,s3,s4;   // 貯留量[mm]
  double q;             // 流出量[mm/h]
  double qobs;          // 実測流出量[mm/h]
  double r;             // 時間tからt+1における有降雨量[mm/h]
  double tmax=94;       // 計算時間[h]

  // 入力データファイルの設定
  ifstream ifs("./data.txt");
  if(!ifs){
    cout<<"cannot open file"<<endl;
    exit(1);
  }
  getline(ifs,ss);

  // 出力データファイルの設定
  ofstream ofs("./out.txt");
  ofs<<"time[h]\train[mm/h]\tQobs[mm/h]\tQsim[mm/h]"<<endl;

  //貯留量の初期値
  s1=0.;
  s2=0.;
  s3=0.;
  s4=0.;

  // 計算開始
  for(int t=1;t<tmax; t++){
    ifs>>ss>>r>>qobs;                   // 時間tからt+1における有降雨量の読み込み
    TankModel4steps(r,q,s1,s2,s3,s4);   // 4段タンクモデルでの計算

    // 計算結果の出力
    ofs<<std::fixed<<std::setprecision(4)<<t<<"\t"<<r<<"\t"<<qobs<<"\t"<<q<<endl;
    cout<<t<<endl;
  }
  return 0;
}


計算結果は下図の通りである．明らかに計算値は実測値に一致していないが，計算プログラムは合っている（はず）．これは用いたパラメータが今回の対象流域に対して適切でないためである．パラメータの最適化は今後の課題とする．


f:id:okome29taberu:20200502194041p:plain:w600


引用文献
  1. 菅原正巳, 渡辺一郎, 尾崎睿子, 勝山ヨシ子 (1986) パーソナル・コンピュータのためのタンク・モデル・プログラムとその使い方, 国立防災科学技術センター研究報告, 37.
  2. （財）北海道河川防災研究センター・研究所（2006）実践流出解析ゼミQ&A編