格子型河道水文モデルの演習問題

河道網

東京大学生産技術研究所の沖大幹氏を中心とするグループによる TRIP (Total Runoff Integrating Pathways; OKI and SUD, 1998; http://hydro.iis.u-tokyo.ac.jp/~taikan/TRIPDATA/TRIPDATA.html) を使う。 OKI and SUD (1998)で 紹介されているのは緯度経度1度格子のものであり、ここではそれを使うが、 ほかに0.5度格子のものも作られている(上記hydro内のページ参照)。

各ます目からの流れは縦・横・斜めのいずれかの隣のます(8つのいずれか)に 行くと考える。 便宜上、流れる先を下のようなコードで表わしている。 河口はコード「9」、水域は「0」となっている。

(なお、沖氏のグループで現在作成中の新しいバージョンでは、 流れる先が隣であるという制約を取り除き、 少し離れたますに行くこともありうるとしているそうだ。そうすると、 1度格子のますの中を複数の大河川が通る場合を無理なく扱うことができる。)

対象地域

メコン川流域およびチャオプラヤ川流域を含む、 北緯8〜34度、東経93〜110度の範囲の陸地。

この領域でのTRIPの流路網は次のようになっている。 (文字幅一定のフォントで表示させてほしい)。 □は海、●は河口として扱われているます目、矢印は東西や南北に向かう流路、 ┛などの罫線記号はとがっている向き(この例では右下つまり南東)に向かう流路を示す。

33.5 ←┛┛┛┛┛→→→┏↓┛↓┗┛┛↓
32.5 ┛┏→┛┛┛┛┛↓↑↓↓┛↓↓→→
31.5 →→┛┛↓↓┛┛┛↓↓┛↓↓┗←┓
30.5 ┛┛↓┛┛┛↓┛↓↓↓↓┛↓┗→→
29.5 →→┗┗↓↓↓↓↓→┛┛→→┓←↓
28.5 ┛┛↓┗┛┛┛┛↓↓→→┓↑↑↑→
27.5 ┗┗←←↓↓↓↓↓┗↑↑↑→┓→┓
26.5 ←←↓┗↓↓↓→→┓↑┛┓↓→→┓
25.5 ↓↓┗↓↓┛┛↓↑↑↑┛┛┛┛→↓
24.5 ←↓┗┗←┗┗↓┛↑→→→→┛→┛
23.5 ↓↓┛←←↓←↓┛┛┛┛→→┛→→
22.5 ┗┛↓↓┗↓┛┛↓┛→┛↓→→→┓
21.5 ↓↓┗←↓↓→↓┗↓→┛↓●●□●
20.5 ●┛↓↓↓┗←┛↓┗┛┛●□□□□
19.5 □┛↓↓↓←↓↓↓↓→→●□□□●
18.5 □●↓↓↓↓↓┗┛→┛→●□□□●
17.5 □●↓●↓┛┛↓→┓→↓←●□□□
16.5 □●●●●↓→↓→→┛┛↓←●□□
15.5 □□□□□┛┛↓→→→→↓↓→●□
14.5 □□□□□●↓↓┓┛↓┗↓┗┗●□
13.5 □□□□□●●●●→┛←↓←←●□
12.5 □□□□□□●□□●→→↓←↑●□
11.5 □□□□□□●□□□●┓┛↓→●□
10.5 □□□□□●□□□□●●●●□□□
 9.5 □□□□□□□□□□□□●□□□□
 8.5 □□□□□●●□□□□□□□□□□

モデルの構成

時間ステップΔt=1日で計算する。

空間は実際は球面で考えるべきだが、 便宜上、緯度経度座標を直交直線座標とみなす。 したがって、ます目の面積はどれも同じとなる。 東西の格子間隔は北緯17°での値で代表させる。 つまり、Δy = 2 π r_E/ 360; Δx = Δy × cos(17°)。

TRIPで陸となっているます目だけ計算する。 海となっているます目では実質計算しないが、計算結果のファイルには欠損の印として 「-1.0」という値を入れておく。

陸面過程

簡単化したバケツモデル (真鍋のバケツと違って蒸発量を土壌水分貯留量や 気象要素に関連させない)を使い、 観測に基づく降水量のデータを入力とし、蒸発量は一定値を仮定して、 流出高を求める。

S_soil(t+Δt) = S_soil(t) + (P - E - R)Δt

ここで、Pは観測データの日降水量を与える。 Eは一定値(プログラム中の変数evap、現在の値は2.5 mm/day)を与える。 Rは、現在の版では、

• S_soil(t) + (P - E)Δt > S_soil,maxのとき、 RΔt = S_soil(t) + (P - E)Δt - S_soil,max
• S_soil(t) + (P - E)Δt ≤ S_soil,maxのとき、 R = 0

とする。ただし、S_soil,maxは「バケツの容量」で、 プログラム中の変数bucket、現在の値は150 mmである。

河川水文過程

各ます目のうち河道は狭い面積比をしめるだけだと考えているのだが、 便宜上、河道の水貯留量S_riverは、 ます目に平均に水が配分された場合の深さ[mm] (または面積あたりの質量[kg/m²])で表わす。 (川の水深と直接比較できないことに注意。)

水は、格子型流路網上を、一定速度 v (プログラム中の変数vel、現在の値は0.4 m/s)で動くと仮定する。 たとえばxの正方向(東)に動くとすると、Δtだけの時間には、 S_river×A×v×Δt/Δxだけの水が東隣に移動し、 S_river×A×(1 - v×Δt/Δx)だけの水がもとのます目に残ることになる (ただし、A = Δx×Δy。 v×Δt/Δxが1を越えないようにΔtを選んで計算することにする。) 同様にy方向に動くとすればv×Δt/Δxのところをv×Δt/Δy、 斜めに動くときは v×Δt/sqrt((Δx)²+(Δy)²)とする(sqrtは平方根)。

[注意] ここでは、Δx、Δyの計算はます目の中心の間の距離をそのまま使った。 OKI & SUD (1998)、 松山ほか(1998)、MIYAOKA et al. (1999)では、 実際の河道の長さが1度格子で表現された河道の長さの約1.4倍になるとしているが、 この演習用プログラムではその1.4という因子を入れていない。 そのぶんだけ流速vの意味が違っていると考えればよいと思う。

降水量データ

総合地球環境学研究所の谷田貝亜紀代氏を代表とするプロジェクト 「アジアの水資源への温暖化影響評価のための日降水グリッドデータの作成」 (APHRODITE's Water Resources, http://www.chikyu.ac.jp/precip/)で作成された 東アジアのデータセットAPHRO_EA_V0804を使用する。 このデータセットは、雨量計観測値を空間内挿したものである。 雨量計の系統誤差の補正は含まれていないが、 地形の効果は地域によって(旧版と同じだとすれば中国とモンゴルについて) 考慮されている。 1980年から2002年まで時間間隔1日で、 空間間隔緯度経度0.25度および0.5度のものが作られている。 なお、このデータセットの旧版はXIE et al. (2007)に紹介されている。

ここでは、期間として、 GAME Phase 2 (TAKAHASHI & AGATA, 2006)のデータが取りこまれている1997-2002年を選び、 空間分解能は緯度経度0.5度から緯度経度1度に変換(空間平均)したものを使うことにする。

計算対象期間、初期条件

水文年を陸水貯留量が最小となるときから始めることが多いので、 ここでは計算開始を1997年4月1日とした (プログラム内でjyearb = 1997; jmonb = 4として設定)。 入力データファイルは1月1日から始まるが、1997年3月までのデータは 読みとばす。 計算の終わりは2002年末とした (プログラム内でjyeare = 2002として設定)。 計算開始時の貯留量は、河道水貯留量が0、 土壌水分貯留量が100 mmとした。 (プログラム内の最初のDO la, DO loのループ内で srivo(lo,la) = 0.0; soilm(lo,la) = 100.0としているところで設定。)

比較する流量データ

メコン川下流4か国で作る国際機関である メコン川委員会(Mekong River Commision, http://www.mrcmekong.org/)がまとめたデータから、 次の地点の日流量を用意する。

• Chiang Saen (チエンセーン): タイ, 北緯20.27度、東経100.08度、上流の流域面積189000平方キロメートル、データ期間1960年〜2000年。
• Pakse (パクセ): ラオス, 北緯15.12度、東経105.80度、上流の流域面積545000平方キロメートル、データ期間1923年〜2001年。

また、タイ王立灌漑局(Royal Irrigation Department)から GAME-Tプロジェクトに提供されたデータから、 チャオプラヤ(Chao Phraya)川の次の地点の日流量を用意する。

• RIDの地点番号C.2、Khai Chira Prawat (カイ・チラ・プラワット; Nakhon Sawan (ナコンサワン)の近郊): タイ, 北緯15.67度、東経100.11度、上流の流域面積110569平方キロメートル、データ期間1956年〜1999年。

ただしデータ期間の一部にデータ欠損の期間がある。 ここに置いたファイルは、データが欠損の日は記録そのものがないようにした。

書き出す変数

プログラム中に(だいたい)次のような行が用意してある。

DATA    ifout/1,0,0,0,0/

この列の左から、河川流量、降水量、流出高、河道水貯留量、土壌水分貯留量に 対応し、値が1のとき結果をファイルに書き出し、0のときは書き出さない。

河川流量については、"flow.txt" というひとつのファイルに、 1時刻(1日)ごとに1行ずつ書き続けていく。 各行は、左から、年、月、(月ごとの)日、年の通しの日、 Chiang Saenの流量、Pakseの流量、C.2の流量とする。 ただし、各地点の流量とは、各地点を含むモデルのます目から出ていく流量を さすものとする。 流量の単位は観測に合わせてm³/sとする。 (本来は時間あたりの質量のほうがよいのだが、標準の密度1 g/cm³を 仮定している。)

その他の変数は、1日ごとに、たとえばpr19970401.txtのような 名前のテキストファイルに2次元の場を書き出すようにした。 ただし、prは降水量(mm/day; 入力値の確認のため書き出し)、 roは流出高 (mm/day)、 rwは河道水貯留量 (mm)、 swは土壌水分貯留量 (mm)である。 そのあとの4けたの数字は日付である。 (pr、roはその1日間の平均値、 rw、swはその日の終わりの値と考える。) ファイル内は、次のように構成している。 つまり、「#」以外の字で始まっている行に、緯度に続いて変数値がならぶ。

• ### 変数名
• #@ 日付
• # 経度 経度 ...
• 緯度 変数値 変数値 ...

プログラムの所在

プログラムrouting1.fはFortranで記述してある。 CygWin上ならば、g77でコンパイルして、./aで 実行できるはず。 ただし、実行前に、流路網データファイルrouting.txtも 用意しておくこと。

作図は各自くふうしてほしい。

gawk -f precont.awk pr19970401.txt

load "cont.gplt"

また、毎日のデータを月ごとに平均するAwkプログラムを用意した。 モデル実行後、次のようにすると、変数名rwのデータを平均して、 結果をrw199704.txtのような形の名前のファイルに書き出す。

gawk -v var=rw -f monav.awk /dev/null

データの所在

以下のデータは、834号室のWindowsサーバー上にあり、 834号室のPCからは「Xドライブのフォルダmasuda」に見える。 CygWinの中からは/cygdrive/x/masuda/の下である。

降水量のデータはmasudaの下のprecipの中に置く。 モデルのプログラムはそこを見にいくので、 データを自分のところにコピーする必要はない。

流量の観測データはmasudaの下のflowの中に置く。 地点ごと、年ごとに、たとえば chiang_saen1997.txt、 c2_1997.txtのような名前にしてある。 内容は1時刻(1日)が1行で、左から、年、月、(月のうちの)日、年の通しの日、 流量(単位はm³/s)である。 gnuplotで作図するときには、たとえば

plot "x:/masuda/flow/chiang_saen1997.txt" using 4:5

データはいずれもテキストファイルだが、 改行コードはLFのみ(Unix用)になっている。 CygWin、emacs、gnuplotではそのまま読めるが、 その他のWindows用ソフトウェアでは、 期待どおり改行されるためにはコードをCR+LFに変える必要があるかもしれない。

課題

モデルを実行し、いずれかの観測地点での流量を観測データと比較して、 まず、大まかな季節変化がどの程度合っているかを見てほしい。 さらに、次のいずれかを論じてほしい。

• 流量の観測値の年々変動の特徴をモデルが再現しているか。
• 流速vの値を変えると流量の特徴はどう変わるか。 最適なvの値は見つかるか。ただし、v×Δt/Δxが1を越えない範囲で。
• バケツの深さbucketを変えると流量の特徴はどう変わるか。

内容的注意

• このモデルでは乾季の流量はあまり現実的でない。 これは、流出が起こるのを「バケツがあふれる」場合に限ったためである。 実際には乾季でも基底流出があるので、 それを表わすプロセスをモデルに組みこむべきだが、次の段階で考えたい。
• チャオプラヤ川は灌漑が発達しているので、 自然プロセスだけを考えた流出の計算は観測値と大きく離れるかもしれない。

作業上の注意

モデルのプログラムがこのままだと、 実行条件を変えても、同じ日付の同じ変数のファイルには同じ名前が使われるので、 上書きになってしまう。 プログラムでデータを書き出す部分を修正してもよいが、 簡単なのは、ディレクトリ(フォルダ)を作って計算結果をそこに移すことである。 ただし、precont.awk, monav.awkなどを実行する際には とくに現在作業中のディレクトリ(current directory)に注意して作業すること。

参考

宮岡健吾氏が東京都立大学地理学科卒業研究で行なった内容を改訂したものが 松山ほか(1998), MIYAOKA et al. (1999)として発表されている。

この前に、Global Soil Wetness Projectの第1期として、 全球1度格子で、さまざまな陸面モデルに共通の気象データを与えて土壌水分の シミュレーションが行なわれていた。 対象期間は1987-1988の2年間。 まとめの報告はDIRMEYER et al. (1999)。 また日本気象学会誌(JMSJ)の特集号(77巻1B号)が組まれた。

宮岡氏はGSWPの結果のうち陸面モデルJMA-SiBで計算された各ます目からの流出高をもとに、 TRIPの格子状河道網を用い、流速は一定値を仮定して、河川流量を計算した。 ただし、流速の値はいろいろ試み、 流量を観測値と比較して季節変化の再現がよいものが適切と判断した。 このモデルの結果から、土壌水分・積雪・河川水を含めた流域水貯留量を求め、 観測データ(大気部分は再解析データ)から大気・流域水収支から求めた結果と 比較した。

文献

• Paul A. DIRMEYER, A.J. DOLMAN & Nobuo SATO, 1999: The pilot phase of the Global Soil Wetness Project. Bulletin of the American Meteorological Society, 80, 851 - 878.
• 池淵 周一, 椎葉 充晴, 宝 馨, 立川 康人, 2006: エース 水文学。 朝倉書店, 201 pp.
• 松山 洋, 宮岡 健吾, 沖 大幹, 1998: 大気・流域水収支法および河川流路網を用いたJMA-SiBの出力の検証 --北米Mackenzie川における事例解析。 気候システム変動の謎に挑む (東京大学気候システム研究センター), 93 - 105.
• K. MIYAOKA, H. MATSUYAMA & T. OKI, 1999: Validation of the output from JMA-SiB using the combined water balance method and a river routing scheme: A case study in the Mackenzie River basin. Journal of Geophysical Research, 104, 31199 - 31206.
• Taikan OKI and Y.C. SUD, 1998: Design of Total Runoff Integrating Pathways (TRIP) A global river channel network, Earth Interactions, Vol. 2, Paper 1. http://earthinteractions.org/からリンクされたところにある。 この論文は無料公開記事になった。
• P. XIE, A. YATAGAI, M. CHEN, T. HAYASAKA, Y. FUKUSHIMA, C. LIU, and S. YANG, 2007: A gauge-based analysis of daily precipitation over East Asia. Journal of Hydrometeorology, 8, 607-627.
• Kiyotoshi TAKAHASHI and Yasushi AGATA eds., 2007: GAME Phase 2 Collected Data. (GAME Data CD No. 12). http://game.suiri.tsukuba.ac.jp/cdroms/CD-ROM.html からCD-ROMイメージが配布されている。 ただし、CD-ROMで不足しているデータ提供もとの各国機関に関する情報については http://hydro.iis.u-tokyo.ac.jp/GAME-T/GAIN-T/routine/と http://hydro.iis.u-tokyo.ac.jp/GAME-T/GAIN-T/link.htmlを参照。