エネルギー収支気候モデル [例題 3] 温室効果の強さが温度とともに変わる

0次元(空間的不均一を省略)のエネルギー収支気候モデルを、 つぎのように、すこしずつ複雑にしていく。

• [[例題 1] 温室効果のないばあい]
• [[例題 2] 温室効果の強さを定数で入れる]
• [例題 3] 温室効果の強さが温度とともに変わる
• [例題 4] 惑星アルベド (太陽放射反射率) が温度とともに変わる

温室効果の強さを数量であらわす方法には、きまった標準はないが、 ここでは、例題2からつづいて、 [温室効果の基本]のページで 説明のためにつかった「ガラス板モデル」で、 大気が理想化したガラス板 n 枚に相当するとき、 大気上端から出ていく地球放射は、地表面温度での黒体放射の 1/(n+1) であることをふまえて、 この n を「温室効果の強さ」と考えることにする。 このように仮定すると、 出ていく地球放射は、[例題1]の式(6)をもとにした、[例題2] のつぎの式できまる。

F_out = σ_B T⁴ / (n+1)

温室効果の強さが温度によってかわる原因としては、 海のある地球のもとでは、水蒸気量の変化がある。 [大気中の水蒸気量の表現]のページを参照してほしい。 温度があがるとともに、飽和水蒸気圧が大きくなる。 大気の総質量はかわらないとみなせるので(こまかいことをいうと水蒸気の質量のぶんだけかわるが)、 平均気圧はかわらないとしてよい。したがって、飽和水蒸気圧に比例して飽和比湿がふえる。 実際の比湿 (空気の質量のうちの水蒸気の質量のわりあい) がどのようにふえるかは自明ではないが、 近似としては、飽和比湿に比例してふえる、つまり、相対湿度が一定である、と、想定できる。

温度から (たとえば、比湿を経由して)、 ここでいう「温室効果の強さ」 n をもとめる計算式をつくることができればよいのだが、 それは簡単ではない。

そこで、温度の1次関数 (関数のグラフでいえば、直線) による近似で考える。 温度 T から 温室効果の強さ n をきめる関係式を n = n(T) とすれば、 その関数 n(T) を T で微分した dn/dT を考えることができる。 ある特定の温度 T_ref での dn/dT の値を (このページかぎりの約束で) b としてみよう。 (この T_refを「参照される温度」という意味で reference temperatureということがある。) すると、温度が T_ref 付近にあるときの n の値を、 T の1次関数 (グラフの接線となる直線) で、次のように近似することができる。

n(T) = n(T_ref) + b × (T - T_ref)

• [注] bのあとにある「×」は省略してよいのだが、 そのつぎにある かっこ が、ひきざんをかけざんよりも優先させるためにいれたものであって、 「b が かっこ内のもの の 関数である」という意味ではないこと (「n(T)」の かっこ とは意味がちがうこと) を明示するために、念のため書いた。

ここでの T の単位は K (ケルビン) でも ℃ でもよいが、一貫していないといけない。 今回のプログラム例では、計算式の中では K 単位だが、人が読み書きする数量は ℃ 単位として、 必要に応じて換算する形で書いてみた。 全部 K でよいならば、換算するところをはずせばよい。

実際にすこしやってみると、例題1や例題2が初期値にあまりよらないふるまいをするのに対して、 例題3のモデルは、初期値にかなり敏感であることがわかった。

とくに、 低い温度ではじめると、温度がますますさがり、n が 負という、モデルの想定外の状態になることもある。 n が負になってもすぐには計算は破綻しないが、n が -1 以下またはそれに近い値になると破綻する。 そこで、プログラムでは、n のありうる最小値 (例文では 0 ) をきめておき、 計算された n がそれ以下になったら、最小値でおきかえることにする。

例題2のモデルで、「太陽定数」と地球の惑星アルベドが現実にちかい値で、 n = 1 (ガラス板モデルのガラス板1枚に相当) としたとき、定常となる温度は 約 30℃である (現実の世界平均地表温度 約 15℃ とは、この程度にちがう)。 そこで、今回の実行例では、T_ref を 30℃とし、n(T_ref) = 1, b = 0.02 としてみた (そうすると、T が -20℃で n = 0 になる。)。 初期値は、T_ref から極端にはちがわない値 (15℃～40℃ぐらい) でやることにする。 このような値をあたえてみると、 現実に近いかどうかはともかく、モデルとしては無理のないふるまいをする。

b (R言語のプログラム例では dghedtemp、表計算のワークシート例ではセル C22 の内容)の 値を 0 にすれば、実質的に 例題 2 と同じ問題の計算になる。

R言語による計算

プログラム例は、 [例題 2]のプログラム、 それはさらに、[[例題 1] R言語による計算]の プログラムを部分的に修正したものなので、 共通の点については例題1、例題2のページを見てほしい。

今回のプログラム例は、日本語文字を 文字コード UTF-8 でふくんでいる。 パソコンのOSやソフトウェアのバージョンによっては、文字コード Shift_JIS に変換してから 読む必要があるかもしれない。

• ex3.R ... 例題2の ex2.R と同様に、1つの条件でシミュレーションをおこない、 結果の温度とエネルギーの出入りのそれぞれの時系列をグラフで表示する。
  • プログラムの実質的な変更箇所はつぎのとおり。
    • 上記の説明の T_refを tempref [単位 K] および tdegcref [単位 ℃]、 b を dghedtemp (変数名は dn/dT のつもり) 、 n のありうる最小値を ghemin として、それぞれ値を代入しておく。 (dghedtempの値を変えて実行することを想定している。)
    • 「時間ステップを deltat ずつ進めて計算」の部分の中で、 (エネルギーの出入りを計算する前に) その時刻の ghe の値を計算する。
    • 計算された ghe の値が ghemin よりも小さいときは、ghemin の値を代入する。
    • グラフの表題に、ghe0, dghedtemp, tdegcref の値をふくめる。
    • グラフの縦軸の範囲の指定をふくめる。 (これは結果の数値の範囲に応じて手動で変えて実行することを想定している。)

実行例では、参照温度を 30 ℃、温度の初期値を 20 ℃ にしてみた。

結果の図はつぎのようになる。

• 計算結果を数字の形でファイルに書き出すプログラムと、 それを読みこんでグラフにするプログラムは、 [例題2] のばあいと同様につくることができるが、まだ用意していない。

表計算ワークシートによる計算

[[例題 1] 表計算ワークシートによる計算] のワークシートを部分的に修正すれば実行できるので、共通の点についてはそのページを見てほしい。

「変わらない量の部分」の変更

表題は「シミュレーションの例題3: 0次元エネルギー収支気候モデル、時間変化あり、温室効果の強さは温度の1次関数 (増田 2021-01-16版)」

モデルパラメータとしてつぎのものを追加する。値はシミュレーション実行ごとに変更できる。 参照温度を K で与えたいならば、セルC23をつかわず、セルC24の式のかわりに値を直接与えてもよい。

例題2 で 第20行 に追加した「温室効果の強さ」 ghe は、今回はつかわないので、消してもよい。 念のため、今回は、数値を配置する場所から A20:D20 を除外することにした。

「変わる量の部分」の変更

セル「L6」の式を変更し、M列を追加する。(L6, M6の両方の式の中にある「I6」の「I」はアイです。)

(「温室効果の強さ」のありうる最小値を計算に入れるのは、L列の計算式でおこなっている。 M列は、そのまえの段階の計算であり、値は負になりうる。)

L6、M6の式をそこから下のL列、M列のセルにコピーする。

複数のシミュレーションを同じワークシート内で実行するやりかた

H列からM列までの6つの列をまとめて、あいている列にコピーしてから、 ちがえるべきところをちがえればよい。 (例題2とは列の数がちがうことに注意。)