数理解析研究所講究録 1413 巻 2005 年 36-44

36

火山噴火の MD シミュレーション
東京大学大学院工学系研究科物理工学専攻
湯川論 (Satoshi Yukawa)* 伊藤伸泰 (Nobuyasu Ito)$\text{、}$Department of Applied Physics, School of Engineering, The University of Tokyo概要に、火山噴火に対してマグマだけではなくマグマ内部に溶けている揮発性のガスが重要な役割をブルカノ式噴火を二成分 Lennard-Jones 粒子

担う。 [2] このガスはほぼ水であり、 マグマに対系を用いてモデル化した。 衝撃波管実験をまねして数質量%含まれていることが知られている。

たシミュレーションによりマグマ中に気泡が成長しながら爆発していく様子を再現できた。 また、 ガスとガスの間の衝撃波やマグマとガスの接触面が伝播する様子、 膨張波など全て再現することができた。 さらに、 火山爆発の圧縮性流
体力学モデルと比較を行い、 振る舞いが定性的に一致していることを示した。

ブルカノ式噴火における噴火のダイナミクス
は次のようなものであるといわれている。 ブル

カノ式噴火を示す火山ては、 マグマが高粘性であるため噴火前には火口の上に溶岩ドームが形成されている。 噴火は、 何らかの原因て火道内のマグマーガス系の圧力が上昇し、 この溶岩ドームを破壊した瞬間から起こる。 通常、 マグマーガ

ス系は地殻圧に支えられガス (水) が飽和してい1

初めに日本は有数の火山国であるため、 火山噴火のる状態てマグマ溜まりにある。 この状態て溶岩

ドームの破壊に伴う減圧が発生すると、 飽和溶

解度の限界を超えガスが過飽和になる。 すると

ダイナミクスを理解することは、学術的に価値が ガス成分の核生成が発生し、 マグマ中にガスのあるだけではなく、 防災上も非常に重要である。 気泡が成長する。気泡を含んだマグマーガス系は、
しかし、 そのダイナミクスの物理的側面の理解 密度が軽くなり上昇を始める。 この上昇に伴い、は十分になされているとはいいがたい。 それは、 圧力が低下しさらに気泡が成長する。 気泡の成

そもそもの火山噴火を詳細に観測することや実長が臨界値を越えると、 気泡が浸透しマグマーガ

験を行うことが非常に難しいためてある。 ここでは火山噴火、 特にブルカノ式噴火に焦点を当

ス系はマグマ気泡流から、 マグマの噴霧流へと遷移する。 このような状態の変化が瞬間的に起てシミュレーションを行った結果を報告する。

こるのが、 ブルカノ式噴火である。 時間的な変ブルカノ式噴火とは、 日本ては桜島の噴火に 化が急であること、 状態の変化が劇的であるこ見られるような火山爆発てあり、 比較的高粘性 とが理論的な解析を難しくしている一因てある。
のマグマが関与する。 (これと対比されるのが低このような変化を標語的に表現すれば次のよ
粘性マグマが重要になるストロンボリ式噴火で うになる。 ブルカノ式噴火のイメージは、 よくあり、 ブルカノ式とは噴火様式が異なる。) 一般
振った炭酸飲料のふたをはすしたときの内部の飲料の噴出と重ねることができる。 つまり炭酸*yukawa@ap.t.u-tokyo.ac.jp

37

ガスの溶け込んだ飲料が、 蓋をはすしたときの急減圧で発泡し、 それと同時に飲料が噴き出すことにする。 つまり、 マグマーガスの二成分系をというものである。火山では、 飲料はマグマ、 炭

離散的な粒子をつかつて表現し、 粒子間の相互作用を決め、 あとはニュートンの運動方程式に

酸ガスはマグマに溶け込んでいるガスである。

従って運動するという単純なモデルを採用する。

このような噴火の過程を物理的側面から研究 このようなニュートン方程式に従って運動する
したいのであるが、実際の火山爆発の解析は、火 粒子系において、 さまざまな流体の振る舞いが山固有の特徴や同じ火山でも噴火ごとに違った 再現されることが知られているので、 連続体の様相を見せるなどの個別性から、 直接物理系と 記述の検証にも用いることができる。 [9]
具体的な構成粒子間の相互作用として、 ここ
して理論モデルを構築し取り扱うのは難しい。 このような状況をふまえ、火山学の方で衝撃波管のでは Lennard-Jones 型のボテンシャノレを採用す実験というブルカノ式噴火のアナログ実験が数 る。 Lennard-Jones 型のポテンシャルを用いるこ多く行われるようになってきた。 [5, 3, 4] そこで とにより、 気相、液相、 固相、 またそれらの共存
は、マグマに見立てた粘弾性体や粉体などといっ など、 マグマダイナミクスに必要な相は全て再たアナログ物質を高温高圧でチューブに閉じこ 現することができる。 また、平衡系の性質がよくめ、 隔壁をはすしたあとの爆発的振る舞いを観 知られているので、 パラメーターを選ぶときに測するということが行われている。 これらの状 無駄なシミュレーションを省くことができる。マ況に対応するとみなせる火山噴火の連続体記述 グマとガスの物性の違いは、 Lennard-Jones 粒子による理論モデルもいくつか提唱され、解析もなされているが、 マグマーガスニ成分系の状態方程式を非常に単純化したモデルで決めるなど理論の正当性などはまだまだ明らかでない。 [6, 7, 8]
このような問題点を考えるため、 ミクロな観点からのモデル化を行いシミュレーションを行った。 そこでは、 マクロな仮定をいつさい排除す
るのて連続体モデルの妥当性を議論できる。 まのパラメーターを変更することで表現する。今のモデル化では質量と相互作用をマグマ粒子、ガス粒子で変更することにより物性の違いを表現している。 実際の系のダイナミクスを支配する
3 次元 $N$ 粒子系のハミルトニアンとして次のような物を使用する。
.

$\mathcal{H}=\sum_{i=1}^{N}\frac{\mathrm{p}_{i}^{2}}{2m_{i}}+\frac{1}{2}\sum_{i,j(i\neq j)}^{N}\alpha_{i}\alpha_{j}\phi(|\mathrm{q}_{i}-\mathrm{q}_{j}|)-$$(1)$

た、 衝撃波管の実験と同じような状況をシミュレートすることにより、 非常に理想的でコント$\phi(r)$

ロールされた衝撃波管の実験と見なすこともできる。 そのためこのようなシミュレーションをり、大規模に行うことはブルカノ式噴火の物理的理解につながると考えることができる。

は Lennard-Jones 12-6 ポテンシャノレてあ$\phi(r)=4\epsilon\{(\frac{\sigma}{r})^{12}-(\frac{\sigma}{r})^{6}\}$と一般的に書かれる。 また、 各粒子の質量としてマグマは $m_{magma}=1$ ’ ガスは $m_{gas}=$$m_{i}$$0\cdot 1\cross m_{magma}$

2モデルとシミュレーション

21

モデル

本研究では、 ミクロな観点からモデル化を行(2):を取ることにする。 さらにエネ

ルギーの次元 \epsilon 、 長さの次元$\sigma$

をともに 1 に取ることにし、 これらて全ての次元を無次元化する。 また、 Boltzmann 定数$k_{B}$

を 1 と取ること

で、 温度とエネルギーは同じ単位で測ることにする。 これより (数) 密度の単位が個う。 ミクロなモデルとは言ってもいろいろな立の単位が場があり得るが、 ここでは分子動力学法を使う

どと決まる。

$\sigma\sqrt{m_{magma}/\epsilon}\text{、}$$/\sigma^{3}$

圧力の単位が、時間$\epsilon/\sigma^{3}$な38

相互作用の前の係数$\alpha_{i}$は、 マグマとガスに対に対して運動方程式は、 次のように変更される。

して相互作用を変えるためであり、 マグマに対しては 1 を、 ガスに対しては 0.1 を取ることに$\dot{\mathrm{q}}_{i}=\frac{\partial \mathcal{H}}{\partial \mathrm{p}_{i}}$(3)$\dot{\mathrm{p}}_{i}=-\frac{\partial \mathcal{H}}{\partial \mathrm{q}_{i}}-\zeta \mathrm{p}_{i}$(4)する。 これでマグマ間の相互作用はガス間の相互作用に対して 100 倍強く働くことになる。 もちろん$\mathrm{p}_{i},$$\mathrm{q}i$は、 粒子の三次元での運動量、 座

標を表す。

2.2

$\dot{\zeta}=\frac{1}{\tau}(\sum_{i\in A}\frac{\mathrm{p}_{i}^{2}}{2m_{i}}-\frac{3}{2}N_{A}T_{A})$

ここで、シミュレーション$\sum_{i\in A}$(5)はマグマ溜まりに入っている粒子に対して和を取ることを表しており、$N_{A},$ $T_{A}$はマグマ溜まりの粒子数、 マグマ溜まりの温度を表す。 また、$\tau$

は温度制御の時定数であり正のパラメーターである。 (が熱浴の効果を現しており、制御したい部分の全運動エネルギーが、設定
したい温度で決まる運動エネルギーに収束する
ように、 一種の 「摩擦」で運動量を変化させる。

この運動方程式を用いて平衡状態に十分緩和したあと、 マグマ溜まりと火道を仕切っている弾性壁を取り除く。 これが衝撃波管の実験での隔壁に穴を開けた瞬間に対応する。 これ以降は熱浴の作用を取り除き、 系の時間発展はニュートン方程式に従って時間発展させる。

初期にマグマ溜まりの温度を比較的高温に保ち、 密度も高密度で準備し、 火道内のガスの温

図 1: 系の形

度、 密度を適当に小さく設定しておけば、 マグ

シミュレーションに際しては衝撃波管の実験を念頭におこなう。 系の形状として直方体の箱($L_{x}\cross L_{y}\cross L_{z^{\text{、}}}L_{i}$考える。 マグマはにする。 $x,$$y$

マ溜まりと火道の隔壁を取り除くと、 自発的に爆発が発生する。

はそれぞれの方向の長さ) を$z$

軸正の方向に爆発すること

方向は周期的境界条件を課し、$z$軸方向の境界条件として, 両端に弾性反射壁をおいておく。 (図 1) この弾性壁は、 Lennard-Jones

3

結果

3.1

物理量の時空プロファイル

ポテンシャルの底から見た斥力部分て表現してシミュレーションの典型的な結果を示す。 全いる。初期状態として、系を「マグマ溜まり」部ての物理量がミクロな量で定義されるので、 密分と、「火道」部分の二つに弾性壁で分け、 マグ

マ溜まり部分に、 マグマ粒子とガス粒子、 火道

度場や圧力場、 温度場など連続体との検証で必要な物は全て計算することができる。 それぞれ部分にガス粒子を入れておく。 それぞれで能勢-の物理量は、 系を

Hoover 熱浴を使って等温の熱平衡状態を準備する。 [10, 11, 12] 例えば、マグマ溜まりにある粒子

1 のスライスの中で計算している。例えば、 密度
場 \rho (z)、 圧力場 供などは次のように定義$z$軸方向に輪切りにした厚み$\beta(Z)$39

されている。

(6)$\rho(z)=\frac{\sum_{i\in z}1}{L_{x}L_{y}}$$\Pi_{\alpha \mathcal{B}}(z)=\frac{1}{L_{x}L_{y}}\sum_{i\in z}\frac{(\mathrm{p}_{i})_{\alpha}(\mathrm{p}_{i})_{\beta}}{m_{i}}$(7)$+ \frac{1}{2L_{x}L_{y}}$$\sum_{or,j\in z,(i<j)}f_{\alpha}^{i,\mathrm{j}}q_{\beta}^{i,j}i\in z$$z$

は一つのスライスを指定するインデックスである。 和$\sum_{i\in z}$

はスライス 2 にふくまれる粒子で和を取ることを表している。 圧力場の定義でのいすれかを取り、 (pi)。は粒子$i$の運動量の $a$ 成分を表す。 また、 \mbox{\boldmath $\xi$}.\dashv は、 粒子$i$$\alpha,$と$\beta$$j$

は $x,$$y$ ) $z$の間に働く力の$\alpha$らの粒子の相対座標の和は、 粒子$i$もしくは成分であり、$\beta$$j$$\mathrm{q}_{\beta}^{i,j}$はそれ成分を表す。 第二項のが考えているスライスの中に入っているときに和を取り、 二つの粒子が同時に入っているときは二重に数える。 これは圧力場を運動量の輸送で定義したことになる。
この定義を用いると、 非平衡状態では圧力場は(b) 圧力 (x)(a) 密度等方的にならす、 一般に流れの方向に依存する。
図 2 には密度場と $xy$ 平面内の圧力場を示し

図 2: 局所的な物理量のプロファイル。 系の大きておいた。 (系のパラメーターは図のキャプショ

さは $L_{x}=40,$ $L_{y}=40,$ $L_{z}=752$ 粒子数は、 マンを参照。 ) 横軸がグマ 57600 粒子、 ガス 118400 粒子。横軸は$z$軸のスライスに対応して$z$

軸

おり、 爆発は左から右に起こる。 また、 左端と方向の座標。 縦軸は時刻である。 時刻 0 で爆発

右端には弾性壁が設けられている。 縦軸は時刻を表し、 時刻 0 が爆発の瞬間になっている。 密
度場において、 座標が正の向きに進行している

しその後の時間発展が示されている。 また、 初$z$期のマグマ溜まりの大きさは $40\cross 40\cross 40$ である。 初期状態ではマグマ溜まりにマグマ 57600

衝撃波が二つ見られる。 初めの速度の速い物は、 粒子と、 ガス 6400 粒子をつめ、 温度 2 で熱平衡
ガスの音速を超えておりこれは噴火によるガス 化させている。 また火道部分 40にガとガスの間の衝撃波である。 このような衝撃波 ス 112000 粒子をつめ温度 08 で熱平衡化させて$\mathrm{x}40\cross 704$

は実際の火山噴火の際にも観測されている。 まある。 直方体をた、 それより遅れてもうひとつ高密度領域と低そのスライスの中で局所的な物理量を計算した。

密度領域を分ける波が伝播しているのが観測される。 この波は、 マグマとガスの間の接触面である。 さらに、 爆発後、 時刻 120 当たりから、
軸負の方向に向かって右の弾性壁から衝撃波が$z$走っているのが見られる。 これは、

$z$

軸の上の境

界条件である弾性壁からの反射波なので、 これ$z$軸方向に厚み 1 で輪切りにし、

40以降の計算は上部の境界の影響が出ている。 また時刻 20 ぐらいから、 左の弾性壁からもう一つ
密度波が出ていることがわかる。 この密度波は、膨張波が、 マグマ溜まりの底の弾性壁に当たって跳ね返ってきた分に対応している。 さらに、時刻 40 ぐらいから、 密度場においてマグマーガス

接触面より後方で内部構造が成長して行くのが観測される。 これは詳細なスナップショットの解析から、 マグマ内部のガスが発泡しながら爆発

図 4: 内部構造、 $t=170$ たったあとの内部構造。

構造が成長し整理され、 いくつかの気泡とマグ
マの液滴が見える。 [14] カラー版ではマグマ粒
子が赤、 ガス粒子が青で示されている。

していく様子であることがわかり、 火山学で言われている発泡しながら爆発するという状況が再現できている。 この内部構造の様子は、 次の小節で詳しく見る。 さらに、 密度場ではわからないが、圧力場のプロファイルにおいて、爆発直
後から

$z$

軸負の向きに、 別の波が伝わっている様子が見られる。 これは、 減圧が伝わっている波であり、 膨張波であると見なすことができる。

噴火後、 40 単位時間たったあとの粒子のスナッ

プショットを図 3 に、 170 単位時間たったあとのスナップショットを図 4 に示しておいた。 図で横方向が 軸方向であり、 噴火は右向きに起こ$z$る。初期状態では、 一様に混ざっていたマグマ粒
子 (カラー版では赤) とガス粒子 (カラー版では青)1 が、 噴火後 40 単位時間たったあとでは、 ガ
ス粒子が析出することにより相分離を起こして

32

内部構造

いることがわかる。 また、 このとき生成してい前小節で、 物理量の時空プロファイルにおいて、 噴火しているマグマに内部構造が成長していく様子が見られた。 ここでは、 この内部構造を詳しく見てみよう。

るガスの気泡の大きさがそれほど均一ではなく、様々なサイズに分布していることも見てわかる。

このあともう少し時間がたつと (図 4)、様々なサイズに分布していた気泡が成長、合体し、比較的大きな気泡になっていることがわかる。 また、成長した大きな気泡の内部にマグマの液滴が存在していることもわかる。 さらに、 マグマ内部にもすこし気泡が存在していることも見て取れ
図 3: 内部構造のスナップショット、 $t=40$ たっる。 このように、 時空間のプロットで内部構造が成長しているように見えたところでは、 実際方向であり、 左に粒子分布を見れば、 内部構造が存在し、 且つがマグマ溜まり、 右が火道方向である。 爆発はその構造が成長していく様子も観測される。 こ

たあとの内部構造。 横方向が$z$

右向きに起こる。 粒子は二種類示しており、 マの振る舞いは、 火山学での噴火のシナリオで噴グマと、 マグマ溜まり由来のガス成分を示して霧流遷移が起きると言うことに対応しているが、いる。 火道由来のガス粒子は示していない。 グ

レースケールでははっきりしないが、 様々な大

現在の計算規模では、 実際の遷移にはまだ至っていない。

きさのガス成分の気泡がみられる。 [14] カラー
版では、 マグマ粒子が赤、 ガス粒子が青で示されている。
1 カラー版の原稿をweb て公開している。 [14]

3.3と表現できる。 また、

連続体モデルとの比較$\gamma_{m}$は、 定圧比熱と定積比

熱の比であり、 今の仮定の下では局所的な状態

火山学の方でブルカノ式噴火に対する連続体のモデルがいくつかある。 その中でも衝撃波管を決めると決まる。 方程式の変数をまとめておであり、 これら以
くと、 場の変数は$w,p_{g},$外の量は定数か、 これらから計算できる物であ$T$$\rho,$の連続体の記述をそのまま火山噴火に応用した
Woods (1995) のモデルが解析的にも取り扱い る。 これで変数の数と方程式の数が一致するのやすい。 [6] ここでは、 Woods のモデ) とシミュで、 ダイナミクスが計算できることになる。
1/– ション結果の比較を行う。 Woods のモデル
この支配方程式から見てわかるとおり、 Woods
は、基本的に一次元の一成分圧縮性流体の方程式 のモデルではマグマに溶けているガスが析出しであり、 以下のような支配方程式で記述される。
てくる効果は全くふくまれていない。 このため、(8) マグマからガスが析出し核生成を行う噴火の初期段階の記述は不完全であると思われる。 また、(9) ガスの気泡の体積が大きくなってきて全系に浸透しきった状態への転移後のダイナミクスの記(10)
ただ元々が圧縮性流体の$\mathrm{s}$$\frac{\partial\rho}{\partial t}+w\frac{\partial\rho}{\partial z}=-\rho\frac{\partial\rho}{\partial z}$$\frac{\partial w}{\partial t}+w\frac{\partial w}{\partial z}=-\frac{1}{\rho}\frac{\partial p_{g}}{\partial z}$$\frac{1-n}{\rho\iota}+\frac{nRT}{p_{g}}=\frac{1}{\rho}$(11)$p_{g}( \frac{\phi}{\rho})^{\gamma_{m}}=const$

述も難しいであろう。
方程式であるため、 シミュレーションや実際の観測で見つかつている衝撃波を定性的に記述することが可能であると思われる。 実際、 この方

ここで初めの二つが一般的な質量保存の式と運に取ってい動方程式であり、 時空間の座標は

程式は理想気体の一次元圧縮性流体の記述と構が造が同じ為、 次のような方程式へと変形できる。

$t,$ $z$

る。

$t$

が時間、

$z$

が噴火方向の座標である。

$\rho$

マグマーガスを合わせた質量密度、 が流体の速がガス成分の圧力であり、度場である。 また、
Woods のモデルてはマグマの圧力と等しいとされる。 三つ目の式が、 マグマーガスの状態方程式
であり、 がガス成分の質量分率を表す。 がマグマの質量密度てあり、 $R$ は気体定数、 が温度である。 Woods のモデル化では質量分率 が$w$$\{\frac{\partial}{\partial t}+(w\pm a(\rho))\frac{\partial}{\partial z}\}$$p_{g}$$n$$\rho\iota$

$T$

$n$

も不変
一定であると仮定し、 マグマ質量密度であるとみなす。 このとき、 ガス成分が理想気体$\rho\iota$

的に振る舞うと仮定すれば、 この状態方程式でマグマーガスの状態を記述できる。 最後の 4 番目の式が、 等エントロピー性を仮定して導かれた式であり、 理想気体の等エントロピー流れの解)析で出てくる物と類似している。 (例えばは、 ガス成分の体積分率を表し、 状態方程式

( At
$w$

ここで、 $a(\rho)$$\int^{\rho}\frac{a(\rho’)}{\rho},d\rho’)=0$(13)

はマグマーガス系の音速であり、という関係がある。 これと、
等エントロピーの条件を連立させると、 音速は$a^{2}(\rho)=\gamma_{m}p_{g}/(\rho\phi)$

密度の関数として一意に求めることができる。 この方程式の形から明らかなように、 $w\pm a$ の速度で移動する観測者から見て$w \pm\int^{\rho}\frac{a(\rho’)}{\rho},d\rho’$

が保存量であり、 これらから膨張波の速度などが議論できる。 また、 適当な境界条件を組み合わせることにより、 衝撃波の速度、 膨張波内部の$[13]_{\text{。}}$$\phi$から$\underline{nRT}$

T(z)、爆発方向の平均速度(12)$\phi=\frac{p_{g}}{\frac{nRT}{p_{g}}+\frac{1-n}{\rho_{l}}}=\frac{1}{1+\frac{1-n}{n}\frac{p_{g}}{\rho_{l}RT}}$状態なども議論することができる。
Woods のモデルで計算できる物理量と対応するものを図 5 に示しておいた。 上から順に温度$(\mathrm{v}(z))_{z}$、爆発方向とでは直交する方向の圧力 供 (z)、質量密度ある。 圧力の定義は、 式 (7) と同じてあり、 質$\tilde{\rho}(z)$

$x$

42

量密度は、 それぞれの成分の数密度から計算することができる。 温度、 平均速度の定義は以下のようにした。 温度は、 一般に非平衡状態での定義はよくわからないので、 ここでは速度分散で定義した。 つまり、 平均的な速度$\mathrm{v}(z)$

を

(14)

$\mathrm{v}(z)=\frac{\sum_{i\in z}\mathrm{p}_{i}/m_{i}}{\sum_{i\in}\sim 1}$,と定義し、 これから

(15)

$T(z)= \frac{1}{3}\frac{1}{\sum_{i\in z}1}\sum_{i\in z}m_{i}|\mathrm{v}_{i}-\mathrm{v}(z)|^{2}$て、 温度を定義した。 一般にマグマとガスの温度がそれぞれで違う可能性があるが、 この定義
ではいっさい考慮していない。
またこの図は、 系の大きさが $L_{x}=L_{y}=$$40,$ $L_{\approx}=848$の計算結果であり、 初期マグマ溜

まりの大きさは 40$\mathrm{x}40\cross 80$である。 この中に、数密度 1 てマグマ粒子とガス粒子が入っており、そのうちの 10% がガス粒子である。 平衡状態の変化。 $t=6$ で空間方向にスライスしてある。時温度は先の計算と同じである。 また、 初期の火道部分のガス粒子の密度が、 先のと比較すると空間のプロファイルから、衝撃波面、 マグマと火

1/5 になってるが、 これは主として計算時間を稼

道内ガスの接触面が同定され、 初期条件からどぐためであり、 定性的な振る舞いは全く変わら

図 5: 温度、速度 (z)、圧力 (x)、 質量密度の空間のような波も通過していない領域が同定される。 ない。
図 5 には膨張波の波頭がマグマ溜まりの底にこれらを図に示しておいた。 グレーの長方形で囲まれている領域が、 左から順に、 初期状態の到達する前の時刻 $t=6$ での物理量のプロファイ

まま、 マグマとガスの接触面、 初期状態のまま、 ルを表示した。 横軸は、 噴火方向の座標 であの領域を表している。 それらに挟まれている領 り、 $z=80$ までが初期のマグマ溜まりの領域で$z$

域が、 左から膨張波領域、 ガスの衝撃波領域である。 火道は $z=848$ まで続いているが、 この

図では平衡状態から何も変わっていない大部分
違い $L_{x}=40,$ $L_{y}=40,$ $L_{z}=848$ である。 また、 のところはのぞいて $z=200$ まてプロットして初期のマグマ溜まりの大きさは、 $40\cross 40\cross 80$ て ある。 初期に準備した平衡状態での温度が、 マある。 それ以外のパラメーターはほぼ先の図の グマ溜まり T=2、 火道 $T=0.8$ であることと、計算と同じであり、 長くなった分粒子数が増え 図の温度プロファイルから、 膨張波の波頭が到ある。計算している系の大きさは, 先ほどのとはている。

達していないところと、 衝撃波の波面が到達していないところがはっきりとわかる。 これらのとし
Equilibriumて表示してある。 また、 時空間のプロットから、
位置 100 程度の所に、 マグマとガスの接触面が領域を図では、$\lceil \mathrm{I}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}$$\mathrm{S}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{e}\rfloor$43

あることがわかっている。 これは、質量密度のプにおいてそうはなっていない。 これは、 衝撃波

ロファイルからもはっきりと読み取れる。 このマ 面の通過による加熱のミクロなダイナミクスをグマーガス接触面と、 膨張波の波頭までの間が膨 反映しているのではないかと考えられ、 もう少張波が存在しているところであり、 温度や速度、 しマグマ溜まり部分が長い系の計算が必要であ
圧力、 密度になめらかな変化が見られる。 さら る。 また定量的な衝撃波面の速度や、 膨張波のに、 マグマーガス接触面と前方の初期状態のまま

波頭速度などの解析はまだ行ってはおらす、 今でいる領域の間は、 ガスの衝撃波が通り抜けた

後の課題である。

あとの高温に加熱されている部分であり、 温度
や、 速度は主としてガス粒子が担っていることが読み取れる。 速度のプロファイルで、 マグマー

ガス接触面直前に比較的高速になっている部分があるが、 これが実際にガス粒子が加速されている領域で、 熱化する前の状態なのか、 たんな

る揺らぎなのかはサンプル平均を取っていない
のでまだ何ともいえない。

4

まとめ
ブルカノ式噴火のミクロな観点からのモデルを構成し、 シミュレーションを行った。 この結果、 ガス間、 ガスとマグマ間の衝撃波の存在や、減圧が伝わる波である膨張波などの波が観測さこの計算結果を定性的に連続体モデルの結果 れた。 これらは、 火山学の方でも存在が認めらと比較してみよう。 連続体モデルは理想気体の れており、 マグマーガス二成分系の相の状態を知圧縮流体モデルであるから、 既にわかつている る上での指標となるので重要である。 また、 定
衝撃波管の解析を流用することができる。 その 性的にではあるが、 噴火直後からガスが核形成
解析によれば上のような図を書いたとき、 プロ を経て、 発泡し複雑な内部構造を生成しているファイルは 5 つの領域に分かれることが知られ ことが確認された。 これらは、 気泡を伴うマグている。 衝撃波が伝わる下流側から見れば、 (1) マ流からマグマの液滴を含んだガス流への転移
初期の低温平衡状態が保たれている領域、 (2) 衝 の際に重要になる発泡の浸透につながる非常に撃波面が通ったあとの低温状態にいたガスが加 重要な現象である。
熱された高温領域 (熱気体の領域) (3) 高温状態
シミュレーションで得られた結果と、 連続体にいたガスが断熱膨張で冷却された領域 (冷気体モデルとの比較も定性的におこない良い 致をの領域) (4) 膨張波が存在する領域、 (5) 初期の見せている。 しかし、 定量的な比較や、 衝撃波$\text{、}$$-\wedge$$\backslash$

高温平衡状態にある領域、 である。 シミュレー
ションの結果をこの領域分割と比較すると、 領

域 (3) をのぞく全ての領域が完全に再現されていることがわかり定性的に良い一致を見せている。 またはっきりとは再現されていない領域 (3)も、 マグマーガス接触面であると見なした部分をこの領域であると見なせないこともなく、 領域(3) が無いとは言い切れない。
定性的には非常に良い一致を見せてはいるが、問題ももちろんある。 理想気体の場合は一般に(4) の領域をのぞくと、 残り四つの領域で物理量は一定てあるが、 シミュレーションでは領域 (2)管の実験結果などとの比較は、 まだ行ってはおらす今後の課題である。 ただ定性的に良い一致を見せているところから定量的にもこのミクロなモデルは良いモデルであろうと思われる。今後、 マグマ流からガス流への完全な転移をしら

ベマグマの破砕の様子を定量的に調べるためには、 もう少し計算規模を大きくする必要がある。
これは、 現在のシミュレーションで系の断面の1/4 の大きさの液滴が見えていることから、横方向の断面積を数倍大きくすれば見えてくるはずの現象であり、 計算機中での火山噴火の完全な再現は、 もう手の届くところにある。

44

参考文献
[1] 特定領域研究 「火山爆発のダイナミクス」研究
集会、年、 火山爆発夏の学校」 テキスト。
$\lceil 2004$[2] 鍵山恒臣編、 東京大学地震研究所編集、 地球科
学の新展開 3 「マグマダイナミクスと火山噴火」
朝倉書房、 2003 年。

[3] M. Ichihara, D. Rittel, and B. Sturtevant: “Fragmentation of a porous viscoelastic
material: Implications to magma fragmentation”, J. Geophys.
. 107(Bl0), 2229,
doi:10.1029/2001JB000591, (2002).
${\rm Res}$[4] O. Spieler, D. B. Dingwell, and M. Alidibirov:“Magma fragmentation speed: an experimendetermination”, J. Volcanol. Geothem.
. $129,109- 123,$ $(2004)$ .
$\mathrm{t}\mathrm{a}\mathrm{l}$${\rm Res}$J.
[5] B. Cagnoli, A. Barmin, O. Melnik, R.
Sparks: “Depressurization of fine powders ina shock tube and dynamics of ffagmented
magma in volcanic conduits”, Eanh Planet.
Sci. Lett. $204,101- 113,$ $(2002)$ .
$\mathrm{s}.$[6] A.
Woods: “A model of vulcanian explosions”, Nucl. $Eng.$ Design, 155, 345-357,(1995).
$\mathrm{w}.$[7] 0. Melnik: “Dynamics of two phase conduit
flow of high-viscosity gas-saturated magma:large variations of sustained explosive eruption
intensity”, Bull. Volcanol. 62, 153-170, (2000).
[8] O. Melnik and R. S. J. Sparks: “Nonlinear dynamics of lava dome extrusion”, Nature $402$ ,37 41, (1999).
$\text{ー}$[9] T. Ishiwata, T. Murakami, S. Yukawa, andN.
: “Particle Dynamics Simulations of theNavier-Stokes Flow with Hard Disks”. Int. .
Mod. Phys. (2004) 印刷中。
$\mathrm{I}\mathrm{t}\mathrm{o}$$J$$\mathrm{C}$[10] S. Nos\’e: “A molecular-dynamics method for
simulations in the canonical ensemble”, $Mol$ .
Phys. 52, 255, (1984).
[11] S. Nos\’e: “A unified formulation of the constant temperature molecular-dynamics meth”, J. Chern. Phys. 81, 511, (1984).
. G. Hoover: “Canonical dynamics: Equi[12]distributions”, Phys. Rev.
librium$\mathrm{o}\mathrm{d}\mathrm{s}$$\mathrm{w}$$\mathrm{p}\mathrm{h}\mathrm{a}\mathrm{s}\triangleright \mathrm{s}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{e}$$\mathrm{A}31,1695,$$(1985)$.

[13] 松尾一泰「圧縮性流体力学」理工学社、 1994 年。
[14] この原稿のカラー版を web ページに期間限定で掲示しておく。 カラーて図を見たい人
や、 動画を見たい人はアクセスしてほしい。

http:$//\mathrm{b}\mathrm{o}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}.\mathrm{t}.$u-tokyo. ac.j$\mathrm{p}/\text{ー}\mathrm{y}\mathrm{u}\mathrm{k}/\mathrm{v}\mathrm{o}\mathrm{l}\mathrm{c}\mathrm{a}\mathrm{n}\mathrm{o}/$