[SSS13-P12] 2016年4月15日熊本の地震(Mj6.4)の運動学的および動力学震源モデル
キーワード：2016年熊本地震、震源過程、動力学シミュレーション
2016年熊本地震の地震群のうち，3番目に大きい2016/4/15 00:03に発生したMj6.4の地震について，まず強震記録のインバージョン解析により震源過程を求め，その結果をもとに，運動学的震源モデルを再現するような動力学震源モデルの構築を試みた．まず，強震記録のマルチプルタイムインバージョン解析により震源過程を求めた．解析には，防災科学技術研究所K-NET, KiK-net，気象庁および自治体震度計の計11観測点を使用した．速度構造モデルは，Yoshida et al. (2017)による余震の波形フィッティングによる構造モデルのチューニングにより求めた．この時，余震の震源位置にはdouble-difference法により再決定されたもの(Uchide et al., 2016)を用いた．本震の断層面は1.5 km×1.5 kmの小断層に，時間方向にはパルス幅0.8秒のスムーズドランプ関数を0.4秒間隔で6個並べ離散化した．各小断層からの理論グリーン関数の計算には，離散化波数法(Bouchon, 1981)と反射・透過係数行列法(Kennett and Kerry, 1979)により点震源の波形を計算して求めた．逆解析にはマルチタイムウィンドウインバージョン法(Hartzell and Heaton, 1983)を用いた．解析対象には，観測加速度記録に0.05－1.0 Hzのバンドパスフィルターを適用し，1回積分して作成した速度波形とした．P波到達から18秒間を切り出し，解析に用いた．なお，KiK-net 観測点では地中波形を逆解析の対象とした．破壊伝播速度VFT は，平滑化係数λを固定して，2.0～2.6 km/sの範囲で0.2 km/sごとに解を求め，それらの残差が最小となった2.4 km/sを採用した．各機関(USGS W-phase, USGS body wave, USGS Regional, Global CMT, JMA, F-net)のCMT解のM0は平均が1.2×1018 Nmであったので，採用したVFTのうちM0がもっともこの値に近くなるようなλの解を採用した．解析により得られた震源モデルでは，破壊開始点より浅部に，最大滑り量は約1 m，平均すべり量が0.6～0.7 mのアスペリティと高すべり速度域（High slip/moment Rate Area, HRA; 吉田・他，2017）が求められた（図左）．アスペリティの上端では，すべり量と比べてすべり速度が小さく，同定されたHRAはアスペリティよりも上端の1メッシュ分小さいものが求められた．この解析で得られた断層面積，アスペリティ面積などはスケーリング則に良く対応した．
次に，この地震のインバージョン解析で得られた震源モデルの応力分布や経験的グリーン関数(EGF)法による解析結果（倉橋・他，2017）などをもとに動力学シミュレーションを行った．広域応力場はUrata et al (2017)などを参考に，σ1=100 MPa，σ3=50 MPaとした．断層面上の応力降下量分布は，インバージョン解析で得られたすべり量分布からOkada (1992)の方法で静的応力降下量分布を計算し，アスペリティやHRAで平均したものをパッチとして与えたモデル，すなわち応力降下量分布を特性化して設定した．設定した静的応力降下量は，それぞれアスペリティモデルでは3.9 MPa, HRAモデルでは4.9 MPa，背景領域は両者とも1.0 MPaである．破壊強度の設定は，設定した応力降下量をもとに，EGFによるSMGAの応力パラメータを実効応力とみなして，実行応力=破壊強度+静的応力降下量となるよう設定した．臨界すべり量Dcは，既往研究を参考にしたほか，KMMH14においてDc’’(Fukuyama and Mikumo, 2007)を求めたところ15 cm程度であったこと，予備検討の結果などから，ここでは断層面全体で10 cmとした．シミュレーションは差分法で行い，格子間隔は50 m×50 m，時間間隔は0.004秒として8秒間分を計算した．破壊が断層面全体に伝わるよう，実効応力がHRAのモデルでは6.3 MPaに，アスペリティでは5.1 MPaとなるよう破壊強度を調整した（両者で実効応力が異なるのは，背景領域の応力降下量が面積の違いにより異なるため，図左下）．シミュレーションの結果は，運動学的震源モデルのモーメントや滑り分布（アスペリティやHRAで0.6 m程度）を良く再現するものが得られた（図右）．ただし，アスペリティ上端部でのすべり量がやや小さいこと，アスペリティやHRAの上端で大きなすべり速度のパルスが現れることがあること，断層面上の一部でスーパーシア破壊（S波速度を超える破壊伝播速度の破壊）が見られること，といった課題が残った．臨界すべり量など，さらに動力学パラメータの検討を進める必要があると考えられる．
謝辞 防災科学技術研究所K-NET, KiK-net, 気象庁震度計，御船町の記録を使った．震源位置は産業技術総合研究所内出崇彦博士に提供いただいた．この研究の一部は原子力規制庁「平成30年度原子力施設等防災対策等委託費（内陸型地震による地震動の評価手法の検討）事業」の成果である．
次に，この地震のインバージョン解析で得られた震源モデルの応力分布や経験的グリーン関数(EGF)法による解析結果（倉橋・他，2017）などをもとに動力学シミュレーションを行った．広域応力場はUrata et al (2017)などを参考に，σ1=100 MPa，σ3=50 MPaとした．断層面上の応力降下量分布は，インバージョン解析で得られたすべり量分布からOkada (1992)の方法で静的応力降下量分布を計算し，アスペリティやHRAで平均したものをパッチとして与えたモデル，すなわち応力降下量分布を特性化して設定した．設定した静的応力降下量は，それぞれアスペリティモデルでは3.9 MPa, HRAモデルでは4.9 MPa，背景領域は両者とも1.0 MPaである．破壊強度の設定は，設定した応力降下量をもとに，EGFによるSMGAの応力パラメータを実効応力とみなして，実行応力=破壊強度+静的応力降下量となるよう設定した．臨界すべり量Dcは，既往研究を参考にしたほか，KMMH14においてDc’’(Fukuyama and Mikumo, 2007)を求めたところ15 cm程度であったこと，予備検討の結果などから，ここでは断層面全体で10 cmとした．シミュレーションは差分法で行い，格子間隔は50 m×50 m，時間間隔は0.004秒として8秒間分を計算した．破壊が断層面全体に伝わるよう，実効応力がHRAのモデルでは6.3 MPaに，アスペリティでは5.1 MPaとなるよう破壊強度を調整した（両者で実効応力が異なるのは，背景領域の応力降下量が面積の違いにより異なるため，図左下）．シミュレーションの結果は，運動学的震源モデルのモーメントや滑り分布（アスペリティやHRAで0.6 m程度）を良く再現するものが得られた（図右）．ただし，アスペリティ上端部でのすべり量がやや小さいこと，アスペリティやHRAの上端で大きなすべり速度のパルスが現れることがあること，断層面上の一部でスーパーシア破壊（S波速度を超える破壊伝播速度の破壊）が見られること，といった課題が残った．臨界すべり量など，さらに動力学パラメータの検討を進める必要があると考えられる．
謝辞 防災科学技術研究所K-NET, KiK-net, 気象庁震度計，御船町の記録を使った．震源位置は産業技術総合研究所内出崇彦博士に提供いただいた．この研究の一部は原子力規制庁「平成30年度原子力施設等防災対策等委託費（内陸型地震による地震動の評価手法の検討）事業」の成果である．