DCRTM 支配方程式とその離散化

{地球流体電脳倶楽部 \atop DCRTM 開発プロジェクト}

平成 28 年 11 月 11 日

目次

第 1 章はじめに 3 1.1 この文書について . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 DCRTM の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

第 2 章座標系 4 2.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 基本グリッド . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.3 物理量の配置 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

第 3 章吸収過程 5 3.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3.2 線吸収 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3.3 連続吸収 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.4 紫外線領域での吸収 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3.4.1 H2O の紫外吸収 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

第 4 章散乱過程 9 4.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 4.2 大気分子による散乱: Rayleigh 散乱 . . . . . . . . . . . . . . . . . . . . . . . . . . 9

4.2.1 空気以外の分子による散乱 . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 4.2.2 空気による散乱 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

4.3 ダスト, エアロゾルによる散乱: Mie 散乱 . . . . . . . . . . . . . . . . . . . . . . . 10

第 5 章光学的厚さの計算 11

第 6 章放射伝達計算 14 6.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 6.2 Two-Stream 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 6.3 Two-Stream Source Function 法 (放射源関数法) . . . . . . . . . . . . . . . . . . . 14

6.3.1 Toon et al. (1989) による Two-Stream Source Function 法 . . . . . . . . . 15 6.3.2 放射源の放射場を self consistent に求める . . . . . . . . . . . . . . . . . . 18

第 7 章相関 k 分布法による吸収断面積, 散乱断面積の計算 22 7.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 7.2 吸収物質が 1 種類の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

7.2.1 吸収断面積の計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 7.2.2 散乱断面積の計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 第 8 章その他 24

8.1 断熱温度勾配 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 8.1.1 乾燥断熱減率 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 8.1.2 湿潤断熱減率 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

8.2 潜熱 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 8.3 H_{2}O の状態方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 8.4 中心星のスペクトル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 8.5 地表面のアルベド . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

第1章はじめに

1.1 この文書について

この文章は, 地球流体電脳倶楽部で開発中の惑星大気放射モデル DCRTM の支配方程式系およびその離散化手法を解説したものである.

1.2 DCRTM の概要

ここでは, DCRTM の概要を示す.

• line-by-line 放射モデル

• GCM のための k 分布法による吸収・散乱断面積計算モデル

第2章座標系

2.1 はじめに

2.2 基本グリッド

DCRTM では, 鉛直レベルとして圧力座標, 波数方向に波数座標を持つ. 圧力座標は, 大気層の境界の値の圧力を与える. 波数座標は, 計算波数領域の最小波数, 最大波数の間を等間隔に切っている. 波数間隔は, 波数解

像度として与える.

2.3 物理量の配置

DCRTM では, 鉛直レベルに関して, 鉛直大気層の中心に配置される物理量と, 大気層の境界に配置される物理量がある.

圧力, 温度, 体積混合比, 大気上端からの光学的厚さ, フラックスは大気境界に配置される. 各大気層の光学的厚さ, 一次散乱アルベド, 加熱率は大気層の中心に配置される.

第3章吸収過程

3.1 はじめに

ここでは, 気体分子による吸収過程について扱う.

3.2 線吸収

気体分子 m の同位体 i による, η から η′ のエネルギー遷移に対応する波数 ν_{ηη′} での吸収断面積 k_{m,i,ηη′} [m^{2}/molecule] は, 圧力 p, 温度 T の関数として以下のようにかける.

k_{m,i,ηη′}(ν, T, p) = S_{m,i,ηη′}(T)f_{m,i,ηη′}(ν, ν_{ηη′}, T, p)  (3.1)

S_{m,i,ηη′}(T)[m^{−1}/(molecule m^{−2})] は吸収線強度, f_{m,i,ηη′}(ν, ν_{ηη′}, T, p)[m] は吸収線の広がり (形) を決める関数である.

吸収線強度

吸収線強度 S_{m,i,ηη′}(T) は温度の関数であり, その温度依存性は以下のようにかける.

S_{m,i,ηη′}(T) = S_{m,i,ηη′}(T_{ref})\frac{Q_{m,i}(T_{ref})}{Q_{m,i}(T)} \frac{exp(−hcE_{η}/k_{B}T)}{exp(−hcE_{η}/k_{B}T_{ref})} \frac{1 − exp(−hcν_{ηη′}/T)}{1 − exp(−hcν_{ηη′}/T_{ref})}  (3.2)

ここで Q_{m,i}(T) は, 温度 T における気体分子 m の同位体 i の分配関数. T_{ref} はある基準の温度である. また E_{η} は基底状態における分子のエネルギーであり, h, c, k_{B} はそれぞれ, プランク定数, 光速, ボルツマン定数である.

S_{m,i,ηη′}(T_{ref}), E_{η} は HITRAN, HITEMP データベースより得ることができる. また, 分配関数 Q_{m,i}(T) は HITRAN2008 に同梱されている parsum.dat ファイル内にテーブルとして格納されている. DCRTM では本テーブルを用いて, 必要な温度の分配関数を線型内挿することで得ている.

吸収線の広がり

吸収線の広がりは主に以下のふたつの過程によって引き起こされる (浅野, 2002).

• 分子同士の衝突

• 分子の視線方向速度によるドップラー効果

これらの過程により引き起こされる吸収線の広がりは, それぞれ Lorentz 線型, Doppler 線型で特徴づけられる. さらに両者の重ね合わせによる吸収線型を Voigt 線型と呼ぶ.

1. 分子同士の衝突: Lorentz 線型 Lorentz 線型で特徴づけられる吸収線の広がりは以下のように表される.

f_{L(m,i,ηη′)}(ν, ν_{ηη′}, T, p) = \frac{1}{π} \frac{γ_{L(m,i,ηη′)}(p, T)}{(γ_{L(m,i,ηη′)}(p, T))^{2} + (ν − ν_{ηη′})^{2}}  (3.3)

但し, γ_{L(m,i,ηη′)} は Lorentz 半値幅で, 以下のように表される.

γ_{L(m,i,ηη′)}(p, T) = γ_{L(m,i,ηη′)}(p_{ref}, T_{ref}) (\frac{p}{p_{ref}} ) (\frac{T_{ref}}{T} )^{n_{L(m,i,ηη′)}}  (3.4)

ここで, p_{ref} は基準の圧力であり, n_{L(m,i,ηη′)} は Lorentz 半値幅の温度依存指数である. Lorentz 半値幅は, 衝突する分子同士が, 同種の分子か異なる分子かによって以下のように表される.

γ_{L(m,i,ηη′)} = p_{s}γ_{L(m,i,ηη′),s} + p_{f}γ_{L(m,i,ηη′),f}  (3.5)

p_{s}, p_{f}, γ_{L(m,i,ηη′),s}, γ_{L(m,i,ηη′),f} はそれぞれ, 分子 m の分圧, それ以外の分子の分圧, 分子 m の self-broadening による半値幅, 分子 m の foreign broadening による半値幅である. さらに圧力により吸収線の中心がシフトする効果を考慮する. シフトした中心波数を ν_{ηη′}^{∗} とすれば,

ν_{ηη′}^{∗} = ν_{ηη′} + δ(p_{ref})p (\frac{T_{ref}}{T} )  (3.6)

となる. ここで, δ(p_{ref}) は圧力シフトを与える係数である. 吸収線の圧力シフトを考慮した線型は以下のようになる.

f_{L(m,i,ηη′)}(ν, ν_{ηη′}, T, p) = \frac{1}{π} \frac{γ_{L(m,i,ηη′)}(p, T)}{(γ_{L(m,i,ηη′)}(p, T))^{2} + (ν − ν_{ηη′}^{∗})^{2}} .  (3.7)

ν_{ηη′}, γ_{L(m,i,ηη′),s}(p_{ref}, T_{ref}), γ_{L(m,i,ηη′),f}(p_{ref}, T_{ref}), n_{L(m,i,ηη′)}, δ(p_{ref}) は HITRAN, HITEMP のデータベースより得ることができる. ※ Rothman et al., 1998 によれば, 圧力シフトによる中心波数は

ν_{ηη′}^{∗} = ν_{ηη′} + δ(p_{ref})p  (3.8)

となっている. 一方で, LBLRTM では, 圧力シフトは p (\frac{T_{ref}}{T} ) に比例するように与えられている. DCRTM では, LBLRTM と同様の方法で圧力シフトを与えている.

2. 分子の視線方向速度によるドップラー効果: Doppler 線型分子のドップラー効果による吸収線の広がりは以下のように表される.

f_{D(m,i,ηη′)}(ν, ν_{ηη′}, T) = \frac{1}{γ_{D}(T)π^{1/2}} exp[−(\frac{ν − ν_{ηη′}}{γ_{D}(T)} )^{2}]  (3.9)
γ_{D}(T) = \frac{ν_{ηη′}}{c} (\frac{2RT}{M_{m,i}})^{1/2}  (3.10)

3. Lorentz 線型と Doppler 線型の重ね合わせ: Voigt 線型実際の大気においては吸収線の広がりは衝突効果とドップラー効果の両方により決まり, Voigt 線型と呼ばれる. 吸収線に及ぼす衝突効果, ドップラー効果は互いに独立であると仮定することで, ふたつの線型のたたみこみにより Voigt 線型を以下のように得ることができる.

f_{V (m,i,ηη′)}(ν, ν_{ηη′}, p, T) = ∫_{0}^{∞} f_{L(m,i,ηη′)}(ν′, ν_{ηη′}, T, p)f_{D(m,i,ν′)}(ν, ν′, T)dν′ (3.11)

= ∫_{0}^{∞} \frac{1}{π} \frac{γ_{L(m,i,ηη′)}(p, T)}{(γ_{L(m,i,ηη′)}(p, T))^{2} + (ν′ − ν_{ηη′}^{∗})^{2}} \frac{1}{γ_{D}(T)π^{1/2}} exp[−(\frac{ν − ν′}{γ_{D}(T)})^{2}]dν′  (3.12)

= \frac{1}{γ_{D}(T)√π} × V oigt(x, y).  (3.13)

ただし,

V oigt(x, y) = \frac{y}{π} ∫_{−∞}^{∞} \frac{1}{y^{2} + (x − t)^{2}} exp[−t^{2}]dt,  (3.14)

x = \frac{ν − ν_{ηη′}^{∗}}{γ_{D}(T)} , y = \frac{γ_{L(m,i,ηη′)}(p, T)}{γ_{D}(T)}, t = \frac{ν − ν′}{γ_{D}(T)}, γ_{D}(T) = \frac{ν_{ηη′}^{∗}}{c} (\frac{2RT}{M_{m,i}})^{1/2}  (3.15)

Voigt 関数を数値的に解く手法はいくつか公開されている. DCRTM では, Humlicek, 1982 による手法を用いている. 尚, LBLRTM でも Humlicek, 1982 を用いている (Clough et al., 2005).

式 (3.1) で得られる吸収線強度を今考えている分子種のすべての遷移エネルギーに対して足し合わせることで, 吸収スペクトルを得ることができる. DCRTM では, 連続吸収を MTCKD モデルを用いて計算するため, 吸収線の広がりは中心波数から 25[cm^{−1}] でカットオフしている.

3.3 連続吸収

大気分子の光の吸収過程は, 分子の振動や回転準位などのエネルギー遷移によるものだけでなく, 連続的な吸収があることが知られている. 特に地球大気においては, 大気の窓領域において水蒸気の連続吸収があることが古くから知られ, その成因などについて今日まで議論がある (例えば, Shine et al., 2012).

近年地球大気の放射計算において最も広く用いられているのは MT CKD 連続吸収モデル (Mlawer et al., 2012) である. MTCKD 連続吸収モデルは, 分子の吸収スペクトルを, 吸収線の中心から両翼 25[cm^{−1}] までの「線吸収」と「それ以外の吸収」にわけ, 「それ以外の吸収」の吸収を連続吸収として扱っている. したがって, MTCKD 連続吸収モデルを使用する場合には, 線吸収の計算において 25[cm^{−1}] のカットオフを設けて線吸収を計算する必要がある. そのように計算された線吸収と MTCKD 連続吸収モデルの連続吸収を合わせることで, 線吸収と連続吸収を考慮した吸収断面積を計算することができる.

3.4 紫外線領域での吸収

HITRAN2008 (Rothman et al., 2009), HITEMP2010 (Rothman et al., 2010) では 30,000 [cm^{−1}] までの吸収線, MTCKD 連続吸収モデルでは, 水蒸気: 20000[cm^{−1}], 二酸化炭素 : 10000[cm^{−1}], 窒素 : 2570[cm^{−1}], 酸素 : 29870[cm^{−1}], オゾン : 54000[cm^{−1}] までの連続吸収を扱っている. そのため, さらに高波数の吸収を取り扱う必要がある場合には, これら以外のデータから吸収断面積を入手する必要がある.

紫外領域の分子の吸収については, 下記の書籍が参考になる.

Joseph Berkowitz, Atomic and Molecular Photoabsorption, Absolute total cross sections, Academic press, 2002.

3.4.1 H2O の紫外吸収

DCRTM では H2O の紫外領域での吸収断面積を下記の文献から得ている.

Chan et al., 1993, The electronic spectrum of water in the discrete and continuum regions. Absolute optical oscillator strengths for photoabsorption (6 - 200 eV), Chemical Phyisics 178, 387 - 400.

第4章散乱過程

4.1 はじめに

4.2 大気分子による散乱: Rayleigh 散乱

入射光の波長に対して小さく, 均質な球形粒子による散乱を考える. この散乱は Rayleigh 散乱と呼ばれ, 大気分子の散乱過程をよく表現する. 入射する電磁場によって粒子内に誘起される微小な電気双極子による 2 次放射を散乱光であるとして, 入射場に対する 2 次放射場の割合を求める. これが散乱断面積 σ_{s}^{R} である.

Rayleigh 散乱による散乱断面積 σ_{s}^{R} を, 分子の分極率 α, 波長 λ でかくと, 以下のようになる (浅野, 2002).

σ_{s}^{R} = \frac{128π^{5}}{3λ^{4}} α^{2}.  (4.1)

一方, 分子の分極率 α は Lorentz-Lrenz の式により, 複素屈折率 m̃ と関係付けられる.

α = \frac{3}{4πN_{s}} \frac{m̃^{2} − 1}{m̃^{2} + 2}  (4.2)

ここで, N_{s} は単位体積あたりの分子の数である. これらにより, 散乱断面積 σ_{s}^{R} は,

σ_{s}^{R} = \frac{24π^{3}}{λ^{4}N_{s}^{2}} \frac{(m̃^{2} − 1)^{2}}{(m̃^{2} + 2)^{2}} f(δ)  (4.3)

となる. ただし, 上式では実際の分子の分極が等方的でないことに対する補正項 f(δ) を付加している. この項は, カバネスの校正因子 (Cabanes correction factor) と呼ばれ,

f(δ) = \frac{6 + 3δ}{6 − 7δ}  (4.4)

で与えられる. ここに, δ は誘導分極の非等方性を表す因子であり, 偏向解消因子 (depolarization factor) と呼ばれる.

4.2.1 空気以外の分子による散乱

空気以外の分子の散乱は, 式 (4.1) に従って計算している. CRC HANDBOOK により, 各分子の分極率を与えることで, 散乱断面積を計算する. 分極率の波数依存性は無視している.

4.2.2 空気による散乱

空気は, 窒素や酸素などの混合気体であるが, 空気による散乱断面積は, Bucholtz, 1995 で与えている, 屈折率, 非等方因子を用いて, 空気による散乱断面積を与える.

Bucholtz, 1995 では, 圧力 101325 [Pa], 温度 15 [℃] における空気の屈折率を与えている. この状態のとき, 単位体積あたりの空気分子の数 N_{s}は, 2.54743 × 10^{25}[m^{−3}] である. また屈折率は実数で与えている. すなわち m̃〜m_{r}(m_{r}は屈折率の実部).

空気の屈折率 m_{r} は以下のように 2 つの波長領域で与えられている. λ > 0.23[μm]:

(m_{r} − 1) × 10^{8} = \frac{5791817}{238.085 − (1/λ)^{2}} + \frac{167909}{57.362 − (1/λ)^{2}}  (4.5)

λ ≤ 0.23[μm]:

(m_{r} − 1) × 10^{8} = 8060.51 + \frac{2480990}{132.274 − (1/λ)^{2}} + \frac{17455.7}{39.32957 − (1/λ)^{2}} .  (4.6)

ここで, 屈折率中の波長 λは [μm] 単位で与える. Buholtz, 1995 では, 非等方因子 δ 値として, Penndorf 1957, Bates 1984, Nicolet 1984 を比較

し, 最終的に Bates 1984 を用いた場合の散乱断面積のフィッティング式 (波長の関数) を得ている:

σ_{s}^{R} = Aλ^{−(B+Cλ+D/λ)}  (4.7)

ここで, 係数 A, B, C, D は以下のように与えられ, λは [μm], 散乱断面積 σ_{s}^{R}は [cm^{2}] 単位である. 0.2 ≤ λ ≤ 0.5[μm]: A = 3.01577 × 10^{−28}, B = 3.55212, C = 1.35579, D = 0.11563. λ > 0.5[μm]: A = 4.01061 × 10^{−28}, B = 3.99668, C = 1.10298 × 10^{−3}, D = 2.71393 × 10^{−2}.

4.3 ダスト, エアロゾルによる散乱: Mie 散乱

第5章光学的厚さの計算

ある大気層の光学的厚さ τ の定式化を記述する。計算しようとする大気層に対して、与えられている物理量は大気層の上端、下端の圧力 p′, 温度

T′, 体積混合比 x′ である。 • p′_{t}: 上端の圧力

• p′_{b}: 下端の圧力

• T′_{t}: 上端の温度

• T′_{b}: 下端の温度

• x′_{m,t} : 大気層上端の分子 m の体積混合比

• x′_{m,b} : 大気層下端の分子 m の体積混合比また、計算しようとしている光学的厚さ τ は大気層に対して与えられる値である。ここでは、大

気層の上端下端で与えられる物理量にはプライム (′) をつけ、大気層を代表する物理量にはプライムをつけない。

• τ_{m} : 大気分子 m の光学的厚さ

• τ : 大気の光学的厚ささらに必要な物理量を以下のように定義する。

• p: 大気層を代表する圧力

• T: 大気層を代表する温度

• x_{m} : 分子 m の、大気層を代表する体積混合比

• k′_{m,t}: 分子 m の、上端の吸収係数

• k′_{m,b}: 分子 m の、下端の吸収係数

• k_{m} : 分子 m の、大気層を代表する吸収係数

• u_{m} : 大気層中の分子 m の分子の数

• g : 重力加速度

• N_{a} : アヴォガドロ数

• M_{ave} : 平均分子量

大気層上端 --------------- p'_t, T'_t, x'_m,t, k'_m,t

           -  -  -  -  -  p, T, x_m, k_m, u_m, tau_m, tau

大気層下端 --------------- p'_b, T'_b, x'_m,b, k'_m,b

1. 大気層の上端、下端の吸収係数の計算大気層中に含まれる分子 m による、波数 ν の吸収係数 k′_{m,j}(ν) は、圧力 p′_{j}、温度 T′_{j}、大気層に

含まれる分子の体積混合比 x′_{m,j} を使って計算することができる (j : t or b, m : 1, 2, ・・・ n(n は今考えている大気層中の大気分子の種類の数))。

k′_{m,j}(ν) = k′_{m,j}(ν, p′_{j}, T′_{j}, x′_{m,j})  (5.1)

※体積混合比は self-broadening, foreign-broadening の大きさを決める。今考えている分子 m と同じ分子がどれだけあるか、違う分子がどれだけあるかが分かればよいので、吸収係数 k′_{m,j}の引数は、分子 m の分圧 x′_{m,j}のみでよい。

2. 大気層を代表する吸収係数の計算分子 m の大気層を代表する吸収係数 k_{m} は、大気上端下端の吸収係数 k′_{m,j} に分子の数の重みを

つけて平均することで求める。大気上端下端における分子 m の数密度 n′_{m,j} は

n′_{m,j} = \frac{p′_{j}}{RT′_{j}} x′_{m,j}  (5.2)

となり (R は気体定数)、分子 m の大気層を代表する吸収係数 k_{m}(ν) は以下のようにかける。

k_{m}(ν) = \frac{n′_{m,b} × k′_{m,b}(ν) + n′_{m,t} × k′_{m,t}(ν)}{n′_{m,b} + n′_{m,t}}  (5.3)

3. 大気層中の分子の数の計算大気層中の分子 m の分子の数は、単位面積当たりのカラム量として下記のようにかける。

u_{m} = \frac{p′_{b} − p′_{t}}{g} \frac{N_{a}}{M_{ave}} x_{m}  (5.4)

ここで、

x_{m} = \frac{n′_{m,b} + n′_{m,t}}{Σ_{m}(n′_{m,b} + n′_{m,t})}.  (5.5)

また大気層の平均分子量は M_{ave} は、以下のように与えている。分子 m の同位体 i の分子量と存在比をそれぞれ M_{m,i}、 x_{m,i} とすれば、分子 m の平均分子量 M_{m}

は

M_{m} = Σ_{i}(x_{m,i} × M_{m,i})  (5.6)

となり、大気層の平均分子量 M_{ave} は、

M_{ave} = \frac{Σ_{m}(n′_{m,b} × M_{m} + n′_{m,t} × M_{m})}{Σ_{m}(n′_{m,b} + n′_{m,t})}  (5.7)

となる。 4. 光学的厚さの計算分子 m の光学的厚さ τ_{m}(ν) は吸収係数と分子の数を掛け合わせて計算することができる。よっ

て、分子 m の光学的厚さ τ_{m}(ν) は

τ_{m}(ν) = u_{m} × k_{m}(ν)  (5.8)

また、全ての分子の寄与を合わせた光学的厚さ τ(ν) は

τ(ν) = Σ_{m=1,n}(τ_{m}(ν))  (5.9)

これを大気の上端、下端の量を使って書き下せば、

τ(ν) = \frac{p′_{b} − p′_{t}}{g} N_{a} \frac{Σ_{j}(p′_{j}/T′_{j}Σ_{m}(x′_{m,j}k′_{m,j}(ν)))}{Σ_{j}(p′_{j}/T′_{j}Σ_{m}(x′_{m,j}M_{m}))}  (5.10)

となる。 ※補足: 大気の上端下端で混合分子の吸収係数を考えて、光学的厚さを計算する場合の定式化層の上端下端における混合分子の吸収係数 k′_{j}(ν) は

k′_{j}(ν) = Σ_{m}(x′_{m,j} × k′_{m,j}(ν))  (5.11)

となる。大気を代表する混合分子の吸収係数 k(ν) は

k(ν) = \frac{n′_{b}k′_{b}(ν) + n′_{t}k′_{t}(ν)}{n′_{b} + n′_{t}}  (5.12)

となる。よって光学的厚さ τ(ν) は

τ(ν) = u × k(ν)
= \frac{p′_{b} − p′_{t}}{g} \frac{N_{a}}{M_{ave}} \frac{n′_{b}k′_{b}(ν) + n′_{t}k′_{t}(ν)}{n′_{b} + n′_{t}}
= \frac{p′_{b} − p′_{t}}{g} N_{a} \frac{Σ_{j}(p′_{j}/T′_{j}Σ_{m}(x′_{m,j}k′_{m,j}(ν)))}{Σ_{j}(p′_{j}/T′_{j}Σ_{m}(x′_{m,j}M_{m}))}.  (5.13)

第6章放射伝達計算

6.1 はじめに

ここでは, 放射伝達計算について述べる. DCRTM の放射伝達計算は, Toon et al. (1989) を基本としている. 中心星放射が惑星大気により散乱されるされることによる放射は, Toon et al., (1989) の Two-Steam 法を用いて計算する. 散乱は δ-Eddington 近似を用いている. 惑星大気が射出することによる放射伝達は, Toon et al., (1989) の Two-Stream Source Function 法を用いて計算する. この計算には惑星放射の散乱も考慮されている. 惑星放射の散乱は, δ-半球平均近似による. この Two-Stream Source Function 法は, Two-Stream 法により計算された放射場を放射源関数として用いて放射伝達を解く方法であり, Toon et al., (1989) によれば, 惑星放射の放射伝達計算については, Two-Stream 法よりも一般に精度がよいとされている. しかし, ある条件のもとでは, Two-Stream Source Function 法による放射場で計算した加熱率が物理的に不自然な値となることがある. これは, Two-Stream 法による放射場と, それを放射源とした Two-Stream Source Function 法による放射場が整合的でない場合に起こる. 本来 2 つの放射場は一致すべきものであり, DCRTM では, Two-Stream Source Function 法による放射場を再度放射源として放射伝達式を解くことを繰り返し, self consistent な解を計算している.

6.2 Two-Stream 法

DCRTM の放射伝達計算モジュールは, DCPAM の放射伝達モジュールをもとにしている. Two-Stream 法の方程式については, DCPAM のドキュメント (DCPAM5 支配方程式とその離散化) を参照のこと.

6.3 Two-Stream Source Function 法 (放射源関数法)

放射伝達方程式 (波数の添え字は省略, 方位角方向には積分済)

μ\frac{∂I(τ, μ)}{∂τ} = I(τ, μ) − S(τ, μ)  (6.1)

より, 第 n 層の上端からの上向き放射 I_{n}^{+}(0, μ) は,

I_{n}^{+}(0, μ) = I_{n}^{+}(τ_{n}, μ) exp[−τ_{n}/μ] + ∫_{0}^{τ_{n}} S_{n}^{+}(t, μ) exp[−t/μ]\frac{dt}{μ}.  (6.2)

同様に, 第 n 層の下端からの下向き放射 I_{n}^{−}(τ_{n}, −μ) は,

I_{n}^{−}(τ_{n}, −μ) = I_{n}^{−}(0, −μ) exp[−τ_{n}/μ] + ∫_{0}^{τ_{n}} S_{n}^{−}(t, −μ) exp[−(τ_{n} − t)/μ]\frac{dt}{μ}.  (6.3)

S(τ, μ) は放射源を表す関数であり,

S(τ, μ) = S_{e}(τ, μ) + \frac{ϖ}{2} ∫_{−1}^{1} P (μ, μ′)I(τ, μ′)dμ′.  (6.4)

右辺第 1 項は, 大気の射出項で,

S_{e}(τ, μ) = 2π(1 − ϖ)B(T (τ)).  (6.5)

右辺第 2 項は, 散乱項である.

6.3.1 Toon et al. (1989) による Two-Stream Source Function 法

Toon et al. (1989) の Two-Stream Source Function 法では, 放射源関数 S(τ, μ) 内の I(τ, μ) に, Two-Stream 法により計算された放射場を与える. Two-Stream 法によるフラックスを F_{T S}(τ) とすると,

I^{±}(τ, μ) = F_{T S}^{±}(τ)/μ_{1},  (6.6)
μ_{1} = \frac{1}{2}.  (6.7)

第 n 層のフラックスは,

F_{T S,n}^{+}(τ) = Y_{1,n} (exp[−λ_{n}(τ_{n} − τ)] + Γ_{n} exp[−λ_{n}τ])  (6.8)
+Y_{2,n} (exp[−λ_{n}(τ_{n} − τ)] − Γ_{n} exp[−λ_{n}τ]) + C_{n}^{+}(τ).
F_{T S,n}^{−}(τ) = Y_{1,n} (Γ_{n} exp[−λ_{n}(τ_{n} − τ)] + exp[−λ_{n}τ])  (6.9)
+Y_{2,n} (Γ_{n} exp[−λ_{n}(τ_{n} − τ)] − exp[−λ_{n}τ]) + C_{n}^{−}(τ).

位相関数 P (μ, μ′) は,

P (μ, μ′) = \{ \begin{matrix}1 + g & (前方散乱) \\ 1 − g & (後方散乱) \end{matrix}  (6.10)

μ > 0 の時,

\frac{ϖ}{2} ∫_{−1}^{1} P (μ, μ′)I(τ, μ′)dμ′ = \frac{ϖ}{2} ∫_{−1}^{0} P (μ, μ′)\frac{F_{T S}^{−}(τ)}{μ_{1}}dμ′ + \frac{ϖ}{2} ∫_{0}^{1} P (μ, μ′)\frac{F_{T S}^{+}(τ)}{μ_{1}}dμ′(6.11)
= \frac{ϖ}{2} (1 − g)\frac{F_{T S}^{−}(τ)}{μ_{1}} + \frac{ϖ}{2} (1 + g)\frac{F_{T S}^{+}(τ)}{μ_{1}}.  (6.12)

μ < 0 の時,

\frac{ϖ}{2} ∫_{−1}^{1} P (μ, μ′)I(τ, μ′)dμ′ = \frac{ϖ}{2} ∫_{−1}^{0} P (μ, μ′)\frac{F_{T S}^{−}(τ)}{μ_{1}}dμ′ + \frac{ϖ}{2} ∫_{0}^{1} P (μ, μ′)\frac{F_{T S}^{+}(τ)}{μ_{1}}dμ′(6.13)
= \frac{ϖ}{2} (1 + g)\frac{F_{T S}^{−}(τ)}{μ_{1}} + \frac{ϖ}{2} (1 − g)\frac{F_{T S}^{+}(τ)}{μ_{1}}.  (6.14)

よって,

\frac{ϖ_{n}}{2} ∫_{−1}^{1} P_{n}(μ, μ′)I_{n}(τ, μ′)dμ′ = \frac{ϖ_{n}}{2μ_{1}} (F_{T S,n}^{+}(τ) + F_{T S,n}^{−}(τ))  (6.15)

+\frac{ϖ_{n}g}{2μ_{1}} (F_{T S,n}^{+}(τ) − F_{T S,n}^{−}(τ))

= ϖ_{n}(F_{T S,n}^{−}(τ) + F_{T S,n}^{+}(τ))  (6.16)

+ϖ_{n}g(F_{T S,n}^{+}(τ) − F_{T S,n}^{−}(τ))

= ϖ_{n}(Y_{1,n} + Y_{2,n})(1 + Γ_{n}) exp[−λ_{n}(τ_{n} − τ)]  (6.17)

+ϖ_{n}(Y_{1,n} − Y_{2,n})(1 + Γ_{n}) exp[−λ_{n}τ]

+ϖ_{n}2π(B_{0,n} + B_{1,n}τ)

±ϖ_{n}g(Y_{1,n} + Y_{2,n})(1 − Γ_{n}) exp[−λ_{n}(τ_{n} − τ)] ∓ϖ_{n}g(Y_{1,n} − Y_{2,n})(1 − Γ_{n}) exp[−λ_{n}τ]

±ϖ_{n}g \frac{2πB_{1,n}}{γ_{1,n} + γ_{2,n}} .

(ここで, ±, ∓ の符号は, 上が μ > 0 下が, μ < 0. 以下同様.) さらに,

S_{e,n}(τ, μ) = 2π(1 − ϖ_{n})B(T (τ)) = 2π(1 − ϖ_{n})(B_{0,n} + B_{1,n}τ)  (6.18)

であるから,

S_{n}^{±}(τ, ±μ) = S_{e,n}(τ, μ) + \frac{ϖ_{n}}{2} ∫_{−1}^{1} P_{n}(μ, μ′)I_{n}(τ, μ′)dμ′  (6.19)
= ϖ_{n}(Y_{1,n} + Y_{2,n})(1 + Γ_{n}) exp[−λ_{n}(τ_{n} − τ)]  (6.20)

+ϖ_{n}(Y_{1,n} − Y_{2,n})(1 + Γ_{n}) exp[−λ_{n}τ]

+2π(B_{0,n} + B_{1,n}τ)

±ϖ_{n}g(Y_{1,n} + Y_{2,n})(1 − Γ_{n}) exp[−λ_{n}(τ_{n} − τ)] ∓ϖ_{n}g(Y_{1,n} − Y_{2,n})(1 − Γ_{n}) exp[−λ_{n}τ]

±ϖ_{n}g \frac{2πB_{1,n}}{γ_{1,n} + γ_{2,n}} .

よって, Two-Stream Source Function 法による第 n 層の上端からの上向き放射 I_{n}^{+}(0, μ), 下端からの下向き放射 I_{n}^{−}(τ_{n}, −μ) は, それぞれ以下のようになる.

I_{n}^{+}(0, μ) = I_{n}^{+}(τ_{n}, μ) exp[−τ_{n}/μ] + ∫_{0}^{τ_{n}} S_{n}^{+}(t, μ) exp[−t/μ]\frac{dt}{μ}  (6.21)
= I_{n}^{+}(τ_{n}, μ) exp[−τ_{n}/μ]  (6.22)

+\frac{G_{n}}{(λ_{n}μ − 1)} (exp[−τ_{n}/μ] − exp[−τ_{n}λ_{n}])

+\frac{H_{n}}{(λ_{n}μ + 1)} (1 − exp[−τ_{n}(λ_{n} + 1/μ)])

+α_{1} (1 − exp[−τ_{n}/μ]) + α_{2} (μ − (τ_{n} + μ) exp[−τ_{n}/μ]) .

I_{n}^{−}(τ_{n}, −μ) = I_{n}^{−}(0, −μ) exp[−τ_{n}/μ] + ∫_{0}^{τ_{n}} S_{n}^{−}(t, −μ) exp[−(τ_{n} − t)/μ]\frac{dt}{μ} (6.23)
= I_{n}^{−}(0, −μ) exp[−τ_{n}/μ]  (6.24)

+\frac{J_{n}}{(λ_{n}μ + 1)} (1 − exp[−τ_{n}(λ_{n} + 1/μ)])

+\frac{K_{n}}{(λ_{n}μ − 1)} (exp[−τ_{n}/μ] − exp[−τ_{n}λ_{n}])

+σ_{1} (1 − exp[−τ_{n}/μ]) + σ_{2} (μ exp[−τ_{n}/μ] + τ_{n} − μ)) .

ここで,

G_{n} = (Y_{1,n} + Y_{2,n})ϖ_{n}(1 + g + Γ_{n}(1 − g))  (6.25)
H_{n} = (Y_{1,n} − Y_{2,n})ϖ_{n}(1 − g + Γ_{n}(1 + g))  (6.26)
J_{n} = (Y_{1,n} + Y_{2,n})ϖ_{n}(1 − g + Γ_{n}(1 + g))  (6.27)
K_{n} = (Y_{1,n} − Y_{2,n})ϖ_{n}(1 + g + Γ_{n}(1 − g))  (6.28)
α_{1} = 2π (B_{0,n} + \frac{ϖ_{n}gB_{1,n}}{γ_{1,n} + γ_{2,n}} )  (6.29)
α_{2} = 2πB_{1,n}  (6.30)
σ_{1} = 2π (B_{0,n} − \frac{ϖ_{n}gB_{1,n}}{γ_{1,n} + γ_{2,n}} )  (6.31)
σ_{2} = 2πB_{1,n}  (6.32)

さらに,

ϖ_{n}(1 + g) = 2 − γ_{1,n}  (6.33)
ϖ_{n}(1 − g) = γ_{2,n}  (6.34)
γ_{2,n}Γ_{n} = γ_{1,n} − λ_{n}  (6.35)
γ_{2,n}/Γ_{n} = γ_{1,n} + λ_{n}  (6.36)
μ_{1} = 1/2  (6.37)

であることに注意すると,

G_{n} = (Y_{1,n} + Y_{2,n})ϖ_{n}(1 + g + Γ_{n}(1 − g))  (6.38)
= (Y_{1,n} + Y_{2,n}) (2 − γ_{1,n} + γ_{2,n}Γ_{n})  (6.39)
= (Y_{1,n} + Y_{2,n}) (2 − γ_{1,n} + γ_{1,n} − λ_{n})  (6.40)
= (Y_{1,n} + Y_{2,n}) (1/μ_{1} − λ_{n})  (6.41)

H_{n} = (Y_{1,n} − Y_{2,n})ϖ_{n}(1 − g + Γ_{n}(1 + g))  (6.42)
= (Y_{1,n} − Y_{2,n}) (γ_{2,n} + (2 − γ_{1,n})Γ_{n})  (6.43)
= (Y_{1,n} − Y_{2,n})Γ_{n} (γ_{2,n}/Γ_{n} + 2 − γ_{1,n})  (6.44)
= (Y_{1,n} − Y_{2,n})Γ_{n} (γ_{1,n} + λ_{n} + 2 − γ_{1,n})  (6.45)
= (Y_{1,n} − Y_{2,n}) (λ_{n} − 1/μ_{1})  (6.46)

同様に,

J_{n} = (Y_{1,n} + Y_{2,n})Γ_{n} (λ_{n} + 1/μ_{1})  (6.47)
K_{n} = (Y_{1,n} − Y_{2,n}) (1/μ_{1} − λ_{n})  (6.48)

また,

α_{1} = 2π (B_{0,n} + \frac{ϖ_{n}gB_{1,n}}{γ_{1,n} + γ_{2,n}} )  (6.49)
= 2π (B_{0,n} + \frac{(1 − (γ_{1,n} + γ_{2,n})/2)B_{1,n}}{γ_{1,n} + γ_{2,n}} )  (6.50)
= 2π (B_{0,n} + B_{1,n} (\frac{1}{γ_{1,n} + γ_{2,n}} − \frac{1}{2}))  (6.51)
= 2π (B_{0,n} + B_{1,n} (\frac{1}{γ_{1,n} + γ_{2,n}} − μ_{1}))  (6.52)

(6.53)

同様に,

σ_{1} = 2π (B_{0,n} − B_{1,n} (\frac{1}{γ_{1,n} + γ_{2,n}} − μ_{1}))  (6.54)

6.3.2 放射源の放射場を self consistent に求める

放射源関数内の I(τ.μ) と放射伝達方程式を解いて得られる放射場は整合的でなければならない. ここでは, self consistent な放射場を得るために, iterative に放射場を求めることを想定して放射伝達方程式を書き下す.

上向き放射 I_{n}^{+}(0, μ)

第 n 層上端から μ の角度に射出する放射 I_{n}^{+}(0, μ) は, 以下のように書ける.

I_{n}^{+}(0, μ) = I_{n}^{+}(τ_{n}, μ) exp[−τ_{n}/μ] + ∫_{0}^{τ_{n}} S_{n}^{+}(t, μ) exp[−t/μ]\frac{dt}{μ}.  (6.55)

S_{n}(τ, μ) = 2π(1 − ϖ_{n})B(T (τ_{n})) + \frac{ϖ_{n}}{2} ∫_{−1}^{1} P_{n}(μ, μ′)I_{n}(τ, μ′)dμ′.  (6.56)

今, 放射源関数内の放射場を I_{s,n}(τ, μ) で与える. 散乱には半球平均近似を用いると, μ > 0 の条件で,

\frac{ϖ_{n}}{2} ∫_{−1}^{1} P_{n}(μ, μ′)I(τ, μ′)dμ′ = \frac{ϖ_{n}}{2} ∫_{−1}^{0} P_{n}(μ, μ′)I_{s,n}^{−}(τ, −μ′)dμ′  (6.57)

+\frac{ϖ_{n}}{2} ∫_{0}^{1} P_{n}(μ, μ′)I_{s,n}^{+}(τ, μ′)dμ′

= \frac{ϖ_{n}}{2} (1 − g) ∫_{−1}^{0} I_{s,n}^{−}(τ, −μ′)dμ′  (6.58)

+\frac{ϖ_{n}}{2} (1 + g) ∫_{0}^{1} I_{s,n}^{+}(τ, μ′)dμ′.

ここで, 天頂角積分の評価にはガウス求積法を用いる. ω_{i} をガウスの重み関数として,

\frac{ϖ_{n}}{2} ∫_{−1}^{1} P_{n}(μ, μ′)I(τ, μ′)dμ′ ∼ \frac{ϖ_{n}}{2} (1 − g) ∑_{i} ω_{i}I_{s,n}^{−}(τ, −μ′_{i})  (6.59)

+\frac{ϖ_{n}}{2} (1 + g) ∑_{i} ω_{i}I_{s,n}^{+}(τ, μ′_{i})

さてここで, 第 n 層内の放射強度を, 層の上端と下端の放射強度の光学的厚さに対する線型な関数で近似すると,

I_{s,n}^{±}(τ, ±μ) = I_{0,n}^{±}(±μ) + I_{1,n}^{±}(±μ)τ  (6.60)

ここで,

I_{0,n}^{±}(±μ) = I_{s,n}^{±}(0, ±μ)  (6.61)
I_{1,n}^{±}(±μ) = (I_{s,n}^{±}(τ_{n}, ±μ) − I_{s,n}^{±}(0, ±μ)) /τ_{n}.  (6.62)

以上から, S_{n}^{+}(τ, μ) は, μに依存しない関数となることが分かる.従って S_{n}^{+}(τ, μ) = S_{n}^{+}(τ) として,

S_{n}^{+}(τ) ∼ 2π(1 − ϖ_{n})(B_{0,n} + B_{1,n}τ)  (6.63)

+\frac{ϖ_{n}}{2} (1 − g) ∑_{i} ω_{i} (I_{0,n}^{−}(−μ′_{i}) + I_{1,n}^{−}(−μ′_{i})τ)

+\frac{ϖ_{n}}{2} (1 + g) ∑_{i} ω_{i} (I_{0,n}^{+}(μ′_{i}) + I_{1,n}^{+}(μ′_{i})τ)

= 2π(1 − ϖ_{n})B_{0,n} + 2π(1 − ϖ_{n})B_{1,n}τ  (6.64)

+\frac{ϖ_{n}}{2} ∑_{i} ((1 − g)ω_{i}I_{0,n}^{−}(−μ′_{i}) + (1 + g)ω_{i}I_{0,n}^{+}(μ′_{i}))

+\frac{ϖ_{n}}{2} ∑_{i} ((1 − g)ω_{i}I_{1,n}^{−}(−μ′_{i}) + (1 + g)ω_{i}I_{1,n}^{+}(μ′_{i})) τ

≡ U_{n} + V_{n}τ + W_{n}^{+} + X_{n}^{+}τ  (6.65)

ここで,

U_{n} = 2π(1 − ϖ_{n})B_{0,n}  (6.66)
V_{n} = 2π(1 − ϖ_{n})B_{1,n}  (6.67)
W_{n}^{±} = \frac{ϖ_{n}}{2} ∑_{i} \{(1 ± g)ω_{i}I_{0,n}^{+}(μ′_{i}\} + (1 ∓ g)ω_{i}I_{0,n}^{−}(−μ′_{i})\}  (6.68)
= \frac{ϖ_{n}}{2} \{(1 ± g)L_{n}^{+}(0) + (1 ∓ g)L_{n}^{−}(0)\}  (6.69)
X_{n}^{±} = \frac{ϖ_{n}}{2} ∑_{i} \{(1 ± g)ω_{i}I_{1,n}^{+}(μ′_{i}) + (1 ∓ g)ω_{i}I_{1,n}^{−}(−μ′_{i})\}  (6.70)
= \frac{ϖ_{n}}{2} \{(1 ± g)(L_{n}^{+}(τ_{n}) − L_{n}^{+}(0)) + (1 ∓ g)(L_{n}^{−}(τ_{n}) − L_{n}^{−}(0))\} /τ_{n}  (6.71)
L_{n}^{±}(τ) = ∑_{i} ω_{i}I_{s,n}^{±}(τ, ±μ_{i})  (6.72)

以上より,

I_{n}^{+}(0, μ) = I_{n}^{+}(τ_{n}, μ) exp[−τ_{n}/μ] + ∫_{0}^{τ_{n}} S_{n}^{+}(t, μ) exp[−t/μ]\frac{dt}{μ}  (6.73)
∼ I_{n}^{+}(τ_{n}, μ) exp[−τ_{n}/μ] + (U_{n} + W_{n}^{+})(1 − exp[−τ/μ])  (6.74)
+(V_{n} + X_{n}^{+}) {μ − (μ + τ_{n}) exp[−τ_{n}/μ]} .

下向き放射 I_{n}^{−}(τ_{n}, −μ)

上向き放射と同様の議論により, 下向き放射 I_{n}^{−}(τ_{n}, −μ) の式を導く.

第 n 層下端から μ の角度に射出する放射 I_{n}^{−}(τ_{n}, −μ) は, 以下のように書ける.

I_{n}^{−}(τ_{n}, −μ) = I_{n}^{−}(0, −μ) exp[−τ_{n}/μ] + ∫_{0}^{τ_{n}} S_{n}^{−}(t, −μ) exp[−(τ_{n} − t)/μ]\frac{dt}{μ}.  (6.75)

S_{n}^{−}(τ, −μ) = 2π(1 − ϖ_{n})B(T (τ_{n})) + \frac{ϖ_{n}}{2} ∫_{−1}^{1} P_{n}(−μ, μ′)I_{n}(τ, μ′)dμ′  (6.76)
= 2π(1 − ϖ_{n})(B_{0,n} + B_{1,n}τ)  (6.77)

+\frac{ϖ_{n}}{2} ∫_{−1}^{0} P_{n}(−μ, μ′)I_{s,n}^{−}(τ, −μ′)dμ′

+\frac{ϖ_{n}}{2} ∫_{0}^{1} P_{n}(−μ, μ′)I_{s,n}^{+}(τ, μ′)dμ′

= 2π(1 − ϖ_{n})(B_{0,n} + B_{1,n}τ)  (6.78)

+\frac{ϖ_{n}}{2} (1 + g) ∫_{−1}^{0} I_{s,n}^{−}(τ, −μ′)dμ′

+\frac{ϖ_{n}}{2} (1 − g) ∫_{0}^{1} I_{s,n}^{+}(τ, μ′)dμ′

∼ 2π(1 − ϖ_{n})(B_{0,n} + B_{1,n}τ)  (6.79)

+\frac{ϖ_{n}}{2} (1 + g) ∑_{i} ω_{i}I_{s,n}^{−}(τ, −μ′_{i})

+\frac{ϖ_{n}}{2} (1 − g) ∑_{i} ω_{i}I_{s,n}^{+}(τ, μ′_{i})

= 2π(1 − ϖ_{n})(B_{0,n} + B_{1,n}τ)  (6.80)

+\frac{ϖ_{n}}{2} (1 + g) ∑_{i} ω_{i} \{I_{0,n}^{−}(−μ′_{i}) + I_{1,n}^{−}(−μ′_{i})τ\}

+\frac{ϖ_{n}}{2} (1 − g) ∑_{i} ω_{i} \{I_{0,n}^{+}(μ′_{i}) + I_{1,n}^{+}(μ′_{i})τ\}

= 2π(1 − ϖ_{n})B_{0,n} + 2π(1 − ϖ_{n})B_{1,n}τ  (6.81)

+\frac{ϖ_{n}}{2} ∑_{i} \{(1 + g)ω_{i}I_{0,n}^{−}(−μ′_{i}) + (1 − g)ω_{i}I_{0,n}^{+}(μ′_{i})\}

+\frac{ϖ_{n}}{2} ∑_{i} \{(1 + g)ω_{i}I_{1,n}^{−}(−μ′_{i}) + (1 − g)ω_{i}I_{1,n}^{+}(μ′_{i})\} τ

= U_{n} + V_{n}τ + W_{n}^{−} + X_{n}^{−}τ  (6.82)

以上より,

I_{n}^{−}(τ_{n}, −μ) = I_{n}^{−}(0, −μ) exp[−τ_{n}/μ] + ∫_{0}^{τ_{n}} S_{n}^{−}(t, −μ) exp[−(τ_{n} − t)/μ]\frac{dt}{μ} (6.83)
∼ I_{n}^{−}(0, −μ) exp[−τ_{n}/μ] + (U_{n} + W_{n}^{−})(1 − exp[−τ/μ])  (6.84)

+(V_{n} + X_{n}^{−}) {τ_{n} − μ(1 − exp[−τ_{n}/μ])} .

iterative に放射場を求める

放射源関数法によって求められた第 n 層からの上向き, 下向き放射 I_{n}^{+}(0, μ), I_{n}^{−}(τ_{n}, −μ) は, 本来放射源関数 S_{n}(τ, μ) と整合的でなければならない. しかし, 放射源関数法で求められた放射場は必ずしもそうなっていない. そこで, I_{n}^{+}(0, μ), I_{n}^{−}(τ_{n}, −μ) を放射源関数とし, 再度放射伝達を解くことを繰り返すことで, 整合的な放射場を得ることを考える.

今, 放射場 I_{s,n,l}(τ, μ) をもとにした放射源関数 S_{n,l}(τ, μ) を使って計算された放射場を, I_{n,l}(τ, μ) とする. すなわち上向き放射を例に取れば,

I_{n,l}^{+}(0, μ) = I_{n,l}^{+}(τ_{n}, μ) exp[−τ_{n}/μ] + ∫_{0}^{τ_{n}} S_{n,l}^{+}(t, μ) exp[−t/μ]\frac{dt}{μ}.  (6.85)

ここで,

I_{s,n,l+1}^{+}(τ, μ) = I_{n,l}^{+}(τ, μ)  (6.86)

として, 再度放射伝達計算を行う. そして,

I_{n,l+1}^{+}(τ, μ) = I_{n,l}^{+}(τ, μ)  (6.87)

となるまで計算を繰り返すことで, 整合的な放射場を得ることができる. 下向き放射についても同様である. ちなみに, 最初の放射源関数は, Two-Stream 法による放射場を与える. すなわち,

I_{s,n,1}^{±}(τ, ±μ) = F_{T S,n}^{±}(τ, ±μ)/μ_{1}.  (6.88)

第7章{相関 k 分布法による吸収断面積, 散乱 \atop 断面積の計算}

7.1 はじめに

GCM での放射計算のための相関 k 分布法による吸収断面積, 散乱断面積の計算について述べる.

7.2 吸収物質が 1 種類の場合

7.2.1 吸収断面積の計算

大気層中の吸収物質が 1 種類の分子種の場合を考える.

今, line-by-line 法により, 波数 ν での吸収断面積 σ(ν) は既知とする. k 分布法のバンド幅 (ν_{min} - ν_{max}) を決め, そのバンド内で吸収断面積を単調増加関数である分布関数 σ_{g}(g) で表す.ただし , g は積算確率である. このときバンドの平均透過率 T̄ は

T̄ = \frac{1}{ν_{min} − ν_{max}} ∫_{ν_{min}}^{ν_{max}} exp[−σ(ν)u]dν  (7.1)
= ∫_{0}^{1} exp[σ_{g}(g)u]dg  (7.2)

ここで, u は大気層中の吸収物質のカラム量である. ここで吸収断面積を積算確率に対して離散化すれば, 近似的に

T̄ ≃ ∑_{0}^{1} exp[σ_{gi}u]Δg_{i}  (7.3)

と表すことができる. ここで σ_{gi}はΔg_{i} での吸収断面積である. σ_{gi} は透過率が変化しないという条件:

∫_{∑_{1}^{i−1} Δg_{i}}^{∑_{1}^{i} Δg_{i}} exp[σ_{g}(g)u]dg = exp[σ_{gi}u]Δg_{i}  (7.4)

が成り立っているとすれば,

σ_{gi} = −\frac{1}{u} log [\frac{1}{Δg_{i}} ∫_{∑_{1}^{i−1} Δg_{i}}^{∑_{1}^{i} Δg_{i}} exp[σ_{g}(g)u]dg]  (7.5)

となる.

この σ_{gi}を適当な圧力 p, 温度 T, 体積混合比 x の範囲でバンドごとに計算し, σ_{gi}(p, T, x) のテーブルを作成する.

ここで σ_{gi} の計算には, u, Δg_{i} を与える必要がある. これらは, 計算しようとしている大気を想定して与える. 今, 想定している大気の典型的なプロファイルを用意し, 圧力レベルが p_{j} で与えられたとする. (p_{a}, T_{a}, x_{a}) での吸収断面積を計算しようとする場合, p_{j} < p_{a} < p_{j−1} となる圧力レベルから Δp_{j} = p_{j−1} − p_{j} を求めれば, このときのカラム量 u_{a} は

u_{a} = x_{a} \frac{Δp_{j}}{g} \frac{N_{a}}{M̄}  (7.6)

となる. ここで, g は重力加速度, N_{a}はアヴォガドロ数, M̄は平均分子量である.

積算確率の離散化は, まず参照圧力, 温度, 体積混合比を決め, そのときの吸収断面積をもとに行う. 今, この吸収断面積をバンド内で吸収断面積が単調増加関数である分布関数 σ_{g}(g) で並べ替え, 吸収断面積を対数等間隔に N 等分する. このとき離散化された積算確率は

σ_{g} (∑_{i=1}^{l} Δg_{i}) = exp [log (\frac{σ_{g}(1)}{σ_{g}(0)}) \frac{l}{N} + log(σ_{g}(0))] .  (7.7)

さらに, 今積算確率を L(L ≤ N) 領域に離散化しようとするとき, Δg_{i} の組み合わせで得ることができる L の場合の数は,_{N−1}C_{L−1} である. これらすべての積算確率の組み合わせに対し, 吸収断面積を計算し, テーブルの作成を行う. 尚, 相関 k 分布法を仮定しているので, 参照圧力, 温度, 体積混合比以外の吸収断面積も離散かされた積算確率は参照条件の場合と同じである. このようにして得られた_{N−1}C_{L−1} 通りの吸収係数セットに対し, 今想定している大気の典型的なプロファイルのもとで line-by-line 放射計算を行った場合の加熱率との比較を各バンドごとに行い, line-by-line 計算との誤差が最も小さいものをバンドごとに選び出す.

errHR_{c,b} = ∑_{j=1}^{jmax} W_{j}|\frac{HR_{c,b,j} − HR_{line,j}}{HR_{line,j}} |  (7.8)

ここで, errHR_{c,b}はバンド b, 積算確率が c の場合の加熱率の誤差, HR_{c,b,j}は k 分布モデルによる加熱率, HR_{line,b,j}はline-by-line モデルによる加熱率, W_{j} は重み関数である. バンドごとに errHR_{c,b} が最も小さくなる積算確率の組み合わせ c を選ぶ.

7.2.2 散乱断面積の計算

分子による散乱を想定し, 可視光領域より長波長領域を考えれば, バンド内 (ν_{min} 〜 ν_{max}) での散乱断面積は単調な関数であると考えられる. 散乱断面積の波数空間での増減は吸収断面積のそれとは無関係であり, またバンド内での散乱断面積の大小の差は, 吸収断面積のそれに比べて小さいと仮定し, バンド内で吸収断面積は 1 つの値 σ_{s} で代表させる. 式 (7.5) を参考にすれば, このときの散乱断面積 σ_{s} は以下のようになる.

σ_{s} = −\frac{1}{u} log [∫_{0}^{1} exp[σ_{g}(g)u]dg]  (7.9)

波数空間での散乱断面積は line-by-line 計算により得られ, 散乱体の大気層中のコラム量 u は, 吸収断面積の計算の場合と同様に与えれば, 上式よりバンドごとの散乱断面積を得ることができる.

第8章その他

8.1 断熱温度勾配

8.1.1 乾燥断熱減率

単位質量あたりの空気塊に ΔQ の熱が加えられた時, ΔW の仕事がされ, Δu の内部エネルギー変化があったとすると, 熱力学第 1 法則から

ΔQ = ΔW + Δu.  (8.1)

なされた仕事 ΔW は圧力 p, 体積変化 (比容の変化) Δα を用いて

ΔW = pΔα.  (8.2)

また理想気体を仮定すれば, 内部エネルギーの変化は温度のみの関数となるので, 定積比熱 C_{v}, 温度変化 ΔT を用いて

Δu = C_{v}ΔT.  (8.3)

よって熱力学第一法則は

ΔQ = pΔα + C_{v}ΔT  (8.4)

となる.

ここで単位質量あたりの理想気体の状態方程式 (pα = RT, R は単位質量あたりの気体定数 (R = R^{∗}/M[JK^{−1}kg^{−1}]) において, 圧力が p + Δp, 比容が α + Δα のとき, 温度が T + ΔT となったとすれば

(p + Δp)(α + Δα) = R(T + ΔT).  (8.5)

2 次の微少量を無視すれば,

\frac{Δp}{p} + \frac{Δα}{α} = \frac{ΔT}{T} .  (8.6)

これを熱力学第一法則に代入すれば

ΔQ = C_{p}ΔT − αΔp  (8.7)

となる. ただし, C_{p} = C_{v} + R である.

今, 断熱温度勾配を考えているので ΔQ = 0 とすれば,

C_{p}ΔT = αΔp.  (8.8)

よって圧力に対する乾燥断熱減率は

\frac{ΔT}{Δp} = \frac{α}{C_{p}}  (8.9)
= \frac{RT}{pC_{p}}  (8.10)

8.1.2 湿潤断熱減率

水蒸気を含む大気の湿潤断熱減率を導出する.

ここでは, 混合比 w を乾燥空気の密度 ρ_{d} に対する水蒸気密度 ρ_{v}(単位容積中に含まれる水蒸気質量) とする. 水蒸気の分圧を e, 湿った空気の圧力を p とすれば, 乾燥空気の分圧は p − e である. 水蒸気も理想気体の状態方程式に従うと仮定すれば,

e = R_{v}ρ_{v}T.  (8.11)

ここで R_{v} = R^{∗}/M_{H_{2}O}. 同様に乾燥空気の状態方程式は

p − e = R_{d}ρ_{d}T  (8.12)

であるから, 混合比 w は

w = \frac{e}{p − e} \frac{R_{d}}{R_{v}} = ϵ\frac{e}{p − e}.  (8.13)

ここで ϵ = R_{d}/R_{v} である. さらに水蒸気が飽和している場合には, 水蒸気の分圧 e は飽和蒸気圧 e_{s} で置き換わる. このとき混合比も w_{s} と書くことにする. すなわち

w_{s} = ϵ\frac{e_{s}}{p − e_{s}} .  (8.14)

今, 乾燥空気 1kg の質量をもち, 水蒸気で飽和した空気塊を断熱的に変化させたとする. 空気塊の乾燥空気分圧を p_{d}(= p − e_{s}) としたとき, 分圧, 温度はそれぞれ p_{d}, T から p_{d} + Δp_{d}, T + ΔT に変化し, 飽和混合比も w_{s}から w_{s} + Δw_{s} に変化したとする. この −Δw_{s} 分の水蒸気が凝結 (あるいは昇華) し, −LΔw_{s} の熱が放出される (L は単位質量あたりの潜熱).

w_{s} + Δw_{s} = ϵ\frac{e_{s} (1 + \frac{Δe_{s}}{e_{s}} )}{p_{d} (1 + \frac{Δp_{d}}{p_{d}} )} ≃ w_{s} (1 + \frac{Δe_{s}}{e_{s}} ) (1 − \frac{Δp_{d}}{p_{d}} )  (8.15)

であるから, 近似的に

Δw_{s} = w_{s} (\frac{Δe_{s}}{e_{s}} − \frac{Δp_{d}}{p_{d}} ) .  (8.16)

ここで, e_{s} は温度のみの関数であるから, 以下のように書ける.

Δw_{s} = w_{s} (\frac{1}{e_{s}} \frac{Δe_{s}}{ΔT} ΔT − \frac{Δp_{d}}{p_{d}} ) .  (8.17)

凝結による潜熱は空気塊の温度上昇のみに寄与すると考えれば (水滴などは温めないと考えれば), 熱力学第一法則で ΔQ = −LΔw_{s} として,

−Lw_{s} (\frac{1}{e_{s}} \frac{Δe_{s}}{ΔT} ΔT − \frac{Δp_{d}}{p_{d}} ) = C_{p}ΔT − αΔp.  (8.18)

これを整理すれば,

\frac{ΔT}{Δp} = \frac{\frac{ϵLe_{s}}{(p−e_{s})^{2}} + α}{C_{p} + \frac{ϵLp}{(p−e_{s})^{2}} \frac{Δe_{s}}{ΔT}} .  (8.19)

ちなみに, Nakajima et al. (1992) 式 (4) は, 上式で飽和水蒸気混合比は全圧に対して小さい (p − e_{s} ∼ p) としたものに対応している.

8.2 潜熱

潜熱 L をクラウジウス=クラペイロンの式から導出する. 今, 液層 (もしくは固層) と気層が層平衡で接している状態を考える. つまり系の圧力は飽和蒸気

圧 e_{s} となっている. 液層 (もしくは固層), 気層の単位質量あたりの水蒸気の体積をそれぞれ, v_{l}, v_{g} とすれば, クラウジウス=クラペイロンの式は以下のようになる.

\frac{de_{s}}{dT} = \frac{L}{T (v_{g} − v_{l})} ∼ \frac{L}{T v_{g}} .  (8.20)

理想気体を仮定すれば v_{g} = R_{v}T/e_{s} なので,

L = \frac{R_{d}T^{2}}{e_{s}} \frac{de_{s}}{dT} .  (8.21)

8.3 H_{2}O の状態方程式

8.4 中心星のスペクトル

8.5 地表面のアルベド