地震波の波線追跡でシャドウゾーンを再現する

大きな地震が起きると、その揺れは地球の裏側にある地震計にも届きます。ところが観測記録を集めてみると、震源からある角度を超えた先ではS波だけがぴたりと途絶えます。届くはずの波が届かない。この欠落こそが、誰も直接見たことのない外核が液体であることの根拠でした。教育用3Dアプリ「地球地下シミュレータ」では、この推論をブラウザ上で追体験できるように、スネルの法則にもとづく波線追跡を実装しました。この記事では、その計算の中身を紹介します。

シャドウゾーンという観測事実

地震波には縦波のP波と横波のS波があり、S波はずれ(せん断)の変形を伝える波なので、ずれに抵抗しない液体の中を伝わりません。20世紀初頭、世界各地に地震計が置かれるようになると、震源からの角度(震央距離)が約103度を超える観測点にS波が届かないことが分かってきました。S波の届かないこの領域がシャドウゾーンです。地球の深部に液体の層があってS波を遮っていると考えれば、この影は説明がつきます。1906年に地震波の解析から核の存在が指摘されてから、1936年に内核が見つかるまで、地球内部の描像はこうした「届く波と届かない波」の分布から組み上げられました。

シミュレータで再現したかったのは、この推論の過程です。断面表示した3D地球に震源を置くと波線が広がり、観測点に到達時刻が記録され、S波の届かない領域が浮かび上がる。目指したのはそういう画面です。

速度勾配で曲がる波線

素朴に考えると、波は震源から直線で広がりそうです。均質な媒質ならそのとおりですが、地球内部の地震波速度はおおむね深さとともに増すため、波線は曲がります。波面のうち深い側が速く進む分だけ、進行方向が浅い側へ傾き続けるからです。

球対称な媒質では、波線の進行方向(単位ベクトルdir)は経路長sに沿って次の式に従います。

d(dir)/ds = -(v'(r)/v) · (r̂ - (r̂·dir)·dir)

vは地震波速度、v'(r)はその半径方向の勾配、r̂は位置の単位ベクトルです。動径と角度の極座標で積分する定式化では、波線が最深点(転回点)で折り返す瞬間に特異性が現れますが、デカルト座標のままこの式を積分すれば転回点の特別扱いが要りません。実装では1ステップ15kmの中点法で積分し、ステップの中間点で曲げを評価し直して精度を確保しています。

// 中点の位置・方向で1ステップ分の曲げを評価する
const [hdx, hdy] = bend(x, y, dx, dy, ds / 2);
const midX = x + hdx * (ds / 2);
const midY = y + hdy * (ds / 2);
[dx, dy] = bend(midX, midY, dx, dy, ds);

計算が正しいかどうかは、保存量で確かめられます。球対称媒質の波線ではレイパラメータ p = r·sinθ/v(θは動径と進行方向のなす角)が経路全体で一定に保たれます。Vitestのテスト(全16ケース)では、斜め下40度に射出したレイのpを経路上の複数点で計算し、3%以内で一致することを検証しました。

層境界の屈折とS波の消滅

速度勾配による屈折だけでは、まだ足りません。地球内部には速度が不連続に跳ぶ層境界があり、本モデルは地殻から内核までを6層で表しています。境界に達した波線には、ベクトル形式のスネルの法則を適用します。

const cos1 = -(nxv * dx + nyv * dy);
const eta = v2 / v1;
const k = 1 - eta * eta * (1 - cos1 * cos1);
if (k < 0) {
  // 全反射
  dx += 2 * cos1 * nxv;
  dy += 2 * cos1 * nyv;
} else {
  const coeff = eta * cos1 - Math.sqrt(k);
  dx = eta * dx + coeff * nxv;
  dy = eta * dy + coeff * nyv;
}

kが負になるのは屈折方向の解が存在しない場合、つまり全反射で、このときは反射して追跡を続けます。そしてS波には、もう一つの分岐があります。外核は液体でS波速度が0のため、核マントル境界(深さ約2891km)で屈折先の速度が0と判定されたレイは、そこで消滅します。

if (v2 <= 0.01) return { points, exit: 'absorbed' }; // S波が液体外核へ

浅い角度で射出されたS波はマントル内で曲がりながら地表へ戻りますが、深い角度のS波は核マントル境界に突き当たって消えます。この分かれ目が、地表ではシャドウゾーンの境界として現れます。

シャドウゾーン境界の解像度

波線は震源から全周へ一様に180本放射しています。ただ、一様サンプリングだけでは、シャドウゾーンの境界がぼやけます。境界付近では、核マントル境界をかすめるかどうかで波線の行き先が急変する(射出角のわずかな差が震央距離の大きな差になる)ため、「境界ぎりぎりまで届くレイ」を取りこぼすからです。そこで、核マントル境界にちょうど接する臨界射出角を先に求め、その近傍に追加のレイを配置しました。

// 臨界のわずかに上(かすめて遠くへ届く)とわずかに下(核へ入る)の両側
const factors = [0.97, 0.99, 0.997, 1.001, 1.003, 1.01, 1.03, 1.1];

臨界射出角の計算にも、さきほどのレイパラメータが使えます。核マントル境界に接する波線のpは境界半径を境界直上の速度で割った値になり、pは保存されるので、震源での射出角はsin i = p·v/rで逆算できます。臨界の0.997倍や1.003倍といった際どい角度のレイを足すことで、シャドウゾーンの境界が滑らかに定まります。

震源設定時の事前計算

波線追跡はアニメーションの再生中には走らせず、震源を設定した瞬間にすべて計算します。P波とS波それぞれで全レイを追跡し、5秒刻みの時刻グリッドへ線形補間で再サンプリングして、Float32Arrayのテーブルに書き込みます。再生時は各フレームでテーブルの該当時刻を引くだけです。マントル対流のように終わりなく時間発展する計算ならWeb Workerに逃がすところですが、波線追跡は有限の計算で完結するので、その必要がありませんでした。事前計算が2波種合わせて500ms未満で終わることも、テストで確認しています。

各レイが地表に達した位置と時刻は走時データとして記録します。断面円周上の観測点ごとにP波とS波の初動到達時刻を求め、走時表と、震央距離を横軸に到達時刻を縦軸にとった走時曲線で表示します。震央距離150度や175度の観測点にはS波が届かず、走時表のS波の欄は空白のまま残ります。20世紀初頭の地震学者が世界中の観測記録を並べて目にしたのと、同じ形の欠落です。シャドウゾーンの範囲も走時データから求めます。S波が届いた最大の震央距離から先をS波シャドウとして断面上にハイライトし、P波についても直接波と核通過波の間に生じる隙間(P波シャドウ帯)を検出して表示します。

実地球とのずれ

出来上がったモデルを動かすと、S波シャドウは約93度から始まります。実地球では約103度です。10度足りません。原因は速度モデルの簡略化にあります。本モデルは標準地球モデルPREMの層境界値を層内で線形補間していますが、実際のPREMは深度の多項式で表され、最下部マントルでは速度勾配が平坦になります。線形補間ではこの平坦化が再現されず、最下部マントルの速度勾配が実際より急になるため、波線が早めに曲がり戻り、実際ほど遠くまで届きません。このずれは、波種変換や表面反射を扱わないことと合わせて、アプリのUI上に注記しています。届かない波から見えない内部を推論するという体験の骨格は、10度ずれても変わらないと判断しました。

人類は地球の中心まで掘って確かめたことは一度もありません。それでも6371km先の核まで層構造を描けているのは、走時曲線の欠落という紙の上のデータを手がかりに、推論を重ねたからでした。波線を1本ずつ追う実装は、その推論を手元の画面でなぞり直すための道具です。なお、同じシミュレータのマントル対流は、終わりのない時間発展計算のためWeb Workerに載せています。その実装は「マントル対流をWeb Workerでリアルタイムに計算する」で紹介しています。