凍結面近傍の不凍水量変化に基づく凍土の水分特性曲線と不飽和透水係数の検討
渡辺晋生 1 · 紀藤哲矢 1 · 坂井 勝 2 · 取出伸夫 1

Evaluation of hydraulic properties of a frozen soil based on observed unfrozen water contents
at the freezing front
Kunio WATANABE1 , Tetsuya KITO1 , Masaru SAKAI2 and Nobuo TORIDE1

不凍水の圧力勾配により低温領域へ液状水が移動する

Abstract: We evaluated hydraulic properties of a frozen
soil based on the observed temperature, water content, andpressure head near the freezing front in a one-dimensional
freezing column experiment for an unsaturated silty loamsoil. Amount of unfrozen water near the freezing front
was more than equilibrium amount of water according to
the retention curve at the equivalent pressure head because
freezing was a slow process. We proposed a modified-θmodel to describe nonequilibrium amount of unfrozen water in the retention curve. Coupled heat and water flow
equations were simulated with a modified-θ model as wellas two other existing hydraulic property models: the classical Harlan model, and the modified-K model. The Harlan model overestimated water flow from the unfrozen region to the frozen region. Although the modified-K model
well simulated water flow in a frozen soil, calculated pressure heads and liquid water content showed unrealisticchanges at the freezing front due to the extremely large
decrease in the hydraulic conductivity. The modified-θmodel agreed well with observed liquid water contents
and pressure heads in a frozen soil. The unsaturated hydraulic conductivity of a frozen soil based on the modifiedθ model would be more physically feasible than other two
hydraulic property models.
Key Words : frozen soil, unfrozen water, water retention
curve, unsaturated hydraulic conductivity, nonequilibrium（Hoekstra, 1966; Fukuda et al., 1980）．こうした土中水の再分布は，春先の農地の水分 · 施肥管理（Baker andSpaans, 1997; Derby and Knighton, 2001）や地耐力，凍土下の土壌微生物の活性（Watanabe and Ito, 2007），融雪水による表土流亡（Singh et al., 2009）などに影響を及ぼすため，凍結 · 融解過程における不飽和土中の水分· 熱移動を予測することは重要である．
凍結にともなう不飽和土中の水分の再分布過程は，水の相転移を含むリチャーズ式と熱移動方程式により予測できる（たとえば，渡辺ら, 2007）．これらの移動式を解くためには，凍土の水分移動特性，すなわち，不凍水量と圧力水頭の関係を与える凍土の水分特性曲線と凍土の不飽和透水係数を適切に表現することが重要である．特に，凍結面近傍の水分移動特性は，未凍土から凍土への水分移動量に大きな影響を及ぼす．しかし，凍土の水分移動特性，とりわけ不飽和透水係数の測定は非常に難し
い．既往の不飽和土の凍結実験（たとえば，Fukuda etal., 1980; Jame and Norum, 1980; 溝口ら, 1986）においても，凍土中の不凍水量や不凍水圧力の測定の難しさから，凍結面近傍における液状水量や土中水の圧力変化の連続測定に成功した例はない．凍土の透水係数の測定も試みられているが（Burt and Williams, 1976; Horiguchi

and Miller, 1986; McCauley et al., 2002），測定可能な温度や水分量範囲は限られている．それゆえ，凍結面近傍

1. はじめに
土中には，0 ◦ C に冷やされても毛管力や表面力の効果
により凍らない水が存在する．この凍土中の不凍水は，さらなる温度低下により氷に相変化して減少する．凍結
面近傍においては，凍土中の不凍水の圧力が未凍土中の土中水の圧力より低いため，未凍土から凍土への水分移動が生じる．氷と不凍水が共存する凍土中においても，
1 Graduate School of Bioresources, Mie University, 1577 Kurima-

Machiya, Tsu, Mie 514-8507, Japan. Corresponding author: 渡辺晋生，の土の水分移動特性の検討は十分に行われていない．
未凍土中の水分 · 熱移動モデルを凍結過程に拡張する
ために，Harlan（1973）は，未凍土の水分量–圧力の関係が凍土の不凍水量–圧力の関係と等しいと仮定し，さらに未凍土の不飽和透水係数–液状水量の関係と凍土の不飽和透水係数–不凍水量の関係を同一の関係式で表現した（Harlan モデル：木下, 1982）．不凍水の圧力は，クラウジウスクラペイロン式により凍土の温度から与えられる．このモデルは，未凍土の移動モデルを比較的単純

三重大学大学院生物資源学研究科に修正することで凍結過程の計算を可能にするが，未凍

2 Dept. Plants, Soils, and Climate, Utah State University, 4820 Old Main

土から凍土への水分移動を過大評価することが指摘され

Hill, Logan, UT 84322-4820, USA2010 年 9 月 23 日受稿 2010 年 10 月 21 日受理
土壌の物理性 116 号，9–18 (2010)ている（Harlan, 1973）．そこで，氷量増加にともなう水
みちの変化や氷による水みちの閉塞を考慮し，未凍土の10

土壌の物理性第 116 号

(2010)

透水係数を抵抗係数と氷量の割合に基づく指数関数を用ある．以下，未凍土，凍土に関わらず，液状水量（不凍いて補正し，凍土の透水係数の低下を表現するモデルが（1）式の Klh , KlT は圧力勾配お水）は θl で表す．また，提案されている（K 補正モデル：Jame and Norum, 1980;

よび温度勾配で移動する液状水の透水係数であり，KlT

Lundin, 1990）．しかし，このモデルは凍土の透水係数をは通常 Klh に土中水の表面張力の温度変化と促進係数を

液状水量の等しい未凍土の透水係数より数オーダー小さ乗じて表現される（Noborio et al., 1996; 渡辺ら, 2007）．

く与えることが多く，物理的に考えにくい（Watanabe,（2）式については，Cp は熱容量，ql は水分フラックス，

2008）. また，抵抗係数の決定法も明確に示されていな

Lf は水の凍結潜熱，λ は熱伝導率である．凍土について

い．抵抗係数は，土質や水分量，さらには凍結面と地下水（1）式と（2）式を連立して解くためには，凍土の水分特面の距離や凍結履歴によっても異なり，氷量のみの関数性曲線（θl − h 関係）と不凍水量曲線（θl − T 関係），透で表現することの限界が指摘されている（Newman and

水係数が必要である．また，熱的性質 C, λ については，Wilson, 1997）．不飽和透水係数は液状水量の影響を最

液状水量 θl だけでなく氷量 θi との関係もそれぞれ必要も強く受けるため，アイスレンズや氷板の生じない不飽となる．

和凍土の透水係数は，氷量ではなく不凍水量の関数で表

氷と水が共存する場合，二相共存系の状態方程式であ現する方が物理的には妥当である（Watanabe and Flury,

るクラウジウスクラペイロン式が成立する．凍上が生2008）．このため，凍土の透水係数モデルの検討には，凍じず，氷圧を大気圧と等しいとみなせる場合，凍土中の結面近傍における不凍水量の把握が重要である．

不凍水の圧力水頭 h と温度 T との関係は次式で与えら

近年，氷量や表面力の効果を考慮した TDR の校正モれる．

デルがいくつか提案され（Bittelli et al., 2004; Watanabe

h=

Lf
T
ln
g Tm

(3)

and Wake, 2009; 渡辺, 2010），凍土中の不凍水量測定の精度が向上している（Stähli and Stadler, 1997; Watanabe
and Wake, 2008）．そこで本研究では，不飽和土カラム

ここで，g は重力加速度，Tm は融点である．そして，土の鉛直一次元凍結実験において，温度，テンシオメータ

形状を同一とみなし（渡辺ら, 2007），凍土の θl − h 関係

による土中水圧力，TDR による液状水量の連続的測定が未凍土の水分特性曲線（θ − h 関係）と等しいと仮定を行った．凍結面近傍の凍土の土中水圧力の変化を温度すると，（3）式用いて任意の温度 T の凍土の不凍水量とクラウジウスクラペイロン式により推定し，凍結過程にある土の水分量と圧力の関係を求めた．得られた凍結

θl を推定できる（Williams, 1964; Koopmans and Miller,1966）．温度平衡状態にある凍土へのクラウジウスクラ

過程の水分量–圧力関係に Durner モデル（Durner, 1994）

ペイロン式の適用の妥当性は，カロリメトリー法や核磁を適用し，未凍土と同様に凍土の不飽和透水係数を表

気共鳴法を用いた実験からも示されている（たとえば，

現するモデルを提案した（θ 補正モデル）．そして，従来の Harlan モデルや K 補正モデルと，ここで提案した

θ 補正モデルを用いた水分 · 熱移動モデルを凍結過程の水分量変化の実測値と比較し，凍結面近傍の土の不飽
和透水係数の検討を行った．なお，凍結過程の不飽和土
中の水分 · 熱移動のそれぞれのモデルの数値計算には修正版 HYDRUS-1D（Hansson et al., 2004; 渡辺ら, 2007;

Šimůnek et al., 2008）を用いた．

2. 凍結過程の土中の水分 · 熱移動モデル

中の水の固液界面と未凍土中の水の気液界面の幾何学的Black and Tice, 1989）．
2.2 水分移動特性と熱的性質
水 分 特 性 曲 線 は ，van Genuchten モ デ ル（vanGenuchten, 1980）で与えた 2 種類の土の水分特性曲線の線形和で表す Durner モデル（Durner , 1994）を用いた．

θ − θrθs − θr= w1 [1 + (α1 h)n1 ]−m1 + w2 [1 + (α2 h)n2 ]−m2

Se =

(4)

ここで，Se は有効飽和度，θr は残留体積含水率，θs は

2.1 水分 · 熱移動モデル

飽和体積含水率，α , n, m = 1 − 1/n は水分特性曲線の傾凍結をともなう土中の水分，熱移動は，非等温下の水きを与えるパラメータ，w1 ，w2 は重み付けパラメータ

分，熱移動式に氷の相変化と潜熱項を加えて表されるであり w1 + w2 = 1 である．Mualem（1976）のモデル（渡辺ら, 2007）．水分量が比較的高く，水蒸気移動が無に基づく閉形式の不飽和透水係数は，次式で与えられる

視できるときは，次式となる．

（Priesack and Durner, 2006）．

()∂ θ l ρi ∂ θ i
∂
∂h
∂T
+
=
Klh
+ Klh + KlT∂tρl ∂ t
∂z
∂z
∂z

(1)

[
]
∂T
∂ θi
∂
∂T
∂T
Cp
λ
− Lf ρi=−Cl ql
∂t
∂t
∂z
∂z
∂z

(2)

Klh (θ ) = Ks (w1 Se1 + w2 Se2 )ℓ
{ (
) }
{ (
) }]2
[
1/m m1
1/m m2
w1 α1 1− 1−Se1 1+w2 α2 1− 1−Se2 2×(w1 α1 + w2 α2 )2(5)

ここで，t は時間，T は温度，h は圧力水頭，z は鉛直位

ここで，Ks は飽和透水係数，ℓ は間隙結合係数である．

置，θ は体積含水率，ρ は密度である．添え字の l，i は

なお，w2 = 0 のとき，Durner モデルは単一な土となり，それぞれ液状水，氷を示し，θl + θ（
i = θ ）は全水分量で通常の van Genuchten モデルとなる．

論文：凍結面近傍の不凍水量変化に基づく凍土の水分特性曲線と不飽和透水係数の検討
液状水（不凍水）量 θl の凍土の透水係数 Klh (θl ) は，間隙氷の成長にともなう水みちの変化などにより， θl と11の初期状態に対して，カラム上下端の温度制御セルに恒温冷媒を循環させることで，試料上端に − 8 ◦ C，下端に等しい水分量の未凍土の透水係数 Klh,unfrozen (θl ) より低

2 ◦ C の温度境界条件を与え，試料を上端から下方へ凍結いと考え，未凍土の透水係数に対して氷量 θi の割合によした．実験中，上下端はカラム内への給排水のない水分り指数関数的に減少する Klh (θl ) を仮定するモデルが提フラックスゼロの境界条件とした．試料の凍結過程にお案されている（Jame and Norum, 1980; Lundin, 1990）．

いて，比誘電率と土中水の圧力水頭，温度の経時変化をθ− θi Ω

Klh (θl ) = Klh,unfrozen (θl ) × 10

(6)

5 分間隔で測定した．TDR 測定には TDR100（Campbell
Scientific, Inc.）を用いた．凍土の比誘電率は，液状水量のみならず，氷量や土粒子表面の吸着水量の影響を受け

ここで，抵抗係数 Ω は経験的に与えられる定数である．

る．そこで，TDR で測定した比誘電率を，核磁気共鳴実（6）式に従えば，不凍水量 θl が等しい凍土であっても，験との比較をもとに作成した凍土用の検量式（Watanabe

全水分量 θ が異なれば氷量 θi（= θ − θl ）も異なるため，

and Wake, 2009; 渡辺ら, 2010）を用いて不凍水量 θl に

Klh − θl 関係は一致しない．また，Ω は氷の発生にとも換算した．実験は三連で行い，凍結開始前（0 h）および

なう不凍水の流路の変化を与えるパラメータとみなせる

凍結開始後異なる時間（24，48 h）に試料を 2.5 cm 間隔

が，砂やシルトでは 8 程度と非常に大きな値の報告もあで切り分け，炉乾により各凍結時間の各深さの全水分量り（Gosink et al., 1988），詳細は不明な点が多い．

θ を求めた．なお，各カラムにおいて温度と比誘電率分

凍土の熱容量 Cp は，各相の熱容量とそれぞれの体積

布の時間変化が一致することを確認した．また，同一高割合の積の総和として表される（渡辺ら, 2010）．一方，さにおけるカラム中心部と壁面近傍の温度差は 0.1 ◦ C 以

凍土の熱伝導率 λ は液状水量と氷量によって非線形に変内であり，土が一次元的に凍結しているとみなした．

化する．本研究では，Hannson et al.（2004）のモデルを用いた．

λ = C1 +C2 (θ + F θi )[]− (C1 −C4 ) exp −{C3 (θ + F θi )}4(7)

ここで，C1 , . . . ,C4 は実験，あるいは固相率などの物性値からそれぞれの土に与えられる定数，F1 , F2 , F =1 + F1 θi F2 は凍土の熱伝導率の非線形性を表現するパラメータである．
（7）式は Campbell（1985）の未凍土の熱伝導率モデルに氷量の項を加えたモデルであり，渡辺ら（2010）はこれを実測値に基づき検討している．

3. 凍結実験

3.2 水分量と圧力の関係
Fig. 1（a）に，25 ◦ C の室温において吸引法（HW），
加圧板法（PP），露点計測式水ポテンシャル計（WP）（Campbell et al., 1973; Gee et al., 1992）で測定した試料の水分特性曲線を示す．Fig. 1 には，カラム凍結実験に
おいて TDR とテンシオメータで測定した各深さの液状
水量（0 ◦ C 以下では不凍水量）と圧力（測定可能範囲：
h ≥ −800 cm）の関係，および凍結後の TDR で測定した不凍水量と（3）式を介して温度から換算した圧力の関係（熱電対の測定誤差が圧力推定に影響しない範囲：

T < −0.08 ◦ C）を示した．Fig. 1（a）の拡大図である
Fig. 1（b）の第 2 横軸には，（3）式に基づく圧力に相当する凍土の温度を示した．圧力，温度軸は，それぞれ絶対値の指数軸を用いた．図は 0 ◦ C 以下の温度低下によ

3.1 カラム実験

り不凍水量が減少すること，またそれにともない圧力も

北海道十勝農業研究センター芽室研究拠点で採土した低下することを表している．また，図から（3）式による

火山灰性のシルトローム（採土時の体積含水率 θ = 0.4

凍結前後の土中水の圧力変化は， −0.08 ◦ C 相当の −800

m3 m−3 ）を，内径 7.8 cm，高さ 35 cm のアクリル鉛直
カラムに，採土時の乾燥密度（ρb = 0.95 g cm−3 ）に等し
く下層より 2 cm ずつ順に充填した．カラム側面に 3 線
式 TDR プローブ（線径 1.5 mm × 間隔 5 mm × 長さ 74mm）とテンシオメータ（直径 6.5 mm × 長さ 10 mm）を
5 cm 間隔で，T 型熱電対（線径 0.08 mm）を 1 cm 間隔で挿入し，カラム壁面を厚さ 1 cm の発泡ゴムシートとグ
ラスウールで断熱した．定容積採土円筒（100 mL）を用

cm 付近においてテンシオメータの測定値と連続的な変化を示すことがわかる．テンシオメータと（3）式のクラ

いた変水位透水試験により求めた試料の飽和透水係数はで測定された凍結面近傍の水分量は同じ圧力の室温の水

18 cm d−1 あった．また，非定常法で求めた試料の熱伝導率変化に適合した（7）式のパラメータは，C1 = 0.286，
C2 = 1.26，C3 = 6.61，C4 = 0.172，F1 = 2.5，F2 = 1 で分特性曲線から得られる水分量より高くなった．こうしあった．

cm）で顕著であった．一方，凍結速度の遅いカラム中央部（12.5 cm）や，カラム上部においても，凍結が進行して温度の変化速度が遅くなる T < −1 ◦ C においては，凍土の水分量 θl と圧力 h の関係は，未凍土の水分特性曲線

カラムを 4 ◦ C の恒温室に 24 h 静置したところ，カラム内の温度分布は T = 4 ◦ C で一定，水分分布は上端で

h = −340 cm，下端で −305 cm の直線分布になった．こウジウスクラペイロン式による温度からの推定を併用す
れば，凍結面近傍の圧力変化を連続的に測定可能である．
カラム実験開始時の未凍結試料の水分量と圧力の関係
は，吸引法や加圧板法で求めた水分特性曲線と一致した．
凍結面への水分移動により未凍土内の水分量が低下して圧力が h = −500 cm 程度以下に減少すると，カラム実験た傾向は，T > −1 ◦ C の温度変化の速い凍結初期における凍土において見られ，とりわけカラム上部（2.5，7.5

12土壌の物理性第 116 号

(2010)
Temperature, T (-oC)
0.1
1

]

0.6

Column experimentinitial condition
0.4

0.2

HW
PP
WP
TDR
(a)

101

102

103

104

Volumetric water content, q (m3m-3)

Volumetric water content, q (m3m-3)

0.4
Harlan model
Modified-K model
Modified-q modelFreezing

2.5 cm depth
0.3
7.5 cm17.5 cm
12.5 cm

0.2

Harlan model
Modified-K model

_

0.1
105
103
Pressure head, h (-cm)

(b)
105

104

(a) シルトロームの水分特性曲線．HW, PP, WP はそれぞれ室温における吸引法，加圧板法，露点計測式水ポテンシャル計の測定値．× 印はカラム凍結実験における TDR による測定値．実線と点線は Table 1 のパラメータ値を用いた Durner モデルによる適合曲線．(b) は (a) の破線枠内の拡大図．カラム実験において TDR で測定した各深さの水分量に対応する圧力は，h > −800
cm はテンシオメータによる測定値（白抜プロット），T < 0◦ C は（3）式による推定値（灰塗プロット）．
Fig. 1(a) Water retention curves for the silt loam soil. HW, PP, and WP represent observed data at the room temperature using a hanging
water column, a pressure plate, and the dew point potensiometer, respectively. Multiple marks are liquid water contents measured with
TDR in a column freezing experiment. Solid and dotted lines are the fitted Durner model with parameter values shown in Table 1. (b)Enlargement of part of Fig. 1(a). Outline marks are observed pressure heads using tensiometers while gray marks are estimated values
based on Eq. (3) for T < 0◦ C.
Table 1 数値計算に用いた水分特性曲線（（4）式），不飽和透水係数（（5）式）のパラメータ．
Soil water retention curve (Eq. (4)) and hydraulic conductivity (Eq. (5)) parameters used in numerical simulation.

Harlan モデル
K 補正モデルθ 補正モデル

θr

θs

a1
cm−1

n1

w2
cm−1

a2

n2

ℓ
cm d−1

Ks

Ω

0.019
0.019
0.085

0.625
0.625
0.625

0.01010.01010.0133

1.35
1.35
2.33

0
0
0.531

−
−
0.0003

−
−
1.85

−1
−3
1.42

18
18
40

0
7
0

に近くなった．温度制御開始後 48 h の凍結面にあたる
17.5 cm については，温度が 0 ◦ C から −0.3 ◦ C へと低下した（Fig. 1（b）第 2 横軸）が，温度の変動も激しかった．これは，熱電対の測定精度が ±0.03 ◦ C 程度あるこ

関係ついては，−300 cm 以上の未凍土の水分特性曲線ととや，潜熱の発生による温度の変動が原因である．その1 に示す．Fig. 1 には，適合した van Genuchten モデル
を実線，Durner モデルを点線で示す．2 種類の土の水分特性曲線を重ね合わせる Durner モデルは，凍結過程の
階段上の水分量–圧力関係を良く表現している．

ため，0 ◦ C から −0.3 ◦ C における 12.5 cm と 17.5 cm の温度から換算される圧力の逆転は，推定値の誤差の範囲とみなすのが妥当と考える．温度の低下速度に氷の成長速度が及ばない非平衡状態が生じると，不凍水量と圧力

カラム凍結実験で測定した 2.5, 7.5, 12.5, 17.5 cm の水分量–圧力関係の全ての測定値に対して Durner モデルを適合した．最適化したそれぞれのパラメータの値を Table

3.3 数値計算の関係は，平衡状態の水分特性曲線より高い液状水量を（1）∼（3）式の水 · 熱移動モデルに対して，次の 3 種持つ．これが，凍結面近傍の凍土が温度平衡状態にある類の水分移動特性モデルを適用し，凍結過程の水分分布

凍土より多くの不凍水を持つ原因と考えられる．一方，の数値計算を行い，それぞれのモデルについて比較検討

不凍水と氷が平衡状態に近づくと，不凍水の圧力と水分した．

量の関係に対して未凍土の水分特性曲線が適用できると（1）Harlan モデル：未凍土の実測値に基づく水分特性考えられる．

Fig. 1 の水分特性曲線と水分量–圧力関係に対して，（4）式の Durner モデルをエクセルのソルバーを利用して非線形最小二乗法により適合した．まず，室温で測定した曲線と（5）式の不飽和透水係数．
（2）K 補正モデル：未凍土の実測値に基づく水分特性
曲線と（6）式の抵抗係数 Ω で補正した凍土の不飽和透水係数．

未凍土の水分特性曲線は，圧力低下に対し水分量が単調（3）θ 補正モデル：凍結初期の温度非平衡の不凍水量に減少しているとみなし，w2 = 0 とした van Genuchten
モデルを適合した．一方，カラム凍結実験の水分量–圧力θl 変化を考慮した非平衡水分特性曲線と（5）式の不飽和透水係数．

論文：凍結面近傍の不凍水量変化に基づく凍土の水分特性曲線と不飽和透水係数の検討

Hydraulic conducitivity, K (cm d-1)0

価した．これは，未凍土から凍結面付近の 0.3 < θ < 0.4

Modified-K model(Unfrozen)

10-2Modified-K
㧔Frozen㧕
7.5 cm12.5 cm-depth

10

-4

10

-6

10

Modified-θ model-8

Harlan model

10

0.1

13

0.6
0.2
0.3
0.4
0.5
Volumetric water content, θ (m3m-3)

3 種類の水分移動特性モデルの不飽和透水係数．K 補
正モデルは凍結時の水分量により，K(θ ) 曲線の経路が異なる（図中 θ = 0.37 周辺）．
Fig. 2

Unsaturated hydraulic conductivities for three hydraulic property
models.

ここで，Harlan モデルと K 補正モデルの未凍土の水分

m3 m−3 の領域の Klh を ℓ のみで再現できないためである．そこで，θs = 0.63 m3 m−3 の飽和透水係数 Ks に仮想的な値を与えることにより，0.3 < θ < 0.4 m3 m−3 の領
域の Klh を変化させた．そして，凍土と未凍土内の両水分分布に最も適合する Ks と ℓ を求めた．Fig. 2 に，それぞれのモデルについて，不飽和透水係数 Klh と水分量
θ （凍土については不凍水量 θl ）の関係を示す．図中，カラム実験の初期水分量（θ = 0.4 m3 m−3 ）を破線で示した．K 補正モデルにおける（6）式の氷量には，数値計算で得られた各深さの氷量を与えた．θ 補正モデルでは Klh − θ 曲線の傾きが未凍土から凍土（θ = 0.4 → 0.3
m3 m−3 ）にかけてなだらかに変化するのに対し，K 補正
モデルでは未凍土と凍土の間で Klh − θ 曲線の傾きが異なり，凍土の Klh の減少割合がいずれのモデルより大きくなった．また，K 補正モデルでは，凍結初期の氷量の違いにより， Klh − θ 曲線が深さ毎で異なった．たとえば，Fig. 2 においては，比較的高い水分量で凍結する 7.5
cm 深に対し，凍結面への水移動により水分量が減少した後に凍結する 12.5 cm 深の透水係数は θ = 0.365 において 1 オーダー小さくなった．

特性曲線は Fig. 1 中の実線の van Genuchten モデル，

θ 補正モデルの非平衡水分特性曲線は点線の Durner モ

4. 結果と考察

デルである．以下の凍結過程の計算においては，初期圧力以下の h < −300 cm の液状水を対象としている．

4.1 温度分布と水分分布

そのため，−300 < h < −80 cm にみられる実線の van
Genuchten モデルと点線の Durner モデルの違いは，計Fig. 3（a）に，試料下端の温度制御開始後 0, 24, 48 hの温度分布を示す．図中，プロットはカラム実験の測定

算結果に影響を及ぼさない．

値，実線は計算値である．カラム実験では，端面の温度

水分移動式，熱移動式における初期条件と境界条件は，制御を開始すると上端より凍結が進行し，24 h で凍結面

カラム実験の実測値を与えた．また，いずれのモデルに（−0.3 < T < 0 ◦ C）が約 11 cm，48 h で約 17 cm まで達

おいても，凍土の土中水圧力は，温度とクラウジウスク

し，凍結面の進行速度は時間の経過とともに遅くなった．

ラペイロンの（3）式により求め，仮定した水分特性曲線

計算結果は，いずれの水分移動特性モデルを用いた場合

により凍土の不凍水量 θl を与えた．

も実測の温度分布と良く一致し，水分特性曲線や凍土の不飽和透水係数 Klh のパラメータの間隙結合係数 ℓ，飽

透水係数モデルの影響は小さかった．

和透水係数 Ks ，抵抗係数 Ω は，水分分布の実測値を出

Fig. 3（b）–（d）に 3 種類の水分移動特性モデルで計

来る限り再現するように試行錯誤で決定した（Table 1）．

算した凍結過程の土中の水分分布を示す．図中，実線とその際，未凍土内の水分分布は，初期水分量の θ = 0.4
m3 m−3 付近の Klh の影響を強く受けること，また，凍土から未凍土への水分移動には，凍結面近傍の 0.3 < θ <0.35 m3 m−3 の Klh の影響が大きいことの 2 点に特に留点線は全水分量 θ と不凍水量 θl を示し，これらの差は
氷量 θi に相当する．白抜プロットは炉乾による全水分量，黒塗プロットは TDR による不凍水量の測定値をそれぞれ示す．カラム実験において，土を上端から一次元

意しながら，それぞれのパラメータ値を試行錯誤により凍結すると，未凍土から凍結面へ上向きに水分が移動し変化させた．また，θ 補正モデルでは，（5）式の間隙結

た．この結果，全水分量は凍土（上方）で増加し，未凍

合係数 ℓ の増大は Klh − θ 曲線の傾きを増加させ，飽和

土（下方）で減少した．また，凍結面は時間の経過にと透水係数 Ks の変化は Klh − θ 曲線の値を平行移動させるもない下方へ進行した．この際，いずれの時間において点を考慮した．

も凍土内（Fig. 3（a）で T < −0.3 ◦ C の領域）の全水分

Harlan モデルは，水分分布の全領域に適合すること

量はほぼ等しく，初期水分量からの増加量はおよそ 0.05ができないため，未凍土内の水分分布と適合するように間隙結合係数 ℓ を求めた．また，Ks には実測値を用い

m3 m−3 であった．また，未凍土内（Fig. 3（a）で T > 0◦ C の領域）の全水分量は全層でほぼ等しく，減少速度もた．K 補正モデルは，凍土内の水分分布と適合するようほぼ一様であった．

に抵抗係数 Ω を決定しつつ，未凍土内の水分分布にも間隙結合係数を ℓ = −1 とした Harlan モデルは，未凍適合するように ℓ を決定した．θ 補正モデルは，ℓ を適

土内の高温端（下端）の水分変化をよく表すが，凍結面合パラメータとするだけでは全領域の透水係数を過小評

近傍の未凍土の水分量を過小に，そして凍土の全水分量

14

土壌の物理性第 116 号

0(2010)

00h

10Depth, z (cm)Depth, z (cm)

6h

20
0h
24h
48h

(a)

30
-8

-6

10

24h
48h

20

(b)

30

-4
-2
0
Temperature, T (oC)

0

2

4

0

Harlan model(Ͱ = -1, Ω = 0)0.2
0.4
Volumetric water content, θ (m3m-3)

0

0h

0h

10

6h

24h

Depth, z (cm)Depth, z (cm)

6h

48h
20

30

(c)
0

Modified-K model( Ͱ = -3, Ω = 7)

10

24h
48h

20

30

0.2
0.4
Volumetric water content, θ (m3m-3)(d)0Modified-θ model(Ͱ = 1.42, Ω = 0)0.2
0.4
Volumetric water content, θ (m3m-3)

Fig. 3 凍結過程の (a)Harlan モデルによる温度分布．プロットは熱電対による測定値．(b) Harlan モデル，(c) K 補正モデル，(d)θ 補正モデルによる水分分布．白抜プロットは炉乾による全水分量の測定値．黒塗は TDR による不凍水量の測定値．実線と点線
は全水分量と不凍水量の計算値．

(a) Temperature profiles at 0, 24 and 48 h after freezing started. Plots are observed data in a freezing soil and solid lines are calculatedwith the Harlan model. (b)(c)(d) Soil water profiles. Solid and dotted lines are total and unfrozen water contents calculated with (b) the
Harlan model, (c) the K-modified model, and (d) the θ -modified model.

0を過大評価した（Fig. 3（b））．これは，凍結面近傍の透Harlan model

水係数を過大評価しているため，未凍土から凍土へ過剰Depth, z (cm)な水分移動が生じたことが原因である．そこで，より大10

きな ℓ を与えることで凍結面近傍の透水係数を低下させたところ，凍土内の全水分量分布はわずかに減少したが，同時に未凍土全体の透水係数も低下するため，未凍土内

20Freezing
frontの水分移動量を過小評価し，高温側の下端周辺の水分分
布を過大評価，凍結面近傍の未凍土の水分量を過小評価

Modified-K model

した．このことから未凍土の水分特性曲線を用い，凍結による不飽和透水係数の低下を考慮しない Harlan モデ

Modified-θ model30

ルでは，凍結過程の水分移動を十分に表現できないと言
える．

-10

-8

-6

-4

-2

10
10
10
10
10
Hydraulic conducitivity, K (cm d-1)

0

10

Fig. 4 48 h 凍結後における 3 種類のモデルの透水係数分布．
Hydraulic conductivity profiles at 48h after freezing for three hydraulic property models.

K 補正モデルと θ 補正モデルについては，θ ≈ 0.4
m3 m−3 の未凍土に対して ℓ = −1 の Harlan モデルとほぼ等しい不飽和透水係数 Klh を与えた上で，凍結面近傍の 0.3 < θ < 0.35 m3 m−3 の凍土の Klh を小さくする分布（Fig. 2）を与えたところ，各時間の全水分量分布の実

論文：凍結面近傍の不凍水量変化に基づく凍土の水分特性曲線と不飽和透水係数の検討

15

22.5
27.5
32.5

0.4

900
7.5

7.5

2.5

0.2

0.1
(a) Modified-K model ( Ͱ = -3, Ω = 7)
22.5
27.5
32.5

0.4

0.3
2.5

0.2

700

17.5 cm
depth
12.5

17.5 cm
depth
12.5

7.5

27.5

22.5
32.5 cm
depth

Pressure head, h (-cm)

Liquid water content, θl (m3m-3)

0.3

12.5 17.5500(b) Modified-K

900
7.5

12.5 17.5

700

27.5

22.5
32.5 cm
depth

500
0.1
Ͱ = 1.42, Ω = 0)(c) Modified-θ model
0(d) Modified-θ
300

0

10

20

30

40

0

10

20

30

40

Elapsed time, t (h)Fig. 5 K 補正モデルと θ 補正モデルによる凍結土中の異なる深さの（a）（c）液状水量と（b）（d）圧力水頭の経時変化．プロッ
トは測定値．↑ と ↓ は 7.5 cm 深と 12.5 cm 深においてそれぞれ −0.3 ◦ C と −1 ◦ C に達した時間．
(a)(c) liquid water contents, and (b)(d) pressure heads as a function of time at different depths in a freezing soil simulated with themodified-K model and the modified-θ model. Plots are observed date in a freezing soil. Arrows indicate times reached −0.3 ◦ C and−1 ◦ C at 7.5 cm and 12.5 cm depths.

測値をよく表した（Fig. 3（c），（d））
．しかし，凍結面近4.2 液状水量と圧力水頭の時間変化

傍の未凍土内における水分減少については，K 補正モデ

未凍土から凍土への水分移動量は，凍結面近傍の圧力

ルの方が少なくなった．また，K 補正モデルの液状水量

変化や液状水量変化にともなう透水係数の変化に最も強は凍結面で急変するのに対し，θ 補正モデルでは凍結面く影響を受ける．そこで，全水分量分布をよく表した K

で未凍土の液状水量と凍土の不凍水量が滑らかに変化し補正モデルおよび θ 補正モデルについて，さらに詳細なた．Fig. 4 に 3 つのモデルで計算した凍結開始後 48 h の検討を行った．Fig. 5（a）–（d）は，異なる深さにおけ

不飽和透水係数の分布を示す．K 補正モデルの透水係数

る土中の液状水量と土中水圧力水頭の経時変化である．

は，凍結面において 10−2 cm d−1 から 10−8 cm d−1 近くが原因である．Fig. 3（c），（d）の −1 ◦ C 以下の凍土内

Fig. 5 中，プロットは TDR による液状水量とテンシオ
メータによる圧力の測定値，実線は計算値である．Fig.
6 は，Fig. 5 に示した凍結面近傍の各深さにおける K 補
正モデルと θ 補正モデルの不飽和透水係数の変化である．凍結面（0 < T < −0.3 ◦ C）は，凍結開始より約 10 h
で 7.5 cm 深，25 h で 12.5 cm 深に達しており（Fig. 5 の上向き矢印），22.5 cm 以深の測点は常に未凍土内に位置の不凍水量分布については，いずれの水分移動特性モデした．液状水量の実測値は，未凍土内では時間とともにルを用いた計算も実測値をよく表した．温度の変化速度

緩やかに減少し，0 ◦ C 近くの凍土内では温度の低下にとの遅い −1 ◦ C 以下の凍土については，水分特性曲線の差もない急激に減少した．また，温度変化が緩やかになっは少なく，未凍土の水分特性曲線と（3）式により温度かた −1 ◦ C（Fig. 5 の下向き矢印）以下の凍土内では液状ら推定した液状水量が，不凍水量を妥当に表していると水量（不凍水量）はほぼ一定となった．圧力水頭は初期言える（Fig. 1（b）の第 2 横軸参照）．

圧力から水分量の減少にともない低下した．K 補正モデまで急激に低下したのに対し，θ 補正モデルの透水係数
は凍結面近傍の未凍土内から凍土内にかけてなだらかに低下し，その低下量も K 補正モデルに比べ少なかった．

Fig. 3（c），（d）に見られる凍結面近傍の水分分布の違い
は，こうした透水係数分布に対応した水分移動量の違い

16

土壌の物理性第 116 号が，Durner モデルにより未凍土の持つ団粒構造を考慮し-1

Hydraulic conducitivity, K (cm d-1)

10た水分特性曲線による検討も必要であろう（坂井 · 取出,2007）．さらに，（1）∼（3）式の温度と圧力が平衡状態

-3

10

-5

7.5 cm
depth

2.5

10

12.5

として与えられる水分移動モデルを，凍結速度を考慮した非平衡モデルへと改良することも今後の課題である．

Modified-θ
model

Modified-K
model

-7

17.5

5. おわりに
凍結過程にある土の水分特性曲線と不飽和透水係数の10検討を目的に，不飽和シルトロームの一次元カラム凍結実験を行い，水分 · 熱移動モデルによる解析を行った．

-910凍結面近傍の圧力水頭は，テンシオメータによる測定と凍土内のクラウジウスクラペイロン式に基づく温度から

-11

10

0

10

20

30

40

Elapsed time, t (h)Fig. 6(2010)

K 補正モデルと θ 補正モデルによる凍結土中の異なる深さの不飽和透水係数の時間変化．

Unsaturated hydraulic conductivities as a function of time at different depths in a freezing soil simulated with the modified-K
model and the modified-θ models.

の推定により，連続的に測定できることを確認した．また，凍結速度が遅い際には，凍土の不凍水の圧力と不凍水量の関係が未凍土の水分特性曲線と一致したが，凍結速度が速い凍結面近傍においては，凍土は水分特性曲線から推定される水分量より過剰の不凍水を持つことが明らかになった．これは，温度の低下に対して氷の成長が瞬時には生じないことによる非平衡状態が原因と考えられた．

ルでは，実測値に見られない液状水量と圧力水頭の極端な変化が凍結面を境に現れた（Fig. 5（a），（b））．また，得られた水分量や圧力の実測値に対して，未凍土の水分特性曲線と不飽和透水係数を凍土にも拡張して用いるの水分や圧力が減少せず，0 ∼ −1 ◦ C の凍結時に不凍水

Harlan モデル，凍土の不飽和透水係数を氷量と抑制係数で指数関数的に減少させる K 補正モデル，凍結初期の量が急激に減少した．これは，氷の発生により不飽和透不凍水量変化を考慮した非平衡の水分特性曲線と不飽和

水係数 Klh (θl ) が，狭い領域で（Fig. 4）短時間に（Fig.

透水係数を凍土と未凍土に用いる θ 補正モデルを用い

6）急減することに起因する非現実的な変化である．また，こうした Klh の急激な変化は，凍結面近傍の数値計て，凍結過程の透水係数の変化を検討した．Harlan モデ

計算結果は，実測値に比べ 17.5 や 22.5 cm 深の未凍結時

ルは，未凍土から凍土への水分移動を過大評価したが，収束に多大な時間を要したり，計算を不能としたりしや

K 補正モデルと θ 補正モデルは凍結過程にある土中の水分分布を概ね表すことができた．この際，K 補正モデル

すい．

については実測値に見られない液状水量と圧力水頭の極算に極端に小さな時間刻みと空間刻みを必要とし，解の一方，θ 補正モデル（Fig. 5（c），（d））では，Fig. 5

端な変化が凍結面を境に現れた．これは K 補正モデル（a），（b）の K 補正モデルに見られた凍結時の液状水量の透水係数が凍結面を境に急激に変化することに起因すと圧力水頭の凍結開始時の急変は現れず，実測値の滑らる．一方，θ 補正モデルは，凍結時の液状水量と圧力水

かな変化をよく再現した．透水係数の変化も，液状水量

頭の滑らかな変化をよく再現し，不飽和透水係数の変化の変化同様，滑らかに変化しており（Fig. 6），K 補正

も滑らかな減少を示した．比較的速い速度で凍結する土モデルに比べてより現実的な変化を示していると思われでは，不凍水量と圧力の関係を表す非平衡状態の水分特る．また，数値的な評価が困難になることは少なく，計性曲線が評価できれば，未凍土と凍土を同一の水分移動

算時間も K 補正モデルに比べて短い利点を持つ．

特性関数を用いることにより，より物理的に妥当な不飽

本研究は，比較的速い速度で凍結する土において，不和透水係数の評価が可能であると考えられる．その際，凍水量と圧力の関係を反映した非平衡状態の水分特性曲

非平衡状態の水分移動特性関数には，自由度の高い曲線線（θ 補正モデル）が評価できれば，より物理的に妥当を与えることができる Durner モデルが有用である．今な不飽和透水係数の評価が可能なことを示していると考後は，凍結過程に加えて融解過程に対しても，異なる土

える．凍土の不凍水量と圧力の関係に対して，土の深さ

質，初期水分量，凍結速度における水分特性曲線や凍結や水分量，凍結速度，すなわち温度の変化速度が及ぼす面近傍の透水係数の変化を検討することが課題である．

効果についてはさらなる検証が必要である．また，本研

謝辞

究では，凍結過程の水分分布の実測値に適合して Durner
モデルの Ks と ℓ を決定したが，たとえば蒸発法やマル

試料の採土にあたって北海道農業研究センター（芽室チステップ法などにより求められる未凍土の不飽和透研究拠点）の岩田幸良氏のご助力を頂いた．本研究の一水係数との比較も必要であると考える．また，未凍土の部は，科学研究費補助金（20589001）による．ここに記

水分特性曲線に対しては van Genuchten モデルを用いた

して謝意を表します．

論文：凍結面近傍の不凍水量変化に基づく凍土の水分特性曲線と不飽和透水係数の検討

17

Koopmans, R.W.R. and Miller, R.D. (1966): Soil freezing and

引用文献

soil water characteristic curves. Soil Sci. Soc. Am. Proc., 30:
680–685.

Baker, J.M. and Spaans, E.J.A. (1997): Mechanics of meltwater

McCauley, C.A., White, D.M., Lilly, M.R. and Nyman, D.M.

movement above and within frozen soi. Symp. Physics, Chem-(2002): A comparison of hydraulic conductivities, permeabil-

istry, and Ecology of Seasonally Frozen Soils, Special Reportities and infiltration rates in frozen and unfrozen soils. Cold

97–10, pp. 31–36, CRREL, Hanover.
Bittelli, M., Flury, M. and Roth, K. (2004): Use of dielectric spectroscopy to estimate ice content in frozen porous media. Water
Resour. Res., 40: W04212.
Black, P.B. and Tice, A.R. (1989): Comparison of soil freezing
and soil water curve data for Windsor sandy loam. Water Resour. Res., 25: 2205–2210.
Burt, T.P. and Williams, P.J. (1976): Hydraulic conductivity infrozen soils. Earth Surf. Processes, 1: 349–360.

Reg. Sci. Technol., 34: 117–125.
溝口 勝, 中野政詩, 白井清恒 (1986): 不飽和土の凍結における水分 · 溶質 · 温度分布の変化. 農土論集, 122: 11–17.

Mualem, Y. (1976): A new model for predicting the hydraulicconductivity of unsaturated porous media. Water Resour. Res.,
12: 513–522.
Newman, G.P. and Wilson, G.W. (1997): Heat and mass transfer
in unsaturated soils during freezing. Can. Geotech. J., 34: 63–
70.

Campbell, E.C., Campbell, G.S. and Barlow, W.K. (1973): A dew-

Noborio, K., McInners, K.J. and Heilman, J.L. (1996): Two-point hygrometer for water potential measurement. Agric. Me-

dimensional model for water, heat and solute transport inteor., 12: 113–121.

furrow-irrigated soil: I. Theory. Soil Sci. Soc. Am. J., 60:

Campbell, G.S. (1985): Soil physics with BASIC, p. 150, Elsevier,
New York.
Derby, N.E. and Knightton, R.E. (2001): Field-scale preferential
transport of water and chloride tracer by depression-focusedrecharge. J. Environ. Qual., 30: 194–199.
Durner, W. (1994): Hydraulic conductivity estimation for soils
with heterogeneous pore structure. Water Resour. Res., 30:
211–223.

1001–1009.
Priesack, E. and Durner, W. (2006): Closed-form expression
for the multi-modal unsaturated conductivity function. Vadosezone J., 5: 121–124.
坂井 勝, 取出伸夫 (2007): 砂質土と団粒土に対する水分移動特性の検討, 土壌の物理性. 107: 63–77.

Šimůnek, J., Šejna, M., Saito, H., Sakai, M., and van Genuchten,
M.Th. (2008): The HYDRUS-1D software package for sim-

Fukuda, M., Orhun, A. and Luthin, J.N. (1980): Experimental

ulating the movement of water, heat, and multiple solutes instudies of coupled heat and moisture transfer in soils during

variably saturated media, Version 4.0, HYDRUS Software Se-

freezing. Cold Reg. Sci. Technol., 3: 223–232.

ries 3, Dep. of Environmental Sciences, Univ. of California

Gee, G.W., Campbell, M.D., Campbell, G.S. and Campbell, J.H.

Riverside, Riverside.

(1992): Rapid measurement of low soil water potentials using

Singh, P., Wu, J.Q, McCool, D.K., Dun, S., Lin, C-H. and Morse,a water activity meter. Soil Sci. Soc. Am. J., 56(4): 1068–

J.R. (2009): Winter hydrologic and erosion processes in the

1070.

U.S. Palouse region: Field experimentation and WEPP Simu-

Gosink, J.P., Kawasaki, K. Osterkamp, T.E. and Holty, J. (1988):

lation. Vadose Zone J., 8: 426–436.

Heat and moisture transport during annual freezing and thaw-Stähli, M and Stadler, D. (1997): Measurement of water and so-ing. Proc. 5th Int. Conf. Permafrost, pp. 355–360, Natl. Acad.

lute dynamics in freezing soil columns with time domain re-

Press, Washington, D.C.

flectometry. J. Hydrol., 195: 352–269.

Hansson, K., Šimůnek, J., Mizoguchi, M., Lundin, L.C. and vanvan Genuchten, M.Th. (1980): A closed-form equation for pre-

Genuchten, M.T. (2004): Water flow and heat transport in

dicting the hydraulic conductivity of unsaturated soils. Soilfrozen soil: Numerical solution and freeze-thaw applications.

Sci. Soc. Am. J., 44: 892–898.

Vadose Zone J., 3: 693–704.
Harlan, R.L. (1973): Analysis of coupled heat-fluid transport inpartially frozen soil. Water Resour. Res., 9: 1314–1323.

Watanabe, K. (2008): Water and heat flow in a directionally frozen
silty soil. Proc. 3rd HYDRUS workshop, pp. 15–22, Tokyo
University of Agriculture & Technology, Tokyo.

Hoekstra, P. (1966): Moisture movement in soils under tempera-

Watanabe, K and Flury, M. (2008): Capillary bundle model ofture gradients with the cold-side temperature below freezing.

hydraulic conductivity for frozen soil. Water Resour. Res., 44:

Water Resour. Res., 2: 241–250.

W12402.

Horiguchi, K. and Miller, R.D. (1983): Hydraulic conductivity

Watanabe, K. and Ito, M. (2007): In situ observation of the dis-of frozen earth materials. Proc. 4th Int. Conf. Permafrost, pp.

tribution and activity of microorganisms in frozen soil. Cold

504–509, Natl. Acad. Press, Washington, D.C.

Reg. Sci. Tech., 54: 1–6.

Jame, Y.W. and Norum, D.I. (1980): Heat and mass transfer in渡辺晋生, 大森陽介, 和気朋己, 坂井 勝 (2010): サーモ TDR に

freezing unsaturated porous media. Water Resour. Res., 16:

よる凍土の不凍水分量 · 熱伝導率の同時測定. 雪氷, 72(3):811–819.
木下誠一 (1982): 凍土の物理学, pp. 51–57, 森北出版, 東京.

157–168.
渡辺晋生, 取出伸夫, 坂井 勝, Šimůnek, J. (2007): 凍結をとも18

土壌の物理性

なう土中の水分 · 熱 · 溶質移動モデル. 土壌の物理性, 106:

21–32.

第 116 号

(2010)

Watanabe, K. and Wake, T. (2009): Measurement of unfrozen
water content and relative permittivity of frozen unsaturated

Watanabe, K. and Wake, T. (2008): Hydraulic conductivity insoil using NMR and TDR. Cold Reg. Sci. Tech., 59: 34–41.

frozen unsaturated soil. Proc. 9th Int. Conf. Permafrost, pp.

Williams, P. (1964): Unfrozen water content of frozen soils and

1927–1932, University of Alaska Fairbanks, Fairbanks.

soil moisture suction. Geotechnique, 14: 231–246.

要

旨

不飽和シルトロームの一次元カラム凍結実験において土中の温度，水分量，圧力水頭の変化を測定し，凍土の水分特性曲線と不飽和透水係数の検討を行った．温度の低下速度が氷の成長速度を上回る凍結面
近傍においては，凍土が未凍土の水分特性曲線より過剰の不凍水を持つ非平衡状態が生じた．この非平衡の水分特性曲線を用いる θ 補正モデルを提案し，水分 · 熱移動モデルによる凍結過程の水分移動を数値計算した．そして，未凍土の水分移動特性を凍土に拡張した Harlan モデル，凍土の不飽和透水係数
を氷量と抑制係数で減少させる K 補正モデルと比較検討をした．Harlan モデルは，未凍土から凍土への水分移動を過大評価した．K 補正モデルは，凍結過程の水分分布を概ね再現したが，凍結面における透水係数の急激な減少のため，実測値とは異なる液状水量と圧力水頭の極端な変化が凍結面に生じた．
一方，θ 補正モデルは，凍結時の液状水量と圧力水頭をよく再現し，より物理的に妥当な凍土の不飽和透水係数を評価することができた．
キーワード：凍土，不凍水，水分特性曲線，不飽和透水係数，非平衡