学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

11-NA13巨大地震発生サイクルシミュレーションの高度化
平原和朗（京都大学大学院理学研究科）
2011 年 3 月 11 日東北地方太平洋沖地震が発生し、未曾有の東日本大震災を引き起こした。宮城沖では M7 クラスの地震が 30 年程度の間隔で繰り返し発生していたが、その海溝側で 50m に達するすべりを生じ巨大津波を発生させ、幅 200km 長さ 500km にわたる広い領域を破壊し、M9.0 の超巨大地震となった。仙台平野の津波堆積物からこういった超巨大地震は 450-800 年間隔で発生していることが示唆されている。今世紀前半にその発生が危惧されている南海トラフの巨大地震などの巨大地震の発生予測に向けた
巨大地震発生サイクルシミュレーションはますます大規模・マルチスケールシミュレーションを必要とするため、昨年に引き続き、階層型行列（H-matrix）法によるコードの省メモリ化・高速化を進める予定であった。ところが、この超巨大地震の発生を受け、現在バージョンの H-matrix 法を用いてこの巨大地震発生サイクルモデルの構築のためシミュレーションを緊急に行なった。沈み込む太平洋プレートの３次元形状を用い、これまでに提唱されているモデルの中で、階層型アスペリティモデルと強パッチアスペリティモデルを検討した。その結果、宮城沖および茨城沖で M7 クラスの地震が発生する
一方、海溝側では固着し続け 800 年程度に一度超巨大地震を発生するという M9 超巨大地震発生サイクルの描像が得られた。

1. 研究の目的と意義
東北日本では東から太平洋プレートが日本海溝

HKNK沿いに、西南日本では南からフィリピン海プレートが南海トラフ沿いに沈み込んでいる。これらのプレートの沈み込みにより、およそ百年間隔でプTH

レート間巨大地震が発生し、人的・物的被害を含NTむ大災害が生じている（図１）。特に、南海トラフ
では今世紀前半に巨大地震の発生が危惧され、その被害総額は国家の財政規模に相当すると試算されている。また、1995 年神戸地震に見られるように、南海トラフ巨大地震発生の前 50 から後 10 年の間に西南日本内陸活断層での地震活動が高まることも報告され、南海トラフ巨大地震による被害に加え、内陸地震による大被害も危惧されている。

図１

プレート収束帯としての日本列島
PAC:太平洋プレート、こういった状況の中、最近、プレートの沈み込

PHS:フィリピン海プレートみによって引き起こされるプレート間巨大地震の

AM：アムールプレート

発生予測を目指して、岩石実験から導かれた摩擦

NA・北米プレート構成則に基づき、プレート運動を原動力とした、カッコ内は収束速度を表す。

プレート間巨大地震発生サイクルシミュレーショ

TH:東北地方太平洋沖地震震源域ンが行われ、東北日本および西南日本の海域でのNT:南海トラフ巨大地震震源域巨大地震発生履歴の再現に成功するなどの成果をHK:北海道東岸 500 年間隔地震上げている。

1

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

しかしながら、現状は、プレート境界に分布させる摩擦パラメータの推定は試行錯誤的に探索している段階であり、次のステップとして、気象海洋分野で開発されたデータ同化手法の導入による地震発生予測システムの構築を考える段階にあると言える。また、特に西南日本の震源域深部で最近観測されている様々なスロースリップの活動は、巨大地震発生前に前兆的変動を示す可能性が指摘されており、巨大地震発生サイクル中での種々のスロースリップ活動変化を含むマルチスケールのプレート間地震サイクルシミュレーションを行う必要がある。こういったデータ同化やマルチスケール地震発生サイクルシミュレーションには長大な計算時間および大容量のメモリを要し、現在のシミュレーション手法では南海トラフといった広領域での現実的な計算は不可能である。このため、本研究では、より高速化・省メモリ化を実現する地震発生サイクルシミュレーションコードを開発する。
また前に述べたように、プレート間地震に加え内陸地震の発生も考慮するには、地震時における図２弾性応答のみならず、地殻下部やマントルの粘弾

陸上の GPS 観測による陸上の地殻変動データと GPS/音響測距による海底地殻変動観性による時間遅れを考慮する必要がある。さらに測データより推定された断層面上でのすべり分

プレートの沈み込みにより日本列島下には大きな布 （ 国 土 地 理 院 ・ 海 上 保 安 庁 ,2011,

不均質構造が作り出されている。現状のシミュレ

http://www.gsi.go.jp/cais/topic110520-index.hーションでは均質半無限弾性体を仮定しており、tml）粘弾性や不均質性の影響は考慮されていない。そこで、内陸地震との相互作用を含むプレート間巨大地震サイクルシミュレーションに向けて、3 次元

方宮城沖では M7〜8 クラスの巨大地震のみが想定

不均質粘弾性媒質中での地震発生サイクルシミュされて、想定外と言える状況にあった。この超巨レーションコードを開発する。

大地震の発生過程を解明する必要があり、予定を上記が申請書に記したもので、今年度も引き続変更して緊急に超大地震サイクルシミュレーショ

き、地震発生サイクルシミュレーションコードのンを行っている。こういった超巨大地震サイクル省メモリ化・高速化を進める予定であった。とこの発生機構の解明は今後の南海トラフ巨大地震のろが、2011 年 3 月 11 日に発生した東北地方太平洋予測に向けても非常に重要な課題と捉え、研究を沖地震（マグニチュード M9.0）が発生し、未曾有進めている。の東日本大震災を引き起こした。こういった数十

宮城沖では M7 クラスの地震が 30 年程度の間隔年から百年に一度のプレート境界地震とは異なる、 で繰り返し発生していたが、その海溝側で 50m に超巨大地震の発生は、南海トラフや北海道東岸の達するすべりを生じ巨大津波を発生させ（図２）、海溝沿いでその発生が予測されていたが、東北地

450x200km の領域を破壊し、M9.0 の超巨大地震と2

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月なった。仙台平野の津波堆積物からこういった超







 i t   i *   i t  ai ln V* / Vi t  V* / Vst

巨大地震は 450-800 年間隔で発生していることが示唆されている。50m を超える巨大なすべりを発生 (2)

 Vpl,i Vpl,i    t  V t  V t  di t  biVpl,iexp   i  exp   i   i  i t   bi ln dt
Libi Li  Vc  Vi t  Vst  した領域は、プレート収束速度（約 9cm/年）を考(3)えると 500 年以上にわたりぴったり固着していたことになる。福島沖では、また、福島沖では 2000を連立させて各プレート境界セルにおけるすべり年頃までは固着していたが、2005 年頃からは地震の発展をシミュレートする。実際には、式(1)と(2）

時すべりより余効すべりの方が大きい地震が発生を時間で微分し連立させ、速度の時間微分を得て、し、また海底 GPS 観測からも固着のはがれが進行

こ れ を 式 (4) と 連 立 さ せて 、 時 間 刻み 幅 可 変 のしていたことが知られていた。これらの特徴を説Runge-Kutta 法により積分している。式(2)のよう明する地震サイクルモデルを構築する必要があるに摩擦がすべり速度 と状態変数 （その時間発展は式(3)に従う）との関数となっており、σは法線

2. 当拠点公募型共同研究として実施した意義

応力で、式(2)の右辺第 2 項は摩擦係数となってい本研究は、京都大学理学研究科・海洋開発研究る。速度と状態に依存する摩擦構成則と言われる機構・東京大学情報基盤センター・京都大学学術所以である。また、ここで、 a および b(または情報メディアセンターからなる、超大規数値計算

A=σa および B=σb)は、各セル i に与える摩擦特系応用分野における、学際大規模情報基盤共同利性をコントロールするパラメータで、２つのパラ用・共同研究拠点公募型共同研究である。

メータの差（a-b または A-B）が正（負）だと速度強化（速度弱化）を示し、すべり速度が増すと3. 研究成果の詳細と当初計画の達成状況

摩擦が増える（減少する）。また、L は特徴的すべ(1) 研究成果の詳細り距離と呼ばれ、大きいほど破壊エネルギーが大以下に示すように、岩石実験から得られたすべきく破壊しにくくなる性質を持っている。これら

り速度と状態に依存する摩擦構成則に基づき、プの摩擦パラメータをプレート境界上に分布させ複レート運動を駆動力とする地震発生サイクルシミ雑な地震発生サイクルを再現する。

ュレーションを行なった。

プレート境界を小さなセルに分割すると書いた

半無限均質弾性体を仮定し、その中にプレート境界が、この分割の影響が地震サイクル計算結果に影を設定する。プレート境界を N 個の小さなセル（小響を及ぼさないようにするには、分割サイズ h が、h < h* = cGL/σ(b-a)

断層）に区切ると、各セル i での準静的運動方程式にとる必要がある(Rice,1993)（ここで、c は分割は次のように書ける。

 するセルの形に依存する定数）。通常用いられてい 2 V t  (1) i t    K ij  j t  V pl , j t j(4)Gるパラメータ範囲では、セルサイズは 1km 以下とiなる。
式(2)と(3)から分かるように、問題の非線形性に

ここで、  i , Vi ,  j , V pl, j , K ij , G,  は、それぞよる繰り返し計算の必要性に加えて、分割セルサれ、セル i でのせん断応力、すべり速度、セルイズの要請から来る、セル数の多さにより、地震ｊでのすべり、プレート収束速度、セル j での発生サイクルシミュレーションは、多大な計算時単位すべりによるセル i での応力変化（すべり間とメモリを要する。式(1)に見られる、すべり応応答関数）、剛性率、S 波速度である。これと、答関数とすべりの行列・ベクトル積

 K ij  j t  V pl , jt 

速度と状態に依存する摩擦則のうち、compositelaw (Kato and Tullis,2006)

j3(5)学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月の部分（またはその時間微分したもの）の省メモリ化・高速化を開発中である。分割したセルの数を N とすると、これには、O(N2)のメモリおよび演算回数を要する。この行列・ベクトル積の演算に H-matrices 法（Börm et al., 2009）を用いて、O(N)〜O(NlogN)となる省メモリ・高速化を開発中である(平原・他,2011a; Ohtani et al., 2011a)。
開発中の H-matrices 法を用いてシミュレーションを行なった(大谷・他,2011；Ohtani et al.,2011b,Ohtani, 2012)。
3-1 モデル
３次元形状を持つプレート境界を 200,704 個のすべり応答関数の H-matrix 表

図３

三角形のセル（小断層）に分割する（N＝200,7004）。

示。もとの行列サイズは（N２、図３にすべり応答関数行列の H-matrices を示す。

N=200,7004）。小行列内の数字は近なお、プレート収束速度 Vpl は 8.4cm/年に設定し似するランク数を表す。ている。また、摩擦分布モデルとしては、階層アスペリティ（HA）モデルと通常のアスペリティ（SA）領域サイズが大きなって、臨界サイズと同程度かモデルを考え、比較した。

大きい A-B<0 の領域は条件付き不安定すべりといHori and Miyazaki(2011)を参考に HA モデル、および Kato and Yoshida(2011)を参考に通常の SAった摩擦特性を持ち、余効すべりと地震すべりを起こすことができる。

モデルを構築した（図４）。彼らのモデルは概念モ

階層アスペリティ（HA）モデルは、M９の震源域デルや２次元モデルであったが、本研究では実際全域のサイズが臨界サイズと同程度になるように的な３次元プレート形状に基づく大規模シミュレ

A-B<0 及び大きな L を設定している。その中で今ーションを実行した。

回 50m 以上のすべり生じた宮城沖浅部の領域でA-B<0 の領域では速度弱化を生じると述べたが、B-A を大きくして地震時応力降下量を大きくし、これは領域サイズが L/(B-A) に比例する臨界サイ

地震の繰り返し間隔が 500 年以上になり、大きなズ Lc より小さい時に生じ、固着―すべりを繰り返すべりが生じるように設定している。さらに数十す。こういった領域をアスペリティと呼んでいる。

年の繰り返し間隔を持つ宮城沖深部、茨城沖で L-1

図４

103

A-B (MPa)log(L/Lo)

L0=6cm0-1

1

A-B (MPa)

HA（左）と SA（右）モデルにおける A-B および L の分布
40log(L/Lo)

2lo=6cm

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月る。

図５階層アスペリティ（HA）モデル（左）と通常のアスペリティ（SA）モデル（右）の場合の各領域におけるすべり蓄積の巨大地震サイクル中での時間変化。
宮城県沖浅部（MYG_S)、宮城県沖深部(MYG_D1,D",D3)、茨城県沖（IBK）
Vpl:太平洋プレート収束速度を小さくして臨界サイズより小さなアスペリティ

HA モデルでは、地震後での余効すべりが終わるを置いている。福島沖深部は繰り返し間隔が明らと以降は M9 震源域全域にわたり強く固着するが、かでないがアスペリティを設定している。これに深部から剥がれていき、アスペリティ部分に応力対し、強パッチを持つ通常の SA モデルでは、L はがかかるようになると、そいて茨城沖（IBK）や宮浅部で大きく深部で小さいといった深さ依存で設城沖深部のアスペリティ（MYG_D1,D2,D3）では地定し、A-B>0 の速度強化領域に、絶対値の大きな震が発生するようになる。宮城沖深部や茨城沖でA-B<０の値を持つ強パッチを宮城沖浅部のすべりは地震は発生するが、蓄積したすべりを全て解放の大きかった領域に設定し、宮城沖深部、福島沖する訳ではなく、MYG_D3 ではサイスミックカップ

深部のアスペリティの設定は HA モデルと同じ考リングは 0.44 となり、従来言われてきた 0.3 に近えで設定しているが、茨城沖のアスペリティはまく、観測が再現されている。実際にも２日前に前だ設定していない。ただ、福島沖浅部では地震時震が発生し、その地震時および余効すべりが本震前のゆっくりすべりと固着を再現するようややをトリガーしたと言われているが、シミュレーシ

B-A<0 に設定している。

ョンでも MYG_D1 の地震後余効すべりで巨大地震はトリガーされ、62 日後に発生している。

3-2 シミュレーション結果

SA モデルでは、巨大地震後深部の A-B>0 領域で図５に各領域における積算すべり量の超巨大地は余効すべりが発生し、そのため宮城沖深部震サイクル内の時間変化を示す。このパラメータ(MYG_D2 や_D3)では地震活動が促進されている。

設定では、HA モデルと SA モデルでは、それぞれモデルパラメータ設定によりこの特徴を薄めるこ

巨大地震の繰り返し間隔は、825 年、853 年で、地とは可能であるが、これは HA モデルと大きく異な震時の最大すべり量は、67.8m、72,7m となっている特徴である。
5

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

図６

宮城県沖と福島県沖浅部（MYG_S 化地震時すべりによると FUKS_S）における積算すべり量のサ

イクル中での時間発展の比較。

陸域 GPS 観測から 2000 年までは福島県沖に強い化に加えて、地震時すべりによる摩擦発熱のため、固着域が観測されていたが、2007 年以降本震前ま断層間隙に含まれる流体が熱せられ間隙流体庄がでこの固着が弱くなり、また 2005 年以降に発生し上昇し（TP: Thermal Pressurization）、有効法線

たこの地域の地震の余効すべりが本震のモーメン

応力、すなわち摩擦力が急減し、特大すべりや大トを上回るといった現象が生じていえた。このよきな応力降下をもたらしたとする考えが出されてうに震源域におけるプレート間固着の剥がれ（弱いる（Mitsui et al., 2012）。こういったモデルまり）は M9 超巨大地震の準備過程を表すとも考えにおいても、強パッチ以外の M９震源域の振る舞られる。

いはここで述べた２つのモデルの振る舞いにようやくできると思われる。

HA モデルでは M9 震源域全体にすべり速度弱化ではるが条件付き安定すべり特性を与えている、(2) 当初計画の達成状況について

そこで福島県沖浅部(FKS_S)では、図６左図に示すように、地震サイクル前半では固着しているが、2011 年東北地方太平洋沖地震の発生を受け、急深部からの剥がれが浅部へ達するサイクル後半に遽その発生要因を探るべく、現在開発中の階層型入ると、ゆっくりすべりー固着をくりかえすよう行列（H-matrices)法を用いて、大規模 M9 超巨大になる。こういったゆっくりすべりの部分が 2000

地震サイクルシミュレーションを緊急に実施した。

年後半のプレートの剥がれに対応している可能性

概念的あるいは２次元モデルで提唱されていた２がある。SA モデルでも福島沖浅部に条件付き安定つのモデルについて、沈み込む太平洋プレートのすべり特性を与えると、図６右図に見られるよう３次元形状を考慮して、実際的なシミュレーショに、ゆっくりすべりー固着を繰り返す。

ンを実行し、M９超巨大地震発生の２つのシナリオをえることができた。

このように、本震前の福島県沖での固着の弱まりは本震破壊の準備過程ではなく、やはり本震は２つのモデルの違いは、地震後に現れる。すな

宮城県沖深部(MYG_D1)での地震発生およびその余

わち HA モデルでは、M９震源域で固着が進み、地効すべりによりトリガーされた可能性が高い。

震サイクルの前半中は、M７地震の発生が抑えられなお、宮城県沖浅部での特大すべり域を生成する。これに対し、SA モデルでは、M9 震源域内のアるメカニズムには、ここで考えた強い固着を持つスペリティ以外の領域では、安定すべり特性を与強パッチを設定する代わりに、高速すべり域で摩えているので、地震後には余効すべりと安定すべ擦力が急減しかつ大きな L を持つ新たな摩擦構成

りにより M７震源域に応力が掛かり、M７地震の発則（Shibazaki et al., 2011）や、摩擦係数の変生が早まりまた発生間隔も短い。やがて余効すべ6

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

りも落ち着き数錠年間隔で M7 クラスの地震が発可能となる。現状では、粘弾性体中のすべり応答生するようになる。このように、M７地震の活動が関数の計算およびメモリ変数を用いる際に階層型大きく違ってくるので、今後の観測が重要となる。

行列法を用いて単純に高速化できないという難点

このように、ある程度のモデルを示せたと言え、がありクリアすべき問題があるが、なんとか解決

M9 超巨大地震の発生を受け急遽変更された当初のして粘弾性計算に進む必要がある。

目的は達成されたと言える。

しかしながら、階層型行列（H-matrices）法を用いた地震サイクル計算の省メモリ化・高速化に

4. 今後の展望おける並列化性能を含む、これらの難点を克服し、階層型行列（H-matrices）法を用いた地震発生

内陸断層とプレート境界地震とがカップルする地サイクルシミュレーションで、並列化性能に問題震発生サイクルシミュレーションへと研究を進めがあることが分かっており、その改善に向けて方る必要がある。

策を取る必要がある。

5. 研究成果リスト

2011 年東北地方太平沖地震については、まずは現行の２つのモデルの詳細を更に調べてより定量的に色々な観測量と照らし合わせて、超巨大地震(1) 学術論文（投稿中のものは「投稿中」と明記）サイクルモデルを構築し、今後の房総沖での地震

Ohtani, M., K. Hirahara, Y. Takahashi, T. Hori, M.
Hyodo, H. Nakashima and T. Iwashita, 「 Fast

発生可能性を検討する必要がある。
更にモデルスペースを抑えるために今回破壊しcomputation of quasi-dynamic earthquake cycle

なかった岩手北・青森沖のアスペリティについて

simulation with Hierarchical Matrices」 Procedia

はモデルに入れていない。H-matrices 法の並列化

Computer Science, 4, 1456–1465, 2011a.
Hori, T. and S. Miyazaki, 「A possible mechanism of

効率の更なる向上を図り更なる大規模シミュレー

M 9 earthquake generation cycles in the area ofションに挑む必要がある。
またここでは弾性媒質での準動的シミュレーシ

repeating M 7 ∼ 8 earthquakes surrounded by

ョンを考えているが、最初に述べたように海溝型

aseismic sliding」, Earth Planets Space, 63, 773,777,

地震と内陸地震との相互作用を考える際に必要と2011.

なることに加えて、500-1000 年と言った長い繰り

Mitsui, Y., N. Kato, Y. Fukahata and K. Hirahara

返し間隔を持つ超巨大地震サイクルシミュレーシ

Megaquake cycle at the Tohoku subduction zone

ョンにはマントルの粘弾性を考慮する必要がある。

with thermal fluid pressurization near the surface,

粘弾性媒質における応力計算には粘弾性すべり応

Earth

答関数とすべり速度の履歴積分（畳み込み積分）

doi:10.1016/j.epsl.2012.01.026, 2012.

PlanetSci.

Lett.,

325-326,

21-26が必要である。そのため、過去の全てのすべり速(2) 国際会議プロシーディングス

度履歴をメモリにおく必要があり大変な計算にななし。

る（FFT を用いた計算も同様）。非弾性媒質中を伝播する地震波の減衰を差分法を用いて時間領域で計算するのに開発された、メモリ変数を用いた有(2) 国際会議発表

効な計算手法を開発中である（平原・他、2011b；

Hirahara, K., M. Ohtani, M. Hyodo, T. Hori. [FastHirahara et al.,2011b）。この方法では、粘弾性

computational methods for large- and multiscale

すべり応答関数を数個の緩和関数で近似し、各緩interplate earthquake cycle simulations,]

和関数に対応する１階の常微分方程式に従うメモIUGG2011, S10S2-3078, Melbourne, 2011.

リ変数を導入して、履歴積分を行わず応力計算が7

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

Hirahara, K., M. Ohtani, Y.

Shikakura,(5) 謝辞[Quasi-dynamic earthquake cycle simulation in a本研究は文部科学省による「東海・東南海・南viscoelastic medium with memory variables,] AGU

海地震の連動性評価研究」および「地震及び火山2011 Fall Meeting、T43I-01, San Francisco, 2011.

噴火予知のための観測研究計画」の支援を受け行Ohtani M., K., Hirahara, T., Hori, and M., Hyodo,

われました。

「Comparative study on quasi-dynamic earthquake
cycle models」for the 2011 giant Tohoku Earthquake,
AGU 2011Fall Meeting, U51B-0044, San Francisco,
2011.
(4) 国内会議発表
鹿倉洋介・深畑幸俊・平原和朗、「南海トラフ沈み込みイベントに伴う西南日本内陸活断層の破壊関数変化」、日本地球惑星科学連合 2011 年大会、SSS035-P01, 幕張、2011a.
平原和朗・大谷真紀子・兵藤守・堀高峰,「大規模マルチスケール地震発生サイクルシミュレーシ
ョンの高速化手法」、日本地球惑星科学連合 2011年大会、SSS035-40, 幕張、2011a.
堀高峰・兵藤守・宮崎真一、「東南海地震後から南海地震までの地殻変動データに対する逐次デー
タ同化の数値実験」日本地球惑星科学連合 2011年大会、SSS035-39, 幕張、2011
鹿倉洋介・深畑幸俊・平原和朗、「フィリピン海プレートと太平洋プレートの沈み込みイベントに伴う西南日本内陸活断層のクーロン破壊関数の変化」、日本地震学会 2011 年秋季大会、P3-45
静岡、2011b
平原和朗・大谷真紀子・鹿倉洋介、「メモリ変数を用いた粘弾性媒質における地震発生サイクルシミュレーション」、日本地震学会 2011 年秋季大会、D21-10, 静岡、2011b.
大谷真紀子・平原和朗・堀高峰・兵藤守、「2011 年東北地方太平洋沖地震発生サイクルモデルの準動的シミュレーションによる比較検証」、日本地震学会 2011 年秋季大会、P1-13, 静岡、2011.
兵藤 守・堀高峰、「東北地方太平洋沖地震に伴う海溝軸付近での地殻変動の解釈」日本地震学会 2011 年秋季大会、P2‒12, 静岡、2011.
8