Python, matplotlibでFlood Routing（洪水時貯水池水位追跡）

はじめに

洪水時の貯水池水位追跡（英語では flood routing）は，洪水が貯水池に流入した時，ダムからの放流量を考慮しながら，水位がどのように変化するかを追跡する解析です．特に湛水面積・容量の大きな貯水池では，洪水流入時の貯水池の貯留効果が大きいため，ある程度の貯水池の水位上昇は伴いますが，流入洪水波形のピークをカットできるとともに，合理的な洪水吐設計を行うことができるため，重要な解析となってきます．また繰り返し計算を伴う解析であるため，プログラム言語による処理が適していると考えています．

筆者も，色々なケースでこの解析を行ってきましたが，プログラムの簡単な解説を行いながら，ここにアップしてみることにしました．

理論

貯水池への貯留量は，流入量と流出量の差に支配され，ある時間間隔$\Delta t$に対し，この関係は，以下のとおり表現できる．

$$
\Delta S=Q_{i}\cdot \Delta t - Q_{o}\cdot \Delta t
$$

• $\Delta t$ : 時間間隔
• $\Delta S$ : $\Delta t$時間での貯水池への貯留量
• $Q_i$ : $\Delta t$時間内の平均流入量
• $Q_o$ : $\Delta t$時間内の平均流出量

上式を解く上で，以下の３つの特性値が必要となる．

• 貯水池の水位ー容量曲線
• 貯水池への流入量時刻歴：通常は洪水のハイドログラフ
• 貯水池からの流出特性：通常は洪水吐の水位ー流出量関係

解析手順は以下のとおり．

(Step-0) 初期値として，初期貯水池水位標高，初期流入量を与える．初期水位標高により，貯水池水位ー容量曲線から初期貯水池貯留量が算定できる．

(Step-1) $\Delta t$時間後の流入量を入力するとともに，水位上昇量（貯水池水位）を仮定し流出量の計算を行う．

(Step-2) $\Delta t$時間前の流入量と流出量を用い，$\Delta t$時間内の流入体積と流出体積を求める．この流入体積と流出体積の差分が，$\Delta t$時間内の貯水池への貯留量増分となる．

(Step-3) 貯留量増分に前回貯水池貯留量を加えることにより，最新の貯水池貯留量が計算され，また貯水池水位ー容量曲線より貯水池水位も計算できる．

(Step-4) (Step-1)において水位上昇量を仮定しているので，(Step-3)で求めた貯水池水位と(Step-1)で仮定した貯水池水位を比較し，これらの差分が許容誤差以内になるまで，(Step-1)から(Step-4)の過程を繰り返す．

なお，通常の場合，洪水初期は貯水池水位および流出量は上昇していき，流入量のピークを過ぎても貯水池水位および流出量は上昇し続けるが，ある時点から水位および流出量は減少し始める．設計上は，この水位および流出量のピークを捉えれば良いわけであるが，解析としてグラフの体裁を考えると水位減少時の挙動も追跡してあるほうがかっこいいので，水位降下時を追跡できるルーチンもプログラムに組み込んでおく．

解析事例

• 貯水池は発電用であり，ゲートを有している．
• 初期水位はEL.130mであり，ゲート部越流頂標高はEL.114mとしている．
• 事例では，ピーク流量12500m3/sの洪水波形が流入した時の最高水位および流出量を追跡している．
• ゲートを全開にした時，水位標高EL.130mだとすれば越流量は約4300m3/sであるため，流入量が4300m3/s以下の場合はゲート操作により貯水位をEL.130mに保ち，それ以上となったらゲート全開の条件で自由越流させることにしている．
• また発電用貯水池なので，洪水ピークが去り貯水位低下が始まってからは，水位標高がEL.130mに達する直前からゲート操作を開始し，貯水位をEL.130mに保つ運用をするように設定している．すなわち貯水位EL.130mとなってからは流入量＝放流量となることとしている．

プログラム

プログラムは，解析用と作図用に分けている．

以下は実行用スクリプト

python3 py_eng_floodrm.py hv1m.txt hydro1m.txt out_flood.txt
python3 py_eng_fig_floodrm.py

解析プログラム

入力データ：貯水池H-Vデータ

• １行目は貯水池容量に乗じる係数（ここでは百万m3）．
• ２行目以降は，１列目が標高，２列めが累計容量を示す．

hv1m.txt

1e6
 60.0         0.000
 61.0         0.045
 62.0         0.187
 63.0         0.437
..........

入力データ：洪水波形（ハイドログラフ）

• １行目１列目は初期貯水池水位標高（ここでは130m）．
• １行目２列目は常時放流量（ここでは０）．
• ２行目以降は，１列目が時刻，２列目が流入量（m3/s）

hydro1m.txt

130.0 0
  0.0           118
  1.0           205
  2.0           289
  3.0           371
  4.0           454
  5.0           537
..........

解析プログラム

関数 def SP(elv,hmax): および def SPQ(elv,qqin,hmax,qref,ELini): は洪水吐きの特性を示すものであり，基本的には色々な貯水池のタイプに合わせて，この２つを書き換えながら使いまわしている．（洪水吐の特性によっては本ルーチンも書き換える必要あり）

py_eng_floodrm.py

# -*- coding: utf-8 -*-
import sys
from math import *

def SP(elv,hmax):
    # Head-Discharge relationship of Spillway
    ELcr=116.0
    b=11.5
    n=4
    x=0.1387
    Cd=1.971+0.498*x+6.63*x*x
    a=(1.6-Cd)/(Cd-3.2)
    h=elv-ELcr
    q=0.0
    if 0.0<h:
        Hd=139.0-ELcr
        c=1.6*(1.0+2.0*a*h/Hd)/(1.0+a*h/Hd)
        kp=0.112703962703963*(h/Hd)**2-0.237885780885781*(h/Hd)+0.126496503496503
        ka=-0.0813160731673404*(h/Hd)**2+0.255143813465017*(h/Hd)
#        kp=0.02
#        ka=10.0*kp
        bx=b*float(n)-2.0*(float(n-1)*kp+ka)*h
        q=c*bx*h**1.5
    return q


def SPQ(elv,qqin,hmax,qref,ELini):
    # To obtain the discharge of overflow crest from the water level
    q=qqin
    qs=SP(elv,hmax)
    if qref<qqin: q=qs
    if qref+0.1<qs: q=qs
    return q

def RET_V(nn,rh,rv,elv):
    # To obtain the cumulative volume from the water level
    for i in range(0,nn-1):
        if rh[i]<=elv and elv<=rh[i+1]: break
    if rh[nn-1]<elv: i=nn-2
    x1=rv[i]
    y1=rh[i]
    x2=rv[i+1]
    y2=rh[i+1]
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    v=(elv-b)/a
    return v

def RET_H(nn,rh,rv,v):
    # To obtain the water level from cumulative volume
    for i in range(0,nn-1):
        if rv[i]<=v and v<=rv[i+1]: break
    if rv[nn-1]<v: i=nn-2
    x1=rv[i]
    y1=rh[i]
    x2=rv[i+1]
    y2=rh[i+1]
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    elv=a*v+b
    return elv


#Main routine
param=sys.argv
fnameR1=param[1] # H-V data of reservoir
fnameR2=param[2] # Hydrograph of design flood
fnameW =param[3] # Output file name

# Input H-V data
rh=[]
rv=[]
fin=open(fnameR1,'r')
text=fin.readline()
text=text.strip()
text=text.split()
vcoef=float(text[0])
while 1:
    text=fin.readline()
    if not text: break
    text=text.strip()
    text=text.split()
    rh=rh+[float(text[0])]
    rv=rv+[float(text[1])*vcoef]
fin.close()
nn=len(rh)

# Input time sequence of inflow
ti=[]
q_in=[]
fin=open(fnameR2,'r')
text=fin.readline()
text=text.strip()
text=text.split()
ELini=float(text[0])
outlet=float(text[1])
while 1:
    text=fin.readline()
    if not text: break
    text=text.strip()
    text=text.split()
    ti=ti+[float(text[0])]
    q_in=q_in+[float(text[1])]
fin.close()
nd=len(ti)

hmax=-99.9
for iii in range(0,10):
    qref=SP(ELini,hmax)
    pEL=[]
    pQo=[]
    fout=open(fnameW,'w')
    # Initial outflow
    elv=ELini
    EL=elv
    VOL=RET_V(nn,rh,rv,elv)
    q_out=SPQ(elv,q_in[0],hmax,qref,ELini)+outlet
    i=0
    iUD=0
    print('{0:>5s} {1:>5s} {2:>10s} {3:>10s} {4:>15s} {5:>15s} {6:>15s} {7:>15s}'.format('i','iUD','time','EL','EL-elv','VOL','q_in','q_out'),file=fout)
    print('{0:5d} {1:5d} {2:10.3f} {3:10.3f} {4:15e} {5:15e} {6:15e} {7:15e}'.format(i,iUD,ti[i],EL,EL-elv,VOL,q_in[i],q_out),file=fout)
    # Iterative calculation
    iUD=1
    dh=0.0001
    eps=0.0001
    itmax=int(1.0/dh)*10
    for i in range(0,nd-1):
        qqin=0.5*(q_in[i]+q_in[i+1])
        tim=0.5*(ti[i+1]-ti[i])
        Qin=0.5*(q_in[i]+q_in[i+1])*(ti[i+1]-ti[i])*3600.0
        hh=0.0
        j=0
        while 1:
            j=j+1;
            hh=hh+float(iUD)*dh
            elv=EL;    q1=SPQ(elv,qqin,hmax,qref,ELini)+outlet
            elv=EL+hh; q2=SPQ(elv,qqin,hmax,qref,ELini)+outlet
            Qout=0.5*(q1+q2)*(ti[i+1]-ti[i])*3600.0
            R=VOL+Qin-Qout
            elv=RET_H(nn,rh,rv,R)
            if abs(q1-qqin)<eps and tim<36.0:
                R=VOL
                hh=0.0
                break
            if iUD==1 and j==10:
                if EL+hh>elv:
                    iUD=-1
                    hh=0.0
                    j=0
            if iUD==-1 and j==10:
                if EL+hh<elv:
                    iUD=1
                    hh=0.0
                    j=0
            if abs(EL+hh-elv)<eps: break
        VOL=R    # Cumulative volume
        EL=EL+hh # Elevation
        q_out=q2 # Outflow
        print('{0:5d} {1:5d} {2:10.3f} {3:10.3f} {4:15e} {5:15e} {6:15e} {7:15e}'.format(i+1,iUD,ti[i+1],EL,EL-elv,VOL,q_in[i+1],q_out),file=fout)
        pEL=pEL+[EL]
        pQo=pQo+[q_out]
        sys.stdout.write('\r%d / %d' % (i+1,nd-1))
        sys.stdout.flush()
    fout.close()
    if abs(ELini+hmax-max(pEL))<0.001: break
    hmax=max(pEL)-ELini
    sys.stdout.write('\n')
    sys.stdout.write('Time: %10.3f\n'% max(ti))
    sys.stdout.write('h   : %10.3f\n'% hmax)
    sys.stdout.write('Qin : %10.3f %10.3f\n'% (min(q_in),max(q_in)))
    sys.stdout.write('Qout: %10.3f %10.3f\n'% (min(pQo),max(pQo)))
    sys.stdout.write('EL  : %10.3f %10.3f\n'% (min(pEL),max(pEL)))
    sys.stdout.write('qref: %10.3f\n'% (qref))
    del pEL,pQo
sys.stdout.write('\n')
sys.stdout.write('Time: %10.3f\n'% max(ti))
sys.stdout.write('h   : %10.3f\n'% hmax)
sys.stdout.write('Qin : %10.3f %10.3f\n'% (min(q_in),max(q_in)))
sys.stdout.write('Qout: %10.3f %10.3f\n'% (min(pQo),max(pQo)))
sys.stdout.write('EL  : %10.3f %10.3f\n'% (min(pEL),max(pEL)))
sys.stdout.write('qref: %10.3f\n'% (qref))

作図プログラム

解析結果出力ファイル

• i：データ番号
• iUD：水位上昇・降下を示す記号（0：初期値，1：水位上昇，2：水位降下から一定）
• time：時刻
• EL-elv：仮定した水位(EL)と貯水池水位ー容量曲線より求めた水位の差分
• VOL：貯水池の貯留量
• q_in：貯水池への流入量（m3/s）
• q_out：貯水池からの流出量（m3/s）

out_flood.txt

    i   iUD       time         EL          EL-elv             VOL            q_in           q_out
    0     0      0.000    130.000    0.000000e+00    1.053350e+09    1.180000e+02    1.180000e+02
    1     1      1.000    130.000    0.000000e+00    1.053350e+09    2.050000e+02    1.615000e+02
    2     1      2.000    130.000    0.000000e+00    1.053350e+09    2.890000e+02    2.470000e+02
    3     1      3.000    130.000    0.000000e+00    1.053350e+09    3.710000e+02    3.300000e+02
    4     1      4.000    130.000    0.000000e+00    1.053350e+09    4.540000e+02    4.125000e+02
    5     1      5.000    130.000    0.000000e+00    1.053350e+09    5.370000e+02    4.955000e+02
..........

作図プログラム

作図プログラムでは，以下のように，numpy配列でデータを読み込み，計算時間内で，流入量と流出量に大きな差がないことを確認するため，Simpson法による数値積分を行っている．

import numpy as np
from scipy import integrate

data=np.loadtxt(fnameR,skiprows=1,usecols=(2,3,6,7))
ti=data[:,0]
EL=data[:,1]
Qi=data[:,2]
Qo=data[:,3]
int_Qi=integrate.simps(Qi*3600,ti)
int_Qo=integrate.simps(Qo*3600,ti)
print('V_in={0:10.3f} million m3'.format(int_Qi/1e6))
print('Vout={0:10.3f} million m3'.format(int_Qo/1e6))

numpy配列では，最大値となる要素のインデックスは以下のように求める．

n1=np.argmax(Qi)
n2=np.argmax(Qo)
n3=np.argmax(EL)

プログラムでは，twinx() を用いて，左縦軸に流量，右縦軸に貯水池水位をとってグラフ化している．また凡例を１つのボックスに収めるため，右縦軸用のダミープロットを行っている．また点線および一点鎖線は，デフォルトはあまりかっこよくないため使わず，別途定義している．

line1=plt.plot(ti,Qi,color='black',lw=2.0,label='Q (inflow)')
line2=plt.plot(ti,Qo,color='black',lw=2.0,label='Q (outflow)')
plt.plot([0],[0],color='black',lw=2.0,label='Water Level') #Dummy for legend
line1[0].set_dashes([5,2])
line2[0].set_dashes([8,4,2,4])

プログラム全文

py_eng_fig_floodrm.py

# -*- coding: utf-8 -*-
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

def DRAWFIG(fnameR,fnameF,ELcr):
    # Parameter setting
    xmin=0.0;xmax=140.0
    qmin=0.0;qmax=20000.0
    hmin=124.0;hmax=144.0
    # Input data
    data=np.loadtxt(fnameR,skiprows=1,usecols=(2,3,6,7))
    ti=data[:,0]
    EL=data[:,1]
    Qi=data[:,2]
    Qo=data[:,3]
    int_Qi=integrate.simps(Qi*3600,ti)
    int_Qo=integrate.simps(Qo*3600,ti)
    print('V_in={0:10.3f} million m3'.format(int_Qi/1e6))
    print('Vout={0:10.3f} million m3'.format(int_Qo/1e6))

    fig=plt.figure()
    line1=plt.plot(ti,Qi,color='black',lw=2.0,label='Q (inflow)')
    line2=plt.plot(ti,Qo,color='black',lw=2.0,label='Q (outflow)')
    plt.plot([0],[0],color='black',lw=2.0,label='Water Level') #Dummy for legend
    line1[0].set_dashes([5,2])
    line2[0].set_dashes([8,4,2,4])
    n1=np.argmax(Qi)
    n2=np.argmax(Qo)
    n3=np.argmax(EL)
    plt.text(ti[n1],Qi[n1],'max:%.0f'%max(Qi),fontsize=10,color='black',ha='left',va='bottom')
    plt.text(ti[n2],Qo[n2],'max:%.0f'%max(Qo),fontsize=10,color='black',ha='left',va='bottom')
    plt.xlabel('Time (hour)')
    plt.ylabel('Discharge (m$^3$/s)')
    plt.xticks(np.arange(xmin,xmax+12,12))
    plt.yticks(np.arange(qmin,qmax+2000,2000))
    plt.grid()
    plt.legend(shadow=True,loc='upper right',handlelength=3)
    plt.twinx()
    plt.plot(ti,EL,color='black',lw=2.0,label='Water Level')
    plt.hlines([ELcr],xmin,xmax,color='black',lw=1.0,linestyles='dashed',label='O.F.crest')
    plt.text(ti[n3],EL[n3],'ELmax:%.3f'%max(EL),fontsize=10,color='black',ha='left',va='bottom')
    plt.text(xmin+1,ELcr,'FSL:EL.%.1f'%ELcr,fontsize=10,color='black',ha='right',va='bottom')
    plt.ylim(hmin,hmax)
    plt.yticks(np.arange(hmin,hmax+2,2))
    plt.ylabel('Water Level (EL.m)')
    plt.savefig(fnameF,dpi=200)
    plt.show()
    plt.close()


# Main routine
fnameR='out_flood.txt'
fnameF='fig_flood.png'
ELcr=130.0
DRAWFIG(fnameR,fnameF,ELcr)

以上