環境システム工学専攻要旨集、平成２７年度

無次元化されたRichards 式を用いた
土壌中及び表層における水分再分布の検討
THE SOIL MOISTURE REDISTRIBUTION IN THE SOIL AND SOIL SURFACE USING
THE NUMERICAL SOLUTION OF DIMENSIONLESS RICHARDS EQUATION

石崎

諒1・陸旻皎2,3

1

Ryo ISHIZAKI and Minjiao LU2,3
1工学士

長岡技術科学大学大学院工学研究科修士課程（〒940-2188 新潟県長岡市上富岡町 1603-1）
2工博 長岡技術科学大学教授（〒940-2188 新潟県長岡市上富岡町 1603-1）
3工博

重慶交通大学兼任教授（〒400074 中国重慶市南岸区学府大道66）
1B. Eng. Student, Nagaoka University of Technology, Niigata, Japan

2Dr. Eng. Professor, Nagaoka University of Technology, Niigata, Japan
3Dr. Eng. Adjunct Professor, Chongqing Jiaotong University, Chongqing, China

流出解析において土壌中の水分移動を考慮することは重要である.土壌浸透は再分布によって形成される直前の土壌水分プロファイルによって大きく変化し,土壌の浸透能によって地下水涵養量に影響するため,高精度の流出解析や洪水予測において考えなければならない物理プロセスである.しかし流域スケールで土壌の不飽和水分移動の支配方程式であるRichards式を解くには計算時間が膨大となってしまう.本研究ではRichards式の無次元化を行い流出解析モデルに導入可能なモデルの作成と共に計算の簡略化を図る.無次元のRichards式を用い,土壌水分の再分布及び表層土壌における水分の推移について検討を行った.
Key Words : 土壌水分再分布, 無次元Richards式, 懸垂水帯水分量, 土壌水分欠損量,地表面水分量

１．はじめに近年の異常気象により極端な降雨や長期的な降雨によって,多くの洪水等の災害が多発している.防災対策として,河川においては流出解析が用いられている.流出解
析は降雨及び地理特性を考慮し,河川での流出量を算出する解析であるが,精度の高い解析を行う際には土壌浸透を考慮しなければならない.土壌浸透は地下水への流出や,蒸発,流出量に影響する.浸透は土壌水分再分布による土壌水分プロファイルによって大きく変化し,高精度の流出解析には考慮しなければならない物理プロセスとなっている.また,地表面水分量は大気に直接接するため,蒸発や後続降雨の浸透,流出に密接に関係し重要な部分となっている.
土壌浸透についての研究は古くから行われており,Horton式（1930）をはじめとする経験式が多く提唱され1),不飽和土壌水分移動の支配方程式であるRichards式（1931）が発表された.1950年代よりRichards式に関する研究が行われ,浸透量,及び浸透速度を用いたPhilip式
や2)湿潤前線を用いるGreen and Ampt式等が提唱された2).
一般的な土壌浸透解析や多くの研究では帯水層から地表面水分・大気との境界までの水分移動を表現できる
Richards式が多く用いられている.しかし,Richards式は強い非線形を有するため,数値計算の負荷は大きく,流域スケールで解析を行う場合,細かい空間分解能で長時間の計算が必要になるため,流出解析モデルへの適用は現在でも難しい.また,検討対象土壌の詳細なパラメータや観測データを十分に把握しきることは時間及びコスト面においても問題が生じるため,十分な検証を用いたモデルへの組み込み事例は少ない.既往の研究3)では,土壌特性を考慮したモデル開発として,鉛直1次元Richards式の数値解に基づき,9種の土壌に対し深さ方向1ｍの土柱内の水分再分布を定量的に表現した.本検討では流出解析モ
デルへの導入を目標とし,Richards式の無次元化を図り,高分解能の数値解を用いて,土壌内の浸透・水分再分布・地表面水分量に注目し,それぞれの水分量変化と土壌パラメータとの関係性を見出した.無次元化を行うことで詳細な観測データを用いることなく解析が可能となり,考慮する土壌パラメータが減るため計算にかかる負荷も減少することが考えられる.

２．研究手法
本検討では深さ方向の鉛直1次元のRichards式を用いて無次元化を図り,数値解を得ることで土壌内水分再分布を表現する.式(1)に本検討で用いたRichards式の基本形を示す.
（1）

ここで は体積含水率, [T]は時間, [L]は上向きを正とする鉛直座標, [LT-1]は土壌の不飽和透水係数, [L]はマトリックスポテンシャルである.式(1)の数値解を得るためには不飽和透水係数 を求める必要があり の値はvan Genuchten(1980) 4)の式にMualem(1976) 5)6)の理論的モデルを用い導かれた次式から求めることができる.
.

図-1 無次元浸透モデル概略図（2）[ LT-1]は飽和透水係数, は有効飽和度, は飽和体積含水率, は残留体積含水率である.またはvan Genuchtenモデルで用いられる各土壌のパラメータである. は空気侵入圧の逆数に密接に関係し, (>1)は空隙系径分布の形状を表す. は一般的に用いられる式(2)を使用する.本研究ではこれらの各パラメータを無次元化した場合を考え,その有用性を検討する.土壌に対し,長さ及び時間の代表値をと取り決め無次元化を図る.ここで無次元化されたパラメータを ∗ , ∗ , ∗ と表現し,以下の式に示す.

∗∗

得られた無次元のパラメータを用いRichards式の無次元化を図る.
∗∗∗∗∗∗∗,,,∗∗

整数レベルでは有効飽和度 を求め,中間レベルでは不飽和透水係数 およびその地点におけるフラックスを定義づけ,解析の初期条件から次式を用いて算出を行う.

∗∗
∗

（3）

式(3)から図-1土壌浸透数値モデルに適用し,数値解析を行う.土注の深さ方向を刻み幅 ∗ により離散化し計算
時間 ∗ における水分移動の推定を行う.

ここで, は計算時間 ∗ を表す.無次元化されたRichards式を用いた解析の初期条件として深さ方向を∗とし,空間分解能 ∗.計算時間を ∗とし時間分解能を ∗とした.また,水分再分布の検討のため降雨直弧後の土壌を想定し解析を行うため, ∗の土壌を湿潤状態とした.上端及び下端はフラックス0として流入・流出は考えないものとして解析を行う.また,本研究では汎用性のある結果を求めるため,SandからClayまでの12種類の土壌に対し解析を行った.

表-1 対象土壌パラメータ
Soil typeSandLoam
Clay loam
Clayθs0.430.430.410.38

θr0.045
0.078
0.095
0.068

α(m -1 )
14.5
3.6
1.9
0.8

Ks(mhr -1 )0.29700.10400.00260.0020ｍ0.627
0.359
0.237
0.083

図-3 SW* およびSMD*の概念図

図-2 Sand土壌における土壌水分プロファイル(実線:Richards式による数値解, 点線:HYDRUS-1Dによる数値解)

12種の土壌パラメータはHYDRUS-1Dの一般的に用いられる値を使用する.本概要では代表土壌として4種のSand, Loam, Clay loam, Clayを対象に結果を掲載し,表-1にそれぞれの土壌パラメータを示す7).初期条件より無次元化された時間における土壌水分プロファイルを得る.
今回の無次元浸透モデル解析による浸透解析による結果とHYDRUS-1Dによる数値解の土壌水分プロファイルの比較図を図-2に示す.結果より,土壌中の水分移動が表現され,計算終了時には安定状態である静水圧平衡状態に近づくことが示される.また,比較を行うとHYDRUS1Dと同様の結果が得られ,無次元化による解析を行った場合でも有次元と同等の結果が得られることが分かった.

図-4 SW* の時間変化

３．土壌水分再分布指標の変化
本検討では,土壌水の指標として扱われている懸垂水帯水分量及び土壌水分欠損量を用いて検討を行う.図-3に各指標の概念図を示す.上層が飽和状態である初期土壌水分プロファイルから徐々に浸透が行われ安定状態に近づく.浸透過程は任意時間での土壌水分プロファイルであり, ∗ における有効飽和度の推移を示す,安定状態は,土壌水分が重力とマトリックスポテンシャルが釣り合っている状態を示し,静水圧平衡状態と呼ぶ.このとき図-3
中の土壌上部に存在する浸透状態と安定状態で囲まれた,時間経過に伴い次第に減少していく領域を ∗ 無次元懸∗は湿潤状態から静水圧平垂水帯水分量と定義する.
衡状態までの推移を表すため,土壌内水分移動の1つのバロメータとして扱うことができる8).また,飽和状態から∗任意時間における地下水面の位置までの総面積を9)無次元土壌水分欠損量として定義する .

図-5 SW* / SMD*の時間変化

各指標における式を以下に示す.
∗∗

∗

∗

∗

∗∗∗

ここで ∗ は安定状態を示す.
は土壌表面から底面までの各層ごとの土壌水分欠損量を足し合わせることで土壌内の乾燥度合を求めることができ,後続降雨の流入や蒸発量の指標となる.
図-4に無次元化されたRichards式によって得られた各土壌の ∗ 時間変化を示す.いずれの土壌も水分量が一∗∗様に減少する様子が見て取れる.
をにて除し,∗が指数的に低減す正規化を行う.結果を図-5に示す.
ることが示され,以下の推定式を提唱することができる.
∗∗∗∗ ∗図-6 aSW,bSWと土壌パラメータｍの関係

図-8 Sand土壌における1mm/dayの無次元蒸発速度を与えた
場合の土壌水分プロファイル

表-2 各土壌の解析結果
Soil Type
SandLoam Sand
Sandy Loam
Loam
SiltSilt Loam
Sandy Clay Loam
Clay Loam
Silt Clay Loam
Sandy Clay
Silt Clay
Clay

a SW
0.001340.001470.00177
0.00226
0.00257
0.00252
0.002410.00262
0.00256
0.002560.001680.00168b SW
0.13820.14170.14400.14640.15000.14900.14760.15180.15480.15480.15950.1595a SW10.004040.004090.004210.00428
0.003540.00370
0.003400.003420.003210.003240.003180.00320b SW11.05971.08581.06610.88600.62860.70230.91960.60510.5201
0.53280.29180.3002a SMD10.00327
0.00638
0.01393
0.02658
0.03172
0.03001
0.014310.03615
0.03646
0.038280.042020.04497b SMD17.01686.48655.57954.3742
3.49243.69024.01733.06492.53162.5005
1.1524
1.1274

図-9 1mm/dayの無次元蒸発速度考慮時のSW*,SMD*時間変化の土壌水分再分布について検討を行った.蒸発速度[ LT-1]は速度の次元を有するため,浸透速度である飽
和透水係数 [ LT-1]によって除すことで無次元蒸発速度 ∗ を得る.
∗

図-7 無次元懸垂水帯水分量比較図（実線:Richards式による蒸発速度を考慮することで表層圧力が低下するが,許とする.
で容最小圧力 を与え, をはフラックスが と等しい境界条件,に達した段階で一定圧力境界条件とする10).

数値解, 点線:推定式による数値解）

ここで,が無次元低減定数,が無次元初期懸垂水∗∗分量を示し,及び ∗ との線形回帰を行い12種の土壌の,を求めた.得られた,と土壌パラメータ との関係性を調べた.得られた結果を図-6および表-2に示す.各係数と推定式を用いて無次元Richards式より得られた数値解と比較を行い,図-7に示す.
比較を行った結果,Sand等の粒径の粗い土壌では高い再現精度を得ることはできなかったが,全体的に懸垂水帯水分量の時間変化は推定式によって再現を行うことができた.無次元化を図ることで,検討する土壌パラメータ∗
は のみとなり解析の簡略化が可能となり,における推定式も提唱することができた.

４．蒸発によるSW*への影響

また,本検討では,現在の日本における平均蒸発散である約2mm/day11)を想定し,1～3mm/dayの低い蒸発速度での土壌水分再分布を想定し検討を行った.また,本概要では1mm/dayの蒸発速度を与えた場合についての結果を示す.
図-8にSand土壌へ1mm/dayの無次元蒸発速度を与えた場
合の土壌水分プロファイルを示す.前章同様に1mm/dayの蒸発速度を与えた場合の ∗ 時間変化を図-9に示す.図4と比較すると蒸発による影響により土壌水分が計算開∗始と共に急激に推移し,が減少することが示された.
∗
また,蒸発による排水が行われることにより,の値は徐々に増加する.図-9において両対数にて表現を行う∗∗
は指数的,は累乗的に変化すると考えことでられ,以下の推定式で表現することができる.
∗∗∗

前章では上端及び下端からのフラックスを0とし,流入・流出がない状態での検討を行う.本章では上端を開放し,蒸発を考慮する.土壌表層へ蒸発速度を与えた場合∗,,,と土得られた蒸発考慮時の壌パラメータ との関係性を得ることで簡易的に蒸発考∗慮時の ∗ ,を得ることが可能となる.前章同様に線形回帰を行い,推定式による再現性を評価する.

図-10 aSW1, bSW1と土壌パラメータｍの関係
図-13 地表面水分量時間変化
表-2 各土壌の解析結果

図-11 aSMD1,bSMD1と土壌パラメータｍの関係

Soil Type
SandLoam Sand
Sandy Loam
Loam
Silt
Silt Loam
Sandy Clay Loam
Clay Loam
Silt Clay Loam
Sandy Clay
Silt Clay
Claya Sr
1.06170.87360.70250.56300.47670.49590.52790.44610.39980.3998
0.2836
0.2836b Sr
0.437480.42660
0.425900.440410.45682
0.45298
0.44669
0.462770.470020.470020.47186
0.47186

図-12 蒸発考慮時SW*及びSMD*比較図（実線:Richards式による数値解, 点線:推定式による数値解）比較を行った結果,閉鎖状態の土壌の場合に比べ再現精度は良好な結果を得られなかった.これらは計算開始直後から ∗ が大きく減少するため,その傾向を推定式では捉えきれていないためと考えられる.SandやClayなどの粒径が極端な土壌以外では傾向は捉えることができ,精度の高い検証が行うことができる可能性が示された.
また,無次元化されたRichards式によって蒸発を考慮した場合での土壌水分再分布が推定可能となり,推定式を用いることで簡易的に求めることができると分かった.

図-14 地表面水分量比較図（実線:Richards式による数値解,点線:推定式による数値解）

地表面水分量は比較的早い ∗ において有効飽和度が推∗に収束すると考えられる.ここで,移しており,∗までの区間に注目し, にて正規化を行う.
,,

５．地表面水分量変化
地表面における土壌水分は直接大気に接する面であり,後続降雨の流入や蒸発等による流出が行われるため,重要な物理量となる.一般的に地表面水分量は降雨,流入後に湿潤となり,浸透・流出が行われることで徐々に減少してゆく.本研究ではモデル第1層 ∗における有効飽和度 の時間変化を示す.図-13に上端及び下端フラックスが無い場合の地表面水分量の推移を示す.各土壌の減少傾向から浸透による減少,一定値に収束が確認できる.一定値は図-3に示される安定状態,静水圧平衡状態の有効飽和度 ∗ と考えられる.

∗!

,各パラメータ , は正規化を行った式を変形し,累乗近似の形を取り,非線形解析を用い算出する.得られた, を表-3に示す.各係数と土壌パラメータ の関係を見ることで簡易的に地表面水分量の推定が可能となる.
, が推定可能となり地表面各係数の回帰式を得る.
における有効飽和度,体積含水率が算出可能となる.図-7
同様に推定式による結果と無次元Richards方程式の数値解から得られた結果を図-14に示す.各土壌において地表面水分量の時間変化をよい精度で再現することができ,簡易的に地表面水分量を推定することができたと言える.

図-15 蒸発考慮時地表面水分量時間変化

６．蒸発による地表面水分量変化
本章では4章同様に蒸発を考慮した場合の地表面水分量について検討を行う.初期条件は4章同様に上端を開放し,無次元蒸発速度を表層へ与える.蒸発を考慮することでより現実に近い推定を行うことができる.図-15に1mm/dayの無次元蒸発速度を与えた場合の地表面水分量の変化を示す.計算開始から有効飽和度は急激に減少し,安定状態である各土壌におけるそれぞれの有効飽和度に収束することがわかる.また,Sand土壌では比較的緩やかに有効飽和度は減少し,粒径の細かい土壌になるほど収束までの時間は短い.粒径の粗い土壌では土壌水分の移動が容易であり比較的深い位置からでも水分が供給できるため,有効飽和度の収束まで時間がかかり,表層は乾燥状態に近づく.一方粒径の細かい土壌では,土壌が密な形状を取るため水分を保持する力が強いことから土柱上層の一部の水分のみが蒸発し,それ以外の土壌水分は保持されているため,収束までの時間が速いと考えられる.1ｍｍ/day以上の蒸発を考慮した場合,各土壌においてより短時間で有効飽和度が一定の値に収束することが確認された.

また,蒸発を考慮した場合の土壌水分プロファイル,懸垂水帯水分量,土壌水分欠損量が算出可能となり,より現実的な解析が可能となった.各指標における推定式を提唱し,各低減係数を土壌パラメータ と関連付けることにより,簡易的に蒸発を考慮した場合の土壌水分再分布が推定可能となった.しかし,現段階の推定式は簡易的なものであり高い再現精度とは言えないため,今後の課題となる.また,本検討では低い蒸発速度のみに焦点を当てているため,蒸発量の大きい場合や植生などによる損失を考慮した場合の検討が今後必要となる.
土壌パラメータ を用いることで土壌水分再分布の指標である懸垂水帯水分量,土壌水分欠損量が推定式により簡易的に推定が可能となった.また,表層における有効飽和度を算出可能となり,後続降雨への流入,蒸発や流出をモデル化し,流出解析へ応用できる可能性が示された.
参考文献
1) Jury,W.A.and Horton,R.(2006):土壌物理学—土中の水·熱·ガス·化学物質移動の基礎と応用—(取出伸夫 監訳: 井上光弘，長裕幸，西村拓，諸泉利嗣，渡辺晋生訳),pp.36–159,築地書店，
東京2) Smith, R. E., K. R. J. Smettem, P. Broadbridge and D. A. Woolhiser:
Infiltration theory for hydrologic applications, American GeophysicalUnion, Vol.15, 2002.
3) 広田浩平(2012)：土壌水分再分布を考慮した土壌浸透モデルの開発. 平成24年長岡技術科学大学修士論文.
4) Mualem,Y. 1976. A new model for predicting the hydraulicconductivity of unsaturated porous media. Water Resour. Res.
12:513-522.
5) 小杉賢一郎(2007)：古典を読む Y.Mualem著「不飽和多孔質体の透水係数を推定する新たなモデルについて」ならびに
M.Th.van.Genuchten著「不飽和土壌の透水係数を推定する閉７．結論

形式解について」. J.Jpn.Soc.Soil Phys. 106:47-60.
6) van Genuchten,M.Th. 1980. A close form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Am.J. 44:892-

本検討では,Richards方程式を基礎方程式として,鉛直方向の無次元Richards式を導出,土壌浸透モデルの構築を行った.無次元Richards式の有用性を示すとともにSandからClayまでの12種の様々な土壌において土壌水分プロ
ファイル,懸垂水帯水分量,土壌水分欠損量の時間変化を表現することができた.
土柱の上端及び下端からの流入,流出が無い閉鎖状態の場合,懸垂水帯水分量は指数的に減少し,静水圧平衡状態へ推移することが明らかになった.推定式を提唱し,低減係数と土壌パラメータ との関係性を見出すことで,簡易的に懸垂水帯水分量の時間変化を表現することができた.また,地表面水分量の推移をモデルによって表現し,地表面水分量の減少を にて正規化を行った.地表面水分量の減少初期における減少過程を,土壌パラメータを用いることで精度よく簡易的に再現することができた.

898.
7) Simunek, J., M. Th. van Genuchten, and M. Sejna, The Hydrus-1D
software package for simulating the one-dimensional movement ofwater, heat, and multiple solutes in variably-saturated media. Version3.0, HYDRUS Software Series 1, Department of Environmental
Sciences, University of California Riverside, Riverside, CA, 270 pp.,2005.
8) 藤原 俊六郎（2013）：図解 土壌の基礎知識.農山漁村文化協会.
9) 金子真郷(2010)：分布型水文モデルへ導入が可能な簡略化した土壌とモデルの開発. 平成22年長岡技術科学大学修士論文.
10) 齋藤広隆・坂井勝・Jiri Simunek・鳥出信夫（2006）：不飽
和土中の水分移動モデルにおける境界条件
11) 国土交通省関東地方整備局:http://www.ktr.mlit.go.jp/(2016/01).
（2016.3.11受付）