土壌水分再分布を考慮した土壌浸透モデルの開発
環境システム工学専攻 水文気象研究室 広田洸平
指導教員 陸 旻皎
1. 背景と目的流出解析において, 土壌中の水分移動を表現するには Richards 式を解く必要がある. しかし Richards 式は強い非線形を有するため, 正確な浸透計算を行うためには, 土壌を細かく刻み, 短い時間間隔で繰り返し計算を行わなくてはならない. そのため分布型水文モデルでRichards 式を用いる場合, 流域スケールでこの数値計算を行う必要があり, 計算時間が膨大なものとなる. 実際にモデルに組み込むことは困難である. 本研究では土壌水分の再分布に着目し, これを数値化することで, 土壌浸透計算の簡略化を目指し, 検討を行った.
2. 研究手法
本研究では Richards 式を数値計算で解くことにより土壌の浸透過程を求めている. 以下に今回使用した Richards 式の基本形を示す.
(1)

図 1

ここで,ここで は体積含水率, t[T]は時間, z[L]は鉛(4)直座標で上向きの土壌深さ, K[LT-1]は土壌の不飽和透水係数,[L]はマトリックポテンシャル水頭である.

モデル計算の概念図[LT-1]は飽和透水係数,は飽和体積含水率,

式(1)にある不飽和透水係数 K の値は vanGenuchten(1980)の式に, Mualem(1976)の理論的モデルを用いることで導かれた次式から(5)は有効飽和度,は残留体積含水率である.

また N, M, αは van Genuchten モデルで用いられる各土壌のパラメータである.
数値計算はこれらの式を用いて行われる. モデルは図 1 に示したように, 土壌を刻み幅⊿z求めることができる.

で細分化し, 地表面から降水を与え, 土壌の i(2)番目の層の体積含水率を計算する構造となっている. このとき透水係数はそれぞれの層の中＋(3)間の位置で計算され, 各層間の水分移動を決定している.

3. 簡略化のための検討θ(0)t=∞

本検討では計算の簡略化のため, 土壌水分はθ(0)t00SWt

計算を開始させた. また境界条件は土壌の上部および下部からの流入や流出は考慮しないよ

Depth [m]態を仮定し, 地表から-0.1m を飽和状態として土壌深さ [m]初め地表面付近に水が溜まっているという状

SMD
初期値
安定状態
t[hr]における水分分布

う設定し, 検討を行った.
3.1.

浸透計算

-1.0-1
図 2

透過程の簡略化について検討を行っている. 図
2 に今回検討を行った土壌水分再分布について

0.20.190.180.170.160.150.140.130.120.11
0.10.03

SW [m]ているため, 地表から-0.1m を飽和状態（θ=土壌水分再分布の概念図0.04の概念図を示す. 今検討では土壌の初期状態を地表面付近に水が溜まっている状態を仮定し0.02θs）としている. 図のように初め上部にあった0.01

0500の出入りがない状態で, 経過時間を無限大にと
図 3っていくと実線で表わされている安定状態にる.

1000

1500

2000

1000
time [hr]

1500

2000

SW の時間変化

1.00
1.0
0.10
10-1SW/SMD

い次第に減少していく領域を SW[L]と定義す

5000

透していく. その後, 土壌表面と底面からの水初期状態と安定状態で囲まれた, 時間経過に伴Loam Sand
Sandy Loam
Loam
Silt
Silt Loam
Sandy Clay Loam
Clay Loam
Silt Clay Loam
Sandy Clay0

土壌水分は, 時間経過に伴い下部に向かって浸なる. このときの, 図中の土壌上部に存在する

θs

体積含水率 θVolumetric Soil Moisture

今研究では土壌水分の再分布に着目して, 浸10-20.01

(6)

0.00
10-30.0010-4

図 3 に計算から得られた, 各土壌の SW と時間との関係を示す. 各土壌によって SW の減少速度は異なるが, いずれの土壌でも SW の値が 0に向かって収束していくのがわかる. 図 4 に,図 3 を無次元化し, 対数をとったグラフを示す.
ここで SW を無次元化するため, 土壌内の乾燥度 合 を示 す 指標 で ある土 壌 水分 欠 損量 (SoilMoisture Deficit, 以下 SMD[L])を用いている.
土壌表面から底面までの各層ごとの水分欠損量を足し合わせることで, 土壌内の乾燥度合を求めることができる.
Σ(7)0

図 4204060time ×Ks / G80100

無次元化した場合の SW の変化

また横軸は飽和透水係数と Capillary Length
Scale を 用 い て 無 次 元 化 す る . Capillary
Length Scale は土壌の水文特性や初期状態に依存している値で, G[L]で表わされる.
(8)
各土壌を無次元化, 対数化することで SW は直線的に減少し, このときの減少式を指数近似を用いて以下のように表わすことが可能とな

る.

41(9)

30.82(10)
ここで得られた各土壌の近似係数

0.9aSW = 5.0134 e-0.918αaSW

ここで漸化式を用いて bSW を消去すると,0.71の値は0.600で, 一定の関係性を見ることができる. 図 5 に
図 5の関係を示す. 図 5 では, α1.0 1は減少している.

なる, つまり土壌水分の減少速度が遅くなる傾向にあることを示している. また図 6 に,Sandy Clay Loam

0.000110-40

図 6いる. このグラフから Richards 式による数値

0.8計算と本検討における近似式は整合性が取れ

204060time ×Ks / G80100

モデル計算値と近似式との比較(SW)

1.01

Se

0.60.4

地表面水分量の計算

浸透モデルにおける土壌の地表面の水分量0.2は, 土壌外部の大気と直接接しているため, 土0壌浸透中の土壌上部からの新たな水分の流入や土壌内水分の大気への蒸発量を決定する上で重要な要素となってくる.

0500

1000
time [hr]1500

図 7地表面の Se の時間変化

0500

図 8地表面の Sr の時間変化

2000

1.01

ここでは式(4)を用いて,地表面の有効飽和度の時間経過に伴う減少の図を示す. 各土壌のは図 7にあるように計算開始直後は急激に減少し、時0.80.6Srの検討を行う. まず図 7 に0.4間の経過とともに次第に一定の値に収束して0.2いく. これは土壌の上部および下部からの水の0出入りが無い状態で, 土壌上部にあった水分が浸透し、土壌全体が安定状態になったためと考えられる. また loam sand や sandy loam といった粗い土壌ほど

Loam
Silt

Sand, Sandy Clay Loam, Loam, Silt を用いて3.2.

15-20.01
10割合で増加させた場合の SW の減少（破線）とていることがわかる.

10土壌パラメータαと aSW の関係-30.00110を一定の比較を示す. ここでは代表土壌として Loam

5Loam Sand

Richards 式を数値計算で解いた場合の SW の減少（実線）と, 式(10)を用いて0.6-1100.1SW/SMD

これはαが大きい土壌ほど図 4 の傾きが小さく0.4

αvan Genuchten パラメータである M と式(10)の値が大きくなるにつれ,0.20

van Genuchten パラメータと対応させることから得られた

Loam Sand
Sandy Loam
LoamSiltSilt Loam
Sandy Clay Loam
Clay LoamSilt Clay Loam
Sandy Clayの減少幅が大きく, 逆

1000time [hr]

1500

2000

3

比較的細かい土壌ほど減少幅が小さくなるこ2とが図から読み取れる. 次に土壌全体の傾向を見るため, 式(11)で地表面の有効飽和度を正規化し, それをプロットしたものを図 8 に示す.

Ln[(1-Sr)/Sr]に silt clay loam や sandy clay といった粒径が10-1-2(11)-30は,き最小となったときの値である. この図からほとんどの土壌で,のときに 1,

図 93式(12)を用いて計算した Sr の変化

1.01のaSr = - 0.045 α + 0.8897

0.8のときに 0 となる性質を式(12) で表現した. これをプロ2ln tの値が一様な減少傾向を示しているとみることができる . またを利用し,1が時間を無限大にとってい0.6aSr

ここで0.4ットすると, 図 9 のように を直線に近い形で0.2表わすことが可能となる.

00(12)

51015α

図 10

土壌パラメータαと aSr の関係

22.0

またこれを漸化式を用いて表わすと,
1.5

ここで

,bSr = 0.1919 lnKs + 1.5345

bSr(13)は図 9 における各土壌の値の減

11.00.5少に関する近似係数である. これらの近似係数はそれぞれの van Genuchten の土壌パラメー

0-8タと対応させると一定の関係性を見ることができる. 図 10 では土壌パラメータのαの増加に伴い

図 11が減少していることがわかる. これはαが大きい土壌ほどく, つまり-6-4lnKs-2

土壌パラメータ Ks と bSr の関係

1.01の減少速度が大き

Silt0.8の減少速度が大きくなることを

Loam0.4

Sandy Clay Loamの 減 少 速 度 が大 き く , つ ま り

0.2

Loam Sandの減少速度が大きくなることを示している.

0の増加に伴い増加することを示しており, これは土壌ほどがが大きいSe0.6示している. 図 11 では0

また図 12 に Richards 式を数値計算で解いた場合の Se の減少（実線）と, 式(13), (11)を用いて計算された Se の減少（破線）との比較を示す.
SW と同様に代表土壌として Loam Sand,0

図 12500

1000
time [hr]15002000

モデル計算値と近似式との比較(Se)20mm/hr beta
0.1mm/hr beta
0.2mm/hr beta

1.80.3mm/hr beta
0.4mm/hr beta

1.60.5mm/hr beta
0.6mm/hr betaSandy Clay Loam, Loam, Silt を用いている.

1.41.01

このグラフから Richards 式による数値計算と1.20.81mm/hr beta

Ep=0mm/hr

Richards modelRichardsmodel

β0.610.40 E =1.0mm/hr612p

182418240.2

4. 蒸発を考慮した土壌水分再分布の推定
ここまでの検討で得られた SW と0.8mm/hr beta
0.9mm/hr beta

本検討における近似式は整合性が取れていることがわかる.

0.7mm/hr betaの減少

02 0式を用いて蒸発を考慮した土壌内の水分の推定を行った. まず蒸発量 E[L] は可能蒸発量[L]を用いて以下のように求め, ここで蒸発612time [hr]0mm/hr beta
0.1mm/hr beta

0.2mm/hr beta
1.80.3mm/hrbeta
Sandy Clay Loam におけるRichards式
0.4mm/hr beta
1.60.5mm/hr beta
を用いて計算したβ(破線)と蒸発能を変化させて

図 130.6mm/hr beta

効率 は地表面の水分量から求めた.
(14)

0.7mm/hr beta
モデル計算から求めたβ(実線)の時間変化
1.40.8mm/hr beta

0.02(15)0.9mm/hr beta

1.2

Ep=0mm/hr

1mm/hr beta
Richards modelRichardsmodelが減少した分だけ SMD が増加する.
(16)
このとき

(17)

1SW [m]

これを各時間の SW から引き, 土壌内の水分量

0.010

Ep=1.0mm/hr

61218240612time [hr]1824は以下のように変化していくもの
0とする.
(18)
ここで得られた値をこの時間ステップにおけ図 14

Sandy Clay Loam における Richards 式る真値とし, 次の時間ステップの計算へと進行を用いて計算した SW(破線)と蒸発能を変化させする. 図 13, 図 14 に代表土壌として Sandyてモデル計算から求めた SW(実線)の時間変化

Clay Loam を 用 い た 場 合 の と SW の ,Richards 式を数値計算で解いた場合の値（蒸発5. 結論

を 0[mm/hr]（蒸発無し）から

本研究では Richards 式を基礎方程式とした1[mm/hr]まで 0.1 ずつ変化させた場合の値の土壌浸透モデルを用いて各土壌の浸透プロフ時間変化を示す. 図からが大きいほど もァイルを求め, そこから得た結果より土壌内のSW も急激に減少しているのがわかる. また水分移動における計算の簡略化を目指した. 結Richards 式で求められた値との場合の果として, 土壌水分再分布の指標である懸垂水モデル計算における値とを比較すると大きな帯水分量は van Genuchten の土壌パラメータ誤差は無く, 同じような傾向で減少しているこで表現できることがわかった. 同様に, 大気-土とがわかる. SW が減少するということは, 初壌間の水分移動を決定する要因である地表面め上部あった水分が下部への浸透および大気水分量も van Genuchten の土壌パラメータを中への蒸発が起こっているということであり,用いて表現することができることがわかった.

この値は土壌内の水分移動を表わす１つの指またこれにより, Richards 式を用いて計算され標と見なすことが出来る. つまりこのモデルをた値と同様な傾向を有する, 蒸発による土壌内用いることで, 土壌内の水分移動を擬似的に求水分の減少を考慮した土壌水分再分布を推定めることが可能である.

することが可能なモデルを構築した.

無し）と,

参考文献
(1) Mualem,Y.

1976.

A

newmodelfor(8) Davies,

J.

A.

andAlle,C.

D.,

predicting the hydraulic conductivity of

1973:Equilibrium, potential

andunsaturated porous media. Water Resour.

evaporation

surfaceRes. 12:513-522

southern Ontario. J. Appl. Met.,12:649-657(2) van Genuchten,M.Th. 1980. A close form
equationforpredicting

thehydraulicfrom

croppedactualin(9) ウィリアム・ジュリー+ロバート・ホートン=著、取手伸夫=監修、井上光弘+長裕幸+西村

conductivity of unsaturated soils. Soil Sci.

拓+諸泉利嗣+渡辺晋生=訳：土壌物理学

Am.J. 44:892-898

中の水・熱・ガス・化学物質移動の基礎と応(3) Simunek, J., M. Th. van Genuchten, andM.

Sejna,

package

TheHydrus-1Dforsoftware

simulating

the用、築地書館、2006(10)

Roger E. Smith with Keith R. J.

Smettem, Philip Broadbridge and D. A.

one-dimensional movement of water, heat,

WoolhiserandHydrologic

multiplesolutes

invariably-saturated media. Version 3.0,
HYDRUS Software Series 1, Department

土:

Infiltration
Applications,

TheoryforAmerican

Geophysical Union Washington, DC, 2009(11)

杉田倫明、田中正編集、筑波大学水of Environmental Sciences, University of

文科学研究室 著： 水文科学、共立出版、

California Riverside, Riverside, CA, 270

2009
(12)

pp., 2005.
(4) 金子真郷(2010)：分布型水文モデルへ導入が可能な簡略化した土壌新とモデルの開発. 平
成 22 年長岡技術科学大学修士論文(5) 小杉賢一郎(2007)：古典を読む

Y.Mualem

著「不飽和多孔質体の透水係数を推定する新た な モ デ ル に つ い て 」 な ら び に
M.Th.van.Genuchten 著「不飽和土壌の透水係数を推定する閉形式解について」.
J.Jpn.Soc.Soil Phys. 106:47-60(6) 小杉賢一郎(1999)：森林土壌の雨水貯留能を評価するための新たな指標の検討. J. Jpn.
For .Soc. 81(3):226-235(7) Kond, J., N. Saigusa and T. Sato, 1990:AParameterization of Evaporation from Bare
Soil Surface. J.Appl. Meteor., 97:739-774

近藤純正 編著：水環境の気象学 ‐地表面の水収支・熱収支‐、朝倉書店、1994