データ同化手法＋地下水モデルによる地下水涵養量の時系列推定手法の確立－石川県手取川扇状地を事例として－
石川県立大学生物資源環境学部 藤原洋一・高瀬恵次・長野峻介
１. はじめに石川県手取川扇状地は北陸を代表する扇状地の一つである。豊富な地下水を有し、私たちに多くの恩恵をもたらしている。地下水の用途としては、例えば、良質な水を用いた産業利用（繊維産業や酒造業、先端産業など）や農業利用が挙げられる。また、上水道や冬季の消雪用としても利用され、私たちの生活に深く結びついている。さらに、地下水は通年で水温変化が小さく、トミヨといった希少生態系の生息環境を形成している。このように、手取川扇状地の地下水は地域水資源として非常に重要な役割を担っている。

図 2 土砂崩壊地の衛星画像（Google、2016）

図 1 手取川流域

図 3 (a)日降水量、(b)濁度（日平均値）の推移

この手取川扇状地上流域において、2014 年 10 月に小規模な土砂崩壊が発生した。そして、2015 年 5月に融雪に伴い、崩壊土砂量 130 万 m3 という大規模な斜面崩壊に拡大した（図 1、図 2）。この土砂が手取川に流入したことで、浮遊土砂濃度 1,000mg/L を超える高濃度濁水が長期間継続した（図 3）。柳井ら（2017）は、機械学習であるランダムフォレストを用いて地下水変動を分析し、地下水位の異常低下時期と濁水発生時期が合致することを指摘している。田中ら（2018）は、実測調査として水田浸透量調査
および河川伏流量調査行い、濁水発生前後の結果を比較したところ、濁水発生後に水田浸透量と河川伏流量が著しく減少したことを確認している。また、吉岡ら（2018）は、扇状地地下水および地表水の酸素・水素安定同位体比の観測を行い、左岸域における河川水の涵養の寄与割合の低下を示唆している。

また、右岸域においては、河川水と田面水の涵養の減少が共に生じたことにより、大きな地下水低下が生じたと推察している。以上のように、高濃度濁水が扇状地地下水環境に及ぼす影響について様々な研究が行われてきた。しかし、これらの研究は扇状地の地下水環境を調べるためには必要不可欠な調査であるが、多大な労力とコストを要し、高頻度に行うことが困難である。そのため、水田および河川からの涵養量の減少がどのくらい続くのか、地下水位の回復がいつどのように生じるのか、浸透量が回復するとなれば何が要因となって回復するのかなど、未だに不明な点が多い。
そこで、本研究では濁水発生前後の水田・河川からの涵養量を時系列推定し、地下水がいつ・どのように回復するのかといった回復プロセスを明らかにすることを試みた。具体的には、手取川扇状地を対象とした地下水モデルを構築し、モデルにデータ同化手法（粒子フィルタ）適用することで水田・河川からの涵養量の時系列変化を１日毎に（日単位で）推定した。
２. 方法
２. １

地下水モデル

手取川扇状地の地下水位変動を再現するために、地表帯・中間帯・地下水帯領域の 3 領域で構成される集中定数型 3 段タンクモデルを構築した（図 4）。手取川上流域からの流量やその他の背後地流量、平野部における降雨などの入力要素に加え、河川からの農業用取水、地下水からの各種揚水取水などの様々な水循環を考慮している。そして、第 3 段目タンク貯留深の変動が扇状地の平均的な地下水位変動を表現している。実測地下水位は、手取川扇状地内 8 ヶ所の観測井の位置からティーセン法によってそれぞれの面積重みを求め、8 ヶ所の加重平均値を実測地下水位としている。
河川からの伏流に関しては、モデル内で計算している。高瀬ら（2019）は、河川流量と伏流量の関係性を考慮したサブモデルを含む集中定数型水循環モデルにより、瀬戸内海沿岸平野部の地下水位を上手く再現している。そこで、手取川扇状地を対象とする本モデルにおいても、類似したサブモデルを用いることにした。サブモデルの概要としては、手取川流量に比例した量が伏流するようになっており、河川流量がある一定値に達すると、伏流量が一定になるような上限値が設けられている（図 5）。水田からの浸透量は、地表帯領域のサブ領域である「水田」において計算される。
本モデルには、各領域タンクからの流出や浸透に関わるパラメータ、河川からの伏流に関わるパラメータなどがある。これらの物理的に決定できないパラメータについては、大域的探索手法の一つである
進化戦略（藤原ら、2003）を用いて同定した。同定期間は濁水発生前の 2013~2014 年とした。地下水位の再現性を確認したところ、濁水発生前の 2013~2014 年における年平均相対誤差は、2 年間の平均で 1.1％となり、本モデルで精度よく地下水位を再現できることが示された。しかし、濁水発生後 2015~2019 年の年平均相対誤差は、5 年間の平均で 29.4 %となり、実測地下水位に比べて過大推定していた（図 6）。
この過大推定は、濁水発生後に水田からの浸透量と河川伏流量が減少したことが原因だと考えられた。
２.２ 粒子フィルタ
モデルを用いて濁水前の地下水変動を上手く再現できていることから、水田および河川涵養量に関わるパラメータの変化を逐次推定できれば、涵養量の時系列変化を把握できると考えられる。そこで、地下水モデルにデータ同化手法を適用し、濁水発生後の地下水位の低下とその後の変動を再現できるパラメータ（影響係数）を逐次推定することにした。本研究では、影響係数の推定にデータ同化手法として粒子フィルタを適用した。

図 4 地下水モデルの概要

図 5 手取川からの伏流量計算

図 6 地下水位の再現結果
粒子フィルタは、ベイズ理論を応用したモンテカルロ法的状態修正法であり、観測値𝑦𝑡 が得られた後の状態量𝑥𝑡 の事後分布𝑝( 𝑥𝑡 | 𝑦𝑡 )を多数のサンプル（粒子）で近似的に表現することが大きな特徴である。非線形・非ガウス型のあらゆる状態空間モデルに適用可能、計算機システムへの実装が極めて容易であることから、気象学・海洋学・経済学など幅広い分野で用いられている。近年では、水工学の分野においても実時間洪水予測（工藤ら、2009）などに粒子フィルタが用いられている。
粒子フィルタで推定する影響係数は、0～1 の範囲で値をとる。影響係数を濁水前の水田・河川涵養量に乗じることで、濁水後の前の水田・河川涵養量を求めることができる。影響係数が 1 に近い場合は、濁水前とほぼ変化がないことを意味し、0 に近い場合には、濁水前に比べ大きく減少したことを意味する。本研究では以下の手順で適用した。
① 水田・河川涵養量の 2 個のパラメータを逐次推定の対象とする。これに 4 個の貯留高（水田、その他、中間帯、地下水帯：図 4）を加えたものを 1 つの粒子とする。それぞれのタンク貯留高を粒子に加えた理由は、前日の貯留高を保存するためであり、逐次推定の対象とはしていない。そして、各パラメータの値に正規乱数を加え、1000 個の初期粒子を生成する。
② 各粒子で地下水位を計算（予測）する。また、各粒子には正規乱数でシステム誤差（各パラメータ範囲の 4%）を加える。

③ 観測地下水位（日単位）を用いて各粒子の尤度を計算する。観測誤差は観測地下水位の 15%とした。さらに、各粒子の尤度を正規化して、各粒子の重みを求める。
④ 重みに応じて、1000 個の粒子のリサンプリングを行う。t＝t＋1 とし、②に戻る。

図 7 粒子フィルタを適用した地下水位再現結果

図 8 影響係数の時系列変化
３. 結果および考察
粒子フィルタを地下水モデルに適用した結果を図 7 に示す。2015 年冬季、2016 年冬季から 2017 年春
季にかけての急激な変動に対して差異は見られるものの、観測地下水位と計算地下水位はほぼ一致していた。 2015~2019 年の年平均相対誤差は 12.4 %、5.2 %、0.9 %、0.8 %、-0.8 %であり、濁水発生後 5 年間の平均は 3.7 %となった。粒子フィルタを適用しない場合の相対誤差が 29.4 %であったことを踏まえると、地下水位を精度よく再現できており、粒子フィルタによって状態（パラメータ）の修正が上手く行われたと考えられる。

図 9 水田浸透量の時系列変化

図 10

河川浸透量の時系列変化

水田および河川の影響係数の時系列変化を図 8 に示す。濁水が発生した 2015 年 5 月以降に、まず水田の影響係数が急激に低下し、遅れて河川の影響係数が急激に低下したことが分かる。2016 年には増加と減少を繰り返し、2 つの影響係数は小さい状態が続いている。2017 年夏季から影響係数は増加し始め、12 月には濁水発生前の状態まで回復していることが分かる。しかし、2018 年に入り再び急激に低下し、夏季から冬季にかけて回復している。2019 年も 2018 年と同様の変化をしている。
扇状地全体に対する水田浸透量の時系列変化の結果を図 9 に示す。農業用水を取水する灌漑期間において浸透量が増加し、非灌漑期間に浸透量が減少していることから、営農活動を反映できていると考えられる。しかし、各年の灌漑期間の平均浸透量を求めたところ、2013 年 8.5 mm/d、2014 年9.1 mm/d、2015 年 3.8 mm/d、2016 年 1.4 mm/d、2017 年 6.0 mm/d、2018 年 8.5 mm/d、2019 年 5.2 mm/dであり、実測調査の結果 2014 年 3.5mm/d、2016 年 3.2 mm/d と比べると過大推定となっていた。
河川伏流量の時系列変化の結果を図 10 に示す。濁水発生前 2013 年 1 月の実測値が 4.19 mm/d、濁水発生後 2018 年 6 月が 5.02 mm/d であったのに対し、図から計算値の最大値が約 1.7 mm/d で打ち止めになっていることが分かる。これは、河川流量がある一定値に達すると、伏流量が一定になるようなサブモデルの上限設定が原因だと考えられる（図 5）。また、河川伏流量を過小推定したことが、水田浸透量の過大推定に繋がったと考えられる。
４. おわりに本研究では、手取川扇状地の地下水位変動を再現するために、地表帯・中間帯・地下水帯領域の 3領域で構成される集中定数型 3 段タンクモデルを構築した。構築したモデルに、データ同化手法として粒子フィルタを適用することで、水田・河川涵養量の時系列変化を日単位で把握することを試みた。以下に本研究により明らかになった事項を列挙する。
（1）モデルに粒子フィルタを適用することで、濁水発生後の観測地下水位と計算地下水位の平均相対誤差は 29.4%から 3.7%となり、精度よく地下水位を再現することができた。したがって、粒子フィルタによって状態（パラメータ）の修正を上手く行うことができたと考えられる。
（2）水田浸透量と河川伏流量に関わるパラメータとして影響係数を逐次推定することで、濁水発生後の涵養量の時系列変化を日単位で推定することができた。しかし、現段階では水田浸透量が過大推定、河川伏流量が過小推定していることが分かった。
（3）今後モデルの改善を行い実測値との整合性が向上した際には、断片的にしか評価できない実測調査とは異なり、連続したデータから濁水の影響評価を行うことができると考えられる。最終的には、地下水位の回復プロセスについて解明できると考えられる。
謝辞：本研究は、北陸地域づくり協会「北陸地域の活性化｣に関する研究助成事業」の支援を受けた。
５. 引用文献藤原洋一・田中丸治哉、畑武志、多田明夫（2003）：進化戦略による流出モデル定数の最適同定、農業土木論文集、227、119-129⚫

Google（2016）：（参照 2016.5.18）：Google Earth、（オンライン）、入手先＜https://earth.google.com/web/@36.88764455、138.42445732、-39409.17061751a、170536.57816239d、35y、62.85905408h、76.84301697t、359.9998r?utm_source=earth7&utm_campaign=vine&hl=ja ＞⚫

工藤亮治・近森秀高・永井明博（2009）：粒子フィルタを用いた河川流域における実時間洪 水予測、
農業農村工学会論文集、259、17-25⚫

高瀬恵次・徳増実（2019）：集中定数型水循環モデルによる沿岸平野部地下水の水収支解析、農業農村工学論文集、309（87-2）、I_211-I_218⚫

田中健二・瀬川学・藤原洋一・高瀬恵次・丸山利輔・長野俊介（2018）：高濃度濁水が扇状地の水田浸透量・河川伏流量に及ぼす影響－手取川流域の大規模土砂崩壊を事例として－、農業農村工学会
論文集、86（1）、I_47-I_54⚫

柳井清治・岡崎正規・高瀬恵次・瀧本裕士・一恩英二・百瀬年彦・藤原洋一・北村俊平・長野峻介・本多裕司（2017）：白山源流で発生した地すべりによる濁水が手取川の流域環境に及ぼす影響とその対策（2）、平成 28 年度石川県立大学年報、9-34⚫

吉岡有美・伊藤真帆・中村公人・瀧本裕士・土原健雄（2018）：酸素・水素安定同位体比からみた手
取川扇状地の河川水－地下水の交流現象と地下水涵養源、地下水学会誌、60（2）、205-221