3-D Numerical Simulation of Volcanic Eruption Clouds
Yujiro J. Suzuki 1, Takehiro Koyaguchi 111

Earthquake Research Institute, The University of Tokyo, Japan

Introduction

During explosive volcanic eruptions, a mixture of hot solid
particles (volcanic ashes) and volcanic gas is released from thevolcanic vent into the atmosphere. Such events are characterized
by the formation of buoyant plumes and/or density currents (Fig.
1). Generally, the ejected material (i.e., volcanic ashes and gas)
have an initial density of several times as large as theatmospheric density since it contains more than 90 wt.%
volcanic ashes at the vent [1]. As the ejected material is mixed
with ambient air, the density of the mixture drastically decreases
and becomes less than the atmospheric density, because theentrained air expands by heating from the hot ashes. When theejected material entrains sufficient air to become buoyant, alarge turbulent plume (“eruption column”) rises up to a height ofseveral tens of kilometers. On the other hand, if the ejectedmaterial does not entrain sufficient air and its vertical velocity
falls to zero before the eruption cloud becomes buoyant, aneruption column collapses and the heavy and hot cloud spreads
radially as a density current (“pyroclastic flow”). The upper
part of the pyroclastic flow entrains the ambient air and forms
another type of buoyant plume (“co-ignimbrite ash cloud”).
Because the impact and type of volcanic hazards are largely
different between the two eruption styles (i.e., eruption column
and pyroclastic flow), it has been a central subject ofvolcanology and hazard science to quantitatively predict thecondition where an eruption column collapses to generate apyroclastic flow.
The plume such as an eruption column and/or co-ignimbrite
ash cloud exhausts its thermal energy and loses its buoyancy
with the stratified atmosphere. At the neutral buoyancy level(NBL) where the plume density is equal to that of theatmosphere, the eruption cloud spreads laterally as a density
current (“umbrella cloud”) [1]. If the plume has much greater
vertical velocity than the speed of cross-wind in the atmosphere,
its trajectory is not significantly controlled by the wind and rises
vertically (Fig. 1a). On the other hand, if the vertical velocity ismuch smaller than the cross-wind speed, the plume is largely
distorted by the cross-wind and shows a bent-over trajectory(Fig. 1b); this effect can modify the efficiency of turbulent
mixing, and hence, the maximum eruption column height.
Because the height of eruption column is one of the few
available data for estimating the eruption conditions, it is of
critical importance to clarify the relationship between theeruption condition and column height under the cross-windconditions.
The dynamics of volcanic eruption clouds are largely
controlled by the efficiency of turbulent mixing between theeruption cloud and ambient air, as described above. Previousstudies using steady one-dimensional (1-D) models of eruption
clouds [e.g., 1] described the behaviour of turbulent mixing in a
simple form based on the theoretical work of pure jet/plumes [2].
Unlike pure jet/plumes, the driving force of eruption clouds
changes with height because of the change in density. The flownear the vent was negative buoyancy and is driven by the initialmomentum. As the eruption column rises, the density of thecloud decreases and the flow is primarily driven by positive
buoyancy. In this study, we aim to develop a numerical modelof eruption cloud which can reproduce the density change of thecloud and the turbulent mixing correctly.

Fig. 1. Schematic illustrations of volcanic eruption
clouds.

2

Model Description

The numerical model of eruption cloud is based on the workof Suzuki et al. [3]. The model is designed to describe theinjection of a gas-ash mixture from a circular vent above a flat
surface in a stratified atmosphere. The physical domain involves
a vertical and horizontal extent of several tens of kilometers. At
the ground boundary, the free-slip condition is assumed for thevelocities of the ejected material and air. At the upper and otherboundaries of computational domain, the fluxes of mass,momentum, and energy are assumed to be continuous. We
assume steady conditions; for each run, vent radius, exit velocity,
initial temperature, and initial mass fraction of gas content arefixed. In this study, we adopt a pseudo gas model; we ignore theseparation of volcanic ashes from the eruption cloud. The
momentum and heat exchanges between the volcanic ashes andgas are assumed to be so rapid that the velocity and temperature
are the same for all phases.
2.1
Equation of StateThe density of eruption clouds varies nonlinearly with themixing ratio between the ejected material and air. We reproduce
the nonlinear feature of mixture density by changing theeffective gas constant and the effective specific heat of themixture asR m  n a Ra  n g R g ,(1)Cvm  na Cva  n g Cvg  ns C vs ,(2)

where R is the gas constant, n is the mass fraction, and Cv
is the specific heat at constant volume. The subscripts m,a, g, and s refer to the mixture, air, volcanic gas, and solid
particles (ashes), respectively. The mass fraction of air(na), volcanic gas (ng), and ashes (ns) satisfy the conditionof na + ng + ns = 1. Using these procedures, the equationof state for the mixture of the ejected material and air canbe approximated by that for an ideal gas asp  R m T ,(3)where p, , T are the pressure, density, and temperature of themixture. On the assumption of ideal gas, the governing equation
is expressed by the Euler equation for compressible flow.

spreads horizontally along the wind direction around the heightof 5-8 km.

2.2
Numerical Procedure
The partial differential equations are solved numerically by theRoe scheme [4] in space, which is a general total variation
diminishing (TVD) scheme for compressible flow and cansimulate a generation of shock waves inside and around thehigh-speed jet correctly. The MUSCL method [5] is applied tointerpolate the fluxes between grid points, and therefore ournumerical model achieves third-order accuracy in space. Theseequations are solved using the time splitting method.
In order to correctly reproduce the general feature of turbulent
mixing that the efficiency of entrainment is independent of theReynolds number [6], it is essential to apply 3-D coordinateswith a sufficiently high spatial resolution [3]. In this study, thecalculations were performed on a 3-D domain with a nonuniform grid. The grid size was set to be sufficiently smaller
than L0/8 near the vent, where L0 is the vent radius. In order toeffectively simulate the turbulent mixing with high spatial
resolutions both far from and close to the vent, the grid size
increases at a constant rate (by a factor of 1.02) up to L0 with thedistance from the vent.
3Fig.2 Numerical simulation result of eruption column
and umbrella cloud in Run 1. Cross-sectional distribution of
the mass fraction of the ejected material in x-z space atResults

We carried out three simulations for explosive eruptions. Table
1 lists the initial conditions and numerical conditions for thesimulations. In all the runs, an initial temperature of T0 = 1000 K
and initial mass fraction of volcanic gas of ng0 = 0.0284.

Fig.3 Numerical simulation result of bent-over eruption
plume in Run 3. Cross-sectional distribution of the massfraction of the ejected material in x-z space in 406 sec isshown.

Table 1 Input Parameters of the Simulations
m0 (kg s-1)

w0 (m s-1)

Simulation

L0 (m)Run 1

45

4.0 x 106

180Run 2189

1.0 x 108115

276115Run 3

2.0 x 10

Our simulations have reproduced the behaviour of eruption
clouds including eruption columns and/or the formation of
pyroclastic flows and the unsteady and multidimensionalfeatures of eruption clouds.
Run 1 shows the typical result of eruption column (Fig. 2). The
eruption cloud is ejected from the vent as a turbulent jet. As it
rises and entrains ambient air, the density of the cloud becomesless than that of air and then a buoyant plume rises. The buoyant
plume is high unstable as it ascends, and undergoes meandering
instability that induces efficient mixing so that its radial scale
gradually increases with height. After the eruption column
reaches up to a height of 14 km, the eruption cloud flows down
to the NBL (~10 km) and spreads laterally as an umbrella cloud.
When the vent radius is larger (189 m in Run 2), a buoyant
plume and pyroclastic flow develop simultaneously (figure isnot shown). Some parts of the mixture remain heavier than air,whereas the others become lighter. The heavier parts collapse tothe ground and spread radially as a pyroclastic flow. The lighterparts of the mixture continue to rise as a buoyant plume andform an umbrella cloud. The upper region of the pyroclastic
flow entrains air and forms co-ignimbrite ash clouds which
subsequently join the buoyant plume.
Simulation of Run 3 reproduces the bent-over plume oferuption cloud (Fig. 3). Near the vent, the eruption cloud rises
vertically. Above 3km, the cloud is largely distorted by thecross-wind. Then, the eruption cloud stops to rise and eventually

4Conclusion

In this study, we have developed a new numerical model ofthe dynamics of eruption clouds. Employing 3-D coordinates, a
3rd order accuracy scheme and a fine grid size, our model cansuccessfully reproduce an eruption column, pyroclastic flow,umbrella cloud, and bent-over plume. These results enable us toestimate the relationship between the eruption conditions andthe observable data.
Acknowledgements
We thank Earth Simulator Center, Japan Agency for
Marine-Earth Science and Technology, and ResearchInstitute for Imformation Technology, Kyushu
University. Part of this study was supported by theERI Cooperative Research Program.
References[1][2]

[3]

[4]
[5]

[6]

A. W. Woods, The fluid dynamics and thermodynamics oferuption columns, Bull. Volcanol., 50, 169-193, 1998.
B. R. Morton, G. I. Taylor, J. S. Turner, Turbulent gravitational
convection from maintained and instantaneous sources, Proc. R.
Soc. London, Ser. A, 234(1196), 1-23, 1956.
Y. J. Suzuki, T. Koyaguchi, M. Ogawa, I. Hachisu, A numerical
study of turbulent mixing in eruption clouds using a threedimensional fluid dynamics model, J. Geophys. Res., 110, B08201,2005.
P. L. Roe, Approximate Riemann Solvers, parameter vectors, anddifference schemes, J. Comput. Phys., 43, 357-372, 1981.
B. van Leer, Towards the ultimate conservative difference scheme
III. Upstream-centered finite-difference scheme for ideal
compressible flow, J. Comput. Phys. 23, 263-275, 1977.
P. E. Dimotakis, H. J. Catrakis, Turbulence, fractals, and mixing,in Mixing: Chaos and Turbulence, NATO ASI Ser., Ser. B, Phys.,
Plenum, New York, 59-143, 1999.