plumes

Journal of Volcanology and Geothermal Research 221–222 (2012) 1–13

Contents lists available at SciVerse ScienceDirect

Journal of Volcanology and Geothermal Research journal homepage: www.elsevier.com/locate/jvolgeores

3-D numerical simulations of eruption column collapse: Effects of vent size on pressure-balanced jet/plumes Yujiro J. Suzuki ⁎, Takehiro Koyaguchi Earthquake Research Institute, University of Tokyo, Yayoi 1-1-1, Bunkyo-ku, Tokyo 113-0032, Japan

a r t i c l e

i n f o

Article history: Received 15 October 2011 Accepted 28 January 2012 Available online 11 February 2012 Keywords: Eruption cloud Column collapse condition Numerical simulation Flow-regime map Turbulent mixing

a b s t r a c t Buoyant columns or pyroclastic flows form during explosive volcanic eruptions. In the transition-state, these two eruption styles can develop simultaneously. We investigated the critical condition that separates the two eruption styles (referred to as “the column collapse condition”) by performing a series of three-dimensional numerical simulations. In the simulation results, we identify two types of flow regime: a turbulent jet that efficiently entrains ambient air (jet-type) and a fountain with a high-concentration of the ejected material (fountain-type). Hence, there are two types of column collapse (jet-type and fountain-type). Which type of collapse occurs at the column collapse condition depends on whether the critical mass discharge rate for column collapse (MDRCC) is larger or smaller than that for the generation of a fountain (MDRJF) for a given exit velocity. Temperature controls the relative magnitude of MDRCC relative to MDRJF, and hence the type of collapse. For given magma properties (e.g., temperature and water content), the column collapse condition is expressed by a critical value of the Richardson number (RiCC ≡ g′0L0/w02, where g′0 is the source buoyancy, L0 is the vent radius, and w0 is the exit velocity). When the jet-type collapse occurs at the column collapse condition, RiCC is independent of the exit velocity. When the fountain type collapse occurs at the column collapse condition, on the other hand, RiCC decreases as the exit velocity increases. As the exit velocity exceeds the sound velocity, a robust flow structure with a series of standing shock waves develops in the fountain, which suppresses entrainment of ambient air and enhances column collapse. © 2012 Elsevier B.V. All rights reserved.

1. Introduction During explosive volcanic eruptions, a mixture of hot ash (pyroclasts) and volcanic gas is released from the vent into the atmosphere. The mixture generally has an initial density several times larger than the atmospheric density at the vent and its ascent is driven only by its momentum. As the ejected material entrains ambient air, the density of the mixture decreases because the entrained air expands by heating from the hot pyroclasts. If the density of the mixture becomes less than the atmospheric density before the eruption cloud loses its upward momentum, a buoyant plume rises to form a Plinian eruption column. On the other hand, if the mixture loses its upward momentum before it becomes buoyant, the eruption column collapses to generate a pyroclastic flow. We refer to the condition that separates these two eruption regimes as “the column collapse condition” (e.g., Sparks and Wilson, 1976). Because the impact and type of volcanic hazards are largely different between the two eruption regimes, predicting the column collapse condition is of critical importance in physical volcanology.

⁎ Corresponding author. Tel.: + 81 3 5841 5680. E-mail addresses: [email protected] (Y.J. Suzuki), [email protected] (T. Koyaguchi). 0377-0273/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.jvolgeores.2012.01.013

The column collapse condition has previously been predicted by one dimensional (1-D) steady eruption column models (e.g., Wilson et al., 1980; Bursik and Woods, 1991; Kaminski and Jaupart, 2001; Carazzo et al., 2008; Koyaguchi et al., 2010). First generation 1-D steady models are based on a classical theory of turbulent jet/plumes from a point source, in which the efficiency of turbulent mixing between the eruption cloud and ambient air is expressed by a single parameter, the “entrainment coefficient” (e.g., Taylor, 1945; Morton et al., 1956). These models captured the basic physics of eruption column dynamics and enabled us to calculate the upward momentum and buoyancy of eruption columns as a function of the downstream distance from the vent, and thus, to predict the maximum heights of the eruption column and column collapse condition for given magma properties (e.g., temperature and water content) and source conditions (e.g., vent radius and velocity). The models of column collapse show a quasi-quantitative agreement with field observations of witnessed eruptions (e.g., Carey et al., 1990). Recently, the agreement between the model predictions and field observations has been improved by a more sophisticated 1-D model in which the variation of entrainment coefficient with height is taken into account (Carazzo et al., 2008). Since the 1980s, 2-D and 3-D models have been developed for eruption column dynamics. They have shown a number of complex fluid dynamical features that cannot be described by 1-D models

2

Y.J. Suzuki, T. Koyaguchi / Journal of Volcanology and Geothermal Research 221–222 (2012) 1–13

(e.g., Wohletz et al., 1984; Valentine and Wohletz, 1989; Neri and Dobran, 1994; Oberhuber et al., 1998; Suzuki et al., 2005; Ogden et al., 2008a,b; Suzuki and Koyaguchi, 2009). Generally, the flow of a gas–pyroclast mixture near the vent is characterized by an annular structure consisting of an outer shear region and dense inner core (referred to as a “potential core” in aerodynamics literature (e.g., Rajaratnam, 1976)). The presence of the dense inner core affects the flow pattern during column collapse (e.g., Suzuki et al., 2005). When the mixture is ejected from a large vent (more than a few hundred meters), the outer shear region cannot reach the central axis before the initial momentum is exhausted, so that the dense inner core generates a fountain structure, called “a radially suspended flow” (Neri and Dobran, 1994). On the other hand, for eruption from a relatively narrow vent, the dense inner core disperses due to turbulent mixing and the eruption cloud collapses without a fountain structure. Ogden et al. (2008a) pointed out that the overpressure at the vent generates another type of annular structure in the flow near the vent, which causes complex flow behavior during column collapse; an oscillatory collapse with a regular periodicity occurs even when a steady condition is applied at the source. In this paper and some accompanying papers, we systematically investigate how these 2-D and 3-D effects modify the column collapse condition on the basis of a number of 3-D numerical simulations of an eruption cloud with a high spatial resolution. We carried out an extensive parameter study to make regime maps of different flow patterns and to determine the column collapse condition for each flow pattern. From the features of the column collapse conditions, as well as those of the boundaries between the different flow patterns, we can infer the governing factors which control the dynamics of column collapse. In this particular paper, we focus on the flow whose pressure balances ambient atmospheric pressure (referred to as “a pressurebalanced jet/plume”). We show that even for this simplest case, the 2-D and 3-D effects of flow (e.g., the presence of a potential core with or without shock waves) substantially modify the column collapse condition.

2. Model description Our main concern is to clarify how the variation in the vent radius influences the global features of eruption clouds such as the formation of an eruption column and/or pyroclastic flow. The most essential physics which governs the dynamics of an eruption cloud is the efficiency of turbulent mixing between the ejected material and ambient air. In order to minimize complexity and to focus on the problem of turbulent mixing, a pseudo-gas model is applied in this study; we ignore the separation of solid pyroclasts from the eruption cloud and treat the eruption cloud as a single gas whose density is calculated from the mixing ratio of the ejected material and entrained air. The fluid dynamics model solves a set of partial differential equations describing the conservation of mass, momentum, and energy, and constitutive equations describing the thermodynamic state of the mixture of pyroclasts, volcanic gas, and air. These equations are solved numerically by a general scheme for compressible flow (Roe, 1981; van Leer, 1977). Details of the numerical procedures used in this study are described in Suzuki et al. (2005). The simulations are designed to describe the injection of a mixture of pyroclasts and volcanic gas from a circular vent above a flat surface of the earth in a stationary atmosphere with a temperature gradient typical of the mid-latitude atmosphere (Table 1). The vent is located in the center of the ground surface. The physical domain extends horizontally and vertically from a few kilometers to several tens of kilometers. At the ground boundary, the free-slip condition is assumed for the velocities of the ejected material and air. At the upper and other boundaries of the computational domain, the fluxes of mass, momentum, and energy are assumed to be continuous, and these

Table 1 List of material properties and values of physical parameters. Variable

Value (unit)

Meaning

g Rg Ra Cvs Cvg Cva Ta0 pa0 ρa0 μ1 μ2

9.81 (m s− 2) 462 (J kg− 1 K− 1) 287 (J kg− 1 K− 1) 1100 (J kg− 1 K− 1) 1348 (J kg− 1 K− 1) 713 (J kg− 1 K− 1) 273 (K) 1.013 × 105 (Pa) 1.29 (kg m− 3) − 6.5 (K km-1) 0.0 (K km− 1)

Gravitational body force per unit mass Gas constant of volcanic gas Gas constant of atmospheric air Specific heat of solid pyroclasts Specific heat of volcanic gas at constant volume Specific heat of air at constant volume Atmospheric temperature at z = 0 km Atmospheric pressure at z = 0 km Atmospheric density at z = 0 km Temperature gradient at z = 0.0–11.0 km Temperature gradient at z = 11.0–20.0 km

boundary conditions correspond to the free outflow and inflow of these quantities. In order to correctly reproduce the general feature of turbulent mixing that the efficiency of entrainment is independent of the Reynolds number (Dimotakis, 2000), it is essential to apply 3-D coordinates with a sufficiently high spatial resolution (Suzuki et al., 2005). In this study, the calculations were performed on a 3-D domain with a non-uniform grid. The grid size was set to be sufficiently smaller than L0/8 near the vent, where L0 is the vent radius. In order to effectively 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.01 for the vertical coordinate and 1.02 for the horizontal coordinate) up to L0 with the distance from the vent, such that the grid size is small enough to resolve the turbulent flow far from the vent (cf. Suzuki and Koyaguchi, 2010). We assume steady conditions; for each run, magma temperature T0, water content ng0, vent radius L0, and exit velocity w0 are fixed. _ 0 is determined from the relationship of The mass discharge rate m _ 0 ¼ πρ0 w0 L20 , where ρ0 is a function of T0 and ng0. In this paper, for m simplicity we assumed that the pressure at the vent is equal to the atmospheric pressure. In general, the pressure of a gas–pyroclast mixture deviates from the atmospheric pressure during explosive eruptions; the mixture is accelerated and/or decelerated because of decompression and/or compression within a short distance from the vent (5–20 vent radius) (Woods and Bower, 1995; Ogden et al., 2008b; Koyaguchi et al., 2010). Another extensive parameter study focusing on this effect is in progress and the results will be reported elsewhere. In order to classify different flow patterns and determine the column collapse condition in pressure-balanced jet/plumes, we performed a parameter study involving 95 numerical simulations. The conditions for these simulations were divided into four groups according to magma temperature (T0) and water content (ng0); group H with a high magma temperature (T0 = 1000 K and ng0 = 2.84 wt.%), group I with an intermediate magma temperature (T0 = 800 K and ng0 = 2.84 wt.%), group L with a low magma temperature (T0 = 550 K and ng0 = 2.84 wt.%), and group N with a small water content (T0 = 1000 K and ng0 = 1.23 wt.%). These magma properties are considered to represent those of typical explosive eruptions (see Carazzo et al., 2008; Koyaguchi et al., 2010). Low temperature mixture of group L can be ejected during eruptions involving the interaction of ground water and/or country rock with magma (e.g., Koyaguchi and Woods, 1996). In order to cover the conditions of typical explosive eruptions, we set the mass discharge rate to be 10 4– 10 9 kg s − 1 and the exit velocity to be 0.5c0–3c0, where c0 is the sound of the gas–pyroclast mixture at the vent, defined as ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi pvelocity c0 ¼ 0:01ng0 Rg T 0 . Generally, exit velocity, magma properties (i.e., T0 and ng0), vent radius, and hence, magma discharge rate are coupled such that they are consistent with the fluid dynamics in the conduit and crater (Wilson et al., 1980; Koyaguchi et al., 2010). In this study, we assessed whether or not the above ranges are physically

Y.J. Suzuki, T. Koyaguchi / Journal of Volcanology and Geothermal Research 221–222 (2012) 1–13 Table 2 (continued)

Table 2 Input parameters and flow regimes. Simulation

_0 m (kg s− 1)

w0 (m s− 1)

Ri

Group H (T0 = 1000 K, ng0 = 2.84 wt.%) Run 24 1.0 × 109 115 Run 25 1.0 × 109 230 9 Run 26 1.0 × 10 173 9 Run 27 1.0 × 10 346 9 Run 28 1.0 × 10 288 230 Run 29 1.0 × 108 Run 30 1.0 × 108 346 Run 31 1.0 × 108 115 Run 32 1.0 × 108 288 8 Run 33 1.0 × 10 173 7 Run 34 1.0 × 10 115 7 230 Run 35 1.0 × 10 Run 36 1.0 × 107 346 Run 37 1.0 × 107 57.6 Run 39 1.0 × 106 57.6 Run 40 1.0 × 106 115 6 Run 41 3.2 × 10 57.6 Run 42 3.2 × 106 115 6 173 Run 43 3.2 × 10 Run 44 1.0 × 106 80.7 Run 45 3.2 × 106 80.7 Run 46 3.2 × 107 196 Run 47 3.2 × 107 253 Run 48 1.0 × 108 150 9 Run 49 1.0 × 10 254 6 Run 50 1.0 × 10 69.1 Run 51 1.0 × 107 80.6 Run 52 1.0 × 107 92.2 Run 53 1.0 × 108 138 Run 54 3.2 × 108 173 Run 55 3.2 × 108 196 8 Run 56 3.2 × 10 219 8 Run 57 3.2 × 10 242 Run 67 3.2 × 106 69.1 Run 68 1.0 × 106 50.0 Run 69 1.0 × 106 34.6 Run 99 1.0 × 107 173 8 Run 101 1.3 × 10 288

598 423 488 345 378 134 109 189 120 154 59.8 42.3 34.5 84.6 26.8 18.9 47.6 33.6 27.5 22.6 40.2 81.6 71.7 166 403 24.4 71.5 66.9 173 275 258 244 232 43.4 28.7 34.5 49.0 136

2.2 × 100 3.9 × 10− 1 8.0 × 10− 1 1.4 × 10− 1 2.2 × 10− 1 1.2 × 10− 1 4.5 × 10− 2 7.0 × 10− 1 7.0 × 10− 2 2.5 × 10− 1 2.2 × 10− 1 3.9 × 10− 2 1.4 × 10− 2 1.2 × 100 3.9 × 10− 1 7.0 × 10− 2 7.0 × 10− 1 1.2 × 10− 1 4.5 × 10− 2 1.7 × 10− 1 3.0 × 10− 1 1.0 × 10− 1 5.5 × 10− 2 3.6 × 10− 1 3.1 × 10− 1 2.5 × 10− 1 5.4 × 10− 1 3.8 × 10− 1 4.4 × 10− 1 4.5 × 10− 1 3.3 × 10− 1 2.5 × 10− 1 1.9 × 10− 1 4.4 × 10− 1 5.6 × 10− 1 1.4 × 100 7.5 × 10− 2 7.5 × 10− 2

1.0 2.0 1.5 3.0 2.5 2.0 3.0 1.0 2.5 1.5 1.0 2.0 3.0 0.5 0.5 1.0 0.5 1.0 1.5 0.7 0.7 1.7 2.2 1.3 2.2 0.6 0.7 0.8 1.2 1.5 1.7 1.9 2.1 0.6 0.4 0.3 1.5 2.5

F-CC F-CC F-CC F-EC F-CC F-EC J-EC F-CC F-EC F-EC F-EC J-EC J-EC F-CC J-EC/F-EC J-EC F-CC J-EC J-EC J-EC J-EC/F-EC J-EC/F-EC J-EC F-EC F-CC J-EC F-CC F-EC F-CC F-CC F-CC F-EC F-EC F-EC F-EC F-CC J-EC F-EC

Group I (T0 = 800 K, ng0 = 2.84 wt.%) Run 120 1.0 × 106 103 Run 121 4.0 × 106 103 Run 122 1.0 × 107 103 Run 123 2.5 × 107 103 Run 124 4.0 × 106 206 7 Run 125 1.0 × 10 206 7 Run 126 2.5 × 10 206 Run 127 4.0 × 105 103 Run 128 6.3 × 107 206 Run 129 1.6 × 108 206 Run 130 1.6 × 106 103 Run 131 2.5 × 106 103

17.9 35.7 56.6 89.7 25.2 40.0 63.4 11.3 101 159 22.5 28.4

1.1 × 10− 1 2.1 × 10− 1 3.4 × 10− 1 5.4 × 10− 1 3.8 × 10− 2 6.0 × 10− 2 9.5 × 10− 2 6.8 × 10− 2 1.5 × 10− 1 2.4 × 10− 1 1.3 × 10− 1 1.7 × 10− 1

1.0 1.0 1.0 1.0 2.0 2.0 2.0 1.0 2.0 2.0 1.0 1.0

J-EC/J-CC F-EC/F-CC F-CC F-CC J-EC J-EC J-CC/F-CC J-EC F-CC F-CC J-CC J-CC/F-CC

7.8 × 10− 2 5.7 × 10− 2 3.8 × 10− 2 6.9 × 10− 2 5.7 × 10− 2 1.2 × 10− 1 2.5 × 10− 1 6.8 × 10− 1 2.2 × 10− 2 4.7 × 10− 2 2.7 × 10− 2 2.4 × 10− 2 3.0 × 10− 2 3.4 × 10− 2 3.5 × 10− 2 2.9 × 10− 2 2.2 × 10− 1 1.2 × 100 5.3 × 10− 1

1.5 1.7 2.0 2.5 2.7 2.0 1.5 1.0 2.5 2.9 2.3 2.4 2.2 2.1 1.3 1.4 1.0 0.5 0.7

F-EC J-CC J-CC J-CC J-CC J-CC F-CC F-CC J-EC J-CC J-EC/J-CC J-EC J-CC J-CC J-CC J-CC J-CC/F-CC F-CC F-CC

13.3 12.5 11.5 32.6 31.4 36.4 42.1 51.5 10.3 30.2 10.7 10.5 11.0 11.2 4.52 4.35 16.3 23.0 19.5

M

Regime

Simulation

_0 m (kg s− 1)

w0 (m s− 1)

Run Run Run Run Run Run Run Run Run Run Run

1.0 × 106 1.0 × 105 1.0 × 105 1.0 × 107 1.0 × 107 1.0 × 107 1.0 × 107 1.0 × 107 1.0 × 107 1.0 × 107 1.0 × 107

68.4 128 137 145 154 162 137 94.0 103 111 120

a

L0 (m)

Group L (T0 = 550 K, ng0 = 2.84 wt.%) Run 58 1.0 × 106 128 Run 59 1.0 × 106 145 Run 60 1.0 × 106 171 Run 61 1.0 × 107 214 Run 62 1.0 × 107 231 7 Run 64 1.0 × 10 171 7 Run 65 1.0 × 10 128 7 Run 66 1.0 × 10 85.5 Run 70 1.0 × 106 214 Run 71 1.0 × 107 248 Run 74 1.0 × 106 197 Run 75 1.0 × 106 205 6 Run 77 1.0 × 10 188 6 Run 78 1.0 × 10 179 5 Run 80 1.0 × 10 111 Run 81 1.0 × 105 120 Run 85 1.0 × 106 85.5 Run 86 1.0 × 106 42.7 Run 87 1.0 × 106 59.8

3

(continued on next page)

88 89 90 91 92 93 94 95 96 97 98

Group N (T0 = 1000 K, ng0 = 1.23 wt.%) 75.6 Run 103 1.0 × 107 Run 104 1.0 × 108 75.6 Run 105 1.0 × 109 75.6 Run 106 1.0 × 107 151 8 Run 107 1.0 × 10 151 Run 108 1.0 × 109 151 6 Run 109 4.0 × 10 75.6 151 Run 110 4.0 × 106 Run 111 4.0 × 105 75.6 Run 112 4.0 × 105 151 Run 115 1.0 × 106 75.6 Run 116 1.6 × 106 75.6 6 Run 117 2.5 × 10 75.6 7 Run 118 1.6 × 10 151 Run 119 2.5 × 107 151

Ri

M

Regimea

18.2 4.21 4.07 39.5 38.4 37.4 40.7 49.1 47.0 45.2 43.5

3.8 × 10− 1 2.5 × 10− 2 2.1 × 10− 2 1.8 × 10− 1 1.6 × 10− 1 1.4 × 10− 1 2.1 × 10− 1 5.4 × 10− 1 4.3 × 10− 1 3.5 × 10− 1 2.9 × 10− 1

0.8 1.5 1.6 1.7 1.8 1.9 1.6 1.1 1.2 1.3 1.4

F-CC J-EC J-EC J-CC/F-CC J-CC/F-CC F-CC/F-CC F-CC F-CC F-CC F-CC F-CC

48.6 154 486 34.4 109 344 30.7 21.7 9.70 6.86 15.4 19.4 24.4 43.3 54.5

1.1 × 100 3.4 × 100 1.1 × 101 1.9 × 10− 1 6.0 × 10− 1 1.9 × 101 6.8 × 10− 1 1.2 × 10− 1 2.1 × 10− 1 3.8 × 10− 2 3.4 × 10− 1 4.3 × 10− 1 5.4 × 10− 1 2.4 × 10− 1 3.0 × 10− 1

1.0 1.0 1.0 2.0 2.0 2.0 1.0 2.0 1.0 2.0 1.0 1.0 1.0 2.0 2.0

F-CC F-CC F-CC F-EC F-CC F-CC F-CC J-EC J-EC J-EC J-EC F-EC/F-CC F-CC F-CC F-CC

L0 (m)

J-EC: jet-type eruption column, J-CC: jet-tyle column collapse. a F-EC: fountain-type eruption column, F-CC: fountain-type column collapse.

realistic on the basis of Koyaguchi et al. (2010). The wide range of exit velocity allows us to investigate features of subsonic, sonic, and supersonic flows. A summary of the simulation parameters for the four groups is listed in Table 2. 3. Results 3.1. Classification of flow patterns The flow patterns in the simulation results are classified into four flow regimes: an eruption column with a jet structure (e.g., run 99), an eruption column with a fountain structure (e.g., run 33), jet collapse (e.g., run 78), and fountain collapse (e.g., run 49). The representative features of each flow regime are described below. 3.1.1. Eruption column regimes _ 0 ¼ 1:0  107 kg s − 1, L0 = 49 m), a stable In run 99 (T0 = 1000 K, m eruption column with a typical jet structure develops (Fig. 1). In this run, the mixture of volcanic gas and pyroclasts is ejected from the vent as a dense, high speed jet. After traveling a short distance from the vent, the shear flow at the boundary between the jet and the ambient air becomes unstable. The jet entrains ambient air by this shear instability; it forms an annular mixing layer which surrounds an unmixed core flow (Fig. 1a). The unmixed core is eroded by the annular mixing layer and disappears at a certain level (at a height of ~ 3 km (60 times L0) in run 99). As the eruption column further ascends, the flow becomes highly unstable and undergoes a meandering instability that induces more efficient mixing (e.g., Suzuki et al., 2005). The stream-wise evolution of the flow results in a complex density distribution in the eruption cloud (Fig. 1b). The unmixed core is denser than the ambient air, whereas the annular mixing layer is lighter than the ambient air and gains buoyancy owing to expansion of the entrained air. After the meandering instability develops, the dense unmixed core disappears before it exhausts its initial momentum. We refer to the eruption column that is characterized by the

3

3

z [km]

Y.J. Suzuki, T. Koyaguchi / Journal of Volcanology and Geothermal Research 221–222 (2012) 1–13

z [km]

4

(b)

a

0

0

(a) 3

0

3

3

0

r [km]

0

0.5

3

r [km]

1

-0 .5

0

0 .5

Fig. 1. Representative numerical results of the eruption column with a jet structure (run 99). Parameters used and conditions at the vent are listed in Tables 1 and 2, respectively. (a) Crosssectional distribution of the mass fraction of the ejected material (ξ) in the r–z space at 600 s. (b) Cross-sectional distribution of the density difference relative to the stratified atmospheric density at the same vertical position (Δρ ¼ ρ=ρa −1) in the r–z space at 600 s.

z [km]

3

z [km]

3

column (run 33). In the jet-type column, the mass fraction of the ejected material, ξ, gradually decreases with height (Fig. 3a). The high-ξ region intermittently rises at a time interval of ~ 10 s; the flow of the unmixed core is highly unsteady. In the fountain-type column, on the other hand, ξ remains at 1.0 below a height of 2.2 km (Fig. 3b) and decreases rapidly at this level, which corresponds to that of the large-scale eddy at the top of the fountain. The height of the rapid change of ξ is constant with time and does not show oscillatory behavior, which suggests that the flow of the unmixed core is almost steady. In Fig. 4, the time-averaged mass fraction of the ejected material, ξ, and normalized upward velocity, w/w0, along the central axis are presented as functions of the normalized vertical position, z*∗ (≡z/(w02/2g)) for the jet-type column (run 99) and fountain-type column (run 33); normalized quantities are used for quantitative comparison between results of different vent conditions in these diagrams. In the jet-type column, ξ gradually decreases with height from 1.0 to 0.8 below z* = 1 and rapidly falls to 0.15 from z* = 1 to 2 (solid curve in Fig. 4a). The averaged upward velocity, w/w0, gradually decreases from 1.0 to 0.5 below z* = 1 and slightly increases above it (solid curve in Fig. 4b).

disappearance of the unmixed core before it exhausts its initial momentum as “the jet-type column”. When the gas–pyroclast mixture is ejected from a larger vent (run 33; L0 = 154 m), a stable eruption column with a fountain structure develops (Fig. 2). In this run, because the vent radius is relatively large, the erosion by the annular mixing layer does not reach the central axis at the height where the unmixed core exhausts its initial momentum (Fig. 2a). The dense unmixed core around the central axis spreads radially at this height (Fig. 2b). Such a radially suspended flow is commonly observed in the fountain which results from the injection of a dense fluid upward into a less dense fluid (e.g., Neri and Dobran, 1994; Lin and Armfield, 2000). The large-scale eddy of this radially suspended flow causes an intensive mixing between the ejected material and ambient air. After this mixing, the resultant mixture becomes buoyant and generates an upward flow from the largescale eddy; this produces another type of stable column. The eruption column that is characterized by the radially suspended flow is referred to as “the fountain-type column”. Fig. 3 shows the time evolutions of the vertical profiles along the central axis for the jet-type column (run 99) and fountain-type

3

0

3

3

r [km]

0

a

0

(b)

0

(a)

0.5

0

3

r [km]

1

-0 .5

0

0 .5

Fig. 2. Representative numerical results of the eruption column with a fountain structure (run 33). Parameters used and conditions at the vent are listed in Tables 1 and 2, respectively. (a) Cross-sectional distribution of the mass fraction of the ejected material (ξ) in the r–z space at 190 s. (b) Cross-sectional distribution of the density difference relative to the stratified atmospheric density at the same vertical position (Δρ ¼ ρ=ρa −1) in the r–z space at 190 s.

5

1

10

Y.J. Suzuki, T. Koyaguchi / Journal of Volcanology and Geothermal Research 221–222 (2012) 1–13

0.5 0

5 0

z [km]

(a) run 99

0

200

400

600

200

300

1

10

Time [sec]

0.5 0

5 0

z [km]

(b) run 33

0

100

Time [sec] Fig. 3. Timewise distributions of the mass fraction of the ejected material (ξ) along the central axis (spatial average within the range of |r| b L0) for (a) run 99 and (b) run 33.

the fountain changes with the exit velocity (Fig. 5). When the exit velocity is sonic (the Mach number is 1.0 in run 34), the time-averaged pressure inside the fountain is uniform and almost equal to the atmospheric pressure (Fig. 5a). When the exit velocity is supersonic (the Mach number is 1.5 in run 33), a spatially periodic pattern of higher and lower pressure is observed in the fountain (Fig. 5b). This periodic pattern results from the generation of standing shock waves in the fountain; like in an overexpanded jet (i.e., a jet with exit pressure lower than the surrounding fluid), compression waves (i.e., oblique shock waves) are produced in the supersonic jet, and these compression waves reflect from the symmetry axis and generate expansion waves (Dobran et al., 1993).

The height of z* = 1 coincides with the level where the unmixed core disappears (z = 3 km). These profiles of ξ and w/w0 are explained by the gradual growth of the annular mixing layer and the erosion of the unmixed core below z* = 1. In the fountain-type column, on the other hand, ξ remains 1.0 below z* = 0.7, rapidly decreases to 0.1 at the fountain top (z* = 0.7), and remains 0.1 above it (dashed curve in Fig. 4a). The averaged upward velocity, w/w0, decreases from 1.0 to − 0.1 in the fountain below z* = 0.7, and rapidly increases above it (dashed curve in Fig. 4b). These profiles of ξ and w/w0 reflect the features of entrainment in the fountain-type column: extensive entrainment occurs at the fountain top, whereas the entrainment is inefficient in the fountain. It is noted that w/w0 has a negative value just above the fountain top (at z* = 0.7–0.9). This negative value of w/w0 is caused by vortical motion associated with the large-scale eddy at the fountain top. The extensive entrainment by the large-scale eddy at the fountain top provides new insight into the mixing process in eruption columns, because the mixing at the vertical edge of column has been assumed in the 1-D models. In the simulation results of the fountain-type column regimes, we observed that the qualitative feature of pressure distribution inside

3.1.2. Column collapse regimes In the present simulations, column collapse occurs when the vent radius is large. For given vent radius and exit velocity, column collapse tends to occur as magma temperature decreases. The flow pattern of the column collapse regimes depends on the magma temperature and vent radius.

2

z*

z*

2

1

(a)

(b)

0 0

0.2

1

0.4

0.6

0.8

1

0 -0.2

0

0.2

0.4

0.6

0.8

1

1.2

w/w0 Fig. 4. Time-averaged vertical profiles of physical quantities along the central axis (spatial average within the range of |r| b L0) for run 99 (solid curves) and run 33 (dashed curves). (a) The mass fraction of the ejected material, ξ. (b) The upward velocity normalized by the exit velocity at the vent, w/w0. The vertical axis represents the vertical position normalized by the characteristic height (z* = z/(w02/2g)).

2

2

z [km]

Y.J. Suzuki, T. Koyaguchi / Journal of Volcanology and Geothermal Research 221–222 (2012) 1–13

z [km]

6

(b) run 33

0

0

(a) run 34 2

0

2

2

0

r [km]

-0 .2

0

2

r [km]

0 .2

-0 .2

0

0 .2

Fig. 5. Time-averaged distribution of the pressure difference relative to the stratified atmospheric pressure at the same vertical position, Δρ ¼ ρ=ρa −1, where ρa is the atmospheric pressure. Cross-sectional distributions for (a) run 34 and (b) run 33. Parameters used and conditions at the vent are listed in Tables 1 and 2, respectively.

referred to as “fountain-type collapse”. When the dense eruption cloud hits the ground and spreads horizontally, it efficiently entrains ambient air. The entrained air expands by heating from the pyroclasts so that the mixture of the ejected material and ambient air forms a buoyant co-ignimbrite ash cloud. Fig. 8 shows the time-averaged vertical profiles of physical quantities along the central axis for the jet-type collapse (run 78) and fountain-type collapse (run 49). In the jet-type collapse, ξ and w/w0 gradually decrease with height and fall to zero at z* ~ 1 (solid curves in Fig. 8). The gradual decreases in ξ and w/w0 represent the dispersion of the unmixed core due to the gradual growth of the annular mixing layer. In the fountain-type collapse (dashed curve in Fig. 8a), ξ remains 1.0 in the fountain (z* b 1.2) and rapidly decreases to 0.1 at the fountain top (z* = 1.2). The value for z* > 1.2 (ξ = 0.1) represents the mass fraction of the ejected material in the co-ignimbrite ash cloud. The upward velocity decreases with height in the fountain and falls to zero at the fountain top (dashed curve in Fig. 8b). In contrast to the case of the jet-type collapse, the unmixed core itself is the main source of pyroclastic flows in the case of the fountain-type collapse.

0.5

3.1.3. Transitional regimes The boundaries of the above four flow regimes (the jet-type column, fountain-type column, jet-type collapse, and fountain-type

z [km]

0.5

z [km]

When a low-temperature gas–pyroclast mixture is ejected from a narrow vent (run 78; T0 = 550 K, L0 = 11 m), the eruption column collapses from the jet (Fig. 6). In this run, the jet entrains ambient air by the shear instability; the annular mixing layer develops and the unmixed core disappears at a height of 0.5 km (Fig. 6a). Because the initial temperature is low, the eruption cloud remains heavier than air even though it entrains a large amount of air. As a result, the low-concentration mixture generates a pyroclastic flow (Fig. 6b). The column collapse regime that is characterized by the disappearance of the unmixed core before it exhausts its initial momentum is referred to as “the jet-type collapse”. In contrast to the jet-type collapse, when a high-temperature gas– pyroclast mixture is ejected from an extremely large vent (run 49; T0 = 1000 K, L0 = 403 m), the eruption column collapses from a fountain (Fig. 7a). As was observed in the case of the fountain-type column (Fig. 2), when the vent radius is large, the unmixed core reaches a height where the initial momentum is exhausted and generates a radially suspended flow which causes mixing by a large-scale eddy. As the vent radius exceeds a certain value, the large-scale eddy of the radially suspended flow cannot entrain sufficient air for the whole mixture to become buoyant (Fig. 7b). As a result, the unmixed part in the radially suspended flow collapses and generates a pyroclastic flow. The column collapse regime that is characterized by the radially suspended flow is

0.5

0

0.5

0.5

r [km]

0

a

0

(b)

0

(a)

0.5

0

0.5

r [km]

1

-0 .5

0

0 .5

Fig. 6. Representative numerical results of the column collapse from the jet structure (run 78). Parameters used and conditions at the vent are listed in Tables 1 and 2, respectively. (a) Cross-sectional distribution of the mass fraction of the ejected material (ξ) in the r–z space at 120 s. (b) Cross-sectional distribution of the density difference relative to the stratified atmospheric density at the same vertical position (Δρ ¼ ρ=ρa −1) in the r–z space at 120 s.

(a)

(b)

a

0

0

7

10

z [km]

10

z [km]

Y.J. Suzuki, T. Koyaguchi / Journal of Volcanology and Geothermal Research 221–222 (2012) 1–13

0

10

10

0

10

r [km]

0

10

r [km]

0.5

1

-0 .5

0

0 .5

Fig. 7. Representative numerical results of the column collapse from the fountain structure (run 49). Parameters used and conditions at the vent are listed in Tables 1 and 2, respectively. (a) Cross-sectional distribution of the mass fraction of the ejected material (ξ) in the r–z space at 500 s. (b) Cross-sectional distribution of the density difference relative to the stratified atmospheric density at the same vertical position (Δρ ¼ ρ=ρa −1) in the r–z space at 500 s.

flow. This pyroclastic flow is less energetic than the pyroclastic flow of the fountain-type collapse regime which falls down directly from the fountain top (cf. Fig. 7). Around the CC conditions for groups L and I, a jet-type column becomes unstable with time and can collapse to generate a dilute pyroclastic flow; for example, in run 120 (T0 = 800 K, L0 = 18 m), the transition from the jet-type column to the jet-type collapse regimes occurs at t ~ 150 s from the beginning of the eruption (the result is not shown).

collapse regimes) are determined by two types of transition. One is the transition between the jet-type and fountain-type regimes; we call the condition for this transition the jet-fountain (JF) condition. The other is the transition between the column convection and collapse regimes; we call the condition for this transition the column collapse (CC) condition. Around the JF condition, some intermediate features between jettype and fountain-type regimes were observed. Run 46 (T0 = 1000 K, L0 = 82 m) shows an intermediate feature between the jet-type column and the fountain-type column regimes (Fig. 9). In this run, the unmixed core is eroded to a considerable extent by the annular mixing layer, while relatively large-scale eddies develop at the fountain top (z ~ 4 km). In run 92 (T0 = 550 K, L0 = 38 m), the type of column collapse changes from the jet-type to the fountain-type with time (Fig. 10). At the initial stage of eruption (t b 50 s), the unmixed core disappears at z = 0.6 km (see the contour of ξ = 0.95), whereas the total height of the negatively buoyant jet reaches 1.4 km high. Subsequently, the height of the unmixed core increases with time (50 s b t b 80 s) and reaches the top of the jet so that the fountain structure (i.e., the radially suspended flow) develops (t > 80 s). Around the CC condition, some intermediate features between column and collapse regimes are observed. In run 27 (T0 = 1000 K, L0 = 345 m), a fountain-type column and pyroclastic flow develop simultaneously (Fig. 11). After mixing due to the large-scale eddy at the fountain top, most of the cloud becomes a buoyant plume, whereas the side of the fountain partially flows down to form a pyroclastic

3.2. Flow-regime map On the basis of extensive parameter studies, we made flow-regime _ 0 ) and exit maps in the parameter space of the mass discharge rate (m velocity (w0) for 550 K b T0 b 1000 K and 1.23 wt.% b ng0 b 2.84 wt.% (Fig. 12). When the temperature is low (group L) or intermediate (group I), the flow changes from the jet-type column to the jet-type collapse regimes, and from the jet-type collapse to the fountain-type collapse regimes as the mass discharge rate increases for a given exit velocity (Fig. 12a,b). The boundary between the jet-type column and the jettype collapse regimes represents the CC condition (solid curves), and that between the jet-type collapse and the fountain-type collapse regimes represents the JF condition (dashed curves). Both the critical mass discharge rate for the CC condition (MDRCC) and for the JF condition (MDRJF) increase as the exit velocity increases. For a low temperature, MDRCC is substantially smaller than MDRJF; as a result, the

2

z*

z*

2

1

1

(a)

(b)

0 0

0.2

0.4

0.6

0.8

1

0 -0.2

0

0.2

0.4

0.6

0.8

1

w/w0

Fig. 8. Time-averaged vertical profiles of physical quantities along the central axis (spatial average within the range of |r| b L0) for run 78 (solid curves) and run 49 (dashed curves). (a) The mass fraction of the ejected material, ξ. (b) The upward velocity normalized by the exit velocity at the vent, w/w0. The vertical axis represents the vertical position normalized by the characteristic height (z* = z/(w02/2g)).

8

Y.J. Suzuki, T. Koyaguchi / Journal of Volcanology and Geothermal Research 221–222 (2012) 1–13

4

0.9

0.8

Ri ¼

g′0 L0 ; w20

ð1Þ

where the source buoyancy, g′0, is defined using the density of the ejected material, ρ0, and the atmospheric density, ρa, as g′0 = (ρ0 − ρa) g/ρ0. The Mach number is defined as the ratio of the exit velocity and the sound velocity of the ejected material, c0:

0

z [km]

In order to understand the physical meanings of the two critical conditions, we redraw Fig. 12 in a non-dimensional space, namely in the space of the Richardson number (Ri) and the Mach number (M) (Fig. 13). The Richardson number is defined by

4

0

4

r [km]

0

0.5

M¼

1

Fig. 9. Cross-sectional distribution of the mass fraction of the ejected material (ξ) in the r–z space for the result of run 46 at 190 s. Parameters used and conditions at the vent are listed in Tables 1 and 2, respectively. The contour levels are ξ = 0.8 and 0.9.

jet-type collapse occurs in a wide range of parameters in this diagram (see Fig. 12a). As the temperature increases, the distance between these two boundaries decreases so that the region of the jet-type collapse regime shrinks (Fig. 12b). When the temperature is high (group H), MDRCC becomes greater than MDRJF for a given exit velocity; as a result, the flow changes from the jet-type column to the fountain-type column regimes, and from the fountain-type column to the fountain-type collapse regimes as the mass discharge rate increases (Fig. 12c,d). In these cases, the CC condition is represented by the boundary between the fountaintype column and the fountain-type collapse regimes, and the JF condition is represented by the boundary between the jet-type column and the fountain-type column regimes. The results of group N (T0 = 1000 K and ng0 = 1.23 wt.%) shows that the distance between these two boundaries decreases as the water content decreases for a given temperature (Fig. 12d).

4. Physical meanings of the two critical conditions

ð2Þ

Fig. 13 indicates that both the two critical conditions are determined mainly by the Richardson number for given magma properties (cf. Valentine and Wohletz, 1989). We define the Richardson number at the CC condition as RiCC and that at the JF condition as RiJF, respectively. For a given Mach number, RiCC increases as T0 and/or ng0 increases (see gray stars in Fig. 13), whereas RiJF is almost independent of T0 and/or ng0 (see white stars in Fig. 13). Whether RiCC b RiJF or RiCC > RiJF determines the type of the intermediate flow regime. In the low-temperature case (Figs. 13a,b), RiCC b RiJF so that the intermediate flow regime is the jettype collapse. In the high-temperature case (Fig. 13c,d), RiCC > RiJF so that the intermediate flow regime is the fountain-type column. The physical mechanisms that cause the dependencies of RiJF and RiCC on magma properties will be discussed below. 4.1. Physical mechanisms to control RiJF The JF condition is determined by whether the unmixed core is eroded by the annular mixing layer or not before the initial momentum is exhausted. This problem has two length scales: the height where the unmixed core disappears (Hcore), and the height where the initial momentum is exhausted (Hmmt). When Hcore > Hmmt, the fountain is generated. According to experimental studies on highspeed jets (e.g., Nagamatsu et al., 1969), Hcore is a function of the vent radius, L0, and M: H core ¼ ða1 M þ a2 ÞL0 ;

ð3Þ

where a1 and a2 are positive constants. On the other hand, Hmmt is estimated from the balance between the kinetic and potential energy: 2

H mmt ¼

w0 : 2g0 0

ð4Þ

1

2

The variation in flow-regime maps suggests that which flow regime can occur is determined by the positional relationship between the JF and CC conditions. This positional relationship depends on the magma properties; for a given exit velocity, MDRrmCC is smaller than MDRJF when the temperature is relatively low, whereas MDRCC is larger than MDRJF when the temperature is high. In this section, we discuss how each critical condition depends on the magma properties.

w0 : c0

0.5

1

0

0.95

0

z [km]

0.9

0

50

100

150

Time [sec] Fig. 10. Timewise distributions of the mass fraction of the ejected material (ξ) along the central axis (spatial average within the range of |r| b L0) for run 92. Parameters used and conditions at the vent are listed in Tables 1 and 2, respectively. The contour levels are ξ = 0.9 and 0.95.

Y.J. Suzuki, T. Koyaguchi / Journal of Volcanology and Geothermal Research 221–222 (2012) 1–13

9

These two equations yield

10

z [km]

Hcore 2ða1 M þ a2 Þg′0 L0 ¼ ¼ 2ða1 M þ a2 ÞRi: H mmt w20

ð5Þ

Consequently, we obtain the JF condition (i.e., Hcore = Hmmt) as

0

RiJF ¼

10

0

10

r [km]

0

0.5

1

Fig. 11. Cross-sectional distribution of the mass fraction of the ejected material (ξ) in the r–z space for the result of run 27 at 590 s. Parameters used and conditions at the vent are listed in Tables 1 and 2, respectively.

Jet-type Column

Exit velocity [m/s]

Exit velocity [m/s]

300

Jet-type Collapse

100

(a) group L 0

200

106

107

0 5 10

108

Jet-type Column Fountain-type Collapse

(c) group H

(T0 =1000 K n0 =0.0248) 106

107

108

109

1010 1011

Mass discharge rate [kg/s]

Exit velocity [m/s]

Exit velocity [m/s]

300

106

200

Fountain-type Column

0 5 10

(b) group I (T0=800 K n0 =0.0248)

107

108

109

Mass discharge rate [kg/s]

400

100

Fountain-type Collapse

Jet-type Collapse

Fountain-type (T0 =550 K n0 =0.0248) Collapse

105

Jet-type Column

100

Mass discharge rate [kg/s]

200

ð6Þ

This equation successfully explains the simulation results in Fig. 13 that RiJF is independent of the magma properties for a fixed M, and that RiJF decreases as M increases. The above result that the transition between the jet and fountain occurs at a certain critical Ri is supported by laboratory experiments of negatively buoyant jet/fountains (e.g., Kaye and Hunt, 2006). In these experiments, when Ri at the source is smaller than a certain value of O(1), the flow behaves as a highly forced jet where the injected material efficiently mixes with the ambient fluid (i.e., the jet-type). When Ri is larger than the critical value, on the other hand, the flow becomes a fountain where there is little mixing between the injected material and the ambient fluid (i.e., the fountain-type). In order to extend the conclusion derived from

300

200

1 : 2ða1 M þ a2 Þ

Fountain-type Column

Jet-type Column

100

Fountain-type Collapse

(d) group N (T0 =1000 K n0=0.0123) 0 105

106

107

108

109

Mass discharge rate [kg/s]

Fig. 12. Flow-regime maps for the results of (a) group L, (b) group I, (c) group H, and (d) group N. Parameters used are listed in Table 1. The conditions at the vent for simulations are listed in Table 2. Pluses represent the jet-type column regime. Circles indicate the fountain-type column regime. Diamonds are the jet-type collapse regime. Triangles represent the fountain-type collapse regime. Solid curves are the column collapse condition. Dashed curves are the jet-fountain condition. Dotted thin curves are the column collapse conditions which are predicted by the 1-D model of Woods (1988) with variable entrainment coefficients. The values at the edge of these curves are the assumed values for the entrainment coefficient.

10

Y.J. Suzuki, T. Koyaguchi / Journal of Volcanology and Geothermal Research 221–222 (2012) 1–13

these experimental results for subsonic flows to the case of supersonic flows, we carried out supplementary simulations (Appendix A).

_ air is proportional to the surface area of the eruption colBecause m umn and the exit velocity, w0, it is expressed by

4.2. Dependency of RiCC on the magma properties

_ air ¼ 2πkρa w0 L0 Hmmt m

The CC condition is determined by whether the eruption cloud becomes buoyant or not at z = Hmmt. When the mass fraction of the ejected material, ξ, is smaller than a certain value, ξcrt, the mixture of ejected material and entrained air is lighter than the ambient air. On the other hand, when ξ > ξcrt, the density of the mixture is larger than the atmospheric density. Therefore, the column collapse occurs when the mass fraction of the ejected material at z = Hmmt, ξmmt, is larger than ξcrt. The value of ξmmt is approximately estimated by

under the assumption that the efficiency of turbulent mixing (i.e., the entrainment coefficient, k (Morton et al., 1956)) is constant. Substituting Eqs. (8) and (9) into Eq. (7), we obtain

ξmmt ¼

_0 m ; _ air _ 0þm m

2

ξmmt ¼

πρ0 w0 L0 : πρ0 w0 L20 þ 2πkρa w0 L0 H mmt

RiCC ¼ k

_ 0 ¼ πρ0 w0 L20 : m

ð8Þ

3

ρa ξcrt : ρ0 1−ξcrt

3

(b) group I (T0=800 K n0 =0.0248)

Mach number

Mach number

(T0 =550 K n0 =0.0248) 2 Fountain-type Collapse

1 Jet-type Column Jet-type Collapse 0.1

2 Jet-type Column

Fountain-type Collapse

1

Jet-type Collapse

1

0 0.01

10

Richardson number

0.1

1

10

Richardson number 3

3

(c) group H (T0=1000 K n0 =0.0248)

2 Fountain-type Collapse

Jet-type Column

1

0 0.01

0.1

1

Richardson number

10

Fountain-type Column

Mach number

Fountain-type Column

Mach number

ð11Þ

This result (i.e., Eq. (11)) successfully explains the simulation result that RiCC depends on magma properties even for a fixed M; because 1/ρ0 and ξcrt increase with T0 and ng0, RiCC increases with T0 and ng0. This result is also consistent with the predictions based on the 1-D model of eruption column by Woods and Caulfield (1992) (Appendix B).

(a) group L

0 0.01

ð10Þ

From Eqs. (1), (4), and (10), we obtain the CC condition (i.e., ξmmt = ξcrt) as

ð7Þ

_ 0 is the mass discharge rate at the vent, and m _ air is the mass where m _ 0 is of air entrained below z = Hmmt per unit time. By definition, m represented by

ð9Þ

(d) group N (T0=1000 K n0=0.0123)

2 Jet-type Column

Fountain-type Collapse

1

0 0.01

0.1

1

10

Richardson number

Fig. 13. Flow-regime maps for the results of (a) group L, (b) group I, (c) group H, and (d) group N in the non-dimensional space. Parameters used are listed in Table 1. The conditions at the vent and non-dimensional parameters for simulations are listed in Table 2. The horizontal axis is the Richardson number (Ri) and the vertical axis is the Mach number (M). Pluses represent the jet-type column regime. Circles indicate the fountain-type column regime. Diamonds are the jet-type collapse regime. Triangles represent the fountain-type collapse regime. Solid curves are the column collapse condition. Dashed curves are the jet-fountain condition. Gray and white stars are the critical points for the column collapse condition and those for the jet-fountain condition for M = 1.2, respectively.

Y.J. Suzuki, T. Koyaguchi / Journal of Volcanology and Geothermal Research 221–222 (2012) 1–13

4.3. Dependency of RiCC on the Mach Number The results of our 3-D simulations show that RiCC is independent of the Mach number when RiCC is smaller than RiJF (Fig. 13a,b). When RiCC > RiJF , on the other hand, RiCC decreases as the Mach number increases (Fig. 13c,d). According to Eq. (11), RiCC is independent of the Mach number under the assumption that the efficiency of turbulent mixing (i.e., the entrainment coefficient) is constant. The variation of RiCC in Fig. 13c,d, therefore, means that the efficiency of turbulent mixing changes with the exit velocity for RiCC > RiJF . The efficiency of turbulent mixing for the results of 3-D simulations (the effective entrainment coefficient) can be estimated by comparing the CC condition based on the 3-D simulations with that predicted by the 1-D model (cf. Suzuki et al., 2005). In Fig. 12, the column collapse conditions predicted by the 1-D model with different assumed entrainment coefficients are also illustrated by the dotted curves. For jet-type collapse (RiCC bRiJF ), the column collapse condition based on the 3-D simulations coincides with that of the 1-D model with a constant entrainment coefficient (~ 0.03 in Fig. 12a,b). For fountain-type collapse (RiCC > RiJF ), the effective entrainment coefficient of the 3-D simulations decreases from 0.10 to 0.03 with the increasing Mach number (Fig. 12c,d). We attribute the systematic decrease in the efficiency of turbulent mixing with the increasing exit velocity to the formation of shock wave structures inside the fountain. When the gas–pyroclast mixture is ejected from the vent at high speed, the standing shock waves are generated in the fountain even if the pressure at the vent is balanced with the atmospheric pressure (Fig. 5). These shock wave structures suppress efficient mixing between the ejected material and ambient air. Recent experimental results by Saffaraval et al. (2010) also indicate that the efficiency of entrainment is reduced when shock wave structures are formed in highspeed jets. As the Mach number increases, the flow with a series of shock waves reaches a higher level, which causes less efficient mixing. Thus, the region of the fountain-type collapse regime expands with the increasing Mach number in the Ri–M space (Fig. 13c,d).

11

in some pyroclastic flow deposits (e.g., Druitt and Sparks, 1982; Walker, 1985). Another implication is concerned with the CC condition predicted by the 1-D models. The efficiency of entrainment estimated in the present 3-D numerical simulations (k ~ 0.03 for the jet-type CC condition and k = 0.03–0.10 for the fountain-type CC condition) is substantially smaller than the value commonly assumed in the 1-D models (e.g., k = 0.09 in Woods (1988)). This means that the eruption column collapse is likely to occur for smaller mass discharge rates than the CC condition predicted by those 1-D models. Our results also suggest that the efficiency of entrainment decreases as the exit velocity increases along the fountain-type CC condition (Fig. 12c,d). Small values of k have been proposed in some 1-D models (Carazzo et al., 2008, 2010; Koyaguchi et al., 2010) on the basis of recent laboratory experiments (Kaminski et al., 2005) and direct measurements in numerical simulations (Suzuki and Koyaguchi, 2010). The effect of the exit velocity on the efficiency of entrainment (e.g., Saffaraval et al., 2010), however, has not been taken into account in any 1-D models. We suggest that this effect may play an important role in column collapse, especially during high-temperature eruptions. On the basis of the 3-D numerical simulations, we have shown that the 3-D structures of flow largely affect the column collapse condition. Our results indicate that the differences between the jet-type and fountain-type collapses are closely related with the unsteady vortical structures (e.g., Suzuki and Koyaguchi, 2010). For better understanding of column collapse condition, further analyses based on 3D flow visualization are required. Acknowledgments The manuscript was improved by comments from L. Wilson, D. E. Ogden and an anonymous reviewer. Computations were carried out in part on Earth Simulator at the Japan Agency for Marine-Earth Science and Technology, Altix 4700 at the Earthquake Research Institute, University of Tokyo, and Primergy RX200S6 at Computer System for Research, Kyushu University. Part of this study was supported by the ERI Cooperative Research Program and by KAKENHI (No. 21340123).

5. Geological implications Appendix A. Heights of negatively buoyant jet and fountain The present 3-D numerical simulations of eruption clouds show two types of column collapse (CC) conditions. The CC condition is characterized by the jet-type collapse when the initial temperature of the ejected material is relatively low, whereas the CC condition is characterized by the fountain-type collapse when the initial temperature of the ejected material is high. The flow patterns in the jet-type and fountain-type collapse regimes are considered to affect depositional structures in pyroclastic flows. When the jet-type collapse occurs, the eruption cloud contains a large amount of air, generating a dilute pyroclastic flow. The column collapse during low-temperature eruptions, therefore, tends to generate a dilute pyroclastic flow. Generally, eruptions with low temperatures are attained by the interaction between the magma and water (e.g., Koyaguchi and Woods, 1996). The present numerical results are consistent with the fact that surge deposits from dilute low-temperature pyroclastic flows are commonly observed in the deposits of phreatomagmatic eruptions (e.g., Moore, 1967). When the fountain-type collapse occurs, on the other hand, high-concentration and high-temperature pyroclastic flows are generated. Such pyroclastic flows account for the occurrence of highly welded pyroclastic flow deposits (e.g., Branney and Kokelaar, 1992). During the fountain-type collapse, a large amount of pyroclasts fall to the ground directly from the fountain top, which generates a chaotic flow around the vent as well as outward pyroclastic flows (see Fig. 6). The chaotic flow around the vent may explain complex features of proximal facies

In Section 4.1, we presented a simple analysis for the JF condition on the assumption that the fountain is generated when Hcore > Hmmt (Eq. (6)). The analysis is based on previous experimental studies on the transition between negatively buoyant jets and fountains for subsonic flows (e.g., Kaye and Hunt, 2006). According to the experimental results, the JF condition of negatively buoyant jet/fountains is closely related to the relationship between the maximum height (Hmax) and Ri. When Ri is sufficiently small so that Hcore b Hmmt (see also Eq. (5)), Hmax of the negatively buoyant jet increases with Ri for a fixed exit velocity. This is because the ratio of the mass of entrained air to that of the ejected material decreases as Ri increases, and hence, Hcore/Hmmt increases. When Ri is large so that Hcore > Hmmt, on the other hand, Hmax of the negatively buoyant fountain remains Hmmt for a fixed exit velocity and is independent of Ri. These results indicate that, at least for subsonic flows under the experimental conditions, the JF condition can be determined by the critical value of Ri at which the Hmax–Ri relationship changes. Here, we attempt to extend the above idea to the JF condition of supersonic flows; we carried out supplementary simulations for negatively buoyant jet/fountains of a gas–pyroclast mixture for various Mach numbers. In order to avoid unnecessary complications (e.g., rapid expansion of the cloud due to heating of the entrained air) and to extract the effect of the Mach number on the JF condition, the temperature of the gas–pyroclast mixture was set to be equal to the ambient temperature (273 K). Under

12

Y.J. Suzuki, T. Koyaguchi / Journal of Volcanology and Geothermal Research 221–222 (2012) 1–13

Height of Negatively Buoyant Jet [m]

104

system can be derived on the basis of the 1-D model with a constant entrainment coefficient (k) as: M=3

 2 _ 5=2 4kα M W0 W0 0 −1 ≥ ; _ 30 W m 2W m 5ðβm −α Þg m

M=2

103

M=1

0.01

0.1

1

10

Richardson Number Fig. A.1. Maximum heights of the negatively buoyant jet/fountains (Hmax) as a function of the Richardson number (Ri) on the basis of the 3-D numerical simulations. White and black points represent the jet-type and fountain-type regimes, respectively. Gray plots are the intermediate regime between the jet-type and fountain-type regimes. Circles, diamonds, and triangles indicate the results for M = 1, 2, and 3, respectively. Arrows represent the JF condition based on the H max –rmRi relationship.

_ (≡ w02L02) is the iniwhere W0 is the initial volume fraction of MEG, M 0 _ 0 (≡ w0L02) is the initial mass flux. It should tial momentum flux, and m be noted that this formula is applicable only to the MEG–water system, because it is based on the fact that the density of the mixture is equal to the ambient density when W = 2Wm(≡ Wcrt) in Eq. (B.1). Here, we extend Eq. (B.2) to the general case and show that it is consistent with our result (i.e., Eq. (11)). In order to obtain the CC condition for the general case, Wm in Eq. (B.2) must be given as a function of general quantities such as the critical volume fraction (Wcrt) or the critical mass fraction (ξcrt) of ejected materials. Considering that Eq. (B.1) and Wcrt = 2Wm are the relationships for the MEG–water system, Eq. (B.2) is rewritten using the general quantities as −Ri≥

this condition, the cloud density is always larger than the ambient density. Fig. A.1 shows the Hmax–Ri relationship of the negatively buoyant jet/fountains for M = 1–3. For sonic flows (M = 1), negatively buoyant jets develop and Hmax increase with Ri when Ri b 0.6. When rmRi > 0:6, negatively buoyant fountains with a constant Hmax develop. The results show that the critical Richardson number for the JF condition (RiJF ) based on the change in the Hmax–Ri relationship decreases with the increasing Mach number (arrows in the figure). From this RiJF –M relationship, we obtain the values of a1 and a2 in Eq. (6); the estimated values are a1 = 1.25 and a2 = 0.42, which are consistent with those estimated from Fig. 13 (a1 = 1.3–3.4 and a2 = 0.16–4.7). Appendix B. Derivation of a formula for the column collapse conditions by 1-D model on the basis of analogue experiments Woods and Caulfield (1992) carried out a series of analogue experiments simulating column collapse and proposed a simple formula for the CC condition. In their experiments, mixtures of methanol and ethylene glycol (MEG) were injected as a downward propagating jet into fresh water. The MEG is lighter than fresh water, whereas the MEG–water mixture can be denser than fresh water. The density difference between the fresh water and MEG–water mixture was approximated by a function of the volume fraction of MEG in the mixture, W, as Δρ≡

ðB:2Þ

  β−α W 2 ¼ 1− 1− ; βm −α Wm

ðB:1Þ

where β is the MEG–water mixture density and α is the water density. The MEG–water mixture has a maximum density of βm when W = Wm. If the MEG entrains sufficient water before its initial momentum is exhausted, the MEG–water mixture becomes denser than the ambient water and continues to flow downward as a negatively buoyant “plume”. If the MEG does not mix with sufficient water, the mixture “collapses” and rises back toward the top of the tank. Although the direction of the flow is upside-down, the experiments reproduce fundamental features of eruption column collapse. Woods and Caulfield (1992) showed that, when the Δρ–W relationship is explicitly given as Eq. (B.1), a simple formula of the CC condition (i.e., the region of collapse regime) for the MEG–water

8k W crt : 5 W 0 −W crt

ðB:3Þ

Here, the definition of Ri is given by Eq. (1) (note that Ri is negative because the direction of the flow is upside-down). Furthermore, because the mass fraction is expressed by ξ¼

Wβ0 ; Wβ0 þ ð1−W Þα

ðB:4Þ

where β0 is the initial density of MEG, Eq. (B.3) is rewritten as −Ri≥Ricrt ≡

8k α ξcrt : 5 β0 1−ξcrt

ðB:5Þ

This has the same form as our formula for the CC condition (i.e., Eq. (11)). These two formulae have different proportionality constants (1.0 in Eq. (11) and 8/5 in Eq. (B.5)). This difference is attributed to the fact that the mass of entrained air is assumed to be constant at any heights in our analysis (see Eq. (9)), whereas it is calculated on the basis of the 1-D model with a constant k in Woods and Caulfield (1992). References Branney, M.J., Kokelaar, P., 1992. A reappraisal of ignimbrite emplacement: progressive aggradation and changes from particulate to non-particulate flow during emplacement of high-grade ignimbrite. Bull. Volcanol. 54, 504–520. Bursik, M.I., Woods, A.W., 1991. Buoyant, superbuoyant and collapsing eruption columns. J. Volcanol. Geotherm. Res. 45, 347–350. Carazzo, G., Kaminski, E., Tait, S., 2008. On the dynamics of volcanic columns: a comparison of field data with a new model of negatively buoyant jets. J. Volcanol. Geotherm. Res. 178, 94–103. Carazzo, G., Kaminski, E., Tait, S., 2010. The rise and fall of turbulent fountains: a new model for improved quantitative predictions. J. Fluid Mech. 657, 265–284. Carey, S., Sigurdsson, H., Gardner, J.E., Criswell, W., 1990. Variations in column height and magma discharge during the May 18, 1980 eruption of Mount St. Helens. J. Volcanol. Geotherm. Res. 43, 99–112. Dimotakis, P.E., 2000. The mixing transition in turbulent flows. J. Fluid Mech. 409, 69–98. Dobran, F., Neri, A., Macedonio, G., 1993. Numerical simulation of collapsing volcanic columns. J. Geophys. Res. 98, 4231–4259. Druitt, T.H., Sparks, R.S.J., 1982. A proximal ignimbrite breccia facies on Santorini, Greece. J. Volcanol. Geotherm. Res. 13, 147–171. Kaminski, E., Jaupart, C., 2001. Marginal stability of atmospheric eruption columns and pyroclastic flow generation. J. Geophys. Res. 106, 21785–21798. Kaminski, E., Tait, S., Carazzo, G., 2005. Turbulent entrainment in jets with arbitrary buoyancy. J. Fluid Mech. 526, 361–376. Kaye, N.B., Hunt, G.R., 2006. Weak fountains. J. Fluid Mech. 558, 319–328. Koyaguchi, T., Woods, A.W., 1996. On the formation of eruption columns following explosive mixing of magma and surface-water. J. Geophys. Res. 101, 5561–5574. Koyaguchi, T., Suzuki, Y.J., Kozono, T., 2010. Effects of the crater on eruption column dynamics. J. Geophys. Res. 115, B07205.

Y.J. Suzuki, T. Koyaguchi / Journal of Volcanology and Geothermal Research 221–222 (2012) 1–13 Lin, W., Armfield, S.W., 2000. Direct simulation of weak axisymmetric fountains in a homogeneous fluid. J. Fluid Mech. 403, 67–88. Moore, J.G., 1967. Base surge in recent volcanic eruptions. Bull. Volcanol. 30, 337–363. Morton, B.R., Taylor, G.I., Turner, J.S., 1956. Turbulent gravitational convection from maintained and instantaneous sources. Proc. R. Soc. Lond. A 234, 1–23. Nagamatsu, H.T., Sheer, R.E.J., Horvay, G., 1969. Supersonic jet noise theory and experiments. Basic Aerodynamics Noise Research. NASA SP-207, pp. 17–51. Neri, A., Dobran, F., 1994. Influence of eruption parameters on the thermofluid dynamics of collapsing volcanic columns. J. Geophys. Res. 99, 11833–11857. Oberhuber, J.M., Herzog, M., Graf, H.F., Schwanke, K., 1998. Volcanic plume simulation on large scales. J. Volcanol. Geotherm. Res. 87, 29–53. Ogden, D., Glatzmaier, G.A., Wohletz, K.H., 2008a. Effects of vent overpressure on buoyant eruption columns: implications for plume stability. Earth Planet. Sci. Lett. 268, 283–292. Ogden, D.E., Wohletz, K.H., Glatzmaier, G.A., Brodsky, E.E., 2008b. Numerical simulations of volcanic jets: importance of vent overpressure. J. Geophys. Res. 113, B02204. Rajaratnam, N., 1976. Turbulent jets. Elsevier. Roe, P.L., 1981. Approximate Riemann solvers, parameter vectors, and difference schemes. J. Comput. Phys. 43, 357–372. Saffaraval, F., Solvitz, S., Mastin, L.G., 2010. Experimental of near-field entrainment of high-speed jets. AGU Fall Meeting 2010, V13C–2379. Sparks, R.S.J., Wilson, L., 1976. A model for the formation of ignimbrite by gravitational column collapse. J. Geol. Soc. Lond. 132, 441–451. Suzuki, Y.J., Koyaguchi, T., 2009. A three-dimensional numerical simulation of spreading umbrella clouds. J. Geophys. Res. 114, B03209. Suzuki, Y.J., Koyaguchi, T., 2010. Numerical determination of the efficiency of entrainment in volcanic eruption columns. Geophys. Res. Lett. 37, L05302.

13

Suzuki, Y.J., Koyaguchi, T., Ogawa, M., Hachisu, I., 2005. A numerical study of turbulent mixing in eruption clouds using a three-dimensional fluid dynamics model. J. Geophys. Res. 110, B08201. Taylor, G.I., 1945. Dynamics of a mass of hot gas rising in air. U. S. Atomic Energy Commission MDDC 919, LADC276. Valentine, G.A., Wohletz, K.H., 1989. Numerical models of Plinian eruption columns and pyroclastic flows. J. Geophys. Res. 94, 1867–1887. van Leer, B., 1977. Towards the ultimate conservative difference scheme, III. Upstreamcentered finate-difference schemes for ideal compressible flow. J. Comput. Phys. 23, 263–275. Walker, G.P.L., 1985. Origin of coarse lithic breccias near ignimbrite source vents. J. Volcanol. Geotherm. Res. 25, 157–171. Wilson, L., Sparks, R.S.J., Walker, G.P.L., 1980. Explosive volcanic eruptions-iV, the control of magma properties and conduit geometry on eruption column behaviour. Geophys. J. R. astron. Soc. 63, 117–148. Wohletz, K.H., McGetchin, T.R., Stanford II, M.T., Jones, E.M., 1984. Hydrodynamic aspects of caldera-forming eruptions: numerical models. J. Geophys. Res. 89, 8269–8285. Woods, A.W., 1988. The fluid dynamics and thermodynamics of eruption columns. Bull. Volcanol. 50, 169–193. Woods, A.W., Bower, S.M., 1995. The decompression of volcanic jets in a crater during explosive volcanic eruptions. Earth Planet. Sci. Lett. 131, 189–205. Woods, A.W., Caulfield, C.-C.P., 1992. A laboratory study of explosive volcanic eruption. J. Geophys. Res. 97, 6699–6712.