— 3 —SENAC Vol. 50, No. 2（2017. 4）［共同研究成果］

数値シミュレーションから見出された
余効すべり伝播速度と摩擦特性の関係- Relationship between afterslip propagation speed and frictional properties 有吉 慶介*1・松澤 暢*2・Roland Burgmann*3・
長谷川 昭*2・日野 亮太*2・堀 高峰*1
*1

独立行政法人 海洋研究開発機構 地震津波・防災研究プロジェクト*2

東北大学大学院理学研究科附属 地震・噴火予知研究観測センター*3

カリフォルニア大学バークレー校
巨大地震発生後に付随する余効すべりの伝播過程を再現するために，これまでの数値シミュレーションでは入力パラメタ―を試行錯誤的に決めていた．本研究では，多数の数値シミュレーション結果から共通の特徴を見出したことで，摩擦パラメタ―との関係式を解析的に導出することに成功した．これにより，本震から余震が発生するまでの時間差を事前に見積もったり，逆に，観測された時間差から摩擦特性を絞り込めることが期待される．

１．はじめに巨大地震が発生して断層にずれが生じた後，その歪を解消するために，その周辺にすべりが広がる現象は，「余効すべり」[1]と考えられている．この「余効すべり」は，地震波伝播速度よりもはるかに遅く，前震→本震や本震→余震の時間差が生じる原因の一つだと考えられている．
例えば，近い将来に巨大地震が発生すると考えられている南海トラフ沿い (図１) では，1707 年の宝永地震では全体が一気に壊れ，1854 年の安政地震では東海→南海で 30 時間の時間差が生じ，1944/1946 年の昭和地震では東南海→南海で 2 年間の時間差が生じているが，このように地震が連動する時間差が毎回異なる理由については未解明のままである．

図１．南海トラフで過去に起きた主要な地震 [2]．

— 4 —

SENAC Vol. 50, No. 2（2017. 4）

一方で，これまでの数値シミュレーション研究では，東北地方太平洋沖地震をはじめとする海溝型巨大地震の発生サイクルを再現すべく，その時間差が実際に観測されたものに近づけるために，入力値となる摩擦パラメタ―を試行錯誤的に多数のモデルを導入してきた [3]．そのため，実際に地震が発生してからモデルが構築されるまでに何年もの時間を要し，即応的な地震発生予測をするのが難しいという問題を抱えている．
そこで本研究では，余効すべりと摩擦特性との関係について，単純化したモデルの数値シミュレーション結果に基づいて近似条件を見出すことで，両者の関係を解析的に導出し，その妥当性を検証すると共に，過去の地震活動に基づいて摩擦特性を考察した．

２．数値シミュレーションに基づく余効すべり伝播速度と摩擦パラメタ―との関係
余効すべりの伝播速度は，摩擦特性だけでなく断層の形状や弾性定数などにも影響を受ける．そこで，他の要因を排除するために，地震サイクル数値シミュレーションとして単純化された数値モデル（図２）を適用した [4]．計算手法等の詳細については，過去の文献を参照されたい [5]．
数値シミュレーションから再現された余効すべり伝播過程の特徴として，図２．適用した３次元沈み込みプレート境界面のバケツリレーのようなかたちで伝播し数値モデル．各区間の指標として記号で示す．ていることが分かった．
これを言い換えると，ある物体（スライダー）がすべると，隣のスライダーとつなぐバネが縮み，その後，隣のスライダーがすべる・・・，という過程が連鎖的に続く，というものである．この近似条件を，岩石実験から得られた摩擦法則（速度・状態依存摩擦構成則）に適用すると，伝播速度 Vprop は，摩擦パラメタ― (�� �� �� )，有効法線応力 (�)，断層形状に基づくばね定数 (�)，スライダーがすべる直前までに負荷する応力擾乱振幅 (�) を関数として，
����� ��� �� �� � �� �� ������ =��� ������
������� ������ +���������と記述することが出来た．ここで，i はスライダーの指標，��� はスライダーの間隔（メッシュサイズ），L は応力擾乱の時間関数（ramp 関数，step 関数）で決まる因子，�は (b, dc,k, σ) で決まる因子，����� は (a, b, dc) で決まる因子，��� はプレート運動速度と (k, σ) で決まる因子である（詳細な記述は省略する）．この導出にあたっては，摩擦構成則の数式を展開しただけでは不十分で，大規模数値計算を共同研究開発によって多数実行することにより，共通した近似条件を見出すことが出来た点が大きく寄与している．

数値シミュレーションから見出された余効すべり伝播速度と摩擦特性の関係— 5 —この関係式で得られる伝播速度と，地震サイクル数値シミュレーションから再現された伝播速度を比較すると，自由表面効果が強い地表付近を除いて，全ての場所で一桁以内の精度で再現することが出来た (図３)．この際，東北大学サイバーサイエンスセンターとの共同研究によるプログラム開発によって多数の数値モデルで検証することになったため，妥当性を客観的に示すことが出来た．

図３．プレート収束速度(Vpl) で規格化した最大すべり速度と，本研究で適用した伝播速度の理論値 (Vana) を数値実験から求められた値 (Vsim) で規格した精度評価との関係．記号は，図１の区分に対応し，青色は Lfunc における ramp 関数，緑色は step 関数をそれぞれ示す．
自由表面の効果が残る△以外は，全て一桁以内の範囲に収まっている．

３．有効法線応力と余効すべり伝播速度との関係
本研究の数値シミュレーションから示された特徴について，ここでは一例として，余効すべり伝播速度と有効法線応力との関係を示す．図４を見ると，応力変化量Δτが大きく，その時間率が大きい（ramp < step = ∞）ほど，有効法線応力に対する伝播速度の依存性が大きいことが分かる．また，(a),(b) と(c),(d) が異なることから，沈み込み方向と走行方向でも伝播速度に違いがあることが示された．灰色の部分は，解析解の伝播速度が負の値となり，物理的に不合理なものとなった．これは，上記で述べた，ばねスライダーモデルにおいて，ばねが縮まらないことを意味する．

— 6 —

SENAC Vol. 50, No. 2（2017. 4）

図４．有効法線応力と余効すべり伝播速度との関係．記号は，図１の区間に対応し，応力変化量Δτは，Hσに相当する．参考として，間隙圧が静水圧と仮定した場合の深さへの変換値も付記した．灰色の部分は，伝播速度が負の値となるため，プロットの対象外とした．
この結果を東北地方太平洋沖地震など，既に既存観測から捉えられたデータに適用することが出来れば，地震発生メカニズムの過程を知る上で重要な情報が得られると期待される．JFAST による東北沖の掘削では，断層が 50m も滑った付近での見かけ摩擦係数は 0.08程度と推定された [6]．これは，有効法線応力が極めて小さく，図４の灰色部分に相当することになる．これをばねスライダーモデルに立ち返って考察すると，ばねが縮まずに棒のような状態のままスライダーと一緒にすべったことを意味し，時間遅れを伴う余効すべり

— 7 —

数値シミュレーションから見出された余効すべり伝播速度と摩擦特性の関係が起こらない環境であったことを示唆するものである．
他にも最近では，スロー地震が南海トラフの走行方向に沿って移動する現象が地震計の観測網から捉えられているが [7]，これが小地震の余効すべりに伴う連鎖破壊として説明出来ると考えると，見かけ摩擦係数や応力降下量を推定することが出来る可能性がある．同様に，2011 年 3 月 11 日の 2 日前に発生した M7 の前震についても，前震の余効すべりが本
震 M9 の震源域まで伝播したとすれば，その間の摩擦特性を推定することが期待される．同様のことは，伝播方向を推定した釜石沖地震 [1] についても当てはまる．今後はこのような現象について，現在整備中の海底観測網 (DONET, S-net) の観測データと併せて検証することで，今後発生する巨大地震が連動する時間差について事前に定量的な評価を示すことを目指す．

謝辞
本研究は，東北大学サイバーサイエンスセンターとの共同研究の一環であり，研究費の一部は科学研究費補助金として，新学術領域（研究領域提案型）「余効すべり過程の理論解
析と数値モデリングに基づく摩擦特性の推定 (15H01148)」および「低速変形から高速すべりまでの地球科学的モデル構築 (16H06477)」から支援を受けたものである．

参考文献[1] 有吉 慶介・内田 直希・松澤 暢・日野 亮太・長谷川 昭・堀 高峰・金田 義行, 2015. 東北地震に伴う固有地震活動の揺らぎから推定された摩擦特性と余効すべり伝播過程,
SENAC 48(1), 33-38.
[2] 地震調査研究推進本部 地震調査委員会, 2013, 南海トラフの地震活動の長期評価（第二版）, pp. 96.
[3] R. Nakata, T. Hori, M. Hyodo, K. Ariyoshi, 2016, Possible scenarios for occurrence of M～7interplate earthquakes prior to and following the 2011 Tohoku-Oki earthquake based on numericalsimulation, Scientific Reports, vol. 6, No. 25704, doi:10.1038/srep25704.
[4] K. Ariyoshi, T. Matsuzawa, A. Hasegawa 2007, The key frictional parameters controlling spatial
variations in the rate of postseismic slip propagation on a subduction plate boundary, Earth andPlanetary Science Letters, vol. 256, 136-146, doi:10.1016/j.epsl.2007.01.019[5] 有吉 慶介・松澤 暢・長谷川 昭, 2005. 沈み込みプレート境界における余効すべり伝播
速度の空間分布, SENAC 38(4), 3-11.
[6] K. Ujiie, et al. 2013, Low Coseismic Shear Stress on the Tohoku-Oki Megathrust Determined
fromLaboratory

Experiments,

Science,

vol.

342,

No.

6153,

1211-1214,

doi:10.1126/science.1243485[7] 有吉 慶介・松澤 暢・日野 亮太・長谷川 昭・堀 高峰・中田 令子・金田 義行, 2014. 超低周波地震の活動に基づくプレート間固着の推定, SENAC 47(1), 29-34.