学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

jh150039-NA23

グリーンランド氷床モデルを用いた 2.5 次元理論地震波形計算
豊国源知（東北大学）

概要 北極グリーンランドでは，2011 年から氷床上での本格的な地震観測がスタートした．申請者も観測隊の一員として，5 年間毎年現地で地震計のメンテナンス作業を続けてきた．現在はリアルタイムでのデータ転送も行われ，氷床の影響を受けた地震波形データが蓄積され始めている．地震波形に現れた氷床の影響は，波形を利用して地下構造を解析する際にはノイズとなる一方，氷床そのものの内部構造や，氷床の状態の時間変化を調べる手がかりにもなる．本研究課題は，効率の良い理論地震波形計算により，様々な構造断面・震源位置・発震機構について，氷床がどのように観測地震波形に影響を与えるかを調べ，その生成メカニズムを明らかにすることを目的とした．本課題によって基本的なプログラム開発が終了し，東北大学サイバーサイエンスセンターの SX-ACE1
ノード 16 プロセッサを用いた自動並列/OpenMP 演算に成功した．震源・構造モデルを様々に変えた複数のシミュレーションを実行し，氷床の影響を受けた地震波形を比較することで，氷床内を多重反射しながら伝播する地震波相を発見し，「Le 波」と命名した．

1. 共同研究に関する情報のことで，グリーンランドと南極にしか存在しな(1) 共同研究を実施した拠点名

い．グリーンランド氷床の厚さは平均で 2 km，最厚部では 3 km にも達し，すべてが融解した場合，東北大学サイバーサイエンスセンター

海水準が 7 m 以上も上昇すると予測されている．

(2) 共同研究分野

近年は地球温暖化の影響により，氷床の融解が指■超大規模数値計算系応用分野

摘されていることから，融解の状況をモニタリン



超大規模データ処理系応用分野

グする方法として，地震観測に期待が寄せられる

超大容量ネットワーク技術分野こととなった．融解に伴い，氷が割れたり滑った

超大規模情報システム関連研究分野りする際に生じる振動は，地震計で即時に検出で(3) 参加研究者の役割分担

研究代表者：豊国源知（東北大・助教）：コードの開発と演算の実行，結果の解釈
共同研究者：竹中博士（岡山大・教授）：
コードの開発共同研究者：江川隆輔（東北大・準教授）：
高速化・最適化支援
共同研究者：小松一彦（東北大・助教）：
高速化・最適化支援

2. 研究の目的と意義
グリーンランドは北極圏に位置する世界最大の島で，表面の 80%以上は氷床に覆われている．氷図 1：GLISN 地震観測点の分布．赤色は日本隊が設

床とは，地表を覆う全面積が 5 万 km2 以上の氷塊置・メンテナンスを行った点．黄色はその他の点．

1

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

きるからである．しかしグリーンランドの内陸部は，人を寄せ付けない過酷な環境にあるため，2008 年の時点では，周辺の島々を合わせても地震観測点は 16 点，うち氷床上はわずか 1 点という状況で，観測網の拡充が望まれていた．
こうした経緯で，2009 年に「グリーンランド氷床モニタリング観測網（GreenLand Ice Sheetmonitoring Network，略称 GLISN）」と呼ばれる
図 2：2013 年 2 月 19 日にグリーンランド南端で発国際プロジェクトが発足し，観測点の増強が進め生した地震（深さ 10 km，mb=4.8）の，3 観測点にられてきた．現在は 11 ヶ国が参加し，氷床上 4おける上下動記録．左に観測点名と氷床上・露岩上点を含む計 33 点の地震観測点が運用されているの別，震央距離（Δ）を示した．観測点の位置は図 1（図 1）．日本はプロジェクト発足時からの参加国

参照．すべての波形に 0.5‐833 s のバンドパスフィで，2011 年以降は毎年アメリカ隊との共同観測隊ルターを適用．群速度 3.1‐2.6 km/s の波群を太線を派遣して，氷床上 3 観測点と沿岸域 3 観測点ので示す．

設置やメンテナンスに携わってきた．氷床上の観測点は維持に多大な困難が伴うが，良質な観測網の構築には均一な観測点分布が不可欠であり，日震源決定，氷床内部の構造とその時間変化の解析

米共同観測隊の貢献はプロジェクト全体でも極め等に利用できる可能性も秘めている．そこで本研て大きいといえる．2014 年には衛星通信を用いた究課題では，グリーンランド氷床によって励起さ氷床上からの即時データ転送にも成功し，極地でれた地震波動場とその励起メカニズムを，現実的観測された地震波形を，世界中ですぐさまダウンな地形，氷床厚，地殻・マントル構造モデルを用ロードできるようになった[例えば，Toyokuni etいた理論地震波形計算で明らかにすることを目的

al., 2014, Antarctic Record]．とする．

GLISN 観測網から得られる地震波形データは，GLISN 観測網では，広帯域 3 成分地震計を用いリアルタイムで国際大学共同機関のデータベースた連続観測を行っているため，得られるデータのにアップロードされ，解析に供されている．しか情報量は極めて多く，氷床融解の検出以外の目的

し氷床を通過した地震波の波形は，氷床による反にも利用することができる．例えば，グリーンラ射波や変換波によって，通常の観測波形と大きくンドは詳細な地下構造が明らかされておらず，乖離した特徴が見られることが知られている．図

GLISN データの解析による研究が待たれている．

2 で示した例はグリーンランド南端の浅部（深さデータ解析の際には，氷床が観測波形に与える影10 km）で発生した地震の波形であり，氷床の影響がノイズとなることが予想されるため，解析に響により (1) P 波や S 波の立ち上がりが不明瞭，先立って氷床の影響を正しく見積もっておく必要(2) 周期数秒で小振幅の鋸歯状の波が現れる，(3)がある．しかし氷床中を伝播する地震波に関する

S 波の後続波部分に大振幅の波群（点線枠部分）

先行研究は，1960 年代に 1 篇が見られるのみでが現れる，といった特徴が見られる[Toyokuni et[Robinson, 1968, JGR]，数値計算等で理論的に取al., 2014, Antarctic Record]．

り扱った研究はいまだ行われていない．

このように氷床の影響を受けた地震波形を地下構造の解析等に利用する場合は，十分に氷床の影3. 当拠点公募型共同研究として実施した意義

響を除去しないと精度良い解析は望めない．一方

本課題で利用する計算手法は，申請者がこれまで，氷床による波形の大きな変形は，浅い地震ので開発を行ってきた，ローカルな地震波伝播を精2

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

度と効率よく計算する手法「円筒座標系 2.5 次元VP, VS)=(0.914g/cm3, 4.0 km/s, 2.0 km/s)を用い差分法」である[Takenaka et al., 2003, GRL;

た．また空気層は(ρ, VP, VS)= (0.00129g/ cm3, 0.0

Toyokuni et al., 2012, InTech; Toyokuni et al.,km/s, 0.0 km/s)とした．空気層を導入すると，地

2015, PEPI, 投稿中]．2.5 次元計算とは，媒質の震波から音波への変換波もシミュレーションでき

2 次元構造断面上で 3 次元の波動場を計算する方るが，音速は岩石中の地震波速度に比べて極端に法であり，通常の 3 次元計算の 1 万～10 万分の遅く，精度良くシミュレーションするためには格1 の計算資源で現実的な波動場を計算できる．本子間隔を細かく切らなければならない．今回の計研究の特色は，手法の効率の良さを活かし，様々

算は音波を対象としたものではないので，空気層な構造・震源モデルに適用し多数の計算を実行すの地震波速度をゼロとすることで，音波の伝播をることで，グリーンランド氷床が観測波形に与え抑制した．

る影響を網羅的に捉えられる点にある．

計算は東北大学サイバーサイエンスセンターの今回の課題で実行した演算は，SX-ACE の 4 ノSX-9 のジョブクラス p8，および SX-ACE のジョード 16 プロセスを利用した自動並列/OpenMP 演

ブクラス p16 による自動並列/OpenMP 演算で，算である．計算規模の例としては，グリーンラン

水平（r）700 km，鉛直（z）150 km の構造断面

ドのある構造断面を切り出した水平 700 km，深を，14000×3000 の空間差分格子に分割して行っ

さ 150 km の領域について，水平 14,000×深さた．格子間隔は水平・鉛直ともに 0.05 km である．

3,000 格子に切り，卓越周期 2 Hz までの地震波伝

時間刻みは 0.0025 s とし，励起後 150 s まで計算

播シミュレーションを行う計算が挙げられる．計した．震源にはピュアな縦ずれ断層型のメカニズ

算結果は全ての空間格子上で書き出さなければな

ムを用い，震源時間関数として幅 2 Hz のベル型らず，ファイル容量は格子数に比例して大きくなパルスを入力した．図 4 は，4 つのシミュレーシるため，CPU 時間とファイル使用量に比例してョンで得られた理論波形を，3 つの震央距離（Δ=50，課金される現在の環境では，目標とする計算は実100，150 km）について並べて表示した図である．

質不可能であり、本制度による使用量の補助が不以下ではこれらの理論波形を用いて，4 つのシミ可欠である．

ュレーションの結果について述べる．

4. 前年度までに得られた研究成果の概要① FLAT2.0-0
氷床が厚さ 2 km 一定の水平成層構造と仮定し，本課題は継続課題ではないので該当しない．

震源を氷床直下に置いた場合，極めて継続時間の長い S 波後続波が現れた．氷床を入れない場合（図5. 今年度の研究成果の詳細
今回は表 1 と図 3 で示す 4 つのシミュレーショ

4 黒線），地殻の中を通過した S 波（Sg）と表面波ンを行った．計算のメインは，実際の地形の起伏による単純なパルスの到着後，振幅はほぼゼロにや氷床厚分布を考慮した TOPO-0 および TOPO-5なっている．一方，氷床を入れた場合（図 4 赤線），であるが，比較対象として，氷床を厚さ一定の水Sg 到着後，数 10 秒以上大振幅が継続する．大振平 成 層 構 造 と し て モ デ ル 化 し た FLAT2.0-0 ，幅の継続時間は震央距離に比例し，Δ=150 km のFLAT0.5-0 でもシミュレーションを行った．氷床地点では 80 s にも達している．
氷床を入れた場合の継続時間の長い S 波後続波の影響を調べるため，氷床部分を岩石で置き換えは，氷床内部に地震波のエネルギーが強くトラッた計算も行った．
構造モデル作成の際，氷床以外の部分は 1 次元プされることで起きる現象である．氷床直下に与地球モデルを用いた．氷床の密度（ρ），P 波速度

えられた震源から放出された地震波は，エネルギ（VP），S 波速度（VS）は，先行研究等を参考に(ρ,ーの大部分が氷床に入射し，薄く低地震波速度の3学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

表 1：4 つのシミュレーションの諸元
識別名

氷床モデル

氷床表面からの震源深さ（km） 氷床基盤からの震源深さ（km）

FLAT2.0-0

水平均一厚（厚さ 2.0 km）

2.20.2

FLAT0.5-0

水平均一厚（厚さ 0.5 km）

0.70.2TOPO-0

図 3(c)

0.60.2

TOPO-5図 3 の A-A’断面

5.65.2では，氷床を意味するドイツ語「Eisdecke」の頭文字から，氷床内を多重反射する S 波を「Le 波」と命名した．同様に，地殻表面を伝播する表面波（レイリー波）は「Rg 波」と呼ばれることから，氷床表面を伝播する表面波を「Re 波」と命名した[Toyokuni et al., 2015, Polar Sci.]．
②FLAT0.5-0
氷床が厚さ 0.5 km 一定の水平成層構造と仮定
し，震源を氷床直下に置いた場合，前節の結果と似た継続時間の長い S 波後続波が現れた．図 4 の赤線が示すように，S 波後続波の継続時間はFLAT2.0-0 とほとんど変わらず，氷床の厚さには依存しないと考えられる．一方，周波数は前節の
図 3：4 つのシミュレーションで用いた構造モデルと震源位場合よりも高周波が卓越している．波がトラップ置の違い．(a)～(c)は計算に用いた構造断面．黄色が空気層，される過程は前節の結果と同様であるが，氷床が青色が氷床，赤色は岩石層．地表付近を鉛直方向に拡大し薄いことで波の成長が阻害され，高周波成分が卓て表示．星印は震源．(a) 水平成層，厚さ 2 km の氷床モデ越したと考えられる．

ルの直下に震源を入力（FLAT2.0-0），(b) 水平成層，厚さ0.5 km の氷床モデルの直下に震源を入力（FLAT0.5-0），(c)② TOPO-0

現実的な地形の起伏と氷床厚分布を使用し，震源は氷床直現実的な地形・氷床厚分布を用い，震源を氷床下（TOPO-0）または氷床下 5 km（TOPO-5）に置いた．(d)直下に置いた場合は，水平成層の場合と大きく異構造断面の A-A’測線の位置を示す平面図．背景色は氷床厚．

なる波形が得られた．S 波後続波の大振幅は Sgの到着からかなり遅れて現れるうえ，継続時間は氷床内部を多重反射しながら伝播する．氷床と短く，波群がコンパクトにまとまっている（図 4

基盤岩との境界が水平な場合，全反射が卓越する赤線）．このケースも震源が浅いため，地震波のエことで，地震波のエネルギーは氷床から漏れ出しネルギーの大部分が氷床内に入射する点は前節まにくく，継続時間の長い S 波後続波が形成される．

でと同じである．しかし地形や基盤岩の凹凸によ低地震波速度の層内で地震波がトラップされる現って P 波から S 波，S 波から P 波への変換が起こ

象としては，地殻内部を多重反射する S 波であるり，エネルギーが分散されることに加え，波の全「Lg 波」がよく知られている．一方，氷床内部の反射が崩れてエネルギーは次々と地殻内に漏れ出多重反射は，P 波の多重反射の観測事例があるもしていく．結果として，氷床内をほぼ水平に伝播のの，理論的研究は本研究が最初である．本研究する波が選択的に残され，コンパクトな Le 波の4

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

図 4：4 つのシミュレーションによる，震央距離Δ=50，100，150 km における上下動成分の理論波形．震央距離とシミュ

図 5：TOPO-0 による地震波伝播の様子．上から励起後 10，30，50，70，90 s のスナップショット．赤色は P 波，緑色

レーションの識別名はパネルの左に表示した．赤線は氷床は S 波を表す．色が濃いほど振幅が大きい．震源の位置はあり，黒線は氷床なし．主要なフェーズ（Pg，Sg など）の星印．Pg，Sg は地殻中を伝播する P 波と S 波，Pn はマン到着時刻を青点線で示した．太線は図 2 に対応する群速度

トル中を伝播する P 波，Lg は地殻内部にトラップされた S

3.1‐2.6 km/s の波群．

波，Le は氷床内部にトラップされた S 波，Re は氷床表面を伝播する表面波．

波群が現れる．氷床から S 波エネルギーが漏れ出していることは，図 5 のスナップショットで，地幅や後続波の継続時間の増大，といった現象に現殻内に S 波を示す緑色が卓越していることから明れている．これは柔らかい地盤で地震動が増幅さ瞭に確認できる．

れる現象と同じメカニズムである．図 4 の太線で示した部分は，群速度 3.1‐2.6 km/s の波群で，
図 2 の観測波形の太線部分と対応している．今回③ TOPO-5
現実的な地形・氷床厚分布を用い，震源を氷床の計算例は構造断面が観測のものと異なるため波の下 5 km に置いた場合は，明瞭な Le 波が現れな形の直接比較はできないが，震源が氷床直下になかった．このケースでは，震源が深いことで，そい場合の波形の特徴は理論波形でよく再現できたもそも氷床内部に入射するエネルギーが少ない．といえる．

また地形・基盤岩の凹凸によって氷床内へのエネ

6. 今年度の進捗状況と今後の展望

ルギーのトラップも弱められている．従って地殻内を伝播する波が卓越し，氷床内を伝播する波は

本課題では，現実的な氷床厚モデルと複数の震不明瞭となる（図 4 赤線）．しかし氷床の影響は，源位置やメカニズムを組み合わせて大量の計算を氷床がない場合（図 4 黒線）に比べた Sg 波の増実行することで，氷床の影響を受けた特殊な波動5

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

場「Le 波」の存在を理論的に示した．今後は観測(4) 国内会議発表

波形との比較をさらに進め，氷床の情報をリアル

豊国 源知・竹中 博士,・金尾 政紀・坪井 誠司・

タイムで抽出する手法の開発・検討を行っていき東野 陽子:たいと考えている．

グリーンランド氷床モデルを用いた地震波伝播モ
デリング, 日本地球惑星科学連合, 5 月 22 日,

一方，当初はプログラムの MPI 化も行う予定で2015.

あったが，地形や非弾性減衰の導入を優先したた
め，年度内には実行できなかった．これは今後の(5) その他（特許，プレス発表，著書等）

課題として取り組んでいきたい．

豊国 源知・竹中 博士,・金尾 政紀・坪井 誠司・

7.

研究成果リスト

東野 陽子:(1) 学術論文

グリーンランド氷床モデルを用いた地震波伝播シ

Toyokuni, G., H. Takenaka, M. Kanao, S. Tsuboi,

ミュレーション, SENAC, 48(4), 10-23, 2015.

and Y. Tono:
Numerical modeling of seismic waves for

豊国 源知・竹中 博士,・金尾 政紀・坪井 誠司・estimating the influence of the Greenland ice

東野 陽子:

sheet on observed seismograms, Polar Sci., 9(1),

氷床によって形成される地震波相「Le」と「Re」

80-93, 2015.の提唱, 月刊地球, 37(11), 463-473.

Toyokuni, G., H. Takenaka, and D. Zhao:
Quasi-cylindrical 2.5-D wave modeling with ageneralpointsourceandtheanelasticattenuation, Phys. Earth Planet. Int., 2016, 投稿中．
(2) 国際会議プロシーディングス
Toyokuni, G., H. Takenaka, M. Kanao, S. Tsuboi,
and Y. Tono:
Seismic waveform modelings for estimating
influenceof

theicesheetonobserved

seismograms, Abstract S23B-2705, 2015 Fall
Meeting, AGU, 2015.
(3) 国際会議発表
Toyokuni, G., H. Takenaka, M. Kanao, S. Tsuboi,
and Y. Tono:
Seismic waveform modelings for estimating
influenceof

theicesheetonobserved

seismograms, 2015 Fall Meeting, AGU, San
Francisco, Calif., 14 Dec., 2015.
6