HYDRUS-1D(Ver.4)を 用 いた 蒸 発 法 による 水 分 移 動 特 性 の 推 定 2008.5.11 HYDRUS-1D には, 土 中 の 圧 力 変 化 などの 測 定 値 に 基 づき 水 分 移 動 特 性 (Soil Hydraulic Parameters)を 逆 解 析 する 機 能 がある ここでは, 土 カラムからの 水 分 蒸 発 過 程 の 圧 力 変 化 から 逆 解 析 する 蒸 発 法 への 適 用 例 を 示 す 例 題 1 土 カラムからの 水 分 蒸 発 と 水 分 移 動 特 性 の 逆 解 析 Evapor ここでは,HYDRUS-1D サンプルプロジェクト Evapor (ワークスペース Inverse 中 )を 用 いて 解 説 する 測 定 項 目 深 さ 10cm 土 層 からの 蒸 発 過 程 における 圧 力 水 頭 ( 深 さ 1,3,5,7,9cm), 蒸 発 速 度, 実 験 終 了 時 平 均 体 積 含 水 率 推 定 項 目 van Genuchten モデルパラメータθ r,θ s,α,n,k s (1) 主 要 な 計 算 条 件 Main Processes Heading: プロジェクトの 見 出 し Simulate: 水 分 移 動 (Water Flow)を 選 択 Inverse Solution: 逆 解 析 (パラメータの 推 定 ) 選 択 (2) 逆 解 析 Inverse Solution Estimate: 水 分 移 動 特 性 (Soil Hydraulic Parameters)の 推 定 Weighting of Inversion Data: 逆 解 析 の 重 み 関 数 通 常 は weighting by mean ratio 選 択 Maximum Number of Iterations: 反 復 回 数 の 上 限 100 程 度 に 設 定 する 0 のときはパラメータの 推 定 は 行 わず, 初 期 値 による 順 解 析 として 計 算 Number of Data Points: 逆 解 析 に 用 いる 総 データ 数 ( 入 力 ファイル Fit.In 中 の NOBB) データは Data for Inverse Solution で 入 力 する データ 入 力 時 にデータ 数 を 変 更 したときは, 自 動 的 に 変 更 される 1
(3) 土 層 に 関 する 情 報 Main Geometry Information Length Units: 長 さ 単 位 Number of Soil Materials: 計 算 に 用 いる 土 の 種 類 の 数 ( 均 一 試 料 に 対 する 逆 解 析 では 1) Number of Layers for Mass Balances: マスバランスを 計 算 する 層 の 数 マスバランスの 計 算 は, 土 の 層 と 一 致 する 必 要 はなく,いくつかの 土 による 構 成 でも 良 い 後 述 のグラフィカルエディターProfile では,マスバランス 領 域 を Subregion とよぶ Decline from Vertical Axes: 流 れ 方 向 角 度 x の cos x 値 1: 垂 直 流 れ; 0: 水 平 流 れ Depth of the Soil Profile: 計 算 する 土 層 の 長 さ (4) 時 間 情 報 Time Information Time Units: 時 間 単 位 Time Discretization: 離 散 化 時 間 の 設 定 Initial Time: 初 期 時 間 Final Time: 計 算 終 了 時 間 Initial Time Step: 初 期 時 間 刻 み 幅 Minimum Time Step: 最 小 時 間 刻 み 幅 計 算 の 非 線 形 性 が 強 いときには, 数 値 計 算 の 収 束 が 困 難 にな る( 水 分 特 性 曲 線 の 傾 きが 急 な 場 合 等 ) この 場 合 は, 初 期 時 間 刻 み 幅, 最 小 時 間 刻 み 幅 を 極 端 に 小 さくする 必 要 がある( 砂 丘 砂 の 蒸 発 法 の 場 合, 数 日 間 の 計 算 に 対 してそれぞれ 10-5,10-8 日 の 設 定 が 必 要 であった) 2
Maximum Time Step: 最 大 時 間 刻 み 幅 Boundary Conditions: 表 面 境 界 条 件 である 蒸 発 速 度 は 時 間 により 変 化 するので Time-Variable Boundary Conditions を 選 択 Number of Time-Variable Boundary Records: 時 間 変 化 の 境 界 条 件 として 入 力 のデータ 数 データは Time Variable Boundary Conditions で 入 力 する データ 入 力 時 にデータ 数 を 変 更 したときは, 自 動 的 に 変 更 される (5) 印 刷 情 報 Print Information Print Options 印 刷 オプション T-Level Information: 選 択 すると, 計 算 時 間 に 対 する 平 均 水 分 量, 平 均 水 分 溶 質 フラックス, 積 算 水 分 フ ラックス, 繰 り 返 し 計 算 等 に 関 する 情 報 が,ファイル(T_LEVEL.OUT)に 出 力 される Every n time steps に 値 を 入 力 することで, 出 力 データを 間 引 くことが 可 能 Print at Regular Time Interval: 入 力 した 時 間 間 隔 で T-Level Information を 出 力 Screen Output: 選 択 すると, 計 算 実 行 中 に 画 面 に 計 算 結 果 が 出 力 される 非 選 択 の 場 合, 反 復 計 算 中 の 水 分 移 動 特 性 パラメータ 値 が 画 面 に 表 示 される Print Fluxes for Observation Nodes: 選 択 すると, 観 測 点 におけるフラックスを 出 力 Hit Enter at End: 選 択 すると, 計 算 終 了 時 に Enter キーで 実 行 画 面 を 閉 じる 計 算 終 了 時 に 表 示 される 計 算 所 要 時 間 を 確 認 することができる Print Times 印 刷 時 間 Number of Print Times: 出 力 ファイル(NOD_INF.OUT)にデータ ( 圧 力 水 頭, 含 水 率, 濃 度, 水 分 フラッ クス, 水 分 収 支, 溶 質 収 支 等 )を 出 力 する 印 刷 時 間 数 データ 数 を 入 力 後,Select Print Times を 選 択 して 印 刷 時 間 を 入 力 する 3
Select Print Times: 印 刷 時 間 の 入 力 Default を 選 択 すると 等 間 隔 時 間 となる Default(log)を 選 択 すると 対 数 時 間 の 等 間 隔 となる 計 算 開 始 時 間 は 含 めず, 印 刷 時 間 の 最 大 値 は 計 算 終 了 時 間 と 同 値 であること (6) 反 復 計 算 Iteration Criteria Iteration Criteria: 反 復 計 算 に 関 する 条 件 設 定 p.1.8 参 照 Maximum Number of Iterations: 各 時 間 ステップでの 最 大 反 復 回 数 10 程 度 Water Content Tolerance: 不 飽 和 流 れの 水 分 量 許 容 誤 差 0.001 程 度 Pressure Head Tolerance: 飽 和 流 れの 圧 力 水 頭 許 容 誤 差 1cm 程 度 Time Step Contorol: 時 間 ステップ 制 御 Lower Optimal Iteration Range: 収 束 回 数 がこの 値 以 下 では 時 間 刻 み 幅 を 増 加 させる 推 奨 値 3 Upper Optimal Iteration Range: 収 束 回 数 がこの 値 以 上 では 時 間 刻 み 幅 を 減 少 させる 推 奨 値 7 Lower Time Step Multiplication Factor: 時 間 刻 み 幅 の 減 少 割 合 推 奨 値 1.3 Upper Time Step Multiplication Factor: 時 間 刻 み 幅 の 増 加 割 合 推 奨 値 0.7 Iternal Interpolation Tables: HYDRUS は, 計 算 を 行 う 前 に, 入 力 された 水 分 移 動 パラメータに 基 づく, 水 分 量, 不 飽 和 透 水 係 数, 比 水 分 容 量 の 表 を 作 成 する そして, 反 復 計 算 の 際 には, 作 成 した 表 の 該 当 部 分 を 線 形 補 完 することで 水 分 移 動 特 性 値 を 決 定 する この 補 間 表 の 利 用 により, 計 算 処 理 を 早 めることができる この 補 間 による 計 算 を 行 いたくない 場 合 は, 補 間 の 圧 力 水 頭 上 下 限 値 の 両 方 に 0 を 入 力 する Lower Limit of the Tension Interval: 1e 006 程 度 Upper Limit of the Tension Interval: 10000 程 度 4
この 反 復 計 算 に 関 する 設 定 は, 計 算 に 問 題 が 生 じない 限 り,デフォルト 設 定 値 を 用 いた 方 がよい (7) 土 の 水 分 移 動 モデル Soil Hydraulic Model Hydraulic model: 水 分 移 動 特 性 モデルの 選 択 ここでは van Genuchten-Mualem を 選 択 (8) 水 分 移 動 特 性 パラメータの 逆 解 析 に 関 する 設 定 Water Flow Parameters-Inverse Solution Initial Estimate: 水 分 移 動 特 性 パラメータ 推 定 の 初 期 値 van Genuchten 式 の 場 合,Qr = θ r, Qs = θ s, Alpha = α,n,k s,l の 値 を 入 力 パラメータ 推 定 の 際 の 初 期 値 は, 適 切 な 推 定 値 を 得 るために 重 要 である 水 分 保 持 曲 線, 及 び 飽 和 透 水 係 数 の 実 測 値 に 基 づく 情 報 がある 場 合 は,その 値 を 用 いる 非 線 形 最 小 2 乗 法 では, 測 定 値 と 計 算 値 に 関 する 残 差 の 2 乗 和 の 極 小 値 を 与 える パラメータの 組 み 合 わせを 推 定 値 として 求 める 極 小 値 がいくつか 存 在 する 場 合 には, 非 現 実 的 なパラメータの 組 み 合 わせが 推 定 結 果 となることも 少 なくない そのため 初 期 値 については, 推 定 結 果 を 確 認 しながら, 試 行 錯 誤 する 必 要 がある なお 測 定 データに 対 して 不 合 理 な 値 を 与 えると, 数 値 的 なエラーが 生 じる 場 合 が 多 い 極 端 に 大 きな( 小 さな)パラメータ 値 (たとえば n = 10)も, 計 算 が 破 綻 しやすい ここではカタログの Clay Loam の 値 を 初 期 値 として 用 いる 5
Minimum Value: パラメータ 推 定 の 最 小 値 Maximum Value: パラメータ 推 定 の 最 大 値 最 小 値 = 最 大 値 を 入 力 すると, 推 定 の 際 のパラメータの 範 囲 は 設 定 されない ( 下 表 の Alpha) Fitted: 推 定 するパラメータを 選 択 ここでは,θ r,θ s,α,n,k s を 選 択 (9) 水 分 流 れの 境 界 条 件 Water Flow Boundary Conditions Upper Boundary Condition: 蒸 発 法 の 上 部 境 界 は 蒸 発 フラックスを 与 える Atmospheric BC with Surface Layer を 選 択 蒸 発 法 では,カラム 重 量 の 測 定 値 から 蒸 発 フラックスを 算 出 する この 条 件 を 選 択 し た 後, 測 定 した 蒸 発 フラックスの 経 時 変 化 は,Time Variable Boundary Conditions で 入 力 する Lower Boundary Condition: 下 端 はゼロフラックスなので,Constant Flux を 選 択 Initial Condition: 初 期 条 件 を 圧 力 水 頭 あるいは 水 分 量 で 与 えるかの 選 択 テンシオメータで 圧 力 を 測 定 している 場 合 は, 初 期 条 件 は 圧 力 で 与 える 一 定 境 界 フラックス Constant Boundary Fluxes Lower Boundary Flux: 下 端 Constant Flux の 条 件 0 を 入 力 6
(10) 変 動 境 界 条 件 Time Variable Boundary Conditions この 画 面 を 通 じて 境 界 条 件 を 入 力 する データの 数 は Time information 中 の Number of Time-Variable Boundary Records で 与 えられている Time: 経 過 時 間 Precip.: 降 水 フラックス( 蒸 発 法 ではゼロ) Evap.: 蒸 発 フラックス( 蒸 発 法 ではカラム 重 量 の 減 少 量 から 求 める) hcrita: 計 算 上 許 容 する 表 面 の 最 小 圧 力 水 頭 の 絶 対 値 十 分 に 乾 燥 した 条 件 として 100000 程 度 を 用 い る (11) 逆 解 析 に 用 いるデータ Data for Inverse Solution 逆 解 析 に 用 いるデータの 入 力 測 定 データの 種 類 に 応 じてデータタイプを 分 類 し, 入 力 する データの 数 は,Inverse Solution 中 の Number of Data Points で 与 えられる Type =0: 境 界 面 での 積 算 フラックス = 1: ある 観 測 点 における 圧 力 水 頭 ( 圧 力 水 頭 の 経 時 変 化 に 対 して 逆 解 析 ) = 2: ある 観 測 点 における 体 積 含 水 率 ( 体 積 含 水 率 の 経 時 変 化 に 対 して 逆 解 析 ) = 3: 境 界 面 でのフラックス = 4: ある 観 測 点 における 濃 度 ( 濃 度 の 経 時 変 化 に 対 して 逆 解 析 ) = 5: 水 分 特 性 曲 線 h(θ )の 測 定 値 7
= 6: 不 飽 和 透 水 係 数 K(h)の 測 定 値 = 7: パラメータα の 予 測 値 = 8: パラメータ n の 予 測 値 = 9: パラメータθ r の 予 測 値 = 10: パラメータθ s の 予 測 値 = 11: パラメータ K s の 予 測 値 = 12: ある 印 刷 時 間 における,ある 位 置 の 圧 力 水 頭 ( 圧 力 水 頭 プロファイルに 対 して 逆 解 析 ) = 13: ある 印 刷 時 間 における,ある 位 置 の 体 積 含 水 率 ( 体 積 含 水 率 プロファイルに 対 して 逆 解 析 ) = 14: ある 印 刷 時 間 における,ある 位 置 の 濃 度 ( 濃 度 プロファイルに 対 して 逆 解 析 ) = 15: ある 印 刷 時 間 における,ある 位 置 の 吸 着 濃 度 ( 吸 着 濃 度 プロファイルに 対 して 逆 解 析 ) X = Type=0,1,2,3,4 に 対 しては 時 間 ( 最 大 値 は 計 算 終 了 時 間 以 下 であること) = Type=5,6 に 対 しては 圧 力 水 頭 = Type=7,8,9,10,11 に 対 してはダミー( 計 算 には 用 いられないので 何 でも 良 い) = Type=12,13,14,15 に 対 しては 測 定 位 置 ( 深 さ) Y = 測 定 値 ただし Type=2 and Position=0 のときは,カラム 全 体 の 平 均 体 積 含 水 率,また Type=4 and Position=0 のときは,カラム 全 体 の 溶 質 量 Position = Type=1,2,4 では 観 測 点 番 号 (soil profile で 指 定 する observation node の 番 号 ) = Type=0,3 では 境 界 のコード 番 号 ( 上 端 :1, 下 端 :2) = Type=5,6,7,8,9,10,11 では 物 質 番 号 (material number) = Type=12,13,14,15 に 対 しては 印 刷 時 間 (Print time)の 番 号 Weight: 逆 解 析 に 対 する 各 データの 重 み 通 常 は 1 とする ここでは,Type1( 観 測 点 における 圧 力 水 頭 変 化 ),Type2Position0( 最 終 時 間 の 平 均 体 積 含 水 率 )を 用 い る ここでのデータ 入 力 では, 別 途 エクセル 等 で 同 様 な 表 を 作 成 し, 入 力 値 を 貼 り 付 ける(Ctrl+V) 方 が 容 易 である *Type1,2,4 を 用 いる 場 合, 後 のグラフィカルエディターによって 先 に 観 測 点 を 設 定 しておく 必 要 がある そのため,まず Cancel を 選 択 し,データの 入 力 は 後 で 行 う 8
HYDRUS-1D guide ここまでの 設 定 に 対 して, 土 層 分 布 の 設 定 を 行 う 場 合 は,OK を 選 択 (12)グラフィカルエディターProfile Profile では, 空 間 刻 み 幅, 土 の 分 布 等 を 設 定 する ツールバーの Conditions( 基 本 的 には 上 から 下 の 項 目 の 順 序 ),あるいはアイコン( 基 本 的 には 左 から 右 の 順 序 )により 入 力 項 目 を 選 択 して 設 定 を 行 う 1. プロファイルの 離 散 化 Profile descritazation Conditions Profile Discretization あるいはアイコン 選 択 a. 節 点 数 (Number of Nodes): Edit Number of Nodes あるいは 画 面 左 の Edit nodes から 選 択 デフォル ト 設 定 値 は 101 である 必 要 に 応 じて 変 更 b. 固 定 点 の 挿 入 (Insert Fixed Point): 何 も 選 択 しないと, 表 面 と 下 端 の 2 点 が 自 動 的 に 選 択 される この 位 置 はマウスで 設 定 したい 位 置 を 選 択 する Option Grid の Height を 0.01 等 小 さくすることで,より 正 確 な 位 置 が 設 定 できる 固 定 ポイントの 削 除 は Delete Fixed Pont により 行 う ここでは, 観 測 点 を 挿 入 す る 深 さ 1,3,5,7,9cm に 固 定 点 を 設 定 する c. 密 度 (density): 各 固 定 ポイントの 上 下 の 密 度 を 与 える ここでは, 表 面 の 節 点 間 隔 を 小 さくするために, 9
最 上 部 の 固 定 点 をクリックし,0.2 を 入 力 する 2. 土 層 の 分 布 Material distribution Conditions Material Distriburion あるいはアイコン 選 択 土 層 の 分 布 を 指 定 する マウスで 該 当 する 節 点 を 選 択 し, 土 の 種 類 を 指 定 する 土 の 種 類 の 数 は,Main Geometry Information 中 の Number of Materials で 入 力 する ただし 均 一 試 料 の 場 合 では, 全 層 が 土 層 1 に 指 定 され,ここでの 変 更 は 必 要 な い 3. 初 期 分 布 Initial distribution Water Flow Boundary Conditions でした 項 目 ( 圧 力 水 頭 あるいは 体 積 含 水 率 )についての 分 布 を 指 定 する 複 数 の 節 点 を 選 択 し, 上 部 と 下 部 に 異 なる 値 を 入 力 すると(Use top value for both のチェックは 外 す),そ の 区 間 では 直 線 的 な 分 布 が 与 えられる 蒸 発 法 では, 初 期 のテンシオメータの 読 みを 初 期 圧 力 水 頭 とし 10
て 与 える ここで 設 定 した 分 布 は,View Graph current quantity で 確 認 可 能 4. マスバランスを 計 算 する 層 Subregion Distribution Conditions Subregions あるいはアイコン 選 択 Main Geometry Information 中 の Number of layers で 指 定 した 数 のマスバランス 領 域 を 設 定 する ここでは, 層 全 体 の 1 領 域 である 11
5. 観 測 点 Observation Points Conditions Observation Points あるいはアイコン 選 択 観 測 位 置 の 節 点 を 選 択 することにより 設 定 する 図 中 には 赤 四 角 のマークで 示 される ここでは 画 面 右 下 の 座 標 を 見 ながら, 位 置 1, 3, 5, 7, 9cm を 選 択 なお, 観 測 点 番 号 は 上 から 順 に 1,2,3 となるので, 前 述 (11) 観 測 点 に 用 いるデータ の Position にはこの 観 測 点 番 号 を 入 力 (13) 土 プロファイルのサマリー Soil Profile Summary グラフィカルエディターProfile を 終 了 すると, 設 定 した 分 布 の 表 が 示 される このデータは Profile.dat に 保 存 されている またこの 表 を 用 いて, 手 入 力 でデータの 修 正 が 行 える z: 位 置,h: 初 期 圧 力 水 頭,Root : 根 の 吸 い 込 み,Axz: 圧 力 に 関 するスケーリングファクター,Bxz: 透 水 係 数 に 関 するスケーリングファクター,Dxz: 含 水 率 に 関 するスケーリングファクター,Mat: 土 の 種 類 (Material) 番 号 である 12
*ここまで 終 了 したらもう 一 度 Data for Inverse Solution に 戻 り,データを 入 力 する( 前 述 (11)) 以 上 の 入 力 に 誤 りがなければ,プログラムを 実 行 する(Calculation あるいはアイコンを 利 用 ) 計 算 結 果 の 多 くはグラフが 作 成 される また 結 果 は 出 力 ファイルが 作 成 される 13
計 算 結 果 ここでは, 逆 解 析 に 関 する 主 な 結 果 のみを 示 す 観 測 位 置 での 圧 力 Observation Points 実 測 値 と 計 算 値 が 同 時 に 表 示 される 推 定 した 水 分 移 動 特 性 Soil Hydraulic prperties 14
逆 解 析 によるパラメータ 推 定 結 果 とその 情 報 Inverse Solution Information 逆 解 析 の 結 果 が 示 される この 結 果 は,Fit.Out にも 出 力 される 15
例 題 2 Andisol1 ここでは, 蒸 発 法 による 熊 本 黒 ボク 土 の 水 分 移 動 特 性 の 推 定 例 を 示 す 測 定 項 目 深 さ 18.39cm 土 層 からの 蒸 発 過 程 における 圧 力 水 頭 ( 深 さ 5.5,10.5,15.5cm), 積 算 蒸 発 量 ( 蒸 発 速 度 ), 実 験 終 了 時 平 均 体 積 含 水 率, 水 分 保 持 曲 線 推 定 項 目 dual porosity モデルパラメータθ r,θ s,α,n,k s,l,w 2,α 2,n 2 以 下 は, 例 題 1 と 大 きく 異 なる 部 分 のみを 示 す (1) 土 の 水 分 移 動 モデル Soil Hydraulic Model ここでは Dual-porosity(Durner, dual van Genuchten-Mualem)を 選 択 (2) 水 分 移 動 特 性 パラメータの 逆 解 析 に 関 する 設 定 Water Flow Parameters-Inverse Solution Initial Estimate: ここでは 水 分 保 持 曲 線 の 実 測 値 から 得 られたパラメータ, 飽 和 透 水 係 数 の 実 測 値 を 初 期 値 として 設 定 Fitted: 推 定 するパラメータを 選 択 ここでは,θ r,θ s,α,n,k s,l,w 2,α 2,n 2 を 選 択 (3) 変 動 境 界 条 件 Time Variable Boundary Conditions Evap.: 測 定 した 積 算 蒸 発 量 をもとに 得 られた 蒸 発 速 度 を 入 力 ここでは, 積 算 蒸 発 量 を 多 項 式 で 近 似 し, その 微 分 値 から 求 めた 蒸 発 速 度 を 入 力 ( 滑 らかな 蒸 発 速 度 の 変 化 を 与 えた 方 が, 計 算 が 安 定 するため) hcrita: HYDRUS では, 表 面 圧 力 水 頭 が-hCritA に 達 した 段 階 で 減 率 蒸 発 となり, 与 えた 蒸 発 速 度 を 満 たさなくなる できる 限 り 蒸 発 速 度 を 満 たすため( 実 験 の 水 収 支 を 満 たすため),hCritA は 安 定 した 計 算 が 得 られる 最 大 値 を 入 力 (ここでは 10 6 ) 16
(11) 逆 解 析 に 用 いるデータ Data for Inverse Solution ここでは,Type1( 深 さ 5.5,15.5cm の 観 測 点 における 圧 力 水 頭 変 化 ),Type0Position1( 積 算 蒸 発 量 ), Type5( 水 分 保 持 曲 線 の 実 測 値 ),Type2Position0( 最 終 時 間 の 平 均 体 積 含 水 率 )を 入 力 実 験 の 水 収 支 をできる 限 り 満 たすように,Type0,2 の 重 みは w = 10 とする また, 実 測 の 水 分 保 持 曲 線 を 信 頼 し,Type5 の 重 みも w = 10 とする Type2 は w = 1 プログラムを 実 行 する(Calculation あるいはアイコンを 利 用 ) 17
計 算 結 果 観 測 位 置 での 圧 力 Observation Points 実 測 値 と 計 算 値 が 同 時 に 表 示 される 積 算 蒸 発 量 Water Flow- Boundary Fluxes and Heads 実 測 値 と 計 算 値 が 同 時 に 表 示 される 18
水 分 移 動 特 性 Soil Hydraulic Properties 実 測 値 と 推 定 値 が 同 時 に 表 示 される( 推 定 値 は-2 < logh < 4 の 範 囲 で 表 示 ) 水 収 支 情 報 Mass Balance Information W-Volume/Area がカラム 内 平 均 体 積 含 水 率 となる ため, 入 力 した 最 終 平 均 体 積 含 水 率 と 比 較 する 計 算 値 : 8.2013/18.390 = 0.446 入 力 値 : 0.461 逆 解 析 によるパラメータ 推 定 結 果 とその 情 報 Inverse Solution Information 上 記 の 実 験 値 と 計 算 結 果 の 比 較, 及 び 推 定 結 果 を 総 合 的 に 評 価 し, 非 現 実 的 であれば, 推 定 パラメータ の 初 期 値,または 水 収 支 にかける 重 みの 大 きさについて 検 討 する 必 要 がある 19