ゎࠉㄝ

J. Jpn. Soc. Soil Phys.
土壌の物理性LecturesNo. 154, p.19 ∼ 27 (2023)

一般化透水モデルによる不飽和透水係数の閉形式解
関 勝寿 1 · 取出伸夫 2

Closed-form equations of unsaturated hydraulic conductivity
with general hydraulic conductivity model
Katsutoshi SEKI1 and Nobuo TORIDE2

1. はじめにを導いた．ここで p は間隙結合係数（tortuosity factor）とよばれるパラメータであり，多くの土壌における平均土壌の不飽和水分移動を解析する際には，マトリック
ポテンシャル水頭 h と体積含水率 θ との関係をあらわす
水分保持関数 θ (h) と，h および θ に依存する土壌の不飽和透水係数 K をそれぞれ関数として表現することが重要である．
水分保持関数 θ (h) は，飽和体積含水率 θs と残留体積
含水率 θr によって定義される飽和度 S によって，S(h) =(1)

圧を h > 0 とする表記（いわゆるサクション）を用いる．

1950 ∼ 1970 年代に水分保持曲線から不飽和透水係数を導出するモデルが多く提唱された．Burdine（1953）
は，S(h) と透水性の間に[ ∫ S(h)∫1[ ∫ S(h)0∫1h(S)−1 dS

0 h(S)

−1

]2

dS

(4)

も示した．Burdine と Mualem モデルそれぞれの物理的な意味合いについては，小杉（2007）に詳しく解説されと表記できる．なお，本稿では表記の簡便さのために負0K(h)Kr (h) == S(h)0.5
Ksている．

θ (h) − θrθs − θr

K(h)Kr (h) == S(h)2
Ks

的な値として p = 0.5 を求め，

h(S)−2 dS]−20 h(S) dS(2)

Hoffmann-Riem et al.（1999）と Kosugi（1999）は，
Burdine モデルと Mualem モデルを一般化した Mualem
and Dagan（1978）に基づき次の一般化透水モデルを提
案した．

K(h)
Kr (h) == S(h) p
Ks[ ∫ S(h)0∫1h(S)−q dS−q
0 h(S) dS]r(5)

ここで， p，q，r はパラメータである．特定の水分保持関数に対してこの一般化透水モデルの積分を解析的にという関係が成立することを土壌間隙径分布のモデルと解いたものは，不飽和透水係数の閉形式解（closed-form

ハーゲン · ポアズイユの法則によって導いた．ここで，solution）とよばれ，Burdine モデル（ p = 2，q = 2，r = 1）
や Mualem モデル（ p = 0.5，q = 1，r = 2）をまとめて

Ks は飽和透水係数，Kr は相対透水係数である．Mualem（1976）は Burdine のモデルをベースとして，異なる間隙径の毛管が連結しているときの等価間隙径（equivalentpore radius）の仮定から，
K(h)= S(h) p
Kr (h) =Ks[ ∫ S(h)0∫1h(S)−1 dS−10 h(S)]2dS

表現できる．本稿では，水分保持曲線−不飽和透水係数
連結モデルの透水モデルを区別する場合は，式（2）はBurdine，式（3）は Mualem，式（5）は General と表記
し，表記のない場合は式（5）の一般化透水モデルとする．
閉形式解による水分保持曲線−不飽和透水係数連結(3)1 Natural Science Laboratory, Toyo University, 5-28-20 Hakusan,

Bunkyo-ku, Tokyo 112-8606, Japan. Corresponding author. 関 勝寿, 東洋大学自然科学研究室

モデルの代表格が Mualem-van Genuchten モデルである（van Genuchten, 1980）．閉形式解の連結モデルは，少ないパラメータで保水性と透水性の両方を表現できるため，特に不飽和水分移動のシミュレーションにおいて有2 Graduate School of Bioresources, Mie University, 1577 Kurimamachiya-

用である．しかし，透水モデルの閉形式解の導出過程は，

cho Tsu, Mie 514-8507, Japan.
2023 年 5 月 8 日受稿 2023 年 6 月 25 日受理

モデルの仮定を確認したり，より利用のしやすいモデルへと改良したりする際に重要であるが，数学的な煩雑さ

土壌の物理性

20

第 154 号もあり，最近では閉形式解の導出を確認する機会は減っ(2023)

これを式（7）に代入して，ていると思われる．

B = A(0) = hb −q

そこで本稿では，van Genuchten モデルを中心としたλ

基本的な水分保持関数 S(h) に対する式（5）の水分保持曲線−不飽和透水係数連結モデルの閉形式解の導出方法(の取り組んだ線型結合モデルを紹介する．

Kr (h) = S(h) p

2． 基本モデル(=

Brooks and Corey（1964）モデル（以下，BC モデル），
van Genuchten モデル（VG モデル），Kosugi（1996）モデル（KO モデル）の 3 つを基本モデルとして, それぞれについて式（5）の閉形式解を導出する．
∫ S(h)

h(S)−q dS

(6)

0

とすると，

S(0) =(=∫ 1

B = A(0) =

h
hb
h
hb

h(S)−q dS(7))−r(λ +q))−pλ −r(λ +q)(13))−(p+r)λ −qr

Kr (S) =0hhb

また，h ≤ hb では Kr (h) = 1 である．これが，General −BC モデルの閉形式解である．式（9）と式（13）から K
を S の関数であらわすと,θ (0) − θr θs − θr==1
θs − θrθs − θrより(12)

したがって，式（8）より h > hb では，を解説する．そして，閉形式解の応用例として，著者ら

A(h) =

)−1
+1

(q

qr
K(S)
= S p+r+ λ
Ks

(14)

Brooks and Corey（1964）では，この式に p = 2，q = 2，
r = 1 を代入した Burdine–BC モデルの Kr = S3+2/λ が導出されている．Mualem–BC モデルでは，q = 1，r = 2
を代入した Kr = S p+2+2/λ となる．
2.2 VG モデルであるので式（5）は，次に，VG モデルの水分保持関数は，]rp A(h)Kr (h) = S(h)
B
[

[

]m
1
S(h) =
1 + (α h)n(8)と書くことができる．

2.1 BC モデル

ここで，α ，m，n はパラメータであり，α > 0，0 < m < 1，n > 0 である．

BC モデルの水分保持関数は，まずは式（15）を変形して h−q の式を得る．S(h) を S

( )−λ h(h > hb )hb
S(h) =
1(h ≤ hb )

と表して，

(9)
1 + (α h)n =ここで，hb ，λ はパラメータであり，hb > 0，λ > 0 である． h > hb において h を S の関数として求めると，

h(S) = hb S− λ

1

(10)

(α h)n =S1/m[(α h)−q =[h−q = α qq
λA(h) = hb −q
S dS
0
(
)−1
q
−q q
= hb
+1
S(h) λ +1
λ()−1 ( )−λ −qqh= hb −q
+1
λ
hb

(11)

1
S1/m

1 − S1/mとなるため，式（6）の A(h) を計算すると，∫ S(h)

(15)

S1/m

]q
n

(16)

1 − S1/m
] nq

1

Sm
1

1−Sm

この式から式（6）の積分を計算する．
ここで，一般に定積分の積分変数を y から x に変換す
るには，y = g(x) として，

解説：一般化透水モデルによる不飽和透水係数の閉形式解∫ g(b)∫ b

f (y)dy =
g(a)

a

dy
f (g(x)) dxdx（ p = 0.5 の式が示されている），m = 1 − 2/n の Burdine–(17)が成立することを使う．S = xm と変数変換をすると式（6）と式（16）から，∫ S(h)

A(h) =
0

= αq

= αq

= αq

∫ S(h) m1

[

α
[

0

∫ S(h) m1

[

0

q

x
1−x

mx

]q
n

VG モデルとして，
[
(
) ]
1 m
Kr (h) = S(h)2 1 − 1 − S(h) m

く使われていることから，Mualem–VG モデルのことをdS

単に VG モデルとよぶことが多い．

VG モデルにおいて，α h ≫ 1 では水分保持関数の式（15）は λ = mn の BC モデルに近似され，また式（21）は，

dS
dx
dx
(18)

]q

Kr (h) ≈ mr α h−(p+r)mn−qr

n

mx

m−1

(24)

の閉形式解を示した．現在では Mualem–VG モデルが広n

1 − S1/mx
1−x∫ S(h) m1

]q

S1/m

21

(25)

dx
と近似される．すなわち．α h ≫ 1 において log h ‒qn +m−10(1 − x)− qn

log K グラフの傾きは −(p + r)mn − qr となって λ = mnとしたときに BC 式と一致する（関，2022）．

dx

2.3 KO モデル

この式はこのままでは解析的に解くことはできないが，間隙径（間隙を毛管とみなしたときの半径）を ra とq/n+m−1 がゼロ，すなわち m = 1−q/n ならば x の指数部がゼロになり容易に積分できることを van Genuchten（1980）は示した．このときに，n > q となる．なお，
m = 1 − q/n の制限は数学的な利便性によるものであって物理的な意味はない（Hoffmann-Riem et al., 1999）．
m = 1 − q/n と定めたことにより，式（18）はして（一般化透水モデルの r と区別するために ra としA(h) = α q∫ S(h) m1
0

1

pore-size distribution function） f (ra ) を
f (ra ) = d θ /dra(26)によって定義する（Mualem, 1976）．すなわち， f (ra )dra
は間隙径が ra 以上 ra + dra 以下の間隙の単位体積あたりm(1 − x)m−1 dx

S(h) m= α q [−(1 − x)m ]0た），水分量変化に基づく間隙径分布関数（volumetricの体積である．KO モデルは，Brutsaert（1966）のように間隙径分布を対数正規分布で表現したモデルである．

(19)

[
]
1 m
= α q 1 − (1 − S(h) m )

すなわち，

[
]
θs − θr[ln(ra /rm )]2√ exp −f (ra ) =2σ 2ra σ 2π(27)

したがって，

B = A(0) = α q

このときに，γ を水と空気の間の表面張力，β を接触(20)

式（8）より General–VG モデルは，

[
(
) ]
1 m r
Kr (h) = S(h) p 1 − 1 − S(h) m(21)角，ρw を水の密度，g を重力加速度として，定数 C を
C = 2γ cos β /ρw g と定めれば，間隙径 ra と間隙水のマト
リックポテンシャル h の間には，毛管上昇のジュレンの式より（宮﨑ら，2005）h=

ここで m = 1 − q/n であることから式（15）の S，式（21）の Kr ともに q の関数である．General–VG の Kr (S) 表
記は

C
ra

(28)

という関係があるため，付録 A で示すように，KO モデ

[
(
) ]
1 m r
Kr (S) = S p 1 − 1 − S m

ルの水分保持関数(22)

van Genuchten（1980）は，m = 1 − 1/n の Mualem–VG
モデルとして

[
(
) ]
1 m 2
Kr (h) = S(h) p 1 − 1 − S(h) m(23)

S(h) = Q

 ln

(

h
hm
σ

) 


(29)

が導かれる．ここで，hm ，σ はパラメータで，hm > 0，

土壌の物理性

22

第 154 号

(2023)

σ > 0，Q(x) は標準正規分布の x から上側の確率（上側確率）z=

Q(x) = 1 − ϕ (x) =

( 2)
1
x
√ exp −
dx
2
2π

∫ ∞
x

(30)

y + qσ√2(38)

とおいて z2 を計算すると，−qσ y −

ただし ϕ (x) は標準正規分布の累積分布関数である（東y2q2 σ 2
=
− z2
2
2

(39)

京大学教養学部統計学教室, 1991）．Q(x) は Microsoft

Excel の norm.dist 関数で Q(x) = 1−norm.dist(x, 0, 1,
TRUE) と計算できる．また，Q(x) を誤差関数
2
erf(x) = √
π

∫ x

( )
exp x2 dx(31)0によってあらわすと，

[
(
)]
1
x
Q(x) =
1 − erf √
2
2

(32)

が得られ，これを式（37）に代入して積分変数を y からz に変換すると，dy/dz =

√

2 より，

hm −q ∫ z(h)A(h) = − √
exp
2π ∞
hm −q
= − √ exp
π

(

(

)
q2 σ 2 2 dy
−z
dz
2
dz

)
(
)
q2 σ 2 ∫ z(h)exp −z2 dy
2
∞

(40)

式（31）と erf(∞) = 1 の関係を用いると，式（6）の積分変数を

(
ln
y=

h
hm

)
(33)

σ

によって S から y に変換すると，S = Q(y) であり．S = 0
において y = ∞ であるので，∫ y(h)

A(h) =∞( 2 2)
q σ
hm −q
exp
[1 − erf(z(h))]22A(h) =

さらに式（32）より，A(h) = hm
h−q

dS
dy
dy

(41)

−q

( 2 2)√q σexp
Q( 2z(h))
2

(42)

(34)
式（38）と式（33）を順に代入すると，式（33）より

h = hm exp(σ y)(35)であることと，式（30）から]( 2 2) [
ln
(h/h
)
q
σ
m
A(h) = hm −q exp
Q
+ qσ
2
σ
したがって，( 2)y
1
dS dQ(y)== − √ exp −
dy
dy
2
2π

(36)

であることから，式（35）と式（36）を式（34）に代入
すると，

B = A(0) = hm−q

( 2 2)
q σ
exp
2

)

1 ∫ y(h) −qy2
A(h) = − √hm exp(−qσ y) exp −
dy
2
2π ∞
(
)
hm −q ∫ y(h)y2=−√exp −qσ y −
dy
2
2π ∞
(37)

(44)

式（8）より General–KO モデルの閉形式解は，

[ [
(

(43)

Kr (h) = S(h)pln (h/hm )
Q
+ qσ
σ

]]r
(45)

また，Kr を S の関数であらわすためには，Q の逆関数

Q−1 が利用できる．式（29）の逆関数は，
(
ln

さらに，
−1

Q (S) =

h
hm

σ

)
(46)

解説：一般化透水モデルによる不飽和透水係数の閉形式解23と変形できることから，

式（45）に代入すると，[ []]r
Kr (S) = S p Q Q−1 (S) + qσ

∫ si (h)(47)

Q−1 は Microsoft Excel の正規分布の累積分布関数の逆 関 数 norm.inv 関 数 を 用 い て 計 算 す る こ と が で き ，
Q−1 (S) = norm.inv(1 − S, 0, 1) で あ る ．な お ，Kosugi（1996）は Mualem モ デ ル と Burdine モ デ ル ，Kosugi（1999）は式（5）の一般化透水モデルに対して間隙径の対数正規分布から式（45），式（47）を導出しているので参照されたい．

Ai (h) =0h(Si )−q dSi∫ 1

Bi = Ai (0) =0h(Si )−q dSi[

3. 線型結合モデル

水分保持関数を足し合わせた線形結合モデルに応用することが可能である．Durner（1994）は，複数の VG モデ
ルの水分保持関数の重み付き線形和(52)としたときに，式（8）より

K(h)Kr (h) == S(h) p
Ks

ここまで示した基本モデルの閉形式解は，それぞれの(51)∑i=1 wi Ai (h)k∑i=1 wi Bi
k

]r
(53)

サブモデルの式（51）と式（52）は，式（6）と式（7）と等しい．そのため，2.2 節の VG モデルと同様に式（49）から Ai (h) の計算ができ，式（53）の閉形式解を得ることができる．また，Priesack and Durner（2006）は，式（53）において r = 1 とすると，[kS(h) = ∑ wi Si (h)

(48)

i=1

p −1

Kr (h) = S(h) ηk]∑ wi Ai (h)

(54)

i=1

によって柔軟な水分保持曲線を表現した．ここで，k は使用する基本モデル（以降サブモデル）の数，i はサブモ
デル番号，wi は重み係数（0 < wi < 1 ∑ wi = 1），Si (h)は i 番目のサブモデルである．サブモデルにおける VGのように，透水性関数がサブモデルの関数 Ai の重み付き線形和となることを指摘した．ここで η = ∑ki=1 wi Bi であり，式（49）では η = ∑ki=1 wi αi である．

モデルの式（15）は，それぞれのパラメータがサブモデ

Seki et al.（2022）は，式（53）を発展させ，式（48）ルで異なるため，それぞれにサブモデル番号 i をつけて，において BC モデル，VG モデル，KO モデルの水分保持[1
Si (h) =1 + (αi h)ni関数の任意の組み合わせの線形モデルに対して，不飽和

]mi
(49)

透水係数の閉形式解を示した．Table 1 にその計算結果をまとめる．それぞれのサブモデルの Ai (h)，Bi は，対応する基本モデルの A(h)，B と等しい．そのため，BCと表記される．ここで，それぞれのサブモデルにおいて

mi = 1 − q/ni である．

モデル，VG モデル，KO モデルのそれぞれのパラメータに対しては，サブモデル番号 i をつけて表示してある．

Priesack and Durner（2006）は，式（48）と式（49）の線形和の水分保持関数に対して，式（5）の一般化透水モ
デルによる不飽和透水係数の閉形式解を導いた．式（6）は積分変数の変換によって，

たとえば，VG モデルと BC モデルの水分保持関数の線形和は，サブモデル番号を下付き文字で示して VG1 BC2モデルと表記する．特に 2 つの同じモデルを重ねる場
合は dual-モデルと称する．たとえば VG1 VG2 モデルは

dual-VG モデルと表記する．
∫ S(h)

A(h) =

h(S)−q dS

4. 一般化透水モデルの有用性

0

∫ h

=

∞

h−q

∫ h

=

∞

h

−q

= ∑i=1 wik=∑

一般化透水モデルは，Burdine モデルと Mualem モデ

dS
dh
dh
(

k
w
i=1 i

ルの積分の式にあらわれる指数を共通の変数で与えるた

dSikw
i=1 i dh

)

∑
∫ h
∞

h−q

∫ si (h)0め，不飽和透水係数の閉形式解を統一して記述できる利

dh

dSi
dh
dh
−q

h(Si ) dSi(50)点がある．一方，Burdine モデルや Mualem モデルにおけるパラメータ値は，それぞれの物理的な考察に基づくものである．そのため，式（5）の p，q，r を変数とする
ことは，本来のパラメータ値の持つ物理的な意味を失う
ことにはなる．しかし，変数が増えることにより透水モ
デルの自由度を上げることでき，より実用的な式となる利点がある．

土壌の物理性

24

第 154 号

(2023)

Table 1 Brooks-Corey（BC），van Genuchten（VG），小杉（KO）モデルのサブ関数 Si (h) と Ai (h) と Bi ．
Sub-retention functions Si (h), and expression for Ai (h) and Bi for the BC, VG and KO models.
サブモデル

Si (h)(

BCh(h > hbi )hbi(h ≤ hbi )[1 + (αi h)ni ]−mimi = 1 − q/ni
[

KO

(

)−λi

1

VG

Ai (h)Q

ln(h/hmi )Bi

Hoffmann-Riem et al.（1999）は，式（23）の Mualem–
VG モデルや式（24）の Burdine–VG モデルのように q，r を固定するよりも，式（21）の General–VG モデルにおいて p，r を自由パラメータとすることで精度良く透水性を表現できるとした．Kosugi（1999）は，General–KO モ
デルにおいて p，q を自由パラメータとする有用性を示している．Seki et al.（2023）は，KO1 BC2 モデルでは pと q，dual-VG モデルでは p と r を自由パラメータとす
れば，多くの土壌の広い圧力範囲における不飽和透水係)−λi −q((h > hbi )

hbi

hbi

(h ≤ hbi )

Bi

[
(
)]mi
Bi 1 − 1 − S(h)1/mi

]

σi

h

Bi

[
Bi Q

ln(h/hmi )

σi

−q

q

λi

(
hmi −q exp



S(h)



S(hb )S′ (h) =



 1


 Kr (h)Kr (hb )Kr′ (h) =



 1

数に対して，どちらのモデルでも精度良くあらわせるこ

+1

αi q

]
+ q σi)−1q 2 σi 2)

2(h > hb )(55)(h ≤ hb )(h > hb )(56)(h ≤ hb )とを示した．なお，実測値からモデルのパラメータを決定するときは，水分保持曲線，不飽和透水係数の順に決修正関数では，h ≤ hb において飽和となり，h > hb にお定することが多い．dual-VG モデルを含む General–VG

ける S′ と Kr′ が h = hb で連続となるように修正前の関数モデルでは，S(h) と K(h) の両者が q の関数であるため，S(h) で q を決めると透水性関数の q を自由パラメータと全体を定数倍して補正する．式（55），式（56）の修正モ

できなくなることに注意が必要である．

et al.（2023）の線型結合モデルの修正モデルを与えることを付録 B に示す．

デルは，Vogel et al.（2001）の Mualem–VG モデル，Seki

5. 修正モデル
6. おわりにq = 1 の Mualem–VG モデルにおいて，n が 1 に近い
とき，飽和近傍で透水性曲線の傾きが極めて大きくなり，

本稿では，水分保持曲線−不飽和透水係数連結モデル

不飽和透水係数が飽和透水係数に対して不自然に小さくの一般化透水モデルに対して，広く用いられている BCなる問題がある．Vogel et al.（2001）は，仮想的な空気

モデル，VG モデル，KO モデルの水分保持関数から不侵入値 hb を与えて，h < hb において一定の飽和透水係数

飽和透水性係数の閉形式解を導出する方法を示した．本を与える修正 VG モデルを提案した．修正モデルでは，稿で示した計算手順にしたがえば，たとえば Excel に

hb に 1 cm 程度の小さな値を与えると，VG モデルの水分保持曲線の形はほぼ変えずに，透水性曲線の h < hb に

よって様々なモデルの水分保持曲線と透水性曲線を計おける急激な変化を抑制することができる．

る SWRC Fit や unsatﬁt といったプログラムにより実測

Seki et al.（2023）は，線型結合モデルにおいても ni が小さいときあるいは σi が大きいときには同様の問題が（Seki et al., 2023）．

算することができる．また，著者が作成して公開していされたデータからパラメータを決定することができる

生じ得るため，この修正モデルを線型結合モデルに適用

現在の不飽和水分移動のシミュレーションにおいてした式を示した．さらに一般的に，任意の水分保持関数

は ，Mualem–VG モデ ルが最 も広 く用 いら れて いる ．

S(h) と透水性関数 Kr (h) に対する修正モデルは，次の修正水分保持関数 S (h) と修正透水性関数 Kr′ (h) で与えら

HYDRUS に代表される汎用水分移動解析プログラム（Šimůnek et al., 2022）が手軽に利用できるようになったれる．

反面，透水モデルの仮定や限界などが検討されることは解説：一般化透水モデルによる不飽和透水係数の閉形式解少なくなったと思われる．本稿で示した基本モデルの閉形式解の導出過程を確認すると同時に，本稿で引用した

x=

文献の議論を再確認することは，不飽和水分移動モデルのさらなる改良において重要と考える．

25

H − ln(hm ) ln (h/hm )
=
σ
σ

(62)

は標準正規分布 N(0,1) に従う（東京大学教養学部統計学
教室, 1991）．

付録 A

さて，式（1）と式（59）から，

KO モデルの間隙径分布の式（27）を導き，さらに式（28）を用いて水分保持関数の式（29) を導く．KO モデ
ルは，間隙径の自然対数 ln(ra ) を R としたときに，R が
平均 ln(rm )，標準偏差 σ の正規分布に従う，すなわちR ∼ N(ln(rm ), σ 2 ) としたモデルである．このときに raの確率分布を対数正規分布（lognormal distribution）と呼び（東京大学教養学部統計学教室, 1991），rm はその中央値で S(rm ) = 0.5 となる．R と ra の確率密度関数（probability density function）をそれぞれ PR (R)，Pr (ra )
とすると，

[

[R − ln rm ]2
1
PR (R) = √ exp −
2σ 2
σ 2π

S(ra ) =∫ ra1
θs − θr0

f (ra )dra

(63)

式（60）を式（63）に代入すると∫ ra

S(ra ) =Pr (ra )dra

0

(64)

この式の右辺は対数正規分布 ra の下側確率（累積分布関数）であるから，対応する正規分布 R ∼ N(ln(rm ), σ 2 ) の下側確率と等しく，∫ R]

S(R) =(57)

0

PR (R)dR(65)

この式の変数 R を H に変換すると，H の確率密度関数
を PH (H) と書けば，Pr (ra ) = PR (R)

dR
dra

∫ H

]
[ln(ra /rm )]2
1√ exp −=2σ 2
ra σ 2π[

S(H) =
(58)

=−
=

H

∫ ra
0

f (ra )dra

∫ H

∞
∫ ∞

ここで，式（26）より θ は積分によって，

θ (ra ) = θr +∞PH (H)

(59)

dR
dH
dH

PH (H)dH(66)

PH (H)dH

さらに変数 H を x に変換すると，x の確率密度関数をPx (x) と書けば，
∫ ∞

∫1
∫1
と計算できる． 0 f (ra )dra = θs − θr かつ 0 Pr (ra )dra =1 であるから，

S(x) =x

Px (x)dx(67)

この式は標準正規分布の上側確率であるから式（30）のf (ra ) = (θs − θr )Pr (ra )(60)

Q(x) の定義と合致し，式（58）を式（60）に代入すると式（27）が得られる．
次に式（28）を用いて式（29) を導く．式（28）の両辺の対数を取ると，H = ln(h) とおいて，

H = ln(C) − R

S(x) = Q(x)(68)

したがって式（62）を代入すれば式（29）が導かれる．こ
こで h の分布，すなわち間隙水圧分布関数（pore capillary(61)

ここで正規分布に従う変数 R の一次関数であらわされる

pressure distribution function）（Kosugi, 1994）は，ra と同様に対数正規分布となり，hm がその中央値で S(hm ) = 0.5となる．

変数 H も正規分布に従い（東京大学教養学部統計学教室，

以上をまとめると，KO モデルは，間隙径の自然対数 R

1991），その平均は ln(C) − ln(rm ) = ln(C/rm ) = ln(hm )が正規分布に従うと定めることにより，間隙水圧の自然（hm は rm に対応する圧力水頭），標準偏差は σ ，すなわ

対数 H は R の一次関数であるため正規分布に従い，そち H ∼ N(ln(hm ), σ 2 ) である（式（57）に式（61）を代の標準化変数 x が標準正規分布に従う，というモデルで入することで確かめられる）．また，H の標準化変数 x

ある．

土壌の物理性

26

第 154 号

(2023)∑i=1 wi Si (h)k

付録 B

S′ (h) =(77)

式（55）と式（56）の修正モデルの一般形に対して

θm − θr k=∑ wi Si (h)θs − θr i=1

Vogel et al.（2001）の Mualem–VG モデルと Seki et al.
（2023）の線型結合モデルの修正モデルを確認する．
まずは，VG モデルの修正モデルでは，h > hb においS(hb )

式（53）を式（56）に代入すれば h > hb において，て，式（1），式（55）から，′

θs − θr
S(h)S(hb )

θ ′ (h) = θr +Kr (h) = S (h)(69)

ここで θm を[′p∑ki=1 wi Ai (h)∑ki=1 wi Ai (hb )

]r
(78)

式（77）と式（78）は Seki et al.（2022）の式と一致する．

θs − θr
S(hb )

θm = θr +(70)すなわち S(hb ) = 1 としたときの仮想的な θs の値として定義すると，式（69），式（70）から，

引用文献

Brooks, R.H. and Corey, A.T. (1964): Hydraulic properties ofporous media. Hydrology Papers 3. Colorado State Uni-θ ′ (h) = θr + (θm − θr )S(h)versity, Fort Collins, CO, USA.

(71)

Brutsaert, W. (1966): Probability laws for pore-size distributions. Soil Science, 101(2): 85–92.

したがって．VG モデルにおいては式（15）から，

θ ′ (h) = θr +

θm − θrm[1 + (α h)n ]Burdine, N.T. (1953): Relative permeability calculations frompore-size distribution data. Petroleum Transactions of the(72)

American Institute of Mining and Metallurgical Engineers,
198: 71–77.
Durner, W. (1994): Hydraulic conductivity estimation for soils

また，

with heterogeneous pore structure. Water Resources Re-

[

]
1 m

search, 30(2): 211–223.

F(h) = 1 − S(h) m(73)

Hoffmann-Riem, H., van Genuchten, M.Th. and Flühler, H.
(1999): A general model of the hydraulic conductivity ofとすると式（21）より，

unsaturated soils. In van Genuchten et al. (ed.) Proceedings of the International Workshop, Characterization andr

Kr (h) = S(h) [1 − F(h)]p(74)

measurement of the hydraulic properties of unsaturated
porous media., Riverside, CA. 22–24 Oct. 1997. Univer-となるから，この式と式（55），式（56）より h > hb において，Kr ′ (h) = S′ (h) p[

1 − F(h)1 − F(hb )sity of California, Riverside.
Kosugi, K. (1994): Three-parameter lognormal distribution
model for soil water retention. Water Resources Research,]r30(4): 891–901.

(75)

Kosugi, K. (1996): Lognormal distribution model for unsaturated soil hydraulic properties. Water Resources Research,

ここで，Kr と F をそれぞれ S の関数とみなすと，

S′ (hb ) = 1 であるため，
Kr ′ (S′ ) = S′p[]
1 − F(S′ ) r
1 − F(1)conductivity for soils with lognormal pore-size distribution. Soil Science Society of America Journal, 63(2): 270–(76)

そして，式（72）と式（76）は Vogel et al.（2001）の式と一致する．
次に線型結合モデルにおいては，式（48），式（70）を式（55）に代入すれば h > hb において，32(9): 2697–2703.
Kosugi, K. (1999): General model for unsaturated hydraulic277.
小杉賢一郎 (2007): 古典を読む Y. Mualem 著 「不飽和多孔質体の透水係数を推定する新たなモデルについて」 なら
びに M.Th. van Genuchten 著 「不飽和土壌の透水係数を推定する閉形式解について」. 土壌の物理性, 106: 47–60.
宮﨑 毅, 長谷川周一, 粕渕辰昭 (2005): 土壌物理学, pp. 19–

20. 朝倉書店, 東京．

解説：一般化透水モデルによる不飽和透水係数の閉形式解27hydraulic properties. Journal of Hydrology and Hydrome-

Mualem, Y. (1976): A new model for predicting the hy-

chanics, 71(1): 22–34.

draulic conductivity of unsaturated porous media. Water
Resources Research, 12(3): 513–522.

Šimůnek, J., Šejna, M., Brunetti, G. and van Genuchten,

Mualem, Y. and Dagan, G. (1978): Hydraulic conductivity

M.Th. (2022): The HYDRUS Software Package for Sim-of soils: Uniﬁed approach to the statistical models. Soil

ulating the One-, Two, and Three-Dimensional Movement

Science Society of America Journal, 42(3):392–395.

of Water, Heat, and Multiple Solutes in Variably Saturated

Priesack, E. and Durner, W. (2006): Closed-form expression

Media, Technical Manual I, Hydrus 1D, Version 5.0, pp.

for the multi-modal unsaturated conductivity function. Va-

334. PC Progress, Prague, Czech Republic.
東京大学教養学部統計学教室編 (1991): 統計学入門, pp.106–

dose Zone Journal, 5(1): 121–124.
関 勝寿 (2022): Burdine 透水モデルの再評価. 東洋大学紀要
自然科学篇, 66: 33–40．

107, p.121, p.128, p.280, 東京大学出版会, 東京．
van Genuchten, M.Th. (1980): A closed-form equation for

Seki, K., Toride, N. and van Genuchten, M.Th. (2022):

predicting the hydraulic conductivity of unsaturated soils.

Closed-form hydraulic conductivity equations for multi-

Soil Science Society of America Journal, 44(5): 892–898.

modal unsaturated soil hydraulic properties. Vadose Zone

Vogel, T., van Genuchten, M.Th. and Cislerova, M. (2001):

Journal, 21(1): e20168.

Effect of the shape of the soil hydraulic functions near sat-

Seki, K., Toride, N. and van Genuchten, M.Th. (2023): Evaluation of a general model for multimodal unsaturated soil

要uration on variably-saturated ﬂow predictions, Advances
in Water Resources, 24(2): 133–144.

旨

水分保持曲線−不飽和透水係数連結モデルの Mualem モデルと Burdine モデルを一般化した一般化透水モデルに対して，Brooks and Corey モデル，van Genuchten モデル，Kosugi モデルの水分保持関数から不飽和透水性係数の閉形式解を導出する方法を解説した．さらに，線型結合モデルに対する不飽和透水係数の閉形式解の導出方法を示した．関連する議論として，一般化透水モデルの有用性を示し，また，飽和近傍で透水性曲線の傾きが極めて大きくなり，不飽和透水係数が飽和透水係数に対して不自然に小さくなることを避けるための修正式の一般式を示した．
キーワード：水分保持関数，透水性関数，閉形式解，一般化透水モデル，線型結合モデル