第 1 章

EBM を用いた太陽定数増減実験

部分凍結解の初期値依存性を調べるために, 南北 1 次元エネルギーバランスモデル (EBM) を用いて実験を行う. EBM は, GCM よりも過程が単純であるため部分凍結解の初期値依存性の解釈がしやすく, 計算量も少ないため様々な実験を行うことができる. 本章でははじめに, 作成した EBM について説明する. 次に, 作成した EBM を用いて部分凍結解の初期値依存性について調べる. ここでは, 水平格子点数や初期条件を変化させ, いくつかの実験を行う.

1.1 EBM の概要

本節では, 作成した EBM について詳細を説明する.

1.1.1 支配方程式系

本研究で作成したエネルギーバランスモデルは Sellers (1969) の定式化に従う. モデルの支配方程式は, 地表面温度の時間発展の式

C \frac{∂T (x, t)}{∂t} = F_{S} − F_{OLR} + \frac{∂}{∂x} D (1 − x^{2}) \frac{∂T (x, t)}{∂x},  (1.1)

\underline{EBM 結果 1}

で表される. ここで, C は熱容量, T は地表面温度, x は緯度の正弦, D は拡散係数である. また, 右辺第一項 F_{S} は太陽放射, 右辺第二項 F_{OLR} は惑星放射, 右辺第三項は水平熱輸送を表す. ここで, 太陽放射 F_{S} と惑星放射 F_{OLR} はそれぞれ

F_{S} = Qs (x) (1 − α (x))  (1.2)
F_{OLR} = A + BT (x, t)  (1.3)

である. ここで, Q は全球平均太陽放射フラックス, α はアルベド, A, B は定数である. 本研究では, Lindzen (1990) の値を参照し, A = 212.05 W m^{−2}, B = 1.55 W m^{−2} K^{−1} とした. s (x) は年平均太陽放射フラックスの子午面分布を表す関数であり,

s (x) = 1 + s_{2}p_{2} (x)  (1.4)

と近似する. ここで, s_{2} = −0.482, p_{2} (x) は 2 次のルジャンドル多項式で, p_{2} (x) = (3x^{2} − 1) /2 である. s_{2} の値は Lindzen (1990) の値を参照した.

1.1.2 境界条件

ここでは年平均の平衡状態について考えるため, 計算領域を半球対称とし, 赤道 (x = 0) から極 (x = 1) までの領域を考える. そして, 計算領域の両端において熱の出入りはないため, 境界条件は,

(\frac{∂T}{∂x} ) = 0, (x = 0) ,  (1.5)
(\frac{∂T}{∂x} ) = 0, (x = 1)  (1.6)

である.

1.1.3 支配方程式の離散化

ここでは, 1.1.1 節に示した支配方程式の離散化について説明する. 図 1.1 に示すように, 赤道から極までを N 個の緯度帯に分割する. (1.1) は下のように離散化する.

\underline{EBM 結果 2}

図 1.1: 離散モデルにおける変数配置

C \frac{T_{k}^{l+1} − T_{k}^{l}}{Δt} = Qs (x) (1 − α (x))−(A + BT_{k}^{l})+\frac{D}{Δx^{2}} [(1 − x_{i+ \frac{1}{2}}^{2} ) (T_{i+1} − T_{i}) − (1 − x_{i− \frac{1}{2}}^{2} ) (T_{i} − T_{i−1})] .

(1.7)

ここで, 上付き添字 l は時間を表し, 下付き添字 k は緯度帯を表す.

1.2 EBM における部分凍結解の初期値依存性

本節では, GCM で見られるような部分凍結解の初期値依存性が EBM においても見られるか否かを確認する. そのため, 解像度や初期条件を変化させ, いくつかの実験を行う.

1.2.1 実験設定

アルベド α (x_{i}) は,

α (x_{i}) = \{ \begin{matrix}α_{i} = 0.6 (T (x) ≤ 271.15 K) , \\ α_{g} = 0.1 (T (x) > 271.15 K) . \end{matrix}  (1.8)

を与える. つまり, 凝結温度は 271.15 K とした. 熱容量 C は 1.0 J K^{−1} とし, 拡散係数 D は 0.2 W m^{−2} K^{−1} とする. モデルの水平格子点数は 16, 50, 500 とする. ここで, 水平格子点数 16 は GCM における T21 に対応する. 与えた全球平均太陽放射フラックスは, 250 Wm^{−2} ≤ Q ≤ 500 Wm^{−2} である.

\underline{EBM 結果 3}

各実験では次に述べる初期条件から平衡に達するまでモデルを積分することで平衡状態を求める. (1.1) の右辺の絶対値が 10^{−5} W m^{−2} より小さくなると平衡状態になったと判断する.

1.2.2 初期条件

初期条件には次の 4 つの温度分布を与える (表 1.2.1, 図 1.2). 一つ目は, 全緯度で 250 K である全球凍結状態の温度分布 (図 1.2a, 青線), 二つ目は, x < 0.40 において 300 K, x ≥ 0.40 において 250 K である部分凍結状態の温度分布 (図 1.2b, 緑線), 三つ目は, x < 0.72 において 300 K, x ≥ 0.72 において 250 K である部分凍結状態の温度分布 (図 1.2b, オレンジ線), 四つ目は, 全緯度で 300 K である氷なし状態の温度分布である (図 1.2a, 赤線). 以下では, 一つ目の初期条件を用いた実験を [GIC250], 二つ目の初期条件を用いた実験を [PIC5_300-250], 三つ目の初期条件を用いた実験を [PIC11_300-250], 四つ目の初期条件を用いた実験を [NIC300] と呼ぶ (表 1.2.1).

表 1.2.1: 実験設定

実験名初期条件
[GIC250]  T (x) = 250 K (0 ≤ x ≤ 1)
[PIC5_300-250] T (x) = \begin{cases}300 K (x < 0.40) \\ 250 K (x ≥ 0.40) \end{cases}
[PIC11_300-250] T (x) = \begin{cases}300 K (x < 0.72) \\ 250 K (x ≥ 0.72) \end{cases}
[NIC300]  T (x) = 300 K (0 ≤ x ≤ 1)

1.2.3 実験結果

まず, 水平格子点数 16, Q = 300 W m^{−2} として, [GIC250], [PIC5_300-250], [PIC11_300-250], [NIC300] の実験を行った.

{(a)  [GIC250] と [NIC300] に用いた \atop 初期温度分布}{(b) [PIC5_300-250] と [PIC11_300-250] \atop に用いた初期温度分布}

図 1.2: 水平格子点数 16 のときの初期温度分布

図 1.3 に水平格子点数 16, Q = 300 W m^{−2} のときの地表面温度の時間発展を示す. 図 1.3a–1.3d を見ると, [GIC250], [PIC5_300-250], [PIC11_300-250], [NIC300] のいずれの場合においても地表面温度はある解に収束している. [GIC250] は全球凍結状態, [PIC5_300-250], [PIC11_300-250] は部分凍結状態, [NIC300] は氷なし状態の平衡状態となっており, 多重解が得られた. さらに, [PIC5_300-250], [PIC11_300-250] はどちらも部分凍結状態だが, 温度は異なっており, 結果として氷線も異なる緯度に位置する.

図 1.4 に水平格子点数 16, Q = 300 W m^{−2} のときの [PIC11_300-250] と [PIC5_300-250] から得られた地表面温度の緯度分布を示す. 2 つの結果を比較すると, 最も低緯度に位置する凍結した格子点は 2 格子点ずれている.

[PIC5_300-250] と [PIC11_300-250] において最も低緯度に位置する凍結した格子点が 3 格子点ずれることを詳しく調べるために, [PIC5_300-250] と [PIC11_300-250] だけでなく, 新たに 3 つの初期値を追加し (図 1.5 の左図の黒線, 赤線, ピンク線) 実験を行った. 図 1.5 の右図は図 1.4 に追加実験の結果を追加したものである. 図 1.5 の右図より, 最も低緯度に位置する凍結した格子点が 1 格子点ずつずれた部分凍結解が存在することが分かる.

次に氷線緯度の太陽定数依存性を調べるために, 250 W m^{−2} ≤ Q ≤ 550 W m^{−2} の範

(a) [GIC250]  (b) [PIC5_300-250]

(c) [PIC11_300-250]  (d) [NIC300]

{図 1.3: 水平格子点数 16, Q = 300 W m^{−2} のときの地表面温度の時間発展 (紫: 0 秒, 青: 5000 秒, \atop 水色: 10000 秒, 緑: 25000 秒, 黄: 50000 秒, オレンジ: 75000 秒, 赤: 平衡時)}

囲で Q を 1 W m^{−2} きざみで変化させて実験を行った.

図 1.6 に水平格子点数 16 のときの氷線緯度と太陽定数の関係を示す. 氷線緯度は, 凝固点の前後の格子点の温度の線形補間で求めた. 全球凍結解, 部分凍結解, 氷なし解の 3 つの解や多重解が得られた. また, 初期条件によって, 2 種類の部分凍結解が存在することが確認でき, [PIC11_300-250] より [PIC5_300-250] から得られた平衡解の氷線の方が低緯度に位置している.

図 1.4: 水平格子点数 16, Q = 300 W m^{−2} のときの平衡状態における地表面温度の緯度分布. 緑色の × 印は [PIC5_300-250] の結果, オレンジ色の × 印は [PIC11_300-250] の結果である.

\begin{matrix}色 & 初期条件 \\ 黒 & T (x) = \begin{cases}300 K (x < 0.47) \\ 250 K (x ≥ 0.47) \end{cases} \\ 赤 & T (x) = \begin{cases}300 K (x < 0.53) \\ 250 K (x ≥ 0.53) \end{cases} \\ ピンク & T (x) = \begin{cases}300 K (x < 0.59) \\ 250 K (x ≥ 0.59) \end{cases} \end{matrix}^{表 1.2.2: 追加した初期温度分布の詳細}

最後に, 氷線緯度と太陽定数の関係の水平格子点数依存性について調べるため, 水平格子点数を 50, 100, 500 として同様の実験を行った. 図 1.7 は水平格子点数 16, 50, 100, 500 における氷線緯度と太陽定数の関係である. 図 1.7 より, 部分凍結解の初期値依存性の幅は水平格子点数が増えると小さくなることがわかる.

図 1.5: 左図: [PIC5_300-250] と [PIC11_300-250] の初期温度分布と追加した初期温度分布 (詳細は表 1.2.2 参照). 右図: 水平格子点数 16, Q = 300 W m^{−2} のときの平衡状態における地表面温度の緯度分布

図 1.6: 水平格子点数 16 のときの太陽定数と氷線緯度の関係. ピンクは [GIC250] の結果, 青は [PIC5_300-250] の結果, 緑は [PIC11_300-250] の結果, 赤は [NIC300] を示す.

(a) 水平格子点数 16 (b) 水平格子点数 50 (c) 水平格子点数 500

図 1.7: 太陽定数と氷線緯度の関係: (a) 水平格子点数 16, (b) 50, (c) 500. ピンクは [GIC250] の結果, 青は [PIC5_300-250] の結果, 緑は [PIC11_300-250] の結果, 赤は [NIC300] の結果を示す. (b), (c), (d) には比較のために (a) の結果を黒点で描画している.

1.2.4 まとめと考察

本節では, GCM で見られた部分凍結解の初期値依存性が EBM でも見られるか調べた. その結果, EBM においても部分凍結解の初期値依存性が存在することや, その部分凍結解の初期値依存性の幅が水平格子点数を増やすと小さくなることが分かった.

図 1.8: 水平格子点数 16, Q = 300 W m^{−2} のとき得られた平衡状態における太陽放射フラックス, 長波放射フラックス, 水平熱輸送フラックスの単位面積当たりの加熱率. 赤線は太陽放射フラックス, 青線は長波放射フラックス, 緑線は水平熱拡散フラックスの加熱率である

従って, 水平格子点数を無限に大きくすると, 部分凍結解は一意に求まるように思われる. つまり, 離散系における表現に部分凍結解の初期値依存性の原因があると想定される. そこで, (1.1) の右辺のうちどの項が特に離散化の影響を受けるのかを調べることにする. 図 1.8 に, 水平格子点数 16, Q = 300 W m^{−2} のとき得られた平衡状態における太陽放射フラックス, 長波放射フラックス, 水平熱輸送の加熱率を示す. 拡散型の水平熱輸送を用いたことで温度は連続的な緯度分布をしているため, 温度の一次関数である長波放射フラックスは連続的な分布となっている. ところが, 吸収される太陽放射フラックスと正味の水平熱輸送の加熱率は氷線近傍で不連続的な分布を示している. 水平熱輸送の加熱率は他の二項と釣り合うように決定されるため, ここでは太陽放射フラックスの加熱率に着目する.

\underline{EBM 結果 10}

(1.2) と (1.4) より, Qs(x) は連続的な分布であるが, α は不連続な分布になっている. つまり, 氷線近傍における太陽放射フラックスの加熱率が不連続な分布になるのはアルベドの分布に起因していると考えられる. そこで, 次節では離散モデルでのアルベドの表現を見直し, それによる解への影響を調べることにする.

1.3 EBM におけるアルベドの表現の影響

1.2 節より, 部分凍結解の初期値依存性は支配方程式の離散化に原因があることがわかった. また, 太陽放射フラックスに着目すると氷線近傍において不連続となっており, これはアルベドの分布に起因していると考えた. そこで本節では, 1.2 節のアルベドの表現 (図 1.10 の黒線) からサブグリッドスケールの氷の面積を考慮し連続系における表現に近づけたアルベドの表現 (図 1.10 の赤線) に変更し, 追加実験を行う. そして, その結果と 1.2 節の実験結果を比較する.

1.3.1 実験設定

本節では, サブグリッドスケールの氷の面積を考慮し, アルベド分布を考慮する. このとき, アルベド α (x_{i}) は,

α (x_{i}) = \frac{1}{Δx} ∫_{Δx} α (x) dx
= \frac{1}{Δx} ∫_{Δx} [α_{i}β (T) + α_{g} (1 − β (T))] dx  (1.9)

とする. ここで,

β (T) = \{ \begin{matrix}1 (T (x) ≤ 271.15 K) , \\ 0 (T (x) > 271.15 K) . \end{matrix}  (1.10)
T (x) = T (x_{i}) + \frac{T (x_{i+1}) − T (x_{i})}{x_{i+1} − x_{i}} (x − x_{i})  (1.11)

である. つまり, (1.11) よりアルベドは格子点の間の温度分布を線形と仮定して, 各格子内の温度が凝固点を下回る面積を推定し (図 1.9 の青い領域), 面積の重みをつけた平均によって格子点のアルベドを求める. その結果, 例えば, 水平格子点数 16, Q = 300 W m^{−2} のときの平衡状態の温度分布 (図 1.9 の赤線) におけるアルベド分布は, 図 1.10 の赤線のようになる. その他の設定は 1.2 と同じである.

図 1.9: 水平格子点数 16, Q = 300 W m^{−2} のと ; きの平衡状態における地表面温度分布. ; 青で塗られた領域はある格子内の温度が ; 凝固点を下回る領域であり, 青線で囲ま ; れ白で塗られた領域はある格子内の温度 ; が凝固点を上回る領域である. ; 図 1.10: 水平格子点数 16, Q = 300 W m^{−2} の ; ときの平衡状態におけるのアルベド分 ; 布. 黒線はアルベドの扱いを変更する ; 前, 赤線はアルベドの扱いを変更した後 ; のアルベド分布である.

1.3.2 初期条件

実験名および初期条件は 1.2 と同じである (表 1.2.1, 図 1.2).

1.3.3 実験結果

まず, 水平格子点数 16, Q = 300 W m^{−2} として, [GIC250], [PIC5_300-250], [PIC11_300-250], [NIC300] の実験を行った.

図 1.11 に水平格子点数 16, Q = 300 W m^{−2} のときの地表面温度の時間発展を示す. 図 1.11a–1.11d を見ると, [GIC250], [PIC5_300-250], [PIC11_300-250], [NIC300] のいずれの場合においてもある解に収束している. [GIC250] は全球凍結状態, [PIC5_300-250], [PIC11_300-250] は部分凍結状態, [NIC300] は氷なし状態の平衡状態となっており, 多重解が得られた. さらに, [PIC5_300-250], [PIC11_300-250] はどちらも部分凍結状態であり同じ温度となっており, 氷線も同じ緯度に位置する.

次に氷線緯度の太陽定数依存性を調べるために, 250 W m^{−2} ≤ Q ≤ 550 W m^{−2} の範囲で Q を 1 W m^{−2} きざみで変化させて実験を行った.

図 1.12 に水平格子点数 16 のときの氷線緯度と太陽定数の関係を示す. 全球凍結解, 部

ebm_conclusion.tex  2018/03/10(松田幸樹)

(a) [GIC250]  (b) [PIC5_300-250]

(c) [PIC11_300-250]  (d) [NIC300]

図 1.11: 水平格子点数 16, Q = 300 W m^{−2} のときの地表面温度の時間発展 (黒: 0 秒, 紫: 5000 秒, 青: 10000 秒, 水色: 50000 秒, 緑: 100000 秒, 黄: 150000 秒, オレンジ: 200000 秒, 赤: 平衡時)

分凍結解, 氷なし解の 3 つの解や多重解が得られた. また, アルベドの表現を変更したことによって, 初期条件によらず部分凍結解は一意に求められている.

次に, 氷線緯度と太陽定数の関係の水平格子点数依存性について調べるために, 水平格子点数を 50, 100, 500 として同様の実験を行った. 図 1.13 は水平格子点数 16, 50, 100, 500 における氷線緯度と太陽定数の関係である. これによると, 部分凍結解は水平格子点数によらずおおよそ一意に求まることが分かる.

図 1.12: 水平格子点数 16 のときの太陽定数と氷線緯度の関係. ピンクは [GIC250] の結果, 青は [PIC5_300-250] の結果, 緑は [PIC11_300-250] の結果, 赤は [NIC300] の結果を示す.

図 1.14 は各解像度ごとの氷線緯度と太陽定数の関係を 1.2 節のアルベドの扱いを変更する前の実験と 1.3 節のアルベドの扱いを変更した後の実験を比較したものである. 図 1.14 より, アルベドの扱い変更後の実験において得られる解はアルベドの扱い変更前の実験において水平格子点数を大きくしたときに収束する解とおおよそ一致する.

(a) 格子点数 16 (b) 格子点数 50 (c) 格子点数 500

図 1.13: 太陽定数と氷線緯度の関係: (a) 水平格子点数 16, (b) 50, (c) 500. ピンクは [GIC250] の結果, 青は [PIC5_300-250] の結果, 緑は [PIC11_300-250] の結果, 赤は [NIC300] の結果を示す. (b), (c), (d) には比較のために (a) の結果を黒点で描画している.

\underline{EBM 結果 16}

(a) 格子点数 16 (b) 格子点数 50 (c) 格子点数 500 図 1.14: 氷線緯度と太陽定数の関係の各解像度ごとの比較: (a) 水平格子点数 16, (b) 50, (c) 100, (d) 500. 赤はアルベドの扱いを変更する前, 青はアルベドの扱いを変更した後の結果で ある.

ebm_conclusion.tex  2018/03/10(松田幸樹)

1.3.4 まとめと考察

本節では, サブグリッドスケールの氷の面積を考慮してアルベドの表現を変更し, 部分凍結解の初期値依存性について調べた. その結果, 部分凍結解の初期値依存性は現れなかった. これは水平格子点数によらない結果である. また, アルベドの表現を変更する前の実験において水平格子点数を大きくしたときに収束した結果とアルベドの表現を変更した実験の結果はおおよそ一致した. 以上より, EBM における部分凍結解の初期値依存性はサブグリッドスケールのアルベド分布を無視したことに起因していることがわかった. そして, 本研究で用いたアルベドの表現は上手く機能しており, それを用いることで部分凍結解の初期値依存性はなくなることがわかった.

1.2 節では, アルベドの表現を変更する前の実験において, 水平格子点数を大きくすると部分凍結解の初期値依存性の幅は小さくなることがわかった. 従って, 1.2 節と同様のアルベドの表現をしている本研究の GCM においても解像度を大きくすると, 部分凍結解の初期値依存性の幅は小さくなることが期待されるが, 容易ではない. しかし, 1.3 と同様に, GCM においてもサブグリッドスケールの氷の面積を考慮したアルベド分布を与えることで第??章で見られた部分凍結解の初期値依存性は小さくなると考えられる.

1.4 初期値依存性の解釈

ここでは, アルベドの表現を変更する前になぜ EBM において部分凍結解の初期値依存性が現れたのかについて改めて考察することにする. それは, モデルの振る舞いの理解に有益であるためである. 本節では, (1.1) の左辺をゼロとした平衡状態におけるエネルギーバランスの式,

F_{S} − F_{OLR} + D \frac{∂^{2}T}{∂x^{2}} = 0  (1.12)

を用いて, この式の解の振る舞いを考える. 任意の全球平均太陽放射 Q に対して, 以下のように (1.12) を解く:

1. 始めに氷線を仮定する. この時仮定した氷線を x_{0} とする.

ebm_conclusion.tex  2018/03/10(松田幸樹)

\underline{EBM 結果 18}

2. 氷線を仮定したことによりアルベド分布が求まり, (1.12) を解くことができる.

3. (1.12) を解いて得られた温度分布から氷線を求める. この時求めた氷線を x とする.

4. このとき, x = x_{0} ならば x は (1.12) の解となる.

ここで, 物理定数の値やアルベドの与え方は 1.2 節と同様である. ただし, 上記の方法では, 小氷冠不安定や大氷冠不安定のような不安定解も解として求まることに注意しておく.

Q = 300 W m^{−2} のときの仮定した x_{0} と x の関係を図 1.15 に示す. 赤線は水平格子点数 16 のときの結果, 青線は水平格子点数 500 のときの結果である. また, 緑線は x = x_{0} の線であり, 緑線と赤線もしくは青線との交点が解となる. まず, 水平格子点数 16 のとき, x_{0} に対応する x が階段状に分布する. これは, 仮定する氷線の位置 x_{0} を変化させても, x_{0} が隣の格子点の位置に至るまではアルベドの分布に変化がなくそれに対応する x の値も変化しないためである. このため, 緑線と赤線の交点となっている部分凍結解が多数の x で実現する. 一方, 水平格子点数 500 では一定の x となる x_{0} の範囲が小さいので, 部分凍結解はおおよそ一意に求まることがわかる. 本節の実験では水平格子点数 500 を用いたが, 水平格子点数を無限にすれば連続系の表現と考えることができ, このとき x は連続な分布となると推測される. 従って, 連続系では部分凍結解は一意に求まると考えられる. また, 1.3 節のようにアルベドの表現を変更すると, x_{0} の変化に対して連続的にアルベドの分布が変化することになる. これによって, 1.3 節では部分凍結解が一意に求められたと考えられる.

GCM では, EBM には含まれていない多数の複雑な過程が取り入れられている. しかし, 本節で述べた離散化された EBM の振る舞いと同様のことは GCM においても起こっているだろう. GCM で見られた異なる氷線を持つ部分凍結解の存在には, ここで述べた機構が重要な役割を果たしていると考えられる.

\underline{EBM 結果 19}

{図 1.15: Q = 300 W m^{−2} のときの x_{0} と x の関係 (赤線は水平格子点数 16 のときの結果, 青線 \atop は水平格子点数 500 のときの結果, 緑線は x = x_{0} の線).}

ebm_conclusion.tex  2018/03/10(松田幸樹)