インヤン格子を用いた三次元球殻マントル対流コードの開発

この記事は、「陰陽格子法を用いた新しい三次元球殻マントル対流数値シミュレーションコードの開発」（吉田晶樹・陰山聡、地球惑星科学関連学会2004年合同大会予稿集、2004年5月）に基づいて書かれています。

■はじめに

　地球のマントルは厚さ約2900kmの球殻の形状をしている（図参考）．全て岩石で構成されるマントルを非常に粘性の高い流体としてレイリー数を見積もると，約３０億年前の太古代から現在に至るまで地球のマントルは対流し続けてきたことはほぼ間違いない．マントル対流の最も簡単な数値モデルは，非圧縮・ブジネスク近似下の流体を対流層の下面（核・マントル境界に相当）で加熱し，上面（大気とマントルの境界）で冷却する，いわゆるレイリー・ベナール対流モデルである．厳密には地球の表面には厚さ100km程度のマントルより冷たく固い岩石からなるリソスフェアで覆われているが，これがマントル対流の上面熱境界層に相当するものとすると，やはりリソスフェアもマントルの一部分であると言える．マントル対流の数値シミュレーションを行って過去や現在のマントル対流の振る舞いを調べることは，究極的には地球変動の歴史を解明するための重要な手段の一つとなる．マントル対流の数値計算はプレートテクトニクス理論の登場直後の1960年代後半から二次元モデルで行われてきた．

図の説明：地球の断面図．黄色い部分が計算対象となるマントルを表す．赤い空洞の部分はコアを表す．

■球殻マントル対流コードの開発

　三次元球殻のモデルを用いたマントル対流の数値シミュレーションは世界では 1980年代半ばから行われてきた．マントル対流を記述する基礎方程式の離散化の方法としては，スペクトル法（Machetel et al., 1986; Glatzmaier, 1988; Bercovici et al., 1989; Harder and Christensen, 1996），有限体積法（Ratcliff et al., 1996; Iwase, 1996），有限要素法（Zhong et al., 2000; Richards et al., 2001）が用いられてきた．このうちスペクトル法は，流体内部で局所的に大きく変化する物理量を取り扱う場合，そのスペクトル展開に莫大な計算時間とメモリを要する．例えば，地球型惑星のマントルを構成する物質（岩石）の粘性率は温度や歪速度，化学組成の変化などによって水平方向に短い空間スケールで何桁も変化する．そのため，マントル対流の数値シミュレーションにおいてスペクトル法は現在あまり用いられていない．一方で，有限差分，有限体積，有限要素法は，粘性率またはその空間微分を各格子点で「局所的に」扱うことが可能である．
　一般的に，差分法系の手法（有限差分法と有限体積法）は地球シミュレータに代表される大規模ベクトル・並列計算機において計算効率が高いことが分かっている．

■「極」の問題の克服

　ところで，差分法系の手法を用いる場合，計算格子として，通常の緯度経度（θ-φ）格子（Ratcliff et al., 1996）を用いると低緯度領域と高緯度領域で格子点間隔が極端に不均一になるため，クーラン条件による時間刻みが高緯度領域の密な格子間隔で決まってしまう．この「極」の問題により，差分法系の手法で現実的なマントルダイナミクスの時間発展問題を地質学的時間スケール（数十億年）にわたって三次元球殻モデルを用いて解くマントル対流のシミュレーションはこれまで行われていなかった．そこで我々は「インヤン格子」を用いることにより新しい三次元球殻マントル数値シミュレーションコードの開発を行った（Yoshida and Kageyama, 2004）．

図の説明：インヤン格子は、球座標格子（移動経度格子）の低緯度領域を二つ組み合わせた格子。

■計算手法

　マントル対流を支配する無次元化された基礎方程式（質量保存式，運動量保存式，エネルギー方程式）を二次精度の有限差分法で離散化する．マントル対流の定常流れ場（速度場と圧力場）は，これまでのマントル対流計算で広く使用されているSIMPLERアルゴリズム（Patankar, 1980）を用いて，質量保存式と運動方程式を同時に満たすように逐次補正しながら陰的に求める．エネルギー方程式（温度場の輸送方程式）の時間差分は，二次精度のクランク・ニコルソン法を用いる．エネルギー方程式の移流項の空間差分は風上差分法を用いる．速度，圧力，温度の変数配置は離散化が簡単で境界条件が扱いやすいコロケート格子法を採用する．陰陽格子法においては，これらの基礎方程式は時間ステップごとに二つの格子系において別々に同時進行で解く．このとき一方の格子系上の境界値データは，SIMPLER計算と温度場計算のループとも，一回の反復計算ごとに他方の格子系上に重なる格子点から線形補間して与える．

■ベンチマークテスト

　開発したコードの妥当性を調べる．そのために，ブジネスク近似，レイリー・ベナール型の熱対流問題を解くことにより，これまでに三次元球殻マントル対流計算から得られている結果とのベンチマークテストを行った．温度場の初期条件は，静水圧下での熱伝導方程式と固定温度の境界条件を満たす温度場に「四面体型」（プルーム状の上昇流が四つ）と「六面体型」（プルーム状の上昇流が六つ）パターンの微小擾乱を加え，対流パターンが幾何学的対称・定常解（図参考）を示す範囲の低いレイリー数（Ra_c =< Ra =< 1.4×10⁴, Ra_cは臨界レイリー数で約712）を与えた．そして，最終的に得られた定常解でのヌッセルト数（対流による熱輸送効率を表す数）と対流速度の二乗和平均値を計算した．その結果，一様粘性率，粘性率が温度に依存する場合の両方において，我々が開発したコードと，上に列挙した従来のコードから得られた結果との間には，離散化法や数値計算手法，格子点数の違いがあるにも関わらず，数%以内の誤差で整合することが確認された．
　一様粘性率でレイリー数が低い場合（Ra = 10³～10⁵）において，詳細なパラメータサーチを行った結果，与えたレイリー数がRa_cの100倍程度で，幾何学的な対称性が崩れ非定常な対流パターンが得られることが分かった．通常の緯度経度格子を用いた有限体積法コード（Ratcliff et al, 1996）では，同じレイリー数で非定常な対流パターンは得られていない．これは，上に述べた格子間隔の極端な不均一性により，特に六面体型対流パターンでは，対流の上昇流が数値計算上の極域で幅広く，赤道域で幅狭くなるという非対称な定常対流のパターンを許してしまうためであると考えられる．

図の説明：低レイリー数における「四面体型」（左）と「六面体型」（右）のマントル対流の温度場．温度異常（各深さでの温度の水平平均からのずれ）の等値面の青色は周囲のマントルより冷たい下降流（ΔT= -0.125，値は無次元），黄色は周囲のマントルより温かい上昇流（ΔT= +0.150，値は無次元）を表す．中心の赤色の球はコアの表面を表す．

■球殻マントル対流コードによるテスト計算：内部発熱の影響

　開発したコードを用いて簡単なテスト計算を行う．
　マントルから地球の外に放出される熱量(44テラワット)は，マントル下面におけるコアからの加熱量とマントル内部に存在する放射性元素（ウラン，トリウム，カリウムなど）の壊変に伴う発熱量が大部分（70～80％程度）を占める（Turcotte and Schubert, 2002）．当然のことながら，この放射性発熱量は地球の歴史の中で徐々に減少する一途にある．これまでの数値シミュレーション結果から，マントル内部に熱源が存在する場合では，レイリー数が低く，対流セルのアスペクト比が小さい場合においても，下面加熱のみの場合の対称・定常的な解とは異なり，対流パターンは時間依存性をもつことが分かっている（例えば，Bercovici et al., 1989）．このように放射性元素による内部発熱の影響は，現実的な地球のマントル対流の非対称性や時間依存性を再現する上で重要な要素であると言える．
　ところが，先に述べた理由により，差分系の手法を用いたマントル対流シミュレーションコードで現実的なレイリー数と内部発熱の効果を考慮した計算はこれまで行われていなかった．そこで我々は，開発したコードを用いて，マントルの内部発熱を考慮したマントル対流シミュレーションを行い，内部発熱量を変化させたときの対流パターンの変化に注目した．
　計算の格子点数はインヤン格子の各格子系（N系とE系）で102(r)×54(θ)×158(φ)とする．マントル対流層上下面の熱的境界層は固定温度（下面で温め，上面で冷やす），力学的境界条件は自由滑り，不透過とした．内部発熱量のみの変化による対流パターンの影響を調べるため，マントルの粘性率やその他の物性値は全て一定とした．内部発熱量は，対流層上下面の温度差で無次元化されたパラメータ（H）をエネルギー方程式の生成項として与えた．簡単のため，放射性熱源はマントル内部に空間的，時間的に一様に分布している（つまり熱源の移動や量の変化なし）と仮定した．対流層上下面の温度差のみで定義されるレイリー数Raは実際の地球マントルのレイリー数に相当する10⁷とした．
　Ra=10⁷で内部発熱がなく加熱が下面のみからの場合（H=0），円筒状の上昇流とシート上の下降流からなる対流パターンが観察される（図参考）．レイリー数が一桁小さいRa=10⁶の場合と比較して対流パターンは短波長構造を示すことが分かる．ここで内部発熱の効果を考慮してみる．内部発熱量は， H = 10（現在の地球マントルの内部発熱量に近い値）， H = 20（約20～30億年前の地球マントルの内部発熱量に近い値）を用いることにする．内部発熱量が大きくなるにつれて，対流パターンはさらに短波長構造になることが分かる．内部発熱がかなり強いH = 20の場合では，数多くの円筒状の下降流とシート状の上昇流が見られる．内部発熱によりマントル内部の平均温度が上昇し，その結果下面の温度境界層が薄くなるために大規模な上昇流が起きにくくなる．このように内部発熱量の変化は，対流セルのアスペクト比のみならず上昇流や下降流の形態を変化させることがわかる．
　地震波トモグラフィーモデルで観測事実として得られているマントル対流パターンは，対流セルがマントルの厚さと同程度かそれ以上の波長が卓越する大規模構造で，円筒状の上昇流とシート状の下降流で構成されていることから，現在の地球のマントル対流パターンは， H=0（もしくは，我々が用いた放射性熱源が一様に分布しているモデルとは異なり，放射性熱源が核－マントル境界に近いマントル深部に集中している）の場合や H=10の場合で得られた状態にほぼ近いと考えられる．また，インヤン格子を用いたマントル対流コードでは， H=20の場合のような非常に短波長な構造をもつマントル対流も解くことが可能であることが分かった．

図の説明：高レイリー数と内部発熱を考慮したマントル対流の温度場． (左上)Ra=10⁶, H=0，(右上)Ra=10⁷, H=0， (左下)Ra=10⁷, H=10, (右下)Ra=10⁷, H=20の場合．カラーコンターは無次元温度を表す（カラーバー参考）．温度異常の等値面の青色は周囲のマントルより冷たい下降流（ΔT= -0.1），黄色は周囲のマントルより温かい上昇流（ΔT= +0.1）を表す．中心の赤色の球はコアの表面を表す．

■おわりに

　インヤン格子を用いることにより，有限差分法を用いて現実的なマントルダイナミクスの時間発展問題を地質学的時間スケールにわたって解くことに成功した．我々が開発したインヤン格子を用いたマントル対流コードは現在のところ，差分法系の手法で現実的なマントル対流問題を三次元球殻モデルで「正確」に解くための唯一のコードであると言える．

■参考文献

Masaki Yoshida and Akira Kageyama, Application of the Yin-Yang grid to a thermal convection of a Boussinesq fluid with infinite Prandtl number in a three-dimensional spherical shell, Geophysical Research Letters, 31(12), L12609, doi:10.1029/2004GL019970, 2004.
Masaki Yoshida and Akira Kageyama, Low-degree mantle convection with strongly temperature- and depth-dependent viscosity in a three-dimensional spherical shell, Journal of Geophysical Research, 111(B3), B03412, doi:10.1029/2005JB003905, 2006.
吉田晶樹, 陰山聡, 陰陽格子法を用いた新しい三次元球殻マントル対流数値シミュレーションコードの開発, 地球惑星科学関連学会2004年合同大会予稿集, 千葉, 2004年5月．
Masaki Yoshida and Akira Kageyama, Application of the Yin-Yang grid to three-dimensional spherical shell convection at infinite Prandtl number, EOS Transactions, American Geophysical Union, 85(47), Fall Meeting Supplement, 2004.
Masaki Yoshida and Akira Kageyama, Development of a new numerical finite difference code for a three-dimensional thermal convection with infinite Prandtl number in a spherical shell using Yin-Yang grid, Abstracts of 9th Symposium on the Study of the Earth's Deep Interior, 67-68, 2004.

トップページ