エントランスへはここをクリック   
地学雑誌 JournalofGeography 116(6)768-7822007

第四紀の氷期サイクルと日射量変動
Insolation Variations and Ice Age Cycles in theQuaternary
TakashiITO*and AyakoABE-OUCHI**+  2007

伊藤孝士* 阿部彩子** *東京天文台、**東京大学理学部
ミランコヴィッチメニューへ戻る

抄録
第四紀の気候は、氷河期の周期によって特徴づけられています。
数万年から10万年の間、長期的な日射量の変動が引き金となって地球の公転運動や年輪運動を計算することができます。長期的には正確に計算できますが氷河期サイクルの本質と起源を理解するためには、地球の複雑な気候システムの中で起こっている物理的・力学的プロセスを知る必要があります。この原稿では、まず、以下のことを簡単におさらいします。また、惑星間の重力相互作用が日射量の変動をどのように引き起こすのかを明らかにしました。次に、氷河期のサイクルがどのようにして地球の気候システムの中で、特に氷の動的モデル化に焦点を当てています。シートを使用しています。
キーワード.第四紀、氷河期サイクル、太陽日射量、摂動理論、気候モデル


Abstract
The climate in the Quaternary is characterized by ice age cycles with periods in the range of tens of thousands to a hundred thousand years, triggered by long-term insolation variations due to the Earth's orbital and precessional motions. Although we can accurately calculate long-term insolation variations at the top of the Earth's atmosphere, we need to know the physical and dynamical processes occurring in the complicated climate system of the Earth in order to understand the true nature and origin of the ice age cycles. In this manuscript, first we briefly review how the gravitational interaction between planets causes insolation variations. Then, we summarize the recent status of large-scale numerical experiments as to how the ice age cycles take place in the climate system of the Earth, with a particular focus on dynamical modeling of ice sheets.

Key words : Quaternary, ice age cycles, solar insolation, perturbation theory, climate model


I. 氷河時代としての第四紀

 微天体の落下と集積により形成した地球の表層は,その初期には極めて高温だったはずである。やがて大陸の形成に伴い大気から二酸化炭素が取り除かれるにつれ,大気による保温効果が弱まって地表付近の温度は次第に下降した。

 現在の地球では炭素循環が気温の変化に対して負のフィードバック効果を持つことで安定な気温状態が保たれている(Berner,1991)。こうした地球の歴史において現在のように高緯度域に巨大な氷の塊(氷床)が存在する時期は氷河時代と呼ばれる。

 氷河時代は地球史の常態であったわけではなく,数億年間隔で断続的に訪れて来た。24-22億年前,約7億年前,約6億年前,約3億年前,そして現在である。約7億年前と約6億年前の氷河時代には氷床が低緯度や赤道付近にも分布し,この時代の地球が全球凍結(スノーボール・アース)に近い状態であったことが判明しつつある(田近,2007)。

 一方で氷河時代と氷河時代の間,例えば約6億年前から約3億年前,また約2億年前から約1億年前にかけては森林が高緯度域に広く分布し,砂漠の広がりも大きく,地球は現在よりずっと温暖だったと考えられている。

 今日の氷河時代は南極大陸に氷床が形成し始めた約3000万年前に始まった。その中で現在から遡ること過去約200万年の時代は第四紀と呼ばれているが,第四紀における気候の最大の特徴は,氷床が拡大する時期(氷期)と縮小する時期(間氷期)が周期的に繰り返されることである(阿部・増田,1996)。



 この繰り返しは氷期・問氷期サイクルまたは単に氷期サイクルと呼ばれ,その存在は19世紀にはすでに知られていた。現在の地球は間氷期の状態にあるが,約2万1千年前の地球は最終氷期と呼ばれる典型的な氷期の状態にあったことがCLIMAPやSPECMAPなどの古気候研究プロジェクトによって明らかにされている(CLIMAPProject,1976,1981;Bergeretal.,1984)。

 海洋底から掘削された堆積物や極域から採取された氷床コアのデータにより,今日では以下のような最終氷期の世界像が得られている。

 最終氷期には北米大陸のローレンタイド氷床とユーラシア大陸のフェノスカンジア氷床が大きく発達しており,現在は緑豊かになっていても当時は氷床に覆われていた寒冷地が広くあった(図1)。全球の平均海面水温は今より3℃以上低く,大気中の二酸化炭素濃度は現在の1/2程度であった。

 中緯度地域には乾燥域が広がり,砂塵を多く含む風が吹き荒れていた。北大西洋の深層循環は現在と比べると不活発であったりと,海洋循環の形態も現在とはだいぶ異なっていた。

 氷期と問氷期が交互に繰り返すことの原因については19世紀後半から20世紀前半にかけてJ.Adhemar,J.Croll,L.Pilgrimらにより議論され,地球軌道の変化と自転軸の歳差による極域への日射量の変動が関与しているという説が形成されて行った(中島,1980)。

 白い氷床が太陽光を反射することで気温が低下し氷床が発達しやすくなる効果(氷床・アルベドフィードバック)も早くから認知され,日射量変動と氷期サイクルが相関する可能性は高いと考えられていた。

II. ミランコビッチ・サイクル

 こうした一連の定性的議論を詳細な計算により初めて定量化したのがセルビアの科学者M.Milankovitchである(木下,1993)。Milankovitchは地球へ入射する日射量の緯度分布と季節変化を当時としては驚くべき高精度で計算し,また,彼の時代の最新の公転軌道変化理論を用いて日射量の長期変動を詳しく計算した(Milankovitch,1920,1930,1941)。

 Milankovitchの研究の大半は1920年代から1930年代にかけて遂行され,1950年代あたりまでは地質データを解読する際の時間目盛としてその結果を使おうとする研究者も存在した。

 当時は放射性同位元素による精密な年代測定法がまだ確立しておらず,氷期に関連する地質データを解釈する上での年代尺度がMilankovitchが計算した日射量の理論値くらいしかなかったからである。当時の年代推定は,地質データから推測される気候条件の変化をMilankovitchの日射量変動理論値と定性的に比較し,データの年代を推算するというものである。

 容易に想像されるとおり,この方法で推定される年代の不確定性は大きい。また1960年代に入り放射性同位元素による年代測定法がその精密さを増して行くと,年代尺度として日射量の理論値を使う必要性が薄れてきた。

 その上,Milankovitchが計算した日射量変動と地質データの示す気候条件の間にある微妙な差異が問題にされたこともあり,氷期サイクルの要因を地球への日射量変動に基づいて議論する風潮はいつの間にか消えてしまった(増田,1993)。

 時は流れ1970年代に入ると海洋底や氷床の掘削試料が豊富に採取され,連続性の良い気候変動指標が手に入るようになった。

 気候変動指標とは気候の状況を示す地質学的な一次データ,例えば海洋底から取得した試料中の酸素同位体比や浮遊性生物化石の種の構成などである(阿部・増田,1996)。近年盛んに研究されているものとしては氷床コア中の気泡に含まれる二酸化炭素濃度や酸素・窒素比なども挙げられる(Kawamuraetal.,2007)。

 同じ頃,時系列データを周波数領域へ変換して議論するスペクトル解析の手法も発達した。すると,地質試料に記録された気候変動指標の変化がMilankovitchの理論的な日射量変化と酷似した周期性を持つことが報告され始めた(Haysetal.,1976)。CLIMAPやSPECMAPなどのプロジェクトにより得られた海底堆積物の数多くから,Milankovitchの理論と調和的な周期を持つ気候変動指標が得られた。

 例えば図2aはSPECMAPから得られた酸素同位体比異常δ180の過去約80万年の時系列だが、,この時系列を周波数分解してみると,後述する日射量変動の特徴的周期である1万9千年,2万3千年,4万1千年,5万4千年といった周期成分が卓越していることがわかりる(図2b)。



 こうした結果は,第四紀の氷期サイクルが日射量変動に駆動されているというミランコヴィッチ(Milankovitch)説の正しさを支持するものである。

 地球の気候システムを理論的にモデル化する研究も同時期に進展し,わずか数%の日射量変化でも氷期サイクルのような激しい気候変動をもたらし得るという認識も次第に広まった。こうして第四紀の氷期サイクルの起源を日射量の変化に求める考え方は復活して一般に定着し、このような周期的な日射量の変化はミランコビッチ・サイクル(MilankovitchCycles)と呼ばれるようになった。

 最近では,第四紀のものに限らず堆積物に記録された周期的な変化を解釈する際にミランコビッチ・サイクルの反映である可能性を考えることはごく自然に行われているように見える(安成・柏谷,1992;増田,1993)。

 もちろん,現時点でもミランコビッチ・サイクルと氷期サイクルの関係について未詳な点は幾つも残されている。図2aに示された氷床変動の時系列が含有する最強の周期成分は約10万年である(図2b)。だが前述した1万9千年,2万3千年,4万1千年といった周期成分が日射量変動の時系列にも共通して存在するのに対し,この10770(a)(b)万年周期は日射量変動にはほとんど含まれていない。

 また,図2aに見られるように氷床の消長は対称的なものではなく,氷床の拡大に要する時間に比べて縮小に要する時間はずっと短い。このような非対称の氷床変動は鋸歯型振動とも称され,10万年周期は日射量変動にはほとんど含まれていない。

 また,図2aに見られるように氷床の消長は対称的なものではなく,氷床の拡大に要する時間に比べて縮小に要する時間はずっと短い。このような非対称の氷床変動は鋸歯型振動とも称され,10万年周期の問題とともに第四紀の氷期サイクル研究における大きな謎として残されている。

 こうした問題を解決するために理論的な気候モデルを用いた数値実験が盛んに行われており,それについては本稿の第V,VI章で述べる。ミランコビッチ・サイクルの原因である地球の軌道要素と自転軸の変化は,惑星間に働く重力に支配された力学過程である。

 次章ではその概略を眺め,惑星間の重力相互作用が地球の日射量変動に至るまでの道筋を復習するが,その前にミランコビッチ・サイクルという用語に関して時折見られる混乱について言及しておく。

 それは「ミランコビッチ・サイクルとは第四紀の氷期サイクルそのものである」という言明である。Milankovitch自身はその著作の中で氷期サイクルの概念にも触れているので,この言明が一概に誤りと言うわけではない。

 しかしMilankovitchが提起した論点のうち後世の研究に最も大きな影響を及ぼしたのは,日射量という外力の変化と気候システムの応答を明確に区別して議論する考え方である。その意味ではあくまで単純に「ミランコビッチ・サイクルとは惑星の公転運動と自転運動に起因する日射量変動の代名詞である」と捉えたほうが,氷期サイクルの研究においては建設的な議論を行うことができると考える。

 本稿でもこのような意味でミランコビッチ・サイクルという用語を使うこととする。

III. 軌道要素の変化を計算する

 地球を始めとする惑星に入射する日射量の時間変化を計算するためには,質点としての公転運動と形状を持つ天体としての自転軸の運動が計算できれば良い。広く知られているように,二天体がお互いの重力のみに従い運動を行う場合(重力二体問題)には厳密な解析解(ケプラー運動)が存在する。

 だが天体が三個以上になるとその運動を表す解析的な解は存在せず,物事は急に複雑になる。惑星の軌道は楕円であるとしばしば表現されるが,これは「惑星軌道は楕円に近い形状をしている」という意味に過ぎない。実際の惑星軌道は互いの重力の影響によって常に揺れ動いており,一定の楕円形状を保っているわけではない。

 しかし太陽系惑星の運動においては太陽重力の影響が他の力に比べて圧倒的に大きい。太陽重力の大きさを1とすると惑星間の重力相互作用の大きさは最大でも1/1000程度である。

 これの意味するところは,太陽系惑星の運動はほぼ太陽重力に支配されており,それに対して惑星間の重力相互作用が微小な摂動として働いているということである。私達は重力二体問題に関してケプラー運動という厳密な解を知っている。

 従って,重力二体問題からさほど外れてはいない惑星運動の第一近似解としてケプラー運動を採用し,それに対して惑星間の重力相互作用による補正を少しずつ積み重ねて行けば(逐次近似),最終的にはかなり正確な惑星運動の解を得ることができる。このような方法で「解けない問題」の近似解を得る手法は摂動論と呼ばれ,惑星運動に限らず今や物理学全般で広く使われている概念である(木下,1998)。

 ところで,惑星の軌道運動はその時間スケールによって二種類に分けられる。まず(1)その年その日の楕円軌道の上で惑星が時々刻々と位置を変える運動がある。これはいわゆる狭い意味での公転運動であり,強大な太陽重力に起因する。この時間スケールは地球であれば一年である。

 次は(2)惑星が日々それに沿って公転する軌道自体の形状や方向が変化する運動で,惑星間の重力相互作用が原因で発生する。氷期サイクルを引き起こす日射量変動をもたらすものはこちらの運動であり,時間スケールは一般に数万年から数十万年と長い。

 一公転周期程度の時間では地球の軌道形状や方向はほぼ変わらず,例えば近日点が現在の場所から大きく移動するには何千年もの時間が必要になるということである。国立天文台が編纂する『理科年表』に記されている地球の近日点方向や離心率の値が昨年も今年も来年もほとんど同じであることがその事実の一端を示している。

 この性質に着目し,惑星の運動方程式を公転周期について平均化するというテクニックが存在する。要するに(1)の効果を無視するのだが,この方法を用いると惑星の運動に関与する項の数が激減し,方程式はかなり簡単な形になる。

 そうして得られた方程式を用いて摂動論的に近似解を求めれば,計算量がとても少なくて済む。この方法もまた古くから知られており,永年摂動論と呼ばれている。最先端の永年摂動論の精度は非常に高く,惑星の運動方程式を計算機で直接解いた結果(数値積分)とも良い一致を示している(Laskar,1988)。

 ちなみにこうした計算結果の不定性がどのくらいであり,得られた解がどの程度の期間にわたり有効であるのかを問われることが頻繁にある。これは実はとても難しい質問と言える。数千万年や数億年という長期の惑星運動に関しては計算結果と比較できる観測事実も存在しないので,計算結果の有効期間を具体的に示すのが難しいのである。

 しかし惑星の運動方程式の解はとても滑らかなので長期間にわたる高精度求解が可能であるし,異なる研究者が異なる方法を用いて求めた解が数億年にわたり良い一致を示すことも知られている(ItoandTanikawa,2002)。従って第四紀の時間スケールでの日射量変動を知るという目的では,こうした計算から得られる解が十分に信頼できると仮定しても良いだろう。

 さて,地球に入射する日射量変動を計算するには公転軌道の変化を知るだけでは十分でない。地球の自転軸は公転軌道面に対して23.4度傾いているが,この値は数万年周期で変動する。また自転軸は慣性系に対して静止しているのではなく,公転軌道面の法線のまわりを約2万6千年で周回する。これが自転軸の歳差運動と呼ばれるものである。

 地球に入射する日射量の変動はこうした自転軸の運動にも大きな影響を受けるので,公転軌道要素の変化に加えて自転軸の歳差運動を詳しく解いて初めて日射量変動の計算が可能になる。一般に惑星自転軸の運動周期はその自転周期や公転周期に比べてかなり長い。

 地球の歳差運動の周期2万6千年に対して公転周期は一年,自転周期は一日,という具合である。また,惑星の形状(正確に言えば重力ポテンシャル)は完全に球対称ではないものの,赤道部分が若干膨らんだ軸対称の回転楕円体としての扱いが良い近似になることが多い。

 従って,軸対称な形状を持つ惑星の自転運動に関する運動方程式を立て,それを公転周期に関して平均して簡略化すれば,公転軌道要素の永年摂動を計算した時と同様にして惑星自転軸の長年変化を計算することが可能になる。

 以上のような手続きを踏み,摂動論や数値積分と言った解法を駆使して惑星の公転運動と自転運動の変化を追うことで,私達は過去や未来の特定の時点での地球軌道の形状と方向,自転軸の傾きと向き,太陽との位置関係を知ることが可能になる。

 それさえわかればあとは幾何学的な関係によって地球上の任意の地点での日射量を計算できる。1930年代にMilankovitchが行ったのはこのような定量的計算であった。

IV. 地球の日射量変動

 惑星の公転運動と自転軸の歳差運動を正確に記述するにはその自由度に応じた多くの変数が必要となる。しかし長期の日射量変動を議論するだけならば関連する変数は数えるほどしかない。

 特に重要なものは軌道の離心率,近日点の方向を表す角度(近日点経度),自転軸の傾きを表す角度(赤道傾角),それに自転軸の方向を表す角度(歳差角)である。これらの量がどのようにして長期の日射量変動に関与するのかを詳細に説明するには本来は数式が不可欠なのだが,本稿ではそれを行わず,簡略な定性的説明を試みる。氷期サイクルを駆動し得る日射量変動には主として以下の二種類がある。

 第一は赤道傾角の変化による日射量変動で,赤道傾角項と呼ばれる(図3a)。順」という呼称が使われるのは,具体的なロ射量変動の数式の各項がこのように呼ばれるからである(中島,1980)。本稿では「項」=「効果」くらいに考えて差し支えない。二番目は歳差角と近日点経度および離心率の変化による日射量変動で,気候的歳差項と呼ばれる(図3b)。



 赤道傾角項の働きは直感的にわかりやすい。地球の赤道傾角は主として公転軌道面の運動に依存して変化するが,赤道傾角が大きい場合,すなわち地球の自転軸が大きく傾いている場合には高緯度の夏の日射量が相対的に大きくなり,低緯度の夏の日射量は小さくなる。

 逆に赤道傾角が小さい場合には高緯度の夏の日射量は減り,低緯度では日射量が増える。つまり赤道傾角項は緯度帯毎の日射量のコントラストを変える効果を持つのである。これに対して気候的歳差項は季節のコントラストに影響を与える。地球が最も太陽に近くなる地点(近ロ点)での自転軸の方向は歳差運動によって刻々と移り変わる。もしもある年,北半球の夏の時点で地球が近日点にあれば,北半球での夏の日射量は非常に大きくなるであろう。

 一方,(a)Obliquity(b)Climaticprecessionその年の冬は地球は遠日点にあり,北半球が受ける日射は少なくなる。逆の場合も同様であり,北半球の冬の時点で地球が近日点にあれば北半球での夏の日射量は比較的小さくなるし,冬は近日点にあたるために日射量は比較的大きくなる。地球の軌道が楕円であることと自転軸が歳差運動することにより,季節毎の日射量のコントラストが長い時間スケールで変動するのである。

 こうした事柄を頭に入れて,地球の離心率,赤道傾角,気候的歳差の時系列とその周波数成分の計算結果を見てみよう(図4)。



 地球の軌道は他の惑星軌道に比べて円に近いと言われるが,その離心率は図4aのように0から0.06程度の値をとりつつ約10万年と約40万年の特徴的周期を持って振動する(図4e)。赤道傾角の現在値は約23.4度だが,惑星摂動による地球の軌道傾斜角変などによりこの角度は約4万1千年と約5万4千年の周期を持って振動を行う(図4b,f)。

 赤道傾角の振幅は小さく,±1度程度である。気候的歳差は離心率×sin[近日点経度]という関数形を持ち,1万9千年付近と2万3千年付近に強い周期性を持つ(図4c,g)。以上のような力学変数の時系列から計算される日射量の長期変動が図4dであり,それを周波数分解したものが図4hである。

 ここではBerger(1978)の定式化を用い,北半球の夏至における北緯65度地点での日平均日射量を計算した。氷期サイクルの要因という観点では氷床分布の代表的緯度として北緯65度付近での夏の日射量変動を議論することが多いので,本稿でもそれに従っている。

 なお夏至一日の代わりに夏季三ヶ月や夏半年の平均日射量を計算すれば,その変動振幅は図4dのものより幾分か小さくなろう。図4e,f,g,hに示した周波数分解の結果を見れば明らかなように,日射量変動はもっぱら気候的歳差と赤道傾角に依存し,1万9千年,2万3千年,4万1千年といった気候的歳差と赤道傾角の典型的周期が日射量変動の周期性(図4h)にそのまま反映されている。

 これらは第四紀の氷期サイクルが持つ周期として図2bのような気候変動指標のデータにも見られるものでもある。一方で,氷期サイクルが持つ最も強い周期性である約10万年の変動成分は図4hに示された理論的な日射量変動にはほとんど見られない。

  前述したように,この事実は氷期サイクル研究における最大の謎として数多くの研究者を悩ませて来た問題である。理論的な日射量変動に10万年の周期成分が全くないわけではなく,約10万年周期で振動する離心率(記号eで表される)に対して(1-e2)-1/2の依存性を持つ項は存在する。

 だが地球の離心率の値は0.01の桁に留まるので,この項の寄与は小さい。図4hにある日射量のスペクトルのうちこの項に起因する10万年周期成分(~97kyrおよび~125kyrと記されたもの)の強度は赤道傾角項の1/10程度,気候的歳差項の1/20程度に過ぎない。

 離心率は気候的歳差項においてもesin[近日点経度]という形で間接的に日射量変動に関与するが,離心率eのみならずsin[近日点経度]も周期的な関数であるため,その積に離心率自身の持つ周期成分は現れない(伊藤,2002)。

 こうした事情は古くから認識されていたものの,氷期サイクルの卓越周期である10万年に近い周期で振動する力学変数が地球軌道の離心率くらいしか見当たらないため,「第四紀の氷期サイクルの直接の原因は離心率の時間変化に起因する日射量変動の10万年周期成分である」という半ば強引な説が唱えられることもあった。

 だがもし(1-e2)-1/2の項に比例する日射量変化が氷期サイクルの10万年周期変動に直接反映されているならば,気候的歳差項や赤道傾角項に対するその相対的な影響度は何らかの原因により数十倍に増幅されていることになる(図2bと図4hにおける10万年周期成分の相対的な強度を比較してみれば良い)。

 このため,日射量変動に見られる微弱な10万年周期変動からいかにして顕著な10万年周期の氷期サイクルが引き起こされるかについての理論的な仮説が山のように提唱されて来た。

 けれども,複雑を極める地球の気候システムの全容がいまだ明らかになっていない現時点では,10万年周期という数値の類似性が明確な物理的根拠を持つものなのか,または単なる偶然の一致なのかは不明と言うべきである。この謎を解くために詳細な観測データの取得と解析が有効であることは当然だが,それとともに重要な手法は第V,VI章で述べられる数値的な気候のモデル化であろう。

 ちなみに火星など他の惑星についても図4と同様な計算を実施することは可能である。そうした結果を見ると,火星らと比べて地球の赤道傾角の変化はとても小さく,結果として長期の日射量変動の振幅が抑制されていることがわかる(伊藤,2004)。自転軸の歳差運動は重力的に扁平な惑星に対して他の天体からトルクが働くことにより発生する(ここで言う「トルク」とは惑星の自転状態を変化させる重力のモーメントである)。

 地球以外の惑星については太陽からのトルクが強く働くが,地球は例外で,月からの重力トルクが775太陽からの二倍も大きい(伊藤,1993;Itoetal.,1995)。つまり地球自転軸の歳差運動は月からの強い重力トルクにより安定化されており,それにより赤道傾角の変動振幅が抑制されているのである(Laskaretal.,1993)。

 月の存在による歳差運動の安定化は地球史の初期から存在していたメカニズムであり,その帰結のひとつと考えられる気候の安定化を高等生命の発生や進化と結び付ける議論も少なくない(井田,2003)。

 日射量変動をつかさどる惑星の公転運動と自転運動の研究は長い歴史を持ち,今や数理科学の中でも最も精緻な理論体系を持つ分野のひとつとなっている。解の精度を上げるための研究は現在も続いているが,第四紀の氷期サイクルを議論するのに必要な精度という意味ではこの分野はすでに完成の域に到達していると言って良い。

 例えば,現在最も広く使われているBerger(1978)の日射量変動の定式化は1940年代から1960年代に構築された方法(BrouwerandvanWoerkom,1950;SharafandBudnikova,1969a,b)を踏襲したものだが,これを使った場合でも,また21世紀の最新の方法(Laskaretal.,2004)を使った場合でも,結果的に計算される日射量変動の時系列は見分けが付かないほど良く似ており,従って氷期サイクルに対する影響もほぼ同等となることが様々な気候モデル計算から示唆されている。

 日射量変動の理論体系構築という観点での関心は,地球と同様な氷期サイクルが存在するのではないかと思われる火星上での日射量の計算や,太陽系外惑星系に存在するだろうまだ見ぬ「第二の地球」での気候変動の予測という方向に向かいつつある(Atobeetal.,2004;Abeetal.,2005)•B

V. 数値的な気候モデル

 前章で述べたように,地球に到達する日射量の分廊や時間変化については今や精密な計算が可能である。

 だが,私達が精密に計算できるのは地球大気の上端に達する太陽からのエネルギー輻射量に過ぎない。大気上端に届いた日射量の変動は大気や海洋を要素とする気候システムを介した後にようやく氷期サイクルとして発現し,その過程は決して単純なものではない。

 日射量変化の時系列には含まれない10万年の周期性が氷床変動の時系列に卓越する事実などは,地球の気候システムの奇怪とも言える複雑性を示す典型的な例であろう。計算機技術が発達した現在では,このような問題を解明するために数値的な気候モデルを用いた研究が盛んに行われている。

 以下の章ではこうした研究に用いられる気候モデルの現状を紹介し,著者ら自身が行っている数値実験の成果に関する簡単な解説を試みたい。第四紀の氷期サイクルのように長い時間スケールの気候変動を解明する上で注目すべき量的変数は,気温,大気の乾燥湿潤状態,氷床量,海洋循環の状態などである。気候モデルを構築する上ではこうした諸量がどのような物理過程に支配されているのかを知っておくことが重要である。

 全球平均的な気温を決定するのは気候系全体としてのエネルギー交換過程,つまり地球に入射する太陽エネルギーと地球から出て行くエネルギーの釣り合いである。局地的な気温や大気の乾燥湿潤状態を決めるものは,大気中および陸面や海面におけるエネルギー・水・運動量の交換過程である。

 エネルギー・水・運動量の交換に関与する気候のサブシステムには大気(輻射,全球規模の循環(大循環),陸地との水循環の過程などがかかわる),海洋(大循環,対流,拡散過程などが関与する),氷床や海氷(相変化や流動過程が関与する)があり,これに加えて陸面自体や植生が水循環や熱交換過程に関与する場合がある。

 海洋表層部分の循環は大気との運動量交換により支配される。海洋深層部分の循環には海水密度の決定要因である水温や塩分が重要であり,海洋内でのエネルギー交換過程と淡水の存在に影響を受けるし,大気との相互作用も関与する。氷床の形状を決定する要因としては降雪に代表される氷床の酒養過程,大気とのエネルギー交換過程,それらに影響を受ける融解などの消耗過程,そして氷床自身の流動が挙げられる。

 氷期サイクルを再現するための数値モデルにはこうした主要な気候構成要素や物理過程が過不足なく取り込まれている必要がある。このような数値モデルの代表的なものには,エネルギーバランスモデル,大気海洋結合大循環モデル,氷床力学モデルがある。

 地球の最も基本的な気候状態は,太陽から受け取る輻射エネルギーと地球から出て行くエネルギーの収支で決まる。この概念に基づき,地表でのエネルギー収支のみを主要な構成要素とした気候モデルをエネルギーバランスモデルと呼ぶ。エネルギーバランスモデルは後述する大循環モデルや氷床力学モデルに比べて単純であり,運動方程式や物質の出入りを表す関係式を含まない。

 こうした簡単なモデルを用いると,エネルギー収支を表現する方程式系の解構造推定や安定性解析といった定性的な見通しを得る作業が容易になるので,気候システム全体を傭轍しての議論を行いやすくなる。一方,エネルギーバランスモデルから得られる解の定量性は拡散係数やアルベドなどモデル内のパラメータとその数学的表現に大きく左右される。

 また,このモデルでは大気・海洋循環や降水量に深く関与する水循環を表現することが困難である。このような困難はエネルギーバランスモデルがその構造上必ず内包するものであり,大循環モデルや氷床力学モデルなどを用いない限り解決は不可能である。

 大気海洋結合大循環モデル(以下では単に大循環モデルと呼ぶ)は,流体の運動方程式を数値的に解くことで大気や海洋の運動を計算機内に再現しようとするモデルである。大循環モデルでは大気や海洋を三次元の格子に細かく区切り,格子問での運動量や熱・水蒸気・塩分のやり取りを計算する。

 大循環モデルを用いて様々な数値実験を試行すれば,どの物理過程がどのような気候状態を支配しており,日射量に代表される外的入力が気温や降水量などの出力にどう影響するのかを事細かに検証することが可能となる。

大循環モデルは主として二種類のサブモデルから構成される。ひとつは大気大循環モデルであり,気温,水蒸気量,風,地上気圧,降水量,積雪量,土壌水分量などを予報変数(解を求めるべき物理量)とする。もうひとつは海洋大循環モデルであり,海流,海水温,塩分などを予報変数とする。この他にも海氷の分布や厚さの変化を計算する海氷モデルや,海洋への陸水供給過程を再現する河川モデルなどが大循環モデル内に含まれることもある。

 これらのサブモデルを計算機の中で結合し,地球表層において物質やエネルギーを循環させてその振舞いを仔細に観察する仕組みこそが大気海洋結合大循環モデルである。大気大循環モデルの基本は天気予報に用いられる数値モデルであり,数値天気予報とともに発展してきたと言って良い。

 だが数値天気予報モデルが数日という短時間スケールの現象を扱う目的で構築されているのに対し,気候研究に使われる大気大循環モデルは初期に与えた気候状態の情報が失われてしまう長時間スケールでの現象を相手にする。言葉を変えれば,日射量や海面水温といった境界条件に対する系の応答が定常状態に移行した時点での平均的な場を知るための道具が大気大循環モデルであると言える。

 大気大循環モデルに含まれる方程式には以下のようなものがある:静水圧近似式,水平方向の運動方程式,連続の式,エネルギーの保存式,水蒸気の保存式,状態方程式,など。モデル内でその時間発展が計算されるのは大気の流れ(いわゆる風),気温の変化,摩擦,加熱,水蒸気発生に関する過程などである。

 ここには雲の発生や水蒸気の凝結,太陽エネルギー輻射などの過程も含まれる。モデル内で設定される格子よりも小さな空間スケールを持つ物理過程(微過程)は直接には計算できないので,経験的な式で表現される。この手法はパラメタリゼーションと呼ばれるが,パラメタリゼーションの違いにより大気大循環モデルが持つ性質は大きく異なるものとなる。

 大気大循環モデルの改良と言えば現在ではもっぱらこうした微過程のパラメタリゼーションを指し,各種の観測結果を参照してより良い経験式を構築する努力が不断に行われている。大循環モデルのもう一方の構成員である海洋大循環モデルは,文字通り海洋の大循環を表現する目的で開発されたものである。

 海洋大循環モデルを構成する方程式系は大気大循環モデルと類似しており,連動方程式,連続の式,熱力学の式,状態方程式,静水圧近似式などである。海洋大循環777モデルと大気大循環モデルの違いは以下のようになろう:状態方程式の形の違い,水蒸気の保存式の代わりに塩分の保存式を用いること,境界条件の設定が海底や海岸で行われること,など。海洋大循環モデルにおいてもパラメタリゼーションは用いられる。

 例えば格子サイズより小さくて直接には計算不能な渦の効果を海水の拡散過程に組み込んで扱う手法などである。大気大循環モデルに比べると海洋大循環モデルはパラメタリゼーションへの依存部分が一般に少なく,構造が幾分か単純であるとは言える。

  ただしパラメタリゼーションの実装の違いが最終的な計算結果に無視できない影響を与えることはあり得るので,計算の際には留意する必要がある。容易に想像されるとおり,大循環モデルを動かして意味のある結果を得るためには膨大な量の計算機資源が必要である。

 大気大循環モデルを稼働させてその定常的な解を得るためには約50年間分の数値積分が必要だが,間隔約200kmという粗い空間格子を設定し,現代最高速度のベクトル計算機を用いてもこの計算には数日を要する。また,気候システムの中で最長の時間スケールを持つプロセスは海洋の深層循環である。

 大気海洋結合大循環モデルを用いて海洋深層循環の定常状態を得るには約2000年間分の数値積分が必要となるが,これは100日掛かりの計算となる。かくして大循環モデルを用いた研究の進展は計算機の能力向上に強く依存していることがわかる。

VI.氷床力学モデルと氷期サイクル

 氷期サイクルを数値実験で再現するには,氷床が直接関与する力学過程の精密なモデル化が欠かせない。前述したように氷床の形状を決定する主な要因には氷床の湧養過程,大気とのエネルギー交換過程,融解などの消耗過程,氷床自身の流動過程があるが,これら一連の過程を計算機内で再現するためのモデルが氷床力学モデルである。氷床力学モデル内で扱われる予報変数は氷床の面積,高度,温度,流動速度,物質分布などである。

 氷床上端の境界条件には氷床表面での質量収支(降雪による氷床の酒養,融解や蒸発による氷床の消耗)が与えられ,これは主として夏季の気温や降水量に影響される。氷床下端での境界条件には基盤岩から氷床への熱流量および基盤岩の地形が与えられる。氷床の流動は連続体の力学に基づいて計算され,氷床の質量分布の変化に応じた基盤岩の変形も考慮される必要がある。氷床と大気の問には様々な相互作用が働くが,以下の四種類のフィードバック効果が特に重要であるとされている。

 まず高度・質量収支フィードバックと呼ばれる効果がある。氷床が成長すると大気温の低い標高(つまり高地)に存在する氷床の面積が大きくなるので,そこでは氷床が融解しにくくなる。すると更に氷床が成長してその高度を増し,ますます融けにくくなるという効果である。

 次は本稿冒頭にも記した雪氷・アルベドフィードバックである。これは白い氷床の存在により局地的にアルベドが増大すると太陽光の吸収率が下がり,周囲の気温が低下して氷床が一層成長し易くなるという効果である。三種目は定在波・温度フィードバックと呼ばれるものであり,氷床という厚さ数千メートルもの冷源の存在により大気の惑星波(定在波)が変調され,大気大循環を介して気温の空間分布が全球的に変わり,氷床の成長を助長し得るという効果である。

 四番目はストーム・トラックフィードバックと呼ばれるもので,氷床の存在により低気圧経路と降水帯が氷床の南縁に局在化し,それが原因して氷床自身の滴養が助けられるという効果である(阿部,2002)。氷床変動に関するこうしたフィードバック効果についてはかつて幾つかのモデル研究がなされて来た。

 だが現実の大気・海洋状態を最も良く模擬すると思われる大循環モデルを氷床力学モデルと組み合わせた研究がほとんどなかったため,どのフィードバック効果が氷床の消長にどのくらい重要な役割を果たしているかは未解明であった。理想的には氷床力学モデルを大循環モデルのサブルーチンのひとつとして直接組み込み,氷床の変動を大気海洋大循環の一過程として逐次追い掛けるのが最善である。

 しかしこれには膨大な計算機資源が必要とされ,現時点で実現可能なやり方と778は言えない。そのような状況の下,筆者らはかつてなく詳細な氷床力学モデルの開発を行い(AbeOuchiandBlatter,1993;Abe-Ouchietal.,1994;SaitoandAbe-Ouchi,2004,2005),それを大循環モデルの計算結果と組み合わせる数値実験を実施して,氷床と大気の相互作用に関する検証を続けている(Abe-Ouchietal.,2007)•B筆者らの方法ではまず,幾つかの典型的な氷床分布を仮定して大循環モデルを用いた計算を行う。

 ここでは様々なパターンで存在する氷床が大気・海洋系にどのような影響を及ぼし,その結果として氷床自身にどのような影響が戻ってくるのかをあらかじめ知ることが目的である(例えばAbe-Ouchietal.(2007)では大循環モデルを用いて18種類の計算を行っている)。

 この計算の結果を踏まえ,高度・質量収支フィードバック効果や雪氷・アルベドフィードバック効果による気温や降水量の変化が氷床の高さや面積にどのように依存するのかを表す経験的な式を構築する。

 また,時刻の経過に伴う日射量および大気中の二酸化炭素濃度の時間変化は大循環モデルに対する境界条件(外力)として扱う。大循環モデルを用いた計算結果を基にして外力の変化による気温や降水量の変化パターンを求めておけば,氷床の形態に応じた気温や降水量の経験的変化式を構築した際と同様にして,様々な外力が働いた場合の気温や降水量変化の経験式を作り上げることができる。

 こうした経験式を用い,時間変化する境界条件の下で氷床力学モデルを動かせば,現実的な地球の気候環境下で氷床が変動する様子を一連の時間発展過程としてつぶさに観察できるのである。こうした計算結果のひとつが図5にある北半球氷床の分布である。この結果は観測データから推定された最終氷期の氷床分布(図1)と良く一致している。



 ローレンタイド氷床に比べてフェノスカンジア氷床が小さめなこや,ラブラドル半島からロッキー山脈に至る地域がすべて氷床で繋がれるといった最終氷期の描像が,一連のモデル計算により再現されたのである。またこうした計算の過程で,氷床の成長に関する四種類のフィードバック効果の中でも高度・質量収支フィードバックと雪氷・アルベドブイードバックが他に比べて特に重要であることも判明した。氷床周辺の気温変化に対する寄与と言う点では,高度・質量収支フィードバックと雪氷・アルベドフィードバックからの影響が他の二種類のフィードバックに比べて10倍以上になることもある。

 なお観測データからの推定(図1)と比べると,モデル計算の結果(図5)ではフェノスカンジア氷床がかなり南に張り出している。これは氷床の成長に対する大気大循環の影響がまだ完全にはモデルに取り込まれておらず,そのためにユーラシァ大陸での降水量が現実よりも過大になったことが原因と考えられる。

 今後のモデル改良点のひとつである。また計算された氷床の高度分布は観測からの推定値に比べてやや低めに見えるが,最終氷期の氷床高度推定については観測データの不定性がまだ大きいので,現時点では氷床高度の比較にさほど大きな意味はないと言える。

 さて第四紀の氷期サイクルに関するモデル計算において氷床の空間分布の再現と同様に重要なのが,図2aにあるような10万年周期の鋸歯型氷床変動の発生要因を解明することである。

 筆者らの数値実験では,10万年という周期の値も含めてこの鋸歯型振動が良く再現された(図6)。氷床モデルに入力される境界条件パラメータ,つまり氷床周辺の気温減率やアルベド・降水量の温度依存性などにより,氷床変動の程度は異なる。

 けれども計算開始時からじわじわと氷床が発達し始め,約10万年で最大化した後には急激に縮小するという鋸歯型の振動は,その振幅の違いこそあれどのようなパラメータを採用しても定性的にはほぼ確実に再現されたのである(図6c)。

 前記したように,この実験では理論的に計算される日射量と南極の氷床コアデータから推定される大気中の二酸化炭素濃度の変化を外力として与えている(図6a,b)。第IV章で述べたように日射量変動には10万年の周期成分がほとんど存在しないが、観測データを基にした二酸化炭素濃度の変化には10万年周期が含まれている。

 よって厳密に言えば,このモデル計算の結果として得られた氷床の10万年周期変動が外力として入力された二酸化炭素濃度の10万年周期変動を単に反映したものである可能性を否定することはできない。この可能性を否定または肯定するには,大気中の二酸化炭素濃度を一定にして日射量のみを時間変化させた数値実験を実施すれば良い。

 本稿に図を付してはいないが,筆者らはこのための予備的な計算をすでに進めている。その結果,同様なモデルを用いつつ二酸化炭素濃度を一定にして日射量のみを時間変化させた場合でも、振輻は大きくないがやはり10万年周期の鋸歯型氷床振動が観察された。これより筆者らは,氷床コアに記録された二酸化炭素濃度の10万年周期変動は氷床変動の原因ではなく,むしろ結果であろうと予測している。

 この点は今後の数値実験で更に検証されるべき最優先課題のひとつだが,その詳細な記載は別の機会に実現したい。ミランコビッチ・サイクルが第四紀の氷期サイクルを駆動するという説が復活して以来,モデル計算により氷期サイクルの再現を試みた研究の数は果多しい(阿部・増田,1993;Paillard,2001)。

本章で述べた筆者らの計算では非現実的な近似や仮定を用いることなく,ごく自然なモデル化により10万年周期の氷床振動を再現することができた。この結果は世界でも類例を見ないものである

 無論このモデルもまだ完全なものではない。例えば上記した大気中の二酸化炭素濃度に関して言えば,本来これはモデル内で自律的にその変化が計算されるべき量である。また,現在のところ筆者らのモデルでは北緯65度での7月平均日射量しか考慮されていない。他の季節や他の緯度帯に注れた二酸化炭素濃度の10万年周期変動を単に反映したものである可能性を否定することはできない。

 この可能性を否定または肯定するには,大気中の二酸化炭素濃度を一定にして日射量のみを時間変化させた数値実験を実施すれば良い。本稿に図を付してはいないが,筆者らはこのための予備的な計算をすでに進めている。

 その結果,同様なモデルを用いつつ二酸化炭素濃度を一定にして日射量のみを時間変化させた場合でも,振輻は大きくないがやはり10万年周期の鋸歯型氷床振動が観察された。これより筆者らは,氷床コアに記録された二酸化炭素濃度の10万年周期変動は氷床変動の原因ではなく,むしろ結果であろうと予測している。

 この点は今後の数値実験で更に検証されるべき最優先課題のひとつだが,その詳細な記載は別の機会に実現したい。ミランコビッチ・サイクルが第四紀の氷期サイクルを駆動するという説が復活して以来,モデル計算により氷期サイクルの再現を試みた研究の数は果多しい(阿部・増田,1993;Paillard,2001)。

 本章で述べた筆者らの計算では非現実的な近似や仮定を用いることなく,ごく自然なモデル化により10万年周期の氷床振動を再現することができた。この結果は世界でも類例を見ないものである。

 無論このモデルもまだ完全なものではない。例えば上記した大気中の二酸化炭素濃度に関して言えば,本来これはモデル内で自律的にその変化が計算されるべき量である。また,現在のところ筆者らのモデルでは北緯65度での7月平均日射量しか考慮されていない。

 他の季節や他の緯度帯に注ぐ日射量をも考慮した場合に氷期サイクルがどうなるのかの検証も喫緊の課題である。本稿では数万年周期の日射量変動に対して氷床がどのように応答するかを考えてきた。そう遠くない未来には,大循環モデルと氷床力学モデルを直接結合した大規模な計算も現実のものとなるだろう。

 また第四紀よりも更に長い時間スケールでの気候変動を考える時には,大陸配置や大気成分の違いが気候へ与える影響を詳しく調べる必要も発生しよう。そうした一般的な気候状況を考察する際,または近い将来発見されるであろう太陽系外にある「地球的」惑星の気候について考える際にも,本稿で述べたような地球の第四紀における氷期サイクルの研究は重要な知識の礎を私達に授けてくれるはずである。

 謝辞本稿に記した氷床力学モデルによる数値実験は瀬川朋紀氏および齋藤冬樹氏との共同研究である。瀬川氏からは本稿の図作製に際しても多大な協力を得た。HeinzBlatter氏との議論は有意義かつ刺激的であり,本稿の執筆にあたり数多くの発想を与えてくれた。山川修治氏および一名の匿名査読者からのコメントにより本稿の質は大いに高まった。各氏には改めて感謝の意を表す。


文 献

阿部彩子(2002):気候システムと地球史、.熊澤峰夫 ・伊藤孝士 ・吉田茂 生編:全地球史解 読 東京大 学 出版会,234-258.   概要スライド(阿部彩 子)

阿部彩子 ・増田耕一(1993)、:氷床と気候感度モデルによる研究のレビュー.気象研究 ノー ト,177,183-222.

阿部彩 子 ・増田耕 一(1996):、第 四紀の気候 変動.住明正 ・平朝 彦 ・鳥 海光弘 ・松井孝典編:気 候 変動論,岩波講座地球惑星科学,No.11.岩 波 書 店,103-156。

Abe, Y., Numaguti, A., Komatsu, G. and Kobayashi, Y.(2005): Four climate regimes on a land planet with wet surface : Effects of obliquity change and implications for ancient Mars. Icarus, 178, 27-39.

Abe-Ouchi, A. and Blatter, H. (1993) : On the initiation of the ice sheets. Ann. Glacial., 18, 200-204.

Abe-Ouchi, A., Blatter, H. and Ohmura, A. (1994):How does the Greenland ice sheet geometry remember the ice age? Glob. Planet. Chang., 9, 133-

Abe-Ouchi, A., Segawa, T. and Saito, F. (2007) : Climatic conditions for modelling the northern hemisphere ice sheets throughout the ice age cycle.Clim. Past, 3, 423-438.

Atobe, K., Ida, S. and Ito, T. (2004) : Obliquity varia.-tions of terrestrial planets in habitable zones. Icarus, 168, 223-236.

Berger, A.L. (1978): Long-term variations of daily insolation and Quaternary climatic changes. J Atmos.Sci., 35, 2362-2367.

Berger, A.L., Imbrie, J., Hays, J., Kukla, G. and Saltzman, B. eds. (1984): Milankovitch and Climate. D. Reidel, Norwell, Mass.Berner, R.A. (1991) : A model for atmospheric CO2 over Phanerozoic time. Am. J. Sci., 291, 339-376.

Brouwer, D. and van Woerkom, A.J.J. (1950): The secular variations of the orbital elements of the principal planets. Astron. Pap. Am. Ephemeris. Naut.Alm., 13 (2) , 81-107.

CLIMAP Project (1976) : The surface of the ice-age Earth. Science, 191, 1131-1137.

CLIMAP Project (1981) : Seasonal Reconstruction of the Earth's Surface at the Last Glacial Maximum.Geological Society of America.

Hays, J., Imbrie, J. and Shackleton, N. (1976) : Variations in the earth's orbit : Pacemaker of the ice ages. Science, 194, 1121-1132.

井 田 茂(2003):異 形 の 惑 星(NHKブ ッ ク ス966).日本 放 送 出版 協 会 、

伊 藤 孝 士(1993):ミ ラ ン コ ビ ッチ サ イ ク ル の 進 化.地球,15,316-323.

伊 藤 孝 士(2002):日 射 量 変 動 の 基 礎 理 論 熊 澤 峰 夫 ・伊 藤 孝 士 ・吉 田 茂 生 編:全 地 球 史 解 読.東 京 大 学 出版 会,137-161.

伊 藤 孝 士(2004):火 星 の 日射 量 変 動 と気 候.遊 ・星 ・人(日 本 惑 星 科 学 会 誌),13,137-144.

Ito, T. and Tanikawa, K. (2002) : Long-term integrations and stability of planetary orbits in our solar system. Mon. Not. R. Astron. Soc., 336, 483-500.

Ito, T., Masuda, K., Hamano, Y. and Matsui, T. (1995) :

Climate friction: A possible cause for secular drift of the earth's obliquity. J. Geophys. Res. , 100,15147-15161.

Kawamura, K., Parrenin, F., Lisiecki, L., Uemura, R.,Vimeux, F., Severinghaus, J.P., Hutterli, M.A., Nakazawa, T., Aoki, S., Jouzel, J., Raymo, M.E., Matsumoto, K., Nakata, H., Motoyama, H., Fujita, S., Goto-Azuma, K., Fujii, Y. and Watanabe, 0. (2007) :

Northern Hemisphere forcing of climatic cycles in Antarctica over the past 360,000 years. Nature,
448, 912-916.

木 下 宙(1993):ミ ラ ン コ ヴ ィ チ周 期 計 算 の 基 礎 に つい て.地 球,15,314-315.

木 下 宙(1998):天 体 と軌 道 の 力 学.東 京 大 学 出 版 会.

Laskar, J. (1988): Secular evolution of the solar system over 10 million years. Astron. Astrophys., 198, 341-362.

Laskar, J., Joutel, F. and Robutel, P. (1993) : Stabilization of the Earth's obliquity by the Moon. Nature, 361, 615-617.

Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A.C.M. and Levrard, B. (2004): A long-term numerical solution for the insolation quantities of the Earth. Astron. Astrophys., 428, 261-285.

増 田 耕 一(1993):氷期 ・問氷期 サイクルと地球の軌道要素.気象研究ノー ト,177,223-248.

Milankovitch, M. (1920) : Theorie Mathematique des Phenomenes Thermique Produis par la Radiation Solaire. Gautier-Villars, Paris.

Milankovitch, M. (1930) : Mathematische Klimalehre und Astronomische Theorie der Klimaschwankungen. in Handbuch der Klimatologie edited by KOpen and Geiger, Band 1. Teil A., Springer-Verlag.

Milankovitch, M. (1941) : Kanon der Erdbestrahlung und seine Anwendung auf das Eiszeitproblem. Vol. 133 of Koniglich Serbische Academie Publication,Koniglich Serbische Academie. (柏 谷健 二 ・山本 淳之 ・大村 誠 ・福 山 薫 ・安成哲 三 訳(1992):ミ ランコビ ッチ気候変動 の天文学 理論 と氷河時代.古 今書 院)
中島映至(1980):地 球軌道 要素の変動 と気候 気象研究 ノー ト,140,81-114.

Paillard, D. (2001) : Glacial cycles : Toward a new paradigm. Rev. Geophys., 39, 325-346.

Peltier, W.R. (2004) : Global glacial isostasy and the surface of the ice-age Earth : The ICE-5G (VM2) model and GRACE. Ann. Rev. Earth Planet. Sci.,
32, 111-149.

Saito, F. and Abe-Ouchi, A. (2004) : Thermal structure of Dome Fuji and east Dronning Maud Land, Antarctica, simulated by a three-dimensional ice-sheet model. Ann. Glaciol., 39, 433-438.

Saito, F. and Abe-Ouchi, A. (2005): Sensitivity of Greenland ice sheet simulation to the numerical procedure employed for ice sheet dynamics. Ann. Glaciol., 42, 331-336.
Sharaf, S.G. and Budnikova, N.A. (1969a) :

On secular perturbations in the elements of the earth's or bit and their influence on the climates in the geological past. Bull. Inst. Theor. Astron., 11, 231-261.
(原 文 は ロ シ ア語 英 訳 は NASA TTE 12, 467, 1-37)

Sharaf, S.G. and Budnikova, N.A. (1969b): Secular perturbations in the elements of the Earth's orbit and the astronomical theory of climate variations. Trudy Inst. Theor. Astron., 14, 48-84. (原文はロシア語)

田 近 英 一(2007):全 球 凍 結 と 生 物 進 化.地 学 雑 誌, 116,79-94.

安 成 哲 三 ・柏 谷 健 二 編(1992):地 球 環 境 変 動 と ミラ ンコ ヴ ィ ッチ ・サ イ クル.古 今 書 院(2007年9月10日

受 付,2007年11月10日 受 理)