— 10 —SENAC Vol. 48, No. 4（2015. 10）［共同研究成果］

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

豊国源知 1・竹中博士 2・金尾政紀 3・坪井誠司 4・東野陽子 5
1 東北大学大学院 理学研究科 地震・噴火予知研究観測センター
2 岡山大学大学院 自然科学研究科 地球生命物質科学専攻
3 国立極地研究所
4 海洋研究開発機構
5 文部科学省
平成 26 年度の共同研究課題では，局地的な地震波伝播を効率良くモデリングできる数値計算手法を開発した．また本手法による計算例として，北極グリーンランド氷床を対象とした地震波
伝播シミュレーションを行った．グリーンランドでは 11 ヶ国の国際協力で，2011 年から本格的な地震観測が進められてきた．特に氷床上の観測点は，日米共同観測隊が毎年現地で行うメンテ
ナンスによって維持されており，得られたデータの解析が鋭意進められている．解析の際には，観測地震波形に与える氷床の影響を正確に考慮する必要があるが，これまで研究はほとんど行われておらず，本研究がその端緒を開くと期待される．本稿では，我々の計算手法とシミュレーシ
ョン成果の概要を紹介する．

1.

はじめに

グリーンランドは北極圏に位置する世界最大の島で，表面の 80%以上は氷床に覆われている．
氷床とは，地表を覆う全面積が 5 万 km2 以上の氷塊のことで，グリーンランドと南極にしか存在しない．グリーンランド氷床の厚さは平均で 2 km，最厚部では 3 km にも達し，すべてが融解した場合，海水準が 7 m 以上も上昇すると予測されている．近年は地球温暖化の影響により，氷床の融解が指摘されていることから，融解の状況をモニタリングする方法として，地震観測に期待が寄せられることとなった．融解に伴い，氷が割れたり滑ったりする際に生じる振動は，地震計で即時に検出できるからである．しかしグリーンランドの内陸部は，人を寄せ付けない過酷な環境にあるため，2008 年の時点では，周辺の島々を合わせても地震観測点は 16 点，うち氷床上はわずか 1 点という状況で，観測網の拡充が望まれていた．
こうした経緯で，2009 年に「グリーンランド氷床モニタリング観測網（GreenLand Ice Sheet
グ リ ス ン
monitoring Network，略称GLISN）」と呼ばれる国際プロジェクトが発足し，観測点の増強が進められてきた．現在は 11 ヶ国が参加し，氷床上 4 点を含む計 33 点の地震観測点が運用されている（図 1）．日本はプロジェクト発足時からの参加国で，2011 年以降は毎年アメリカ隊との共同観測隊を派遣して，氷床上 3 観測点と沿岸域 3 観測点の設置やメンテナンスに携わってきた．氷床上の観測点は維持に多大な困難が伴うが，良質な観測網の構築には均一な観測点分布が不可欠であり，日米共同観測隊の貢献はプロジェクト全体でも極めて大きいといえる．2014 年には衛星通信を用いた氷床上からの即時データ転送にも成功し，極地で観測された地震波形を，世界中ですぐさまダウンロードできるようになった[1-5]．
GLISN 観測網では，広帯域 3 成分地震計を用いた連続観測を行っているため，得られるデー
タの情報量は極めて多く，氷床融解の検出以外の目的にも利用することができる．例えば，グリーンランドは詳細な地下構造が明らかされておらず，GLISN データの解析による研究が待たれている．データ解析の際には，氷床が観測波形に与える影響がノイズとなることが予想されるため，解析に先立って氷床の影響を正しく見積もっておく必要がある．しかし氷床中を伝播する地震波に関する先行研究は，1960 年代に 1 篇が見られるのみで[6]，数値計算等で理論的に取り扱った研究はいまだ行われていない．

グリーンランド氷床モデルを用いた地震波伝播シミュレーション— 11 —

図 1：（左）GLISN 地震観測点の分布．
赤色は日本隊が設置・メンテナンスを行った点．黄色はその他の点．（右）日本隊が設置を行った ICESG 観測点の外観（2011 年）．

平成 26 年度の共同研究課題では，局地的な地震波伝播を効率よくシミュレーションすることができる数値計算手法を開発した．またこの手法に氷床を含むグリーンランドの地下構造モデルを入力してシミュレーションを行い，観測地震波形に氷床が与える影響を調べた．以降では本年出版した論文[7]の流れに沿いながら，手法と結果について述べていく．

2.

円筒座標系 2.5 次元モデリング

地震波伝播シミュレーションでは，ある地下構造モデルと震源モデルを与えて，地震波伝播の支配方程式を数値的に解くことで，任意の地点で観測される地震波形を理論的に予測する．モデリング手法には様々なものが存在するが，計算精度と効率の良さを併せ持っている手法は少ない．
本誌 46 巻 3 号で紹介した「軸対称 2.5 次元モデリング」は，その代表的な手法である[8]．この手法では，震源を通る鉛直軸の周りに構造の軸対称性を仮定することで，3 次元的な空間の広がりの効果を 2 次元断面に落とし込んだのち，その断面内において，面内運動 2 成分と，面外運動
1 成分の，3 次元の波動伝播を計算する．
「波動場 3 次元，構造 2 次元」という意味で，「2.5 次元」という呼称が用いられている．計算は断面上だけで行うので計算資源を節約でき，計算時間やメ
モリは 3 次元でフルに計算する場合の数千～数万分の 1 で済む．また 3 次元の波動場が得られる
ので，観測地震波形と直接比較できるというメリットを持つ．計算では最初に回転対称軸を与える必要があることから，座標系としては球座標系と円筒座標系が利用される．
我々はこれまで，地球や他の天体の全球地震波伝播をシミュレーションするため，球座標系を使ったプログラム開発を行ってきた[8-13]．しかしグリーンランドのような局地的な地震波伝播をシミュレーションする場合には，地球の中心まで含む計算を行うとメモリや計算時間が余分にかかってしまうことから，地球表層に特化したプログラム開発が必要であった．
円筒座標系��� �� ��において，�軸を鉛直下方，�軸を水平方向にとって，方位角を�で表すと，� � 0を対称軸とした��断面上で 2.5 次元計算が可能となる．��座標は，�軸が� � 0で反転することを除けば 2 次元デカルト座標系と同じ扱いができるので，地表面を� � 0とおいた矩形断面上で計算が行える．矩形断面では地表の曲率を考慮できないので，シミュレーションの対象としては，地表を平面として取り扱える，局地的（伝播距離≲ 1000 km）な地震波伝播が適している．円筒座標系を用いた 2.5 次元計算は，1960 年代から局地的な地震波伝播シミュレーションに活用されてきた[14-17]．
しかし上述の�軸が反転する問題は，従来のシミュレーションに深刻な影響を与えてきた．すな
わち，ある構造断面を��平面上で表そうとすると，� � 0の軸を挟んで，� � 0と� � �における 2

— 12 —

SENAC Vol. 48, No. 4（2015. 10）

図 2：（a）通常の円筒座標領域と（b）準円筒
座標領域[18]．星印は震源．

枚の矩形を用いなければならない．このため任意の構造モデルを߶ ൌ Ͳの面に入力すると，軸対称性によって߶ ൌ ߨの面における構造は߶ ൌ Ͳのものと対称にならざるを得ず，現実的な地下構造を取り扱うことができない（図 2a）．また対称軸である‫ ݎ‬ൌ Ͳが，剛体の壁と同じ働きをして，そこで人工的な反射波を生じさせる問題もあった．これらのデメリットは，円筒座標領域が通常（െλ ൏ ‫ ݖ‬൏ λǡ Ͳ ൑ ‫ ݎ‬൏ λǡ െߨ ൑ ߶ ൑ ߨ）という範囲を用いていたことに起因している．
これに対して文献[18]では，準円筒座標領域（െλ ൏ ‫ ݖ‬൏ λǡ െλ ൏ ‫ ݎ‬൏ λǡ െߨȀʹ ൑ ߶ ൑ ߨȀʹ）が提案された．このような範囲を用いると，構造断面は߶ ൌ Ͳにおける 1 枚の矩形のみで表され，この面には任意の不均質構造モデルが入力できる（図 2b）．座標軸の取り方を変えるだけなので，従来の軸対称モデリングのメリットは保存したまま，非対称な構造モデルが取り扱える．文献[18]は爆破地震探査を想定した研究であったことから，計算では軸対称な震源のみを扱っていた．我々はこの手法をさらに拡張し，断層を含む任意の震源メカニズムと任意形状の地形の導入も行った[7]．これにより，現在のプログラムでは極めて現実に近い地震波伝播をシミュレーションできる．
手法の詳細については，文献[19]も参照されたい．
円筒座標系を使った地震波伝播シミュレーションで基本となる方程式は，3 次元の弾性波の支配方程式である．ここで‫ݐ‬は時間，ߩは密度，ߨとߤは弾性定数，‫ݒ‬௜ ，݂௜ ，ߪ௜௝ ，‫ܯ‬௜௝ はそれぞれ粒子速度，体積力，応力テンソル，モーメントテンソルの成分である（݅ǡ ݆ ‫ א‬ሼ‫ݖ‬ǡ ‫ݎ‬ǡ ߶ሽ）．変数の上についたドットは時間微分を表す．弾性定数としてはラメ定数ߣǡ ߤを用いるのが一般的であるが，ߨ ‫ ߣ ؠ‬൅ ʹߤを導入すると，P 波の性質を単独の弾性定数で表せるメリットがある．
我々の開発したプログラムでは，この方程式に現れるすべての物理量を߶方向にフーリエ級数
展開したのち，差分法で離散化して計算を行う．フーリエ級数展開によって，震源メカニズムの߶— 13 —

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

図 3：計算で用いた差分格子の模式
図．黒が粒子速度，
白が応力の格子を表す．星印は震源．

依存性が考慮でき，任意のメカニズムによる波動場の励起が可能となる．図 3 に差分法で用いる格子配置の模式図を示した．気体・液体を含む媒質でも安定に計算できるよう，粒子速度と応力の別々の格子に配置される「スタガード格子」を採用している．震源は必ず� = �の軸上に置かなければならないので，水平方向の震源位置を変えると�座標も変化する．

3.

グリーンランド氷床モデルによるシミュレーション

本節では円筒座標系 2.5 次元モデリングによるグリーンランドの地震波伝播シミュレーションの成果について述べる．
図 4 はグリーンランドの南端，深さ 10 km で発生した地震（実体波マグニチュード�� = 4.8）の上下動記録である．GLISN 観測網の露岩上 1 観測点（NRS），氷床上 2 観測点（DY2G ，ICESG）の波形を示している．氷床上の観測波形には，S 波の後続波部分に群速度 3.1‐2.6 km/s の特徴
的な振幅増幅が見られることがわかる．氷床に関連すると思われるこのような特徴的な波形の成因について調べるため，今回は表 1 と図 5 で示す 4 つのシミュレーションを行った．
図 4：2013 年 2 月 19 日にグリーンランド南端で発生した地震（深さ 10 km，�� = 4.8）の，3 観測点における上下動記
録 [7] ． 左 に 観 測 点 名 と 氷 床上・露岩上の別，震央距離（Δ）
を示した．観測点の位置は図 1
参照．すべての波形に 0.5‐833 s のバンドパスフィルタ
ー を 適 用 ． 群 速 度 3.1 ‐ 2.6
km/s の波群を太線で示す．
表 1：4 つのシミュレーションの諸元
識別名

氷床モデル

FLAT2.0-0

水平均一厚（厚さ 2.0 km）

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

0.2

FLAT0.5-0

水平均一厚（厚さ 0.5 km）

0.7

0.2

TOPO-0

図 4 の A-A’断面

0.6

0.2

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

5.6

5.2

— 14 —

SENAC Vol. 48, No. 4（2015. 10）

図 5：4 つのシミュレーションで用いた構造モデルと震源位置の違い（文献[7]の図を修正）．(a)～(c)は計算に用いた構造断面．黄色が空気層，青色が氷床，赤色は岩石層．地表付近を鉛直方向に拡大して表示．星印は震源．(a) 水平成層，厚さ 2 km の氷床モデルの直下に震源を入力（FLAT2.0-0），(b) 水平成層，厚さ 0.5 km の氷床モデルの直下に震源を入力（FLAT0.5-0），(c) 現実的な地形の起伏と氷床厚分布を使用し，震源は氷床直下（TOPO-0）または氷床下 5 km
（TOPO-5）に置いた．(d) 構造断面の A-A’測線の位置を示す平面図．背景色は氷床厚．
計算のメインは，実際の地形の起伏や氷床厚分布を考慮した TOPO-0 および TOPO-5 であるが，比較対象として，氷床を厚さ一定の水平成層構造としてモデル化した FLAT2.0-0，FLAT0.5-0
でもシミュレーションを行った．氷床の影響を調べるため，氷床部分を岩石で置き換えた計算も行ったので，トータルの計算回数は 8 回である．構造モデル作成の際，氷床以外の部分は 1 次元
地球モデルである PREM[20]を用いた．TOPO-0，TOPO-5 で用いた地形や氷床厚分布は，文献，S 波速度（�� ）は，先行研究[23-25]等を参[21,22]に基づく．氷床の密度（�），P 波速度（�� ）�
考に(�, �� , �� ) = (0.914g/cm , 4.0 km/s, 2.0 km/s)を用いた．また空気層は(�, �� , �� ) = (0.00129g/
cm� , 0.0 km/s, 0.0 km/s)とした．空気層を導入すると，地震波から音波への変換波もシミュレー
ションできるが，音速は岩石中の地震波速度に比べて極端に遅く，精度良くシミュレーションするためには格子間隔を細かく切らなければならない．今回の計算は音波を対象としたものではな
いので，空気層の地震波速度をゼロとすることで，音波の伝播を抑制した[26-28]．
計算は SX-9 のジョブクラス p8 による自動並列演算で，水平（�）700 km，鉛直（�）150 kmの構造断面を，14000×3000 の空間差分格子に分割して行った．格子間隔は水平・鉛直ともに0.05 km である．時間刻みは 0.0025 s とし，励起後 150 s まで計算した．震源にはピュアな縦ず
れ断層型のメカニズムを用い，震源時間関数として幅 2 Hz のベル型パルスを入力した．図 6 は，
4 つのシミュレーションで得られた理論波形を，3 つの震央距離（Δ =50，100，150 km）について並べて表示した図である．また図 7～11 では，励起後 10，30，50，70，90 s における構造断面上の地震波伝播のスナップショットを並べている．赤色が P 波，緑色が S 波で，色が濃いほど振幅が大きい．以下ではこれらの理論波形とスナップショットを用いて，4 つのシミュレーショ
ンの結果について述べる．

グリーンランド氷床モデルを用いた地震波伝播シミュレーション— 15 —

図 6：4 つのシミュレーションによる，震央距離Δ =50，100，150 km における理論波形[7]．
（左）上下動，（右）水平動．震央距離とシミュレーションの識別名はパネルの左に表示した．
赤線は氷床あり，黒線は氷床なし．主要なフェーズ（Pg，Sg など）の到着時刻を青点線で示した．太線は図 4 に対応する群速度 3.1‐2.6 km/s の波群．

3.1 FLAT2.0-0
氷床が厚さ 2 km 一定の水平成層構造と仮定し，震源を氷床直下に置いた場合，極めて継続時
間の長い S 波後続波が現れた．氷床を入れない場合（図 6 黒線），地殻の中を通過した S 波（Sg）と表面波による単純なパルスの到着後，振幅はほぼゼロになっている．一方，氷床を入れた場合（図 6 赤線），Sg 到着後，数 10 秒以上大振幅が継続する．大振幅の継続時間は震央距離に比例し，Δ = 150 kmの地点では 80 s にも達している．
図 7 は氷床を入れない場合，図 8 は入れた場合の地震波伝播のスナップショットである．図 7と図 8 の比較から示されるように，氷床を入れた場合の継続時間の長い S 波後続波は，氷床内部に地震波のエネルギーが強くトラップされることで起きる現象である．氷床直下に与えられた震源から放出された地震波は，エネルギーの大部分が氷床に入射し，薄く低地震波速度の氷床内部を多重反射しながら伝播する．氷床と基盤岩との境界が水平な場合，全反射が卓越することで，地震波のエネルギーは氷床から漏れ出しにくく，継続時間の長い S 波後続波が形成される．
低地震波速度の層内で地震波がトラップされる現象としては，地殻内部を多重反射する S 波である「Lg 波」がよく知られている．一方，氷床内部の多重反射は，P 波の多重反射の観測事例が文献[6]で示されているものの，理論的研究は本研究が最初である．本研究では，氷床を意味するドイツ語「Eisdecke」の頭文字から，氷床内を多重反射する S 波を「Le 波」と命名した．同様に，地殻表面を伝播する表面波（レイリー波）は「Rg 波」と呼ばれることから，氷床表面を伝播する表面波を「Re 波」と命名した[7]．

— 16 —

SENAC Vol. 48, No. 4（2015. 10）

図 7：FLAT2.0-0 で氷床がない場合の地震波伝播の様子．図 5(a)の氷床モデルを直下の地殻構
造で置き換えて計算した．上から励起後 10，30，50，70，90 s のスナップショット．赤色は P
波，緑色は S 波を表す．色が濃いほど振幅が大きい．構造断面のうち，深さ 30 km 以浅の部分を拡大して表示した．震源の位置は星印．Pg，Sg は地殻中を伝播する P 波と S 波，Pn，Snはマントル中を伝播する P 波と S 波，Lg は地殻内部にトラップされた S 波，Rg は地殻表面を伝播する表面波．

グリーンランド氷床モデルを用いた地震波伝播シミュレーション— 17 —

図 8：FLAT2.0-0 による地震波伝播の様子（文献[7]の図を修正）．上から励起後 10，30，50，
70，90 s のスナップショット．赤色は P 波，緑色は S 波を表す．色が濃いほど振幅が大きい．
構造断面のうち，深さ 30 km 以浅の部分を拡大して表示した．震源の位置は星印．Pg，Sg は地殻中を伝播する P 波と S 波，Pn はマントル中を伝播する P 波，Lg は地殻内部にトラップされた S 波，Le は氷床内部にトラップされた S 波，Re は氷床表面を伝播する表面波．図 7 と比
較すると，トラップ波によって形成された地表付近の大振幅が，長い水平距離にわたって続いていることがわかる．

— 18 —

SENAC Vol. 48, No. 4（2015. 10）

図 9：FLAT0.5-0 による地震波伝播の様子（文献[7]の図を修正）．上から励起後 10，30，50，
70，90 s のスナップショット．赤色は P 波，緑色は S 波を表す．色が濃いほど振幅が大きい．
構造断面のうち，深さ 30 km 以浅の部分を拡大して表示した．震源の位置は星印．Pg，Sg は地殻中を伝播する P 波と S 波，Pn，Sn はマントル中を伝播する P 波と S 波，Lg は地殻内部にトラップされた S 波，Le は氷床内部にトラップされた S 波，Re は氷床表面を伝播する表面波．氷床が薄いため，トラップされた波が十分成長できず，図 8 に比べて細かな振動が卓越している．大振幅の部分の水平距離は，図 8 のものとほぼ同じである．

グリーンランド氷床モデルを用いた地震波伝播シミュレーション— 19 —

図 10：TOPO-0 による地震波伝播の様子（文献[7]の図を修正）．上から励起後 10，30，50，
70，90 s のスナップショット．赤色は P 波，緑色は S 波を表す．色が濃いほど振幅が大きい．
構造断面のうち，深さ 30 km 以浅の部分を拡大して表示した．震源の位置は星印．Pg，Sg は地殻中を伝播する P 波と S 波，Pn はマントル中を伝播する P 波，Lg は地殻内部にトラップされた S 波，Le は氷床内部にトラップされた S 波，Re は氷床表面を伝播する表面波．地形や氷床基盤の起伏によってトラップ波のエネルギーが漏出するため，図 8，9 と比べて地殻内部に緑色部分が卓越している．

— 20 —

SENAC Vol. 48, No. 4（2015. 10）

図 11：TOPO-5 による地震波伝播の様子（文献[7]の図を修正）．上から励起後 10，30，50，
70，90 s のスナップショット．赤色は P 波，緑色は S 波を表す．色が濃いほど振幅が大きい．
構造断面のうち，深さ 30 km 以浅の部分を拡大して表示した．震源の位置は星印．Pg，Sg は地殻中を伝播する P 波と S 波，Pn はマントル中を伝播する P 波，Lg は地殻内部にトラップされた S 波．氷床に入射するエネルギーが少なく，トラップ機構も弱いため，氷床内の振幅は図8～10 と比べて小さく，氷床がない場合の波動場（図 7）に近い．

グリーンランド氷床モデルを用いた地震波伝播シミュレーション— 21 —

3.2 FLAT0.5-0
氷床が厚さ 0.5 km 一定の水平成層構造と仮定し，震源を氷床直下に置いた場合，前節の結果と似た継続時間の長い S 波後続波が現れた．図 6 の赤線が示すように，S 波後続波の継続時間は
FLAT2.0-0 とほとんど変わらず，氷床の厚さには依存しないと考えられる．一方，周波数は前節の場合よりも高周波が卓越している．図 9 からは，波がトラップされる過程は前節の結果と同様
であるが，氷床が薄いことで波の成長が阻害され，高周波成分が卓越したことがわかる．
3.3 TOPO-0
現実的な地形・氷床厚分布を用い，震源を氷床直下に置いた場合は，水平成層の場合と大きく異なる波形が得られた．S 波後続波の大振幅は Sg の到着からかなり遅れて現れるうえ，継続時
間は短く，波群がコンパクトにまとまっている（図 6 赤線）．このケースも震源が浅いため，地震波のエネルギーの大部分が氷床内に入射する点は前節までと同じである．しかし地形や基盤岩の凹凸によって P 波から S 波，S 波から P 波への変換が起こり，エネルギーが分散されることに加え，波の全反射が崩れてエネルギーは次々と地殻内に漏れ出していく．結果として，氷床内をほぼ水平に伝播する波が選択的に残され，コンパクトな Le 波の波群が現れる．氷床から S 波エ
ネルギーが漏れ出していることは，図 10 のスナップショットで，地殻内に S 波を示す緑色が卓越していることから明瞭に確認できる．
3.4 TOPO-5
現実的な地形・氷床厚分布を用い，震源を氷床の下 5 km に置いた場合は，明瞭な Le 波が現れなかった．このケースでは，震源が深いことで，そもそも氷床内部に入射するエネルギーが少ない．また地形・基盤岩の凹凸によって氷床内へのエネルギーのトラップも弱められている．従って地殻内を伝播する波が卓越し，氷床内を伝播する波は不明瞭となる（図 6 赤線，図 11）．しかし氷床の影響は，氷床がない場合（図 6 黒線）に比べた Sg 波の増幅や後続波の継続時間の増大，といった現象に現れている．これは柔らかい地盤で地震動が増幅される現象と同じメカニズムである．図 6 の太線で示した部分は，群速度 3.1‐2.6 km/s の波群で，図 4 の観測波形の太線部分と対応している．今回の計算例は構造断面が観測のものと異なるため波形の直接比較はできないが，震源が氷床直下にない場合の波形の特徴は理論波形でよく再現できたといえる．

4.

まとめ

本稿では平成 26 年度の共同研究の成果として，局地的な地震波伝播を精度と効率よくモデリングする手法の開発と，グリーンランド氷床の現実的な地形・氷床厚分布を考慮して実行した地
震波伝播シミュレーションの成果を報告した．南極では以前から氷床上での地震観測が続けられているうえ，グリーンランドでも 2011 年から本格的な観測がスタートした．航空機によるメンテナンスや，衛星通信によるデータ転送によって，氷床上で得られた観測地震波形は急速に身近になりつつある．しかし氷床が観測地震波形に与える影響についての本格的な研究はほとんどな
く，試算が行われる場合も水平成層構造が用いられてきた．本研究では，水平成層構造モデルと現実的な地形・氷床厚モデルによるシミュレーションを比較して，構造や震源位置のわずかな違
いが，大きく異なる波形を生み出すことを明らかにした．シミュレーションの結果は次のようにまとめられる：(1) 震源が氷床直下にある場合，氷床の影響が波形に強く表れる；(2) 氷床内部での多重反射による地震波エネルギーのトラップ機構と，地形・基盤岩の起伏によるトラップ崩壊
機構の複合度合いによって，特徴の異なる波動場が形成される．さらに本研究では，氷床内にトラップされた S 波を「Le 波」，氷床表面を伝播する表面波を「Re 波」と命名した．氷床上で得られたデータの普及に伴い，氷床そのものの構造解析や氷床の影響の除去を目的として，観測波

— 22 —

SENAC Vol. 48, No. 4（2015. 10）形に氷床が与える影響の研究は重要な課題となりつつある．本研究は，世界に先駆けてその端緒を開く成果を挙げた．
なお本研究の計算は弾性波のみを扱っているが，実際の氷床や岩石は完全弾性体からずれた性質を持つため，波動場を減衰させる（＝非弾性減衰）．我々は本年，円筒座標系 2.5 次元差分法のプログラムに，非弾性減衰を導入した[29]．今後は氷床の非弾性減衰も考慮し，さらに多くの構造断面について大規模な計算を実行し，観測波形と理論波形の直接比較を行っていきたい．

謝辞
サイバーサイエンスセンターの小林広明センター長，江川隆輔准教授，小松一彦助教，大泉健治氏，山下毅氏には，本共同研究の機会をお与えいただいたうえ，プログラム開発やシステム利用について様々なご助言・ご助力をいただいた．GLISN 日本隊の活動は，科研費（課
題番号 24403006）により運営されている．記して感謝する．
参考文献[1] 豊国源知, グリーンランド氷床上での地震観測点設置, 極地, 94, 30-39, 2012.
[2] 豊国源知, GLISN 日本隊：2012 年の活動紹介, 極地, 96, 16-26, 2013.
[3] Toyokuni, G., Kanao, M., Tono, Y., Himeno, T., Tsuboi, S., Childs, D., Anderson, K.,
Takenaka, H., Monitoring of the Greenland ice sheet using a broadband seismometer
network: the GLISN project, Antarctic Record, 58(1), 1-18, 2014.
[4] 豊国源知・金尾政紀・東野陽子・姫野哲人・坪井誠司, GLISN 計画による日米共同地震観測隊の活動報告(2011～2014 年), 月刊地球, 37(3), 83-91, 2015.
[5] 豊国源知, 地震で調べる北極グリーンランドの環境変動, 科学フォーラム, 370, 28-33, 2015.
[6] Robinson, E.S., Seismic wave propagation on a heterogeneous polar ice sheet, J. Geophys.
Res., 73(2), 739-753, 1968.
[7] Toyokuni, G., Takenaka, H., Kanao, M., Tsuboi, S., Tono, Y., Numerical modeling ofseismic waves for estimating the influence of the Greenland ice sheet on observed
seismograms, Polar Sci., 9(1), 80-93, doi:10.1016/j.polar.2014.12.001, 2015.
[8] 豊国源知・竹中博士・趙大鵬・石原吉明, 月と火星の全球地震波伝播シミュレーション,
SENAC, 46(3), 1-12, 2013.
[9] Toyokuni, G., Takenaka, H., Wang, Y., Kennett, B.L.N., Quasi-spherical approach forseismic wave modeling in a 2-D slice of a global Earth model with lateral heterogeneity,
Geophys. Res. Lett., 32, L09305, doi:10.1029/2004GL022180, 2005.
[10] Toyokuni, G., Takenaka, H., FDM computation of seismic wavefield for an axisymmetricearth with a moment tensor point source, Earth Planets Space, 58(8), e29–e32, 2006.
[11] Toyokuni, G., Takenaka, H., ACE–A FORTRAN subroutine for analytical computation ofeffective grid parameters for finite-difference seismic waveform modeling with standard
Earth models, Comput. Geosci., 35(3), 635–643, doi:10.1016/j.cageo.2008.05.005, 2009.
[12] Toyokuni, G., Takenaka, H., Accurate and efficient modeling of global seismic wave
propagation for an attenuative Earth model including the center, Phys. Earth Planet. Int.,
200-201, 45-55, doi:10.1016/j.pepi.2012.03.010, 2012.
[13] Toyokuni, G., Takenaka, H., Kanao, M., Wiens, D., Nyblade, A., Comparison of global
synthetic seismograms calculated by the spherical 2.5-D finite-difference method withobserved long-period waveforms including data from intra-Antarctic region, Polar Sci., 6,155-164, doi:10.1016/j.polar.2012.06.001, 2012.
[14] Alterman, Z., Karal, F.C., Propagation of elastic waves in layered media by finite
difference methods, Bull. Seism. Soc. Am., 58(1), 367–398, 1968.

グリーンランド氷床モデルを用いた地震波伝播シミュレーション— 23 —[15] Stephen, R.A., A comparison of finite difference and reflectivity seismograms for marine
models, Geophys. J. R. Astr. Soc., 72(1), 39–57, 1983.
[16] Stephen, R.A., A review of finite difference methods for seismo-acoustics problems at theseafloor, Rev. Geophys., 26(3), 445–458, 1988.
[17] Igel, H., Djikpéssé, H., Tarantola, A., Waveform inversion of marine reflection
seismograms for P impedance and Poisson's ratio, Geophys. J. Int., 124(2), 363–371, 1996.
[18] Takenaka, H., Tanaka, H., Okamoto, T., Kennett, B.L.N., Quasi-cylindrical 2.5D wave
modeling for large-scale seismic surveys, Geophys. Res. Lett., 30(21), 2086,
doi:10.1029/2003GL018068, 2003.
[19] Toyokuni, G., Takenaka, H., Kanao, M., Quasi-axisymmetric finite-difference method forrealistic modeling of regional and global seismic wavefield-review and application. In: Kanao,
M. (Ed.), Seismic Waves, Research and Analysis. InTech, ISBN 978-953-307-944-8, 85-112,
2012.
[20] Dziewonski, A.M., Anderson, D.L., Preliminary reference Earth model, Phys.Earth
Planet. Int., 25(4), 297–356, 1981.
[21] Bamber, J.L., Ekholm, S., Krabill, W.B., A new, high-resolution digital elevation model ofGreenland fully validated with airborne laser altimeter data, J. Geophys. Res., 106(B4),6733–6745, 2001.
[22] Bamber, J.L., Layberry, R.L., Gogineni, S.O., A new ice thickness and 400 bed data set
for the Greenland ice sheet, 1. Measurement, data reduction, and errors, J. geophys. Res.,
106(D24), 33773–33780, 2001.
[23] Livingston, C.W., Explosions in ice, Technical Report 75, U.S. Army Snow, Ice, andPermafrost Research Establishment, Engineer Research and Development Center, pp 100,
1960.
[24] Robinson, E.S., Seismic wave propagation on a heterogeneous polar ice sheet, J. Geophys.
Res., 73(2), 739–753, 1968.
[25] Ishizawa, K., The measurement of the velocities of P and S waves propagating in thesurface layer of ice sheet at Mizuho Station, East Antarctica, Antarctic Record, 73, 147–160,
1981.
[26] Takenaka, H., Nakamura, T., Okamoto, T., and Kaneda, Y., A unified approach
implementing land and ocean-bottom topographies in the staggered-grid finite-difference
method for seismic wave modeling, Proc. 460 9th SEGJ Int. Symp., Sapporo,doi:10.1190/SEGJ092009-001.13, 2009.
[27] Nakamura, T., Takenaka, H., Okamoto, T., and Kaneda, Y., FDM simulation ofseismic-wave propagation for an aftershock of the 2009 Suruga Bay earthquake: Effects 440
of ocean-bottom topography and seawater layer, Bull. Seism. Soc. Am., 102(6), 2420–2435,
doi:10.1785/0120110356, 2012.
[28] Nakamura, T., Takenaka, H., Okamoto, T., and Kaneda, Y., Seismic wavefields in the
deep seafloor area from a submarine landslide source, Pure Appl. Geophys., 171, 1153–1167,
doi:10.1007/s00024-013-0717-3, 2014.
[29] Toyokuni, G., Takenaka, H., Zhao, D., Quasi-cylindrical 2.5-D wave modeling with ageneral point source and the anelastic attenuation, Phys. Earth Planet. Int., submitted, 2015.