はじめに

この文章では, Weidenschilling and Lewis (1973), Atreya and Romani (1985) で開発された ECCM (Equilibrium Cloud Condensed Model) と, 我々の用いている RAND 法とを比較検討し, それぞれの利点と欠点を評価する.

地球での経験から明らかなように, 気液平衡を考慮するだけなら ECCM のやり方は簡潔で使いやすい. しかし, 溶液のような「混じり合ったもの」や化学反応も含めて扱うには RAND 法の方が向いているだろう. 以下の文章ではその点について検証していきたい.

第1章 ECCM の概要

エネルギー保存則 dS = 0(dQ = 0) を陽に解く. 相変化と化学反応を考慮した場合, エネルギー保存則は以下のように書ける.

c̄_{p}dT − vdp + ∑ L_{k}dX_{k} + ∑ L_{r}dX_{r} = 0.  (1.1)

ここで c̄_{p} は定圧モル比熱の平均値, T は温度, v はモル体積, p は圧力, L_{k} はモル当たりの潜熱, dX_{k} は凝結に伴うモル数変化, L_{r} はモル当たりの化学反応の反応熱, dX_{r} は化学反応に伴うモル数変化である. 気塊に含まれる全モル数を 1 mol とみなせば, dX_{k}, dX_{r} をモル分率の変化と読み換えることができる.

エネルギー保存則 (1.1) 式の変数は, dT, v, dp, dX_{k}, dX_{r} である. 理想気体, 静水圧平衡, 化学平衡を仮定した時に成立する関係式を用いて, これらの変数を温度変化 dT と高度変化 dZ の項で表現することを目指す.

(1.1) 式を見てわかる通り, エネルギー保存の式を立てるには, 大気中でどのような凝縮, 化学反応が生じているか予め把握しておく必要がある.

1.1 エネルギー保存則の各項の表現

1.1.1 仕事 −vdp の表現

仕事 dW に対応する項, −vdp の変形を行う. 静水圧平衡の式と理想気体の状態方程式を用いると,

dp = ρgdz = \frac{1}{v} gdz  (1.2)

である. それを用いて仕事項を変形すると,

−vdp = gdz  (1.3)

となる.

1.1.2 相変化に伴う潜熱の表現

相変化に伴う潜熱項の変形を行う. これはモル分率変化を分圧の形に変形した上で, クラウジウス– クラペイロンの式を用いて分圧を温度の関数として書き換える. モル分率は X_{k} = p_{k}/p と書けるので, それを微分すると,

dX_{k} = \frac{dp_{k}}{p} − (\frac{p_{k}dp}{p^{2}})  (1.4)

となる. 上式中の dp を静水圧平衡で, dp_{k} をクラウジウス–クラペイロンの式で dT の式に変形する.

クラウジウス–クラペイロンの式は

dp_{k} = \frac{p_{k}L_{k}dT}{RT^{2}}  (1.5)

のように書ける. クラウジウス–クラペイロンの式は, 気液平衡が成立するときに成り立つ関係式で, 気相と液相の化学ポテンシャルの値が等しいという条件から求まる.

1.1.3 化学反応に伴う反応熱の表現

化学反応に伴う反応熱項の変形を行う. これを行うためには, 具体的な化学反応式とその反応の平衡定数を予め把握しておく必要がある. 平衡定数は一般に反応に関係する化学種の分圧と温度の関数で書かれるので, それを微分することでモル分率の変化を温度と高度 (圧力) の関数で書く事ができる.

例えば Areya and Romani (1985) によると, 硫化水素とアンモニアから硫化アンモニウムの生成される反応では,

log_{10} (p_{NH_{3}} · p_{H_{2}S}) = log_{10}(K_{p}) = 6.467 − \frac{10834}{T}  (1.6)

の関係式が成立する. 但し単位系は cgs 単位系である. (1.6) 式を微分することで, モル分率変化を温度変化と圧力変化で表現することができる.

異なる化学反応を考えるたびに, 反応式を扱うためのルーチンを追加する必要がある. 大気中で生じる化学反応を予め掌握する必要があると言える.

1.1.4 水溶液での潜熱解放の表現

水溶液の生成される際の潜熱解放について考える. 問題となるのは, 水溶液と気相が平衡にあるときの飽和蒸気圧の表現である. 純物質の飽和蒸気圧のように, 温度だけの関数で書けるなら問題ないが, 水溶液の場合は必ずしもそうはならない.

アンモニア水溶液の場合, 温度と水溶液濃度によって水溶液に溶け込む気相の飽和蒸気圧が変化する. Atreya and Romani (1985) では, 潜熱は温度に対して線形であると仮定し, クラウジウス– クラペイロンの式を積分して求めた式,

ln (p_{k}) = A + \frac{B}{T} + C ln (T).  (1.7)

を用いて水溶液と平衡状態にある気相の飽和蒸気圧を表現している. 但し A, B, C は濃度の関数となる. 彼らはこの A, B, C を実際のデータにフィットするよう, 3 次のスプライン関数として表現している. 結果として飽和蒸気圧を計算するためには, 水溶液濃度と温度を知る必要がある. ある高度での温度は与えられるので, 実際の計算は以下の手順で行う.

(1) 水溶液濃度を仮定し, その濃度での飽和蒸気圧を計算する

(2) (1) で求めた飽和蒸気圧を用いて, 水とアンモニアが凝縮するか判定する. 両者が凝縮しない場合は (1) に戻る (仮定する水溶液濃度を変えてやり直す).

(3) (2) で水とアンモニアの両方が凝結するならば, その時の水溶液濃度を計算する. 得られた水溶液濃度が (1) で仮定した水溶液濃度に等しいか否かを判断する. 等しくない場合は (1) に戻る.

Atreya and Romani (1985) では, アンモニア水溶液に対する硫化水素の溶解も考慮している. 彼らは Leyko (1964) によって開発された経験式を用いている. この経験式は硫化水素の分圧を温度, 溶液中のアンモニア濃度と硫化水素濃度によって表現したものである.

p_{H_{2}S}^{∗} = exp (22.221 − \frac{5388}{T})C_{H_{2}S}^{1.130} (\frac{C_{H_{2}S}}{C_{NH_{3}}} )^{1.8953}  (1.8)

ここで p^{∗} は飽和蒸気圧, C は濃度である. すでに温度とアンモニア水溶液でのアンモニア濃度は分かっているので, 硫化水素の濃度をパラメタとして硫化水素の飽和蒸気圧を計算することができる. 仮定した硫化水素濃度と, 硫化水素の飽和蒸気圧と分圧から求まる硫化水素濃度が等しくなるまで計算を反復的に繰り返す.

アンモニア水溶液の潜熱は平均値 L_{s} を用いて,

L_{s} \frac{dX_{H_{2}O}}{1 − C_{NH_{3}}}  (1.9)

と表現する. Atreya and Romani (1985) では, 水溶液に伴う潜熱解放はだいたい一定に保たれると述べている. また硫化水素のアンモニア水溶液への溶解に伴う潜熱に関しては, Atreya and Romani (1985) では触れられていない. 考慮していないのかもしれない.

Weidenschilling and Lewis (1973) は Atreya and Romani (1985) に比べて, 水溶液の雲の扱いが 2 つの点で異なる.

(1) 飽和蒸気圧の表現 Weidenschilling and Lewis (1973) では水溶液と平衡にある気相の飽和蒸気圧を, 同じ温度における純粋な得気相の蒸気圧と水溶液濃度に関する 5 次の方程式としてフィッティングした. Atreya and Romani (1985) はこの方法では求まる低温での蒸気圧の値は現実的ではないと批判している.

(2) 水溶液濃度の計算 Weidenschilling and Lewis (1973) と Atreya and Romani (1985) は, 水溶液雲の雲底での水溶液濃度を決める方法は同じである. Atreya and Romani (1985) は雲底以外の水溶液濃度の決め方も雲底での方法に同じであるが, Weidenschilling and Lewis (1973) は雲底以外の水溶液濃度を別のやり方で求めている. 彼らはある濃度 C の溶液が存在する平衡な系が微小量だけ断熱膨張するときの濃度変化の式を用いて水溶液濃度の変化を計算している.

1.2 数値解法

ECCM の計算では, 凝結するかどうかの判定式を別途持つ必要がある. 凝結と硫化アンモニウムの生成反応を考慮した場合 (溶液の凝結は添字 k に含まれる), 温度変化の式は以下の様に書ける.

\frac{dT}{dz} = −\frac{m̄g}{c̄_{p}} \frac{[1 + \frac{1}{RT} \{∑_{k} L_{k}X_{k} + \frac{2(X_{H_{2}S}·X_{NH_{3}} )L_{r}}{(X_{H_{2}S}+X_{NH_{3}} )} \}]}{[1 + \frac{1}{c̄_{p}T^{2}} \{∑_{k} \frac{L_{k}^{2}X_{k}}{R} + \frac{(X_{H_{2}S}·X_{NH_{3}} )}{(X_{H_{2}S}+X_{NH_{3}} )} · L_{r} · 10834\}]}.  (1.10)

この式で既知の量は, 高度変化 dz (Weidenschilling and Lewis, 1973 では 1 km), 平均分子量 m̄, 重力加速度 g, 平均比熱 c̄_{p}, 気体定数 R, 潜熱 L_{k}, 反応熱 L_{r} である. 潜熱と反応熱は定数とみなしているようだ. 変数は, 温度変化 dT, 温度 T, モル分率 X_{k}, X_{H_{2}O}, X_{NH_{3}} である. 温度とモル分率を何らかの方法で決めることで, 高度変化したときの温度変化を求めることができる.

ある高度での温度とモル分率を与えれば, そこでの温度減率が計算できる. その温度減率が dz の区間で維持されるとみなせば, 高度が dz だけ離れた点での温度を求めることができる. その意味で温度とモル分率は 1 つ前のステップでの値を用いる. Atreya and Romani (1985) ではクッタ-シンプソンの 4 次式を用いると述べているので, 実際の計算では前のステップ数点の値を用いているようだ.

モル分率は, 純物質の飽和蒸気圧の経験式, 化学反応の平衡定数に関する経験式, 水溶液上の飽和蒸気圧の数値的な経験式, といった判定式から判断する. これらの経験式は蒸気圧の最大値を表すものである. 全圧とモル分率の積から蒸気圧を計算し, それが経験式から求まる飽和蒸気圧を越えないか否かを判定する. 蒸気圧がその最大値を越えなければモル分率の値はそのまま保管され, 蒸気圧がその最大値を越えた場合には最大値に固定するという操作を行う. この方法では, 凝結量が判定式を適用する順番に依存する可能性がある.

凝結物質は「落下」するかどうかはきちんと明記されていない. しかし, 結果の絵から判断するに, 気塊のみが上昇する過程を考えているものと想像される.

第2章ギブス自由エネルギー最小化法

ギブス自由エネルギー最小化法は, 平衡状態は熱力学関数が最小化された状態として表現できるという化学熱力学の定義を用いた方法である. 熱力学関数の中からギブス自由エネルギーを選択したのは, 大気の状態を温度, 圧力, 大気組成として特徴づけたいからである.

ギブス自由エネルギー最小化法を用いると, ある温度, 圧力における平衡組成を計算することができる. 温度と圧力をパラメタとして変化させることで, 温度圧力空間における平衡組成のコンター図を書く事ができる. 平衡組成から解析的にエントロピーを計算することができるので, 温度圧力空間においてエントロピーの等値線をプロットし, エントロピーの等値線を結んだものが求めるべき大気の断熱状態となる.

擬断熱過程を考える場合には, 多少やっかいになる. 擬断熱過程とは, 気体が凝縮した際に凝縮物が系から離脱すると考えた場合に実現する断熱減率で, 離脱した凝縮物質も含めた全エントロピーが保存すると考える. この場合には系に含まれる元素量が変化するので温度圧力空間におけるエントロピーのコンターマップを引くのは困難である. ある温度, 圧力, 元素存在度にあった空気塊を追いかけていく必要がある. 圧力を dp だけ変化させた際に空気塊で凝縮が生じたならば, 次に圧力を変化させるときに保存される元素量は, 気相のそれである.

2.1 最適化法による解法

熱力学変数として, 温度, 圧力, 組成を選択する. このとき大気の熱力学状態を与える適切な熱力学関数はギブス自由エネルギーとなる. 大気の平衡状態は, 元素数保存の条件の下でギブス自由エネルギーが最小化した状態であるとする. 温度と圧力を与えた時, ギブスの自由エネルギーと元素数保存式は以下のように書ける.

G(T, p, n_{i}^{ϕ}) = ∑_{ϕ} ∑_{i} μ_{i}^{ϕ}(T, p, n_{i}^{ϕ})n_{i}^{ϕ}
= ∑_{ϕ} ∑_{i} \{μ_{i}^{∘}^{ϕ}(T) + RT ln \frac{n_{i}^{ϕ}}{∑ n_{i}^{ϕ}} + α^{ϕ}RT ln \frac{p}{p_{0}} \} n_{i}^{ϕ}  (2.1)
∑_{ϕ} ∑_{i} a_{ie}^{ϕ}n_{i}^{ϕ} = B_{e}  (2.2)

元素数保存の条件下でギブスの自由エネルギーを求めるという, 条件付き極値問題に帰結することができる. 要はラグランジュの未定乗数法を用いて解けば良い.

実際に数値的にギブス自由エネルギーを最小化するための最適化法として, White et al (1958) によって開発された RAND 法を用いる. RAND 法の骨子は, 逐次近似法によってギブス自由エネルギーが最小となる物質存在量を求めることにある. ギブス自由エネルギーを適当な物質存在量のまわりで物質量に関する 2 次の関数に近似し, ラグランジュの未定乗数法を用いて元素数保存条件の下で物質存在量の極値を計算する. 得られた物質存在量が収束するまで同様の操作を繰り返す.

この RAND 法を適用できない状況がいくつか存在する.

• ギブス自由エネルギーを組成に関する 2 次の関数に近似できない場合には, RAND 法の前提条件が成り立たない.

• 物質量が 0 に近付く場合には化学ポテンシャルが発散する.

また,

• 初期近似解 (計算に与える初期組成) が十分に最適解 (平衡組成) に近くないと, 最適解に収束しないことがある. そのため, 初期近似解の選び方が重要となる.

ということも言える. 例えば同じ RAND 法を利用している Wood and Hashimoto はなんでも混ざるという近似をしているので (凝縮物質の相をすべてまとめて 1 つとしているので), 上記の適用条件が当てはまらない.

2.2 化学ポテンシャル中の各項の表現

2.2.1 エントロピーの表現

エントロピーは Maxwell の関係式から求めることができる. ギブス自由エネルギーとエントロピーの関係は,

S = − (\frac{∂G}{∂T})  (2.3)

なので, G が温度, 圧力, 組成の関数として書けているのであれば, 解析的に解くことができる. 実際, 扱った全ての物質に関して解析的に計算している.

2.2.2 物性値の表現

必要となる物性値は, 系に含まれる化学種の標準化学ポテンシャルである. 標準ポテンシャルを求める方法としては, (1) 比熱および基準状態での温度・圧力を用いて解析的に計算, (2) 気液平衡時に気相と液相の化学ポテンシャルが等しいことを利用して飽和蒸気圧の式から計算, (3) 化学反応式の圧平衡定数と反応時の化学ポテンシャルが等しい, ことから求まる.

比熱は NIST JANAF の熱力学テーブルから情報を得た. 離散的なデータとして与えられているので, スプライン関数で連続的に表現している.

飽和蒸気圧の式は, 化学便覧, 化学工学便覧から得た. 適用範囲外に外挿する際は, 温度 0 K でエントロピーが 0 になるように補間している.

2.2.3 溶液の表現

理想溶液を仮定するなら, その化学ポテンシャルは,

μ = μ^{∘} + RT ln \frac{n_{i}^{ϕ}}{∑ n_{i}^{ϕ}}  (2.4)

と書けるので, 何らかの方法で標準状態の化学ポテンシャルさえ分かれば, あとは機械的に扱うことができる.

第3章議論

ECCM の方法とギブス自由エネルギー最小化法は, 化学平衡状態を計算しているという点で全く等価である. それぞれの化学種の物性や反応の情報を熱の形 (温度変化 dT) で導入するか, 化学ポテンシャルの形で導入するかが異なるに過ぎない. しかし, 化学反応や溶液を扱うには, 化学ポテンシャルの形で物性を取り扱う方が圧倒的に簡単である. 化学ポテンシャルを計算するまでは面倒かもしれないが, その後の操作を一元的に扱えることは大きなメリットであると考える.

経験式を用いるという点では, ECCM も RAND 法も同じである. RAND 法では化学ポテンシャルを計算するために, ECCM では飽和蒸気圧を計算するために, 経験式を用いている.

ECCM に比べてギブス自由エネルギー最小化には以下のような利点がある.

• 大気中でどのような化学反応, 水溶液の雲が生じているか考える必要がない. ECCM では大気中で生じる凝結・化学反応を予め把握しておかないと基礎方程式自体が立てられない.

• 水溶液の雲の扱いが簡潔である. RAND 法では水溶液の雲も純粋な凝縮相と同様に扱うことができる. ECCM では, 変数が濃度, 温度, 圧力にも関わらず圧力しか与えることができないので, 濃度を仮定してシューティングする必要があった.