透水係数の入替え処

PHREEQC を Python で動かしてみました。

IPhreeqc と PhreeqPy を試してみましたが、残念ながら動いたのは前者のみ。

で、動かして気づいたのですが、これ、２,３の species に着目するだけではダメですね。輸送部分を FiPy に任せるつもりだったのですが、できた solution をそのまま次のステップに持っていく必要があります。例えば、calcite の溶解を扱うにしても pH、CaCo3、Ca2+、CO2 のみではダメ。HCO3-、 CaHCO3+、CO3-2はもちろん、その他の elements, species のみならず redox や温度などの条件も次のステップに持っていかなければなりません。先の文献のように、簡単ではありませんでした。

そうなるとこの種の問題は PHAST などで輸送まで計算しておいて、ステップ間にて透水係数を変化させた方が簡単そうですね。FiPy での物理量の入れ替えまで確認できていたのですが、見通しが甘かったようです。

ま、バカなりに手を動かして、見通しを得られただけでも良しとしましょう。

OGS#IPhreeqc

注文していた図書が来ました。

OpenGeoSys Tutorial: Computational Hydrology III: OGS#IPhreeqc Coupled Reactive Transport Modeling (SpringerBriefs in Earth System Sciences)

PHREEQC‎ 関連の図書としては久しぶりの発刊です。ここ数年のモデリング例や表流水との連成例でも掲載されているか？と期待して注文しました。が、残念ながら外れ。

簡単な2例でのチュートリアルをベースに、機能の説明と使用法が記載されています。

OGS#IPhreeqc の操作性に言及すれば、今までの Reactive Transport Modelling に対応したソフト（PHAST など）より、退化しています。ジオメトリ作成は GINA。確認しましたが8月の試行時点より Ver. は上がっていません。残りはエディタで手打ちです。

https://phreeqc.blogspot.jp/2017/08/ginaogs.html

機能は同等、インターフェースは劣る、といった評価になるため、表流水との連成例などが出てこない限り使用しないでしょう。

このような RTM に関するソフトの利用は、日本では全く流行っていません。化学物質による汚染の拡散を数値解析で解く場合でも、国家プロジェクトを除けば、3次元移流分散レベルで止まっています。吸着計算の代わりに遅延係数などフィッティングパラメータを入れて「答えがあえば良いでしょう」状態から脱却していません。ま、簡単でビジネスとして成り立つ価格になりますし、化学計算にもフィッテイングパラメータは入るので、流行らないだけなのかもしれませんが。

さらにその先、熱-水-応力-化学の連成は、日本では程遠い世界のようです。

Asの含有形態

土のAs含有形態について考えています。

熱力学データを用いた反応計算でAsの含有形態を推定するには、土と水（間隙水）の両者のデータが必要です（定量的にも、定性的にも）。ところが、実務で実施される含有量試験は湿潤状態の試料であり、間隙水に含まれるAsも一緒に測定してしまいます。その影響がどの程度かは現場によって異なるでしょうし、区分して測らない限りは良く分からないデータが得られることになります。つまり、計算には使えません。
特に吸着の場合はひと手間必要です。純水で土を洗うと間隙水の影響は抜けますが、吸着しているAsも一緒に洗い流してしまいます。水を使用してはダメということです（海外の参考書では詳細が述べられています）。

一方、国交省「建設工事における自然由来重金属等含有岩石・土壌への対応マニュアル（暫定版）平成22年3月」には資料集に「3.岩石・土壌に含まれる重金属等の起源を識別する試験法」が掲載されています。
3段階抽出法が紹介されていますが、評価の仕方は書かれていません。溶出量は土1kg当たりに換算するのでしょうか？文献を読みなさいということでしょうか？
3段階のうち、2段階（溶出量試験・含有量試験）は通常行いますので、あとはXRFで分析すれば良いだけです。先に述べたように水と土の分離があいまいである以上、このような大雑把な試験法や比較でも実務的には価値があると判断されているのでしょう。
5段階も紹介されていますが、こちらも十分対応できるレベルだと思います。コチラのほうが説明しやすいですね。

実務で含有形態を知りたい場合では、まずマニュアルの5段階抽出法を実施し、重要案件では化学計算を行うのが良いのかもしれません。
含有形態が分かり、それを計算で再現ができたなら、コンセプトモデルは出来上がりです。あとは搬出先の条件に合わせて、水-反応連成計算を行い、敷地境界での濃度を推定すれば良いですね。こちらは3次元まで対応可能でしょう。MATLAB+COMSOL+IPhreeqc がスマートですが、HYDRUS, PHREEQC, PHASTの組み合わせでも対応できますね。
http://phreeqc.blogspot.com/2011/05/matlab-comsol-iphreeqc.html

調査計画、特に採水方法や保管方法など、また岩石・鉱物・土の判定方法や欲しい精度を、事前に良く考えておかないと問題は解けませんね。以前にも書きましたが、こういった環境問題は地質屋さんの力量に大きく依存していると思います。
http://phreeqc.blogspot.com/2011/05/limestone4.html

実務で含有形態の考慮やコンセプトモデルの作成は、まだまだ実施されていません。法的調査に必要ないからでしょう。業界の意図も含まれている気がしますね。

今年の目標

今年も半分終わりました。

短期目標として、今年は以下の点をこなそうと思っていました。進捗は50%程度と言ったところでしょうか？

1. Dtransuの高速化　OpenMP, GPGPU
2. 飽和・不飽和輸送と反応の連成
3. 変形解析

Dtransuは半分クリアーです。あとはGPGPUだけですね。
ただし、Dtransuのコード読みが終わっていないので、CUDAコードと併せて読み込む必要がありそうです。目標にはしていませんが、使用している技術者として当然のことなので、やっておく必要があります。でも、難しいんですよね。

輸送と反応は1次元のみHP1でクリアーです。3次元はMATLAB+COMSOL+IPhreeqc と連携コードで可能ですが、ソフトがありません。COMSOLのセミナーには出席する予定ですが、ペイできなければ買ってくれないでしょう。
しかし、お金を出せばできることは分かりました。この辺りでクリアーにしましょう。

変形解析は、まだまだですね。
分かったと思ったら、次に分からないことがドンと出てくる段階です。山の裾野のほうでしょう。早く連続体から不連続体、粒子法へ行きたいのですが、先は長そうです。学ぶべきことは、まだまだたくさんあります。

残り半年、焦らずに時間をかけてやりましょう。

MATLAB + COMSOL + IPhreeqc

概要は理解できました。COMSOLで移流分散、IPhreeqcで反応、それらの結果の橋渡しをMATLABで行うというものでした。

最初に MATLAB が 初期化を行ったあと、輸送計算を COMSOL で行います。 次にMATLABが comsol2phreeqc.m で結果を変換し、COMモジュール（IPhreeqc）を呼び出して反応計算を行います。反応の結果を phreeqc2comsol.m で変換し、COMSOLで次のステップの計算を行うといった流れです。というわけで、計算方法は non‐iterative です。HP1と同じですね。

早速、見積もりを取りました。両方合わせて300万を超えました。ムリです。買えません。
大企業か、行動力のある小企業しか無理でしょう。

しかし、お金を出せば3次元飽和・不飽和輸送＋化学反応連成計算ができることが分かりました。
ペイできるだけの業務が取れる見込みもありませんし、、当分は HP1 ＋ PHAST で我慢しましょう。

IPhreeqc

PHREEQC 2.18と同時に公開された IPhreeqc、いいですね。
http://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc/

C、C++、Fortranやなどで利用できるよう、モジュールとして配布されています。

COMモジュールをインストールし、EXCEL2007のVBAで試しましたが、動きますね。びっくりです。ただ、EXCELではメリットないですけどね。MATLAB　＋　COMSOL　Multiphysics　で3Dの飽和・不飽和移流分散と組み合わせた例があるようですので、文献を取り寄せてみましょう。

Fortran で利用できるため、理屈では Dtransu に組み込むことができます。アイデアはありますが、組んでくれるプログラマー、いないでしょうか。

しかし、USGSはすごいですね。こんなパワフルなコードを、多くの言語で使えるようにモジュール化して、しかも無償で公開しているのですから。お金の流れはどういう仕組みなんでしょう。
日本も見習ってほしいものです。