数理解析研究所講究録
第 1800 巻 2012 年 195-202

195

火山噴煙挙動の数値的研究

東京大学地震研究所 鈴木 雄治郎 (Yujiro Suzuki)
Earthquake Research Institute,
The University of Tokyo1.

はじめに火山噴火は，溶岩流出から噴煙放出まで多岐に渡る地中のマグマが地上に噴出する現象である．一般に，地中で高温高圧状態にある液体マグマ (溶融した岩石) は，主に水からなる揮発性成分を含有している．マグマが上昇して減圧すると，揮発性成分は液体マグマ中に気泡として析出する．更に減圧が進むと，ガスの体積分率が 80% を超え，ガスが主体の高速流れとなる．この時，細かく分断したマグマは，火山灰や軽石として高速のガスに加速され上昇する．火山
灰と火山ガスからなるこの混合物が火口に達すると，混合物は爆発的に大気中へと放出される (Fig. 1). このような噴火様式を爆発的噴火と呼ぶ．一方で，揮発性成分が液体マグマから効率的に抜けると，火口からは液体マグマが溶岩として流出することとなる (非爆発的噴火). 火山噴火時に爆発的噴火となるか非爆発的噴火となるかは，火山学上最重要課題の一つとして長年研究され，マグマの物性や火道 (マグマ溜まりから火口への通路) の径，ガス抜けの浸透現象，マグマの破砕条件などが複合的に絡み合って決まっていることが分かってきている [1, 2].

Fig. 1 Eruption cloud in Sakura-jima volcano, Feb. 82010.

196
爆発的噴火では，噴出物 (火山灰 火山ガス) が大気と混合し，噴煙として大気中に拡大する．噴煙は固体粒子である火山灰を含むため，噴出直後は周囲の大気よりも重い状態にある．しかし，取込んだ大気を火山灰の熱によって急$+$激に膨張させると，周囲の大気よりも軽くなって浮力で上昇する．噴煙柱と呼ばれるこの上昇流は高度数百 m$\sim$

数 km, 巨大噴火では $30km$ まで達することがある．一方，噴出物が浮力を得るために十分な大気を取り込まないと，噴煙は周囲の大気よりも重いままで山体を流れ下る．この下降流は火砕流と呼ばれ，火山帯周辺を荒廃させることとなる．火砕流はその上面から周囲の大気を取込み，軽くなって「灰かぐら」と呼ばれる上昇流を形成する．噴煙柱や灰かぐらといった上昇流は，その密度が成層大気の密度と釣り合うと，それ以上は上昇することはできず，傘型噴煙と呼ばれる水平方向に拡大する噴煙形態となる．
一般に，噴出率 (単位時間当たりの噴出密度体積) が大きくなるほど火砕流が発生しやすく，また，より噴煙高度や傘型噴煙半径の拡大速度が大きくなることが分かっている [3]しかし，噴出率とそれらの定量的な関係性は未だ明らかになっておらず，防災上も重要な問題である．

本稿では，火山噴煙の数値モデル化とこれまでに得られた結果を紹介し，その発展性について議論する．

2. 火山噴煙の数値モデル化
噴煙ダイナミクスは，成層大気中の固気混相の高速流れとして捉えられる．
その支配要因として以下が挙げられる :a) 大気との乱流混合，b) 密度の非線形変化，c) 噴煙内外の衝撃波形成，d) 大気の成層構造，e) 大気の風，f) 火山
灰の分離，g) 水の相変化，h) 火山灰の凝集，i) 火山の地形，j) 噴火の非定常性．筆者はまず，なるべく単純な系で噴煙挙動を再現することを目的に，a-d の効果を入れた噴煙の 3 次元モデルを構築した [4].
火山灰が噴煙から分離する効果 ( を無視し，噴出物をひとつの流体として捉$\mathfrak{d}$えると，噴煙挙動は質量・運動量．エネルギー保存式からなる Euler 方程式と，色付け関数として噴出物のみの質量保存式で与えられる．
$\frac{\partial}{\partial t}\{\begin{array}{l}\rho\rho u\rho e\rho\epsilon\end{array}\}+\nabla\cdot\{\begin{array}{l}\rho u\rho uu+pI(\rho e+p)u\rho\epsilon u\end{array}\}=\{\begin{array}{l}0-\rho g-\rho g\cdot u0\end{array}\}$ここで，は密度，は速度，は全エネルギー，は噴出物の質量分率，は圧力，は単位行列， は重力加速度である．また，火山灰，火山ガス，大気の混合物$u$$\rho$$I$$e$$\epsilon$$g$

である噴煙の状態方程式は以下で表される．
$\frac{1}{\rho}=\frac{n_{s}}{\rho_{s}}+\frac{(n_{g}R_{g}+n_{a}R_{a})T}{p}$

$p$

197
$n_{s},n_{g},$ $n_{a}$はそれぞれ火山灰，火山ガス，大気の質量分率， は火山灰密度，$R_{g},$ $R_{a}$$\rho_{s}$はそれぞれ火山ガス，大気の気体定数である．右辺第一項は火山灰の占める体積を表しており，大気圧中では噴煙全体の体積から比べて無視できるぐらい小さい ( $<1$ vol.%). したがって，噴煙の状態方程式は，$p\approx\rho(n_{g}R_{g}+n_{a}R_{a})T$と，理想気体の状態方程式で近似できる．ただし，各地点における噴煙全体の気体定数 (ngRg$+$naRa) を，混合比によって随時変化させる必要がある．

理想気体近似を用いた噴煙の支配方程式は，一般的な圧縮性流体の計算手法で解析することができる．数値計算では，流束計算に Roe 法 [5] を用 , MUSCL法 [6] を適用することで 3 次精度化している．時間積分に関しては，時間分割法を用いている．数値計算は，海洋研究開発機構の地球シミュレータで 64-256$\mathfrak{h}\backslash$

CPU を用いて実行した．
3. 数値計算結果
数値シミュレーションの結果，噴煙の特徴的な挙動を再現することができた．
中規模噴火に属する噴出率が $10^{7}kg/s$ の計算では，噴煙柱と傘型噴煙が再現された (Fig. 2). 平坦な地表にある円形の火口からの噴出物は，高度約 $15km$ までの$15km$ ま噴煙柱を形成した．噴煙の密度は高度約 $12km$ で大気密度と釣り合い，で達した噴煙は $12km$ まで流れ落ち，水平方向へと拡大して傘型噴煙を形成した．一方，大規模噴火に属する噴出率が $10^{9}kgls$ の計算では，火砕流の発生と，火砕流からの上昇流 (灰かぐら) が再現された (Fig. 3). これらの結果は，噴出率が大きいほど火砕流が発生しやすいという観測事実とも整合的である．

$v\iota m**\S\infty$meX 評鴫$\oint$Fig. 2 Numerical simulation result of eruption column and umbrella cloud.

198

Fig. 3 Numerical simulation result of pyroclastic flow and $co$ -ignimbrite ash clouds.

$x$

[km]

$x$

[km]

Fig. 4 Flow patterns of eruption clouds near the volcanic vent. Upper left: Eruption

column with ajet stmcture. Upper right: Eruption column with a fountain. Lower
left: Column collapse with ajet structure. Lower right: Column collapse with afountain.

199
火口付近の流れに注目し，噴出温度，噴出率，噴出速度のパラメータスタデ

イを行なったところ，4 っの流れのレジームが存在することが分かった (Fig. 4).
噴出率が小さい時には，中心軸の高濃度領域が上昇とともに消滅するジェットのような構造を持つ噴煙柱となった (Jet タイプ噴煙柱 ;Fig. 4 左上). 噴出率を大きくすると，高濃度領域が消滅する前に運動量を失うため，噴水のような構造を持つ噴煙柱となった (Fountain タイプ噴煙柱 ;Fig. 4 右上). 更に噴出率を大きくすると，噴水構造は十分な大気と混合することができず，噴煙柱は崩壊し火砕流を発生させた (Fountain タイプ崩壊 ;Fig.4 右下).

一方で，噴出率が小さくて噴出温度が低い場合には，ジェット構造を持ち，周囲の大気とよく混合するが，取り込んだ大気を十分に膨張させるだけの熱量がないために噴煙柱は崩壊する (Jet タイプ崩壊 ;Fig. 4 左下).

以上 4 つのレジームについて詳細なレジームマップを作成した (Fig. 5). 噴出率が小さく，噴出速度が大きい場合には，Jet タイプ噴煙柱となった．逆に，噴出率が大きく，噴出速度が小さい場合には，Fountain タイプ崩壊レジームとなった．それらの間の遷移レジームは噴出温度によって変化する．噴出温度が低い場合には，Jet タイプ崩壊が遷移レジームとなった (Flg. 5 左). 一方，噴出温度が高い場合には，Fountain タイプ噴煙柱が遷移レジームとなった (Fig. 5 右).
レジーム境界のパラメータ依存性を明らかにするために，リチャードソン数とマッハ数の無次元空間でのレジームマップを作成した (Fig. 6). その結果，火砕流の発生条件は温度によってマッハ数依存性が変化することが分かった．噴出温度が低い場合，火砕流発生条件は Fountain 形成条件よりも低リチャードソン数側にあり，その境界値はマッハ数に依存しない．一方，噴出温度が高い場
合，火砕流発生条件は Fountain 形成条件よりも高リチャードソン数側にあり，マッハ数が増加するほどその境界値は減少する．更なる解析を行なったところ，噴出率$\zeta k\mathfrak{g}/s]$

噴出率 [夏$glsJ$

噴出率 $[|cg/s]$

Fig. 5 Flow regime maps of eruption clouds in dimensional spaces. Left: lowtemperature case. Center: intermediate temperature case. Right: high temperaturecase.

200Richardson number$R$に hardson number$Rid|ard\infty n\mathfrak{n}umber$

Fig. 6 Flow regime maps of eruption clouds in non dimensional spaces (Richardsonand Mach numbers). Left: low temperature case. Center: intermediate temperature$-$case. Right: high temperature case.
このマッハ数依存性は，Fountain 内に衝撃波構造が形成することに起因することが分かった [8]. Fountain が形成する場合，バレルショックやマッハディスクと呼ばれる定常衝撃波が形成しやすく，それらが噴煙と大気の混合を阻害する．
結果として，噴煙は十分な浮力を得ることができずに火砕流が発生しやすくなる．

4. モデル発展の方向性

これまでに，固気混相流である噴煙の流れを理想気体近似することで，噴煙柱や火砕流といった現象を再現し，火砕流の発生条件を定量的に求めることができた．火山灰がどこまで拡大し，降灰するかを理解するためには，2 章で挙げた噴煙ダイナミクスの支配要因 $e-j$ についても考慮したモデルを構築しなければならない．
風の効果 (e) については，既にモデルを構築し，風のプロファイルについてのパラメータスタディに入りつつある (Fig. 7). 風の存在によって噴煙はたなびくだけではなく，大気との混合率も変化することで，最終的に到達する高度にも影響する可能性がこれまでの結果が示唆している．この風の存在と混合率の関係は，火山学上のみならず，流体力学上も非常に重要な問題として解析を進める必要がある．

火山灰粒子の分離 (0 については，これまでの理想気体近似モデルにマーカー粒子を入れるモデルを構築し，予備計算を行なっている (Fig. 8). 実際の噴煙では，火山灰が分離することによって，噴煙の質量と熱源が失われることとなり，噴煙のダイナミクスは大きく変化することが予想される．この火山灰から噴煙への影響を含めたモデルを構築するためには，無数に存在する火山灰粒子をどのように数値的に扱うかが問題となる．ひとつの方向性としては，計算機で扱

201

X 家1$[$Fig. 7 Numerical simulation result of eruption cloud in a cross-wind field.

Fig. 8 Preliminary result of eruption cloud with marker particles.

うことが出来る範囲の粒子数で全火山灰を代表し，衝突・合体などを実際の現象を正しく表すように確率的に与える超水滴法 [7] などを用いることが考えられる．超水滴法は降雨シミュレーションで成功を収めており，火山噴煙の問題に適用することはチャレンジングであると考えられる．

202
5. 参考文献[1] T. Kozono, T. Koyaguchi, “Effects of relative motion between gas andliquid on 1-dimensional steady flow in silicic volcanic conduits: 1. An$180,21-36$ , 2009.
analytical method”, J. Volcanol. Geotherm.
[2] T. Kozono, T. Koyaguchi, “Effects of relative motion between gas andliquid on 1-dimensional steady flow in silicic volcanic conduits: 2. Origin of$180,37-49$ , 2009.
diversity of eruption styles, J. Volcanol. Geotherm.
[3] A. W Woods, “The fluid dynamics and thermodynamics of eruption
columns”, Roll. Volcanol., 50, 169-193, 1988.
[4] Y. J. Suzuki, T. Koyaguchi, M. Ogawa, I. Hachisu, ‘A numerical study ofturbulent mixing in eruption clouds using a three dimensional fluid, 110, B08201, 2005.
dynamics model”, J. Geophys.
[5] P. L. Roe, “Approximate Riemann solvers, parameter vectors, anddifference schemes”, J. Comput. Phys., 43, 357-372, 1981.
[6] B. van Leer, “Towards the ultimate conservative difference scheme III.
Upstream-centered finite-difference schemes for ideal compressible flow“, J.
${\rm Res}.,$${\rm Res}.,$$-$${\rm Res}.$Comput. Phys., 23, 263-275, 1977.
[71 S. Shima, K. Kusano, A. Kawano, T. Sugiyama, S. Kawahara, “Thesuper-droplet method for the numerical simulation of clouds andprecipitation: a particle-based microphysics model coupled with
non-hydrostatic model“, Q. J. R. Meteorol. Soc., 135, 1307-1320, 2009.
[8] Y. J. Suzuki, T. Koyaguchi, 3-D numerical simulations of eruption
column collapse: Effects of vent size on pressure-balanced $jet/plumes$ ”, J.
Volcanol. Geotherm. Res., accepted.