[共同研究成果]
相似地震震源域における非相似地震発生の可能性
有吉 慶介，松澤 暢，長谷川 昭
東北大学大学院理学研究科地震・噴火予知研究観測センター

相似地震震源域であっても，場合によっては通常とは異なる地震（非相似地震）が発生するであろうと一般に考えられているが，そのメカニズムについては不明であった．本研究ではそのメカニズムの有力な原因の一つとして，他の震源域で発生した大規模な余効すべりの通過に伴うすべり速度場の揺らぎであることを，数値シミュレーションから示すことに成功した．これにより，相似地震解析に基づいて積算すべり量の時系列を定量的に推定する上で，大きな指針を与えると考えられる．
１．はじめに日本の沈み込みプレート境界に着目すると，大規模な地震[1]から小規模な地震[4]までの広いマグニチュード範囲で，同じような規模の地震が同じ場所で繰り返し起こっていることが，近年明らかになってきた．このことは，（断層が一気にずれるために地震を引き起こす）アスペリティ領域と（断層がゆっくりすべるため地震が起こらない）非アスペリティ領域とが棲み分けるような強度分布[10]が，図１に示すように少なくとも第一近似的には時空間的に大規模〜小規模で存在していることを意味する．特に，地震の規模に比べてその震源域近傍での固着期間内で生じる速度・応力場の揺らぎが小さい場合には，地震波形で表わされる破壊様式が似ている（相似地震）だけでなく，発生間隔まで規則的になることも，観測から見出されてきている[4]．また，地震の規模と地震一回当たりのすべり量は比例関係にあることが知られており[4]，相似地震の発生頻度からその震源域周辺でのすべり速度場をモニタリングする手法も確立されつつある[4][6]．
これらの現象を別の観点からみると，個々のアスペリティは，地震の破壊様式も発生間隔も常に同一の地震が規則的に繰り返し発生している性質を本来持っているものの，アスペリティが近接しているなどで相互作用が生じることによって，発生間隔が異なり，相互作用の影響が強い場合には破壊様式まで異なり複雑になるという考えも可能である．これに相当する可能性のある観測事例として，釜石沖の相似地震活動を挙げると，M4.8±0.1 の規模の地震が 5.52±0.62 年の間隔で発生していることが知られている[9]．この相似地震は一般に知られている地震に比べて，変動係数（相対標準偏差；標準偏差と平均値との比）が他の地震（例えば宮城県沖地震では，37.1±6.61 年）に比べて小さい．このことは，図１に示すように釜石沖の相似地震震源域周辺では，他のアスペリティから孤立しており，速度や応力の揺らぎが小さいためだと考えられている[8]．釜石沖地震の発生間隔や規模が揺らぐ具体的な原因としては，前回 1995 年の地震発生間隔が短いのは，1992 年に起きた最大 M6.9 の群発地震に伴う余効すべりの加速によるものだという考えが挙げられているが，今のところ他の手段から独立に確かめられていないのが現状である．
本稿では，シナジーセンターとの共同研究によって可能となった，極めて大規模な数値シミュレーションを用いて，1992 年に発生した三陸はるか沖地震のような大規模なアスペリティによって生じるすべり速度場・応力場の揺らぎによって，相似地震活動がどのように変動するのか定量的に考察することを目的とする．

図４−１．本章でモデル化の対象となる，釜石沖地震周辺におけるアスペリティの空間分布の概念図[8]．

２．数値シミュレーションの計算手法
基本的には均質半無限弾性体中にある３次元プレート境界面にすべり速度/状態依
存摩擦構成則を用いた数値シミュレーション[11]であり，以下その概略について述べる．

図２．本章で扱う共存アスペリティモデルの概要．（ａ）沈み込みプレート境界におけるセルの空間配置の概略図．ここでは，実際のメッシュサイズより 4 倍粗く表示している．（ｂ）プレート境界面上における摩擦パラメター γ = ( a − b) の空間分布図．
青色がアスペリティ領域，桃色が弱安定すべり領域，赤色が強安定すべり領域に相当する．

本稿では，図２(a)に示すような離散化された各計算セルにおいてすべり速度の時系列を逐次的に求めることで，プレート運動を表現することになる．すべりの成分は逆断層型すべりのみを考慮し，横ずれ成分は計算に含まないものとする．また，地表を除く対象領域外側では，定常的プレート収束速度 V pl で一様にすべっているものとする．
計算領域内のプレート境界面では，dip slip により生じるせん断応力と摩擦力が準静的につり合っていると仮定し，摩擦係数は速度/状態依存摩擦構成則[12][13]のslowness law[14]に従うものとする．これらの条件のもとに，プレート境界面上での力の釣り合いは以下のように記述することができる．
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 は分割セルの指標， N はセルの総数，(1)式の K ij は半無限均質媒質での弾性論から求められたグリーン関数[15]， u は変位， t は時間，式の右辺第２項は地震波放射による減衰項[16]， σ は有効法線応力，(2)式の z は深さ，(6)式は食い違いによるせん断応力 τ dislocation と摩擦力 τ force の準静的なつり合い条件をそれぞれ示す．応
力計算については，(1)式に示すように対象領域内での変位量と定常的なすべり量 V pl tとの差分のみとし，定常的なすべりによる応力変化成分は高速すべりによる応力変化に比べて十分小さいものとして無視する[17]．(4)(5)式は速度/状態依存摩擦構成則を表わし， a, b, d c は，摩擦特性を示すパラメターを示し[12][13]， ( a − b) < 0 ではすべり速度と共に摩擦係数が低下し，不安定すべりを生じることが知られている．逆に(a − b) > 0 では，粘性のようにすべり速度の増加と共に摩擦が増加する性質をもつ．
(1) -(6)から導出される dV/dt 及び dτ/dt に関する微分方程式について解析的に求め，Runge-Kutta 法[18]を用いて数値的に解いていく．
定数については次のように与えた．剛性率 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（等方均質弾性体を仮定）．
N

本稿では計算の高速化を図るために，(1)式中にある∑ K (u (t ) − V t ) の積和に
j =1ijjpl

ついて，畳み込み計算と同様の形式をしている[16]ことから，海溝軸方向（[ X ]軸方向）について FFT を適用することで，計算コストを約 1/30 まで大幅に軽減することに成功した．一般に，FFT の計算ではベクトル化が難しいとされているが，シナジーセンターで導入されている SX-7 ではベクトル化が可能であるため，FFT を導入してもベクトル化率が 98%と高いことが大きく関係している．このため，計算セルの総数が数十万のオーダーであっても，半日〜数日程度で計算を実行することが可能となった．
摩擦パラメターの空間分布モデルについては，図２(b)に示すように，計算対象領域を γ = ( a − b ) の値によって，①大アスペリティ領域，②小アスペリティ領域，③弱安定すべり領域，④強安定すべり領域の４つに分けることにした．
３．数値シミュレーションの計算結果から示された相似地震の揺らぎ
図３(a)に印された（▲，☆，◆，〇，■）の５点と小アスペリティ領域中心点での積算すべり履歴などを図５(b)に示す．☆は大アスペリティ領域の中心点に相当する．
大アスペリティの地震規模と発生間隔はそれぞれ M7.3 と 49.3 年のまま単一周期の地震が繰り返し発生しており，これは小アスペリティ領域を弱安定すべり領域に置き換えて単一のアスペリティとした場合でも同一であった．このことから，小アスペリティ領域から生じる揺らぎは大アスペリティでの地震活動に対してほとんど影響を及ぼさないことがわかる．〇は小アスペリティ領域中心点と対称の位置にあり，ここでの積算すべり履歴を図５(b)の赤と青の線で示すが，これらの曲線も小アスペリティの有無に関わらず同じ時系列を示すことから，これらの曲線は大アスペリティ領域が小アスペリティ領域に及ぼすすべり速度場の揺らぎに相当するものだと考えることが出来る．この二つの曲線に挟まれたものが小アスペリティ領域中心点でのすべり履歴として示しているが，両者はよく一致していることがわかる．相似地震解析に基づく

積算すべり量の推定は，第１章で述べたように，地震の規模と１回当たりのすべり量が比例関係にあり，その比例係数を San Andreas 断層間に設置したクリープメーターで観測された地殻変動量から求めている[19]．図５の結果は，この関係が沈み込みプレート境界面上のすべりにおいて相互作用が働く複雑な場合でも成り立つことを意味し，別の見地から証明することが困難であった相似地震解析に基づく積算すべり量の推定の妥当性について，速度/状態依存摩擦構成則に基づく数値シミュレーションから示したことになる．また，逆に言えば，小アスペリティ領域で発生する相似地震の発生間隔及び規模の揺らぎは，周囲のすべり速度場と対応していることを示している．従来の地震発生時刻や規模の予測には，過去の地震カタログを元に確率論的に基づいて推定したり，すべり速度場が時間変化せずに一定とした前提で見積もったりしていたが，この結果によれば，GPS 観測等からアスペリティ領域周辺でのすべり速度場を精度良く捉えることが出来れば，決定論的に推定できる可能性がある．
地震性すべりを抽出するために，小アスペリティ領域の中心点においてすべり速度が 1mm/sec を超える高速すべりの成分を地震性すべりとして V>1mm/sec に示した．
この履歴に注目すると，大アスペリティ領域で発生した地震に伴う余効すべりが小アスペリティ領域を通過する前後で，地震性すべりが発生していないことがわかる．地震性すべりを含めた総すべり量では大きなすべり量が解放されていることから，ここでは非地震性のスロースリップイベント[20]が発生していることがわかる．このことは，相似地震震源域において，余効すべりに伴うすべり速度場の揺らぎによって非相似性の地震が発生することを示唆するものである．言い換えると，相似地震解析において，非相似地震と認識されるために積算すべり量の過小評価が為される可能性があることを意味する．
この観測事例として，十勝沖地震の余効すべり解析結果を図４に示す．図４(a)に示した８つの相似地震グループ領域について，相似地震解析から推定された平均すべり量（赤線）と GPS 解析から推定されたすべり量（青線）を比較する[21]と，両者が一致しておらず，全体的に相似地震解析から推定されたすべり量について過小評価の傾向がみられる．中でも特に顕著なものとして，図４(a)中の★で印された十勝沖地震
震源及びその周囲に黒色の曲線で囲まれた地震性すべり領域（アスペリティ領域）から数十 km しか離れていない，相似地震グループ３，４に注目すると，ここでは GPS解析から示されているように，震源域から近いため余効すべりが卓越すると考えられるが，相似地震解析ではほとんど捉えられていない．これは，相似地震震源域で非相似地震が起きたために相似地震解析ではすべり量として積算されないためだと考えられる．それに対して相似地震グループ５，７，８では，両者が一致していることから，十勝沖震源域から十分離れたところでは，余効すべりによる速度場の揺らぎの影響が小さいため，相似性が失われないと考えられる．
４．まとめ「相似地震震源域では地震波形がそっくりな地震が毎回必ず起きるのか？」という疑問や問題点は，多くの研究者が気づいていたと考えられるが，その可能性及びメカニズムについては，これまで全くの未知のものであった．本研究では，速度/状態依存摩擦構成則に基づいた数値シミュレーションによって，巨大地震震源域近傍に相似地震震源域が共存する場合，大規模な余効すべり通過前後において相似地震震源域において非相似であるスロースリップイベントが発生することが示された．このことは十勝沖地震の余効すべりにおいて，GPS 解析結果と相似地震解析結果が定量的に一致していない現象を説明出来る可能性があると考えられる．
近年観測から捉えられているＭ４クラスの相似地震について，M７クラスの大アスペリティとの相互作用を数値シミュレーションで表現することは，二つのアスペリティのスケールの差が大き過ぎるために，計算コストが非常に高く困難であった．本研
究では，FFT を用いた計算の高速化に成功したことで，それを実現することが十分可能なものとなった．実際の沈み込みプレート境界面は，図１に示すように，東北日本であれば三陸はるか沖地震や十勝沖地震や宮城県沖地震など巨大なアスペリティが複数存在する一方で，釜石沖のような小規模のもの，発生間隔がそれほど規則的でないものも含めれば，それよりさらに数桁のスケールで規模の小さい相似地震震源域が多数存在し，それらが相互作用の影響を受けて複雑な振る舞いをしているのが実際の場であると考えられる．従って，実際の地震発生過程について数値シミュレーションからさらに踏み込んだかたちで説明するためには，今後さらに大規模な数値シミュレーションが出来るように計算環境を整える必要が迫られている．また，今回は計算コストの軽減のために，すべりの方向を固定し，有効法線応力の時間変化を無視した．
これらもいずれ取り入れることで，より高精度の下で定量的な議論が出来ると期待される．

(a)
▲

☆
〇◆
■

(b)
▲

☆
◆■

図３．(a)で印した点[あるいは座標（ X , W ）]における(b)積算すべり履歴図．赤と青の実線は，〇の位置での履歴に相当し，赤青の間にある黒色の実線は，小アスペリティ中心点での履歴を示す．灰色の実線は，小アスペリティ領域のみが単独で存在する場合のものを示し，V>1mm/secは，共存する場合の小アスペリティ領域中心点でのすべり履歴のうち，すべり速度が1mm/secを超えた地震性すべりのものを示す．

図４．（ａ）（左）相似地震解析によって得られた 2003 年十勝沖地震（M8.0)後から 2004 年釧路沖の地震（M7.1)までの期間の準静的滑り速度分布．カラースケー
ルは相似地震解析によって得られた準静的滑り速度分布，青コンターは GPS から推定した余効滑り速度分布．赤の数字はそれぞれの相似地震グループである．
（ｂ）（上）GPS と相似地震解析から推定されたプレート面すべり履歴．青が GPSの推定値，赤が相似地震解析による推定値を示す．左の番号は(a)のそれぞれの相似地震グループの領域番号を示す．[21]参考文献[1] 永井理子・菊地正幸・山中佳子，三陸はるか沖における再来大地震の震源過程の比較研究―1968 年十勝沖地震と 1994 年三陸はるか沖地震の比較，地震 2, 54,
267-280, 2001.
[2] Yamanaka Y., and M. Kikuchi, Source process of the recurrent Tokachi-okiearthquake on September 26, 2003, inferred from teleseismic body waves,

Earth Planets Space, 55, e21-e24, 2003.
[3] Yamanaka Y., and M. Kikuchi, Asperity map along the subduction zone innortheastern Japan inferred from regional seismic data, J. Geophys. Res., 109,10.1029/2003JB002683, 2004.
[4] Igarashi, T., T. Matsuzawa, A. Hasegawa, Repeating earthquakes andinterpolate aseismic slip in the northeastern Japan subduction zone, J.

Geophys. Res., 108, 10.1029/2002JB001920, 2003.
[5] Matsuzawa,T.,

T.

Igarashi,andA.

Hasegawa,

Characteristic

small-earthquake sequence off Sanriku, northeastern Honshu, Japan,Geophys. Res. Lett., 29, doi:10.1029/2001GL014632, 2002.
[6] Uchida, N., T. Matsuzawa, T. Igarashi, and A. Hasegawa, Interplate
quasistatic slip off Sanriku, NE Japan, estimated from repeating earthquakes,Geophys. Res. Lett., 30, doi:10.1029/2003GL017452, 2003.
[7] Uchida, N., A. Hasegawa, T. Matsuzawa, and T. Igarashi, Pre- andpost-seismic slip on the plate boundary off Sanriku, NE Japan associatedwith three interplate earthquakes as estimated from small repeatingearthquake data, Tectonophysics, 385, 1-15, 2004.
[8] Uchida, N., T. Matsuzawa, A. Hasegawa, and T. Igarashi, Recurrence
intervals of characteristic M4.8+/-0.1 earthquakes off Kamaishi, NE Japan Comparison with creep rate estimated from small repeating earthquake data,Earth Planet. Sci. Lett., 233, 155-165, 2005.
[9] Okada T., T. Matsuzawa, and A. Hasegawa, Comparison of source areas ofM4.8+-0.1 earthquakes off Kamaishi, NE Japan - Are asperities persistentfeatures?, Earth Planet. Sci. Lett., 213, 361-374, 2003.
[10] Boatwright, J., and M. Cocco, Frictional constraints on crustal faulting, J.

Geophys. Res., 101, 13895-13909, 1996.
[11] Liu, Y., and J. R. Rice, Aseismic slip transients emerge spontaneously inthree-dimensional rate and state modeling of subduction earthquake
sequences, J. Geophys. Res., 110, B08307, doi:10.1029/2004JB003424, 2005.
[12] Dieterich, J. H., Modeling of rock friction, 1, Experimental results andconstitutive equations, J. Geophys. Res., 84, 2161–2168, 1979.
[13] Ruina, A. L., Slip instability and state variable friction laws, J. Geophys. Res.,88, 10359-10370, 1983.
[14] Beeler, N. M., T. E. Tullis, and J. D. Weeks, The roles of time anddisplacement in the evolution effect in rock friction, Geophys. Res. Lett., 21,
1987-1990, 1994.
[15] Okada, Y., Internal deformation due to shear and tensile faults in a halfspace,Bull. Seism. Soc. Am., 82, 1018–1040, 1992.
[16] Rice, J. R., Spatio-temporal complexity of slip on a fault, J. Geophys. Res., 98,9885– 9907, 1993.
[17] Savage, J. C., A dislocation model of strain accumulation and release at asubduction zone, J.Geophys. Res., 88, 4984-4996, 1983.
[18] Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical

Recipes, 2nd. ed., Cambridge Univ. Press, New York, 1992.
[19] Nadeau, R. M., and L.R. Johnson, Seismological studies at Parkfield: VI.
Moment release rates and estimates of source parameters for small repeating
earthquakes, Bull. Seismol. Soc. Am., 88, 790–814, 1998.
[20] Kawasaki, I., Silent earthquakes occurring in a stable-unstable transitionzone and implications for earthquake prediction, Earth Planet. Sci., 56,813-821, 2004.
[21] 油井智史・三浦哲・長谷川昭・八木勇治, 2003 年十勝沖地震(M8.0)前後の北海道
地方 のプレート間カップリング状況, 日本地震学会 2005 年秋季大会, 札幌, 北海道大学, 2005 年 10 月.