主成分凝結大気における 2 次元非静力学モデルの定式化

山下達也

YAMASHITA Tatsuya 北海道大学理学院宇宙理学専攻

2012 年 3 月 27 日

\underline{2 大気主成分凝結対流支配方程式}

目次

第 1 章支配方程式 3 1.1 大気主成分の凝結を考慮した 2 次元準圧縮方程式系 . . . . . . . . . 3 1.2 乱流パラメタリゼーション . . . . . . . . . . . . . . . . . . . . . . . 6

1.2.1 乱流拡散 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2.2 地表フラックス . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.3 雲物理パラメタリゼーション . . . . . . . . . . . . . . . . . . . . . 9 1.3.1 凝結量の定式化 . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3.2 雲粒落下量の定式化 . . . . . . . . . . . . . . . . . . . . . . 11

付録 B 圧力方程式の導出 19 付録 C 主成分凝結系における乱流パラメタリゼーション 22 付録 D (??) の導出 26 付録 E 下層で飽和比一定, 上層で等温な基本場温度分布 29

\underline{大気主成分凝結対流支配方程式 3}

第1章支配方程式

第 1 章では雲対流モデル deepconv/arare4 で用いられている支配方程式について説明する.

1.1 大気主成分の凝結を考慮した 2 次元準圧縮方程式系

本モデルでは支配方程式系として準圧縮方程式系 (Klemp and Wilhelmson, 1978) に大気主成分の凝結効果を考慮したもの (Odaka et al., 2005) を用いる. 準圧縮方程式系とは弾性方程式系のうち, 基本場を導入して線形近似を行うものを指す. 弾性方程式系とは大気の圧縮性を考慮する方程式系を指す. 準圧縮方程式は更に線形近似の適用の有無により, 準圧縮方程式系と完全圧縮方程式系に分類される. 準圧縮方程式では時間積分にモード別時間分割法 (Klemp and Wilhelmson, 1978) を用いることが出来, 計算のコストを下げることが出来る^{1 )}. 一方, 完全圧縮方程式は準圧縮方程式よりも精度の点で優れているものの, モード別時間分割法を適用すると安定に計算することが出来ないことが指摘されている (MacDonald et al., 2000). そこで計算の効率性を重視して準圧縮方程式系を採用することとする.

また本モデルでは水平 1 次元, 鉛直 1 次元の 2 次元直線直交座標系を採用した. 本研究では主成分凝結を考慮した対流の性質を理解することを目的としている為, 一般に 3 次元よりも現象の解析を行いやすい 2 次元で議論を行うこととする.

以下, 主成分凝結を考慮した準圧縮方程式系について説明する. 各変数は基本場成分と擾乱成分の和で表現できると仮定する. 擾乱成分は基本場成分に比べて十分小さいと仮定し, 各物理量の時間発展の方程式における 2 次の微小項を無視する. 但し本研究では対流に着目した計算を行なうので, 移流項については 2 次の微小量を無視しない. 線形化を行なっている準圧縮方程式系では質量保存及びエネルギー保存が成立しないことに注意しなければならない^{2 )}. 準圧縮方程式系では水平流速 u, 鉛直流速 w, 温位 θ, エクスナー関数 Π, 雲密度 ρ_{s} を予報し, 圧力 p, 密度 ρ, 温度 T などの他の物理量については予報変数を用いて診断的に算出する. 温位, エクスナー関数は以下のように定義される.

1 )モード別時間分割法の詳細については離散化のドキュメントで述べる.

2 )主成分凝結を考慮した準圧縮方程式系でのエネルギー方程式の導出及びエネルギー非保存の議論については, 山下 (2009) の付録 A を参照されたい.

θ ≡ T (\frac{p_{0}}{p} )^{\frac{R}{cpv}} ,  (1.1)
Π ≡ (\frac{p}{p_{0}} )^{\frac{R}{cpv}} = \frac{T}{θ} .  (1.2)

ここで p_{0}, R, c_{pv} はそれぞれ基準圧力, 単位質量あたりの気体定数, 定圧比熱を表し, R ≃ 188.9 [J K^{−1} kg^{−1}], c_{pv} ≃ 734.1 [J K^{−1} kg^{−1}] である. 本モデルでは地表面圧力を基準圧力にとることにし, p_{0} = 700 [Pa] とする.

本モデルにおいては静止状態にある乾燥した水平一様な基本場を与える. 物理量 ϕ の基本場成分を\overline{ϕ}, 擾乱成分を ϕ′ と表すことにすると, 各予報変数は以下のように表される.

u(x, z, t) = u′(x, z, t),  (1.3)
w(x, z, t) = w′(x, z, t),  (1.4)
θ(x, z, t) = θ′(x, z, t) + \overline{θ}(z),  (1.5)
Π(x, z, t) = Π′(x, z, t) + \overline{Π}(z),  (1.6)
ρ_{s}(x, z, t) = ρ_{s}′(x, z, t).  (1.7)

流速及び雲密度の基本場成分はゼロであるので, 以後流速及び雲密度の擾乱成分については簡単の為にプライム記号 “ ′ ” を省略する. 基本場は静止状態にあるので, 静水圧平衡の式

\frac{∂\overline{p}}{∂z} = −\overline{ρ}g  (1.8)

が成り立つ. 但し g は重力加速度である. (1.8) を\overline{Π}, \overline{θ} で表現すると

\frac{∂\overline{Π}}{∂z} = −\frac{g}{c_{pv}\overline{θ}}  (1.9)

となる.

簡単の為大気は CO_{2} のみから成るものとし, 理想気体とみなせるものとする. また惑星の回転効果を無視する. このとき主成分凝結を考慮した準圧縮方程式は以下のように記述される^{3 )}.

\frac{∂u}{∂t} = −u\frac{∂u}{∂x} − w\frac{∂u}{∂z} − c_{pv}\overline{θ}\frac{∂Π′}{∂x} + D_{u} + D_{u}^{sfc},  (1.10)
\frac{∂w}{∂t} = −u\frac{∂w}{∂x} − w\frac{∂w}{∂z} − c_{pv}\overline{θ}\frac{∂Π′}{∂z} + D_{w} + g (\frac{θ′}{\overline{θ}} − \frac{R}{p_{0}} \frac{\overline{θ}}{\overline{Π}^{c_{vv}/R}} ρ_{s}) , (1.11)

3 )(1.10), (1.11) の導出については付録 A, (1.13) の導出については付録 B を参照されたい.

\overline{2012/03/27(山下達也)}

\frac{∂θ′}{∂t} = −u\frac{∂θ′}{∂x} − w\frac{∂θ′}{∂z} − w\frac{∂\overline{θ}}{∂z}
+\frac{1}{\overline{Π}}(Q_{cond} + Q_{rad} + Q_{dis}) + D_{θ} + D_{θ}^{sfc},  (1.12)
\frac{∂Π′}{∂t} = −\frac{\overline{c}^{2}}{c_{pv}\overline{ρ}\overline{θ}^{2}} [\frac{∂(\overline{ρ}\overline{θ}u)}{∂x} + \frac{∂(\overline{ρ}\overline{θ}w)}{∂z}]
+\frac{\overline{c}^{2}}{c_{pv}\overline{θ}^{2}} \frac{1}{\overline{Π}}(Q_{cond} + Q_{rad} + Q_{dis}) − \frac{\overline{c}^{2}}{c_{pv}\overline{ρ}\overline{θ}^{2}} M_{cond},  (1.13)
\frac{∂ρ_{s}}{∂t} = −\frac{∂(ρ_{s}u)}{∂x} − \frac{∂(ρ_{s}w)}{∂z} + M_{cond} + M_{fall} + D_{ρ_{s}}.  (1.14)

ここで\overline{c} は音速の基本場成分を表し, 以下の関係式により定められる.

\overline{c}^{2} = \frac{c_{pv}}{c_{vv}} R\overline{T} = \frac{c_{pv}}{c_{vv}} R\overline{θ}\overline{Π}.  (1.15)

但し c_{vv}, c_{pv} はそれぞれ CO_{2} 気体の定積比熱, 定圧比熱を表し, c_{vv} ≃ 545.2 [J K^{−1} kg^{−1}], c_{pv} ≃ 734.1 [J K^{−1} kg^{−1}], である. D_{u}, D_{w}, D_{θ}, D_{ρ_{s}} はそれぞれ水平流速, 鉛直流速, 温位, 雲密度に関する乱流拡散項である. D_{u}^{sfc}, D_{θ}^{sfc} はそれぞれ地表面からの運動量, 熱のフラックスの輸送を表す項である. Q_{cond}, Q_{rad}, Q_{dis} はそれぞれ温度に対する潜熱加熱項, 放射加熱項, 散逸加熱項である. また M_{cond}, M_{fall} はそれぞれ単位時間単位体積当たりの凝結量, 単位時間単位体積当たりの雲粒の落下量である. D_{u}, D_{w}, D_{θ}, D_{ρ_{s}} の具体的な表式については 1.2 節で述べる. また M_{cond}, M_{fall} の具体的な表式については 1.3 節で述べる.

熱力学の式 (1.12), 圧力方程式 (1.13) の中の非断熱加熱項 Q_{cond}, Q_{rad}, Q_{dis} について述べる. 潜熱加熱項については

Q_{cond} = \frac{LM_{cond}}{c_{pv}\overline{ρ}}  (1.16)

と表す. ここで L は単位質量当たりの CO_{2} の潜熱であり, 本モデルでは一定とみなして L ≃ 5.86 × 10^{5} [J kg^{−1}] とする. Q_{dis} の具体的な表式については 1.2 節で述べる. Q_{rad} は陽に時間発展を解くことはせずに, 実験設定に応じて高度と時間に依存するパラメータとして与える.

圧力方程式 (1.13) について注意すべき点をいくつか述べる. (1.13) の右辺第 2 項は非断熱加熱に伴う熱膨張の効果, (1.13) の右辺第 3 項は凝結による大気質量の減少に伴う圧力低下の効果を表している. Klemp and Wilhelmson(1978) では, (1.13) の右辺第 2 項及び第 3 項は十分小さいものとして無視しているが, 本モデルでは無視しない. (1.16) を用いると, (1.13) は以下のように書き換えられる^{4 )}.

\frac{∂Π′}{∂t} = −\frac{\overline{c}^{2}}{c_{pv}\overline{ρ}\overline{θ}^{2}} [\frac{∂(\overline{ρ}\overline{θ}u)}{∂x} + \frac{∂(\overline{ρ}\overline{θ}w)}{∂z}]

4 )(1.17) の導出については付録 B を参照されたい.

+\frac{\overline{c}^{2}}{c_{pv}\overline{ρ}\overline{θ}} (\frac{L}{c_{pv}\overline{T}} − 1) M_{cond} + \frac{\overline{c}^{2}}{c_{pv}\overline{θ}^{2}\overline{Π}}(Q_{rad} + Q_{dis}).  (1.17)

火星極域の場合 O(\overline{T}) ∼ 150 [K] であるので, c_{pv} ≃ 734.1 [J K^{−1} kg^{−1}], L ≃ 5.86 × 10^{5} [J kg^{−1}] より L/(c_{pv}\overline{T}) > 1 が成り立つ. 即ち火星極域で CO_{2} の凝結が生じた場合, 圧力は増加することになる.

以下, 境界条件について述べる. 水平方向の境界では周期境界条件を適用し, 鉛直方向の境界では応力無し条件を適用する. 計算領域の水平方向の長さを L とすると, 任意の物理量 ϕ に対して課されるべき水平方向の境界条件は

ϕ(x, z, t) = ϕ(x + L, z, t)  (1.18)

と表される. また計算領域の鉛直方向の長さを D, w 以外の予報変数を ψ とすると, 鉛直方向の境界条件は

w = 0,  (z = 0, D)  (1.19)
\frac{∂ψ}{∂z} = 0,  (z = 0, D)  (1.20)

と表される.

1.2 乱流パラメタリゼーション

一般に大気の運動を記述する場合, 乱流運動 (小規模運動) も考慮する必要がある. しかし数値モデルの場合有限個の格子点を配置して計算を行う為, 格子点間隔よりも小さなスケール (サブグリッドスケール) の現象を取り扱うことが出来ない. そこでサブグリッドスケールの乱流運動に関する物理量を格子点上の物理量と結びつけて表現する (乱流パラメタリゼーション).

1.2.1 乱流拡散

本モデルでは Klemp and Wilhelmson(1978), 坪木・榊原 (2001) に従って乱流拡散項, 乱流拡散係数, 散逸加熱項を表現するのに 1.5 次のクロージャモデルを用いる. 1.5 次のクロージャモデルを用いると D_{u}, D_{w}, D_{θ}, D_{ρ_{s}} は以下のように表現される.

D_{u} = −\frac{∂}{∂x} [−2K_{m} \frac{∂u}{∂x} + \frac{2}{3}E_{turb}] − \frac{∂}{∂z} [−K_{m} (\frac{∂u}{∂x} + \frac{∂w}{∂z})] , (1.21)
D_{w} = −\frac{∂}{∂x} [−K_{m} (\frac{∂w}{∂x} + \frac{∂u}{∂z})] − \frac{∂}{∂z} [−2K_{m} \frac{∂w}{∂z} + \frac{2}{3}E_{turb}] , (1.22)

D_{θ} = −\frac{∂}{∂x} [K_{h} \frac{∂θ′}{∂x} ] − \frac{∂}{∂z} [K_{h} \frac{∂(θ′ + \overline{θ})}{∂z}] ,  (1.23)
D_{ρ_{s}} = −\frac{∂}{∂x} [K_{h} \frac{∂ρ_{s}}{∂x} ] − \frac{∂}{∂z} [K_{h} \frac{∂ρ_{s}}{∂z} ]  (1.24)

但し K_{m}, K_{h} はそれぞれ運動量, 熱に対する乱流拡散係数であり, E_{turb} はサブグリッドスケールの乱流運動エネルギーである. 次元解析に基づき, E_{turb} を K_{m} で表現すると

E_{turb} = (\frac{K_{m}}{C_{m}l})^{2}  (1.25)

となる. ここで l は混合距離であり, 水平方向の格子間隔 Δx 及び鉛直方向の格子間隔 Δz を用いて l = (ΔxΔz)^{1/2} と表現する. また C_{m} は無次元定数である. K_{h} については Deardorff(1972) に従って

K_{h} = 3K_{m}  (1.26)

と与える. 1.5 次のクロージャモデルでは以下に示す K_{m} の時間発展を解き, 各乱流拡散項を定める^{5 )}.

\frac{∂K_{m}}{∂t} = − (u\frac{∂K_{m}}{∂x} + w\frac{∂K_{m}}{∂z} ) − \frac{3gC_{m}^{2}l^{2}}{2} (\frac{1}{\overline{θ}} \frac{∂θ}{∂z} − \frac{R\overline{θ}}{p_{0}\overline{Π}^{c_{v}/R}} \frac{∂ρ_{s}}{∂z} )

+C_{m}^{2}l^{2} [(\frac{∂u}{∂x})^{2} + (\frac{∂w}{∂z})^{2}] + \frac{C_{m}^{2}l^{2}}{2} (\frac{∂u}{∂z} + \frac{∂w}{∂x})^{2}

−\frac{K_{m}}{3} (\frac{∂u}{∂x} + \frac{∂w}{∂z}) + \frac{1}{2} (\frac{∂^{2}K_{m}^{2}}{∂x^{2}} + \frac{∂^{2}K_{m}^{2}}{∂z^{2}} )

+ (\frac{∂K_{m}}{∂x} )^{2} + (\frac{∂K_{m}}{∂z} )^{2} − \frac{C_{ε}}{2C_{m}l^{2}} K_{m}^{2}.  (1.27)

ここで C_{ε} は無次元定数である. 無次元定数の値については Deardorff(1975) に従って C_{m} = C_{ε} = 0.2 とする.

散逸加熱項 Q_{dis} については

Q_{dis} = \frac{C_{ε}}{c_{pv}l}E_{turb}^{3/2} = \frac{C_{ε}}{c_{pv}l} (\frac{K_{m}}{C_{m}l})^{3}  (1.28)

と与える.

5 )式 (1.27) の導出については, 付録 C を参照されたい.

1.2.2 地表フラックス

地表面からの運動量, 熱のフラックスをそれぞれ F_{u}, F_{θ} とすると, 地表面からの運動量, 熱の輸送量 D_{u}^{sfc}, D_{θ}^{sfc} はそれぞれ

D_{u}^{sfc} = −\frac{1}{ρ} \frac{∂F_{u}}{∂z} ,  (1.29)
D_{θ}^{sfc} = −\frac{1}{Πρc_{p}} \frac{∂F_{θ}}{∂z}  (1.30)

と表される. 本モデルでは, 地表面からの運動量, 熱のフラックスの表現として, Louis(1979) を簡単化した小高他 (2001) によるバルク法の定式化を採用している. 以下では, 火星大気に対するバルク法の表式を示す.

小高他 (2001) では D_{u}^{sfc}, D_{θ}^{sfc} を以下のように表す.

F_{u} = −C_{D}V_{sfc}ρu_{z_{1}},  (1.31)
F_{θ} = −c_{p}C_{D}V_{sfc}ρ(T_{z_{1}} − T_{sfc}).  (1.32)

ここで C_{D} はバルク係数, z_{1} はモデル最下層高度, T_{sfc} は地表面温度, V_{sfc} は地表面風速である. 地表面におけるエクスナー関数, 密度がモデル最下層高度での値と等しいとし,

D_{u}^{sfc} ≈ \frac{∂}{∂z} (C_{D}V_{sfc}u_{z_{1}}),  (1.33)
D_{θ}^{sfc} ≈ \frac{∂}{∂z} [C_{D}V_{sfc} (θ_{z_{1}} − \frac{T_{sfc}}{Π_{z_{1}}} )]  (1.34)

と近似する. V_{sfc} は

V_{sfc} = \sqrt{u_{z_{1}}^{2} + v_{0}^{2}}  (1.35)

で表される. v_{0} は風が吹いていなくても日射による加熱や蒸発といった現実的な現象が起きることを保証する為の風速の最低値であり, 本モデルでは v_{0} = 0 m/s としている. C_{D} はバルクリチャードソン数 Ri_{B} の関数であり,

C_{D} = \begin{cases}C_{DN} (1 − \frac{A_{1}Ri_{B}}{1+A_{3}|Ri_{B}|^{1/2}} ) , (Ri_{B} < 0) \\ C_{DN} \frac{1}{1+A_{2}Ri_{B}} , (Ri_{B} > 0) \end{cases}  (1.36)

Ri_{B} = \frac{gz_{1}(θ′_{1} + \overline{θ}_{1} − T_{sfc})}{\overline{θ_{1}}^{x}max(u_{1}^{2}, u_{min}^{2})}  (1.37)

と表される. ここで u_{min} はバルクリチャードソン数に対する風速の最低値であり, u_{min} = 0.1 m/s としている. また C_{DN} は中立成層でのバルク係数, A_{1}, A_{2}, A_{3} は定数係数であり,

C_{DN} = [\frac{κ}{ln(z_{1}/z_{0})}]^{2} ,
A_{1} = 9.4,
A_{2} = 4.7,
A_{3} = 7.4C_{DN}A_{1} (\frac{z_{1}}{z_{0}} )^{1/2}  (1.38)

である. ここで κ = 0.35 はカルマン定数, z_{0} は粗度長であり, Sutton et al.(1978) に従い, 10^{−2} m とする.

1.3 雲物理パラメタリゼーション

1.3 節では単位時間単位体積当たりの凝結量 M_{cond} 及び単位時間単位体積当たりの雲粒の落下量 M_{fall} の表し方について述べる.

地球大気の雲対流モデルでは, 過飽和状態は瞬時に解消されるものとみなし, 飽和状態となるように凝結量を調整することが多い (例えば Soong and Ogura, 1973 ). しかし火星大気の場合, 過飽和状態が有限の時間存在しうることが観測・実験・理論によって示唆されている (例えば Colaprete et al., 2003; Glandorf et al., 2002; Määttänen et al., 2005). そこで本モデルでは雲粒が有限の時間をかけて成長するものとして定式化を行なう.

火星大気での温度・圧力環境下では, CO_{2} は気相或いは固相のいずれかの状態で存在する. 従って本モデルで考えるべき雲物理過程は CO_{2} 大気と CO_{2} 氷雲の変換のみである. 地球大気の場合, 氷晶核の形成のメカニズムは大気中に存在する微粒子の影響を受ける不均一核形成と, その影響を受けない均一核形成の 2 種類に分類されることが知られている (例えば水野 (2000)). 地球上の雲や霧などの大部分は不均一核形成によって生じると考えられている. 火星大気中にはダスト粒子や H_{2}O 氷でコーティングされたダストが浮遊しており, 不均一核形成が生じると考えられる (Gooding, 1986). そこで本モデルでは核形成として不均一核形成のみを考える.

雲粒の成長過程は拡散成長と併合成長の 2 種類に分類される. 拡散成長とは気相中の分子が拡散によって雲粒表面に移動・凝結することによって雲粒が成長することを指す. 併合成長とは雲粒同士が衝突し, 合体することによって雲粒が成長することを指す. Tobie et al.(2003) によると, 凝結核となるダストの平均半径が約 0.5 [μm] であり (Ockert-Bell et al., 1997), 可視光に対する典型的な光学的厚さがおおよそ 0.1 となる火星大気の場合, 半径 50 [μm] 程度の雲粒が併合成長で出来るのに要する時間は拡散成長で出来るのに要する時間よりも十分長い. 北守 (2006) は流れの存在しない火星大気での拡散成長についての数値計算を行ない, 拡散成長の時間スケールが 1 – 20 [s] 程度であること, 雲粒半径が 5 – 25 [μm] となることを見出した. 雲粒半径が小さいほど拡散成長は併合成長に比べて優勢となるので, 火星大気環境では拡散成長が優勢となると考えられる. そこで本モデルでは併合成長による雲粒成長は無視できると仮定し, 拡散成長のみを考える.

1.3.1 凝結量の定式化

以下, 単位時間体積当たりの凝結量 M_{cond} の取り扱いについて述べる. 本モデルでは単位質量の気相に含まれる凝結核の個数及び半径は空間的・時間的に一様と仮定する. また雲粒の半径は各格子内において空間的に一定であると仮定する. 更に雲粒は球形の凝結核を核として形成され, 雲粒自身も球形となると仮定する. このとき, 単位時間単位体積当たりの凝結量 M_{cond} は Tobie et al. (2003) 同様, 以下のように表される^{6 )}.

M_{cond} = \frac{4πr_{d}\overline{ρ}N^{∗}}{R_{h}} (S − 1).  (1.39)

R_{h}, S, r_{d} はそれぞれ熱輸送に関する係数, 飽和比, 雲粒半径であり,

R_{h} = \frac{L^{2}}{kRT^{2}} = \frac{L^{2}}{kRθ^{2}Π^{2}} ,  (1.40)
S = \frac{p}{p_{∗}} = \frac{p_{0}^{\frac{cpv}{R}} Π}{p_{∗}}  (1.41)
r_{d} = (r_{aero}^{3} + \frac{3ρ_{s}}{4ρ_{I}π\overline{ρ}N^{∗}} )^{1/3}  (1.42)

と表される. 但し k, p_{∗} はそれぞれ熱拡散係数, CO_{2} の飽和蒸気圧, ρ_{I} は CO_{2} 氷の密度, r_{d} は雲粒半径, r_{aero} は凝結核の半径, N^{∗} は単位質量の気相に含まれる凝結核の数である. 本モデルでは ρ_{I} = 1.565 × 10^{3} [kg/m^{3}] と与え, r_{aero}, N^{∗} は実験に応じて与える. 図 1.1 に雲密度の関数として雲粒半径をプロットしたものを示す. CO_{2} の飽和蒸気圧については半経験式である Antoine の式 (Antoine, 1888) に基づき,

p_{∗} ≃ exp (A_{ant} − \frac{B_{ant}}{T} )  (1.43)

とする. ここで A_{ant}, B_{ant} は実験により定まる係数であり, CO_{2} の場合 A_{ant} = 27.4, B_{ant} = 3103 である (化学工学会, 1999).

6 )(1.39) の導出については付録 D を参照されたい.

\overline{2012/03/27(山下達也)}

1.3.2 雲粒落下量の定式化

以下, 単位時間体積当たりの雲粒落下量 M_{fall} の取り扱いについて述べる. M_{fall} は Kessler(1969) と同様に, 雲粒の終端速度 V_{term} での移流として表現する. 即ち

M_{fall} = \frac{∂}{∂z} (ρ_{s}V_{term})  (1.44)

と表す. 火星の大気密度は小さく, 微小な雲粒・ダスト粒子の終端速度は Stokes 則で評価される値よりも有意に大きくなることが指摘されている. そこで雲粒の終端速度については, Conrath(1975) 以来多くの火星大気研究で行なわれているのと同様に, Cunningham 補正を Stokes 則に適用して

V_{term} = C_{sc} \frac{2r_{d}^{2}gρ_{I}}{9η}  (1.45)

と表す (Stokes, 1845; Cunningham, 1910). η は粘性係数であり, Sutherland の公式に従い,

η = η_{ref} (\frac{T_{ref} + C_{CO_{2}}}{\overline{T} + C_{CO_{2}}} ) (\frac{\overline{T}}{T_{ref}} )^{3/2} = η_{ref} (\frac{T_{ref} + C_{CO_{2}}}{\overline{θ}\overline{Π} + C_{CO_{2}}} ) (\frac{\overline{θ}\overline{Π}}{T_{ref}} )^{3/2} (1.46)

で表現する (Sutherland, 1893). η_{ref}, T_{ref}, C_{CO_{2}} はそれぞれ粘性係数の基準値, 温度の基準値, CO_{2} に関する Sutherland 定数であり, η_{ref} = 1.47 × 10^{−5} [Pa · s], T_{ref} = 293 [K], C_{CO_{2}} = 240 [K] と与える (国立天文台, 2004). C_{sc} は微小な粒子における Stokes 則からのずれを補正する係数 (Cunningham 補正係数) であり,

C_{sc} = 1 + \frac{4}{3}K_{n}  (1.47)

と与える (Rossow, 1978: Tobie et al., 2003). 但し K_{n} ≡ λ/r_{d} は雲粒に対する Knudsen 数である^{7 )}. λ は CO_{2} の平均自由行程であり,

λ = \frac{k_{B}\overline{T}}{\sqrt{2}πσ^{2}\overline{p}} = \frac{k_{B}\overline{θ}}{\sqrt{2}πσ^{2}p_{0}\overline{Π}^{c_{vv}/R}}  (1.48)

と与える (Chapman and Cowling, 1970). k_{B} は Boltzmann 定数, σ は CO_{2} 分子の直径であり, k_{B} = 1.38 × 10^{−23} [m^{2} kg s^{−2} K^{−1}], σ = 3.3 × 10^{−10} [m] である (Golden and Sircar, 1994). 図 1.2 に火星極冠周縁の凝結高度付近 (T ∼ 150[K], p ∼ 400[Pa]), 湿潤断熱層上端 (T ∼ 135[K], p ∼ 90[Pa]) での雲粒の終端速度を示す.

C_{sc} = 1 + K_{n} [1.257 + 0.400 exp (−\frac{1.10}{K_{n}})] .

7 )Cunnigham(1910) によって初めて与えられた補正係数の表式は

C_{sc} = 1 + 1.255K_{n} である. その後より精度の良い Cunningham 補正の表式がいくつか提案された. 例えば Davies(1945)

では以下の表式が提案された.

図 1.1: 雲粒半径を雲密度の関数としてプロットしたもの. 実線は火星極冠周縁の凝結高度付近の特徴的な温度圧力条件 (T = 150[K], p = 400[Pa]) を想定したもの, 破線は火星極冠周縁の湿潤断熱層上端の特徴的な温度圧力条件 (T = 135[K], p = 90[Pa]) を想定したものである. 凝結核数密度は 5.0 × 10^{8}[kg^{−1}] としている.

図 1.2: 雲粒の終端速度を雲粒半径の関数としてプロットしたもの. 実線は火星極冠周縁の凝結高度付近の特徴的な温度圧力条件 (T = 150[K], p = 400[Pa]) を想定したもの, 破線は火星極冠周縁の湿潤断熱層上端の特徴的な温度圧力条件 (T = 135[K], p = 90[Pa]) を想定したものである.

参考文献

Antoine, C., 1888: Tensions des vapeurs: nouvelle relation entre les tensions et les températures. Les Comptes Rendus de l′ Acadèmie des Sciences, 107, 681 – 684

Colaprete, A., Haberle, R. M., Toon, O. B., 2003: Formation of convective carbon dioxide clouds near the south pole of Mars. Journal of Geophysical Research, 108(E7), 5081, doi:10.1029/2003JE002053

Chapman, S., Cowling, T. G., 1970: The Mathematical theory of non-uniform gases. Cambridge University Press, 423 pp

Conrath, B. J., 1975: Thermal structure of the Martian atmosphere during the dissipation of the dust storm of 1971. Icarus, 24, 36 – 46

Cunningham, E., 1910: On the velocity of steady fall of spherical particles through fluid medium. Proceedings of the Royal Society of London, Series A, 83, 357 – 365

Davies, C. N., 1945: Definitive equations for the fluid resistance of spheres. The Proceedings of the Physical Society, 57, 259 – 270

Deardorff, J. W., 1972: Numerical investigation of neutral and unstable planetary boundary layers. Journal of the Atmospheric Sciences, 29, 91 – 115

Deardorff, J. W., 1975: The development of boundary-layer turbulence models for use in studying the severe storm environment. Proceeding of the SESAME meeting, Boulder, NOAA-ERL, 251 – 261

Gerber, H., 1991 : Supersaturation and droplet spectral evolution in fog. Journal of the Atmospheric Sciences, 45, 1651 – 1664

Glandorf, D. L., Colaprete, A., Tolbert, M. A., Toon, O. B., 2002: CO_{2} snow on Mars and early Earth : Experimental constraints. Icarus, 160, 66 – 72

Golden, T. C., Sircar, S., 1994: Gas adsorption on silicate. Journal of Colloid and Interface Science, 162, 182 – 188 Gooding, J. L., 1986: Martian dust particles as condensation cuclei — A preliminary assessment of mineralogical factors. Icarus, 66, 56 – 74

Kessler, E., 1969: On the distribution and continuity of water substance in atmospheric circulation. Meteorological Monographs, 10, 1 – 84

Klemp, J. B., Wilhelmson, R. B., 1978: The simulation of three-dimensional convective storm dynamics. Journal of the Atmospheric Sciences, 35, 1070 – 1096

Louis, J. 1979: A parametric model of vertical eddy fluxes in the atmosphere. Bound.- Layer Meteor., 17, 187 – 202

MacDonald, A. E., Lee, J. L., Xie, Y., 2000: The use of quasi-nonhydrostatic models for mesoscale weather prediction. Journal of the Atmospheric Sciences, 57, 2493 – 2517

Määttänen, A., Vehkamäki, H., Lauri, A., Merikallio, S., Kauhanen, J., Savijärvi, H., Kulmala, M., 2005: Nucleation studies in the Martian atmosphere, Journal of Geophysical Research, 110, E02002, doi:10.1029/2004JE002308

Ockert-Bell, M. B., Bell, III, J. F., Pollack, J. B., McKay, C. P., Forget, F., 1997: Absorption and scattering properties of the martian dust in the solar wavelength. Journal of Geophysical Research, 102, 9039 – 9050

Odaka, M., Kitamori, T., Sugiyama, K., Nakajima, K., Takahashi, Y. O., Ishiwatari M., Hayashi, Y.-Y., 2005: A formulation of non-hydrostatic model for moist convection in the Martian atmosphere. Proceeding of the 38 th ISAS Lunar and Planetary Symposium, The InInstitute of Space and Astronautical Science., 173 – 175

Rossow, W. B., 1978: Cloud microphysics: Analysis of the clouds of Earth, Venus Mars, and Jupiter. Icarus, 36, 1 – 50

Soong, S-T., Ogura, Y., 1973: A comparison between axi-symmetric and slab-symmetric cumulus cloud models. Journal of the Atmospheric Sciences, 30, 879 – 893

Stokes, G. G., 1845: On the effect of internal friction on the motion of pendulum, Transactions of the Cambridge Philosophical Society, 9, part2, 8 – 14

Sutherland, W., 1893: The viscosity of gases and molecular force. Philosophical Magazine, S. 5, 36, 507 – 531

Sutton, J. L., C. B. Leovy, and J. E. Tillman, 1978: Diurnal variations of the Martian surface layer meteorological parameters during the first 45 sols at two Viking lander sites, J. Atmos. Sci., 35, 2346 – 2355 Tobie, G., Forget, F., Lott, F., 2003: Numerical simulation of winter polar wave clouds observed by Mars Global Surveyor Mars Orbiter Laser Altimeter. Icarus, 164, 33 – 49

小高正嗣, 中島健介, 石渡正樹, 林祥介, 2001: 2 次元非弾性系を用いた火星大気放射対流の数値計算, ながれ, 20, ながれマルチメディア 2001, http://www.nagare.or.jp/mm/2001/odaka/

化学工学会編, 1999: 化学工学便覧, 丸善, 1339 pp

北守太一, 2006: 大気主成分の凝結を考慮した火星大気非静力学モデルの開発. 北海道大学大学院地球惑星科学専攻修士論文, 99 pp http://www.gfd-dennou.org/arch/prepri/2006/hokudai/kitamo/paper/pub/main.pdf

国立天文台編, 2004: 理科年表, 丸善, 1015 pp

坪木和久, 榊原篤志, 2001: CReSS ユーザーズガイド第 2 版. 253 pp http://cf.tokyo.rist.or.jp/CReSS.top.html

水野量, 2000: 雲と雨の気象学. 朝倉書店, 196 pp

山下達也, 2009: 大気主成分の凝結を考慮した 2 次元湿潤対流— 火星大気中の CO_{2} 氷雲を伴う対流に関する基礎的数値実験 — . 北海道大学大学院宇宙理学専攻修士論文, 179 pp \underline{16 大気主成分凝結対流運動方程式の導出}

付録A 運動方程式の導出

本付録では運動方程式 (1.10), (1.11) の導出を行なう. 固相密度 ρ_{s} が気相密度 ρ_{v} に比べて十分小さい場合, 一般的な圧縮性流体の方程式は

\frac{Du}{Dt} = −\frac{1}{ρ} \frac{∂p}{∂x} + D_{u},  (A.1)
\frac{Dw}{Dt} = −\frac{1}{ρ} \frac{∂p}{∂z} − g + D_{w}  (A.2)

と表される. 但し ρ ≡ ρ_{v} + ρ_{s} は全密度である. 圧力傾度力を p の代わりに θ, Π で表すと,

−\frac{1}{ρ} \frac{∂p}{∂x_{i}} = −\frac{ρ_{v}}{ρ} \frac{1}{ρ_{v}} \frac{∂p}{∂x_{i}}
= −\frac{ρ_{v}}{ρ} \frac{RθΠ}{p} \frac{∂}{∂x_{i}} [p_{0}Π^{c_{pv}/R}]
= −\frac{ρ_{v}}{ρ} \frac{RθΠ}{p} \frac{p_{0}c_{pv}}{R} Π^{(c_{pv}/R)−1} \frac{∂Π}{∂x_{i}}
= −\frac{ρ_{v}}{ρ} \frac{RθΠ}{p}p_{0}Π^{c_{pv}/R} \frac{c_{pv}}{RΠ} \frac{∂Π}{∂x_{i}}
= −\frac{ρ_{v}}{ρ} \frac{RθΠ}{p}p\frac{c_{pv}}{RΠ} \frac{∂Π}{∂x_{i}}
= −\frac{ρ_{v}}{ρ} c_{pv}θ\frac{∂Π}{∂x_{i}}  (A.3)

となる. 但し式変形の途中で理想気体の状態方程式 p = ρ_{v}RT = ρ_{v}RθΠ を用いた. (A.3) を (A.1), (A.2) に代入すると,

\frac{Du}{Dt} = −\frac{ρ_{v}}{ρ} c_{pv}θ\frac{∂Π}{∂x} + D_{u},  (A.4)
\frac{Dw}{Dt} = −\frac{ρ_{v}}{ρ} c_{pv}θ\frac{∂Π}{∂z} − g + D_{w}  (A.5)

となる.

次に (A.4) を線形化する. 基本場には雲は存在しないと想定しているので, \overline{ρ} = \overline{ρ_{v}}, ρ′ = ρ′_{v} + ρ′_{s}, ρ = \overline{ρ_{v}} + ρ′_{v} + ρ′_{s} となる. (A.4) の各変数を基本場成分と擾乱成分の和 \underline{大気主成分凝結対流運動方程式の導出 17} で表し, 擾乱成分の 2 次以上の量は十分小さいものとして無視すると,

\frac{∂u}{∂t} + u\frac{∂u}{∂x} + w\frac{∂u}{∂z} = −\frac{\overline{ρ_{v}} + ρ′_{v}}{\overline{ρ} + ρ′} c_{pv}(\overline{θ} + θ′)\frac{∂Π′}{∂x} + D_{u}
= −\frac{\overline{ρ_{v}} + ρ′_{v}}{\overline{ρ}} \frac{1}{1 + (ρ′/\overline{ρ})}c_{pv}(\overline{θ} + θ′)\frac{∂Π′}{∂x} + D_{u}
≃ −\frac{\overline{ρ_{v}} + ρ′_{v}}{\overline{ρ}} (1 − \frac{ρ′}{\overline{ρ}} ) c_{pv}(\overline{θ} + θ′)\frac{∂Π′}{∂x} + D_{u}
= −\frac{\overline{ρ_{v}} + ρ′_{v}}{\overline{ρ_{v}}} (1 − \frac{ρ′}{\overline{ρ}} ) c_{pv}(\overline{θ} + θ′)\frac{∂Π′}{∂x} + D_{u}
= − (1 + \frac{ρ′_{v}}{\overline{ρ_{v}}} ) (1 − \frac{ρ′}{\overline{ρ}} ) c_{pv}(\overline{θ} + θ′)\frac{∂Π′}{∂x} + D_{u}
≃ −c_{pv}\overline{θ}\frac{∂Π′}{∂x} + D_{u},
\frac{∂u}{∂t} = −u\frac{∂u}{∂x} − w\frac{∂u}{∂z} − c_{pv}\overline{θ}\frac{∂Π′}{∂x} + D_{u}  (A.6)

となり, (1.10) が得られる. 次に (A.5) を線形化する. 基本場において静水圧平衡が成り立つとき, (A.5) より

0 = −c_{pv}\overline{θ}\frac{∂\overline{Π}}{∂z} − g  (A.7)

となる. (A.5) を線形化し, (A.7) を用いて書き換えると,

\frac{∂w}{∂t} + u\frac{∂w}{∂x} + w\frac{∂w}{∂z} = −\frac{\overline{ρ_{v}} + ρ′_{v}}{\overline{ρ} + ρ′} c_{pv}(\overline{θ} + θ′) (\frac{∂\overline{Π}}{∂z} + \frac{∂Π′}{∂z} ) − g + D_{w}
≃ − (1 + \frac{ρ′_{v}}{\overline{ρ_{v}}} ) (1 − \frac{ρ′}{\overline{ρ}} ) c_{pv}(\overline{θ} + θ′) (\frac{∂\overline{Π}}{∂z} + \frac{∂Π′}{∂z} )

−g + D_{w}

≃ −c_{pv}\overline{θ}\frac{∂\overline{Π}}{∂z} − c_{pv}\overline{θ}\frac{∂Π′}{∂z}

−\frac{ρ′_{v}}{\overline{ρ_{v}}} c_{pv}\overline{θ}\frac{∂\overline{Π}}{∂z} + \frac{ρ′}{\overline{ρ}} c_{pv}\overline{θ}\frac{∂\overline{Π}}{∂z} − c_{pv}θ′ \frac{∂\overline{Π}}{∂z} − g + D_{w}

= −c_{pv}\overline{θ}\frac{∂Π′}{∂z} − \frac{ρ′_{v}}{\overline{ρ}} c_{pv}\overline{θ}\frac{∂\overline{Π}}{∂z} + \frac{ρ′}{\overline{ρ}} c_{pv}\overline{θ}\frac{∂\overline{Π}}{∂z} − c_{pv}θ′ \frac{∂\overline{Π}}{∂z} + D_{w}
= −c_{pv}\overline{θ}\frac{∂Π′}{∂z} + \frac{ρ′_{s}}{\overline{ρ}} c_{pv}\overline{θ}\frac{∂\overline{Π}}{∂z} − c_{pv}θ′ \frac{∂\overline{Π}}{∂z} + D_{w}
= −c_{pv}\overline{θ}\frac{∂Π′}{∂z} − \frac{ρ′_{s}}{\overline{ρ}} g + \frac{θ′}{\overline{θ}} g + D_{w}
= −c_{pv}\overline{θ}\frac{∂Π′}{∂z} + g (\frac{θ′}{\overline{θ}} − \frac{R\overline{θ}}{p_{0}\overline{Π}^{c_{vv}/R}} ρ′_{s}) + D_{w},

\overline{2012/03/27(山下達也)}

\underline{18 大気主成分凝結対流運動方程式の導出}

\frac{∂w}{∂t} = −u\frac{∂w}{∂x} − w\frac{∂w}{∂z}
−c_{pv}\overline{θ}\frac{∂Π′}{∂z} + g (\frac{θ′}{\overline{θ}} − \frac{R\overline{θ}}{p_{0}\overline{Π}^{c_{vv}/R}} ρ′_{s}) + D_{w}  (A.8)

となり, (1.11) が得られる.

\overline{2012/03/27(山下達也)}

\underline{大気主成分凝結対流圧力方程式の導出 19}

付録B 圧力方程式の導出

本付録では圧力方程式 (1.13), (1.17) の導出を行なう. 状態方程式を θ, Π で表現すると,

ρ_{v} = \frac{p}{RT} = \frac{p_{0}Π^{c_{pv}/R}}{RΠθ} = \frac{p_{0}}{R} \frac{Π^{c_{vv}/R}}{θ}  (B.1)

となる. (B.1) の Lagrange 微分をとると,

\frac{Dρ_{v}}{Dt} = \frac{D}{Dt} (\frac{p_{0}}{R} \frac{Π^{c_{vv}/R}}{θ} )
= \frac{p_{0}c_{vv}}{R^{2}} \frac{Π^{c_{vv}/R−1}}{θ} \frac{DΠ}{Dt} − \frac{p_{0}}{R} \frac{Π^{c_{vv}/R}}{θ^{2}} \frac{Dθ}{Dt}
= \frac{c_{vv}}{RΠ} (\frac{p_{0}}{R} \frac{Π^{c_{vv}/R}}{θ} ) \frac{DΠ}{Dt} − \frac{1}{θ} (\frac{p_{0}}{R} \frac{Π^{c_{vv}/R}}{θ} ) \frac{Dθ}{Dt}
= \frac{c_{vv}ρ_{v}}{RΠ} \frac{DΠ}{Dt} − \frac{ρ_{v}}{θ} \frac{Dθ}{Dt}  (B.2)

ここで連続の式

\frac{Dρ_{v}}{Dt} = −ρ_{v}∇ · u − M_{cond}  (B.3)

及び熱力学の式

\frac{Dθ}{Dt} = \frac{1}{Π} (Q_{cond} + Q_{dis} + Q_{rad})  (B.4)

を (B.2) に代入すると,

\frac{DΠ}{Dt} = \frac{RΠ}{c_{vv}} (\frac{1}{ρ_{v}} \frac{Dρ_{v}}{Dt} + \frac{1}{θ} \frac{Dθ}{Dt})
= \frac{RΠ}{c_{vv}} [−∇ · u − \frac{1}{ρ_{v}} M_{cond} + \frac{1}{θΠ} (Q_{cond} + Q_{dis} + Q_{rad})]
= \frac{c^{2}}{c_{pv}θ} [−∇ · u − \frac{1}{ρ_{v}} M_{cond} + \frac{1}{θΠ} (Q_{cond} + Q_{dis} + Q_{rad})]
= −\frac{c^{2}}{c_{pv}θ}∇ · u − \frac{c^{2}}{c_{pv}ρ_{v}θ}M_{cond} + \frac{c^{2}}{c_{pv}θ^{2}Π} (Q_{cond} + Q_{dis} + Q_{rad}) (B.5)
\overline{2012/03/27(山下達也)}

\underline{20 大気主成分凝結対流圧力方程式の導出} となる. 但し熱力学の式において, 乱流拡散項は他の項に比べて十分小さいとして無視した. また c は音速であり,

c^{2} = \frac{c_{pv}}{c_{vv}} RθΠ  (B.6)

である. M_{cond}, Q_{∗} が擾乱成分であることに注意して, (B.5) の線形化を行なうと,

\frac{∂Π′}{∂t} = −w′ \frac{∂\overline{Π}}{∂z} − \frac{\overline{c^{2}}}{c_{pv}\overline{θ}}∇ · u′ − \frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}}M_{cond}
+\frac{\overline{c^{2}}}{c_{pv}\overline{θ}^{2}\overline{Π}} (Q_{cond+}Q_{dis} + Q_{rad})  (B.7)

となる. ここで

\overline{ρ}\overline{θ} = \frac{p_{0}}{R} \overline{Π}^{c_{vv}/R}  (B.8)
\overline{c^{2}} = \frac{c_{pv}}{c_{vv}} R\overline{θ}\overline{Π}  (B.9)

となることに着目すると,

\frac{∂}{∂z} (\overline{ρ}\overline{θ}) = \frac{∂}{∂z} (\frac{p_{0}}{R} \overline{Π}^{c_{vv}/R})
= \frac{c_{vv}}{R} \frac{p_{0}}{R} \overline{Π}^{c_{vv}/R−1} \frac{∂\overline{Π}}{∂z}
= \frac{c_{vv}}{R} \frac{\overline{ρ}\overline{θ}}{\overline{Π}} \frac{∂\overline{Π}}{∂z}
= \frac{c_{pv}\overline{ρ}\overline{θ}^{2}}{\overline{c^{2}}} \frac{∂\overline{Π}}{∂z}  (B.10)

となる. (B.10) より

−w′ \frac{∂\overline{Π}}{∂z} = −\frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}^{2}} w′ \frac{∂}{∂z} (\overline{ρ}\overline{θ})
= −\frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}^{2}} u′ · ∇ (\overline{ρ}\overline{θ})  (B.11)

となる. (B.11) を (B.7) に代入すると,

\frac{∂Π′}{∂t} = −\frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}^{2}} u′ · ∇ (\overline{ρ}\overline{θ}) − \frac{\overline{c^{2}}}{c_{pv}\overline{θ}}∇ · u′

+\frac{\overline{c^{2}}}{c_{pv}^{2}\overline{θ}^{2}\overline{Π}} (Q_{cond} + Q_{dis} + Q_{rad}) − \frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}}M_{cond}

\underline{大気主成分凝結対流圧力方程式の導出 21}

= −\frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}^{2}} [u′ · ∇ (\overline{ρ}\overline{θ}) + \overline{ρ}\overline{θ}∇ · u′]

+\frac{\overline{c^{2}}}{c_{pv}^{2}\overline{θ}^{2}\overline{Π}} (Q_{cond} + Q_{dis} + Q_{rad}) − \frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}}M_{cond}

= −\frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}^{2}} ∇ · (\overline{ρ}\overline{θ}u′) + \frac{\overline{c^{2}}}{c_{pv}\overline{θ}^{2}\overline{Π}} (Q_{cond} + Q_{dis} + Q_{rad}) − \frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}}M_{cond}

(B.12)

となり, (1.13) が得られる. 更に (B.12) に

Q_{cond} = \frac{LM_{cond}}{c_{pv}\overline{ρ}}  (B.13)

を代入して Q_{cond} を消去すると,

\frac{∂Π′}{∂t} = −\frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}^{2}} ∇ · (\overline{ρ}\overline{θ}u′) + \frac{\overline{c^{2}}}{c_{pv}\overline{θ}^{2}\overline{Π}} (\frac{LM_{cond}}{c_{pv}\overline{ρ}} + Q_{dis} + Q_{rad})

−\frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}}M_{cond}

= −\frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}^{2}} ∇ · (\overline{ρ}\overline{θ}u′) + \frac{\overline{c^{2}}}{c_{pv}\overline{θ}^{2}\overline{Π}} · \frac{LM_{cond}}{c_{pv}\overline{ρ}} − \frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}}M_{cond}

+\frac{\overline{c^{2}}}{c_{pv}\overline{θ}^{2}\overline{Π}} (Q_{dis} + Q_{rad})

= −\frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}^{2}} ∇ · (\overline{ρ}\overline{θ}u′) + \frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}} (\frac{L}{c_{pv}\overline{θ}\overline{Π}} − 1) M_{cond}

+\frac{\overline{c^{2}}}{c_{pv}\overline{θ}^{2}\overline{Π}} (Q_{dis} + Q_{rad})

= −\frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}^{2}} ∇ · (\overline{ρ}\overline{θ}u′) + \frac{\overline{c^{2}}}{c_{pv}\overline{ρ}\overline{θ}} (\frac{L}{c_{pv}\overline{T}} − 1) M_{cond}
+\frac{\overline{c^{2}}}{c_{pv}\overline{θ}^{2}\overline{Π}} (Q_{dis} + Q_{rad})  (B.14)

となり, (1.17) が得られる.

\underline{22 大気主成分凝結対流乱流パラメタリゼーション}

{付録C 主成分凝結系における乱流パ \atop ラメタリゼーション}

乱流パラメタリゼーション

Klemp and Wilhelmson (1978) および CReSS で用いられている 1.5 次のクロージャーを用いる. このとき乱流運動エネルギーの時間発展方程式は,

\frac{DE_{turb}}{Dt} = B + S + D_{E} − (\frac{C_{ε}}{l} ) Eturb^{\frac{3}{2}}  (C.1)

と与えられる. l は混合距離で, l = (ΔxΔz)^{1/2} とする. B と S はそれぞれ浮力と流れの変形速度による乱流エネルギー生成項, D_{E} は乱流エネルギー拡散項, 第 4 項は乱流エネルギーの消散項であり,

B = \frac{g_{j}}{\overline{θ}} \overline{u′_{j}θ′_{v}},  (C.2)
S = −\overline{(u′_{i}u′_{j})}\frac{∂u_{i}}{∂x_{j}} ,  (C.3)
D_{E} = \frac{∂}{∂x_{j}} (K_{m} \frac{∂E_{turb}}{∂x_{j}} )  (C.4)

である. 主成分凝結系においては

θ′_{v} = θ′ − \frac{\overline{θ}}{\overline{ρ}}ρ′_{s}
= θ′ − \frac{R\overline{θ}^{2}}{p_{0}\overline{Π}^{c_{v}/R}} ρ′_{s}  (C.5)

であるので,

B = \frac{g_{j}}{\overline{θ}} \overline{u′_{j}θ′_{v}}
= \frac{g_{j}}{\overline{θ}} \overline{u′_{j} (θ′ − \frac{R\overline{θ}^{2}}{p_{0}\overline{Π}^{c_{v}/R}} ρ′_{s})}

\underline{大気主成分凝結対流乱流パラメタリゼーション 23}

= \frac{g_{j}}{\overline{θ}} \overline{u′_{j}θ′} − g_{j} \frac{R\overline{θ}}{p_{0}\overline{Π}^{c_{v}/R}} \overline{u′_{j}ρ_{s}}

(C.6)

と表される. 1.5 次のクロージャーでは, レイノルズ応力を以下のように定義する.

\overline{(u′_{i}u′_{j})} = −K_{m} (\frac{∂u_{i}}{∂x_{j}} + \frac{∂u_{j}}{∂x_{i}} ) + \frac{2}{3}δ_{ij}E_{turb},  (C.7)
\overline{u′_{j}ϕ′} = −K_{h} \frac{∂ϕ}{∂x_{j}} .  (C.8)

ここで ϕ′ は任意のスカラー場の擾乱成分, K_{m} は運動量に対する渦粘性係数であり, E_{turb} はサブグリッドスケールの乱流運動エネルギー, K_{h} は渦拡散係数である. K_{m}, K_{h} は E_{turb} を用いて以下のように与えられる.

K_{m} = C_{m}Eturb^{\frac{1}{2}}l,  (C.9)
K_{h} = 3K_{m}.  (C.10)

パラメータ C_{ε}, C_{m} はともに 0.2 である.

2 次元の場合の表現

2 次元の場合の (C.1) 式の各項を書き下す. 浮力による乱流エネルギー生成項は,

B = \frac{g_{j}}{\overline{θ}} \overline{u′_{j}θ′} − g_{j} \frac{R\overline{θ}}{p_{0}\overline{Π}^{c_{v}/R}} \overline{u′_{j}ρ_{s}}
= \frac{g}{\overline{θ}} \overline{w′θ′} − g \frac{R\overline{θ}}{p_{0}\overline{Π}^{c_{v}/R}} \overline{w′ρ_{s}}
= −\frac{g}{\overline{θ}} (K_{h} \frac{∂θ}{∂z} ) + g \frac{R\overline{θ}}{p_{0}\overline{Π}^{c_{v}/R}} (K_{h} \frac{∂ρ_{s}}{∂z} )  (C.11)

である. 次に流れの変形速度による乱流エネルギー生成項 S は,

S = −\overline{(u′_{i}u′_{j})}\frac{∂u_{i}}{∂x_{j}}
= − \{−K_{m} (\frac{∂u_{i}}{∂x_{j}} + \frac{∂u_{j}}{∂x_{i}} ) + \frac{2}{3}δ_{ij}E_{turb}\} \frac{∂u_{i}}{∂x_{j}}
= \{K_{m} (\frac{∂u_{i}}{∂x_{j}} + \frac{∂u_{j}}{∂x_{i}} ) − \frac{2}{3}δ_{ij}E_{turb}\} \frac{∂u_{i}}{∂x_{j}}
= \{K_{m} (\frac{∂u}{∂x_{j}} + \frac{∂u_{j}}{∂x} ) − \frac{2}{3}δ_{1j}E_{turb}\} \frac{∂u}{∂x_{j}}

+ \{K_{m} (\frac{∂w}{∂x_{j}} + \frac{∂u_{j}}{∂z} ) − \frac{2}{3}δ_{3j}E_{turb}\} \frac{∂w}{∂x_{j}}

= \{2K_{m} (\frac{∂u}{∂x}) − \frac{2}{3}E_{turb}\} \frac{∂u}{∂x} + K_{m} (\frac{∂w}{∂x} + \frac{∂u}{∂z}) \frac{∂u}{∂z}

+K_{m} (\frac{∂w}{∂x} + \frac{∂u}{∂z}) \frac{∂w}{∂x} + \{2K_{m} (\frac{∂w}{∂z}) − \frac{2}{3}E_{turb}\} \frac{∂w}{∂z}

= 2K_{m} \{(\frac{∂u}{∂x})^{2} + (\frac{∂w}{∂z})^{2}\} + K_{m} (\frac{∂u}{∂z} + \frac{∂w}{∂x})^{2}
−\frac{2}{3}E_{turb} (\frac{∂u}{∂x} + \frac{∂w}{∂z})  (C.12)

である. 乱流エネルギー拡散項 D_{E} は,

D_{E} = \frac{∂}{∂x_{j}} (K_{m} \frac{∂E_{turb}}{∂x_{j}} ) ,
= \frac{∂}{∂x} (K_{m} \frac{∂E_{turb}}{∂x} ) + \frac{∂}{∂z} (K_{m} \frac{∂E_{turb}}{∂z} )  (C.13)

である. 以上の (C.11), (C.12), (C.13) 式を (C.1) 式に代入することで以下の式を得る.

\frac{DE_{turb}}{Dt} = −\frac{g}{\overline{θ}} (K_{h} \frac{∂θ}{∂z} ) + g \frac{R\overline{θ}}{p_{0}\overline{Π}^{c_{v}/R}} (K_{h} \frac{∂ρ_{s}}{∂z} )

+2K_{m} \{(\frac{∂u}{∂x})^{2} + (\frac{∂w}{∂z})^{2}\} + K_{m} (\frac{∂u}{∂z} + \frac{∂w}{∂x})^{2} − \frac{2}{3}E_{turb} (\frac{∂u}{∂x} + \frac{∂w}{∂z})

+\frac{∂}{∂x} (K_{m} \frac{∂E_{turb}}{∂x} ) + \frac{∂}{∂z} (K_{m} \frac{∂E_{turb}}{∂z} ) − (\frac{C_{ε}}{l} ) Eturb^{\frac{3}{2}}.  (C.14)

乱流拡散係数を用いた表現

(C.14) 式を (C.9) 式を用いて K_{m} に関する式に変形する. 右辺第 6, 7 項の乱流エネルギー拡散項を書き下すと,

\frac{∂}{∂x} (K_{m} \frac{∂E_{turb}}{∂x} ) + \frac{∂}{∂z} (K_{m} \frac{∂E_{turb}}{∂z} )
= \frac{1}{C_{m}^{2}l^{2}} \{\frac{∂}{∂x} (K_{m} \frac{∂K_{m}^{2}}{∂x} ) + \frac{∂}{∂z} (K_{m} \frac{∂K_{m}^{2}}{∂z} )\}
= \frac{1}{C_{m}^{2}l^{2}} \{K_{m} \frac{∂^{2}K_{m}^{2}}{∂x^{2}} + \frac{∂K_{m}}{∂x} \frac{∂K_{m}^{2}}{∂x} + K_{m} \frac{∂^{2}K_{m}^{2}}{∂z^{2}} + \frac{∂K_{m}}{∂z} \frac{∂K_{m}^{2}}{∂z} \}

= \frac{K_{m}}{C_{m}^{2}l^{2}} (\frac{∂^{2}K_{m}^{2}}{∂x^{2}} + \frac{∂^{2}K_{m}^{2}}{∂z^{2}} ) + \frac{2K_{m}}{C_{m}^{2}l^{2}} \{(\frac{∂K_{m}}{∂x} )^{2} + (\frac{∂K_{m}}{∂z} )^{2}\}

(C.15)

となるので, (C.14) 式を変形すると,

\frac{2K_{m}}{C_{m}^{2}l^{2}} \frac{DK_{m}}{Dt} = −\frac{g}{\overline{θ}} (K_{h} \frac{∂θ}{∂z} ) + g \frac{R\overline{θ}}{p_{0}\overline{Π}^{c_{v}/R}} (K_{h} \frac{∂ρ_{s}}{∂z} ) + 2K_{m} \{(\frac{∂u}{∂x})^{2} + (\frac{∂w}{∂z})^{2}\}

+K_{m} (\frac{∂u}{∂z} + \frac{∂w}{∂z})^{2} − \frac{2}{3} \frac{K_{m}^{2}}{C_{m}^{2}l^{2}} (\frac{∂u}{∂x} + \frac{∂w}{∂z})

+\frac{K_{m}}{C_{m}^{2}l^{2}} (\frac{∂^{2}K_{m}^{2}}{∂x^{2}} + \frac{∂^{2}K_{m}^{2}}{∂z^{2}} ) + \frac{2K_{m}}{C_{m}^{2}l^{2}} \{(\frac{∂K_{m}}{∂x} )^{2} + (\frac{∂K_{m}}{∂z} )^{2}\}

−\frac{C_{ε}}{C_{m}^{3}l^{4}} K_{m}^{3}.  (C.16)

係数を整理すると,

\frac{∂K_{m}}{∂t} = − (u\frac{∂K_{m}}{∂x} + w\frac{∂K_{m}}{∂z} ) − \frac{gC_{m}^{2}l^{2}}{2} \frac{K_{h}}{K_{m}} (\frac{1}{\overline{θ}} \frac{∂θ}{∂z} − \frac{R\overline{θ}}{p_{0}\overline{Π}^{c_{v}/R}} \frac{∂ρ_{s}}{∂z} )

+ (C_{m}^{2}l^{2}) \{(\frac{∂u}{∂x})^{2} + (\frac{∂w}{∂z})^{2}\}

+\frac{C_{m}^{2}l^{2}}{2} (\frac{∂u}{∂z} + \frac{∂w}{∂z})^{2} − \frac{K_{m}}{3} (\frac{∂u}{∂x} + \frac{∂w}{∂z})

+\frac{1}{2} (\frac{∂^{2}K_{m}^{2}}{∂x^{2}} + \frac{∂^{2}K_{m}^{2}}{∂z^{2}} ) + (\frac{∂K_{m}}{∂x} )^{2} + (\frac{∂K_{m}}{∂z} )^{2}

−\frac{C_{ε}}{2C_{m}l^{2}} K_{m}^{2}

となる. ここで C_{m} = C_{ε} = 0.2 と K_{h} = 3K_{m} という関係を用いると,

\frac{∂K_{m}}{∂t} = − (u\frac{∂K_{m}}{∂x} + w\frac{∂K_{m}}{∂z} ) − \frac{3gC_{m}^{2}l^{2}}{2} (\frac{1}{\overline{θ}} \frac{∂θ}{∂z} − \frac{R\overline{θ}}{p_{0}\overline{Π}^{c_{v}/R}} \frac{∂ρ_{s}}{∂z} )

+ (C_{m}^{2}l^{2}) \{(\frac{∂u}{∂x})^{2} + (\frac{∂w}{∂z})^{2}\}

+\frac{C_{m}^{2}l^{2}}{2} (\frac{∂u}{∂z} + \frac{∂w}{∂x})^{2} − \frac{K_{m}}{3} (\frac{∂u}{∂x} + \frac{∂w}{∂z})

+\frac{1}{2} (\frac{∂^{2}K_{m}^{2}}{∂x^{2}} + \frac{∂^{2}K_{m}^{2}}{∂z^{2}} ) + (\frac{∂K_{m}}{∂x} )^{2} + (\frac{∂K_{m}}{∂z} )^{2}

−\frac{1}{2l^{2}} K_{m}^{2}  (C.17)

が得られる.

\overline{2012/03/27(山下達也)}

\underline{26  大気主成分凝結対流  (1.39) の導出}

付録D (1.39) の導出

質量 m をもつ雲粒がまわりの大気の拡散によって成長する過程を考える. 拡散による雲粒の成長は

\frac{dm}{dt} = 4πr_{d}^{2}D\frac{∂ρ}{∂r}  (D.1)

で表される. r は雲粒の中心からの距離, D は大気の分子拡散係数である. 大気の状態は定常かつ等方であると仮定してこの式を r について解くと

\frac{dm}{dt} = −4πr_{d}D(ρ_{r_{d}} − ρ_{∞})  (D.2)

となる. ただし, 境界条件として r = r_{d} で ρ = ρ_{r_{d}}, r → ∞ で ρ → ρ_{∞} を用いた. ここで理想気体の状態方程式 p = ρRT を代入すると

\frac{dm}{dt} = −4πr_{d}D (\frac{p_{∗}(T_{r_{d}})}{RT_{r_{d}}} − \frac{p_{∞}}{RT_{∞}} )  (D.3)

となる. この式の両辺を p ∗ (T_{∞}) で割ると

\frac{1}{p_{∗}(T_{∞})} \frac{dm}{dt} = −\frac{4πr_{d}D}{R} (\frac{p_{∗}(T_{r_{d}})}{p_{∗}(T_{r_{d}})T_{r_{d}}} − \frac{p_{∞}}{p_{∗}(T_{r_{d}})T_{∞}} )
= −\frac{4πr_{d}D}{RT_{∞}} (\frac{p_{∗}(T_{r_{d}})}{p_{∗}(T_{∞})} − \frac{p_{∞}}{p_{∗}(T_{∞})})  (D.4)

となる. ただし, 最後の変形には (T_{r_{d}} − T_{∞})/T_{∞} ≪ 1 を用いた.

凝結が起きた時潜熱が解放される. この潜熱が熱伝導によって輸送されると仮定すると

L\frac{dm}{dt} = −4πr^{2}k \frac{dT}{dr}  (D.5)

が成り立つ. ここで k は大気の熱拡散係数である. 大気密度の拡散方程式と同様にこの式を解くと

L\frac{dm}{dt} = 4πrk(T_{r_{d}} − T_{∞})  (D.6)

\underline{大気主成分凝結対流  (1.39) の導出  27}

となる. ただし, 境界条件として r = r_{d} で T = T_{r_{d}}, r → ∞ で T → T_{∞} を用いた. ここでクラウジウス-クラペイロンの式

\frac{dp_{∗}}{dT} = \frac{Lp_{∗}}{RT^{2}}

を積分すると

ln \frac{p_{∗}(T_{r_{d}})}{p_{∗}(T_{∞})} = \frac{L}{R} \frac{T_{r_{d}} − T_{∞}}{T_{r_{d}}T_{∞}}
∼ \frac{L}{R} \frac{T_{r_{d}} − T_{∞}}{T_{∞}^{2}}  (D.7)

となる. よって

\frac{p_{∗}(T_{r_{d}})}{p_{∗}(T_{∞})} = exp (\frac{L}{R} \frac{T_{r_{d}} − T_{∞}}{T_{r_{d}}T_{∞}} )
∼ 1 + \frac{L}{R} \frac{T_{r_{d}} − T_{∞}}{T_{∞}^{2}}
∼ 1 + \frac{L^{2}}{4πr_{d}kRT_{∞}^{2}} \frac{dm}{dt}  (D.8)

となる. 最後の変形には式 (D.6) を用いた.

式 (D.4), (D.8) より

\frac{dm}{dt} = −\frac{4πrDp_{∗}(T_{∞})}{RT_{∞}} (1 + \frac{L^{2}}{4πr_{d}kRT_{∞}^{2}} \frac{dm}{dt} − \frac{p_{∞}}{p_{∗}(T_{∞})})  (D.9)

となる. この式を整理すると

(\frac{RT_{∞}}{Dp_{∗}(T_{∞})} + \frac{L^{2}}{kRT_{∞}^{2}} ) \frac{dm}{dt} = 4πr_{d} (\frac{p_{∞}}{p_{∗}(T_{∞})} − 1)  (D.10)

となる. ここで

R_{m} = \frac{RT_{∞}}{Dp_{∗}(T_{∞})}
R_{h} = \frac{L^{2}}{kRT_{∞}^{2}}
S = \frac{p_{∞}}{p_{∗}(T_{∞})}

とおくと

\frac{dm}{dt} = \frac{4πr_{d}}{R_{m} + R_{h}} (S − 1)  (D.11)

となる. R_{m}, R_{h} はそれぞれ質量輸送, 熱輸送に関係する係数である. S は飽和比である.

単位体積当たりの雲粒の個数を N とし, 雲粒の大きさが全て同じであると仮定すれば単位体積当たりの凝結量 M_{cond} は

M_{cond} = N \frac{dm}{dt} = \frac{4πr_{d}N}{R_{m} + R_{h}} (S − 1)  (D.12)

で与えられる. 火星極冠周縁での温度・圧力条件を想定すると, R_{m} ≪ R_{h} である. 従って

M_{cond} ∼ \frac{4πr_{d}N}{R_{h}} (S − 1)  (D.13)

となる.

\underline{大気主成分凝結対流下層で飽和比一定, 上層で等温な基本場温度分布 29}

{付録E 下層で飽和比一定, 上層で等 \atop 温である基本場温度分布について}

Colaprete et al.(2003) の数値実験では初期に飽和比が一定である温度分布を与えている. ここでは下層で飽和比一定, 上層で等温である基本場温度分布\overline{T} を与える場合を考える. 上層での温度 T_{iso}, 下層での飽和比の値 S_{0}, 地表面圧力 p_{s} は既知であるとする. 先ず飽和比一定の層での温度分布 T_{S_{0}}(z) を求める. 基本場においては静水圧平衡が成り立つので,

\frac{dp}{dz} = −ρg  (E.1)

となる. 飽和蒸気圧を p_{∗} とすると,

p = S_{0}p_{∗}  (E.2)

が成り立つ. 式 (E.2) を式 (E.1) に代入すると,

\frac{dp}{dz} = S_{0} \frac{dp_{∗}}{dz} = S_{0} \frac{dp_{∗}}{dT} \frac{dT}{dz} = −ρg  (E.3)

となる. ここでクラウジウス・クラペイロンの式

\frac{dp_{∗}}{dT} = \frac{Lp_{∗}}{RT^{2}}  (E.4)

及び理想気体の状態方程式 p = ρRT を用いると,

S_{0} \frac{Lp_{∗}}{RT^{2}} \frac{dT}{dz} = −\frac{S_{0}p_{∗}}{RT} g,
\frac{dT}{dz} = −\frac{g}{L}T  (E.5)

となる. z = 0 での温度を T_{0} とし, 式 (E.5) を z = 0 から z まで積分すると,

T_{S_{0}}(z) = T_{0} exp (−\frac{g}{L}z)  (E.6)

となる. Antoine の式

p_{∗} = exp (A_{ant} − \frac{B_{ant}}{T} )  (E.7)

\underline{30 大気主成分凝結対流下層で飽和比一定, 上層で等温な基本場温度分布}

において z = 0 の場合を考えると,

\frac{p_{s}}{S_{0}} = exp (A_{ant} − \frac{B_{ant}}{T_{0}} ) ,
T_{0} = \frac{B_{ant}}{A_{ant} − ln(p_{s}/S_{0})}  (E.8)

となる. 式 (E.8) を式 (E.7) に代入すると,

T_{S_{0}}(z) = \frac{B_{ant}}{A_{ant} − ln(p_{s}/S_{0})} exp (−\frac{g}{L}z)  (E.9)

となる. 更に温度一定の分布 T = T_{iso} と飽和比一定の分布 T = T_{S_{0}} が交わる高度を z_{T} とすると, 式 (E.9) より

T_{iso} = \frac{B_{ant}}{A_{ant} − ln(p_{s}/S_{0})} exp (−\frac{g}{L}z_{T} ) ,
z_{T} = −\frac{L}{g} ln [\frac{A − ln(p_{s}/S_{0})}{B}T_{iso}]  (E.10)

となる. 以上より,

\overline{T} (z) = \begin{cases}\frac{B_{ant}}{A_{ant}−ln(p_{s}/S_{0})} exp (−\frac{g}{L}z) , (z ≤ z_{T}) \\ T_{iso}  (z > z_{T}) \end{cases}  (E.11)

となる. 但し温度分布を式 (E.11) と与える場合, z = z_{T} に温度勾配の不連続が存在し, 局所的に大きな温度の拡散が生じてしまうことになる. そこで実際の数値計算では, 以下のように重み関数 W (z) を導入し, 温度一定の分布と飽和比一定の温度分布をなめらかにつないだものを与えている.

\overline{T} (z) = [1 − W (z)]\frac{B_{ant}}{A_{ant} − ln(p_{s}/S_{0})} exp (−\frac{g}{L}z) + W (z)T_{iso}, (E.12)
W (z) = \frac{1}{2} [tanh (\frac{z − z_{T}}{Δh} ) + 1] .  (E.13)

但し Δh は温度分布の遷移の特徴的な長さスケールである.