第7章 「地殻活動シミュレーション手法」研究計画
１．はじめに
地殻活動シミュレーション研究の目標は，隣接するプレート同士が複雑に相互作用する日本列島域の地殻活動のシミュレーションモデルを構築し，広域GPS観測網や地震観測網等からの膨大な地殻活動データをリアルタイムで解析・同化することにより，プレ−ト相対運動によって駆動されるテクトニック応力の蓄積から準静的な破壊核の形成を経て動的破壊の開始・伝播・停止に至る大地震発生過程の定量的予測を行うことにある．上記の目標を達成するためには，全国の大学及び関係諸機関が適切な役割分担の下に連携・協力し，複数の要素モデルをシステム結合した日本列島域を対象とするプロトタイプの地殻活動統合シミュレーションモデルを構築する一方，大学等の研究グループが中心となって，以下に挙げるようなモデリング及びシミュレーション手法の高度化のための基礎研究を重点的に推進し，その成果をプロトタイプモデルに逐次組み込むことで地殻活動予測シミュレーションモデルを継続的に改良・発展させていく必要がある．
大学等の研究グループが重点的に推進すべき基礎研究項目
１）断層破砕帯の素過程
２）断層間相互作用
３）内陸活断層の地震発生過程
４）地殻活動データの解析・同化
５）特定地域に於ける地震発生サイクルモデルの開発
６）日本列島域の広域変形・応力場のシミュレーション
地震予知のための新たな観測研究計画（第１次：平成11年度〜15年度）においては，建議の遅れもあり，新研究項目である地殻活動シミュレーション手法の開発に関して，建議の内容を実施する体制が充分整わないままスタートせざるを得なかったが，その後は段階的に研究体制を充実し，最終的には，第１次の地震予知研究事業として，以下の５研究課題を実施した．
第１次計画（平成11年度〜15年度）で実施した研究課題
0127）地殻内流体の挙動とその地震発生に対する力学的効果に関する研究
0128）断層間相互作用による断層成熟度の変化についての研究
0703）下部地殻流動特性とプレート内応力の蓄積・解放過程のシミュレーション研究
0908）海溝型巨大地震の地震サイクルモデリング研究
0120）地殻応力・歪変化シミュレーション手法に関する研究
また，これらのモデリング及びシミュレーション手法高度化のための基礎研究と並行して，全国のシミュレーション研究者の連携・協力により，複数の要素モデルをシステム結合した日本列島域を対象とするプロトタイプの地殻活動統合シミュレーションモデルの開発が進められた（Matsu’ura, 2001; 松浦, 2003）．
２．第１次計画（平成11年度〜15年度）の研究成果
「地震予知のための新たな観測研究」の第１次計画（平成11年度〜15年度）では，将来の定量的地震発生予測の実現に向けた，地殻活動統合シミュレーションモデルの開発とその高度化のための多様な基礎研究を実施し，以下のような成果が得られた．
2.1 断層破砕帯の素過程に関する研究
地殻内流体の挙動とその地震発生に対する力学的効果に関する研究（課題番号0127）を実施し，流体移動を考慮した地震活動の数値シミュレーションにより，破壊成長と流体移動の相互作用を考えることで，一見多様に見える地震破壊現象を統一的に説明できることを示した（Yamashita, 2003）．また，熱的・水力学的効果を考慮した動的破壊の数値シミュレーションにより，断層帯内の流体がすべり時間関数や破壊伝播の特徴的な振る舞いに関与している可能性を明らかにした（鈴木・山下，2003）．
破壊成長と流体移動の相互作用を考慮した数値シミュレーションでは，断層破壊が生じるとそこでの凝着強度は低下し透水係数は増大すると考えれば，余震活動に関して，大森公式とグーテンベルグ・リヒターのマグニチュードと発生頻度の関係式が統一的に説明できることが分かった．また，グーテンベルグ・リヒターの関係式を満たすものは繰り返しすべりを起こしている断層の破壊であること，余震時系列の初期に比較的大きなイベントが起きる傾向にあることなどが定量的に示された．一方，熱多孔質弾性体中の動的破壊伝播の数値シミュレーションでは，液相の熱膨張率が固相のそれに比べ非常に大きい場合は，破壊開始直後に短時間のパルス状の断層すべり（self-healing slip）が生じるが，熱膨張率の比が小さい場合には，最終的に定常すべりが実現されることが明らかにされた（図１）．
2.2 断層間相互作用に関する研究
断層間相互作用による断層成熟度の変化についての研究（課題番号0128）を実施し，断層形状や複数断層の相互作用が地震破壊の伝播過程や停止過程に及ぼす影響の重要性を明らかにするとともに，動的破壊による断層の成長・合体により断層形状や断層破砕帯が時間発展していく過程を数値シミュレーションによって示した（Kame & Yamashita, 2003; 安藤・山下, 2003）．
亀裂間の動的相互作用を考慮したシミュレーションでは，隣接する亀裂が互いに接近・合体するか否かは主に初期亀裂の空間配置によること，亀裂成長速度が大きくなると亀裂端の屈曲の程度は大きくなり，隣接する亀裂は合体しやすくなることを明らかにした．また，断層破砕帯のシミュレーションでは，断層破砕帯（主断層と多数の微小弱面で表現）の幅が主断層の破壊の進展に伴って増大すること，主断層の破壊に伴い微小弱面の破壊がS波によって励起されることを明らかにした（図２）．また，主断層破壊の初期段階において一旦加速した破壊が減速したり，最終的な破壊伝播速度がレーリー波速度以下に抑えられたりするのは，微小弱面の破壊により生じたstress shadowによる負の相互作用の影響であることを明らかにした．
また，より巨視的な観点からは，すべり速度と状態に依存する摩擦構成則（Kato & Tullis, 2003）を用いた地震サイクルシミュレーションを行い，アスペリティの破壊過程の多様性や複数のアスペリティ間の相互作用を明らかにした．まず，エピソディックな非地震性すべりに関する数値シミュレーションでは，様々な時定数をもつ非地震性すべりイベントや地震発生層よりも深部で発生する先駆的すべりを理解するのに役立つ結果を得た（Kato, 2003a, 2003b）．また，アスペリティ間の相互作用に関するシミュレーションでは，非地震性すべりがトリガーする遅れ破壊や複雑な地震サイクルを再現することに成功した（加藤，2003a）．更に，このモデルを用いて，1968年十勝沖地震と1994年三陸はるか沖地震が発生した三陸沖でのプレート境界地震発生サイクルシミュレーションを行い，２つのアスペリティが同時に破壊される巨大地震と１つのアスペリティだけが破壊されるやや小さめの地震とが交互に発生することやアスペリティの外側では顕著な余効すべりが発生することなど，三陸沖での地震発生サイクルについての幾つかの特徴を再現することができた（図３; 加藤，2003b）．また，2003年十勝沖地震の余効すべりの時空間変化をGPSデータから推定し，その結果を用いてプレート境界面上の応力とすべり速度の時間変化を調べることにより，プレート境界面上の摩擦構成則パラメータの推定を試みた（Miyazaki et al., 2004）．
2.3 内陸活断層の地震発生過程に関する研究
下部地殻流動特性とプレート内応力の蓄積・解放過程のシミュレーション研究（課題番号0703）を実施し，異方的な流動特性を持つ粘弾性物体の力学的応答の定式化を試みる一方，デタッチメントから分岐する断層の地質学的時間スケールでの挙動を数値シミュレーションにより明らかにした．また，現実的な３次元プレート境界モデルを用いて，太平洋プレート及びフィリピン海プレートの定常沈み込みに伴う日本列島域の地殻変形運動を再現するとともに（Hashimoto et al., 2004），地殻内応力の蓄積メカニズムについても考察し，プレート収束運動の１割が地殻内変形で解消される (衝突率10％) とすると，東北日本弧の応力蓄積パターンが合理的に説明されることを示した（図４；Matsu’ura et al., 2003）．
2.4 地殻活動データの解析・同化手法に関する研究
地殻応力・歪変化シミュレーション手法に関する研究（課題番号0120）を実施し，GPS等の測地測量データから地殻応力の変化を推定するためのインバージョン解析手法を開発した（Hori, 2000a, 2000b）．この解析手法は，観測された地表変位データからプレート境界での固着に起因する広域的変動を分離するための逆解析部分と局所的変位データから歪み変化を計算して応力変化を推定する逆解析部分とから成るが，前段の部分に関してはグリーン関数のスペクトル分解に基づく逆解析理論を，後段の部分については応力成分を生成する応力関数推定の逆解析理論を構築した．
この他にも，別の視点から，ABIC (Akaike’s Bayesian Information Criterion) を用いて測地測量データから大地震１サイクル間のプレート境界面のすべり履歴を推定するインバージョン解析手法（Fukahata et al., 2004）や観測された地震波形データから震源での破壊過程を偏り無く推定するインバージョン解析手法（Fukahata et al., 2003）が開発された．
2.5 特定地域に於ける地震発生サイクルモデルの開発に関する研究
海溝型巨大地震の地震サイクルモデリング研究（課題番号0908）を実施し，日本海溝沿いの東北日本と南海トラフ沿いの西南日本を対象として，プレート形状を考慮した２次元及び３次元の運動学的地震サイクルモデルを開発する一方，簡単なシステムですべり速度と状態に依存する摩擦構成則を用いた準静的地震サイクルシミュレーションを行った．
＜運動学的地震サイクルモデリング＞
大規模汎用有限要素法ソフトウエアGeoFEMの地震サイクルモジュール開発（Hirahara, 2002a, 2002b; Hyodo & Hirahara, 2004; 平原他，2003）の一環として，東北日本の３次元粘弾性FEMモデルを構築し，地震サイクルに伴う地殻変形・応力場の時間変化を運動学的に再現し，東北日本の最近100年間の水平歪み速度との比較から，1900年頃三陸沖でモーメントマグニチュード8.4に相当する大規模な非地震性すべりが発生していた可能性を見出した（Suito, et al., 2002）．また，西南日本についても３次元粘弾性FEMモデルを構築し，内陸部での東西圧縮場およびGPS観測で発見されたひずみ集中帯（NKTZ:新潟—神戸変動帯）の成因について定量的解析を行った結果，沈み込む太平洋プレートが系の緩和時間よりも長期にわたり固着していればプレート沈み込みによる圧縮力を内陸部へ効率的に伝達できること，またNKTZに歪みを集中させるには地殻の不均一が必要であることを定量的に示した（図５；Hyodo, & Hirahara, 2003）．
＜準静的地震サイクルシミュレーション＞
ブロックーバネモデルによる南海トラフ沿いの巨大地震の発生シミュレーションを行い，ブロックA〜Dがほぼ同時にすべる場合が2/3程度の頻度で発生するのに対し，交互にすべる場合が1/3程度の頻度で発生することが明らかになった．このことは，発生パターンの大きく異なる２つの地震サイクル状態が存在し得ることを示している（Mitsui & Hirahara,, 2004）．また，海溝軸と平行な方向に均質な摩擦パラメータ分布を持つ沈み込みプレート境界での３次元地震発生サイクルシミュレーションを行い，海溝に沿った地震発生帯の長さHと地震発生様式の関係を調べた結果， H < 300 kmの場合は走向方向の中央部で固有地震が周期的に発生するのに対し，H = 1000 kmとなると色々な場所で色々な大きさを持つ通常地震とゆっくり地震が発生するようになり，中間の 400 km < H < 800 kmでは遷移的な傾向を示すことが分かった（図６）．こうした地震発生様式の多様性は自発的に形成される不均質な応力分布に起因しているものと理解される（Hirose & Hirahara, 2002）．
2.6 日本列島域の広域変形・応力場のシミュレーションに関する研究
地殻応力・歪変化シミュレーション手法に関する研究（課題番号0120）を実施し，GPS観測から得られる変位データから逆解析手法（Hori, 2000a, 2000b）を用いて応力変化をモニタリングするシステムを構築した(Hori et al., 2000)．この解析手法を日本列島域に適用して応力変化の空間分布を推定し，それを歪変化の空間分布と比較することで，剛性率の小さな領域で地震活動度が高いという結果を得た（Hori et al., 2001）．しかし，この結果の妥当性は，データ解析の際に用いた仮定の妥当性に強く依存するので，多少問題があった．そこで，今度は，変位増分の時系列データから構成法則のパラメータ値を推定する手法（弾性係数逆解析手法）を開発し，GPSデータから日本列島域の各地域のポアソン比分布の推定を試みた（図７）．但し，データに含まれる誤差によりポアソン比が有意に決まらない場合は，ポアソン比を０としている．
2.7 日本列島域の地殻活動予測シミュレーションモデルの開発に関する研究
平成10年度に５ヶ年計画としてスタートした固体地球シミュレータ計画（Matsu’ura, 2004）の一環として，平成12年度末までには，強度回復メカニズムを内包する断層構成則（Aochi & Matsu’ura, 2002）を用いた横ずれプレート境界での３次元準静的応力蓄積モデルを完成させ（Hashimoto & Matsu’ura, 2000, 2002），全長1200kmのサン・アンドレアス断層系に適用できるように拡張し，地球シミュレータ上での最適化を行った（Nakajima, 2003）．また，それを動的破壊伝播モデル（Aochi et al., 2000a, b）とシステム結合することにより，テクトニック応力の蓄積から破壊核の形成を経て動的破壊の開始・伝播・停止に至る地震発生サイクルの全過程のシミュレーションを可能にした（Fukuyama et al., 2002）．平成13年度後半から日本列島域の地殻活動シミュレーションモデルの開発に着手し，平成14年度にはモデリングのベースとなる日本列島周辺域の３次元プレート境界形状モデルを完成させ（Hashimoto et al., 2004），GPSデータのインバージョン解析により，日本列島周辺域のプレート境界におけるすべり遅れ分布（固着域の分布）を推定した．平成15年度には，シミュレーション・コードのベクトル化及び並列化を進め，太平洋プレートの沈み込みに伴う1968年十勝沖地震の震源域での準静的応力蓄積過程のシミュレーションと３次元屈曲断層での破壊伝播シミュレーションをシステム結合することにより，2003年時点で次の十勝沖地震が発生する可能性の検証を行った（図８；Matsu’ura et al., 2003）．
３．まとめ
平成11年度にわずか１つの研究課題（0120）でスタートした本研究計画は，その後平成12年度・13年度に４つの新たな研究課題（0127, 0128, 0703,0908）を加えることで，最終的にはモデリングやシミュレーション手法の高度化のための基礎研究を相当に活性化することができた．一方，固体地球シミュレータ計画（平成10年度〜14年度）の一環として，全国のシミュレーション研究者の連携・協力の下に進められた日本列島域の地殻活動シミュレーションモデルの開発は，平成13年度末の地球シミュレータの完成により加速され，平成15年度末にはそのプロトタイプがほぼ完成した（橋本他，2003）．これらの研究成果は，平成16年度からスタートする「地震予知のための新たな観測研究計画（第２次）」の項目２-(1) 地殻活動予測シミュレーションモデルの構築に関する５ヶ年の研究計画を推進していく上で極めて重要な基礎となるものである．
参考文献
安藤亮輔・山下輝夫，断層間の動力学的相互作用と断層形状の生成，地震，56, 1-10, 2003．
Aochi, H., E. Fukuyama, and M. Matsu’ura, Spontaneous rupture propagation on a non-planar fault in 3D elastic medium, Pure Appl. Geophys., 157, 2003-2027, 2000a.
Aochi, H.,
Aochi, H. and M. Matsu’ura, Slip- and time-dependent fault constitutie law and its significance in earthquake generation cycles, Pure and Applied Geophysics, 159, 2029-2047, 2002.
Fukahata, Y., Y. Yagi, and M. Matsu’ura, Waveform inversion of earthquake rupture processes using ABIC with two sorts of prior constraints: Comparison between proper and improper formulations, Geophys. Res. Lett., 30, 1305, doi:10.1029/2002GL016293, 2003.
Fukahata, Y., A. Nishitani, and M. Matsu’ura, Geodetic data inversion using ABIC to estimate slip history during one earthquake cycle with viscoelastic slip-response functions, Geophys. J. Int., 156, 140-153, 2004.
Hashimoto, C., K. Fukui, and M. Matsu’ura, 3-D
modelling of plate interfaces and numerical simulation of long-term crustal
deformation in and around
Hashimoto, C. and M. Matsu’ura, 3-D physical modelling of stress accumulation processes at transcurrent plate boundaries, Pure Appl. Geophys., 157, 2125-2147, 2000.
Hashimoto, C. and M. Matsu’ura, 3-D simulation of earthquake generation cycles and evolution of fault constitutive properties, Pure Appl. Geophys., 159, 2175-2199, 2002.
橋本千尋・中島研吾・福井健史・佐藤利典・岩崎貴哉・松浦充宏，日本列島域の地殻活動予測シミュレーションモデルの開発，月刊地球，Vol.25, No.9, 675-681, 2003．
Hirahara, K., Interplate earthquake fault slip during periodic earthquake cycle in a viscoelastic medium at a subduction zone, Pure Appl. Geophys., 159, 2201-2220, 2002a.
Hirahara, K., A simple review on the simulation of earthquake cycle at subduction zones, in Seismotectonics in Convergent Plate Boundary, edited by Y. Fujinawa and A. Yoshida, 273-282, Terra Scientific Publishing Company, Tokyo, 2002b.
平原和朗・水藤尚・兵藤守・光井能麻・飯塚幹夫・殷峻・加藤尚之・宮武隆・堀高峰，GeoFEM地震サイクルシミュレーション，月刊地球，Vol.25, No.9, 694-698, 2003.
Hirose, H. and K. Hirahara, A model for complex slip behavior on a large asperity at subduction zones, Geophys. Res. Lett., 29, 22, 2968, 10.1029/2002GL015825, 2002.
Hori, M., Inversion of stress and constitutive relations using strain data for Japanese Island, Int. Symposium on Inverse Problems in Eng. Mech. II, 349-358, 2000a.
Hori, M., Inversion method using spectral decomposition of Green’s function, IUTAM-Symposium on Field Analyses for Determination of Material Parameters: Experimental and Numerical Aspects, IUTAM, 2000b.
Hori, M., T. Kameda, and T. Kato, Prediction of stress
field in
Hori, M., T. Kameda, and T. Kato, Application of the
inversion method to a GPS network for estimating the stress increment in
Hyodo, M. and K. Hirahara, A viscoelastic model of interseismic strain accumulation in Niigata-Kobe Tectonic Zone of central Japan, Earth Planets Space, 55, 667-675, 2003.
Hyodo, M. and K. Hirahara, GeoFEM kinematic earthquake
cycle simulation in southwest
Kame, N. and T. Yamashita, Dynamic branching, arresting of rupture and the seismic wave radiation in self-chosen crack path modelling, Geophys. J. Int., 155, 1042-1050, 2003.
Kato, N., Repeating slip events at a circular asperity:
Numerical simulation with a rate- and state-dependent friction law, Bull.
Earthq. Res. Inst., Univ.
Kato, N., A possible model for large preseismic slip on a deeper extension of a seismic rupture plane, Earth Planet. Sci. Lett., 216, 17-25, 2003b.
加藤尚之，アスペリティの相互作用に関する数値シミュレーション，月刊地球，25, 699-703, 2003a．
加藤尚之，プレート境界面の摩擦パラメターの推定について —三陸沖のアスペリティを例として—，地学雑誌，112, 857-868, 2003b．
Kato, N. and T. E. Tullis, Numerical simulation of seismic cycles with a composite rate- and state-dependent friction law, Bull. Seismol. Soc. Am., 93, 841-853, 2003.
Matsu’ura, M., The crustal activity modelling program: Progress toward scientific forecast of earthquake generation, Proceedings of the 2nd ACES Workshop, 23-26, 2001.
松浦充宏，地殻活動の予測シミュレーションとモニタリング，月刊地球，Vol.25, No.10, 767-772, 2003．
Matsu’ura, M., C. Hashimoto, and
Mitsui, N. and K. Hirahara, Simple spring-mass model
simulation of earthquake cycle along the Nankai trough, southwest
Nakajima, K., C. Hashimoto and M. Matsu’ura, Parallel performance of tectonic loading process model at transcurrent plate boundaries, Proceedings of the 3rd ACES Workshop, 243-248, 2003.
Suito, H., M. Iizuka and K. Hirahara, 3-D viscoelastic FEM modeling of crustal deformation in northeast Japan, Pure Appl. Geophys., 159, 2239-2260, 2002.
鈴木岳人・山下輝夫，地震の初期破壊に対する温度・流体圧・空隙率の非線形な相互作用の効果，日本地震学会秋季大会講演予稿集，B003, 2003．
Yamashita,T., Regularity and complexity of aftershock
occurrence due to mechanical interactions between fault slip and fluid flow,
Geophys. J. Int., 152, 20-33, 2003.