Determination of the mass eruption rate for the 2014 Mount Kelud eruption using three-dimensional numerical simulations of volcanic plumes

Accepted Manuscript Determination of the mass eruption rate for the 2014 Mount Kelud eruption using three-dimensional numerical simulations of volcanic plumes

Y.J. Suzuki, M. Iguchi PII: DOI: Reference:

S0377-0273(17)30365-7 doi: 10.1016/j.jvolgeores.2017.06.011 VOLGEO 6132

To appear in:

Journal of Volcanology and Geothermal Research

Received date: Revised date: Accepted date:

24 August 2016 14 February 2017 5 June 2017

Please cite this article as: Y.J. Suzuki, M. Iguchi , Determination of the mass eruption rate for the 2014 Mount Kelud eruption using three-dimensional numerical simulations of volcanic plumes, Journal of Volcanology and Geothermal Research (2017), doi: 10.1016/ j.jvolgeores.2017.06.011

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Determination of the mass eruption rate for the 2014 Mount Kelud eruption

using three-dimensional numerical simulations of volcanic plumes

Earthquake Research Institute, The University of Tokyo

SC

1-1-1 Yayoi, Bunkyo-ku, Tokyo, 113-0032, Japan

RI

a

PT

Y. J. Suzukia, M. Iguchib

Sakurajima Volcano Research Center, Disaster Prevention Institute, Kyoto

MA

b

NU

[email protected] (corresponding author)

D

University

PT E

1722-19 Sakurajima-Yokoyama, Kagoshima, 891-1419, Japan

AC

CE

[email protected]

1

ACCEPTED MANUSCRIPT

Abstract

Using a three-dimensional fluid dynamic model for eruption cloud dynamics,

we present numerical simulations of the development of a volcanic eruption column

PT

and umbrella cloud formed during the 13 February 2014 eruptions of Mount Kelud,

SC

RI

Indonesia. The model used in this study quantitatively reproduces the observed data

NU

and, in particular, the plume height and the horizontal expansion level and rate of the

umbrella cloud. The simulation results suggest that on the basis of the plume height

MA

and structures and the horizontal expansion level of the umbrella cloud the mass

D

eruption rate (MER) for the 2014 Kelud eruptions was 3 × 107 to 4 × 107 kg s–1. On

PT E

the basis of the horizontal expansion rate of the umbrella cloud, on the other hand, the

CE

estimated MER was 1 × 108 kg s–1. The difference between these two estimates

AC

implies the strong diffusion in the umbrella cloud or the efficient entrainment of

ambient air by eruption column.

Keywords: Mount Kelud; Explosive volcanism; Eruptive plume dynamics; Fluid

dynamic models; Three-dimensional simulations.

2

ACCEPTED MANUSCRIPT

1. Introduction

Explosive volcanic eruptions are characterized by the formation of an eruption

column and umbrella cloud. A mixture of solid pyroclasts and volcanic gases is

PT

ejected from the volcanic vent into the atmosphere and rises as a buoyant plume.

SC

RI

Subsequently, the mixture laterally spreads as an umbrella cloud at the level where its

NU

density becomes equal to the atmospheric density (i.e., the neutral buoyancy level;

NBL). These eruption columns and umbrella clouds transport a large amount of solid

MA

pyroclasts and volcanic gases into the atmosphere. The accumulation of pyroclasts in

D

the atmosphere can affect aviation traffic. Emission of sulfur dioxide gas can be a

PT E

cause of climate change (Robock, 2000). The deposition of pyroclasts on the ground

CE

can impact agriculture, infrastructure, and public health.

AC

In order to reconstruct or predict the progression of volcanic eruptions,

estimations of eruption conditions at the vent from observable quantities such as

plume heights have long been a major focus of research in volcanology (e.g., Costa et

al., 2016). Previous studies using one-dimensional (1D) models (e.g., Woods, 1988;

Bursik, 2001; Woodhouse et al., 2013), two-dimensional (2D) models (e.g., Wohletz

et al., 1984; Neri and Dobran, 1994; Herzog et al., 1998), and 3D models (e.g., Neri et 3

ACCEPTED MANUSCRIPT

al., 2007; Van Eaton et al., 2015; Esposti Ongaro and Cerminara et al., 2016; Suzuki

et al., 2016) investigated eruption cloud dynamics in a stationary environment and in a

wind field. Suzuki and Koyaguchi (2009) performed three-dimensional (3D)

PT

numerical simulations of gigantic eruption clouds observed during the 1991 Pinatubo

RI

eruption to estimate its mass eruption rate (MER). Given that the intensity of this

NU

SC

eruption was strong (MER ~ 109 kg s–1) and the effects of wind on the eruption cloud

dynamics was small, the eruption column rose almost vertically and the umbrella

MA

cloud spread axisymmetrically as a gravity current (Holasek et al., 1996; Costa et al.,

D

2013). The MER was estimated from the maximum height of the eruption column and

PT E

from the height and level of the axisymmetrical spreading umbrella cloud. Suzuki and

CE

Koyaguchi (2013) performed a 3D numerical simulation of a relatively small eruption

AC

cloud observed during the 2011 Shinmoe-dake eruption. In contrast to the 1991 Pinatubo eruption, this eruption was relatively small (MER ~ 106 kg s–1) and the

atmospheric wind was strong. Therefore, the eruption cloud was strongly distorted by

the wind in a single direction. The MER was estimated by comparing the simulated

level of horizontally moving cloud with the observations.

4

ACCEPTED MANUSCRIPT

The eruption conditions at the vent are also estimated from the spreading rate

of umbrella cloud. The spreading of umbrella cloud was expressed by the analytical

models of gravity current (Woods and Kienle, 1994; Holasek et al., 1996; Sparks et al.,

PT

1997). The MER at the vent and/or the volume flux into the umbrella cloud were

SC

RI

estimated from the satellite images, using the analytical model of gravity current

NU

(Pouget et al., 2013) and shallow water model (Pouget et al., 2016). In addition, some

previous studies tried to explain the dynamics of umbrella cloud spreading and

MA

estimate the eruption conditions, using a combination of plume dynamics model,

D

gravity current model, plus sedimentation model (e.g., Bursik et al., 1992; Koyaguchi

PT E

and Ohno, 1992) and a combination of gravity current model plus advection-

CE

diffusion-sedimentation model (Costa et al., 2013).

AC

When the eruption intensity is intermediate between the 1991 Pinatubo

eruption and the 2011 Shinmoe-dake eruption, eruption cloud dynamics are expected

to be complex. The horizontally spreading clouds can be controlled by both the

gravity current of the umbrella cloud and advection due to wind. The 13 February 2014 eruption of Mount Kelud, Indonesia, had an estimated MER of ~107 kg s–1

(Maeno et al., this issue) and is considered to be a typical case of an intermediate 5

ACCEPTED MANUSCRIPT

eruption. On February 13 2014, a major Plinian eruption began at approximately

16:00 UT (23:00 local time). Nakashima et al. (2016) investigated seismic records and

GNSS (Global Navigation Satellite System) data, and interpreted the initial explosion

PT

to have occurred at 15:46 UT, and the Plinian phase to have started at ~16:00 UT and

SC

RI

terminated at ~18:00 UT. They also found that an event occurred at ~16:14 UT that

NU

generated a Rayleigh wave. Images obtained by the MTSAT-1R satellite (Global

Volcanism Program (GVP), 2014) indicate that the plume entered the troposphere at

MA

~16:10 UT. Therefore, the climactic phase of the eruption is considered to start

D

around 16:10–16:15. The eruption column reached ~26 km above sea level (asl), and

PT E

the umbrella cloud laterally spread as a gravity current at a height of 18–20 km asl. At

CE

the same time, the umbrella cloud was affected by an easterly wind and drifted. After

AC

the climactic eruption phase, the eruption cloud moved toward the west or west–

southwest.

In this study, we present 3D numerical simulations for the 2014 Kelud

eruption. We aim to comprehensively explain the data obtained by various

observations and to constrain the eruption conditions at the vent and, in particular, the

mass eruption rate (MER). 6

ACCEPTED MANUSCRIPT

2. Method and simulation inputs

The 3D numerical simulations were designed to reproduce the injection of a

PT

mixture of pyroclasts and volcanic gas from a circular vent located at 1500 m asl.

SC

RI

Atmospheric conditions were based on the meteorological reanalysis data provided by

NU

the Japan Meteorological Agency’s Non-Hydrostatic Model at 16:00 UT on 13

February 2014 (Fig. 1). The temperature profile was typical of tropical atmosphere;

MA

the mean lapse rate was –6.7 K km−1 below 16.5 km, and 2.6 K km−1 above it. These

D

profiles are similar to those obtained by radiosonde measurements.

PT E

We used a combination of a pseudo-gas model for fluid motion and a

CE

Lagrangian model for particle motion (Suzuki and Koyaguchi, 2013). In the pseudo-

AC

gas model, the mixture of the ejected material (solid pyroclasts and volcanic gases)

and the entrained air is treated as a single fluid so that the disequilibrium between the

solid pyroclast and the gas phases is neglected. The fluid dynamic model solves a set

of partial differential equations that describe the conservation of mass, momentum,

and energy, and constitutive equations describing the thermodynamic state of the

mixture of pyroclasts, volcanic gas, and air. Details of the numerical procedures of the 7

ACCEPTED MANUSCRIPT

fluid dynamic model are described in Suzuki et al. (2005). The possible effects of

disequilibrium between the solid pyroclasts and the gas phases are discussed in

Cerminara et al. (2016).

PT

The Lagrangian model for particle motion was employed to calculate ash

SC

RI

particle transport. The particle shape was assumed to be an ideal sphere, and 100–200

NU

marker particles were ejected from the vent every 10 s at the same velocity as the

pseudo-gas. Particle grain sizes are randomly selected within a range of

MA

3.9 × 10−3 mm to 256 mm (i.e., 8  to –8 ). The particle density was assumed to be

D

1500 kg m−3 for particles of 8  to 2 , 1250 kg m−3 for particles of 2  to –2 , and

PT E

1000 kg m−3 for particles of –2  to –8 . The marker particles moved with a relative

CE

velocity to the ambient fluid. The particles’ terminal velocities were dependent on

AC

their Reynolds number. Each particle settled as sediment when it reached the ground

surface. The terminal velocity and particle Reynolds number are described in Suzuki

and Koyaguchi (2013).

Numerical calculations were performed on a generalized coordinate system.

The vent diameter consists of 20 grid points. The grid size increased with distance

from the vent at a constant rate (by a factor of 1.02 for vertical and horizontal 8

ACCEPTED MANUSCRIPT

coordinates) up to the D0, where D0 is the vent diameter. The time-step was

constrained by the CFL (Courant–Friedrich–Lévy) condition using the smallest ratio

of the grid size divided by the characteristic velocity (fluid velocity plus sound

PT

velocity) for the whole domain. For simulations with a background wind, the vent

SC

RI

location was placed off-center to allow for a larger model domain on the downwind

NU

side of the vent. A free-slip condition was applied at the ground surface boundary,

whereas inflow/outflow conditions were implemented at the upper and other

MA

boundaries of the computational domain. The grid numbers and spatial domain for

D

each run are reported in Table 1. For all the simulations, 512 CPUs were used.

, ranged from 5.0 × 105 to 1.0 × 108 kg s−1. The other parameters

CE

1–8). The MER,

PT E

We carried out eight simulations of eruption plumes with variable MER (Runs

AC

were kept fixed in all of the simulations. Magmatic temperature and water content

were assumed to be 1273 K and 5.0 wt.%, respectively. The pressure at the vent was

assumed to be in equilibrium with the atmospheric conditions. According to the

equation of state, the initial density of the ejected material, 2.89 kg m−3 at the vent. The exit velocity,

, was estimated to be

, was 173 m s−1, corresponding to a

Mach number of 1.0. Magmatic temperature, water content, pressure, density, and exit 9

ACCEPTED MANUSCRIPT

velocity were assumed to be constant during eruption. The vent radius was calculated

from the relationship:

SC

RI

1, and the input values for the runs are reported in Table 2.

PT

where R0 ranged from 80 to 252 m. The common parameters used are listed in Table

NU

3. Simulation results

Results of all the 3D simulations are shown in Fig. 2. The eruption columns

MA

exceed 20 km high. The umbrella clouds show asymmetrical spreading; the radial

D

lengths of the downwind umbrella clouds are longer than those of the upwind clouds.

PT E

The downwind umbrella clouds are slightly lower and thicker than the upwind clouds.

CE

We can identify different features in the different MER cases. In Run 1, the

AC

upper cloud edge was almost flat, and the height of the eruption column was at

approximately the same spreading level as the umbrella cloud (Fig. 2a). When the MER is larger (2 × 107 kg s−1; Run 2), the eruption column continued upward to

overshoot the umbrella cloud, and reached a maximum height of ~28 km (Fig. 2b). When the MER is still larger (3.0 ×107 and 3.6 × 107 kg s−1; Runs 3 and 4,

respectively), a parcel of eruption clouds drifted downwind from the overshooting top. 10

ACCEPTED MANUSCRIPT

When the MER is greater than 5 × 107 kg s−1 (Runs 5–7), the overshooting tops

exceeded 30 km. Several layers were developed above 18 km and the cloud became

thicker than 20 km.

PT

Figure 3 illustrates the time evolution of the mass fraction of the erupted

SC

RI

material () along the central axis estimated from the spatial average within 10 grid

points from the central axis. The total height oscillates with time. The maximum

NU

height increases as the MER increases (see also Table 3). The value of  around the

MA

plume top also increases as the MER increases.

D

Figure 4 illustrates how the area of umbrella clouds evolves. The simulation

PT E

results show the asymptotic behavior after 1000 s. For a fixed time, the area of the

CE

umbrella cloud increases as the MER increases. However, the areas of umbrella

AC

clouds are almost indistinguishable, especially between Runs 2 to 5. As shown in the

right panels of Fig. 2, in the range of this parametric study, the depth of umbrella

cloud increases as the MER increases, whereas the area does not increases so much.

In addition, the area of umbrella cloud is controlled not only by the MER but also the

strength of atmospheric stratification and the drag force along the interfaces of the

current (Pouget et al., 2016). 11

ACCEPTED MANUSCRIPT

In the previous studies, the growth of the umbrella cloud was modeled as an

axisymmetrical gravity current, spreading out to form a uniform layer at the NBL

(Woods and Kienle, 1994; Holasek et al., 1996; Sparks et al., 1997). These gravity

PT

current models suggest that the area of umbrella clouds is proportional to t4/3; the

SC

RI

power-low exponent is 1.33. Pouget et al. (2016) reported that the best power-law fit

NU

is ~1.38 at the first stage of eruption (t < 3600 s) and ~0.80 at the later stage (5400 < t

< 9000 s). Using the data of Run 3, we estimated the power-law exponents for

MA

spreading time. Power-law exponents in our simulations are 2.14 from 200 to 1000 s

D

and 1.30 after 1000s. This indicates that the power-low exponent in our simulations is

PT E

substantially greater than the simple gravity model or the analyses by Pouget et al.

CE

(2016) before 1000 s, whereas it is consistent with 1.33 but greater than 1.38 after

AC

1000 s. These results mean that although the simulated umbrella cloud develops

asymmetrically, the asymptotic behavior of its lateral evolution can be approximately described by the simple gravity current model (i.e., the power law of t4/3).

4. Comparison with the observations

12

ACCEPTED MANUSCRIPT

We have presented 3D numerical simulation results of the eruption column

and umbrella cloud development under the same atmospheric conditions as the 2014

Kelud eruption. In this section, we compare the simulation results with the data

PT

obtained from the observation of the 2014 Kelud eruption, and estimate the MER that

SC

RI

can comprehensively explain the observations.

NU

The level of lateral spreading of the umbrella cloud was measured by CALIOP

LIDAR from onboard the CALIPSO satellite (GVP, 2014). The CALIPSO orbit at

MA

this time was ~30 km west of the volcanic vent. The CALIOP measured its strongest

D

signal in altitudes around 17–18 km asl., with the top of the umbrella cloud at 18–

PT E

19 km asl. (Kristiansen et al., 2015). The results of Runs 1–4 show the horizontal

CE

spreading of the umbrella cloud at approximately 18 km, whereas those of Runs 5–7

AC

show the thick layer of the umbrella cloud spreading above 20 km (Fig. 2).

The CALIOP measurements also indicated some higher spreading layers at

approximately 22 and 25 km (Kristiansen et al., 2015). The results of Runs 3–7 show

multiple layers of eruption cloud, whereas those of Runs 1 and 2 indicate only one

layer at approximately 18 km (left panels in Fig. 2).

13

ACCEPTED MANUSCRIPT

The CALIOP data revealed that some clouds reached a maximum height of

~26 km asl. Because the CALIPSO did not pass directly over the vent, this

comparison between the simulation results and the CALIOP data is subsidiary. The

PT

maximum height in Run 1 is lower than the observed height, and those in Runs 5–7

SC

RI

exceed 30 km and are significantly higher than the observed height. The maximum

NU

heights in Runs 2–4 are ~25 km asl, which are consistent with the observations.

The area of the umbrella cloud can be estimated from the observations of the

MA

MTSAT satellite (Fig. 4). The areas of the umbrella clouds obtained from all the runs

D

show the similar behavior as the observations. In detail, the result of Run 7 seems

PT E

more consistent with the observations, whereas those of Runs 1–6 are slightly smaller

CE

than the observations.

AC

Through comparison, we conclude that the results of Runs 3 and 4 (MERs of 3.0 × 107 and 3.6 × 107 kg s−1, respectively) explain the observed plume height, and

the maximum spreading level of the umbrella cloud, whereas those of Run 7 (MER of 1.0 × 108 kg s−1) explain the observed expansion rate of the umbrella cloud. We also

qualitatively compare the dispersal of marker particles with the field observations.

Figure 5 indicates that the distribution of marker particles obtained by the simulation 14

ACCEPTED MANUSCRIPT

of Run 3 is similar to the satellite images. The cloud expanded almost

axisymmetrically in the 60 min after the onset of eruption. After that, the cloud was

elongated by the easterly wind. The cloud expanded ~40 km upwind forming a

PT

stagnation point. At this point, the radial expansion velocity was considered to be

SC

RI

equal to the wind velocity.

NU

Figure 6 shows the depositional pattern of the marker particles at the ground.

The main axis of the deposits is west, which roughly coincides with the field

MA

observations (Maeno et al., this issue). However, our simulations did not reproduce

D

the north and northeast deposits from the volcano. In our simulations, the coarse

PT E

particles separated from the eruption column and settled near the volcano, whereas the

CE

fine particles were affected by both the umbrella cloud spreading and the easterly

AC

wind, and settled west of the volcano. If the eruption is weak and the fine particles

separated from the small eruption column, then the southerly and southwesterly wind

at lower levels can cause them to drift and settle north and northeast of the volcano (cf. Tanaka et al., 2016). Even when the MER is small (5.0 × 105 kg s−1; Run 8; results not

shown), the fine particles were transported up to ~10 km asl and drifted toward the

west. The deposits observed in the north and northeast area of the volcano may have 15

ACCEPTED MANUSCRIPT

been formed by particles separating from relatively low clouds, such as a co-

ignimbrite ash cloud generated by pyroclastic flows. Indeed, some pyroclastic flows

were observed during the late stages of the Kelud eruption.

PT

Our simulation results suggest that a MER of 3 × 107 to 4 × 107 kg s−1 can

SC

RI

explain the observed data from the Kelud 2014 eruption from the viewpoint of the

plume height, the expansion level of umbrella cloud, and structure of cloud, whereas a

NU

MER of 1 × 108 kg s−1 can explain the observed data from the viewpoint of the

MA

expansion rate of umbrella. These estimates of MER are consistent with the

D

sedimentological studies; Maeno et al. (this issue) estimated an MER of

PT E

6.5 ± 2.8 × 107 kg s−1 by dividing the total mass of the deposits (3.3–6.6 × 1011 kg) by

CE

the duration of eruption (2.0–2.5 h). There are two possible explanations for the

AC

difference between our two estimations (i.e., 3–4 × 107 and 1 × 108 kg s−1). One is the

diffusion in the atmosphere. In the present 3D model, the small-scale turbulence into

the atmosphere is not taken into account. The strong diffusivity may lead to a rapid expansion of umbrella cloud even when the MER is 3–4 × 107 kg s−1. Another

possible explanation for the difference is related to the entrainment of ambient air by

eruption column. From the analytical study (Morton et al., 1956), the characteristic 16

ACCEPTED MANUSCRIPT

height of the eruption column (H) and the volumetric flow rate into the umbrella ． cloud (V) are approximately give as a function of the efficiency of entrainment (k), the

(1)

RI

PT

． MER (m0), and the stratification of the atmosphere (N) as

SC

where C1 and C2 are the proportionality constants (see also Suzuki and Koyaguchi,

NU

． 2009). Equation (1) indicates that H decreases as k increases, whereas V increases

MA

with k. The small-scale disturbance that cannot be captured by the reanalysis data of

atmospheric conditions may cause the efficient mixing between eruption cloud and

PT E

D

ambient air. In addition, the disequilibrium with the atmospheric pressure and the

supersonic near the vent can change the flow structures in a lower region of eruption

CE

column (Ogden et al., 2008; Carcano et al., 2013) and the efficiency of entrainment in

AC

a whole region of column (Koyaguchi et al., in prep).

5. Concluding remarks

We carried out 3D numerical simulations of the eruption cloud observed

during the 2014 Kelud eruption. The simulation results indicate that the mass eruption rate (MER) of this eruption was 3 × 107 to 4 × 107 kg s−1 from the viewpoint of the

17

ACCEPTED MANUSCRIPT

plume height, level of umbrella cloud spreading, and structure of cloud, whereas the MER was 1 × 108 kg s−1 from the viewpoint of the expansion rate of umbrella cloud.

These values are consistent with sedimentological studies. The difference between

PT

two estimations implies the strong diffusion in the atmosphere or the efficient

SC

RI

entrainment of ambient air by eruption cloud. To clarify this difference, a systematic

NU

investigation using the model in which the small scale deviation in atmosphere is

taken into account is required. Moreover, to quantitatively reproduce depositional

MA

patterns, advection–diffusion simulations of tephra combined with a three-

PT E

AC

CE

in future studies.

D

dimensional fluid dynamic model, such as our simulation model, should be performed

18

ACCEPTED MANUSCRIPT

Acknowledgements

This study was supported by SATREPS of JST and JICA, the Ministry of Education

PT

Culture, Sports, Science and Technology (MEXT) of Japan, under its Observation and

RI

Research Program for the Prediction of Earthquakes and Volcanic Eruptions. The

SC

computations were carried out in part on the Earth Simulator at the Japan Agency for

NU

Marine-Earth Science and Technology (JAMSTEC) and also on the Primergy

MA

RX200S6 at the Research Computer System, Kyushu University. We thank the

Meteorological Research Institute and Toshiki Shimbori for providing the satellite

PT E

D

images from the Earthquake Research Institute cooperative research program: “Development of volcanic ash transport model integrating meteorological model and

CE

three-dimensional volcanic eruption dynamics model”. The original data for the

AC

isopach in Fig. 6 were provided by Dr. F. Maeno. We wish to thank A. Van Eaton, A.

Costa, T. Koyaguchi, and anonymous reviewers for constructive suggestions that

improved the quality of the manuscript.

19

ACCEPTED MANUSCRIPT

References

Bursik, M., 2001. Effect of wind on the rise height of volcanic plumes. Geophys. Res.

Lett. 28(18), 3621–3624. http://dx.doi.org/10.1029/2001GL013393.

PT

Bursik, M.I., Carey, S.N., Sparks, R.S.J., A gravity current model for the May 18,

SC

RI

1980 Mount St Helens plume. Geophys. Res. Lett. 19(16), 1663–1666.

NU

Carcano, S., Bonaventura, L., Esposti Ongaro, T., Neri, A., A semi-implicit, second-

order-accurate numerical model for multiphase underexpanded volcanic jets.

MA

Geosci. Model Dev. 6, 1905–1924. http://dx.doi.org/10.5194/gmd-6-1905-2013.

D

Cerminara, M., Esposti Ongaro, T., Berselli, L.C., 2016. ASHEE 1.0: a compressible,

PT E

equilibrium-Eulerian model for volcanic ash plumes. Geosci. Model Dev. 9, 697–

CE

730. http://dx.doi.org/10.5194/gmd-9-697-2016.

AC

Cerminara, M., Esposti Ongaro, T., Neri A., 2016. Large eddy simulation of gas-

particle kinematic decoupling and turbulent entrainment in volcanic plumes. J.

Volcanol. Geotherm. Res. 326, 143–171.

http://dx.doi.org/10.1016/j.jvolgeores.2016.06.018.

20

ACCEPTED MANUSCRIPT

Costa, A., Folch, A., Macedonio, G., 2013. Density-driven transport in the umbrella

region of volcanic clouds: Implications for tephra dispersion models. Geophys.

Res. Lett. 40(10), 4823–4827. http://dx.doi.org/10.1002/grl.50492.

PT

Costa, A., Suzuki, Y.J., Cerminara, M., Devenish, B., Esposti Ongaro, T., Herzog, M.,

SC

RI

Van Eaton, A.R., Denby, L.C., Bursik, M., de’Michieli Vitturi, M., Engwell, S.,

NU

Neri, A., Barsotti, S., Folch, A., Macedonio, G., Girault, F., Carazzo, G., Tait, S.,

Kaminski, E., Mastin, L.G., Woodhouse, M.J., Phillips, J.C., Hogg, A.J.,

MA

Degruyter, W., Bonadonna, C., 2016. Results of the eruption column model inter-

D

comparison exercise. J. Volcanol. Geotherm. Res. 326, 2–25.

PT E

http://dx.doi.org/10.1016/j.jvolgeores.2016.01.017.

CE

Esposti Ongaro, T., Cerminara, M., 2016. Non-equilibrium processes in ash-laden

AC

volcanic plumes: new insights from 3D multiphase flow simulations. J. Volcanol.

Geotherm. Res. 326, 127–142. http://dx.doi.org/10.1016/j.jvolgeores.2016.04.004.

Global Volcanism Program, 2014. Smithsonian Institution, Kelud Eruption Page.

http://www.volcano.si.edu/volcano.cfm?vn=263280.

Herzog, M., Graf, H.-F., 2010. Applying the three-dimensional model ATHAM to

volcanic plumes: Dynamics of large co-ignimbrite eruptions and associated 21

ACCEPTED MANUSCRIPT

injection heights for volcanic gases. Geophys. Res. Lett. 37(18), L19807.

http://dx.doi.org/10.1029/2010GL044986.

Herzog, M., Graf, H.-F., Textor, C., Oberhuber, J.M., 1998. The effect of phase

PT

changes of water on the development of volcanic plumes. J. Volcanol. Geotherm.

SC

RI

Res. 87, 55–74. http://dx.doi.org/10.1016/S0377-0273(98)00100-0.

NU

Holasek, R., Self, S., Woods, A.W., 1996. Satellite observations and interpretation of

the 1991 Mount Pinatubo eruption plumes. J. Geophys. Res. 101(B12), 27635–

MA

27655. http://dx.doi.org/10.1029/96JB01179.

D

Koyaguchi, T., Ohno, M., 2001. Reconstruction of eruption column dynamics on the

PT E

basis of grain size of tephra fall deposits 2. Application to the Pinatubo 1991

CE

eruption. J. Geophys. Res. 106(B4), 6513 – 6533.

AC

Kristiansen, N.I., Prata, A.J., Stohl, A., Carn, S.A., 2015. Stratospheric volcanic ash

emissions from the 13 February 2014 Kelud eruption. Geophys. Res. Lett. 42(2),

588–596. http://dx.doi.org/10.1002/2014GL062307.

Maeno, F., Nakada, S., Yoshimoto, M., Shimano, T., Hokanishi, N., Zaennudin, A.,

Iguchi, M., this issue. A sequence of a plinian eruption preceded by dome

destruction at Kelud volcano, Indonesia, on February 13, 2014, revealed from 22

ACCEPTED MANUSCRIPT

tephra fallout and pyroclastic density current deposits. J. Volcanol. Geotherm.

Res.

Morton, B.R., Taylor, G.I., Turner, J.S., 1956. Turbulent gravitational convection

PT

from maintained and instantaneous sources. Proc. R. Soc. London, Ser. A,

SC

RI

234(1196), 1–23. http://dx.doi.org/10.1098/rspa.1956.0011.

NU

Nakashima, Y., Heki, K., Takeo, A., Cahyadi, M.N., Aditiya, A., Yoshizawa, K.,

2016. Atmospheric resonant oscillations by the 2014 eruption of the Kelud

MA

volcano, Indonesia, observed with the ionospheric total electron contents and

D

seismic signals. Earth Planet. Sci. Lett. 434, 112–116.

PT E

http://dx.doi.org/10.1016/j.epsl.2015.11.029.

CE

Neri, A., Dobran, F., 1994. Influence of eruption parameters on the thermofluid

AC

dynamics of collapsing volcanic columns. J. Geophys. Res. 99 (B6), 11833–

11857. http://dx.doi.org/10.1029/94JB00471. Neri, A., Esposti Ongaro, T., Menconi, G., De’Michieli Vitturi, M., Cavazzoni, C.,

Erbacci, G., Baxter, P.J., 2007. 4D simulation of explosive eruption dynamics at

Vesuvius. Geophys. Res. Lett. 34(4), L04309.

http://dx.doi.org/10.1029/2006GL028597. 23

ACCEPTED MANUSCRIPT

Ogden, D., Glatzmaier, G.A., Wohletz, K.H., 2008. Effects of vent overpressure on

buoyant eruption columns: implications for plume stability. Earth Planet. Sci. Lett.

268, 283–292. http://dx.doi.org/10.1016/j.epsl.2008.01.014.

PT

Pouget, S., Bursik, M., Webley, P., Dehn, J., Pavolonis, M., 2013. Estimation of

NU

J. Volcanol. Geotherm. Res. 258, 100–112.

SC

RI

eruption source parameters from umbrella cloud or downwind plume growth rate.

http://dx.doi.org/10.1016/j.jvolgeores.2013.04.002.

MA

Pouget, S., Bursik, M., Johnson, C.G., Hogg, A.J., Phillips, J.C., Sparks, R.S.J., 2016.

D

Interpretation of umbrella cloud growth and morphology: implications for flow

PT E

regimes of short-lived and long-lived eruptions. Bull. Volcanol. 71:1.

CE

http://dx.doi.org/10.1007/s00445-015-0993-0.

AC

Robock, A., 2000. Volcanic eruptions and climate. Rev. Geophys. 38:2, 191–219.

http://dx.doi.org/10.1029/1998RG000054.

Sparks, R.S.J., Bursik, M.I., Carey, S.N., Gilbert, J.S., Glaze, L.S., Sigurdsson, H.,

Woods, A.W., 1997. Volcanic Plumes. John Wiley & Sons, Chichester. 574 pp.

24

ACCEPTED MANUSCRIPT

Suzuki, Y.J., Koyaguchi, T., 2009. A three-dimensional numerical simulation of

spreading umbrella clouds. J. Geophys. Res. 114, B03209.

http://dx.doi.org/10.1029/2007JB005369.

PT

Suzuki, Y.J., Koyaguchi, T., 2013. 3D numerical simulation of volcanic eruption

SC

RI

clouds during the 2011 Shinmoe-dake eruptions. Earth Planets Space 65, 581–

NU

589. http://dx.doi.org/10.5047/eps.2013.03.009.

Suzuki, Y.J., Koyaguchi, T., Ogawa, M, Hachisu, I., 2005. A numerical study of

MA

turbulent mixing in eruption clouds using a three-dimensional fluid dynamics

D

model. J. Geophys. Res. 110, B08201. http://dx.doi.org/10.1029/2004JB003460.

PT E

Suzuki, Y.J., Costa, A., Koyaguchi, T., 2016. On the relationship between eruption

CE

intensity and volcanic plume height: Insights from three-dimensional numerical

AC

simulations. J. Volcanol. Geotherm. Res. 326, 120–126.

http://dx.doi.org/10.1016/j.jvolgeores.2016.04.016.

Tanaka, H.L., Iguchi, M., Nakada, S., 2016. Numerical simulations of volcanic ash

plume dispersal from Kelud volcano in Indonesia on February 13, 2014. J.

Disaster Res. 11, 1, 31–42.

25

ACCEPTED MANUSCRIPT

Van Eaton, A.R., Mastin, L.G., Herzog, M., Schwaiger, H.-F., Schneider, D.J.,

Wallace, K.L.W., Clarke, A.B., 2015. Hail formation triggers rapid ash

aggregation in volcanic plumes. Nature Communications 6:7860.

PT

http://dx.doi.org/10.1038/ncomms8860.

SC

RI

Wohletz, K.H., McGetchin, T.R., Stanford II, M.T., Jones, E.M., 1984.

Geophys. Res. 89 (B10), 8269–8285.

NU

Hydrodynamic aspects of caldera-forming eruptions: Numerical models. J.

MA

Woodhouse, M.J., Hogg, A.J., Phillips, J.C., Sparks, R.S.J., 2013. Interaction between

D

volcanic plumes and wind during the 2010 Eyjafjallajökull eruption, Iceland. J.

PT E

Geophys. Res. Solid Earth 118, 92–109.

CE

http://dx.doi.org/10.1029/2012JB009592.

AC

Woods, A.W., 1988. The fluid dynamics and thermodynamics of eruption columns.

Bull. Volcanol. 50, 169–193. http://dx.doi.org/10.1007/BF01079681.

Woods, A.W., Kienle, J., 1994. The dynamics and thermodynamics of volcanic

clouds: Theory and observations from the April 15 and April 21, 1990 eruptions

of Redoubt Volcano, Alaska. J. Volcanol. Geotherm. Res. 62, 273–299.

http://dx.doi.org/10.1016/0377-0273(94)90037-X. 26

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

27

ACCEPTED MANUSCRIPT

Tables

Table 1. Common input parameters and constants for simulations. Value

Vent elevation

PT

Variable

1500 m

1273 K

Exit water fraction

0.05

SC

Exit temperature

RI

173 m s−1

Exit velocity

2.89 kg m−3

Exit density

9.81 m s−2

NU

Gravity body force

462 J kg−1 K−1

Gas constant of atmospheric air

287 J kg−1 K−1

Specific heat of solid pyroclasts

1100 J kg−1 K−1

MA

Gas constant of volcanic gas

1348 J kg−1 K−1

Specific heat of volcanic gas

D

at constant volume

717 J kg−1 K−1

PT E

Specific heat of air

AC

CE

at constant volume

28

ACCEPTED MANUSCRIPT

Table 2. Input parameters for simulations. MER

Vent radius

(kg s−1)

(m)

Grid number

Domain

Final time

x×y×z

(km3)

(sec)

1.0 × 107

80

1520 × 1040 × 510

121 × 83 × 35

1800

2

2.0 × 107

113

1040 × 1040 × 320

117 × 117 × 29

2057

3

3.0 × 107

138

2960 × 2000 × 382

408 × 275 × 42

7200

4

3.6 × 107

151

1000 × 1000 × 312

151 × 151 × 37

3403

5

5.0 × 10

7

178

1040 × 1040 × 312

185 × 185 × 44

3228

6

7.0 × 107

211

1000 × 1000 × 280

210 × 210 × 45

3674

7

1.0 × 108

252

1456 × 1456 × 320

366 × 366 × 64

3050

8

5.0 × 105

18

1016 × 1016 × 1016

18 × 18 × 17

862

AC

CE

PT E

D

MA

NU

SC

RI

1

PT

Run

29

ACCEPTED MANUSCRIPT

Table 3. Comparison between the simulation results and observed data. Maximum

heights of eruption columns were estimated from the level where the mass fraction of

the ejected material is 0.01, using the time window 500–1800 s. The levels of

PT

umbrella cloud were determined as center of the cloud at x = –30 km for direct

SC

RI

comparison with the CALIPSO pass (see also the right panels in Fig. 2). The results

consistent with the observations are highlighted in grey. Level of

(kg s-1)

umbrella (km)

Maximum

Multiple

Area of

height (km)

layer

umbrella

20.7–23.6

–

< obs.

23.6–27.8

–

< obs.

25.6–31.7

✔

< obs.

25.6–33.0

✔

< obs.

NU

MER

~18

2

2.0 × 10

7

~18

3

3.0 × 107

~18

4

3.6 × 107

~18

5

5.0 × 107

> 20

30.8–40.0

✔

< obs.

6

7.0 × 107

> 20

38.4–47.0

✔

< obs.

7

1.0 × 108

> 20

46.9–56.1

✔

~ obs.

8

5.0 × 105

~8

~11

–

< obs.

CE

AC

MA

1.0 × 107

PT E

1

D

Run

30

ACCEPTED MANUSCRIPT

Figure captions

Fig. 1. Initial atmospheric conditions used for simulations: (a) temperature, (b) wind

velocity from west to east, and (c) wind velocity from south to north. Solid curves: the

PT

reanalysis data calculated by the Global Spectral Model of the Japan Meteorological

SC

RI

Agency. Dashed curves: the radiosonde data in Surabaya. The altitude of the volcanic

NU

vent is 1.5 km asl, as illustrated by the horizontal lines.

MA

Fig. 2. Simulation results of (a) Run 1, (b) Run 2, (c) Run 3, (d) Run 4, (e) Run 5, (f)

D

Run 6, and (g) Run 7. Left column figures show iso-surfaces of magma mass fraction

PT E

 = 0.02 at 1800 s from eruption initiation and their vertical cross-sections at x = –

CE

30 km. Right column figures show the vertical cross-sections of the mass fraction of

AC

the erupted material at y = 0 km. Parameters used and simulation settings are given in

Tables 1 and 2.

Fig. 3. Timewise distributions of the mass fraction of the erupted material for Runs 1–

7.

31

ACCEPTED MANUSCRIPT

Fig. 4. Umbrella cloud evolution for Runs 1–7 (curves). The area of umbrella cloud are calculated from the region where  > 10-4 and z > 10 km. The observed data from

the Kelud 2014 eruption are also plotted. To plot the observed data, we determined

PT

the area of umbrella cloud obtained from the MTSAT-1R images. The images were

SC

RI

converted and expressed by 256 gradation (brightness). We counted the number of

NU

pixels where the brightness is larger than the threshold, and calculated the area from

the number of pixels. We tentatively assumed that the threshold values are 200 and

MA

256, which represent the minimum and maximum values on the vertical error bars,

D

respectively. To determine the time from the eruption beginning, we assume that the

PT E

climactic phase of eruption started at 16:10. The horizontal error bars represent ±600

CE

s. Numbers next to curves are the best power-law best fits using the data of Run 3

AC

(2.14 from 200 to 1000 s and 1.30 from 1000 to 7200 s).

Fig. 5. Horizontal expansion of umbrella clouds obtained from the simulation of Run

3. Top views of the marker particles are shown for (a-1) 0 min, (a-2) 10 min, (a-3)

30 min, (a-4) 50 min, (a-5) 60 min, (a-6) 70 min, (a-7) 80 min, (a-8) 90 min, (a-9)

100 min, and (a-10) 120 min after eruption initiation. For comparison, satellite images 32

ACCEPTED MANUSCRIPT

on 13 February 2014 at (b-1) 16:13, (b-2) 16:24, (b-3) 16:47, (b-4) 17:03, (b-5) 17:13,

(b-6) 17:23, (b-7) 17:33, (b-8) 17:44, (b-9) 17:56, and (b-10) 18:14 UTC are shown.

Note that the time represented in panels (b) are not exact because the duration of

SC

RI

PT

capture was several minutes.

NU

Fig. 6. Depositional pattern of marker particles at 7200 s obtained from simulation

Run 3. Curves represent isopachs of 0.1, 1.0, 4.0, and 15.0 cm for the Kelud 2014

AC

CE

PT E

D

MA

eruption, as observed in the field (Maeno et al., this issue).

33

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

Figure 1

34

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Figure 2

35

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Figure 3

36

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC

CE

Figure 4

37

AC

Figure 5

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

38

RI

PT

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

NU

SC

Figure 6

39

ACCEPTED MANUSCRIPT

Highlights

We performed 3D numerical simulations for the 2014 Kelud eruption.

-

Simulation results agree well with the various observations.

-

The mass eruption rate for this eruption is estimated to be 0.3 – 1.0 × 108 kg s-1.

AC

CE

PT E

D

MA

NU

SC

RI

PT

-

40