4D simulation of explosive eruption dynamics at Vesuvius

Click Here GEOPHYSICAL RESEARCH LETTERS, VOL. 34, L04309, doi:10.1029/2006GL028597, 2007 for Full Article 4D simulation of explosive eruption dynamics at Vesuvius Augusto Neri,1 Tomaso Esposti Ongaro,1 Gianluca Menconi,1 Mattia De’Michieli Vitturi,1 Carlo Cavazzoni,2 Giovanni Erbacci,2 and Peter J. Baxter3 Received 27 October 2006; revised 23 December 2006; accepted 10 January 2007; published 24 February 2007. [1] We applied a new simulation model, based on multiphase transport laws, to describe the 4D (3D spatial coordinates plus time) dynamics of explosive eruptions. Numerical experiments, carried out on a parallel supercomputer, describe the collapse of the volcanic eruption column and the propagation of pyroclastic density currents (PDCs), for selected medium scale (subPlinian) eruptive scenarios at Vesuvius, Italy. Simulations provide crucial insights into the effects of the generation mechanism of the flows - partial collapse vs boiling-over on their evolution and hazard potential, the unstable dynamics of the fountain, and the influence of Mount Somma on the propagation of PDCs into the circumVesuvian area, one of the world’s most hazardous volcanic settings. Results also show that it is possible to characterize the volcanic column behavior in terms of percentage of the mass of pyroclasts collapsed to the ground and how this parameter strongly influences the dynamics and hazard of the associated PDCs. Citation: Neri, A., T. Esposti Ongaro, G. Menconi, M. De’Michieli Vitturi, C. Cavazzoni, G. Erbacci, and P. J. Baxter (2007), 4D simulation of explosive eruption dynamics at Vesuvius, Geophys. Res. Lett., 34, L04309, doi:10.1029/2006GL028597. 1. Introduction [2] A transient 3D numerical simulation model is an essential advance in increasing our understanding of the dynamics of PDCs in explosive eruptions, since they are the most destructive and deadliest of volcanic phenomena [Simkin et al., 2001], as attested by the thousands of deaths and devastation caused in 79AD and 1631 by the largest eruptions of Vesuvius in its last two millennia of activity. Explosive eruptions are characterized by the injection of fragmented magma into the atmosphere where it disperses as a multi-particle and multi-component mixture of gases and pyroclasts in a process that can occur on a wide range of spatial and temporal scales as well as with different mechanisms. In Plinian and sub-Plinian eruptions the buoyant mixture lofts in an erupting column to tens of kilometers in height, but a loss of efficiency of atmospheric mixing or a structural collapse of the volcanic edifice can precipitate a partial or total collapse of the column to form PDCs, or pyroclastic flows (s.l.), in which the hot mixture propagates at high velocity along the volcano flanks [Sparks et al., 1 Istituto Nazionale di Geofisica e Vulcanologia, Section of Pisa, Pisa, Italy. 2 High Performance Computing Group, Consorzio Interuniversitario per il Calcolo Automatico dell’Italia Nord Orientale, Bologna, Italy. 3 Institute of Public Health, University of Cambridge, Cambridge, UK. Copyright 2007 by the American Geophysical Union. 0094-8276/07/2006GL028597$05.00 1997] with potentially devastating consequences. Although several different types of PDCs have been recognized, in most volcanological models they comprise a basal concentrated dense flow underlying a more dilute surge-like ash cloud of fine particles [Druitt, 1998; Branney and Kokelaar, 2002]. [3] In recent decades, numerous authors have used different numerical models based both on 1D steady-state homogeneous flow [Sparks et al., 1978; Malin and Sheridan, 1982; Bursik and Woods, 1996; Dellino et al., 2004] and 2D transient homo/multiphase flow assumptions [Wohletz et al., 1984; Valentine and Wohletz, 1989; Dobran et al., 1993; Neri et al., 2003; Dartevelle et al., 2004; Ishimine, 2005] to describe the mechanics of column collapse and to capture the dynamics and hazards of PDCs. Although these models elucidate several of the first-order characteristics of explosive eruptions, many of them ignore the description of the column collapse and all neglect crucial 3D features of the system, such as the volcano topography, flow instabilities, and atmospheric conditions, making their application to real examples problematic. Only recently has the availability of supercomputers allowed the analysis of some 3D features of buoyant plumes by adopting homogeneous flow models [Oberhuber et al., 1998; Suzuki et al., 2005]. [4] Here we present, for the first time, 4D multiphase flow simulations of the formation of PDCs by the collapse of the volcanic column and their propagation over the actual volcano topography. The model, named PDAC (Pyroclastic Dispersal Analysis Code) (T. Esposti Ongaro et al., A parallel multiphase flow code for the 3D simulation of explosive volcanic eruptions, submitted to Parallel Computing, 2006, hereinafter referred to as Esposti Ongaro et al., submitted manuscript, 2006), represents the development of the PDAC2D model [Neri et al., 2003] already used for the analysis of the 2D dynamics of these phenomena [Dobran et al., 1994; Clarke et al., 2002; Todesco et al., 2002]. The first successful application of PDAC was carried out in the reconstruction of the main 3D features of the PDC in the Boxing Day volcanic blast that occurred at the Soufrière Hills volcano, Montserrat in 1997 [Esposti Ongaro et al., 2005]. In particular, the model reproduced the reconstructed distribution of products and damages according to independent observations, when the simulation incorporated the available data on the lava dome’s physical characteristics. In the following, we present the application of the code to explosive events at Vesuvius. [5] It is worth mentioning that the initial and boundary conditions describing the eruptive scenarios reported below are only meaningful in the context of a parametric study and that major epistemic and aleatoric uncertainties affect several of them [Sparks and Aspinall, 2004]. Thus the simu- L04309 1 of 7 L04309 NERI ET AL.: 4D EXPLOSIVE ERUPTION DYNAMICS AT VESUVIUS lated scenarios presented here need to be interpreted as ideal representations of likely events at Vesuvius. 2. Model Formulation and Numerics [6] In the applied model the fundamental processes governing explosive eruptions are described in terms of first principles whereas the constitutive equations describe the physical properties of the eruptive fluids [Valentine and Wohletz, 1989; Dobran et al., 1993; Gidaspow, 1994; Dartevelle et al., 2004]. The PDAC model describes the Eulerian transport equations of mass, momentum and enthalpy for a multiphase mixture formed by a continuous multi-component gas phase and n solid particulate phases representative of magma fragments (such as ash, crystals and lapilli). In more detail, the atmospheric dispersal of pyroclasts is computed over the 3D spatial domain and over time by solving a set of mass, momentum, and energy balance equations for each considered particulate phase. The model represents each pyroclast category (ash, crystals, etc.) by a solid particulate phase characterized by a diameter, density, specific heat, thermal conductivity, etc. and provides as output its concentration, velocity, temperature over the 3D domain at each instant of time. For simplicity, we assumed no change in particle size due to the effects of secondary fragmentation or aggregation processes. Pyroclast sedimentation and elutriation, as well as the air entrainment and heating, are explicitly computed through the solution of the non-equilibrium multiphase flow equations. The deposition process (implying a particle loss from the flow bottom) was not taken into account in the present formulation. [7] The solution of the model equations is based on a second-order accurate, finite-volume discretization scheme and a semi-implicit time-advancing scheme. The new 3D code has been parallelized following an approach based on the domain-decomposition strategy and the MPI protocol (Esposti Ongaro et al., submitted manuscript, 2006). According to this, the physical domain was subdivided in a number of smaller sub-domains, with each one assigned to a single processor of the parallel supercomputer. A novel immersed-boundary technique suited to compressible multiphase flows has been adopted in some of the simulations in order to accurately describe the no-slip flow condition at the interface between the fluid flow and the irregular 3D topography, even with relatively coarse meshes (M. De’Michieli Vitturi et al., An immersed boundary method for compressible multiphase flow: Application to the dynamics of pyroclastic density currents, submitted to Computational Geosciences, 2006). [8] The transport equations were solved on a Cartesian non-uniform structured mesh, with minimum horizontal and vertical spacing of 20 m and a maximum cell size of 100m. The topography was acquired from an accurate 10 mresolution Digital Elevation Model (DEM). The typical time-step was of 0.01 s. The adopted domain extended 12 and 14 Km in the two planar directions and up to 8 km in the vertical. The simulations required about 14,000 hrs of total CPU time per 1,000 s of simulation running on the IBM SP5 (512 Power 5 processors at 1.9 GHz) available at CINECA (Bologna, Italy). [9] Sensitivity tests were performed in order to assess the influence of several numerical parameters, including grid size, time step, domain extension, boundary conditions and L04309 numerical accuracy (Esposti Ongaro et al., submitted manuscript, 2006). The results showed that the investigated numerical parameters do not significantly affect the main simulation outputs (fountain height, column behavior, percentage of mass collapsed, flow propagation, etc.). On the other hand, the adopted vertical mesh resolution along the terrain prevents an accurate description of the basal coarserich granular flow underlying the more dilute fine-rich ash cloud, thus introducing a bias in the estimate of the absolute values of flow runout and hazardous variables. Nevertheless, parametric studies performed with different grid resolutions have shown that the propagation of the more dilute surge-like component of the flow is weakly affected by this effect [Neri et al., 2003; this study]. 3. Scenarios Investigated [10] Vesuvius, a strato-volcano in central Italy, has a history of alternating explosive and effusive periods of activity. Since the 79 AD eruption, several other catastrophic eruptions have occurred with one of the largest being the sub-Plinian event in 1631, when over 4,000 people were killed and extensive damage was caused by pyroclastic flows in the surrounding towns [Rosi et al., 1993]. A cessation of activity since 1944 has created conditions that may favour a further sub-Plinian event in the next eruption, when today about 550,000 people live in the designated evacuation area (Red Zone) around the volcano. The dynamics of the column collapse, the mechanisms of generation of pyroclastic flows, as well as the influence of the volcano’s topography - particularly that of Mount Somma, a caldera rim about 200 m high on the N flank - on the propagation of the flows has until now been only poorly understood. [11] Several scenarios of collapsing column and pyroclastic flow propagation along the Vesuvius flanks were simulated with reference to a sub-Plinian scale event, similar to the Verdoline (16 ka BP), 472 and 1631 AD sub-Plinian eruptions, with an estimated intensity range of 2– 8  107 kg/s [Rosi et al., 1993; Cioni et al., 2003; Sulpizio et al., 2005]. [12] Here we present two distinct simulations representative of two different collapsing mechanisms: (1) SimA, characterized by the partial collapse of the volcanic column, and (2) SimB, characterized by a total collapse or ‘‘boilingover’’ type injection of the mixture in the atmosphere. For both cases, flow conditions at the crater exit were assumed constant in time, pressure-balanced, and the intensity was assumed equal to 5  107 kg/s. Therefore the different dynamics of the columns are due uniquely to the different vent conditions assumed which, in turn, affect the regime of the volcanic fountain. SimA adopts a crater diameter of 250 m, exit velocity of 175 m/s, and a mixture density of 6.5 kg/m3. These conditions are representative of a quasi steady-state sustained column and were computed by simulating the decompression of the mixture in a 2D axisymmetric crater by assuming vent conditions (base of crater) derived by magma ascent modelling using magmatic properties typical of Vesuvius [Todesco et al., 2002]. In contrast, SimB adopts a crater diameter of 350 m (i.e. the crater surface area is about twice that of SimA), exit velocity of 90 m/s and the same exit density of SimA; such 2 of 7 L04309 NERI ET AL.: 4D EXPLOSIVE ERUPTION DYNAMICS AT VESUVIUS L04309 Figure 1. Large-scale dynamics of the partial column collapse scenario at Vesuvius (SimA). The two color isosurfaces represent the total volume particle concentration in the plume and PDCs and correspond to 10 4 (dark brown) and 10 6 (light brown). (a) 30, (b) 300, (c) 600, (d) 900, and (e) 1300 s from the injection of the mixture into the atmosphere taken from a SW view. (f) 1300 s from the injection of the mixture into the atmosphere taken from N. The topographic relief of Mount Somma, on the N slope of the volcano, is clearly visible, as is its major influence on the propagation of the flows that are unable to overcome it. crater conditions were assumed representative of those produced during vent disruption and widening, as supposed to have occurred during the 1631 eruption [Rosi et al., 1993]. In both simulations, a water content of 2 wt% and magma temperature of 950°C were adopted [Todesco et al., 2002] to favour collapsing conditions of the column, although similar dynamics could be reproduced by using different combinations of vent flow conditions. A threephase formulation of the mixture, with two particle size classes of 30 and 500 micron and a two-component gas phase - formed of atmospheric air and water vapour leaving the vent - was adopted based on available grain-size data [Rosi et al., 1993; Cioni et al., 2003; Sulpizio et al., 2005] and previous studies [Todesco et al., 2002]. Particles of 30 and 500 micron were assumed equally present in weight, with densities of 2,800 and 1,000 kg/m3, respectively, and were assumed representative of the fine ash and coarse pumiceous pyroclasts of the mixture. Although a more realistic description of the grain-size distribution could be implemented in the model [Clarke et al., 2002; Neri et al., 2003] the three-phase formulation here adopted represents a good trade-off between an accurate description of the column dynamics and the computational resources required. 4. Results [13] Figure 1 illustrates the large-scale dynamics of SimA through a series of snapshots of ash concentration at different times from the injection of the mixture into the atmosphere. At about 30 s (Figure 1a) the gas-particle jet has lost its vertical momentum and the head of plume is close to collapse at about 2000 m above the vent. At 300 s (Figure 1b) the column has collapsed to the ground and pyroclastic flows are beginning to form all around the fountain. Simultaneously, a number of thermals above the PDCs are pulled into the rising central plume. By 600 s (Figure 1c) three main pyroclastic flows are moving in the W, SW, and SE directions. From inspection of results (see 3 of 7 L04309 NERI ET AL.: 4D EXPLOSIVE ERUPTION DYNAMICS AT VESUVIUS L04309 Figure 2. Dynamics of the fountain of the partial column collapse scenario (SimA). (a) 600, (b) 630, and (c) 660 s are taken from SW and show the evolution of the isosurface of the total particle volume concentration corresponding to 10 3. (d) The oscillations of the vertical velocity of the flow at various heights along the fountain axis. The largest oscillations are predicted close to the fountain top, with dominant frequencies in the range 0.02– 0.05 Hz. Animation S1 in the auxiliary material1 for details) it is evident how the W and S flows are produced by the diversion of the flow path caused by Mount Somma. The SW flow is the most advanced with a front width of about 500 m and, at this time, has reached a distance of about 3.5 km from the vent. By 900 s (Figure 1d) the three main flows start producing secondary branches which further propagated down slope with the most advanced reaching a distance of about 5 km from the vent. At 1300 s (Figure 1e) the flows moving in the SW direction have already reached the coastline of the Tyrrhenian sea whereas those produced by the diversion of Mount Somma are slowly moving in proximity to the border of the Red Zone area, at about 8 km from the vent. As discussed above a precise determination of the flow-runout is not possible because the density stratification of the flow occurs on a length scale comparable to the vertical grid resolution (20 m, in this case) and roughness length-scale (about 10– 30 m in the urbanized areas). However, from the distribution computed at 1500 s (not shown), it appears that large sectors of the south flank of the volcano would be affected by surge-like PDCs in a scenario of this type. Of great importance for hazard assessment is the finding that the flows were unable to significantly overcome the Mount Somma rim despite the remarkable collapse height, as is clearly visible in Figure 1f (N view). Only minor, dilute, low-temperature, short-runout flows were able to overcome Mount Somma mostly in the W direction where the rim altitude is lower. 1 Auxiliary materials are available in the HTML. doi:10.1029/ 2006GL028597. [14] In addition, we have made an attempt at quantifying the hazardous impacts of the flows using the simulation outcomes. Dynamic pressure and temperature values inferred from the simulation are below 1 kPa and 100°C, respectively, for distances from the vent greater than about 4 km, although dynamic pressures were corrected by a conservative factor of five to account for the effect of the vertical mesh size on the value of the mixture density in the basal cells, as deduced by 2D simulations performed at higher resolution. At these values the PDCs would cause relatively minor to moderate damage in the impacted populated area [Spence et al., 2004; Baxter et al., 2005] and they lend further weight to our view that eventually 3D modeling could be used in mitigation methodologies. [15] The outcomes of SimA also clearly show the strongly unsteady and chaotic behavior of the collapsing column that feeds the flows during a partial collapse (see Animation S2 in the auxiliary material). For example, Figure 2 illustrates three subsequent snapshots of the fountain ash concentration at 600 (Figure 2a), 630 (Figure 2b), and 660 (Figure 2c) s. In this case, the depicted isosurface displays the inner portion of the fountain including the single batches of mixture that intermittently and anisotropically leave the fountain and feed the pyroclastic flows. The transient dynamics produce main oscillations of the fountain height – defined as the quota where the average vertical velocity is nil – with variations up to 50%. Analysis of the oscillations of the vertical mixture velocity along the jet axis (Figure 2d), indicates that their amplitude and frequency depend mainly on the region of the fountain investigated. The largest oscillations are predicted close to the fountain top, with dominant frequencies in the range 0.02–0.05 Hz. Finally, integration of raw data within 4 of 7 L04309 NERI ET AL.: 4D EXPLOSIVE ERUPTION DYNAMICS AT VESUVIUS L04309 Figure 3. Comparison between (a and b) partial collapse (SimA) and (c and d) boiling-over (SimB) scenarios at Vesuvius, 500 s (Figures 3a and 3c) and 1000 s (Figures 3b and 3d) from the injection of the mixture in the atmosphere. Colorcontours on the 2D vertical plane represent the total volume particle concentration in the plume and PDCs for the two simulations. Colors on the topography indicate the computed total volume concentration of particles in the first cell along the ground. the fountain indicates that about 50% of both particle sizes ejected from the vent collapse to the ground. In contrast, at the end of the simulation, less than 10 wt% of 500 mm particles and 5 wt% of 30 mm particles ejected from the vent form the density currents, reflecting an effective elutriation by natural convection of both particle types during the proximal-median propagation of the currents [Druitt, 1998; Branney and Kokelaar, 2002; Neri et al., 2003]. [16] Figure 3 illustrates the boiling-over scenario (SimB) through the evolution of the total particle volume fraction in the atmosphere and in comparison to the partial collapse scenario (SimA) (see Animation S3 in the auxiliary material). In detail, Figure 3 shows the distribution of particles on the ground and along a vertical plane passing through the fountain axis and oriented in the SW direction for SimA (Figures 3a and 3b) and SimB (Figures 3c and 3d) at 500 and 1000 s, respectively. At 30 s the volcanic jet has lost its momentum at about 450 m above the crater and the collapsing stream of the fountain has already reached the ground. By 90 s, pyroclastic flows are spreading almost symmetrically and a dilute, rising plume starts to develop above the fountain. At 500 s (Figure 3c) a main pyroclastic flow is propagating in the SW direction, while two more flows propagate in the W and SE directions again due to the diversionary effect of Mount Somma. Similarly, Figure 3d shows the particle concentration at the ground and over the 2D plane at 1000 s. At this time, all the main flows have reached a distance of about 8 km and are entering the Tyrrhenian Sea or leaving the Red Zone area. It is worth noting the control of topography on the denser part of the flows, as shown by the red color along the ground, as well as the most dispersed pattern of the dilute portions of the flow (green color). In comparison to SimA, SimB shows much faster dynamics as well as more critical conditions associated with the flows. Thus, dynamic pressures of the order of 1 – 3 kPa and temperatures over 250°C are predicted for SimB at distances of about 7.5 km from the crater, and so much more extensive devastation of the Red Zone would be likely, which can be attributed to the greater mass collapsed from the column. Indeed about 90 wt% of the mass erupted has collapsed to the ground, but only about 20% forms the flows at the end of the simulation (see the effective elutriation of particles from the upper part of the 5 of 7 L04309 NERI ET AL.: 4D EXPLOSIVE ERUPTION DYNAMICS AT VESUVIUS flows in Figures 3c and 3d). Nevertheless, for both SimA and SimB, Mount Somma presents an insurmountable barrier for the flows that are, in effect, redirected in the W and SE directions. This effect of Mount Somma was also observed in another simulation performed with an intensity of 8  107 kg/s and similar crater conditions (see Animation S4 in the auxiliary material). [17] Due to the uncertainty associated with the reconstruction of past explosive events at Vesuvius a rigourous validation of the above described modeling results is not possible. However, it is notable that the aerial extension of the regions affected by the flows during the 1631 eruption is comparable to that predicted by the model and also that Mount Somma relief was not overcome by the flows in that event [Rosi et al., 1993]. Some discrepancies with the evidence from other sub-Plinian eruptions can be explained in terms of different vent location and topography with respect to the values assumed in the present study [Cioni et al., 2003; Sulpizio et al., 2005]. 5. Conclusions [18] Our study shows that 4D multiphase numerical models can illuminate the non-intuitive and internal dynamics of explosive eruptions that cannot otherwise be studied by direct observation or using previous models. In particular, from the simulation performed, it is clear that the collapse mechanism - partial collapse vs boiling-over strongly controls the dynamics, evolution, and hazard potential of the associated pyroclastic flows. In the range of vent conditions investigated, the different outcomes are attributable to the amount, or percentage, of erupted mass that collapse to the ground rather than to the collapse height of the column, suggesting that, at constant any other parameters, low fountaining boiling-over events likely represent the most hazardous scenarios. Simulation results also indicate that, for the eruptive scale and vent conditions investigated, Mount Somma represents an almost insurmountable barrier for the flows. [19] Further studies using different eruptive inputs (e.g., vent location, crater shape, vent and atmospheric conditions) should lead to an improved understanding of the column collapse and generation and propagation of PDCs at Vesuvius in determining spatial risk as well as for enumerating their most hazardous characteristics in future eruptions. Similarly, investigations aimed at model validation need to be further developed by comparing model results to ad hoc lab experiments and well-known previous eruptive events. This work is another step towards the long-term goal of developing a multi-disciplinary approach to disaster reduction based on simulation modeling of the main eruption hazards and quantifying their potential impacts in highrisk volcanic regions. [20] Acknowledgments. This work was funded by the EU-project EXPLORIS and by the Dipartimento di Protezione Civile, Italy. We warmly thank our colleagues of the EXPLORIS project for their assistance and contributions during the carrying out of this work. Maria Teresa Pareschi and Marina Bisson (INGV Pisa) and Roberto Gori, Tiziano Diamanti, and Antonella Guidazzoli (CINECA Bologna) are thanked for providing the DEM of Vesuvius as well as for contributing to the development of the graphics tools used to produce Figure S4 and Animation S1 of the auxiliary material. L04309 References Baxter, P. J., R. Boyle, P. Cole, A. Neri, R. Spence, and G. Zuccaro (2005), The impact of pyroclastic surges on buildings at the eruption of the Soufrière Hills volcano, Montserrat, Bull. Volcanol., 67, 292 – 317. Branney, M. J., and P. Kokelaar (2002), Pyroclastic Density Currents and the Sedimentation of Ignimbrites, Geol. Soc. Mem., vol. 27, 143 pp., Geol. Soc., London. Bursik, M. I., and A. W. Woods (1996), The dynamics and thermodynamics of large ash flows, Bull. Volcanol., 58, 175 – 193. Cioni, R., R. Sulpizio, and N. Garruccio (2003), Variability of the eruption dynamics during a sub-Plinian event: The Greenish Pumice eruption of Somma-Vesuvius (Italy), J. Volcanol. Geotherm. Res., 124, 89 – 114. Clarke, A. B., B. Voight, A. Neri, and G. Macedonio (2002), Transient dynamics of Vulcanian explosions and column collapse, Nature, 415, 897 – 901. Dartevelle, S., W. I. Rose, J. Stix, K. Kelfoun, and J. W. Vallance (2004), Numerical modeling of geophysical granular flows: 2. Computer simulations of plinian clouds and pyroclastic flows and surges, Geochem. Geophys. Geosyst., 5, Q08004, doi:10.1029/2003GC000637. Dellino, P., R. Isaia, and M. Veneruso (2004), Turbulent boundary layer shear flows as an approximation of base surges at Campi Flegrei (southern Italy), J. Volcanol. Geotherm. Res., 133, 211 – 228. Dobran, F., A. Neri, and G. Macedonio (1993), Numerical simulation of collapsing volcanic columns, J. Geophys. Res., 98, 4231 – 4259. Dobran, F., A. Neri, and M. Todesco (1994), Assessing pyroclastic flow hazard at Vesuvius, Nature, 367, 551 – 554. Druitt, T. H. (1998), Pyroclastic density currents, in The Physics of Explosive Volcanic Eruptions, edited by J. S. Gilbert and R. S. J. Sparks, Geol. Soc. Spec. Publ. London, 145, 27 – 50. Esposti Ongaro, T., A. Clarke, A. Neri, B. Voight, and C. Widiwijayanti (2005), A new high-performance 3D multiphase flow code to simulate volcanic blasts and pyroclastic density currents: Example from the Boxing Day event, Montserrat, Eos Trans. AGU, 86(52), Fall Meet. Suppl., Abstract V31D-0645. Gidaspow, D. (1994), Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions, 467 pp., Elsevier, New York. Ishimine, Y. (2005), A numerical study of pyroclastic surges, J. Volcanol. Geotherm. Res., 139, 33 – 57. Malin, M. C., and M. F. Sheridan (1982), Computer-assisted mapping of pyroclastic surges, Science, 217, 637 – 640. Neri, A., T. Esposti Ongaro, G. Macedonio, and D. Gidaspow (2003), Multiparticle simulation of collapsing volcanic columns and pyroclastic flows, J. Geophys. Res., 108(B4), 2202, doi:10.1029/2001JB000508. Oberhuber, J. M., M. Herzog, H. F. Graf, and K. Schwanke (1998), Volcanic plumes simulation on large scales, J. Volcanol. Geotherm. Res., 87, 29 – 53. Rosi, M., C. Principe, and R. Vecci (1993), The 1631 eruption of Vesuvius reconstructed from the review of chronicles and study of deposits, J. Volcanol. Geotherm. Res., 58, 151 – 182. Simkin, T., L. Siebert, and R. Blong (2001), Volcano fatalities: Lessons from the historical record, Science, 291, 255. Sparks, R. S. J., and W. P. Aspinall (2004), Volcanic activity: Frontiers and challenges in forecasting, prediction and risk assessment, in The State of the Planet: Frontiers and Challenges in Geophysics, Geophys. Monogr. Ser., vol. 150, edited by R. S. J. Sparks, and C. J. Hawkesworth, pp. 359 – 373, AGU, Washington, D. C. Sparks, R. S. J., L. Wilson, and G. Hulme (1978), Theoretical modeling of the generation, movement, and emplacement of pyroclastic flows by column collapse, J. Geophys. Res., 83, 1727 – 1739. Sparks, R. S. J., et al. (1997), Volcanic Plumes, 557 pp., John Wiley, Hoboken, N. J. Spence, R. J. S., G. Zuccaro, S. Petrazzuoli, and P. J. Baxter (2004), The resistance of buildings to pyroclastic flows: Analytical and experimental studies and their application to Vesuvius, Nat. Haz. Rev., 5, 48 – 59. Sulpizio, R., D. Mele, P. Dellino, and L. La Volpe (2005), A complex Subplinian-type eruption from low-viscosity, phonolitic to tephri-phonolitic magma: The AD 472 (Pollena) eruption of Somma-Vesuvius, Italy, Bull. Volcanol., 67, 743 – 767. Suzuki, Y. J., T. Koyaguchi, M. Ogawa, and I. Hachisu (2005), A numerical study of turbulent mixing in eruption clouds using three-dimensional fluid dynamics model, J. Geophys. Res., 110, B08201, doi:10.1029/ 2004JB003460. Todesco, M., A. Neri, T. Esposti Ongaro, P. Papale, G. Macedonio, R. Santacroce, and A. Longo (2002), Pyroclastic flow hazard assessment at Vesuvius (Italy) by using numerical modelling: I. Large-scale dynamics, Bull. Volcanol., 64, 155 – 177. Valentine, G. A., and K. H. Wohletz (1989), Numerical models of Plinian eruption columns and pyroclastic flows, J. Geophys. Res., 94, 1867 – 1887. 6 of 7 L04309 NERI ET AL.: 4D EXPLOSIVE ERUPTION DYNAMICS AT VESUVIUS Wohletz, K. H., T. R. McGetchin, M. T. Sandford, and E. M. Jones (1984), Hydrodynamic aspects of caldera-forming eruptions: Numerical models, J. Geophys. Res., 89, 8269 – 8285. P. J. Baxter, Institute of Public Health, University of Cambridge, Robinson Way, Cambridge CB2 2SR, UK. L04309 C. Cavazzoni and G. Erbacci, High Performance Computing Group, Consorzio Interuniversitario per il Calcolo Automatico dell’Italia Nord Orientale, Casalecchio di Reno, Bologna I-040033, Italy. A. Neri, T. Esposti Ongaro, G. Menconi, and M. De’Michieli Vitturi, Istituto Nazionale di Geofisica e Vulcanologia, Section of Pisa, via della Faggiola 32, Pisa I-56126, Italy. ([email protected]) 7 of 7