木星型惑星大気を念頭においた

熱力学計算手法

A thermodynamic calculation scheme for jovian planetary atmospheres

北海道大学大学院理学研究科地球惑星科学専攻

杉山耕一朗

SUGIYAMA Ko-ichiro 2001 年 2 月 17 日

\underline{木星型惑星大気の熱力学計算目次 1}

目次

1 はじめに 1

2 計算手法の検討 3 2.1 全体の枠組 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 大気の平衡状態の記述 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.3 大気の断熱変化の記述 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3 大気の平衡状態の記述 6

3.1 大気の平衡状態の計算方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3.1.1 大気が理想気体のみから成る場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3.1.2 多相系への拡張 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.1.3 多相系への拡張 2 (相の扱いを物質の違いとする場合) . . . . . . . . . . . . . . . . . 12

3.2 数値計算上のテクニック . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2.1 相の扱い (その 1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2.2 対数項 (ln(m_{i}^{ϕ}/m^{ϕ})) の扱い . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2.3 負のモル数の処理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2.4 計算を停止させるための条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

4 大気の断熱変化の記述 16 4.1 大気の断熱変化の計算方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 4.2 数値計算上のテクニック . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

4.2.1 相の扱い (その 2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.2.2 熱力学関数 S の計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

5 化学ポテンシャルの与え方 20 5.1 標準化学ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 5.2 気体の標準化学ポテンシャルの計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 5.3 固体・液体の標準化学ポテンシャルの計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

5.3.1 飽和蒸気圧が知られている場合の計算 . . . . . . . . . . . . . . . . . . . . . . . . . 21 5.3.2 化学反応の圧平衡定数が知られている場合 . . . . . . . . . . . . . . . . . . . . . . . 23

6 結果 26

6.1 モデルの信頼性の検証 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 6.1.1 G を用いた他の物質量の計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 6.1.2 地球大気の乾燥断熱減率の計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 \underline{木星型惑星大気の熱力学計算目次 2}

6.2 地球大気の断熱温度減率と凝縮物質の鉛直分布 . . . . . . . . . . . . . . . . . . . . . . . . . 29 6.3 木星条件での計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

7 まとめ 38 8 参考文献 39

Abstract

木星表層はアンモニアの雲に覆われており, 雲の発生・消滅が観測されている. しかしアンモニア雲層高度以下の領域は光学的にほとんど観測できないために, 木星の対流がどのような構造をしているかはよく分かっていない.

木星大気対流を実際にシミュレーションするためには, 大気の平均的な状態として断熱温度減率と凝縮物質の鉛直分布を把握することが重要である. 本研究では木星型惑星大気の断熱温度減率と凝縮物質の鉛直分布とを計算するための手法を検討する.

木星大気の相変化や化学反応を計算した例として Atreya and Romani (1985) がある. 彼らは熱力学関数 S の保存式を理想気体の状態方程式と潜熱・反応熱を用いて熱力学変数 T, p, n_{i} の式として定式化した. この定式化では空気塊において生じうる全ての化学反応を知る必要がある. しかし木星型惑星大気における化学反応の全容は明らかになっていない. 組成の違う大気を考える場合には, 生じうる化学反応が変わるので, 数値コードを書き換えねばならない. そもそも平衡状態を仮定しているので, 化学反応に関する情報がなくとも断熱温度減率と凝縮物質の鉛直分布は計算できるはずである.

本研究では大気中で生じうる反応式を考えずに済み, さらに大気組成を簡単に変更することのできる手法を検討する. それを行うためには (1) 大気の平衡状態の記述, (2) 大気の断熱変化の記述, を検討することが必要となる.

(1) 大気の平衡状態の記述であるが, 熱力学変数として温度 T, 圧力 p, 物質存在量 n_{i} を選択する. このとき大気の状態を与える適切な熱力学関数は S ではなく G である. 大気の平衡状態は, 熱力学関数 G が最小化された状態であるとする. 温度・圧力を与えると, 熱力学関数 G は,

G = ∑ μ_{i}n_{i},

と書ける. ここで μ_{i} は化学ポテンシャルである. 上式より熱力学関数 G の最小値は, 物質存在量 n_{i} をパラメータとして変化させることによって求めることができる. 大気の状態を熱力学関数 G によって記述したので, 潜熱・反応熱の代わりに化学ポテンシャル μ_{i} のリストを与えればよい. 化学ポテンシャルは既存の物性データから求めることができる.

(2) 大気の断熱変化の記述であるが, 大気の断熱変化は熱力学関数 S の保存として記述することができる. 熱力学関数 S は G を用いて以下のように与えられる.

(\frac{∂G}{∂T})_{p,n_{i}} = −S.

温度・圧力空間で dS = 0 の曲線の通る軌跡を順にたどれば, 断熱温度減率と凝結物質の存在量を求めることができる.

大気成分気体やその凝結成分を理想気体/液体/固体/溶液と仮定し, 上記の計算手法に則って実際に木星大気における断熱温度減率と凝結物質の鉛直分布を計算した. 大気組成はボイジャーの観測結果から予想されるものとした. 凝結物質の種類と各凝結物質量のピーク値の圧力は, Atreya and Romani (1985) の結果にほぼ等しい. したがって本研究によって開発された計算手法を用いることによって, 他の木星型惑星に対しても妥当な断熱温度減率と凝結物質の鉛直分布が得られると期待される.

\underline{木星型惑星大気の熱力学計算第 1 章はじめに 1}

第 1 章

はじめに

木星型惑星 (木星・土星・天王星・海王星) は水素・ヘリウムを主成分とする分厚い大気をまとっている. その表層は厚い雲によって覆われ, 外部からの観測によって雲層の下をうかがい知ることは困難であり, 未だ成功していない. 外惑星のうち, 最も観測の進んでいる木星も例外ではない. 木星の表層はアンモニアの雲に覆われており, 東西に延びる縞と帯, 大赤斑や白斑と呼ばれる渦等の特徴的な構造が見られる. そのような大規模の構造の中に, 小スケールの雲塊が発生・消滅していることがガリレオ探査機の観測によって明らかにされている. しかしあまねく存在するアンモニア雲層のために雲層高度以下の観測は断片的にならざるをえない. その雲層の下の温度構造, 物質の分布, 大気循環がどうなっているのかについては, 依然謎なままである.

木星型惑星の大気の循環構造を解明するには, (1) 木星を地上観測もしくは探査機による観測, (2) 理論的なシミュレーションからの予想, という方法がある. しかし前述したように, 地上観測やボイジャー・ガリレオの観測ではアンモニア雲層以下の領域の大気循環は明らかにされていない. 複数の探査機によるアンモニア雲層以下の領域の観測が行われない限り, 大気の循環構造は観測的に明らかにされないであろう. そのため現在は (2) のシミュレーションからの予想を行うことで木星大気の循環構造を予想することが最も現実的である. シミュレーションの予想をもとに, 新たなる観測計画を立案できる可能性もある.

木星大気のシミュレーションを行うためには, 大気の平均的な状態として断熱温度減率と凝縮物質の鉛直分布とを把握することが重要である. 過去の木星型惑星大気の断熱温度減率と凝縮物質の鉛直分布を計算した例として Atreya and Romani (1985) がある. 彼らは熱力学関数 S の保存式を理想気体の状態方程式と潜熱・反応熱を用いて熱力学変数 T, p, n_{i} の式として定式化した. この定式化では空気塊において生じうる全ての化学反応を知る必要がある. しかし木星型惑星大気における化学反応の全容は明らかになっていない. 組成の違う大気を考える場合には, 生じうる化学反応が変わるので, 数値コードを書き換えねばならない. そもそも平衡状態を仮定しているので, 化学反応に関する情報がなくとも断熱温度減率と凝縮物質の鉛直分布は計算できるはずである.

そこで本研究では大気中で生じる反応式を考えずに済み, さらに大気組成を簡単に変更することのできる熱力学計算手法を検討する. そしてその手法を用いて大気の断熱温度減率と凝縮物質の鉛直分布を計算する.

熱力学計算手法を開発するためには以下の点に関して検討を行う必要がある.

• 大気の平衡状態をどのように記述するか.

• 大気の断熱変化をどのように記述するか.

それをふまえ, 以下の章では上記の 2 点について検討を行い, それぞれを数値的に扱うための方法を開発する.

最後に本計算手法を用いて計算した木星大気の断熱温度減率と凝縮物質の鉛直分布を Atreya and Romani (1985) の計算結果や, ガリレオプローブの観測結果と比較する. それにより本手法の妥当性を検討する.

\underline{木星型惑星大気の熱力学計算第 2 章計算手法の検討 3}

第 2 章

計算手法の検討

本章では大気の断熱温度減率と凝縮物質の鉛直分布を計算するための熱力学計算手法に関して検討する.

2.1 全体の枠組

図 2.1 のように, 圧力 p において温度 T, 圧力 p, 物質存在量 n_{i} の空気塊を考える. その空気塊において, • 空気塊は平衡状態にある.

• 空気塊に含まれる物質は全て理想気体/ 液体/ 固体/ 溶液である.

• 空気塊は運動していない. 即ち上昇流や下降流といった流れは存在しない. の条件が成立するものとする. 但し 2 番目の条件は以下の 2.2 節, 2.3 節の議論とは本質的に関係ない. 2 番目の条件は物質の化学ポテンシャルを表現するときに必要となる条件である.

図 2.1: 対流熱力学モデルの系の設定. 空気塊の圧力を断熱的に dp づつ変化させる.

空気塊の断熱変化を仮定し, 空気塊の圧力を dp づつ変化させた時の温度と物質存在量を計算する. 圧力に対する温度の変化より大気の断熱温度減率が求まる. 圧力に対する凝縮物質の存在量より凝縮物質の鉛直分布が求まる.

2.2 大気の平衡状態の記述

大気の状態を表現するための熱力学変数として温度 T, 圧力 p, 物質存在量 n_{i} を選択する. 大気の断熱温度減率と凝縮物質の鉛直分布を計算したいためである. このとき大気の状態を与える適切な熱力学関数は G である.

大気の状態を熱力学関数 G で表現することに決めたので, 大気の平衡状態は熱力学関数 G の最小化された状態である. 温度・圧力を与えると熱力学関数 G は,

G = ∑ μ_{i}n_{i},  (2.1)

と書ける. ここで添字 i は各物質を意味する添字, μ_{i} は化学ポテンシャルである. (2.1) 式より熱力学関数 G の最小値は, 物質存在量 n_{i} をパラメータとして変化させることによって求めることができる. したがって大気の平衡状態は

• 予め与えられた大気の温度 T ・圧力 p

• (2.1) 式の G を最小化する物質存在量 n_{i}

によって記述される. また同時に平衡時の熱力学関数 G の値も決まる.

大気の状態を熱力学関数 G によって表現することにしたので, 潜熱・反応熱の代わりに化学ポテンシャル μ_{i} のリストを与えればよい. 化学ポテンシャルは既存の物性データから求めることができる. 今や Atreya and Romani (1985) のように, 大気中で生じる全ての化学反応を知っている必要がなくなった.

2.3 大気の断熱変化の記述

大気の断熱変化は熱力学関数 S の保存として記述することができる. 熱力学関数 S は G を用いて以下のように与えられる.

(\frac{∂G}{∂T})_{p,n_{i}} = −S  (2.2)

である. 大気の状態を熱力学関数 G によって記述したために, 温度・圧力空間における G の分布が求まる. G の分布より温度・圧力空間における熱力学関数 S の分布が計算できる. 熱力学関数 S の等値線の通る軌跡を dp づつ順にたどれば, 圧力に対する温度と物質存在量を求めることができる. それらは断熱温度減率と凝縮物質の鉛直分布に相当する.

大気が乾燥 (= 凝縮物質が存在しない) している場合は温度・圧力空間の S の等値線を順にたどれば良い. しかし大気が湿潤 (= 凝縮物質が存在する) になった場合は, たどるべき S の等値線を変化させる必要がある. 大気では凝縮物質は降水過程を通じて気塊から除去されるので, その過程を組み込む必要があるからだ. この過程は「擬断熱 (pseudo moist adiabat) 過程」と呼ばれる.

• 擬断熱温度 (pseudo moist adiabat) 過程: 空気塊に含まれる可凝縮気体が凝縮した場合, その凝縮物質は空気塊からただちに離脱する.

擬断熱過程にある場合, 大気の平衡は保たれるが空気塊の熱力学関数 S 及び凝結性成分に関する元素量は保存されない. 擬断熱過程では, 大気 (気相のみ) の S の値と離脱した凝縮物質の S の値が保存されているものとみなす. すなわち凝縮物質が生成された場合, 大気の熱力学関数 S は減少する.

図 2.1 を例にして考える. 番号 2, 3 の空気塊において凝縮物質が生じているとする. それぞれの空気塊の熱力学関数 S の値をそれぞれ S_{1}, S_{2}, S_{3} とする. また番号 2, 3 で生じた凝縮物質のみの S の値をそれぞれ S_{2}^{c}, S_{3}^{c} とする. このとき,

S_{1} = S_{2} + S_{2}^{c}
S_{2} = S_{3} + S_{3}^{c}
S_{1} = S_{2} + S_{2}^{c} = S_{3} + S_{2}^{c} + S_{3}^{c}

が成立する.

\underline{木星型惑星大気の熱力学計算第 3 章大気の平衡状態の記述 6}

第 3 章

大気の平衡状態の記述

3.1 大気の平衡状態の計算方法

大気の温度・圧力を与えると, 熱力学関数 G は以下のようにかける.

G = ∑ μ_{i}n_{i}.  (3.1)

と書ける. ここで添字 i は大気に含まれる i 番目 (0 ≤ i ≤ N) の物質を意味し, μ_{i} は i 番目の物質の化学ポテンシャルを意味する. また平衡状態を計算する際に, 大気中に含まれる元素量は変化しないという条件を付ける.

∑_{i} a_{ie}n_{i} = B_{e}.  (3.2)
n_{i} ≥ 0  (3.3)

ここで添字 e は系に含まれる e 番目 (0 ≤ e ≤ M) の元素を意味し, また a_{ie} は i 番目の物質に含まれる e 番目の元素の数, B_{e} は系全体に含まれる e 番目の元素の個数を意味する. 例えば元素 H (e = 1), O (e = 2) と物質 H_{2}O (i = 1) を考えると a_{1,1} = 2, a_{1,2} = 1 となる. (3.3) 式はモル数は常に正の数であることを意味する.

本章では (3.1) 式と (3.2) 式を連立し, G の最小値を与える物質存在量の数値的計算方法を議論する.

3.1.1 大気が理想気体のみから成る場合

大気が理想気体のみから成る場合を考える. この場合, White et al. (1958) らによって開発された「RAND 法」を用いる (Van Zeggerren and Storey, 1970).

この計算方法の骨子は, 逐次近似法によって熱力学関数 G が最小となる物質存在量を求めることにある. 熱力学関数 G を適当な物質存在量のまわりで 2 次の関数に近似し, その最小値計算する. 得られた物質存在量が収束するまで同様の操作を繰り返す.

適当な物質存在量を m_{i} とする. m_{i} よりも熱化学平衡に近い物質存在度を n_{i}(= m_{i} + δm_{i}) とする. 平衡時の熱力学ポテンシャル G を m_{i} のまわりでテーラー展開し, δm_{i} に関する 2 次の微小量まで近似式 Q(n_{i}) を作る.

\underline{木星型惑星大気の熱力学計算第 3 章大気の平衡状態の記述 7} 図 3.1: 逐次近似法の概要. 適当な物質存在度のまわりで熱力学関数 G を 2 次関数に近似し極値を求める. 得られた極値を使って同様の操作を繰り返し, 平衡時の物質存在量を求める.

Q(n_{i}) = G(m_{i}) + ∑_{i} (\frac{∂G}{∂m_{i}} ) δm_{i} + \frac{1}{2} ∑_{i} ∑_{j} (\frac{∂^{2}G}{∂m_{i}∂m_{l}} ) δm_{i}δm_{j},  (3.4)

と書ける. ここで理想気体の化学ポテンシャルの定義式を用いると,

(\frac{∂G}{∂m_{i}} ) = μ_{i} = μ_{i}^{∘}(T) + RT ln (\frac{m_{i}}{m} )p.  (3.5)

と書ける. 但し

m = ∑_{i} m_{i}  (3.6)

という関係式が成立する. (3.5) 式を用いて (3.4) 式の右辺第 2 項のカッコ内を書き下すと,

(\frac{∂^{2}G}{∂m_{i}∂m_{l}} ) = RT (\frac{δ_{il}}{m_{i}} − \frac{1}{m}) ,  (3.7)

となる. (3.5), (3.7) 式を (3.4) 式に代入し, RT で割ることによって,

\frac{Q(n_{i})}{RT} = \frac{G(m_{i})}{RT} + ∑_{i} [\frac{μ_{i}^{∘}(T)}{RT} + ln p + ln \frac{m_{i}}{m} ] δm_{i} + \frac{1}{2} [∑_{i} \frac{(δm_{i})^{2}}{m_{i}} − \frac{(δm)^{2}}{m} ] ,  (3.8)

が導かれる.

\overline{main.tex 1991/02/17(杉山耕一朗)}

(3.8) 式の Q(n_{i})/RT と質量保存則から, ラグランジュの未定係数法を用いて極値の条件を与える. L は以下のように定義する.

L = \frac{G(m_{i})}{RT} + ∑_{i} [\frac{μ_{i}^{∘}(T)}{RT} + ln p + ln \frac{m_{i}}{m} ] δm_{i} + \frac{1}{2} [∑_{i} \frac{(δm_{i})^{2}}{m_{i}} − \frac{(δm)^{2}}{m} ]
− ∑_{e} χ_{e}[∑_{i} a_{i,e}n_{i} − B_{e}].

よって極値を取る際に

\frac{∂L}{∂n_{i}} = [\frac{μ_{i}^{∘}(T)}{RT} + ln p + ln \frac{m_{i}}{m} ] + [\frac{n_{i}}{m_{i}} − \frac{n}{m}] − ∑_{e} χ_{e}a_{ie} = 0,  (3.9)

を満足する. (3.9) 式と (3.2) 式, n = ∑ n_{i} を連立することで, 変数 n_{i}, n, χ_{e}(ラグランジュの未定係数) に関する N + M + 1 個の線形方程式が得られ, 方程式系が閉じる.

計算量を減らすために Γ_{i} という量を導入する. Γ_{i} は物質 i の熱力学関数 G = m_{i}μ_{i}(T, p) を RT で規格化した値, すなわち

Γ_{i} ≡ \frac{m_{i}μ_{i}}{RT} = m_{i} [\frac{μ_{i}^{∘}}{RT} + ln p + ln \frac{m_{i}}{m} ]  (3.10)

として定義する.

Γ_{i} を用いて熱力学関数 G が極値をとる時の条件式 (3.9) 式を n_{i} に関する式に書き直す.

n_{i} = −Γ_{i} + (\frac{m_{i}}{m} ) n + ∑_{e} χ_{e}a_{ie}m_{i}.  (3.11)

この (3.11) 式を i に関して 1 から N までの和をとると,

∑_{e} B_{e}χ_{e} = ∑_{i} Γ_{i},  (3.12)

となる. また (3.11) 式の両辺に a_{if} をかけて i に関して和を取る.

∑_{e} \{∑_{i} a_{ie}a_{if}m_{i}\} χ_{e} + B_{f} (\frac{n}{m} − 1) = ∑_{i} a_{if}Γ_{i},  (3.13)

である. a_{if} の添字 f (0 ≤ f ≤ N) はダミーサフィックスで, 実質的には e と同等である.

Γ_{i} の導入によって得られた (3.12), (3.42) 式は M + 1 個の方程式から成り, それらを連立することで χ_{e} と (n/m − 1) の値が求まる. 各物質のモル数を求めるには, 得られた χ_{e} と (n/m − 1) の値を (3.11) 式に代入する必要がある.

Γ_{i} を導入した利点は, N + M + 1 個の方程式を連立すべきところを M + 1 個の方程式の連立で済ませられる点にある. この変形によって逆変換する行列の大きさが M + 1 となるので, 計算量や計算時間を短縮できるものと期待される.

この数値解法の具体的な計算サイクルは以下のようになる (図 3.2 参照).

[figure]

(1) 適当な物質存在量 m_{i} を与える. この値は物質保存則を満足するように与える.

(2) (3.42) 式の係数や Γ_{i} を計算する.

(3) (3.12), (3.42) 式を連立し, χ_{e} と n について解く. 一般的な連立一次方程式の数値解法を用いることができる.

(4) (3.11) 式より, m_{i} よりも平衡に近い物質存在量が見積もられる.

(5) (4) の操作で得られた n_{i} を m_{i} とし, (2)–(4) の操作を繰り返す.

(6) (4) の操作で得られた n_{i} が m_{i} に十分近くなった場合, n_{i} を平衡時の物質存在量とする.

[figure]

図 3.2: 大気の平衡計算の流れ図. 化学ポテンシャルは比熱や飽和蒸気圧等の既存のデータから予め計算しておく.

最終的に解くべき式をまとめる.

[figure]

温度・圧力一定, 初期状態での物質のモル数 m_{i}, 理想気体を仮定した場合, 平衡状態での物質のモル数は以下の式から得られる.

n_{i} = −Γ_{i} + (\frac{n}{m}) m_{i} + ∑_{e} χ_{e}a_{ie}m_{i}.

但し, 上式に含まれる変数 χ_{e}, n は以下の 2 式を連立させることで得られる.

∑_{e} B_{e}χ_{e} = ∑_{i} Γ_{i}.

∑_{e} \{∑_{i} a_{ie}a_{if}m_{i}\} χ_{e} + B_{f} (\frac{n}{m} − 1) = ∑_{i} a_{if}Γ_{i}.

尚, 定数 Γ_{i} は以下に定義される.

Γ_{i} ≡ m_{i} [\frac{μ_{i}^{∘}}{RT} + ln p + ln \frac{m_{i}}{m} ] .

[figure]

3.1.2 多相系への拡張

RAND 法を気相系から多相系へ拡張する (Van Zeggerren and Storey, 1970). 各相を独立に扱い, 多相系へ拡張する. 相の区別を添字 ϕ (1 ≤ ϕ ≤ Φ) で表し, 前章で導出した式を書き換える.

まず物質の保存則は以下のように書き換えることができる.

∑_{ϕ} ∑_{i} a_{ie}^{ϕ}n_{i}^{ϕ} = B_{e}  (3.32)

尚, (3.32) 式の B_{e} は各元素の全相における総和に拡張されている. 多相系で保存されているのは, 全相にわたる元素の存在量である.

次に理想系での化学ポテンシャルの表現を拡張する.

μ_{i}^{ϕ}(T, p) = μ_{i}^{∘}^{ϕ}(T) + RT ln (\frac{m_{i}^{ϕ}}{m^{ϕ}} ) + αRT ln p  (3.33)

但し, m^{ϕ} = ∑_{i} m_{i}^{ϕ} は相 ϕ に存在する物質量 (モル数) の和である. また係数 α は気相ならば 1, 液相・固相ならば 0 になる.

前章の方法を多相系に拡張した方程式に対して適用する. 適当な物質存在量 m_{i}^{ϕ} のまわりで 2 次のテーラー展開を行い, その極値を求める. (m_{i}^{ϕ} においても質量保存を満足していると仮定する.)

Q(n_{i}^{ϕ}) = G(m_{i}^{ϕ}) + ∑_{ϕ} ∑_{i} \underbrace{(\frac{∂G}{∂m_{i}^{ϕ}} )}_{A} δm_{i}^{ϕ} + \frac{1}{2} ∑_{ϕ} ∑_{i} ∑_{j} \underbrace{(\frac{∂^{2}G}{∂m_{i}^{ϕ}∂m_{l}^{ϕ}} )}_{B} δm_{i}^{ϕ}δm_{j}^{ϕ},  (3.34)

\underline{木星型惑星大気の熱力学計算第 3 章大気の平衡状態の記述 11} と書ける. ここで (3.34) 式の A 項は (3.33) 式を用いることによって,

(\frac{∂G}{∂m_{i}^{ϕ}} ) ≡ μ_{i}^{ϕ} = μ_{i}^{∘}^{ϕ}(T) + RT ln (\frac{m_{i}^{ϕ}}{m^{ϕ}} ) + αRT ln p.  (3.35)

と書ける. また (3.34) 式の B 項を (3.33) 式を用いて書き下すと,

(\frac{∂^{2}G}{∂m_{i}^{ϕ}∂m_{l}^{ϕ}} ) = RT (\frac{δ_{il}}{m_{i}^{ϕ}} − \frac{1}{m^{ϕ}} ) ,  (3.36)

となる. (3.35), (3.36) 式を (3.34) 式に代入し, RT で割ることによって,

\frac{Q(n_{i}^{ϕ})}{RT} = ∑_{ϕ} \frac{G(m_{i}^{ϕ})}{RT} + ∑_{ϕ} ∑_{i} [\frac{μ_{i}^{∘}^{ϕ}(T)}{RT} + α ln p + ln \frac{m_{i}^{ϕ}}{m^{ϕ}} ] δm_{i} + \frac{1}{2} ∑_{ϕ} [∑_{i} \frac{(δm_{i}^{ϕ})^{2}}{m_{i}^{ϕ}} − \frac{(δm^{ϕ})^{2}}{m^{ϕ}} ] ,(3.37)

が導かれる. (3.37) 式の Q(n_{i})/RT と質量保存則から, ラグランジュの未定係数法を用いて極値の条件を与える. L を以下のように定義する.

L = ∑_{ϕ} \frac{G(m_{i}^{ϕ})}{RT} + ∑_{ϕ} ∑_{i} [\frac{μ_{i}^{∘}^{ϕ}(T)}{RT} + α ln p + ln \frac{m_{i}^{ϕ}}{m^{ϕ}} ] δm_{i} + \frac{1}{2} ∑_{ϕ} [∑_{i} \frac{(δm_{i}^{ϕ})^{2}}{m_{i}^{ϕ}} − \frac{(δm^{ϕ})^{2}}{m^{ϕ}} ]
− ∑_{e} χ_{e} ∑_{ϕ} [∑_{i} a_{i,e}^{ϕ}n_{i}^{ϕ} − B_{e}^{ϕ}] .

極値をとる条件は ∂L/∂n_{i}^{ϕ} = 0 なので,

(\frac{μ_{i}^{∘}^{ϕ}(T)}{RT} + ln p + ln \frac{m_{i}^{ϕ}}{m^{ϕ}} ) + (\frac{n_{i}^{ϕ}}{m_{i}^{ϕ}} − \frac{n^{ϕ}}{m^{ϕ}} ) − ∑_{e} χ_{e}a_{ie}^{ϕ} = 0,  (3.38)

である. (3.38) 式と (3.32) 式, n^{ϕ} = ∑ n_{i}^{ϕ} を連立することで, 変数 n_{i}^{ϕ}, χ_{e} (ラグランジュの未定係数) に関する ΦN + M 個の線形方程式が得られ, 方程式系が閉じる. ここで Γ_{i}^{ϕ} という量を導入する. これは (3.10) 式の Γ_{i} を多相系に拡張したものである.

Γ_{i}^{ϕ} ≡ m_{i}^{ϕ} \frac{μ_{i}^{ϕ}}{RT} = m_{i}^{ϕ} (\frac{μ_{i}^{∘}^{ϕ}}{RT} + ln \frac{m_{i}^{ϕ}}{m^{ϕ}} + α ln p)  (3.39)

として定義する. (3.38) 式を n_{i}^{ϕ} に関して整理すると,

n_{i}^{ϕ} = −Γ_{i}^{ϕ} + (\frac{n^{ϕ}}{m^{ϕ}} ) m_{i}^{ϕ} + ∑_{e} χ_{e}a_{ie}^{ϕ}m_{i}^{ϕ}.  (3.40)

この (3.40) 式を i に関して 1 から N までの和をとると,

∑_{e} a_{i,e}^{ϕ}m_{i}^{ϕ}χ_{e} = ∑_{i} Γ_{i}^{ϕ}  (3.41)

となる. また (3.41) 式の両辺に a_{if}^{ϕ} をかけて i と ϕ に関して和を取る.

∑_{e} \{∑_{i} ∑_{ϕ} a_{ie}^{ϕ}a_{if}^{ϕ}m_{i}^{ϕ}\} χ_{e} + ∑_{ϕ} ∑_{i} a_{if}^{ϕ}m_{i}^{ϕ} (\frac{n^{ϕ}}{m^{ϕ}} ) = ∑_{i} ∑_{ϕ} a_{if}^{ϕ}Γ_{i}^{ϕ} + B_{f}  (3.42)

である. a_{if} の添字 f (0 ≤ f ≤ M) はダミーサフィックスで, 実質的には e と同等である.

(3.41), (3.42) 式は M + Φ の方程式であり, また変数の数は χ_{e} (M 個) と n^{ϕ} (Φ 個) なので, この 2 式を連立することによって方程式が閉じる.

最後に解くべき式は以下のようになる.

[figure]

温度・圧力一定, 初期状態での物質のモル数 m_{i}, 理想気体を仮定した場合, 平衡状態での物質のモル数は以下の式から得られる.

n_{i}^{ϕ} = −Γ_{i}^{ϕ} + (\frac{n^{ϕ}}{m^{ϕ}} ) m_{i}^{ϕ} + ∑_{e} χ_{e}a_{ie}^{ϕ}m_{i}^{ϕ}.

但し, 上式に含まれる変数 χ_{e}, n^{ϕ} は以下の 2 式を連立させることで得られる.

∑_{e} a_{ie}^{ϕ}m_{i}^{ϕ}χ_{e} = ∑_{i} Γ_{i}^{ϕ}

∑_{e} \{∑_{i} ∑_{ϕ} a_{ie}^{ϕ}a_{if}^{ϕ}m_{i}^{ϕ}\} χ_{e} + ∑_{ϕ} ∑_{i} a_{if}^{ϕ}m_{i}^{ϕ} (\frac{n^{ϕ}}{m^{ϕ}} ) = ∑_{i} ∑_{ϕ} a_{if}^{ϕ}Γ_{i}^{ϕ} + B_{f}

尚, 定数 Γ_{i}^{ϕ} は以下に定義される.

Γ_{i}^{ϕ} = m_{i}^{ϕ} [\frac{μ_{i}^{∘}^{ϕ}}{RT} + ln \frac{m_{i}^{ϕ}}{m^{ϕ}} + α ln p] .

[figure]

3.1.3 多相系への拡張 2 (相の扱いを物質の違いとする場合)

(3.40), (3.41), (3.42) 式では相の区別を添字 ϕ によって表現した. この方法では添字の数が多くなり繁雑であり, また空気塊に存在する物質の全部の相を考えねばならず無駄が多い. 例えば H_{2} の液相や固相は通常考える必要がない. そこで各物質の相の区別を i に押しつけることにする. つまり今までは氷と水と水蒸気は同じ物質 (i の値は同じ) で ϕ によって区別していたが, これからは氷と水と水蒸気は別々の物質 (i の値が異なる) と考える. ここでは, 気相な物質の番号を 1 ≤ i ≤ N1, 液相な物質の番号を N1 + 1 ≤ i ≤ N2, 固相な物質の番号を N2 + 1 ≤ i ≤ N とする.

まず (3.40) 式を書き換える.

1 ≤ i ≤ N1,  n_{i} = −Γ_{i} + (\frac{n_{g}}{m_{g}} ) m_{i} + ∑_{e} χ_{e}a_{ie}m_{i},  (3.43)
N1 + 1 ≤ i ≤ N2,  n_{i} = −Γ_{i} + (\frac{n_{l}}{m_{l}} ) m_{i} + ∑_{e} χ_{e}a_{ie}m_{i},  (3.44)
N2 + 1 ≤ i ≤ N,  n_{i} = −Γ_{i} + (\frac{n_{s}}{m_{s}} ) m_{i} + ∑_{e} χ_{e}a_{ie}m_{i},  (3.45)

ここで, m_{g}, m_{l}, m_{s} は以下のように定義する.

m_{g} = ∑_{i=1}^{N1} m_{i},  m_{l} = ∑_{i=N1+1}^{N2} m_{i},  m_{s} = ∑_{i=N2+1}^{N} m_{i}  (3.46)

次に (3.41) 式を書き換える. この式は ϕ 本の方程式より成っていたので, 気相, 液相, 固相の 3 本の式が必要となる.

∑_{e} {∑_{i=1}^{N1} a_{ie}m_{i}} χ_{e} = ∑_{i=1}^{N1} Γ_{i}.  (3.47)
∑_{e} { ∑_{i=N1+1}^{N2} a_{ie}m_{i}} χ_{e} = ∑_{i=N1+1}^{N2} Γ_{i}.  (3.48)
∑_{e} { ∑_{i=N2+1}^{N} a_{ie}m_{i}} χ_{e} = ∑_{i=N2+1}^{N} Γ_{i}.  (3.49)

最後に (3.42) 式を書き換える.

∑_{e} \{∑_{i=1}^{N} a_{ie}a_{if}m_{i}\} χ_{e} + ∑_{i=1}^{N1} a_{if}m_{i} (\frac{n_{g}}{m_{g}} ) + ∑_{i=N1+1}^{N2} a_{if}m_{i} (\frac{n_{l}}{m_{l}} )

+ ∑_{i=N2+1}^{N} a_{if}m_{i} (\frac{n_{c}}{m_{c}} ) = ∑_{i=1}^{N} a_{if}Γ_{i} − B_{e}.  (3.50)

定数 Γ_{i} は以下のように書ける.

1 ≤ i ≤ N1,  Γ_{i} = m_{i} [\frac{μ_{i}^{∘}}{RT} + ln \frac{m_{i}}{m_{g}} + α ln p]  (3.51)
N1 + 1 ≤ i ≤ N2,  Γ_{i} = m_{i} [\frac{μ_{i}^{∘}}{RT} + ln \frac{m_{i}}{m_{l}} + α ln p]  (3.52)
N2 + 1 ≤ i ≤ N,  Γ_{i} = m_{i} [\frac{μ_{i}^{∘}}{RT} + ln \frac{m_{i}}{m_{c}} + α ln p]  (3.53)

となる.

3.2 数値計算上のテクニック

数値計算を安定して稼働させるためにいくつかテクニックを使っている. 本節では平衡解を求めるためのテクニックを紹介する.

3.2.1 相の扱い (その 1)

ある相に存在する物質存在量が全て 0 に近付くと, (3.41) 式は

0 · χ_{e} = 0  (3.54)

となり, χ_{e} は不定となる. それを避けるために, 相に存在する全ての物質の存在量が 0 に近付いた場合, その相が存在しないものとして計算を進める.

現在は相を削除するかどうか判断するためのしきい値として 1.0 × 10^{−10} という値を用いている. このしきい値はかなりナイーブで, 例えば 1.0 × 10^{−18} にすると, 計算が不安定になる. 行列のコンディションが悪くなるためだと想像される.

3.2.2 対数項 (ln(m_{i}^{ϕ}/m^{ϕ})) の扱い

(3.39) 式の対数項を考える.

Γ_{i}^{ϕ} = m_{i}^{ϕ} (\frac{μ_{i}^{∘}^{ϕ}}{RT} + ln \frac{m_{i}^{ϕ}}{m^{ϕ}} + α ln p)

物質存在量 m_{i} が 0 に近付いても, Γ_{i}^{ϕ} が −∞ に発散しないことが必要となる. そのため実際の計算では, (3.39), (3.40), (3.41), (3.42) 式の Γ_{i}^{ϕ} に関する項を, m^{ϕ} で割った形式で置き換えている.

\frac{Γ_{i}^{ϕ}}{m^{ϕ}} = \frac{m_{i}^{ϕ}}{m^{ϕ}} (\frac{μ_{i}^{∘}^{ϕ}}{RT} + ln \frac{m_{i}^{ϕ}}{m^{ϕ}} + α ln p) .  (3.55)

(3.55) 式の対数項は

\frac{m_{i}^{ϕ}}{m^{ϕ}} · ln \frac{m_{i}^{ϕ}}{m^{ϕ}}  (3.56)

となっているので, 物質存在量 m_{i} が 0 に近付く場合には対数項は 0 収束する.

3.2.3 負のモル数の処理

(3.40), (3.41), (3.42) 式を連立することによって, より平衡状態に近い物質存在量 n_{i} が求まる. しかし得られた n_{i} は負になる可能性がある. もしも計算から得られた存在量 n_{i} が負の場合, 以下の式より得られる n_{i}^{∗} をより平衡状態に近い物質存在量とみなす.

n_{i}^{∗} = m_{i} + λ{n_{i} − m_{i}}.  (3.57)

ここで λ(> 0) は定数である. (3.57) 式の n_{i}^{∗} は元素保存の式 ((3.2) 式) を満たす. n_{i}^{∗} は必ず正なので λ は以下のようになる.

λ < \frac{m_{i}}{m_{i} − n_{i}} .  (3.58)
λ = \frac{m_{i}}{m_{i} − n_{i}} × 0.99  (3.59)

(3.58) 式の λ は全ての i に関して成立するものとする. 但し (3.58) 式の値を具体的に与える必要があるので, 実際には (3.59) 式のようにしている.

もしも全ての n_{i} が負でなければ λ = 1 である.

3.2.4 計算を停止させるための条件

本手法では逐次近似を反復的に行うために, 計算を停止させるかどうかの判定条件を設定する必要がある. 本手法では熱力学関数 G の変化率から計算を停止するかどうかの判断を行う.

• 計算前と後の熱力学関数 G を比較し, 前の値が後の値よりも小さい場合.

• 計算前と後の熱力学関数 G の変化率が 10^{−15} 以下になった場合.

以上の条件のどちらか一方が満たされた場合, 空気塊は平衡状態に達したとみなし計算を停止させる.

\underline{木星型惑星大気の熱力学計算第 4 章大気の断熱変化の記述 16}

第 4 章

大気の断熱変化の記述

4.1 大気の断熱変化の計算方法

第 3 章の議論より, 温度・圧力空間における熱力学関数 G の分布も計算できる. 熱力学関数 S は

(\frac{∂G}{∂T})_{p,n_{i}} = −S  (4.1)

なので, G の分布より温度・圧力空間における S の分布も求まる. しかし興味の対象は温度・圧力空間における S の分布ではない. 興味の対象は, 大気が任意の S の等値線に沿って変化するときの温度変化と物質存在度変化である.

任意の S の等値線に沿った大気の断熱温度減率と凝縮物質の鉛直分布を以下の方法によって計算することができる (図 4.1 参照). この方法ではある温度・圧力条件にある空気塊を出発点とする. 出発点から圧力を dp づつ変化させていき, 擬断熱条件を満たす温度・物質存在量を順次計算する.

[figure]

1. 計算開始時の温度 T = T0, 圧力 p = p0 における平衡状態, エントロピー S(T,p) を求める. 凝縮が生じた場合は凝縮物質を除去し, 空気塊のエントロピーは気相だけのエントロピーに等しいものとする.

2. 圧力を p + dp にする.

3. 温度を T から dT づつ k 回ずらした点 (T + k dT) でのエントロピー S(T + k dT, p + dp) を計算し, S(T, p) の値と一致する温度を求める.

4. 凝縮物質を除去し, 空気塊のエントロピーは気相だけのエントロピーに等しいものとする (擬断熱の条件). 気相だけのエントロピーを新たに S(T,p) として定義する.

5. 手順 2–5 を規定回数繰り返す.

[figure]

4.2 数値計算上のテクニック

数値計算を安定して稼働させるためにいくつかテクニックを使っている. 本節では擬断熱変化を求めるためのテクニックを紹介する.

4.2.1 相の扱い (その 2)

擬断熱を仮定しているので, 凝縮物質は空気塊から離脱する (存在量を 0 にする). その状態でさらに 1 ステップ進めて空気塊を p = p + dp に持ち上げると, 凝縮物質の存在量が全て 0 なので凝縮相に関する方程式が除かれる. (3.2.1 節「相の扱い (その 1)」参照). 熱力学的に凝結する高度でも, 凝縮相の方程式が除去されてしまったので, 凝縮が生じないということが起きる. そのような状態を回避するために, 圧力 p を dp 変化させるたびに, 凝縮相に凝縮するための「タネ」を入れている. 例えば水蒸気の存在量の 1/10 を液体の水の存在量に振り分けるという作業を行っている (図 4.2 参照).

4.2.2 熱力学関数 S の計算

熱力学関数 S を数値的に求めるためには, G の差分を求める必要がある. 現在, 差分形式として中央差分を採用している. 中央差分の取り方は図 4.3 に示される. 圧力 p を固定し, 温度を ±dT だけずらした点 T ± dT での熱力学関数 G をそれぞれ G(T − dT), G(T + dT) とする. そのとき温度 T, 圧力 p での G の差分は以下のようになる.

(\frac{∂G}{∂T})_{p,n_{i}} = (\frac{ΔG}{ΔT})_{p,n_{i}} ≃ (\frac{G(T + dT) − G(T − dT)}{2dT})_{p,n_{i}} .  (4.2)

この差分の誤差は (dT)^{2} 程度である.

(4.2) 式は物質量は変化しないという条件で行わねばならない. そのため温度が ±dT だけずれた点同士での各物質存在量を比較し, その変化率が差分の誤差程度に収まるようにチェックしている.

\underline{木星型惑星大気の熱力学計算第 4 章大気の断熱変化の記述 18} 図 4.1: 大気の擬断熱変化の計算手順. {図 4.2: 相の扱い (その 2) の概要. 圧力を dp だけ変化させるたびに可凝縮気体の存在量を減らし, その分 \atop 凝縮物質の存在量を増やしている.}

\overline{main.tex 1991/02/17(杉山耕一朗)}

\underline{木星型惑星大気の熱力学計算第 4 章大気の断熱変化の記述 19} 図 4.3: 中央差分の模式図. ×印の点での差分を求めるために, その前方の点と後方の点の値を用いる.

\overline{main.tex 1991/02/17(杉山耕一朗)}

\underline{木星型惑星大気の熱力学計算第 5 章化学ポテンシャルの与え方 20}

第 5 章

化学ポテンシャルの与え方

大気の状態を表現するための熱力学関数として G を選択したので, 物性値として化学ポテンシャルを与えれば良くなった. 本章では化学ポテンシャルの計算方法を議論する.

5.1 標準化学ポテンシャル

大気中の物質が全て理想気体/ 液体/ 固体/ 溶液であるとすると, その化学ポテンシャルは以下のように書ける.

μ_{i}^{ϕ}(T, p) = μ_{i}^{∘}^{ϕ}(T) + RT ln (\frac{m_{i}^{ϕ}}{m^{ϕ}} ) + αRT ln p  (5.1)

但し, μ_{i}^{∘}^{ϕ}(T) は標準化学ポテンシャル, m^{ϕ} = ∑_{i} m_{i}^{ϕ} は相 ϕ に存在する物質量 (モル数) の和である. また係数 α は気相ならば 1, 液相・固相ならば 0 になる.

(5.1) 式の右辺第 2 項, 第 3 項は空気塊の温度・圧力・物質存在量が決まれば計算できる. しかし右辺第 1 項は温度・圧力・物質存在量が決まっても, その関数の形がわからないと計算できない. そのためデータとして知らねばならないのは, 各物質の標準化学ポテンシャルとなる.

1 モル当たりの標準化学ポテンシャルは基準となる温度・圧力における熱力学関数 S, H のモル当量と比熱を用いて以下のように書ける.

μ_{i}^{∘}(T) = h(T_{0}) − T s(T_{0}, p_{0}) + ∫_{T_{0}}^{T} c_{p}dT − T ∫_{T_{0}}^{T} \frac{c_{p}}{T} dT  (5.2)

但し h は熱力学関数 H のモル当量, s は熱力学関数 S のモル当量, c_{p} はモル定圧比熱である. 尚, 全ての相において同様に書けるので, 相を区別するための添字 ϕ を外した. 基準となる温度・圧力であるが, 気体ならば標準状態 (T_{0} = 298.15 K, p_{0} = 1.0 × 10^{5}Pa), 液体・固体ならば 3 重点での温度・圧力とする.

(5.2) 式より, 標準ポテンシャルを与えるためには h(T_{0}), s(T_{0}, p_{0}), c_{p}(T) を既存のデータ集から取得すれば良いことが分かる.

5.2 気体の標準化学ポテンシャルの計算

(5.2) 式より, 標準ポテンシャルを与えるためには h(T_{0}), s(T_{0}, p_{0}), c_{p}(T) を既存のデータ集から取得すれば良いことが示された. 本節では具体的な気体の標準化学ポテンシャルの与え方を議論する.

気体の h(T_{0}), s(T_{0}, p_{0}), c_{p}(T) は以下のように与えた.

• h(T_{0})

JANAF Thermochemical Tables から与えた.

• s(T_{0}, p_{0}) JANAF Thermochemical Tables から与えた.

• c_{p}(T)

JANAF Thermochemical Tables から与えた. データとして与えられていない温度のモル定圧比熱は 3 次の Spline 補完を用いて求めた.

比熱のデータとその補完状況の例を図 5.1 に示す. 黒丸は JANAF Thermochemical Tables から引用した実測データ, 実線は実測データを 3 次の Spline 補完した値である. 但し希ガスだけは Spline 補完を用いていない. 希ガスの比熱は 20.786 J K^{−1} mol^{−1} の一定値なので Spline 補完するまでもないからである. しかも Spline 補完では一定値をうまく表現できない.

5.3 固体・液体の標準化学ポテンシャルの計算

固体・液体の化学ポテンシャルの計算も, 比熱のデータと 3 重点での h, s が分かれば計算できる. しかし液相・固相の c_{p} の実測値の表はあまり見出すことができなかった. そのため液相・固相の標準化学ポテンシャルは飽和蒸気圧, もしくは圧平衡定数を用いて計算した. 本節ではその方法を議論する.

5.3.1 飽和蒸気圧が知られている場合の計算

多くの物質に関して飽和蒸気圧の経験式が与えられている. それらは NIST Chemical Web Book, 化学工学便覧, 化学便覧から引用可能である.

1 成分 2 相系からなる系を考える. 相平衡が成立しているならば液相 (固相) と気相の化学ポテンシャルは等しいので,

μ_{l}(T, p^{∗}(T)) = μ_{g}^{∘}(T) + RT ln{p^{∗}(T)}  (5.3)

となる. ここで μ_{l}(T, p^{∗}(T)) は液相の化学ポテンシャル, μ_{g}^{∘} は気相の標準化学ポテンシャル, p^{∗} は飽和蒸気圧である. ギブスの相律より 1 成分 2 相系の自由度は 1 なので, 飽和蒸気圧は温度のみの関数として書ける. 飽和蒸気圧がわかれば (5.3) 式より, 液相 (固相) の化学ポテンシャルを計算することができる.

飽和蒸気圧の経験式として有名なのは Antone の式である. Antone の式は,

ln\{p^{∗}(T)\} = A − \frac{B}{C + T}  (5.4)

図 5.1: N_{2}, O_{2}, Ar, CO_{2} の比熱の補完状況. 黒丸は JANAF Thermochemical Tables から引用した実測データ, 実線は実測データを Spline 補完した値である. 但し Ar だけは Spline 補完を用いていない. 横軸は温度, 縦軸は c_{p} の値である.

の形式で与えられる. A, B, C は定数であり, T は温度 (K) である. モデルで使用した物質に対する A, B, C は表 5.1 に示す.

また, NH_{3} や H_{2}O に関しては, Briggs and Sackett (1989) の経験式を使うこともできる. かれらの経験式は,

ln\{p^{∗}(T)\} = \frac{a_{1}}{T} + a_{2} + a_{3} ln T + a_{4}T + a_{5}T^{2}  (5.5)

で与えられる. a_{1}, a_{2}, a_{3}, a_{4}, a_{5} は定数であり, 表 5.2 によって与えられる. 但し cgs 単位系での値なので, 使用する際には SI 単位系に変換する必要がある.

図 5.2, 図 5.3 は Antonie の経験式と Briggs and Sacket (1989) の経験式を用いて NH_{3}(l) と H_{2}O(l) の蒸気圧をプロットしたものである. 図 5.2 に示されるように, NH_{3} の飽和蒸気圧は Antonie の経験式と Briggs and Sacket (1989) の経験式とで差はほとんどない. しかし図 5.3 に示されるように, H_{2}O の飽和蒸気圧に関しては Briggs and Sacket (1989) の経験式の方がより良い近似となっている. そこで実際の計算には Briggs and Sacket (1989) の経験式を用いた.

5.3.2 化学反応の圧平衡定数が知られている場合

Lewis (1979) によって, 木星大気中に NH_{4}SH の雲が生じる可能性が示唆されている. この反応に関しては圧平衡定数が知られているので, それを用いて標準化学ポテンシャルを求める.

NH_{4}SH の生成反応は以下のように書ける.

NH_{3} + H_{2}S → NH_{4}SH  (5.6)

Briggs and Sackett (1989) によれば, (5.6) 式の圧平衡定数は以下のように書くことができる.

K_{p} = ln (p_{NH3} · p_{H_{2}S}) = 61.781 − \frac{10834}{T}  (5.7)

但し, p_{NH3}, p_{H_{2}S} はそれぞれ, アンモニア気体の分圧, 硫化水素の分圧である. 単位系は cgs 単位系である (圧力の単位は dynes/cm^{2}).

相平衡が成立するならば, NH_{4}SH の化学ポテンシャルは以下のように書ける.

μ_{NH_{4}SH,s}(T) = μ_{NH_{3},g}(T, p) + μ_{H_{2}S,g}(T, p)
= \{μ_{NH_{3},g}^{∘}(T) + RT ln (\frac{p_{NH3}}{p_{0}} )\} + \{μ_{H_{2}S,g}^{∘}(T) + RT ln (\frac{p_{H_{2}S}}{p_{0}} )\}
= μ_{NH_{3},g}^{∘}(T) + μ_{H_{2}S,g}^{∘}(T) + RT ln (\frac{p_{NH3} · p_{H_{2}S}}{p_{0}^{2}} )
= μ_{NH_{3},g}^{∘}(T) + μ_{H_{2}S,g}^{∘}(T) + RT ln (p_{NH3} · p_{H_{2}S}) − RT ln p_{0}^{2}
= μ_{NH_{3},g}^{∘}(T) + μ_{H_{2}S,g}^{∘}(T) + RT (61.781 − \frac{10834}{T}) − RT ln (1.0 × 10^{12}) (5.8)

但し, μ_{NH_{3},g}^{∘}(T) はアンモニア気体の標準化学ポテンシャル, μ_{H_{2}S,g}^{∘}(T) は硫化水素気体の標準化学ポテンシャルである. また 1 atm = 1.0 × 10^{6} dyn cm^{−2} であることを用いた. (5.8) 式より, NH_{4}SH の標準化学ポテンシャルはアンモニア気体と硫化水素気体の標準化学ポテンシャルと系の温度から求まる.

図 5.2: NH_{3}(l) の飽和蒸気圧. × は実測データ (化学便覧による). 青線は Antonie の経験式, 赤線は Briggs and Sacket (1989) の経験式.

図 5.3: H_{2}O(l) の飽和蒸気圧. × は実測データ (化学便覧による). 青線は Antonie の経験式, 赤線は Briggs and Sacket (1989) の経験式.

表 5.1: Antoine の式の定数 A, B, C の値 (化学工学便覧より)

物質名  A  B  C 適用温度 (K)
H_{2}O(l) 23.1964 3816.44 -46.13  284–441
H_{2}S(l) 20.9968 1768.69 -26.06  190–230
NH_{3}(l) 21.8409 2132.50 -32.98  179–261

表 5.2: Briggs and Sacket (1989) の経験式の定数 a_{1}, a_{2}, a_{3}, a_{4}, a_{5} の値 (Table 4 of Briggs and Sacket, 1989, ICARUS, 80, p98.)

物質名  a_{1}  a_{2}  a_{3}  a_{4}  a_{5}
NH_{3}(s)  -4122  41.67871  -1.81630  0  0
NH_{3}(l) -4409.3512 76.864252 -8.4598340  5.51029 × 10^{−3} 6.804632 × 10^{−6}
H_{2}O(s) -5631.1206 -8.363602 8.2312 -3.861449 × 10^{−2} 2.77494 × 10^{−5}
H_{2}O(l) -2313.0338 -164.03307 38.053682 -1.3844344 × 10^{−1} 7.4465367 × 10^{−5}

\underline{木星型惑星大気の熱力学計算第 6 章結果 26}

第 6 章

結果

本章ではモデル計算の結果を示す.

6.1 モデルの信頼性の検証

外惑星条件でモデルを走らせる前に, よく知られている問題を用いてモデルの信頼性を確認する. 検証内容としては, (1) 熱力学関数 G から得られた物理量がもっともらしいかの検証, (2) 地球の乾燥大気を想定した場合に得られる乾燥断熱減率がもっともらしいかの検証, を行う.

6.1.1 G を用いた他の物質量の計算

熱力学関数 G に関して以下の式が成立する.

(\frac{∂G}{∂T})_{p,n_{i}} = −S  (6.1)
(\frac{∂G}{∂p})_{T,n_{i}} = V  (6.2)

(6.1) 式, (6.2) 式から得られた S, V を, 解析的に計算して得られた S, V と比較する.

モデル大気として表 6.2 を設定した. モデル大気 A は地球大気を想定し, 小倉 (1999) から引用した. この検証ではモデル大気 A を用いて計算を行う.

S の比較を行う. 大気は理想気体であると仮定しているので, S は解析的に

s = ∫_{T_{0}}^{T} c_{p}dT + R ln p  (6.3)

である. 標準状態において (6.1) 式, (6.3) 式それぞれから得られた S を比較すると,

[figure]

{解析的に求めた値:21763.16796255409 [J] \atop G から求めた値: 21762.95779971406 [J]}

[figure]

という結果が得られた. 解析解は平衡時の物質存在量と個々の物質の比熱, 標準エントロピー, 標準エンタルピーを用いて計算した. 値としては 10^{−}4 の桁でずれる程度で, なかなか良い結果であると考える.

V の比較を行う. 理想気体であれば 1 モル当たり 22.4 l (= 2.24 × 10^{−2} m^{3}) になる. 標準状態において理想気体の体積 V と (6.2) 式から得られた V を比較すると,

[figure]

{理想気体の体積: 2.463992857238437 [m^{3}] \atop G から求めた値: 2.496777103049681 [m^{3}]}

[figure]

という結果が得られた. 解析解の方は平衡時の全物質量 (モル) に 22.4 l を乗じたものである. この結果は 10^{−2} の桁まで一致している.

表 6.1: モデル大気 A に含まれる元素とその初期存在量. これは地球の乾燥大気を念頭においている. 添字 g は気相, 添字 l は液相, 添字 s は固相を意味する. 小倉 (1999) を元に作成. 但し H_{2}O の存在量は地球大気の上限値を与えている.

物質名初期存在量
N_{2}(g)  78.088
O_{2}(g)  20.949
Ar(g)  0.93
CO_{2}(g)  0.03
CO(g) 1.0×10^{−5}
Ne(g) 1.8×10^{−3}
He(g) 5.24×10^{−4}
CH_{4}(g) 1.4×10^{−4}
Kr(g) 1.14×10^{−4}
N_{2}O(g) 5.0×10^{−5}
H_{2}(g) 5.0×10^{−5}
H_{2}O(g)  1.0

6.1.2 地球大気の乾燥断熱減率の計算

大気に含まれる物質が全て気相とした時の断熱温度減率 (乾燥断熱減率) を計算する. モデル大気 A の乾燥断熱温度減率は地球の経験から良く知られているので, モデルの動作を確認する都合良い.

モデルの設定は以下の通り.

[figure]

初期温度: 288.15 K

初期圧力: 1.0×10^{5} Pa (1 気圧)

鉛直変化 (dp): -5.0 × 10^{2} Pa 温度変化 (dT): 1.0 × 10^{−3} T 鉛直プロット数: 100 点

[figure]

モデルの出力結果を図 6.3 に示す. 図 6.3 は圧力と温度の関係をプロットしたものである. 500 hPa 変化すると 4.2 K 温度が下がることが示される.

図 6.1: モデル大気 A (地球大気, 凝縮成分無し) での断熱温度減率.

理想気体と静水圧平衡を仮定した場合, 乾燥断熱減率はおおざっぱに計算して,

\frac{∂T}{∂p} = −\frac{R_{d}T}{pc_{p}} ∼ \frac{300 · 300}{10^{5} · 10^{4}} ∼ −1.0 × 10^{−4}  (6.4)

となり, 500 hPa 下降すると 5 K 下降することになる. 但し, R_{d} は地球大気の気体定数で 300 J K^{−1} kmol^{−1}, C_{p} は地球大気の比熱で 1000 J K^{−1} kg である. この値はモデルの結果と整合的である.

6.2 地球大気の断熱温度減率と凝縮物質の鉛直分布

大気成分気体やその凝縮成分を全て理想気体/ 液体と仮定し, 本手法に則って地球大気の断熱温度減率と凝縮物質の鉛直分布を計算する. 地球を想定してモデル大気 B を設定した. これは小倉 (1999) を元に作成した. 但し H_{2}O の存在量は地球大気の上限値として与えてある. この手法の特徴は, 大気に含まれる物質の化学ポテンシャルさえ与えられれば木星大気でも地球大気でも同様に計算することができることにある.

モデルの設定は以下の通り.

[figure]

初期温度: 288.15 K

初期圧力: 1.0×10^{5} Pa (1 気圧)

鉛直変化 (dp): -5.0 × 10^{2} Pa

温度変化 (dT): 1.0 × 10^{−3} T

鉛直プロット数: 100 点

[figure]

米国標準大気モデルでは, 1013 hPa から 540 hPa までの間で温度が 32.47 K 変化する. 本手法を用いて計算した結果では, 同一区間での温度変化は約 36 K であった. 計算結果は実際の地球大気と整合的であるといえる. また計算結果によると, 凝縮が生じていない領域での温度減率は -9 K/ 100 hPa, 凝縮が生じている領域での温度減率は -6 k / 100 hPa であった. 地球大気の乾燥断熱減率と湿潤断熱減率の比は約 10 : 6 である. そのため計算結果は凝縮が生じるような領域でも正しく計算できているといえる.

表 6.2: モデル大気 B に含まれる元素とその初期存在量. これは地球の乾燥大気を念頭においている. 添字 g は気相, 添字 l は液相, 添字 s は固相を意味する. 小倉 (1999) を元に作成. 但し H_{2}O の存在量は地球大気の上限値を与えている.

物質名初期存在量
N_{2}(g)  78.088
O_{2}(g)  20.949
Ar(g)  0.93
CO_{2}(g)  0.03
CO(g) 1.0×10^{−5}
Ne(g) 1.8×10^{−3}
He(g) 5.24×10^{−4}
CH_{4}(g) 1.4×10^{−4}
Kr(g) 1.14×10^{−4}
N_{2}O(g) 5.0×10^{−5}
H_{2}(g) 5.0×10^{−5}
H_{2}O(g)  1.0
H_{2}O(l) 1.0×10^{−8}

\overline{main.tex 1991/02/17(杉山耕一朗)}

6.3 木星条件での計算

大気成分気体やその凝縮成分を全て理想気体/ 液体と仮定し, 本手法に則って木星大気の断熱温度減率と凝縮物質の鉛直分布を計算する. 木星大気を想定してモデル大気 C を設定した. この大気組成はボイジャーの観測記録から予想されるもので, Weidenschilling and Lewis (1973) が計算に用いた値を引用した.

モデルの設定は以下の通り.

[figure]

初期温度: 340.0 K, 350.0 K (2 回分計算)

初期圧力: 1.0×10^{6} Pa (10 気圧)

鉛直変化 (dp): -9.7 × 10^{3} Pa

温度変化 (dT): 1.0 × 10^{−3} T

鉛直プロット数: 100 点

[figure]

木星深部の温度構造は未だ 1 点でしか観測されていないので, 計算開始温度を 340.0 K と 350.0 K にした時の 2 回分計算した. 初期温度 340 K はガリレオプローブの観測記録をもとに与えた. 初期温度 350 K は Atreya and Romani (1985) の計算結果をもとに与えた. 本研究の計算結果を図 6.5 から図 6.8 までに示す. 但し計算結果をガリレオの観測と Atreya and Romani (1985) の結果と比較しやすくするために, 図の鉛直軸の幅を変化させている. 図 6.5, 図 6.6 は 10 bar 面高度から 0.3 bar 面高度までプロットしてある. 図 6.7, 図 6.8 は 5 bar 面高度から 0.3 bar 面高度までプロットしてある.

ガリレオプローブによって観測された木星大気の鉛直温度分布を図 6.4 に示す (Seiff et al., 1998). ガリレオプローブによって 10 bar 面高度で気温が 340 K であることが示された. ガリレオプローブの落下点は雲の少ない乾燥領域であったので, 落下点の大気組成はモデル大気 C に比べて可凝縮気体の存在量が極端に少ない. そのため凝縮の生じていない圧力領域において図 6.4 と図 6.5 を比較することは意味がある. 図 6.4 と図 6.5 を比較すると, 凝縮の生じない 10 bar 面高度から 5 bar 面高度の温度は等しい. 例えば 5 bar 面高度の温度は 275 K と等しい. したがって木星大気の乾燥断熱温度減率は正しく計算されていることが示唆される.

過去の木星大気の断熱温度減率と凝縮物質の鉛直分布を計算した代表例として Atreya and Romani (1985) の結果 (図 6.9 参照) を示す. 彼らも木星の大気組としてボイジャーの観測記録から予想されるものを用いた. 図 6.9 と図 6.8 を比較すると, 凝縮物質の種類とそのピークの存在高度は等しい. また大気上部 (1 bar 付近) の温度も等しい. したがって本手法を用いた計算結果は, 今までの木星大気断熱温度減率と物質存在量の鉛直分布に関する計算結果と整合的である.

\underline{木星型惑星大気の熱力学計算第 6 章結果 31}

図 6.2: モデル大気 B (地球大気, 凝縮成分有り) での断熱温度減率. 横軸は温度 (K), 縦軸は圧力 (Pa) 図 6.3: モデル大気 B (地球大気, 凝縮成分有り) での凝縮物質の鉛直分布. 凝縮物質は H_{2}O のみ. 横軸は温度 (K), 縦軸は圧力 (Pa)

\overline{main.tex 1991/02/17(杉山耕一朗)}

表 6.3: モデル大気 C に含まれる元素とその初期存在量. これは木星大気を念頭においている. Weidenschilling and Lewis (1973) より作成. 添字 g は気相, 添字 l は液相, 添字 s は固相を意味する.

物質名初期存在量
H_{2}(g) 8.86×10^{−1}
Ar(g) 1.12×10^{−1}
H_{2}O(g) 1.05×10^{−3}
CH_{4}(g) 6.30×10^{−4}
NH_{3}(g) 1.52×10^{−4}
H_{2}S(g) 2.90×10^{−5}
H_{2}O(l) 1.0×10^{−8}
CH_{4}(l) 1.0×10^{−8}
NH_{3}(l) 1.0×10^{−8}
H_{2}S(l) 1.0×10^{−8}
H_{2}O(s) 1.0×10^{−8}
NH_{3}(s) 1.0×10^{−8}
NH_{4}SH(s) 1.0×10^{−8}

図 6.4: ガリレオプローブによって観測された木星大気の鉛直温度分布. (Fig. 34 of Seiff et al. 1998, J. Geophys. Res., 103, p.22878)． \underline{木星型惑星大気の熱力学計算第 6 章結果 33} 図 6.5: モデル大気 C (木星大気) での断熱温度減率. 横軸は温度 (K), 縦軸は圧力 (Pa). 計算開始温度は 340 K. 実験の初期値は Seiff et al (1998) の観測結果から与えた.

\overline{main.tex 1991/02/17(杉山耕一朗)}

\overline{main.tex 1991/02/17(杉山耕一朗)}

図 6.6: モデル大気 C (木星大気) での凝縮物質の鉛直分布. 実線は H_{2}O(s), 点線は NH_{4}SH(s), 破線は NH_{3}(s). 横軸はモル分率, 縦軸は圧力 (Pa). 計算開始温度は 340 K. 実験の初期値は Seiff et al (1998) の観測結果から与えた.

\underline{木星型惑星大気の熱力学計算第 6 章結果 35} 図 6.7: モデル大気 C (木星大気) での断熱温度減率. 横軸は温度 (K), 縦軸は圧力 (Pa). 計算開始温度は 350 K. 実験の初期値は Atreya and Romani (1985) より与えた.

\overline{main.tex 1991/02/17(杉山耕一朗)}

図 6.8: モデル大気 C (木星大気) での凝縮物質の鉛直分布. 実線は H_{2}O(s), 点線は NH_{4}SH(s), 破線は NH_{3}(s). 横軸はモル分率, 縦軸は圧力 (Pa). 計算開始温度は 350 K. 実験の初期値は Atreya and Romani (1985) より与えた. 横軸が Atreya and Romani (1985) の図と異なることに注意.

\underline{木星型惑星大気の熱力学計算第 6 章結果 37} 図 6.9: 木星大気の凝縮物質の鉛直分布 (Fig.14 of Atreya and Romani 1973, In Planetary Meteorology,Cambridge Univ. Press, Cambridge, UK. 17, p.53).

\overline{main.tex 1991/02/17(杉山耕一朗)}

\underline{木星型惑星大気の熱力学計算第 7 章まとめ 38}

第 7 章

まとめ

木星型惑星大気の断熱温度減率と凝縮物質の鉛直分布を計算するための手法を検討し, 実際に計算を行った. 本手法の特徴は大気中の物質の化学ポテンシャルの表さえあれば, 「どのような組成の大気でも計算できる」ことにある. 第 6 章で示したように, 木星大気であろうと, 地球大気であろうと, 数値計算コードを変更することなく計算することができた.

本研究では木星以外の惑星大気を計算していない. それは低温 (T < 100 K) での化学ポテンシャルの表を持っていないためである. 土星・天王星・海王星では, 大気温度が 100 K 以下になるので, 現在計算できない. この温度領域の化学ポテンシャルを求めることは, 以後の課題である.

木星大気の計算結果は過去の代表的な研究 Atreya and Romani (1985) の結果にほぼ等しい. またガリレオプローブの観測にも一致する. したがって本研究によって開発された計算手法を用いることによって, 他の木星型惑星に対しても妥当な断熱温度減率と凝縮物質の鉛直分布が得られると期待される.

謝辞

本研究を行うにあたり多くの助言をして頂いた林祥介先生, 倉本圭先生に深く感謝いたします. 本論文を書くにあたり多くのフリーソフトウェアを利用しました. 特に Debian GNU/Linux を開発した Debian プロジェクトの皆様, 電脳倶楽部ライブラリ (DCL) を開発した地球流体電脳倶楽部 (http://www.gfd-dennou.org) の皆様には深く感謝致します.

\underline{木星型惑星大気の熱力学計算第 8 章参考文献 39}

第 8 章

参考文献

American Chemical Soc, American Institute of Physics, National Bureau of Standards, 1986: JANAF

Thermochemical Tables (3rd. Edition), Part I, II. American Chemical Society, Atreya, S.K., Romani. P.N., 1985: Photochemistry and clouds of Jupiter, Saturn and Uranus. In:

Hunt. G.E. (Ed.), Planetary Meteorology. Cambridge University Press, pp. 17–68. Briggs, F. H., Sackett, P. D., 1998: Radio Observations of Saturn as a Probe of its Atmosphere and

Cloud Structure. ICARUS, 80, p77–103. Lewis, J.S., 1969: Observability of Spectroscopically Active Compounds in the Atmosphere of Jupiter.

ICARUS, 10, 393–409 National Institute of Standards and Technology, 2000: NIST Chemical Web Book. http://webbook.nist.gov/chemistry/. Seiff, A., Kirk, D.B., Knight, T.C.D., Young, R.E., Mihalov, J.D., Young, L.A., Milos. F.S., Schubert,

G., Blanchard, R.C., Atkinson, D., 1998: Thermal structure of Jupiter’s atmosphere near the edge of a 5-μm hot spot in the north equatorial belt. J. Geophys. Res. 103, 22,857–22,889.

Van Zeggerren, F., Storey, S.H., 1970: The computation of chemical equilibria, Cambridge Univ. Press. Weidenschilling, S.J., Lewis, J.S., 1973: Atmospheric and cloud structure of the Jovian planets. I-

CARUS, 20, 465–476. White, W.B., Johnson, S.M., Dantzig, G.B., 1958: Chemical equilibrium in complex mixture. J. Chem.

Phys. 28, 751–755 化学工学会編, 1999: 化学工学便覧, 改訂 6 版. 丸善. 日本化学会編, 1999: 化学便覧, 改訂 4 版. 丸善. 久保亮五, 1998: 熱学・統計力学 (修正版). 裳華房. プリゴジーヌ, デフェイ, 1966: 化学熱力学 I, II. みすず書房.