J. Jpn. Soc. Soil Phys.
土壌の物理性No. 153 p.33 ∼ 39 (2023)

汎用流域水循環シミュレーションツールにおける陸面から地下までの水 · 化学物質移動の数値的取扱いについて
田原康博 1

Numerical treatment of surface water, groundwater ﬂow and solute transport processes
in general-purpose, physically-based watershed modeling tools
Yasuhiro TAWARA1

1. はじめにやエネルギーの移動自体は考慮されない．陸面は地形 ·土地利用 · 土地被覆 · 植生の不均一分布をもつ 2 次元空間，地下は土壌や岩石の広がる 3 次元空間として表現さ陸面から地下までの水や化学物質の移動を定量的に捉れるのが一般的である．この技術を利用することによっえていくことは，流域における水循環の実態を理解してて，気象条件や土地 · 水利用等の様々な自然 · 人工環境いく上での助けとなり得るものである．しかしながら，流速が大きく異なる地表水と地下水を同時に取扱うことや，不飽和帯における毛細管現象など水と空気との相互の変化が流域内の物質 · エネルギー挙動にどのような影響を及ぼすかの定量評価を与えることができる．
菱谷 · 田原（2021）によれば，流域水循環シミュレー作用，さらに，降雨や気温，植生状態，地表水の変化なションツールは採用しているモデリング手法の違いにどによって時々刻々と変化する水分状態など，非常に複よって大きく 2 つのタイプに分類される．

雑で多くのプロセスを考慮していく必要がある．

1 つ目は，流域における水の流れを，河川などの水路陸面から地下までの水や化学物質の移動を定量的に捉

流れや斜面などの地表流れ，地下の不飽和流れ，地下飽えていくための技術として，流域水循環シミュレーショ

和流れのように複数のコンポーネントに分割して，コンンが挙げられる．流域水循環シミュレーションを可能とポーネント間でデータをやり取りしながら独立に計算するツールは近年ではいくつか存在し，中にはパッケーを行っていくもので，GSFLOW（例えば，Markstrom et

ジ化 · 汎用化され，高度なものであるにも拘わらず非常に利用しやすいものも普及するようになってきている．
ここでは，流域水循環シミュレーションにおいて対象al., 2008）や MIKE SHE（例えば，DHI, 2017）などが代表的なものとして挙げられる（菱谷 · 田原, 2021）．この手法では，コンポーネント間でデータをやり取りする際とする場である，陸面から地下までの水や化学物質の移に必要となる経験的パラメータの評価が極めて重要とな動に対する数値的な取扱いについて解説する．

る．また，コンポーネント毎に独立に計算していくため，次に示す 2 つ目のモデリング手法と比較して，多くの場

2. 流域水循環シミュレーション

合，計算時間が短くなる．しかしながら，解析者が対象とする流域に対してコンポーネントを決定し，コンポー

流域水循環シミュレーションは，陸域における様々なネント間のデータのやり取りを設計していく必要がある物質やエネルギーの輸送過程を計算対象とし，地表面とため，コンポーネントの選び方や，コンポーネント間のいう陸面と地下間のインターフェースで分断することつながりが実際の流域における水などの流れに対応付いなく，流域内の一体的な輸送過程を捉えようとするアプていないと現実とは大きく異なる結果となる場合がある

ローチである．流域水循環シミュレーションの対象とすことに留意する必要がある．

る場は，大気，海洋，陸面，地下に大別され，その中で2 つ目は，地表面から地下までを含めた流域全体を多は，流体（主に水，空気），エネルギー（熱），化学物質，数の四面体や六面体などによって直接表現し，流域にお土砂が考慮される．大気，海洋は降水や気温等の気象変

ける水の流れを 1 つのモデル（模型）の中で計算して化や潮位変化を反映させた境界場として取扱われ，物質いくもので，ParFlow（例えば，Maxwell et al., 2017）やHydroGeoSphere（例えば，Aquanty Inc., 2015），GETFLOWS（例えば，Tosaka et al., 2010; Mori et al., 2015）などが該当すると考えられる（菱谷 · 田原, 2021）．この1 Geosphere Environmental Technology Corp. 2-1, Kanda-Awajicho,
Chiyoda-ku, Tokyo 101-0063, Japan. Corresponding Author: 田原康博,地圏環境テクノロジー．
2022 年 12 月 5 日受稿 2023 年 2 月 3 日受理

土壌の物理性

34

第 153 号（2023）

Fig. 1 GETFLOWS が取扱うプロセスの概念図（所ら（2022）より引用）．

手法では，対象とする領域を丸ごとシミュレーションす

熱の移流，熱伝導，熱の発生，相変化 · 相間物質移動がるため，水などの流れはシミュレーションの結果として主なプロセスとして考慮される．これらの取扱い方法に現れてくる．反面，地下水とは流動する速度が大きく異ついても Table 2 に示した．

なる地表水流動を同時に取扱い，加えて，非線形性の強3.2 支配方程式

い地下不飽和帯も同時に計算するため，本質的に計算が水，空気と水に溶けて移動する化学物質に焦点を当て複雑で，計算時間が長くなる傾向にある．
以降では，2 つ目のモデリング手法を採用している汎用
流域水循環シミュレーションツールである GETFLOWSた場合，等温状態における水 · 空気 2 相 3 成分流体系の支配方程式は，以下に示す水相，空気相，水相中の溶存物質の質量収支式により記述される．

を例に，考慮するプロセスの数値的な取扱いの詳細について示す．

−∇ · (ρw υw ) + ρw qw =∂(ρw ϕ Sw )∂t(1)−∇ · (ρa υa ) + ρa qa =∂(ρa ϕ Sa )∂t

(2)

3. 流域水循環シミュレーションツールの数値的取扱い（ GETFLOWS を例に）
3.1 考慮するプロセス
GETFLOWS が取扱うプロセスの概念図を Fig. 1 に示す．

−∇ · (ρw υwCw ) + ∇ · Dw ∇(ρwCw ) + ρw qwCw=∂(ρw ϕ SwCw )∂t(3)

陸面で考慮する主なプロセスとしては，降水，積雪 ·融雪，蒸発散，地表水流動，浮遊砂輸送，雨滴衝撃によここに，∇は微分作用素，ρp は流体相 p の密度（kg m−3 ），る土砂生産，水流による浸食 · 堆積，地形変化，河川取の相互干渉，土壌蒸発，固相圧縮，物質の移流，水力学

vp は流体相 p の流速（m s−1），Cw は水相中の溶存物質の濃度（m3 m−3 ），Dw は水相中の溶存物質の水力学的分散
2−1係数（m s ），qp は流体相 p の生成 · 消滅量（m3 m−3−1s ），ϕ は間隙率（–），Sp は流体相 p の飽和率（–），t は時間（s），下付き文字 p は水相（w），空気相（a）を示す．

的分散，固相への吸着 · 脱着，崩壊 · 分解，物質の生成，陸面と地下，陸面–地下間で流体相 p の流速 vp は以下

排水，農地かんがい，ダム操作，雨水排水が挙げられる．
これらの取扱い方法を Table 1 にまとめた．
加えて，地下では，地下流体流動，毛細管効果，流体解説：汎用流域水循環シミュレーションツールにおける水 · 化学物質移動の数値的取扱いについて

35

Table 1 GETFLOWS が陸面で考慮する主なプロセスと取扱い方法．
プロセス

取扱い方法

降水

地表面へのソース項として考慮

積雪 · 融雪

積雪 · 融雪タンクが設けられ，降水の一時貯留（積雪），融雪モデルにより融雪量を推定しソース項として考慮

蒸発散

蒸発散モデル（単層熱収支法，ハーモン法などの経験法）により蒸発散量を推定し，シンク項として考慮

地表水流動

浅水波近似式に線形化拡散波近似を適用

浮遊砂輸送

乱流状態の平均的な振る舞いを記述した移流拡散方程式に従う

雨滴衝撃による土砂生産

直達雨量による生産と樹冠からの雨滴による生産，なお，地表水の水深による生産の減衰を考慮

水流による浸食 · 堆積

水中の浮遊砂濃度と流速によって浸食量 · 堆積量を計算

地形変化

土砂生産，浸食 · 堆積に伴う標高変化を考慮

河川取排水

河川におけるシンク · ソース項として考慮

農地かんがい農地へのソース項として考慮

ダム操作ダムにおける一時貯留や放流をシンク · ソース項として考慮

雨水排水

降雨量から差し引くことで考慮

Table 2 GETFLOWS が地下で考慮する主なプロセスと取扱い方法．
プロセス

取扱い方法

地下流体流動一般化ダルシー則に従う多相流れ

毛細管効果

毛細管圧力を水相飽和率の関数によって設定

流体の相互干渉相対浸透率を水相飽和率の関数によって設定

土壌蒸発蒸発散モデルから推定された土壌蒸発量に対して土壌蒸発効率を掛けることでシンク項として考慮

固相圧縮

流体圧力の線形関数に従う

物質の移流

一般化ダルシー則に従う地下流体流動による移動

水力学的分散分子拡散（フィックの法則に従う）と機械的分散固相への吸着 · 脱着
崩壊 · 分解

Freundlich 型 ·Langmuir 型の吸着モデル（Freundlich, 1907; Langmuir, 1918）によってシンク · ソース項として考慮（平衡，非平衡）
一次反応則に従い，シンク · ソース項として考慮

物質の生成

ソース項として考慮

熱の移流

一般化ダルシー則に従う地下流体流動による移動熱伝導

フーリエ則に従う

熱の発生ソース項として考慮

相変化 · 相間物質移動

Higbie の浸透理論（Higbie, 1935）によりモデル化，l は地
ここに，R は径深（m），nl は粗度係数（m−1/3 s）のように整理される．

表水流動方向の距離（m），h は水深（m），z は基準面か

陸面：

s2R3
υw = −nlυa = −

らの標高（m），K は絶対浸透率（m2 ），krp は流体相 p の∂h ∂z∂h ∂z+sgn+∂l ∂l∂l ∂l(4)Kkra∇ · (Pa + ρa gz)µa(5)相対浸透率（–），µp は流体相 p の粘性係数（Pa·s），Pp は流体相 p の圧力（Pa），g は重力加速度（m s−2 ）を示す．
なお，式（4）の sgn は符号関数であり，（ ∂∂ hl + ∂∂ zl ）の値が正の場合は 1，0 の場合は 0，負の場合は −1 となる．
陸面における地表水流動の流速 vw は，式（4）に示す
マニング則を仮定した開水路流れの浅水波近似式から拡散波近似したものによって評価する．一方，陸面におけ

地下，陸面–地下間：

υp = −Kkrp∇ · (Pp + gµp

Zる空気の流速 va については，式（5）に示す一般化ダルρp dz)(6)

シー則によって評価し，式中の絶対浸透率に 1.0 × 1030

土壌の物理性36

第 153 号（2023）

Fig. 2 GETFLOWS の三次元格子モデルの概念（Mori et al.,（2017）より引用し一部改変）．

m2 などの非常に大きな値を設定することで，地表水流子と呼ばれる六面体形状を用いることができ，それぞれ動に影響を及ぼさないものとして取扱う．また，地下なの格子体積および隣接格子間の絶対浸透率を正確に評価

らびに陸面–地下間の流体相 p の流速 vp については，式する．本手法は，格子毎の厳密な質量収支を保存するこ（6）に示す一般化ダルシー則によって評価する．とができる．

以下のパラメータは等温状態を仮定して以下のような関数として取扱う．

流体の圧縮性や，不飽和帯における相対浸透率，毛細管圧力などによる強い非線形性を有する三次元流体問
題を解く必要性から，時間離散化として，支配方程式

Pc = Pc [Sw ] = Pa − Pw

(7)を一次風上法によって完全陰的に有限差分展開を行い，Newton-Raphson 法によって反復的に解く手法を採用し

ϕ = ϕ [Pa ] = ϕ0 (1 +Cr (Pa − P0 ))(8)

krp = krp [Sw ](9)ている．
これらの離散化によって 7 重対角行列を係数行列とした連立方程式を得ることができ，行列の各成分は，空気相圧力 Pa ，水相飽和率 Sw ，水相中の溶存物質の濃度 Cw
を未知数とした 3 × 3 小行列となり，システム全体の自 ρp = ρp Pp

(10)

 
µp = µp Pp(11)

由度は格子数×未知数となる．これを，前処理付き Gbi-CGStab（s, L）法（Tanio and Sugihara, 2010）などの線形ソルバーによって解く．前処理には Nested Factorization（Appleyard et al., 1981）などを適用している．また，大容量計算を高速処理するため，非線形反復過程の中で，設定した許容誤差の数オーダー小さい値以下に収斂しここに，Pc は水–空気系の毛細管圧力（Pa），Cr は固相圧た格子を係数行列から自動的に除外する逐次陽化処理

縮率（Pa−1 ），ϕ0 は標準状態における間隙率（–），P0 は（Successive Locking Process）（Tosaka et al., 1987）や領標準状態における圧力（Pa）を示す．式中の krp [Sw ] な域分割法（Domain Decomposition Method）（Toselli and

どは krp が Sw の関数であることを意味する．なお，支配Widlund, 2005）による高速並列計算を採用し，大規模三方程式中の主変数（未知変数）は空気相圧力 Pa ，水相飽次元問題をより効率的に解くことができる．

和率 Sw ，水相中の溶存物質の濃度 Cw の 3 つである．

三次元格子モデルの概念を Fig. 2 に示す．格子システ

3.3 数値解法

ムは大別して役割の異なる 3 つの層（大気層，地表層，空間離散化は，複雑な地形起伏や地盤物性分布を三次地下地盤層）から構成される．各層の概要を以下に示す．

元格子によって表現可能な積分型有限差分法（Integral

地上の大気層は格子モデルの第一層として表現されFinite Difference Method: IFDM）による．構造格子を前る．ここは，物理的には毛細管圧力は 0（多孔質媒体でな提としており，格子形状にはコーナーポイント型差分格い空間）で，非常に大きな容量をもつ大気空間であるこ解説：汎用流域水循環シミュレーションツールにおける水 · 化学物質移動の数値的取扱いについて

37

Fig. 3 地表面における浸透 · 湧出の計算の模式図（菱谷 · 田原（2021）より引用し一部改変）．とから間隙率は 1.0 × 1030 のような極端に大きな値を設の数値的取扱いについて説明する．

定する．絶対浸透率については 1.0 × 1030 m2 など非常

Fig. 3 に示すように，地表面に高さ H の地上（地表層）に大きな値とし，水相飽和率は 0 として設定する．容量格子 1 があり，隣接する地下最上層（地下地盤層）格子の大きい大気空間では圧力や飽和率などの状態量の変化

2 との間の流動量を考える．ここでは簡単のために地表面に傾斜は無く，水相の密度 ρw や粘性係数 µw は格子間で同じとする．地上格子 1 に水深 h の状態を設定した
とき，地表面には水圧 Pa1 +ρw gh が与えられ，地下最上層格子 2 との水理ポテンシャル差で浸透率 K2 の場を浸透または湧出することになるため，その流動量 Qz （m3s−1 ）は以下の式で表現される．

が生じない．大気層格子と地表層格子の境界では式（6）を適用し，圧力勾配に従う空気の流動のみが計算される．
地表層格子では河川や斜面を流れる地表水，湖沼 · 海
洋の停留水を表現する．地表層格子における水平方向の水と空気の流動は，それぞれ，式（4），（5）によって評価される．前述したように，空気については絶対浸透率に 1.0 × 1030 m2 などの非常に大きな値を設定することで，水の流動に影響を及ぼさないものとして取扱う．場Qz

所による地表水の移動し易さの相違は土地利用や被覆状=−態などに対応付けた粗度係数として考慮され，個々の格子に与えられる．また，地表層格子–地下地盤層格子間の移動に対しては，式（6）で評価され，水については後述する疑似的な毛細管圧力（Tosaka et al., 2000）を考慮することによって浸透 · 湧出が整合的に表現される．なお，間隙率は多孔質媒体（固相）がないことから 1.0 を設定する．
地下地盤層を構成する格子は，式（6）に示した一般化ダルシー則に従った多相多成分流体流動を表現する．

K2 A12 Pa1 + ρw gh + ρw g(z1 − dz1 ) − (Pw2 + ρw gz2 )µw
dz2
K2 A12 Pa1 + ρw gh − ρw gdz1 + ρw g(z1 − z2 ) − Pw2=−µwdz2(12)
ここに，A12 は格子間面積（m2 ）を示す．その他の記号の説明については Fig. 3 を参照されたい．
ところで，Fig. 3 に示す格子間流動量を地下での取扱い（一般化ダルシー則）と同形式で記述すると次式のようになる．

個々の格子には地層の水理物性（間隙率，絶対浸透率，毛細管圧力，相対浸透率，固相圧縮率）を与える．なお，Qz

地下における飽和–不飽和浸透流シミュレーションと比K12 A12 Pw1 + ρw gz1 − (Pw2 + ρw gz2 )µw
z1 − z2
K12 A12 Pa1 − Pc1 + ρw g(z1 − z2 ) − Pw2=−µw
dz1 + dz2

較した場合，空気の流動を解いているか否かが大きく異なる点である．

3.4 地表面での浸透 · 湧出の取扱い=−(13)

ここでは，流域水循環における水や化学物質の移動を定量的に捉えていくために重要となる要素の 1 つである

K12 は調和平均によって求められる格子間の平均浸透率

地表面での浸透および湧出に着目して，GETFLOWS で（m2 ）である．地表層格子には，絶対浸透率が非常に大土壌の物理性

38

第 153 号（2023）

Fig. 4 GETFLOWS よって計算された 3 次元的な水の流れの可視化（所ら（2022）より引用）．

きい空間と捉え，K1 に 1.0 × 1030 m2 などの極端に大き間 · 空間スケールの幅の大きさから，土壌物理学の分野な値を設定するため，以下に示す式が導かれる．

だけでなく，様々な分野の問題に適用可能なものであると考えられる．今後は，化学反応モデルや生態系モデル

K2 A12 Pa1 − Pc1 + ρw g(z1 − z2 ）− Pw2
Qz = −µw
dz2など異分野モデルとの連成 · 連携によってさらにその適(14)地上格子 1 の毛細管圧力 Pc1 を ρw gH(0.5 − Sw1 ) に置き換えると，式（12）と式（14）は同一のものとなる．式（13）は地下における格子間流動量を計算するための一般的なダルシー則であり，地上格子の絶対浸透率に非常に大きな値を，疑似的な毛細管圧力として前記の式を与えることによって地表面における流動量が地下の流れと全く同じ手続きで取扱うことができるようになる．

Fig. 4 に GETFLOWS によって計算された水の流れのシミュレーション結果を例示した．これは，地表に降り注いだ雨を粒子に見立て，その流動軌跡を三次元的に表示したもので，白い粒子とラインが地上の水を，黒が地用範囲の拡大が期待される．一方で，考慮しているプロセスの数値的な取扱いを正しく理解した上で，解析者が対象とする問題に対して適用して良いのかを都度判断していくことが必要である．また，出力された結果が現実にはあり得ない状態を示していないかを良く確認することも同時に重要である．
常日頃から研究や実務などでシミュレーションツールを用いて様々な土壌物理学の問題に向き合っておられる方，これから初めて使用していきたいと考えておられる方に対しても，本稿が研究の進展や速やかな実務進行の助けになればうれしく思う．

引用文献

下水を表している．この図に示したように，流域内の水の流れがシームレスに表現されており，流域水循環におAppleyard J.R., Cheshire I.M. and Pollard R.K. (1981): Spe-ける重要な要素の 1 つである浸透と湧出を整合的に取扱

cial techniques for fully-implicit simulators. Proceedings

うことができていることを示している．

of the European Symposium on Enhanced Oil Recovery,
pp. 395–408. Bournemouth, England.

4. おわりにAquanty Inc. (2015): HydroGeoSphere User Manual.
DHI (2017): MIKE SHE Volume 1: User Guide.
Freundlich H. (1907): Über die Adsorption in Lösungen.

汎用流域水循環シミュレーションツールにおいて，陸面から地下までの水や化学物質の移動に着目し，考慮するプロセスが数値的にどのように取扱われているかについて，GETFLOWS を例に解説した．
流域水循環シミュレーションツールは，対象とする時

Zeitschrift für Physikalische Chemie – Stöchiometrie undVerwandschaftslehre, 57(4): 385–470.
Higbie R. (1935): The rate of absorption of a pure gas into astill liquid during short periods of exposure. Transactionsof the American Institute of Chemical Engineers, 31: 365–解説：汎用流域水循環シミュレーションツールにおける水 · 化学物質移動の数値的取扱いについて377.

and transport processes derived from the Fukushima ac-

菱 谷 智 幸, 田 原 康 博 (2021): 地 下 水 流 動 シ ミ ュ レ ー ションツールにおける降雨などの地表面からの浸透に関 す る 取 扱 い に つ い て. 地下 水 学会 誌, 63(2): 57–66,

doi:10.5917/jagh.63.57.

cidents. Environmental Modelling & Software, 72: 126–146, doi:10.1016/j.envsoft.2015.06.012.
Tanio M. and Sugihara M. (2010): GBi-CGSTAB (s, L):
IDR(s) with higher-order stabilization polynomials. Jour-Langmuir I. (1918): The adsorption of gases on planesurfaces of glass, mica, and platinum. Journal of
the American Chemical Society, 40(9):391361–1403,

nal of Computational and Applied Mathematics, 235(3):765–784, doi:10.1016/j.cam.2010.07.003.
所眞理雄, 高橋桂子, 石川幹子 (2022): 水大循環と暮らしⅢ
持続可能な水大循環を実現する, pp. 152–155. 丸善プラ

doi:10.1021/ja02242a004.
Markstrom, S.L., Niswonger, R.G., Regan, R.S., Prudic, D.E.

ネット, 東京.

and Barlow, P.M. (2008): GSFLOW — Coupled ground-Tosaka H. and Matsumoto Y. (1987): An efﬁcient reservoir

water and surface-water ﬂow model based on the Integra-

simulation by the successive explicitization process. Jour-tion of the precipitation-runoff modeling system (PRMS)

nal of the Japanese Association for Petroleum Technology,

and the modular ground-water ﬂow model (MODFLOW-

52 (4): 307–313, doi:10.3720/japt.52.307.

2005). U.S. Geological Survey Techniques and Methods
6–D1.

Tosaka H., Itoh K. and Furuno T. (2000): Fully coupled formulation of surface ﬂow with 2-phase sub-Maxwell, R.M., Kollet, S.J., Smith, S.G., Woodward, C.S. ,surface ﬂow for hydrological simulation. Hydrologi-Falgout, R.D., Ferguson, I.M., Engdahl, N., Condon, L.E.,

cal Processes, 14(3): 449–464, doi:10.1002/(SICI)1099-Hector, B., Lopez, S.R., Gilbert, J., Bearup, L., Jeffer-1085(20000228)14:3<449::AID-HYP948>3.0.CO;2-9.

son, J., Collins, C., de Graaf, I., Prubilick, C. Baldwin, C.,

Tosaka H., Mori K., Tada K., Tawara Y. and Yamashita K.

Bosl, W.J., Hornung, R. and Ashby, S. (2017): PARFLOW(2010): A general-purpose terrestrial ﬂuids/heat ﬂow sim-user’s manual. Integrated Ground-Water Modeling Center

ulator for watershed system management. Proceedings ofReport GWMI 2016-01.

IAHR International Groundwater Symposium 2010, Va-Mori K., Tada K., Tawara Y., Ohno K., Asami M., Kosaka

lencia, Spain.

K. and Tosaka H. (2015): Integrated watershed mod-Toselli A. and Widlund O.B. (2005): Domain decomposi-eling for simulation of spatiotemporal redistribution oftion methods — algorithms and theory. Springer Series in
Computational Mathematics.

post-fallout radionuclides: Application in radiocesium fate

要旨

陸面から地下までの水や化学物質の移動を定量的に捉えていくことは，流域における水循環の実態を理解していく上での助けとなり得るものであるが，陸面や地下で起こる複雑で様々なプロセスを同時に考慮していく必要がある．本稿では，陸面から地下までの水や化学物質の移動に焦点を当て，汎用流域水循環シミュレーションツールにおいてどのように数値的に取扱われているかについて，GETFLOWS を例に解説した．
キーワード：流域水循環シミュレーション，陸面，地下，浸透 · 湧出，GETFLOWS