博士論文
土中の不飽和水分移動と水分移動特性関数に関する研究
Soil water hydraulic functions forwater flow in variable saturated soils

平 成 19 年 9 月

三重大学大学院生物資源学研究科

坂井 勝

目次
第 1 章

序章
1.1 はじめに1.2 不 飽 和 土 中 の水 分 移 動 式 と HYDRUS
1.3 不 飽 和 水 分 移 動 特 性 関 数
1.3.1 Mualem の間 隙 径 分 布 モデル
1.3.2 van Genuchten モデル
1.3.3 van Genuchten モデルの問 題 点 と修 正 モデル(1) 高 水 分 領 域(2) 低 水 分 領 域
1.4 不 飽 和 水 分 移 動 特 性 の測 定 と蒸 発 法
1.5 本 研 究 の目 的 と方 法

1
2
3
4
5

8
9

第 2 章

蒸 発 法 における水 分 移 動 特 性 推 定 の最 適 条 件 の検 討
2.1 はじめに
12
2.2 数 値 計 算
132.3 応 答 局 面
152.4 蒸 発 法
2.4.1 蒸 発 実 験
16
2 . 4 . 2 パラ メ ー タ の 推 定
17
2.5 結 果 と考 察
2 . 5 . 1 パラ メ ー タ の 推 定
182.5.2 応 答 局 面
20(1) 観 測 位 置(2) 蒸 発 速 度(3) 試 料 長
2.5.3 蒸 発 法 によるパラメータの推 定
26(1) θ r が不 飽 和 透 水 係 数 に及 ぼす影 響(2) 試 料 長 がパラメータの推 定 に与 える影 響(3) 観 測 位 置 がパラメータの推 定 に与 える影 響
2.6 まとめ
31

第 3 章

砂 質 土 と団 粒 土 に対 する水 分 移 動 特 性 関 数 の検 討
3.1 はじめに3.2 水 分 移 動 特 性 関 数
3.2.1 Mualem モデル
3.2.2 水 分 保 持 曲 線 ・不 飽 和 透 水 係 数 連 結 モデル(1) van Genuchten モデル(2) Fayer モデル(3) Durner モデル
3.3 蒸 発 法
3.3.1 蒸 発 実 験
3.3.2 逆 解 析(1) 水 分 移 動 式(2) パラメータの推 定

33
35
35

38
40

3.4

3.5
付録

結 果 と考 察
3.4.1 水 分 保 持 曲 線 関 数
3.4.2 鳥 取 砂 丘 砂(1) 砂 丘 砂 に対 する水 分 移 動 特 性 関 数(2) 水 蒸 気 移 動 を考 慮 した推 定(3) 高 水 分 領 域 の不 飽 和 透 水 係 数 の推 定
3.4.3 熊 本 黒 ボク土
3.4.4 不 飽 和 透 水 係 数 の信 頼 性
3.4.5 van Genuchten モデルの適 用 範 囲
まとめ

43
44

50
52
54
55
57

第 4 章

温 度 勾 配 下 における凝 縮 ・蒸 発 を伴 う砂 中 の水 分 移 動
4.1 はじめに
4.2 試 料 と方 法
4.2.1 凝 縮 実 験
4.2.2 数 値 計 算(1) 不 飽 和 水 分 移 動 式(2) 不 飽 和 水 分 移 動 特 性(3) 熱 移 動 式(4) 溶 質 移 動 式
4.2.3 初 期 ・境 界 条 件
4.2.4 不 飽 和 水 分 移 動 特 性 パラメータ
4.3 結 果 と考 察
4.3.1 不 飽 和 透 水 係 数
4.3.2 温 度 分 布
4.3.3 液 状 水 ・水 蒸 気 フラックス
4.3.4 表 面 境 界 条 件
4.3.5 凝 縮 ・蒸 発 速 度
4.3.6 水 蒸 気 促 進 係 数
4.3.7 溶 質 分 布
4.4 まとめ
第 5 章
結論
記号
参考文献
謝辞
要旨
Abstract

60
61
62

66
67
69
71
72
72
75
77
79
79
81

第1章

1.1

序 章

はじめに

近年，地球規模の水循環の把握に注目が集まっており，不飽和帯における土中水分移動の重要性が認識されてきている。不飽和帯は，降雨浸透，地表面蒸発，地下水涵
養，植物根による吸水，土壌浸食等に対して非常に重要な役割を果たしている。不飽和
土中の水分移動については，畑地において根圏へ効率よく水を供給し，農作物の収量を増加させるために，古くから多くの研究が行なわれてきている。近年では，工場からの産業廃棄物による土壌汚染，それに続く地下水汚染，さらには畑地における過剰な肥料や農薬の投入による地下水の富栄養化，及び汚染という環境問題の視点からも不飽和
土中の水分・溶質移動に対する研究が進められてきた。
これまでの研究の成果により，不飽和土中の水分移動及び溶質移動についての偏微分方程式が構築され，それを解くことによって土中水分・溶質移動の予測が可能となった。さらに，数値計算ソフトの汎用化により，強い非線形性を持つ偏微分方程式を短時間で手軽に 解く こ とが でき る ように なっ て い る 。 不飽和水分移動汎用プログ ラム
HYDRUS-1D（Šimůnek et al., 2005）では，1 次元土中の水分・溶質・熱移動のシミュレー
ションが可能である。また，コンピュータの高速化に伴い，2 次元，3 次元の土中水分移動（HYDRUS2D3D, Šimůnek et al., 2006）や，土中の温度変化に伴う熱・液状水・水蒸気
同時移動（Philip and de Vries, 1957; Saito et al., 2006），複雑な化学反応を伴う溶質移動（Parkhurst and Appelo, 1999）等，より詳細で複雑な計算が短時間で可能となっている。
土中水分移動を表現する微分方程式を解くためには，土の保水性，透水性を特徴付
けるパラメータが必要である。土の保水性は体積含水率と土中水圧力の関係である水分保持曲線で表され，土の透水性は体積含水率や土中水圧力の関数である不飽和透水係数で表される。この両者を合せて，soil water hydraulic function（Jury and Horton,
2006）等とよばれ，本研究では，「不飽和水分移動特性」と定義する。また，地盤工学の分野では不飽和浸透特性と定義されている。
土の不飽和水分移動特性を測定する方法について，これまで様々な研究が行なわれてきた。しかし，不飽和透水係数の測定には多くの労力と時間を要し，特に低水分領域
1の測定は非常に困難である。また，不飽和透水係数を表現する関数型には，土性や土の構造に関わらず，後述する van Genuchten-Mualem モデルに代表される単一の関数のみが用いられる場合がほとんどである。そこで本研究では，より正確に土中水分移動予測が可能となるように，不飽和水分移動特性の測定法，及びその関数型に焦点を当てた。

1.2

不飽和土中の水分移動式と HYDRUS

鉛直 1 次元の不飽和土中の水分移動は，土中水の全水頭（h + z）の勾配により生じ，
次式のバッキンガム‐ダルシー則で表すことができる（e.g., Jury and Horton, 2006）。
qz = − K ( h )∂ (h + z)⎛ ∂h ⎞= − K ( h ) ⎜ + 1⎟∂z
⎝ ∂z ⎠(1.1)

h は土中水の圧力水頭（L），z は地表面を 0 とする上向き正の空間座標（L）および土中
水の重力水頭（L），K（h）は不飽和透水係数（L T−1），qz 鉛直方向の水分フラックス（LT−1）である。ここで水頭とは，単位重量当たりの水のエネルギーであり，長さの次元（L）を持つ。また，バッキンガム‐ダルシー式を水保存則と結合させることで，リチャーズ式が導かれる。
∂θ ∂ ⎡⎛ ∂h ⎞ ⎤= ⎢ K (h) ⎜ + 1⎟ ⎥ − S∂t ∂z ⎣⎝ ∂z ⎠ ⎦(1.2)

ここで，θ は体積含水率（L3 L−3），t は時間（T），S は植物根による吸水などを示す吸水項である。拡散方程式であるリチャーズ式に，適切な初期条件・境界条件（斎藤ら，2006）
を与えて解くことで，体積含水率や土中水圧力の分布など，不飽和土中の水分移動を予測することができる。また吸水項 S について，植物根の吸水モデル（e.g., Feddes et al.,
1978; van Genuchten, 1987）と組み合わせて解くことで，根の吸水による水分移動を予測することも可能である。さらに，次式に示す水蒸気移動を考慮したより一般的なリチャー
ズ式（Philip & de Vries, 1957）を解くことで，日々の温度変化に伴う土中の水分移動を予測することができる。
∂θ ∂ ⎡∂T∂h∂T ⎤⎛ ∂h ⎞= ⎢ K Lh ( h ) ⎜ + 1⎟ + K LT
+ K vh+ K vT−S∂t ∂z ⎣∂z∂z∂z ⎥⎦⎝ ∂z ⎠(1.3)

ここで，KLh は圧力勾配による不飽和透水係数（L T−1），KLT は温度勾配による不飽和透2水係数（L2K−1T−1），Kvh は圧力勾配による水蒸気透過係数（L T−1），KvT は温度勾配による水蒸気透過係数（L2K−1T−1）である。
ところが，リチャーズ式は不飽和透水係数 K(h)により，強い非線形性を示す。このため，特殊な条件（例えば，雨が常に一定の強度で降り，根による吸水はないなど）を除いて，解析解は得られない。自然界における土中水分移動は降雨や蒸発を繰り返すように，一定な条件では表現できない。このような自然界の変動的な条件のもとで土中水分移動を予測するためには，リチャーズ式を数値的に解く必要がある。U.S. Salinity Laboratry（米国農務省塩類研究所）において開発された HYDRUS は，リチャーズ式の数値解を導く不飽和水分移動汎用プログラムである。近年では，1 次元水分移動（HYDRUS-1D,Šimůnek et al., 2005）に加え，2 次元や 3 次元の数値計算プログラムも開発されている（HYDRUS2D3D, Šimůnek et al., 2006）。また，安定した計算に加え，インターフェイスが整備されているため，世界中で研究・教育・コンサルタント等の広い分野で用いられている。
本研究では HYDRUS-1D を用いた数値計算により，不飽和土中の水分移動の数値解を導いた。また，カリフォルニア大学リバーサイド校の Šimůnek 博士との共同研究として，HYDRUS-1D を改良して用いた。

1.3

不飽和水分移動特性関数

1.3.1 Mualem の間隙径分布モデル

不飽和水分移動特性関数は，水分保持曲線関数と不飽和透水係数関数で与えられるが，特に不飽和透水係数は測定が困難である。そのため古くから，比較的測定が容易な水分保持曲線に基づき，不飽和透水係数を推定する試みが行なわれてきた。その主流である間隙径分布モデルは，土中の間隙径分布に基づいた毛管モデルにより保水性を表現し，毛管内のポワズイユ流れを仮定することで異なる水分量の不飽和透水係数を評価する（Childs and Collis-George, 1950; Burdine, 1953; Mualem, 1976; Jury and Horton,
2006）。つまり，間隙径分布モデルに水分保持曲線関数を代入することにより，物理的根拠を持つ不飽和透水係数関数を導くことができる（水分保持曲線‐不飽和透水係数連結3モデル）。少ないパラメータ数で土の保水性と透水性を表す水分移動特性関数を表現できるため，間隙径分布モデルは非常に有用性が高い（小杉，2007）。また近年注目が集まっている，土中の水分移動過程における水分量や土中水圧力変化から水分移動特性関数のパラメータを推定する，マルチステップ法や蒸発法等の逆解析手法においても，少ないパラメータ数で水分移動特性を表現できることは非常に重要である（Eching andHopmans，1993; Šimůnek et al, 1998; 坂井・取出，2007）。
間隙径分布モデルの中でも広く用いられるのが，Mualem モデルである（Mualem,1976）。
⎡ Se 1
K ( Se ) = K s K r = K s Se ⎢ ∫
dSe
0 h
⎣

⎤
1
dS
e
∫0 h ⎥⎦

2

1

l

(1.4)

ここで，K は不飽和透水係数（L T−1），Ks は飽和透水係数（L T−1），Kr は比透水係数（−），
Se は有効飽和度（−）である。l は間隙結合係数（pore-connectivity coefficient）で，土中間隙の屈曲や連結による影響を補正するための係数であり，Mualem は 45 種類もの土試
料の平均値として l = 0.5 を推奨している。等価毛管半径に基づく場合，本来ならば（1.4）
式右辺の 2 乗の位置は，h2 となるべきである。小杉（2007）は，Mualem が（1.4）式を示した理由として，物理的根拠に基づくものではなく，積分が容易となる数学的な便宜性を強く意識した結果ではないかと考察している。（1.4）式に関数化した h と Se の関係，すな
わち水分保持曲線関数を与えることで，不飽和透水係数関数が導かれる。

1.3.2 van Genuchten モデル

Mualem モデルに基づく代表的な水分保持曲線‐不飽和透水係数連結モデルは，vanGenuchten（VG）モデルである（van Genuchten, 1980）。van Genuchten（1980）は，（1.4）式に代入したときに積分可能な水分保持曲線関数として次式を提案した。

()

θ − θrn −m= Se = 1 + α hθs − θr(1.5)

ここで，θ r は残留体積含水率（L3 L−3），θ s は飽和体積含水率（L3 L−3），α（L−1），n（−），

m（= 1−1 / n）は水分保持曲線の形状を与えるパラメータである。（1.5）式を（1.4）式に代入し積分することで，不飽和透水係数関数を解析的な表現（閉形式解，closed-form）が
導かれる。
4

m
K ( Se ) = K s Sel ⎡⎢1 − (1 − Se1 m ) ⎤⎥
⎣
⎦

2

(1.6)

VG モデルは，水分保持曲線の実測値に対し適合性が高く，また，簡単で滑らかな関数型で表されるため数値計算との相性がよいことから，非常に広く用いられている（小杉，2007）。

1.3.3 van Genuchten モデルの問題点と修正モデル(1) 高水分領域

VG モデルは非常に汎用性が高く，有用な不飽和水分移動特性関数であるが，その一方で問題点も指摘されている。粘質土のような細粒土が示す，傾きの緩やかな水分保持曲線に対して VG モデルの水分保持曲線関数（（1.5）式）を適合した場合，（1.6）式で表される不飽和透水係数が飽和近傍で大きく変化する（Vogel and Cislerova, 1988; 坂
井・取出，2007）。この非現実的な不飽和透水係数の変化は，数値計算を不安定にし，計算結果に誤りをもたらす（Rassam et al., 1998）。Vogel and Cislerova（1988）は，この飽和近傍の問題を解決するために，VG モデルを改良した不飽和水分移動特性関数を示
している。

(

)

n −m
⎧
⎪θ a + (θ m − θ a ) 1 + α h
θ =⎨
⎪⎩
θs

Ks Kr
⎧
⎪
( h − hk )( K s − K K )⎪K ( Se ) = ⎨ K k +hs − hk⎪⎪⎩Ks

h < hsh ≥ hsh ≤ hkhk < h < hs

5

(1.8)

h ≥ hsK ⎛ S ⎞ ⎡ F (θ r ) − F (θ ) ⎤Kr = k ⎜ e ⎟ ⎢⎥K s ⎝ Sek ⎠ ⎣⎢ F (θ r ) − F (θ k ) ⎦⎥⎡ ⎛ θ − θ ⎞1 m ⎤a
F (θ ) = ⎢1 − ⎜⎟ ⎥⎢⎣ ⎝ θ m − θ a ⎠ ⎥⎦

(1.7)

2

(1.9)

m

(1.10)

Sek =

θk − θrθs − θr(1.11)

ここで，θk，Kk はそれぞれ任意の土中水圧力 hk における体積含水率と不飽和透水係数の実測値である。このモデルでは，現実の飽和体積含水率θ s より大きな仮想の飽和体
積含水率θ m（L3 L−3）を用いて，空気侵入値（AEV）hs（L）を明確に定義することで，VG
モデルにおける飽和近傍の問題を解決している。HYDRUS においては特に，仮想的に

hs = −2 cm と設定し，問題を解決している（AEV=−2 cm オプション，Rassam et al., 1998）。
また，Vogel and Cislerova モデルと同様に明確な空気侵入値を示す Brooks and Corey（1964）モデルでは，飽和近傍の問題は生じない。

(2) 低水分領域
また，低水分領域についても問題点が指摘されている。VG モデルは無限小の土中
水圧力で水分量が残留体積含水率θr に漸近する水分保持曲線関数であるため，土中水圧力が非常に小さくなったとき体積含水率θ が 0 に近づく現実の水分保持曲線を十分に表現できない（Rossi and Nimmo, 1994; Fayer and Simmons, 1995）。Campbell andShiozawa（1992）は，低水分領域における土粒子表面の吸着保水量が，土中水圧力の対数軸に対して直線的に減少するモデルを示した。

⎡

θ = ⎢1 −⎣⎢ln ( − h ) ⎤⎥ θaln ( −hm ) ⎦⎥(1.12)

ここで，θa は h = 1 cm における吸着保水量であり，hm はθ = 0 となる土中水圧力であり，

Rossi and Nimmo（1994）は−10７ cm としている。Fayer and Simmons（1995）は，（1.12）式の吸着保水モデルを用いて VG モデルの残留体積含水率θr を修正し，−10７ cm の土中
水圧力で水分量が 0 となる水分保持曲線関数を示している。

θ = χ ( h )θ a + (θ s − χ ( h )θ a ) ⎡1 + α h ⎤
n

⎣

⎦

−m

(1.13)

ここで，χ (h) = 1 – ln (−h)/ ln (−hm)である。また，Fayer and Simmons（1995）は，（1.13）をMualem モデルに代入して得られる不飽和透水係数関数の閉形式解を示している（Fayer モデル）。

K ( Se ) = K s Sel ⎡⎣ Γ s ( h ) Γ max ⎤⎦62(1.14)

黒ボク土に代表される団粒構造を持つ土は，団粒間間隙と団粒内間隙における 2 種
類の保水形態を持ち，階段状の水分保持曲線を示すため，VG モデルでは低水分領域を十分に表現できない（徳本ら，2005）。Durner（1994）は，2 つ（もしくはそれ以上）の VG
モデルを結合することで，階段状の水分保持曲線関数を示している（Durner モデル）。

Se = Se1 + Se 2 = w1 ⎡1 + α h ⎤
⎣
⎦
n

−m

+ w2 ⎡1 + α 2 h 2 ⎤
⎣
⎦
n

− m2

(1.15)

ここで，w1 + w2 = 1 である。また，Priesack and Durner（2006）は（1.20）式を Mualem モデルに代入して得られる閉形式解を示している。

( w1Se1 + w2 Se 2 ) ⎛⎜ w1α ⎡⎣⎢1 − (1 − Se11m ) ⎤⎦⎥ + w2α 2 ⎡⎣⎢1 − (1 − Se12m ) ⎤⎦⎥ ⎞⎟ml

K ( Se ) = K s⎝( w1α + w2α 2 )

2

m2

2

⎠

2

(1.16)
さらに小杉（2007）は，（1.6）式で表される不飽和透水係数はθr の定め方に依存して大きく変化することを示している。また Tuller and Or（2001）は，低水分領域では間隙流による透水に比べ膜流による透水が卓越するため，Mualem の間隙径分布モデルにより導かれた不飽和透水係数（（1.6）式）は，θr に近い低水分領域で現実を過小評価することを示している。特に，比較的高土中水圧力で水分が急激に低下する砂質土に対して VG
モデルを適用した場合，速やかにθ r に漸近するため，不飽和透水係数を正しく評価できる土中水圧力範囲は限られる。

VG モデルに代表される閉形式解は，水分保持曲線関数から不飽和透水係数関数を導くことができるため非常に利便性が高い。しかし多くの場合，水分保持曲線にのみ注
目し，導かれた不飽和透水係数について十分な検討が行なわれないまま数値計算に用いられている。特に，低水分領域の不飽和透水係数について検討した研究は少なく，土中水圧力や水分量等の数値計算結果と実測値を比較して検証する必要がある。また，

Mualem の間隙径分布モデルに含まれる間隙結合係数 l は，水分保持曲線から独立して不飽和透水係数を決める重要なパラメータである。それにも関わらず，多くの場合

Mualem の推奨値である 0.5 に固定して用いられ，十分な検討がなされていない。

7

1.4

不飽和水分移動特性の測定と蒸発法

水分保持曲線については，吸引法や加圧板法（中野ら，1995），蒸気圧法（藤井・中野，1984）等により，幅広い圧力範囲での直接測定が可能である。一方，不飽和透水係数の測定は，不飽和流れや土中水圧力の制御が難しく，多くの時間を要する。不飽和
透水係数の直接測定法には，圧力制御法や流速制御法等の定常法（中野ら， 1995；

Klute and Dirksen，1986）が挙げられるが，飽和に比較的近い限られた圧力領域の測定に限られる。
近年，数値計算プログラムの高速化と汎用化に伴い，不飽和水分移動特性関数のパラメータを推定する逆解析法に注目が集まっている。逆解析法の 1 つに，土カラムに対して，試料内の空気圧を段階的に増加させ，積算流出量や試料内部の土中水圧力の測定値を用いて不飽和水分移動特性パラメータを求めるマルチステップ流出法（Echingand Hopmans，1993；van Dam et al.，1994）がある。マルチステップ流出法は，段階的な空気圧の調整を多くの試料に対して同時に適用して測定することが可能であるために広く用いられているが，正確なパラメータ推定には，試料下端に用いる多孔質フィルターの透水係数の評価が必要である（井上，2004）。また，圧力制御装置に加えて，テンシオメ
ーターによる土中水圧力測定部も含めた実験装置の耐圧性が必要である。
同様の逆解析法である蒸発法は，蒸発過程における土カラム中の土中水圧力変化の測定値から不飽和水分移動特性パラメータを求める方法である（Šimůnek et al，1998）。
蒸発法は，マルチステップ流出法に比べて実験装置は単純であり，飽和近傍から−1000

cm 程度の比較的低圧力水頭領域までの水分移動特性の推定が可能である。また，試料下端は水分フラックスが生じない条件であるため，下端のフィルターに関する評価が不要である。このため，蒸発法は不飽和水分移動特性を推定する上で，非常に有用性の高い方法だと考えられる。

Hopmans et al. （2002）は，わずかなパラメータの値の変化が地表面付近の土中水圧力の計算値に大きな影響を与えることから，地表面付近の土中水圧力データが推定に有効であることを示す一方，その測定の困難さも指摘している。Šimůnek et al.（1998）は，

Wendroth ら（1993）の 2 段階の蒸発強度を用いた蒸発実験は，実験時間を短縮させるが，パラメータの推定に対しては利点がないことを示した。Šimůnek et al.（1998）は，土カラムからの水分蒸発を想定した数値実験を用いて，VG モデル（（1.5）式，（1.6））の各パラメ
8

ータ（θ r，θ s，α，n，Ks）の変化が土中水圧力に及ぼす影響である感度とパラメータ間の相関を，様々な組み合わせの 2 個のパラメータに対する応答曲面（response surface）を用いて調べている。しかし，蒸発法の適用事例や実験条件である蒸発速度，試料長，推定に用いる圧力の観測地点の位置と数などの研究例は，マルチステップ流出法に比べて少ないのが現状である。このため，不飽和水分移動特性関数のパラメータをより正確に推定するためには，蒸発法における適切な実験条件について検討を行なうことは，重要な課題であると考えられる。

1.5

本研究の目的と方法

これまで，土中水分移動を評価する上で必要である不飽和水分移動特性の測定法，
関数型について数多くの研究が行なわれてきた。しかし，土中水分移動過程の水分量
や土中水圧力変化といった実測データに対して，数値計算の結果を比較して検証した例は少ない。また，水分移動特性関数に含まれる間隙結合係数 l は，不飽和透水係数を表現する上で重要なパラメータであるにも関わらず，十分な評価がなされて来なかった。そこで本研究では，実際の不飽和水分移動過程の土中水圧力変化に基づき，水分保持曲線‐不飽和透水係数連結モデルで表される不飽和水分移動特性を，逆解析手法
により推定する蒸発法に注目した。蒸発法によってより正確に不飽和透水係数を推定するために，数値実験及び室内実験を行ない，蒸発実験の実験条件が間隙結合係数 l の推定に与える影響を明らかにした。また，蒸発法に基づいた蒸発実験と逆解析より，砂質土や日本に広く分布する団粒構造を持つ黒ボク土に対して，より広水分範囲に適用できる不飽和水分移動特性関数について検討を行なった。図 1-1 に論文構成の流れ図を示す。
第 1 章では，不飽和土中の水分移動を予測する上で，不飽和水分移動特性を把握することの重要性を述べた。その上で，不飽和水分移動特性関数の既存のモデルとその問題点を整理した。
次章の第 2 章では，蒸発法において，試料長，可能蒸発速度，観測地点という実験条件の違いが VG モデルの各種パラメータの推定に与える影響について，特に不飽和透水係数のパラメータである間隙結合係数 l の推定を中心に整理する。まず，Šimůnek
9et al.（1998）と同様に，シルト質土カラムからの水分蒸発を想定した数値実験を行ない，間隙結合係数 l と他のパラメータ（θ r，θ s，α，n，Ks）の組み合わせによる土中水圧力変化に対する 2 次元応答曲面を示す。それにより，それぞれのパラメータの感度やパラメータ間の相関を評価する。そして，各実験条件がパラメータの推定の信頼性へ及ぼす影響について検討する。次に，藤森シルトを用いた蒸発実験を行ない，異なる条件でパラメータの推定を行なう。そして，蒸発法における試料長，観測点の位置と数が推定パラメータの信頼性に与える影響について考察する。
第 3 章では，広く用いられる VG モデルでは表現できない砂質土と団粒土を対象に不飽和水分移動特性の関数型，特に不飽和透水係数について検討する。VG モデ
ルの水分保持曲線関数を低水分領域の吸着保水を考慮して修正した Fayer モデルと，2 つの VG モデルを線形結合して階段状の水分保持曲線関数を示した Durner モデルを用いて，砂質土の鳥取砂丘砂と団粒構造の発達した黒ボク土を対象に蒸発法を行ない，パラメータの推定を行なう。そして，蒸発過程の土中水圧力の時間変化の実測値と計算値を比較，また水分保持曲線の実測値と推定値を比較することで，Mualem の間隙分布モデルから導かれる不飽和透水係数関数の適用性について明らかにする。
第 4 章では，非等温でより複雑な条件下における水分移動の実測値と計算値の比較により，不飽和水分移動特性関数の適用性について検証する。風乾した砂丘砂カラ
ムの上端を高温で湿った大気に開放し，大気からの温度勾配による砂中への水蒸気の浸入，凝縮を測定した水蒸気凝縮実験（宮崎，1979）の水分分布の実測値と，Fayer モ
デルを用いた計算値を比較することで，低水分量の砂質土に対する Fayer の不飽和透水係数関数の適用性について明らかにする。そして，正しく評価された不飽和水分移動
特性関数を用いた数値計算により，試料中の水分・熱・溶質移動を定量的に評価する。

10

序論
第1章

理論（モデル）

第2章

第3章

第4章

蒸発法の最適実験条件

砂丘砂の水分移動特性

黒ボク土の水分移動特性

凝縮実験
への適用

第5章

結論

図 1-1 本論文の構成
Fig1-1 Logical sequence of this study.

11

第２章

2.1

蒸発法における水分移動特性推定の最適条件の検討

はじめに蒸発法は，蒸発過程における土カラム中の土中水圧力変化の測定値から不飽和水
分移動特性パラメータを求める逆解析手法の１つである（Šimůnek et al，1998）。蒸発法
は，マルチステップ流出法に比べて実験装置は単純であり，飽和近傍から−1000 cm 程
度の比較的低圧力水頭領域までの水分移動特性の推定が可能である。また，試料下端は水分フラックスが生じない条件であるため，下端のフィルターに関する評価が不要である。Hopmans et al. （2002）は，わずかなパラメータの値の変化が地表面付近の土中水
圧力の計算値に大きな影響を与えることから，地表面付近の土中水圧力データが推定に有効であることを示す一方，その測定の困難さも指摘している。Šimůnek et al.（1998）
は，Wendroth ら（1993）の 2 段階の蒸発実験は，実験時間を短縮させるが，パラメータの推定に対しては利点がないことを示した。

Šimůnek et al.（1998）は，土カラムからの水分蒸発を想定した数値実験を用いて，VG
モデルの各パラメータ（θ r，θ s，α，n，Ks）の変化が土中水圧力に及ぼす影響である感度とパラメータ間の相関を，様々な組み合わせの 2 個のパラメータに対する応答曲面（response surface）を用いて調べた。そして，n，θ s は土中水圧力変化に対する感度が高
いこと，それに比べてθ r， α，Ks の感度は相対的に低いこと，また，θ r と n の間には強い正の相関があること，また，蒸発法の測定圧力水頭範囲（0 ≤ h ≤ −800 cm）を超えた範囲の水分保持曲線や不飽和透水係数は信頼性が低いことを示した。
一方，不飽和透水係数を表す Mualem モデルに含まれる間隙結合係数 l は，屈曲度の影響を示す土固有なパラメータである（Rassam et al., 2004）。しかし，多くの逆解析の報告では，Mualem（1976）が推奨した l = 0.5 が用いられ，推定パラメータに含まれていない（Inoue et al.，1998；Šimůnek et al.，1998）。蒸発法は，アメリカ土壌学会の測定法を示した Methods of soil analysis Part 4（Hopmans et al., 2002）においてもマルチステップ
流出法と並び紹介されている。しかし，蒸発法の適用事例や実験条件である蒸発速度，
試料長，推定に用いる圧力の観測地点の位置と数などの研究例は，マルチステップ流出法に比べて少ないのが現状である。特に，不飽和透水係数に関する間隙結合係数 lの推定に焦点を当てて検討した研究例はない。
12

そこで本章では，蒸発法において試料長，可能蒸発速度，観測地点の違いが水分移
動特性関数の各種パラメータの推定に与える影響について，特に不飽和透水係数のパラメータである間隙結合係数 l を中心に検討を行なった。まず，Šimůnek et al.（1998）と同様に，シルト質土カラムからの水分蒸発を想定した数値実験を行ない，間隙結合係数 lと他のパラメータの組み合わせによる土中水圧力変化に対する 2 次元応答曲面を示した。
それにより，それぞれのパラメータの感度やパラメータ間の相関を評価した。そして，各実験条件がパラメータの推定の信頼性へ及ぼす影響について検討した。次に，藤森シ
ルトを用いた蒸発実験を行ない，異なる条件でパラメータの推定を行なった。そして，蒸発法における試料長，観測点の位置と数が推定パラメータの信頼性に与える影響について考察した。

2.2

数値計算

水分飽和した長さ L の土カラムに対して，地表面からの水分蒸発過程の数値計算を行なった。鉛直 1 次元で水蒸気移動を考慮しない不飽和土中の水分移動は，リチャーズ
式で表される。

∂θ ∂ ⎡⎛ ∂h ⎞ ⎤= ⎢ K ( h) ⎜ + 1⎟ ⎥∂t ∂z ⎣⎝ ∂z ⎠ ⎦(2.1)

ここで，θ は体積含水率（L3 L−3），h は土中水の圧力水頭（L），K は不飽和透水係数（L T−1），t は時間（T），z は上向き正の空間座標（L）であり地表面で 0 である。
初期条件は，蒸発開始時の初期土中水圧力 hi（z）を与えた。
h ( z , 0 ) = hi ( z )(2.2)

カラム下端は水分フラックスが生じないように閉じているので，下端境界はゼロフラック
ス境界条件とした。

⎛ ∂h ⎞− K ( h ) ⎜ + 1⎟=0
⎝ ∂z ⎠ z =− L

(2.3)

ここで，L は試料長（L）である。
地表面蒸発を表す上端境界条件には，次の 2 種類の境界条件を用いた（斎藤ら，2006）。
13⎛ ∂h ⎞− K ( h ) ⎜ + 1⎟ = qevap ( t )⎝ ∂z ⎠ z =0⎧
⎛ ∂h ⎞⎪ − K ( h ) ⎜ + 1⎟ = E⎝ ∂z ⎠ z =0⎨⎪h(0, t ) = hA
⎩

(2.4)

( h ( 0, t ) > h )
A

(2.5)

（2.4）式は変動フラックス境界条件であり，各時間の蒸発速度 qevap（t）（L T−１）を上端境界に与えた。（2.5）式は土表面の乾燥状態に依存する系依存型境界条件であり，地表面
圧力 h（0, t）が許容最小圧力 hA（L）より大きいときは可能蒸発速度 E（L T−１）の一定フラックス境界条件を与え（恒率蒸発過程），また，地表面の乾燥が進行して h（0, t）が hA に達した段階で h（0, t） = hA の一定圧力境界条件を与える（減率蒸発過程）。蒸発法によるパラメータの推定においては，地表面の境界フラックスを，実験条件に忠実に再現することが重要である。そのため，本解析の蒸発実験の実測データに対するパラメータの推定には，（2.4）式の変動フラックス境界条件を用いた。また，数値実験に対しては（2.5）式の系依存型境界条件を用いた。
水分保持曲線と不飽和透水係数の水分移動特性関数には，van Genuchten（VG）モデル（van Genuchten，1980）を用いた。

Se =

(

)

θ − θrn −m= 1+ α hθs − θrm
K ( h ) = K s Sel ⎡⎢1 − (1 − Se1 m ) ⎤⎥
⎣
⎦

(2.6)
2

(2.7)

ここで，Se は有効飽和度，θ r は残留体積含水率（L3 L−3），θ s は飽和体積含水率（L3L−3），α（L−1），n（−），m（= 1−1 / n）は水分移動特性の形状を与えるパラメータ，Ks は飽
和透水係数（L T−1）である。l は間隙結合係数（pore-connecting coefficient）であり，水分保持曲線とは独立に不飽和透水係数の圧力依存性を与えるパラメータであり，Mualem（1976）の推奨値は l = 0.5 である。本研究における水分移動の数値計算と後述のパラメータの推定は，不飽和水分移動汎用プログラム HYDRUS-1D（Šimůnek et al., 2005）を
用いた。

14

2.3

応答曲面

試料長，可能蒸発速度，観測点の位置や数の違いが水分移動特性関数の推定に及ぼす影響を調べるために，土カラムからの水分蒸発を想定した数値実験を行なった。試料長は 10 cm と 20 cm の 2 種類とし，VG モデルのパラメータとしてθ r = 0.034，θ s = 0.46，

α = 0.016 cm−1，n = 1.4, Ks = 6 cm d−1，l = 0.5 を持つシルト質土（Carsel and Parrish,1988）を対象とした。初期条件として，試料下端の土中水圧力を 0 cm とし，鉛直方向の平衡圧力分布として h (z, 0) = −L − z を与えた（（2.2）式）。上端境界条件には（2.5）式を用い，可能蒸発速度として，通常の実験条件において大きな値である E = 1 cm d−1 を最大に，E = 0.5，0.1 cm d−1 の 3 条件を与えた。なお，許容最小圧力は hA = −105 cm とした（Rassam et al., 2004）。土中水圧力の観測位置は，技術的に測定可能な最も浅い深さとして深さ 0.5 cm，地表面に近い深さ 1 cm，表層の影響が小さい深さ 9 cm の 3 地点とし，計算は，深さ 0.5 cm の土中水圧力が−800 cm に達するまで行なった。
測定した土中水圧力変化に基づきパラメータ値を逆解析で推定する場合，与えたパラメータ値に対して計算を行なう順解析において，パラメータ値の変化が土中水圧力に及ぼす感度を把握しておくことが重要である（Šimůnek et al.,1998）。そこで，計算に用いたシルト質土の水分移動特性パラメータ値からそれぞれのパラメータの値を変化させ，土中水圧力変化の計算結果に及ぼす感度とパラメータ間の相関を調べた。まず，任意に選択した 2 個のパラメータ値のみを変化させた（それ以外のパラメータ値はシルト質土の値で固定）仮想土を考える。そして，同じ初期，境界条件による計算を行ない，観測地
点 j における仮想土とシルト質土の土中水圧力の平方残差φ j を求めた。
nj

φ j = ∑ ⎡⎣ p*j (ti ) − p j (ti ) ⎤⎦

2

(2.8)

i =1

ここで，nj は観測地点 j における土中水圧力データの数，p*j (ti)，pj (ti)はそれぞれ観測地
点 j の時刻 ti における仮想土とシルト質土の土中水圧力を示す。そして，変化させた 2個の水分移動特性パラメータに対し，平方残差φ j の等高線図，すなわち 2 次元の応答曲面（Response Surface）を作成した。VG モデルは 6 個のパラメータがあるため，合計 15
種類の応答曲面が存在するが，本研究では，特に不飽和透水係数を決定する間隙結合係数 l の応答曲面（l-α，l-Ks，l-n，l-θ r，l-θ s）に注目した。

15

2.4

蒸発法

2.4.1 蒸発実験

試料には藤森シルトを用いた。図 2-1 は，吸引法，蒸気圧法で測定した水分保持曲線である。実線は，非線形最小 2 乗法で実測値に適合した VG モデル（ θ r = 0.02，θ s =0.574，α = 0.00235，n = 1.4）である（van Genuchten et al., 1991）。藤森シルトは，比較的明瞭で小さな空気侵入圧（約−100 cm）を持つのが特徴である。変水位透水試験による飽和透水係数は Ks = 5.9 cm d−1 であった。
内径 3.9 cm で高さ 8.3 cm，16 cm の 2 種類のアクリルカラムに対して，含水比 0.041 kgkg−1 の藤森シルトを乾燥密度ρ b = 1.16 Mg m−3 で充填した（以降，それぞれ 8.3 cm カラム，16 cm カラム）。そして，マリオット管を用いて地下水位を徐々に上げ，3 日間かけて毛管飽和した。毛管飽和後，マリオット管を切り離して試料下端を閉鎖し，試料表面に小型ファンを用いて一定の風をあてながら蒸発を促進させた。8.3 cm カラムには深さ 0.7，

3.1，5.2 cm，16 cm カラムには深さ 1，5，13 cm の位置にテンシオメーターを水平方向に

Water content, θ (cm3cm-3)

0.6

0.4

0.2
observed
fitted
0
100

101

102

103

104

105

106

Pressure head, h (-cm)
図 2-1 藤森シルトの水分保持曲線の実測値と VG モデルによる適合曲線。
Fig. 2-1 Observed and fitted water retention curve for Fujinomori silty soil. Fitted line is the
VG model.
16

挿入し，圧力変換器を用いて土中水圧力を測定した。また，カラム全体を電子天秤の上に設置し，重量変化から積算蒸発量を求めた。圧力変換器，電子天秤は Campbell 社のデータロガーCR1000 に接続し，試料内の土中水圧力，積算蒸発量を自動計測した。
試料の上部が乾燥し，最上部のテンシオメーターが−700 ~ −800 cm 程度に達して測定不能になった段階で，実験を終了した。その後，カラムを解体し，炉乾法で試料の平均
体積含水率θ ave を測定した。

2.4.2 パラメータの推定

蒸発実験の実測データ対して，（2.6），（2.7）式のθ r，θ s，α，n，Ks，l の 6 個のパラメー
タを逆解析により推定した。本研究においては，不飽和透水係数の推定が主な目的であるため，間隙結合係数 l も推定パラメータに含めた。間隙結合係数 l が不飽和透水係数の圧力依存性に関わるのに対して，飽和透水係数 Ks は，不飽和透水係数の最大値を与えるパラメータである。藤森シルトの Ks は変水位透水試験により得られているが，Ks
は飽和土中の封入空気量に依存して変化する。封入空気量は，試料長，飽和時間，飽
和方法によって異なり，Ks は実験条件により変動することが予測される。そのため，本解析では，飽和水分量θ s に加えて Ks の推定も行なった。また，残留体積含水率θ r も，

2.5.3（1）節で述べるように不飽和透水係数の形状に影響を及ぼすパラメータであるため，
推定パラメータとした。
計算は，蒸発実験を再現する条件を与えた。まず，試料長，観測地点を実験と等しく与え，（2.2）式の初期圧力には，各深さの測定値を線形補間して与えた。（2.4）式の上端境界条件に対しては，実測した積算蒸発量の 5 時間毎の差に基づいた蒸発速度 qevap (t)を与えた。
逆解析には，土中水圧力の変化，実験終了時の平均体積含水率（8.3 cm カラムは 56
時間後でθ ave = 0.40 cm3 cm−3，16 cm カラムは 82 時間後でθ ave = 0.42 cm3 cm−3）のデータ用いて目的関数Φを定義した。なお，最終平均体積含水率θ ave は，マスバランスに
より水分保持曲線のθ の絶対値を決める情報であるため，必ず目的関数に入れる必要がある（Šimůnek et al.,1998）。
m

nj

j =1

i =1

*
Φ = ∑ v j ∑ ⎡⎣ p*j (ti ) − p j (ti ) ⎤⎦ + ⎡⎣θ ave− θ ave ⎤⎦
2

17

2

(2.9)

ここで，上付*の有無はそれぞれ実測値と計算値を示し，pj (ti)は観測地点 j の時刻 ti における土中水圧力，m は観測地点の数，nj は観測地点 j におけるデータの数，θ ave は実験終了時の平均体積含水率である。vj は異なる深さにおける土中水圧力の絶対値の違いを正規化するする重み係数であり，次式で定義した。
vj =1n jσ 2j(2.10)

ここで，σ j 2 は観測地点 j における土中水圧力の測定値の標本分散である。そして，非線形最小 2 乗法に基づく Levenberg-Marquardt 法（Marquardt, 1963）により，目的関数Φを最小とするパラメータの最適化を行なった。
逆解析によるパラメータの推定には，適切な初期値を与えることが不可欠である（Rassam et al., 2004）。そこで，図 2-1 の水分保持曲線の実測値に VG モデルを適合して得られたθ r，θ s，α，n のパラメータ値（θ r = 0.02，θ s = 0.574，α = 0.00235，n = 1.4）を初期値として用いた。また，飽和透水係数 Ks には変水位透水試験で得られた Ks = 5.9 cm

d−1，間隙結合係数 l には Mualem（1976）の推奨値 l = 0.5 を用いた。なお，初期値の組
み合わせ次第で，目的関数が他の極小値（local minimum）に収束することが非線形最
小 2 乗法の問題点として知られている。そのため，推定値が適当でないと考えられる場
合は，他の条件での収束値を初期値としてパラメータの再推定を行ない，推定結果を比較検討した。

2.5

結果と考察

2.5.1 パラメータの推定

図 2-2 は，試料長 L = 10 cm のシルト質土に対して可能蒸発速度 E = 1，0.5，0.1 cm

d−1，また L = 20 cm に対して E = 0.5 cm d−1 の条件を与えた実蒸発速度の時間変化である。L = 10 cm，E = 1 cm d−1 の条件では，蒸発速度が一定となる恒率蒸発期間は 1.5 日と短く，その後は蒸発速度が減少する減率蒸発期間となった。一方，可能蒸発速度が非常に小さい E = 0.1 cm d−1 のときは，深さ 0.5 cm の土中水圧力が−800 cm に達した段階（22 日）でも恒率蒸発が継続した。また，中間の E = 0.5 cm d−1 の場合，減率蒸発となるのは，L = 10 cm では 3.3 日，L = 20 cm では 5 日と，試料長が長いほど恒率蒸発期間は18

長いが，試料長が 2 倍になっても，恒率蒸発期間は 1.5 倍程度の増加であった。
図 2-3（a）は，図 2-2 に示した L = 10 cm，E = 1，0.1 cm d−1 の条件における土中水圧力分布を示す。図には，地表面の土中水圧力が−200 cm に達した時間と深さ 0.5 cm の土中
水圧力が−800 cm に達した時間として，E = 1 cm d−1 では 1 日と 2.6 日，E = 0.1 cm d−1 では 14.5 日と 22 日の分布を示した。恒率蒸発期間が 1.5 日と短い E = 1 cm d−1 の条件では，短期間に表面付近の土中水圧力が著しく低下して大きな圧力勾配が形成されるが，下層部は 2.6 日においても−200 cm 程度までしか低下しない。一方，22 日においても恒率蒸発が継続する E = 0.1 cm d−1 の条件では，蒸発初期から 14.5 日程度までカラム全体にほぼ一様な上方への水分フラックスが生じ，深さ方向に一様に土中水圧力が低下した。
また，小さな蒸発速度では，長時間にわたり下層からの水分供給により蒸発が進行する
ため，22 日において，試料下端の土中水圧力は−350 cm と，E = 1 cm d−1 の場合に比べて低下が大きい。同様に図 2-3（b）は，E = 0.5 cm d−1 を与えた L = 10 cm と 20 cm の土カ
ラムの土中水圧力分布を示す。地表面の土中水圧力が−200 cm に達した時間と深さ 0.5

cm の土中水圧力が−800 cm に達した時間として，L = 10 cm では 2.4 日と 4 日，L = 20cm では 3.7 日と 6.2 日の分布を示した。なお，L = 20 cm の分布は，上層 10 cm の土中
水圧力のみの表示である。L = 10 cm の 4 日における下端の土中水圧力は−200 cm であるのに対し，L = 20 cm の 6.2 日における深さ 10 cm の土中水圧力は−150 cm であり低下が小さい。試料長が 2 倍になった時の恒率蒸発期間の増加が 1.5 倍であることから，これは，L = 20 cm では相対的に下層部の乾燥の進行が抑えられた結果である。
Evaporation rate (cm/d)

1.2
L (cm) / E (cm/d)
cm / 1
10
cm / 0.1
10cm / 0.5
10
20
cm / 0.5

0.8

0.4

0
0

2

4

6

8

10

Time, t (day)
図 2-2 異なる試料長 L と可能蒸発速度 E におけるシルト質土の実蒸発速度。
Fig.2-2 Actural evaporation rates from a silty soil with different column length L and thepotential evaporation rate E.
19

0

1d

2.6d

Depth, -z (cm)

2

6.2d

(a)

(b)

L=10cm

E=0.5cm/d

22d

3.7d

14.5d

4

6
L=20cm
L=10cm

E=1cm/dE=0.1cm/d

8

2.4d
4d

10
-1000

-800

-600

-400

-200

-1000

-800

-600

-400

-200

0

Pressure head, h (cm)

図 2-3 異なる試料長 L と可能蒸発速度 E におけるシルト質土の土中水圧力分布：（a） L = 10
cm における E = 1，0.1 cm d−1，（b）E = 0.5 cm d−1 における L = 10, 20 cm。
Fig. 2-3 Pressure head profiles in a silty soil with different column length L and the potential
evaporation rate E: （a）E = 1，0.1 cm d−1 for L = 10 cm, （b）L = 10, 20 cm for E = 0.5 cm d−1.

2.5.2 応答曲面

VG モデルの 5 個のパラメータ x（α，Ks，n，θ r，θ s）と間隙結合係数 l の組み合わせ（x,l）について，蒸発過程の土中水圧力変化に対するパラメータの感度と，それぞれの組み合わせのパラメータ間の相関について応答曲面を用いて調べた。図 2-4 は，L = 10 cm，E = 1 cm d−1 の条件の深さ 0.5 cm の土中水圧力変化に対する 5 個の組み合わせ（x, l）
についての平方残差φ （（2.8）式）の等高線図（応答曲面）である。各パラメータは，それぞれ等間隔な 30 個の離散データとして与えた。そのため 1 つの応答曲面は，合計 900
個のパラメータ値の組み合わせである格子点から成る。パラメータがそれぞれシルト質
土の値 x0（θ r = 0.034，θ s = 0.46，α = 0.016 cm−1，n = 1.4, Ks = 6 cm d−1），l0 = 0.5 と等しいとき，（2.8）式の平方残差φ は 0 となり，応答曲面上で極小値をとる。図 2-4 では，l の範
囲を 0.2 < l < 0.8 とし，（x0，l0）が中央に位置し，また異なるパラメータの平均残差を相対的に比較するために，x の範囲は，グラフ左端中央部の（x0 – ∆x，l0）においてφ ＝ 105程度となる x0 – ∆x < x < x0 + ∆x とした。
平方残差φ は，異なるパラメータ値を持つ仮想土とシルト質土の土中水圧力が異なるほど大きい。そして，いずれかのパラメータを固定した断面は，もう一方のパラメータのφ20に対する感度を示す。たとえば，l0 = 0.5 の縦軸中央の水平断面のφ の分布は，他のパラメータ値をシルト質土の値に固定したときのパラメータ x の感度である。また，x = x0 の
横軸中央の垂直断面のφ の分布は l の感度である。φ が急な勾配を持つほど，パラメー
タの小さな変化が土中水圧力の計算結果に大きな違いを与えること，つまりパラメータの感度が高いことを示す。
さらに，等高線の広がりの傾きはパラメータ間の相関を示す。Ks，θ r は l と正の相関を持ち（図 2-4（b），（d）），またα，n，θ s は l と負の相関を持つ（図 2-4（a），（c），（e））。また，
等高線の分布が細長くなるほど相関が強い。φ の値の小さい白色の領域が細長いα と lの間には強い負の相関があり（図 2-4（a）），Ks と l の間には強い正の相関がある（図 2-4（b））。このような強い相関が存在するとき，（x, l）の異なる組み合わせが等しいφ の値を与えるため，φ が小さい図中の白色の領域において，パラメータ値の決定が難しい。さら

0.8

0.8

(a)

0.7

0.7

0.6

0.6

0.5

0.5

0.5

l

0.6

l

l

0.7

0.8

(b)

0.4

0.4

0.4

0.3

0.3

0.3

0.015

0.016

α (cm−1)

0.017

5

6

7

1.35

Ks (cmd−1)

0.8

0.8

(d)
0.7

0.7

0.6

0.6

0.5

0.5

(c)

1.36

1.37

1.38

1.39

n

(e)
500000100000
50000

l

l

10000

0.4

0.4

0.3

0.3

0

0.02

0.04

0.06

5000
1000
500
100

0.43

0.45

θr

0.47

0.49

θs

図 2-4 試料長 L = 10 cm，可能蒸発速度 E = 1 cm d−1 におけるシルト質土の深さ 0.5 cm の土中水圧力変化に対する平方残差φ の応答曲面：（a）l−α，（b）l−Ks，（c）l−n，（d）l−θ r，（e）l−θ s。
Fig. 2-4 Response surfaces of square residual φ for pressure heads at 0.5 cm depth for a silty
soil with column length L = 10 cm silt and the potential evaporation rate E = 1 cm d−1: (a)l-α , (b) l-Ks, (c) l-n, (d) l-θ r, (e) l-θ s.
21

に細かく応答曲面を見ると，図 2-4（a）のα − l の応答曲面では複数の極小値が現れている。複数の極小値の存在は，非線形最小 2 乗法によるパラメータ推定の際に，収束値がパラメータの初期値に依存する可能性が高いことを示す。図 2-4（a）では，l の推定を行なう場合，α の初期値に関して慎重な検討が必要であることを示している。
一方，l と n，θ r，θ s の相関は弱く，等高線中心部のφ の値の小さい白色の領域は，ほぼシルト質土の値（x0，l0）を中心にした小さな円状の分布となっている（図 2-4（c），（d），（e））。このように感度が高く，また相関が弱い場合，それぞれのパラメータが独立して収束する可能性が高く，パラメータの決定が容易である。Šimůnek et al.（1998）は，l を除く

VG モデルのパラメータ（α，Ks，n，θ r，θ s ）間の相関を応答曲面を用いて示し，θ s は他のパラメータとの相関が弱く，決定しやすいパラメータであること，また，n とθ r の間には強い正の相関があることを示した。また，Ks を含む応答曲面は，Ks 軸方向に長い等高線の分布となることを示し，Ks の初期値の重要性を述べている。

(1) 観測位置
図 2-5（a），（b）は，図 2-4（e）と同じ L = 10 cm，E = 1 cm d−1 の条件下で，観測位置 1，

9 cm における飽和体積含水率θ s と l の応答曲面である。比較のため，パラメータの範囲
は図 2-4 と等しくした。土中水圧力変化が大きい深さ 0.5 cm では（図 2-4（e）），深い位置の 1 cm，9 cm に比べて，l = 0.5，θ s = 0.46 の極小値付近においてφ の等高線が局所的に密な分布，すなわちφ が急な勾配を持つ。θ s = 0.46 における l の感度に注目すると，深さ 0.5 cm におけるφ < 1000 の範囲は 0.48 ≤ l ≤ 0.53，であり，深さ 1 cm の場合は 0.45

≤ l ≤ 0.58（図 2-5（a）），深さ 9 cm の場合はさらに大きくなり，等高線の中心部には複数の極小値が現れた（図 2-5（b））。そのため，φ が急な勾配を持つ表面に近い観測位置ほど lの感度が高く，収束しやすい条件である。
図 2-5（c），（d）は，同じ条件の観測位置 9 cm におけるα − l，Ks − l の応答曲面である。

θ s − l の応答曲面と同様に，観測位置が深い 9 cm では，φ が小さい領域が広がり，複数の極小値が現れた。図 2-2（a）に示したように，地表面に近い観測地点ほど乾燥が進行し，土中水圧力が大きく低下する。変化の大きい地表面に近い観測地点の土中水圧力を用いるほどパラメータの感度が増加するため，信頼度の高い推定に有効となる。しかし，現実に土中水圧力変化の非常に大きい地表面付近にテンシオメーターを設置して測定する場合，テンシオメーターの大きさ，設置位置の正確な確定等の実験上の制約から，
22

0.5〜1cm 程度の深さが，最も地表面に近い観測地点となると考えられる。
複数の深さの観測位置のデータを用いた場合の応答曲面は，（8）式に示すように，図2-4（e），図 2-5（a），（b）に示した観測地点 j の平方残差φ j を足し合わせたφ = Σφ j で表される。そのため，φ 分布の勾配は大きくなり，パラメータの感度が増加して収束のしやすい条件となる。これは，より多くの情報を与えると，逆解析の信頼度が向上することを意味する。そのため，一番浅い観測位置に加えて，複数の観測位置のデータを加えて逆解析を行なうことが望ましいと考えられる。

0.8

0.8

(b)9cm

0.7

0.7

0.6

0.6

0.5

0.5

l

l

(a)1cm

0.4

0.4

0.3

0.3
0.43

0.45

θs

0.47

0.49

0.43

θs

0.47

0.49

0.8

0.8

(c)9cm(d)9cm

0.7

0.7

0.6

0.6

0.5

0.5

l

l

0.45

0.4

0.4

0.3

0.3
0.015

0.016

α (cm−1)

500000100000
50000
10000
5000
1000
500
100
50
10

5

0.017

6

Ks (cmd−1)7

図 2-5 試料長 L = 10 cm，可能蒸発速度 E = 1 cm d−1 におけるシルト質土の土中水圧力変化に対する平方残差φ の応答曲面：（a）l−θ s（深さ 1 cm），（b）l−θ s（深さ 9 cm），（c）l−α（深さ 9 cm），（d）l−Ks（深さ 9 cm）。
Fig. 2-5 Response surfaces of square residual φ for pressure heads for a silty soil with column
length L = 10 cm silt and the potential evaporation rate E = 1 cm d−1: (a) l−θ s at 1 cm
depth, (b) l−θ s at 9 cm depth, (c) l−α at 9 cm depth, (d) l− Ks at 9 cm depth.

23

(2) 蒸発速度
図 2-6（a），（b）は，L = 10 cm，E = 0.1 cm d−1 の条件の観測位置 0.5 cm におけるθ s − lと Ks − l の応答曲面である。可能蒸発速度の大きい E = 1 cm d−1 の対応する応答曲面は，それぞれ図 2-4（e），（b）である。可能蒸発速度の小さい E = 0.1 cm d−1 では，恒率蒸発期が継続し（図 2-2），カラム全体の土中水圧力が低下した（図 2-3（a））。θ s − l の応答曲
面は，E = 1 cm d−1 の l に関して感度の高い等高線（図 2-4（e））とは対照的に，E = 0.1 cm

d−1 では，平方残差φ の小さい範囲が，θ s = 0.46 付近において l 軸方向に長く広がり，lに対して感度が非常に低いことを示した（図 2-6（a））。E = 1 cm d−1 の Ks − l の応答曲面
は，強い正の相関を示すものの，φ が小さい白色の領域は比較的小さい（図 2-4（b））。一方，図 2-6（b）の E = 0.1 cm d−1 では，φ が小さい白色の領域が右斜め方向に大きく広がり，
複数の極小値が現れた。このように Ks と l の間の正の相関が非常に強く，複数の極小値が存在するため，Ks と l の推定は困難である。
応答曲面を使った検討から，可能蒸発速度 E が大きく，短期間で恒率蒸発から減率蒸発へ移行するときの土中水圧力変化は，パラメータが収束しやすい条件であるため，推定される不飽和透水係数の信頼度は高いと考えられる。蒸発実験において，ファン等を使って蒸発をある程度促進することは，実験時間を短縮するのみならず，より信頼性

0.8

0.8

(a)

(b)

0.7

0.7

0.6

0.6

0.5

0.5

500000100000
50000

l

l

10000

0.4

0.4

0.3

0.3
0.43

0.45

θs

0.47

0.49

5000
1000
500
100

5

6

Ks (cmd−1)7

図 2-6 試料長 L = 10 cm，可能蒸発速度 E = 0.1 cm d−1 におけるシルト質土の深さ 0.5 cmの土中水圧力変化に対する平方残差φ の応答曲面：（a）l−θ s，（b）l−Ks。
Fig. 2-6 Response surfaces of square residual φ for pressure heads at 0.5 cm depth for a silty
soil with column length L = 10 cm silt and the potential evaporation rate E = 0.1 cm d−1: (a)l−θ s, (b) l−Ks.

24の高い不飽和透水係数，特にパラメータ l を推定するために必要である。ただし，あまり極端に蒸発速度が大きいと，測定期間の極端な減少や，また本研究では考慮していな
い表層における水蒸気移動の卓越など，推定精度を低下させる要因が生じることが予測される。

(3) 試料長
図 2-7 は，E = 0.5 cmd−1 に対する L = 10 cm，20 cm の深さ 0.5 cm の土中水圧力に対する応答曲面である。L = 10 cm の場合（図 2-7（a）），θ s = 0.46 における平方残差φ <10000 の範囲は 0.38 ≤ l ≤ 0.62 であるのに対し，L = 20 cm の場合（図 2-7（b）），0.42 ≤ l ≤0.60 と狭くなり，応答曲面は比較的急なφ の勾配を示す分布となった。これは，試料長 Lが長いほど，l に関する感度が高いことを示す。E = 0.5 cm d−1 の場合，E = 1 cm d−1 に比べて恒率蒸発期間が長く（図 2-2），地表面付近の土中水圧力低下も遅い。しかし，試料
長 L を長くすると，試料下層部の土中水圧力が高く保たれる分布となり，地表面と下層部の圧力差が大きくなる（図 2-3（b））。このような条件では，l の感度が相対的に増加して収束が容易になると考えられる。しかし，前節の蒸発速度の影響に比べて応答曲面の変化が小さいことから，試料長を長くするよりは蒸発を促進する方が，不飽和透水係数の推定精度の向上には有効と考えられる。
0.8

0.8
0.7

(a)
L=10cm0.7(b)
L=20cm

500000100000

0.6

0.6

0.5

0.5

50000

l

l

10000

0.4

0.4

0.3

0.3

0.43

0.45

θs

0.47

0.49

5000
1000
500
100

0.43

0.45

θs

0.47

0.49

図 2-7 可能蒸発速度 E = 0.5 cm d−1 を与えた異なる試料長 L のシルト質土における深さ
0.5 cm の土中水圧力変化に対する平方残差φ の l−θ s 応答曲面：（a）L = 10 cm，（b）L =20 cm。
Fig. 2-7 Response surfaces of square residual φ for pressure heads at 0.5 cm depth for a silty
soil with different column lengths L and the potential evaporation rate E = 0.5 cm d−1: (a) L= 10 cm, (b) L = 20 cm.
25

2.5.3 蒸発法によるパラメータの推定(1) θ r が不飽和透水係数に及ぼす影響
本研究で用いた藤森シルトは，図 2-1 に示すように比較的明瞭な空気侵入圧を持ち，またその値が−100 cm と小さい特徴を持つ。VG モデルを図 2-1 の水分保持曲線に適用
すると，勾配がなだらかな形状を持つため，小さな n の値（ = 1.4）が得られる。ところが，
小さな n の値を持つ VG モデルは，飽和近傍の不飽和透水係数が非現実的に大きく変化する問題が知られている（Rassam et al., 2004）。一方，VG モデルの n はθ r と正の相関を持つ（Šimůnek et al., 1998）。このため，土中水圧力の測定範囲（−800 < h < 0 cm）において，大きな n に対してθ r を増加させることにより，異なるパラメータ値の組み合わせに対してほぼ等しい水分保持曲線を表現することが可能である。
図 2-8（a）の破線は，図 2-1 の実測値に適合させた VG モデル（実線）に対して，他のパラメータは変化させずに，n = 2 に対してθ r = 0.25 と増加させたときの水分保持曲線である。点線の水分保持曲線が空気侵入圧付近でわずかにθ が大きい点を除き，両者はほぼ同一な水分保持曲線を示す。一方，図 2-8（b）に示す飽和付近の不飽和透水係数

0.6

101

Water content, θ (cm3cm-3)

10-1

0.4
101

10-3

0.2

100

10-5

n=1.4, θr =0.01
n=2, θr =0.25
0
100

0.4

101

102

103

104

Pressure head, h (-cm)

0.3

10-1
0.6

0.5

0.4

0.5

Water content, θ (cm3cm-3)

10-7
0.6

図 2-8 異なるθ r と n の組み合わせによる VG モデル（θ s = 0.574，α = 0.00235，Ks = 5 cm
d−1，l = 0.5）による（a）水分保持曲線，（b）不飽和透水係数と体積含水率の関係。
Fig. 2-8 (a) Volumetric water content vs. pressure head and (b) unsaturated hydraulicconductivity vs. volumetric water content for different sets of parameter values of θ r and n
for the VG model (θ s = 0.574, α = 0.00235, Ks = 5 cm d−1, l = 0.5).

26

Hydraulic conductivity, K (cmd-1)

(b)

(a)

は，n = 1.4 では飽和付近での極端な低下が生じるのに対し，n = 2 ではこの問題が回避される。また，このようにθ r は不飽和透水係数の形状に影響を及ぼすパラメータであるため，低水分領域の水分移動に関わらない水分量としての物理的な意味よりも，推定パラメータとしての性格が強いことも指摘されている（小杉，2007）。藤森シルトの蒸発実験の圧力測定範囲の−800 < h < 0 cm においては，図 2-8（a）の実線と破線の水分保持曲線が十分に良く一致すること考慮して，本研究ではθ r を推定パラメータとして用いた。

(2) 試料長がパラメータの推定に与える影響
本研究の計算に用いた HYDRUS では，逆解析の際に得られる推定パラメータの共分
散行列中の対角成分である自己相関成分に基づいた標準偏差σ と，σ を用いた 95 ％の信頼区間が結果と併せて与えられる。この信頼区間は，パラメータの推定精度の指標であり，2.5.2 節で示したパラメータの感度と相関のみならず，目的関数に与えるデータの数が不十分な場合や，誤差が大きく不適切である場合にも，信頼区間は大きくなる。
しかし，応答曲面で示したような収束しやすい条件ほど信頼区間は小さく，また逆に収束が困難な条件ほど信頼区間は大きい傾向は維持されると考えられる。そのため，以下，蒸発法で得られたパラメータ値に関して，得られた信頼区間と 2.5.2 節の応答曲面の解
析結果を用いて考察する。
表 2-1 は，2 種類の試料長の蒸発実験に対するパラメータの推定値と標準偏差である。

8.3 cm カラム（L = 8.3 cm）では，深さ 0.7，3.1，5.2 cm の 3 点で測定した土中水圧力を（2.9）式の目的関数に用い，16 cm カラム（L = 16 cm）では，深さ 1，5，13 cm の 3 点の組
み合わせを目的関数に用いた。
図 2-9 は，推定したパラメータ値に基づく水分保持曲線である。どちらの試料長も，すべての深さの土中水圧力データを用いた場合の結果である。8.3 cm カラム，16 cm カラムから推定される水分保持曲線は，−200 < h < −800 cm の範囲でほぼ同一の形状を示した。図 2-1 の水分保持曲線に対しては n = 1.4 とθ r = 0.02 が適合値であったが，蒸発法の推定では，比較的大きな n とθ r が得られた（8.3 cm カラムの場合で n = 1.90，θ r =0.222，16 cm カラムの場合で n = 3.33，θ r = 0.367）。これは，2.5.3（1）節で述べたように，
n が小さい場合，飽和付近における非現実的な不飽和透水係数の低下が生じ，測定データとの適合ができないためである。本解析では，θ r を推定パラメータとして用いたため，いずれの試料長においても，低圧力領域の水分保持曲線の実測値を大きく過大評価し27

表 2-1 蒸発法による VG モデルのパラメータの推定値と標準偏差（下段）。
Table 2-1 Optimized parameter values and standard deviations (lower row) for the VG model
based on the evaporation method.
Length tensiometers
L =16cm
1,5,13cm
1cm
5cm
13cm
1,5cm
1,13cm5,13cm
L =8.3cm 0.7,3.1,5.2cm

θrθsα (cm−1)

0.367

0.530

0.00232

0.006

0.008

0.389

0.528

0.005

0.008

2.14×10−5

0.398

0.532

0.003

2.39×10−5n
K s (cmd−1)
3.33
1.26

l
8.12×10−5

0.07

0.13

1.34×10−4

3.82

15.9

0.002

0.08

2.47

0.013

0.002512.68

0.00980.004

1.87×10−5

4.13
0.10

0.62

0.027

0.409

0.539

0.002614.00

0.0089

0.008

0.011

3.09×10−5

4.31
0.14

6.38

0.037

0.372

0.530

0.00234

3.58

1.69

3.56×10−4

0.005

0.007

2.10×10

0.07

0.26

3.47×10−3

0.349

0.532

0.00223.03

0.84

1.31×10−6

0.010

0.010

3.58×10

0.08

0.09

0.401

0.533

0.00255

4.17

3.33

0.020

0.003

1.67×10

0.09

0.36

2.80×10−4

0.222

0.588

1.90

1.73

0.050

0.019

0.008

0.05

1.01

0.960

0.00236

−5

−5

−5

0.00302

8.04×10−5

1.75×10−4

3.32×10−4

Water content, θ (cm3cm-3)

0.6

0.5

0.4

0.3

observed
estimated
L=16cm
L=8.3cm

0.2
10

100

1000

Pressure head, h (-cm)
図 2-9 蒸発法で推定した藤森シルトの水分保持曲線：プロットは実測値，実線は 16 cm カラムによる推定値，点線は 8.3 cm カラムによる推定値。
Fig. 2-9 Water retention curves for a Fujinomori silty soil estimated with the evaporation
method: plots are observed data, the solid line is for L = 16 cm, and dashed line is for L =8.3 cm.

28

た結果となっている。そのため，今回の推定は，あくまで蒸発法の適用領域（h > −800cm）に対する結果であることに注意する必要がある（Šimůnek et al.,1998）。
また，試料長により飽和体積含水率θ s の推定値に違いが見られた（8.3 cm カラムでθ s= 0.588，16 cm カラムでθ s = 0.530）。試料長が長い 16 cm カラムではθ s が小さいのは，試料が長いほど毛管飽和による水分飽和度が低いためと考えられる。この飽和度の違
いは，飽和透水係数 Ks の推定値にも現れた（8.3 cm カラムで Ks = 1.73 cm d−1，16 cm カラムで Ks = 1.26 cm d−1）。なお，変水位透水試験による飽和透水係数は Ks = 5.9 cm d−1である。
通常，水分保持曲線の測定データを目的関数に入れると，より信頼度の高い水分保持曲線を推定でき，また計算の収束も速くなることが期待できる（Šimůnek et al.,1998）。しかし，測定データの飽和度が異なる場合，Ks や飽和付近の水分保持曲線の実測データを目的関数に入れても推定の信頼度は向上しない点には注意が必要であると考えられる。今回の結果は，別途測定したθ s，Ks に関する実測データがある場合においても，実測値と推定値との比較を行ない，妥当な値について慎重に検討する必要があることを示している。

8.3 cm カラム，16 cm カラムによる推定において，水分保持曲線を決めるパラメータθ r，θ s，α，n の標準偏差は十分に小さく（表 1），試料長 L の影響は小さい（Hopmans et al.，2002）。一方で，不飽和透水係数を決めるパラメータ Ks と l の標準偏差は，8.3 cm カラム
ではそれぞれ 1.01，0.960 と大きな値を持つのに対し，16 cm カラムではそれぞれ 0.13，

1.34×10−4 と小さな値となった。推定値 le に対して le − 2σ ≤ l ≤ le + 2σ と定義される 95％の信頼区間は，8.3 cm カラムでは−1.87 ≤ l ≤ 1.97，16 cm カラムでは−1.86 ×10−4 ≤ l ≤3.49 ×10−4 であった。図 2-10 には，推定した不飽和透水係数に加えて，l の 95％の信頼区間として与えられる最大値と最小値を用いた計算値を示す。8.3 cm カラムの不飽和透水係数の 95％信頼区間に対応する範囲は，たとえば h = −800 cm において，2 ×10−3 ≤K ≤ 5 ×10−2 cm d−1 と広い。一方，16 cm カラムでは図上で確認できないほどこの信頼区間の幅が非常に狭い。これらの結果は，試料長 L が長いと不飽和透水係数，特にパラメータ l の推定精度を向上させることを示しており，数値実験による応答曲面（2.5.2（3）節）で得られた結果とも一致する。

29

Hydraulic conductivity, K (cmd-1)

101

100

10-1

10-2

10-3

L=16cm
L=8.3cm

10-4
10

100

1000

Pressure head, h (-cm)
図 2-10 蒸発法で推定した藤森シルトの不飽和透水係数と土中水圧力の関係。実線は推定値，点線は間隙結合係数 l の 95 %信頼区間に対応した信頼範囲。
Fig. 2-10 Unsaturated hydraulic conductivity as a function of soil water pressure head for a
Fujinomori silty soil estimated with the evaporation method: solid lines are optimized
conductivity functions and dashed lines represent the confidence ranges based on the 95%
confidence limits for the pore-connecting parameter l.

(3) 観測位置がパラメータの推定に与える影響
数値実験による感度解析により，パラメータの推定では，表層に近い複数の深さの土
中水圧力の測定値を用いることが望ましいことを示した。そこで，蒸発実験の土中水圧力測定データを用いて，異なる深さのデータを組み合わせ，パラメータの推定について，
表 1 の推定結果と標準偏差σ に基づき検討した。
深さ 1 cm のみ，5 cm のみ，13 cm のみの場合，θ r，θ s，α の 推定値は良く一致し，また，それぞれのσ も十分に小さい値が得られた。一方 n，l については，地表面に近い観測点ほどσ が小さくなり，信頼度が高い。これは，前述の 2.5.2（1）節の応答曲面で得られた収束条件の考察と一致する。

2 地点の土中水圧力を組み合わせてパラメータを推定した場合，l の標準偏差はσ =10−4~10−3 と極めて小さくなった。深さ 1 cm のみではσ = 0.013 であったことから，2 地点の土中水圧力を用いることは，l の推定の信頼度を大きく向上させている。また，Ks の標
準偏差もσ = 0.09 ~ 0.36 であり，深さ 1 cm のみ（σ = 2.47）と比較しても十分に小さくなっ30

た。しかし， 3 地点の土中水圧力を用いた場合， l と Ks の標準偏差はそれぞれ σ =1.34×10−4，0.13 であり，2 地点の場合に比べての改善は見られなかった。以上より，地表面付近を含めた 2 地点の土中水圧力データを用いると，不飽和透水係数の推定における信頼度を向上できると考えられる。
なお，前述の通り，パラメータの推定には，原則として図 2-1 の適合値と Ks = 5.9 cm d−1，

l = 0.5 を初期値に用いたが，深さ 1 cm のみの場合，推定結果はθ r = 0.001，θ s = 0.515，

α = 0.00199，n = 2.15，Ks = 0.004 cm d−1，l = 0.001 となり，他の推定と極端に異なる値，
特に Ks は非現実的に小さな値を示した。このとき，得られた n と Ks の相関係数は−0.8132 であり，強い負の相関を示した。図 2-4（a），（b）に示したように，パラメータ間の相関が強い場合，真の値とは異なるパラメータの組み合わせが推定される可能性がある。
そこで，深さ 1 cm のみの推定は，表 1 に示すすべての深さを用いた（1，5，13 cm）推定
結果を初期値として，再度推定を行なった。推定結果が他の深さと大きく異なるときや非現実的なときには，異なる初期値を用いた検討が必要である。

2.6

おわりに

本章では，蒸発法による水分移動特性の推定における適切な実験条件や推定方法について，不飽和透水係数の推定に焦点を当てて検討した。まず，水分飽和したシルト
質土カラムからの蒸発過程を対象とした数値実験を行ない，水分移動特性モデルの van

Genuchten-Mualem（VG）モデルのパラメータについて，不飽和透水係数を決める間隙
結合係数 l と他のパラメータの組み合わせによる土中水圧力に関する 2 次元応答曲面を示した。そして，それぞれのパラメータの感度やパラメータ間の相関に基づき，蒸発法の最適条件について検討した。その上で藤森シルトを用いた蒸発実験を行ない，測定した
土中水圧力を用いた逆解析により水分移動特性パラメータの推定を行なった。そして，実験条件の違いが蒸発法による水分移動特性パラメータ，特に l の推定に与える影響についての検討を行なった。
土中水圧力変化に対する 2 次元の応答曲面は，l とα の強い負の相関，また l と Ks の強い正の相関を示した。また応答曲面上で複数の極小値を示すことから，異なるパラメ
ータ値の組み合わせに収束する可能性が示された。そのため，l，α，Ks の推定に関して
31

は，複数の初期値の組み合わせについての推定結果の比較が必要である。一方，n，θ s，

θ r と l との間の相関は弱く，パラメータの決定がα，Ks に比べて容易であることが示された。
シルト質土からの水分蒸発の数値実験において，間隙結合係数 l と飽和体積含水率

θ s が異なる観測位置における土中水圧力変化に及ぼす感度について応答曲面をもちいて検討したところ，l の感度は地表面に近い位置の土中水圧力に対してほど高かった。
これは，土中水圧力の時間変化が大きい表層のデータほど収束しやすく，信頼度の高い推定が可能なことを意味する。さらに，複数の観測位置のデータを加えて逆解析を行なうことは，収束の条件が改善されるので望ましいと考えられた。
同様に，蒸発速度と試料長についても検討した。蒸発速度を促進して恒率蒸発期間を短くすると，大きな圧力勾配が試料内に生じた。また，試料長を長くすると，下層土の土中水圧力が高く保たれ，試料内の土中水圧力差が増加した。地表面付近で土中水圧力が低下する分布が形成されると，応答曲面における不飽和水分移動特性パラメータの感度が高くなり，パラメータの収束が容易になる。特に蒸発をある程度促進することは，
実験時間の短縮のみならず，推定した不飽和透水係数の信頼性を高めるために望まし
いと考えられた。
藤森シルトを用いた蒸発法では，異なる試料長と観測点の位置と数について検討した。
異なる試料長の推定したパラメータの標準偏差σ を比較したところ，特に Ks と l についての違いが大きく，試料長が長いカラムのσ は小さかった。また，1 地点の観測点の土中水圧力を用いて推定を行なった場合，θ r，θ s，α の 結果に違いは見られなかったが，n，lについては表面に近い観測点ほど，σ が小さくなった。こうしたσ の結果は，応答曲面による収束条件の検討結果と一致している。また，異なる 2 地点の土中水圧力データを組
み合わせてパラメータを推定した場合，1 地点のみの場合に比べて，Ks と l の推定の信頼度が大きく向上した。以上より，感度の高い地表面付近を含めた 2 地点の土中水圧力
データを用いることは，不飽和透水係数の推定において望ましいと考えられた。

32

第3章

3.1

砂質土と団粒土に対する水分移動特性関数の検討

はじめに

土の水分移動特性は，体積含水率θ と土中水圧力 h（以下，圧力）の関係である水分保持曲線と，θ や h の関数である不飽和透水係数で与えられる。リチャーズ式に初期条
件，境界条件を与え，数値計算によって土中の不飽和水分移動を予測する場合，水分保持曲線と不飽和透水係数をそれぞれ連続的な関数で表現する必要がある。この水分保持曲線と不飽和透水係数の両者を数理モデルにより表現した関数は，水分移動特性関数（soil water hydraulic function）とよばれている（Jury and Horton, 2006）。
水分保持曲線については，多くのモデルが提案されている。広く知られたモデルである Brooks and Corey（1964）の水分保持曲線関数は，明確な空気侵入圧を示す。また，

van Genuchten（1980）のモデルは，S 字型のなめらかな形状を示すため，実測値に対する適合性が高く，最も広く用いられている。Kosugi（1996）は，対数正規分布で表現した間隙径分布の毛管保水に基づく物理的なモデルを示している。
水分保持曲線が比較的容易に測定できるのに対し，不飽和透水係数の測定は，多くの労力と時間を要する。特に，低水分領域の不飽和透水係数の測定は非常に困難である。そこで，古くから水分保持曲線に基づいて不飽和透水係数を推定する試みが行なわれてきた。その多くは，土の保水性を土の間隙径分布と毛管保水で表現し，さらに毛管内のポワズイユ流れを仮定することにより不飽和透水係数を推定する間隙径分布モ
デルである（Childs and Collis-George, 1950; Burdine, 1953; Mualem, 1976; Jury andHorton, 2006）。このような水分保持曲線‐不飽和透水係数連結モデルは，不飽和透水係数の情報が限られている場合においても物理的根拠に基づいた不飽和透水係数の推定を可能にする上，水分保持曲線関数と不飽和透水係数関数を別個に与える場合に比べて少ないパラメータ数で土の水分移動特性関数を表現できるため，有用性が高い（小杉，2007）。特に，近年広く用いられているマルチステップ法や蒸発法といった水分
移動特性関数の逆解析手法においても，少ないパラメータ数で水分移動特性を表現することは，解の収束性という点で重要である（Eching and Hopmans，1993; Šimůnek et al,

1998; 坂井・取出，2007）。
Mualem（1976）の提案した間隙径分布モデルは，モデル中の積分から解析解が得や33

すいように数学的利便性を求めて工夫された関数型で表現されている（小杉，2007）。

van Genuchten（1980）は，水分保持曲線関数を Mualem モデルに代入し，不飽和透水係
数関数を解析的な閉形式解（closed-form）により表現した。この van Genuchten の導いた
水分保持曲線-不飽和透水係数連結モデル（以下 VG モデル）は，他のモデルに比べて単純な関数であり，また，数値計算に相性のよいなめらかな関数であるため，最も広く用いられている（小杉，2007）。
しかし，VG モデルの不飽和透水係数に関してもいくつかの問題点が指摘されており，修正モデルが提案されている。VG モデルの水分保持曲線関数は，残留体積含水率θ rを定義し，圧力の低下に伴いθ r に漸近する関数型で表現される。θ r は水分移動に関与しない不動水と仮定しているため，VG モデルはθ r に近づいた低水分領域において，不飽和透水係数を過小評価することが指摘されている（Tuller and Or, 2001; 小杉，2007）。
特に，高い土中水圧力で水分が急激に低下してθ r に漸近する砂質土では，VG モデル
が不飽和透水係数を正しく評価できる圧力範囲は非常に狭い。 Fayer and Simmons（1995）は，乾燥領域においてθ r に漸近するのではなく，圧力が−107 cm のときに水分量が 0 となるように，VG モデルの残留体積含水率θ r を修正した水分移動特性関数を示した。
また，水分保持曲線が S 字型の形状を示す VG モデルでは，保水性と透水性を適切に表現できない土も存在する。たとえば，団粒構造が発達した土では，団粒間間隙と団粒内間隙における水分保持特性が異なるため，2 段の階段状の水分保持曲線（徳本ら，2005）を示すことが多い。Durner（1994）は，団粒土のように異なる土性から構成される土を対象に，複数の VG モデルを足し合わせることにより，階段状を示す水分保持曲線関数を表現した。

Fayer and Simmons のモデル，および Durner のモデルについても，Mualem モデルを適用した水分保持曲線‐不飽和透水係数連結モデルが示されている（以下，それぞれの連結モデルを Fayer モデル，Durner モデルとよぶ）。しかし，比較的最近提案されたため適用事例が少ない Fayer モデルと Durner モデルに関しては，不飽和透水係数の検討がまだ十分に行なわれていない。また，実際の土中の水分移動現象を対象に Fayer モ
デルや Durner モデル適用して数値計算を行なった事例も少ない。
そこで本研究では，砂質土と団粒構造の発達した黒ボク土に対して水分蒸発実験を行ない，蒸発法により Fayer モデル，Durner モデルのパラメータを推定した。そして，水34

分保持曲線と蒸発過程の圧力変化に対して，実測値と逆解析で得たパラメータによる計算値を比較することにより，それぞれのモデルの不飽和透水係数関数の砂質土および団粒土への適用性について検討した。

3.2

水分移動特性関数

3.2.1 Mualem モデル

Mualem（1976）は，不飽和透水係数関数を推定する間隙径分布モデルとして次式を与えた。

⎡ Se 1
K ( Se ) = K s K r = K s Se ⎢ ∫
dSe
0 h S
(
)
e
⎣

⎤
1
∫0 h ( Se )dSe ⎥⎦
1

l

2

(3.1)

ここで，K は不飽和透水係数（L T−1），Ks は飽和透水係数（L T−1），Kr は Ks に対する比透水係数（−），Se は有効飽和度（−），l は間隙結合係数（pore-connectivity coefficient）（−），

h は土中水圧力（L）である。Mualem（1976）は，45 種類もの土試料に対して適用し，その平均値として l = 0.5 を示している。（3.1）式の Mualem モデルに対して水分保持曲線関数 h(Se)を与えると，K を h あるいは Se の関数として表すことができる。この水分保持曲線‐不飽和透水係数連結モデルにより，保水性と透水性を同じパラメータを用いた関数で表現することができる。

3.2.2 水分保持曲線・不飽和透水係数連結モデル(1) van Genuchten モデル
現在，Mualem モデルに基づく水分保持曲線−不飽和透水係数連結モデルの中で最も広く利用されているものが，van Genuchten（1980）が提案した van Genuchten（VG）モデルである。VG モデルは，水分保持曲線関数を次式で与える。

()

θ − θrn −m= Se = 1 + α hθs − θr(3.2)

ここで，θ は体積含水率（L3 L−3），θ r は残留体積含水率（L3 L−3），θ s は飽和体積含水率（L3 L−3），α（L−1），n（−），m（= 1 − 1/n）は水分保持曲線の形状を与えるパラメータであ35

る。（3.2）式を（3.1）式の Mualem モデルに代入すると，次の不飽和透水係数関数の閉形式解が得られる。
m
K ( Se ) = K s Sel ⎡⎢1 − (1 − Se1 m ) ⎤⎥
⎣
⎦

2

(3.3)

(2) Fayer モデル

VG モデルでは，残留体積含水率θ r として表現される水分量が水分移動に関与しな
いと仮定するため，低水分領域の不飽和透水係数を過小評価する（Tuller and Or, 2001;
小杉，2007）。Fayer and Simmons（1995）は，土粒子表面の吸着保水を考慮し（Campbell

and Shiozawa, 1992），炉乾に対応する土中水圧力で水分量が 0 となるように VG モデルのθ r を修正したモデルを提案した。

θ θ s = Se = (1 − χ ( h )θ a θ s ) ⎡1 + α h ⎤
n

⎣

−m

⎦

+ χ ( h )θ a θ s(3.4)

ここで，θ a は水分保持曲線の形状を与えるパラメータであり，χ (h) = 1 – ln (−h)/ ln (−hm)
である。hm はθ = 0 となる炉乾に対応する土中水圧力を示す。Rossi and Nimmo（1994）は，

105 ~ 110 °C で相対湿度 50 %の乾燥炉に対応する圧力として hm = −107 cm としている
が，本研究では鳥取砂丘砂の実測値に対する適合性を考慮して，hm = −106 cm を用いた（後述，3.4.1 節）。さらに（3.4）式を書き直すと次式が得られる。

θ θ s = Se = ⎡1 + α h ⎤
n

⎣

−m

⎦

{

n
+ 1 − ⎡1 + α h ⎤⎣⎦−m} χ ( )θ θ
h

a

s

(3.5)

図 3-1（a）に，後述の表 3-1 に示す鳥取砂丘砂のパラメータを用いた（3.5）式の水分保持曲線を示す。また図には，（3.5）式右辺の第 1 項と第 2 項についても示した。右辺第 1 項は，θ r = 0 とした VG モデル（（3.2）式）であり，高土中水圧力領域において排水する間隙
保水に相当する成分を表す。一方，右辺第 2 項は，h < −40 cm の低圧力領域において hの対数軸に対して Se が直線的に減少する形状を表現している。すなわち，第 1 項の VG
モデルで不動水と仮定していた低土中水圧力領域の保水成分を補正する項である。

Fayer and Simmons の水分保持曲線関数（（3.4）式）を Mualem モデル（（3.1）式）に代入すると次の不飽和透水係数関数が得られる。

K ( Se ) = K s Sel ⎡⎣ Γ s ( h ) Γ max ⎤⎦
36

2

(3.6)

1.2

Effective saturation, Se (-)(a) Fayer model(b) Durner model1
1st term
2nd term
total

0.8
0.6
0.4
0.2
0 0
10

101

102

103

104

105

106 100

101

102

103

104

105

106

Pressure head, h (-cm)

図 3-1 （a）Fayer モデル（（3.4）式）と（b）Durner モデル（（3.6）式）の水分保持曲線関数が表す 2
つの保水形態：点線は右辺第 1 項，一点鎖線は右辺第 2 項，実線は 2 項の和（左辺）。
Fig. 3-1 Two different regions for water retention of (a) the Fayer model (eq. (3.4)) and (b) the
Durner model (eq. (3.6)): dashed and chain lines indicate the first and second terms of theright-hand side of the equations, respectively, and solid lines indicate the sum of the two terms(left-hand side).

ここで，Γs(h)とΓmax の定義は付録 1 に示す。以下，（3.4）式と（3.6）式で示される水分保持
曲線‐不飽和透水係数連結モデルを Fayer モデルとする。

(3) Durner モデル

VG モデルの水分保持曲線は S 字型の形状を示すため，団粒構造の発達した土にみられる 2 段の階段状の水分保持曲線は表現できない。Durner（1994）は，複数の VG モ
デルを線形結合することで，多段の水分保持曲線を表現した。本研究では，団粒の内外の間隙における 2 つの保水形態を想定し，2 つの VG モデルを結合したモデルを対象とする。このとき水分保持曲線関数は，次式で表される。

θ − θrn −mn −m= Se = w1 ⎡1 + α h ⎤ + w2 ⎡1 + α 2 h ⎤⎣⎦⎣⎦θs − θr
2

2

(3.7)

ここで，w は重み係数で，w1 + w2 = 1 である。図 3-1（b）には，後述の表 3-1 に示す熊本
黒ボク土のパラメータを用いた（3.7）式の水分保持曲線を示す。また図には，図 3-1（a）の Fayer モデルと同様に，（3.7）式右辺の第 1 項と第 2 項についても示した。右辺第 137

項は，団粒間の大きな隙間における間隙保水を表し，およそ h = −1000 cm までにすべての水分が排水する。一方，右辺第 2 項は，小さな隙間の団粒内保水を表し，空気侵入
圧が h = −1000 cm 程度の保水特性を示す。（3.5）式の Fayer モデルにおいては，圧力範囲によってそれぞれの項が卓越するのに対して，Durner の水分保持曲線関数は，飽和の水分量を 2 分割し，2 つの異なる保水特性の足し合わせとして表現するため，2 つの項の水分保持曲線に対する役割はそれぞれのモデルで異なる。しかし，両モデルとも異なる 2 つの保水形態の和として土の保水性を表現したモデルである点は共通である。さらに，パラメータ数の多い関数である Durner のモデルは様々な水分保持曲線の形状を表すことが可能であり，Fayer モデルの低水分領域を直線で表現した形状を表すこともできる。

Durner の水分保持曲線関数（（3.7）式）を Mualem モデル（（3.1）式）に代入すると，次の不飽和透水係数関数の閉形式解が得られる（Priesack and Durner, 2006）。

(

)

m2
m
l
( w1Se1 + w2 Se 2 ) ⎛⎜ w1α ⎡⎣⎢1 − (1 − Se11m ) ⎤⎦⎥ + w2α 2 ⎡⎢⎣1 − 1 − Se12m2 ⎤⎥⎦ ⎞⎟⎝⎠K ( Se ) = K s2( w1α + w2α 2 )

2

(3.8)
以下，（ 3.7 ）式と（ 3.8 ）式で示される水分保持曲線‐不飽和透水係数連結モデルを

Durner モデルとする。

3.3

蒸発法

3.3.1 蒸発実験
試料には鳥取砂丘砂と熊本の黒ボク土を用いた。図 3-2 のプロットは，吸引法，加圧板法，蒸気圧法で測定した水分保持曲線である。砂丘砂の水分保持曲線は，−20 ≥ h ≥−50 cm で水分量が急激に低下し，h < −50 cm では緩やかに低下する特徴を示す。低土
中水圧力領域においてやや階段状の形状が見られるが，この主な原因は，h ≥ −1000

cm は加圧板法（中野ら，1995）で測定し，h < −1000 cm は蒸気圧法（藤井・中野，1984）
で測定したため，測定時に完全な平衡状態でなかった可能性も含めて測定法の違いに起因すると考えられる。一方，黒ボク土は明瞭な 2 段の水分保持曲線を示す。黒ボク土
においては，h < −150 cm の範囲はすべて加圧板法で測定した。
38

0.8

0.4

(b) Andisol

Water content, θ (cm3cm-3)(a) dune sandobserved
Fayer model
Durner model
VG model

0.3

0.6

0.2

0.4

0.1

0.2

0
100

101

102

103

104

105

106

0
100

101

102

103

104

105

106

Pressure head, h (-cm)

図 3-2 （a）鳥取砂丘砂，（b）熊本黒ボク土の水分保持曲線：プロットは実測値，実線，点線，一点鎖線はそれぞれ Fayer モデル，Durner モデル，VG モデルによる適合曲線。
Fig. 3-2 Water retention curves for (a) Tottori dune sand and (b) Kumamoto Andisol: open
circles are the observed data, and solid, dashed, and chain lines describe the fitted curves
described with the Fayer, Durner, and VG model, respectively.

蒸発実験では，砂丘砂（乾燥密度ρb = 1.62 g cm−3）と黒ボク土（ρb = 0.48 g cm−3）を，それぞれ内径 3.8 cm で高さ 8 cm，内径 7.8 cm で高さ 18.6 cm のアクリルカラムに充填し，マリオット管を用いて地下水位を徐々に上げ，3 日間かけて毛管飽和した。毛管飽和後，マリオット管を切り離して試料下端を閉鎖し，試料表面に小型ファンを用いて一定の風をあてながら蒸発を促進させた。砂丘砂には深さ 0.5，1，2，3，5 cm，黒ボク土には深さ 5.5，

10.5，15.5 cm の位置にテンシオメーターを水平方向に挿入し，圧力変換器を用いて土中水圧力を測定した。同時に，カラム全体を電子天秤の上に設置し，重量変化から積算蒸発量を求めた。圧力変換器，電子天秤は Campbell 社のデータロガーCR1000 に接続
し，試料内の土中水圧力，積算蒸発量を自動計測した。試料の上部が乾燥し，最上部のテンシオメーターが測定不能になった段階で，実験を終了し，カラムを解体して炉乾法で試料の平均体積含水率θ ave を測定した。なお，黒ボク土の表層付近の圧力は，テンシオメーターの不良により，今回は解析の対象から外した。

39

3.3.2 逆解析(1) 水分移動式
水分飽和した長さ L の土カラムの地表面からの水分蒸発過程を対象に，以下の水分移動式を解析に用いた。鉛直 1 次元で水蒸気移動を考慮しない不飽和土中の水分移動は，下記のリチャーズ式で表される。

∂θ ∂ ⎡⎛ ∂h ⎞ ⎤= ⎢ K ( h) ⎜ + 1⎟ ⎥∂t ∂z ⎣⎝ ∂z ⎠ ⎦(3.9)

ここで，t は時間（T），z は上向き正の空間座標（L）であり，地表面で 0 である。また，水蒸
気移動を含むリチャーズ式は次式で表される（Philip and de Vries, 1957; 坂井・取出，

2006）。
∂θ ∂ ⎡∂h ⎤⎛ ∂h ⎞= ⎢ K (h) ⎜ + 1⎟ + K vh ⎥∂t ∂z ⎣∂z ⎦⎝ ∂z ⎠(3.10)

Kvh は圧力勾配による水蒸気透過係数（L T−1）であり，詳細は付録 2 に示す。なお本解析
では，蒸発実験における温度は一定とみなし，温度勾配による水蒸気フラックスは無視した。
数値計算の初期条件は，圧力分布 hi（z）を与えた。
h ( z , 0 ) = hi ( z )(3.11)

カラム下端は水分フラックスが生じないように閉じているので，下端境界条件にはゼロフ
ラックス境界条件を与えた。

∂h⎛ ∂h ⎞− K ( h ) ⎜ + 1⎟
− K vh=0∂z z =− L⎝ ∂z ⎠ z =− L

(3.12)

ここで，L は試料長（L）である。水蒸気移動を考慮しない計算においては，左辺第 2 項は0 である。
試料表面からの蒸発を表すために，上端境界条件には次の地表面蒸発モデルを用いた（Rassam et al., 2004; 斎藤ら，2006）。

⎧∂h⎛ ∂h ⎞= E (t )⎪− K ( h ) ⎜ + 1⎟ − K vh∂z z =0⎝ ∂z ⎠ z =0⎨⎪h(0, t ) = hA⎩40

( h ( 0, t ) > h )A(3.13)

ここで，水蒸気移動を考慮しない計算においては，左辺第 2 項は 0 である。（13）式は，地表面の土中水圧力 h（0, t）が許容最小圧力 hA（L）より大きいときは，可能蒸発速度 E（L T−１）の一定フラックス境界条件を与えて恒率蒸発段階を表現し，地表面の乾燥が進行して h（0, t）が hA に達した段階で h（0, t） = hA の一定圧力境界条件を与えて減率蒸発段階を表現するモデルであり，可能蒸発速度 E は一定値で与えるのが一般的である。しかし，蒸発法によるパラメータの推定では，地表面の境界フラックスを，蒸発実験に忠実に再現することが重要である。そこで本解析では実測の蒸発速度 E を時間変動条件として与えた。許容最小圧力 hA には，安定した計算が行なえる範囲のできる限り小さい値として，hA = −106 cm を与えた。そのため，計算期間の大半で h（0, t） > hA であり，変動フラックス E（t）によって現実の恒率蒸発，減率蒸発段階を表現した。

(2) パラメータの推定
逆解析において最小化する目的関数Φは，土中水圧力の時間変化，実験終了時の平均体積含水率θ ave（砂丘砂は 63 時間後で 0.069 cm3 cm−3，黒ボク土は 208 時間後で0.46 cm3 cm−3），積算蒸発量の時間変化 Ecum(ti)のデータ，水分保持曲線の実測値（図3-2）を用いて定義した。
m

nj

j =1

i =1

*
Φ = ∑ v p , j ∑ wp ⎡⎣ p*j (ti ) − p j (ti ) ⎤⎦ + wθ ⎡⎣θ ave− θ ave ⎤⎦2nE+ vE ∑ wE ⎡⎣ E
i =1

*
cum

nθ

2

( ti ) − Ecum ( ti ) ⎤⎦ + vθ ( h) ∑ wθ ( h) ⎡⎣θ ( h ) − θ ( h )⎤⎦
2

*

2

i =1

(3.14)
ここで，上付*の有無はそれぞれ実測値と計算値を示し，pj (ti)は観測地点 j の時刻 ti における土中水圧力，m は観測地点の数，nj は観測地点 j におけるデータの数，Ecum(ti)は時刻 ti における積算蒸発量である。wp，wθ，wE，wθ (h)はそれぞれ，圧力，平均体積含水率，積算蒸発量，水分保持曲線のデータに与える重み成分である。異なる 2 地点の圧力データを組み合わせてパラメータを推定することで，パラメータの信頼度は向上する（坂井・取出，2007）。そこで本解析では，最上部と最下部の 2 点の測定点（鳥取砂丘砂では深さ 0.5 cm と 5 cm，黒ボク土では深さ 5.5 cm と 15.5 cm）を用いた。またθ ave は，推定された水分保持曲線のθ 軸上の位置を決める情報であるため，必ず目的関数に入れる必要がある（Šimůnek et al.,1998）。また，前節で述べたように，計算の過程において地41

表面圧力 h（0, t）が許容最小圧力 hA に達し，蒸発速度の計算値が実測値に基づく蒸発速度 E（t）よりも小さくなり，マスバランスを満たせなくなることがある。この問題を最小限にするため，積算蒸発量の時間変化 Ecum(ti)のデータを目的関数に加えて，可能な限り蒸発速度を再現するようにした。本解析では wp を 1 とし，水分保持曲線の実測値と蒸発実験における水収支を特に信頼できるデータとして重視するために，wθ (h)，wθ，wE をそれ
ぞれ 10 とした。また，v は絶対値やデータ数が大きく違う，異なる種類のデータを正規化する重み係数である。たとえば，砂丘砂の深さ 0.5 cm の土中水圧力は 0 > h > −700 cm
で変化するのに対し，水分保持曲線の実測値は 0 < θ < 0.35 であるため，変動の大きさが異なる。v は各データの変動の大きさを表す標本分散σ 2，およびそのデータ数 n を用いて次式で定義される。

v=

1
nσ 2

(3.15)

たとえば，（3.14）式右辺第 3 項に含まれる vE は，積算蒸発量の実測値の標本分散σE2 とそのデータの数 nE を用いて vE = 1/(nEσE2)で定義される。非線形最小 2 乗法に基づくLevenberg-Marquardt 法（Marquardt, 1963）により，目的関数Φを最小とするパラメータの最適化を行なった。
砂丘砂，黒ボク土に対し，それぞれ Fayer モデル（θ a，θ s，α，n，Ks，l），Durner モデル（θ r，θ s，α，n，w2，α2，n2，Ks，l）のパラメータを逆解析により推定した。このとき，飽和体
積含水率θ s は，封入空気量に依存する飽和度の違いにより変化するため（坂井・取出，

2007），水分保持曲線の実測値に基づく値（後述，3.4.1 節）で固定せず，推定パラメータとした。このため，（3.14）式の目的関数Φにおける水分保持曲線の実測値も，飽和近傍のデータは除外し，空気侵入圧付近以下のデータのみ（砂丘砂，黒ボク土ともに h ≤ −30cm）を用いた。同様に，飽和透水係数 Ks も飽和度に依存して変化するため，変水位透水試験による実測値（後述，3.4.1 節，表 3-1）で固定せず，推定パラメータとした。間隙結合係数 l は，水分保持曲線の形状とは独立して不飽和透水係数の圧力依存性を与える
重要なパラメータであるため，推定パラメータとした。また，VG モデルについても同様にパラメータ（θ r，θ s，α，n，Ks，l）の推定を行なった。
逆解析によるパラメータの推定では，初期値の組み合わせ方次第で，目的関数の最小値ではなく極小値に収束することがある。このため，適切な初期値を与えることが不可欠である（Rassam et al., 2004）。本解析では，推定パラメータの初期値には，水分保持曲42表 3-1 水分保持曲線の実測値に対する水分移動特性パラメータの適合値，および飽和
透水係数 Ks の実測値と Mualem（1976）による間隙結合係数 l の提唱値。
Table 3-1 Soil water hydraulic parameters of the VG, Fayer, and Durner models fit to theobserved soil water retention curves, the observed saturated hydraulic conductivity Ks, andthe pore-connectivity coefficient l presented by Mualem (1976).

Tottori dune sand
VG modelFayer model
Durner model
Andisol
VG modelFayer model
Durner modelθsα (cm−1)w 2 α 2 (cm− n 2
−
−
−
−
−
−
−4
1.5
0.02 9.96×10

θr

θa

0.030
−
0.002

−
0.341
0.085 0.336
0.337
−

0.035
0.036
0.036

n
5.20
7.06
6.63

0
−
0

0.769
−
0.754 0.724
−
0.744

0.192
0.040
0.041

1.14
−
−
−
14.92 −−−
2.44 0.59 1.01×10−4 1.47

K s (cmd−1)
550
550
550

l
0.5
0.5
0.5

200
200
200

0.5
0.5
0.5

線の実測値に各水分保持曲線関数をあてはめて得たパラメータ値を用いた（後述，3.4.1
節，表 3-1）。また，飽和透水係数 Ks の初期値には変水位透水試験で得た値（砂丘砂はKs = 550 cm d−1，黒ボク土は Ks = 200 cm d−1），間隙結合係数 l の初期値には Mualem（1976）の推奨値（l = 0.5）を用いた。
試料長，観測地点の深さを蒸発実験と等しく与え，（3.11）式の初期圧力には，各深さの測定値を線形補間して与えた。（3.13）式の上端境界条件に対しては，実測した積算
蒸発量の変化に基づき蒸発速度 E (t)を与えた。なお，水分移動の数値計算と逆解析は，不飽和水分移動汎用プログラム HYDRUS-1D（Šimůnek et al., 2005）を用いて行なった。

3.4

結果と考察

3.4.1 水分保持曲線関数

図 3-2（a）に鳥取砂丘砂の，図 3-2（b）に熊本の黒ボク土の水分保持曲線の実測値に対する Fayer モデル（実線），Durner モデル（点線），VG モデル（一点鎖線）の適合曲線を示す。実測値に対するモデルの適合は非線形最小 2 乗法で行なった。黒ボク土では，
低水分領域の水分保持曲線の実測値数が少ないことから，非線形最小 2 乗法が収束し43

なかったり，階段状の形状が得られなかったりした。そのため，Durner モデルを適合する
場合は，パラメータにθ r > 0，α 2 < 10−3，w2 > 0 といった制限範囲を設けて，最適化を行なった。適合によって得られたパラメータを表 3-1 に示す。砂丘砂に Fayer モデルを適合した場合，低圧力領域でθ が h の対数軸に対して直線的に減少する形状であるため，全圧力領域で実測値とよく適合した。また Durner モデルは，低圧力領域においてやや階段状を示す実測値の傾向もほぼ完全に再現した。黒ボク土に Fayer モデルを適合した
場合，高圧力領域（h ≥ −200 cm）では実測値とよく一致するが，低圧力領域（h < −200cm ）では水分量が直線的に減少するため，階段状の実測値を表現できなかった。
Durner モデルを用いた場合は，全圧力領域において階段状の水分保持曲線に非常によく適合した。一方，砂丘砂に VG モデルを適合した場合，高圧力領域の水分量の急激な低下をよく表現するが，一定値である残留体積含水率θ r（= 0.03）に漸近する形状のため，−80 ≥ h ≥ −1000 cm では実測値を過小評価し，h ≤ −10000 cm では過大評価した。また，黒ボク土に対しては，階段状を示す実測値を平均化するような緩やかな曲線となり，適合度は低かった。この表 3-1 のパラメータ値を，蒸発法の逆解析における初期値として用いた。

3.4.2 鳥取砂丘砂図 3-3 に鳥取砂丘砂の蒸発過程の圧力変化の実測値を示す。地表面に近い深さ 0.5

cm の実測値は，蒸発による乾燥に伴いおよそ 40 時間後から大きく低下し，h = −700 cm
程度にまで達した。一方で，深さ 5 cm の圧力低下は緩やかであり，h = −60 cm 程度しか低下しなかった。このように，表層付近で乾燥が進行するのが砂質土の蒸発過程の特徴である。

(1) 砂丘砂に対する水分移動特性関数
ここでは，まず水蒸気移動を含まない水分移動式を用いて逆解析を行ない，それぞれのモデルを比較検討した。図 3-3 に Fayer モデル，Durner モデル，VG モデルをそれぞれ用いて逆解析したときの蒸発過程圧力変化の計算結果を示す。また，図 3-4 には得られた水分保持曲線と不飽和透水係数を示し，それぞれのパラメータの推定値は表 3-2
に示した。
44

Pressure head, h (cm)

0

-200

0

observed
0.5cm
5cm

-400
-50

calculated
Fayer model
Durner model
VG model

-600
-100
0

20

40

60

80

-800
0

20

40

60

80

Time, t (hour)
図 3-3 蒸発過程における鳥取砂丘砂の土中水圧力変化：プロットは実測値，実線，点線，
一点鎖線はそれぞれ Fayer モデル，Durner モデル，VG モデルによる計算値。
Fig. 3-3 Pressure heads change during the evaporation process for Tottori dune sand: plots are
the observed data, and solid, dashed, and chain lines indicate predicted pressure headchanges using the Fayer, Durner, and VG model, respectively.

表 3-2 蒸発法による鳥取砂丘砂，熊本黒ボク土の水分移動特性パラメータの推定結果。
Table 3-2 Soil water hydraulic parameters for Tottori dune sand and Kumamoto Andisol
obtained inversely with the evaporation method.
Tottori dune sandθr
VG model0.0290
Fayer model−
Durner model0
Fayer (with vapor)−
Fayer (with K data)−
Andisol
VG model0.0004Durner model
0.0108θaθs− 0.330
0.089 0.317
0.333−
0.089 0.319
0.089 0.313
−
−

0.686
0.678

α (cm−1)
0.034
0.035
0.035
0.035
0.035

n
4.72
7.44
6.42
7.34
7.60

w 2 α 2 (cm− n 2 K s (cmd−1)
l
−
−
−
74.88
2.59×10−3
98.88
2.83×10−4−−−−4
0.17 6.73×10 1.42
65.24
3.30×10−3−−−
100.460.15
781.94
1.18
−
−
−

0.086
0.030

1.14
2.29

−
−
−
0.63 1.11×10−4 1.42

843.6
49.5

−0.092
0.87

Fayer モデルの場合，圧力変化の計算値は，表層付近のみで大きく低下する実測値の傾向をよく再現した（図 3-3）。また，水分保持曲線の推定値も実測値をよく再現した（図 3-4（a））。不飽和透水係数が h の対数軸に対して直線的に単調減少する VG モデル
に対し，Fayer モデルは，h < −80 cm で折れ曲がり，減少が緩やかになる形状を示す（図3-4（b））。Mualem モデルでは，有効飽和度 Se の関数として不飽和透水係数が導かれる（（3.1）式）。残留体積含水率θ r に漸近する VG モデルに対して，hm = −106 cm でθ = 0 と45

なる Fayer モデルでは，低水分領域の有効飽和度 Se が大きい。このため，不飽和透水係数の低下の割合が小さく，VG モデルに比べて大きな値となり，不飽和透水係数を過小評価する VG モデルの問題点が改善されている。折れ曲がった形状の不飽和透水係数は，物理的には，高水分領域では間隙流，低水分領域では土粒子表面に沿った膜流と水分量に応じて流れの形態が異なると考えられるため，不飽和透水係数の土中水圧力
依存性が圧力範囲によって異なることを示していると言える（Tuller and Or, 2001）。Fayer
モデルでは，全圧力範囲に対して水分保持曲線−不飽和透水係数連結モデルとして毛管保水の間隙流に基づく Mualem モデルを用いており，不飽和透水係数の推定では流れの形態の違いを考慮していない。しかし，この圧力変化の計算結果は（図 3-3），Fayer
モデルにより低水分領域の水分保持特性を適切に表現し，さらに間隙結合係数 l を推定
すると，蒸発過程の圧力変化をよく再現する低水分領域の不飽和透水係数が推定できることを示している。

0.4

Water content, θ (cm3cm-3)

observed
Fayer model
Durner model
VG model

0.3

100

0.2
10-4
0.1

0
100

101

102

103

104

105

106

101

102

103

104

105

10-8
106

Pressure head, h (-cm)

図 3-4 蒸発法で推定した鳥取砂丘砂の(a)水分保持曲線，(b)不飽和透水係数と土中水圧力の関係：プロットは実測値，実線，点線，一点鎖線はそれぞれ Fayer モデル，Durner
モデル，van Genuchten モデルによる推定値。
Fig. 3-4 (a) Water retention curve and (b) unsaturated hydraulic conductivity as a function of
the pressure head for Tottori dune sand estimated with the evaporation method: open circles
are the observed retention data, and solid, dashed, and chain lines are those predicted with
the Fayer, Durner, and VG model, respectively.
46

Hydraulic conductivity, K (cmd-1)

(b)

(a)

104

また，Durner モデルを適用した場合も，土中水圧力の時間変化（図 3-3），水分保持曲線の実測値（図 3-4（a））をよく再現した。Durner モデルの不飽和透水係数は，その関数型から 2 段の階段状を示すが，Fayer モデルとほぼ等しい低水分領域で低下が緩やかになる形状を示した（図 3-4（b））。Durner モデルは，2 つの VG モデルを結合して階段状の水分保持曲線を表現するモデルである。パラメータ数が多く自由度の高いため，Fayer モデルをほぼ包括することができ，図 3-4 の推定結果がよく一致したと考えられる。
しかし，推定パラメータ数が多い場合，逆解析のよる解の一意性が低くなる可能性がある。Fayer モデルは，VG モデルの低水分領域の問題を解消する上で，パラメータ数が少ない点で有利である。なお，VG モデルを用いて最適化を行なった場合，深さ 0.5 cm の圧力の計算値は 60 時間後においても−70 cm 程度までしか低下せず，実測値を再現することができなかった（図 3-3）。これは，VG モデルが砂丘砂の低水分領域の不飽和透水係数を正しく表現できないことが原因である。

(2) 水蒸気移動を考慮した推定
砂丘砂の蒸発過程では，地表面付近の乾燥が進行すると水蒸気移動が卓越する。このため，水蒸気移動を考慮すると，より精度の高い不飽和透水係数の推定が可能となると考えられる。図 3-5 に，（3.10）式により水蒸気移動を考慮して推定した Fayer モデルの不飽和透水係数（灰色実線）を，水蒸気移動を考慮しない場合の結果（点線）と併せて示す。図中では重なっているため判別できないが，両結果とも h ≥ −20 cm の高水分領域
ではおよそ K = 100 cm d−1 の一定値，そして−20 > h > −70 cm において 0.01 cm d−1 まで低下と一致した。パラメータの推定結果は表 3-2 に示す。水分保持曲線を表すパラメータθ a，θ s，α，n の推定値は，水蒸気移動の有無に関わらず，ほぼ同じ値となった。また，
水蒸気を考慮した圧力変化の計算値も，図 3-3 の液状水のみの計算結果とほぼ一致した。一方，間隙結合係数 l の推定値は，水蒸気移動を考慮しない場合は l = 2.83×10−4であるのに対して，水蒸気を考慮した場合はやや大きい l = 0.15 であった。この l の推定値の違いにより，h ≤ −60 cm の範囲で水蒸気移動を考慮した場合の方が不飽和透水係
数はやや小さな値となった。しかしその違いは小さく，蒸発法によりパラメータを推定する上で，水蒸気移動の影響は小さいと言える。
図 3-6 に 63 時間後の地表面付近の水分分布を示す。水蒸気移動を考慮した計算においては，地表面から深さ 0.2 cm 付近に乾燥層が形成されている。土中の水分フラック47スは，液状水フラックスと水蒸気フラックスの和として表されるが，水分量が非常に低い
乾燥層では，水蒸気移動が支配的に生じる。一方，水蒸気移動を考慮しない計算においては，表面付近で大きく水分量が減少する分布となっている。これは，（3.13）式の条
件では，液状水のみの移動で地表面蒸発速度を満たすために，仮想的な大きな圧力勾配を地表面付近に与えるためである（斎藤ら，2006）。このように表層付近の水分と圧力
分布は，水蒸気移動を考慮するかどうかにより異なるが，乾燥層の下部（深さ 0.5 cm 以下）の圧力変化から推定する水分移動特性が，水蒸気の有無に関わらずほぼ一致する

Hydraulic conductivity, K (cmd-1)

104
observed
liquid without K
liquid and vapor
without K
liquid with K

100

10-4

10-8 1
10

101

102

103

104

105

106

Pressure head, h (-cm)
図 3-5 水蒸気移動，および高水分領域における不飽和透水係数の実測値が，鳥取砂丘
砂の不飽和透水係数の推定に与える影響：プロットは実測値，点線，灰色実線，黒実線
はそれぞれ，水蒸気移動，実測値を考慮しない場合，水蒸気移動を考慮した場合，実測値を考慮した場合の Fayer モデルの推定値。
Fig. 3-5 Impacts of accounting for vapor flow and the observed unsaturated hydraulicconductivity data for higher water content on the estimation of hydraulic conductivity of
Tottori dune sand: open circles are the observed hydraulic conductivity data, and dashed,gray, and black solid lines describe predicted hydraulic conductivity using the Fayer model
without neither vapor flow nor the observed data, with vapor flow, and with the observed
data, respectively.

48

0
63h

Depth (cm)

0.2

0.4

0.6
without vapor
with vapor

0.8

1
0

0.02

0.04

Water content, θ (cm3cm-3)0.06

図 3-6 地表面付近における水分分布の計算値：点線，実線はそれぞれ水蒸気移動を考慮しない場合と考慮した場合の計算値。
Fig. 3-6 Calculated water content profiles near the soil surface with (solid line) and without(dashed line) considering vapor flow.

点は興味深い。この結果は，水分移動特性の推定に対しては，水蒸気移動の考慮よりも，
地表面水分境界フラックスを正しく与えることの重要性を示していると考えられる。また，
砂丘砂の蒸発過程では，土中水分移動に対する水蒸気移動の影響は他の土性に比べ大きい。そのため，他の多くの土性においても，蒸発法による不飽和透水係数の推定に対して，（3.13）式の境界条件を用いて蒸発速度を与えれば，水蒸気移動を考慮することによる効果は小さいと考えられる。

(3) 高水分領域の不飽和透水係数の推定
砂丘砂の蒸発実験において，深さ 5 cm の土中水圧力の実測値は，実験開始後 3 時
間で h = − 22 cm まで低下する。水分飽和状態から始まる蒸発実験では，飽和近傍の高土中水圧力領域の不飽和透水係数が影響するのは，実験開始後数時間の間だけであり，そのため，蒸発法を砂質土に適用する際，飽和透水係数 Ks は他のパラメータに比べて推定が難しいパラメータになる。Ks が他のパラメータと相関が高く推定が難しいことは，
49

蒸発過程の圧力変化に対する 2 次元応答曲面を用いたパラメータの感度解析からも示されている（Šimůnek et al.,1998; 坂井・取出，2007）。本解析において推定した鳥取砂
丘砂の Ks も，実測値（Ks = 550 cm d−1）に比べて小さな値である（たとえば，Fayer モデルの推定値 Ks = 98.88 cm d−1）。そこで，定常法（中野ら，1995）によって測定した不飽和透水係数の実測値（h = − 27.7 cm における K = 152.16 cm d−1，h = − 35.5 cm における K =8.23 cm d−1）を目的関数（（3.14）式）に代入し，逆解析により Fayer モデルのパラメータの推定を行なった。不飽和透水係数の実測値にかける重み成分は w = 10 とし，計算は水蒸気移動を考慮せずに行なった。図 3-5 に，不飽和透水係数の実測値を与えて逆解析した推定結果（実線）を実測値とともに示し，パラメータの推定結果を表 3-2 に示す。不飽
和透水係数の実測データを加えた場合，飽和透水係数の推定結果は Ks = 843.6 cm d−1となり，実測値に近づいた。このとき，不飽和透水係数の実測値を加えても，蒸発過程の圧力変化の計算値に違いは見られず，また，h < −40 cm では，両条件の不飽和透水係
数はほぼ一致している。このことは，砂質土の飽和付近の透水性が蒸発過程に及ぼす影響が小さいことを示唆している。これは，飽和透水係数の大きい砂質土における蒸発法で，飽和領域の不飽和透水係数を精度良く推定するためには，飽和近傍の不飽和透水係数の情報を与えるとよいことを意味している。また，飽和透水係数の実測値に十分な信頼が持てる場合，その実測値を（3.14）式の目的関数に重みをつけて代入することによっても，飽和近傍の不飽和透水係数推定の精度を向上させることができると考えられる。

3.4.3 熊本黒ボク土

図 3-7 に黒ボク土の蒸発過程の土中水圧力変化の実測値を示す。深さ 5.5 cm，15.5

cm の土中水圧力変化の測定値は，蒸発による試料の乾燥により，実験開始から約 170
時間後までは同時に緩やかに低下した。170 時間以降は，表面に近い深さ 5.5 cm の測定値が−400 cm 程度まで急激に低下し，一方で深さ 15.5 cm の測定値は−150 cm 程度まで緩やかに低下した。
図 3-7 に Durner モデルと VG モデルを水蒸気移動を含まない計算によって，最適化
した土中水圧力変化の計算結果を示す。また，図 3-8 に水分保持曲線と不飽和透水係数の推定結果，表 3-2 にパラメータの推定値を示す。なお，Fayer モデルを用いた逆解50

Pressure head, h (cm)0-100

observed
5.5cm
15.5cm-200

calculated
Durner model
VG model

-300
-400
0

50

100

150

200

250

Time, t (hour)
図 3-7 蒸発過程における熊本黒ボク土の土中水圧力変化：プロットは実測値，実線，一点鎖線はそれぞれ Durner モデル，VG モデルにより最適化した計算値。
Fig. 3-7 Pressure heads change during the evaporation process for Kumamoto Andisol: plots
are the observed data, and solid and chain lines indicate predicted pressure head changes
using the Durner and VG models, respectively.

103

0.8

Water content, θ (cm3cm-3)

100

0.6
10-3

0.4
10-6

0.2

0 0
10

observed
Durner model
VG model
101

102

103

10-9

104

105

106

101

102

103

104

105

10-12
106

Pressure head, h (-cm)

図 3-8 蒸発法で推定した黒ボク土の(a)水分保持曲線，(b)不飽和透水係数と土中水圧力の関係：プロットは実測値，実線，一点鎖線はそれぞれ Durner モデル，VG モデルによる推定値。
Fig. 3-8 (a) Water retention curve and (b) unsaturated hydraulic conductivity as a function of
the pressure head for Kumamoto Andisol estimated with the evaporation method: open
circles are the observed retention data, and solid and chain lines indicate those predictedwith the Durner and VG model, respectively.

51

Hydraulic conductivity, K (cmd-1)

(b)

(a)

析では，計算が不安定となり，パラメータが推定不能となった。これは，低土中水圧力領域の水分保持曲線を直線で表す Fayer モデルが，2 段の階段状を示す黒ボク土の水分保持曲線を適切に表現できないことが原因であると考えられる。

Durner モデルの場合，圧力変化の計算値は実測値とよく一致し（図 3-7），また，階段状の水分保持曲線の実測値をよく再現した（図 3-8（a））。Durner モデルの不飽和透水係
数は，水分保持曲線に対応して階段状を示し，h = −200 cm 付近で折れ曲がる形状を示した。これは，団粒間間隙と団粒内間隙の透水性の圧力依存性が異なることを反映している。すなわち，h > −100 cm の高水分領域の急激な不飽和透水係数の低下は，団粒間の粗間隙における水分の減少に伴う透水性の変化を示しており，h ≤ −100 cm の低水分領域での不飽和透水係数の緩やかな低下は，団粒内間隙の透水形態が支配的になったことを示していると解釈できる。黒ボク土の階段状の水分保持曲線によく適合し，団粒内外の透水性を反映した現実的な不飽和透水係数を表現できる Durner モデルは，
黒ボク土中の水分移動を計算する上で適した水分移動特性関数であると言える。
一方，VG モデルを用いて最適化を行なった場合，土中水圧力変化の実測値をある程度再現するものの，深さ 5.5 cm の 170 時間以降の急激な低下は再現できない（図3-7）。水分保持曲線は階段状の実測値をなだらかにした形状を示し，不飽和透水係数もそれに対応した単調に減少する形状を示した（図 3-8）。このように水分移動特性関数の適合が悪いため，VG モデルは黒ボク土中の圧力変化を適切に再現できなかったと考えられる。

3.4.4 不飽和透水係数の信頼性

Mualem モデルは，毛管保水と毛管内のポワズイユ流れを仮定して，土の不飽和透水性を表現したモデルである。砂丘砂においては，低水分領域の土粒子表面に沿った膜
流は，高水分領域の間隙を満たした流れとは透水形態が異なる。そのため厳密には，

Mualem モデルは，全水分範囲の流れの抵抗則を統一しては表現できないと考えられる。
また，団粒土においても，団粒内外で厳密には流れの抵抗則が異なる可能性がある。しかし，本研究において，2 種類の保水形態を反映した水分保持曲線‐不飽和透水係数
連結モデルである Fayer モデルと Durner モデルは，砂丘砂や黒ボク土の蒸発過程を良く表現することができ，それぞれの土の不飽和透水係数を適切に表現していると考えら52

れた。
一方，Mualem モデルにおける間隙結合係数 l は，屈曲度の影響を表すパラメータとして解釈されるが，Mualem 自身も 45 種類の試料の平均値として l = 0.5 を示しているに過ぎない（Mualem, 1976）。そのため l には，物理的な性格以上に，水分保持曲線とは独立して不飽和透水係数と土中水圧力の関係を与えるフィッティングパラメータとしての性格が強い。広い圧力範囲の水分保持曲線に対して適合性の高い Fayer モデルや

Durner モデルにおいて，さらに間隙結合係数 l を推定パラメータとして用いることにより，Mualem モデルによる水分保持曲線‐不飽和透水係数連結モデルの利便性を失うことな
く，飽和から乾燥領域までの広い圧力範囲の不飽和透水係数が表現できたと考えられる。
蒸発法による水分移動特性の推定の信頼範囲は，土中水圧力の測定範囲（本解析における砂丘砂では h > −700 cm）とされている（Šimůnek et al.,1998）。しかし，圧力の測定範囲を超えた領域の不飽和透水係数の信頼性は，用いた水分移動特性関数，言い換えると，Mualem モデルによる水分保持曲線‐不飽和透水係数連結モデルの信頼性に依存する。蒸発法に代表される逆解析手法での水分移動特性関数の推定において，推定結果の信頼性を厳密に論じることは難しい。しかし，本研究においては，砂丘砂の水分保持曲線によく適合する Fayer モデルは，蒸発過程の圧力変化の実測値をよく再現し，吸着保水が卓越する低水分領域の不飽和透水係数もよく表現していると考えられた。同様に，団粒構造を反映した水分保持曲線によく適合する Durner モデルは，団粒内間隙保水に対応する低水分領域の不飽和透水係数も適切に表現していると考えられた。こうした結果は，水分保持曲線を広範囲にわたり表現できるモデルを用いることの重要性を示している。
現在までの研究において，毛管内のポワズイユ流れ以外の仮定を用いた水分保持曲線‐不飽和透水係数連結モデルは，実用上利用できるものは存在しない。前述の本解析の結果に基づけば，Mualem モデルにおける間隙結合係数 l は，流れの抵抗則の不完全さを補正するパラメータとして位置づけることが可能であろう。このように Mualem モ
デルの実用上の利便性に加えて，低水分領域における流れの抵抗則に対する妥当性が裏付けられると，蒸発実験の測定範囲よりも低い圧力領域においても，水分保持曲線を正確に表現できれば，信頼度の高い不飽和透水係数の推定が可能であると考えられる。そのため，VG モデルの低水分領域を修正した Fayer モデルと Durner モデルは，実53

用上，非常に有用なモデルであると考えられる。

3.4.5 van Genuchten モデルの適用範囲

ここまでで VG モデルは，残留体積含水率θ r 付近の低水分領域の不飽和透水係数を表現できないことを示してきた。しかし，VG モデルは比較的な簡単な関数で，また，数値計算に相性のよいなめらかな関数であるため，非常に有用性が高い。そこで，水分保持曲線のみならず，不飽和透水係数を正しく表現できる，VG モデルの適用基準があれば便利であると考えられる。

van Genuchten（1978）は，S 字型を示す VG モデルの水分保持曲線の形状を特徴付けるものとして，θ s とθ r の中点，h−θ 関数（（3.2）式）の変曲点，log h−θ 関数の変曲点の3 点を挙げている。これらの 3 点における接線と直線θ = θ r の交点の土中水圧力 hB を下限値とし，VG モデルの適用範囲（hB ≤ h ≤ 0 cm）とする。ここでは，3 点について検討し，最も単純な関数で表される h−θ 関数の変曲点を VG モデルを特徴づける点として用いた。
図 3-9 に VG モデルの水分保持曲線を示す。点 A は変曲点を示し，その土中水圧力 hAは次式で表される（付録 3 参照）。

hA = − m1− m α(3.16)

変曲点における接線と直線θ = θ r の交点（図 3-9，点 B）の土中水圧力 hB は VG モデルのパラメータα と m を用いて，

hB = − m−1− m α(3.17)となり，hB を下限値とし，VG モデルの適用可能範囲とした。

−m−1− m α ≤ h ≤ 0(3.18)（3.18）式による VG モデルの適用可能範囲は，本研究で用いた鳥取砂丘砂（α = 0.036，

m = 0.857）では− 37.02 ≤ h ≤ 0 cm，熊本の黒ボク土（α = 0.0418，m = 0.561）では− 58.93≤ h ≤ 0 cm となる。蒸発法の適用範囲（h > − 1000 cm）に比べて非常に狭く，VG モデル
では蒸発法における全土中水圧力範囲の不飽和透水係数を表すことができないことを示している。ここで示した VG モデルの適用可能範囲は便宜的なものであるが，VG モデルが正しく不飽和透水係数を表すことができる範囲の目安として有効であり，この適用可能範囲を把握した上で，VG モデルを数値計算に用いることで，土中水分移動予測の信頼性が向上すると考えられる。
54

Water content, θ (cm3cm-3)

0.4

applicable range

0.3

A

0.2

0.1
B
0
0

20

40

60

Pressure head, h (-cm)

図 3-9 VG モデルの適用範囲。点 A は変曲点，点 B は変曲点における接線（点線）と直線θ = θr（一点鎖線）の交点。
Fig. 3-9 The applicable range for the VG model.

3.5

まとめ

本研究では，VG モデルの残留体積含水率θ r を修正した Fayer モデルと，2 つの VG
モデルを結合することで 2 段の階段状の水分保持曲線を表現した Durner モデルをMalem モデルに代入して導いた不飽和透水係数関数について，砂質土および団粒土への適用性を検討した。そのために，鳥取砂丘砂と熊本黒ボク土に対する水分蒸発実
験，蒸発法による Fayer モデル，Durner モデルのパラメータの推定を行ない，水分保持曲線と蒸発過程の圧力変化について，実測値と各モデルの計算値との比較を行なった。
鳥取砂丘砂に対して Fayer モデルを適用した場合，計算値は水分保持曲線の実測値によく適合するとともに，蒸発過程の圧力変化をよく再現し，その不飽和透水係数関数は，h < −80 cm で不飽和透水係数の低下が緩やかになる形状を示した。hm = −106 cm
で体積含水率θ = 0 となる Fayer モデルでは，低圧力領域でも有効飽和度 Se が大きな値を持つことで，θ r 付近で不飽和透水係数を過小評価する VG モデルの問題点が修正さ55れた。折れ曲がった形状の不飽和透水係数は，高水分領域の間隙流と低水分領域の土粒子表面における膜流の 2 つの異なる透水形態を反映していると考えられた。Durner
モデルは，Fayer モデルを包括するが，パラメータ数が多いため解の一意性が低いこと
からも，砂丘砂に対しては Fayer モデルが有利である。
砂丘砂については，水蒸気移動を考慮したリチャーズ式を用いたパラメータの推定も行なったが，水蒸気移動の考慮による不飽和透水係数の推定結果の違いは小さかった。
これは，水蒸気フラックスを無視し，液状水フラックスのみで境界フラックスを表現した場
合でも，蒸発実験の蒸発速度を正しく与えることの重要性を示している。乾燥層が形成される砂丘砂の蒸発過程は，他の土性に比べて水蒸気移動の影響は大きい。そのため，他の土を用いた蒸発法による不飽和透水係数の推定に対しても，境界条件としての蒸発速度を正しく与えれば，水蒸気移動を考慮する効果は小さいと考えられた。
黒ボク土に対して Durner モデルを適用した場合，水分保持曲線の実測値によく適合するとともに，蒸発過程の圧力変化をよく再現した。Durner モデルの不飽和透水係数関数は階段状を示し，h = −200 cm 程度で折れ曲がる形状を示した。これは，団粒間間隙と団粒内間隙の異なる 2 つの透水形態を反映していると考えられた。水分保持曲線によく適合し，より現実的な不飽和透水係数を表現できる Durner モデルは，黒ボク土中の水分移動の計算に適した水分移動特性関数であることが示された。
広い圧力範囲の水分保持曲線に対して適合性の高い Fayer モデルや Durner モデル
において，さらに間隙結合係数 l を推定パラメータとして用いることにより，Mualem モデ
ルによる水分保持曲線‐不飽和透水係数連結モデルの利便性を失うことなく，飽和から乾燥領域までの広い圧力範囲の不飽和透水係数が表現できた。また，蒸発法による水分移動特性推定の信頼範囲は，測定された圧力範囲とされている。しかし，砂質土に対する Fayer モデルや黒ボク土に対する Durner モデルのように，低水分領域の実測値を正確に表現できていれば，蒸発法で得られる不飽和透水係数は，測定した圧力範囲を超えても信頼度が高いと考えられた。そのため，VG モデルの低水分領域を修正した

Fayer モデルと Durner モデルは，実用上，非常に有用なモデルであると考えられた。

56

付録 1
（3.6）式中のΓs (h)，およびΓmax の定義を以下に示す（Fayer and Simmons, 1995）。まず，（3.4）式は次式に書き直すことができる。
nSe = ( X + γ ) Sa + ⎡⎣1 − ( X + γ ) Sa ⎤⎦ ⎡1 + ( −α h ) ⎤
⎣
⎦

−m

(A1)

ここで，Sa = θ a / θ s，X = 1−ln(−αh) / ln(−hm)，γ = ln(α) / ln(−hm)である。Mualem モデル（（3.1）式）の分子部分の積分は次式で表される。
Γ s ( h ) = I1 ( w ) α (1 − γ S a − S a ) +I3 ( w) ⎤⎡⎢ I 2 ( h, w ) +⎥n ⎦ln ( − hm ) ⎣

α Sa

(A2)

ここで，w = [1+(−αh)n]−1 であり，I1(w)，I2(h, w)，I3(w)の定義は以下に示す通りである。

I1 ( w ) = (1 − wm ) − (1 − w )m

I 2 ( h, w ) =m

1 ⎛1 1 ⎞( m −1)( m −1)− (1 − w )⎜ − ⎟ + (1 − wm )α ⎝ h hm ⎠I 3 ( w ) = I 3a ( w ) + I 3b ( w0 ) w > w0
I 3 ( w ) = I 3b ( w )w ≤ w0

(A3)
(A4)

(A5)

ここで，wm は土中水圧力 hm に対する w の値であり，w0 は 10−10 である。（A5）式中の I3a(w)，I3b (w)は次式で表される。
⎛ 1 − w0 ⎞ 1 ⎤m ⎡m ⎡⎛ 1− w ⎞ 1 ⎤I 3a ( w ) = (1 − w0 ) ⎢ln ⎜⎟ − ⎥ − (1 − w ) ⎢ ln ⎜⎟ − ⎥ + G (1 − w ) − G (1 − w0 )⎣ ⎝ w ⎠ m⎦
⎣ ⎝ w0 ⎠ m ⎦(A6)
I 3b ( w ) = m ⎡⎣( m − wm ) + wm ln ( wm ) − w ln ( w ) ⎤⎦(A7)

また，（A6）式中の G(w)は次式で表される。
G ( w ) = Gs (W0 ) + H (1 − W0 ) − H (1 − w ) w > W0
G ( w ) = Gs ( w )w ≤ W0(A8)

ここで，W0=0.9 であり，H(w)，Gs(w)はそれぞれ，H ( w ) = ln ( w ) − mw + m ( m − 1)w2
w3
− m ( m − 1)( m − 2 )
4
18

100

wk
k =0 m + 1 + k

Gs ( w ) = wm +1 ∑である。また，Mualem モデル（（3.1）式）の分母部分の積分は次式で表される。
57

(A9)

(A10)

Γ max = Γ s ( h0 ) + Γ 0 c(A11)

ここで，h0 は（A2）式を評価することができる最大の土中水圧力を示し，およそ 0.025/αで表すことができる。またΓ0c (h)は，⎛⎞Sa( n −1)( n −1)
Γ 0 c = α ⎜⎜1 − γ Sa − Sa +− ( −α hc ) ⎤⎟⎟ ⎡⎣( −α h0 )⎦n ln ( − hm ) ⎠⎝α Sa ⎡
1 ⎞
1 ⎞⎤
n −1
( n −1) ⎛+( −α h0 )( ) ⎛⎜ ln ( −α h0 ) −⎟ − ( −α hc )⎜ ln ( −α hc ) −⎟⎢ln ( −hm ) ⎣n −1 ⎠
n − 1 ⎠ ⎦⎥
⎝
⎝
(A12)
である。hc は（A2）式が数値的に 1 を示す最小の土中水圧力を示し，本研究では−10−6

cm を用いた。

付録 2
圧力勾配による水蒸気透過係数 Kvh の定義を以下に示す。

K vh =

D

ρw

ρvs

Mg
Hr
RTabs

(B1)

ここで，D は土中の水蒸気拡散係数(L2T−1)，ρw は液状水密度（ML−3 = 1 − 7.37 × 10−6 ×(T−4)2 + 3.79 × 10−8 × (T−4)3 kg m−3），T は相対温度（K，本研究では 20 C で一定），ρvs
は飽和水蒸気密度（ML−3 = T−1 × 10−3 × exp(31.37− 6014.79T−1 − 7.92 × 10−3T) kg m−3），

M は水の分子量（Mmol−1, = 0.018015 kg mol−1），g は重力加速度（LT−2, = 9.81 m s−2），
R は気体定数（ML2T−2mol−1K−1, = 8.341 J mol−1 K−1），Tabs は絶対温度（K），Hr は相対湿度である。土中の水蒸気拡散係数 D は次式で導かれる。

D = τ aθ air Da(B2)

ここで，θ air は気相率（L3L−3），Da は大気中の水蒸気拡散係数（L2T−1 = 2.12 × 10−5
（Tabs/273.15）2 m2 s−1）である。τ a は気相の屈曲度で，次式で示される（Millington and

Quirk, 1961）。

τa =θ air 7/3
θ s2

(B3)

相対湿度は次式で導かれる（Philip and de Vries, 1957）。
⎡ hMg ⎤
H r = exp ⎢⎥⎣ RTabs ⎦
58

(B4)

付録 3

VG モデルの変曲点は，（3.2）式の 2 階微分が 0 となる点で定義される。
d 2 Se
=0
dh 2

（C1）

m = 1 − 1 / m の関係を用いて（3.2）式から n を消去すると，{}

1 1− m − m

Se = 1 + ( −α h )(C2)となる。（C2）式を（C1）式に代入すると，

−

mα 2

(1 − m )−α h )2 (2 m −1 1− m{

} i{m − ( −α h ) } = 0

1 1− m − m − 2

i 1 + ( −α h )1 1− m(C3)となる。h = 0，h → ∞以外で（C3）式を満たす変曲点の土中水圧力 hA は，

hA = − m1− m α(C4)となり，本文中の（3.16）式が得られる。変曲点における接線は，hA を用いて次式となる。

θ=(θ s − θ r ) m2−m
( h − hA ) + {(1 + m ) (θ s − θ r ) + θ r }1+ m(1 + m ) ( m − 1) hA(C5)(C5)式と直線θ =θ r の交点の土中水圧力 hB は，

hB = − hA m 2となり，本文中の（3.17）式が得られる。

59

(C6)

第4章

4.1

温度勾配下における凝縮・蒸発を伴う砂中の水分移動

はじめに乾燥地や半乾燥地において，水循環や溶質移動を予測する上で，液状水による水分移動のみではなく，水蒸気移動について把握することは重要である。1 日の気温変化による地表面温度変化は，土中の水蒸気移動に大きな影響を与え，時にそれは土中の水分移動において支配的になる（Milly, 1984）。さらに水蒸気移動は，蒸発潜熱として多量の熱量を運ぶため，熱移動に対しても大きな影響を持つ（Cahill and Parlange, 1998;

Saito et al., 2006）。
土中の水蒸気移動理論は，Fick の拡散則に基づき，古くは Penman（1940）により提唱された。その後，Philip and de Vries（1957）は，非等温条件下の液状水と水蒸気の同時移動について理論化し，全水分フラックスは，水分量勾配と温度勾配による液状水フラックスと水蒸気フラックスの 4 成分に分離できることを示した。Milly（1982）は，水分量勾配を土中水圧力勾配に置き換え，Philip and de Vries モデルを修正した。この修正により，
ヒステリシスや不均一な土層を考慮することが可能となった。Nassar and Horton （1989）
や Nassar et al.（1992b）は，水蒸気，液状水移動に対する浸透圧の効果を考慮すること
で，Philip and de Vries モデルをさらに展開した。Philip and de Vries（1957）はさらに，温度勾配による水蒸気フラックスに対して，水蒸気促進係数（enhancement factor）を導入した。土中の水蒸気は，土中間隙に存在する「液島」を，一方で凝縮し，一方で蒸発することで通過するため，水蒸気フラックスは促進される。また，気相の局所的な温度勾配は土全体の平均温度勾配に比べ著しく大きいことを，水蒸気フラックスが促進されるもう 1 つの理由として挙げている。水蒸気促進係数を直接測定することは難しい。 Cass et al.
（1984）は，熱伝導率の圧力依存性に基づいて，水蒸気促進係数を間接的に推定している。
一方で，乾燥した土中の液状水フラックスを予測する上で，低水分領域の不飽和水分移動特性を評価することが必要不可欠となる。不飽和水分移動特性を測定する方法は多く提案されているが（Klute and Dirksen, 1986），低水分領域の不飽和透水係数を測定することは非常に困難である。近年，マルチステップ流出法（Eching and Hopmans, 1993;

Inoue et al., 1998）や蒸発法（Kool et al., 1987; Šimůnek et al., 1998）といった，土中の水60分移動過程のデータをもとに不飽和水分移動特性を推定する逆解析手法が提案されている。不飽和透水係数の関数は，間隙径分布に基づき Burdine （ 1953 ）や Mualem（1976）により提案されている。Mualem のモデルは，間隙の連続性や屈曲度を考慮した間隙結合係数を持ち，Mualem は 45 種類もの土試料の平均値として 0.5 を提唱している（Mualem モデルを Brooks and Corey（1964）のモデルに適用したとき，間隙結合係数の値は 1 となる）。間隙結合係数は，土中水圧力に対する不飽和透水係数の傾きに影響を与え重要なパラメータであるにも関わらず，0.5 に固定して用いられることが多い。
そこで本研究の目的は，まず，水蒸気凝縮過程における水分分布の実測値に基づき，
低水分領域の砂の不飽和透水係数を推定することである。温度勾配下の乾燥した砂カ
ラム中の，水蒸気凝縮・蒸発を伴う，液状水，水蒸気，熱移動を測定した宮﨑（1976）の実験に対して，Philip and de Vries モデルを含む HYDRUS-1D（Šimůnek et al., 2005）を用いて数値計算を行なった。そして，間隙結合係数が数値計算結果に与える影響について調べた。2 つ目の目的は，全水分フラックスの 4 成分，土中水圧力勾配と温度勾配による液状水フラックスと水蒸気フラックス，を定量的に評価することである。また水収支に基づき，凝縮速度と蒸発速度を評価する。さらに，凝縮・蒸発が生じる過程における水蒸気促進係数の影響について調べる。

4.2

試料と方法

4.2.1 凝縮実験

宮﨑（1976）は，1 次元の室内カラム実験により，高温で多湿な大気から低温で乾燥した土への水蒸気移動について測定を行なった。試料には風乾した浜岡砂丘砂（体積含
水率 0.0045 cm3 cm−3）を用い，高さ 10 cm，内径 10 cm のアクリルカラムに乾燥密度 1.5 g

cm−3 で充填している。作製した砂カラムを，37 ºC で相対湿度 85 ∼ 90 %に設定した恒温
恒室チャンバーに設置し，カラム上端は湿った大気に開放としている。また，カラム下端は閉鎖系とし，20 ºC の水をポンプを用いて循環させることで，一定温度に保っている。
カラム側面を断熱することで熱の流れを 1 次元に制御し，熱電対を用いて試料中の温度分布を測定している。カラム総重量を測定し，その変化量から積算水蒸気侵入量を測定している。また，カラムは複数用意し，任意の時間にカラムを 1 cm または 2 cm 幅で切断
61

し，炉乾法により体積含水率分布の時間変化を測定している。
浜岡砂丘砂の排水過程の水分保持曲線は吸引法で測定した。また，吸水過程の圧力変化と水分量変化をそれぞれ，圧力センサーと TDR プローブで同時測定することで，吸水過程の水分保持曲線を測定した。また，変水位透水試験で測定した飽和透水係数は，3456 cm d−1 であった。

4.2.2 数値計算

(1) 不飽和水分移動式

土中水の表面張力は温度に依存して変化するため，温度変化が生じる条件下では，土中の液状水フラックス qL を表すダルシー則は，土中水圧力勾配のみではなく，温度勾配も考慮する必要がある（Philip and de Vries, 1957）。また，水蒸気密度は土中水圧力と温度の関数であるため，フィックの拡散則で表される水蒸気フラックス qv も土中水圧力勾
配と温度勾配の 2 成分で表される。そして，不飽和土中の全水分フラックス qTotal は液状
水フラックスと水蒸気フラックスの和で表される。

qTotal = qL + qv = qLh + qLT + qvh + qvT =∂T∂h∂T⎛ ∂h ⎞= − K Lh ( h ) ⎜ + 1⎟ − K LT
− K vh− K vT∂z∂z∂z⎝ ∂z ⎠(4.1)

ここで qLh は土中水圧力勾配による液状水フラックス（LT−1），qLT は温度勾配による液状
水フラックス（LT−1），qvh は土中水圧力勾配による水蒸気フラックス（LT−1），qvT は温度勾配による水蒸気フラックス（LT−1），KLh（LT−1）と KLT（L2K−1T−1）はそれぞれ土中水圧力勾
配，温度勾配による液状水移動における不飽和透水係数，Kvh（LT−1）と KvT（L2K−1T−1）はそれぞれ土中水圧力勾配，温度勾配による水蒸気移動における透過係数，h は土中水
圧力（L），T は温度（K），z は上向きを正とする空間座標(L)を示す。（4.1）式を水の保存則に代入することで，非等温条件下のリチャーズ式が導かれる（ Nassar and Horton,

1997; Noborio et al., 1996a）。
∂q∂θ= − Total =
∂t
∂z
∂ ⎡
∂T
∂h
∂T ⎤
⎛ ∂h ⎞= ⎢ K Lh ( h ) ⎜ + 1⎟ + K LT
+ K vh+ K vT
∂z ⎣
∂z
∂z
∂z ⎦⎥
⎝ ∂z ⎠
62(4.2)

ここで，θ は液状水量と水蒸気量を合せた全体積含水率（L3L−3），t は時間（T）を示す。
（4.2）式に示す水の保存則は，液状水と水蒸気に関する 2 つの保存則に分割することができる（Nassar and Horton, 1992）。

∂θ L∂q=− L −E
∂t
∂z

(4.3)

∂θ v
∂q
=− v +E
∂t
∂z

(4.4)

ここで，θL は液状水の体積分率（L3L−3），θv は液状水等量の水蒸気の体積分率（L3L−3），

E は蒸発速度（T−1）を表す（負の E は凝縮速度を表す）。

(2) 不飽和水分移動特性
水分保持曲線関数には，Fayer and Simmons（1995）により提案されたモデル（Fayer モ
デル）を用いた。このモデルは，低水分領域における土粒子表面における吸着保水を考慮し（Campbell and Shiozawa, 1992），数値計算に広く用いられる van Genuchten モデル（van Genuchten, 1980）を修正したものである。

θ L = χθ a + (θ s − χθ a ) ⎡1 + ( −α h ) ⎤
n

⎣

−m

⎦

(4.5)

ここで，θs は飽和体積含水率（L3L−3），α（L−1），n（−），m （= 1 − 1/n），θa （L3L−3）はそれ
ぞれ水分保持曲線の形状を示すパラメータである。（4.5）式の右辺第 1 項は吸着保水項
を示しており，χ は次式で表される。

χ ( h ) = 1 − ln ( −h ) ln ( − hm )(4.6)

ここで，hm は体積含水率が 0 となる炉乾時の土中水圧力を示し，一般的に−107 cm が用いられる。
土中水圧力勾配による液状水フラックスにおける不飽和透水係数 KLh は，（4.5）式をMualem（1976）の間隙径分布モデルに代入することで導かれる。
⎡ Se dSe
K Lh = K s Se ⎢ ∫
0
h
⎣
l

dSe ⎤
∫0 h ⎥⎦
1

2

(4.7)

（4.7）式の積分の評価は Fayer and Simmons（1995）に詳しい。Brooks and Corey モデル（Brooks and Corey, 1964）は，間隙径分布に基づく，広く用いられる不飽和水分移動特性関数の 1 つである。
63⎧⎪(θ s − θ r ) α h − n + θ r
θL = ⎨θs⎪⎩h < −1/ α(4.8)

h ≥ −1/ α

K Lh ( h ) = K s Se 2 / n +l + 2(4.9)

ここで，θr は残留体積含水率（L3L−3），α（L−1），n（−）はそれぞれ水分保持曲線の形状を示すパラメータ，Se（= （θ − θr）/（θs − θr））は有効飽和度である。Brooks and Corey（1964）
は，（4.9）式の間隙結合係数 l の推奨値を 1 としている。
温度勾配による液状水フラックスにおける不飽和透水係数 KLT は，次式で定義される（e.g., Noborio et al., 1996a）。
⎛
1 dγ ⎞K LT = K Lh ⎜ hGwT⎟
γ 0 dT ⎠⎝(4.10)

ここで，γ は土中水の表面張力（MT−2）で，その温度依存性はγ = 75.6 − 0.1425Τ − 2.38

× 10−4T 2 (g s−2)と表される。ここでγ 0は 25 °C における表面張力である（= 71.89 g s−2）。
GwT は表面張力の温度依存性を修正する促進係数であり，本研究では 7 に固定して用いた（Nimmo and Miller, 1986）。
土中水圧力勾配，および温度勾配による水蒸気フラックスの透過係数 Kvh，KvT はそれ
ぞれ次式で定義される。

K vh =Dρw

K vT =

ρvs

D

ρw

Mg
Hr
RTabs

(4.11)d ρ vsdT(4.12)

η Hr

ここで，D は土中の水蒸気拡散係数(L2T−1)である。ρw は液状水密度（ML−3）でρw = 1− 7.37 × 10−6(T−4)2 + 3.79 × 10−8(T−4)3 （kg m−3）と定義される。ρvs は飽和水蒸気密度で（ML−3），ρvs = T−1 × 10−3exp(31.37− 6014.79T−1 − 7.92 × 10−3T) （kg m−3）と定義される。

M は水の分子量（Mmol−1, = 0.018015 kg mol−1），g は重力加速度（LT−2, = 9.81 m s−2），
R は気体定数（ML2T−2mol−1K−1, = 8.341 J mol−1 K−1），Tabs は絶対温度（K），Hr は相対湿度である。また，ηは水蒸気促進係数であり，Cass et al.（1984），Campbell（1985）は次式で定義している。
⎧

4⎫θ
2.6 ⎞ θ L ⎤ ⎪⎪ ⎡⎛η = a + 3 L − ( a − 1) exp ⎨− ⎢⎜1 +
⎟ ⎥ ⎬
⎜
θs
f c ⎟⎠ θ s ⎥⎦ ⎪
⎪ ⎢⎣⎝
⎩

64

⎭

(4.13)

ここで，fc は粘土分率（−）で本研究では 0.02 に固定した。a は定数で，Cass et al.（1984）
は，9.5 としている。土中の水蒸気拡散係数 D は，大気中の水蒸気拡散係数 Da（L2T−1）に気相の屈曲度τa（Millington and Quirk, 1961）を掛け合わせることで導かれる。

D = τ aθ a Da =

θ a7 / 3
θ D
θ s2 a a(4.14)

ここで，θa は気相率（L3L−3）である。
相対湿度は次式で表される（Philip and de Vries, 1957）。
⎡ hMg ⎤
H r = exp ⎢⎥⎣ RTabs ⎦(4.15)(3) 熱移動式
熱の保存則の微分方程式は次式で表される（de Vries, 1958）。
∂C pT
∂t

+ L0

∂θ v ∂ ⎡
∂q T
∂q
∂T ⎤
∂q T
= ⎢ λ (θ L ) ⎥ − Cw L − Cv v − L0 v
∂t
∂z ⎣
∂z ⎦
∂z
∂z
∂z

(4.16)

ここで，L0（= Lwρw）は単位体積当たりの水の蒸発潜熱（J m−3）で，Lw（= 2.501 × 106 –2369.2 T ）は 単位重量当た り の 蒸発潜熱（ J kg−1 ） であ る 。 λ(θL) は 土の 熱伝導率
（MLT−3K−1）で次式で表される（Chung and Horton, 1987）。

λ (θ L ) = b1 + b2θ L + b3θ L 0.5(4.17)

ここで，b1，b2，b3 は定数で，砂質土では b1 = 0.228 W m−1 K−1，b2 = −2.406 W m−1 K−1，

and b3 = 4.909 W m−1 K−1 である。（4.16）式の左辺第 1 項と第 2 項はそれぞれ，土中の顕熱変化量と水蒸気による潜熱変化量を示す。（4.16）式の右辺は順に土中の熱伝導項，
液状水移動による顕熱輸送項，水蒸気拡散による顕熱輸送項，水蒸気拡散による潜熱輸送項をそれぞれ示す。土の体積熱容量 Cp（ML−1T−2K−1）は，固相の体積熱容量 Cs，
液相の体積熱容量 Cw，気相の熱容量 Cv にそれぞれの体積分率θ を掛けて，足し合わせることで表される。

C p = Csθ solid + Cwθ L + Cvθ v(4.18)(4) 溶質移動式
土中の溶質移動は移流分散式（the convection-dispersion equation: CDE）で表される。

65

∂θ L c ∂ ⎡∂c ⎤ ∂q c= ⎢(θ L Dwτ w + λL | qL |) ⎥ − L∂t∂z ⎣∂z ⎦ ∂z(4.19)

ここで，c は溶液濃度（ML−3），Dw は水中の分子拡散係数（L2T−1），λL は分散長（L）である。τw は液相の屈曲度で，次式で表される（Millington and Quirk, 1961）。

τw =θ L7 / 3θ s2(4.20)

本研究では，Dw = 1.27 cm2 d−1, and λl = 0.5 cm とした。

4.2.3 初期・境界条件

初期条件，境界条件は，水蒸気凝縮実験に合せて設定した。初期体積含水率は，砂カラム中で一定に設定した。

θ L ( z, 0 ) = 0.0045( −10 ≤ z ≤ 0 )(4.21)

水分移動の上端境界条件は，全水分浸入量で表した。

qTotal (0, t ) = ( qLh + qLT + qvh + qvT ) z =0(t > 0)(4.22)

全表面フラックス qTotal (0, t)は，砂試料中の水蒸気凝縮量を反映する，カラムの総重量の変化量とした。カラム総重量の変化（積算凝縮量）を 2 次関数で近似し，その時間微分から qTotal (0, t)を求めた。
カラムの下端は閉じてあるので，下端境界条件にはゼロフラックス境界条件を用いた（Nassar and Horton, 1992）。

qTotal (−10, t ) = ( qLh + qLT + qvh + qvT ) z =−10 = 0(t > 0)(4.23)

液状水フラックス（qLh と qLT の和）と水蒸気フラックス（qvh と qvT の和）が 0 となるように設定した。
初期温度は，T ( z , 0 ) = 23.5

( −10 ≤ z ≤ 0 )(4.24)とした。熱移動の計算には，ディリクレ型の境界条件を適用した。実験上は，試料境界の周りの温度を設定するため，試料境界の温度はその設定温度とはならない（ Farlow,1993）。このため，上端，下端の境界温度は，境界近辺で測定された試料の平均温度とした。
T ( 0, t ) = 36.3
66

(t > 0)(4.25)

T ( −10, t ) = 23.5

(t > 0)(4.26)( −10 ≤ z ≤ 0 )(4.27)

初期溶液濃度は，c ( z , 0 ) = 0.93とした。上端と下端の溶質境界条件には，溶液フラックス 0 の境界条件を与えた。

∂c ( 0, t )∂z=0∂c ( −10, t )∂z=0(t > 0)(4.28)(t > 0)(4.29)（4.2）式，（4.16）式の偏微分方程式で表される，液状水，熱，水蒸気同時移動は，修正版の HYDRUS-1D（Šimůnek et al., 2005）を用いて，数値計算により求めた。

4.2.4 不飽和水分移動特性パラメータ

図 4-1（a）は吸引法と加圧板法で測定した，浜岡砂丘砂の排水過程の水分保持曲線を示す。Fayer モデル（（4.5）式）を用いて，最小 2 乗法により適合し，θ a = 0.0354，θ s =0.318，α = 0.0316 cm−1，n = 7.59 が得られた。砂質土は排水過程と吸水過程の水分保
持曲線の間に，大きなヒステリシスを持つことが知られている。凝縮実験は吸水過程で生じるため，計算には吸水曲線を用いた。吸水過程の水分保持曲線は，Kool and Parker（1987）の過程に基づき，排水過程の水分保持曲線から得られたαを 2 倍して用いることで推定した（α = 0.0632 cm−1）。
図 4-1（b）は，Fayer モデルを Mualem モデルに代入することにより導かれた不飽和透水係数と土中水圧力の関係を示す。ここで，飽和透水係数 Ks には実測値（3456 cm d−1）を用いた。間隙結合係数 l が計算結果に与える影響を調べるために，Mualem の推奨値
（l = 0.5）の代わりに，様々な値を用いた（例，l = 5.4，図 4-1（b））。l の値を大きくすること
で，高土中水圧力における不飽和透水係数の傾き dK/dh が大きくなり，特に，低水分領域の不飽和透水係数を小さくする。また，図 4-1（b）に Brooks and Corey モデルにより導かれた不飽和透水係数を示す。Fayer モデルの不飽和透水係数が折れ曲がる形状を示すのに対し，Brooks and Corey モデルの不飽和透水係数は直線的に単調減少する形状を示す。本研究では，比較として，Brooks and Corey モデルを用いた計算も行なった。

67

104

(b)

observeddrainage
imbibition

(a)

0.3

10-4

F&S
l = 0.5

0.2
F&S
l = 5.40

10-12

0.1
B&C
l=5
0
100

101

102

103

104

105

100

101

102

103

104

105

Hydraulic conductivity, K (cmd-1)

Water content, θL (cm3cm-3)

0.4

10-20
106

Pressure head, h (-cm)

図 4-1 （a）浜岡砂丘砂の水分保持曲線の実測値と Fayer モデルによる適合曲線。（b）
Fayer モデルと Brooks and Corey モデルによる不飽和透水係数。
Fig. 4-1 Water retention curve (a) and unsaturated hydraulic conductivity functions (b)described using the Fayer model and the Brooks and Corey model.

Enhancement factor, η20

16

a = 15 (Case 6)12
a = 9.5 (Case 1, 2, 3)
8
a = 5 (Case 5)
4
η = 1 (Case 4)
0
0

0.1

0.2

Water content, θL (cm3cm-3)0.3

図 4-2 異なるパラメータ a を用いた場合の，水蒸気促進係数ηと体積含水率θL の関係。
Fig. 4-2 Enhancement factor for different values of parameters a and b.

68さらに，水蒸気促進係数ηが計算結果に与える影響を調べるために，（4.13）式におけるパラメータ a を変化させて用いた（図 4-2）。また，水蒸気促進係数を考慮しない場合（η= 1）の計算も行なった。

4.3

結果と考察

4.3.1 不飽和透水係数

まず，不飽和水分移動特性関数に Fayer モデルを用い，間隙結合係数 l には Mualem（1976）の推奨値である 0.5 を用いた結果を示す。体積含水率分布の実測値と計算値を
図 4-3 に示す。実測値と同様に計算値は，水蒸気の凝縮により，カラム下端から水分量が増加する過程を再現した。しかし，計算の体積含水率分布の形状は，実測値の形状と大きく異なった。計算の体積含水率は，カラム下端付近で実測値を過小評価し，一方で，カラム中央付近では過大評価した。この実測と計算の相違は，液状水の毛管上昇の過大評価，もしくは下方への水蒸気移動の過小評価が原因だと考えられる。間隙結合係数 l を大きくすることにより，不飽和透水係数 K を小さくすることができる（図 4-1（b））。そこでまず，l を大きな値に設定することで K を小さくし，液状水の上方移動を減少させた。
図 4-4（a）は，l = 5.4 を用いて計算した場合の体積含水率分布を示す。液状水の上方移動を減少させることで，計算値は実測値をよく再現することができた。
また，図 4-4（b）に，広く用いられている不飽和水分移動特性関数である Brooks andCorey モデルを用いて計算した結果を示す。間隙結合係数は Brooks and Corey（1954）の推奨値である l = 1 の代わりに，l = 5 を用いることで，不飽和透水係数を小さくした（図4-1（b））。計算結果の水分量増加前線の位置は，実測値のものをよく再現した。しかし，実測値の水分分布が緩やかな勾配の形状を示すのに対し，計算の結果は急な勾配の分布を示した。土への水の浸潤を考えた場合，湿潤域と乾燥域の遷移領域は，細粒土に比べ粗粒土のほうが著しく狭くなる（Jury and Horton, 2006）。これは，土中水圧力に対する不飽和透水係数の傾き dK/dh が細粒土に比べ，粗粒土の方が大きくなるためである。
この水蒸気凝縮過程における液状水移動を試料下端からの浸潤と仮定した場合，l = 5で表された Brooks and Corey モデルの不飽和透水係数は急な dK/dh の勾配を持つため，水分量増加前線の形状が急になったと考えられる。一方で，Fayer モデルの不飽和透水690observed
5d
10 d
16 d

Depth (cm)2

4calculated
initial
2d
10 d
16 d

6

8

10
0

0.05

0.1

Water content, θL (cm cm )
3

0.15

-3

図 4-3 Fayer モデルと間隙結合係数 l = 0.5 を用いた場合の水分分布の計算値と実測値の比較。
Fig. 4-3 A comparison of observed and simulated volumetric water content profiles calculatedusing the Fayer model with the pore-connectivity coefficient l = 0.5.

0

(b)

(a)
observed
5d
10 d
16 d

Depth (cm)2

4calculated
initial
5d
10 d
16 d

6

8

10
0

0.05

0.1

0

0.05

Water content, θL (cm3cm-3)

0.1

0.15

図 4-4 （a）Fayer モデルと間隙結合係数 l = 5.4 を用いた場合の水分分布の計算値と実測値の比較，（b）Brooks and Corey モデルと間隙結合係数 l = 5 を用いた場合の水分分布の計算値と実測値の比較。
Fig. 4-4 A comparison of observed and simulated volumetric water content profiles calculatedwith the Fayer model (a) and calculated with the Brooks and Corey model.

70係数は折れ曲がる形状となるため，低土中水圧力領域で，緩やかな dK/dh の勾配を示す。この不飽和透水係数の形状により，Fayer モデルは水蒸気凝縮過程の水分量変化をよく再現することができたと考えられ，また，Fayer モデルは低水分領域も含む，広い水分領域の不飽和透水係数をよく表していると言える。以下については，不飽和水分移動
特性関数に Fayer モデルを用い，間隙結合係数 l を 5.4 として計算した結果について示す。

4.3.2 温度分布

図 4-5 に，温度分布の実測値と計算値を示す。計算値は，上端，下端境界温度（（4.25）式，（4.26）式）により，実験開始から 0.1 日間の短時間で直線分布となった。6 日目と 16 日目の計算温度分布は，ほぼ直線分布を示したが，6 日目では深さ 9 cm で，16日目では深さ 7 cm でやや上に凸の形状を示した。この直線分布からのわずかなずれは，水分量増加前線における液状水の蒸発（後述）と，水分量が多いことにより下層の方が熱伝導率が大きいことが原因である。6 日目，16 日目ともに，計算値は実測値をよく再現した。

0initial2Depth (cm)

0.01
4

16 d
0.02

6

0.1
6

8

10
20

24

28

observed
6d
16 d
32

36

40

Temperature (°C)
図 4-5 実測温度分布と計算温度分布の比較。
Fig. 4-5 A comparison of observed and simulated soil temperature profiles.
714.3.3 液状水・水蒸気フラックス

土中水圧力分布と温度分布が得られたら，土中水圧力勾配，温度勾配による液状水，
水蒸気フラックスを定量的に計算することができる。図 4-6 に，6 日目，16 日目の液状水
フラックス分布と水蒸気フラックス分布を，2 つのフラックスを足し合わせた全フラックス分
布と共に示す。また図 4-7 に，液状水フラックス，水蒸気フラックスをさらにそれぞれ土中
水圧力勾配，温度勾配による 2 つの成分に分けた，4 成分のフラックス分布を示す。正の値は上向きの水分移動を，負の値は下向きの水分移動を示す。
既往の報告と同様に（Taylor and Cavazza, 1954; Nassar et al., 1992a），水蒸気は高温
側から低温側へ移動し，逆に液状水は低温側から高温側へ移動した。水蒸気フラックス

qv は，温度勾配による水蒸気フラックス qvT が主で(Scanlon et al., 2003)，土中水圧力勾配による水蒸気フラックスは，比較的大きな圧力勾配が生じる蒸発前線（後述，6 日目では深さ 9 cm 付近）と地表面を除いては，小さな値となった（図 4-7，6 日目）。また，低温のカラム下端からの液状水の上昇は，土中水圧力勾配による液状水移動 qLh で生じている（図 4-7）。温度勾配による液状水移動は qLT は，液状水量が多い下端付近で大きな値を持つ（図 4-7）。しかし，その値は qLh に比べ小さく，qLh と qLT の和で表される液状水フラッ
クス qL は上向き（正の値）を示す。
水蒸気フラックス qv と液状水フラックス qL は共に，水分量増加前線（図 4-4（a））において最大値を持つ。液状水量がカラム下端から徐々に増加し，前線が上方に移動するに従い，両フラックスの最大値の位置も上昇し，また両フラックスの最大値も増加した（図4-6）。また，全水分フラックス qTotal は常に下向き（負の値）を示した。つまり正味の水分移
動は高温側の上端から，低温側の下端へ向かって起こっていることが分かる（Jones andKohnke, 1952）。

4.3.4 表面境界条件

Milly（1984）や Noborio et al.（1996b）が，表面の水分フラックスをそれぞれ気象条件や熱収支を用いて計算しているのに対し，本計算では，カラム総重量の変化から求められた浸入フラックスを qTotal(0,t)を表面フラックスとして与えた。ここで，（4.22）式で表されているように，qTotal(0,t)は表面を通過する液状水フラックスと水蒸気フラックスの和で与えら
72

0
16 d

6d

Depth (cm)

2

4

6
qv
8

10
-0.04

qv

qL

qTotal
-0.02

qL

qTotal

0

-0.04

0.02

-0.02

0

0.02

0.04

Flux (cmd-1)

図 4-6 6 日後と 16 日後の全水分フラックス（qTotal），液状水フラックス（qL），水蒸気フラックス（qv）分布の計算値。
Fig. 4-6 Simulated profiles of the total water flux (qTotal), liquid water flux (qL), and vapor flux(qv) after 6 (left) and 16 (right) days.

0
16 d

6d

Depth (cm)

2
qLh
qLT
qvh
qvT
qTotal

4

6

8

10
-0.04

-0.02

0

-0.4

0.02

-0.2

0

0.2

0.4

Flux (cmd-1)

図 4-7 6 日後と 16 日後の全水分フラックス（qTotal），土中水圧力勾配による液状水フラックス（qLh），温度勾配による液状水フラックス（qLT），土中水圧力勾配による水蒸気フラックス（qvh），温度勾配による水蒸気フラックス（qvT）分布の計算値。
Fig. 4-7 Simulated profiles of the total water flux (qTotal), liquid water flux due to pressure
head (qLh) and temperature (qLT) gradients, and vapor flux due to pressure head (qvh) andtemperature (qvT) gradients after 6 (left) and 16 (right) days.

73

れる。土中水圧力勾配，温度勾配による液状水フラックス，水蒸気フラックスで表される
全水分フラックスの 4 成分は，（4.1）式を用いて求めることができる。4 成分のフラックスは，試料表面の土中水圧力勾配，温度勾配により，それぞれ違う向きで違う大きさを持つ。
土カラム中の温度勾配は実験開始後すぐに定常状態に達するため（図 4-5 では 0.1 日後），計算において，qTotal を水分浸入量の実測値に合せることができるのは，試料表面の土中水圧力変化のみである。
図 4-8（a）は，計算期間中の表面土中水圧力変化を示す。実験開始時の−15000000

cm から，2 日目までに急激に増加し，16 日目には−200000cm に達した。図 4-8（b）に，
表面フラックス qTotal(0,t)の 4 成分のフラックスの時間変化を示す。液状水フラックスの 2成分（qLh と qLT）は，水蒸気フラックスの 2 成分（qvh と qvT）に比べ，無視できるほど小さい値となった。これは液状水移動に対する透水係数（KLh ≈ 10−10 cm d−1 と KLT ≈ 10−8 cm2

K−1 d−1）が水蒸気移動に対する透過係数（Kvh ≈ 10−5 cm d−1 と KvT ≈ 1 cm2 K−1 d−1）に比べ，非常に小さいことによる。この水蒸気凝縮実験においては，下向きを示す温度勾配
による水蒸気移動と，上向きを示す土中水圧力勾配による水蒸気フラックスが支配的であると言える。

0.01

0

(b)
qvh

-400000

qLT

0

qLh

-0.01

-800000

qTotal-1200000

-0.02

qvT

Surface flux (cmd-1)

Surface pressure head (cm)(a)

observed inflow-1600000
0

4

8

0

12

4

8

12

-0.03
16

Time (d)図 4-8 （a）地表面の土中水圧力変化の計算値，（b）地表面の水分フラックス変化の計算値。
Fig. 4-8 Calculated pressure heads (a) and surface water fluxes (b) at the soil surface as afunction of time.

74

4.3.5 凝縮・蒸発速度

内部格子点における蒸発速度は，次に示すように（4.3）式を差分化することで求めることができる。

Ei = −jθ Lij − θ Lij −1∆t+qLij−−11 22 − qLij+−11 22∆z(4.30)

ここで，下付きの i は格子点の位置を示し，上付きの j および j−1 は，現在の，また 1 つ
前の時間レベルを示す。i = 1+1/4 と i = N-1/4 における蒸発速度は（4.31）式と同様に，表現することができる。ここで i = 1 と N はそれぞれ，試料の下端境界と上端境界を示す。
j1+1/ 4E=−

ENj −1 4 = −θ L1j+1 4 − θ L1j+−11 4∆tθ L Nj −1 4 − θ L Nj −−11 4∆t+

qL1j −1 2 − qL1j+−11 22+∆z / 2
qL Nj −1 2 − qL Nj −−1122∆z / 2(4.31)(4.32)

θL1+1/4 のような i = 1+1/4 および i = N−1/4 における値は，隣接する 2 つの格子点の値を線形補間すること求めることができる。ここで，qL1 と qLN は，式（4.22），式（4.23）に示されるように全フラックスが 0 である境界条件であるため，必ずしも 0 である必要はない。表面境界における蒸発速度 E1 は，E2 を基準として E1+1/4 を線形外挿することで求めることができる。また，同様に EN−1 を基準に EN−1/4 を外挿することで，下端境界における蒸発速度

EN を求めることができる。
図 4-9（a）は，（4.31）式から（4.32）式を用いて計算した蒸発速度分布である。ここで，負の蒸発速度は水蒸気凝縮速度を示す。蒸発速度の最大値は，6 日目では深さ 9 cmに生じ，16 日目には深さ 7 cm の位置に上昇した。この局所的に大きな蒸発速度により，水分量増加前線における最大の温度勾配による水蒸気フラックス qvT が生じる（図 4-7）。
図 4-9（b）には，（4.15）式から求めた 6 日目と 16 日目の相対湿度分布を示す。相対湿度
は，試料下端から水分量増加前線（図 4-3（a））の間はほぼ 1 を示す。前線より上方では，乾燥した領域であるため，地表面に近くなるほど相対湿度は小さくなる。
水蒸気凝縮は，蒸発速度が最大値を示す深さの直下から，つまり相対湿度が 1 となる深さから生じる。土中の気相に含まれる水蒸気量は相対湿度，つまり水蒸気密度によって一定量にきまっているため，上方からの水蒸気の凝縮が生じる。図 4-6 に示すように，
水蒸気はさらに下方へ降下し，カラム下端に達した時に，凝縮する。図 4-9（a）に示すよ
75

うに，最大の凝縮（最小の蒸発）はカラム下端で生じ，6 日目では−0.349 d−1，16 日目では−0.177 d−1 を示した。
土カラム中の総凝縮量は，凝縮速度分布（蒸発速度分布の負の領域）を深さ方向に積分することで得られる。図 4-10 は，積算凝縮量の時間変化を示す。土カラム中の総凝縮量は，カラム下端における凝縮量のみではなく，その他の部分における凝縮量も含む。
実験開始直後においては，カラム下端における凝縮が支配的である。その後は，下端以外のカラム内部における凝縮量が増加し，16 日目においては全凝縮量の 60 ％が内部で生じた（下端では残りの 40 %）。地表面からの水蒸気浸入量の測定値を，図 4-10 に同時に示す。総凝縮量が水蒸気浸入量の実測値よりも大きくなるのは，カラム内部で蒸発が生じ，生じた水蒸気がさらに凝縮するためである（図 4-9（a））。言い換えれば，総凝縮量と積算水蒸気浸入量の差が，カラム内部における蒸発量となる。カラム内部の蒸発量は，時間に関わらず，総凝縮量の 55 から 58 ％であった。

0

(b)

(a)

Depth (cm)

2

6d
16 d

6d
16 d

4

6

8

10
-0.15

-0.1

-0.05

0

0.05

0.8

Evaporation rate (d-1)

0.85

0.9

0.95

1

Relative humidity

図 4-9 （a）蒸発速度分布の計算値，（b）相対湿度分布。
Fig. 4-9 Simulated profiles of the evaporation rate (a) and relative humidity (b) after 6 and 16
days.

76

Cumulative condensation (cm)

0.8
Total
BottomInflow

0.6

0.4

0.2

0
0

4

8

12

16

Time (d)図4-10 土カラム中の総凝縮量（Total），カラム下端における凝縮量（Bottom）と地表面からの積算水蒸気浸入量（Inflow）の比較。
Fig. 4-10 Calculated cumulative condensation amounts in the entire soil column (Total) and at
the bottom part of the soil column (Bottom), compared with the cumulative water inflow
from the surface (Inflow) as a function of time.

4.3.6 水蒸気促進係数

水蒸気促進係数ηが，計算結果に与える影響については，液状水，水蒸気フラックス
分布（図 4-11（a））と（図 4-11（b））により示す。ここに結果は示さないが，水蒸気促進係数を無視した場合（η = 1），下向きの水蒸気フラックスを大きく過小評価し，カラムの上部において水蒸気の凝縮が生じた。a ≥ 5 の結果については，計算の水分分布には大きな違
いは見られなかったが，水蒸気促進係数が大きくなるほど下向きの水蒸気フラックスは大きくなった（図 4-11（a））。また，最大の蒸発速度と凝縮速度も，蒸気促進係数が大きくなるほど増加した（図 4-11（b））。本研究の水蒸気凝縮過程においては，蒸発と凝縮という 2 つの過程が競合し，打消しあうことで，水蒸気促進係数の水分分布に対する感度を小さくしたものと考えられる。また，砂質土の水蒸気促進係数がおよそ 10 であるというCass et al.（1984）の提案は，妥当であると言える。

77

0

(a)
2

Depth (cm)

(b)

16 d
a = 9.5
a = 15

4

6

qvT

qLh

8

10
-0.06

-0.03

0

-0.2

0.03

-0.1

0

0.1

0.2

Evaporation rate (d-1)

Flux (cmd-1)

図 4-11 （a）水蒸気促進係数のパラメータ a = 9.5 と a = 15 の液状水フラックスと水蒸気フラックスの比較，（b）蒸発速度分布の比較。
Fig. 4-11 A comparison of profiles of liquid water and vapor fluxes (a) and evaporation rates (b)simulated using different values of the parameter a in the enhancement factor, a = 9.5 (solidline) and a = 15 (dashed line).

0

(b)

(a)

observed
6d
16 d

Depth (cm)

2

4

6
16 d

16 d
8
6

6

initialinitial

10
0

0.4

0.8

1.2

0

Solute concentration, c (mmol cm-3)

0.002

0.004

0.006

0.008

Amount of solute, θc (mmol cm-3soil)

0.01

図 4-12 （a）溶液濃度分布の計算値，（b）単位土体積当たりの溶質量分布の計算値と実測値。
Fig. 4-12 Simulted solution concentration profiles (a) and a comparison of observed and simulated
solute concentration per unit volume of soil (b) after 6 and 16 days.

784.3.7 溶質分布

溶質は，液状水と共に移動するため，土中の溶質分布は液状水，水蒸気移動を反映
する（Gurr et al., 1952）。図 4-12 は，6 日目と 16 日目の溶質濃度分布の実測値と計算値を示す。図 4-12（a）は溶液濃度 c 分布を示し，図 4-12（b）は，単位土体積当たりの溶質
濃度分布θc を示す。初期溶液濃度は 0.93 mmol cm−3(water) で，初期溶質濃度は0.0042 mmol cm−3(soil)である。溶液濃度は，試料表面除いて初期値より低下した。これ
は，水蒸気凝縮により液状水量が増加した結果である（図 4-4（a））。特にカラム下端では，
水蒸気凝縮により溶質を含まない水が大量に生じるため，溶液濃度はほぼ 0 となった
（図 4-9（a））。また，水分量増加前線より上方においても，初期の溶液濃度よりも小さくなった。これは，実験開始初期に生じる凝縮による水分量増加によるものである（例えば深さ 5 cm において，θL = 0.0045 から 0.009）。溶液濃度分布は 2 ヶ所でピークを持った。1つはカラム上端で，液状水の上方への移動と試料表面における蒸発によるものである。
もう 1 つは，最大の蒸発が生じる深さ（図 4-9（a）），つまり水分量増加前線で，6 日目，16
日目でそれぞれ深さ 9 cm，6.5 cm である。
単位土体積当たりの溶質濃度θ c も，最大の蒸発が生じる深さにおいてピークを示した（図 4-12（b））。また，溶質は液状水と共に動くため，ピークの深さは時間と共に上方へ移動した。電気伝導度（EC）測定に基づく溶質濃度の実測値を図 4-12（b）に示す。異なる時間の EC 測定は，異なる土カラムにおいて行なわれたため，初期の溶質量は各土カラムで多少異なっていた。しかし，計算値は実測値のピーク位置をよく再現した。

4.4

まとめ

温度勾配下の砂カラム中の水蒸気，液状水移動（宮﨑，1976）を Philip and de Vries（1957）モデルを用いて定量的に評価した。水蒸気，液状水，熱，溶質同時移動について，修正した HYDRUS-1D を用いて数値計算を行なった（Philip and de Vries, 1957）。

Philip and de Vries モデルは砂カラム中の凝縮・蒸発過程をよく再現した。
不飽和透水係数を，実測の体積含水率分布に基づき評価した。不飽和水分移動特性関数に Fayer and Simmons（1995）のモデルを用い，間隙結合係数 l = 5.4 を用いたとき，計算の水分分布は実測値と一致した。また，Brooks and Corey（1964）のモデルを用79

いた場合は，小さな l の値を用いても，実測の水分分布を十分に再現することができなかった。このことから，凝縮過程においても，土中水圧力に対する不飽和透水係数が曲がる形状を示す Fayer モデルは，低水分領域の砂の不飽和透水係数を表現するのに適していることが裏付けられた。
全水分フラックスの 4 つの成分，土中水圧力勾配と温度勾配による液状水フラックスと水蒸気フラックス，について評価した。また，これらのフラックスを用いることで，蒸発速度および凝縮速度を評価した。水蒸気は，高温の地表面から砂カラム内へ浸入し，温度勾配により下方へ移動した。そして，低温のカラム下端において凝縮した。一方，凝縮により生じた液状水は，土中水圧力勾配により上方へ移動した。そして，上方へ移動した液状水は，相対湿度が 1 から減少する水分量増加前線において蒸発した。蒸発により生じた水蒸気は，試料表面から浸入した水蒸気と共に，再び下方へ移動した。水蒸気の一部は，相対湿度がおよそ 1 となる，蒸発速度がピークを示す深さの直下において凝縮した。最大の凝縮は低温のカラム下端で生じ，実験初期は総凝縮量の大部分が，16 日目には 40 %がカラム下端で生じた。このように液状水と水蒸気は，凝縮と蒸発を繰り返しな
がら，カラム下端と水分量増加前線の間を循環していることが明らかとなった。温度勾配
による液状水フラックスは，水分量が大きいカラム下部において大きな値を持ち，下向きを示した。しかし，その絶対値は上向きを示す土中水圧力勾配による液状水フラックスよりも小さく，正味の液状水フラックスは上向きを示した。また，土中水圧力勾配による水蒸気フラックスは，前線付近で大きな値を持ち，上向きを示したが，その絶対値は温度勾配による水蒸気フラックスに比べて小さかった。
さらに，Cass et al.（1984）のモデルで示される水蒸気促進係数が，計算結果に与える影響についても調べた。Cass モデルのパラメータ a が 5 より大きい時，計算の水分分布に大きな違いは見られなかった。大きな水蒸気促進係数を用いると，下向きの水蒸気フラックスと凝縮速度が増加した。しかし，同時に上向きの液状水フラックスと，前線における蒸発速度も増加した。これらの 2 つの逆の過程が競合することで，水蒸気促進係数の水分分布に対する感度が小さくなったと考えられた。砂質土の水蒸気促進係数がおよそ10 であるという Cass et al.（1984）の提案は，妥当であると言えた。

80

第5章

結論

本研究では，Mualem の間隙径分布モデルにより水分保持曲線‐不飽和透水係数連結モデルで表される不飽和水分移動特性の評価と，その関数型の適用性，特に低水分領域の不飽和透水係数について検討することを主目的とした。そこで，実際の不飽和水分移動過程の土中水圧力変化に基づき，逆解析手法によって不飽和水分移動特性を推定する蒸発法に注目した。まず，仮想試料を用いた数値実験，及びシルト質土を用いた室内実験により，蒸発法により不飽和透水係数を推定する上での最適な実験
条件について明らかにした。次に，蒸発法に基づいた蒸発実験と逆解析を行ない，Fayer モデル，Durner モデルの砂質土と団粒構造を持つ黒ボク土への適用性を明らかにした。さらに，低水分領域の砂質土への水蒸気凝縮実験について数値解析を行ない，Fayer モデルの適用性について検証を行なった。
第 1 章では，数値計算による不飽和土中の水分移動予測における，不飽和水分移
動特性関数の重要性について述べ，その既往の研究について，Mualem（1980）の間隙径分布モデルを中心に整理した。そして，閉形式解を用いて水分保持曲線から
不飽和透水係数を導く方法は非常に簡便である一方，多くの場合，導かれた不飽和透水係数について十分な検討がなされないまま，数値計算に用いられている現状を指摘した。
第 2 章では，不飽和水分移動特性を評価する方法として蒸発法に注目し，水分保持曲線‐不飽和透水係数連結モデルの各種パラメータを推定するために有効な蒸発実
験の実験条件について明らかにした。特に不飽和透水係数のパラメータである間隙結合係数 l の推定を中心に，数値実験及び室内実験から整理した。シルト質土カラムからの蒸発過程を対象とした数値実験を行ない，VG モデルのパラメータについて，土中水
圧力に関する 2 次元応答曲面を示し，それぞれのパラメータの感度やパラメータ間の相関に基づき，蒸発法の最適条件について検討した。その結果，不飽和透水係数を決める間隙結合係数 l の感度は地表面に近い位置の土中水圧力に対してほど高いことから，
土中水圧力の時間変化が大きい表層のデータほどパラメータが収束しやすく，不飽和
透水係数の信頼度の高い推定が可能なことが示唆された。さらに，複数の観測位置のデータを加えて逆解析を行なうことは，収束の条件が改善されるので望ましいと考えられ81た。また，蒸発を促進するほど，また試料長を長くするほど，試料内の土中水圧力差が増大することから，不飽和水分移動特性パラメータの感度が高くなり，パラメータの収束が容易になる，すなわち不飽和水分移動特性の推定の信頼性を高めることができることを明らかにした。
さらに，藤森シルトを用いて蒸発法を行ない，異なる試料長と観測点の位置と数について検討した。試料長が長いほど，また地表面に近い観測点ほど，推定したパラメータの標準偏差σ が小さくなり，不飽和水分移動特性の推定の信頼度が向上することが示された。また 2 地点の土中水圧力データを用いることは，不飽和透水係数の推定において望ましいことが示された。
第 3 章では，広く用いられる VG モデルでは表現できないとされる砂質土と団粒土を対象に不飽和水分移動特性の関数型，特に不飽和透水係数について検討した。

VG モデルの水分保持曲線関数を低水分領域の吸着保水を考慮して修正した Fayer モ
デルと，2 つの VG モデルを線形結合して階段状の水分保持曲線関数を示した Durner
モデルを用いて，砂質土の鳥取砂丘砂と団粒構造の発達した黒ボク土を対象に蒸発法を行ない，パラメータの推定を行なった。そして，蒸発過程の土中水圧力の時間変化の実測値と計算値を比較，また水分保持曲線の実測値と推定値を比較し，砂丘砂には

Fayer モデルと Durner モデルを，黒ボク土には Durner モデルを用いることで，低水分領域の不飽和透水係数を表現できることを明らかにした。
広い圧力範囲の水分保持曲線に対して適合性の高い Fayer モデルや Durner モデル
において，さらに間隙結合係数 l を推定パラメータとして用いることにより，Mualem モデ
ルによる水分保持曲線‐不飽和透水係数連結モデルの利便性を失うことなく，飽和から乾燥領域までの広い圧力範囲の不飽和透水係数が表現できた。また，蒸発法による水分移動特性推定の信頼範囲は，測定された圧力範囲とされている。しかし，砂質土に対する Fayer モデルや黒ボク土に対する Durner モデルのように，低水分領域の実測値を正確に表現できていれば，蒸発法で得られる不飽和透水係数は，測定した圧力範囲を超えても信頼度が高いと考えられた。そのため，VG モデルの低水分領域を修正した

Fayer モデルと Durner モデルは，実用上，非常に有用なモデルであると考えられた。
第 4 章では，水蒸気凝縮，蒸発が生じる非等温で乾燥した砂質土中の水分移動実験の水分量変化の実測値と，Fayer モデルを用いた計算値を比較することで，第 3 章で示された低水分量の砂質土に対する Fayer モデルの不飽和透水係数の適用性について82

検証を行なった。Fayer モデルを用いた計算値は，高温で湿った大気から，水蒸気が温度勾配により乾燥した砂中へ浸入し，試料下端で凝縮する過程を表し，実測の水分分
布をよく再現した。このことからも，Fayer モデルにより，砂質土の低水分領域の不飽和透水係数を表現することができることが検証された。
さらに，Fayer モデルにより正しく評価された不飽和水分移動特性を用いて，計算を行なうことで，試料中の熱移動，溶質移動，4 成分の水分フラックス（土中水圧力勾配，温度勾配による液状水，水蒸気フラックス），蒸発および凝縮速度を定量的に評価することができ，また，試料中を液状水と水蒸気が蒸発，凝縮の相変化を繰り返しながら循環する過程を明らかにした。

以上のように，本研究では，不飽和水分移動特性について，室内実験による不飽和
水分移動の実測値と，数値計算による予測値の比較を行なうことで，大きさを評価
し，砂質土や黒ボク土への適用性を明らかにした。不飽和水分移動特性を評価する方法として，蒸発法の適切な実験条件を整理することで，より正確な不飽和水分移
動特性の把握が可能となったと考える。また，水分保持曲線の実測値に対し広水分量範囲で適合する Fayer モデルや Durner モデルといった水分保持曲線関数を用いることで，Mualem モデルの閉形式解により，低水分領域に対してもある程度正し
い不飽和透水係数を表現できることは，数値計算において非常に有益な知見であると言える。

これらの本研究で得られた知見により，不飽和土中の水分移動をはじめ，熱移動，汚染物質などの溶質移動の数値計算によるより正確な予測が可能になると考える。
本研究では第 4 章において，砂質土中の水分移動に対する Fayer モデルの検証しか行なわなかったが，同様に黒ボク土中の水分移動に対する Durner モデルの検証を行なうことは今後の課題として残す。また，降雨・蒸発といった再分布過程を伴うより現実に近い水分移動を正しく予測するためにも，ヒステリシスモデルを含む
不飽和水分移動特性関数についても検討していく必要がある。

83

記号

SymbolDescription

Dimension

a
b
c
Cs
Cp
Cv
Cw
D
Da
Dw
E
E
Ecum
fc
g
GwT
h
hA
hi
hk
hm
hs
Hr
K, KLh
Kk
KLT
Kr
Ks
Kvh
KvT
l
L
L0
Lw
m
m
M
n
nj
p*j(ti)pj(ti)
qevap
qL
qLh
qLT

水蒸気促進係数のパラメータ
熱伝導率のパラメータ（b1，b2，b3）
溶液濃度
固相の体積熱容量
土の体積熱容量
気相の熱容量
液相の体積熱容量
土中の水蒸気拡散係数
大気中の水蒸気拡散係数
水中の分子拡散係数
可能蒸発速度
蒸発速度（第 4 章）
積算蒸発量
粘土分率
重力加速度（= 9.81 m s−2）
表面張力の温度依存性を修正する促進係数
土中水の圧力水頭
許容最小圧力水頭
初期土中水圧力水頭
任意の土中水の圧力水頭
炉乾時の土中水の圧力水頭
空気侵入値
相対湿度
不飽和透水係数
任意の土中水圧力 hk における不飽和透水係数
温度勾配による液状水移動の不飽和透水係数
比透水係数
飽和透水係数圧力勾配による水蒸気透過係数
温度勾配による水蒸気透過係数
間隙結合係数
試料長
単位体積当たりの水の蒸発潜熱
単位重量当たりの蒸発潜熱（J kg−1）
水分保持曲線の形状を与えるパラメータ（m2）
観測地点の数（目的関数）
水の分子量（= 0.018015 kg mol−1）
水分保持曲線の形状を与えるパラメータ（n2）観測地点 j における土中水圧力データの数
観測地点 j の ti における土中水圧力の計算値
観測地点 j の ti における土中水圧力の実測値
蒸発速度
液状水フラックス
土中水圧力勾配による液状水フラックス
温度勾配による液状水フラックス
84

M L−3
ML T K−1
M L−1 T−2 K−1
M L−1 T−2 K−1
M L−1 T−2 K−1
L2 T−1L2 T−1L2 T−1
L T−１
T−1
L
−1 −2

L T−2
L
L
L
L
L
L
L T−1
L T−1
L2 K−1 T−1
L T−1
L T−1
L2 K−1 T−1
L
M L−1 T−2
L2 T−2
M mol−1
L
L
L T−１
L T−１
L T−１
L T−１

qTotal
qv
qvh
qvT
R
Se
t
T
Tabs
ti
w1
w
z

α
χ(h)
φ
Φ

γ
γ0
Γ

η
λL
λ(θL)
v

θ, θL
θa
θa
θave
θk
θm
θr
θs
θsolid
θv
ρb
ρvs
ρw
σ2
τa
τw

全水分フラックス
水蒸気フラックス
土中水圧力勾配による水蒸気フラックス
温度勾配による水蒸気フラックス
気体定数（= 8.341 J mol−1 K−1）
有効飽和度（Se1, Se2）
時間
相対温度（ C）
絶対温度
測定時刻
Durner モデルの重み係数（w1 + w2 = 1）
重み係数（wp，wθ，wE，wθ (h)）
上向き正の空間座標
水分保持曲線の形状を与えるパラメータ低水分領域における水分保持曲線の補間関数
平方残差
目的関数
土中水の表面張力
25 °C における表面張力（= 71.89 g s−2）
Mualem モデルの積分（Γs(h)，Γmax）
水蒸気促進係数
分散長
土の熱伝導率
重み係数
体積含水率
水分保持曲線の形状を与えるパラメータ
土中の気相率（第 4 章）
蒸発実験終了時の平均体積含水率
任意の土中水圧力 hk における体積含水率
仮想の飽和体積含水率
残留体積含水率
飽和体積含水率
固相率
水蒸気の体積分率（液状水等量）
乾燥密度
飽和水蒸気密度
液状水密度
標本分散（σj2）
気相の屈曲度
液相の屈曲度

85

L T−１
L T−１
L T−１
L T−１
ML2T−2mol−1K−1
T
K
K
T
L
L−1

M T−2
M T−2
L
MLT−3K−1
L3 L−3L3 L−3L3 L−3L3 L−3L3 L−3L3 L−3L3 L−3L3 L−3L3 L−3
L3L−3
M L−3
M L−3
M L−3

参考文献
（1）
（2）
（3）
（4）
（5）

（6）
（7）
（8）
（9）
（10）
（11）
（12）
（13）
（14）
（15）
（16）
（17）
（18）
（19）

（20）

Brooks, R.H., and Corey, A.T. (1964): Properties of porous media affecting fluid flow. J.
Irrig. Drainage Div., ASCE Proc. 72(IR2): 61-88.
Burdine, N.T. (1953): Relative permeability calculations from pore-size distribution data.
Trans. Am. Inst. Min. Metall. Pet. Eng., 198: 71-77.
Cahill, A.T. and Parlange, M.B. (1998): On water vapor transport in field soils. Water
Resour. Res., 34: 731-739.
Campbell, G.S. (1985): Soil Physics with BASIC. Elsevier, New York.
Campbell, G. S. and Shiozawa, S. (1992): Prediction of hydraulic properties of soils using
particle size distribution and bulk density data, in International workshop on indirect
methods for estimating the hydraulic properties of unsaturated soils, University ofCalifornia Press, Berkeley, 1992.
Carsel, R.F. and Parrish, R.S. (1988): Developing joint probability distributions of soil
water retention characteristics. Water Resour. Res, 24: 755-769.
Cass, A., Campbell, G.S., and Jones, T.L. (1984): Enhancement of thermal water vapor
diffusion in soil. Soil Sci. Soc. Am. J., 48: 25-32.
Celia, M., Bououtas, E.T., and Zebra, R.L. (1990): A General mass-conservative numericalsolution for the unsaturated flow equation. Water Resour. Res., 26: 1483-1496.
Childs, E. C. and Collis-George, N. (1950): The permeability of porous materials. Proc. Roy.
Soc. London, Ser. A, 201: 392-405.
Chung, S.O. and Horton, R. (1987): Soil heat and water flow with a partial surface mulch.
Water Resour. Res., 23: 2175-2186.
de Vries, D.A. (1958): Simultaneous transfer of heat and moisture in porous media. Trans.
Am. Geophys. Union, 39: 909-916.
Durner, W. (1994): Hydraulic conductivity estimation for soils with heterogeneous pore
structure. Water Resour. Res., 30: 211-223.
Eching, S.O. and Hopmans, J.W. (1993): Optimization of hydraulic functions from transient
outflow and soil water pressure data. Soil Sci. Soc. Am. J., 57: 1167-1175.
Farlow, S.J. (1993): Partial differential equations for scientists and engineers. pp. 19-26,
Dover, New York.
Fayer, M.H. and Simmons, C.S. (1995): Modified soil water retention functions for all
matric suctions. Water Resour. Res., 31: 1233-1238.
藤井克己・中野政詩（1984）：ベントナイト吸着水の化学ポテンシャルについて．農業土
木学会論文集，112：43-53．
Fujimaki, H. and Inoue, M. (2003): A transient evaporation method for determining soil
hydraulic properties at low pressure, Vadose zone J., 2: 400-408.
Gurr, C.G., Marshall, T.J., and Hutton, J.T. (1952): Movement of water in soil due to atemperature gradient. Soil Sci., 74: 335-345.
Hopmans, J.W., Šimůnek, J., Romano, N., and Durner, W. (2002): Inverse
method-Evaporation method, Methods of soil analysis Part 4. pp. 978-981, ASA and SSSA.,
Madison.
井上光弘（2004）：不飽和地盤の挙動と評価．pp．202-204，地盤工学会，東京．
86

（21）

（22）
（23）

（24）
（25）
（26）
（27）
（28）

（29）
（30）
（31）

（32）
（33）
（34）
（35）
（36）
（37）
（38）

（39）

Inoue, M., Šimůnek, J., Hopmans, J.W., and Clausnitzer, V. (1998): In situ estimation of soil
hydraulic functions using a multistep soil-water extraction technique. Water Resour. Res.,
34: 1035-1050.
Jones, H.E. and Kohnke, H. (1952): The influence of soil moisture tension on vapor
movement of soil water. Soil Sci. Soc. Am. Proc., 16: 245-248.
Jury, W.A. and Horton, R. (2006) 土壌物理学−土中の水・熱・ガス・化学物質移動の基
礎と応用−(取出伸夫 監訳: 井上光弘・長裕幸・西村拓・諸泉利嗣・渡辺晋生訳),築地書店,東京.
Klute, A. and Dirksen, C. (1986): Hydraulic conductivity and diffusivity-Laboratory
methods. Methods of soil analysis Part 1. pp, 687-734, ASA and SSSA., Madison.
Kool, J.B., Parker, J.C. (1987): Development and evaluation of closed-form expressions forhysteretic soil hydraulic properties. Water Resour. Res., 23(1): 105-114.
Kool, J.B., Parker, J.C., and van Genuchten, M.Th. (1987): Parameter estimation forunsaturated flowand transport models-A review. J. Hydrol., 91: 255-293.
Kosugi, K. (1996): Lognormal distribution model for unsaturated soil hydraulic properties.
Water Resour. Res., 32: 2697-2703.
小杉賢一朗（2007）：古典を読む：Y. Mualem 著「不飽和多孔質体の透水係数を推定する
新たなモデルについて」ならびに M. Th. van Genuchten 著「不飽和土壌の透水係数を推定する閉形式解について」，土壌の物理性，106: 105-112.
Marquardt, D.W. (1963): An algorithm for least-squares estimation of nonlinear parameters.
SIAM J. Appl. Math, 11: 431-441.
Millington, R.J. and Quirk, J.M. (1961): Permeability of porous solids. Trans. Faraday Soc.,
57: 1200-1207.
Milly, P.C.D. (1982): Moisture and heat transport in hysteric, inhomogeneous porous
media: A matric head-based formulation and a numerical model. Water Resour. Res., 18:489-498.
Milly, P.C.D. (1984): A simulation analysis of thermal effects on evaporation. Water Resour.
Res., 20: 1087-1098.
宮﨑 毅（1976）：温度勾配下の砂における水蒸気の凝縮と移動．農業土木学会論文集，
61: 1-8.
Mualem, Y. (1976): A new model for predicting the hydraulic conductivity of unsaturated
porous media. Water Resour. Res., 12: 513-522.
中野政詩・宮﨑 毅・塩沢 昌・西村 拓（1995）：土壌物理環境測定法．pp．115-136，東京大学出版会，東京．
Nassar, I.N. and Horton, R. (1989): Water transport in unsaturated nonisothermal salty soil:
II. Theoretical development. Soil Sci. Soc. Am. J., 53: 1330-1337.
Nassar, I.N. and Horton, R. (1992): Simultaneous transfer of heat, water, and solute inporous media: I. Theoretical development. Soil Sci. Soc. Am. J., 56: 1350-1356.
Nassar, I.N. and Horton, R. (1997): Heat, water, and solute transfer in unsaturated porous
media: I. Theory development and transport coefficient evaluation. Transport in Porous
media., 27: 17-38.
Nassar, I.N., Globus, A.M., and Horton, R. (1992)a: Simultaneous soil heat and watertransfer. Soil Sci., 154: 465-472.
87

（40）

（41）
（42）

（43）

（44）

（45）
（46）
（47）
（48）

（49）
（50）
（51）
（52）
（53）

（54）

（55）

（56）

Nassar, I.N., Horton, R., and Globus, A.M. (1992)b: Simultaneous transfer of heat, water,and solute in porous media: II. Experiment and Analysis. Soil Sci. Soc. Am. J., 56:
1357-1365.
Nimmo, J.R. and Miller, E.E. (1986): The temperature dependence of isothermal moisturevs. potential characteristics of soils. Soil Sci. Soc. Am. J., 50: 1105-1113.
Noborio, K., McInnes, K.J., and Heilman, J.L. (1996)a: Two-dimensional model for water,
heat, and solute transport in furrow-irrigated soil: I. Theory. Soil Sci. Soc. Am. J., 60:
1001-1009.
Noborio, K., McInnes, K.J., and Heilman, J.L. (1996)b: Two-dimensional model for water,
heat, and solute transport in furrow-irrigated soil: II. Field Evaluation. Soil Sci. Soc. Am. J.,
60: 1010-1021.
Parkhurst, D.L. and C.A.J.Appelo (1999): The PHREEQC Software Package forSimulating Chemical Reactions and Transport Processes in Natural or Polluted water,
Version 2.0. U.S. Department of the Interior Bruce Babbitt, Secretary.
Penman, H.L. (1940): Gas and vapor movement in soil. I. J. Agr. Sci., 30: 437-462.
Philip, J.R., and de Vries, D.A. (1957): Moisture movement in porous media undertemperature gradients. Eos Trans. AGU, 38(2): 222-232.
Priesack, E. and Durner, W. (2006): Closed-form expression for the multi-modal
unsaturated conductivity function. Vadose zone J., 5: 892-898.
Rassam, D., Šimůnek, J., and van Genuchten, M.Th. (2004): HYDRUS-2D による土中の不飽和流れ計算 (取出伸夫・井上光弘 監訳), 農業土木学会土壌物理研究部会 HYDRUS
グループ, 東京．
Rossi, C. and Nimmo, J.R. (1994): Modeling of soil water retention from saturation to ovendryness. Water Resour. Res., 30: 701-708.
斎藤広隆・坂井 勝・Jiri Šimůnek・取出伸夫（2006）：不飽和土中の水分移動モデルにおける境界条件，土壌の物理性，104: 63-73.
Saito, H., Šimůnek, J., and Mohanty, B.P. (2006): Numerical analysis of coupled water,
vapor, and heat transport in the vadose zone. Vadose Zone J., 5: 784-800.
坂井 勝・取出伸夫（2007）：蒸発法における不飽和水分移動特性推定の最適条件の検討，土壌の物理性，106: 33-46.
Šimůnek, J., van Genuchten, M.Th., and Šejna, M. (2006): The HYDRUS software package
for simulating two- and three dimensional movement of water, heat, and multiple solutes invariably saturated media. Version 1.0. HYDRUS Software Series 1, Dep. of Environmental
Sciences, Univ. of California Riverside, Riverside, CA.
Scanlon, B., Keese, K., Reedy, R.C., Šimůnek, J. and Andraski, B.J. (2003): Variations inflow and transport in thick desert vadose zones in response to paleoclimatic forcing (0−90kyr): Field measurements, modeling, and uncertainties. Water Resour. Res., 7, 1179, doi:
10.1029/2002WR001604. 3.1-3.7.
Šimůnek, J., van Genuchten, M.Th., and Šejna, M. (2005): The HYDRUS-1D software
package for simulating the movement of water, heat, and multiple solutes in variably saturated
media. Version 3.0. Technical manual, PC Progress, Prague, Czech Republic.
Šimůnek, J., Wendroth, O., and van Genuchten, M.Th. (1998): Parameter estimation
analysis of the evaporation method for determining soil hydraulic properties. Soil Sci. Soc.
88

（57）
（58）
（59）
（60）

（61）
（62）
（63）

（64）
（65）

Am. J., 62: 894-905.
Taylor, S.A. and Cavazza, L. (1954): The movement of soil moisture in response totemperature gradients. Soil Sci. Soc. Am. Proc., 18: 351-358.
徳本 家康・取出 伸夫・井上 光弘（2005）：団粒構造を持つ黒ボク土の溶質分散について，水文・水資源学会誌，18 (4): 401-410.
Tuller, M. and Or, D. (2001): Hydraulic conductivity of variably saturted porous media.
Water Resour. Res., 37: 1257-1276.
van Dam, J.C., Stricker, J.N.M., and Droogers, P. (1994): Inverse method for determining
soil hydraulic functions from multi-step outflow experiments. Soil Soc. Am. J., 58:
647-652.
van Genuchten, R. (1978): Calculating the unsaturated hydraulic conductivity with a newclosed-form analytical model. Res. Rep. 78-WR-08. Princeton Univ., Princeton, USA.
van Genuchten, M.Th. (1980): A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Soc. Am. J., 44: 892-898.
van Genuchten, M.Th., Leij, F., and Yates, S. (1991): The RETC code for quantifying thehydraulic functions of unsaturated soils. Technical Report EPA / 600/ 2-91/065, US
Environmental Protection Agency.
Vogel, T. and Císlerová, M. (1988): On the reliability of unsaturated hydraulic conductivity
calculated from the moisture retention curve. Transport in porous media, 3: 1-15.
Wendroth, O., Ehlers, W., Hopmans, J.W., Kage, H., Halbertsma, J., and Wosten, J.H.M.
(1993): Reevaluation of the evaporation method for determining hydraulic functions in
unsaturated soils. Soil Sci. Soc. Am. J., 57: 1436-1443.

89謝辞
本論文は，平成 16 年 4 月に三重大学生物資源学研究科に入学してから，平成 19 年

9 月に卒業するまでの間に，著者が土壌圏循環学研究室で行なった研究成果をまとめたものである。その間，多くの方々のご指導やご支援を頂いたので，ここに感謝を表したい。
取出伸夫教授には，三重大学に誘って頂いたことから始まり，現象に対する鋭い視点
からのご指導，海外生活を踏まえた幅広い人生観，さらにはマラソンや水泳など，多くの貴重な経験を頂いた。渡辺晋生准教授には，様々なアドバイスを頂くとともに，研究者としてのお手本を示して頂いた。小畑仁教授，亀岡孝治教授，清沢秀樹准教授，また，東京大学農学生命科学研究科の宮﨑毅教授には，お忙しい中審査員として適切なご意
見，ご指導を頂いた。
陳代文君とは，同じ博士後期課程の学生として，お互いに励まし合いながら研究を進めることができた。野口淳平君，鈴木康生君，和気朋己君をはじめ，多くの後輩達とは，研究室の設立に始まり，実験の手伝いをして頂いたり，また日々の運動をしたりと，刺激を受けながら共に楽しい時間を過ごす事ができた。
カリフォルニア大学リバーサイド校の Jirka Šimůnek 教授には，HYDRUS を用いた数値計算において，多くのアドバイスを頂いた。徳本家康君には，黒ボク土についての貴重なデータを提供して頂いた。また，武藤由子さん，中野恵子さん，斎藤広隆さん，辻田
あずささんには，多くの暖かい励ましの言葉を頂いた。
最後に，私のわがままを寛大に受け入れてくれ，いつでも理解を示してくれた両親に感謝の意を表したい。

90

要旨
土中の水分移動の定量的な予測のために，土壌物理や水文学の分野では，水分不飽和土中の水分移動式としてのリチャーズ式に対する数値解法が長年にわたり研究されており，その集大成として不飽和土中の水分・溶質移動に対する汎用プログラムが開発されている。特に本研究で用いた HYDRUS は，最も広く利用されている汎用プログラ
ムのひとつである。汎用プログラムによる精度の高い予測には，土の水分保持曲線と不飽和透水係数を与える不飽和水分移動特性の正確な評価が最も重要である。水分移動
特性は，Mualem の間隙分布モデルにより水分保持曲線と不飽和透水係数を同じパラメータで表現する水分移動特性関数が主流である。水分移動特性関数では，多くの土に対して比較的適合が良く滑らかな関数型を与える van Genuchten（VG）モデルが広く用いられているが，低水分領域で不飽和透水係数を過小評価すること，また粗粒な砂や団粒構造を持つ土などは十分に表現できないことが指摘されている。そのため，低水分領域の吸着保水を考慮した Fayer モデルや，団粒内外の水分保持を異なる土性の重ね合
わせと見なす Durner モデルなど，VG モデルを修正したモデルが提唱されている。しかし，そうした修正モデルを不飽和透水係数の観点からの検証した研究例は少ない。一方，
汎用プログラムの発展に伴い，非定常な水分移動過程における土中水圧力などの実測値から，水分移動特性関数のパラメータを逆解析する手法が注目されている。マルチス
テップ法は代表的な逆解析法であるが，下端フィルターの評価が困難なことが指摘されており，より単純な装置で測定が可能な蒸発法が提唱されている。しかし蒸発法に対しては，まだ十分な検討が行なわれていないのが現状である。
本研究では，不飽和透水係数の推定法として蒸発法に注目し，信頼性の高い推定のための最適な実験条件を検討した。その上で，水分移動特性関数としての VG モデル，

Fayer モデル，Durner モデルを蒸発法に適用し，主に不飽和透水係数の観点から比較検討を行ない，それぞれのモデルの特性について検討した。さらに，砂丘砂の水分凝縮実験に対し，液状水，水蒸気，熱，溶質の移動式を適用することにより，砂丘砂に対する Fayer モデルの適用性を検討すると共に，水・溶質・熱の同時移動現象の解明を試みた。なお，本研究は，HYDRUS の作者である U.C. Riverside の Šimůnek 博士の協力を得ながら，水分移動汎用プログラム HYDRUS の改良の一環として行なった。
まず，蒸発法によるパラメータ推定のための適切な実験条件について検討するため
91

に，数値実験及び，実測のデータに対するパラメータの推定を行ない，特に不飽和透水係数の間隙結合係数 l の推定を中心に整理した。数値実験では，蒸発強度，試料長，土中水圧力測定の位置と数を変化させ，パラメータの相関と蒸発過程の土中水圧力変化に対する感度を調べた。また，藤森シルトに対して蒸発法を適用し，実際にパラメータを推定した。蒸発を促進すること，試料長を長くすること，表面付近を含む 2 地点の圧力を測定することで，推定の信頼度を向上できることを明らかにした。
続いて，砂質土の鳥取砂丘砂と団粒構造の発達した黒ボク土を対象に蒸発法を行な
い，VG モデル，Fayer モデル，Durner モデルのパラメータを推定した。そして，蒸発過程の土中水圧力変化の実測値と計算値を比較，また水分保持曲線の実測値と推定値を比較し，砂丘砂には Fayer モデルと Durner モデルを，黒ボク土には Durner モデルを用いることで，高水分領域と低水分領域の 2 つの異なる透水形態（砂丘砂は間隙流と膜流，
黒ボク土は団粒間と団粒内の間隙流）による不飽和透水係数を表現できることを明らかにした。そして，水分保持曲線を適切に表現したモデルを用いた上で，間隙結合係数 lを蒸発法により推定することで，不飽和透水係数の推定に大きな改善が図れることを示した。また，VG モデルについては，不飽和透水係数を表現できる適用可能圧力範囲を便宜的に提案した。
さらに，水蒸気凝縮，蒸発が生じる非等温で乾燥した砂質土中の水分移動実験に対して数値計算を行ない，水分量変化の実測値と，Fayer モデルを用いた計算値を比較することで，砂質土に対する Fayer モデルの不飽和透水係数の適用性について検証を行なった。そして，Fayer モデルにより評価した水分移動特性を用いて，水分凝縮過程の熱移動，溶質移動，4 成分の水分フラックス（土中水圧力勾配，温度勾配による液状水，水
蒸気フラックス），蒸発および凝縮速度を定量的に評価し，試料中を液状水と水蒸気が蒸発，凝縮の相変化を繰り返しながら循環する過程，またそれに伴う熱，溶質移動を明らかにした。
本研究では，水分保持曲線の実測値に対し広水分量範囲で適合するモデルを用い，蒸発法により間隙結合係数 l と共にパラメータを推定することで，Mualem モデルにより不飽和透水係数も広範囲でより正しく表現できることを明らかにした。これは，不飽和土中の物質移動予測の実用面の上からも意義深く，予測精度の向上へ寄与できるものと考
えられる。

92

Soil water hydraulic functions for water flow in variable saturated soils
Masaru Sakai
Graduate school of Bioresources, Mie University
Abstract
Analyzing water flow in unsaturated (vadose) zone plays important roles in assessing
movement of contaminants, fertilizers and pesticides to groundwater aquifer and planning
irrigation schedules for arid or semi-arid regions. In these decades, computer programs
numerically solving the highly nonlinear equations have been developed for predicting
water and solute movement in the vadose zone and for analyzing laboratory or field
experiments. For simulating water flow in the unsaturated soil accurately, evaluating
unsaturated hydraulic property consisting of unsaturated hydraulic conductivity and soil
water retention curve is necessary. In this study, we focus on unsaturated hydraulic property
in specific, and evolve the HYDRUS-1D code simulating water flow in unsaturated soils.
Since most methods measuring unsaturated hydraulic conductivity remain relatively time
consuming and costly, and are limited to relatively narrow range of water content,
unsaturated hydraulic conductivity function is often estimated from soil water retentionfunction based on statistical pore-size distribution model of Mualem. One simple laboratory
method for both retention curve and unsaturated hydraulic conductivity is evaporation
method, optimizing parameters of soil water hydraulic function as numerically generated
data fit to experimentally obtained pressure heads in evaporation process. In this study,
optimum experimental conditions for predicting unsaturated hydraulic properties using theevaporation method were investigated. Then applicability of new unsaturated hydraulicproperties to dune sand and aggregated Andisol was investigated using the evaporation
method.
Sensitivities and correlations for two selected parameters including the pore-connectivity
coefficient, l, which describes the tortuosity factor for the unsaturated hydraulicconductivity of the van Genuchten-Mualem model were evaluated with the response
surface in terms of the pore-water pressure heads inside the soil column. Higher maximum
potential rate of evaporation, longer column length, and pressure heads measured close to
the soil surface were found to be more valuable for the parameter estimation.
The Fayer model well described unsaturated hydraulic conductivity for the dune sandreflecting pore water retention for higher water contents and absorbed water retention on
the surface of soil particles for lower water contents. The Durner model was suitable to theAndisol whose hydraulic properties were strongly affected with the distinct the aggregatedstructure. The Mualem model can predict unsaturated hydraulic conductivity for a wide93

range of soil water pressure when an appropriate model is used for the water retention curve
and the pore-connectivity coefficient, l is estimated. Furthermore, it was also suggested that
we need to use an appropriate soil water retention function to predict unsaturated hydraulicconductivity for much smaller pressure heads than the minimum pressure head observed
during the evaporation experiment.
Furthermore, water vapor condensation process to dry sand was investigated with using
Fayer model, which describes unsaturated hydraulic conductivity of sand well, andsimultaneous mavement of liquid water, water vapor, heat, and solute movement were
evaluated quantitatively.

94