(57)【要約】      （修正有）【課題】土粒子の間隙に水、空気、非水溶性流体（NAPL）が介在する多相系地盤での間隙流体の浸透現象を精度よくシミュレートし、土壌汚染の拡がりや程度の予測や、土壌汚染の浄化対策の効果をより正確に評価・検証する不飽和土壌中の非水溶性液体の挙動解析方法を提供する。【解決手段】水−空気−NAPLの３相流の支配方程式を設定し、支配方程式に対し、空間をガラーキン有限要素法、時間を後退差分法で離散化し、修正Picard法を適用することにより、非水溶性液体の挙動を解析する。【選択図】図５
2015225494.tif 000002
【特許請求の範囲】
【請求項１】
不飽和土壌中での非水溶性液体の挙動を解析する方法であって、
下記の式（１）、式（２）のように、土粒子の間隙に水、空気、非水溶性流体が介在する多相系地盤での間隙流体の浸透現象を表すモデルを構成するとともに、
下記の式（３）、式（４）、式（５）のように水−空気−NAPLの３相流の支配方程式を設定し、
式（３）、式（４）、式（５）に対し、空間をガラーキン有限要素法、時間を後退差分法で離散化し、修正Picard法を適用することによって得られる下記の式（６）、式（７）、式（８）を用いて非水溶性液体の挙動を解析することを特徴とする不飽和土壌中の非水溶性液体の挙動解析方法。
【数１】
2015225494.tif 000050
【数２】
2015225494.tif 000051
【数３】
2015225494.tif 000052
【数４】
2015225494.tif 000053
【数５】
2015225494.tif 000054
【数６】
2015225494.tif 000055
【数７】
2015225494.tif 000056
【数８】
2015225494.tif 000057
ここで、
ｎ：地盤の間隙率、Ｓ：飽和度、ｐ：圧力、ｋｒ：比透水係数、Ｋ：飽和透過係数、γ：単位体積重量、Ｋｖ：体積圧縮係数、ｑ：流量、ρ：密度、ｂ：物体力、μｒ：非粘性係数、Ｅｊｉ：∂Ｓｉ／∂ｐｊ、
［上付き／下付きのw, a, n, l］：［それぞれ水、空気、非水溶性液体（NAPL）、液体（水＋NAPL］、［上付きのｍ、ｋ］：［それぞれ時間ステップ、繰り返しステップ］、［・］：時間微分、［−（上付き）］：要素内の平均値、［−（下付き）］：境界における値
を示す。
【請求項２】
請求項１記載の不飽和土壌中の非水溶性液体の挙動解析方法において、
土粒子の間隙中での水、空気、非水溶性流体の任意の２相と３相の移り変わりを、下記の式（９）で表される中間流体圧係数μを用い、下記の式（１０）、式（１１）、式（１２）のように０から１の変数で表し、
且つ、下記の式（１３）を満たすβｗ、βｌとμの関係をベジエ曲線により与え、
水−空気のサクションｓａｗに、μにより定まる水相と、水及び非水溶性流体の液体相のスケール変数βｗ、βｌを乗じて、下記の式（１４）に代入することにより、水相と液体相の飽和度を示す特性曲線を得るようにしたことを特徴とする不飽和土壌中の非水溶性液体の挙動解析方法。
【数９】
2015225494.tif 000058
【数１０】
2015225494.tif 000059
【数１１】
2015225494.tif 000060
【数１２】
2015225494.tif 000061
【数１３】
2015225494.tif 000062
【数１４】
2015225494.tif 000063
ここで、Ｓは飽和度、Ｓｍａｘ、Ｓｍｉｎはそれぞれ最大，最小飽和度である。ｓはサクションで、ｐを圧力とするとｓｉｊ＝ｐｉ−ｐｊである（添え字ｉ、ｊは２相のうちそれぞれ濡れ性の低い流体、高い流体を示す）。
【請求項３】
請求項２記載の不飽和土壌中の非水溶性液体の挙動解析方法において、
水相と液体相の流出入履歴を反映した状態変数Ｉｗ、Ｉｌを、水及び液体の流入過程で１、排出過程で−１にそれぞれ近づく関係を満たす下記の式（１５）で与え、
下記の式（１６）のように、水及び液体の流出入履歴を考慮したＩｗ、Ｉｌでμをμｗ＊、μｌ＊に修正するようにしたことを特徴とする不飽和土壌中の非水溶性液体の挙動解析方法。
【数１５】
2015225494.tif 000064
【数１６】
2015225494.tif 000065
ここで、ｉには水相ｗ及び液体相ｌが対応する。ζｗ、ζｌは水と液体の飽和度変化に対するＩｗ、Ｉｌの変化率を決めるパラメータである。ｍｉは流出入時のヒステリシスの影響の強さを決めるパラメータである。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、不飽和土壌中の非水溶性液体の挙動解析方法に関する。
【背景技術】
【０００２】
１９８２年に実施された環境庁による全国規模の地下水汚染調査を契機に、トリクロロエチレン、テトラクロロエチレン、１，１，１−トリクロロエタンなどの揮発性有機塩素化合物や、ガソリン、灯油、軽油、重油（油類）などによる土壌・地下水汚染が都市部を中心に全国的に存在することが確認されている。
【０００３】
一方、このような汚染土壌・地下水を効果的且つ効率的に浄化、修復するためには、汚染物質の挙動、すなわち拡散域、濃度分布などをできるだけ正確に捉えることが重要になる。特に、トリクロロエチレンやテトラクロロエチレン、ガソリン、軽油、灯油などのNAPL（非水溶性液体：Non-Aqueous Phase Liquid）は、長期にわたって地盤中に滞留し、深刻な汚染を引き起こすおそれがあるため、その挙動を正確に捉えることが重要である。
【０００４】
そして、従来、地盤中での汚染物質の挙動を予測する技術手法として、１）飽和地盤中の水の浸透現象を解析する手法（単相）、２）任意の２相の浸透現象を解析する手法（２相）、３）水−NAPL−空気の浸透現象を解析する手法（３相）が用いられている。
【０００５】
１）飽和地盤中の水の浸透現象を解析する手法（単相）は、地下水面以深など土粒子間隙が水で満たされた飽和地盤中の間隙水（地下水）の浸透（移動）現象をシミュレートする手法である。土粒子間隙が単一の流体で満たされた状態での浸透現象をDarcy則と連続式（間隙水の質量保存則）を支配方程式とし、有限要素法や差分法あるいは有限体積法により空間離散化して解くことにより、地下水の挙動、ひいては地下水に含まれた汚染物質の挙動を精度よくシミュレートすることができる。
【０００６】
２）任意の２相の浸透現象を解析する手法（２相）は、例えば、地下水面以浅など、土粒子の間隙に水と空気が介在する不飽和地盤中の間隙水及び間隙空気の浸透（移動）現象をシミュレートする手法である。水−空気の浸透現象を解析する手法としては、間隙水と間隙空気それぞれのDarcy則と質量保存則に加え、水相及び空気相の圧力と飽和度の関係を与える水分特性曲線、各相の飽和度と透水・透気係数の関係、空気相の状態方程式（理想気体を仮定）を支配方程式とし、有限要素法や差分法により解く水−空気２相連成浸透解析手法が提案、実用化されている。
なお、水相及び空気相の圧力と飽和度の関係を与える水分特性曲線としてはvan Genuchtenの式（非特許文献１参照）やBrooks-Coreyの式（非特許文献２参照）がよく用いられている。また、各相の飽和度と透水・透気係数の関係にはMualemの式（非特許文献３参照）などを用いる手法が幾つか提案されている。
また、間隙流体の組み合わせが異なる２相系の水−NAPLあるいはNAPL−空気の浸透現象は、間隙流体の界面張力の違いに考慮して水−空気の水分特性曲線をスケーリングした間隙流体の圧力と飽和度の関係を用いて、水‐空気の浸透現象とほぼ同様の手法により解析されており、水−NAPLの解析は飽和地盤中、NAPL−空気の解析は乾燥地盤中のNAPLの浸透（移動）の解析に用いられている。
【０００７】
多相の浸透現象の数値解析手法には、有限要素法や有限差分法、有限体積法が多用されているが、間隙中に複数の流体が存在する地盤の浸透現象は飽和度変化に伴う相対透過度の変化が大きく、非線形性が強い問題となるため、支配方程式の離散化手法や数値計算技術に関する検討も続けられてきた。
【０００８】
そして、１９９０年代には２相系の浸透流の解析手法の開発が進められ、Celia and Bouloutas（非特許文献４参照）は水−空気の不飽和浸透の解析について、空間離散化する支配方程式を圧力で表示して解くと質量保存則（マスバランス）が満たされず数値振動が発生すること、圧力と飽和度で表示した混合型方程式で解くと質量保存則は満足されるものの数値振動が発生することを示すとともに、混合型方程式にPicard法を適用する修正Picard法を用いれば質量保存則が満たされ、圧力の振動も起こらないことを示した。さらに、日比ら（非特許文献５参照）は水−NAPL２相系の浸透解析手法について、Fully upwind法（非特許文献６参照）などの解析手法の適用性の比較を行い、修正Picard法による解法は、水−NAPL２相系の浸透解析に対しても質量保存則が満たし、圧力の振動も起こさない合理的な手法であることを示した。
【０００９】
３）水−NAPL−空気の浸透現象を解析する手法（３相）は、NAPLが地表付近で漏洩し、水や空気と混ざり合うことなく地盤中に浸透し、水相、空気相と独立した流れを形成するような場合において、土粒子の間隙に水と空気、NAPLが介在する多相地盤中の水、NAPL及び空気の多相浸透（移動）現象をシミュレートする手法である。
【００１０】
そして、この水−NAPL−空気の浸透現象を解析する手法では、水相、空気相とNAPL相それぞれのDarcy則と質量保存則に加え、３相の間隙流体の圧力と飽和度の関係、及び各相の飽和度と透水・透NAPL（油）・透気係数（あるいは相対透過度）の関係、並びに空気相の状態方程式（理想気体の圧縮特性）を支配方程式として有限要素法などにより解く必要がある。
【００１１】
水−NAPL−空気３相系地盤では、水−NAPL界面、NAPL−空気界面のみが出現し、水−空気界面は出現しないとするLeverettの仮定（非特許文献７参照）を用い、Lenhard and Parkerによる水−NAPL２相系及びNAPL−空気２相系のモデルを組み合わせた式（非特許文献８参照）が用いられている。
また、Stone（非特許文献９参照）やParker and Lenhard（非特許文献１０、非特許文献１１、非特許文献１２参照）は水−NAPL−空気３相系地盤の各相の飽和度と相対透過度の関係式を提案している。
【００１２】
水−NAPL−空気３相系地盤の浸透解析においても、２相系と同様に、圧力を未知数とした支配方程式を離散化するとマスバランスが維持されない（質量保存則が満たされない）という課題が指摘されたが、Oostrom et al. （非特許文献１３参照）はFully upwind型の有限差分法による解析手法により、日比・藤縄（非特許文献１４参照）は修正Picard法を３相系に適用することにより、それぞれマスバランスを維持した解析を実現している。
【先行技術文献】
【非特許文献】
【００１３】
【非特許文献１】van Genuchten, M.Th.(1980) A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc.Am.J.,44,pp.892-898.
【非特許文献２】Brooks, R.J and Corey, A.T.(1964) Hydraulic properties of porous media, Hydrol.Pap.3, Colo. State Univ., Fort Collins.
【非特許文献３】Mualem,Y. (1976) A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res.,12,pp.513-522
【非特許文献４】Celia, M. and Bouloutas, E. T.（1990）A general mass-conservative numerical solution for the unsaturated flow equation, Water Resources Research,26(7),p.1483-1496.
【非特許文献５】日比義彦，藤縄克之，藤原幸彦（2002）多孔体中における水−NAPL２相流の数値モデルの比較，土木学会論文集，No.720/VII-25，pp.27-38．
【非特許文献６】Thomson, N. R., Graham, D. N., and Farquhar, G. J.（1992）One-dimensional immiscible displacement experiments, J. Contamination Hydrology, 10,p.197-223.
【非特許文献７】Leverett, M.C. (1941) Capillary behavior in porous solids, Trans. Am. Inst. Min. Metal.Eng.,142,pp.152-169.
【非特許文献８】Lenhard, R.J. and Parker, J.C. (1987) A model for hysteretic constitutive relations governing multiphase flow 2. Permeability- saturation relations, Water Resour. Res., 23(12),pp.2197-2206.
【非特許文献９】Stone, H.I. (1973) Estimation of three-phase relative permeability and residual oil data,J. Can. Pet. Tech. 12(4),pp.53-61.
【非特許文献１０】Parker, J.C. and Lenhard, R.J. (1990) Determining three-phase permeability-saturation-pressure relations from two-phase system measurements, J. Pet. Sci. Eng. 4,pp.57-65,1990.
【非特許文献１１】R. J. Lenhard and Parker, J. C. : Water Resour. Res. Vol.24(3), pp.373-380,1988.
【非特許文献１２】Parker, J.C., R.J. Lenhard, and T. Kuppusamy : Water Resour. Res. Vol.23(4), pp.618-624,1987.
【非特許文献１３】Oostrom, M., Hofstee, C., Lenhard, R.J. and Wietsma, T.W.（2003）Flow behavior and residual saturation formation of liquid carbon tetrachloride in unsaturated heterogeneous porous media, Journal of Contaminant Hydrology, 64(1-2),p.93-112.
【非特許文献１４】日比、藤縄 (2005) 多孔体中における空気−水−NAPL３相流の数値モデルの比較，土木学会論文集，79781-94．
【発明の概要】
【発明が解決しようとする課題】
【００１４】
しかしながら、上記従来の水−NAPL−空気３相系の浸透解析では、各相の質量保存則とDarcy則及び圧力−飽和度（ｐ−Ｓ）関係、飽和度−相対透過度（Ｓ−ｋ）関係から得られる支配方程式を有限要素法や有限体積法などにより空間離散化するが、圧力−飽和度関係を与える既往のモデルではNAPL（非水溶性流体）の圧力が比較的小さい状態、すなわちNAPLの飽和度が低い条件下における３相の圧力−飽和度関係を正しく捉えることができないという問題があった。
【００１５】
そして、従来の水−NAPL−空気３相系の解析手法では、NAPLの圧力が比較的低い条件下でNAPLの飽和度を極端に過小評価するため、NAPL相の相対透過度も過小評価することになる。このため、NAPLが地中に侵入する際、あるいはNAPLを地盤中から排出（浄化など）する際の圧力−飽和度関係及び飽和度−固有透過度関係を正しく捉えることができないという問題があった。
【００１６】
言い換えれば、上記従来の水−NAPL−空気３相連成浸透解析手法においては、NAPLが浸透する際に飽和度が低い状況で圧力−飽和度関係を適切に表せるモデルがなく、また、NAPLの圧力が比較的低い条件でNAPLの侵入・排出時のヒステリシスを考慮することができない。このため、不飽和地盤中にNAPLが侵入する際、あるいはNAPLで汚染された土壌を浄化する際の水−NAPL−空気３相浸透現象を精度よくシミュレートすることができないという問題があった。
【００１７】
本発明は、上記事情に鑑み、土粒子の間隙に水、空気、非水溶性流体（NAPL）が介在する多相系地盤での間隙流体の浸透現象を精度よくシミュレートすることができ、土壌汚染の拡がりや程度の予測や、土壌汚染の浄化対策の効果をより正確に評価・検証することを可能にする不飽和土壌中の非水溶性液体の挙動解析方法を提供することを目的とする。
【課題を解決するための手段】
【００１８】
上記の目的を達するために、この発明は以下の手段を提供している。
【００１９】
本発明の不飽和土壌中の非水溶性液体の挙動解析方法は、不飽和土壌中での非水溶性液体の挙動を解析する方法であって、下記の式（１）、式（２）のように、土粒子の間隙に水、空気、非水溶性流体が介在する多相系地盤での間隙流体の浸透現象を表すモデルを構成するとともに、下記の式（３）、式（４）、式（５）のように水−空気−NAPLの３相流の支配方程式を設定し、式（３）、式（４）、式（５）に対し、空間をガラーキン有限要素法、時間を後退差分法で離散化し、修正Picard法を適用することによって得られる下記の式（６）、式（７）、式（８）を用いて非水溶性液体の挙動を解析することを特徴とする。
【００２０】
【数１】
2015225494.tif 000003
【数２】
2015225494.tif 000004
【数３】
2015225494.tif 000005
【数４】
2015225494.tif 000006
【数５】
2015225494.tif 000007
【数６】
2015225494.tif 000008
【数７】
2015225494.tif 000009
【数８】
2015225494.tif 000010
【００２１】
ここで、
ｎ：地盤の間隙率、Ｓ：飽和度、ｐ：圧力、ｋｒ：比透水係数、Ｋ：飽和透過係数、γ：単位体積重量、Ｋｖ：体積圧縮係数、ｑ：流量、ρ：密度、ｂ：物体力、μｒ：非粘性係数、Ｅｊｉ：∂Ｓｉ／∂ｐｊ、
［上付き／下付きのw, a, n, l］：［それぞれ水、空気、非水溶性液体（NAPL）、液体（水＋NAPL］、［上付きのｍ、ｋ］：［それぞれ時間ステップ、繰り返しステップ］、［・］：時間微分、［−（上付き）］：要素内の平均値、［−（下付き）］：境界における値
を示す。
【００２２】
また、本発明の不飽和土壌中の非水溶性液体の挙動解析方法においては、土粒子の間隙中での水、空気、非水溶性流体の任意の２相と３相の移り変わりを、下記の式（９）で表される中間流体圧係数μを用い、下記の式（１０）、式（１１）、式（１２）のように０から１の変数で表し、且つ、下記の式（１３）を満足するβｗ、βｌとμの関係をベジエ曲線により与え、水−空気のサクションｓａｗに、μにより定まる水相と、水及び非水溶性流体の液体相のスケール変数βｗ、βｌを乗じて、下記の式（１４）のvan Genuchtenの式に代入することにより、水相と液体相の飽和度を示す特性曲線を得ることが望ましい。
【００２３】
【数９】
2015225494.tif 000011
【数１０】
2015225494.tif 000012
【数１１】
2015225494.tif 000013
【数１２】
2015225494.tif 000014
【数１３】
2015225494.tif 000015
【数１４】
2015225494.tif 000016
【００２４】
ここで、Ｓは飽和度、Ｓｍａｘ、Ｓｍｉｎはそれぞれ最大，最小飽和度である。ｓはサクションで、ｐを圧力とするとｓｉｊ＝ｐｉ−ｐｊである（添え字ｉ、ｊは２相のうちそれぞれ濡れ性の低い流体、高い流体を示す）。
【００２５】
さらに、本発明の不飽和土壌中の非水溶性液体の挙動解析方法においては、水相と液体相の流出入履歴を反映した状態変数Ｉｗ、Ｉｌを、水及び液体の流入過程で１、排出過程で−１にそれぞれ近づく関係を満たす下記の式（１５）で与え、下記の式（１６）のように、水及び液体の流出入履歴を考慮したＩｗ、Ｉｌでμをμｗ＊、μｌ＊に修正するようにしたことを特徴とする。
【００２６】
【数１５】
2015225494.tif 000017
【数１６】
2015225494.tif 000018
【００２７】
ここで、ｉには水相ｗ及び液体相ｌが対応する。ζｗ、ζｌは水と液体の飽和度変化に対するＩｗ、Ｉｌの変化率を決めるパラメータである。ｍｉは流出入時のヒステリシスの影響の強さを決めるパラメータである。
【発明の効果】
【００２８】
本発明の不飽和土壌中の非水溶性液体の挙動解析方法においては、土粒子の間隙に水、空気、非水溶性流体（NAPL）が介在する多相系地盤での間隙流体の浸透現象を精度よくシミュレートすることが可能になる。
【図面の簡単な説明】
【００２９】
【図１】LNAPL、DNAPLの浸透現象を示す模式図である。
【図２】本発明の第１実施形態に係る不飽和土壌中の非水溶性液体の挙動解析方法によるシミュレーション例で用いた解析モデルを示す図である。
【図３】本発明の第１実施形態に係る不飽和土壌中の非水溶性液体の挙動解析方法を用いて行ったLNAPLの浸透解析結果を示す図である。
【図４】本発明の第１実施形態に係る不飽和土壌中の非水溶性液体の挙動解析方法を用いて行ったDNAPLの浸透解析結果を示す図である。
【図５】水−非水溶性液体（NAPL）−空気の３相流体が土粒子間隙中において濡れ性が高い水、NAPL、空気の順に土粒子の接触点近傍を満たすというLeverettの仮定を示す図である。
【図６】本発明の第２実施形態に係る不飽和土壌中の非水溶性液体の挙動解析方法の説明に用いた本発明のモデルのパラメータβと中間圧係数μの関係を示す図である。
【図７】従来のＬＰモデル（図７（ａ））と本発明のモデル（図７（ｂ））のパラメータβと中間圧係数μの関係を示す図である。
【図８】棚橋らによる油層厚実験（カラム実験）の模式図と圧力仮定を示す図である。
【図９】本発明の第２実施形態に係る不飽和土壌中の非水溶性液体の挙動解析方法の優位性を確認するために行ったシミュレーション結果であり、棚橋らによる油層厚実験（カラム実験）の結果（図９（ａ））と、従来のＬＰモデルの結果（図９（ｂ））と、本発明のモデルの結果（図９（ｃ））を比較した図である。
【図１０】本発明の第３実施形態に係る不飽和土壌中の非水溶性液体の挙動解析方法の説明に用いた中間流体圧係数μと、履歴を考慮した中間流体圧係数μ＊と、飽和度Ｓの関係を示す図である。
【図１１】ＬＰモデルと、第２実施形態のモデルと、第３実施形態のモデルのシミュレーション結果を比較した図である。
【発明を実施するための形態】
【００３０】
以下、図１から図４を参照し、本発明の第１実施形態に係る不飽和土壌中の非水溶性液体の挙動解析方法について説明する。ここで、本実施形態は、土粒子の間隙に水、空気、NAPL（非水溶性流体）が介在する多相系地盤での間隙流体の浸透現象を精度よくシミュレートすることができ、土壌汚染の拡がりや程度の予測や、土壌汚染の浄化対策の効果をより正確に評価・検証することを可能にする不飽和土壌中の非水溶性液体の挙動解析方法に関するものである。
【００３１】
はじめに、NAPLによる土壌汚染現象やその後の浄化問題を解くには、地盤の浸透特性や流体保持特性を適切に考慮した数値解析手法が必要である。特に、浅い不飽和地盤のNAPL汚染では、土中に混在する水−NAPL−空気３相の圧力と飽和度の関係を精度よく反映した特性曲線とそれに基づく多相流解析手法が有効である。
従来の３相系特性曲線としては、水−NAPL及びNAPL−空気の２相特性曲線を用いたLenhard and Parkerのモデル（非特許文献１０、１１、１２：以下、この従来のモデルを「ＬＰモデル」という）が簡便で、初期値問題の解析にもよく用いられてきた。同モデルは、次の式（１７）、式（１８）に示すように、水飽和度は水圧とNAPL圧、液体飽和度（水飽和度＋NAPL飽和度）はNAPL圧と空気圧を変数として与えている。
【００３２】
ここで、本実施形態において、使用する記号は以下のように定義する。
ｎ：間隙率、Ｓ：飽和度、ｐ：圧力、ｋｒ：比透水係数、Ｋ：飽和透過係数、γ：単位体積重量、Ｋｖ：体積圧縮係数、ｑ：流量、ρ：密度、ｂ：物体力、μｒ：非粘性係数、Ｅｊｉ：∂Ｓｉ／∂ｐｊ、
［上付き／下付きのw, a, n, l］：［それぞれ水、空気、NAPL、液体（水＋NAPL］、［上付きのｍ、ｋ］：［それぞれ時間ステップ、繰り返しステップ］、［・］：時間微分、［−（上付き）］：要素内の平均値、［−（下付き）］：境界における値
【００３３】
【数１７】
2015225494.tif 000019
【数１８】
2015225494.tif 000020
【００３４】
この式（１７）、式（１８）に示すように、ＬＰモデルでは、水飽和度及び液体飽和度がそれぞれ空気圧及び水圧と無関係に与えられており、各相の相互作用を考慮していない。すなわち、水−空気の相互作用を無視して、水−NAPL、及びNAPL−空気の２相特性曲線で簡便に表現している。このため、NAPLの飽和度が低く、水−空気の２相系に近い状態にある地盤中の飽和度分布（３相地盤の圧力−飽和度関係）を適切に表現することができず、NAPLによる不飽和地盤の汚染メカニズムやその浄化過程を正確に解くことができない（後述する図９参照）。
なお、式（１７）、式（１８）における４つの偏微分係数Ｅは、例えば後述の式（３２）で算出することができる。
【００３５】
これに対し、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法では、次の式（１９）、式（２０）を用いてモデルを構成し、各相の圧力の相互作用を考慮した３相系特性曲線によって、土粒子の間隙に水、空気、非水溶性流体（NAPL）が介在する多相系地盤での間隙流体の浸透現象を、精度よくシミュレートできるようにする。
【００３６】
【数１９】
2015225494.tif 000021
【数２０】
2015225494.tif 000022
【００３７】
具体的に、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法では、地盤内の多相流れの支配方程式として、質量保存則と間隙流体の飽和度変化を記述する特性曲線（式（１９）、式（２０））を考えるとともに、各相の浸透現象がそれぞれDarcy則に従うものとする。
ここで、飽和度と比透過係数の関係はMualemの式（非特許文献３）を用いる。そして、水、NAPLは非圧縮流体、空気は理想気体と仮定すると、水−空気−NAPL３相流の支配方程式は次の式（２１）、式（２２）、式（２３）のようになる。
【００３８】
【数２１】
2015225494.tif 000023
【数２２】
2015225494.tif 000024
【数２３】
2015225494.tif 000025
【００３９】
そして、式（２１）、式（２２）、式（２３）に対し、空間をガラーキン有限要素法、時間を後退差分法で離散化して、修正Picard法を適用することにより、次の式（２４）、式（２５）、式（２６）を得る。
【００４０】
【数２４】
2015225494.tif 000026
【数２５】
2015225494.tif 000027
【数２６】
2015225494.tif 000028
なお、式（２４）、式（２５）、式（２６）の｛  ｝は節点における値を示すベクトルである。
【００４１】
次に、Ｓｗｍ＋１，ｋ＋１、Ｓｎｍ＋１，ｋ＋１、Ｓａｍ＋１，ｋ＋１をｐｗ、ｐｎ、ｐａでテーラー展開して整理すると、次の式（２７）、式（２８）、式（２９）が得られる。
【００４２】
【数２７】
2015225494.tif 000029
【数２８】
2015225494.tif 000030
【数２９】
2015225494.tif 000031
【００４３】
ここで、δｍ＋１，ｋ＋１＝ｐｍ＋１，ｋ＋１−ｐｍ＋１，ｋであり、δを未知変数としてδ≒（ニアリーイコール）０ となるまで繰り返し計算を行う。
【００４４】
なお、式（２７）〜式（２９）中の行列は下記の式（３０）で与えられる。
【数３０】
2015225494.tif 000032
【００４５】
次に、上記のように修正Picard法を用いて定式化した本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法による計算例（シミュレーション結果）について説明する。
【００４６】
ここでは、登坂らの２次元解析（登坂ら：地下水学会誌，38(3)，pp.167-180，1996．）を参考にし、地下水面を持つ不飽和地盤中に地表からNAPLが流入する現象のシミュレーションを行っている。
具体的に、解析モデルは、図１、図２に示すような２次元断面であり、地下水面以浅には２箇所、３相いずれも浸透しない不透過層を設けている。
【００４７】
また、主な物性値を表１に示す。
【００４８】
【表１】
2015225494.tif 000033
【００４９】
そして、本シミュレーションでは、表１に示す水より軽いLNAPLと、水より重いＤNAPLの２種類について計算を行っている。また、NAPLの侵入は流量制御で与え、その流量を単位奥行当たり２．０×１０−４ｍ２／ｓｅｃで一定とした。さらに、境界条件として、側方境界では地下水面を一定に保ち、NAPLについては侵入箇所を除いて流量０の位置を境界とした。
【００５０】
図３、図４はシミュレーション結果であり、NAPLの飽和度分布の変化を示している。
図３及び図４に示すように、NAPLは地表面から流入すると不透過層に滞留しつつ徐々に下方へ浸透し、且つ、水より軽いLNAPLに比べ、水より重いDNAPLが早く流下する結果が得られた。また、NAPLの比重による差異は地下水面に到達してからより明確になり、LNAPLは地下水面付近に滞留して水平方向に広がるのに対し、DNAPLは地下水面上に滞留せず、水との密度差により地下水面下に侵入して底部に滞留する結果が得られた。さらに、LNAPLは地下水面上に滞留する際に、地下水面をやや押し下げる様子も確認された。
すなわち、図１に示すLNAPL、DNAPLの浸透現象をよく捉えていた。
【００５１】
このように、従来、地表付近でDNAPL（水より重い非水溶性流体、トリクロロエチレンなど）やLNAPL（水より軽い非水溶性流体、重油・軽油・灯油など）が漏洩した際の汚染流体の拡がりは、ボーリング調査などにより離散的にしか把握することができなかったが、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法によれば、既往の調査法と組み合わせることにより、精度よく地盤中での汚染の拡がりを特定することができる。
【００５２】
なお、本実施形態では２次元の結果を計算例として示したが、本発明の解析方法は基本的に１次元でも３次元でも解くことが可能である。
【００５３】
したがって、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法においては、式（２１）、式（２２）、式（２３）のように水−空気−NAPLの３相流の支配方程式を設定し、空間をガラーキン有限要素法、時間を後退差分法で離散化して、修正Picard法を適用し、式（２４）、式（２５）、式（２６）のように定式化したことで、土粒子の間隙に水、空気、非水溶性流体（NAPL）が介在する多相系地盤での間隙流体の浸透現象を精度よくシミュレートすることが可能になる。
【００５４】
また、NAPLで汚染された地盤を浄化処理する際に、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法を用いてその浄化効果等を予め検証することにより、より合理的な浄化対策法の選定が可能になる。
【００５５】
次に、図５から図９を参照し、本発明の第２実施形態に係る不飽和土壌中の非水溶性液体の挙動解析方法について説明する。
ここで、本実施形態は、第１実施形態と同様、土粒子の間隙に水、空気、非水溶性流体（NAPL）が介在する多相系地盤での間隙流体の浸透現象を精度よくシミュレートすることができ、土壌汚染の拡がりや程度の予測や、土壌汚染の浄化対策の効果をより正確に評価・検証することを可能にする不飽和土壌中の非水溶性液体の挙動解析方法に関し、特に、従来の解析手法では捉えることができなかった水−NAPL−空気３相系の圧力−飽和度関係から、水−空気２相系の圧力-飽和度関係まで連続的に精度よく捉えることを可能にする不飽和土壌中の非水溶性液体の挙動解析方法に関するものである。
なお、本実施形態では、第１実施形態と同様の事項に対し、その詳細な説明を省略する。
【００５６】
はじめに、van Genuchtenの式を任意の２相系流体に拡張させた特性曲線（非特許文献１２）は、次の式（３１）で与えられる。
【００５７】
【数３１】
2015225494.tif 000034
ここで、Ｓは飽和度、Ｓｍａｘ、Ｓｍｉｎはそれぞれ最大，最小飽和度である。ｓはサクションで、ｐを圧力とするとｓｉｊ＝ｐｉ−ｐｊである。また、αは材料パラメータ、βは２相の間隙流体の組み合わせで決まるパラメータでβａｗ＝１である。添え字ｉ、ｊは２相のうちそれぞれ濡れ性の低い流体、高い流体を示す。
【００５８】
そして、Lenhard and Parker（非特許文献１１）は、図５に示すように、間隙中の３相流体について、濡れ性が高い水、NAPL、空気の順に土粒子の接触点近傍を満たすというLeverettの仮定に則って、それぞれの境界に注目し、次の式（３２）で表されるように、水飽和度を水−NAPLの２相特性曲線、水とNAPLの飽和度の和（液体飽和度）をNAPL−空気の２相特性曲線で与えている。
【００５９】
【数３２】
2015225494.tif 000035
なお、第１実施形態と同様、添え字のｗ、ａ、ｎ、ｌはそれぞれ水、空気、NAPL、液体（水＋NAPL）を表し、下付き文字は流体の組み合わせ、上付き文字は飽和度を与える流体を示す。
【００６０】
そして、この式（３２）から分かるように、ＬＰモデルは２相特性曲線を用いる点で簡便であるが、Ｓｗａｎｗ＞ Ｓｌａｎｗとなり、NAPL飽和度が負になるときは、便宜的に場合分けし、NAPLが消失してＳｗａｎｗ＝ Ｓｌａｎｗとみなし、ＳｗａｎｗとＳｌａｎｗを水−空気２相特性曲線で与えるようにしている。
【００６１】
また、ＬＰモデルでは、式（３１）に着目すると、NAPL圧が水圧と一致するとき、NAPLは水に支配されて消失し、水−空気の２相系になる。さらに、NAPL圧がｐｗ＜ｐｎ＜ｐａを満たすとき、水−空気−NAPLの３相系になる。さらに、NAPL圧が空気圧に一致するとき、NAPLは空気を消失させて、水−NAPLの２相系になる。
【００６２】
一方、本願の発明者らは、このような任意の２相と３相の移り変わりを、新たに次の式（３３）で定義し、次の式（３４）、式（３５）、式（３６）のように中間流体圧係数μによって０から１の変数で表すことができるようにした。すなわち、μを用いて、前述した間隙の状態を示すと次の式（３４）、式（３５）、式（３６）となる。
【００６３】
【数３３】
2015225494.tif 000036
【００６４】
【数３４】
2015225494.tif 000037
【数３５】
2015225494.tif 000038
【数３６】
2015225494.tif 000039
【００６５】
ここで、式（３４）〜（３６）の関係を得るためには、関数β（μ）は式（３７）を満足する必要がある。
【００６６】
【数３７】
2015225494.tif 000040
【００６７】
本実施形態では、式（３７）を満たすβｗ、βｌはそれぞれ３次のベジエ曲線で与える。そして、ベジエ曲線を描くには３つの制御点ｒ１、ｒ２、ｒ３が必要になるが、βｗ、βｌそれぞれに、図６に示す（μ，β）平面上の制御点ｒ１ｗ、ｒ２ｗ、ｒ３ｗ、及びｒ１ｌ、ｒ２ｌ、ｒ３ｌを用いる。ベジエ曲線は両端でｒ１とｒ３を通り、その間の制御点ｒ２で曲率が決まる。また、曲線はｒ１とｒ３での接線がそれぞれｒ１−ｒ２及びｒ２−ｒ３と平行になり、滑らかに接続する。ベジエ曲線上の任意の点ｒは次の式（３８）のような線形和で与えられる。
【００６８】
【数３８】
2015225494.tif 000041
【００６９】
ｋ１、ｋ２、ｋ３は媒介変数ｔ（０≦ｔ≦１）で次の式（３９）のように与えられる。
【００７０】
【数３９】
2015225494.tif 000042
【００７１】
なお、ｔ＝０でｒはｒ１、ｔ＝１でｒはｒ３となる。ここではベジエ曲線の曲率を調整できるように、曲率のパラメータａ（−２≦ａ≦２）を導入し、式（３９）の代わりに次の式（４０）を用いる。
【００７２】
【数４０】
2015225494.tif 000043
【００７３】
また、ｋ０からｋ３は次の式（４１）で与える。
【００７４】
【数４１】
2015225494.tif 000044
【００７５】
そして、任意の２相と３相の移り変わりを新たな式（３３）で定義するようにし、式（３３）〜式（４１）で得られるβｌ、βｗとｓａｗを次の式（４２）に代入することにより、式（１９）、式（２０）の６つの偏微分係数Ｅを与えることができる。
【００７６】
【数４２】
2015225494.tif 000045
【００７７】
なお、μを通してＬＰモデルの特徴を考察すると、まず、式（３３）を式（３２）に代入して変形すると次の式（４３）を得る。
さらに、式（４３）は式（４４）のように簡潔に表現できる。
【００７８】
【数４３】
2015225494.tif 000046
【数４４】
2015225494.tif 000047
【００７９】
そして、式（４４）から、この図７（ａ）に示す通り、ＬＰモデルは水−空気の２相系と水−空気−NAPLの３相系において別々の式で規定しているため、両者の移行する点において必ずしも連続とならない。また、μ＝０のときではなくβｌ＞βｗのときに水−空気の２相系となる。
【００８０】
これに対し、本実施形態のように、式（３７）を満足するβｗ、βｌとμの関係を、ベジエ曲線を用いて図７（ｂ）のように与えることで、３相系と水−空気２相系の移行を滑らかに表すことができる。また、図７（ｂ）に示すように、μが０に近づくときβｗとβｌはともに１に近づき、水飽和度と液体飽和度が接近してNAPL飽和度が０に向かう。
したがって、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法を用いることにより、土粒子の間隙に水、空気、非水溶性流体（NAPL）が介在する多相系地盤での間隙流体の浸透現象を精度よくシミュレートすることが可能になる。
【００８１】
次に、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法を用いて行ったシミュレーションを、棚橋らの３相系のカラム実験（棚橋秀行,佐藤健,小西純一：土木学会論文集C Vol.62 No.2, pp. 320-334, 2006.）と比較し、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法の優位性を確認した結果について説明する。
【００８２】
棚橋らによるカラム実験では、模型地盤内に設けた井戸の水面上に油（NAPL）を溜めた際の長時間経過後の地盤内の３相の飽和度分布を観察している。地中の圧力分布は長時間経過後の定常状態を想定し、密度を勾配として図８に示すような線形分布を仮定している。
【００８３】
図９は、カラム実験の実験値と、ＬＰモデル、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法による解析値を比較した結果を示している。
【００８４】
まず、図９（ａ）に示すように、カラム実験では油層厚が小さい場合でも油面上に油が拡がることが確認されている。このことから、汚染サイトの原位置の観測井では数ｃｍオーダーの油層厚の検出が推測され、油面以浅での土壌の汚染が懸念させる。
【００８５】
これに対し、図９（ｂ）に示すように、ＬＰモデルでは油面以浅の油の存在が算出されず、油飽和度が過小評価され、実際の土壌汚染に対して危険側の解を与えることが確認された。他方、図９（ｃ）に示すように、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法では油面上の地中の油の存在を算出しており、カラム実験の実験値をよく捉えている。
【００８６】
したがって、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法においては、中間流体圧係数μに基づいたモデルを用いることで、棚橋らのカラム実験を精度よく再現できることが実証された。
【００８７】
よって、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法によれば、土粒子の間隙に水、空気、非水溶性流体（NAPL）が介在する多相系地盤での間隙流体の浸透現象を、より精度よくシミュレートすることができ、土壌汚染の拡がりや程度の予測や、土壌汚染の浄化対策の効果を、より正確に評価・検証することが可能になる。
【００８８】
次に、図１０、図１１を参照し、本発明の第３実施形態に係る不飽和土壌中の非水溶性液体の挙動解析方法について説明する。なお、本実施形態では、第１、第２実施形態と同様の事項について、その詳細な説明を省略する。
【００８９】
ここで、第２実施形態では、中間流体圧係数μに基づいたモデルを用いることで、棚橋らのカラム実験をよく再現できることを実証した。一方で、第２実施形態では、流体間のサクションｓと各相の飽和度Ｓに一義的な関係を与えたことで、NAPLの流入や排出に伴うヒステリシス現象を考慮していない。
このため、本実施形態では、水相及び液体相（水相とNAPL相の和）の流出入を規定する状態変数を新たに導入し、ヒステリシス現象も考慮できるようにした。
【００９０】
具体的に、第２実施形態に示した通り、中間流体圧係数μ（＝（ｐｎ−ｐｗ）／（ｐａ−ｐｗ））は３相流体の圧力の相対的関係を表す変数で、０から１への変化が水−空気２相系から水−NAPL−空気３相系、水−NAPL２相系への移り変わりに対応する。そして、第２実施形態では、水−空気のサクションｓａｗに、μにより定まる水相と液体相のスケール変数βｗ、βｌを乗じて、van Genuchtenの式に代入することにより、水相と液体相の飽和度を得る特性曲線モデルにすることを説明した。すなわち、この第２実施形態のモデルでは３相の圧力の変数であるｓａｗとμに対して飽和度を唯一に決めるので、圧力履歴に依存した応答が反映されていない。
【００９１】
これに対し、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法では、水相と液体相の流出入履歴を反映した状態変数Ｉｗ、Ｉｌも考慮してβｗ、βｌを与えることでヒステリシス現象を表現する。
【００９２】
まず、Ｉｗ、Ｉｌは水及び液体の流入過程で１、排出過程で−１にそれぞれ近づくとする。このような関係を満たす発展則として次の式（４５）を与える。
【００９３】
【数４５】
2015225494.tif 000048
ここで、ｉには水相ｗ及び液体相ｌが対応する。ζｗ、ζｌは水と液体の飽和度変化に対するＩｗ、Ｉｌの変化率を決めるパラメータである。
【００９４】
そして、履歴依存性は、水及び液体の流出入履歴を考慮したＩｗ、Ｉｌでμをμｗ＊、μｌ＊に修正することにより、次の式（４６）のように表現する。
【００９５】
【数４６】
2015225494.tif 000049
【００９６】
ここで、ｍｉは流出入時のヒステリシスの影響の強さを決めるパラメータである。つまり、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法では３相の圧力に対して唯一に決まるμに代わり、状態変数Ｉｗ、Ｉｌで補正した水及び液体の係数μｗ＊、μｌ＊を用いることで履歴依存性を考慮している。
【００９７】
図１０は、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法による計算過程を示している。この図１０の左上、左下の図は、μに対して履歴依存性を反映した状態変数Ｉｗ、Ｉｌによりヒステリシスを考慮したμｗ＊、μｌ＊を得る過程である。Ｉｉ＝０のときμｉ＝μｉ＊（４５度の破線）となり、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法のモデルはヒステリシスを考慮しない第１実施形態のモデルと同じ解を与える（図１０の黒点線）。
また、図１０の中段左図は、μｗ＊、μｌ＊と水相及び液体相のスケール変数βｗ、βｌの関係を示している。μｗ＊、μｌ＊には履歴依存性の影響が既に反映されている。
最後に、図１０の中段右図のようにβｗ、βｌをｓａｗに乗じてvan Genuchten式に代入し、水と液体の飽和度Ｓｗ、Ｓｌを求める。
【００９８】
これにより、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法では、第２実施形態の長所をそのままに、状態変数Ｉｉを媒介変数としてヒステリシス現象を確実に反映させることができる。
【００９９】
次に、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法を用いたシミュレーションについて説明する。
【０１００】
ここでは、地下水面の上に薄くNAPL（Ｌｎ＝０．０６７ｍ）を張った井戸のカラム実験を初期状態として、油面を井戸の上端（Ｌｎ＝１．６ｍ）まで上昇させた後、再び元の位置に戻すことで、NAPLを流入出させた。
【０１０１】
そして、この本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法によるモデル（第３実施形態のモデル）と、Lenhard & Parkerのモデル（ＬＰモデル）及びヒステリシスを考慮しないモデル（第２実施形態のモデル）でそれぞれシミュレーションを行い、その結果を比較した。なお、圧力分布は長時間経過後の定常状態を想定して静水圧、静油圧分布としている。
【０１０２】
図１１は、上段がLenhard & Parkerのモデル（ＬＰモデル）、中段がヒステリシスを考慮しないモデル（第２実施形態のモデル）、下段が本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法のモデルによる解析結果の一例を示している。また、左列が流入過程、右列が流出過程のＬｎ＝０．８ｍでの地中の飽和度分布を示している。
【０１０３】
そして、図１１に示すように、Lenhard & Parkerのモデル（上段）は、同じ油層厚に対して、ヒステリシスによる水飽和度の分布の違いを算出するが、NAPL飽和度が極端に過小評価されている。
【０１０４】
また、ヒステリシスを考慮しないモデル（中段）は、既往実験でも報告された薄いNAPL層厚下でのNAPL飽和度の分布をよく捉えているが、履歴特性を考慮していないため、流入・流出過程で全く同じ解になってしまう。
【０１０５】
これに対し、実施形態の不飽和土壌中の非水溶性液体の挙動解析方法によるモデル（下段）は、NAPLの流入過程よりもNAPLの流出過程でNAPL飽和度が高い様子を算出しており、NAPL侵入時にはNAPL飽和度が上昇しにくく、排出時にはNAPL飽和度が下降しにくいという土壌中の流体の流出入時のヒステリシス現象をよく捉えている。
【０１０６】
したがって、本実施形態の不飽和土壌中の非水溶性液体の挙動解析方法においては、ヒステリシス現象をよく捉え、土粒子の間隙に水、空気、非水溶性流体（NAPL）が介在する多相系地盤での間隙流体の浸透現象を、さらに精度よくシミュレートすることができる。よって、土壌汚染の拡がりや程度の予測や、土壌汚染の浄化対策の効果を、さらに正確に評価・検証することが可能になる。
【０１０７】
以上、本発明に係る不飽和土壌中の非水溶性液体の挙動解析方法の実施形態について説明したが、本発明は上記の第１、第２、第３実施形態に限定されるものではなく、その趣旨を逸脱しない範囲で適宜変更可能である。
【図６】
2015225494.tif 000066
【図７】
2015225494.tif 000067
【図１０】
2015225494.tif 000068
【図１】
2015225494.tif 000069
【図２】
2015225494.tif 000070
【図３】
2015225494.tif 000071
【図４】
2015225494.tif 000072
【図５】
2015225494.tif 000073
【図８】
2015225494.tif 000074
【図９】
2015225494.tif 000075
【図１１】
2015225494.tif 000076