火山 第 58 巻 ( 2013)
第 4 号 551-555 頁

寄 書

個別要素法を用いた中-高粘性マグマの噴出シミュレーション
酒井望＊・楠本成寿＊＊・清水賀之＊＊＊（2012 年 3 月 8 日受付，2013 年 9 月 2 日受理）

Numerical Simulations of Middle-high Viscous Magma Extrusion
by Means of Discrete Element Modeling
Nozomi SAKAI＊, Shigekazu KUSUMOTO＊＊ and Yoshiyuki SHIMIZU＊＊＊
We developed a discrete element model to describe extrusion of middle-high viscous lava at the surface. We set the710viscosity of magma at 10 to 10 Pa s for the numerical simulations and arbitrarily set the rising velocity to be constantto the bottom of the conduit, as an extrusion model of magma. The conditions under which the lava dome or lava flow
will form depend on not only the viscosity of the magma but also its rising velocity in the conduit. The velocity at which
the lava dome is formed is dependent on the viscosity of magma, with magma of lower viscosity requiring higher risingvelocity in order to form a dome. The inside structure of the lava dome is a concentric circle and is similar to thoseformed in previous analogue experiments.
Key words : discrete element modeling, lava flow, lava dome, PFC2D

1．は じ め に成がドームの成長に重要な役割を果たすことを，Hale

マグマの緩やかな噴出は，地表に溶岩ドームや溶岩流

and Wadge (2008) は，ドーム底部と火道の分岐点に形成を形成する．
溶岩ドームを形成するマグマは，玄武岩質から流紋岩質まで幅広い（例えば，Siebert et al., 2010）．溶岩ドームされるシアーバンドの発達がドームの内成的成長から外成的成長への転換をもたらすことを，それぞれ明らかにしている．の形成・発達には，内成的と外成的の 2 つの様式がある．

一方，溶岩流についても野外調査やアナログ実験，数内成的成長は，新鮮な溶岩が表面に出ることなくドーム

値 シ ミ ュ レ ー シ ョ ン が 行 わ れ て き て い る（例 え ば，を膨張させる成長様式のことであり，外成的成長は，上Walker et al., 1973; Blake and Bruno, 2000; Filippucci et al.,

昇してきたマグマがドームの殻を突き破り，新鮮な溶岩2010）．Dragoni et al. (2005) は，ニュートン流体よりもがドーム表面に付け加わることで成長する様式のことでビンガム流体の方が溶岩流先端部の勾配が急になることある（例えば，中田，1996）．溶岩ドームの形成や成長過や，粘性が高いほど溶岩流が厚くなることを明らかにし程を解明するために，これまで多くの野外調査やアナロている．

グ実験，数値シミュレーションが行われてきた（例えば，

このように，これまで，それぞれの構造や形成プロセ

Blake, 1990; Nakada et al., 1995; Balmforth et al., 2000;

スについて多くの研究成果が挙げられてきているが，溶

Zavada et al., 2009）．Fink and Griffiths (1990) は，殻の形

岩流と溶岩ドームの形状は連続的なものであり，このよ＊〒930-8555 富山市五福 3190
富山大学理学部地球科学科
Department of Earth Science, Faculty of Science, University of Toyama, 3190 Gofuku, Toyama-shi, Toyama, 9308555, Japan
＊＊〒930-8555 富山市五福 3190
富山大学大学院理工学研究部（理学）Graduate School of Science and Engineering for Research(Science), University of Toyama, 3190 Gofuku, Toyama-shi, Toyama, 930-8555, Japan
〒424-8610 静岡県静岡市清水区折戸 3-20-1
東海大学海洋学部
School of Marine Science and Technology, Tokai University, 3-20-1 Orido, Shimizu-ku, Shizuoka, 424-8610,
Japan

＊＊＊

Corresponding author: Shigekazu Kusumoto
e-mail: kusu@sci.u-toyama.ac.jp

552

酒井望・楠本成寿・清水賀之Fig. 1. Schema of magma flow and numerical experimentparameters, including viscosity. Length of slope is 400 m.

Fig. 2. Schema of magma extrusion, showing the numericalexperiment parameters.

うな構造の違いは，溶岩の物性，噴出率，地形によって支配されることが知られている（例えば，中田，1996）．

たマグマのことを，モデルマグマとよぶ）し，この集合

溶岩の物性や噴出率は，温度に依存するため，地表に噴体を流体として振舞わせる．そのため，粒子の半径をモ出した溶岩形状の議論をこれらのパラメータで行なうこ

デルの規模 (Fig. 1 and 2) に対して十分小さく (0.5-0.6 m)とは難しい．しかしながら，このことを理解した上で，とり，ボンドを 0 とした．

溶岩の物性や噴出率と溶岩形状の関係を議論することは，後のより進んだ研究に有用であろう．

一般に，間隙ばねのばね定数が大きくなると，外力に対し，粒子の運動が鈍くなる．本シミュレーションの場

そこで本研究では，重要ではあるが，条件を複雑にす合，粘性が大きくなる．また，せん断方向の間隙ばねのる温度の効果にはひとまず目をつぶり，マグマの粘性やばね定数を，垂直方向のそれより小さくすると，流動し上昇速度が，地表に噴出する溶岩形状に与える影響を数やすくなる．これらの特性を考慮し，間隙ばねのばね定値シミュレーションにより議論した．数値シミュレー数（垂直方向：kn, せん断方向：ks）の組合せを変化させションには，媒質の大変形や破壊，流動の動的解析が可ることで，モデルマグマの粘性を変化させた．なお本研

能な個別要素法に基づく 2 次元粒状体挙動解析ソフト

究では，粒子の移動速度が小さいため，加速運動をする

ウェア PFC2D ver. 4.0 (Itasca, 2008) を用いた．

粒子のみ減衰させる（収束が速い）Local damping（例えば，Itasca, 2008）を採用している．したがって，間隙ば2．パラメータ設定

ねにダッシュポット要素は含まれていない．

個別要素法とは，粒子同士の接触と力と変位の関係かモデルマグマの粘性ηは，傾斜角θの坂を流れ下るモら接触力を求め，個々の粒子に作用する力からその運動デルマグマの表層速度 v を計測する (Fig. 1) ことで，式を，陽解法を用いて逐次追跡していく数値解析手法であ(1) より推定される（例えば，Turcotte and Schubert, 2002）．

る（例えば，Cundall and Strack, 1979）．2 次元個別要素法
では，物体を弾性ばね（間隙ばね）で繋がれた剛体円盤

η=ρgh sin θ2ν(1)の集合で表現する．本研究で用いる PFC2D では，物体ここで g は重力加速度 (9.81 m/s2) であり，ρはモデルの力学的性質は，剛体円盤を繋ぐ垂直方向とせん断方向マグマの密度である．また h は流れ下るモデルマグマのの間隙ばねのばね定数，剛体円盤間の摩擦係数，さらに厚さである．

固体材料の強度を制御するボンドとよばれるパラメータ

本研究では，剛体円盤間の摩擦係数を 0.6，ρ=2300により表現される．固体材料をモデル化する際には，解kg/m3とし，θ=30°の坂を流れ下るモデルマグマの厚さ

析者が材料に設定したい弾性定数や強度となるよう，計h と表層速度 v を計測した．表層速度 v は，Fig. 1 中の四算機中で 2 軸圧縮試験を行い，これらのパラメータを試角領域を通過する要素の速度を計測することで求めた．

行錯誤的に決定する．

Table 1 に作成したモデルマグマ (Model 1-4) のパラ

本研究ではマグマを剛体円盤粒子の集合体としてモデメータと粘性を示す．粘性より，Model 1 は安山岩質マル化（以降，本シミュレーションで用いるモデル化されグマ，Model 2 と 3 はデイサイト質マグマ，Model 4 は流個別要素法を用いた中-高粘性マグマの噴出シミュレーション紋岩質マグマに相当する．
なお，いずれのモデルマグマにおいても，ポアソン比が 0.5 という流体の基本的性質は維持されている．

553

4．結果と考察モデルマグマが，1×10−3 m/s という速度で地表に押し出された時の溶岩の形状を Fig. 3 に示す．
Fig. 3a と Fig. 3b は，それぞれ，Model 1 と Model 2 の3．シミュレーション方法

モデルマグマが地表に噴出したときの形状であり，Fig.

本研究では，水平な地表面に一定速度でモデルマグマ

3c と Fig. 3d は Model 3 と Model 4 のモデルマグマが噴を噴出させる (Fig. 2)．シミュレーションは二次元で行出したときの形状である．Model 1 と Model 2 は，いずなった．

れも地表に噴出した溶岩は薄く広く分布しているが，火道は地表面に垂直であり，その幅と長さ（深さ）はそれぞれ 20 m と 200 m である．したがって，モデルマ

Model 3 と Model 4 では，モデルマグマは，噴出した場所周辺に留まっている．

グマの噴出は，火道上端から 200 m 分のモデルマグマが溶岩流と溶岩ドームの形状は連続的なものであるが，上方に押し出されることによって実現される．実際のシ溶岩の噴出した高さと広がりの比（アスペクト比）がミュレーションでは，火道底部に水平な壁を作成し，そ1/8，すなわち 0.125 以上となるものを溶岩ドームとするの壁を速度 vupで持ち上げることにより，モデルマグマ区分法がある（例えば，横山・他，1979; 中田，2003）．

を地表に押し出した．なお，重力加速度 (9.81 m/s2) は，この区分法によると，Model 1 と Model 2 は溶岩流を形シミュレーションを行っている間，常に下向きに作用さ成し (Fig. 3a and 3b)，Model 3 と Model 4 は，溶岩ドームせた．

を形成していることになる (Fig. 3c and 3d)．
本シミュレーションで再現された溶岩ドームの形状
は，一定の傾斜をもつのではなく，ドーム頂部付近と裾野

Table 1. Viscosities of magma models and combinationsof normal and shear stiffness of springs between
particles for modeling the viscosities shown.の傾斜が急で，その間の部分は緩やかである (Fig. 3c and3d)．また溶岩ドーム内部は，噴出初期の溶岩が外側に，新しい溶岩が内側に同心円状に分布している．これは，地下からのマグマの上昇によって，当初，上方に向かって成長したドームが，側方に崩れ落ちた痕跡であり，Reyer(1888) によるアナログ実験結果とよく似ている．したがって，本シミュレーションは，溶岩ドーム形成の基本的な部分は再現できていると考えてよいと判断できる．

Fig. 3. Lava dome formation and lava flow at the surface with rising velocity set to vup=1×10−3 m/s. The ratio of length
to breadth is 1 : 1 in all figures. (a) Lava flow, η=4.9×107 Pa s. (b) Lava flow, η=2.0×108 Pa s. (c) Lava dome, η=2.1×109 Pa s. (d) Lava dome, η=1.2×1010 Pa s.

554

酒井望・楠本成寿・清水賀之Fig. 4. Lava dome can be formed by an increase in therising velocity of magma. The ratio of length tobreadth is 1 : 1 in all figures. (a) Lava dome formed
when the rising velocity of magma is 4.5×10−2 m/s(Model 1). (b) Lava dome formed when the risingvelocity of magma is 1×10−2 m/s (Model 2).

Fig. 5. Relationship between viscosity and rising velocityof magma required for forming a lava dome (aspectratio, more than 0.125).

Fig. 3 に示された結果より，モデルマグマの上昇速度が1×10−3 m/s である場合，粘性の高い Model 3 や Model 4は溶岩ドームを形成し，これらより粘性の低い Model 1や Model 2 のモデルマグマは溶岩流を形成することが分かった．しかしながら，実際には Model 1 の粘性に相当する安山岩質マグマの溶岩ドームが存在している．
幾つかの観測事例より，溶岩ドーム形成時のマグマの性の関係は，最小二乗法により，log vup=−1.15 log η+16.8(2)と推定された．
本シミュレーションでは，噴出したマグマが冷えて殻を形成するプロセスが考慮されていない．したがって，噴出率には，2 桁程度の幅のあることが知られている（例本シミュレーションで示された溶岩ドームは，形成途中

えば，中田，1996）．中田 (1996) の第 2 表に示される噴のいわゆるスナップショットを見ていることになる．そ出率 m3 /s を，同表記載の火道径を用いて上昇速度 m/sのため，実際に観測される溶岩ドームの形状との直接比に換算したところ，概ね，1.3×10−3 m/s−6.4×10−1 m/s

較は難しい．実際，中程度の粘性のマグマが，上昇速度であり，SiO2の比率が低い程，上昇速度が大きくなる傾が大きくなることで溶岩ドームが形成され得る可能性が向がみられた．

出てくるというのは直感と合わない．これは，マグマのそこで，マグマの上昇速度として，Model 1 に 4.5×冷却固化の効果が含まれていないため，マグマの上昇が10−2 m/s，Model 2 に 1×10−2 m/s を与えて再計算を行っ早く，見かけ上，溶岩ドームのようになったと考えられたところ，Fig. 4 に示すような溶岩ドームが形成された．

る．

78

これらの結果より，10 Pa s や 10 Pa s 程度の粘性のマグ

本シミュレーション結果は，地表に噴出したマグマのマであっても，上昇速度が大きくなると溶岩ドームが形挙動や，溶岩ドーム形成の初期プロセスの議論に有用で成され得る可能性のあることが示された．

あると思われる．今後は，モデルマグマの温度に応じてそこで，Model 1-4 について，上昇速度を変化させて，間隙ばね定数を変化させ，粘性の変化や殻の形成をモデ

アスペクト比がおよそ 0.125 となるときのモデルマグマルに組み込んでいくことが課題である．本研究は，今後の上昇速度と粘性の関係を調べた．その結果，溶岩ドーの進んだシミュレーションを行う際の基本パラメータ設

ムを形成するためには，低粘性のモデルマグマほど大き定等で，幾らかの指標となり得るであろうと考えている．

な上昇速度が必要とされることが分かった (Fig. 5)．Fig.
5 に示される結果は，粘性が低い程（SiO2の比率が低い5．ま程），上昇速度が大きくなるという，先述の傾向と調和的二次元個別要素法を用いて 107-1010 Pa s オーダーのモとめである．Fig. 5 に示されるデータから，溶岩ドームを形デルマグマの地表への噴出シミュレーションを行った．

成するために必要な上昇速度 (vup) とモデルマグマの粘本シミュレーションではマグマの冷却固化が考慮されて個別要素法を用いた中-高粘性マグマの噴出シミュレーションいないため，得られた結果に直感と合わない部分もあるが，概ね高粘性マグマは溶岩ドームを作りやすく，その内部はこれまでのアナログ実験結果と調和的であった．
また，溶岩ドームと溶岩流の境界の目安である，噴出溶岩のアスペクト比が 1/8 以上となるマグマの上昇速度を調べたところ，溶岩ドームを形成するためには，粘性が低いほど大きな上昇速度が要求され，その速度はマグマの粘性の関数になっていることが示された．
謝辞

査読いただいた，防災科学技術研究所藤田英輔 博士，および産業技術総合研究所 下司信夫博士から，非常に有益な助言やコメントをいただき，本稿は大幅に改善されました．心より感謝致します．また，本稿をご担当下さいました，編集委員青木陽介 博士にも深く感謝致します．
引用文献
Balmforth, N. J., Burbidge, A. S., Craster, R. V., Salzig J. andShen, A. (2000) Visco-plastic models of isothermal lava
domes. Jour. Fluid. Mech., 403, 37-65.
Blake, S. (1990) Viscoplastic models of lava domes, In Lava
flows and domes. Emplacement mechanisms and hazard
implications (Frink, J. H., ed.), Springer-Verlag, New York,88-126.
Blake, S., and Bruno, B. C. (2000) Modelling the emplacement of compound lava flows. Earth Planet. Sci. Lett., 184,181-197.
Cundall, P. A., and Strack, O. D. L. (1979) A discrete numerical model for granular assemblies. Géotechnique, 29, 4765.
Doragoni, M., Borsari, I. and Tallarico, A. (2005) A model forthe shape of lava flow fronts. Jour. Geophys. Res., 110, doi :10.1029/2004JB003523.

555

Filippucci, M., Tallarico, A. and Dragoni, M. (2010) A threedimensional dynamical model for channeled lava flow withnonlinear rheology. Jour. Geophys. Res., 115, doi : 10.1029/2009JB006335.
Fink, J. H., and Griffiths, R. W. (1990) Radial spreading ofviscous gravity currents with solidifying crust. Jour. Fluid.
Mech., 221, 485-509.
Hale, A. J., and Wadge, G. (2008) The transition from endogenous to exogenous growth of lava domes with the development of shear bands. Jour. Volcanol. Geotherm. Res.,
171, 237-257.
Itasca. (2008) PFC2D - Particle Flow Code in 2 Dimensions,ver. 4.0, Minneapolis, Minnesota.
中田節也 (1996) 溶岩ドーム噴火の特徴と普賢岳ドームの成長モデル．地質学論集，46, 139-148.
中田節也 (2003) 噴出物から読みとる噴火．地球科学の新展開 3 マグマダイナミクスと火山噴火（東京大学地震研究所編集 鍵山恒臣 編），朝倉書店，156-179.
Nakada, S., Miyake, Y., Sato, H., Oshima, O., and Fujinawa,A. (1995) Endogenous growth of dacite dome at Unzen
volcano (Japan), 1993-1994. Geology, 23, 157-160.
Reyer, E. (1888) Theoretische Geologie. E. SchweizerbartʼscheVerlagshandlung, Stuttgart, 867p.
Siebert, L., Simkin, T., and Kimberly, P. (2010) Volcanoes ofthe world. 3rd edn., University of California Press, Berkeley
and Los Angeles, California, 551p.
Turcotte, D. L., and Schubert, G. (2002) Geodynamics. 2ndedn., Cambridge University Press, Cambridge, 456p.
Walker, G. P. L., Huntingdon, A. T., Sanders, A. T., andDinsdale, J. L. (1973) Length of lava flows [and discussion]. Phil. Trans. R. Soc. Lond. A, 274, 107-118.
横山 泉・荒巻重雄・中村一明 編 (1979) 火山．岩波講座地球科学，7, 294 p.
Zavada, P., Kratinova, Z., Kusbach, V., and Schulmann, K.
(2009) Internal fabric development in complex lava domes.
Tectonophys., 466, 101-113.
（編集担当青木陽介）