A semi-implicit, second-order-accurate numerical model for multiphase underexpanded volcanic jets

Abstract. An improved version of the PDAC (Pyroclastic Dispersal Analysis Code, Esposti Ongaro et al., 2007) numerical model for the simulation of multiphase volcanic flows is presented and validated for the simulation of multiphase volcanic jets in supersonic regimes. The present version of PDAC includes second-order time- and space discretizations and fully multidimensional advection discretizations in order to reduce numerical diffusion and enhance the accuracy of the original model. The model is tested on the problem of jet decompression in both two and three dimensions. For homogeneous jets, numerical results are consistent with experimental results at the laboratory scale (Lewis and Carlson, 1964). For nonequilibrium gas–particle jets, we consider monodisperse and bidisperse mixtures, and we quantify nonequilibrium effects in terms of the ratio between the particle relaxation time and a characteristic jet timescale. For coarse particles and low particle load, numerical simulations well reproduce laboratory experiments and numerical simulations carried out with an Eulerian–Lagrangian model (Sommerfeld, 1993). At the volcanic scale, we consider steady-state conditions associated with the development of Vulcanian and sub-Plinian eruptions. For the finest particles produced in these regimes, we demonstrate that the solid phase is in mechanical and thermal equilibrium with the gas phase and that the jet decompression structure is well described by a pseudogas model (Ogden et al., 2008). Coarse particles, on the other hand, display significant nonequilibrium effects, which associated with their larger relaxation time. Deviations from the equilibrium regime, with maximum velocity and temperature differences on the order of 150 m s−1 and 80 K across shock waves, occur especially during the rapid acceleration phases, and are able to modify substantially the jet dynamics with respect to the homogeneous case

Il progetto EPLORIS: La ricostruzione virtuale dell'eruzione del Vesuvio

The main objective of the Exploris project consists in the quantitative analysis of explosive eruption risk in densely populated EU volcanic regions and the evaluation of the likely effectiveness of possible mitigation measures through the development of volcanic risk facilities (such as supercomputer models, vulnerability databases, and probabilistic risk assessment protocols) and their application to high-risk European volcanoes. Exploris’ main ambition is to make a significant step forward in the assessment of explosive eruption risk in highly populated EU cities and islands. For this project, a new simulation model, based on fundamental transport laws to describe the 4D (3D spatial co-ordinates plus time) multiphase flow dynamics of explosive eruptions has been developed and parallelized in INGV and CINECA. Moreover, CINECA developed specific tools to efficiently visualise the results of simulations. This article presents the results of the large numerical simulations, carred out with CINECA’s Supercomputers, to describe the collapse of the volcanic eruption column and the propagation of pyroclastic density currents, for selected medium scale (sub-Plinian) eruptive scenarios at Vesuvius

Insights into the formation and dynamics of coignimbrite plumes from one-dimensional models

Coignimbrite plumes provide a common and effective mechanism by which large volumes of fine-grained ash are injected into the atmosphere. Nevertheless, controls on formation of these plumes as a function of eruptive conditions are still poorly constrained. Herein, two 1-D axysymmetric steady state models were coupled, the first describing the parent pyroclastic density current and the second describing plume rise. Global sensitivity analysis is applied to investigate controls on coignimbrite plume formation and describe coignimbrite source and the maximum plume height attained. For a range of initial mass flow rates between 108 and 1010 kg/s, modeled liftoff distance (the distance at which neutral buoyancy is attained), assuming radial supercritical flow, is controlled by the initial flow radius, gas mass fraction, flow thickness, and temperature. The predicted decrease in median grain size between flow initiation and plume liftoff is negligible. Calculated initial plume vertical velocities, assuming uniform liftoff velocity over the pyroclastic density current invasion area, are much greater (several tens of m/s) than those previously used in modeling coignimbrite plumes (1 m/s). Such velocities are inconsistent with the fine grain size of particles lofted into coignimbrite plumes, highlighting an unavailability of large clasts, possibly due to particle segregation within the flow, prior to plume formation. Source radius and initial vertical velocity have the largest effect on maximum plume height, closely followed by initial temperature. Modeled plume heights are between 25 and 47 km, comparable with Plinian eruption columns, highlighting the potential of such events for distributing fine-grained ash over significant areas

Large-eddy simulation of pyroclastic density currents

We investigate the dynamics of turbulent pyroclastic density currents (PDCs) by adopting a 3D, Eulerian-Eulerian multiphase flow model, in which solid particles are treated as a continuum and the grain-size distribution is simplified by assuming two particulate phases. The turbulent sub-grid stress of the gas phase is modelled within the framework of Large-Eddy Simulation (LES) by means of a eddy-viscosity model together with a wall closure. Despite the significant numerical diffusion associated to the upwind method adopted for the Finite-Volume discretization, numerical simulations demonstrate the need of adopting a Sub-Grid Scale (SGS) model, while revealing the complex interplay between the grid and the SGS filter sizes. We also analyse the relationship between the averaged flow dynamic pressure and the action exerted by the PDC on a cubic obstacle, to evaluate the impact of a PDC on a building. Numerical results suggest that the average flow dynamic pressure can be used as a proxy for the force per unit surface acting on the building envelope (Fig. 5), even for such steeply stratified flows. However, it is not possible to express such roportionality as a constant coefficient such as the drag coefficient in a steady-state current. The present results indeed indicate that the large epistemic and aleatory uncertainty on initial and boundary conditions has an impact on the numerical redictions which is comparable to that of grid resolution

An application of parallel computing to the simulation of volcanic eruptions

A parallel code for the simulation of the transient 3D dispersal of volcanic particles produced by explosive eruptions is presented. The model transport equations, based on the multiphase flow theory, describe the atmospheric dynamics of the gas-particle mixture ejected through the volcanic crater. The numerics is based on a finite-volume discretization scheme and a pressure-based iterative non-linear solver suited to compressible multiphase flows. The code has been parallelized by adopting an ad hoc domain partitioning scheme that enforces the load balancing. An optimized communication layer has been built over the Message-Passing Interface. The code proved to be remarkably efficient on several high-performance platforms and makes it possible to simulate fully 3D eruptive scenarios on realistic volcano topography

Ash plume properties retrieved from infrared images: a forward and inverse modeling approach

We present a coupled fluid-dynamic and electromagnetic model for volcanic ash plumes. In a forward approach, the model is able to simulate the plume dynamics from prescribed input flow conditions and generate the corresponding synthetic thermal infrared (TIR) image, allowing a comparison with field-based observations. An inversion procedure is then developed to retrieve ash plume properties from TIR images. The adopted fluid-dynamic model is based on a one-dimensional, stationary description of a self-similar (top-hat) turbulent plume, for which an asymptotic analytical solution is obtained. The electromagnetic emission/absorption model is based on the Schwarzschild's equation and on Mie's theory for disperse particles, assuming that particles are coarser than the radiation wavelength and neglecting scattering. In the inversion procedure, model parameters space is sampled to find the optimal set of input conditions which minimizes the difference between the experimental and the synthetic image. Two complementary methods are discussed: the first is based on a fully two-dimensional fit of the TIR image, while the second only inverts axial data. Due to the top-hat assumption (which overestimates density and temperature at the plume margins), the one-dimensional fit results to be more accurate. However, it cannot be used to estimate the average plume opening angle. Therefore, the entrainment coefficient can only be derived from the two-dimensional fit. Application of the inversion procedure to an ash plume at Santiaguito volcano (Guatemala) has allowed us to retrieve the main plume input parameters, namely the initial radius $b_0$, velocity $U_0$, temperature $T_0$, gas mass ratio $n_0$, entrainment coefficient $k$ and their related uncertainty. Moreover, coupling with the electromagnetic model, we have been able to obtain a reliable estimate of the equivalent Sauter diameter $d_s$ of the total particle size distribution. The presented method is general and, in principle, can be applied to the spatial distribution of particle concentration and temperature obtained by any fluid-dynamic model, either integral or multidimensional, stationary or time-dependent, single or multiphase. The method discussed here is fast and robust, thus indicating potential for applications to real-time estimation of ash mass flux and particle size distribution, which is crucial for model-based forecasts of the volcanic ash dispersal process

Pyroclastic flow dynamics and hazard in a caldera setting: application to Phlegrean Fields

Numerical simulation of pyroclastic density currents has developed significantly in recent years and is increasingly applied to volcanological research. Results from physical modeling are commonly taken into account in volcanic hazard assessment and in the definition of hazard mitigation strategies. In this work, we modeled pyroclastic density currents in the Phlegrean Fields caldera, where flows propagating along the flat ground could be confined by the old crater rims that separate downtown Naples from the caldera. The different eruptive scenarios (mass eruption rates, magma compositions, and water contents) were based on available knowledge of this volcanic system, and appropriate vent conditions were calculated for each scenario. Simulations were performed along different topographic profiles to evaluate the effects of topographic barriers on flow propagation. Simulations highlighted interesting features associated with the presence of obstacles such as the development of backflows. Complex interaction between outward moving fronts and backflows can affect flow propagation; if backflows reach the vent, they can even interfere with fountain dynamics and induce a more collapsing behavior. Results show that in the case of large events ( 108 kg/s), obstacles affect flow propagation by reducing flow velocity and hence dynamic pressure in distal regions, but they cannot stop the advancement of flows. Deadly conditions (in terms of temperature and ash concentration) characterize the entire region invaded by pyroclastic flows. In the case of small events (2.5 107 kg/s), flows are confined by distal topographic barriers which provide valuable protection to the region beyond

Fluid-dynamics of the 1997 Boxing Day volcanic blast on Montserrat, W.I.

Directed volcanic blasts are powerful explosions with a significant laterally¬directed component, which can generate devastating, high-energy pyroclastic density currents (PDCs). Such blasts are an important class of eruptive phenomena, but quantified understanding of their dynamics and effects is still incomplete. Here we use 2D and 3D multiparticle thermofluid dynamic flow codes to examine a powerful volcanic blast that occurred on Montserrat in December 1997. Based on the simulations, we divide the blast into three phases; an initial burst phase lasts roughly 5 s and involves rapid expansion of the gas-pyroclast mixture, a gravitational collapse phase which occurs when the erupted material fails to mix with sufficient air to form a buoyant column and thus collapses asymmetrically, and a PDC phase which is dominated by motion parallel to the ground surface and is influenced by topography. We vary key input parameters such as total gas energy and total solid mass to understand their influence on simulations, and compare the simulations with independent field observations of damage and deposits, demonstrating that the models generally capture important large-scale features of the natural phenomenon. We also examine the 2D and 3D model results to estimate the flow Mach number and conclude that the range of damage sustained at villages on Montserrat can be reasonably explained by the spatial and temporal distribution of the dynamic pressure associated with subsonic PDCs

Inter-comparison of three-dimensional models of volcanic plumes

We performed an inter-comparison study of three-dimensional models of volcanic plumes. A set of common volcanological input parameters and meteorological conditions were provided for two kinds of eruptions, representing a weak and a strong eruption column. From the different models, we compared the maximum plume height, neutral buoyancy level (where plume density equals that of the atmosphere), and level of maximum radial spreading of the umbrella cloud. We also compared the vertical profiles of eruption column properties, integrated across cross-sections of the plume (integral variables). Although the models use different numerical procedures and treatments of subgrid turbulence and particle dynamics, the inter-comparison shows qualitatively consistent results. In the weak plume case (mass eruption rate 1.5 × 106 kg s− 1), the vertical profiles of plume properties (e.g., vertical velocity, temperature) are similar among models, especially in the buoyant plume region. Variability among the simulated maximum heights is ~ 20%, whereas neutral buoyancy level and level of maximum radial spreading vary by ~ 10%. Time-averaging of the three-dimensional (3D) flow fields indicates an effective entrainment coefficient around 0.1 in the buoyant plume region, with much lower values in the jet region, which is consistent with findings of small-scale laboratory experiments. On the other hand, the strong plume case (mass eruption rate 1.5 × 109 kg s− 1) shows greater variability in the vertical plume profiles predicted by the different models. Our analysis suggests that the unstable flow dynamics in the strong plume enhances differences in the formulation and numerical solution of the models. This is especially evident in the overshooting top of the plume, which extends a significant portion (~ 1/8) of the maximum plume height. Nonetheless, overall variability in the spreading level and neutral buoyancy level is ~ 20%, whereas that of maximum height is ~ 10%. This inter-comparison study has highlighted the different capabilities of 3D volcanic plume models, and identified key features of weak and strong plumes, including the roles of jet stability, entrainment efficiency, and particle non-equilibrium, which deserve future investigation in field, laboratory, and numerical studies.YJS was partially supported by the ERI Cooperative Research Program and KAKENHI (25750142). The computations of SK-3D were carried out in part on the Earth Simulator at the JAMSTEC and also on the Primergy RX200S6 at the Research Computer System, Kyushu University. AC was partially supported by a grant of the International Research Promotion Office Earthquake Research Institute, the University of Tokyo. AC, TEO and MC were partially supported by the EU-funded project MEDiterranean Supersite Volcanoes (MEDSUV; grant no. 308665). MC acknowledges CINECA award N. HP10BKFD9F (2013) for high performance computing resources and support. AVE acknowledges NSF Postdoctoral Fellowship EAR1250029, a U.S. Geological Survey Mendenhall fellowship, and grant GID 61233 from NASA Ames Supercomputing Center