Computational modelling of the transient dynamics of the August 1997 Vulcanian explosions at Soufriere Hills Volcano, Montserrat: influence of initial conduit conditions on near-vent pyroclastic dispersal A. B. CLARKE 1 , A. NERI 2 , B. VOIGHT1, G. MACEDONIO3 & T. H. DRUITT4 1 Department of Geosciences, Penn State University, University Park, PA 16802, USA (e-mail:
[email protected]) 2 CNR-CSGSDA, Department of Earth Sciences, Pisa, Italy 3 Osservatorio Vesuviano, Napoli, Italy 4 Laboratoire Magmas et Volcans, Universite Blaise Pascal et CNRS, Clermont-Ferrand 63038, France Abstract: This paper presents numerical models of the Vulcanian explosions that occurred in 1997 at Soufriere Hills Volcano. Plume evolution and velocities were calculated for the well-documented and typical explosions of 6 and 7 August 1997, and these data and other observations were compared to transient, axisymmetric, multiphase flow simulations of coupled conduit evacuation and pyroclastic dispersal. Pre-explosion conduit conditions were estimated from Montserrat data, using a simple gas solubility law and assuming that conduit magma flow had stagnated with a constant overpressure prior to the explosions. Reference simulation input parameters include conduit diameter of 30m, crater diameter of 300m, meltwater content of 4.3±0.5wt%, grain sizes of 30, 2000 and 5000um, and conduit overpressure of l0MPa. The numerical simulations of the explosions resolved highly unsteady vent exit conditions such as velocity, pressure and mass flux, and the spatial and temporal dispersal of pyroclasts during the initial few minutes was investigated using one gas phase and two or three solid phases representing pyroclasts of different size. Our simulations produced transitional eruptive regime behaviour, dividing the erupted mass into a portion that generated a radial pyroclastic current fed by a collapsing column, and a convective portion that generated a buoyant plume. This behaviour generally mimicked the observed explosions. The movement of different particle sizes was tracked, with fine particles dominantly influencing the convective behaviour of the central plume and ash plume thermals generated above the pyroclastic currents. Simulated initial vent velocities ranged from 85 to 120 m s - 1 , collapse heights ranged from 450 to 1370m above the vent, initial pyroclastic current velocities ranged from 40 to 6 0 m s - 1 with surge runouts to 1.8km, drawdown depths in the conduit were a few hundred metres, and simulated pyroclastic current deposit temperatures ranged between 135 and 430°C. Subsets of these results are in reasonable agreement with observed and measured parameters of the 1997 explosions. The best match was intermediate between our reference simulation, which assumed no loss of volatiles from the conduit during rise from the magma reservoir and which appeared too energetic, and another simulation in which much volatile leakage was assumed. The results suggest that volatile depletion in the conduit was an important factor in influencing the dynamic behaviour of the Vulcanian explosions on Montserrat.
A total of 88 short-duration Vulcanian explosions, nearly all accompanied by radial fountain collapse, occurred at Soufriere Hills Volcano, Montserrat in 1997 (see Fig. 1). Thirteen occurred in August and 75 more occurred during September and October. These explosions provided an unprecedented opportunity for repeated observation and monitoring (Fig. 2; Druitt et al 2002), as only a few events of this type have been observed closely (Nairn & Self 1978; Sparks & Wilson 1982; Hoblitt 1986). Because of the intensive real-time monitoring of precursory seismicity and tilt, it was possible over the short term to forecast (on the timescale of several hours with an uncertainty of tens of minutes) the onset of many of these explosions (Voight et al 1998, 1999; Druitt et al 2002), thus facilitating scientific preparations for impending explosions. Therefore, many explosions were documented in unusual detail by video and repetitive still photography, by theodolite surveying of eruption plumes, by tephra sampling, and by monitoring of broadband seismicity and deformation. Numerical models of short-duration Vulcanian explosions with time-varying vent flux have not received much attention, and no eruption model to date has tried to combine highly unsteady vent dynamics with explosive dispersal of pyroclastics. Recent improvements in numerical code development provided new opportunities for our study. Therefore, we have used an axisymmetric, multipleparticle-size, numerical code (Neri 1998; Neri et al. 2001b) to calculate the unsteady and transient vent flux following vent cap failure and resulting pyroclast dispersal (explosion) for a number of assumed pre-explosion conduit profiles. The observational data fall into two categories: (1) those which constrained initial conduit conditions and input parameters for our numerical models of explosive pyroclast dispersal; and (2) those which describe the explosions and provide a basis for comparison with the results of the numerical models.
In this paper, we first present a general description of the 1997 Montserrat explosions, and follow this by a brief review of eruption modelling, focusing on what is known about fountain-collapse pyroclastic current generation and Vulcanian explosions. Next, we summarize the pyroclast dispersal model we used (Neri 1998; Neri et al. 2001b) and our methods for determining the pre-explosion conduit conditions. The solutions focus on the near-vent region (within a few kilometres) over the first 90 to 150 s of the explosions. Then we use several numerical simulations to investigate the influence of key input parameters on results, namely conduit overpressure, mass fraction of water, and fines fraction. Finally, we compare our results against observations and inferred parameters of the 1997 explosions at Soufriere Hills Volcano in order to test the pyroclast dispersal (explosion) model and its assumptions. Characteristics of the cyclic Vulcanian explosions An overview of the 1997 Vulcanian explosions and their physical parameters is given in Druitt et al. (2002). The main features relevant to our numerical simulations are summarized here. The Vulcanian explosions of 1997 occurred in two episodes, the first between 4 and 12 August, and the second between 22 September and 21 October. The explosions began with an initial phase of ballistic throwout followed by a phase of powerful fountaining of a pyroclastic mixture, which lasted 10-20s. Condensation of atmospheric moisture caused by shock waves was noted on video footage of some explosions. The condensation patterns propagated radially away from the vent ahead of the expanding fountain of gas and pyroclasts during the first few seconds of the explosions. The vertical jets collapsed to generate fountains and high-velocity radial pyroclastic surges that moved to distances of about 2 km (Fig. 3),
DRUITT, T. H. & KOKELAAR, B. P. (eds) 2002. The Eruption of Soufriere Hills Volcano, Montserrat, from 1995 to 1999. Geological Society, London, Memoirs, 21, 319-348. 0435-4052/02/$15 © The Geological Society of London 2002.
319
320
A. B. C L A R K E ET AL.
Fig. 1. Map of southern Montserrat and Soufriere Hills Volcano, showing the active lava dome (vent) location inside English's Crater, the principal drainages (ghauts) about the dome, and the area affected by pyroclastic surges and flows during 1995-1999 (after Druitt et al. 2002). Photograph and video observation points include Montserrat Volcano Observatory (MVO South) and Fleming.
and also high-concentration pumice-and-ash flows that were generally confined to channels and ran out to distances of 3 to 6 km. Subsequently, buoyant-convecting plumes developed from the fountains and rose to heights of 3 to 15 km, where they generally spread out as umbrella clouds. The explosions were dominated by the buoyant plumes for roughly 500 s. Buoyant ash plume thermals also formed above the surges and channelled flows and tended to rise and move inward to join the main plume. The explosions ended with approximately an hour of waning exhalations characterized by a bent-over plume. Collapse of the fountain generally occurred 10 to 20s after explosion initiation, from heights of 300-650 m above the crater rim (400-750 m above the vent). The individual explosions expelled on average 3.0 x 10 5 m 3 of magma 1.1 x 10 5 m 3 as fallout and 1.9 x 105 m3 as surges and pumice-and-ash flows, evacuating the conduit to depths of 500 to 2000 m. Vent exit velocities ranged from 40 to 140ms - 1 . The pyroclastic surges developed slope-parallel velocities of 30-60 m s - 1 , whereas the channelled pumice-and-ash flows typically had velocities of 10 m s-1. The emplacement temperatures of pumice-and-ash flows ranged from 180 to 2200C (Druitt et al. 2002; Cole et al. 2002), with air entrainment cooling the mixture with respect to its eruption temperature of about 8600C (Barclay et al. 1998). The explosions occurred from a conduit 30m ( 5m) in diameter (Voight et al. 1999) into a flared crater, approxi-
mately 300m (±20 m) in diameter. Our paper focuses mainly on two well-documented events at 14:35 local time (LT) on 6 August and at 12:05 LT on 7 August 1997 (all times given are local times), with simulations attempting to reproduce the near-vent behaviour over the first 90 to 150s. Background on explosion modelling and overview of our model Fountain collapse occurs when an ejected pyroclast-gas jet does not have enough momentum to continue rising, and, having failed to become positively buoyant, falls back toward and across the ground surface as a dispersion of hot particles (fragmented magma and lithic clasts) and gas (Sparks & Wilson 1976; Sparks et al. 1978; Neri & Macedonio 1996a). Early understanding of plume behaviour was derived from numerical work on turbulent gravitational convection (Morton et al. 1956). Experimental and steady-state, pseudogas models Sparks et al. (1978) described the physics of column collapse and pyroclastic flows by applying steady-state conservation of mass
MODELLING DYNAMICS OF VULCANIAN EXPLOSIONS
321
Fig. 2. Photographs of the explosion at 14:35 LT on 6 August 1997. Times after the low-frequency onset of the explosion signal at 14:35:12 LT: (a) 22s (b) 37s (c) 55 s (d) 91 s (e) 96 s. In (a) the column was about 800m high and collapse had begun, but pyroclastic currents had not clearly penetrated the veil of fallout. Gages Mountain (Fig. 1) is in foreground below plume. In (b), pyroclastic currents were descending Mosquito Ghaut (left) and Gages valley (right), and other major drainages around the volcano. In (c), the pyroclastic currents had thickened by development of dilute ash plumes, and the bulbous central plume, depleted in coarse particles, was rising buoyantly. By (d), active tan-coloured ash-rich thermals that developed over the pyroclastic current, and rose upward and inward, were being incorporated into the stalk of the buoyant central plume. A full view of the central plume is shown in (e), which rose ultimately to about 12km above sea level, where it spread to form an umbrella cap that was subsequently sheared from the main column and transported NE by high-altitude winds. Photographs taken from 7 km NW at MVO South by B. Voight.
and momentum to a gas-particle mixture in which the particles and gas are in thermal and kinetic equilibrium. Bursik & Woods (1996) added the conservation of energy expression to the description of the ash flows but did not develop the fountain model further. In all these studies the gas-particle mixture, or pseudogas, is treated as a single-phase fluid with bulk properties determined by the volumetric proportion, size and temperature of the particles. The main parameters controlling fountain formation were found to be the vent radius, gas content and initial vent velocity (or a combination of exsolved gas content and vent velocity, i.e. vent mass flux), whereas flow runout is also controlled by air entrainment and particle sedimentation (Bursik & Woods 1996). Reasonable agreement has been claimed between theory and observations (Wilson et al 1978; Turner 1979; Sparks & Wilson 1982; Sparks 1986); however, as noted by the original authors, these relationships are valid only within the context of one-dimensional steady-state pseudogas theory, where pyroclasts are very fine-
grained, and they were generally developed for steady Plinian-type eruptions. They have been used to estimate first-order time-averaged column behaviour, but the required simplifications impede detailed comparison with the transient, unsteady, multidimensional, and multiphase nature of real eruptions. Shock-tube and gas-particle experiments also suggest that the pseudogas, steady-state, onedimensional approximation may not be a good assumption for highvelocity two-phase volcanic flows (Anilkumar et al. 1993; Sparks et al. 1997; Neri & Gidaspow 2000).
Steady discharge, multiphase models Computer codes that were originally designed for atomic explosion simulation or fluidization studies have been adapted and further developed for the study of volcanic eruptions (Harlow & Amsden
322
A. B. C L A R K E ET AL.
Fig. 2. (continued)
1971, 1975; Amsden & Harlow 1974; Valentine & Wohletz 1989; Valentine et al. 1992; Dobran et al 1993; Gidaspow 1994). These codes enabled solutions to time-dependent, two-phase (one solid, one gas) compressible Navier-Stokes equations. They addressed thermal and kinetic disequilibrium between the solid and gas phases and time-dependent behaviour of the plume and the flows, while holding vent conditions steady. In addition to the previous findings of Sparks et al (1978) and Woods (1988) these results showed that plume behaviour is sensitive to the ratio of conduit pressure to atmospheric pressure at the vent (Valentine & Wohletz 1989). The codes also led to recognition of important relationships between plume behaviour and simulation particle size. Under otherwise identical vent conditions, simulations with larger particle sizes (larger Rouse numbers) are more likely than those with smaller particle sizes to fall back and form pyroclastic currents (Valentine & Wohletz 1989). This is caused by two phenomena, namely slower heating of the surrounding gas by the larger particles and their higher settling velocity.
Neri & Macedonio (1996b) made a first attempt to account for the grain-size distribution of the eruptive mixture by adding a second particle phase to the numerical solution of Dobran et al. (1993). Different drag terms between gas and particles were included for different particle sizes, and a drag term was added to account for collisions between particles of different sizes. Their model results indicate that the particles of different sizes have considerably different dynamics and affect one another's behaviour, thus changing the resulting plume behaviour and runout dynamics of the pyroclastic currents. Although the aforementioned numerical models provide insight into many of the factors that control eruption column behaviour, like the pseudogas models described above, they specifically address Plinian-style eruptions where the vent conditions are quasi-steady. None directly treats the sudden decompression of a pressurized conduit or chamber, or resulting short-pulse explosions, and therefore they do not directly apply to the Vulcanian explosions witnessed at Soufriere Hills Volcano in 1997.
MODELLING DYNAMICS OF VULCANIAN EXPLOSIONS
Fig. 2. (continued)3
323
324
A. B. CLARKE ET AL.
Fig. 3. Schematic representation of selected features of Vulcanian explosions.
Vulcanian models Several studies have explored the physics of short-duration Vulcanian explosions, with important contributions by Self et al. (1979), Kieffer (1981), Turcotte et al. (1990), Fagents & Wilson (1993) and Woods (1995). Self et al. (1979) treated the mixture as a pseudogas, and compared results to Vulcanian explosions observed at Ngauruhoe, New Zealand. Kieffer (1981) developed equations to describe the sudden decompression of a high-pressure magma chamber due to the failure of the north flank of Mount St Helens on 18 May 1980. Fagents & Wilson (1993) examined ranges of pyroclasts ejected. Turcotte et al. (1990) developed a set of equations based on the one-dimensional shock-tube problem, using unsteady conservation laws for mass and momentum and treating the mixture as a pseudogas. Woods (1995) extended the Turcotte et al. model by allowing the mixture to cool adiabatically, and by including a factor to account for particles not in thermal equilibrium with the gas. In addition, Wohletz et al. (1984) used the KACHINA numerical code (Amsden & Harlow 1974) to simulate a caldera-forming eruption through a sudden decompression of a non-homogeneous magma chamber. The code treats one solid phase and one gas phase, and focuses on the initial unsteady blast phase, and on shock wave phenomena.
Our model The previous Vulcanian models addressed transient vent conditions, and related the initial conduit pressure to exit conditions. In this paper, we advance the work of our predecessors by (1) applying a two-dimensional model with multiple particle sizes, which solves for unsteady conduit flow and resulting pyroclast dispersal, (2) constraining pre-eruptive conduit conditions with analyses of observational data, (3) comparing explosion model simulations with well-documented explosions, and (4) varying some key conduit parameters in order to better understand their effects on explosion simulation results.
The pyroclast dispersal model, named PDAC2D (Neri 1998; Neri et al. 200la,b), is an extension of the three-phase flow code used by Neri & Macedonio (I996b) and solves a set of equations expressing the conservation of mass, momentum and energy for one gas phase, and a number of solid phases representative of particles of different sizes. The gas phase consists of a mixture of water vapour and atmospheric air. The fundamental transport equations are solved on an axisymmetric computational grid (5m to 40m grid spacing), as discussed in detail in the Appendix.
Input conditions for explosion model For our simulations of short-duration Vulcanian explosions, several input parameters were required. The most important of these were the conduit, crater and plugging cap geometry, the twodimensional (axisymmetric) topography of the region surrounding the vent, the initial conduit gas pressure, gas mass fraction and gas volume fraction as functions of depth in the conduit, and the sizes and densities of solid particles. Information acquired on Montserrat was used to constrain these input parameters, as described in the following sections.
Conduit and crater diameters, cap thickness and ground topography The conduit diameter of approximately 30m (±5m) was constrained by spine dimensions, magma ascent rates and volume extrusion rates (Voight et al. 1999). The crater diameter of roughly 300m (±20 m) was measured using dual-position theodolite measurements approximately one half hour prior to the 12:05 event on 7 August. The crater depth was approximately 100 m. The cap thickness was assumed to be 20m, which is a size roughly consistent with the estimated total volume of dense clasts per explosion (roughly 1-5% by volume of pumice-and-ash flow deposits).
Fig. 4. Topography used in the simulations. The grid is defined for an axisymmetric solution, thus the crater and conduit dimensions are radii.
MODELLING DYNAMICS OF VULCANIAN EXPLOSIONS
We used the ground topography representative of the north side of the volcano (Fig. 4). The north sector is characterized by a relatively smooth, low-slope fan with average slope approximately 22°, bordered by channels which, by mid-summer 1997, were nearly filled by pyroclastic debris. The model sloping-fan topography thus was more or less realistic for a broad sector of the near-vent region. The model assumption of axisymmetric flow is a reasonable approximation to reality, because the Vulcanian events in August of 1997 are described as axisymmetric (Druitt et al 1998, 2002), with fountain collapses and radial pyroclastic surges occurring simultaneously in all sectors and with pyroclastic current runouts similar in all directions.
Pressure and gas fraction profiles along the pre-explosion conduit For this study, the pre-explosion conduit was modelled quite simply, mostly on the basis of Sparks (1997). We assumed that magma flow in the conduit immediately prior to the explosions had stagnated due to viscous plugging at the vent. Specifically we assumed constant overpressure with depth, since for magma systems with large vertical viscosity gradients, almost all excess pressure drop from the chamber to the surface is concentrated near the top of the conduit (due to the higher viscosity of the magma near the surface), with excess pressure being roughly constant below the upper reaches of the conduit. There is some basis for assuming stagnated flow, because the Vulcanian explosions occurred in association with oscillatory flow (Voight et al. 1999), and modelled conduit flow rates were relatively small just prior to explosions (Wylie et al. 1999; Denlinger & Hoblitt 1999). The actual flow conditions and overpressure distributions were undoubtedly more complicated (Melnik & Sparks 2002a), but our assumptions were practical for a first approximation. We estimated gas mass and gas volume fractions as functions of depth in the conduit using the following solubility and hydrostatic laws. Symbols are summarized in Table 1. The mass fraction, ne, of exsolved water vapour in the bulk magma (melt + crystals + vesicles) at a given conduit depth was taken according to Henry's law as: (1)
where n0 is the mass fraction of water in the melt, z is depth in the conduit, Pg(z) is the total gas pressure in the conduit at depth z, s and (3 are experimentally determined constants for the appropriate melt composition, and 0 is the total volume fraction of crystals in the upper conduit magma (excluding vesicles). The water dissolved in the melt in the chamber, n0 was taken as 4.3 ± 0.5 wt% (Devine et al. 1998; Barclay et al. 1998) and the total crystal volume fraction, was taken to a constant 0.65 (Murphy et al. 2000). The chemistry of the melt phase is rhyolitic (Murphy et al. 1998), for which s = 4.1 x l0 -6 N1/2m -1 and 0 = 0.5 (Wilson et al. 1980). The pressure distribution, Pg(z), is given by: Jo
(2)
where AP is the overpressure at a given depth, which we assume to be constant with depth for each of the simulations presented in this paper, p is the bulk density of the overlying gas-magma mixture, and g is acceleration due to gravity. The water vapour is treated as an ideal gas with density, pg(z), as follows:
325
Table 1. Symbols and variables Symbol Description Cross-sectional area of ballistic block Acoustic velocity in atmospheric air Acceleration of ballistic block in j direction Acoustic velocity in water vapour Drag coefficient for kth solid particle size Specific heats of gas and kth solid particle size Smagorinsky's constant (assumed to be 0.1) Particle diameter of kth and jth solid particle sizes Gas-solid drag coefficient Drag coefficient between kth andy'th solid particle sizes Restitution coefficient for a particle-particle collision Specific expansion energy of conduit water vapour Acceleration due to gravity Solid elastic modulus Fountain collapse height (above the vent) Gas enthalpy Thermal conductivity of kth solid particle size Mass of ballistic block Mach number of shock wave Mass fraction of exsolved water vapour in the bulk magma Mass fraction of water in the melt Nusselt number for kth particle size Atmospheric pressure Total gas pressure Total gas pressure in the conduit (at depth z) Pressure on high-pressure side of shock wave (see Fig. 5) Prandtl number Heat transfer coefficient between the gas and the kth particle size Ideal gas constant for water vapour Radius of conduit Reynolds number of kth particle size Solubility constant for Henry's law Time from onset of explosion Viscous stress tensor for the gas and kth solid particle size DRE volume available in the conduit DRE volume ejected from the conduit Maximum vent velocity for a given simulation Velocity of ballistic in j direction Velocity of water vapour Velocity of kth and j3th solid particle sizes Depth in conduit Restitution coefficient for non-head-on collisions Exponent for Henry's law Gas overpressure Gas volume fraction E Volume fraction of kth solid particle size k Ekj Maximum solids volume fraction for particles of size j and k ES Solid volume fraction in conduit Specific heat ratio for atmospheric air atm m Specific heat ratio for water vapour-magma mixture w Specific heat ratio for water vapour Gas molecular viscosity g Effective gas viscosity ge Effective gas turbulent viscosity gt k Viscosity of kth solid particle size Particle volume fraction of kth or jth size at maximum packing k.j p Bulk density of gas-magma mixture pa Atmospheric air density Pg Gas density Pkj Density of kth and jth solid particle sizes ps Density of melt-crystal mixture (solid) in the pre-explosion conduit Gas deformation tensor g Coulombic repulsive component among solid particles Solid stress tensor for kth particle size
Ah aatm aj aw CD..k Cpgik Cs dk,j Dgk Dkj e Ev g G(Eg) hc hg Kk M Ms ne n0 Nuk Patm Pg Pg(z) P2 Pr Qk R R Rek s t Tg,k VDREA VDREE Vmax Vj vg vkj z a (3 P £g
(3)
where R = 4621kg -1 K - 1 i s the universal gas constant divided by the molecular weight of water, and T is the temperature of the water vapour, which we assume to be isothermal with the meltcrystal mixture (from here on referred to as 'solid') at 1133K
(860°C). If the density of this solid phase, ps, is constant with depth, then the bulk density of the gas-solid mixture is: (4)
326
The solid volume fraction
A. B. CLARKE ET AL. s
at a given depth is: (5)
And the gas volume fraction is simply: (6)
At the depth where the solid fraction s — 0.70 (30% vesicularity), which represents a bulk density slightly below 2000 kg m - 3 , a solid boundary is assumed in our pyroclastic dispersal model. Although this assumption is arbitrary, it is necessary to specify a solid boundary at the base of the conduit in order to effectively simulate the dispersal of the particles. This does not appear to be an unreasonable assumption because our pumice density measurements show that only 13% of randomly sampled clasts have vesicularities less than 30%. Because magma viscosity was very high, bubble expansion after fragmentation was probably minimal (Thomas et al. 1994; Druitt et al, 2002). Therefore, we assume that the vesicularity of pumice represents essentially the pre-explosion state of magma in the conduit, allowing us to claim that roughly 87% of material ejected had pre-explosion solid volume fraction, s < 0.70 (>30% vesicularity). The depth at which s = 0.70 (30% vesicularity) is referred to as effective conduit depth, DE. The dense rock equivalent (DRE) volume of solid material above this boundary is called DRE volume available, VDREA. The DRE volume of solid material ejected permanently from the conduit, that is, that which does not fall back into the conduit during the simulation, is termed DRE volume ejected, VDREE .
Pre-explosion specific expansion energy It should be noted that changing a single conduit parameter (e.g. overpressure or bulk H2O mass fraction), changes several characteristics of the modelled pre-explosion conduit, including VDREA and the effective depth of the simulation conduit. Because of this, it is difficult to attribute differences in simulated explosion behaviour to a single changed conventional conduit parameter. Therefore, to assist interpretation, we combined total pressure, exsolved gas mass fraction, DRE volume available, and conduit depth into a single parameter. This parameter represents the expansion energy due to gas overpressure per unit mass of solid material and is called specific expansion energy, Ev, which provides a convenient way to compare simulation results. Due to the very fast decompression of the system, we assumed adiabatic conditions (Woods 1995) and used the following equation to calculate this energy (Self et al. 1979; Neri et al. 1998):
(7)
where is the specific heat ratio for water vapour (the specific heat at constant pressure divided by the specific heat at constant volume: . = 1.25), Patm is atmospheric pressure, Pg(z) is the total gas pressure, DE is the effective depth of the conduit, r is the radius of the conduit, and is the mass of the material for a segment of the conduit of depth dz.
Ballistic analysis Numerous ballistics were observed during the initial seconds of the Vulcanian explosions of early August 1997. These ballistics were not affected by the motion of the plume because, in general, they were very large and dense, with fall velocities close to the initial explosion gas velocity, making them Type I ballistics according to Self et al. (1980). Documentation of the resulting crater size and position, and clast diameters, estimated from a helicopter, resulted
in a set of 24 ballistic range and size data (Druitt et al. 1998, 2002). We have used these data, along with basic dynamics equations and drag relationships, to estimate the initial velocities of the ballistics (Wilson 1972; Fagents & Wilson 1993; Waitt et al. 1995). The basic equation adopted for the ballistic analysis is dv j /dt = aj where aj is the acceleration (or deceleration) of the ballistic due to components of gravity and drag in the j direction. This acceleration is approximated by the methods of Wilson (1972) and Waitt et al. (1995) as follows: (8) and (9)
where vx and vv are the ballistic velocities in the horizontal and vertical directions respectively, pa is the density of the surrounding atmosphere, v is the total velocity of the ballistic, Ab is the ballistic cross-sectional area, CD is the drag coefficient for the ballistic, M is the ballistic mass, and g is acceleration due to gravity. Calculations for M and Ab assume a spherical clast and dome rock density of 2250 kg m-3 (Sparks et al. 1998). The coefficient of drag, CD, for smooth spheres was assumed. CD varies from 0.07 to 0.20 and is a function of Reynolds number, Re, such that CD = 0.11 log(Re) - 0.55 over a range 4 x 105 < Re < 6 x 106, estimated from Achenbach (1972) and Waitt et al. (1995). The drag coefficient is poorly constrained because the shapes of the ballistics vary greatly. Others, such as Wilson (1972), Fagents & Wilson (1993), Fudali & Melson (1972) and Self et al. (1979), have assumed that irregularly shaped ballistics have CD ranging from 0.7 to >1.0. However, in a more recent analysis, Waitt et al. (1995) suggest that these high values of drag coefficient result in unrealistically high launch velocities, particularly for smaller ballistics. As theoretical justification, they offer the evidence that a dimpled golf ball achieves a greater range than a smooth one, illustrating that the rough surface of ballistics can contribute to reduced drag. In general, although roughness increases skin-friction drag, it delays flow separation and reduces form drag, which is the largest component of total drag. Equations 8 and 9 were solved for initial velocities (vx.i and vy.i) by a fourth-order Runge-Kutta numerical scheme as done by Wilson (1972) and Waitt et al. (1995). We used an average ground slope of 22°. Drag was neglected in the first 1 0 - 2 s because the initial velocity and therefore the initial drag were not known. There was no significant difference in calculated initial velocity when drag was neglected during the first 1 0 - 2 s compared to when drag was neglected during the first 10 - 4 s, indicating that neglecting drag for the first 1 0 - 2 s did not significantly affect results. Total initial ballistic velocity is referred to as Vi. Launch angle is a significant factor. The angles for any specific ballistic clasts and associated craters are unknown. Therefore we executed the method above using angles at 5C intervals, from 30° to 70°. Values of initial velocity ranged from 93 to 1 1 6 m s - 1 for a 30C launch angle and from 131 to 1 6 9 m s - 1 for a 70 0 launch angle. The optimum launch angle was between 30 and 350.Druitt et al. (2002) carried out similar calculations using the ballistic model of Self et al. (1980), which adopts higher drag coefficients than those used here. For an optimum launch angle of about 350, they estimated launch velocities of up to 1 6 0 m s - 1 . Analysis of video footage gives exit velocities of between 40 and 1 4 0 m s - 1 (our work, and Druitt et al. 2002), suggesting that the Self et al. (1980) model overestimates exit velocity and that the values calculated in this paper are probably more realistic.
Estimation of overpressure Estimation of the gas overpressure beneath the conduit cap has been approached in several ways. Robertson et al. (1998) suggested explosion pressures of 10-27MPa for the sustained explosion of
MODELLING DYNAMICS OF VULCANIAN EXPLOSIONS
327
17 September 1996. Voight et al (1999) suggest that tilt data could be used to calculate pre-explosion overpressures, but simple halfspace elastic models resulted in only upper-bound estimates of several tens of megapascals. In the same paper, an extrusion model led to an estimate of overpressure approximately 11 to 25 MPa. The conduit flow models of Melnik & Sparks (2002a) predict overpressures up to 10 MPa. A lower-bound estimate of overpressure is the rock tensile strength, which is on the order of 4 MPa (Voight et al 1999). Druitt et al. (2002) use pumice vesicularities (55-75%) to place total confining pressures between 5 and 15 MPa. The highest of these vesicularities probably represent the upper parts of the conduit immediately beneath the cap. Our assumed 20m thick cap represents less than 0.5 MPa of overburden; therefore the calculated overpressure falls between 4.5 and 14.5 MPa, which is in reasonable agreement with the aforementioned values determined by independent methods. We have therefore chosen the middle of the range of estimates, roughly 10 MPa, as a reasonable estimate of overpressure for our reference simulation. In one of our variations, we used an overpressure of 7 MPa to illustrate the effect of reduced overpressure on results.
Fragmentation description Models of magma flow unsteadiness developed by Melnik & Sparks (2002b) test the significance of different fragmentation assumptions, such as fragmentation at a fixed volume concentration of bubbles (Sparks 1978; Wilson et al 1980), at a bubble overpressure threshold in excess of the magma tensile strength (Melnik 1999), or at a critical elongation strain rate (Papale 1999). The field evidence at Montserrat includes angular platy pumice that lacks bubble elongation texture (Druitt et al 2002). These observations suggest brittle fragmentation produced by a fragmentation wave (Alidibirov & Dingwell 1996) and support the overpressure threshold fragmentation criterion of Melnik (1999). At Soufriere Hills Volcano the sudden decompression of the conduit, due to fracture and disruption of the conduit cap, occurred when the strength of the cap rock was exceeded by the overpressure in the conduit. This disruption of the cap greatly increased the pressure difference between the vesicles near the top of the conduit and the surrounding environment, which was suddenly reduced to atmospheric pressure (Fig. 5). This new pressure state exceeded the fragmentation threshold, thus initiating the fragmentation wave. The following paragraphs describe how our model represents this condition and enumerate the associated assumptions. In our numerical model, the pre-explosion conduit at time t=0 was represented by a two-phase mixture whose properties (pressure, gas mass fraction and gas volume fraction) were described as a function of depth by the equations presented above. The 30m diameter conduit was discretized in a number of cells 10 m deep and 7.5m wide and no radial variation in the initial conditions was considered. At t = t1 > 0, the pressure disequilibrium between the overpressured conduit and the confining pressure produced a decompression wave which travelled down the conduit. As the decompression wave reached a particular depth, the mixture of gas and particles at that depth began to flow out of the conduit, whereas the undisturbed portion of the mixture below the decompression wave remained in magmastatic equilibrium. Such a representation of the conduit embodies three key assumptions. We first assumed that the melt phase at any depth quenched and fragmented nearly instantaneously when the decompression wave reached it, allowing the mixture of gas, crystals and melt to be represented as a gas-particle flow. This is supported by the results of Melnik & Sparks (2002b) which indicate that fragmentation occurs immediately after the decompression wave has reached a given depth. The second assumption is that the fragmentation process involved no energy loss, which in any case was probably relatively small (1-5% according to Alidibirov 1994). Third, we assumed that no further water exsolution occurred once
Fig. 5. Schematic representation of the numerical conduit, which defines the initial conditions for explosion simulations. The conduit is subdivided into grid cells, with input parameters that include solid volume fraction, S(Z), water vapour volume fraction, g(z), diameters of particles, and total pressure of water vapour, Pg(z}. Input values are discussed in the text and summarized in Table 2. The figure on the left also represents a shock tube at time zero, representing the pre-explosion, capped conduit, and on the right at time t1, representing the condition after disruption of the caprock.
disruption of the cap and decompression began. The propagation of a fragmentation front, and acceleration of the fragmented gasparticle mixture to the surface, is too rapid for significant diffusive mass transfer to occur during Vulcanian explosions (Gardner et al. 1999; Melnik & Sparks 2002b), and thus only gas that had exsolved before the cap was disrupted participated in the explosion.
Particle sizes The sizes of solid particles chosen for our simulations (up to three sizes) have been constrained by our grain-size analyses of the pumice-and-ash flow deposits, which considered the full range of particle sizes observed in the outcrops. Our estimates of the clasts ranged from a maximum size of 50cm to a minimum size of 0.9 m. The mean grain size of the pumice-and-ash flow deposits, accounting for a 30 wt% loss of fines (
