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

地下水涵養過程の数値解析再訪
斎藤広隆 1 · 取出伸夫 2

Numerical analysis of groundwater recharge revisited
Hirotaka SAITO1 and Nobuo TORIDE2

1. はじめに
降雨後の不飽和帯の下方浸透による水分移動は，地下水面に到達すると地下水涵養流となるが，定量的な測定には一定の困難さが伴い，モデルを用いた解析が有用である．塩沢 · 中野（1981）は，降雨による地下水涵養過程では圧力勾配が小さくなり重力が支配的となるため，
重力流れのみを考慮した伝達モデルを用い地下水涵養過程を解析し，その有用性を示した．しかし自然条件下の土は，植物根を含む土中生物由来の粗間隙や亀裂などの構造に由来する間隙に富んでおり，土中水はこれら間隙構造を水みちとして流れ，局所的な水の流れが発生することは古くから知られている（Jury and Horton, 2006）．
多孔質体中を進む流れをマトリックス流と呼ぶのに対して，このような局所的な流れを総称して選択流（または Preferential Flow）とよび，マトリックス流と比べてリチャーズ式に基づく汎用一次元不飽和土中水分 · 溶
質移動モデル HYDRUS を用いた計算結果と比較する．

HYDRUS のような汎用モデルの多くは，熱移動や溶質移動との同時解析も可能で，地下水涵養過程のような水文物理現象を理解する上での有用性に加えて，その応用の広さについても考察する．本報ではその一例として，
浸透する降雨をトレーサーとみなした移流分散式による溶質移動解析により，土に侵入した新しい水と土層内の古い水を区別した解析が可能であることを示す．また，選択流の評価には， HYDRUS に実装されている DP モ
デルを用いて，降雨浸透下で選択流が水分と溶質の地下
水涵養過程に与える影響のモデル化について計算例を示
す．なお，本報のすべての解析は，オープンソースバー
ジョンである汎用一次元飽和 · 不飽和水分移動解析プロ
グラム HYDRUS-1D ver4.17（Šimůnek et al., 2013）を用いた．

水分移動速度が有意に速い．地下水涵養過程において，
涵養量の半分以上が選択流によることもあるともいわ

2. 地下水涵養流に重力の果たす役割：

れ（Huang et al., 2019 など），また選択流による水分移伝達モデル

動に伴い汚染物質は化学的 · 生物学的な反応性の高い
表層をバイパスし地下水に到達することが知られてい塩沢 · 中野（1981）は，地下水涵養流れにおいて地下る（Allaire et al, 2009 など）．そのほかにも封入空気に水が十分に深いとき土中水圧力水頭勾配が水分移動に与よるマトリックス流の阻害（Wang et al., 1998 など）なえる影響は重力と比較して十分小さく，重力のみを考慮

ど，地下水涵養過程は複雑な現象となる．しかし，複雑した伝達モデルを使うことで，近似的な扱いが可能であな構造をもたない多孔質を仮定し，飽和 · 不飽和土中水るとした．この伝達モデルを使うことで，地表面水分フ

分移動の基礎方程式であるリチャーズ式にもとづいて，ラックスの変化が地下水涵養流れとしての土中水フラッ

降雨浸透による地下水涵養過程を計算することは，地下クスの変化を予測することができ，非定常不飽和浸透流

水涵養フラックスの動態を不飽和帯の水文物理現象としとしての地下水涵養流れを特徴付けることができる．こて理解する上で重要である．加えて，様々な土中水分移

こではまず，塩沢 · 中野（1981）による伝達モデルを用い動に関するモデルの開発が進み，選択流の影響を考慮すて，そこからわかる地下水涵養流れの特徴について示す．

るため，土を水分フラックスが遅い領域と速い領域に分

まず，リチャーズ式による降雨による浸潤と排水の地下割して，それぞれの領域にリチャーズ式を適用する Dual

水涵養過程の計算結果に伝達モデルを適用し，重力流れ

Permeability（DP）モデルによる補正も可能である．

の役割とモデルの有用性について論じる．なお，計算に本報では，降雨による地下水涵養過程について，まず用いた土性は，シルト（Silt）および砂質ローム（Sandy

塩沢 · 中野（1981）の先駆的な伝達モデルを振り返り，

Loam）である．Fig. 1 に van Genuchgten – Mualem モデル（VGM モデル）によるそれぞれの水分保持曲線 θ (h)および不飽和透水係数 K(h) を示す．ここで θ は体積含水率，h は土中水圧力水頭，K は不飽和透水係数である．

1 東京農工大学大学院農学研究院．
2 三重大学大学院生物資源学研究科．

2023 年 2 月 2 日受稿

2023 年 2 月 28 日受理

土壌の物理性

14

第 153 号（2023）

Fig. 1 シルト，砂質ローム，および仮想的な亀裂の（a）水分保持曲線 θ (h) と（b）不飽和透水係数 K(h)．

2.1 伝達モデル

力勾配は小さくなり負となる（深くなるにつれて圧力がここでは，下方浸透流（基底流）が生じている土層に大きくなる）．この負の勾配は正の勾配を持つ重力と打対して降雨により浸潤前線が地下水に到達し，定常流れ

ち消しあって小さな全水頭勾配を作り，透水係数が大きとなった後に降雨が止み，再び基底流による定常流れにくなっても一定の下方水分フラックスを実現する．

戻る排水過程を考える．基底流のフラックスは，東京の降雨が止み土層からの排水が始まると，土中水圧力

年間降水量が約 1500 mm で，可能蒸発量が 800 mm 程が初期値より低下する排水前線が上端より形成され下度であり（新井, 1989），仮に余剰水がすべて浸透すると方へ移動し，それに伴って土層全体で土中水圧力が低約 0.2 cm d−1 の下方浸透流となることを考慮して，0.1

下する．L = 500 cm のとき，上端の土中水圧力水頭がcm d−1 とした．Fig. 2（a）–（c）は，シルトを対象として，地下水面の深さ L が 500 cm のとき，下向き 0.1 cm
d−1 の定常基底水分フラックスが生じている土層に，計算開始時間（t = 0 d）に 1.0 cm d−1 の降雨が生じ，そq2 = −K(h2 ) に対応する h = h2 となる前に排水前線は地れが計算期間を通して継続して起こるときに，浸潤前線

布をみると，地下水面近傍では地下水面の影響が少ない

が下端の地下水に到達して定常水分流れとなるまでの土深さと比較して，体積含水率の変化は地下水面に向かっ中水圧力水頭 h，体積含水率 θ ，水分フラックス q 分布て小さくなる．一方で Fig. 2（f）に示す水分フラックスの計算結果である（今後座標は地表面を z = 0 cm とし，分布からは，浸潤過程では地下水面近傍でも水分フラッ

フラックスは上向きを正として設定する）．また Fig. 2

クス変化は大きく，排水過程では地下水面に向かって水（d）–（e）は，この後 1.0 cm d−1 の降雨が止んで，基底流

分フラックスの変化が大きくなる．このことを説明する

である 0.1 cm d−1 に戻るまでの排水過程における土中水

ために，塩沢 · 中野（1981）は連続式を変形し得た次式

圧力水頭，体積含水率，水分フラックス分布である．基を用いた（詳細な導出は付録に示す）．ここで，下付き文

底流や降雨による定常水分流れが生じる条件では，地下字は固定する変数を表し，z は任意のフラックス q を与水の影響が少ない土層上部では均一な土中水圧力分布をえる座標を表す．

下水面に到達し，土層上部より徐々に土槽全体の圧力水
頭が h2 に近づく．このとき，土中水圧力の勾配は全体に緩やかなものとなる． Fig. 2（e）に示す体積含水率分

示し，圧力勾配がゼロとなる重力流れが生じる．このと



∂z
−
∂t

きの土中水圧力は q = −K(h) となる h の値を示し，シルトにおいて大きさ 1.0 cm d−1 の降雨では，土中水圧力は




=
q

∂q
∂θ


(1)
z

−16.5 cm となった．一方，地下水面の影響が及ぶ土層下部では，土中水圧力は地下水面に向かって増加し，透この式は，任意の水分フラックス q となる座標 z の移動水係数が大きくなる．そのため，定常水分流れが続いて速度が，その地点の体積含水率 θ に対する水分フラック

いるときは圧力勾配が重力を打ち消す方向に形成され，

ス q の変化率に等しいことを示している．Fig. 2（b）（e）

圧力分布は水分移動が静止する平衡分布に漸近する．

に示すように座標 z が地下水面に近づくにしたがって体Fig. 2（a）–（c）より，浸潤過程では初期の土中水圧

積含水率 θ の変化が小さくなるため，地下水面に向かっ力（および体積含水率）分布が一定の深さにおいて，明て降雨の水分フラックス q の伝達が速くなることを示し

瞭な伝達領域と浸潤前線が形成され（取出ら, 2009），そている．

の形が維持されながら下方へ移動する．地下水の影響がまた，塩沢 · 中野（1981）は，Fig. 2 の排水過程では大きくなる深さにおいては，徐々に浸潤前線先端部の圧圧力勾配が小さく，水分動態は重力による要素が重要で

解説：地下水涵養過程の数値解析再訪

15Fig. 2 地下水面が深さ L = 500 cm（z = −500 cm）に位置するときの，HYDRUS を用いた数値計算によるシルトの浸潤過程（境界で下向き 1.0 cm d−1 ）とその後の排水過程（境界で下向き 0.1 cm d−1 ）における（a）（e）土中水圧力水頭分布，（b）（d）体積含水率分布，（c）（f）水分フラックス分布．

あるとして，圧力勾配による水分移動を無視した伝達モ
水分移動の基礎方程式であるリチャーズ式で重力項のみ

∂z
dK
=−
∂t
dθ

を考慮した伝達モデルでは，不飽和土中水分移動基礎方式 (3) はある θ を与える z の降下が，dK/dθ （> 0）でデルを使って，排水過程の水分動態を解析した．不飽和与えられることを示しており，dK/dθ は t によらず一程式としては次式があたえられる．

∂θ∂ K(h) ∂ θ dK
=
=
∂t
∂z
∂ z dθ

(3)

定となるため，この降下速度は一定となることを示して(2)

ここでは，θ は z と t の関数であるが，z を θ と t の関数

いる．また，K(θ ) は通常下に凸な単調増加関数であり（d2 K/dθ 2 > 0），θ が大きくなると dK/dθ は大きくなる．伝達モデルから，θ が大きいところから排水が素早とすると，次式のように書き換えることができる（詳細く進行し，θ が小さいところはゆっくりと排水が進むとな導出は付録に示す）．

いう性質が示される．Fig. 2（e）に示す排水過程の土中

土壌の物理性

16

第 153 号（2023）はできないので，微分し，dK/dθ = L/t となる θ を数値的に求め，q(L,t) を得ることができる．ここで，Ks は飽
和透水係数，l および m は経験パラメータ，θs および θrはそれぞれ飽和および残留体積含水率である．このようにして求められる地下水涵養フラックス q(L,t) は深さ Lにのみ依存し，基本排水曲線とよばれる．なお，塩沢 · 中
野（1981）では地下水面が深さ L にある場合，その影響がおよぶ範囲では重力流れを仮定できないため，その影
響範囲 r を除いた L (= L − r) におけるフラックスで評価しており，本報ではシルト，砂質ロームともに r = 50

cm として解析を行った．
Fig. 3 に異なる地下水面深さに対して，シルトと砂質
ロームを対象として伝達モデルより求められた地下水面
深さ（z = −L）におけるフラックスとリチャーズ式に基づく数値解析によるフラックスを示す．数値解析による
排水フラックスと伝達モデルによるフラックスは土性やL によらずよく一致しており，排水過程，つまり地下水
Fig. 3 （a）シルトと（b）砂質ロームの排水過程における，
地下水面深さ（z = −L）における水分フラックスの数値解（実涵養過程においては重力による水分移動が支配的である

線）と伝達モデルによる予測（破線）．

とが確認できた．

ことについて，HYDRUS の計算結果でも再現できるこ

2.2 降雨による地下水涵養フラックス
水分分布の計算結果から，飽和に近い θ の場合は 5 d でほぼ排水が完了しているにもかかわらず，0.4 cm3 cm−3付近の θ の場合は 60 d でも地下水に到達しておらず，
上記の性質を確認することができる．
このように地下水面近傍と地下水面の影響のない位置
では，土中水の動態は大きく異なる．そこで，地下水面を無限に深いとし，地下水面の影響がない状態で，重力項のみを考慮した伝達モデルを使うことで，深さ L におけるフラックスを評価することができる．フラックス qとなる z が距離 L 移動するのに要する時間 t は，L を移動速度で割ってL(∂ z/∂ t)q

分フラックスに達することなく，降雨が止む．したがって水分供給としては，一定量の水分が有限時間に与えられるパルス入力となる．このとき，降雨浸透過程では浸潤前線が形成され浸潤は進行するが，上部からの水分供給停止後は，浸潤前線上部では排水過程となる．Fig. 4（a）に，地下水面が深さ 500 cm に位置するシルトと砂質ロームの土層に対して豪雨を想定して 10 cm の降雨を10 cm d−1 の強度（つまり 1 日かけて 10 cm の降雨量）で与えたときの土中水圧力変化を示す．降雨浸透過程でするものの，土中水圧力が急激に変化する浸潤前線が下(4)で与えられる．式（3）および（4）より L/t = dK/dθ となる θ を得ることで，深さ L，時間 t における水分フラッ
クスの大きさが |q(L,t)| = K(θ ) より得られる．塩沢 · 中
野（1981）では，パラメータ α と p を用いて，不飽和透水係数を K(θ ) = αθ p と単純な体積含水率の指数関数で与えることで，伝達モデルを解析的に評価した．本解析に対して伝達モデルを検証するためには，次式の VGM
モデルの体積含水率による微分が必要である．

K (θ ) = Ks

うに降雨時に浸潤前線が地下水面に到達し，一定の水は水分供給が持続されないためピーク圧力は徐々に低下t=

実際の降雨下の地下水涵養過程では，前節の例題のよ" !m #2
θ − θr lθ − θr 1/m
1− 1−
θs − θrθs − θr(5)方へ移動する．一方で，浸潤前線より上部では排水により全体として緩やかな土中水圧力分布が形成される．このように，パルス入力下の水分移動でも前節に示した排水過程が進行し，この過程は土中水の「再分布」として知られている．

Fig. 5 にシルトについて様々な降雨強度 i に対して，その強度が 1 日続くと仮定し，地下水面深さが 500 cmと 200 cm のときの，地下水面深さ（z = L）における排水フラックスの時間変化および 2.2 節で求めた伝達モデルに基づく排水フラックス q(L,t) を示す．浸潤前線が地下水面深さに到達すると水分フラックスが大きくなり，あるピークに達したのち排水過程となりフラックスは時間とともに小さくなる．フラックスの立ち上がり時間は
降雨強度が小さくなるにつれて遅くなり，ピークフラッ
クス値は小さくなる．一方で排水過程では，地下水面深さや降雨強度 i によらず，伝達モデルによる基本排水曲

しかし，VGM の場合 dK/dθ = L/t を解析的に解くこと線に漸近し，パルス降雨においても排水過程の水分移動解説：地下水涵養過程の数値解析再訪

17Fig. 4 シルトと砂質ロームの均一土層に 10 cm d−1 の強度（i = −10 cm d−1 ）で 10 cm の降雨を与えたときの，（a）土中
水圧力水頭分布と（b）地下水面深さ（z = −500 cm）における水分フラックスの時間変化．

Fig. 5 シルトにおいて，地下水面深さが（a）500 cm と（b）200 cm のときに，降雨強度 i が異なる条件のもと HYDRUS
で求めた地下水面フラックスの時間変化と伝達モデルによる予測．

は重力が支配的な流れとなる．

よるピークフラックスの時間が遅くなるが，これは数値

降雨による地下水涵養を評価する際に，水分フラック

解析では圧力勾配を駆動力として考慮するため，L が大スのピーク時間を知ることも重要である．伝達モデルできくなるにつれて圧力の拡散が起こりピーク到達時間が求まる基本排水曲線において，時間 t 以降の基本排水曲

遅れるためである．

線の積分が総降雨量と一致するような時間 t が得られ

2.3 古い水と新しい水：移流分散式による解析

ると，伝達モデルでは重力のみを考慮した水分移動であ

ここまで，降雨浸透による地下水涵養過程の特性は，り，その時間に地下水面においてピークフラックスが発

重力項のみを考慮した伝達モデルによって把握すること生する． Fig. 6 に，リチャーズ式に基づく数値解析によを示した．一方，リチャーズ式に基づく数値解析では地る降雨浸透過程の地下水涵養フラックスのピーク発生時

下水涵養過程の詳細をさらに知ることができる．古くか間と伝達モデルによるピークフラックス発生時間の比較

ら，地下水に到達する水が，土に侵入した「新しい水」でを示す．おおよそ 1 : 1 直線にのるものの，地下水面が深あるか，土層内に存在していた「古い水」であるかにつ

くなるにつれて（L が大きくなるにつれて）数値解析にいて議論されてきた．しかし，水分移動モデルで水分量

土壌の物理性

18

第 153 号（2023）の流れである選択流が発生する．この選択流による水分移動に伴う汚染物質は，特に地下水汚染評価の上で極めて重要な役割を示す．ここでは計算領域を二分割してリチャーズ式を適用する DP モデルをとりあげ，その考え方および計算例を紹介する．

3.1 Dual Permeability モデル
選択流を含む土中水分 · 溶質移動の解析に用いられるモデルの一つに DP モデルがある（Šimůnek and van

Genuchten, 2008）．DP モデルでは，Fig. 8 に示すように
土をフラックスの小さい遅い領域（マトリックス部分）とフラックスの大きい速い領域（亀裂部分）にわけ，それぞれの領域に対して水分移動にはリチャーズ式を，溶質移動には CDE を適用する．なお，本報では遅い領域
をマトリックス，速い領域を亀裂として表記するが，DPFig. 6 シルトにおいて HYDRUS と伝達モデルを用いて求めた降雨下の地下水涵養フラックスのピーク強度発生時間の比較．

モデルの本質は二領域に分割することであり，マトリッ
クスと亀裂の適用に限定されるものではない．DP モデ
ルでは亀裂とマトリックスの体積割合を w，1 − w としてリチャーズ式を適用する．また，領域間の水分交換はを求めることでは，この新しい水と古い水の区別はでき

圧力差に比例する吸込み · 湧き出し項として一次水分移ない．そこで，浸透する新しい水を土粒子と反応しない

動速度 Γw により与える．ここで，w = 0 の場合はマトトレーサーとみなし，移流分散式（CDE）を仮定した溶リックスのみの均一モデルとなる．

質移動問題として検討した．

分フラックスのピークより大幅に遅れ，水分フラックス

"
!#
∂ θf
∂
∂ hf
Γw
=
Kf (hf )
+1
−
∂t
∂z
∂z
w
"
!#
∂ θm
∂
∂ hf
Γw
=
Km (hm )
+1
+
∂t
∂z
∂z
1−w

がピークとなるフラックスにはトレーサーは一切含まれΓw = αw (hf − hm )ない．基底流が −0.1 cm d−1 のとき，平均間隙流速は約

θ = wθf + (1 − w) θm = θF + θM

Fig. 7 に，L = 500 cm で 10 cm の降雨を土粒子に吸着しないトレーサーとみなし，分散長を 10 cm とした CDEを適用したときの下端の水分フラックスとトレーサー濃度の時間変化と鉛直分布を示す．トレーサーの到達は水0.27 cm d−1 であり，500 cm の移動にはおよそ 1850 d(6)かかる．この時間は，CDE に基づく数値解析による溶

ここで，下付き文字 f と m はそれぞれ亀裂とマトリッ

質移動のピーク到達時間とほぼ一致する（Fig. 7）．CDE

クスを指し，αw は一次水分移動係数，θf と θm はそれによる溶質移動では，浸透したトレーサーは断面方向にぞれ亀裂とマトリックスにおける体積含水率，θF と θM

完全に混合しながら下方へ移動する（Jury and Horton,

はそれぞれ土全体積に対する亀裂とマトリックス領域の2006）．そのため，CDE に基づく地下水涵養過程は，浸体積含水率である．亀裂とマトリックス間の水分移動速

透する新しい水が古い水を押し出す現象であることを示

度 Γw は，二領域間の土中水圧力水頭の差を駆動力とすしている．一方，現実の地下水涵養過程は，局所的な流る一次式で与える．なお αw は [T−1 ] の次元を持ち，それが生じていて，断面方向の溶質の混合は不十分であり，の逆数 1/αw は亀裂とマトリックスの水分交換の時間ス

CDE の予測よりは速やかに地下水に到達すると考えら

ケールを与える．溶質移動ではそれぞれの領域に以下のれる．この溶質混合の程度は溶質の移動速度に影響する

CDE を適用する．ただし，ここでは吸着のないトレーので，特に汚染物質の移動を考える上で重要である．こサーを対象とする．

の点については，次節の DP モデルでさらに検討する．

3. 選択流の影響
2 章ではマトリックス流のみ考慮した均一土層にて，
地下水涵養過程について伝達モデルによる現象の解析!
∂ θf cf∂∂ cf∂ qf cf Γs
θf Df
=
−
−
∂t
∂z
∂z
∂z
w
!
Γs
∂ θm cm∂∂ cm∂ qm cm
=
θm Dm
−
+
∂t
∂z
∂z
∂z
1−w

(7)

とリチャーズ式に基づく数値解析の有用性について示した．一方で，自然条件下の土は，植物根を含む土中生

Γs は二領域間の溶質の交換速度を表し，水分交換と同様

物由来の粗間隙や亀裂などの間隙構造に富んでおり，土に濃度差に比例する交換項と水分交換による移流項の和

中水はこれら間隙構造を水みちとして流れ，局所的な水で与える．

解説：地下水涵養過程の数値解析再訪

19Fig. 7 シルトにおいて，地下水面深さが 500 cm のときに，100 mm の降雨を与えたときの（a）地下水面における水分フ
ラックスおよび降雨に含まれるトレーサー溶質フラックスの時間変化と（b）土中水圧力水頭分布とトレーサー濃度分布．

3.2 DP モデルの計算例
本来，亀裂内に水分流れが生じるのは強度の強い雨のときで，弱い雨の時は水分流れが生じない．そのため，解析の際に境界であたえる水分フラックスの大きさは，水分移動の計算に極めて重要な役割を果たす．しかし，
ここでは Fig 4 の例題と同じ境界条件とし，境界条件が亀裂内水分フラックスの違いに与える影響については検討していない． Fig. 9 にシルトを対象として Fig. 4と同じ境界条件で w = 0.05，ωdp = 0（有効拡散移動なし）のとき，αw を変化させ亀裂とマトリックスそれぞれの地下水面深さ（L = 200 cm）における水分フラックスを示す．ここで，亀裂の水分移動特性は，Fig. 1 に示す

θs = 0.8 で粗間隙が大半を占め保水力が弱く飽和領域での透水係数が大きい水分保持関数と不飽和透水係数を用

Fig. 8 Dual Permeability モデルの概念図．

いた．二領域間の水分移動を考慮しない場合（αw = 0），Γs = ωdp (1 − w) θm (cf − cm ) + Γw c∗

亀裂では浸潤前線の到達が速く，ピークフラックスも大(8)きくなる．一方マトリックスでは亀裂よりも大きく遅れて浸潤前線が地下水面に到達し，水分フラックスのピー

ここで，c∗ は Γs が正の場合（亀裂からマトリックスへ

クは Fig. 5(b) と同様およそ 10 d 前後に現れる．αw をの水分移動発生時）は cf ，負の場合は cm となる．また0.0001，0.001，0.01 と大きくし亀裂からマトリックスへ

c は溶質濃度，D は分散係数，q は水分フラックス，ωdpの水分移動が増えると，亀裂を流れる水分フラックスのは二領域間の濃度差を駆動力とする溶質移動係数である

地下水到達時時間がわずかに遅れ，フラックスのピーク（Gerke and van Genuchten, 1993）．次に DP モデルによが減少する．亀裂とマトリックス間での水分交換を伴うる地下水涵養過程の水分 · 溶質移動について，ωdp = 0 と場合，亀裂で浸潤前線が進行すると同時に，マトリック

して有効拡散による移動がない単純な場合を用いて，前スへの水分移動が生じ，その結果としてマトリックスか

節の均一モデルに対して αw による補正の効果を示す．

らの水分フラックスのピークが大きくなり，ピーク到達

土壌の物理性

20

第 153 号（2023）

Fig. 9 シルトで地下水面深さ L = 200 cm のとき Dual Permeability モデルにより求められた，異なる一次水分移
動係数 αw における亀裂とマトリックスの地下水面深さにおける水分フラックスの変化．一部 αw = 0.01 についての結果を拡大して示す．

Fig. 10 シルトで地下水面深さ L = 500 cm のとき Dual Permeability モデルにより求められた，地下水面
深さにおけるトレーサの溶質フラックスの変化．初期の溶質フラックスの立ち上がりについては拡大して示す．

解説：地下水涵養過程の数値解析再訪

21Fig. 11 シルトで地下水面深さ L = 200 cm で αw = 0.01 のとき Dual Permeability モデルにより求められた，亀裂内の（a）土中水圧力水頭分布，（b）亀裂とマトリックス間の溶質移動（正が亀裂からフラックス）分布，（c）マトリックス内の溶質濃度分布．

4. まとめ

時間が早くなる．αw = 0.01 のとき，1/αw = 100 d となり領域間の水分移動の時間スケールが移動時間のスケー
ルと近くなり，下方フラックスの移動と比較して亀裂と本報では，まず一定の深さに地下水面が存在する均一

マトリックス間の水分移動が十分に速く進行するため，な土層を想定し，リチャーズ式による降雨による浸潤と結果として亀裂内の水分移動が遅れ，マトリックス内の排水の地下水涵養過程の数値計算に対して，塩沢 · 中野

水分移動が速くなり，ピークが発生する時間がほぼ一致（1981）による伝達モデルを適用し，重力流れの役割とする（Fig. 9）．

モデルの有用性について示した．伝達モデルにより，排Fig. 10 に，地下水面深さ（z = −L）での全溶質フラッ

水過程において，地下水の影響が及ばない土層上部ではクスの時間変化を示す．亀裂とマトリックスで水分移動ある水分フラックスを与える座標の移動速度は一定であ，ωdp = 0 と拡散を考慮していないがない場合（αw = 0）

るのに対し，地下水面に向けて水分量が増加する領域でため領域間の溶質移動もなく， 亀裂からの極めて速い溶は速くなる特性を確認できた．また，地下水面における質の流出とマトリックスからの極めて遅い流出に二分さ

水分フラックス，すなわち地下水涵養フラックスの変化

れ，亀裂を通過する溶質の地下水面への到達時時間は，は，伝達モデルと数値計算の結果はよく一致し，伝達モ

マトリックスのみを考慮したときの溶質の地下水への到デルによる現象理解の有用性を示した．伝達モデルは地

達時間は約三十分の一程度まで速くなる．一方，αw が下水涵養過程における水分移動の特徴を理解する上では大きくなるにつれて，溶質移動はマトリックスで速く，有用であるが，限られた条件のもとでの解析に限定され，亀裂で遅くなるため，亀裂とマトリックスからの溶質フ汎用性は必ずしも高くない．リチャーズ式に基づく数値

ラックスのピークは近づく．Fig. 11 に αw ＝ 0.01 のと解析では，伝達モデルによる涵養フラックス変化を再現

きの亀裂内の土中水分圧力水頭の変化と，そのときの亀することができ，数値解析モデルは，溶質や熱移動モデ

裂 · マトリックス間の溶質交換分布（正は亀裂からマト

ルとの同時解析や，様々な水分特性関数（水分保持関数

リックスへの移動）およびマトリックス内の溶質分布をと不飽和透水係数）を用いることが可能で，現象の理解示す．Fig. 7 のようにマトリックス内の水分と完全混合

のみならず実問題への適用など汎用性が高く極めて有用する場合と比較して，亀裂内を速い速度で下方へ移動すである．

る水分によって，亀裂からマトリックスへの溶質移動もたとえば降雨として土層へ侵入する「新しい」水と，大きく，その結果として 50 d 後の表層近くのマトリック

もともと土層に存在する「古い水」を区別するために，

ス内の溶質分布に加えて，深さ 150 cm より深いところ
まで溶質が到達する．

新しい水を土粒子と反応しないトレーサーの溶質移動問
題として扱い，CDE に基づく数値解析を行うことができ土壌の物理性

22

第 153 号（2023）

る．その結果，新しい水の移動が圧力伝播としての浸潤

Gerke, H.H. and van Genuchten, M.Th. (1993): A dual-

前線の移動に対して大幅に遅れることを示した．数値解

porosity model for simulating the preferential movement

析手法を適切に使うことで，地下水涵養過程において新of water and solutes in structured porous media. Water Re-

しい水と古い水を区別した解析も可能となる．

sources Research, 29(2): 305–319.

また，遅いマトリックス流と速い亀裂流れを考慮可能

Huang, Y., Evaristo, J. and Li, Z. (2019): Multiple tracers

な DP モデルが実装されている数値解析プログラムを用

reveal different groundwater recharge mechanisms in deep

いることで，選択流の発生が問題となる自然条件下の土

loess deposit. Geoderma, 353: 204–212.

に対して，均一な間隙構造をもつマトリックスを仮定しKlaus, J., Zehe, E., Elsner, M., Kulls, C. and McDonnell,

た水分 · 溶質移動解析結果の補正が可能である．本報でJ.J. (2013): Macropore ﬂow of old water revisited: experi-は，DP モデルの概要と適用事例について，特にマトリックスと亀裂の間の水分および溶質の交換が地下水涵養に与える影響に焦点を絞り解説した．亀裂に伴う選択流によって，地下水に素早く浸潤前線が到達し，マトリック
スと亀裂の間の水分と溶質の交換の程度によって，特にマトリックス中の水分 · 溶質移動が影響を受けることが示された．なお有効拡散を含む例（ωdp ̸= 0）について

mental insights from a tile-drained hillslope. Hydrological
Earth System Sciences, 17: 103–118.
Jury, W.A. and Horton, R. (2006): 土壌物理学 — 土中の水 ·熱 · ガス · 化学物質移動の基礎と応用 —（取出伸夫 監訳:
井上光弘 · 長 裕幸 · 西村 拓 · 諸泉利嗣 · 渡辺晋生訳）, 築地書店, 東京．
塩沢 昌, 中野政詩 (1981): 不飽和層の排水過程としての地下
水涵養機構. 農業土木学会論文集, 92: 35–42.

は，本報では扱わないが，パラメータ ωdp の与える影響Šimůnek, J., Šejna, M., Saito, H., Sakai, M. and vanや適用事例の検討などは今後の課題とする．また亀裂にGenuchten, M.Th. (2013): The Hydrus-1D software pack-

よる選択流が発生する場合でも，新しい水，古い水の動

age for simulating the movement of water, heat, and multi-

態を知ることは溶質移動の観点から極めて重要である．

ple solutes in variably saturated media, version 4.17, HY-

Fig. 12 に，Klaus et al. （2013）による，亀裂を伴う土中

DRUS Software Series 3, pp. 342. Department of Environ-

での水分移動の概念図を示す． 亀裂を含む土層に降雨がmental Sciences, University of California Riverside, River-

あると新しい水が亀裂を流れながらマトリックス内へとside, California, USA.

移動するが，マトリックスが十分に湿ってくると古い水Šimůnek, J. and van Genuchten, M.Th. (2008): Modelingがマトリックスへ移動し，亀裂を浸潤する新しい水と混

nonequilibrium ﬂow and transport processes using HY-

合しながら流出するとしており，DP モデルを活用する
ことで，Fig. 12 に示す概念を実際に数値計算することがDRUS, Vadose Zone Journal. 7(2): 782–797.
取出伸夫, 渡辺晋生, 坂井 勝 (2009): 土中への水の浸潤

1．

フラックス境界と圧力境界条件. 土壌の物理性, 31–41.

できる．

Wang, Z., Feyen, J., van Genuchten, M.Th. and Nielsen, D.R.
(1998): Air entrapment effects on inﬁltration rate and ﬂow
instability, Water Resources Research, 34(2): 213–222.

付録
ここでは，塩沢・中野（1981）に基づいて，本文中の式（1）および式（3）を導出する．式（1）の水分フラッ
クス q となる座標 z の移動速度は水分フラックス q の変化率に等しい関係は，水の連続式の別の表現形である．
不飽和土中の連続式は，

Fig. 12 亀裂とマトリックス間の物質移動の概念図（Klaus
et al.（2013）を改変）．白矢印が新しい水，灰色矢印が古い
水と新しい水の混合物を表し，色が濃いほど古い水の割合が
高い．

∂θ
∂t



 
∂q
=−
∂z t
z

(A1)

ここで，q は上向き正の水分フラックス，z は上向き正の座標，θ は体積含水率，t は時間を表し，下付き添字（t引用文献と z）は固定する変数を示す．ここでは q と θ はともにz と t の関数であり，式（A1）の連続式は，位置 z における θ の時間変化は，時刻 t における q 分布の勾配にマ

Allaire, S.E., Roulier, S. and Cessna, A.J. (2009): Quantifying

イナスをかけたものに等しいことを示す．水分フラック

preferential ﬂow in soils: A review of different techniques.

ス q が一定となる z の移動速度は t の関数であるため，

Journal of Hydrology, 378: 179–204.
新井 正 (1989): 東京の小河川の変化と水収支. 水質汚濁研
究, 12(7): 7–10.

q(t, z) = q(t, g(t)) = const と書ける．このとき，q の時間微分はゼロであり，陰関数の定理を用いると次式となる．

解説：地下水涵養過程の数値解析再訪

dq
=
dt




=

∂q
∂t
∂q
∂t



なお，塩沢 · 中野（1981）では q を下向き正としている
ため，式（A1）と式（A9）の符号が逆である．

 dt∂ q dz+∂ z t dt
z dt




+
z

∂q
∂z

式（3）は，重力のみを考慮した水分移動式である本文

(A2)



dz
=0
t dt

dz
=−
dt

を偏微分として表すと，

∂z
∂t



!

q


=

z

(A4)

dθ
=
dt

dq
=
dt

∂q
∂z



dz
∂q
+
∂θ
θ dt




また，q(z(t), θ (t)) の時間微分は次式となる．






=

t

=



dθ
z dt

(A5)

∂θ
∂t




=
z

∂q
∂z

t

dK
dθ


(A10)
t

∂θ
∂t



 
∂ θ dzdt+∂ z t dt
z dt





∂θ
+
∂z
z

ここで，q と θ は z と t の関数であるので，ある位置 z

∂q
∂t



  dz∂q∂θ+∂θ z ∂t z
θ dt

dz
=−
dt



(A6)

(A11)



dz
=0
t dt

したがって，を固定した式（A5）は偏微分を用いて表せる．



∂θ
∂z

θ 分布の勾配に透水係数の θ に関する勾配をかけたものと等しいことを示す．この関係は，式（A1）の右辺と係数は異なるが q を θ に置き換えた関係と同じである．θ が一定となる z の移動速度は t の関数であるため，
θ (t, z) = θ (t, g(t)) = const と書ける．このとき θ の時間微分はゼロであり，式（A2）と同様に陰関数の定理を用t∂ K(h)∂zいると次式となる．

!z
∂q
∂z

=−



は，位置 z における θ の時間変化は，時刻 t における

(A3)

ここで，z を q と t の関数 z(q,t) として式（A3）の左辺



∂θ
∂t

ここで，θ は z と t の関数であり，式（A10）の移動式

t

∂q
∂t

中の式（2）より得られる．



したがって，

!
∂q
∂t
!z
∂q
∂z

23

∂θ
∂t
∂θ
∂z

!
!z

(A12)

t

ここで，z を θ と t の関数として式（A12）の左辺を偏微分によって表すと，また，式（A6）は z を固定したときの時間微分であるた
め dz/dt = 0 となり，右辺第一項はゼロである．



∂q
∂t




=
z

∂q
∂θ


z

∂θ
∂t




(A7)
z

式（A7）を式（A4） に代入すると，



∂z
∂t


=−

∂q
∂θ

q

!

∂θ
∂t
z
!
∂q
∂z


θ

=−

∂θ
∂z

!
!z

(A13)

t

式（A10）を式（A13）に代入すると，本文中の式（3）が

!
z

∂z
∂t

∂θ
∂t

得られる．

(A8)

dK
∂z
=−
∂t
dθ

(A14)

t

さらに式（A1）を用いると，本文中の式（1）が得られる．



∂z
∂t




=
q

∂q
∂θ


(A9)
z

引用文献
塩沢 昌, 中野政詩 (1981): 不飽和層の排水過程としての地下
水涵養機構. 農業土木学会論文集, 92: 35–42.

24

土壌の物理性

第 153 号（2023）

要旨

降雨後の不飽和帯の下方水分移動は，地下水面に到達すると地下水涵養流となるが，定量的な測定は困難で，リチャーズ式に基づく数値解析による評価が有効である．本報では，まず降雨浸透後の地下水涵
養過程では土中水圧力勾配が小さくなり重力が主な駆動力となるため，重力のみを考慮した塩沢 · 中野
（1981）による伝達モデルについて解説した．そして，HYDRUS-1D による数値解析の結果を伝達モデ
ルと比較し，数値解析手法の有用性について考察した．一方，自然条件下の土に多く観察される局所的な選択流が，地下水涵養過程にあたえる影響について，遅いマトリックス流と速い選択流を考慮可能なDual Permeability モデルを適用し，その特徴を示した．
キーワード：リチャーズ式，重力流れ，Dual Permeability モデル，数値計算