2．（1）「地殻活動予測シミュレーションモデルの構築」研究計画

　「地殻活動予測シミュレーションモデルの構築」研究計画の目標は、北米、太平洋、フィリピン海、ユーラシアの四つのプレートが複雑に相互作用する日本列島域を一つのシステムとしてモデル化し、常時観測網からの膨大な地殻活動データをリアルタイムで解析・同化することで、プレート相対運動によって駆動される広域応力の増加から準静的な震源核の形成を経て動的破壊の開始・伝播・停止に至る大地震発生過程の定量的な予測を行うことである。本研究計画では、上記の目標を達成するために、「ア．日本列島域」を対象とした地殻活動予測シミュレーションモデルを開発した。また、稠密（ちゅうみつ）な観測が行われている「イ．特定の地域」においては、より詳細な地域モデルを開発し、地殻変動データや地震活動データをリアルタイムで取り込んだ予測シミュレーションを行った。さらに、「ウ．予測シミュレーションモデルの高度化」のために、地震発生の物理・化学過程に関する基礎的シミュレーション研究を推進した。

ア．日本列島域

　（日本列島域の地殻活動シミュレーションモデルのプロトタイプ（CAMPモデル）の構築と準静的応力蓄積過程‐動的破壊‐波動伝播の連成シミュレーションの試み）
　プレート運動による準静的歪蓄積モデルと動的地震破壊伝播モデルをシステム結合することにより、地殻マントルの弾性‐粘弾性構造、プレート境界の三次元形状、断層構成則の環境依存性等を考慮した、日本列島域の地殻活動シミュレーションモデルのプロトタイプ（CAMPモデル）を構築した（東京大学理学系研究科［課題番号：1502］）。このモデルを用いて、1968年（昭和43年）十勝沖地震（M7.9）の震源域における巨大地震の発生予測シミュレーションを試み、震源域の応力が充分に大きくなっていないと何らかの要因で動的破壊が開始されたとしても大地震には発展しないこと、逆に震源域の応力が臨界状態に達していると動的破壊は加速され大地震に発展することが再現できた（東京大学理学系研究科［課題番号：1502］、防災科学技術研究所［課題番号：3016］）。また、2003年十勝沖地震について、準静的応力増加‐動的破壊伝播‐地震波動伝播の連成シミュレーションを行い、理論的に予測された地震波形と実際に観測された地震波形が概ね一致することを示した（図32）（東京大学理学系研究科［課題番号：1502］、防災科学技術研究所［課題番号：3016］、Fukuyama et al., 2009）。

（GPSデータ逆解析による北海道・東北地域のプレート境界面の固着状態の推定）
　直接的及び間接的先験情報を観測データと結合した統計モデルを構築し、地殻変動データ解析の新しい逆解析手法を定式化することに成功した。この新しい定式化に基づく、北海道‐東北地域のGPSデータ解析により見出された、プレート境界面上での大きな滑り遅れ領域は、過去に発生した大地震の震源域とほぼ一致し、過去の大地震の震源域は現在固着していることが分かった（図33）（東京大学理学系研究科 [課題番号：1502]、Hashimoto et al., 2009）。

（CMTデータ逆解析による日本列島全域の三次元地殻応力パターンの推定）
　地震破壊は震源域周辺の応力場を反映するという考えに基づき、地震のCMTデータから地震発生応力場を推定する逆解析手法を開発した（Terakawa and Matsu’ura, 2008）。この新しい応力逆解析手法を防災科学技術研究所の15,000個の地震のCMTデータに適用し、日本列島域の三次元地殻応力分布を求めた（図34）（東京大学理学系研究科 [課題番号：1502]）。

イ．特定の地域

（三陸沖プレート境界型地震発生サイクル・シミュレーションモデルの構築）
　三陸沖のプレート境界では、1968年十勝沖地震の際に破壊された二つのアスペリティのうち一つだけが1994年三陸はるか沖地震（M7.6）で破壊された。このような地震サイクルは二つのアスペリティの摩擦特性の違いにより説明可能である。摩擦特性の不均一性を考慮したシミュレーションにより、1994年三陸はるか沖地震の余効滑りと最大余震（M7.2）を含めて、三陸沖のプレート境界地震の繰り返しを説明するアスペリティモデルを構築した。これにより、摩擦特性の大まかな分布について知見が得られた（図35、36）（東京大学地震研究所［課題番号：1411］、Kato,2008）。 また、速度・状態依存摩擦則と媒質の粘弾性をともに考慮したモデルを構築した。これにより、粘弾性による応力変化を考慮した地震サイクルモデルや、余効滑りと媒質の粘弾性変形をともに考慮した地殻変動のモデルに関する研究を進めることができるようになった（東京大学地震研究所［課題番号：1411］、海洋研究開発機構[課題番号：4001]）。

（南海トラフ沿い巨大地震発生サイクル・シミュレーションモデルの構築）
　地球シミュレータを用いた半無限弾性媒質中の平面断層モデルで、南海トラフに沿って変化するプレート相対運動速度、プレート境界面の三次元形状、沈み込む海山等による摩擦パラメータの不均質空間分布等を考慮した数値シミュレーションを行い、南海トラフ沿いの過去の巨大地震発生系列の主要な特徴を再現することができた。1944年東南海地震や1946年（昭和21年）南海地震（M8.0）のように、フィリピン海プレートが高角度で沈み込む紀伊半島沖で破壊が開始する可能性が高いことがシミュレーション結果から示唆された（名古屋大学[課題番号：1704]、東京大学地震研究所［課題番号：1411］、海洋研究開発機構[課題番号：4001]、気象庁[課題番号：7008]、高山・他、2008）。
　しかしながら、シミュレーションにプレート三次元形状を導入すると歴史地震で生じている連動・非連動が再現できなかった。その他にも、平面モデルでも再来間隔が平均値の半分程度から二倍程度までの範囲で大きく変化することが再現できないなどの課題があった。摩擦特性のうち速度依存性を規定するパラメータではなく、特徴的滑り量に不均一を与えることで、上記の問題が解決可能であることが分かった（名古屋大学 [課題番号：1704]、海洋研究開発機構[課題番号：4001]）。
　東海地域に発生したゆっくり滑りと想定東海地震の関連について、東海地域におけるゆっくり滑りの領域、滑り量などの観測結果と整合するようモデルを構築し、ゆっくり滑りが繰り返した後に東海地震に至る可能性を指摘した（気象庁[課題番号：7008]、弘瀬・他,2008）。さらに、ゆっくり滑りの発生要因について、プレートの形状に注目した研究（名古屋大学 [課題番号：1704]）や、遷移領域の摩擦構成則を考慮したモデル化が行われた。ゆっくり滑りについてのシミュレーション研究について、これまでに提案されているモデルを整理した （東京大学地震研究所［課題番号：1411］、芝崎，2009）。

ウ．予測シミュレーションモデルの高度化

　摩擦パラメータの推定や、現在の地殻活動予測シミュレーションでは考慮されていない以下の物理過程（間隙流体・非弾性）の組み込みなど、予測シミュレーションの高度化を行った。

（GPSデータ時間依存逆解析（余効変動）による摩擦パラメータ推定）
　GPSデータを用いて地殻応力変化やプレート境界面上の摩擦パラメータを推定する手法の開発を進め、実際のデータに適用して、その手法の有効性を検証した。2003年十勝沖地震の余効滑りと2000年から始まった東海地方でのゆっくり滑り等について、GPSデータからプレート境界面上の滑りと応力の時空間変化を復元し、摩擦構成則パラメータの推定を行った（東京大学地震研究所［課題番号：1412］、Fukuda et al., 2008）。

（熱多孔質媒質中での破壊伝播シミュレーション）
　熱多孔性媒質を仮定して自発成長する二次元動的破壊の理論的研究を行った。数値シミュレーションから、摩擦発熱を考慮した場合、空隙率の非弾性的な増加率とせん断変形帯幅の積に比例する量である無次元パラメータを用いて多様な動的断層滑りを統一的に説明できることが分かった（図37）。いくつかの観測事実と照らし合わせることにより、このパラメータのとるべき範囲についての推定を行った。この範囲内においては、ほぼ、一般的に滑り強化が起きていることが確認された（東京大学地震研究所［課題番号：1412］）。

（離散要素法による粉体シミュレーション）
　断層の摩擦構成則の研究として、断層ガウジ（断層帯内の細粒破砕物）を模擬した粉体層の摩擦強度回復過程をモデル化した計算コードを使用して不安定滑りの動的過程の数値シミュレーションを行い、薄い粉体層では地震的な高速滑りと非地震的な低速滑りの両方が発生することが発見された（図38）。低速滑りは、粉体層の速度強化的な摩擦法則に起因することが分かった。粉体層が滑りの全過程において弱化を続けるため、臨界滑り長が全滑り量に等しくなることも発見された（東京大学地震研究所［課題番号：1412］、Hatano, 2008）。

（内陸地震のモデル化に向けた粘弾塑性有限要素法（FEM）コードの開発）
　東北日本脊梁山脈を対象に、非線形粘弾性と塑性を考慮した有限要素法により、変形特性の空間的不均一性を考慮した三次元の変形と断層形成過程のモデル化を行った（図39）（Shibazaki et al., 2008）。シミュレーション結果から、北上低地西縁断層帯と横手盆地東縁断層帯に相当する断層が形成され、北部と南部に存在する高温地帯で断層帯の走向が高温地帯に向かうように変化することがわかった（図40）。2008年岩手・宮城内陸地震は、北上低地西縁断層帯の南部延長から、火山地帯に近づくような走向で破壊が生じたが、これはシミュレーションで得られた断層の走向と調和的であった（東京大学地震研究所［課題番号：1412］）。

（有限要素法およびその拡張によるプレート境界域周辺の地殻変動シミュレーション）
　2004年9月5日に発生した紀伊半島南東沖の地震の余効変動に関して、沈み込み帯の三次元構造を考慮して、粘性緩和による変動を有限要素法による数値計算により推定した。東海および紀伊半島周辺のGPS連続観測の時系列データから、粘性率が1018～1019Pa・sであると推定された。地震発生後三年間で粘性緩和による変動量は、志摩半島周辺で南方向に1cm弱、東海地方周辺で南西方向に5mm弱の大きさと推定され、2000年秋以降浜名湖直下で発生しているゆっくり滑りによる変動の推定には、余効変動の影響を取り除く必要があることが分かった（国土地理院［課題番号：6024］、水藤・小沢,2009）。

（不均質媒質における破壊伝播シミュレーションモデルの高度化）
　FEM‐β法により、不均質媒質中における地震断層の準静的成長に関する定量的シミュレーションを行い、地表堆積物の層構造での横ずれ断層成長停止機構を調べた（Kame et al., 2008）。地表堆積層の剛性率と初期せん断応力を深部の半無限媒質のそれらの1/10とし、破壊伝播による応力を計算して、境界面に向かう破壊伝播のシミュレーションを行った。破壊が境界を横切ると歪が連続であっても弾性定数が不連続に1/10になることにより応力レベルが低下して、断層成長が層の途中で停止する可能性がある。また、不均質媒質中において非平面地震断層の動的破壊解析を行うため、境界積分法と有限差分法を組み合わせた新たな数値解析手法の開発を行った（Kame and Aochi, 2009）（東京大学地震研究所［課題番号：1412］）。

課題と展望

　これまで、地殻・マントルの弾性‐粘弾性構造、プレート境界の三次元形状、断層構成則の環境依存性を考慮した、日本列島域の地殻活動シミュレーションモデルの原型を地球シミュレータ上に構築し、プレート境界での準静的応力増加‐動的破壊伝播‐地震波動伝播の連成シミュレーションに成功した。しかし、このモデルによるシミュレーションが真に予測としての意味を持つためには、広域地震／地殻変動観測データから予測シミュレーションに必要な情報を様々な方法で抽出する必要がある。地震発生サイクルシミュレーションなどを通じて、プレート境界面の摩擦構成則パラメータの空間分布についての知識は次第に充実したものになってきている。また、地殻変動データ逆解析により、地殻変動データから日本列島周辺域のプレート境界での滑り分布の変動履歴を復元することが可能になりつつある。また、地震のCMTデータを用いた応力逆解析手法を広域地震活動データに適用する研究も進められている。今後は、これまでの研究成果を集約した標準モデルを設定し、各種データとシミュレーションを統合した研究をより充実させる必要がある。
　現在の地殻活動予測シミュレーションでは考慮されていない物理過程として、流体の移動および空隙の非弾性的生成を考慮した破壊伝播のシミュレーションなどが行われた。また、粉体のシミュレーションも行われ、微視的なものも含め、破壊や摩擦の物理過程に関する理解は深まった。しかしながら、これらを大規模な地殻活動予測シミュレーションモデルに組み込む努力は不十分であり、今後の重要な課題である。また、不均質媒質中での破壊伝播を扱う新たな数値解析手法が開発され、今後の応用に期待できる。地殻活動予測シミュレーションでは、主としてプレート境界型地震を対象としているが、有限要素法により内陸地震発生機構のためにモデルも開発されている。これについても今後の進展に期待がもてる。

参考文献

　Fukuda, J., S. Miyazaki, T. Higuchi, and T. Kato, Geodetic inversion for space  time distribution of fault slip with time‐varying smoothing regularization, Geophys. J. Int., 173, 25  48, 2008.
　Fukuyama, E., R. Ando, C. Hashimoto, S. Aoi and M. Matsu'ura, Simulation of the 2003 Tokachi‐oki, Japan, earthquake to predict strong ground motions, Bull. Seismol. Soc. Am., 99（6）, doi:10.1785/0120080040, 2009 （in press）.
　Hashimoto, C., A. Noda, T. Sagiya, and M. Matsu'ura, Interplate seismogenic zones along the Kuril‐Japan trench inferred from GPS data inversion, Nature Geoscience, 2, 141‐144, 2009.
　Hatano, T., Scaling Properties of Granular Rheology near the Jamming Transition, J. Phys. Soc. Jpn., 77（12）, 123002, 2008.
　弘瀬冬樹・前田憲二・高山博之, 東海地域の長期的スロースリップイベントの再現‐その3‐, 日本地球惑星科学連合2008年大会予稿集, S142‐P010, 2008.
　Kame, N., S. Saito, and K. Oguni, Quasi‐static analysis of strike fault growth in layered media,　Geophys. J. Int., 173, 309‐314, 2008.
　Kame, N. and H. Aochi, A hybrid FDM‐BIEM approach for earthquake dynamic rupture simulation, 12th International Conference on Fracture Proceedings, Ottawa, Canada, 2009, accepted.
　Kato, N., Numerical simulation of recurrence of asperity rupture in the Sanriku region, northeastern Japan, J. Geophys. Res., 113, B06302, doi:10.1029/2007JB005515, 2008.
　芝崎文一郎，沈み込み帯深部で発生するスロースリップイベントのモデル化，地震2，受理，2009.
　Shibazaki, B., K. Garatani, T. Iwasaki, A. Tanaka, and Y. Iio, Faulting processes controlled by the nonlinear flow　in the deeper crust and upper mantle beneath the　northeastern Japanese island arc, J. Geophys. Res., 113, B08415, doi:10.1029/2007JB005361, 2008.
　水藤尚・小沢慎三郎，東海地方の非定常地殻変動―東海スロースリップと2004年紀伊半島南東沖の地震の余効変動，地震，61, 113‐135, 2009.
　高山博之・前田憲二・弘瀬冬樹，南海トラフ大地震の開始位置に与えるプレート境界の形状の効果，地震，60, 279‐284, 2008.
　Terakawa, T. and M. Matsu’ura, CMT data inversion using a Bayesian information criterion to estimate seismogenic stress fields, Geophys. J. Int., 172, 674‐685, 2008.

図32：2003年十勝沖地震における準静的応力蓄積‐動的破壊‐波動伝播の練成シミュレーション。（左図）2003年十勝沖地震の震源域と観測点分布。（中図）震源域における滑り速度の時間発展と最終滑り変位。（右図）各観測点における波動伝播シミュレーションにより合成波形（赤）と観測波形（黒）（防災科学技術研究所［課題番号：3016］、Fukuyama et al., 2009）。

図33：GPSデータの逆解析で推定した北海道・東北地域の北米  太平洋プレート境界面の固着状態。（左図）解析に用いたGPS水平速度データと三角網。（右図）北米  太平洋プレート境界面での滑り遅れ（青）と滑り過剰（赤）速度分布。コンター間隔は3cm/yr。緑の星印と楕円、過去に発生したプレート間地震の震央位置および津波波源域を示す（東京大学理学系研究科 [課題番号：1502]、Hashimoto et al., 2009）。

図34：CMTデータ逆解析により推定された日本列島全域の三次元地殻応力パターン（東京大学理学系研究科 [課題番号：1502]、Terakawa and Matsu’ura, 2008）。（上図）震源球で表された応力場（青：圧縮場・逆断層解、緑：圧縮場・横ずれ断層解、伸張場・正断層解）、（下図a）琉球断面図：伸張場が卓越、（下図b）西南日本断面図：横ずれ断層が卓越、（下図c）東北日本断面図：逆断層が卓越。

図35：三陸沖地震発生サイクルシミュレーション（東京大学地震研究所［課題番号：1411］、Kato,2008）
　（左図）1968年十勝沖（中図のLA1）、1994年三陸はるか沖地震（中図LA1＆LA2）とモデル設定。（中図）摩擦パラメータ（A‐B）分布。SA1は三陸はるか沖地震の最大余震、SA2は青森沖地震。（右図）摩擦パラメータ（特徴的長さL）分布。

図36：三陸沖地震発生サイクルシミュレーション（東京大学地震研究所［課題番号：1411］、Kato,2008）
シミュレーション結果：プレート速度Vplで規格化した滑り速度のスナップショット。1994年三陸はるか沖に対応するLA1の単独破壊（d）、余効変動により引き起こされた余震に対応数SA1の破壊（i）、青森沖の地震に対応するSA2の破壊（n）、1968年十勝沖地震に対応するLA1とLA2の連動破壊（t,u）。

図37：熱多孔質媒質中で二次元動的破壊シミュレーション（東京大学地震研究所［課題番号：1412］）。
断層面上の滑りとせん断応力の関係。
kは透水率、Suは空隙率の非弾性的な増加率とせん断変形体幅の積に比例する無次元量である。Suの値により、滑り開始後に滑り弱化特性を示す場合と滑り強化特性を示す場合があることが分かった。

図38：離散要素シミュレーションで再現された粉体層の滑りの様子。黒線は滑り量、赤線は滑り速度を表す。急加速する初期過程と非常にゆっくりした減速を示す後期過程に分かれている（東京大学地震研究所［課題番号：1412］）。

図39：（上図）三次元のモデル。（下図）地温勾配の分布。赤色の部分で温度が高く、青色の部分で温度が低い（東京大学地震研究所［課題番号：1412］、Shibazaki et al., 2008）。

図40：5000年間に生じる（上図）全歪と（下図)隆起量（cm）（東京大学地震研究所［課題番号：1412］）。