第66回(2022年度)

北海道開発技術研究発表会論文

泥炭土を基盤とする農地における水移動シミュレーション―蒸発散時の事例紹介―
国立研究開発法人 土木研究所 寒地土木研究所 資源保全チーム ○奥田 涼太
長竹 新
清水 真理子暗渠が設置された泥炭土を基盤とする大豆畑の無降雨期間を対象に、二次元の水分移動予測
プログラムHYDRUS-2Dで水分移動のシミュレーション（再現計算）を行った。実測の土壌水分変化に基づき逆解析で求めたパラメータを入力することで、計算値は実測値をよく再現した。
なお、鉛直方向の水分移動が支配的な蒸発散時には、無償で利用可能な1次元のHYDRUS-1Dでの計算が暗渠設置圃場において有効であると示唆された。
キーワード：農地、水分移動、シミュレーション

1. はじめに

0m

70 m

220 m

調査地点

暗渠管

図-1 調査圃場の模式図
土中水圧力・土壌水分量測定
0①作土

2. 現地調査
調査は北海道今金町の転作田（約250 m×60 m）の暗渠から水平方向に5 m離れた地点で行った（図-1）。
(1) 試料採取と断面調査
2020/8/18に試料採取と断面調査を行った。断面調査の結果、0～20 cm深は作土層、20～30 cmは圧密された作土
OKUDA Ryota, NAGATAKE Arata, SHIMIZU Mariko-154-

水閘

5 m（計算領域）

排水路取水桝

近年HYDRUS1)のような、土中の水・熱・ガス・化学物質の移動予測汎用ソフトウェアで水分移動を手軽に計算できるようになった。しかし、HYDRUSなどこれらのソフトウェアは、海外の製品であることから、技術者からの認知度が低いのが現状である。適切な水分移動に関するパラメータと、初期・境界条件をソフトウェアに入力し、水分移動のシミュレーションを行うことで、土中水圧力の時間変化や、任意時刻の土中水圧力の分布、土中水の移動速度と方向を表す水フラックス等の結果を得ることが出来る。こうした結果は、日消費水量の算出、地下水位制御システムによる地下灌漑時の水分変化の可視化、暗渠管埋設間隔の検討などに有用であろう。
本報では、今金町の水田転作畑における大豆畑の蒸発散を例に、HYDRUS-2Dを用いた水分移動シミュレーシ
ョンの事例を紹介する。

層、30～40 cmは作土と泥炭土が混じった層、40～60 cm
は地下水位変動の影響で乾燥収縮や分解が生じている泥炭土層、60～100 cmは一年の大部分で地下水面下にある泥炭土層だった（図-2）。層別に100 cm3の不撹乱土を3連で採取した。

⑥暗渠
疎水材
暗渠

②作土
③作土 泥炭土混在
④泥炭土
⑤泥炭土

500cm
図-2 調査圃場断面模式図

20406080100cm

3. 水分移動特性パラメータのフィッティング
水分移動のシミュレーションを行うためには、飽和透水係数や、土中水圧力と土壌水分量の関係といった土の水分移動特性を数値で表す必要がある。そのために、まずは現地調査で採取した不撹乱土を用いて室内試験を行い、試験結果に基づいて水分移動特性モデルにフィッティングしパラメータを決定した。
(1) 飽和透水係数の決定現地調査で採取した不撹乱土の変水位飽和透水試験を行い、不撹乱土の飽和透水係数Ks（cm day –1）を得た。
(2) 土中水圧力と土壌水分量の関係
飽和透水試験後、吸引法と加圧板法で不撹乱土の水分を調整の上、脱水過程における土中水圧力hと土壌水分量θを測定し相互の関係を得た。
(3) モデルへのフィッティング
不撹乱土、テンシオメータと土壌水分計の実測値から得られたhとθの関係を基にvan Genuchten式(1a)のパラメータを決定した。

rs

K h

Se

1hnmr

K SeleKs S 1

1 S

1m me2

3章で不撹乱土を用いて水分移動特性モデルのパラメータを決定したが、100 cm3の不撹乱土が圃場を正確に表わす事は少ない。より正確に圃場のシミュレーションを行うためには、土中水圧力hや土壌水分量θの実測値を用いて逆解析をすることが有効である。今回のような蒸発散時には上向きの水分移動が支配的であることから、暗渠から水平方向に5 m離れた地点では横方向の水移動が無視できると仮定し、鉛直1次元方向のみの水分移動を計算するHYDRUS-1Dの逆解析機能を用い、土中水圧力hの実測値に基づいてパラメータを決定した。HYDRUS1Dの計算条件は以下のとおりである。
(1) 水分移動特性パラメータ3(3)で求めた式(1a)のパラメータを逆解析の初期値として入力した。推定するパラメータは、①、②の作土層はα、n、Ksで、③の作土と泥炭土が混在した層、④、⑤の泥炭土層はKsとした。
(2) 計算領域の作成
HYDRUS-1D上でz軸（鉛直）方向に100 cmの計算領域を作成し，断面調査の結果から計算土層を5つに分けた。
(3) 初期圧力条件
2021/7/19 0:00における土中水圧力hの実測値に基づき、測定深度間の値を直線で補完して初期条件を入力した。
(4) 上端境界条件
上端の境界条件には蒸発速度、蒸散速度を入力した。
30分間隔の土壌水分変化量は小さかったため、まず土壌水分減少法2)によって土壌水分量θの値から日消費水量求めた。ただし、土中水圧力hと位置水頭zの和である全水頭Hの傾きと地下水位の値から、何cm深までの土壌水分減少が蒸発散に寄与しているかを考慮した。日消費水量は1.0～3.7 mmの間で推移し、期間中の平均値は2.3 mmだった。太陽から圃場に供給される正味の熱エネルギーで3-40(1a)–1

ここで、K(h)は透水係数（cm day ）、Seは有効飽和度、
θr 、θs は残留および飽和土壌水分量（cm3 cm–3 ）、α（cm–1）、n、m（= 1 – 1 / n）、lは適合パラメータである。
ただしlは0.5で一定とし、暗渠疎水材については過年度に求めたパラメータを用いた。

2.5-50地下水位

21.5
1-60-70

蒸発散速度
蒸散速度
蒸発速度-800.507/19-90-1007/247/29

図-3 蒸発散速度と地下水位
OKUDA Ryota, NAGATAKE Arata, SHIMIZU Mariko-155-

地下水位, WL (cm)(3) ドローンによる空中写真撮影
2021/7/18と2021/8/11にドローン（DJI, Phantom 4 RTK）
を用いて高度50 mから空中写真を撮影した。

4. 水分移動特性パラメータの逆解析

蒸発散速度, ETa (cm day–1)(2) 土中水圧力、土壌水分量、地下水位、放射収支測定
同一地点で大豆栽培中の無降雨期間である2021/7/19～30にテンシオメータで土中水圧力h（cm）、TDR土壌水分計で土壌水分量θ（cm3 cm– 3）、地下水位計で地下水位（cm）、放射収支計で放射収支（W m 2）を30分間隔で測定した（図-2）。hは5、15、25、35、45、65、85 cm深を、θは0～10、10～20、20～30、30～40、40～70 cm深毎の平均値を測定した。

ある純放射の値と蒸発散速度に相関があると仮定し、30分間隔の放射収支計の値に応じて日消費水量を按分したものを、30分間隔の蒸発散速度とした。
HYDRUSユーザーマニュアル1)に従い、圃場の表面被覆率を用いて蒸発散速度を蒸散速度と蒸発速度に分離した。図-3に蒸発散速度、分離した蒸散速度と蒸発速度を示す。ここで表面被覆率とは、圃場を真上から見たときの（作物の面積）/（圃場の面積）である。ドローンで撮影した写真をソフトウェア（Agisoft, Metashape）を用いて1枚の画像データにし、プログラミング言語（Python）の画像処理・解析ライブラリであるOpen-CV
を 用 い て 圃 場 の表 面 被 覆率 を 求 め た 。 2021/7/18 と2021/8/11の表面被覆率をそれぞれ56.9、93.4%だった。ドローンでの撮影は2日のみだったため、期間中の表面被覆率は線形で増加したと仮定した。

5. 数値計算条件の入力
HYDRUS-1Dの逆解析機能で決定した圃場を表すパラメータを用いて、HYDRUS-2Dで水分移動をシミュレーションした時の計算条件を以下に示す。

(1) 計算領域の作成
HYDRUS-2D上でx軸（水平）方向に500 cm、z軸（鉛直）方向に100 cmの長方形を計算領域として作成し，断面調査の結果から計算土層を5つに分けた（図-2）。領域左端の25～80 cm深に暗渠疎水材を設定し，直径10 cmの暗渠を管心が85 cm深になるように設定した。
(2) 水分移動特性
逆解析で求めたパラメータ（表-1）を入力した。

(5) 下端境界条件
下端の境界条件には地下水位（図-3）から変換した
100 cm深における土中水圧力を入力した。今回は計算深度が全領域で100 cmなので、100 cmから地下水位（cm）を引いたものが100 cm深における土中水圧力である。地下水位は期間中–44 cmから–64 cmにほぼ線形に減少し、傾きは–1.6 cm dayだった。

(3) 初期圧力条件、上下端境界条件、根の分布
土中水圧力の測定は1地点のみ行ったので、横方向の分布は一定であると仮定して、逆解析時と同様の値を入力した（図-4）。上下端の境界条件と根の分布も逆解析時と同様の値を入力した

6) 根の分布
今回根密度調査は行わなかったため、根の分布は地表面から作土層の下端である30 cm深まで一定密度で存在したと仮定して与えた。

(4) その他の境界条件
左右端の境界条件は圃場の対称性より流入出なしとし、暗渠の境界条件は浸出面（土中水圧力hが正の時、境界の圧力を0にする）とした。

① 0～20 cm② 20～30 cm③ 30～40 cm④ 40～60 cm⑤ 60～80 cm⑥ 暗渠疎水材

表-1 逆解析から求めた van Genuchten 式のパラメータθrθsαn(cm3 cm–3)(cm3 cm–3)(cm–1)(-)
0.00141
0.475
1.39×10-3
1.57
0.001
0.540
3.88×10-3
1.17
0.001
0.826
6.00×10-3
1.63
0.01404
0.8465
4.83×10-3
1.38
0.01404
0.8465
4.83×10-3
1.38
0.092
0.7725
0.585
1.26-2500 -2250 -2000 -1750 -1500 -1250 -1000 -750 -500

-250

0Ks
(cm day1)
1.7×10-2
6.5×10-3
2.5×10-1
1.2×10-2
1.1×10-3
4.2×104

250（cm）

020

暗渠疎水材

地下水面

4060

暗渠

80
100

0100

200300図-4 初期土中水圧力分布図（7/19 0:00）

OKUDA Ryota, NAGATAKE Arata, SHIMIZU Mariko-156-400500（cm）

(3) 水フラックスの分布
図-7に7/25正午の水フラックスの分布を示す。矢印の方向は水の向きを、色は水の移動速度を表している。蒸HYDRUS-2Dで初期からの水分移動を計算することで、 発散時は主に地下水面より上の層で水分移動が生じてい土中水圧力の時間変化や、任意時刻の土中水圧力の分布、 た。暗渠からおおよそ50 cm以上離れた領域では鉛直方向の水分移動のみが生じていた。このように、鉛直方向土中水の移動速度と方向を表す水フラックスを求めるこの水分移動が支配的な蒸発散時には、1次元のHYDRUSとが出来る。蒸発散時の水分移動をシミュレーションし1Dでの計算が暗渠設置圃場においても有効であると考た結果と考察を以下に示す。
えられた。今回はHYDRUS-1Dを用いて逆解析を行ったが，大雨後の排水や地下灌漑時の水分移動をシミュレーション(1) 土中水圧力の時間変化
する際には暗渠からの排水量や給水量が蒸発散時と比較
図-5に土中水圧力の実測値（点線）と計算値（実線）して大きくなり，計算領域全体で水フラックスの横向きを示す。決定係数R2は0.7で、計算値は実測値をよく再成分が大きくなると考えられる。そのため、排水や地下現した。計算期間中、土中水圧力hはいずれの深度でも灌漑時の水分移動をシミュレーションするための逆解析減少傾向だった。一番乾燥した7/31の値を見ると、5、35cm深でおおよそh = –800、–100 cmだった。これらは、植はHYDRUS-2Dで行う必要 があるだろう。
物が有効に利用できる水分量の指標である有効水分量の0範囲内（pF = 1.5～3、h = –32～–1000 cm）であり、乾燥や35湿潤による蒸発散速度の制限は受けていないと考えられ-20025た。
-40015(2) 土中水圧力分布-600
図-6に7/25正午の土中水圧力hの分布を示す。地表面で5 cm深蒸発が、0～20 cm深で根から吸水が生じるため、上端ほ-800
どhが小さく乾燥しており、下端に向かうにつれてhが大-1000
きく湿潤となる分布になった。また、暗渠から50 cm程
7/197/227/25
7/287/31度の領域を除き、hの分布は横方向に一定で、地下水位
図-5 土中水圧力の経時変化
も61 cm深でほとんど横方向に一定だった。
土中水圧力, h (cm)6. 数値計算の結果と考察-2500 -2250 -2000 -1750 -1500 -1250 -1000 -750 -500

-2500250（cm）

020406080
1000100200300

400500（cm）

図-6 土中水圧力分布図（7/25 12:00）
00.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.911.1（cm d-1）

0204060801000100200300

図-7 水フラックス分布図（7/25 12:00）
OKUDA Ryota, NAGATAKE Arata, SHIMIZU Mariko-157-400500（cm）

謝辞：本研究の現地試料採取にあたり、土地所有者、関係機関の各位にご協力を賜った。ここに記して謝意を表する。

7. まとめHYDRUS-2Dの水分移動シミュレーション機能の事例紹介として、今金町の水田転作畑における大豆の蒸発散を例に示した。HYDRUS-2D上で圃場を模した計算領域を作成し、逆解析で求めた水分特性曲線のパラメータ、初期圧力条件、境界条件、根の情報を入力することで、水分移動の結果である土中水圧力分布や、水分移動そのものである水フラックスの分布図が得られた。土中水圧力の結果からは、計算期間中は有効水分量の範囲内であり、乾燥や湿潤による蒸発散速度の制限は受けていないと考えられた。水フラックスの結果からは、鉛直方向の水分移動が支配的な蒸発散時には、HYDRUS-1Dでの計算が暗渠設置圃場においても有効であることが示唆された。

参考文献
1) Riverside, California: The HYDRUS-1D Software Package for Simulating
the One-Dimensional Movement of Water, Heat, and Multiple Solute in
Variably-Saturated Media
2) 農林水産省：8. 計画日消費水量等の算定、土地改良事業設計
基準 計画「農業用水（畑）」技術書、176～196.

OKUDA Ryota, NAGATAKE Arata, SHIMIZU Mariko-158-