JP5154918B2 - 地質環境の長期変動が地下水流動に与える影響を評価するための解析方法 - Google Patents

Description

本発明は、地質環境の長期変動が地下水流動に与える影響を評価するための解析方法に関する。

地下深部の地下水流動場は地形や地質、気侯の変化など自然現象の長期的な変遷の影響を受けており、現在の流動場を理解し、将来の流動場を予測するには、これらの影響を考慮する必要がある。特に高レベル放射性廃棄物の地層処分の安全評価では、１０万年以上の長期評価が必要となることから、地形・地質構造及び地層の透水パラメータなど、地質環境の長期変動を考慮した地下水流動への影響評価が必要になる。

しかしながら、従来の地下水流動解析では、多くの場合、地形・地層構造などのモデル形状や透水パラメータの解析過程での変化を解析に反映させることはできなかった。

ここで、従来の地下水流動解析方法を例示すれば、特許文献１には、解析対象流域の水循環モデルをコンピュータ上に構築し、この水循環モデルにより流域の水循環状態を再現し、汚染物質が排出される地点あるいは地域から、流下地域までの流動経路を地表水・地下水系の連続した流線として再現し、陸水汚染の危険度を解析的に評価する方法が記載されている。

また特許文献２には、任意点から地下水の流動経路および任意時間経過後の到達点に関する期待値や分散を示す確率的情報を求めるための数値解析方法が記載されており、この数値解析方法では、地下水流動特性の統計量を用いて確率有限要素法により浸透流解析し、各有限要素の流速ベクトル毎の期待値と分散を算出し、算出した流速ベクトル毎の期待値と分散に従うようなランダムな要素流速ベクトルをモンテ・カルロ法を適用して設定し、設定したランダムな要素流速ベクトルから粒子追跡法によって流動経路および到達点の試行結果の１つを算出し、同様な試行により算出した多数の結果から流動経路や到達点に関する確率的情報を求めている。

以上の特許文献１及び２は、いずれも地質環境の長期変動を考慮した地下水流動への影響評価が可能なものではない。

これ以外では、非特許文献３に、沿岸域の塩淡境界挙動の将来予測に際して、海水準の挙動、地盤の堆積構造の変化、地形変化など塩淡境界挙動に影響を与える要因を考慮することができる塩淡境界挙動モデルの構築を目的とした研究が記載されている。
しかしながら、この非特許文献３でも、地層の追加や削除により地質構造の変化を表現するにとどまり、地層の形状変化や断層の発生・進展、透水パラメータの時間的変化の考慮のような課題が残されており、地形・地層構造などのモデル形状や透水パラメータの解析過程での変化を十分に解析に反映させることはできなかった。
特開２００７−０７２７５３号公報 特開２００６−２４２８９３号公報 與田敏昭、徳楠充宏、新田昭、西垣誠『１５．沿岸域の地形変化を考慮した塩淡境界挙動モデルの構築』、２００３年１０月２３〜２４日、日本地下水学会２００３年秋季講演会講演要旨、第５８〜６３頁

以上のような現状を鑑みて本発明の課題は、地形・地質構造及び地層の透水パラメータなど、地質環境の長期変動を考慮した地下水流動の解析方法を提供することである。

上記課題を解決するために、本発明では下記（１）〜（５）の手段が提供される。

（１）地質環境の長期変動が地下水流動に与える影響を評価するための解析方法であって、地下水流動の解析対象流域において、現在から過去所定時まで遡る解析対象期間の複数時点における地質構造断面図を現在の地質構造断面図から推定し、各地質構造断面図を所定の座標軸に基づいて数値化して地質構造数値データを作成し、前記解析対象期間を任意の個数のタイムステップに分割し、各タイムステップの所定時点に対応した地質構造断面図を前記地質構造数値データから線形補間によりそれぞれ作成し、当該線形補間により作成された各地質構造断面図から水平方向及び垂直方向のメッシュからなるＦＥＭモデルをそれぞれ作成し、最も古いタイムステップのＦＥＭモデルに所定の境界条件を設定して定常地下水流動解析を実施し、当該定常地下水流動解析の結果を初期値として最も古いタイムステップのＦＥＭモデルにより非定常の飽和・不飽和地下水流動解析を実施し、当該地下水流動解析の結果を次に古いタイムステップのＦＥＭモデルに初期値として入力し、当該タイムステップにおいて地形・地質構造は変化しないと仮定すると共に、所定の境界条件を入力して非定常の飽和・不飽和地下水流動解析を実施し、これ以降、前のタイムステップにおける地下水流動解析の結果を、次のタイムステップのＦＥＭモデルに初期値として入力し、当該タイムステップにおいて地形・地質構造は変化しないと仮定すると共に、所定の境界条件を入力して非定常の飽和・不飽和地下水流動解析を実施する工程をタイムステップＮまで繰り返すことを特徴とする地下水流動の解析方法。

（２）地形及び地質構造の特徴部がそれぞれ対応するように、前記線形補間により作成された各タイムステップの地質構造断面図を複数の区分（「ＦＥＭメッシュデータ」、又は「ジオモデル」とも言う）にそれぞれ分割し、所定のタイムステップのＦＥＭモデルにおける節点に初期値を入力するに際して、当該節点が含まれる区分に対応する区分を、前のタイムステップの地質構造断面図において求め、当該対応区分の対応位置に在るＦＥＭメッシュの要素から補間計算によって初期値を求めることを特徴とする前記（１）に記載の地下水流動の解析方法。

（３）前記境界条件が、海水準及び涵養量である前記（１）に記載の地下水流動の解析方法。

（４）前記解析対象の各地質構造ごとに、透水に関するパラメータとして時間を変数とする関数をそれぞれ設定し、当該パラメータを各タイムステップのＦＥＭモデルに入力することを特徴とする前記（１）に記載の地下水流動の解析方法。

（５）前記透水に関するパラメータが、透水係数、間隙率及び比貯留係数である前記（４）に記載の地下水流動の解析方法。

本発明の解析方法では、現在から過去所定時まで遡る解析対象期間を任意の個数のタイムステップに分割し、タイムステップごとに地質構造断面のＦＥＭモデルをそれぞれ作成し、前のタイムステップのＦＥＭモデルにより得られた解析結果を、順次、次のタイムステップのＦＥＭモデルに初期値として入力し、タイムステップごとに所定の境界条件や透水に関するパラメータを入力して非定常の飽和・不飽和地下水流動解析を実施する工程をタイムステップＮまで繰り返すことにより地下水流動場を求めるものである。
したがって、従来の地下水流動の解析では、地形・地層構造等のモデル形状・構造や透水に関するパラメータの長期変動の影響を解析過程で考慮することは不可能であったが、本発明の解析方法では、タイムステップごとに該当する地質環境の長期変動を反映させると共に、各タイムステップ間で解析結果を引き継ぎながら解析を行うものであるため、地質環境の長期変動が地下水流動に及ぼす影響を評価することができる。

以下、本発明の実施の形態を説明するが、本発明はこれに限定されるものではない。
本実施形態では、実際の地下水流動場、すなわち、北海道幌延地域の地下水流動場に本発明を適用して１５０万年前から現在までの地質環境の長期的変遷が地下水流動に与える影響を評価した例について説明する。

最初に、現在から過去所定時まで遡る解析対象期間の複数時点における地質構造断面図を、現在の地質構造断面図を基に、既存データを参考にして復元する。図１は、このように復元した北海道幌延地域の１５０万年前から現在までの地質構造断面図であり、現在の地質環境から想定される地下水流動方向に沿って、ほぼ東西方向で作成したものである。海岸線は各年代における海水準から推定し、１５０万年前は全域が海水面下である。隆起・沈降速度は１０００年間に数１０ｃｍオーダである。

地質構造は時間経過と共に変化するものであるため、ＦＥＭモデルの形状も変化させる必要がある。そのため、本実施形態では、図１のような地質構造断面図を所定の座標軸に基づいて数値化して地質構造数値データを作成し、図１に示した時点以外の地質構造断面図は、地質構造数値データから線形補間によりそれぞれ作成した。そして、各地質構造断面図から水平方向及び垂直方向にメッシュを切ってＦＥＭモデルを作成した。図２は、現在の地質構造断面図から作成したＦＥＭモデルの例であり、図３は、これに各地層、断層等の名称を付したものである。

次に、地質構造変遷の解析上の扱い方について説明する。
本実施形態では、１５０万年前から現在まで解析対象時間をＮ個（Ｎ＝１０）のタイムステップに分割した。各タイムステップ（期間）ではＦＥＭモデルの形状は一定として、各タイムステップ内において海水準や涵養量が変化する非定常解析を実施した。ＦＥＭモデルの形状は、各タイムステップの中央時点に該当するモデルを用いた。各タイムステップの最終時間ステップの解析結果（圧力及び濃度）を、次タイムステップの初期値として入力し連続的な計算を実施した。
なお、このタイムステップの分割数は任意に増やすことが可能であり、これにより地質構造の変化をより詳細に表現することもできる。

各タイムステップのＦＥＭモデル間で初期値を引き継ぐ際には、地質構造の変遷を考慮した値の受け渡しを実施する。受け渡しの一般的なフロー（Ａ１）〜（Ａ６）の概要を図４（ａ）（ｂ）を参照しながら説明する。本実施形態の解析では、地層構造は主に鉛直方向に変形し、水平方向の移動はほぼ無視しうると考えられるため、地質構造の水平方向（ｘ）を固定し、鉛直方向の変化だけで各ステップ間の初期値を引き継ぐこととした。
（Ａ１）初期値の引き継ぎは、新タイムステップのＦＥＭモデルの全節点に関して実施する。
（Ａ２）次に、解析用のＦＥＭモデルとは別に、図４（ａ）に示したように、地質構造断面における特徴部、例えば、断層の交点、背斜軸及び向斜軸などの地層境界形状を離散化したＦＥＭメッシュデータ(以下、ジオモデルという)に準じたデータセットを作成する。
（Ａ３）対象とする新ステップのＦＥＭモデルのメッシュの節点が新ステップのジオモデルにおいて含まれる要素を同定し、さらに、その要素内での局所座標(要素内の相対位置)を算定する。
（Ａ４）上記で同定された要素と算定された局所座標を前タイムステップのジオモデルに適用し、前タイムステップのＦＥＭモデルにおける座標値を算定する。
（Ａ５）算定された座標値に対応する前タイムステップの最終時間ステップの圧力及び濃度を算定する。
（Ａ６）算定された圧力及び濃度が新タイムステップのＦＥＭモデルの初期値となる。

次に、解析のフローについて説明する。解析は（Ｂ１）〜（Ｂ７）の順に実行する。
（Ｂ１）最も古い第一のタイムステップの地質構造断面図及びＦＥＭモデルを作成する。
（Ｂ２）第一のタイムステップのＦＥＭモデルで海水位及び涵養量を設定し、定常解析を実施する。
（Ｂ３）第一のタイムステップのＦＥＭモデルで非定常解析を実施する。
（Ｂ４）第二のタイムステップの地質構造断面図及びＦＥＭモデルを作成する。
（Ｂ５）第一のタイムステップの非定常解析による最終時点の圧力水頭及び濃度から、第二のタイムステップのＦＥＭモデルに入力する初期値を算定する。
（Ｂ６）算定された初期値を第二のタイムステップのＦＥＭモデルに入力して非定常解析を実施する。
（Ｂ７）これ以降、前のタイムステップのＦＥＭモデルで非定常解析して求めた最終時点の圧力水頭及び濃度から、次のタイムステップで用いる初期値を算定し、この初期値を次のタイムステップのＦＥＭモデルに入力して非定常解析する工程を最終タイムステップ（すなわち、第十のタイムステップ）まで繰り返す。

次に、解析手法について説明する。解析は、地下水流動解析、塩分濃度に関する物質移行解析を実施する。地下水流動解析では塩分濃度よりも、地形の地下水流動への影響が大きいと考え、密度の影響は考慮せずに二次元地下水流動解析を実施し、この流速を利用して塩分濃度の非定常変化を移流・分散を考慮した物質移行解析により算定した。
解析は現状での地下水流動状況と想定される条件で地下水に関する定常解析を実施し、この結果を初期条件として海水準変動、涵養量変動の影響を境界条件として非定常の地下水流動解析を実施するとともに、並行して物質移行解析により塩分濃度の挙動に関しても推定した。
解析には、既存のコンピュータソフトウェアであるＦＥＭ地下水流動解析用プログラムDtransu・3D・ELを基本コードとし、部分的に本実施形態の解析目的に合うように透水係数設定、境界条件制御などについてソースコードを修正して使用した。
Dtransuは、岡山大学西垣誠教授、三菱マテリアル株式会社、ダイヤコンサルタント株式会社の三者共同で開発されたプログラムである。Dtransuは、密度勾配を考慮した飽和・不飽和浸透流及び移流・分散問題を対象とした解析コードで、移流・分散解析においてオイラリアン・ラグランジアン手法を用い、高ペクレ数から低ペクレ数の問題に対して安定した解析が可能である。

次に、解析条件について説明する。
地下水流動の境界条件は、図５に示したように、解析領域側部及び底部を不透水とし、海底部は海水から算出される水位を固定水位として設定する。また陸部は涵養量とともに自由浸出条件を設定する。海水準変動は既存の調査結果から図６（ａ）に示した海水準変動曲線を設定する。涵養量の変化は、現在の値を３６５ｍｍ／年として、海水準と涵養量が線形の相関性を有すると仮定して推定する。既存の調査結果から凍土厚の変遷を推定し、凍土が存在する場合は涵養量が０となるように、図６（ｂ）の曲線を設定する。
物質移行解析の境界条件を図７に示す。モデル側部、底部は不透水境界でゼロフラックスとする。陸部表面では濃度Ｃ＝０の天水が涵養される。海部では濃度Ｃ＝１.０として海水であることを示す。厳密には流出箇所では濃度を有した地下水の湧出も想定されるが、設定した初期濃度分布から濃度を有した湧出の可能性が小さいこと、計算の収束性を確保する観点から、この境界条件とした。他の境界部に関しては特に境界条件を設置していない。海進、海退に伴い沿岸海底部は図８に示したように境界条件が変化する。
現在の海水準と涵養量を与えた条件で地下水の定常解析を実施し、この結果を地下水解
析の初期条件とした。また物質移行の初期条件は全領域濃度１.０を与えた。分散長は、縦分散長αＬ＝１００ｍ、横分散長αＴ＝２０ｍと設定した。

本実施の形態では、透水に関するパラメータとして、図９に示したように透水係数、間隙率及び比貯留係数を設定した。図９中のDepth＿Yt、Depth＿Kt、Depth＿Wkは各層の深度依存性モデルを示す。これらは、透水係数ｋ(ｍ／ｓ)、深度Ｚ（ｍ）とするとそれぞれ次式で表され、上限と下限は下記の通り設定した。
Depth＿Yt：log10(k)＝-0.0034Ｚ-8.3665
＜上限：1×10^-8ｍ／ｓ，下限：1×10^-11ｍ／ｓ＞
Depth＿Kt：log10(k)＝-0.OO39Ｚ-7.5935
＜上限：1×10^-7ｍ／ｓ，下限：1×10^-11ｍ／ｓ＞
Depth＿Wk：log10(k)＝-0.0061Ｚ-5.5626
＜上限：1×10^-6ｍ／ｓ，下限：1×10^-11ｍ／ｓ＞
図９中では、透水係数、間隙率及び比貯留係数を定数、または各層の深度を変数とする関数で示したが、これら透水に関するパラメータは地層厚さの時間的変化が与えられると、時間を変数とする関数でそれぞれ設定することも可能であり、例えば、透水係数はKozenyの式である下記の式（１）、間隙率は下記の式（２）、比貯留係数は下記の式（３）のような関数を用いることができる。

ｋ＝Ｃ・ｎ³／（１−ｎ）²・・・・・・・式（１）
ｋ：透水係数
ｎ：間隙率
ｃ：定数

ｎ＝（Ｈ／Ｈ₀−１＋ｎ₀）／（Ｈ／Ｈ₀）・・・・・・・式（２）
ｎ：間隙率（ｎはＨの関数で時間の関数である）
Ｈ：地層の厚さで時間とともに変化
Ｈ₀：初期の地層厚さ
ｎ₀：初期間隙率

Ｓｓ＝γ_w・ｍ_v・・・・・・・式（３）
Ｓｓ：比貯留係数
γ_w：水の単位体積重量
ｍ_v：体積圧縮係数
ここで、ｍ_v＝（ΔＨ／Ｈ）／Δｐ
すなわち、ｍ_vは体積圧縮係数で堆積に伴う荷重の変化（Δｐ）と地層厚の変化（ΔＨ）
、地層厚（Ｈ）より与えられる。

解析は、地質構造の変遷の有無、透水係数の変遷の有無を組み合わせて、下記Case1〜3について実施した。
Case1：地形・地質構造の変化無し
Case2：地形・地質構造が変化有り
Case3：地形・地質構造が変化有り＋透水係数も変化有り
Case1は地形・地質構造に変化のない従来の解析方法によるケース、Case2は本発明を適用した地形・地質構造の変化を考慮した解析方法によるケース、Case3は地形・地質構造の変化に加え、地層の圧密過程で変化することが想定される透水係数の低下も考慮したケースである。

以上のようにCase1〜3について解析し、１５０万年経過時点（現在）で得られた結果を図１０〜図１３に示した。図１０は現在の解析領域右側４０ｋｍでの全水頭の分布図である。図１１は、図１０の分布図において矢印で示した（ａ）海域、（ｂ）沿岸、（ｃ）丘陵部における鉛直方向分布である。図１２は現在の解析領域右側４０ｋｍでの塩分濃度の分布図である。図１３は、図１２の分布図において矢印で示した（ａ）海域、（ｂ）沿岸、（ｃ）丘陵部における鉛直方向分布である。
（ａ）海域、（ｂ）沿岸では、地形・地質構造の変化を考慮したCase2、Case3で全水頭が深部で小さくなる傾向が見られ、沈降の影響が現れていると考えられる。本実施形態である地形・地質構造の変化を考慮したCase2、Case3の解析では、地形・地質の変化するタイムスチップ間で圧力水頭が引き継がれる。これに対して、沈降している次タイムスチップでは沈降に伴う低い位置水頭が加算されるため全水頭は低下する。
計算開始時点の透水係数を大きく設定したCase3では、洗い出しによる淡水化が深部まで及んでいる。Case1とCase2の比較では、Case2よりCase1で淡水化部が深部まで及んでいる。これは隆起に伴い下方から塩分の高い地層が上昇する効果である。

前記Case1〜3の解析結果を比較すると、下記のような知見が得られた。
（１）地形・地質構造の変化を考慮したCase2、Case3では、それらを考慮しないCase1よりも全水頭が小さな値を示す傾向にある。これは、１５０万年前以降の期間では、海水準が上昇傾向にあることを反映した結果である。
（２）地形・地質構造の変化に加え透水係数の変化を考慮したCase3では、塩分の希釈速度は早くなる。これは、解析の初期段階において１桁大きな透水係数を設定し、それが圧密の進行とともに小さくなると設定したことによる影響である。
（３）地形・地質構造の変化を考慮したCase2、Case3では、それらを考慮しないCase1よりも塩分の希釈が遅い。これは、隆起部において、比較的塩分濃度の高い地下水を含む地層が上昇するとともに、洗い出しによって比較的塩分濃度の低くなった地下水を含む地層が、削剥されるためである。なお、海域で沈降が生じた場合には、海水の存在により、淡水による上方からの希釈は生じない。

以上の解析結果から得られる全体的な傾向から、過去所定時からの地質環境の長期的変遷が地下水流動や塩分濃度に与える影響を評価できる可能性が確認できた。

現在の地質環境から復元した北海道幌延地域の１５０万年前から現在までの地質構造断面図である。 現在の地質構造断面図から作成したＦＥＭモデルの例である。 図２に各地層、断層等の名称を付したものである。 （ａ）（ｂ）は各タイムステップのＦＥＭモデル間で初期値を受け渡す方法の概要を示した図である。 地下水流動解析の境界条件を示した図である。 （ａ）は地下水流動解析の境界条件である海水準の変動曲線であり、（ｂ）は同じく涵養量の変動曲線である。 物質移行解析の境界条件を示した図である。 海進、海退に伴う沿岸海底部の境界条件の変化を示した図である。 透水に関するパラメータである透水係数、間隙率及び比貯留係数を示した表である。 １５０万年経過時点の解析結果であって、全水頭の分布図である。 図１０の分布図における（ａ）海域、（ｂ）沿岸、（ｃ）丘陵部の鉛直方向分布である。 １５０万年経過時点の解析結果であって、塩分濃度の分布図である。 図１２の分布図における（ａ）海域、（ｂ）沿岸、（ｃ）丘陵部の鉛直方向分布である。

Claims (5)

1. 地質環境の長期変動が地下水流動に与える影響を評価するための解析方法であって、
   地下水流動の解析対象流域において、現在から過去所定時まで遡る解析対象期間の複数時点における地質構造断面図を現在の地質構造断面図から推定し、各地質構造断面図を所定の座標軸に基づいて数値化して地質構造数値データを作成し、
   前記解析対象期間を任意の個数のタイムステップに分割し、各タイムステップの所定時点に対応した地質構造断面図を前記地質構造数値データから線形補間によりそれぞれ作成し、当該線形補間により作成された各地質構造断面図から水平方向及び垂直方向のメッシュからなるＦＥＭモデルをそれぞれ作成し、
   最も古いタイムステップのＦＥＭモデルに所定の境界条件を設定して定常地下水流動解析を実施し、当該定常地下水流動解析の結果を初期値として最も古いタイムステップのＦＥＭモデルにより非定常の飽和・不飽和地下水流動解析を実施し、
   当該非定常の地下水流動解析の結果を次に古いタイムステップのＦＥＭモデルに初期値として入力し、当該タイムステップにおいて地形・地質構造は変化しないと仮定すると共に、所定の境界条件を入力して非定常の飽和・不飽和地下水流動解析を実施し、
   これ以降、前のタイムステップにおける地下水流動解析の結果を、次のタイムステップのＦＥＭモデルに初期値として入力し、当該タイムステップにおいて地形・地質構造は変化しないと仮定すると共に、所定の境界条件を入力して非定常の飽和・不飽和地下水流動解析を実施する工程をタイムステップＮまで繰り返すことを特徴とする地下水流動の解析方法。
2. 地形及び地質構造の特徴部がそれぞれ対応するように、前記線形補間により作成された各タイムステップの地質構造断面図を複数の区分にそれぞれ分割し、
   所定のタイムステップのＦＥＭモデルにおける節点に初期値を入力するに際して、当該節点が含まれる区分に対応する区分を、前のタイムステップの地質構造断面図において求め、当該対応区分の対応位置に在るＦＥＭメッシュの要素から補間計算によって初期値を求めることを特徴とする請求項１に記載の地下水流動の解析方法。
3. 前記境界条件が、海水準及び涵養量である請求項１に記載の地下水流動の解析方法。
4. 前記解析対象の各地質構造ごとに、透水に関するパラメータとして時間を変数とする関数をそれぞれ設定し、当該パラメータを各タイムステップのＦＥＭモデルに入力することを特徴とする請求項１に記載の地下水流動の解析方法。
5. 前記透水に関するパラメータが、透水係数、間隙率及び比貯留係数である請求項４に記載の地下水流動の解析方法。