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

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 timeand 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 ms−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.


Introduction
During explosive volcanic eruptions, a mixture of gases, magma fragments, crystals and eroded rocks is injected into the atmosphere at high velocity, pressure and temperature.The diverse and unpredictable variability of eruptive styles depends mostly on the complex rheology of magma and the nonlinear processes leading to the fragmentation of the viscous melt into a polydisperse mixture of gases and particles (Gonnermann and Manga, 2007).Nonetheless, the explosive character of an eruption is always associated with the rapid decompression and the consequent abrupt expansion of gases in the magma (the exsolved magmatic volatiles in magmatic eruptions, vaporized free water or hydrothermal fluids in hydromagmatic and phreatomagmatic eruptions) (Parfitt and Wilson, 2008).Under such conditions, in the proximity of the volcanic vent, the erupted underexpanded multiphase mixture can manifest the features of supersonic flows (Kieffer, 1984;Woods and Bower, 1995;Esposti Ongaro et al., 2008;Ogden et al., 2008b;Orescanin et al., 2010), which, in turn, directly affect the plume source conditions.After decompression, column behavior is controlled by the balance between its negative buoyancy, associated with the load of solid particles, and the positive buoyancy due to air heating and expansion.Plume dynamics are therefore mainly influenced by (subsonic) turbulent mixing and mass and thermal exchange between the eruptive mixture and the atmosphere.Depending upon the efficiency of the turbulent entrainment, the gas-particle mixture can form a buoyant plume in the atmosphere or collapse under its particle load forming pyroclastic density currents (Valentine, 1998).
A general understanding of the transport dynamics of pyroclasts in the atmosphere was first achieved by describing the eruptive mixture as homogeneous, i.e., by assuming kinetic and thermal equilibrium between gas and particles and by solving the resulting transport equations under simplified conditions (e.g., one-dimensional and steady-state approximations) (Wilson, 1976;Woods, 1988;Sparks et al., 1997).Such an approach has also been extended to two and three dimensions and a transient regime (Oberhuber et al., 1998;Suzuki et al., 2005;Ogden et al., 2008a) in order to carry out numerical simulations of volcanic processes at the large scale, highlighting the key roles of environmental atmospheric conditions (Graf et al., 1999), large-eddy turbulence (Suzuki and Koyaguchi, 2010), vent overpressure (Ogden et al., 2008b) and boundary-layer processes (Doronzo et al., 2012).
However, the detailed reconstruction of well-documented eruptions and the growing need to quantify and map the hazards associated with future explosive events require the simulation of full eruptive scenarios.To this aim, eruption models able to incorporate the main dynamic processes and more realistic input conditions are needed.The problem, in its general multidimensional and unsteady formulation, is extremely challenging due to the multiphase nature of the flow as well as its multiscale features.
Mathematical models based on multiphase flow formulation have been proposed starting from the late 1980s (Valentine and Wohletz, 1989;Dobran et al., 1993), but have become more popular in the last decade (Dartevelle et al., 2004;Pelanti and LeVeque, 2006;Dufek and Bergantz, 2007;Esposti Ongaro et al., 2007) also thanks to the impressive development of computational techniques that allow for the solution of the complicated set of transport equations on modern high-performance parallel computers.However, further work is still necessary to achieve their systematic validation and to test their adequacy in simulating volcanic multiphase flows, especially as far as complex phenomena such as underexpanded jets are concerned.
This work is focused on the enhancement and validation of the PDAC model (Pyroclastic Dispersal Analysis Code), which is described in Sect. 2 and in more detail in Neri et al. (2003) andEsposti Ongaro et al. (2007).This model is able to solve the multiphase flow equations for a mixture of volcanic gases and pyroclasts in nonequilibrium conditions and exchanging momentum and heat.The PDAC numerical solution procedure is based on the original algorithm by Harlow and Amsden (1975), in which a first-order semi-implicit treatment for multiphase flows was combined with simple, donor-cell-based finite-volume conservative advection schemes.Typically, such first-order techniques introduce large amounts of numerical diffusion.Esposti Ongaro et al. (2007) extended the first-order spatial discretization to second order in each separate spatial direction by adopting the one-dimensional MUSCL scheme (Sweby, 1984), as is common practice in other multiphase flow codes (e.g., MFIX in Syamlal et al., 1993;Syamlal, 1998), but did not modify the semi-implicit time-advancement scheme based on the backward Euler discretization of time derivatives.The resulting numerical method was effective, but its results still display significant numerical diffusion, especially in multidimensional simulations, which implies the need for very high spatial resolution and small time steps to achieve an accurate simulation.
The first objective of this work is to modify the numerical algorithm in order to increase the accuracy in the simulation of the near-vent decompression dynamics, potentially involving supersonic regimes and shock waves, and the three-dimensional, transient dynamics of turbulent eddies that control, for example, the atmospheric air entrainment.We propose an improvement of the original PDAC numerical model to achieve higher accuracy and robustness in the resolution of compressible regimes, while reducing the numerical diffusion that may significantly damp turbulent eddies in subsonic regimes.More specifically, a second-order Crank-Nicolson-type time discretization (Crank and Nicolson, 1947) and a more accurate and fully multidimensional advection scheme (LeVeque, 1996) are introduced in the framework of the semi-implicit approach proposed by Harlow and Amsden (1975).
The resulting model is applied to the simulation of the decompression structures that form in the lower portion of volcanic columns, and is validated for the first time against laboratory experiments and by comparison with other model results.Although in the context of volcanic eruption simulations a rigorous model verification or validation is not possible (Oreskes et al., 1994), in recent years, three-dimensional multiphase flow models have demonstrated the potential for providing a good representation of the actual processes occurring in the real system (Dufek and Bergantz, 2007;Esposti Ongaro et al., 2012).The validation of numerical results against empirical observations of well-documented eruptions together with the congruence of numerical benchmarks with experimental and theoretical results are at present the only available instruments to assess the "empirical adequacy" (Oreskes et al., 1994) of models to simulate eruptive scenarios.As will be shown below, the results of the present model are in good qualitative and quantitative agreement with a number of experimental and numerical results available in the literature.
Finally, we investigate the capability of the enhanced multiphase model to catch nonequilibrium dynamics between gas-phase and solid particles in transonic regimes.Indeed, the assessment of the interactions between shock waves and solid particles with different sizes is still an almost unexplored subject in the multiphase flow literature, and its deep understanding could also have further impact beyond volcanological applications.
In Sect.2, we describe briefly the PDAC model and the simplifying assumption adopted in the present work.In Sects. 3 and 4, the present, improved version of the PDAC numerical method is described.Numerical benchmarks and comparison with analogous results obtained with other models are presented in Sect. 5. Finally, some conclusions and perspectives for future work are presented in Sect.6.

Multiphase flow equations
In this work we employ the same model equations as in the original PDAC model proposed by Neri et al. (2003).They are appropriate to describe the injection and dispersal of a hot and high-velocity gas-pyroclast mixture in a standard reference atmosphere.The model is based on the following main hypotheses: -the solid particles and the gas are considered as interpenetrating continua, following an Eulerian-Eulerian approach; -the gas phase is compressible and obeys the ideal gas law; -mass transfer processes due to phase changes and chemical reactions are neglected; -solid particles are assumed to be spherical, and each class is assumed to consist of particles of equal radius and density; -particles are assumed to maintain their original size, thus neglecting the effect of any secondary fragmentation or aggregation process on the large-scale dispersal dynamics; and -the heat transfer between different solid phases, as well as the viscous dissipation effects, are neglected due to their second-order effect in comparison with advection, conduction and gas-particle heat exchange.
While in general turbulence and other dissipative effects cannot be neglected, in this study, for the sake of simplicity, we only focus on inviscid equations and regimes.As reported in detail in Carcano et al. (2012), a scaling analysis of the model equations shows that the typical values of the Reynolds number Re = ρU L µ and the Péclet number P e = c p ρU L k (based on diameter and velocity at the vent and on average mixture properties) in a volcanic jet vary from 10 6 to 10 11 in the regimes of interest.Consequently, the only dissipative terms retained in the following are those representing interphase exchange processes between the gas and the solid phase.However, all the physical processes neglected in this study are actually accounted for in the complete PDAC model in the same way as in the original model proposed in Neri et al. (2003), to which we refer for a more complete description of these terms.Numerical results, not shown here, confirm the validity of our approximation in the explosive regimes.
The gas phase is composed of different chemical components leaving the crater, such as water vapor and carbon dioxide, and atmospheric air, considered as a single chemical component.The pyroclasts are described by N classes of solid particles, each one characterized by a diameter, density, specific heat and thermal conductivity.In the following sections, we will denote with the subscript s = 1 . . .N the classes of solid particles and with l = 1 . . .M the chemical components of the gas phase g.The model variables can be defined as follows: g , s = volumetric fractions of gas and solid particles; if V is the representative volume and V g and V s are the volumes occupied by gas and particles, respectively, then the gas and solid volume fractions are defined by ρ g , ρ s = microscopic densities of gas and solid particles; y l = mass fractions of the gas components; v g , v s = velocities of gas and solid particles; p g = gas pressure; h g , h s = enthalpies of gas and solid particles; and -T g , T s = temperatures of gas and solid particles.
The model consists of 5(N + 1) + M coupled partial differential equations for the independent variables p g , ρ g , s , v s , h s (or T s ), y l , with s = g, 1 . . .N and l = 1 . . .M. The mass conservation equations for the gas phase g, the s-th solid phase and the l-th gas chemical component are

S. Carcano et al.: A numerical model for volcanic jets
The momentum balance equations for the gas phase and the s-th solid phase, for all s = 1 . . .N are written as (with p = s).Here, D p,s represents the drag coefficient describing the interaction between the phase p and the phase s, and g denotes the gravitational acceleration vector.In the present formulation, we adopted model A by Gidaspow (1994), in which the gas pressure gradient is present in both the gas and particle momentum equations.Numerical experiments assuming a pressureless particulate phase (model B) do not show significant differences, at least in the dilute regime under investigation.The energy balance equations for the gas phase and the solid phases s = 1 . . .N are written in terms of their enthalpies: Here, Q s is the volumetric heat transfer rate between the gas and the s-th solid phase.For the gas phase, we have considered the reversible rate of enthalpy change due to compression or expansion, which is important in transient, compressible flows.Heat transfer between different solid phases is negligible, and also radiative heat transfer has not been considered.As remarked before, viscous dissipation has been neglected for the applications considered in this paper based on the results of the scale analysis.By definition of the volumetric and mass fractions, one also has the relations The gas phase is compressible, and we suppose that thermodynamic quantities are related by the ideal gas law: where R is the gas constant of the mixture of gaseous components.Particulate solid phases are considered incompressible.Consequently, their microscopic density is assumed to be constant and denoted by ρ s , s = 1 . . .N .The temperature of each phase is derived from its enthalpy as where particle specific heats c p,s are assumed to be constant and to correspond to average values due to their minor sensitivity on temperature.The specific heat of the gas phase c p,g depends on temperature, and it is computed as a weighted average of the specific heats of the M chemical components: Interphase drag coefficients and heat transfer rates are derived from semi-empirical correlations for dilute and dense regimes.
The initial values of all field variables must be specified for the entire computational domain.Usually, a standard atmosphere, vertically stratified in pressure, temperature and density, is considered throughout the domain.The atmosphere is composed of dry air at rest, and no particle of any size is considered present in the computational domain.Appropriate boundary conditions will be described later for each specific test case.

The numerical method: semi-implicit time discretization
The model equations described in the previous section are discretized in time by a second-order extension of the original first-order implicit multifield (IMF) algorithm proposed in Harlow and Amsden (1975).We will describe the time discretization method in the simpler case of a single solid phase s.We employ a semi-implicit time discretization based on a Crank-Nicolson-type time averaging (also known as the θ method) with averaging parameter θ ∈ [0, 1] (Crank and Nicolson, 1947).It is well known (see, e.g., Quarteroni et al., 2002) that, for unconditional linear stability, one has to choose θ ≥ 1/2, while full second-order accuracy is only granted for the limit value θ = 1/2.In general, the value θ = 0.55 is employed in most numerical simulations in order to guarantee stability also in the nonlinear case.
The continuity equation for the phase s is discretized as The momentum equation for the phase s, with p = s is discretized as Notice that the pressure, the gravity and the drag terms are discretized in time by the θ method, while flux terms are treated explicitly.The enthalpy equations for the gas phase g and the solid phase s are solved after continuity and momentum equations.Only the interphase exchange terms are treated semi-implicitly by the θ method, while convective terms are treated explicitly by using the updated densities and velocities.
where the E terms include all the explicit terms.
For each time step t n+1 , Eq. ( 16) is solved by this secondorder extension of the IMF algorithm as follows: -temperature-dependent coefficients of the gas phase are computed; -the interphase coefficients D g,s and Q s and the explicit E terms are computed; -the coupled continuity and momentum equations are solved iteratively by the approximate Newton method to update velocity fields, pressure and volumetric fractions; -gas mass fractions y l are computed by solving the linear transport equations; and -the energy equations are linear in the temperatures and decoupled from the continuity and momentum equations and can finally be solved directly.
We observe that, since the enthalpy equations are solved explicitly after the solution of the momentum and continuity equations, the temperature is kept constant during the solution procedure.The effect of the temperature variation on the gas pressure and density are deferred to the next time-step computation.

The numerical method: space discretization
The model equations are discretized by a finite-volume approach on an orthogonal, nonuniform mesh.A staggered discretization grid with N x × N y × N z Cartesian cells is introduced, along the lines of popular discretization methods such as the MAC (marker and cell) approach, introduced in Harlow and Welch (1965), or the Arakawa C grid (see, e.g., Arakawa and Lamb, 1981).The mesh is composed by rectangular control volumes obtained from the cartesian product of three 1-D discretization intervals along each axis.Each cell is numbered at its center with indices i, j and k, for the x, y and z directions, respectively.If we denote by i = 1 . . .N x , j = 1 . . .N y and k = 1 . . .N z the intervals along the x, y and z axis, respectively, the center of each cell can be identified by a triplet of indices (i, j, k).The length of the cell sides in each directions are denoted by x i , y j and z k , and they are assumed to vary in their respective direction only.The cell volume is given by V i,j,k = x i y j z k , and staggered spacings x i+ 1 2 are defined as the arithmetic average of the neighboring integer index values.
The discrete u velocity is defined at half-integer i and integers j and k; v is defined at integers i, k and half-integer j ; and w is defined at integers i, j and half integer k.Finally, p and all other three-dimensional scalar variables (i.e., pressure, densities, volumetric fractions and enthalpies) are defined at integers i, j and k.Therefore mass and enthalpy equations are solved on the cell centers, whereas the momentum equations are solved at the staggered locations.At points where they are not defined, the discrete variables are generally computed by linear interpolation of the nearest values.Averaged quantities will usually be denoted by an overbar.On a uniform grid, for example, ūi,j,k = If we denote with brackets . . . the discretization of the advective fluxes and we adopt the staggered approach described above, we obtain, for both the gas phase and the solid phase, for each cell (i, j, k) of the mesh, the following system of discretized equations: for all s, p = g, 1 . . .N and p = s, where the discretization of the E terms is specified in Appendix B. The E terms are computed explicitly before the resolution of the system.Finally, we introduce an appropriate discretization technique for the advective fluxes.From the nondimensional analysis (Carcano et al., 2012), we know that advection is one of the dominating phenomena in the process, so we expect that a proper numerical treatment of the advection terms should be necessary in order to obtain an accurate numerical solution.Therefore, one of the first possible modifications of the donor-cell scheme is to introduce in the upwind discretization the so-called transverse fluxes (Colella, 1990;Saltzman, 1993;LeVeque, 1996).In the standard donor-cell upwind method, the advective flux through one single cell boundary is split into independent fluxes along the x, y and z directions by using the velocities u, v and w in the directions normal to each interface.More accurate methods can be obtained by considering the flux with the proper speed v = (u, v, w) without splitting it along the three space directions.Unsplit upwind schemes are more accurate (in partic-ular, less diffusive) and more stable than the original donorcell upwind method.
This improved version of the upwind method is called the corner transport upwind (CTU) method (Colella, 1990).A hierarchy of methods for the numerical solution of advective transport in conservation equations in several space dimensions based on CTU was proposed by LeVeque (1996).In the present model, we employ one of the second-order versions of the algorithm described in LeVeque (1996), including minmod flux limiting (Roe, 1986) to avoid the creation of spurious extrema in the solution.
As an example, in the simpler two-dimensional case, the advective flux Qu ik of the scalar quantity Q along the x direction in the computational cell (i, k) is computed as For where lim represents the flux limiter.Analogous expressions are written for u ik− 1 2 < 0 and/or w ik− 1 2 < 0. The first term on the right-hand side represents the donor-cell upwind flux, the second term represents the CTU correction and the last term represents its second-order extension.

Model validation: axisymmetric underexpanded jet
The proposed numerical method was tested on two-and three-dimensional simulations of underexpanded jets, and the numerical results obtained were compared with both experimental and numerical results available in the literature.In order to simulate a cylindrical underexpanded jet, we assume that each phase enters the domain through a fixed inlet where the volume fraction, velocity and temperature of each phase and the gas pressure are imposed.Mechanical and thermal equilibrium between the phases at the vent are assumed.In two-dimensional tests, we solve the model equations in cylindrical coordinates, and we impose symmetry conditions at the left lateral boundary.At the bottom boundary, no mass and heat transfer are allowed, and free-slip conditions are assumed for the velocity of each phase.At the upper boundary, free-outflow/inflow conditions are assumed, whereas at the lateral boundaries it is possible to assume either free-slip or free-outflow/inflow conditions.In particular, at the outflow boundaries, the mass and momentum equations of the mixture are solved for pressure, assuming a null velocity gradient along the boundary.At the lateral inflow boundaries, incoming air is assumed to be free of particles and to have pressure and temperature characteristics corresponding to those of the standard reference atmosphere.
Here and in the following sections, the chosen time step is such that the maximum Courant number based on the gas velocity (CFL max u = max(|v g | t/ x)) reached in the domain is equal to 0.2 in order to guarantee the stability of the explicit part of the numerical scheme.All the computations were carried out with a parallel version of the improved PDAC model based on parallel implementation described and tested extensively in Esposti Ongaro et al. (2007).The most intensive computations for the three-dimensional tests reported in section 5.3 were carried out in parallel on 32 processors on a Linux cluster with 128 CPUs at 2.4 GHz, with total peak performance of 580 Gflops s −1 , and required about 44 h to reach the final time (T = 2 s).

Comparison with laboratory results and empirical laws
We present here a set of numerical tests aimed at the simulation of pure gas and gas-particle jets at the laboratory scale.It has been proven theoretically and experimentally that vents with supersonic or sonic vertical velocity and gas pressure greater than the atmospheric one result in a rapid expansion and acceleration of the fluid to a high Mach number (Lewis and Carlson, 1964).A series of expansion waves form at the vent exit (Prandtl-Meyer expansion), which are reflected as compression waves at the jet flow boundary.The compression waves coalesce to form a barrel shock and a standing normal shock wave (Mach disk), across which the vertical velocity is reduced and the pressure in the core of the jet increases.The fluid that crosses the Mack disk is rapidly compressed and decelerated to subsonic speeds.Above the Mach disk, the fluid moves slowly in the core of the jet and is surrounded by a supersonic moving shell, with a slip line or a shear layer dividing these regions, as shown in Fig. 1.One of the important parameters describing these supersonic jets is the distance between the vent and the normal shock wave.Experimental results reported in Lewis and Carlson (1964) show that the height of the Mach disk h d in a pure gas jet depends on the vent diameter D v , the exit Mach number Ma v and the ratio of specific heats γ as well as the ratio K of the exit static pressure P v and the atmospheric pressure P atm .The empirical relationship is In the case of multiphase gas-particle underexpanded jets, the location of the Mach disk depends also on the particle loading η = s ρ s g ρ g at the inlet.Even if different empirical relationships between Mach disk height and particle loading are proposed in the literature (e.g., Lewis and Carlson, 1964;Jarvinen and Draper, 1967;Sommerfeld, 1994), all of them predict an upstream movement of the Mach disk and a reduction of the Mach disk distance from the inlet.
In this section, the aim is to verify whether the proposed multiphase model is able to reproduce correctly the wave pattern that forms above an overpressured vent by evaluating the Mach disk location -first in the case of a homogeneous gas, and then for a gas-particle mixture.

Homogeneous jet
We consider a homogeneous fluid (dry air with standard chemical components), and we impose underexpanded sonic or supersonic conditions at the inlet (Table 1); that is, the gas pressure at the inlet is larger than the atmospheric and the Mach number Ma = |v|/c s ≥ 1, where c s is the speed of sound.
The computational domain is a box of size 0.1 × 0.2 m 2 , whose left side coincides with the axis of the vent.The side and the bottom boundaries of the axisymmetric domain are impermeable and stress-free.Two uniform meshes of 160 × 320 ( x = z = 6.25 × 10 −4 m) and 500 × 1000 ( x = z = 2 × 10 −4 m) cells have been employed, with time steps of t = 10 −7 s and t = 5 × 10 −8 s, respectively.
We consider different values of overpressure levels K, and we evaluate the height of the Mach disk h d .We obtain a good agreement between experimental results and numerical simulations, as shown in Fig. 2. In Figs. 3 and 4, the results obtained in terms of vertical velocity and temperature are shown.The improved version of the PDAC code has a better fit with experimental results and is able to describe the shear layer instability above the Mach disk.In general, the  first-order version of the model tends to underestimate the Mach disk height.Moreover, for small values of overpressure K, using first-order methods we do not see the formation of the Mach disk.Second-order methods are able to capture the sharp discontinuity in the flow, as shown in Fig. 5, and allows for a better estimate of the empirical law in Eq. ( 25) to be obtained.The numerical simulation with the improved numerical scheme have been repeated with different meshes in order to check the dependence of the results on the grid resolution.Results in Figs. 6 and 7 show the axial profiles of gas pressure obtained with x = 1, 0.625, 0.5, 0.25, 0.2 mm and the estimates of the Mach disk height, respectively.We observe that the Mach disk height estimate improves with grid reso- lution and that the second-order method is able to provide a reliable estimate of the Mach disk position even with coarse meshes.

Nonhomogeneous gas-particle jet
When solid particles are added to the gas flow, new phenomena associated with kinetic and thermal nonequilibrium between the gas and particulate phases arise.Such effects are controlled by drag and energy exchange terms in the momentum and energy equations, which are recalled in Appendix A.
To quantify the importance of nonequilibrium regimes, dilute gas-particle flow can be characterized by a timescale (the particle relaxation time) determined by the balance between particle inertia and gas-particle viscous drag (Marble (1970) and Burgisser (2002)).Its expression for monodisperse mixtures can be derived from the momentum balance equation for the solid particles (Eq.5) by neglecting all the terms except the drag and inertial terms: where we considered only the component along z, since it is much greater than the horizontal components in the decompression region.The relaxation time is thus defined from Eq. (26) as A simple analysis, e.g., that in Marble (1970), suggests that the timescale for thermal relaxation has the same order of magnitude.For dilute mixtures ( g ≈ 1) and low gas-particle Reynolds number, the particle relaxation time approximates that of a single particle in a laminar flow (Stokes' regime): However, in general, the drag coefficient in Eq. ( 27) is a complex function of the particle concentration and the gasparticle Reynolds number, which is defined as where d s is the particle diameter and µ g is the dynamic viscosity of the gas phase.In the underexpanded jet under investigation, the flow is always in a dilute regime, with g > 0.8.In these conditions, Eq. (A1) for the drag coefficient (Wen and Yu, 1966) can be adopted, and the relaxation time becomes where the coefficient C d,s depends on the gas-particle Reynolds number Re s , as reported in Eq. (A2).
The relaxation time τ s gives an order of magnitude of the time delay with which a particle equilibrates to a timevarying gas flow.In a supersonic jet, such delay may occur in the rapid expansion region above the vent, where a difference between gas and particle velocity w may be expected.We should then compare the particle relaxation time with the formation time of the Mach disk (Orescanin et al., 2010), estimated as where c s,mix is the mixture speed of sound, as defined, e.g., by Pelanti and LeVeque (2006): In this section, we consider a mixture of dry air and fine solid particles with diameter equal to 10 µm and density equal to 2500 kg m −3 .Inlet flow parameters correspond to experimental and simulation conditions investigated by Sommerfeld (1994).Gas and particle velocities are both equal to the speed of sound in the pure gas, whereas the overpressure of the gas phase is K = 31, producing supersonic underexpanded conditions at the inlet (Table 2).According to the simple scaling analysis discussed above, the ratio between the particle relaxation time and jet timescale is about τ s /τ Ma > 10 2 .Therefore, we expect that particles will be loosely coupled to the gas phase and that they will not have the time to equilibrate to the expanding gas flow near the vent.
The computational domain is a box of size 0.15 × 0.225 m 2 and, as in the previous test cases, the left side coincides with the axis of the vent, whereas the side and the bottom boundaries of the axisymmetric domain are impermeable and stress-free.A nonuniform mesh of 500 × 750 computational cells have been employed, with time step of t = 2×10 −8 s.The maximum resolution is imposed above the inlet, where x = z = 10 −4 m.We consider different values of particle volume fractions s at the inlet, and we evaluate the height of the Mach disk h d .
In Fig. 8 we report the results of four different simulations of particle-laden underexpanded jets with different particle concentrations.The gas phase expands radially as in the homogeneous case, thus increasing the final jet radius up to three times in correspondence to the Mach disk location.
On the other hand, as expected from the scaling analysis and also observed in laboratory (Sommerfeld, 1994), particles are almost unaffected by the rapid gas expansion.Particle trajectories remain nearly vertical, with some radial spreading that is almost independent of particle concentration.In the expansion region, particles are only slightly and gradually accelerated by the gas phase (at most by 20 m s −1 , up to 368 m s −1 ) where the Mach disk is located, and then slowly decelerated in the subsonic region, but they never reach an equilibrium condition with the gas phase.However, solid particles tend to deform the Mach disk, moving it towards the vent and causing it to become more concave as the particle loading increases.For initial particle volume fraction equal to 0.0005, the normal shock is located 14.5 mm from the inlet, and the distance is reduced to 12.2 mm when s = 0.004.Moreover, when the particle loading is increased, the expansion, the acceleration and the cooling of the gas phase are reduced, as shown in Fig. 9.As regards particle distribution, the mixture density profile along the jet axis is not affected by the presence of the shock wave but remains almost constant and displays only a small reduction due to the radial spreading, which tends to increase with increasing particle loading, as observed also by Sommerfeld (1994).

Pseudogas regime
When particle relaxation time is much lower than the jet timescale, particles are tightly coupled to the gas phase.
Under such assumption, Ogden et al. (2008b) assumed perfect kinematic and thermal equilibrium between the phases and described the eruptive mixture as a homogeneous pseudogas characterized by average thermodynamic and rheologic properties.Two-dimensional numerical simulations of  (kg m −3 ) 2500 underexpanded volcanic jets were performed with CFDLib, a computational fluid dynamics library developed at Los Alamos National Laboratory that uses a finite-volume computational scheme with cell-centered state variables.CFDLib applies a variation of the implicit continuous-fluid Eulerian (ICE) method, proposed in Harlow and Amsden (1968) and Harlow and Amsden (1975), and a modified Godunov method (Godunov, 1999) to solve shock waves.Following Ogden et al. (2008b), numerical simulations presented in this section are performed in absence of gravity in order to focus on compressibility and nonequilibrium multiphase effects.We assume choked flow conditions at the vent; that is, the inflow velocity is equal to the speed of sound of the mixture, defined by Eq. ( 32).
Inlet conditions for numerical simulations are specified in Table 3.The mixture is composed by water vapor and solid particles that are injected in a standard atmosphere composed by dry air.We first consider a single solid dispersed phase with particle diameter equal to 10 µm.Two different inlet pressure ratios of K = 20 (case A) and K = 5 (case B) were adopted.A third run (case C) is performed with K = 5 and a mixture of two particle phases of 10 and 1000 µm equally distributed in weight.
In case A, the computational domain is a box of size 800 × 2400 m 2 , and we use a uniform 200 × 600 mesh, with x = z = 4 m and a time step t = 10 −3 s. Figure 10 shows the vertical velocity field of the gas phase and the particle distribution above the vent after 20 s, when quasisteady-state conditions are reached.The simulation reproduces the expected behavior of a supersonic underexpanded jet, displaying the barrel shock with a convex Mach disk, at about 320 m above the vent, which decelerates the mixture down to subsonic velocities (see Fig. 10a) and compresses the gas phase so that the particle volumetric fraction increases by one order of magnitude across the discontinuity, as shown in Fig. 10b.
To better analyze the jet dynamics, and to quantitatively compare our results with those of Ogden et al. (2008b), we study the time-averaged vertical profiles along the axis of pressure, mixture density, gas vertical velocity and gas temperature, shown in Fig. 11.The gas phase undergoes a rapid expansion from the initial pressure of 2.02 × 10 6 Pa to pressure values below atmospheric pressure.The minimum of the pressure is 9.1 × 10 3 Pa, and it is reached at the height of 324 m above the vent.The ratio between Mach disk height and vent radius is 8.1, and the difference with respect to the corresponding result of Ogden et al. (2008b) is around 1 %.Through the normal compression shock, the gas phase returns to atmospheric value.During the expansion, as expected in supersonic flows, the gas phase accelerates up to 482 m s −1 , and then through the shock it abruptly decelerates to a subsonic regime, with a vertical velocity around 33 m s −1 .During the expansion and acceleration phase, the gas decreases its temperature down to 1104 K and then warms up again by about 70 K through the Mach disk.Mixture density, defined as ρ m = g ρ g + s ρ s , decreases by two orders of magnitude above the vent, and then it increases by one order of magnitude through the shock.The difference in the Mach disk position with respect to the results reported in Ogden et al. (2008b) is around 4 %, whereas the difference in the maximum vertical velocity is around 2 %.
In case B we consider an inlet pressure ratio of K = 5 and a vent diameter of 20 m (Table 3) in order to maintain the sonic conditions at the vent.The computational domain is a box of size 200 × 400 m 2 , and we use a uniform 200 × 1000 mesh, with x = 1 m, z = 0.4 m and a time step t = 5×10 −5 s. Figure 12 shows the gas vertical velocity and the particle volume fraction when the quasi-steadystate configuration of the normal shock is achieved.The twodimensional jet pattern and shape closely fit the results presented in Fig. 3b of Ogden et al. (2008b).The results obtained with the multiphase model are thus in quantitative agreement with the result obtained by Ogden et al. (2008b), demonstrating that the pseudogas approximation is appropriate for the description of the underexpanded jet regime when particle sizes are lower than 10 µm.
Finally, in case C we consider the same configuration as in case B and we change the gas and particle mixture by introducing a second class of solid particles with diameter equal to 1000 µm.The vertical, axial profiles of case B and C, displayed in Fig. 13, report a Mach disk height around 39 m from the vent and the ratio between the Mach disk height and (kg m −3 ) 1000 1000 1000 the vent radius equal to 3.9, in agreement with the numerical and experimental results reported in Ogden et al. (2008b).In case B, the maximum vertical velocity of the gas phase is 413 m s −1 , while the minimum temperature along the axis is 1132 K, showing a temperature decrease of about 5 % with respect to the vent temperature.
In case C, the Mach disk position is unchanged, whereas we observe a peak velocity about 30 m s −1 lower.Above the normal shock, gas velocity is 56 m s −1 in case B and 87 m s −1 in case C. Flow density is also considerably higher in case C. Such differences between case B and C are associated with the presence of a coarser particulate phase, whose effects are here analyzed in more detail.

Assessment of multiphase effects
The assessment of the influence of fine and coarse particles on the structure of shock waves is still an open issue in the multiphase flow literature, and its deep understanding could also have further impact beyond volcanological applications.To assess the influence of nonequilibrium effects on the jet dynamics, we adopt the scaling analysis presented above.To estimate the magnitude of the relaxation time, we first estimate from numerical results the maximum relative Reynolds number, as defined in Eq. ( 29) in order to estimate the drag coefficient D g,s .We consider that the maximum disequilibrium is achieved across the Mach disk, where the gas phase is decelerated almost instantaneously, while particles cross the discontinuity undisturbed before slowing down by the effect of viscous drag.By using the gas velocity jump across the shock as a proxy of the velocity difference between the gas and the particles and the gas density after the shock, we estimate the maximum relative gas-particle Reynolds number on the order of 24 (case A) and 19 (case B).Therefore, we can use Eq. ( 28) to estimate the particle relaxation time, obtaining τ s 1.5×10 −4 s (by assuming water vapor viscosity at 1000 K equal to µ g = 3.7 × 10 −5 kg m −1 s −1 ).We can compare the particle relaxation time with the formation time of the Mach disk defined by Eq. (31).Taking T = 1000 K and considering a dilute mixture of water vapor and solid particle with ρ s = 1000 kg m −3 , s = 0.01 and ρ g = 0.2 kg m −3 , we obtain c s,mix 300 m s −1 ; in case A τ Ma 0.1 s, and in case B τ Ma 0.03 s.Therefore in both case A and case B, the particle relaxation time is much smaller than the formation time of the Mach disk (τ s τ Ma ), thus meaning that fine particles dynamics are strongly coupled with the gas dynamics.Multiphase effects are negligible and the pseudogas approximation is appropriate.
In case C, we can proceed in an analogous way by assuming that coarse particles move in a fluid composed by water vapor plus fine particles in mechanical and thermal equilibrium, described as a pseudogas (as verified for cases A and B), and by adopting the same Eqs.(A1) and (A2).We therefore compute the average properties of the pseudogas ρ ps , µ ps , and use them to estimate the relaxation time for the coarser particles in a bidisperse mixture.
The particles volume fraction of fine particles decreases by one order of magnitude in the first 10 m above the vent, and further down to 5 × 10 −4 before the shock.The pseudogas density before the shock is ρ ps = g ρ g + s 1 ρ s 1 = 1.2 kg m −3 .Its mean viscosity can be computed as µ ps µ g (1 + 2.5 s 1 ) = 3.71 × 10 −5 Pa s (Einstein, 1906), not very different from the value for pure water vapor.The maximum gas-particles Reynolds number, computed from the gas velocity jump across the shock |w b g − w a g | = 298 m s −1 , is Re s 9 × 10 3 , thus implying that the low-Re approximation in Eq. ( 28) for the relaxation time is not applicable.In regimes where Re s > 1000, we can estimate the relaxation time with the Reynolds number correction, as reported in Eq. (A2); that is,    In case C the particle relaxation time and the formation time of the Mach disk are comparable, and thus it is worth investigating in detail the nonequilibrium effects for the coarsest particles.
We first estimate the magnitude of the velocity difference between gas and particle below the shock as w est ≈ ατ s w b g , where α = dw g dz is the gas velocity vertical gradient below the shock and w b g is the gas vertical velocity below the shock, as computed in the equilibrium gas-particle flow.Based on the results discussed in the previous section (see Figs. 11 and 13), α is taken as constant.In the expansion region, the velocity gradient α is approximatively equal to 6 s −1 , the gas velocity is w b g = 386 m s −1 , and we obtain w est C 23 m s −1 .Across the normal shock, on the other hand, particles will equilibrate to the gas flow within a distance l, also known as shock relaxation (Marble, 1970).We can derive an estimate of the relaxation distance above the shock as l est = |w b g −w a g |•τ s , where |w b g −w a g | is the gas velocity jump across the normal shock.Across the Mach disk, the jump in vertical gas velocity is 299 m s −1 .The distance to which particles equilibrate to the gas flow above the shock can be computed by assuming an initial disequilibrium velocity of the same order of magnitude and the relaxation time, obtaining 14 shows the differences between gas and particles velocity and temperature in case C as they result from the numerical code.We observe that smaller particles are essentially in thermal and mechanical equilibrium with the gas phase, as expected from theoretical results.Larger particles just below the normal shock are slower than the gas by about 25 m s −1 , and they cross the shock with a vertical velocity that is about 140 m s −1 larger than the gas vertical velocity before they finally reach an equilibrium velocity close to the gas velocity.Larger particles tends to cool slower than the gas and the smaller particles, and as such their temperature is about 70 K higher than the gas temperature when they reach the Mach disk.The estimated relaxation length is larger than the vertical grid size (0.4 m), and as a result numerical resolution appears adequate to quantitatively resolve shock relaxation.In particular, after the normal shock, particles slow down to the gas velocity within 6 computational cells, from 40.4 m to 42.8 m.The numerical estimate of the relaxation distance l C 2.4 m is thus comparable with the theoretical one l est C 3 m.Numerical results are therefore consistent with estimates derived from the pseudogas solution and a simple dimensional analysis based on the particle relaxation time.

Three-dimensional simulations
We repeated some of the simulations presented in the previous sections in a three-dimensional configuration in order to compare the results with those obtained in the twodimensional axisymmetric tests.In this section, we present the results obtained with the vent conditions of case B (see Table 3).The computational domain is a box of size 400 × 400 × 400 m 3 , and we use a nonuniform mesh of 120 × 120 × 120 cells and a time step of t = 10 −4 s, with maximum grid resolution of 1 m (equal to that employed in 2-D) in a subdomain of 50 × 50 × 50 m 3 above the vent.
Figure 15 shows the isosurfaces of the gas vertical velocity at t = 2 s. Figure 16 shows the 3-D vertical velocity and the logarithm to the base 10 of total particle volumetric fraction averaged along the aximuthal angle.With respect to the 2-D simulation in cylindrical symmetry (Fig. 12), 3-D simulation displays a more diffused jet boundary likely associated with the effect of the noncircular inlet.Indeed, with Cartesian mesh discretization, the circular vent is approximated by squared cells.The flow density have been opportunely corrected in cells cut by the inlet rim by proportionally reducing the particle concentration in order to impose the correct mass flow rate.However, boundary conditions do not describe the curved inlet rim.This produces some axial switching of the jet cross section in the subsonic region above the Mach disk, analogous to that observed in noncircular subsonic jets (e.g., Gutmark and Grinstein, 1999).
Nonetheless, the shock wave pattern (location and shape of the Mach disk and slip lines) is analogous to the 2-D case.Figure 17 shows that the time-averaged axial profiles of pressure, velocity and mixture density are consistent with the results obtained with two-dimensional axisymmetric simulations.Temperature difference probably reflects some differences in the average axial distribution of solid particles associated with the approximate vent geometry.

Conclusions
The PDAC multiphase flow model has been improved in several aspects of the numerical algorithm to modified the temporal and spatial accuracy of the simulation of explosive volcanic eruptions.In the new model version, a secondorder Crank-Nicolson-type time discretization has been introduced and the fully multidimensional advection schemes proposed by LeVeque (1996) have been employed.The new model has been tested against the complex problem of volcanic jet decompression in both two and three dimensions.Since a proper validation with volcanic jet data is not yet possible, due to the large scale, the dangerous nature of the phenomenon and the difficulty of remote measurements, we have verified that the numerical results adequately reproduce some similar phenomenology (i.e., an underexpanded, supersonic gas jet) as measured at the laboratory scale, where the new numerical scheme demonstrates better performance (in terms of accuracy and reduced numerical diffusion) with respect to previous model versions at all regimes.
However, several aspects of the dynamics of volcanic jets make them different from their laboratory analogues: volcanic jets involve the explosive decompression of a multiphase gas-particle mixture at high temperature, with a wide spectrum of particle grain sizes.To account for the presence of solid particles in supersonic volcanic jets, previous works have described the eruptive mixture as a homogeneous pseudogas (e.g., Kieffer, 1984;Woods and Bower, 1995;Ogden et al., 2008b;Orescanin et al., 2010).In the limiting case of fine particles (having a diameter on the order of 10 microns and relaxation time on the order of 10 −4 s, much smaller than the characteristic time for decompression), we have shown that the multiphase PDAC model consistently reproduces predictions of the pseudogas model.However, in the case of coarse particles and polydisperse mixtures, multiphase ef-fects become more important and also affect the average jet dynamics.
The new numerical code appears suited for the multiphase flow simulation of explosive regimes characterized by rapid decompression of the eruptive mixture and possible transition to the supersonic regime, including the development of impulsive Vulcanian eruptions and volcanic blasts.This application will be addressed in future works, which will be aimed at quantifying multiphase effects also for polydisperse mixtures, where particle-particle drag might play a key role in the nonequilibrium dynamics of gas-particle flows.
The present rigorous verification of a multiphase flow model for volcanological applications demonstrates the need for (at least) second-order numerical discretization schemes.The use of numerical approximations of an order higher than two may not be a major concern, if only volcanological applications are considered, due to the large uncertainties in the available data and especially in the initial conditions.However, in order to improve the description of different portions of the volcanic column, where turbulence effects become more important, and in view of different applications of the multiphase model to industrial problems on the laboratory scale, we are planning to further improve the accuracy by adopting spatial discretizations based on discontinuous Galerkin methods, which allow for increasing the order of the approximation while retaining a compact stencil.Concerning computational efficiency of the time-advancing scheme, we are also considering explicit time discretizations, such as Runge-Kutta methods, that have been widely applied in the literature to solve gas-dynamics problems at high Mach numbers.
Concerning the model formulation, the most critical aspect pertains to particle thermodynamics, since particle-particle collisions are nonnegligible for volume concentrations above about 10 −3 (Gidaspow, 1994).To improve model reliability in such regimes, we are moving from the present formulation (employing a semi-empirical description of solid pressure and equation of state) to a more rigorous closure based on the kinetic theory of dispersed granular materials.However, in the present application, such improvement is not critical since the dynamics are strongly driven by gas pressure terms and gas-particle drag, whose form is well established.The expressions of the discrete explicit terms of the energy equations of the gas and the solid phase are the following:

Fig. 2 .
Fig. 2. Comparison between experimental and numerical results in terms of Mach disk height for different values of the vent overpressure K.The results inLewis and Carlson (1964) are compared with numerical simulation applying the first-order donor-cell upwind method (FOU), upwind method with second-order MUSCL fluxes (MUSCL) and second-order corner transport upwind method with θ method time discretization (CNCTU).

Fig. 5 .
Fig.5.Gas pressure at t = 2 × 10 −3 s.Axial profile for different values of the vent overpressure K computed on a coarse mesh.Comparison between first-order (FOU) and second-order methods (MUSCL and CNCTU) on a 160 × 320 mesh.

Fig. 6 .
Fig.6.Axial profiles of gas pressure at t = 2×10 −3 s obtained with second-order methods and different grid resolutions.

Fig. 7 .
Fig. 7. Estimates of the Mach disk height obtained with secondorder methods with different grid resolution.Comparison with the theoretical estimate.

Fig. 9 .
Fig.9.Axial profiles of gas pressure, mixture density, gas vertical velocity and gas temperature at t = 3 × 10 −4 s.Comparison between the homogeneous jet's profile and results obtained for different values of initial particle volume fraction s = 0.0005, 0.001, 0.002, 0.004.

Fig. 10 .
Fig. 10.Case A. (a) Gas vertical velocity [m s −1 ] and (b) logarithm to the base 10 of particle volume fraction.Snapshots at t = 20 s.

Fig. 11 .
Fig. 11.Case A. Time-averaged axial profiles computed over the interval [16, 20] s of gas pressure, mixture density, gas vertical velocity and gas temperature.

Fig. 12 .
Fig. 12. Case B. (a) Gas vertical velocity [m s −1 ] and (b) logarithm to the base 10 of particle volume fraction.Snapshots at t = 2 s.

Fig. 13 .
Fig. 13.Case B and C. Time-averaged axial profiles computed over the interval [1.6, 2.0] s of gas pressure, mixture density, gas vertical velocity and gas temperature.

Fig. 14 .
Fig.14.Case C. Time-averaged axial profiles of gas-particle mechanical and thermal nonequilibrium computed over the time interval (1.6, 2.0) s.Difference between gas vertical velocity and particle vertical velocities (left) and difference between gas temperature and particle temperatures (right).

Fig. 16 .
Fig. 16.(a) Gas vertical velocity [m s −1 ] and (b) logarithm to the base 10 of particle volume fraction.Snapshots at t = 2 s of a vertical section of quantities averaged along the azimuthal angle.

Table 1 .
Inlet conditions for a homogeneous underexpanded jet at the laboratory scale.

Table 2 .
Inlet conditions for a inhomogeneous underexpanded jet at the laboratory scale.

Table 3 .
Inlet conditions of the inhomogeneous underexpanded jets.
k T s − T g