緒論

第2章では、 現地の地下水水質の形成を考察するために、 水質化学的、 地球化学的な検討 を行った。 その結果、 Ca、 Mg、 Na、 Kの陽イオン濃度の変動を解明する上で、 陽イオンの 輸送過程を把握することが重要であることが示唆された。

地下環境中では、 固相と液相聞における化学反応は主として陽イオン交換反応である。 カ リウムなどのような陽イオン交換を伴う物質が土壌に負荷された場合には、 他の陽イオンと の間に化学反応を考慮した地下水中での物質輸送モデルの開発が必要となる。 陽イオン交換 反応は瞬時に化学平衡に達する反応であるため、 モデル化にあたっては反応速度論的な取り 扱いは不要となる。

本章では、 地下水中の主な化学反応として陽イオン交換反応を取り上げ、 複数の化学種の 化学的な相互作用を考慮した数値モデルについて述べる。 多成分化学反応項を含む移流分散 方程式を数値的に解くために、 移流分散と化学反応の2つの過程を同時に満足する数値モデ ルを提案する。すなわち、 移流分散方程式中の分散項は差分法を、 移流項は粒子移動法30) (特 性曲線法)を適用し、 さらに未定の化学反応項の数値解を得るためにLevenberg-Marq uard t

法3 1)による非線形最小二乗推定を採用している。 一方、 陰イオンの輸送では、 土壌表面で 吸着し、 輸送に遅れを生じる物質を考慮した物質輸送モデルを適用する。

-40-第3章 数値計算モデル 第2節 陽イオン交換反応を考慮した物質輸送に関する基礎式

地下水の流れや拡散によって移動する陽イオン(可動陽イオン) として、 Naヘr、Mg針、Ca2+

を、 また対応する陰イオンとして土壌に吸着しない 化学的に不活性な陰イオンとしてC1ーを考 える。 次に、 吸着サイトに吸着して移動しない不動の陽イオンとして、 Na+、 K+、 Mg2+、Ca2+を 考える。 すなわち、 移動可能な化学種として5種と移動しないが可動陽イオンとイオン交換 する吸着サイト上の陽イオン4種の合計9種を対象とする。 なおFに関し、 本研究ではpHの 変動は無視されるため対象化学種から除外した。

いま、 空間座擦をzとすると、 5種の可動イオンに対する一次元移流分散方程式は、 次工1 で与えられる。

d

(

^O

If

^]

^L

^L

^{(

^e附]

^)}

^{= sα}

θ

(

^OI

7

²¹

¹

⁺

^バ(

^8[Mg2+

^功

^=ぬ

θ

(

⁰

7

^]

¹

⁺

^ペ(

^8[Na+]

)}

^=九

θ

(17

^+]

^L

_L

_{(

e[K+l

)}

^{= SK}

d

(士川

⁺

^バ(

^8[Cr]

^)}

^=。

、、，，， -ZE-­ • 内《JV fa‘、 、、，，，，nノι• 内《JM fEK

、、IJ qd 内4d rtk

(3.4)

、、，，，Fhd

-q‘u fEK

ここに、[ J:液相中の化学種の濃度(mmol'l・1)、 t :時間、 SCa、 SlIg、 SNa、 Sx:複数の化

学種間で起こる化学反応項、 。:体積含水率である。 また、 式の表記を簡単にするために、

次に示す演算子記号:

、，IIri--J Z 一 Z C 一 d 、σ一 D d 一 h rill〈ilL QU a 一 z r 一 ょ Q U 一 司σ一 日v 可lEL，llJ +

QU a c

rfi，、11K

L

を用いている。 ここに、v:実流速、 D:縦方向分散定数αLと実流速vとの積で与えられる流

---aA d斗A

速依存型分散係数である。 また、 土壌間隙の捻れ率(屈曲度)はDの中に含ませている。

一方、 吸着サイトの不動陽イオンに対しては、 次式が成立する。

a [Ca2+]一

(3. 6)

= SCa

at

a[Mg2+]

�

(3.7)

= SMg

at

a[Na+]

⁼

s

^{Na (3.8)}

at θ[K+]

^=JSK

(3. 9)

at

ここに、 上付きバーは土壌に吸着された陽イオンの濃度(mmol.l・1)を表す。 また右辺は化学 反応項を示す。 式(3. 1) '"'-' (3. 4)と式(3. 6) '"'-' (3. 9)の右辺の化学反応項には次の 関係がある。

SCa = -SCa (3.10)

SMg = -SMg 内《u，，E‘、 - .•••• 、、，，，，.••••

SNa = -SNa 内ed，，目、、 - .•••• 、1Jnノι

SK = -SK 内《JM，，t、 • 1 、‘，，，内4d

ワムA斗A

第3章 数値計算モデル 第3節 化学反応式

固相表面と液相聞での陽イオン交換反応は、 例えばCaれ濃度が液相中で増加した場合に、

以前から吸着していたrが固相から離脱して Caがとイオン交換する。 この反応は瞬時に化学 平衡に達する反応であり、 反応速度論的な取り扱いは不要となる。

いま、 固相をrとすると上記の反応は次式で表される。

2 r ^-K+ + C ^a 2 + = r ^-C a 2 + + 2 K+ (3.14)

この反応は固相と液相間の反応で、 この反応が平衡状態に達した時には以下の関係が成立す る。

K"

^{"f =}

fr-α勺(K+f (r -K+ f(αつ

ここに、 KCa/宜は平衡定数、 ( )は成分の活動度である。 なお、 液相中化学種の活動度につい ては後述する液相中の化学種組成から半経験式を用いて計算できる。

一方、 固相上の陽イオンの活動度を測定あるいは計算することは困難であるため、 陽イオ ン交換平衡にある系の液相中の化学種組成と固相上の陽イオン組成から平衡定数の値を計算 するのは難しい。 従って、 固相上の陽イオンに対しては活動度のかわりに、 固相上の2T­

K十およびT-C _a 2+のモル分率(XK、 Xω;固相上における交換性陽イオンのモル分率に等 しい)を用いて

vV

Xca(K+ r

n ー

.I.l>..Ca/ Kーが(ω μ )

^{(3. 15)}

という量を定義する。 ここに、 J(r'Ca/Kは rを基準にしたCa2+のVanselowの選択係数3 2)であ

る。

同様にNゲとCa2+、 Mg2千とCa2+�こ対して

-43

-ど ω/仙 XNa2(ωサ ーん(地+f

^{，IK qd - - nhu 、、，，，}

vV

xcca( Mg2+)

�""Ca/ 均 一XMg(ω勺

^(3.17)

のように表すことができる。 ここに、J(fCa/Na'

J(fl伽皆

はNa+、 Mg2+を基準にしたCa2+のVanselow の選択係数、

ま

^た、

XN

^a、

X

^Mgは固相上に

お

ける交換性陽イオンのモル分率

で

ある。 ( ) は次 式に

示

^す

よ

^{うにそれぞれの化学種}

の

液相中濃度に対する活動度

で

、 モル濃度と活動度係数の 積で与えられる。 すなわち、 イオン強度を

I

として ^:

q，-

z t

c

マム t 1一2 一一

rt

(3. 18)

(iに関する総和は陽・陰イオンの全てが対象) を算定し、 次いでDaviesの式3 3)に基づく次式:

ゅ←叶占-叶

^{fak qd • --- nud 、1J}

により活動度係数Gjの算定を行なう。 式(3. 18) のCJ はiイオンの濃度、 Zj はiイオンの 電荷、 式(3. 19)のAは温度に依存した定数で、 例えば250Cのとき0.5116の値をとる。 こ

こでは、 添字i=1、 2、 3、 4、 5であり、 それぞれ、 CaZ+、 Mg計、 Naヘr、 C1ーを表している。

さて、 固相上の陽イオン濃度を当量濃度で表現するためは、Vanselowの選択係数の定義で 用いられているモル分率のかわりに当量分率(E胸、 Ex、 EIIg、 Eca)で表現する。 当量分率と は吸着イオン量を、 当量の単位で表したもので、 ナトリウム、 カリウム、 マグネシウム、 カ ルシウムの4成分の場合には

E ι - - qNa

… q Na

⁺

q

^K +

2 q Mg

+

2q Ca

^(3.20)

-44-第3章 数値計算モデル

E，.，=_.. ___^.Lk

q Na + q K + 2q Mg + 2q Ca

、】/-nノ』• 内《uft、

E 'Mg I- = ^2qMg

qNa +qk+2qMg+2qca

、、，，，nJι 円/』• 内4dJ'目、、

E"，_ =

2qCa

ca - qNa+qk+2qMg+2qca

、、，Jqd nノ』• 内《JU，，EE、、

で定義される。 ここで、 qNa、 qx、 qJlg、 qCaは、 それぞれの吸着量Cmmol'lOOg.l)である。

当量分率を用いた選択係数 は式(3.14)の交換平衡式で は:

vGT Eca

(

^K+

f

...Ca/K一EK2

(

^Ca2+

)

^(3.24)

同様にNa吃Ca!+、崎2+とCa!+に対して

vGT Eω

(

^Na+

f

n ー

.L ....Ca/陥 一 品2

(

^ω

勺

vGT Eca

(

^Mg2+

)

...Ca/均 一 九

(

^ω

つ

、、，Jphd nノι

内《d，，t‘、 、、，JFhu nノι内4dfzk

であり、 これ はGaines-Thomasの選択係数34)である。

Vanselowの選択係数およびGaines江、homas の選択係数のいずれも古典熱力学的に は等価 であり35)、 この両者の選択係数も互いに換算することができる。 以下、 本章での選択係数と はGaines

次に、 当量分率CENa、 Ex、 E段、 Eω)と、 吸着されている化学種のモル濃度は、 次の関係 を満足している。

F同ud斗A

一一一�

F'r:'F' Ph ECa

[Ca2+] = CECユーァ

n 乙

九一2

A一nc E C

一00一M

(3.27) [.M〆]=CEC企ENa

n

日=CECPb EK

n

ここに、 CEC:陽イオン交換容量、 ρb • 土のみかけの乾燥密度、 n:空隙率である。 さらに、

回相上の当量分率は、 国相上に交換吸着している陽イオンの総量に等しいとすると、 次式が 成り立つ。

ENa +EK+EMg +ECa ⁼ 1 (3.28)

また、 対象化学種は電気的中性条件として次式を満足しなければならない。

2(Ca2+)+2の1�+)+(Na+)+(K+)=(CI-) 内4d，，，‘、 nノι 、、IJn『Jw

-46

-第3章 数値計算モデル 第4節 数値計算の方法

4. 1 陽イオン交換反応を考慮した数値計算の手順

数値 計算では、 式(3.1) '"" (3.9)を満足する合計9種の化学種の濃度を数値的に求める ことになる。 この場合、 対象とする化学反応過程(式(3.24) '"" (3.26)、 式(3.28)、(3.29)) を満足する化学反応項SCa、Slig、SNa、Sxを求めることが必要となる。 図-3.1には化学反応

を考慮した物質輸送解析のフローチャートを示している。 図中の記号の添字iは化学種 、 添 字jは空間座標を表す差分格子点、 添字nは時間ステップ、 および添字νは繰返し回数 を表 している。 また、 Initial Conditions中のnとαは、van Genuchtenが提案した式3 6)の 不飽和パラメータである。 以下に、 その解法の概略を説明する。

S tep 1)未定の化学反応項OSu IMl を仮定する。 ここでは0九州 =SLj nとした。

また、OSIJ I =Oとした。

Síep 2)仮定値νSI.} 州に対する移流分散方程式(式(3.1) '"" (3.5))の解を特性曲線法 37)

により求める。

S t ^ep 3)仮定値νSME4Iに対 する式(3.6) '"" (3.9)の不動陽イオンの解を差分法により求

める。

Síep 4)仮定値νSLjE#Iに対して求めた可動陽イオンνCJ.J 0+1および不動陽イオンνJ川oflの

値を、 式(3.18)、(3.19)、(3.27)に代入する。

Step 5)得られた可動イオン および不動陽イオンが、 対象とする化学反応過程(式(3.24)

'"" (3.26)、 式(3.28))を満足するように、 化学反応項の推定を行う 。 ここでは、

計算値と予め与えられている値との誤差の二乗和を最小にする化学反応項を求める ことにする。 ここに、 計算 値とは式(3.24) '"" (3.26)では右辺 、 式(3.28)では 左辺であり、予め与えられている値とは式(3.24)では当量分率 (Ex)と 選択係数 (j(GTCtJ/K)聞の実験近似式から得られる選択係数、 式(3.25)、(3.26)では実験値を

平均した選択係数、 式(3.28)では不動陽イオンの当量分率の和lのことである。

また、 誤差の評価は、 図中のStep 4における式(a) '"" (e)によって、 無次元表示 した残差11'""んにより行う。 図中の(a)は式(3.28)に、 式(b)、(C)、(d)、(e) は式(3.24)、(3.25)、(3.26)、(3.29)に対応しており、 式(b)、(C)、(d)の右辺

第2項の j(GUCa/JIg、 j(GUCtJ/袖、J(GT'CtJ/Kは、 可動イオンの活動度と 不動陽イオンの当量分 率を式(3.24)、(3.25)、(3.26)の右辺に代入して算定される値である。 残差 二乗

弓，，A斗A

和が最小となるような未定の非線形パラメータである化学反応項 νS，・j ulを繰返し 計算により推定する。 ここでは、 Levenberg - Marquardt法3 1)を適用する。 化学 反応項 νslj EHの最適値が得られた場合には次の格子点に移動し、 全格子点終了後 には、 以上の算定ルーチンを最終計算時間まで繰り返す。

S tep 3においては、 前の時間ステップにおける吸着サイトの濃度に比べ、化学反応項νS人j aH の推定値と差分時間間隔L1tとの積が過大な値となることがある。 すなわち、 推定値に制約を 課していないので、 νC"j nf-Jが負値として算定される場合もある。 この場合は、 化学反応項の 仮定値を現在の推定値より小さい値に再度設定し、Step 2とStep 3 の計算をやり直すこと によって処理している。

4. 2 化学反応項の推定方法

化学反応項の推定には、 式(3.24) ""' (3.26) および式(3.28)、(3.29) を満足するよう に次の5つの無次元化された残差関数を用いている。

f1=W1(1-

2，

_V

EC

₁₎

12 = w2{log(悶�a) -log( �Na)}

13 = w3{log(悶�)-log(�K)}

九 = w4{log( Kg;Mg) -log( �Mg)}

人 =ws{lー やvc:;1叫;1叫1ぺ�l) /り

ここに、 cは濃度、 記号の添字fは化学種、 添字jは空間座標を表す差分格子点、 添字nは 時間ステップ、 Wk (k=1""'5)は重み係数、 および添字νは繰返し回数を表している。

数値計算では、 残差関数11 ""' 15の2乗和が最小となる化学反応項の推定を行うために、

非線形最小2乗法の解法のlつであるLevenberg - Marquardt法31)を用いている。

ここで重み係数の値に関し、 後述する第4章では wt=l (k=l ""'5)を用いて計算を行ってい る。 また後述する第5章のように最終計算時間に到達する聞に濃度の境界条件が変化する数 値シミュレーションを実行する際には、 重み係数wt=l (k=l ""'5)では数値計算上、 化学反応 項の推定値に誤差が蓄積され精度が落ちることがある。 従って第5章の数値計算では、 重み 係数による残差関数I} ""' 15の誤差評価を行い、 重み係数wt=13 (k=l)、wt=l (k=2�5)を用

いて計算することで化学反応項の推定に伴う誤差の蓄積を回避している。

ドキュメント内 農耕地域地下水系の水質形成機構に関する研究 (ページ 47-57)