「地殻活動予測シミュレーションモデルの構築」研究計画の目標は、北米、太平洋、フィリピン海、ユーラシアの四つのプレートが複雑に相互作用する日本列島域を一つのシステムとしてモデル化し、常時観測網からの膨大な地殻活動データをリアルタイムで解析・同化することで、プレート相対運動によって駆動される広域応力の増加から準静的な震源核の形成を経て動的破壊の開始・伝播・停止に至る大地震発生過程の定量的な予測を行うことにある。本研究計画では、上記の目標を達成するために、「ア．日本列島域」を対象とした地殻活動予測シミュレーションモデルを開発する。また、稠密（ちゅうみつ）な観測が行われている「イ．特定の地域」においては、より詳細な地域モデルを開発し、地殻変動データや地震活動データをリアルタイムで取り込んだ予測シミュレーションを行う。さらに、「ウ．予測シミュレーションモデルの高度化」のために、地震発生の物理・化学過程に関する基礎的シミュレーション研究を推進する。
地殻活動シミュレーションモデルの原型（プロトタイプ）を高度化し、1968年（昭和43年）十勝沖地震（マグニチュード7.9）及び2003年（平成15年）十勝沖地震（マグニチュード8.0）を例として、プレート境界での準静的応力増加-動的破壊伝播-地震波動伝播の連成シミュレーションを行った（東京大学理学系研究科［課題番号：1502］，防災科学技術研究所［課題番号：3016］）。また、シミュレーションと地殻活動観測の融合に向け、直接的及び間接的に得られている既知の情報を併用した新しい逆解析手法を発展させ、地震のCMT（セントロイド・モーメント・テンソル：地震波形から求められた地震の位置・規模・発震機構のこと）解をデータとして日本列島域の地殻応力状態を推定する解析プログラムを開発し、北海道-東北地域の三次元地震発生応力場の推定を行った（東京大学理学系研究科［課題番号：1502］）。
太平洋-北米プレート境界面上の1968年十勝沖地震の震源域に大小二つのアスペリティを設定し、太平洋プレートの沈み込み運動に伴う震源域での準静的応力増加-動的破壊伝播の連成シミュレーションを行った（図40）。震源域のせん断応力は地震発生後120年で臨界状態に達した。これを初期状態として応力集中域に小さな擾乱を与え、動的破壊伝播過程のシミュレーションを行った。開始した動的破壊は急激に加速され、大小二つのアスペリティを破壊する大地震に発展した。これに対し、60年後の未だ応力が充分に増加していない状態で大きな擾乱を与え強制的に動的破壊を開始させると、開始した動的破壊は殆ど加速されることなく停止してしまった。このことは、何らかの擾乱がトリガーになって開始した動的破壊が大地震まで発展するか否は、震源域の応力状態に強く依存することを示している（Hashimoto et al., 2006）。また、2003年十勝沖地震については、準静的応力増加-動的破壊伝播-地震波動伝播の連成シミュレーションを行い、理論的に予測された最大地動速度分布と実際に観測された最大地動速度分布が概ね一致することを示した（Fukuyama et al., 2006）。
地震間の滑り遅れ分布（アスペリティ分布）を正しく推定するためには、直接的及び間接的に得られている既知の情報を観測データと結合して、地殻変動データを解析する新しい逆解析手法を開発した。この新しい方法で、日本列島全域を対象とする解析プログラムの開発を進めるとともに、平成17年度に開発したCMTデータ逆解析手法を北海道-東北地域の地震データに適用し、同地域の三次元地震発生応力場の推定を行った（Terakawa et al., 2006）。解析結果（図41）は、沈み込む太平洋プレート内と陸側プレート内の応力場のパターンの違いを明瞭に示している。また、千島海溝と日本海溝の接合点近くでは、応力パターンの複雑な変化が見られる。
地球シミュレータを活用して、南海トラフ沿い巨大地震発生サイクルの準静的及び動的シミュレーションモデルを開発した（Hori, 2006）。西南日本の下に沈み込むフィリピン海プレートの三次元形状モデルに構造探査結果に基づく摩擦特性の不均質分布を与えて地震発生サイクル・シミュレーションを行い、地表変位の時空間変化を計算した。その結果、固着の剥がれに伴って隆起域と沈降域の境界が15年間で数キロメートル程度移動することが分かった（図3）。また長期間のデータがある四国での水準測量データとの比較から、これまでシミュレーションで設定していた地震発生帯の下限が深すぎること、余効滑りの時間スケールが短すぎることが分かった（名古屋大学［課題番号：1704］、海洋研究開発機構［課題番号：4001］）。
余効滑り・余震域拡大の数値シミュレーションにより、余効滑りや余震域の拡大速度は定常的摩擦の速度依存性を支配するパラメータ（A-B）に依存することを明らかにした（Kato, 2007）。また、三陸沖と十勝沖での摩擦特性の違いについて考察し、シミュレーション結果は2003年十勝沖地震の余効滑り速度と応力の関係から推定したA-Bの値と調和的であることを示した。
中国鮮水河断層を例にとり、多数の断層セグメント間の相互作用を考慮した地震発生サイクル・シミュレーションを行った（Kato et al., 2007）。断層モデルは、過去の地震活動記録を参考に、全長350キロメートルで9つのセグメントから成るとした。地震が発生するセグメントでは速度弱化の摩擦特性を仮定し、セグメント間の摩擦特性は（a）強い速度強化と小さな特徴的滑り量、（b）弱い速度強化と大きな特徴的滑り量の二つの場合を考えた。いずれの場合もセグメント間で余効滑りが発生するが、（a）の場合は余効滑りの時間変化は対数関数で記述できるのに対して、（b）の場合は余効滑りの立ち上がり時間が短く、比較的早く余効滑りが減衰する。その結果、（a）では、隣接するセグメントでの地震発生の時間差は比較的長く、ばらつきが大きくなる。これに対して、（b）では、地震発生の時間差は比較的短くなる。各セグメントでの地震発生間隔の統計的性質は、（a）ではBPT（Brownian Passage Time）分布または対数正規分布でほぼ説明できるが、（b）ではそのような分布では説明できない（図42）。また、滑り量の時間変化については、時間予測モデルでも規模予測モデルでも上手く説明できなかった。（東京大学地震研究所［課題番号：1411］）
三陸沖の構造探査データに基づく粘弾性有限要素モデルを構築し、1896年（明治29年）陸羽地震（マグニチュード7.2）の余効変動を計算して水準点での余効変動の時系列や空間パターンと比較した。その結果、構造探査に基づくプレート境界形状を用いた場合、先行研究で指摘されている沈み込んだ太平洋プレートがマントル内の流動を妨げるような幾何学的効果は小さく、地表で観測される地殻変動の時定数を支配している主要因は地殻・上部マントルの弾性構造不均質があることを確認した。（海洋研究開発機構［課題番号：4001］）
滑り速度の減速率が時間的に急変する余効滑りの場合、従来用いてきた「時間に関する平滑化の度合いが解析期間にわたって一定」という先験的拘束条件を課すと、地震直後の速度は小さめに推定され、逆に充分減速した後の推定値にはノイズの影響が大きくでてしまうため、滑り速度の履歴やそれに基づく摩擦パラメータの適切な推定が困難であった。この問題を克服するため、新しい時間依存逆解析手法を開発し、2003年十勝沖地震発生後のGPSデータに適用し、余効滑り域の摩擦パラメータがより正確に推定できるようになった（図43）（東京大学地震研究所［課題番号：1412］）。
動的断層滑り過程における摩擦発熱、流体圧変化、空隙発生及びそれらの間の非線形相互作用についてシミュレーション研究を行った。空隙は断層の動的滑りに伴う衝撃によって生成するものと仮定し、摩擦発熱による流体圧変化を考慮すると、滑り弱化や滑り強化などの多様な動的断層滑りを、単一の無次元パラメータを用いて統一的に説明することができた。滑りが極めて短時間で終了するため流体の流出が無視できるとして、流体圧、温度、断層滑りの時間変化の解析解を求めた。求められた滑り弱化距離は、地震学的に推定されたものとよく一致した。さらに、上述のパラメータがある臨界値より大きい場合には、パルス状の断層滑りが生じることも分かった（東京大学地震研究所［課題番号：1412］、Suzuki and Yamashita, 2006）。
離散要素法を用いて、ガウジ（断層破砕物）を挟んだ断層の運動について岩石の溶融を考慮した大変位高速滑りシミュレーションを行った。平成18年度は、溶融物の発生を空隙率の上昇とみなした簡単なモデルで、岩石の溶融が破壊特性に及ぼす影響を調べた。その結果、ガウジ層の空隙率がある臨界値（36パーセント）よりも大きくなるとせん断応力が大幅に低下すること、また摩擦強度の速度依存性が臨界値を境にして定性的に異なることを明らかにした。（東京大学地震研究所［課題番号：1412］）
平成9年に実施された東北日本横断地殻構造探査による地殻・上部マントル構造と温度の空間分布を考慮して、応力集中過程のモデル化を行った。下部地殻及び上部マントルの温度条件下では岩石の変形は流動則に従うと考え、非線形粘弾塑性有限要素法を用いてシミュレーションを行った。モデル領域の側面から短縮変形を加えて圧縮応力場を作り、地殻がどのように変形するか調べた。このモデル化では応力が高まると非線形流動が生じ、さらに応力が増加すると塑性変形が開始する。シミュレーションでは、奥羽脊梁山地及び出羽山地付近で高温のためにリソスフェアが薄くなり、上部地殻の下部に応力が集中する様子が再現された（図4）（東京大学地震研究所［課題番号：1412］）。
また、下部地殻での非線形流動による短縮変形により、3本の断層帯が形成された。変位の垂直成分からは、奥羽脊梁山地及び出羽山地を含む大規模なポップアップ構造の形成が確認でき、特に脊梁山地付近では、大きな隆起運動が生じている。このシミュレーションは、上平、千屋及び北由利の三つの断層帯の運動を上手く説明している（東京大学地震研究所［課題番号：1412］、Shibazaki et al., 2007）。
有限要素法を用いて、2000年（平成12年）の伊豆諸島の火山活動と東海地域の長期的ゆっくり滑りの関係を調べる数値シミュレーションを行った。先ず、東海地域の長期的ゆっくり滑りの周期、滑り速度が再現できるように、モデルパラメータの最適値を求めた。次に、この地震モデルと伊豆諸島の火山活動による変形量を用いて、長期的ゆっくり滑りの再現を試みた。その結果、幾つかのケースで長期的ゆっくり滑りが再現されたが、観測された特性とは少し異なった。（国土地理院［課題番号：6024］）
これまで、地殻・マントルの弾性-粘弾性構造、プレート境界の三次元形状、断層構成則の環境依存性を考慮した、日本列島域の地殻活動シミュレーションモデルの原型を「地球シミュレータ」上に構築し、プレート境界での準静的応力増加-動的破壊伝播-地震波動伝播の連成シミュレーションに成功した。しかし、このモデルによるシミュレーションが真に予測としての意味を持つためには、広域地震／地殻変動観測データから予測シミュレーションに必要な情報を様々な方法で抽出する必要がある。地震発生サイクル・シミュレーションなどを通じて、プレート境界面の摩擦構成則パラメータの空間分布についての知識は、次第に充実したものになってきつつある。また、地殻変動データ逆解析プログラムが完成すれば、地殻変動データから日本列島周辺域のプレート境界での滑り分布の変動履歴を復元することが可能になり、プレート境界及び周辺域の地殻応力状態の変化をシミュレーションにより予測することができる。他方、CMTデータ逆解析手法を広域地震活動データに適用すれば、地殻内の応力状態をモニターすることも充分可能となるであろう。従って、本研究計画を更に促進するためには、これまでの研究成果を集約した標準モデルを設定すべき段階にきている。
図40
1968年十勝沖地震の震源域での準静的応力蓄積-動的破壊伝播の連成シミュレーション（Hashimoto et. al., 2006）。1968年の大地震発生から60年後，120年後の応力蓄積状態（上段）における動的破壊伝播のスナップショット（下段）（東京大学地震研究所［課題番号：1502］、防災科学技術研究所［課題番号：3016］）。
図41
CMTデータインバージョンで推定した北海道-東北地域の三次元地震発生応力場（Terakawa et al., 2006）。左：震源球で表した地震発生応力場のパターン。右：最大主圧縮軸の方向。いずれも深さ10キロメートルでの水平断面を示す。図の背景色は，逆解析結果の推定誤差をカラースケールで表示したもので，青-緑-赤の順で推定誤差が大きくなる（東京大学理学系研究科［課題番号：1502］）。
図42
シミュレーションで得られた地震の繰り返し間隔の分布関数（Kato et al., 2007）。上段：セグメント境界域でA-B
0の摩擦パラメターを仮定。下段：セグメント境界域で非常に大きいL（
1メートル）を仮定。3つのセグメントS2，S4，S8での地震について結果を示した。赤線，青線，緑線は，それぞれ，BPT分布，対数正規分布，ワイブル分布。COV（α）は標準偏差/平均値を表す（東京大学地震研究所［課題番号：1411］）。
図43
新しい時間依存逆解析法による2003年十勝沖地震の余効滑りの解析。左：90日間の余効滑りによる累積滑り量。右：逆解析結果から推定されたせん断応力の速度依存性は、A-B0.07MPa（メガパスカル）となった（東京大学地震研究所［課題番号：1412］）。