数値シミュレーションの例として、 [第四紀の氷期サイクルと地球の軌道要素]に関連する Imbrie & Imbrie (1980)のモデルを、パソコンで動かしてみることにする。 このモデルについては、増田(1993)の第5節「10万年周期を出すための数理的モデル」の うしろのほうで紹介している。
気候システムの数理モデルは、たいてい、 システムの状態が、いくつかの変数で表現され、 その変数の時間変化(時間による微分値)が、 システムの現在の状態と、システム外からの条件とによって決まるという構造になっている。
dX / dt = f(X, 外部条件)
ある時刻(「初期時刻」という)でのシステムの状態をあらわす値(「初期値」という)を与え、 時間変化を求め、時間について積分することによって システムの状態変数の時系列が得られる。
微分 dX / dt は、{X(t+Δt) - X(t)}/Δt の Δt → 0 の極限である。 微分方程式が積分によってとけないときは、 極限をとるのをやめて、有限の時間ステップ Δt のあいだの変化の式
{X(t+Δt) - X(t)}/Δt = f( X(t), 外部条件(t) )
を考える。これを変形して、時刻 t での X から 時刻 t+Δt での X の値をもとめる式をつくることができる。
X(t+Δt) = X(t) + Δt × f( X(t), 外部条件(t) )
そしてt+Δt を あらためて t とみなして次のステップに進む。
システム外からの条件が、システム内の状態とは基本的に無関係に変わるとき、 それはシステムにとって、外部からの強制(forcing)となる。 (力学でいう力ではないが、「外力」「強制力」ともいう。増田(1993)では「外力」と書いた。)
Imbrie & Imbrie (1980, 以下 II80と略す)のモデルは、 地球の軌道要素という強制を与えられて、氷期サイクルというシステムのふるまいが生じる 数理的なしくみをデモンストレーションしようとしたものであり、 必ずしも、実際の気候システムで働いている物理的なしくみを表現しているわけではない。 しかし、質量保存やエネルギー保存などの物理的なしくみを表現したモデルは もっと複雑になるので、最初の例題としては、この簡単なモデルをあつかうことにする。
しかし、残念ながら、(この論文が短く書かれたもので説明が少ないせいでもあるが)、 わたしは、この論文の結果をよく再現することが、まだできていない。 この論文の結果の大まかな特徴を再現するところまでで、例題とすることにした。
II80のモデルの時間発展する状態変数は1個で、 世界の氷床の氷の合計量に対応するものである。 論文では y、増田(1993)では I と書かれている。 Rプログラム例文での変数名は ice_mass とした。 【この変数は必ずしも現実の氷の量に比例するものではなく、 むしろ、氷の量の1次関数と思ったほうがよさそうだ。 したがって、計算結果が負になってもあまり気にしないで、増減だけを見ることにする。】
モデルの外から与えられる強制は、 地球の軌道要素変化にともなう(緯度別・季節別の)大気上端にはいる太陽放射量 (「日射量」)が想定されているが、それを1個の変数で代表できると仮定する。 論文では x、増田(1993)では Ieqと書かれている。 Rプログラム例文での変数名は forcing とした。 この変数は、上記の「氷の合計量」と同じ物理量の次元をもつ。 モデルの構造からは、 「仮に日射量がその一定値でずっと続いた場合に達するだろう定常状態の氷の量」という 意味づけがされている。 しかし、現実に大陸氷床がとりうる大きさよりも大きな変化をすることもある。 モデルを単純化したので、いくらか非現実的なところがあるのだ。
y の時間変化は、次のように書かれる。
dy/dt = (1/T) (x - y)
ここで T は時間という次元をもつ定数で、システムの時定数と呼ばれる。 増田 1993 では、τ [タウ]と書いた。 Rプログラム例文での変数名は time_const とした。
Tがほんとうに一定値ならば、このシステムは数学の意味で線形システムで、 入力を正弦波に分ければ、出力は入力と同じ振動数成分をもつ。 (ただし、時間に関する積分がはいるので、振動数ごとの重みは かわってくる。)
大陸氷床には、雪がつもって氷になるのは遅いが、とけるのは速いという性質がある。 そこで、T を、x が y よりも大きいとき(氷がふえるとき)には大きい値 Tc、 x が y よりも小さいとき(氷が減るとき)には小さい値 Tw をとると仮定してみる。 II80での具体的な値は、Tc = 42500 年、Tw = 10600 年。 Rプログラム例文での変数名は time_const_accum、time_const_melt とした。
すると、システムは非線形になり、 出力には入力とちがった振動数成分があらわれることがある。 このようにして、「日射量」の時系列には10万年周期成分はわずかしかふくまれていないが、 氷の量の時系列は10万年周期が卓越している、ということを説明可能である。
入力 x は、緯度別・季節別の「日射量」だが、 緯度別・季節別にどこがどれだけきくかは、簡単には決められない。 しかし、いずれにせよ、次の地球軌道要素で表現できるはずだ。 (ギリシャ文字などの記号はII80論文で使われているもの、 そのあとの文字列は Berger & Loutre 1991のデータで使われているもの。)
ただし、ωは、e sin ω と e cos ω の組み合わせ、 書きかえると、e sin (ω - φ) のような形であらわれる。 (φ[ファイ]は近日点黄経の初期値に関連する定数。)
そこで II80 では、強制項を、
x = ε + α e sin(ω - φ)
のように書き、α = -2、φは 2000年にあたる値を与えている。 (いろいろためして、結果が最近50万年間の海底コアの酸素同位体比の時系列に 似ているものをえらんだようだ。) ただし、この式で、εそのほかは、平均値と標準偏差で規格化してから使うとあり、 εについてはわかるのだが、 そのほかの変数の扱いについて、論文での記述はじゅうぶんでないので、 わたしは論文に出た結果のくわしい再現ができていない。
そこで、今回の例題での強制の与えかたは、次のふたとおりとする。
例題2でこの緯度・季節をえらんだのは、 Milankovitchの1930年の著作以来、北半球高緯度の夏の日射量が重要と考えられているからである。 しかし、計算結果を見ると、10万年周期の氷期サイクルが明確に出てこない。 II80で10万年周期がうまく出たばあいの入力は、これとよくにたものではないようだ。
同じモデルについて、もうすこしくわしく検討した Snieder (1985) の論文がある。 しかし、わたしはまだその内容をよく検討できていない。
計算は、時間による微分を、時間区間Δtの差分で近似しておこなう。 時間区間 Δtは(例題2の入力データにあわせて) 1000年とする。 計算の対象期間は、1000ステップ、つまり百万年間とする。 例題2については、百万年前から現在までの期間をあつかう。
