3.9 データ活用予測研究(1) 業務の内容(a) 業務題目「データ活用予測研究」(b) 担当者
所属機関

役職

氏名

国立大学法人東北大学災害科学国際研究所

教授

日野亮太

国立大学法人東北大学大学院理学研究科

助教

太田雄策

助教

内田直希

教授

平原和朗

准教授

宮崎真一

准教授

西村卓也

准教授

伊藤喜宏

教授

鷺谷威

国立大学法人名古屋大学大学院

教授

古本宗充

環境学研究科附属地震火山研究センター

助教

伊藤武男

独立行政法人産業技術総合研究所

研究員

安藤亮輔

独立行政法人海洋研究開発機構

サブリーダー

堀高峰

地震津波・防災研究プロジェクト

特任技術研究副主任

中田令子

特任研究員

兵藤守

技術研究副主任

有吉慶介

国立大学法人京都大学大学院理学研究科

国立大学法人京都大学防災研究所

国立大学法人名古屋大学
減災連携研究センター

活断層・地震研究センター(c) 業務の目的
・地殻変動観測データのコンパイルと解析、データベース化
京 コ ン ピ ュ ー タ 等 で 計 算 さ れ る 多 数 の 地 震 シ ナ リ オ と そ の 前 後 の ゆ っ く り す べりのシミュレーション結果のデータベースと、GEONET・DONET・水準測量や三辺・三角 測 量 等 の 地 殻 変 動 デ ー タ や 相 似 地 震 の デ ー タ を コ ン パ イ ル し た 観 測 デ ー タ ベ ースを構築する。
・プレート境界すべりの推移予測の妥当性検証と予測の試行
過去の地震やゆっくりすべりの観測データを逐次入力し、シミュレーションデータ
ベースと比較する同化の模擬テストを行なう。また、次の豊後水道のゆっくりすべりや東北地方太平洋沖地震後の余効すべり等に対して予測の試行実験を行う。これらを通して同化手法や予測の問題点を洗い出し、改善に必要な知見を得る。
・逐次データ同化手法の改良
精 度 の 異 な る 複 数 デ ー タ を 用 い る 場 合 の 尤 度 評 価 や リ ア ル タ イ ム で デ ー タ を 取り
込 み な が ら モ デ ル パ ラ メ タ や 初 期 値 を 更 新 で き る よ う に 同 化 手 法 を 改 善 す る と ともに、地殻変動以外の観測データ（特に地震活動データ）の同化手法を開発する。
293

(d) ８か年の年次実施業務の要約

平成 25 年度:
これまでに整備した地殻変動データのコンパイル結果に最新データを加えて、追加解析を行うとともに、過去の豊後水道や八重山のゆっくりすべりで逐次同化の模擬テストを行う準備を進めた。また、観測データベースの仕様を検討し、東北地方太平洋沖地震前後のデータのコンパイルを始めた。複数データの尤度評価手法の開発を進めた。
平成 26 年度:
前年度までに整備した地殻変動データのコンパイル結果に最新データを加えて、さらに追加解析を行うとともに、過去の豊後水道や八重山のゆっくりすべりで逐 次同化の模擬テストを行う。また、観測データベースの仕様を決め、東北
地方太平洋沖地震前後のデータを引き続きコンパイルする。複数データの尤度評
価手法を開発する。
平成 27～28 年度：
東北地方太平洋沖地震前後のデータを解析し、同化手法の模擬テストを行う。
モデルパラメタを逐次更新する手法を開発する。
平成 29～30 年度：
模擬テスト結果にもとづいて同化手法を改善するとともに、この期間中に発生が期待される豊後水道のゆっくりすべりでリアルタイム予測試行実験を行う。地
震活動データの同化手法を提案するとともに、観測データベースに地震活動を加える。
平成 31～32 年度：
模擬テスト結果や手法改良研究にもとづいて同化手法の改善を進めるとともに、南海トラフならびに日本海溝において、リアルタイム予測試行実験を行い、予測の問題点を整理する。また、地殻変動・相似地震と地震活動データを用いた
同化システムのプロトタイプを構築する。
(e) 平成 25 年度業務目的
余効すべり・ゆっくり滑りデータ同化手法を開発し、八重山や豊後水道ゆっくり滑りの観測データを整理して、観測データをシミュレーションデータベースと比較するための研究開発を行う。また、東北地方の海陸地殻変動観測のデータを用いて、2011
年東北地方太平洋沖地震の発生前後に起こった、プレート間固着強度の時空間ゆらぎ（余効変動を含むゆっくり滑りイベント、前震活動、余震活動）の同化の模擬テストに着手する。さらに、ゆっくりすべりの観測事例や発生機構に関する海外の研究動向の情報収集を行う。

(2) 平成 25 年度成果
① 余効すべり・過去のゆっくりすべりデータ同化手法の開発
294

(a) 業務の要約
1) 余効すべりのデータ同化手法の開発
大気海洋分野で開発されたデータ同化手法の一手法であるアジョイント法（４次元変分法）を用いて、余効変動データからプレート境界における摩擦パラメータ分布を推定する手法の開発を行い、2003 年十勝沖地震の余効変動実データに適用し、プレー
ト境界における余効すべり域の摩擦パラメータ分布の推定を試みた。
最初に摩擦パラメータの推定可能性を検討するために、単純な系を用いた数値実験を行なった。すべり速度と状態に依存する摩擦則と、状態変数の発展則として、スローネス則（Dieterich, 1979）を用いて、モデル変数（すべり速度 V および状態変数 θ）の初期値や地震発生域の摩擦パラメータは既知とし、余効すべり域の摩擦パラメータ(a, a－b, L)のみを推定した。まず、断層面上のすべり速度を観測データとした場合、観測データが余効すべりの加速過程の情報を含めば、全ての摩擦パラメータが真値か
ら±10%の範囲内で求まることが示された。また、地表面における変位を観測データ
として同様の実験を行なった場合も結果はほぼ同様で、実際の GPS データを用いた摩擦パラメータ推定が可能であることが示された（Kano et al., 2013）。
この結果を受けて、開発した手法を 2003 年十勝沖地震後の余効すべりに適用した。
観測データとして国土地理院により公表されている GPS 解析結果を用いた。また、震源域の摩擦パラメータは、シミュレーションにおいて現実的な地震サイクルが得られる値に固定し、すべり速度の初期値は複数の地震サイクルを経て発生した地震直後の値に固定した。また、初期時刻では定常状態にあるとした。断層面の摩擦パラメータの空間分布を６通り仮定し、これらに対して実際の GPS データを同化した。得られた赤池情報量規準（AIC）の値から最適なモデルを選択し、余効すべり域の摩擦パラメータ分布を推定した。その結果、本震震源域の東と西とで摩擦特性が有意に異なることや、余効すべり域の摩擦特性が有意に速度強化を示すことが示唆された。また、a－b の値は 10 -3 程度で先行研究の結果と整合的であり、さらに、新たに L や a の値の推定にも成功した。ただ、求まった解は初期条件に依存することが明らかになり、初期条件も推定する必要があるが、状態変数 θ の初期条件の推定に問題が残る（Kano,2014）。
2) 八重山スロースリップのデータ同化a) スロースリップ検出法の開発と八重山スロースリップへの適用
AIC を用いて GPS の時系列データのオフセットを自動検出し、スロースリップイ
ベント（SSE）を検出する手法（Nishimura et al., 2013）に改良を加え、南琉球域の八重山地域の GPS 時系列データに適用し SSE を検出した。適用において、AIC の差を取る際、継続時間を含む指数関数付きの直線を用いることで、SSE 検出の精度が向上した。Heki and Kataoka, 2008 により既に八重山 SSE の検出は行われているが、新手法により検出された SSE を彼らの研究と比較すると、彼らの解析期間より後に発生した SSE では再来間隔が大きくなる傾向があることが分かった。

295b) データ同化を目指した八重山 SSE の再現シミュレーション
八重山 SSE を起こしているプレート境界面上の摩擦パラメータを推定するデータ
同化手法開発のため、フォワード計算により、観測された約半年の平均的な再来間隔を再現するシミュレーションを行い、摩擦パラメータを見積もった。また、八重山地域では 1771 年八重山津波が発生しており、琉球海溝に近い浅部域で津波地震が発生した可能性が指摘されている（Nakamura, 2009）。また津波石の年代測定から繰り返しこのような津波地震が発生しており、琉球海溝に近い浅部域では固着している可能性が指摘されている。そこで、SSE 領域より浅部での固着の影響が、SSE に及ぼす影響（SSE で解放されるすべり速度への影響）を調べた。

(b) 業務の実施方法
1) 余効すべりのデータ同化手法の開発
すべり速度と状態に依存する摩擦則（RS 則）、そして例えば状態変数 θ がスローネ
ス則（Dieterich, 1979）に従うとして、プレート境界面上の分割したセル上に摩擦パラメータを分布させ、すべり速度 V と状態変数θの初期値を与えて、プレート運動を駆動力とする準動的地震サイクルシミュレーション（Rice, 1993）が行われ、摩擦パラメータ分布を試行錯誤的にチューニングして、再来間隔などの地震発生履歴を再現している。本 研究 では、 デ ータ同化 の一 手法で あ るアジョ イン ト法（ 例 えば、淡 路ほ か,2009）を用いて、地震後の余効変動をデータとして、余効すべり域での摩擦パラメー
タ分布を推定する手法を開発した。

アジョイント法による推定
各セル i のすべり速度 V i と状態変数 θ i をまとめて次元 N（セル数）のベクトルで表し、V と θ とする。また摩擦パラメータ A、B、L を推定する領域の次元を持つベクトル A、B、L とする。時刻 n でのすべり速度、状態変数および摩擦パラメータを表す変数ベクトル X(n)

をX( n ) = ( V( n ) , θ( n ) , A − B, B, L)T（１）と定義する。時刻 n=０でのすべり速度と状態変数の初期値を設定し、準動的運動方程式を時間積分して、各セル内でのすべり速度と状態変数の時間発展を追う、時刻 n から n＋１へ進めるフォワード計算を

X ( n +1) = M ( n ) X ( n )（２）と表す（具体的な表式は例えば Kano (2014)参照のこと）。従って、X ( n ) = M ( n −1) X ( n −1) = M ( n −1) M ( n −2) X ( n −2) = M ( n −1) M ( n −2) L M (1) X (0)

（３）

と表せる。
通常観測データとして GPS で観測された地表変位を用いることになるが、観測値d(n)に対応する計算された地表変位速度ベクトルはd cal ( n ) = H ' V ( n ) = HX ( n )（４）と書ける。ここで H ' は断層での変位速度から地表変位速度を求める行列で、各要素
296

は例えば Okada, 1992 に従い計算できる。
次に、以下の評価関数

J ( X (0) ) =

−1
1 (0)
1 n =nT( X − X b )T B −1 ( X (0) − X b ) + ∑ ( H ( X ( n ) ) − d ( n ) ) T R ( n ) ( H ( X ( n ) ) − d ( n ) )2
2 n =1（５）bを考える。ここで、 X は背景ベクトルと呼ばれ初期値の初期推定値ベクトルである。
この評価関数を最小化するために、勾配を計算し変形すると以下のように書ける。

∂J / ∂X (0) = B −1 ( X (0) − X b ) + M (1)T ( M ( 2)T (L ( M ( nT −1)T r ( nT ) + r ( nT −1) ) L) + r ( 2) ) + r (1) )
（６）
ここで、

M ( n ) = ∂M ( n ) / ∂X ( n )−1r ( n ) = H T R ( n ) d ( n ) − H ( X ( n ) ) （７）である。すなわち、勾配を求めるには、データのミスフィット r(n)にアジョイント演算

子 M ( n −1)T を作用させ、 r ( n−1) を加えるという計算を、初期値を r T として n = nT から１(n )まで繰り返し、最後に背景誤差に関する勾配を加えれば良いことになる。このように勾配を求めた後は、降下法のプログラムを用いて X (0) の値をより良い（評価関数がより小さくなる）推定値に更新し、再びフォーワードシミュレーションを行いデータミスフ
ィット（残差）を計算し、アジョイントモデルにより勾配を求めるという手順を、推定値が収束するまで繰り返す。

数値モデル実験
すべり速度をデータとする場合と地表地殻変動データを用いる場合とで、摩擦パラメ
ータ推定に関する考察を行うため、単純化された状況でプレート境界における余効すべりモデルを考える（Kano et al., 2013）。図３－９－①－１（左図）に示すように、傾斜角 20°で沈み込むプレート境界面上に一様なすべり強化摩擦特性（A－B >０）を与え、その中に速度弱化（A－B <０）摩擦特性を持つ２つのアスペリティ（LA と SA）を配置し、LA は約 100 年、SA は平均 60 年の繰り返し間隔で破壊する地震サイクル計算を行い、あるタイミングで LA が破壊したあとの、余効すべりの伝播シミュレーションを行う。その際プレート境界面上の点 A～E でのすべり速度をデータとして、アジョイン
ト法で余効すべり域の摩擦パラメータが推定できるかという問題設定をした（実験１）。
次に図３－９－①－１（右図）に示すように、陸側地表に 01～50 の GPS 観測点を、海底に圧力計観測点を配置する。シミュレーションにより得られた余効すべりから Okada(1992)に従い、これらの観測点でのデータを作成するが、GPS 点では３成分変位を、海底観測点では、上下動データのみをノイズを加えて模擬データを作成する。

297

図３－９－①－1（左図）傾斜角 20°で沈み込むプレート境界面上での摩擦分布。LA とSA の四角内のみ速度弱化（A－B <０）。残りの余効すべり域（afterslip region）には速度強化(A－B >０)摩擦特性を与える。黒丸 A－E はプレート境界面上ですべり速度を観測として与える点を示す。（右図）摩擦モデルは同じで、黒丸は、陸域 GPS 観測点 01～50 の位置、黒四角は海域の水圧計海底観測点を示す。

2003 年十勝沖余効変動データへの適用
次に開発されたデータ同化手法の 2003 年十勝沖地震の余効変動データへの適用を示す。データは図３－９－①－２に示す GPS 観測点の水平２成分変位で、地震後２日
目から 15 日目までのデータを用いる。

図３－９－①－２

2003 年十勝沖地震後２日目から 15 日目までの間の累積余効変動の水平ベクトル。
（緑：観測値、赤：推定された摩擦パラメータ、すべり速度から計算される値）

298

図３－９－①－３

十勝沖地震断層モデルとその領域内では摩擦パラメータは一定として推定する６つの余効すべり領域（１～６）をもつ推定モデル６。
西の赤い矩形領域は、2003 年十勝沖地震震源域、東の赤い矩形領域は、2004 年釧路沖地震
震源域を表す。

図３－９－①－４

八重山地域における GPS 観測点配置（右図）と琉球狐の GPS ベクト

ル図および八重山地域の位置（Heki and Kataoka, 2008）
2) 八重山スロースリップのデータ同化a) SSE 検出法の開発と八重山スロースリップへの適用
Nishimura et al. (2013)では、GPS 時系列データから以下の方法で SSE を検出している。
i) GPS データの 360 日移動平均を用いて短周期成分の除去。
ii) 180 日のタイムウィンドウを取って１日ずつ移動。
iii) 時系列にフィットする直線と中央にオフセットの入った直線の AICの差 ∆AIC を計算。
iv) ∆AIC が負の値を取る場合、SSE が発生している可能性有りと判定。
299

以上に加えて、∆AIC を計算する際、Heki and Kataoka (2008)で得られている以下の座標値 x ( t ) の時間発展式n−(t − Ti ) x (t ) = at + b + ∑ X i 1 − expτ i 
i =1
(T1 < T2 < L < Tn < t )（８）を用いている。ここで、 b はオフセット、a は長期間にわたる背景変位速度で地震間の地殻変動を表す。i 番目の SSE の発生時刻は Ti で、時定数 τ i は 0.1 年～0.15 年の値を取る。
図３－９－①－４に示す GPS 観測点で得られた時系列データに上記 SSE 検出法を適用した。

b) データ同化を目指した八重山 SSE の再現シミュレーション
Heki and Kataoka , 2008 では、八重山域でのプレート収束速度は 12.5 cm/年で、石垣島下に傾斜角 15°で沈み込むフィリピン海プレート境界面上の深さ 20 km から 40
km に最大 70 km x 90 km の大きさを持つ SSE 断層モデルを提出している。これを参考に、図３－９－①－５に示すプレート境界モデルを作成（深さ 32 km 以深の速度強化（A－B>０）域内に SSE を発生させるパッチを設定）し、速度状態依存摩擦則・スローネス則を用いて、SSE 発生域の摩擦パラメータ分布を推定した。
アスペリティは速度弱化摩擦特性(A－B<０)を持つが、そのサイズを R 、摩擦パラ
メータで決まる臨界半径（臨界核形成サイズ）を Rc =π

GBL
（ G ：剛性率）と4 ( B − A)2

すると、すべりが進展しサイズが Rc を超えると高速破壊を起こす。逆に言えばアス
ペリティサイズが臨界核サイズより小さい（ R < Rc ）と高速すべりを起こせないで、
スロースリップで終わってしまう。このように速度弱化摩擦特性を持つが、そのサイズが臨界核形成サイズより小さい状態を、条件付き安定すべりという。本研究では、この条件付き安定すべり摩擦特性を利用して SSE を発生させる。

図３－７－１－①－５

八重山 SSE を起こすアスペリティのモデル設定。

300

(c) 業務の成果
1) 余効すべりのデータ同化手法の開発
数値モデル実験
図３－９－①－１に示したモデル実験は、全余効すべり域が同じ摩擦パラメータを持つ場合のモデル実験で、初期値は模擬データを作成したフォワード計算の時と同じ値を用いている。すなわち、すべり速度データまたは地表変位データから摩擦パラメ
ータだけの推定を行なっている。
これらの実験から、以下のことが分かった。１）すべり速度データでも地表変位データでも、真値より１桁大きなおよび小さな範囲の背景値からでも摩擦パラメータを推定できる。２）プレート境界上の点 A－E における合成模擬すべり速度データを図３－９－①－６に予測値と併せて示しているが、(B)のように減速フェイズに加えて、速度の加速フェイズを含むデータが１点でもあると、全ての摩擦パラメータが真値か
ら±10%の範囲内で求まる。３）一方、地表変位データから直接データ同化により摩擦パラメータを推定する実験では、数点の観測点で、ノイズより大きなシグナルを持
つ大きな振幅を持っていれば全ての摩擦パラメータが推定可能である。４）海底観測
点データに大きなノイズを加えたため、摩擦パラメータ推定には大きな寄与はなかったが、ノイズを抑えることができれば、摩擦パラメータ推定に大きな役割を果たすと期待される。

図３－９－①－６
図３－９－①－１の左図に示した各点 A－E における模擬すべり速度データ（丸印）とデータ同化に用いた初期背景値（赤線）とデータで同化により得られた予測値（緑線）を示す。
301

2003 年十勝沖地震余効変動データへの適用
数値モデル実験では、余効すべり域内の摩擦パラメータは一定としていたが、ここでは領域に分けて推定を行なった。Miyazaki et al., 2004 は GPS 観測データから、余効
すべり域でのすべりの発展をインバージョンにより求めている。そのすべり速度分布を参考に領域分けを行い、図３－９－①－３にあるように１～６の領域を設定し、その全てを推定するモデル６に対し、一様なモデル１、領域１～５を一様とするモデル
２、領域１～３、４～５、６の３領域の摩擦パラメータを推定するモデル３、領域１～２、３、４～５、６の４領域のパラメータを推定するモデル５で推定して、AIC でモデルを判定した。その結果、モデル１、２とそれ以降のモデルでは AIC に大きな差が出た。これは領域１～３と領域４～５では、摩擦パラメータが有意に異なっている
ことを示している。すなわち、本震震源域の東と西とで摩擦特性が有意に異なることを示している。更に、余効すべり域の摩擦特性が有意に速度強化を示すことが示唆された。
色々なデータやモデル設定でデータ同化を行なったが、最終的には、図３－９－①－２に示した GPS 観測点での地震後２日から 15 日目までの水平２成分の時系列デー
タを用い、すべり速度の初期値としては、Miyazaki et al., 2004 によるすべり速度を、定常状態を仮定して状態変数 θ を θ＝L/V とした。また、最初はデータおよびすべり
速度方向は海溝に直交する方向にとっていたが、余効すべりの時間発展が、Miyazaki etal., 2004 によるものと大きく異なっていたため、太平洋プレートの斜め沈み込みを考慮してすべり速度方向を決め、摩擦パラメータのデータ同化を行なった。そのデータ同化による予測ベクトルは、図３－９－①－２に観測値とともに示してある。
推定された摩擦パラメータの値を、表３－９－①－１に示す。ただし、有効法線応力は一定で 20 MPa としている。

表

３－９－①－１

モデル６による生成摩擦パラメータの値

モデル６

領域１

領域２

領域３

領域４

領域５

領域６

L (mm)

22.74

71.95

2.45

25.09

0.08

62.46

-4

a (x10 )

3.87

32.85

21.02

25.04

12.06

31.51

-4

8.08

3.04

8.09

5.96

5.89

3.21

b (x10 )

最後に、モデル６により推定されたパラメータで得られた積算すべりの時間発展を、
図３－９－①－７に示す。これらのすべり分布は Miyazaki et al., 2004 の Fig.４.１の結果とよく似ている。まだ相違が見られるが、運動学的インバージョンにおける時間および空間方向での平滑化の影響によるものかも知れない。本研究で開発したデータ同化手法では、こういった時空間方向での平滑化は行なっていない。
Kano et al., 2010 では、３自由度の簡単な断層モデルで、初期条件と摩擦パラメータの推定を行い、両者の同時推定は、摩擦パラメータ推定に大きな影響を与えないとしている。彼らは、速度の初期値はよく拘束されるが、状態変数の初期値はあまり拘束302できないと述べている。確かにすべり速度の初期値は地表データから推定可能であろうが、状態変数の初期値推定には困難が伴う。従って、状態変数の背景値は大きな推定誤差を与えて求める必要がある。

図３－９－①－７

モデル６による積算すべり量の時間発展。

2) 八重山スロースリップのデータ同化a) SSE 検出法の開発と八重山スロースリップへの適用
直線とオフセット方式(A)と直線と時間発展（式８）(B)による SSE 検出率を先行研
究（Heki and Kataoka, 2008）の結果と比較すると以下のようになる。

上表の、本研究 A と B の精度良く検出された数の比較から分かるように、自動でSSE を検出するアルゴリズムに式（８）の時間発展を入れると、SSE 検出の精度が上昇し、Mw や検出時期も丹念に調べた先行研究に近づいてくるのが分かる。
Heki and Kataoka (2008)では、1997～2007 年までの解析であるに対し、それ以降 2013
年までの解析を行い、９個の SSE を検出することができた。図３－９－①－８に宮303

古島に対する波照間観測点での変動量の N20°W 成分の時間変化と検出 SSE を示す。
対応する 18 個の SSE に対する再来間隔は、Heki and Kataoka (2008)の 6.31 ヶ月に対し、本研究では、6.26 ヶ月となった。また、本研究の全期間(1997～2013 年)では、28 個の SSE に対して、6.46 ヶ月となり、2008 年以降の SSE の再来間隔が大きくなっている。

図３－９－①－８

宮古島に対する波照間の変動量の N20°W 成分の時間変化と検出 SSE。

b) データ同化を目指した八重山 SSE の再現シミュレーション
図３－９－①－５に示したモデル設定で、約６ケ月の再来間隔を持つ SSE を発生させる摩擦パラメータは、プレート収束速度を 12.5cm/年として、A=50 kPa、B＝56 kPa、
L=2.2 mm となった。この場合の臨界半径は Rc = 80.63km となり、 R / Rc = 0.37 となる。
Heki and Kataoka , 2008 によれば、SSE により解放されるすべり速度は、11 cm/年
で、12.5 cm/年より小さい。１個のアスペリティをおいたシミュレーションでは、SSEにより解放されるすべり速度は 12.5 cm/年となる。そこで、琉球海溝よりの浅部に固着域を設定し、固着率を変えて、解放速度が 11.5 cm/年となる固着率を探すと、すぐ近傍に置いた場合、固着率が 40％程度であると、説明可能である。
(d) 結論ならびに今後の課題
アジョイント法を用いて、余効すべり速度や地表での余効変動変位をデータとする、
摩擦パラメータ推定のデータ同化手法の開発を行い、2003 年十勝沖地震余効変動デー
タに適用して、摩擦パラメータ推定を試みた。先行研究の運動学的インバージョンで得られた余効すべりの時間発展と同じような結果が得られたが、細部においてはまだ異なっている。運動学的インバージョンでは、時間的空間的平滑化が行われているのに対し、304

本手法ではそういった平滑化は行われていなくて、摩擦構成則に従ってすべりが発展するという物理的拘束条件が課せられているという違いがある。また、本手法では、初期値を必要とし、状態変数の初期値推定に問題が残されている。今後解決すべき課題であろう。
八重山スロースリップのデータ同化研究では、GPS データから SSE を自動検出する手法を八重山 SSE に適用し、先行研究との比較を試みた。オフセットではなく時間発展式を加えることで精度良く検出された SSE の自動検出数が向上した。2008 年以降の
SSE の再来間隔はやや大きくなっている可能性がある。また、データ同化を目指した
SSE の再現シミュレーションにより、平均再来間隔を説明するアスペリティモデルを構築した。今後は再来間隔のばらつきや SSE による解放すべり速度を説明するモデル構
築と実際のデータ同化が課題となろう。

(e) 引用文献
1) 淡 路 敏 之 , 蒲 地 政 文 , 池 田 元 美 , 石 川 洋 一 （ 編 ） , デ ー タ 同 化 , 京 都 大 学 学 術 出 版会,pp.284,2009.
2) Dieterich, J. H., Modeling of rock friction 1. Experimental results and constitutive
equations, J. Geophys. Res., 84, 2161-2168, 1979.
3) Heki, K., and T. Kataoka, On the biannually repeating slow-slip events at the Ryukyu
Trench, southwestern Japan, J. Geophys. Res., 113, B11402, doi:10.1029/2008JB005739,
2008.
4) Kano, M.,

Estimation of frictional parameters in afterslip areas by assimilating GPS

data: Application to the 2003 Tokachi-oki earthquake, Ph.D Thesis, Earth and Planetary
Sciences, Graduate School of Science, Kyoto University, 2014.
5) Kano, M., S. Miyazaki, K. Ito, and K. Hirahara, Estimation of Frictional Parameters andInitial Values of Simulation Variables Using an Adjoint Data Assimilation Method with
Synthetic Afterslip Data (in Japanese with English abstract), Zisin 2, 63, 57–69, 2010.
6) Kano, M., S. Miyazaki, K. Ito, and K. Hirahara, An adjoint data assimilation method for
optimizing frictional parameters on the afterslip area, Earth Planets Space, 65, 1575-1580,
2013.
7) Miyazaki, S., P. Segall, J. Fukuda, and T. Kato, Space time distribution of afterslip
following the 2003 Tokachi-oki earthquake: Implications for variations in fault zone
frictional properties, Geophys. Res. Lett., 31, L06623, doi:10.1029/2003 GL019410,
2004.
8) Nakamura, M., Fault model of the 1771 Yaeyama earthquake along the Ryukyu Trench
estimated

from

the

devastatingtsunami,Geophys.

Res.

Lett.,

36,

doi:

10.1029/2009GL039730, 2009.
9) Nishimura, T., T. Matsuzawa, and K. Obara, Detection of short-term slow slip events
along the Nankai Trough, southwest Japan, using GNSS data, J. Geophys. Res., 118,
3112–3125, doi:10.1002/jgrb.50222, 2013.
30510) Okada, Y., Internal deformation due to shear and tensile faults in a half

space, Bull.

Seismol. Soc. Am., 82, 1018–1040, 1992.
11) Rice, J. R., Spatio-temporal complexity of slip on a fault, J. Geophys. Res., 98,
9885–9907, 1993.

② 予測システムの検証・予測試行実験(a) 業務の要約
次の豊後水道スロースリップイベントや 2011 年東北地方太平洋沖地震前後のゆっくりすべりに対する予測の試行実験に向けて、南海トラフ地域における観測データベ
ース構築用のシステムを日本海溝地域へ拡大するとともに、シミュレーションデータ
ベースの源となるモデルの構築を行った。モデル構築は、シミュレーションによるこれらの現象の再現性や、予測精度向上のための課題を把握するために、地震発生サイ
クルの数値シミュレーションを行った。豊後水道スロースリップイベントについては、
文部科学省の受託研究「東海・東南海・南海地震の連動性評価のための調査観測・研究」（平成 20～24 年度）の成果で得られた、南海トラフ巨大地震震源域西縁のモデルを基に、パラメタスタディを行った。得られた結果を詳細に検討し、当該モデルの妥当性を検証した。また、日本海溝地域においてモデル化に着手し、南海トラフ地域で用いているものと同様の仮定に基づくモデルを用いて、2011 年東北地方太平洋沖地震
前後の観測を、定性的に再現できることが分かった。両地域のモデルにおける課題も明確になった。

(b) 業務の実施方法
プレート間固着強度の時空間ゆらぎの同化の模擬テストを行うために必要な、観測データベースとシミュレーションデータベースを構築する。観測データベースの構築
として、GEONET の地殻変動データの整理と可視化によるモニタリングを行う。シミュレーションデータベースの構築は、地球シミュレータを用いた地震発生サイクルのフォワード計算によって行う。
フォワード計算では、プレート境界で発生する地震の繰り返しを、プレート相対運動からのずれの蓄積と解放過程としてモデル化する（e.g., Rice, 1993）。準動的仮定に基づくプレートの相対運動に起因するすべりによって生じるせん断応力と、強度の変化に伴うすべりの時空間変化を計算する。プレート境界面の摩擦は、すべり速度・状態依存摩擦則（Dieterich, 1979）に従うと仮定し、断層の構成則は Nakatani, 2001 に従うと仮定した。強度の時間発展則には、composite law（Kato and Tullis, 2001）または
slowness (aging) law（Dieterich, 1979; Ruina, 1983）を用いた。摩擦パラメータ A（=aσ）,
B（=bσ）, L（特徴的すべり量）について様々な空間分布を仮定して、フォワード計算を行った。モデル領域は、豊後水道スロースリップイベントを同化模擬テストの対象とした南海トラフ巨大地震震源域西縁（日向灘－紀伊半島沖）、および 2011 年東北地
方太平洋沖地震前後のゆっくりすべりを対象とした日本海溝地域（宮城県沖－茨城県沖）である。
306

さらに、ゆっくりすべりの観測事例や発生機構に関する海外の研究動向の情報収集を行うために、学会等に参加した。

(c) 業務の成果
2011 年東北地方太平洋沖地震の発生前後に観測されたゆっくりすべりによるプレー
ト間固着強度の時空間ゆらぎ（Ito et al., 2013; Uchida and Matsuzawa, 2013）の同化の模擬テストを行うために必要となる観測データベースの構築を、日本海溝地域について行った。海洋研究開発機構において、前プロジェクトから引き続きモニタリングを行っている南海トラフ地域と同様に、図３－９－②－１のような図を作成した。この図を、週１回、最新のデータを加えて自動更新するシステムに組み込み、地殻変動のモニタリン
グを開始した。
南海トラフ巨大地震震源域西縁（日向灘－紀伊半島沖）について、1946 年南海地震・
1968 年日向灘地震の震源域（八木ほか, 1998）と、豊後水道スロースリップイベントの震源域（e.g., Hirose et al., 1999）に周囲と異なる摩擦パラメータを図３－９－②－２のように仮定して地震発生サイクルの数値計算を行ったところ、南海地震震源域では、Mw8.2 の昭和タイプの地震が平均間隔 150 年で繰り返し発生した。そのうち、500～600
年に１度は、浅部や日向灘へすべり域が広がり、Mw8.7 の最大級の地震が発生した（図３－９－②－３）。これらの規模と繰り返しは、観測で得られている知見を概ね再現できている。豊後水道スロースリップイベントは、Mw6.9 のイベントが、繰り返し間隔
約８年で発生しており（図３－９－②－４）、観測された規模・繰り返し間隔で、再現できることを確認した。豊後水道では、規模の大きな(Mw8.7 の)南海地震発生後の約200 年間、スロースリップイベントは発生していなかった可能性を示唆する結果も得られた。このモデルでは、日向灘地震は Mw7.5 の地震が繰り返し発生したが、発生間隔
は、48 年から 209 年まで大きな変化を示していた。数通りのパラメタスタディを行った結果、日向灘地震の発生間隔（1707 年に日向灘地震の震源域が破壊していたとして、1968 年までの 260 年間）については、観測の知見を再現できるパラメタセットは得られなかった。
日本海溝地域（宮城県沖－茨城県沖）について、図３－９－②－５に示した摩擦パ
ラメータ分布を仮定して計算を行ったところ、スロースリップイベントが発生し、そのすべりが伝播して、Mw7.1 の地震（前震）が発生した。さらに、その余効すべりが伝播し、1.6 年後に Mw8.6 の地震（本震）をトリガーする、というシナリオが得られた（図３－９－②－６）。このシナリオは、定性的ではあるが、観測された事象を再現している。しかし、スロースリップイベントの継続期間やすべり速度、および前震と本震との時間差については、観測値との差が大きい。なお、M９地震のサイクル中に、宮城県沖
や茨城県沖で M７クラスの地震が 30～40 年の繰り返し間隔で発生する様子は再現できていた。
さらに、「スロー地震」研究集会（2013 年９月）や、AGU Fall Meeting（2013 年 12月）の「Interplay of Slow Slip, Tremor, and Earthquakes」などのセッションへ参加し、南海トラフと類似した条件の沈み込み帯である北米のカスカディア地域等で発生してい307

るゆっくりすべりの観測事例や発生機構に関する国内外の研究動向の情報収集を行った。

図３－９－②－１（左）日本海溝沿いの GEONET 観測点における上下動成分。（右）水

平成分の１年間の変位ベクトル。

図３－９－②－２

南海トラフ巨大地震震源域西縁（日向灘－紀伊半島沖）のモデルで仮定した摩擦パラメータ分布。NK、HG、BG は、それぞれ、図３－９－②－３
で黒実線、点線、灰色の線で示す時間変化を計算した位置を示す。

308

(a)

図３－９－②－３(a)南 海 地 震 震 源 域 と 日 向 灘 地 震 震 源 域 の す べ り 量 の 積 算 の 時間

変化（黒実線と点線）。灰色の線は豊後水道スロースリップイベント震源域におけるすべり速度の時間変化を示す。(b)Mw8.7 の南海地震のすべり分布。(c)日向灘地震のすべり分布。(d)Mw8.2 の南海地震のすべり分布。

図３－９－②－４

シミュレーションで得られた豊後水道スロースリップイベントの(a)すべり速度と累積すべり量の時間変化。(b)すべり分布。

309

図３－９－②－５

日本海溝地域（宮城県沖－茨城県沖）のモデルで仮定した摩擦パ

ラメータ分布。

図３－９－②－６

すべり速度分布のスナップショット。

(d) 結論ならびに今後の課題
2011 年東北地方太平洋沖地震前後のプレート境界すべりの同化の模擬テストを行うために必要な、地殻変動の観測データベースを構築した。シミュレーションデータベ
ース構築に向けて、南海トラフ巨大地震震源域西縁、および日本海溝地域において、地震発生サイクルの数値シミュレーションを行った。南海トラフ巨大地震震源域西縁で発生している豊後水道スロースリップイベントについては、現行のモデルで、予測試行実験に利用可能な精度で再現できた。今後は、過去のスロースリップイベント（1997 年，2003 年，2010 年）について GEONET で観測された地殻変動データと比較し、予測の試行実験を行うために、パラメタスタディを行い、シミュレーションデー
タベースを増やしていく。並行して、日向灘地震の繰り返しを歴史記録に基づく知見
310

と整合するモデルを得るために、パラメータ探索範囲を広げて、摩擦パラメータのチューニングを行うことが必要である。日本海溝地域については、2011 年東北地方太平洋沖地震の前に観測された、スロースリップイベント発生から本震に至る過程が、定性的にではあるが、再現できた。しかし、それぞれのイベントの規模や時間発展に関して、定量的に再現することが、予測精度の向上に向けた今後の課題である。

(e) 引用文献
1) Dieterich, J. H., Modeling of rock friction, 1. Experimental results and constitutive
equations, J. Geophys. Res., 84, B5, 2161-2168, 1979.
2) Hirose, H., K. Hirahara, F. Kimata, N. Fujii, and S. Miyazaki, A slow thrust slip event
following the two 1996 Hyuganada earthquakes beneath the Bungo Channel, southwest
Japan, Geophys. Res. Lett., 26(21), 3237-3240, 1999.
3) Ito Y, R. Hino, M. Kido, H. Fujimoto, Y. Osada, D. Inazu, Y. Ohta, T. Iinuma, M. Ohzono,
S. Miura, M. Mishina, K. Suzuki, T. Tsuji, and J. Ashi, Episodic slow slip events in theJapan subduction zone before the 2011 Tohoku-Oki earthquake, Tectonophysics, 600,
14-26, 2013.
4) Kato, N., and T. E. Tullis, A composite rate- and state-dependent law for rock friction,
Geophys. Res. Lett., 28(6), 1103-1106, 2001.
5) Nakatani, M., Conceptual and physical clarification of rate and state friction: Frictionalsliding as a thermally activated rheology, J. Geophys. Res., 106(B7), 13347-13380, 2001.
6) Rice, J. R., Spatio-temporal complexity of slip on a fault, J. Geophys. Res., 98 (B6),9885-9907, 1993.
7) Ruina, A., Slip instability and state variable friction laws, J. Geophys. Res. 88(B12),10359-10370, 1983.
8) Uchida, N., and T. Matsuzawa, Pre- and postseismic slow slip surrounding the 2011
Tohoku-oki earthquake rupture, Earth and Planetary Science Letters, 374, 81-91, 2013.
9) 八木勇治, 菊池正幸, 吉田真吾, 山中佳子, 1968 年 4 月 1 日, 日向灘地震(M JMA 7.5)の震源過程とその後の地震活動との比較, 地震, 2, 51, 139-148, 1998.

③ 地殻変動データ整備および東北地方太平洋沖地震の余効変動解析(a) 業務の要約
1) 東北地方太平洋沖地震の余効変動解析
2011 年東北地方太平洋沖地震（以下、「 東北沖地震」という。）時と地震後の GEONET
や海底地殻変動観測によって観測された余効地殻変動データを用いてプレート境界
面上の余効すべりの時空間分布を求めた。この余効変動解析は粘弾性応答を同時に考慮した余効すべりインバージョン解析を実施することで、粘弾性緩和に起因する余効変動を除いたプレート境界面上の余効すべりを推定し、余震活動や粘弾性の影響を評価した。

311

2) 地殻変動データの整備
a) 西南日本
過去の西南日本の水準測量データを整理し、地震サイクル全体にわたる時空間的な地殻変動の概要を整理することにより、西南日本のプレート間カップリングの時空間変化を推定するためのデータセットの準備をする。
b) 南西諸島
南西諸島沿いの琉球海溝からフィリピン海プレートが沈み込んでいるが、プレート間固着は弱いとされてきた。しかしながら、測地学的なアプローチからのプ
レート間の固着は明らかにされていないため、地震の長期評価もされていない。
そのため、喜界島での水準測量による傾斜速度の測定と横当島による GPS 観測を実施することで、プレート境界に直交する喜界島･奄美大島･横当島の測線にてプ
レート間の固着を評価するための準備を行った。

(b) 業務の実施方法
1) 東北地方太平洋沖地震の余効変動解析
2011 年３月 11 日に東北地方太平洋沖地震（Mw9.0）が発生した。この地震による地殻変動が GEONET や海底地殻変動観測点により観測されており、地震時には日本
列島が東向きに変動した（図３－９－③－１(d)）。地震後の変動は陸上の観測点では東向きの変動であるが、「宮城沖１」などの海底地殻変動観測点では西向きの変動を示している。これらの地殻変動時系列を説明するために、余効すべりと粘弾性応答の２つの現象を同時に考慮し、GEONET と海底地殻変動観測のデータを時空間
インバージョン解析することで、2011 年東北沖地震の地震時、地震後のすべりの時
空間分布を推定した。粘弾性応答関数は、３次元的なプレートの沈み込みの形状や
地殻・マントルなどの物性値を考慮した上で有限要素法を用いて作成した。本研究では東日本の 2600 × 1500 × 400 km の領域をモデル化し、上部マントルの粘性率はOhzono et al. (2012)を参考に 4.0 × 10 18 Pa・s を採用した（図３－９－③－１(a-c)）。

312図３－９－③−１(a) 赤い矩形は有限要素法にてモデル化した領域。(b) プ

レート境界面上に設定した小断層。(c) 図３－９－③－１(a)の AB ライン上の断面図。(d) 地震時に観測された地殻変動｡(e) 地震後に観測された陸上での地殻変動｡

2) 地殻変動データの整備
a) 西南日本
長期的地殻変動の重要な測量データの１つとして、全国繰り返し水準測量がある。この測量は、全国の主要道路に沿って約２km 毎に１等から３等の水準点を設置し、約 21000 ヶ所を 100 年以上前から繰り返している水準測量である。西南日
本においては、1886 年から 1899 年にかけて第１回水準測量が実施され、第２回水313

準測量は西南日本では 1927 年から開始され、南海･東南海地
震以前である 1943 年までに西南日本の大半の測量が終わっている。また、第３回水準測量は南海 ･東南 海地震 を契機として開始され、その後は約７年周期で繰り返し測量が実施されている。水準測量は 0.1
mm の 単 位 ま で 測 量 記 録 が 残されており、それらの網平均から求めた１km あたりの標準
偏差は年代によって異なるが、

図３－９－③−２

1.3～3.5 mm である。それらの図 ３ － ９ － ③ − ６ (a) に 示 さ れ た 路 線 図 の 南

水準測量結果の精査を行った。

海･東南海地震の変動 を 含む上下変動図。 (a)

図３－９－③−２は水準測量

オリジナルの記録。(b)異常データを取り除き結果の精読例を示す。

修正した記録。

b) 南西諸島
南西諸島の北部地方は大隅 ･ト カ ラ 列 島 と 奄 美 ･沖 縄 諸
島群から構成され、西側の沖縄トラフでは拡大が指摘されている。一方、南西諸島の東側の琉球海溝からフィリピン
海プレートが沈み込んでおり、複雑なテクトニック背景を持っている（図３－９－③−３）。
ほぼ一列に並ぶ南西諸島内において琉球海溝に直交した約100 km の測線を構築出来る場
所 は 横 当 島 ･奄 美 大 島 ･喜 界 島である。よって、横当島（無人島）に GPS 観測点を 2013 年９月に設置し、測量すること
で、沖縄トラフの拡大と琉球

図３ － ９－ ③−３

海溝からのフィリピン海プレック 背 景。 星 印は 横当 島（ 無 人島 ）、 黒 矢印 はートの沈み込みの影響の分離

フィリピン海プレートの速度、赤矢印は沖縄トを試みる。また、喜界島は南ラフの拡大を示す。
314

南 西 諸島 に おけ る テク トニ

西諸島内でも琉球海溝に最も近い島である。そのため、喜界島内の水準測量を実施することで、喜界島内の傾斜速度を測量し、プレート間カップリングを評価す
る基礎情報とする。

(c) 業務の成果
1)東北地方太平洋沖地震の余効変動解析
地殻変動データは GEONET と海底地殻変動観測点から得られた結果を用いた。地
震時のすべりは宮城沖の海溝付近に位置しており、最大すべり量は海溝付近で約 60m となった。粘弾性と弾性の両方の応答関数を用いたもので同様な結果が得られた（図３－９－③−４(a)と(a’)）。余効す
べりの時空間分布は、粘弾性を考慮する場合としない場
合では大きな違いが生じた（図３－９－ ③ − ４ (b-f) と(b’-f’)）。余効すべりの大きさは粘弾性を考慮した解析の場合は地震後 2.5
年間で最大で２mに達し、粘弾性を考慮しない場合では最大で４m に達した。また、余効すべりの分布は粘弾性を考慮しない場合に比べ、福島沖の余
効すべりは顕著に小さくなっており、
海底地殻変動観測によって観測され
図３－９－③−４
た、「宮城沖１」の弾性応答関数と粘弾性応答関数を用いたインバージョン結果の比較。(a－f)弾性応答関数、(a’－f ’)

西向きの変動も粘粘弾性応答関数を用いた結果。(a)地震時のすべり分布。(b弾性応答で説明可－f)は積算された余効すべり分布。(b)地震後 0.5 年。(c)地震能であった（図３－

後１年。(d)地震後 1.5 年。(e)地震後２年。(f)地震後 2.5 年。

９－③−５）。この事から上部マント315

ルの粘弾性応答に起因する地殻変動の影響は大きく、逆解析にて考慮する必要がある。 粘弾性を考慮した余効すべりの領域は地震時の破壊領域とは相補的な関係があり、過去のプレート境界での大地震の破壊領域とも相補的な関係がある（図３－９－③−６(a)）。なお、プレート境界で発生した地震の分布と余効すべりの分布を比較すると、余効すべりが発生している周りで地震が発生していると考えられる（図３－９－③−６(b)）。

図３－９－③−５海底地殻変動観

図３－９－③−６(a)地震時のすべり分測で得られた余効変動の 3 成分の地布と余効すべりの分布。主な地震時のすべ

殻変動時系列（青色）と粘弾性応答

りをハイライト系の色で示している。(b)を考慮した余効すべりモデルから計

地震後 0.5 年間の余効すべり分布とプレー

算された地殻変動時系列（赤色）｡ト境界の地震との比較。プレート境界地震

それぞれの観測点の位置は図３－９

は Asano et al. (2011)を用いた。

－③−１(e)を参照｡

316

2) 地殻変動データの整備
a) 西南日本
水準測量データを整理したことにより、既存の研究（Kobayashi and Yamamoto,2011; 小林,2012）で指摘されているスロースリップに起因する地殻変動を確認し、適切なデータの整理が出来ているかを確認した。図３－９－③−７は四国東部
における水準路線の上下動変動の時空間変動を示している。南海･東南海地震発生
後の余効すべりとスロースリップが測量されている
（図３－９－③−７(d,e)）。また、図３－９－③−８は西日本の水準測量
による比高差を空間分布として示している。これらの測量結果は広域にわたるため、全ての路線を網羅することは難しいのが現状である。また、測量路線も地域毎に分断されており、相互に路線を接続することができない。そのため、以前の研究は地域的な領域に限られ、個別研究として実施されてきた。広域にモデルで接続され、異なる種類のデータと異なる
図３－９－③−７

水 準 測量結果から各繰り 返 し測量間地域の解析を実施する必要がある。

の比高差。 (a) 水準路線図。 (b-f) それぞれの測量期間の比高差を表示。(c)以外は上下変位速度として表示（左側スケール）。(c) 南海･東南海地震の影響が含まれる測量期間であるため、変位として表示（右側スケール）。
317b) 南西諸島
南西諸島の横当島･奄美諸
島･喜界島近辺では GEONETにて地殻変動が観測されている。特に喜界島の上下変動については奄美大島に対して約２mm/yr の速度で高速隆起しており、地質学的に求められている隆起速度とほぼ一致す
る（図３－９－③−９）。一方、喜界島は傾動していることが地質学的手法により指摘され、長期的には傾動しているが、地震間でどのような傾動をしているのかを把握することは、プレート間のカップリングを調べる上でも重要である。そこで、我々は喜界島内の水準測量を実施し、1997 年９月に
国土地理院によって測量された２等水準路線を再測量することで、喜界島内の上下変動パターンを明らかにした（図３－９－③−１０）。その結果、
約 10 -7 /yr の速度で海溝側に傾動し沈降していることを示した。これらの結果は、プレー
ト間カップリングがある事を示唆している。

図３－９－③−８

繰り返し水準測量による、

広域的な上下変動パターン。それぞれの測量
期間の比高差の空間パターンを表示。測量期間についてはそれぞれの図の中に記載｡

318

図３－９－③−９

喜界島･奄美大島･横当島近郊の GEONET による 2002 年９月から 2012 年９月までの 10 年間の地殻変動。固定点は上対馬。星印と灰色の点は横当島と地震分布を示している。

図３－９－③−１０

喜界島内の水準測量結果。1997 年９月に実施された国

土地理院の水準測量と 2014 年３月に実施した水準測量結果の比高差。

319

(d) 結論ならびに今後の課題
1) 東北地方太平洋沖地震の余効変動解析
GEONET および海底地殻変動観測にて観測された 2011 年東北沖地震の余効変動か
ら、３次元構造を考慮した上でのアセノスフェアーの粘弾性応答を考慮して余効す
べりの時空間分布を求めた。その結果、粘弾性応答を考慮する場合としない場合では大きな余効すべりの時空間分布に大きな違いが確認できた。また、海底地殻変動
観測によって、観測された西向きの海底地殻変動も説明可能である。今後は注意深
く、余効すべりの推移の時空間分布を把握し、プレート境界面の摩擦特性などについて言及していく必要がある。

2)地殻変動データの整備
a) 西南日本
広域の水準測量結果をコンパイルし、その結果、地震サイクルにおいて南海･東南海地震の余効すべりだけでなく、スロースリップに起因するプレート間カップ
リングの揺らぎが確認できる。しかしながら、実際の水準測量の路線は分岐路線も多く有り、全てをコンパイルするにはかなりの時間を要する。そのため、コンパイルするためのアルゴリズムの必要がある。
b) 南西諸島
喜界島内での水準測量により、その結果、約 10 -7 /yr の速度で海溝側に傾動し沈降していることを示した。これらの結果は、定性的ではあるが、南西諸島の喜界
島近辺においてプレート間カップリングがある事を示唆している。今後、定量的な評価を実施する必要がある。

(e) 引用文献
1) Asano, Y., T. Saito, Y. Ito, K. Shiomi, H. Hirose, T. Matsumoto, S. Aoi, S. Hori, and S.
Sekiguchi, Spatial distribution and focal mechanisms of aftershocks of the 2011 off the
Pacific coast of Tohoku Earthquake, Earth Planets Space, 63, 669–673, 2011.
2) Kobayashi, A., and T. Yamamoto, Repetitive long-term slow slip events beneath the
Bungo Channel, southwestern Japan, identified from leveling and sea level data from1979 to 2008, J. Geophys. Res., 116, B04406, doi:10.1029/2010JB007822, 2011.
3) 小林昭夫, 高知市付近で 1978〜1980 年頃に発生した長期的スロースリップ, 地震
2, 64, 63-73, 2012.
4) Ohzono M., Y. Ohta, T. Iinuma, S. Miura, and J. Muto, Geodetic evidence of viscoelastic
relaxation after the 2008 Iwate-Miyagi Nairiku earthquake, Earth Planets Space, 64,
759-764, doi:10.5047/eps.2012.04.001, 2012.

④ 東北地方太平洋沖地震前後の地殻活動の特徴抽出(a) 業務の要約
東北地方の海陸地殻変動観測のデータを用いて、2011 年東北地方太平洋沖地震（以下、「東320

北地震」という。）の発生以前に発生したプレート間固着強度の時空間ゆらぎ（余効すべりを含むゆっくりすべりイベント、前震活動、余震活動）の予測実験に着手した。今年度は、過去の地震やゆっくりすべりの観測データの整理を進め、なゐふる 97 号にまとめた。また、観測データをシミュレーションデータベースと比較するための方法を検討した。本サブテーマに参加する研究者間で、観測データを共有するための VPN を用いたサーバーを構築した。さらに、
ゆっくりすべりの観測事例や発生機構に関する海外の研究動向の情報収集を行った。

(b) 業務の実施方法
東北地方の海陸地殻変動観測データおよび地震波形データを用いて、2011 年東北地震発生前に発生したプレート間固着強度の時間揺らぎの時空間的分布を精査することによって、シミュレーションによる再現実験で必要となる実観測結果のコンパイルを進めた。また、得られた観測データをどのようにシミュレーションデータベースと比較するのが適切かの検討を進めた。さらに本サブテーマ実施担当者の所属する東北大学および京都大学間で VPN を用いた回線を経由したデータミラーリングシステムの構築を行い、観測データの円滑な交換を可能とした。また米国地球物理学連合秋季大会等への参加により、ゆっくりすべりの観測事例や発生機構に関する海外の研究動向の情報収集を行った。

(c) 業務の成果
1) 数十年間のタイムスパンでのすべりの変化
東北地震前の数十年間の地震活動についてまとめた結果をここで述べる。この期間については、小繰り返し地震と GPS データから調べられている。図３−９−④−１に、東北地震のすべり域の内外について、小繰り返し地震から推定したスロースリップの推移を示す（Uchida and Matsuzawa, 2013）。将来地震が起こることになる大すべり
域内では、間欠的なすべりが 1998 年ころから見え、すべりの増加は 2008 年ころから顕著であった。長期間の GPS データからも、2004 年からいくつかの比較的大きな
地震とスロースリップが発生していたことが指摘されていた（Suito et al., 2011;
Ozawa et al., 2012）。これらのことは、東北地震のすべり域周辺で、固着が数年かけて緩んでいったことを示す。

321

図３−９−④−１ 小繰り返し地震データから得られた、東北地震の地震時大すべり域（左図
灰色；すべり量 10 m 以上）の内外でのプレート境界のスロースリップの推移（右図）（ Uchida
and Matsuzawa, 2013 を改変）。地震時すべり域内（右上）では、東北地震前にすべりが増加した。地震時すべり域外では定常的。東北地震後も対照的な時間変化を示す。

2) 数年間のタイムスパンでのすべりの変化
東北大学は、2011 年東北地方太平洋沖地震の震源域となった宮城県沖の海域で、2008 年からスロースリップの観測を行っていた（図３−９−④−２）。この観測では、海底に圧力計を設置し、海底面の隆起や沈降を測定する。観測は、2011 年東北地方
太平洋沖地震時まで実施され、本震発生前に特徴的な地殻変動を観測することができた（Ito et al., 2013; Hino et al., 2013）。この地殻変動は、2008 年 12 月と 2011 年１
月の２回観測された。2011 年の地殻変動は、2011 年１月下旬から観測され始めた。
2011 年の地殻変動は、１月 25 日頃から海溝軸に近い海底圧力計（GJT3-TJT1、
P08-P09）で変動が観測されはじめ、３月９日の最大前震直前まで継続した（図３−
９−④−３a）。また、同時期に陸上の体積ひずみ計（KNK、場所は図３−９−④−２を参照）でも変動が観測された（図３−９−④−３b）。これらの地殻変動を詳しく調べたと
ころ、Mw7.0 に相当するスロースリップが宮城県沖で発生していたことが分かった（図３−９−④−２）。2011 年のスロースリップが発生している間には、スロースリッ
プ域内で M５クラスのプレート境界型の地震のほか、小繰り返し地震も相次いで発生した（図３−９−④−２、図３−９−④−３）。海溝軸から最も近い観測点ペア（GJT3-TJT1）の記録から、スロースリップ域が、２月半ばから海溝側に範囲を拡大していたことを示している。そして、最終的にはスロースリップ域の西側の固着域、すなわち 2011 年３月９日 11 時 45 分 （日本時間）の M7.3 のプレート境界地震（最大前震）を誘発した可能性が高い。

322

図３−９−④−２ スロースリップ観測網（黄ひし形：海底圧力計、黄四角：体積ひずみ計）とそれにより得られた 2011 年のスロースリップ域（赤）。同時に、本震時の大すべり域（オレンジ）、本震（橙色星）と最大前震（濃青星）の震央の位置を示す。水色の星印は、2011
年２月に発生した M５以上の地震の震央の位置を示す。

323

図３−９−④−３ 海底圧力計記録(a)、体積ひずみ計記録(b)と海底圧力計周辺の地震活動度(c)。
海底圧力計記録は、隣り合った２つの観測点同士の圧力計記録の差（相対上下変動）を示す（Ito et al., 2013 を改変）。地震活動は、2010 年 12 月 17 日から 2011 年３月 11 日の本震まで M2.5 以上について示す。黒矢印はスロースリップの開始および海側にその範囲を開始した時期を示す。

3) １ヶ月～数日間のタイムスパンでのすべりの変化
Kato et al. (2012)は、東北地震前約１ヶ月間に発生した前震活動の詳細な時空間発展を推定し、その結果、本震発生の約１ヶ月前の２月中旬と、約２日前の最大前震（M7.3）の発生後の２度、本震の破壊開始点へ向かう震源移動現象がほぼ同じ領域（図３−９−④−４中の EMZ）で起きていたことを示した。それぞれの震源の移動速度
は１度目が２～５ km/日、２度目は平均約 10 km/日であった。最大前震後の地震活動度を見てみると、EMZ の北側では経過時間に反比例して地震活動度が単調に減少するのに対して（図３−９−④−４中の挿絵）、EMZ 内では１日程遅れて活動のピーク324を示し、EMZ 内外で異なる時間変化が観測された。
上記の震源移動を伴う前震活動には、小繰り返し地震が含まれる。小繰り返し地震に基づいてスロースリップの推移を分析することで、２度の震源移動に対応した
間欠的なスロースリップの伝播が認められた。これは、２度の震源移動は本震の破壊開始点へ向かうスロースリップの伝播を意味する。スロースリップの伝播が、本震の破壊開始点へ応力の集中を引き起こし、本震の発生を促した可能性が高い。

図３−９−④−４ 東北地震前の前震活動の時空間発展の様子（Kato et al., 2012 を改変）。青色の丸印は震源を表し、その大きさはマグニチュードに比例する。横軸は日付、縦軸は海溝軸に沿う距離を示す。赤い破線は震源移動のフロントの位置。黒色の星印：本震 M9.0の震源、黄色の星印：最大前震 M7.3 の震源、赤と緑の星印：小繰り返し地震の震源、EMZ：
震源の移動現象が見られた領域。左上の図中の破線は太平洋プレートの等深線、ピンクと青のコンターはそれぞれ 1978 年宮城県沖地震（Yamanaka and Kikuchi, 2004）、2005 年宮城県沖の地震（Yaginuma et al., 2006）のすべり域を示す。
4) 数日間のタイムスパンでのすべりの変化
2011 年３月９日の前震の震源は、1978 年宮城県沖地震（M7.4）、2005 年宮城県沖の地震（M7.2）の震源域（図３−９−④−４左上を参照）よりも沈み込む太平洋プレー
トの浅い部分（深さ 20km 程度）に位置する。最大前震の地殻変動は、海底水圧計および陸上 GPS 観測網で観測された。すなわち、2011 年３月９日の最大前震とそれに引き続く余効変動、そして M9.0 本震に伴う上下地殻変動を連続時系列データとして325

得られた（Ohta et al., 2012; Hino et al., 2013）。
図３−９−④−５に海底水圧計および陸上の GPS 観測点から推定された３月９日の最大前震の地震時すべり分布および、地震後から M9.0 本震までの間に発生した余効
すべり（スロースリップの一種で地震後に見られるもの）の分布を示す。図中に示した白い矢印が最大前震にともなう地殻変動を、同様に黒い矢印が地震後余効変動にともなう地殻変動をそれぞれ示す。これらの海底水圧計による上下変動場に着目すると、例えば GJT3 という観測点では地震時にはほとんど地殻変動が確認できない
ものの、地震後の余効変動では対照的に顕著な隆起が観測された。つまり、地震時と地震後でその空間パターンが大きく異なることが分かった。これはすなわち、地震時のすべりと地震後の余効すべりの発生している領域が異なることを示唆する。
これらの地殻変動データから推定された地震時（水色で塗り潰した部分）、地震後
余効すべり分布（緑色で塗り潰した部分）を見ると、両者のすべった場所が異なり、あまり重なり合っていないことが分かった。地震時すべりのマグニチュードは Mw7.2、
地震後余効すべり（スロースリップ）によって解放されたマグニチュードは Mw6.8であった。また、推定された地震後余効すべり分布の南端は東北地震の震央（赤星）のすぐ北側まで到達し、先に述べた震源分布から推定された非地震性すべりの伝播
とも調和的な結果となった。

326

図３−９−④−５ 海陸測地データから推定された 2011 年３月９日最大前震の地震時すべり
分布および地震後余効すべり分布（Ohta et al., 2012 を改変）。黄色のひし型が海底水圧計の場所を示し、上向きの矢印が隆起の地殻変動、下向きの矢印が沈降の地殻変動を示す。
白色の矢印は地震時の地殻変動、黒色の矢印は地震後から３月 11 日東北地方太平洋沖地震までに生じた地震後余効変動をそれぞれ示す。両矢印それぞれのスケールを図中左上に示す。

(d) 結論ならびに今後の課題
平成 25 年度は、東北地震後の３年間に行われた研究をレビューし、地震のおよそ 10年前から、直前にいたるまでの様々なスロースリップ現象についてまとめた（図３−９−④−６）。本震の 10 年程度前から、図の南側の細い赤四角の領域などで、すべりの増加が見られ始めた。その後 2008 年には、やや海溝に近い太い赤長方形の近く、および南と北の２つの細い赤四角でスロースリップが発生した。2011 年に入ると、太い赤長方形の場所で再びすべりが見られるとともに、２回の震源移動のうちの最初の移動が２月に見られた（白丸）。その後、2011 年３月９日の前震が水色の領域で発生し、その後の余効すべりが、薄緑色の領域や北側の細い赤四角内で発生した。また、この時期には本震の震央（橙色星）に向かうすべりの伝播を示すと考えられる２度目の震源の移動（黄丸）も観測された。そして最後に白コンターの場所を破壊した 2011 年３月 11 日の東北地震が発生した。
327

発生頻度が数百年から千年に１度といわれる巨大プレート境界地震に対して、これらの解析期間は非常に短く、スロースリップの長期的な振る舞いや大地震との関係など、地震の準備過程の全貌については、未だ解明されていない。しかし、今回、稀にしか発生しない巨大地震の前にスロースリップが起きていたことを捉えた点は、今後の地震発
生過程の研究にとって重要な知見である。地震の実態をより深く理解し、現実的な地震発生シミュレーションに貢献するために、スロースリップと通常の地震との相互作用など、次年度以降、更に詳しい分析を継続する。

図３−９−④−６ 東北地震前に震源域周辺でみられた様々なすべり。赤細線の四角領域でのすべりは小繰り返し地震によるもの、赤太線矩形領域でのすべりは、海底圧力、陸上ひずみデータによるもの、白と黄色の丸は、３月９日前と後の地震活動、青と緑の領域は、水圧計のデータおよび陸上 GPS データによる３月９日の前震の地震時すべりとその後の余効すべりを示す。図中の日付はそれぞれのイベントの発生時、橙色星は本震の震央、白いコンターは Iinuma et al., (2012)による本震のすべり分布を示す。
(e)引用文献
1) Hino , R., D. Inazu, Y. Ohta, Y. Ito, S. Suzuki, T. Iinuma, Y. Osada, M. Kido, H. Fujimoto,and Y. Kaneda, Was the 2011 Tohoku-Oki earthquake preceded by aseismic preslip ?
328

Examination of seafloor vertical deformation data near the epicenter, Mar. Geophys. Res.,
doi:10.1007/s11001-013-9208-2, 2013.
2) Iinuma, T., R. Hino, M. Kido, D. Inazu, Y. Osada, Y. Ito, M. Ohzono, H. Tsushima, S.
Suzuki, H. Fujimoto, and S. Miura, Coseismic slip distribution of the 2011 off the Pacific
Coast of Tohoku Earthquake (M9.0) refined by means of seafloor geodetic data, J.
Geophys. Res., 117, B07409, doi:10.1029/2012JB009186, 2012.
3) Ito, Y., R. Hino, M. Kido, H. Fujimoto, Y. Osada, D. Inazu, Y. Ohta, T. Iinuma, M.
Ohzono, S. Miura, M. Mishina, K. Suzuki, T. Tsuji, and J. Ashi, Episodic slow slip
events in the Japan subduction zone before the 2011 Tohoku-Oki earthquake,
Tectonophysics, 600(0), 14-26, 2013.
4) Kato, A., K. Obara, T. Igarashi，H. Tsuruoka, S. Nakagawa, and N. Hirata, Propagation of
Slow Slip Leading Up to the 2011 Mw 9.0 Tohoku-Oki Earthquake, Science, 335(6069),
705-708, 2012.
5) Ohta, Y., R. Hino, D. Inazu, M. Ohzono, Y. Ito, M. Mishina, T. Iinuma, J. Nakajima, Y.
Osada, K. Suzuki, H. Fujimoto, K. Tachibana, T. Demachi, and S. Miura, Geodetic
constraints on afterslip characteristics following the March 9, 2011, Sanriku-oki
earthquake, Japan, Geophys. Res. Lett, 39, L16304, doi:10.1029/2012GL052430, 2012.
6) Ozawa, S., T. Nishimura, H. Munekane, H. Suito, T. Kobayashi, M. Tobita, and T.
Imakiire, Preceding, coseismic, and post seismic slips of the 2011 Tohoku earthquake,Japan, J. Geophys. Res., 117, B07404, doi.org/10.1029/2011JB009120, 2012.
7) Suito, H., T. Nishimura, M. Tobita, T. Imakiire, and S. Ozawa, Interplatefaultslip along
theJapanTrench before the occurrence of the 2011off the Pacific coast of Tohoku
Earthquake as inferred from GPS data, Earth Planets Space, 63, 615–619, 2011.
8) Uchida, N., and T. Matsuzawa, Pre- and postseismic slow slip surrounding the 2011
Tohoku-oki earthquake rupture, Earth and Planetary Science Letters, 374(0), 81-91,
2013.
9) Yaginuma, T., T. Okada, Y. Yagi, T. Matsuzawa, N. Umino, and A. Hasegawa, Co-seismic
slip distribution of the 2005 off Miyagi earthquake (M7.2) estimated by inversion ofteleseismic and regional seismograms, Earth Planets Space, 58, 1549-1554, 2006.
10) Yamanaka, Y., and M. Kikuchi, Asperity map along the subduction zone in northeasternJapan inferred from regional seismic data, J. Geophys. Res., 109, B07307, doi:
10.1029/2003JB002683, 2004.
(3) 平成 26 年度業務計画案
前年度に引き続き、余効すべり・ゆっくり滑りデータ同化手法を開発し、八重山や豊後水道ゆっくり滑りの観測データを整理して、観測データをシミュレーションデータベースと比較するための研究開発を行う。また、東北地方の海陸地殻変動観測のデータを用いて、
平成２３年東北地方太平洋沖地震の発生前後に起こった、プレート間固着強度の時空間ゆらぎ（余効変動を含むゆっくり滑りイベント、前震活動、余震活動）の同化の模擬テスト
329

に着手する。さらに、ゆっくりすべりの観測事例や発生機構に関する海外の研究動向の情報収集を行う。

330

3.10

震源モデル構築・シナリオ研究(1) 業務の内容(a) 業務題目「震源モデル構築・シナリオ研究」(b) 担当者
所属機関

役職

氏名

東京大学情報学環教授

古村孝志

総合防災情報センター

特任助教

原田智也

東京大学地震研究所

准教授

市村強

教授

平原和朗

准教授

橋本千尋

独立行政法人防災科学技術研究所

総括主任研究員

福山英一

独立行政法人海洋研究開発機構

主任研究員

堀高峰

独立行政法人産業技術総合研究所

研究員

安藤亮輔

京都大学

教授

岩田知孝

防災研究所

准教授

関口春子

助教

浅野公之

東北大学

教授

今村文彦

災害科学国際研究所

助教

今井健太郎

関西大学

教授

高橋智幸

巨大地震津波災害予測研究センター
京都大学
大学院理学研究科
名古屋大学大学院
環境学研究科

社会安全学部(c) 業務の目的

南海トラフの過去地震の震源再解析に加え、国内外を含む他の地震発生帯での巨大地震の震源解析、並びに地震発生シミュレーションの結果等を参考にして、南海
トラフで発生する巨大地震の震源特性と連動様式の一般化（レシピ）を図る。また、
津波堆積物や津波石の移動を考慮した津波シミュレーション法を開発し、周辺諸国を含めて古文書等の記 録を集めることにより データが少ない過去地 震の震源過程を評価する。
強震動・津波シミュレーション法の大規模並列化を進め、高分解能・広帯域化するとともに、プレート詳細形状・物性モデルと高分解能地殻・堆積層モデルを結合した高分解能地下構造モデルを構築して、巨大地震シナリオの高度なハザード評価を行う。震源や地下構造モデルの不確定性と地震シナリオの不確実性（多様性）に伴う短周期強震動と長周期地震動の予測のバラツキを適切に評価し、防災に資することのできる実用的なハザード評価を行う。南海トラフ地震と南西諸島海溝地震の連動可能性や、相模トラフの地震や日本海溝の地震との最大連動の可能性、こうし331

た地震津波による広域津波について評価も行う。
(d) ８か年の年次実施業務の要約

平成 25～26 年度:
日本列島広域構造モデルを開発し、データに整合するようパラメータを推定した。過去の地震とその後の粘弾性応答を考慮した応力・強度分布の推定手法を構築した。津波に伴う土砂移動を評価できる津波浸水計算コードを開発した。
プレート形状と３次元地殻・堆積層構造モデルを整備して、近年の大地震の強震動・長周期地震動生成メカニズムを検証するとともに、地震動・津波シミュ
レーションを高精度化・広帯域化した。
平成 27～28 年度：
広域構造モデルを大地震発生前後の内陸地震の活発化問題に適用する。高精度化・広帯域化した強震動・津波シミュレーションと、震度・地殻変動・津波
データにより南海トラフの過去地震の震源過程を再評価するとともに、その結果を考慮した応力・強度分布推定結果にもとづく地震発生サイクルシミュレー
ションを行い、連動の多様性を含む地震発生シナリオを提案する。上記結果等を用いて南海トラフ地震の震源モデルのレシピ化を行い、平均及び最大クラスの震源モデルを明確化する。トラフ軸付近での地震発生域調査観測研究の結果を用いて、津波地震の発生と超巨大地震化可能性について、発生間隔（確率）を含めた現実的なハザード評価を行う。
平成 29～30 年度：
南西諸島海溝から南海トラフまでの３次元プレートモデルを構築して日本
列島広域構造モデルを完成させるとともに、南西諸島海溝の巨大地震の震源モデルを構築する。また、本地域での観測データと強震動・津波シミュレーショ
ンに基づき地震津波ハザードを評価する。また、他課題の成果や本研究が提供するシナリオに対するフィードバックに基づいて地震発生シナリオを改訂する。粘弾性を考慮した構造モデルでのデータ同化手法を開発する。
平成 31～32 年度：
本研究により整備された地下構造モデルとシナリオに基づき、南海トラフ地震および西南諸島海溝地震による全国の地震津波ハザード評価を行う。南海ト
ラフ地震と相模トラフ、日本海溝の地震との巨大連動や、西南諸島海溝地震との巨大連動の可能性と、これによる強震動と津波のハザード評価を行う。
(e) 平成 25 年度業務目的
日本列島広域構造モデルを開発し、データに整合するようパラメータを推定する。
過去の地震とその後の粘弾性応答を考慮した応力・強度分布の推定手法を構築する。津
波に伴う土砂移動を評価できる津波浸水計算コードを開発する。プレート形状と３次元
地殻・堆積層構造モデルを整備して、近年の大地震の強震動・長周期地震動生成メカニ
ズムを検証するとともに、地震動・津波シミュレーションを高精度化・広帯域化する。
332

(2) 平成 25 年度成果
① 日本列島粘弾性モデル構築とその応用(a) 業務の要約
最近の観測データを反映した高詳細な地殻モデルを用いた沈み込み帯の応力・強度分布の推定とそれに基づく地震発生予測を行うための、粘弾性を考慮可能な高詳細な日
本列島広域構造モデルの一次モデルを開発した。構造探査結果を随時取り入れて三次元
不均質構造モデルの改良を図ることができる仕様としている。なお、日本列島広域構造モデルを構築するための基礎検討として、西南日本・東北日本における高詳細な地殻モ
デルを用いた地殻変動解析を行い、その有効性及び列島モデル構築の上で留意点を確認した上で、上記のモデル開発を行っている。

(b) 業務の実施方法
本年度は、以下の 1）及び 2）の業務を行った。1）最新の観測データに基づき構築された高詳細な西南日本地殻モデルを用いて南海地震の想定断層滑りに対する地殻変
動解析を行った。また、高詳細な東北日本地殻モデルを用いた逆解析を行った。これらの解析から、高詳細な地殻モデルを用いた地殻変動解析の有効性と列島モデル構築のための留意点を確認した。これらを踏まえ、2）粘弾性を考慮可能な高詳細な日本列島広域構造モデルの一次モデルを開発した。

(c) 業務の成果
本年度の業務の成果を、(b) 業務の実施方法に示した 1）及び 2）にあわせて示す。
1）高詳細な地殻モデルを用いた地殻変動解析観測された地殻変動データと地殻変動解析を組み合わせた逆解析により、沈み込
み帯におけるプレート境界面の固着状態等の地殻の状態推定を試みる研究が多くなされている。本業務においても同様のアプローチにより、地殻の状態推定を行っている。
国土地理院の GEONET による陸域での密な地殻変動データ「国土地理院，2013」
や海洋研究開発機構の DONET 「海洋研究開発機構，2013」等海底での地殻変動データ等の観測データが整備されている。これらのデータは GEONET が mm から cm
オーダー、DONET が cm オーダーの誤差となっている。地殻の状態推定には、地上と海底で観測されたこれらの高精度な地殻変動観測データが今後用いられていくと考えられる。また、観測技術の発達により地殻データの蓄積も進んでいる。具体的には日 本列 島の 地殻 構 造物性 値デ ータ や分 解 能１km の地殻 形状 デ ータの 構築 が進められている「地震調査研究推進本部地震調査委員会，2013」、「 K. Koketsu et al., 2008」。
一方、地殻変動解析には、半無限均質（弾性）モデル等単純化された地殻構造モデ
ルでの解析解が主に使用されている。いくつかの最近の研究では日本列島規模の三次元有限要素モデルを用いた地殻変動解析が試みられ、地殻構造の形状平坦化や物333

性値の均質化等が解析結果に有意な影響を及ぼすことが示されてきた。そのため、できるだけ単純化のない地殻構造モデルを構築し、GEONET や DONET の地殻変動
観測データに見合うように、より高分解能・高精度での地殻変動解析を行うことが必要である。地殻構造の形状や物性値の不均質性を、解析解を用いた地殻変動解析
により扱うのは困難である。 これに代わる手法として、地殻構造の数値解析モデルを構築し、地殻変動の大規模な数値解析を行うことが考えられる。上述のように、数値解析に際しては近似や単純化を極力行わず、地殻構造の地表面・層同士の境界面の複雑な形状を、入手可能な日本列島の地殻構造データとできるだけ同じ分解能で再現した数値解析モデルを用いることが望ましい。ここではこのようにして生成される数値解析用モデルを「高詳細」なモデルと呼ぶことにする。既往研究「K. Hughes
et al., 2010」で用いられている三次元有限要素モデルは前述の日本列島地殻データと比べて分解能が低く、自由度は 10 6 オーダー程度である。一方、高詳細なデータを用いて三次元有限要素モデル生成を試みる場合、モデルの自由度は 10 8 オーダーとなるので、モデル生成・計算ともに大きな計算負荷を克服するための新たな手法が必要となる。そのため、高詳細モデルを用いた解析を行うとともに、必要に応じてハイパフォーマンスコンピューティングの発達を踏まえた解析環境の整備も行う必要がある。なお、下記及び付随する検討から、地殻形状と物性の不均質性の両者が解析結果に影響すること、その影響はモデルやモデルに入力する境界条件により異なることが分かっている。地殻変動解析の信頼性を高めるためには、高詳細な地殻モデルを用いることが望ましいことが示唆されている。
高詳細モデルの有効性を確認するため、高詳細な西南日本地殻モデルを用いて、南海地震の想定断層すべりによる地殻変動解析を行い、その結果を従来の手法による結果と比較する。具体的には、高詳細な西南日本地殻モデルデータを用いた三次元
有限要素モデルを生成し、これに南海トラフでの地震サイクルシミュレーションから得られた南海地震の想定断層すべりを入力し、西南日本における地殻変動解析を行った。モデル化範囲は東西 976 km、南北 784 km、鉛直 400 km である。地殻構造
モデルは４ 層構造とし、同様の問題設定を扱っている既往研究を参考に定めた各層の物性値を用いた。有限要素モデル生成パラメータには数値検証において収束が確認された値を用いた。生成したモデルを図３－10－①－１ に示す。複雑な地表面・地層境界面形状が細かく再現されていることがわかる。

334

976km

784km

NS

4 UD

400km

EW

図３－10－①－１ 高詳細な西南日本地殻有限要素モデル。全体図（左図）とその左下端上端部の拡大図（右図）。拡大図では、有限要素モデルの要素も併せて可視化している。

図３−10−①−２ 高詳細地殻モデルと従来モデルを用いた地殻変動解析結果の比較。左図は高詳細モデルを用いた地殻変動（単位は m）、右図は高詳細モデルによる結果と従来モデルによる結果の相対差。

モ デ ル の 自 由 度 は 134,000,160 、 四 面 体 要 素 と 立 方 体 要 素 の 要 素 数 は そ れ ぞ れ
44,713,627、36,405,531 となった。比較対象として、地表面形状も物性の不均質性も考慮に入れないモデル（従来用いられる半無限媒体での解析解による評価に相当）も生成した。これら二つのモデルで地殻変動解析を行い、地表面変位を比較する。
335

高 詳 細 モ デ ル と 従 来 法 に よ る 地 表 面 変 位 の 絶 対 値 の 比 較 を 図 ３ － 10－ ① － ２ に 示す。地表面変位分布域のうち、トラフ軸近くの縁部において相対差が大きくなっていることがわかる。図３－10－①－２右図中の赤点は現在稼働している DONET 観測点を示しており、DONET が相対差の大きい地点の地殻変動データを観測することがわかる。DONET の観測点のうち、相対差が大きい領域にある観測点での相対差は45％にも達する。このことは、DONET で観測した高精度な地殻変動データと地殻変動解析を用いた逆解析により地殻の状態推定を行う場合、地殻変動解析に高詳細な地殻構造モデルを用いることで、推定結果が大きく変わる可能性があることを示唆する。半無限弾性体等の単純化された地殻構造モデルに基づく従来法ではなく、本手法を用いることで、地殻の状態推定の信頼性改善が期待される。以上の成果は、引用文献「縣亮一郎・他，2013a」にまとめられている。
高詳細モデルの有効性を確認することを目的として、高詳細な東北日本地殻モデルを用いた断層滑りの逆解析も行った。まず、高詳細な東北日本地殻モデルから図３－10－①－３左図に示す東西 784 km、南北 976 km、鉛直 400 km の大きさの有限要素モデルを構築した。

(NS)

(UD)

(EW)

図３－10－①－３ 高詳細な東北日本地殻有限要素モデル（左図）とこのモデル領域の地表面に仮想的に配置した観測点（右図）。

数値検証で検証された解が収束すると期待される大きさの有限要素で離散化している。高詳細な地殻モデルの有効性を検証するため、現実的な設定を反映した仮想的な問題設定を行う。即ち、図３－10－①－３右図に示す様に観測点配置計画を踏
まえて海域にも観測点を配置し、「W. Suzuki et al., 2011」の断層滑りを入力し観測点で得られた地殻変動量を用いて断層滑りの逆解析を行った。比較として、地表面形
状及び物性の不均質性を無視したモデルを用いた逆解析（従来の半無限弾性界を用いた逆解析に相当）も行う。図３－10－①－４左図に逆解析結果を示す。この結果
336

は入力とほぼ変わらない結果となっている。仮想的に配置した観測点が多数断層面上を覆っているため、逆問題としては、悪くない性質となったため、良好に推定されている。図３－10－①－４右図に従来法との差を示す。西南日本での解析で示したように、地表面形状の無視及び物性の均質化は順解析結果にも大きな影響が生じる。順解析で得られるグリーン関数を用いて行う逆解析では、その影響を受けることとなり、結果として断層滑りの評価に大きな差を生じる可能性が示された。最新の観測データに基づく高詳細な地殻モデルに基づく解析の有用性が示唆された。以
上の成果は引用文献「縣亮一郎・他，2013b」、「R. Agata et al., 2013a」、「R. Agata et al.,2013b」、「R. Agata et al., 2014」にまとめられている。

図３－10－①－４ 高詳細地殻モデルを用いた逆解析結果（左図）と高詳細モデルによる結果と従来モデルによる結果の相対差（右図）。

上記の高詳細地殻モデルを用いた地殻変動解析では、複雑な地殻の幾何形状を取り入れることが可能な点、自由境界面を解析的に取り扱える点から、有限要素法を用いている。有限要素法による coseismic な地殻変動解析で十分解が収束する高詳細な西南日本地殻有限要素モデルや東北日本地殻有限要素モデルの自由度は１億を超える。このような大規模有限要素モデルの構築及びその求解には特段の手法が必要となる。本業務では、「T. Ichimura et al., 2013」に示す方法を東北日本、南西日本地
殻モデル用に修正して用いた。以下に、その概要を示す。数日程度の短い時間スケールを対象とするため、地殻を線形弾性体とみなし、食い違い断層に対する地殻変
動の静弾性解析を行う。境界条件を考慮して整理すると、問題は Ku = f という一次337

連立方程式に帰着する。ここで、K、u、f はそれぞれ全体剛性マトリクス、変位ベ
クトル、外力ベクトルである。断層面上での変位の食い違いは、節点分割法「H. J.
Melosh et al., 1981」によりモデル化する。節点分割法では断層面上の節点を仮想的に上下に分割し、節点の上下に逆向きの仮想変位に相当する節点荷重をそれぞれ与えることにより、食い違い断層を簡便に表現する。モデルの側面と底面には無限要素「O. C. Zienkiewicz et al., 1983」を生成し、半無限性を導入する。従来の CAD モデルを介した三次元有限要素モデル生成手法では、手作業等による調整を施すことでモデルのメッシュ品質が保証される。だが対象領域が大きく、自由度が 10 8 オーダ
ーの高分解能なモデル生成を行わなければならない場合、この手法ではコストが膨大となり、多くの場合モデル生成自体が困難となる。本手法は複雑な地表面・境界
面形状を持った成層構造の地殻データ（数値標高モデルデータと物性値）を入力とし、バックグラウンド構造格子を用いた三次元有限要素モデルの自動生成を行う。
これにより、地殻の複雑な形状や物性の不均質性を考慮した高詳細な地殻構造モデ
ル を 、 メ ッ シ ュ 品 質 を 確 保 し つ つ 高 速 に 生 成 す る こ と が 可 能 に な る 。 こ れ は 「 T.
Ichimura et al., 2009」の手法の拡張である。以下にモデル生成の手順を簡潔に述べる。
(i) 要求される分解能 ds に格子間隔を設定したバックグラウンド構造格子を、モデル化対象領域をカバーするように被せる。構造格子の生成領域は x、y、z 方向についてそれぞれ lx × 2kds、ly ×2kds、lz × 2kds と設定する。 ここで lx、ly、lz は対象領
域がすべてカバーされるように適当な自然数を選ぶ。k は八分木の深さである。(ii)
生成した各構造格子単位で独立に要素を生成する。地表面・境界面が格子と交差する場合、地表面・境界面と格子の辺の４つの交点とその格子の８つの頂点を用いて四面体要素を生成する。要素生成にはデローニ分割を用いる。交差しない場合は、格子はそのまま立方体要素として扱う。この手法により整合性のとれた要素が格子単位で独立に生成される。４つの交点と 8 つの頂点の距離がある閾値より小さい場
合は、アスペクト比の悪い要素が生成される可能性があるので、地表面・境界面の形状の近似を行う。上述のように、四面体要素と立方体要素は近似された地表面・境界面とバックグラウンド構造格子に基づいて格子単位で生成される。(iii) 物性が均質な部分（立方体要素が生成されている部分）では、立方体要素のサイズを大きくし、数値解が収束する範囲で自由度を少なくすることを考える。(iv) 側面・底面付
近では要素を大きくしても解析精度が受ける影響は小さいので、(iii) と同様にして立方体要素を融合する。立方体要素を八分木構造に従って s 倍まで大きくする。(v)
サイズが変わった立方体要素の周辺ではメッシュが不整合となっているので、サイズの違う立方体要素同士をつなげる四面体要素をデローニ分割により生成する。また 半 無 限 性 の モ デ ル 化 の た め 、 無 限 要 素 を モ デ ル の 側 面 と 底 面 に 生 成 す る 。「 T.
Ichimura et al., 2009」からの拡張点は地殻形状の近似精度の改善、サイズの違う立方体要素の生成とその制御、無限要素による半無限性の導入である。なお用いる要素の種類は四面体一次要素、立方体一次要素、三角柱型無限要素の二次要素である。
本手法の機能はあくまでも成層構造に対するメッシュ生成に限定される。だが入力の地殻構造データが成層構造として与えられていれば、要素が格子単位で生成され338

るために、対象領域の大きさに関わらず堅牢なメッシュ生成を可能とする。また節点の座標が格子点として明示的に与えられるので、メッシュの品質の制御が容易である。これにより任意の領域、分解能に対して高品質なメッシュを生成することができる。以上の有限要素モデル生成の手順を図３－10－①－５に示す。次に求解手
法について説明する。本手法によって生成されたモデルは、地表面・境界面形状を高い分解能で詳細に再現するため、自由度が多くなっている。計算時間を実用可能な範囲にするためには計算の並列化が不可欠である。

図３－10－①－５ 有限要素モデル生成の概念図。

また将来的に断層の状態推定や地殻データのあいまいさを定量的に扱う解析を行う場合には、グリーン関数計算やモンテカルロ的アプローチ等により多数ケースの計算が必要となる。このことを踏まえ本手法では、複数の計算ノードで多数回計算を行うことを前提とする。その際、分散メモリ型並列化により「複数の計算ノード
を用いた１ケースの計算」を繰り返すのではなく、共有メモリ型並列化により「個々の計算ノードを用いた多数ケースの独立な解析」が可能な手法を開発する。この手法では計算ノード間の通信によるオーバーヘッドがほとんどなく、トータルの計算時間は短くなると期待される。並列化には OpenMP「OpenMP API specification forparallel programming, 2013」を用いている。対象問題の自由度が大きい場合、計算資源の制約上 K マトリクスをそのまま計算メモリ上に記憶させることは難しい。そこ
で K マトリクスが疎行列であることを利用して使用メモリ量を抑える
Element-by-Element（EBE）法「T. J. R. Hughes et al., 1983」と反復法に基づく求解法を用いる。なお、求解法において、可変前処理、精度混合演算、マルチグリッド等の改良により、単なる EBE を用いた反復法に比べ、７倍程度高速に地殻変動解析が可能となっている。なお、参考文献「Okada, 1985」による半無限解を用いて数値検証を行っており、十分収束した解が得られていることも確認されている。

2) 粘弾性を考慮可能な日本列島広域構造一次モデル
1）において高詳細な地殻モデルを用いた地殻変動解析の有効性を示しつつ、留意
339

すべき事項を洗い出した。これらを踏まえて、日本列島全域を含む、東西 1700 km、
南北 2600 km、鉛直 400 km の領域に対し高詳細モデルを生成する（対象領域は図３－10－①－６参照）。
地表面形状のモデルとしては日本水路協会発刊の日本近海 30 秒（900 m）グリッ
ド分解能の水深データ「JTOPO30，2003」を用い、プレート間境界の形状には CAMP
スタンダードモデル「C. Hashimoto et al., 2004」を用いた。粘弾性層はユーラシアプ
レート及び北アメリカプレートの深さ 30 km から 60 km 部分に設定し、他の領域は全て弾性層としている。

図３－10－①－６ 日本列島広域構造モデル構築の対象領域。

Viscoelastic
layer

Crust layer

Philippine
Sea Plate

Pacific Plate

図３－10－①－７ 構築された日本列島広域構造一次モデルを用いて生成した３次元有限要素モデル。左図は鳥瞰図、右図は各層の形状。

340

狭い領域を対象とする場合、地球の曲率の影響は小さいが、今回のような大領域
でモデルを構築する場合には、曲率の影響は大きくなる。元のデータは球殻座標系
により定義されているため、これをそのまま直接直交座標系に変換し、モデル生成のための標高データとした。すなわち、直交座標系で、地表面が曲がったモデルが生成されることとなる。図３－10－①－７に生成された高詳細な日本列島広域構造
一次モデルを用いて生成した三次元有限要素モデルを示す。coseismic だけでなく、
postseismic な地殻変動を解析するために構築したモデルである。
(d) 結論ならびに今後の課題
本年度は、高詳細な地殻モデルを用いた地殻変動解析の有効性を示しつつ、留意すべき事項を洗い出した。また、これらを踏まえて日本列島広域構造一次モデルを構築した。
来年度以降、日本列島広域構造一次モデルを用いた coseismic と postseismic な地殻変動解析を行う予定である。その際、対象となるシミュレーションの自由度が膨大となるため、ハイパフォーマンスコンピューティング的な検討を加えつつ、解析環境の整備も併せて行うこととする。

(e) 引用文献
1) 縣亮一郎，市村強，兵藤守，堀高峰，平原和朗，堀宗朗, 震源断層の想定高度化に 関 す る 基 礎 研 究 -高 詳 細 な 三 次 元 地 殻 構 造 モ デ ル を 用 い た 地 殻 変 形 解 析 手 法の開発-, 土木学会論文集 A1（構造・地震工学），69，I767-I776, 2013a.
2) 縣亮一郎，市村強，兵藤守，堀高峰，平原和朗，堀宗朗 高 詳 細 な 三 次 元 モ デ ル
を用いた高速な地殻変形解析手法の開発とその東北地方太平洋沖地震断層すべり
推定への適用，計算工学講演会論文集，18，C-11-3

2013b.

3) Agata R., T. Ichimura, K. Hirahara, T. Hori, M. Hyodo and M. Hori, Geodetic Inversion
Analysis Method of Coseismic Slip Distribution, Using a Three-dimensional Finite
Element High-fidelity Model,

American Geophysical Union Fall 2013 Meeting,

DI31A-2202, San Francisco, 2013a.
4) Agata R., T. Ichimura, M. Hyodo, T. Hori, K. Hirahara, M. Hori Inversion

Analysis

Method of Coseismic Slip Distribution Using Three-Dimensional High-fidelity Model of
Crustal Structure, 5th Asia Pacific Congress on Computational Mechanics & 4th
International Symposium on Computational Mechanics, 1428, Singapore, 2013b.
5) Agata R., T. Ichimura, K. Hirahara, M. Hyodo, T. Hori and M. Hori, Several hundred
finite element analyses of an inversion of earthquake fault slip distribution using ahigh-fidelity model of the crustal structure, Procedia Computer Science (2014 Accepted).
6) Hashimoto, C., K. Fukui and M. Matsu’ura, 3-Dmodelling of plate interfaces andnumerical simulation of long-term crustal deformation in and around Japan, Pure appl.
Geophys., 161, 2053–2067, 2004
7) Hughes K., T. Masterlark and W. Mooney, Poroelastic stress-triggering of the 2005 M8.7
Nias earthquake by the 2004 M9.2 Sumatra-Andaman earthquake, Earth and Planetary341

Science Letters, 293, 289-299, 2010.
8) Hughes, T. J. R., I. Levit and J. Winget, An element-byelement solution algorithm forproblems of structural and solid mechanics, Computer Methods in Applied Mechanics
and Engineering, 36, 241-254, 1983.
9) Ichimura T., R. Agata, T. Hori, K. Hirahara, and M. Hori, Fast numerical simulation ofcrustal deformation using a Three-dimensional high-fidelity model, Geophysical JournalInternational, 195, 1730–1744, 2013.
10) Ichimura, T., M. Hori and J. Bielak, A hybrid multiresolution meshing technique for
finite element three-dimensional earthquake ground motion modeling in basins including
topography, Geophysical Journal International, 177, 1221-1232, 2009.
11) 地震調査研究推進本部地震調査委員会，「長周期地震動予測地図」2012 年試作版,
http://www.jishin.go.jp/main/chousa/12_choshuki/index.htm. （2013 年 3 月 1 日閲覧）．
12)

JTOPO30-

日 本 近 海30

秒 グ リ ッ ド 水 深 デ ー タ ，

http://www.mirc.jha.jp/products/finished/JTOPO30/，2003.
13) 海洋研究開発機構，DONET, http://www.jamstec.go.jp/donet/e/. （2013 年 3 月 1 日
閲覧）．
14) Koketsu, K., H. Miyake, H. Fujiwara and T. Hashimoto, Progress towards a Japan
integrated velocity structure model and long-period ground motion hazard map, The 14th
World Conference on Earthquake Engineering, China, 2008.
15) 国土 地理院，電子 基準点データ提供サ ー ビス, http://terras.gsi.go.jp/ja/index.html.
（2013 年 3 月 1 日閲覧）．
16) Melosh, H. J. and A. Rafesky, A simple and efficient method for introducing faults into
finite element computations, Bull. seism. Soc. Am., 71, 1391-1400, 1981.
17) OpenMP API specification for parallel programming, http://openmp.org/wp/. （2013 年
3 月 1 日閲覧）18) Okada, Y., Surface deformation due to shear and tensile faults in a half-space,
Bull.Seism.Soc.Am., 75, 1135-1154, 1985.
19) Suzuki, W., S. Aoi, H. Sekiguchi, and T. Kunugi, Rupture process of the 2011
Tohoku-Oki

mega-thrustearthquake(M9.0)

invertedfrom strong-motion

data.

Geophysical Research Letters, doi:10.1029/2011GL049136, L00G16, 2011.
20) Zienkiewicz, O. C., C. EMson and P. Bettess, A novel boundary infinite element,International Journal for Numerical Methods in Engineering, 19, 393-404, 1983.

② 地震発生サイクル計算手法の高度化(a) 業務の要約
本研究では、現状の地震サイクル計算で取り入れられていない要素の組み込みや計算の効率化に関する地震サイクル計算手法の高度化を行う。現状の地震サイクル計算では、均質半無限弾性体を仮定しているが、現実には地表は例えば海溝を含む海底地形など起伏に富んでおり、この形状が地震サイクルに及ぼす影響を調べる必要がある。そこ
342

で、今年度はまず、海底地形など任意の形状を有する自由表面（地表面）を持つ均質弾性体中での地震サイクル計算手法の開発を行なった。
全無限媒質中でのすべり応答関数を用いる、「S. Hok and E. Fukuyama，2011」に従う
方法で、すべり応答関数（ある断層セル i での単位すべりによる他の断層セル j における応力変化で通常 Kij と記す）の計算精度を保つには、より大きな領域の地表面をできるだけ小さなセルに分割して計算する必要があるが、計算に要するメモリ量や計算時間が膨大になる。そこで、精度と計算効率のバランスを取って、断層セルの深さにスケールして地表面領域サイズおよび分割セルサイズを決定する方法を提案した。
また、開発した手法を用いて、深さ７km に達する日本海溝に沈みこむ太平洋プレー
ト境界面で発生する地震について、海底地形が地震サイクルに及ぼす影響を調べた。海
溝の深さをフラットな地表面とする半無限媒質中での解析的すべり応答関数と比較す
ると、プレート境界からの距離が大きくなるためすべり応答関数が小さくなり、同じ摩擦特性を持つアスペリティを設定しても、海溝を有する海底地形を持つ場合の方が、地震間隔が伸びて全無限媒質の場合に近くなることが分かった「大谷真紀子・平原和朗,2013」。
なお、すべり応答関数の計算には通常の半無限媒質の解析解（ここでは屈曲したプ
レート境界を扱うのに適した三角形断層セルに対するすべり応答関数を「M. Comninou
and J. Dundurs，1975」に従い計算している）を用いるよりは計算時間がかかるが、一度階層型行列（H-matrices）の形「M. Ohtani et al., 2011」で計算しておくと、地震サイ
クル計算ではこれまでの計算の形をとりメモリ量や計算量は増加しないので、任意形状を有する地表面の場合でもこれまで同様、地震サイクル計算の省メモリ化・高速化が可能である。
次に、階層型行列（H-matrices）法を用いた動的破壊を含む地震サイクル計算法の開発中であるが、これについては、P 波や S 波などの位相に応じて時間依存して振幅が時間的に変化する動的すべり応答関数に、効率的に H-matrices 法を適用するかの概念設計を行なった。しかしながら、実際のコーディングには至っていないため、本報告ではこの業務の要約のみにとどめる。

(b) 業務の実施方法「S. Hok and E. Fukuyama，2011」は動的破壊シミュレーションの際、半無限均質弾性体中でのすべり応答関数の表現が複雑で計算効率の上で問題があるとして、全無限媒質中に断層とその面に働く応力がゼロとなる自由表面（地表面）を導入し、断層面セル
でのすべりにより生じる地表面での変位を計算し、断層セルのすべりによる応力と地表面で発生する変位による応力を、全無限媒質中でのすべり応答関数を用いて表すことにより、半無限均質弾性媒質中でのすべり応答関数を計算し、動的破壊過程をシミュレー
ションする新手法を提案している。
彼らの方法を適用して、任意形状を有する地表面を有する均質弾性媒質中での地震サイクル計算法を開発した。彼らの方法は、任意形状を有する地表面を持つ均質無限媒質でも適用可能であり、実際彼らは論文中にその指摘をしてはいるが、フラットな
343

地表面しか扱っていない。またこの方法では、有限のサイズを持つ地表面領域を小さなセルに分割して各地表面セルでの変位を計算する必要があるが、領域を広く取りセ
ルサイズをできるだけ小さくすると計算精度が上がる。しかしながら、計算に必要なメモリ量や計算時間を抑えるには、計算精度を保ちつつ、できるだけ小さなサイズの
地表面領域をできるだけ大きなセルサイズに分割する必要がある。彼らの論文ではこれらのことに触れられていない。そこで、図３－10－②－１に示すモデルを用いて、解析的表現のある平面自由境界の場合の半無限均質弾性体のすべり応答関数と比較して、計算精度を保ちつつ計算効率を上げるように、断層セルの深さおよび傾斜角に依存する地表面領域サイズと分割地表面セルサイズの断層セルの深さへのスケーリングを見出した。地震サイクル計算にはこのようにして一度滑り応答関数を計算しておけ
ば、通常の半無限媒質での場合と同様に、階層型行列（H-matrices）法による、すべり応答行列とすべり（速度）ベクトルの積の計算の省メモリ化・高速化が可能である。
以上開発した手法を、宮城県沖で７km の深さを持つ日本海溝に沈みこむ太平洋プレート境界面で発生する地震サイクルに適用して、海溝などの海底地形が地震サイクルに及ぼす影響を評価した。プレート境界面に働く摩擦力は、すべり速度と状態に依存
する摩擦則に従い、状態変数の時間発展は Composite 則「N. Kato and T. E. Tullis, 2001」に従うとして、放射ダンピング項「J. R. Rice, 1993」を用いる準動的地震サイクル計算を行なった。実際の宮城沖での海底地形を用いているが、問題を簡単化するため、海溝軸方向には、海底地形およびプレートの沈み込み形状は一様としている（図３－10－②－２）。このプレート境界面上に、図３－10－②－３に示す摩擦パラメータ（A-B、L）分布（すべり弱化摩擦特性）を持つ半径 18km のアスペリティを、深さ 20 km と40 km（本シミュレーションでは、深さ７km を Z=０としているので Z=13 km と 33 km）に設定した。このような設置条件のもとで、地震サイクルシミュレーションを実行し、海溝の深さを基準面とする半無限均質弾性体と海底地形を有する地表面の場合に生じる地震サイクルの違いを論じる。

図３－10－②－１

半無限媒質中の深さ 10 km に傾斜角 15°サイズ１km x １km の断層セルに１m の逆断層すべりを与えた場合。ここで、tsize は計算に使用する地表面領域サイズ、
csize は分割した地表面セルサイズを表す。
344

図３－10－②－２

宮城沖で深さ７km の日本海溝に沈み込む太平洋プレートの図。海溝軸方向（Y 軸）方向には海底地形およびプレート形状も一様とする簡単を行なっている。
青（real topography）は実際の海底地形を、緑(flat structure)は海溝深さを地表とする平面地
表面を表す。

図３－10－②－３

プレート境面上での摩擦パラメータ分布（アスペリティ形状）（左図）A-B 分布（MPa）（右図）L（特徴的長さ：これは一様としている）(c) 業務の成果
まず、傾斜角 15°サイズ１km x １km の矩形断層セルを深さ 10 km に置き、１m の逆断層すべりを与えた場合に、解析解と比較して、分割セルサイズ（csize）と地表面領域サイズ（tsize）をそれぞれ、３km および 90 km にとれば、十分な精度が保てることが分かった。
次に、深さを５km、10 km、20 km と変えた場合の地表変位の傾斜方向の変位と上下変位を、傾斜方向に沿って測った距離（X）を深さ（depth）でノーマライズして、（X/depth）
を横軸にして描いた図３－10－②－４を見ると、tsize の幅を取れば類似の形をしており、最大振幅でノーマライズすると同じ精度の変位となることが分かる。
以上から、csize/depth=0.3、tsize/depth=30 ととれば、精度が保てることが分かる。また傾斜角が小さい５°の場合でもこのスケーリングで十分な精度が保てることも分かっている。従って、このような断層セル深さによるスケーリングを行なって計算したすべ
り応答関数を解析解と比較すると十分な精度を持つことも分かっている。
345

次に、図３－10－②－２の海底地形を有する場合において、図３－10－②－３のア
スペリティを設定した場合の、地震サイクルに及ぼす影響について述べる。図３－10－②－５は、アスペリティの中心の深さが Z=13 km の場合の、中心におけるすべり速度の時間変化を比べたものである。赤線と緑線は、それぞれ全無限媒質(inf)と海溝深さに平面地表を持つ半無限媒質（half）における解析的すべり応答関数を用いた計算である。また、青線と桃線は、それぞれ新手法により計算された、平面地表(flat)および海底地形を含む地表を持つ均質媒質(real)による計算結果を示している。緑線（half）と青線（flat）はほとんど重なっており区別できず、ここで新手法での計算が解析解を用いたものと良い一致をしていることを示している。実際の海底地形を持つ場合（桃線：
real）は、海溝深さを地表面に設定した平面地表の場合(half、flat)に比べ、短くなり、全無限媒質の場合(inf)と平面地表面（half,flat）の場合の間に来ている。

図３－10－②－４

深さ 5 km、10 km、15 km での傾斜角 15°1 km x 1 km サイズの断層

セルでの１ｍの逆断層すべりによる地表変位（左：傾斜角方向の水平変位、右：上下変位）を傾斜方向の距離 X を断層深さ（depth）でノーマライズした X/depth を横軸にプロット。tsize の幅を取れば十分な精度で変位が表されている。

346

図３－10－②－５

アスペリティの深さ Z=13km の時の、アスペリティ中心セルにおける、すべり速度の時間変化。赤線と緑線は、それぞれ全無限媒質(inf)と海溝深さに平面地表を持つ半無限媒質（half）における解析的すべり応答関数を用いた計算である。
また、青線と桃線は、それぞれここで開発された手法により計算された、平面地表(flat)
および海底地形を含む地表を持つ均質媒質(real)における計算結果を示す。
それぞれの繰り返し間隔を、T inf 、T half 、T flat 、T real とすると、T inf =34.86 年、T half =38.65
年、T flat =38.55 年 、T real =36.86 年と なる 。これ を無限 媒質 からの ずれ の比で 見る と、(T flat -T inf )/(T real -T inf )=0.97 と な る 。 ま た 、 ア ス ペ リ テ ィ の 深 さ Z=33 km の 場 合 は 、T inf =33.33 年、T half =33.98 年、T flat =33.98 年、T real=33.73 年で、(T flat -T inf )/(T real -T inf )=0.98となる。
最後に、すべり応答関数 Kij（ソース断層セル j に単位すべりを与えた時のレシーバ
断層セル i に於ける）応力変化）の深さ依存性を示す。すなわち、ソース断層セルｊの深さを Z=0.02 km から Z=50 km まで変えた場合の、傾斜方向のプレート境界に沿ってレシーバセル i を取ったときの K ij の変化を示す。図３－10－②－６の赤線、緑線、および青線がそれぞれ半無限解析解（Kij half ）、新手法による平面地表面の場合（Kij flat ）、
および実際の海底地形の場合（K ij real ）の全無限媒質の解析解からの差（K ij inf ）を示す。
ソースセルが浅い場合のレシーバセルが浅い場合を除いて、Kij half と Kij flat はよく一致していて、新手法によるすべり応答関数は精度を持っていることを示している。また、ソースセルが深いとき（得に深さ 50km の時）、レシーバセルが浅いと精度が悪くなっている。これは地表近くのセルで地表の分割セルサイズが大きすぎて精度が落ちていることによる可能性がある。
ソースセルが地表に近づくにつれすべり応答関数の値は大きくなり、海底地形の効果も大きくなるのが見て取れる。このように、プレート境界の断層セルと地表の距離がすべり応答関数に影響を与え、地震サイクルに影響してくることになる。

347

図３－10－②－６

ソース断層セルの深さに対するプレート境界に沿ってのレシーバ

セル点でのすべり応答関数（全無限媒質の場合からの差）。縦軸は、すべり応答関数の差(MPa)、横軸は X(km)を示している。
348

(d) 結論ならびに今後の課題
海底地形などの任意形状を有する地表面を持つ均質弾性媒質中での、準動的地震サ
イクル計算の新手法を開発した。浅い部分での精度向上の問題は残されているが、計算効率を保ち、おおむね良好な精度を出せる地表領域サイズおよぶ分割セルサイズの深さへ依存を織り込んだ決定法を提案した。この新手法を日本海溝に沈み込む太平洋プレー
ト境界での地震サイクル計算に適用し、海底地形が地震サイクルに及ぼす影響を調べた
ところ、同じ摩擦パラメータでも地震の再来間隔が異なってくることが確認された。
今後は南海トラフなどの海底地形を考慮した実際的な地震サイクル計算へ適用して、その影響を見積もる必要があろう。また、より大きな領域、例えば、日本海溝沿いおよび南海トラフ沿いのプレート間地震サイクルを同時に解析領域に加える場合、これまでのやり方では、海溝の深さの違いが問題となったが、新手法を用いることで、日本海溝域と南海トラフ域に発生する巨大地震の相互作用を含む大規模地震サイクル計算も可能になると期待される。

(e) 引用文献
1) Comninou, M., and J. Dundurs, The angular dislocation in a half space, J. of Elasticity,5(3–4), 203–216, doi:10.1007/BF00126985, 1975.
2) Hok, S. and E. Fukuyama, A new BIEM for rupture dynamics in half-space and its
application to the 2008 Iwate-Miyagi Nairiku earthquake, Geophys. J. Int., 184, 301-324,
doi: 10.1111/j.1365-246X.2010.04835.x, 2011.
3) Kato, N., and T. E. Tullis, A composite rate- and state-dependent law for rock friction,
Geophys. Res. Lett., 28, 1103–1106, doi:10.1029/2000GL012060, 2001.
4) Ohtani M., K. Hirahara, Y. Takahashi, T. Hori, M. Hyodo, H. Nakashima, and T. Iwashita,
Fast computation of quasi-dynamic earthquake cycle simulation with HierarchicalMatrices, Procedia Com. Sci., 4, 1456–1465, 2011.
5) 大谷真紀子・平原和朗，地表面形状が地震発生サイクルに与える影響，日本地震学会 2013 年秋季大会，P1-45,2013.
6) Rice, J.R., Spatio-temporal complexity of slip on a fault, J. Geophys. Res., 98,
9885–9907, doi: 10.1029/93JB00191, 1993.

③ 南海トラフ地震発生サイクルの再現・モデル検証(a) 業務の要約
南海トラフ地震を対象とした地震発生サイクルの再現及びモデルの検証を目的とする研究を進める上で、数値シミュレーション・システムの大規模化・高解像度化と共に、多様な観測・データ解析を通したモニタリング情報をシミュレーションに取り込むことが重要である。今年度は先ず、「C. Hashimoto et al., 2013」の地震発生サイクルシミュレーション・システムを大規模汎用計算機に適用して、高解像度・広範囲のモデル領域で実行する為、計算コードのチューニングを行った。更に、高解像度モデルを用いて準静349

的テクトニックローディングと動的破壊伝播の連成シミュレーションのテスト計算を実施した。これらの業務を行なうに当たり必要な計算機環境を整備した。動的地震破壊
伝播シミュレーションは、テクトニックローディングを通して形成された剪断応力分布と断層構成関係パラメータ分布を用いて行なった。破壊は、応力集中域の周囲から開始
した後、徐々に伝播し、最終的に応力集中域全域に及んだ。この計算で得られたすべり分布を用いて、次のサイクルのテクトニックローディング・シミュレーションを継続することができる。
地震発生シナリオを生成する為には、過去のすべり履歴を整合的に再現し得る現実的な断層構成関係を設定することが重要な課題となる。そこで次に、その基本的な情報を得る為に、西南日本域の GPS 測地データについて、モデル領域を拡張して再解析を行ない、プレート境界面上のすべり遅れレートの分布を検討した。解析結果は、「橋本
千尋・他，2009」の結果と基本的に一致し、プレート境界面のすべり遅れが南海トラフに沿って帯状に分布することを示す。推定されたすべり遅れレート分布は、室戸岬付近に顕著なピークを持つ南海・東南海域に及ぶ大きなすべり遅れ域と、御前崎周辺に顕著なピークを持つ東海域のすべり遅れ域に分かれる。

(b) 業務の実施方法
南海トラフ地震を対象とした地震発生サイクルの再現及びモデルの検証を目的とする研究を進める上で、数値シミュレーション・システムの大規模化・高解像度化と共に、多様な観測・データ解析を通したモニタリング情報をシミュレーションに取り込むことが重要である。
「C. Hashimoto et al., 2013」が構築した地震発生サイクルシミュレーションシステム
は、三次元プレート境界面形状モデルを共通の基盤とした、準静的テクトニックローディングモデルと動的地震破壊伝播モデルから成る。プレート境界面上のすべり運動を規定する断層構成則は、「H. Aochi and M. Matsu'ura, 2002」のすべりと時間に依存する断層構成則を用いている。「C. Hashimoto et al., 2013」は、一例として、1968 年十勝沖地
震（M=8.2）の震源域を想定した二つのアスペリティに於ける地震発生サイクルシミュ
レーションを行い、地震間のテクトニックローディングを通して実現した応力状態が、それに続く動的破壊過程を支配していることを示した。この結果は、或る時点の応力状態を適切に推定することによって次ステップの地震発生シナリオを生成することが可能であることを示している。或る時点の応力状態を知る為には、過去のすべり履歴を再現する必要がある。従って、地殻変動データのインバージョン解析「C. Hashimoto et al.,2009; C. Hashimoto et al., 2012」のデータ解析を通してシミュレーション結果を検証することにより、過去のすべり履歴を整合的に再現し得る現実的な断層構成関係を設定することが重要な課題となる。更に、「E. Fukuyama et al., 2009」は、地震発生直前の応力状態と断層構成関係が与えられれば、その後の動的破壊伝播が計算され、そこから生じる地震波動も同時に計算できる事を示している。
以上を踏まえて、今年度は、名大では、先ず、「C. Hashimoto et al., 2013」の地震発生
サイクルシミュレーション・システムを大規模汎用計算機に適用する為のチューニング
350

を行なった。この際、作業を効率的に進める為に必要な計算機環境の整備を行なった。
次に、プレート境界面の摩擦特性の基本的な情報を得る為に、西南日本域の GPS 測地
データについて、「橋本千尋・他，2009」と同様の手法でモデル領域を拡張して再解析を行なった。更に、名古屋大学（以下、「名大」という。）と防災科学技術研究所（以下、「防災科研」という。）との共同で、高解像度モデルを用いて、準静的テクトニックロ
ーディングと動的破壊伝播の連成シミュレーションのテスト計算を行なった。防災科研
では、名大によって計算されたプレート境界上での剪断応力分布と断層構成関係パラメータ分布を用いた動的破壊伝播のシミュレーションが実行可能となるよう、計算機環境の整備を行った。動的破壊伝播の計算は、「S. Hok and E. Fukuyama, 2011」の方法により行ない、構築された計算機環境が十分機能している事を確認した。この際、断層破壊の nucleation zone については、アスペリティ領域の周囲から開始すると仮定し、その領域の剪断応力を降伏応力まで持ち上げる事により破壊伝播を開始させて計算を行なっている。

(c) 業務の成果
今年度の業務実施の結果、以下の成果を得た。
１）地震発生サイクルシミュレーション・システムを大規模汎用計算機に適用して、高解像度・広範囲のモデル領域で実行する為、整備された計算機環境の下で、計算
コードのチューニングを行なった。これにより、南海・東南海地震の震源域に対応する領域で、スプライン関数節点間隔 4 km，計算点間隔 2 km の高解像度の準静的テ
クトニックローディングのシミュレーションを実現した。

２）西南日本域の GPS 測地データについて、モデル領域を拡張して再解析を行ない、プレート境界面上のすべり遅れレートの分布を検討した（図３−10−③−１）。解析結果は、「橋本千尋・他，2009」の結果と基本的に一致し、プレート境界面のすべり遅れが南海トラフに沿って帯状に分布することを示す。推定されたすべり遅れレート
分布は、室戸岬付近に顕著なピークを持つ南海・東南海域に及ぶ大きなすべり遅れ
域と、御前崎周辺に顕著なピークを持つ東海域のすべり遅れ域に分かれる。また、「K.
Obara, 2006」によって得られた深部低周波微動の震源分布は、推定されたすべり遅れ
域の下限付近とよい一致を示す。

351

図３−10−③−１

西南日本域の GPS 測地データから推定したプレート境界のすべり遅れレート分布。
「橋本千尋・他，2009」と同様の手法でモデル領域を拡張して再解析を行なった。
青と赤のコンターは、それぞれ、すべり遅れとすべり過ぎを表わす。コンター間隔は、1.5cm/yr。緑色の星印は、1944 年東南海及び 1946 年南海地震の震央を表わす。緑色の領域は、1944 年東南海及び 1946 年南海地震の津波波源域を表わす「羽鳥徳太郎, 1974」。オレンジ
色の点は、深部低周波微動の震源を表わす「K. Obara, 2006」。

３）動的破壊伝播の計算を効率的に進める為に，計算機環境の整備を行なった。構築した計算機環境上で、地震発生直前の剪断応力分布と摩擦構成関係を用いて、南海トラフに於ける地震の動的破壊伝播のテスト計算を行なった。その結果、地震時の破壊伝播過程の高精度計算が可能であり、次のサイクル計算に必要なすべり分布を得ることが可能であることが分かった。

４）スプライン関数節点間隔４km，計算点間隔２km の高解像度モデルを用いて、準静的テクトニックローディングと動的破壊伝播の連成シミュレーションのテスト計算を実施した。テスト計算のモデル領域は、東南海地震の震源域付近とした（図３−10−③−２）。三次元プレート境界面形状は共通のモデルを用いている。準静的テク
トニックローディングモデルでは、モデル領域内のすべり分布を双三次スプライン
関数の重ね合わせで表現することにより離散化し、応力蓄積過程のシミュレーショ
ンを行なう。プレートの沈み込みに伴い周辺領域ではすべりが進行し、固着域では
352

すべり遅れによる応力蓄積が安定的に進行する。蓄積した応力が臨界状態になると、不安定な破壊が開始することになる。この時点で、準静的テクトニックローディン
グモデルから得られた剪断応力分布（図３−10−③−２）と断層構成関係パラメータ分
布を引き継ぎ、動的破壊伝播モデルによるシミュレーションを開始する。動的破壊
伝播モデルでは、モデル領域を計算点に対応した三角要素で離散化し（図３−10−③−３）、急激な応力解放過程のシミュレーションを行なう。断層破壊の nucleation zoneについては、アスペリティ領域の周囲から開始すると仮定し、その領域の剪断応力
を降伏応力まで持ち上げる事により破壊伝播を開始させた。動的破壊は、開始点から急激に加速して、応力集中域全域に及ぶ大地震に発展する（図３−10−③−４）。

図３−10−③−２

準静的テクトニックローディング－動的破壊伝播シミュレーションのテスト計算に用いたモデル領域とテクトニックローディングモデルにより計算した剪断応力
分布。準静的テクトニックローディングの計算に用いたプレート境界は青色で示されており、その中で動的破壊シミュレーションの計算に用いた領域は、赤線で囲まれた領域である。絶対位置の参考とするため日本列島の位置を重ね書きしてある。図は、北西から見た図になっており、下側がプレートが沈み込む方向となっている。

353

図３−10−③−３

動的破壊伝播の計算に用いた計算メッシュ。三角形要素のつなぎ合わせにより、プレート境界の形状を表現し、その上に、自由表面メッシュをおき、自由表面の効果を考慮した計算モデルとなっている。

図３−10−③−４

東南海地震の震源域付近をモデル領域とした動的破壊伝播シミュレーシ

ョン。上段はプレート境界における大地震によるすべり量の時間発展を示しており、下段
は、すべり速度の時間発展を示している。また、断層上部の水平面は、自由表面を表しており、自由表面上における変位及び変位速度を示している。
354

(d) 結論ならびに今後の課題
南海トラフに於ける過去の地震発生サイクルの再現、モデルの検証に向けて、シミュレーション・システムの大規模化・高解像度化の為のチューニングとテスト計算を行ない、今回構築した計算機環境が十分機能することを確認した。今後は、利用可能な計
算資源を最大限利用することのできる最適な計算規模を見積もり、計算資源を最大限に生かし、最大限の成果を生み出すような計算結果の導出をアレンジすることが課題である。地震発生サイクルシミュレーションでは、大地震直前の応力状態を推定し、その時点に於ける剪断応力分布と断層構成関係を用いて動的破壊伝播の計算を行ない、そこで得られた最終すべり分布を次のサイクルの計算の初期値として再び準静的テクトニッ
クローディングの計算を始めることになる。今後は、このサイクルの計算を進める事に
より、多様な地震発生シナリオの構築を目指す。更に、西南日本域の GPS 測地データについて、モデル領域を拡張して再解析を行ない、プレート境界面上のすべり遅れレー
トの分布を検討した。今後は、すべり遅れレート分布などの解析結果に基づき、より現
実的なプレート境界面の摩擦特性を推定し、それを用いた地震発生サイクルシミュレー
ションを実現することが課題となる。

(e) 引用文献
1) Aochi, H. and M. Matsu'ura, Slip-and time-dependent fault constitutive law and its
significance in earthquake generation cycles, Pure Appl. Geophys., 159, 2029-2046,
2002.
2) Fukuyama, E., R. Ando, C. Hashimoto, S. Aoi, and M. Matsu'ura, A physics-based
simulation of the 2003 Tokachi-oki, Japan, earthquake to predict strong ground motions,
Bull. Seismol. Soc. Am., 99, 3150-3171, doi:10.1785/0120080040, 2009.
3) Hashimoto, C., E. Fukuyama, and M. Matsu'ura, Physics-based 3-D simulation forearthquake generation cycles at plate interfaces in subduction zones, Pure Appl. Geophys.,doi: 10.1007/s00024-013-0716-4, 2013.
4) 橋本千尋・鷺谷威・松浦充宏，ＧＰＳデータインバージョンによる西南日本のプレート間カップリングの推定, 日本地震学会 2009 年度秋季大会予稿集，A32-08，
2009．
5) Hashimoto, C., A. Noda, T. Sagiya, and M. Matsu'ura, Interplate seismogenic zones along
the Kuril-Japan trench inferred from GPS data inversion, Nature Geoscience, 2, 141-144,
2009.
6) Hashimoto, C., A. Noda, and M. Matsu'ura, , The Mw9.0 northeast Japan earthquake:
total rupture of a basement asperity, Geophys. J. Int., 189, 1-5, 2012.
7) 羽鳥徳太郎，東海・南海道沖における大津波の波源―1944 年東南海、1946 年南海道津波波源の再検討と宝永・安政大津波の規模と波源域の推定―，地震 2，27，10-24,
1974.
8) Hok, S. and E. Fukuyama, A new BIEM for rupture dynamics in half-space and its
application to the 2008 Iwate-Miyagi Nairiku earthquake, Geophys. J. Int., 184, 301-324,355

doi: 10.1111/j.1365-246X.2010.04835.x, 2011.
9) Obara, K. and H. Hirose, Non-volcanic deep low-frequency tremors accompanying slow
slips in the southwest Japan subduction zone, Tectonophysics, 417, 33-55, 2006.

④ 地震発生シナリオ評価(a) 業務の要約
南海トラフにおける、過去地震の発生サイクルの再現とモデルの検証のための数値シミュレーションを実施する。平成 25 年度は、準動的地震発生サイクルの数値シミュ
レーションから予測される南海トラフ地震のシナリオの部分すべりを初期条件として、
津波伝播シミュレーションを実施し、津波シミュレーションによる予測値と、過去発生した最大規模の南海トラフ地震である、宝永地震について利用可能な太平洋沿岸・瀬戸
内海の被害データとの整合性を検証した。これによって、地震サイクルシミュレーショ
ンで再現すべき、最大規模の南海トラフ地震の規模・すべり分布を制約した。

(b) 業務の実施方法
地震発生サイクルの数値シミュレーションから予測される南海トラフ地震の様々なシナリオを初期条件として、津波伝播シミュレーションを実施することにより過去発生
した南海トラフ地震について利用可能なデータとの整合性を評価、過去地震の発生サイ
クルを再現するようなシナリオの絞り込みを行う。平成 25 年度はまず、実現可能性が高い南海トラフ地震の規模の上限・すべり分布をおさえることを目標として、準動的地
震発生サイクルシミュレーションから予測された最大規模の南海トラフ地震シナリオ（Mw〜9.0）「M. Hyodo and T. Hori, 2013」を対象とし、そのシナリオに起因する沿岸での津波高の予測値を津波伝播シミュレーションコード JAGURS 「T. Baba et al., 2014」
を利用し見積もった。これらの予測値と、歴史的に最大規模と考えられている 1707 年の宝永地震津波での被害との比較から、被害に整合する地震シナリオの規模・すべり分
布を導出する。津波高の比較には、従来から存在する九州―東海地方にかけての太平洋岸の津波高の観測データに加え、瀬戸内海の津波高データ「北原他, 2012」を考慮するなど、新たに収集されたり、最新の調査によって改訂されたデータも加え整合性を評価する。

(c) 業務の成果「M. Hyodo and T. Hori, 2013」による地震発生サイクルシミュレーションから予測された最大規模の南海トラフ地震シナリオ（Mw〜9.0、図３－10－④－１a）と、過去地震との整合性を検証するため、そのシナリオ自身及び、そのシナリオから抽出した部分
すべり３ケース（図３－10－④－２）に対する地震時地殻変動を初期条件として非線形長波近似による津波伝播シミュレーションコード JAGURS「T. Baba et al., 2014」を利用した津波伝播シミュレーションを実施し、各シナリオから予測される西南日本沿岸での津波高を見積もった。比較のために、従来考えられてきた最大規模（８クラス）南海
地震シナリオ（３－10－④－１b）についても津波高を見積もった。これら予測値と過356

去地震との整合性を検証するため、南海トラフにおける、最大規模の過去地震である
1707 年の宝永地震津波の被害データと津波伝播シミュレーションによる予測値との比較を行い、南海トラフにおける過去地震の規模の上限やすべり分布を制約した。
(b) Smaller-earthquake scenario (Mw = 8.66)
Vertical crustal deformation
at the surface(a) Larger-earthquake scenario (Mw = 9.03)

図３－10－④－1 （a）地震サイクルシミュレーションにより得られた最大規模（Mw = ９）の地震シナリオのすべり分布（緑線、コンター間隔は２m）と、そのすべりに起因する地殻上下変動（赤の領域が隆起、青の領域が沈降）。
（b）M８クラスの地震シナリオのすべり
分布と、そのすべりに起因する地殻上下変動。

(a) Sub-slip 1 (Seismogenic)(b) Sub-slip 2 (Up-dip + Seismogenic)

Vertical crustal deformation
at the surface(c) Sub-slip 3 (Seismogenic + Down-dip)

図３－10－④－２（a）最大規模（Mw = ９）の地震シナリオから日向灘と浅部・深部す

べりを除去した部分すべりの分布（緑線、コンター間隔は２m）と、その部分すべりに起因する地殻上下変動（赤の領域が隆起、青の領域が沈降）。
（b）最大規模（Mw = ９）の地
震シナリオから日向灘と深部すべりを除去した部分すべりの分布と、その部分すべりに起因する地殻上下変動。（c）最大規模（Mw = ９）の地震シナリオから日向灘と浅部すべりを除去した部分すべりの分布と、その部分すべりに起因する地殻上下変動。

まず、トラフ軸に沿う、西側への震源域の広がりを制約するため、九州地方の太平洋沿岸での被害（図３－10－④－３の橙の丸印の位置）と、５つの地震シナリオ（図３357

－10－④－１、図３－10－④－２）による予測値とを比較した （図３－10－④－４）。
その結果、四国より西側の日向灘まで 10 m 以上のすべりを伴うシナリオ（図３－10－④－１a）でのみ、被害データと同等の津波高が予測された。沿岸の津波高は近傍のプレート境界でのすべりの有無により大きな影響を受け、四国より西側でのすべりが 10
m に満たなかったり（図３－10－④－１b）、四国側のみが様々なすべり分布で 10 m
以上すべる場合（図３－10－④－２の３つのシナリオ）では、九州沿岸で予測される津波高は被害データに及ばず、観測と整合しない。つまり宝永地震では九州側まで大きなすべりが発生した可能性が高く、これは 「T. Furumura et al., 2011」と調和的な結果である。

図３－10－④－３

宝永地震津波の被害の分布。丸印は太平洋岸での被害地の位置、白抜

きの四角は瀬戸内海での被害地の位置を示す。

358

図３－10－④－４

九州地方太平洋岸での、地震シナリオを初期条件とした最大津波高の予測値（実線）と宝永地震津波被害の比較（橙色の丸印）。
（a）地震サイクルシミュレーシ
ョンから得られた M９クラス（Larger）と M８クラス（Smaller）の地震シナリオによる最大津波高の予測値（M９クラスが黒線で M８クラスが灰色線）と被害の比較。（b）M９ク
ラス地震シナリオの部分すべりによる最大津波高の予測値と被害の比較。

同様に、四国の太平洋沿岸での被害（図３－10－④－３のピンクの丸印の位置）と、5 つの地震シナリオによる予測値とを比較した。図３－10－④－５a、図３－10－④－
５b に示すように、トラフ軸までが 10 m 以上すべり、マグニチュードが９に達するシ
ナリオでは、いくつかの地点で、被害データを数倍程度上回るような予測となっている
ものの、そのシナリオから等深線に沿ってトラフ軸近傍の部分すべり域を除去したシナ
リオでは、被害データと同等の津波高となっている（図３－10－④－５b の赤線・緑線）。
このことは宝永地震時に四国沖ではトラフ軸まで達するような大きなすべりは発生しなかったことを示唆している。

図３－10－④－５

四国地方太平洋岸での、地震シナリオを初期条件とした最大津波高の予測値（実線）と宝永地震津波被害（ピンク色の丸印）の比較。（a）地震サイクルシミュ
レーションから得られた M９と M８クラスの地震シナリオによる最大津波高の予測値（M９クラスが黒線で M８クラスが灰色線）と被害の比較。（b）M９クラス地震シナリオの部
分すべりによる最大津波高の予測値と被害の比較。
359

これに対し、近畿地方の太平洋沿岸での被害（図３－10－④－３の緑の丸印の位置）と、５つの地震シナリオによる予測値とを比較すると、トラフ軸までが 10 m 以上すべり、マグニチュードが９に達するシナリオと被害データはかなりよく整合する（図３－10－④－６a）。部分すべりの予測値からは、近畿地方沿岸の津波高に対しトラフ軸近傍での大すべりが及ぼす影響は小さく、大すべりが発生しても被害データを逸脱する津波高を引き起こさないことが分かる（図３－10－④－６b）。このことは、『ちきゅう』に よ る 熊 野 灘 沖 の 掘 削 コ ア か ら 見 つ か っ た ト ラ フ 軸 近 傍 の 大 す べ り の 痕 跡 「 A.
Sakaguchi et al., 2011」と矛盾しない。また、すべりの深部延長の有無によって生じる
沿岸での地震時地殻変動によって、太平洋岸の津波高の予測値は大きく影響を受ける。
つまり、深部すべりの水平投影位置と近畿地方太平洋岸の位置関係により、深部すべりが生じない場合には、地震時すべりは太平洋岸の沈降を促進し、予測される津波高は被害データを超越してしまう。これに対し、深部すべりによる地震時地殻変動は、太平洋岸の沈降を抑制し、予測津波高を被害データと整合させるのに寄与している。

図３－10－④－６

近畿地方太平洋岸での、地震シナリオを初期条件とした最大津波高の予測値（実線）と宝永地震津波被害（緑色の丸印）の比較。（a）地震サイクルシミュレー
ションから得られた M９と M８クラスの地震シナリオによる最大津波高の予測値（M９ク
ラスが黒線で M８クラスが灰色線）と被害の比較。（b）M９クラス地震シナリオの部分すべりによる最大津波高の予測値と被害の比較。

伊勢湾から駿河湾にかけての被害データ（図３－10－④－３の青の丸印）について
も、近畿地方と同様に浅部すべりよりも、深部延長でのすべりの寄与が大きい（図３－10－④－７）。つまり、近畿地方以東に対しては、従来言われている地震発生帯（深さ方向の下限〜25 km）より深部（SSE の発生深度程）にまですべりが広がることが被害データとの整合性を高めるために必要である。

360

図３－10－④－７

東海地方太平洋岸での、地震シナリオを初期条件とした最大津波高の予測値（実線）と宝永地震津波被害（青色の丸印）の比較。（a）地震サイクルシミュレー
ションから得られた M９と M８クラスの地震シナリオによる最大津波高の予測値（M９ク
ラスが黒線で M８クラスが灰色線）と被害の比較。（b）M９クラス地震シナリオの部分すべりによる最大津波高の予測値と被害の比較。

近畿以西の被害データに対する、すべりの深部延長の影響は、太平洋岸の被害に対しては小さい（図３－10－④－５）。しかし、瀬戸内での被害（図３－10－④－３の白抜きの四角の位置）に対しては、被害地点とすべりの深部延長との水平位置関係の結果として、すべりの深部延長の有無により地震時地殻変動の上下変動パターンが大きく変化する。つまり、図３－10－④－８に示すように、深部延長でのすべりに起因する地震
時地殻変動によって瀬戸内海ではメートルオーダーの沈降を生じ、そこでの津波高の予
測値は、深部延長ですべりが発生しない場合より有意に大きくなる。このように、深部
延長でのすべりを含むシナリオは、中央部を除けば、瀬戸内海の宝永地震津波の被害データをよく説明することができる。

361

図３－10－④－８

瀬戸内海での、地震シナリオを初期条件とした最大津波高の予測値（実線）と宝永地震津波被害（黒色の丸印）の比較。（a）瀬戸内海において宝永地震津波によって津波被害が報告されている地点。
（b）地震サイクルシミュレーションから得られた M
９と M８クラスの地震シナリオによる最大津波高の予測値と被害の比較。（b）M９クラス
地震シナリオの部分すべりによる最大津波高の予測値と被害の比較。

(d) 結論ならびに今後の課題
平成 25 年度は、南海トラフにおける、最大規模の過去地震である宝永地震の規模・すべり分布を制約するため、準動的地震サイクルシミュレーションから得られた地震シ
ナリオの部分すべりを初期条件とした津波伝播シミュレーションを実施し、その予測値と宝永地震時の津波被害との整合性を評価した。その結果、宝永地震は、以下の特徴を持つことが分かった。
①震源域の西端は日向灘に達し、そこでは全体的に 10 m 以上のすべりが生じた。 ②
四国沖ではトラフ軸近傍での大すべりは被害に矛盾するものの、熊野灘以東ではトラフ軸近傍で大すべりが発生した可能性を否定できない。③瀬戸内や紀伊半島以東の太平洋岸の被害を説明するには SSE の震源付近にまで地震時すべりが生じた可能性が高い。
ただし、上で述べた宝永地震に整合するシナリオは、地震サイクルシミュレーショ
ン自体によって再現された地震シナリオでなく、地震シナリオから部分すべりを抽出す
362

ることによって構築したものである。このため、今後の課題として、地震サイクルシミュレーション自体により直接、今回得られた宝永地震に整合するすべり分布をもつ地震
シナリオを再現するとともに、安政・昭和の過去地震にも整合する地震サイクルを再現する必要がある。

(e) 引用文献
1) Baba, T., N. Takahash, Y. Kaneda, Y. Inazawa, & M. Kikkajin, Tsunami inundation
modeling of the 2011 Tohoku earthquake using three-dimensional building data for
Sendai, Miyagi Prefecture, Japan. In: Kontar, Y.A., -Fandino, V.S., Takahashi, T. (eds.),
Tsunami Events and Lessons Learned; Ecological and Social Significance, 89–98.
SPRINGER, 35, Dordrecht Heidelberg New York London, 2014.
2) Furumura, T., K. Imai, & T. Maeda, A revised tsunami source model for the 1707 Hoei
earthquake and simulation of tsunami inundation of Ryujin Lake, Kyushu, Japan, J.
Geophys. Res., 116, 1–17, 2011.
3) Hyodo, M., and T. Hori, Re-examination of possible great interplate earthquake scenariosin the Nankai Trough, southwest Japan, based on recent findings and numericalsimulations, Tectonophysics, 600, 175–186, 2013.
4) 北原糸子・松浦律子・木村玲欧(編)，日本歴史災害辞典，吉川弘文館，2012．
5) Sakaguchi, A., Chester, F., Curewitz, D., Fabbri, O., Goldsby, D., Kimura, G., Li, C.F.,
Masaki, Y., Screaton, E.J., Tsutsumi, A., Ujiie, K., and Yamaguchi, A., Seismic slip
propagation to the up-dip end of plate boundary subduction interface faults: vitrinitereflectance geothermometry on Integrated Ocean Drilling Program NanTroSEIZE cores.
Geology, 39, 395–398, 2011.

⑤ 巨大地震の震源モデル及び地殻・地盤モデル開発(a) 業務の要約
巨大地震時の強震動予測の高度化には、強震動予測のための巨大地震の震源モデル像の構築と、震源域から強震動予測を行う対象地点を含む地殻・地盤の地震波速度構造モデルの構築が必要である。強震動予測は地震被害の予測や軽減のために行うものである
ので、地震被害を受ける構造物の周期帯域を意識する必要がある。近年の構造物の巨大化や複雑化を受け、広帯域（具体的には周期 0.1 秒程度から 10 秒程度）の強震動予測が必要となってきている。この観点に基づいた震源モデルの複雑さや、地下構造モデルの複雑さを取り入れた検討を進めていく必要がある。
震源モデル及び地下構造モデルの構築にあたってはゼロからスタートする必要はない。政府の地震調査研究推進本部によって進められてきた地震動予測地図作成において、震源モデルや地下構造モデルは、広帯域の時刻歴強震動予測手法とともにプロトタイプが提案されて活用されている。また、震源モデルの広帯域化や各地域の、例えば大規模堆積盆地モデルといった地下速度構造モデルについては、個々の研究者による高精度化が行われている。本業務にあたっては、これらの震源モデル、地殻・地盤速度構造モデ
363

ルに関する既往研究に基づき、それらの広帯域化や信頼性を上げていくための方法論の検討と、改良モデルを提案していくことを目的としている。

(b) 業務の実施方法
強震動予測の高度化のための、巨大地震に対する震源モデルと地盤・地殻速度構造モデルの高度化を行う。今年度は、震源モデルの高度化に関しては、2011 年東北地方
太平洋沖地震の震源像について得られた知見に照らして既往の手法の改良の方向について検討を行い、地下速度構造モデルの高度化に関しては、既往速度構造モデルのひとつである全国一次地下構造モデル「K. Koketsu et al., 2008，2012」の妥当性の検討として、2011 年東北地方太平洋沖地震余震の震源域から大阪堆積盆地に至る長周期地震動
シミュレーションを行った。

(c) 業務の成果
1) 強震動予測のための巨大地震震源モデルの構築
海溝型巨大地震の広帯域震源モデル作成手法には、「関口春子・吉見雅行，2006」
や「N. Pulido et al., 2012」などがある。いずれも、長波長のすべりの不均質は、過去の地震の震源像、プレート境界面の幾何形状や固着域分布など、大きなスケールでしかわからないデータをベースに作成し、そこに、地震の統計的性質に則ってランダムな小さいスケールの不均質を付加して、広帯域で地震規模に見合う量の地震波を発生するようにしたものである。今後、「関口春子・吉見雅行、2006」をベースに改良していくことを念頭に、本節では、2011 年東北地方太平洋沖地震の震源像について得られた知見に照らして「関口春子・吉見雅行，2006」の手法を検討する。検討項目は、広帯域震源モデル作成手法自体、および、震源パラメータの変動による地震動のばらつきの見積り方法について、である。

・広帯域震源モデル作成手法の検討
まず、「関口春子・吉見雅行，2006」の広帯域震源モデル作成手法を概説する。
「関口春子・吉見雅行，2006」では、震源の長波長不均質モデルとして、既往のモデルを用いている。南海トラフ沿いの地震については、大きさは様々だが同じ
応力降下量を持つ複数のアスペリティと均一の背景領域からなる、中央防災会議の想定南海地震の震源モデル「中央防災会議，2003」を選んだ。長波長不均質震源モデルの広帯域化は、長波長不均質モデルの最小のアスペリティより小さいパッチを断層面上に均一かつランダムに配置し、そのパッチ内のすべり量および破壊伝播速度を上下させることによって行う。パッチのサイズとその数、およびすべり量に与える変動量は、フラクタル的な分布になるように与え、すべり分布は、最後にその波数スペクトルを過去の地震のすべり分布モデルの経験的波数スペク
トル形状 k-1.75 「M. P. Mai and G. C. Beroza, 2002」を持つよう調整する（図３−10−⑤−１）。破壊伝播速度の変動量は、震源インバージョン解から推定される変動量「宮腰研・A. Petukhin, 2005」と同等になるような定数値を与えている。なお、す364

べり速度時間関数は、Kostrov 型の形状を持つ「中村洋光・宮武隆，2000」を用いている。このような広帯域化作業により、震源モデルの震源スペクトルは ω -2 形状に近い形に改良される。

図３−10−⑤−１「中央防災会議，2003」の震源モデルを元に作成した南海地震広帯域震源モデルのすべり分布「関口春子・吉見雅行, 2006」。

図３−10－⑤－２

2011 東北地方太平洋沖地震の 0.1～10 秒の地震動に基づく強震動生成域（青枠）と、10 秒以上の地震動に基づくすべり分布（断層面上のカラートーン）「 K. Asanoet al., 2014」。
365以下、2011 年東北地方太平洋沖地震で得られた知見から、この手法を見直してみる。
東北地方太平洋沖地震については、解析対象周期の異なる多数の解析より、プレート境界面上で深さに応じて破壊過程の性格が大きく異なっていた可能性が示唆されている。10 秒以下の周期の地震波形では破壊開始点あたりから断層面最下部の深さ 50km あたりまでに主破壊域が複数個点在する様子が見られ（「K. Asano and
T. Iwata, 2012」など）、10 秒以上の周期の地震波形や津波波形からは破壊開始点あたりから断層面最上端にかけて広大な１つの主破壊域が広がっている様子が見られる（「W. Suzuki et al., 2011」など）（図３−10－⑤－２）。このような知見に基づき、プレート境界面上の深さ方向にすべりの特徴の異なる領域に分けられるという考えが提案されている（例えば、「T. Lay et al., 2012」など）。
地震動被害に関わる周期を用いて求められる強震動生成域（SMGA）も 10 秒以下の周期の解析に含まれるが、「K. Asano and T. Iwata, 2012」の強震動生成域モデ
ルによると、破壊開始点付近の強震動生成域は長周期の大すべり域に近接、もしくは重なるが、南部の強震動生成域は長周期で見たすべり量のあまり大きくない部分にあり、また、応力降下量とライズタイムの強震動生成域のサイズに対する比は、破壊開始点付近の強震動生成域が南部のものに比べ大きい傾向が見られる。
２つの領域で強震動生成域の破壊様式が大きく異なっていたことが示唆されるが、この違いを深さだけで区分するのも難しい。強震動生成域の特徴に幅があることがわかるが、これがプレート境界面のローカルな性質による固定した特徴なのか、
地震毎に替わりうる変動のうちなのかを判定するのは困難と思われる。なお、「K.
Asano and T. Iwata，2012」の強震動生成域モデルの、各強震動生成域の地震モーメ
ント量から見積もられるすべり量は、 「W. Suzuki et al., 2011」の震源モデルにおけるすべり量と同程度の値である。
アスペリティ像については、これまでの想定海溝型巨大震源モデルで想定されていたものと、東北地方太平洋沖地震のそれとは大きく異なる。これまでの想定
海溝型巨大震源モデルのアスペリティは、すべりの大きい場所でかつ強震動の発生源でもあった。すべり量で言うと、東北地方太平洋沖地震のアスペリティは、長周期の解析で見られた破壊開始点～海溝沿いの広大なひとかたまりの大すべり域ということになるが、ここは強震動生成域とは（重なっているとしても）一部
しか重なっておらず、そこでのすべりの大半は津波地震型のすべりであったと考えられる。大すべりを起こす場所としてアスペリティを設置し、そこを強震動の主要発生源とする従来の想定海溝型巨大震源モデルの組み立ては、この地震には当てはまらない。これまで海溝沿い浅部は津波地震を起こす場所として、別個に津波地震が設定され津波予測のみ行われこともあったが「中央防災会議, 2005」、
その考え方と同様、海溝沿い浅部の大すべりを別扱いするという考え方もできるかもしれない。しかし、今回の大すべり域はその津波地震型すべり域のみに収まらない可能性が高く、確認が必要である。ライズタイムの長い大すべりのエリアが典型的なプレート境界地震にも広がっている場合は、これまでのいわゆる、ア366

スペリティ-背景領域の枠組みの、背景領域の一部が地震によっては時定数の長い
大すべり域になると考えることもできるかもしれない。
強震動生成の観点からは、「K. Asano and T. Iwata, 2012」が示したように、複数のパッチからなる強震動生成域で東北地方太平洋沖地震の 0.1-10 秒の強震動はほぼ表現されている。この強震動生成域は、サイズや応力降下量の点で、これまでの海溝型巨大地震震源モデルのアスペリティと似通っている。しかし、前述のようにライズタイムの点でこれまでのアスペリティと違っている可能性もある。アスペリティ（強震動生成域）のこのバラエティを、基本的なスケーリングは変わらずにこれまでの考え方のアスペリティパラメータの変動幅の範疇としてよいか、
もしくは、今回のような津波地震型破壊域との同時発生時の付加的な性質とすべ
きか、はたまた、スケーリングの異なるアスペリティがあると考えるべきか、さらに研究を詰める必要があると思われる。
短周期強震動発生源の分布の観点で、「関口春子・吉見雅行，2006」の震源モデ
ルの広帯域化では、長波長震源モデルに無いより小さいスケールの不均質を、断層面上に均等にランダムに配しているが、この作業が海溝型地震の震源像に合っているのかどうか、高密度の観測網で観測された東北地方太平洋沖地震の震源像の解析により、詳細に確認してみる必要がある。
また、すべり時間関数について、「関口春子・吉見雅行，2006」では、「中村洋
光・宮武隆，2000」を用いているが、この関数は M7 クラスの震源像に基づいた検討であり、海溝型巨大地震にそのまま適用するには疑問が残る。難しい問題ではあるが、強震動生成に関与する部分のすべり関数は、予測地震動に大きな影響を与えるため、東北地方太平洋沖地震の観測に基づいて検討すべき課題の１つと考えられる。

・震源パラメータの変動による地震動のばらつきの見積りについて「関口春子・他，2006」では、「関口春子・吉見雅行，2006」の手法を用いた想定南海地震の地震動予測において、微視的震源パラメータに経験的に得られている程度の変動量を与え、地震動のばらつきの生じ方を調べた。変動を与えた微視的震源パラメータは、平均破壊伝播速度、アスペリティの応力降下量、破壊開始点、および、マルチスケールの不均質化で様々な大きさのパッチをランダムに配する際に用いられる乱数である。アスペリティの応力降下量の増加，破壊伝播速度の増加は，地震動のレベルをおよそ全体的に上げるので、変動による影響が予測しやすい。しかし、破壊伝播方向を変化させるようなパラメータの変動は、影響の出る範囲や影響のレベルが大きく変わり、複雑である。
「関口春子・他，2006」
では、破壊開始点を数十 km ずらしただけでも地震動が何倍も変わる地点も生じた。
そのため、破壊開始点位置など破壊伝播方向を左右するパラメータについては、極端に違いのある設定（例えば、破壊開始点を西端に置いた場合と東端に置いた場合など）だけでなく、たくさんのケースを実施する必要がある。また、「関口・他，2006」では、マルチスケールの不均質を与える際に用いる乱数が、地震動の367

見積もりに少なからぬ幅を与えることを示した。東北地方太平洋沖地震の短周期成分発生域と、強震動生成域や大すべり域の位置関係を再検討し、小スケール不均質の分布についてより狭い拘束を与えることが可能になれば、地震動の見積り幅を小さくすることができるかもしれない。

2) 2011 年東北地方太平洋沖地震・余震時の大阪堆積盆地の長周期地震動シミュレーション
2011 年３月 11 日の東北地方太平洋沖地震の本震（Mw9.0）とその最大余震である
2011 年茨城県沖を震源とする地震(Mw7.7)では、それぞれの震源域から 500km 以上も離れた大阪堆積盆地の湾岸部において大振幅の長周期地震動が観測され、その長周期地震動の卓越周期に近い固有周期をもつ超高層ビルにおいて非構造部材の破損やエレベータの不具合が生じた。南海トラフの巨大地震はその震源域が大阪に更に近く、より大振幅の長周期地震動が一般の震度と関係する周期帯域の強震動と相まって襲うことは確実であり、この長周期地震動特性を把握しておくことは、南海トラ
フの巨大地震時の大都市圏の持つ地震ハザードを正確に知り，減災対策に活かす手始めとなる。
東北地方太平洋沖地震本震と最大余震の大阪堆積盆地内外の地震動記録を比較す
ると、長周期地震動の特徴が類似していることがわかった。周期１～10 秒程度の周期帯域においての擬似速度応答スペクトルの空間分布は、本震時と最大余震時では卓越周期と最大応答値は似通っていた（周期 10 秒以上では異なる。これは本震が周期 10 秒以上の波動にもエネルギーを持つが最大余震はそれほどでもないから）。周期 10 秒以下の周期帯域の本震の震源過程は複雑であることから（例えば「K. Asanoand T. Iwata, 2012」）、震源過程が比較的簡単な最大余震である茨城県沖の地震を検討
対象として、波動伝播シミュレーションを行った。地震波動の計算手法には、不連続格子を用いた３次元差分法「S. Aoi and H. Fujiwara, 1999」を用いた。
与えた震源モデルと地下構造モデルは以下の通りである。最大余震について、計算したい周期帯域を再現できる適切な震源モデルがなかったため、ここでは震源に近い CHBH14（KiK-net 銚子)の周期５～10 秒の波形を再現する点震源モデルによる
震源時間関数を推定し、その震源時間関数と点震源を用いた波動シミュレーションを行った。
地下構造モデルは全国一次地下構造モデル 2009 年版「K. Koketsu et al., 2008」及び2012 年版「K. Koketsu et al., 2012」を用いた。本研究対象の大阪堆積盆地付近の地下
速度構造モデルは 2009 年版から 2012 年版では一部が変更されている。2012 年版作
成にあたり、大阪堆積盆地周辺の岩盤観測点においての地震動の再現性が上がるよ
うに地下構造モデルの修正が行われたためである「地震調査委員会, 2012」。図３−10−⑤−３には２つの地下構造モデルの、大阪湾岸の観測点 SKS(咲洲)を通る緯度での東西断面を示している。2009 年度版に比べ 2012 年版では大阪堆積盆地の東側の盆地外に Vs2.4 km/s およびそれより S 波速度が遅い層が厚く設定されていることがわかる。
盆地内外の速度コントラストは、2012 年版が 2009 年版に比べ小さくなっている。こ368

のため、2012 年版モデルによる計算結果は、盆地内の地震動がより小さくなっていることが確認されたため、ここでは 2009 年版モデルの結果を用いることとした。なお、最大余震の震源は、2009 年版モデルの範囲（東経 141.1 度以西）よりも東側にあるため、東経 141.1 度以東は 2012 年度版モデルを利用し、東経 141.1 以西は 2009年度版モデルを用い、両モデルの境界は線形補間でつないだ。海洋性地殻上面以深は 2012 年版のモデルを全域に用い、層毎の物性値は 2012 年版全国一次地下構造モ
デルのパラメータに統一した。速度構造モデルの最小 S 波速度は 350m/s で、数値的分散の影響を抑えるため、１波長あたり最小６グリッドを与え、計算有効周期３秒以上を確保する gridding を行った。図３−10−⑤−４に計算領域を、図３−10−⑤−５には震源と SKS を通る測線での地下速度構造モデルの断面を示す。
図 ３ −10− ⑤ − ６ に 震 源 か ら SKS に 向 か う 強 震 観 測 点 で の ト ラ ン ス バ ー ス 成 分(N344E 方向)の観測とシミュレーション結果を示す。波形例は、図３−10−⑤−４の地
図に赤で示した観測点の、周期４～10 秒の速度波形を比較している。関東平野内（震央距離約 200km 以内）、濃尾平野内（約 400～450km）、大阪堆積盆地内（約 550km）
において、震幅の増大と長い震動継続時間が再現されている。一番下の SKS（咲洲）
観測点の地震動には２つの波群が存在する点などは再現されているが、振幅自体は過小評価であることがわかる。水平動の擬似速度応答スペクトル値（7 秒、減衰 5％）の計算値／観測値を距離別に示したものが図３−10−⑤−７である。ほとんどの観測点において、この値が 1/2～２倍の範囲に収まっていることから、現在のモデルの妥当性はこの程度あることがわかる。ただし、大阪堆積盆地内、及び大阪堆積盆地よりやや東の地点では、それより過小評価の点が多くみられることから、この付近のモデルの検証と向上が今後の課題と考えられる。

図３−10−⑤−３

大阪堆積盆地の東西方向における全国一次地下構造モデル 2009 年版と2012 年版の違い。左が 2009 年版、右が 2012 年版。オレンジ色のところが SKS 観測点の経度に対応する。

369

CHBH14SKS

図３−10−⑤−４

最 大 余震から大 阪堆積盆 地 に至る観測 記録の観 測 点とモデル 計算領域

（赤波線の枠）。

図３−10−⑤−５（左）震源から SKS を通る測線断面のモデルの S 波速度構造。▲が SKS
で、★が震源位置を表す。（右）深さ５km までの S 波速度構造。

370

CHBH14SKS

図３−10−⑤−６

N344E 方向の震源から大阪咲洲（SKS）の観測点での観測波形（黒）とシミュレーション波形（赤）の比較。周期４～10 秒の速度波形。一番上が CHBH14 で一番下が SKS。

図３−10−⑤－７

震源から大阪堆積盆地に至る観測点における、擬似速度応答スペクト

ル（７秒、５％）の合成波形／観測波形を距離で並べた。青、黄土、緑はそれぞれ関東
平野、濃尾平野、大阪盆地内の観測点に対応。
371

(d) 結論ならびに今後の課題
震源モデルの高度化に関しては、2011 年東北地方太平洋沖地震の震源像について得られた知見に照らして既往の手法で採られている考え方や設定を見直し、今後の改良の方向について検討を行った。これまでの予測震源モデルで広く用いられてきたアスペリ
ティ-背景領域の考え方やアスペリティに設定していたスケーリング則について、見直す必要があることがわかったが、その際、この目的に注目して東北地方太平洋沖地震の震源を調べ直す必要がある。
地下速度構造モデルの高度化に関しては、2011 年東北地方太平洋沖地震余震の震源域から大阪堆積盆地に至る長周期地震動シミュレーションを行い、モデルの妥当性検証を行った。周期４～10 秒の長周期地震動が、関東平野、濃尾平野、大阪盆地で増幅され、震動継続時間が伸長する特徴や卓越周期といった特徴はある程度は再現された。擬似応答スペクトルで比較すると、ほとんどの観測点においては 1/2 から２倍の、いわゆ
る倍半分の中で評価できているが、一部に過小評価の点などがあることがわかった。今後、これらの定量的評価や、他の研究グループで提案されている大阪堆積盆地の地下構造モデルなどのモデル妥当性検証を継続し、速度構造モデルの高度化を進めていく。

(e) 引用文献
1) Aoi, S. and H. Fujiwara, 3-D finite difference method using discontinuous grids, Bull.
Seism. Soc. Am., 89, 918–930, 1999.
2) Asano, K. and T. Iwata, Source model for strong ground motion generation in thefrequency range 0.1-10 Hz during the 2011 Tohoku earthquake, Earth Planets Space, 64,
1111–1123, 2012.
3) Asano, K., H. Sekiguchi, T. Iwata, W. Suzuki, S. Aoi, and T. Kunugi, Source Process ofthe 2011 Off the Pacific Coast of Tohoku Earthquake, H. Kawase (ed.), in Studies on the2011 Off the Pacific Coast of Tohoku Earthquake, Natural Disaster Science and
Mitigation Engineering: DPRI reports, Springer, Japan, 17–36, 2014.
4) 中央防災会議，東南海，南海地震等に関する専門調査会 (第 16 回)資料，2003．
5) 中央防災会議，日本海溝・千島海溝周辺海溝型地震に関する専門調査会（第 10 回）
資料，2005.
6) Koketsu, K., H. Miyake, H. Fujiwara, and T. Hashimoto, Progress towards a Japan
Integrated Velocity Structure Model and Long-Period Ground Motion Hazard Map, Proc.
14th World Conf. Earthq. Eng., Paper No. S10-038, 2008.
7) Koketsu, K., H. Miyake, and H. Suzuki, Japan Integrated Velocity Structure Model
Version 1, Proc. 15th World Conf. Earthq. Eng., Paper No. 1773, 2012.
8) 地震調査委員会,南海地震（昭和型）の長周期地震動予測について（2012 年試作版）
http://www.jishin.go.jp/main/chousa/12_choshuki/index.htm, 2012.
9) Lay, T., H. Kanamori, C.J. Ammon, K. D. Koper, A. R. Hutko, L. Ye, H. Yue, T. M.
Rushing, Depth-varying rupture properties of subduction zone megathrust earthquakes, J.
Geophys., Res., 117(B4), B04311, doi:10.1029/2011JB009133, 2012.
37210) Mai, M. P. and G. C. Beroza, A spatial random field model to characterize complexity inearthquake slip，J. Geophys. Res., doi:10.1029/2001JB000588, 2002.
11) 宮腰研・A. Petukhin, 内陸地震の震源インバージョン結果に基づいた破壊伝播速度 の 不 均 質 性 に 関 す る 検 討 ， 地 球 惑 星 科 学 関 連 学 会 2005 年 合 同 大 会 予 稿 集 ，
S046-P002，2005．
12) 中村洋光・宮武隆，断層近傍強震動シミュレーションのための滑り速度時間関数の近似式，地震 2，53，1–9，2000.
13) Pulido, N., H. Tavera, Z. Aguilar, D. Calderon, M. Chlieh, T. Sekiguchi, S. Nakai and F.
Yamazaki, Mega-earthquakes rupture scenarios and strong motion simulations for Lima,
Peru, Proceedings of The International Symposium for CISMID 25th Anniversary
Technological Advances and Learned Lessons from Last Great Earthquakes and Tsunamis
in the World, Paper No. TS-6-2, 2012.
14) 関口春子・吉見雅行，広帯域地震動予測のための海溝型巨大地震のアスペリティ
モデルのマルチスケール不均質化，巨大地震による長周期地震動―予測と今後の対応策－，月刊地球，号外 No.55，103–109, 2006．
15) 関口春子・吉見雅行・吉田邦一・堀川晴央，海溝型巨大地震の微視的震源パラメータの変動がもたらす予測地震動のばらつき，第 34 回地盤震動シンポジウム，
83-88，2006.
16) Suzuki, W., S. Aoi, H. Sekiguchi, T. Kunugi, Rupture process of the 2011 Tohoku-Oki
megathrust earthquake (M9.0) inverted from strong-motion data, Geophys. Res. Lett., 38,
L00G16, doi:10.1029/2011GL049136, 2011.

⑥ シミュレーションに基づく南海トラフ地震津波の検証・評価(a) 業務の要約
1707 年宝永地震以前に発生した南海トラフ沿いの巨大地震については、史料等の残存状況の悪さから、破壊域の大きさなど、その震源象がはっきりしない地震が含まれる。
これらの地震の震源象を歴史資料の調査と津波・地震動のシミュレーションによって解明することが、近い将来発生する南海トラフ巨大地震の地震像を予測する上で重要となる。慶長九年十二月一六日（1605 年２月３日）の地震（以下、「慶長地震」）は、京都を含む西南日本では地震動の記録がないが、九州、四国、近畿、東海から外房、八丈島にかけての広範囲で大津波の記録が残っているので、南海トラフ沿いの津波地震と考えられてきた。しかし、同時史料における関東の大地震の記述と、外房における地震動と大津波の記録は、南海トラフ沿いの津波地震では説明できない。慶長地震の特徴は、2010
年 12 月 22 日に父島近海で発生したＭ7.4 のアウターライズ地震の特徴（西南日本は無感であったが関東地方以北に震度１～２の異常震域を生じ、父島、八丈島、南関東以西
で 50 cm 以下の津波を観測した「気象庁，2012」）と類似している（図３－10－⑥－１）。
したがって、慶長地震は南海トラフ沿いの津波地震ではなく、伊豆-小笠原海溝沿いの巨大地震であった可能性もあると考えられるので、歴史資料の再調査と津波の数値シミュレーションから検討を行った。
373

検討では、伊豆−小笠原海溝沿いのプレート間地震とアウターライズ地震の双方の断層モデルを、断層の大きさとすべり量を複数仮定して津波の数値シミュレーションを行い、慶長地震の津波高を最もよく説明できる断層モデルを探した。その結果、北緯 30.2°、
東経 141.8°付近のプレート間における 100×100 km の断層面に 10 m のすべり量を与えた断層モデルが慶長地震の津波高を最もよく説明することができた。また、アウターラ
イズ地震でも、北緯 30.2°、東経 142.5°付近の 100×50 km の高角な断層面に 10ｍのすべり量を与えた断層モデルが、慶長地震の津波高分布を良く説明した。しかしながら、プ
レート間地震の断層モデルより、計算された津波高と慶長地震の津波高との一致はよくなかった。
したがって、慶長地震は従来考えられてきた南海トラフ沿いの津波地震ではなく、
伊豆-小笠原海溝沿いの北緯 30.2°、東経 141.8°付近において、100×100 km 程度の領域（アスペリティ）が 10ｍ程度のすべり量で破壊するようなプレート間巨大地震であった可能性が非常に高いと結論できる。あるいは、アウターライズ地震であった場合でも、
北緯 30.2°、東経 142.5°付近で 10 m 程度の破壊を生じさせるような巨大地震であった可能性が高いことも結論された。しかしながら、それぞれの地震の海溝軸方向の広がりについては、少なくとも、10 m 程度のすべりの破壊が長さ 100 km 以上にわたっては生じていないことが津波シミュレーションから明らかになったが、すべり量が 5 m 程度かそれ以下のあまり大きくない破壊がどこまで広がっていたのかは、慶長地震の津波高の記録が少なくて精度もよくないので、明確な結論を出すには至っておらず、今後の課題である。
次に、津波・地震動シミュレーション結果と比較するためのデータである過去の南海トラフ巨大地震の津波高、震度分布、被害分布の再検討と再整理、未発掘の文献を探し出すために、徳島県、高知県、愛媛県、三重県、奈良県の公立図書館や文書館をまわり、歴史資料と地方史研究会や郷土史研究会などの論文集に発表された過去の南海巨大地震の津波・地震動や被害に関する文献の複写収集を行った。さらに、1854 年安政東海地震の断層モデルの改善に向けて、三重県沿岸での詳細な津波高分布図を作るため、三重県度会郡南伊勢町の五ヶ所湾における津波堆積物の野外調査を行った。また、南海
トラフとプレート間地震の繰り返しの時空間パターンが異なる日本海溝におけるプレート間巨大地震の連動性を明らかにし、南海トラフのプレート間巨大地震の連動性と比較して次の南海巨大地震の連動性を評価するために、日本海構における巨大地震の発生
履歴を長期にわたって解明する目的で、青森県沿岸（おいらせ町、階上町）、岩手県沿岸（洋野町、普代村、田野畑村、宮古市）の沿岸部において、津波堆積物の野外調査を
行った。

374

図３−10−⑥—１

2012 年 12 月 22 日の父島近海地震（Ｍ7.4）によって観測された津波高分布と震度分布「気象庁，2012」(b) 業務の実施方法
慶長地震の断層モデルの検討を行うために、まず、慶長地震で記録された津波高の再検討を行った「石橋克彦・原田智也，2013」。各地の津波高について、歴史史料の再検討と津波高の推定に関する既往研究（「羽鳥徳太郎，1975」、「 山本武夫･萩原尊礼，1995」、「村上仁士・他，1996」）を比較検討した。八丈島の津波高については「古村孝志・他，
2012」の調査結果を参照し、高知県土佐清水市の三崎における津波高は本研究で独自に推定した。次に、伊豆−小笠原海溝沿いのプレート間地震とアウターライズ地震の双方の断層モデルを仮定し、慶長地震の津波高と比較した。プレート間地震については、海溝軸に沿った 100×100 km の 6 枚の小断層を仮定し、複数の小断層の組み合わせによってプレート間巨大地震を表現した。各小断層の幅は、スラブ上面の 20 km 等深線と断層面の下端の深さが一致するようにした。断層面の傾斜角は 25 度である。すべり角は90°を仮定した。アウターライズ地震に関しては、アウターライズ領域に、海溝軸とほぼ平行になるように 100 km×50 km の小断層を 6 枚仮定し、複数の小断層の組み合わせ
でアウターライズ巨大地震を表現した。各小断層の傾斜角は 70°を仮定した。津波の伝播計算は、26〜36°N、130〜145°E の範囲において、コリオリ力を考慮した非線形長波式「K. Satake, 1995」を用いて行った。地震の発生後４時間の津波を計算し、沿岸における最大津波高を比較した。津波堆積物の野外調査は、簡易掘削装置を用いて１〜２ｍの深さの範囲で地中の試料の採取を行った。

(c) 業務の成果
1）史料の再調査と津波の数値シミュレーションによる慶長地震の断層モデルの再検討
a）慶長地震の津波高の再検討
慶長地震の津波高の再検討結果を表３−10−⑥−１と図３−10−⑥−２に示す。慶長
地震では四国の太平洋岸における津波高が非常に高いことが特徴である。特に、375

足摺岬付近、室戸付近や徳島南東岸では 10 m 程度の津波高になる。それに対して、九州の南東岸である鹿児島県の大隅半島では３〜４m 程度の津波高であり、和歌
山県有田市付近でも南の田辺市付近でも被害がない（有田市付近：被害記録は『和歌山県有田群地震津浪の記事』によるが、これは明治以降に書かれたものである
ようで信憑性が低い。そして、「天正年間地大に震ひ、紀州沿岸に津浪襲来し云々」という記事において、「地大に震ひ」という常套句も疑わしい「石橋克彦，2014」。
田辺市付近：「山本武夫・萩原尊禮，1995」は、宝永地震・安政南海地震の大被害を受けた田辺市の史料を検討し、慶長九年時は地震動・津波ともに記録されるほどでは無かったとしている。田辺市は有田市より約 40 km 南に位置するので、このことから有田市付近での被害津波も否定される）ことから大津波は襲わなかったと結論された。また、史料の記述の再検討から鹿児島湾奥の加治木における３m以上の津波高は否定された。さらに、八丈島、房総半島における津波高は、それぞれ、７〜８m、３〜５m となった。

表３−10−⑥−１

慶長地震の既往研究による津波高の比較と本研究で再検討された津波高「石橋克彦・原田智也，2013」（最右列黄色の部分）。
「？」は津波高が推定できなかったこと、「×」は津波が到達していないか、到達しても津波高ががほとんど高くないことを示す。
県地

名｢羽 鳥 ,1975｣｢山 本 ･荻原 ,1995｣｢村 上 ･他 ,1996｣｢石 橋 ･原田 ,2013｣

千葉県

外房

5〜 7

＜3

3〜 5

東京都

八丈島

10〜 20

＜ 10

7〜 8

静岡県

下 田 ・ 田牛2〜 3

西 伊 豆 ・仁 科

3〜 4

舞阪

4〜 5

3〜 4
1〜 2

1〜 2

新居

5〜 6

5〜 6

旧 白 須 賀宿

6〜 7

5〜 6

3〜 4

2〜 3

大 湊 ,被 害無？

愛知県

伊 良 湖 岬・堀切

三重県

伊勢

4〜 5

和歌山

辰 ヶ 浜 (有 田 )

4〜 5

広

4〜 5

徳島県

3〜 4

牟岐

×
被害無

高知県

5〜 6

6

5

9

9

4〜 5

4〜 5

5〜 6

5〜 6

8.6〜 13.0

10〜 13

8〜 10

8〜 10

宍喰

5〜 6

5〜 6

野根

4〜 5

被害なし

佐喜浜

8〜 10

10〜 13

室津
奈半利

3〜 4

6

6

安芸

3〜 4

3〜 4

3〜 4

佐賀

4〜 5

4〜 5

4〜 5
6〜 8

三崎浦
鹿児島

＞6

4〜 5

浅川
鞆浦

×

南 の 田 辺 ,被 害 無

薩 摩 ・ 大隅

1〜 2

3〜 4
＞3

鹿 児 島 〜加 治

376

×

図３−10－⑥－２

再検討された慶長地震の津波高分布b）津波シミュレーションによる慶長地震の断層面の推定
伊豆—小笠原海溝のプレート間とアウターライズに仮定した断層の大きさとすべり量の組合せを色々変えて津波シミュレーションを行い、慶長地震の津波高と比較した。その結果、慶長地震の津波高分布を最もよく説明できる断層モデルは、
北緯 30.2°、東経 141.8°付近のプレート間の 100×100 km の断層面に 10ｍのすべり量を与えたモデル（以下、プレート間地震モデル）であった（図３−10−⑥－３、
図３−10−⑥－４）。また、アウターライズ地震では、北緯 30.2°、東経 142.5°付近の100×50 km の断層面に 10 m のすべり量を与えたモデル（以下、アウターライズ地震モデル）が、慶長地震の津波高分布を最もよく説明した（図３−10−⑥－３、図３−10−⑥－４）。

図３−10−⑥−３

慶長 津波の 津波 高分 布を 最 も良く 説明 する 伊豆−小 笠原海 溝沿 いの プレート間地震モデル、アウターライズ地震モデルと「安中正・他，2003」による南海巨大地震モデルの断層面の地表投影。断層面の付近の数値はすべり量を示す。
377

図３−10−⑥−４

慶長地震の津波を最も良く説明する伊豆−小笠原海溝沿いのプレート間地震モデルとアウターライズ地震モデル、「安中正・他，2003」による南海巨大地震モデル
により計算された津波高（青）と慶長地震の津波高（赤）との比較。

次に、八丈島−房総半島、大阪湾−和歌山において、伊豆−小笠原海溝のプレート間地震モデルによる津波高と「安中正・他，2003」による南海巨大地震の断層モ
デル（以下、「南海巨大地震モデル」）による津波高との比較を行った（図３−10−⑥−５）。八丈島では、南海巨大地震モデルによって計算された津波高は慶長地震の津波高の半分程度にしかならないが、プレート間地震モデルによって計算された津波高は慶長地震の津波高を良く説明する。また、外房における慶長地震の津波高もプレート間地震モデルの方が良く説明する。大阪湾周辺、和歌山県和歌山市付近においては、南海地震モデルによって計算された津波高は３〜４ m 程度になり、和歌山市付近の津波はほとんどなかったとする本研究で再検討した結果（表３−10−⑥−１、図３−10−⑥−２）と矛盾するが、プレート間地震モデルによって計算された津波高は１ｍ程度以下になり、再検討結果をよく説明することが分かった。
したがって、プレート間地震モデルの方が、南海巨大地震モデルよりも慶長地震による津波をより良く説明する事が分かった。

378

図３−10−⑥−５

八丈島−房総、大阪−和歌山におけるプレート間地震モデルと南海巨大地震モデルによって計算された津波高の比較

2）徳島県、高知県、愛媛県、奈良県における南海巨大地震の津波に関する文献・史
料調査
津波・地震動シミュレーション結果と比較するためのデータである過去の南海ト
ラフ巨大地震の津波高、震度分布、被害分布の再検討と再整理、未発掘の文献を探し出すために、徳島県、高知県、愛媛県、奈良県の公立図書館や文書館をまわり、歴史資料と地方史研究会や郷土史研究会などの論文集に発表された過去の南海巨大地震の津波・地震動や被害に関する文献の複写収集を行った。例えば、高知県では『土佐史談』や『須崎史談』などの地方史研究論文集に収録された宝永地震、安政南海地震に関する調査・研究部分を、奈良県大和郡山市の柳沢文庫では史料『楽只堂年録』の宝永地震発生前後の部分を複写収集した。

3）1854 年安政東海地震の断層モデルの改善に向けた、三重県度会郡南伊勢町の五ヶ所湾周辺における津波高の詳細分布作成のための文献・史料調査と津波堆積物調査
a）文献史料調査
三重県庁県史編さん班において古地図の調査を行い、明治初期の地籍地図を撮影記録した。同時代に関東地方について作成された迅速測図ほどの精度は無いとはいえ、近代的な測量のもとに作成された三重県下の地図としては最古級のもの
379

である。また、三重県玉城町教育委員会において 1854 年安政東海地震津波に関する史料調査を行い、田丸藩の公文書に該当記録を発見し、複写収集した。当史料は未刊行の新史料である。同時代史料であり、度会郡慥柄組（現、南伊勢町とほぼ同じ領域）に所属する各村落について被害規模がまとめられており、その精度も含め非常に貴重な史料である。

b）津波堆積物調査
三重県度会郡南伊勢町神津佐区、泉区、宿区において、簡易掘削機を用いて７ヶ所、最大深度 6.0ｍまでの掘削を行った。それぞれの内訳は神津佐区後河内地区（U-001: コア長 0.8ｍ、U-002: コア長 1.3ｍ）、神津佐区藤原地区（U-003: コア長1.4ｍ）、泉区（U-004: コア長 2.3ｍ、U-005: コア長 1.7ｍ、U-006: コア長 1.5m）である。コアが得られた地点を図３−10−⑥−６に示す。これら全ての地点において、
表層 0.4～0.8ｍは人為的に造成された砂礫層からなる（盛土）。その下位は比較均質で有機質な泥層、もしくは礫層（角礫）、礫混じり砂層からなり、砂層および礫層は津波による堆積物の可能性があり、今後の調査・分析によって検討を進める必要がある。

図３−10−⑥−６ 南伊勢町津佐区、泉区で行われた津波堆積物調査地点

4）日本海溝沿いのプレート間巨大地震の発生履歴の解明に向けた岩手県−青森県の太平洋岸における津波堆積物調査
南海トラフとプレート間地震の繰り返しの時空間パターンが異なる日本海溝におけるプレート間巨大地震の連動性を明らかにし、南海トラフのプレート間巨大地震の連動性と比較して次の南海巨大地震の連動性を評価するために、日本海溝沿いで発生する巨大地震の発生履歴を長期にわたって解明する目的で、2013 年 11 月に岩手
県宮古市田老地区（図３−10−⑥−７）において東京大学地震研究所によって行われた
380

津波堆積物調査に参加し、地質試料の採取を行った。また、青森県沿岸と岩手県北部の沿岸地域において行われた、次年度の掘削予定地選定のための東京大学地震研究所による予備調査にも参加した。予備調査では、青森県沿岸（おいらせ町、階上町）、岩手県沿岸（洋野町、普代村、田野畑村、宮古市）の沿岸部において、簡易掘削機を用いた地質試料採取を行った。

図３−10−⑥−７

岩手県宮古市における津波堆積物調査地点(d) 結論ならびに今後の課題
慶長九年十二月一六日（1605 年２月３日）の地震（「慶長地震」）は南海トラフ沿いの津波地震ではなく、伊豆-小笠原海溝沿いの巨大地震であった可能性もあると考えられるので、歴史資料の再調査と津波の数値シミュレーションから検討を行った。検討の結果、慶長地震は従来考えられてきた南海トラフ沿いの津波地震ではなく、伊豆-小笠
原海溝沿いの北緯 30.2°、東経 141.8°付近において、100×100 km 程度の領域（アスペリティ）が 10 m 程度のすべり量で破壊するようなプレート間巨大地震であった可能性が非常に高いと結論された。あるいは、アウターライズ地震であった場合でも、北緯 30.2°、
東経 142.5°付近で 10 m 程度の破壊を生じさせるような巨大地震であった可能性が高い
ことも結論された。しかしながら、それぞれの地震の海溝軸方向の広がりについては、
少なくとも、10 m 程度のすべりの破壊が長さ 100 km 以上にわたっては生じていないことが津波シミュレーションから明らかになったが、あまり大きくない破壊がどこまで広がっていたのかは、慶長地震の津波高の記録が少なくて精度もよくないので、明確な結論を出すには至っておらず、今後の課題である。文献・史料調査は、本年度では徳島県、高知県、愛媛県、奈良県で行ったが、次年度以降は九州の大分県、宮崎県、東海地方の三重県、愛知県、静岡県でも行わなければならない。津波堆積物の野外調査に関しては、
本年度は予備調査的な試掘削であったので、次年度以降は本格的な調査に移行していかなければならない。

(e) 引用文献
1） 安中正・稲垣和男・田中寛好・柳沢賢，2003，津波シミュレーションに基づく南381

海トラフ沿いの大地震の特徴，土木学会地震工学論文集 (CD-ROM)．
2） 古村孝志・今井健太郎・前田拓人・原田智也，2012，1605 年慶長地震における八丈島の津波痕跡高の再検討，日本地球惑星科学連合 2012 年大会，SSS38-P1
3） 羽鳥徳太郎，1975，明応 7 年・慶長 9 年の房総および東海南海道大津波の波源，
地震研究所彙報，50，171-185.
4） 石橋克彦・原田智也，2013，1605（慶長九）年伊豆-小笠原海溝巨大地震と 1614
（慶長十九）年南海トラフ地震という作業仮説，日本地震学会 2013 年度秋季大会，D21-03．
5） 石橋克彦，2014，『南海トラフ巨大地震—歴史・科学・社会』，叢書-震災と社会，
岩波書店．
6） 気象庁，2012，地震年報（CD-ROM）
7） 村上仁士・島田富美男・伊藤禎彦・山本尚明・石塚淳一，1996，四国における歴史津波（1605 慶長・1707 宝永・1854 安政）の津波高の再検討，自然災害科学，
15，39-52.
8） Satake K., 1995, Linear and nonlinear computations of the 1992 Nicaragua earthquake
tsunami, Pure Appl. Geophys. 144, 455–470.
9） 山本武夫・萩原尊禮，1995，慶長九年（一六〇五）十二月十六日地震について―
東海・南海沖の津波地震か，萩原尊禮（編著）『古地震探求―海洋地震へのアプ
ローチ』，東京大学出版会，160-251．

⑦ 実用的な津波土砂移動シミュレーション手法の開発(a) 業務の要約
本業務では、これまでに開発されてきた津波土砂移動シミュレーションモデルの網羅的なレビューを行い、それらモデルのモデリング手法や物理過程について検討を行い、実用的な津波土砂移動シミュレーション手法の開発に向けた問題点を抽出した。さらに、
津波土砂移動モデルの検証を行うために、2011 年東北地方太平洋沖地震による津波で大規模な地形変化が生じた宮城県石巻市新北上川河口部を対象とした今次津波前後の地形データの作成を行った。

(b) 業務の成果
1) 津波による土砂移動モデルに関するレビュー
津波土砂移動モデルに関する研究の目的は、津波による地形変化・土砂災害の予測と、津波堆積物を用いた古い津波の規模や波源の推定、の２つの面がある。防災上の理由から、近年益々注目され、重要度の高い研究分野である。特に古津波に関しては、湿地や湖沼に津波で流入した砂は識別と認定が比較的容易であるので、これを使った津波規模の研究は有望であると考えられている（例えば，「行谷佑一・他，
2010」）。2011 年の東北地方太平洋沖地震津波の浸水域は、堆積物に基づいて推定されていた貞観津波の浸水域と類似しており「菅原大助・他，2011」、津波堆積物の有用性が注目されるようになった。日本では、地震津波の履歴と規模、リスクを評価
382

するために、津波堆積物の積極的活用が進められようとしている。
津波堆積物の堆積形態は非常に多様であり、堆積物の分布や特徴は津波の流況や
地形、植生等の影響を受けるほか、土砂供給源の分布と土砂粒子構成によっても大きく変わる。たとえば、津波の砂は地形の凹凸を埋めるように堆積することが知られている。供給源から離れるに従って細粒化を見せるが、遡上途中に別の供給源があれば、再び粗粒になることもある。また、波の繰り返しや押し引き等、流れ状況の変動を反映して複数の層が積み重なることも知られている「J. Goff et al., 2012」。
しかし、このような堆積の特徴と津波流況がどのように関連しているかはそれほど明らかでない。津波後の現地調査で得られる情報は、通常は浸水深（あるいは痕跡高）のみである。最近では津波来襲時の模様がビデオ映像に記録され、流速を推定する貴重な資料となっているが、数はまだまだ少ない「林里美・越村俊一，2012」。
しかも、これらのデータが津波堆積物を調査した地点と同じ場所で得られるとは限らない。
津波堆積物に内在する情報に、津波の流況がどのように寄与しているかも明確とはいえない。土砂移動に関する数理モデル構築に当たっては、複雑な物理過程をモ
デリングする必要があり，堆積物がどのようなプロセスを経て形成されたのか、また形成に至る時間スケールなども含めて仮定する。津波による土砂移動については，
水理実験による物理的な考察とともに，数理モデルの検証が可能であるが，設備の空間スケールや波・流れの強さ・時間的スケールに制約があり、実際の津波による土砂移動の状況を断片的にしかに再現できない。
このような背景から、数値モデルにより津波の流れ場を再現し、その状況下で起こる土砂移動の解析を行うことで、津波堆積物の分布や性質が流れの状況とどのように関連するかが研究されている。ここでは、津波土砂移動の既往研究を概観し、これまでに指摘されてきた、あるいは新たに見えてきた問題点を検討するとともに、今後の研究の展望について述べる。

a) 土砂移動モデルの基本的な構成
津波による土砂移動モデルは氾濫解析モデルと土砂移動に関する数理モデルを
カップリングさせて解く方法が一般的であり，土砂移動に関しては掃流砂と浮遊砂を区別して取り扱うことが多い．その構成は図３－10－⑦－１のようになっている「T. Takahashi et al., 2000」。まず、津波流況により算出される掃流力を用いて、掃流砂量、浮遊砂の巻き上げ量を求める（図３－10－⑦－１中の Bed load layer）。
次に、水深と流速を用いて掃流砂フラックス、浮遊砂の濃度分布 C とその移流拡散を計算する（図３－10－⑦－１中の Suspended load layer）。最後に、質量保存則に基づいて浸食・堆積による地形変化が計算される。掃流砂量 q B=f q (τ * ，τ *c ，…)
は無次元掃流力（シールズ数）の関数であると考えられ、既往の提案式やそれぞれの論文著者らが独自に実験を行って決定した式を用いる事もある。浮遊砂巻き上げ（浸食）量 E S =f E (τ * ，τ *c ，…)も掃流砂量と同様にシールズ数の関数として表す
方法、基準面浮遊砂濃度を求め、その拡散を計算して求める方法があり，掃流砂383

量式と同じく様々な提案式が存在する。

図３－10－⑦－１

土砂移動形態の構成図「T. Takahashi et al., 2000」。土砂移動形

態の構成図（Takahashi et al., 2000）．図中、h B は砂層厚，h s は水深，z B は砂層を含む地盤高，ρ s は砂の比重を示す。
浮遊砂の移動は移流拡散方程式が用いられる。次式は、「T. Yoshii et al., 2010」
による、水深平均された平面 2 次元の浮遊砂連続式である。

上式は第 1 項が浮遊砂濃度の時間変化を表し、第 2・3 項は浮遊砂濃度の水平移流、第 4・5 項は水平拡散項，第 6・7 項はそれぞれ浮遊砂巻き上げと沈降を表す項である。ここで、C は水深平均された浮遊砂濃度、D は全水深、M（=UD）および N（=VD）は流量フラックス、U およびは断面水深平均流速、K h は水平方向の渦動粘性係数である。浮遊砂沈降（堆積）量 D S =w S C は、上式で算定された平均または基準面の浮遊砂濃度と、粒子の沈降速度 w S を用いて計算する。沈降速度の算定
式は様々あるが、例えば Rubey の式「W. W. Rubey, 1933」で計算する。
3 次元計算の場合、静水圧近似で 3 次元化する。z 方向の離散化には、海底から水表面までをある割合で分割する σ 座標系が使われる。次の式は、「 A. Apotsos et al.,2011a」、「A. Apotsos et al., 2011b」で用いられている浮遊砂の連続式である。

ここで、c は浮遊砂の質量濃度、u 、v 、w はそれぞれ x、y、z 方向の流速、ε s,x 、ε s,y 、384

ε s,z はそれぞれ x, y, z 方向の渦拡散係数である。上式において、第 4 項の w S C= D Sが堆積量に、第 5 項から第 7 項が浸食量に相当する。
最終的には、質量保存則に基づいて浸食・堆積による地盤高の変化が計算される。次式は「T. Takahashi et al., 2000」による地形変化算定式の例である。

ここで、 z b は地盤高、 λ は空隙率である。地形変化は流れ場の解析に使う海底・陸上地形に反映される。
掃流力、掃流砂フラックス、浮遊砂の巻き上げ量、浮遊砂濃度とその分布が土砂移動モデルの核となる要素である。それぞれの算定式は、物理過程の考慮の度
合い、経験パラメータの決定方法が異なる様々なものが存在しており、どの式・パラメータを採用するかが、モデルの振る舞いに大きく影響する。

b) 代表的な津波による土砂移動モデル
以下では、各モデルの特徴と適用例を簡単に解説する。

i) Delft3D モデル
Delft3D は、洪水流れ、津波、跳水、段波の解析に適用することを目的に WL/Delft
Hydraulics（現 Deltares）がデルフト工科大と共同で開発したモデルで、非線形長波理論式と土砂移動・地形変化解析モデルを組み合わせた１～３次元の解析ができる。3 次元解析では、静水圧近似を仮定し、σ 座標系を用いている。波・流れ場の解析の検証は、ダムブレイク等の水理実験のほか、インド洋大津波を例に行われ、再現性が確かめられているほか、津波氾濫に対する海岸林や防波堤、サンゴ礁の影響も調べられている。土砂移動については混合粒径の掃流・浮遊による移動、土砂混入による流体の密度成層、干渉沈降など多様な物理プロセスを取り扱うことが可能で、津波に対してはこれまでに平面 2 次元解析（2DH）、断面2 次元解析（2DV）の解析が行われている「A. Apotsos et al., 2011a」、「A. Apotsos
et al., 2011b」、「A. Apotsos et al., 2011c」。Delft3D では、掃流砂フラックス、巻き
上 げ 量 、 浮 遊 砂 濃 度 分 布 は 既 往 の 様 々 な 提 案 式 か ら 選 択 で き る 。 例 え ば 、「 A.
Apotsos et al., 2011a」では、掃流砂は Van Rijn の式「L. C. Van Rijn et al., 2004」
で、浮遊砂巻き上げ量は k-ε 乱流クロージャモデルで決定される渦拡散係数をパラメータとして求めている。実例に対する適用としては、インドネシア国 Kuala
Meurisi 沿岸を例に 2DV 解析が行われ、倍半分のオーダーの概ね良好な結果を得ている「A. Apotsos et al., 2011a」。浸水域端の崖付近で細粒物が厚く堆積した現象は再現できていない。これは、単一粒径の計算のためであると説明されている。
「A. Apotsos et al., 2011c」では詳細な物理過程を考慮した計算結果が示されており、モデルのパフォーマンスを確認できる。しかし、計算条件に関して不確定な点があるため、現地砂分布データによる検証は定性的な比較にとどまっている。
385

シミュレートされた津波堆積物の粒度の内部構造とその形成プロセスは興味深
い結果であるが、津波堆積物形成過程の実スケールの観測データは存在せず、検証方法が課題である。2DH 解析は、アメリカ領サモア Fagafue Bay 沿岸を例に粗度の空間変化や潮位などを考慮して行われている「A. Apotsos et al., 2011c」。空間解像度は２～５m、0.2mm の単一粒径を仮定している。現地で見られた堆積物のパッチ状の分布の再現など、定性的には良好な結果を得ている。

ii) XBeach
XBeach は、The consortium of UNESCO-IHE, Deltares (Delft Hydraulics), Delft
University of Technology and the University of Miami が共同して開発したモデルで、非線形長波理論式と土砂移動モデルの 2DH 解析を同時に行うことができる。このモデルは、様々な解析解、大型の水槽実験、幾つかの現地データを対象に検証されている。混合粒径の砂の移動を、掃流砂・浮遊砂の両方を集約した全流砂量を用いて解析できる。これは、平衡土砂濃度と実際の土砂濃度の差から巻き上
げ・沈降量を求める方法である。平衡土砂濃度は Van Rijn の式「L. C. Van Rijn etal., 1993」で求める。密度成層や干渉沈降は土砂移動の緩和時間を用いて陰的に評価されると説明されている。砂の移動は、「G. Galappatti and C. B. Vreugdenhil,
1985」による移流拡散式を用いるが、水平拡散項は省略されている。混合粒径の砂の移動を解析でき、津波堆積物の水平・鉛直の粒度構造のシミュレートが可能である。
「L. Li et al., 2012」では、次に述べる津波の伝播・遡上モデル COMCOTと連成され、インドネシア国 Lhok Nga 沿岸を対象に解析を行っている。衛星画
像をもとに解析領域をカバーする土地条件図を作成し、土地利用に応じて粗度係数を与えていることが特色である。堆積物は、厚さ３cm の層が最大 20 層重なったものとして計算される。「R. Paris et al., 2010」による現地砂層厚データとの比較が行われ、砂の分布パターンを再現するなど定性的には良い結果を得ているが、場所によってはかなり過小評価になっている。粒度の鉛直構造の解析も行っているが、現地データに見られるような複数の上方細粒化等の構造は再現できていない。これは、用いた砂粒径が２階級（0.2 mm、0.8 mm）、計算で再現できる砂層の鉛直分解能が３m と低かったためと説明されている。

iii) Ontowirjo et al.のモデル
コーネル大学による津波伝播・遡上モデルである COMCOT（Cornell Multi-gridCoupled Tsunami Model）と、Ontowirjo らによる土砂移動モデルを連成したモデルである。COMCOT の信頼性と精度は 1992 フローレス、2003 アルジェリア、
2004 インド洋大津波を例に検証されている。
「B. Ontowirjo et al., 2012」においては、ブジネスク方程式の分散関係を考慮するため、方程式中に数値分散項が導入されている。また、浮遊砂巻き上げと基準面浮遊砂濃度に対する砕波の影響を考慮するため、底面剪断応力項と乱流粘性項を導入したと説明されている。掃流砂
量式は Ribberink の式「J. S. Ribberink, 1998」で求めている。浮遊砂の移流拡散
386

は Delft3D で用いられているものとほぼ同じ構成で、水平方向の移流拡散、鉛直方向の拡散を考慮している。浮遊砂巻き上げは渦粘性係数と乱流シュミット数をパラメータとして求める。基準面浮遊砂濃度は Van Rijn の式「L. C. Van Rijn et al.,1993」、「L. C. Van Rijn et al., 2004」で計算する。実験や実例による土砂移動モデ
ルの検証については触れられていないが、一般的な流砂量式、移流拡散式を用いており、構成は妥当であると考えられる。Lhok Nga を対象に、堆積物を単一粒径と仮定して解析が行われ、「R. Paris et al., 2009」のデータと比較されている。
定性的には良く一致する結果を得たと述べられ、土砂厚さのデータに対してかな
り良い整合性が示されている。しかし、浸水深の現地データに不確実性が見られる上、そのデータに対する計算値の偏差が大きく、土砂移動計算の精度に影響する流れ場の再現に疑念が残る。「R. Paris et al., 2009」では、現地調査データに基づく砂層厚の平面分布推定結果が示されているが、これを用いた計算結果の検証は行われていない。

iv) C-HYDRO3D「N. Kihara and M. Matsuyama, 2010」、「N. Kihara et al., 2012」は、非線形長波モデルと独自の土砂移動モデルを連成した数値モデル C-HYDRO3D を提案している。Delft3D と同じく、σ 座標を導入した 3 次元解析を行うことができる。このモデルでは土砂移動は掃流と浮遊の両方を考慮しており、掃流砂量式と浮遊砂
量式は Van Rijn の式「L. C. Van Rijn et al., 1993」、「L. C. Van Rijn et al., 2004」を用いる。浮遊砂量式は「B. Ontowirjo et al., 2012」のモデルと同じであるが、渦拡散係数の算定では「G. L. Mellor and T. Yamada, 1982」のクロージャモデルが使用されている。浮遊砂連続式は Delft3D や「B. Ontowirjo et al., 2012」のモデルと同じく、水平方向の移流拡散、鉛直方向の拡散を考慮している。C-HYDRO3D
は、インド洋大津波におけるスリランカ国 Kirinda 漁港と港湾を想定した水理実
験で検証が行われている「N. Kihara and M. Matsuyama, 2010」、「N. Kihara et al.,2012」。水理実験では、後述の STM「T. Takahashi et al., 2000」では再現できなかった、港内中央における局所的な堆積の再現に成功している。著者らは、このことは３次元解析で渦の形成にともなう２次流を解像することが出来たためと分析し、渦が生じる地形状況における３次元計算の必要性を指摘している。一方、同じく港湾地形の Kirinda の解析では、港湾周辺の津波の挙動の再現するため、インド洋地震津波に関する２つの波源モデルを比較して、津波到達時刻や高さが良く合うモデルを使って解析している。防波堤が設置された位置よりも陸上で堆積が起こった状況が再現されたが、海底の浸食・堆積の状況は再現できていない。

v) 「T. Takahashi et al., 2000」のモデルとその発展モデル「T. Takahashi et al., 2000」は、従来の平面２次元津波伝播遡上モデルと連成
して地形変化を解析する土砂移動モデルを提案した。「T. Nakamura et al., 2009」はこのモデルをシンプルに STM と呼んでいる。土砂移動は掃流砂・浮遊砂を分387

離して考慮されるが、堆積物は単一粒径で取り扱われる（図３－10－⑦－１）。
掃流砂量・浮遊砂巻き上げ量は著者らが独自に行った津波を模擬した高速な流れを再現できる閉管路実験で定めた式を用いている。Delft3D や Ontowirjo らのモデルとは異なり、浮遊砂の連続式で水平拡散項が省略されている。
「T. Nakamura
et al., 2009」は、津波を想定した構造物周辺の洗掘実験を対象に STM に水―砂
界面での水の浸潤と滲出速度の影響を考慮した修正シールズ数導入モデル MSM、STM に砂層表面の応力変動を考慮した有効応力モデル ESM の検証を行っている。
なお、彼らは、非線形長波方程式の代わりに、NWT-SWM と呼ばれる３次元数
値波動水路と砂骨格間隙水解析モデルを使用している。それによれば、STM による浸食量・堆積量は、水理実験データや MSM・ESM の予測に対して基本的に過小評価になることが示されている。
「T. Yoshii et al., 2010」は、STM を参考に、
同じく非線形長波理論と掃流・浮遊による土砂移動を考慮したモデルを提案し、STM との比較を行っている。STM との違いは、掃流砂量算定に芦田・道上の式「芦田和男・道上正規，1972」を用いている事、浮遊砂巻き上げの算定に独自の実験で定めた式を用いている事、浮遊砂の鉛直濃度分布と水平拡散を考慮していることである。「T. Yoshii et al., 2010」では、港湾を想定した大型の水理実験結果で検証を行っている。それによれば、STM の浸食量予測が過大であったのに対し、「T. Yoshii et al., 2010」のモデルではより実際に近くなったことが示されている。実例への適用は、気仙沼「T. Takahashi et al., 2000」、Kirinda「高橋潤・他，2008」・Lhok Nga「A. Gusman et al., 2012」を対象に行われている。気仙沼の解析では、湾内の浸食卓越傾向は再現されたが、浸食量が実測で 2,618,000 m3であったのに対し、計算では 458,000 m3 の過小評価になっている。Kirinda の解析では堆積量・浸食量のオーダーでは妥当な結果が得られているようであるが、実際の浸食・堆積分布とは大きな差が見られる。Lhok Nga への適用事例も、堆積物厚さのオーダーは妥当であるものの、計算による堆積物分布パターンは現地データとかなり異なるようである。

c) 津波による土砂移動モデルの問題点
上記の土砂移動モデルの研究を通じて、様々な問題が提起されてきた。大きくまとめると、土砂移動に関する物理過程、すなわち土砂移動モデルの構成、計算
条件、モデルの検証の３つに分類できる。それぞれの問題の解決には、今後の更なる研究の進展が必要である。ここでは、物理過程に関する検討を行う。
津波の挙動・土砂移動のプロセスは複雑であり、様々な計算条件・物理過程をどの程度詳細に考慮するかによって解析の精度が変わる。堆積物の浸食にとって最も重要なパラメータは、土砂の移動量を決める掃流力である。掃流力の算定では、しばしば Manning 則を仮定した粗度係数が用いられるが、これは一様定常流の仮定に基づいて決定されている。津波の流れ場は時間的に変動し、特に遡上先端や波の押し引きが顕著な場合、非一様性・非定常性は大きいと考えられる。このとき、掃流力と流速の時間的変動には位相差があることが知られており、流れ
388

場の流速と Manning 則の粗度係数を用いて掃流力を計算した場合、土砂の移動方
向や量を正しく求める事ができないことが知られている「M. B. Adityawan and H.
Tanaka, 2012」。さらに、Manning 則による掃流力は、水深が浅い時には過大に、深い時には過小に評価されることも指摘されている。Manning の粗度係数が表す植生など地表の状態は、その空間的分布を考慮すべきであることや時間発展的に破壊
さ れ る と し て 時 間 的 に も 変 化 す る と 考 え る の が 妥 当 で あ る と い わ れ て い る 「 A.
Apotsos et al., 2011a」、「A. Apotsos et al., 2011c」。
掃流力を精度良く求めたとしても、それを用いて算定される流砂量が正しいことを確認することは容易ではない。これまでに、様々な流砂量式が提案されてきているが、その中からどの式を使うのが良いか、選択基準は自明でない。しかも、河川の洪水や風波を想定したそれらの式が、流速１～20 m/s、周期 10 min 以上の
津 波 に 適 用 可 能 か ど う か に つ い て も 疑 念 が 指 摘 さ れ て い る （「 A. Apotsos et al.,2011b」、「N. Kihara and M. Matsuyama, 2010」）。混合粒径での解析も重要である事が指摘されている（「A. Apotsos et al., 2011c」、「L. Li et al., 2012」）。いくつかの解析例では土砂を単一粒径で取り扱っているが、この場合、津波堆積物に見られる重要な特徴である水平・鉛直方向の粒度構造を再現できない上、土砂の最終的な分布範囲に影響してしまうことから、モデルの検証を行う上でも不都合である。
海岸地形が複雑、あるいは津波の浸入方向が海岸線に垂直でない場合、断面計
算で津波の挙動を解析することの妥当性には疑問が指摘されている「A. Apotsos etal., 2011b」。他にも、遡上に伴うがれきの増加が流れの力学を変えてしまうこと「B.
Ontowirjo et al., 2012」、砕波による乱流の発達を考慮すること「A. Apotsos et al.,2011b」、間隙水圧勾配・堆積物への海水浸透の影響を考慮すること（「A. Apotsos etal., 2011b」、「T. Nakamura et al., 2009」）、土砂移動の３次元性を考慮すること等が指摘されている「N. Kihara et al., 2012」。しかし、これらのプロセスをモデルに取り込む こと は、 解析 の 複雑さ ・不 確実 さを 増 加させ てし まう 「A. Apotsos et al.,2011b」。考慮するプロセスが多いほど、それを表す物理過程モデルの選択基準、任意性のあるパラメータの選択基準に苦慮することになる。詳細なプロセスを導入する場合、パラメータに対するモデルの応答の確認が必要である。古津波堆積物への適用を念頭に置いた場合、後述の通りデータの不確実性の影響が大きいので、数々の複雑なプロセスをモデルに取り入れることの意義は良く検討する必要があるだろう。現地調査や実験を通じて解析すべき対象が何であるかを明確にし、それを解像できるモデルを構成することが重要である。

2) 新北上川河口部における津波前後の地形データおよび土地利用図の作成
上述のとおり、津波による土砂移動モデルの検証には、津波前と後の地形変化を詳細に知る必要がある。本業務では、宮城県石巻市新北上川河口部における地形変化についてのデータを作成した。
表３－10－⑦－１に用いたデータを示す。これらのデータを用いて、被災前後の地形データ（図３－10－⑦－２、図３－10－⑦－３）を作成した。被災後の地形デ389

ータについては、大規模な洗掘を受けた長面地区の多くのエリアでは、測量データを入手することができなった。そのため、国土地理院による地殻変動データや、宮城県による長面地区での現況地形図などの情報を用いて空間的に補間を行った。
土砂移動解析を行うためには、土砂供給が可能なエリアなどの選別が必要になる。
そのため、空中写真などを用いて、当該領域における土地利用図についても作成を行った。図３－10－⑦－４に新北上川周辺における土地利用図を示す。この図から、おおよその土砂供給源を推定することができる。

表３－10－⑦－１

地形作成に用いたデータ

データ

作成機関

発行年

備考

基 盤 地 図情 報 （公 共 測量 ）国 土 地 理院

2009

地 形 デ ータ 作 成に 使 用(陸 域 )

基 盤 地 図情 報 （基 本 測量 ）国 土 地 理院

2012

地 形 デ ータ 作 成に 使 用(陸 域 )

北 上 川 低水 路 平均 河 川縦 断図

国 土 交 通省

2010

地 形 デ ータ 作 成に 使 用(北 上 川 )中 央 防 災会 議 50m メ ッ シュ デー タ

内 閣 府「 日本 海
溝・千 島 海溝 周
辺 海 溝 型地 震に 関 す る専 門
調査会」

2006

地 形 デ ータ 作 成に 使 用(海 域 ほ か )55m DEM

国 土 地 理院

2011（震 災 後）

地 形 デ ータ 作 成に 使 用(陸 域 )

6

現 況 地 形図 （ 宮城 県 によ る長面 周 辺の 農 地 測量 デ ータ ）宮城県

2011（震 災 後）

地 形 デ ータ 作 成に 使 用(長 面 周 辺 圃 場 )7

基 盤 地 図情 報

国 土 地 理院

825000 分 の 1 地 図 (地 図画 像 )国 土 地 理院

被災前

9オ ル ソ 画像平 成 23 年 (2011 年 )東 北 地方 太平 洋 沖
地 震 正 射画 像 デー タ （オ ルソ画 像 ）宮 城 県 三陸 海 岸地 区 （ 3 月 19 日 撮 影分）

国 土 地 理院

2011 年 3 月 19日

平 成 23 年 (2011 年 )東 北 地方 太平 洋 沖
地 震 正 射画 像 デー タ （オ ルソ画 像 ）南 三 陸・女川 地 区（ 6 月 5 日～8 日 撮
影）

国 土 地 理院

2011 年 6 月 5
～8 日

No.
1
2
3

4

390

地 目 デ ータ 作 成に 使 用(道 路 )
地 目 デ ータ 作 成に 使 用

図３－10－⑦－２

新北上川周辺における被災前の地形データ

図３－10－⑦－３

新北上川周辺における被災後の地形データ

391

図３－10－⑦－４

新北上川周辺における土地利用図(c) 結論ならびに今後の課題
実用的な土砂移動シミュレーション手法の開発を目指し、これまでに開発されてきた津波土砂移動シミュレーションモデルの網羅的なレビューを行った。津波土砂移動モデ
ルの実用的かつ高度化に関する課題として、物理過程、すなわち土砂移動モデルの構成、計算条件、モデルの検証の３つに分類でき、それぞれの問題の解決には、今後の更なる研究の進展が必要であることを示し、実用的な津波土砂移動シミュレーション手法の開発に向けた問題点を抽出した。さらに、津波土砂移動モデルの検証を行うために，2011
年東北地方太平洋沖地震による津波で大規模な地形変化が生じた宮城県石巻市新北上
川河口部を対象とした今次津波前後の地形データの作成を行った。
これらの知見や離散化手法に関する検討を行い、土砂移動量に関する精度検証やその安定性などについて検討を行う必要がある。2011 年東北地方太平洋沖地震による宮城県気仙沼市気仙沼湾沿岸、同県石巻市新北上川河口部において、土砂移動に関する精度検証を実施する。さらに、1707 年宝永地震における沿岸湖沼の津波堆積物に関する再現解析を試みる。

(d) 引用文献
1) Adityawan, M.B., Tanaka, H., Bed stress assessment under solitary wave run-up. Earth
Planets Space 64, 945-954, 2012.
2) Apotsos, A., Buckley, M., Gelfenbaum, G., Jaffe, B., Vatvani, D., 2011a. Nearshore
Tsunami Inundation Model Validation: Toward Sediment Transport Applications. Pureand Applied Geophysics 168, 2097-2119, 2011a.
3) Apotsos, A., Gelfenbaum, G., Jaffe, B., Process-based modeling of tsunami inundation
392

and

sediment transport. J. Geophys. Res. 116, F01006. doi:10.1029/2010JF001797,

2011b.
4) Apotsos, A., Gelfenbaum, G., Jaffe, B., Watt, S., Peck, B., Buckley, M., Stevens, A.,
Tsunami inundation and sediment transport in a sediment-limited embayment on
American Samoa. Earth-Sci. Rev. 107, 1-11. doi:10.1016/j.earscirev.2010.11.001, 2011c5) 芦田和男，道上正規，移動床流れの抵抗と掃流砂量に関する基礎的研究，土木学会報告集，206，59-69，1972．
6) GALAPPATTI G., and VREUGDENHIL C. B., A depth-integrated model for suspended
sediment transport, J. Hydraul. Res., 23, 4, 359-377, 1985.
7) Goff, J., Chagué-Goff, C., Nichol, S.L., Jaffe, B. Dominey-Howes, D., Progress inpalaeotsunami research. Sedimentary Geology 243-244, 70–88, 2012.
8) Gusman, A. R., Tanioka, Y., Takahashi, T., Numerical experiment and a case study of
sediment transport simulation of the 2004 Indian Ocean tsunami in Lhok Nga, Banda
Aceh, Indonesia. Earth Planet Space 64, 817-827, 2012.
9) 林里美，越村俊一，映像解析による 2011 年東北地方太平洋沖地震津波の流速測定，土木学会論文集 B2（海岸工学），68，2-(1)，366-370，2012．
10) Kihara, N., Matsuyama, M.,

Numerical simulations of sediment transport induced by

the 2004 Indian Ocean tsunami near Kirinda port in Sri Lanka, Proceedings of 32nd
Conference on Coastal Engineering, Shanghai, China, 6p, 2010.
11) Kihara, N., Fujii, N., Matsuyama, M., Three-dimensional sediment transport process ontsunami-induced topography changes in a harbor, Earth Planets Space 64, 787-797, 2012.
12) Li, L., Qiu, Q., Huang, Z., Numerical modeling of the morphological change in Lhok
Nga, west Banda Aceh, during the 2004 Indian Ocean tsunami: understanding tsunami
depositsusingaforwardmodeling

method.

Nat.

Hazards,

DOI

10.1007/s/11069-012-0325-z, 2012.
13) Mellor, G. L., and T. Yamada, Development of a turbulence closure model for
geophysical fluid problems, Rev. Geophys. Space Phys, 20, 851-875, 1982.
14) Nakamura T., Mizutani, N., Yim, S.C., A Three-Dimensional Coupled Fluid-SedimentInteraction Model With Bed-Load/Suspended-Load Transport for Scour Analysis Arounda Fixed Structure, Journal of Offshore Mechanics and Arctic Engineering 131, 031104,
2009.
15) Ontowirjo, B., Paris, R., Mano, A., 2012. Modeling of coastal erosion and sedimentdeposition during the 2004 Indian Ocean tsunami in Lhok Nga, Sumatra, Indonesia, Nat.
Hazards, DOI 10.1007/s11069-012-0455-3
16) Paris, R., Fournier, J., Poizot, E., Etienne, S., Morin, J., Lavigne, F., Wassmer, P.,
Boulder and fine sediment transport and deposition by the 2004 tsunami in Lhok Nga(western Banda Aceh, Sumatra, Indonesia) - a coupled offshore-onshore model. Marine
Geology 268, 43–54, 2010.
17) Paris, R., Wassmer, P., Sartohadi, J., Lavigne, F., Barthomeuf, B., Desgages, E.,
393

Grancher, D., Baumert, P., Vautier, F., Brunstein, D., Gomez, C., Tsunamis as geomorphic
crises – lessons from the December 26, 2004 tsunami in Lhok Nga, West Banda Aceh
(Sumatra, Indonesia). Geomorphology 104, 59–72, 2009.
18) Ribberink, J.S., 1998. Bed-load transport for steady flows and unsteady oscillatory
flows, Coast .Eng. 34, 59-82.
19) Rubey, W.W., Settling velocities of gravels, sand and silt

particles, American Journalof Science, Vol.25, pp.325-338, 1933.
20) 菅原大助，今村文彦，松本秀明，後藤和久，箕浦幸治，地質学的データを用いた
西暦 869 年貞観地震津波の復元について，自然災害科学，29，4，501−516，2011.
21) Takahashi, T., Shuto, N., Imamura, F., Asai, D., Modeling sediment transport due totsunamis with exchange rate between bedload layer and suspended load layer. Proc. Int.
Conf. Coast. Eng. 2000, 1508-1519, 2000.
22) 高橋潤，後藤和久，大家隆行，柳澤英明，今村文彦，スリランカ・キリンダ漁港を対象とした 2004 年インド洋大津波による土砂移動過程の解析，日本土木学会海岸工学論文集， 55，251-255，2008．
23) Yoshii, T., Ikeno, M., Matsuyama, M., Fujii, N., Pick-up rate of suspended sand due totsunami. Proc. Int. Conf. Coast. Eng. 15p, 2010.
24) 行谷佑一，佐竹健治，山木滋，宮城県石巻・仙台平野および福島県請戸川河口

低地における 869 年貞観津波の数値シミュレーション，活断層・古地震研究報告，
No. 10, p. 9-29, 2010.
25) Van Rijn, L.C., Walstra, D.J.R., van Ormondt, M., Description of TRANSPOR 2004
(TR2004) and implementation in DELFT3D-online, Rep. Z3748, Delft Hydraul., Delft,
Netherlands, 2004.
26) Van Rijn, L.C., Principles of sediment transport in rivers, estuaries and coastal areas.
Aqua Publications, Amsterdam, 1993.
(3) 平成 26 年度業務計画案
日本列島広域構造モデルを継続開発し、地震地殻変動観測データに整合するようプレー
ト境界や地殻・マントルの摩擦・粘性パラメータ等の推定を行う。過去の地震とその後の粘弾性応答を考慮した、応力・強度分布の推定手法を構築する。津波に伴う土砂移動を評価できる津波浸水計算コードを継続開発し、東北地方太平洋沖地震の浸水・土砂移動シミュレーションを実施して検証を行うとともに、南海トラフ地震における龍神池（大分県佐伯市）等の津波池の津波堆積物と津波高の関係を検討する。南海トラフ地震の地震津波の史料調査と地震津波シミュレーションに基づき、過去の南海トラフ地震（慶長地震、宝永、
安政地震など）の震源像を明確化する。フィリピンプレート形状と３次元地殻・堆積層構
造モデルを整備して、近年の大地震の強震動・長周期地震動生成メカニズムを検証するとともに、地下構造の不均質性が地震動に与える影響に感度解析を行い、地震動・津波シミュレーションの高精度化・広帯域化を進める。

394