[共同研究成果]
沈み込みプレート境界における余効すべり伝播速度の空間分布
有吉 慶介，松澤 暢，長谷川 昭
東北大学大学院理学研究科地震噴火予知・研究観測センター

本研究では，地震発生後にすべりがどのように震源域（アスペリティ）から周辺へと伝播していくのかを数値シミュレーションによって再現し，伝播過程の不均質性の主な原因がプレートの重みによる実効法線応力によるものだということを示した．
１．はじめに地震の相互作用・トリガー作用の主な要因として，地震時に高速にすべったアスペリティの周囲で余効すべりが生じ，それが S 波速度よりもはるかに低速で伝播して他のアスペリティの破壊を促進する，ということが挙げられる．一例として，三陸はるか沖地震における余効すべりの伝播速度は陸域に近いプレート境界面深部では GPS 観
測によって 10km/month のオーダーを示し，海溝近くのプレート境界面深部では相似
地震解析によって 10km/day のオーダーと見積もられている[1]．このような伝播速度の違いが何で決まるのかを調べることは，今後，アスペリティの相互作用を理解し，プレート境界面の摩擦パラメターを推定する上で重要だと考えられる．
２．数値シミュレーションの計算手法
基本的には，均質半無限弾性体中にあるプレート境界面にすべり速度/状態依存摩擦構成則を用いた数値シミュレーション[2]を３次元に拡張したものなので，以下その概略のみ述べる．
本稿では，M7 クラスの地震が繰り返し発生する一つのアスペリティを含む３次元モデルを想定した（図１）．すべりの成分は逆断層型の dip slip のみを考慮し，地表を除く対象領域外側では，定常的プレート収束速度 V pl で一様にすべっているものとする．

図 1

３次元プレート境界面における計算対象領域．

計算領域内のプレート境界面では，dip slip により生じるせん断応力と摩擦力が準静的につり合っていると仮定し，摩擦係数は速度/状態依存摩擦構成則の slownesslaw[3]に従うものとする．これらの条件のもとに，プレート境界面上での力の釣り合いは以下のように記述することができる．
Nτ idislocation = ∑ K ij (u j (t ) − V pl t ) −j =1G dui
2 β dt(1)

σ i ( z ) = ( ρ rock − ρ w ) gz

(2)

τ i force = µiσ i(3)µi = µ0 + ai log(Vi / V0 ) + bi log(Viθi / d c )(4)dθ i / dt = 1 − Viθi / d ci

(5)

τ idislocation = τ i force(6)i添え字 i, j は分割セルの指標，(1)式の K ij は半無限均質媒質での弾性論から求められたグリーン関数[4]，u は変位，t は時間，式の右辺第２項は地震波放射による減衰項[5]，σ は有効法線応力，(6)式は食い違いによるせん断応力 τdislocationと摩擦力 τforceの準静的なつり合い条件をそれぞれ示す．応力計算については，(1)式に示すように対象領域内での変位量と定常的なすべり量 V pl t との差分のみとし，定常的なすべりによる応力変化成分は高速すべりによる応力変化に比べて十分小さいものとして無視する[6]．(4)(5)式は速度/状態依存摩擦構成則を表わし， a, b, d c は，摩擦特性を示すパラメターを示し[7][8]， ( a − b) < 0 ではすべり速度と共に摩擦係数が低下し，不安定すべりを生じることが知られている．逆に ( a − b) > 0 では，粘性のようにすべり速度の増加と共に摩擦が増加する性質をもつ．(1) -(6)から導出される微分方程式を，Runge-Kutta 法を用いて数値的に解いていく．
定数については次のように与えた．剛性率 G =30GPa．S 波速度 β =3.75km/s，岩石の密度 ρ rock =2.7g/cm ，水の密度 ρ w =1.0g/cm ，重力加速度 g =9.8m/s ，基準速度332V0 = 1µ m/s，基準速度時の摩擦係数 µ0 =0.6，ポアッソン比 0.25．
本研究で用いる数値モデルでは，特に地震が発生する領域内では，極めて高速なすべりが生じるので，メッシュサイズを十分小さく設定する必要がある．そのため，メッシュ数は数万〜十万程度と膨大になり，計算負荷が極めて大きくなるが，大型計算機センターとの共同研究によって，当初のものと比べてバンクコンフリクトが 500 分の 1 にまで減少し，ベクトル化率も 98.8%から 99.4%へと効率化することで，実質的計算速度を３倍ほど上昇させることができた．

３． 摩擦パラメターの空間分布モデル
摩擦パラメター ( a, b, d c ) の空間分布については，できるだけ不均質性による影響を排除するために， d c は 4cm として全空間で一様とし， ( a − b) については，ここでは①アスペリティ領域，②安定すべり領域の２種類の領域に分け，境界ではステップ的に変化するとした．図１において，メッシュ間隔が 0.5km 以下となっているところが①に相当し，他は②の領域となる．①では ( a − b) = −0.00072 < 0 として，摩擦力の不安定性によって地震性の高速すべりが起こる領域，②では ( a − b) = 0.0001 > 0として，安定化させる摩擦特性によって余効すべりが生じ，遠くに伝播するに従って減衰していく安定的な領域とした．

４．シミュレーション結果
１．で述べたようなプレート境界面の浅部と深部での余効すべり伝播速度の違いは，摩擦パラメターが等しいことから，(2)式で示されている有効法線応力の深さ依存によるものだと考えられる．そこで，基準モデル（以下，モデル A と呼ぶ）として(2)式を使い，比較対象のモデル（以下，モデル B と呼ぶ）として，有効法線応力の空間分布を，モデル A におけるアスペリティの中心位置での値σ=480MPa として，深さに関わらず一定値を与えた．両者を比較した結果を図２に示す．図２(a)(c),(b)(d)は，それぞれモデル A と B の地震発生から約 0.1 年後と 3 年後におけるすべり速度の空間分布を表わし，定常的なプレートの収束速度で規格化された常用対数の値を色で示している．両者の図を見比べてわかるように，図２(a)(c)では浅部で速く，深部で遅く余効すべりが伝わっているのに対し，図２(b)(d)ではほぼ同心円状に伝播していく様子がわかる．これらの結果より，実効法線応力が余効すべりの伝播過程に大きな影響を及ぼすことが示された．

図 2(a)(b)モデル A,B それぞれでの地震発生から約 0.1 年後におけるすべり速度の空間分布．(c)(d)は同様に約 3 年後のもの．どちらもアスペリティ領域では地震発生後の強度回復（固着）過程に入っている．アスペリティの周囲にある暖色系の領域が余効すべり過程に相当する．

ここで実例として，1994 年 12 月に起きた M7.6 の三陸はるか沖地震（図３）とモデル A における余効すべりの伝播過程（図４）を見比べると，図３の A,B,C 領域では，地震発生後 100 日ほどで余効すべりが顕著であるのに対し，図３の D,E,F 領域
では地震発生後から数日以内に顕著な余効すべりが生じていることがわかる．これらの傾向を図４では，時間スケールは数倍程度長くなっているが，よく再現されていることがわかる．

図 3 相似地震解析によって求められた 1994 年 12 月 28 日に起きた三陸はるか沖地震（M7.6）前後での積算すべり量の時系列．A,B,C の領域は陸域に近く，プレート境界面は深さ 40-60km にあるのに対し，D,E,F は海溝に近く，数 km の深さにあると推定される．

図 4 本研究によって再現された積算すべり量の時系列．図１の中央（X=0）地表面付近からプレート境界面に沿って 20km 置きに並べたものであり，上３つが深部，下３つが浅部安定すべり領域，中３つがアスペリティ領域に相当する．

５．まとめ
図３と図４を比べると，浅部と深部で余効すべりの伝播速度が異なる様子がよく再現され，その原因が図２で示されたように，実行法線応力の影響によるものだということがわかった．これらの結果から，実際の観測で見積もられている余効すべりの伝播速度を定量的に再現するためには，どのような条件を与えればよいのかを考えると，実際の実効法線応力が静岩圧で見積もられているものよりも数分の１〜数十分の１程度小さい可能性が挙げられる．これは，間隙水圧が静岩圧に近い値になっていることを示唆するものである．
今後は実効法線応力だけでなく，摩擦パラメター a, b, d c や地震の規模などが及ぼす余効すべり伝播過程への影響なども調べていく予定である．

謝辞
本研究は東北大学情報シナジーセンターとの共同研究として行われました．スーパーコンピューターSX-7 上でのプログラミング最適化や利用環境などについては皆様より多大なご協力を頂きました．ここに記して感謝致します．

参考文献[1] Matsuzawa, T., N. Uchida, T. Igarashi, T. Okada, and A. Hasegawa,
Repeating earthquakes and quasi-static slip on the plate boundary east off
northern Honshu, Japan, Earth Planets Space, 56, 803–811, 2004.
[2] Kato, N. and Hirasawa, T. A Numerical study on seismic coupling along
subduction zones using laboratory-derived friction law, Phys. Earth Planet.
Inter. 102, 51-68, 1997.
[3] Beeler, N. M., Tullis, T. E., and Weeks, J. D., (1994): The roles of time anddisplacement in the evolution effect in rock friction, Geophys. Res. Lett. 21,1987-1990.
[4] Okada, Y., Internal deformation due to shear and tensile faults in a halfspace,Bull. Seism. Soc. Am., 82, 1018–1040, 1992.
[5] Rice, J. R., Spatio-temporal complexity of slip on a fault, J. Geophys. Res., 98,9885– 9907, 1993.
[6] Savage, J. C., A dislocation model of strain accumulation and release at asubduction zone, J.Geophys. Res., 88, 4984-4996, 1983.
[7] Dieterich, J. H., Modeling of rock friction, 1, Experimental results andconstitutive equations, J. Geophys. Res., 84, 2161–2168, 1979.
[8] Ruina, A. L., Slip instability and state variable friction laws, J. Geophys. Res.,88, 10359-10370, 1983.