ASHEE: a compressible, equilibrium-Eulerian model for volcanic ash plumes

A new fluid-dynamic model is developed to numerically simulate the non-equilibrium dynamics of polydisperse gas-particle mixtures forming volcanic plumes. Starting from the three-dimensional N-phase Eulerian transport equations for a mixture of gases and solid particles, we adopt an asymptotic expansion strategy to derive a compressible version of the first-order non-equilibrium model, valid for low concentration regimes and small particles Stokes $St<0.2$. When $St<0.001$ the model reduces to the dusty-gas one. The new model is significantly faster than the Eulerian model while retaining the capability to describe gas-particle non-equilibrium. Direct numerical simulation accurately reproduce the dynamics of isotropic turbulence in subsonic regime. For gas-particle mixtures, it describes the main features of density fluctuations and the preferential concentration of particles by turbulence, verifying the model reliability and suitability for the simulation of high-Reynolds number and high-temperature regimes. On the other hand, Large-Eddy Numerical Simulations of forced plumes are able to reproduce their observed averaged and instantaneous properties. The self-similar radial profile and the development of large-scale structures are reproduced, including the rate of entrainment of atmospheric air. Application to the Large-Eddy Simulation of the injection of the eruptive mixture in a stratified atmosphere describes some of important features of turbulent volcanic plumes, including air entrainment, buoyancy reversal, and maximum plume height. Coarse particles partially decouple from the gas within eddies, modifying the turbulent structure, and preferentially concentrate at the eddy periphery, eventually being lost from the plume margins due to the gravity. By these mechanisms, gas-particle non-equilibrium is able to influence the large-scale behavior of volcanic plumes.


Introduction
Explosive volcanic eruptions are characterized by the injection from a vent into the atmosphere of a mixture of gases, liquid droplets and solid particles, at high velocity and temperature. In typical magmatic eruptions, solid particles constitute more than 95% of the erupted mass and are mostly produced by the brittle fragmentation of a highly viscous magma during its rapid ascent in a narrow conduit (Wilson, 1976;Sparks, 1978), with particle sizes and densities spanning over a wide range, depending on the overall character and intensity of the eruption (Kaminski and Jaupart, 1998;Kueppers et al., 2006). The order of magnitude of the plume mixture volumetric concentration very rarely exceed s ∼ 3 * 10 −3 , because the order of magnitude of the ejected fragments density isρ s ∼ 10 3 kg/m 3 . Thus, the plume mixture con be considered mainly as a dilute suspension in the sense of Elghobashi (1991Elghobashi ( , 1994. This threshold for s is overcome in the dense layer forming in presence of pyroclastic density currents (see e.g. Orsucci, 2014). In the literature, collisions between ash particles are usually disregarded when looking at the dynamics of volcanic ash plume, because this dilute character of the plume mixture (cf. Morton et al., 1956;Woods, 2010).
After injection in the atmosphere, this multiphase eruptive mixture can rise convectively in the atmosphere, forming ei-ther a buoyant volcanic plume or collapse catastrophically forming pyroclastic density currents. Since these two endmembers have different spatial and temporal scales and different impacts on the surrounding of a volcano, understanding the dynamics of volcanic columns and the mechanism of this bifurcation is one of the topical aims of volcanology and one of the main motivations for this work.
The term volcanic column will be adopted in this paper to generically indicate the eruptive character (e.g. convective/collapsing column). Following the fluid-dynamic nomenclature, we will term jet the inertial regime of the volcanic column and plume the buoyancy-driven regime. A forced plume is characterized by an initial momentum-driven jet stage, transitioning into a plume.
In this work, we present a new computational fluiddynamic model to simulate turbulent gas-particle forced plumes in the atmosphere. Although the focus of the paper is on multiphase turbulence in subsonic regimes, the model is also suited to be applied to transonic and supersonic flows. In many cases, indeed, the eruptive mixture is injected in the atmosphere at pressure higher than atmospheric, so that the flow is initially driven by a rapid, transonic decompression stage. This is suggested by numerical models predicting choked flow conditions at the volcanic vent (Wilson, 1980;Wilson et al., 1980), implying a supersonic transition above the vent or in the crater (Kieffer, 1984;Woods and Bower, 1995;Koyaguchi et al., 2010) and it is supported by field evidences of the emission of shock waves during the initial stages of an eruptions (Morrissey, 1997). Despite the importance of the decompression stage for the subsequent development of the volcanic plume (Pelanti and LeVeque, 2006;Ogden et al., 2008b;Orescanin et al., 2010;Carcano et al., 2013) and for the stability of the eruptive column (Ogden et al., 2008a), our analysis is limited to the plume region where flow pressure is equilibrated to the atmospheric pressure. From laboratory experiments, this is expected to occur within less than 20 inlet diameters above the ground (Yüceil andÖtügen, 2002).

Dusty gas modeling of volcanic plumes
Starting from the assumption that the large-scale behavior of volcanic columns is controlled by the bulk properties of the eruptive mixture, most of the previous models of volcanic plumes have considered the eruptive mixture as homogeneous (i.e., they assume that particles are perfectly coupled to the gas phase). Under such hypothesis, the multiphase transport equations can be largely simplified and reduce to a set of mass, momentum and energy balance equations for a single fluid (named dusty-gas or pseudo-gas) having average thermo-fluid dynamic properties (mixture density, velocity and temperature) and equation of states accounting for the incompressibility of the particulate phase and gas covolume (Marble, 1970).
By adopting the dusty gas approximation, volcanic plumes have been studied in the framework of jet (Prandtl, 1963) and plume theory (Morton et al., 1956;Morton, 1959). Onedimensional, steady-state pseudo-gas models of volcanic plumes have thus had a formidable role in volcanology to identify the main processes controlling their dynamics and scaling properties (Wilson, 1976;Woods, 1988;Sparks et al., 1997).
Accordingly, volcanic plume dynamics is schematically subdivided into two main stages. The lower, jet phase is driven by the initial flow momentum. Mixture buoyancy is initially negative (the bulk density is larger than atmospheric) but the mixture progressively expands adiabatically thanks to atmospheric air entrainment and heating, eventually undergoing a buoyancy reversal. When buoyancy reversal does not occur, partial or total collapse of the jet from its maximum thrust height (where the jet has lost all its initial momentum) and generation of pyroclastic density currents are expected.
Above the jet thrust region, the rise of volcanic plumes is driven by buoyancy and it is controlled by turbulent mixing until, in the stratified atmosphere, a level of neutral buoyancy is reached. Above that height, the plume starts to spread out achieving its maximum height and forming an umbrella ash cloud, dispersing in the atmosphere and slowly falling-out.
In one-dimensional, time-averaged models, entrainment of atmospheric air is described by one empirical coefficient (the entrainment coefficient) relating the influx of atmospheric air to the local, vertical plume velocity. The entrainment coefficient also determines the plume shape (Ishimine, 2006) and can be empirically determined by means of direct field observations or ad-hoc laboratory measurements.
Further development of volcanic plume models have included the influence of atmospheric stratification and humidity (Woods, 1993;Glaze and Baloga, 1996), the effect of cross wind (Bursik, 2001), loss and reentrainment of solid particles from plume margins (Woods and Bursik, 1991;Veitch and Woods, 2002), and transient effects (Scase, 2009;Woodhouse et al., 2015). However, one-dimensional models strongly rely on the self-similarity hypothesis, whose validity cannot be experimentally ascertained for volcanic eruptions.
To overcome the limitations of one-dimensional models, three-dimensional dusty-gas models have been developed to simulate volcanic plumes. Suzuki (2005) have developed a three-dimensional dusty gas model (SK-3D) able to accurately resolve the relevant turbulent scales of a volcanic plume, allowing a first, theoretical determination of the entrainment coefficient , without the need of an empirical calibration.
To simulate the three-dimensional large-scale dynamics of volcanic plumes including particle settling and the complex microphysics of water in volcanic plumes, the ATHAM (Active Tracer High Resolution Atmospheric Model) code has been designed (Oberhuber et al., 1998;Graf et al., 1999;Van Eaton et al., 2015). ATHAM describes the dynamics of gasparticle mixtures by assuming that particles are in kinetic equilibrium with the gas phase only in the horizontal component, whereas along the vertical direction they are allowed to have a differential velocity. Thermal equilibrium is assumed. In this sense, ATHAM relaxes the dusty-gas approximation (while maintaining its fundamental structure and the same momentum transport equations) by describing the settling of particles with respect to the gas.

Multiphase flow models of volcanic plumes
Notwithstanding all the above advantages, dusty-gas models are still limited by the equilibrium assumption, which can be questionable at least for the coarsest part of the granulometric spectrum in a plume. Turbulence is indeed a non-linear, multiscale process and the time and space scales of gas-particle interaction may be comparable with some relevant turbulent scales, thus influencing the large-scale behavior of volcanic plumes.
To model non-equilibrium processes, Eulerian multiphase flow models have been developed, which solve the full set of mass, momentum, and energy transport equations for a mixture of gas and dispersed particles, treated as interpenetrating fluids. Valentine and Wohletz (1989) and Dobran et al. (1993); Neri and Dobran (1994) have first analyzed the influence of erupting parameters on the column behavior to identify: By means of two-dimensional numerical simulations, they individuate a threshold from collapsing to convective columns. Lately, two-dimensional (Di Muro et al., 2004;Dartevelle et al., 2004) and three-dimensional numerical simulations (Esposti Ongaro et al., 2008) has contributed to modify the view of a sharp transition between convecting and collapsing columns in favor of that of a transitional regime, characterized by a progressively increasing fraction of mass collapsing. However, previous works could not investigate in detail the non-equilibrium effects in volcanic plumes, mainly because of their averaged description of turbulence: a detailed resolution of the relevant turbulent scales in three dimensions would indeed be computationally prohibitive for N-phase systems.
The main objective of the present work is therefore to develop a new physical model and a fast three-dimensional numerical code able to resolve the spatial and temporal scales of the interaction between gas and particles in turbulent regime and to describe the kinetic non-equilibrium dynamics and their influence on the observable features of volcanic plumes. To this aim, a development of the so-called equilibrium-Eulerian approach (Ferry and Balachandar, 2001;Balachandar and Eaton, 2010) has been adopted. It is a generalization of the dusty-gas model keeping the kinematic nonequilibrium as a first order correction of the Marble (1970) model with respect to the Stokes number of the solid particles/bubbles in the mixture.
The derivation of the fluid dynamic model describing the non-equilibrium gas-particle mixture is described in detail in Section 2. The computational solution procedure and the nu-merical code development are reported in Section 3. Section 4 focuses on verification and validation issues in the context of applications to turbulent volcanic plumes. In particular, here we discuss: three-dimensional numerical simulations of compressible, isotropic turbulence (with and without particles); experimental scale forced plumes; Sod's shock tube problem. Finally, Section 5 presents numerical simulations of volcanic plumes and discusses some aspects related to numerical grid resolution in practical cases.

The multiphase flow model
To derive an appropriate multiphase flow model to describe gas-particle volcanic plumes, we here introduce the nondimensional scaling parameters characterizing gas particle and particle particle interactions.
The drag force between gas and particles introduces in the system a time scale τ s , the particle relaxation time, which is the time needed to a particle to equilibrate to a change of gas velocity. Gas-particle drag is a non-linear function of the local flow variables and, in particular, it depends strongly on the relative Reynolds number, defined as: here d s is the particle diameter,ρ g is the gas density, µ is the gas dynamic viscosity coefficient and u g(s) is the gas (solid) phase velocity field. Beingρ g(s) the gaseous (solid) phase density and s = V s /V the volumetric concentration of the solid phase, it is useful to define the gas bulk density ρ g ≡ (1 − s )ρ g ρ g and the solid bulk density ρ s ≡ sρs (even though in our applications s is order 10 −3 , ρ s is nonnegligible sinceρ s /ρ g is of order 10 3 ).
For an individual point-like particle (i.e., having diameter d s much smaller than the scale of the problem under analysis), at Re s < 1000, the drag force per volume unity can be given by the Stokes' law: where is the characteristic time of particle velocity relaxation with respect the gas,ρ s is the particle density, ν is the gas kinematic viscosity and φ c = 1 + 0.15Re 0.687 s is a correction factor (obtained from the Schiller-Naumann correlation) for finite particle Reynolds number (cf. Clift et al., 1978;Balachandar, 2009;Balachandar and Eaton, 2010;Cerminara, 2015b). In Eq. (2) we disregard all the effects due to the pressure gradient, the added mass, the Basset history and the Saffman terms, because we are considering heavy particles: ρ s /ρ g 1 (cf. Ferry and Balachandar, 2001;Bagheri et al., 2013). Equation (2) has a linear dependence on the fluidparticle relative velocity only when Re s 1, so that φ c 1 and the classic Stokes drag expression is recovered. On the other hand, if the relative Reynolds number Re s grows, nonlinear effects become much more important in Eq. (3). The Clift et al. (1978) empirical relationship used in this work has been used and tested in a number of papers (e.g., Balachandar and Eaton, 2010;Wang and Maxey, 1993;Bonadonna et al., 2002), and it is equivalent to assuming the following gas-particle drag coefficient: Wang and Maxey (1993) discussed nonlinear effects due to this correction on the dynamics of point-like particles falling under gravity in an homogeneous and isotropic turbulent surrounding. We recall here the terminal velocity that can be found by setting u g = 0 in Eq. (2) is: As previously pointed out, correction used in Eq. (4) is valid if Re s < 10 3 , the regime addressed in this work for ash particles much denser then the surrounding fluid and smaller than 1mm. As shown by Balachandar (2009), maximum values of Re s are associated to particle gravitational settling (not to turbulence). Using formula (4) and (5), it is thus possible to estimate Re s of a falling particle with diameter d s . We obtain that Re s is always smaller than 10 3 for ash particles finer than 1 mm in air. If regimes with a bigger decoupling needs to be explored, more complex empirical correction has to be used for φ c (Neri et al., 2003;Bürger and Wendland, 2001). The same reasoning can be applied to estimate the thermal relaxation time between gas and particles. In terms of the solid phase specific heat capacity C s and its thermal conductivity k g , we have: where Nu s = Nu s (Re s ,Pr) is the Nusselt number, usually function of the relative Reynolds number and of the Prandtl number of the carrier fluid (Neri et al., 2003). In terms of τ T , the heat exchange between a particle at temperature T s and the surrounding gas at temperature T g per unit volume is: Comparing the kinetic and thermal relaxation times we get: In order to estimate this number, firstly we notice that factor 2φ c /Nu s tends to 1 if Re s → 0, and it remains smaller than 2 if Re s < 10 3 (Neri et al., 2003;Cerminara, 2015b). Then, in the case of ash particles in air, we have (in SI units) µ 10 −5 , C s 10 3 , k g 10 −2 . Thus we have that τ T /τ s 1, meaning that the thermal equilibrium time is typically of the same order of magnitude of the kinematic one. This bound is very useful when we write the equilibrium-Eulerian and the dusty gas models, because it tells us that the thermal Stokes number is of the same order of the kinematic one, at least for volcanic ash finer than 1 mm.
The non-dimensional Stokes number (St) is defined as the ratio between the kinetic (or thermal) relaxation time and a characteristic time of the flow under investigation τ F , namely St s = τ s /τ F . The definition of the flow time-scale can be problematic for high-Reynolds number flows (typical of volcanic plumes), which are characterized by a wide range of interacting length-and time-scales, a distinctive feature of the turbulent regime. For volcanic plumes, the more energetic time-scale would be of the order τ L = L/U , where L and U are the plume diameter and velocity at the vent, which gives the characteristic turnover time of the largest eddies in a turbulent plume (e.g., Zhou et al., 2001). On the other hand, the smallest time-scale (largest St s ) can be defined by the Kolmogorov similarity law by τ η ∼ τ L Re −1/2 L , where the macroscopic Reynolds number is defined, at a first instance, by Re L = U L/ν, ν being the kinematic viscosity of the gas phase alone. For numerical models, it is also useful introducing the Large-Eddy Simulation (LES) time-scale τ ξ , with respect to the resolved scales ξ, related to the numerical grid resolution, size of the explicit filter, and discretization accuracy (Lesieur et al., 2005;Garnier et al., 2009;Balachandar and Eaton, 2010;Cerminara et al., 2015). At LES scale ξ, St s is not as large as at the Kolmogorov scale, thus the decoupling between particles and the carrier fluid is mitigated by the LES filtering operation. We found that St s 0.2 for LES of volcanic ash finer than 1 mm.
The model presented here is conceived for resolving dilute suspensions, namely mixtures of gases and particles with volumetric concentration Vs V ≡ s 10 −3 . We here use the definition of dilute suspension by Elghobashi (1991Elghobashi ( , 1994 and Balachandar (2009), corresponding to regimes in which particle-particle collisions can be disregarded. This can also be justified by analyzing the time scale of particle-particle collisions. In the dilute regime, in which we can assume an equilibrium Maxwell distribution of particle velocities, the mean free path of solid particles is given by (Gidaspow, 1994): Consequently, particle-particle collisions are relatively infrequent (λ p-p ∼ 0.1 m d s ), so that we can neglect, as a first approximation, particle-particle collisions and consider the particulate fluid as pressure-less, inviscid and nonconductive. In volcanic plumes the particle volumetric concentration can exceed of one order of magnitude the threshold s 10 −3 only near the vent (see, e.g., Sparks et al., 1997;?). However, the region of the plume where the dilute suspension requirement is not fulfilled remains small with respect the size of the entire plume, weakly influencing its global dynamics. Indeed, as we will show in Sec. 5, air entrainment and particle fallout induce a rapid decrease of the volumetric concentration. On the contrary, the mass fraction of the solid particles can not be considered small, because particles are heavy: s * ρ s ≡ ρ s ρ g . Thus, particles inertia will be considered in the present model: in other words, we will consider the two way coupling between dispersed particles and the carrier gas phase.
Summarizing, our multiphase model focuses and carefully takes advantage of the hypotheses characterizing the following regimes: heavy particles (ρ s /ρ g 1) in dilute suspension ( s 10 −3 ) with dynamical length scales much larger than the particles diameter (point-particle approach) and relative Reynolds number smaller than 10 3 .

Eulerian-Eulerian multiphase flow model
When the Stokes number is smaller than one, and the number of particles is very large, it is convenient to use an Eulerian approach, where the carrier and the dispersed phase are modeled as interpenetrating continua, and their dynamics is described by the laws of fluid mechanics (Balachandar and Eaton, 2010).
Here we want to model a polydisperse mixture of i ∈ [1,2,...,I] ≡ I gaseous phases and j ∈ [1,2,...,J] ≡ J solid phases. From now on, we will use the subscript (·) j instead of (·) s for the jth solid phase. Solid phases represent the discretization of a virtually continuous grain-size distribution into discrete bins, as usually done in volcanological studies (cf. Cioni et al., 2003;Neri et al., 2003). Another possible approach is the method of moments, in which the evolution of the moments of the grain size distribution is described. This has recently been applied in volcanology to integral plume models by de' Michieli Vitturi et al. (2015). In the present work we opted for the classical discretization of the grain size distribution (cf. Neri et al., 2003). In (Cerminara, 2015b) we analyze the Eulerian-Eulerian model under the barotropic regime to show the existence of weak solutions of the corresponding partial differential equations problem.
In the regime described above, the Eulerian-Eulerian equations for a mixture of a gas and a solid dispersed phase are (Feireisl, 2004;Marble, 1970;Neri et al., 2003;Gidaspow, 1994;Garnier et al., 2009;Berselli et al., 2015;?): with the following constitutive equations (g is the gravitational acceleration): -Given y i(j) the mass fractions of the gaseous (solid) phases and ρ m the bulk density of the mixture, the bulk density of the gas phase is ρ g = I ρ i = I y i ρ m , while the mass fraction of the solid phases ρ s = J ρ j = J y j ρ m . Consequently, ρ m = ρ g + ρ s .
-The volumetric concentration of the ith(jth) phase is given by i = ρ i /ρ i .
-Perfect gas: p = Iρ i R i T g , with R i the gas constant of the ith gas phase. This law can be simplified by nothing that s 1, thus i 1 andρ i ρ i (cf. Suzuki, 2005). Anyway, in this work we use the complete version of the perfect gas law. It can be written in convenient form for a poly-disperse mixture as: -Newtonian gas stress tensor: where µ(T ) = I i µ i (T ) is the gas dynamic viscosity and µ i is that of the ith gas component.
-Enthalpy per unit of mass of the gas (solid) phase: h g = I ρ i C p,i T g /ρ g h j = C j T j , with C p,i C j the specific heat at constant pressure of the ith (jth) phase.
-The Fourier law for the heat transfer in the gas: q = −k g ∇T , where k g = I i k i and k i is the conductivity of the ith gas component.
-Q j and f j refer to Q T and f d of the jth solid phase; S j is the source or sink term (when needed) of the jth phase. K i = |u i | 2 /2 is the kinetic energy per unit of mass of the ith gas phase (K j for the jth solid phase).

Equilibrium-Eulerian model
In the limit St j 1, the drag terms f j and the thermal exchange terms Q j can be calculated by knowing u g and T g , and the Eulerian-Eulerian model can be largely simplified by considering the dusty-gas (also known as pseudo-gas) approximation (Marble, 1970). A refinement of this approximation (valid if St j 0.2), has been developed by Maxey (1987), as a first-order approximation of the Lagrangian particle momentum balance (see Eq. (10)d): By using the Stokes law and a perturbation method, and by defining a ≡ D t u g (with D t = ∂ ∂t + u · ∇ ), we obtain a correction to particle velocity up to first order leading to the so-called equilibrium-Eulerian model developed by Ferry and Balachandar (2001) and Balachandar and Eaton (2010) for incompressible multiphase flows. It is worth noting that at the zeroth order we recover To write the compressible version of that model, we define the relative jth particle velocity field v j so that u j = u g +v j .
Recalling the definitions of the mass fraction and the mixture density given above, we define: By summing up the gas momentum equation with the solid momentum equations in Eq. (10), we thus obtain: This momentum balance equation is equivalent to the compressible Navier-Stokes equation with the substitution u g → u m and the addition of the term ∇ · (ρ m T r ) which takes into account the first order effects of particle decoupling on momentum (two-way coupling). We keep this term because of the presence of the settling velocity w j in v j which is at the leading order.
Moving to the mass conservation, summing up over i and j the continuity equations in (10), we obtain the continuity equation for the mixture: while for the phases we have: It is worth noting that the mixture density follows the classical continuity equation with velocity field u m . We refer to u m as the mixture velocity field. As pointed out in Eq. (8) and below, in our physical regime the thermal Stokes time is of the same order of magnitude of the kinematic one. However, this regime has been thoroughly analyzed in the incompressible case by Ferry and Balachandar (2005), demonstrating that the error made by assuming thermal equilibrium is at least one order of magnitude smaller than that on the momentum equation (at equal Stokes number), thus justifying the limit T j → T g = T as done for the thermal equation in the dusty gas model.
By summing up all enthalpy equations in (10), and by defining h m = I y i h i + J y j h j = C m T and K m = take into account the combined effect of the kinematic decoupling and the difference between the specific heat (v C ) and kinetic energy (v K ) of the mixture and of the jth specie. Summarizing, the compressible equilibrium-Eulerian model is: (25) The first equation is redundant, because it is contained in the second and third set of continuity equations. Note that we have not used the explicit form of v j for deriving Eqs. (25), which therefore can be used for any multiphase flow model with I phases moving with velocity u g and temperature T , and J phases each moving with velocity u j = u g + v j and temperature T . However, in what follows we will use Eq. (14) when referring to the compressible Equilibrium-Eulerian model.
It is also worth noting that in the Navier-Stokes equations it is critical to accurately take into account the nonlinear terms contained by the conservative derivative ∂ t ψ + ∇ · (ψu) because they are the origin of the major difficulties in turbulence modeling. A large advantage of the dusty gas and equilibrium-Eulerian models is that in both models the the most relevant part of the drag ( J f j ) and heat exchange ( J Q j ) terms have been absorbed into the conservative derivatives for the mixture. This fact allows the numerical solver to implicitly and accurately solve the particles contribution on mixture momentum and energy (two-way coupling), using the same numerical techniques developed in Computational Fluid Dynamics for the Navier-Stokes equations. The dusty gas and Equilibrium-Eulerian models are best suited for solving multiphase system in which the particles are strongly coupled with the carrier fluid and the bulk density of the particles is not negligible with respect to that of the fluid.
The equilibrium-Eulerian model thus reduces to a set of mass, momentum, and energy balance equations for the gasparticle mixture plus one equation for the mass transport of the particulate phase. In this respect, it is similar to the dustygas equations, to which it reduces for τ s ≡ 0. With respect to the dusty-gas model, here we solve for the mixture velocity u m , which is slightly different from the carrier gas velocity u g . Moreover, here kinematic decoupling is taken into account by moving the I gas phases and the J solid phases with different velocity fields, respectively u g and u j . Thus, we are accounting for the imperfect coupling of the particles to the gas flow, leading to preferential concentration and settling phenomena (the vector v j includes a convective and a gravity accelerations terms).
The equilibrium-Eulerian method becomes even more efficient (relative to the standard Eulerian) for the polydisperse case (J > 1). For each species of particle tracked, the standard Eulerian method requires four scalar fields, the fast method require one. Furthermore, the computation of the correction to v j needs only to be done for one particle species. The correction has the form −τ j a, so once the term a is computed, velocities for all species of particles may be obtained simply by scaling the correction factor based on the species' response times τ j . To be more precise, the standard Eulerian method needs I + 5J + 4 scalar partial differential equations, while the equilibrium-Eulerian model needs just I + J + 4, i.e. 4J equations less.

Sub-grid scale models
The spectrum of the density, velocity and temperature fluctuations of turbulent flows at high Reynolds number typically span over many orders of magnitude. In the cases where the turbulent spectrum extend beyond the numerical grid resolution, it is necessary to model the effects of the high-frequency fluctuations (those that cannot be resolved by the numerical grid) on the resolved flow. This leads to the so-called Large-Eddy Simulation (LES) techniques, in which a low-pass filter is applied to the model equations to filter out the small scales of the solution. In the incompressible case the theory is welldeveloped (see Berselli et al., 2005;Sagaut, 2006), but LES for compressible flows is still a new and open research field. In our case, we apply a spatial filter, denoted by an overbar (δ is the filter scale): Some example of LES filters G(x;δ) used in compressible turbulence are reviewed in Garnier et al. (2009). In compressible turbulence it is also useful to introduce the so-called Favre filter:ψ First, we apply this filter to the Equilibrium-Eulerian model fundamental equation (14) modified as follows: moving all the new second order terms into O(τ 2 j ), using ∂ t y j + u j · ∇y j = 0 and defining: Multiplying the new expression for u j by ρ m and Favrefiltering, at the first order we obtain: where we have usedτ j = τ j and consequentlyw j = w j because the Stokes time changes only at the large scale and it can be considered constant at the filter scale. Moreover, we have defined the subgrid-scale Reynolds stress tensor: As discussed and tested in Shotorban and Balachandar (2007), the subgrid terms can be considered O(τ j ) and neglected when multiplied by first order terms. Another form of Eq. (30) can be recovered by noting that at the leading orderũ m ũ g −w r : (32) We recall here the Boussinesq eddy viscosity hypothesis: where the deviatoric part of the subgrid stress tensor can be modeled with an eddy viscosity µ t times the rate-of-shear tensorS m = sym(∇ũ m ) − 1 3 ∇ ·ũ m I. The first term on the right hand side of Eq. (33) is the isotropic part of the subgridscale tensor, proportional to the subgrid-scale kinetic energy K t . While in incompressible turbulence the latter term is absorbed into the pressure, it must be modeled for compressible flows (cf. Moin et al. (1991) and Yoshizawa (1986)). Ducros et al. (1995) showed another way to treat this term by absorbing it into a new macro-pressure and macro-temperature (cf. also Lesieur et al. (2005) and Lodato et al. (2009)). We recall here also the eddy diffusivity viscosity model (cf. Moin et al. (1991)): any scalar ψ transported by u m generates a subgrid-scale vector that can be modeled with the large eddy variables. We have: where Pr t is the subgrid-scale turbulent Prandtl number. We apply the Favre filter defined in Eq. (27) to Eqs. (25) (for the application of the Favre filter to the compressible Navier-Stokes equations cf. Garnier et al. (2009), Moin et al. (1991) and Erlebacher et al. (1990)), obtaining: where are: the subgrid eddy diffusivity vector of the ith phase; of the jth phase; the subgrid-scale stress tensor; the diffusivity vector of the temperature; respectively. Other approximations have been used to derive the former LES model: the viscous terms in momentum and energy equations, and the pressure-dilatation and conduction terms in the energy equations are all non-linear terms and we here treat them as done by Erlebacher et al. (1990) and Moin et al. (1991). The subgrid terms corresponding to the former non-linear terms could be neglected so that, for example, p∇ · u g p∇ ·ũ g .
In particular, this term has been neglected also in presence of shocks (cf. Garnier et al. (2002)). We refer to Vreman (1995) for an a priori and a posteriori analysis of all the neglected terms of the compressible Navier-Stokes equations.
Moreover, in our model the mixture specific heat C m and the mixture gas constant R m vary in the domain because y i and y j vary. Thus, also the following approximations should be done, coherently with the other approximations used:h m = C m T C mT and R m T R mT . In order to close the system, terms µ t , K t and Pr t must be chosen on the basis of LES models, either static or dynamic (see Moin et al., 1991;Bardina et al., 1980;Germano et al., 1991). In the present model, we implemented several sub-grid scale (SGS) models to compute the SGS viscosity, kinetic energy and Prandtl number (Cerminara, 2015b). Currently, the code offers the possibility of choosing between: 1) the compressible Smagorinsky model, both static and dynamic (see Fureby, 1996;Yoshizawa, 1993;Pope, 2000;Chai and Mahesh, 2012;Garnier et al., 2009); 2) the subgrid scale K-equation model, both static and dynamic (see Chacón-Rebollo and Lewandowski, 2013;Fureby, 1996;Yoshizawa, 1993;Chai and Mahesh, 2012); 3) the dynamical Smagorinsky model in the form by Moin et al. (1991); 4) the WALE model, both static and dynamic (see Nicoud and Ducros, 1999;Lodato et al., 2009;Piscaglia et al., 2013).
All through this paper, we present results obtained with the dynamic WALE model (see Fig. 5 and the corresponding section for a study on the accuracy of this LES model). A detailed analysis of the influence of subgrid-scale models to simulation results is beyond the scopes of this paper and will be addressed in future works.

Numerical solver
The Eulerian model described in Section 2, is solved numerically to obtain a time-dependent description of all independent flow fields in a three-dimensional domain with prescribed initial and boundary conditions. We have chosen to adopt an open-source approach to the code development in order to guarantee control on the numerical solution procedure and to share scientific knowledge. We hope that this will help building a wider computational volcanology community. As a platform for developing our solver, we have chosen the unstructured, finite volume (FV) method based open source C++ library OpenFOAM (version 2.1.1). OpenFOAM, released under the Gnu Public License (GPL), has gained a vast popularity during the recent years. The readily existing solvers and tutorials provide a quick start to using the code also to inexperienced users. Thanks to a high level of abstraction in the programming of the source code, the existing solvers can be freely and easily modified in order to create new solvers (e.g., to solve a different set of equations) and/or to implement new numerical schemes. OpenFOAM is well integrated with advanced tools for preprocessing (including meshing) and post-processing (including visualization). The support of the OpenCFD Ltd, of the OpenFOAM foundation and of a wide developers and users community guarantees ease of implementation, maintenance and extension, suited for satisfying the needs of both volcanology researchers and of potential users, e.g. in volcano observatories. Finally, all solvers can be run in parallel on distributed memory architectures, which makes Open-FOAM suited for envisaging large-scale, three-dimensional volcanological problems.

Finite Volume discretization strategy
In the FV method (Ferziger and Perić, 1996), the governing partial differential equations are integrated over a computational cell, and the Gauss theorem is applied to convert the volume integrals into surface integrals, involving surface fluxes. Reconstruction of scalar and vector fields (which are defined in the cell centroid) on the cell interface is a key step in the FV method, controlling both the accuracy and the stability properties of the numerical method.
OpenFOAM implements a wide choice of discretization schemes. In all our test cases, the temporal discretization is based on the second-order Crank-Nicolson scheme (Ferziger and Perić, 1996), with a blending factor of 0.5 (0 meaning a first-order Euler scheme, 1 a second-order, bounded implicit scheme) and an adaptive time stepping based on the maximum initial residual of the previous time step (Kay et al., 2010), and on a threshold that depends on the Courant number (C < 0.2). All advection terms of the model are treated implicitly to enforce stability. Diffusion terms are also discretized implicit in time, with the exception of those representing subgrid turbulence. The pressure and gravity terms in the momentum equations and the continuity equations are solved explicitly. However, as discussed below, the PISO (Pressure Implicit with Splitting of Operators, Issa (1986)) solution procedure based on a pressure correction algorithm makes such a coupling implicit. Similarly, the pressure advection terms in the enthalpy equation and the relative veloc-ity v j are made implicit when the PIMPLE (mixed SIMPLE and PISO algorithm, Ferziger and Perić (1996)) procedure is adopted. The same PIMPLE scheme is applied treating all source terms and the additional terms deriving from the equilibrium-Eulerian expansion.
In all described test cases, the spatial gradients are discretized by adopting an unlimited centered linear scheme (which is second-order accurate and has low numerical diffusion -Ferziger and Perić, 1996). Analogously, implicit advective fluxes at the control volume interfaces are reconstructed by using a centered linear interpolation scheme (also second order accurate). The only exception is for pressure fluxes in the pressure correction equation, for which we adopt a TVD (Total Variation Diminishing) limited linear scheme (in the subsonic regimes) to enforce stability and non-oscillatory behavior of the solution. We refer to Jasak (1996) for a complete description of the discretization strategy adopted in OpenFOAM.

Solution procedure
Instead of solving the set of algebraic equations deriving from the discretization procedure as a whole, most of the existing solvers in OpenFOAM are based on a segregated solution strategy, in which partial differential equations are solved sequentially and their coupling is resolved by iterating the solution procedure. In particular, for Eulerian fluid equations, momentum and continuity equation (coupled through the pressure gradient term and the gas equation of state) are solved by adopting the PISO algorithm. The PISO algorithm consists of one predictor step, where an intermediate velocity field is solved using pressure from the previous time-step, and of a number of PISO corrector steps, where intermediate and final velocity and pressure fields are obtained iteratively. The number of corrector steps used affects the solution accuracy and usually at least two steps are used. Additionally, coupling of the energy (or enthalpy) equation can be achieved in OpenFOAM through additional PIMPLE iterations (which derives from the SIMPLE algorithm by Patankar, 1980). For each transport equation, the linearized system deriving from the implicit treatment of the advectiondiffusion terms is solved by using the PbiCG solver (Preconditioned bi-Conjugate Gradient solver for asymmetric matrices) and the PCG (Preconditioned Conjugate Gradient solver for symmetric matrices), respectively, preconditioned by a Diagonal Incomplete Lower Upper decomposition (DILU) and a Diagonal Incomplete Cholesky (DIC) decomposition. The segregated system is iteratively solved until a global tolerance threshold PIMPLE is achieved. In our simulations, we typically use PIMPLE < 10 −7 for this threshold.
The numerical solution algorithm is designed as follows: 1. Solve the (explicit) continuity equation (35) for mixture density ρ m (predictor stage: uses fluxes from previous iteration).
3. Solve the (semi-implicit) momentum equation to obtain u m (predictor stage: uses the pressure field from previous iteration).
4. Solve the (semi-implicit) enthalpy equation to update the temperature field T and the compressibility ρ m /p (pressure from previous iteration).
5. Solve the (implicit) pressure equation to update the pressure (uses predicted values of fluxes).
6. Correct density, velocity and fluxes with the new pressure field (keeping T and ρ m /p fixed).
7. Iterate from 5 evaluating the continuity error as the difference between the kinematic and thermodynamic calculation of the density (PISO loop).
9. Evaluate the numerical error PIMPLE and iterate from 2 if prescribed (PIMPLE loop).

Compute LES subgrid terms.
With respect to the standard solvers implemented in OpenFOAM (v2.1.1) for compressible fluid flows (e.g. sonicFoam or rhoPimpleFoam), the main modification required are the following: 1. The mixture density and velocity replaces the fluid ones.

2.
A new scalar transport equation is introduced for the mass fraction of each particulate and gas species.
3. The equations of state are modified as described in Eq.(11).
4. First-order terms from the equilibrium-Eulerian model are added in the mass, momentum and enthalpy equations.
5. Equations are added to compute flow acceleration and velocity disequilibrium.
6. Gravity terms and ambient fluid stratification are added.
7. New SGS models are implemented.
Concerning point 5, it is worth remarking that, accordingly to Ferry et al. (2003), the first-order term in τ j in Eq.(14) must be limited to avoid the divergence of preferential concentration in a turbulent flow field (and to keep the effective Stokes number below 0.2). In other word, we impose at each time step that |τ j (a + w j · ∇u g )| ≤ 0.2|u g + w j |. We tested the effect of this limiter on preferential concentration in Sec. 4.2 below.

Verification and validation study
A wide set of numerical tests has been performed to assess the adequacy of the ASHEE model for the intended volcanological application and the reliability of the numerical solution method. Validation tests are focused on the dynamics of gas (Section 4.1) and multiphase (Section 4.2) turbulence and on the mixing properties of buoyant plumes (Section 4.3). Compressibility likely exerts a controlling role to the near-vent dynamics during explosive eruptions (e.g., Carcano et al., 2013). Although this is not the focus of this work, we briefly discuss in Section 4.4 the performance of the model on a standard one-dimensional shock wave numerical test.

Compressible decaying homogeneous and isotropic turbulence
Turbulence is a key process controlling the dynamics of volcanic plumes since it controls the rate of mixing and air entrainment. To assess the capability of the developed model to resolve turbulence (which requires low numerical diffusion and controlled numerical error Geurts and Fröhlich, 2002), we have tested the numerical algorithm against different configurations of decaying homogeneous and isotropic turbulence (DHIT). In this configuration, the flow is initialized in a domain Ω which is a box with side L = 2π with periodic boundary conditions. As described in Lesieur et al. (2005) with peak initially in k = k 0 and so that the initial kinetic energy and enstrophy are: As reviewed by Pope (2000), the Taylor microscale can be written as a function of the dissipation = 2νH: thus in our configuration, the initial Taylor micro scale is: We have chosen the non-dimensionalization keeping the root mean square of the magnitude of velocity fluctuations (u ) equal to u rms : We also chose to make the system dimensionless by fixing ρ m,0 = 1, T 0 = 1, Pr = 1, so that the ideal gas law becomes: and the initial Mach number of the mixture based on the velocity fluctuations reads: This means that Ma rms can be modified keeping fixed u rms and modifying p. Following Honein and Moin (2004), we define the eddy turnover time: The initial compressibility ratio C 0 is defined as the ratio between the kinetic energy and its compressible component K c : Here, u c is the compressible part of the velocity fluctuations, so that ∇ · u = ∇ · u c and ∇ ∧ u c = 0. The last parameter, i.e. the dynamical viscosity, can be given both by fixing the Reynolds number based on λ T,0 or k 0 : It is useful to define the maximum resolved wavenumber on the selected N -cells grid and the Kolmogorov length scale based on Re k0 . They are, respectively: In order to have a DNS, the smallest spatial scale δ should be chosen in order to have k max η K > 2 ( Pirozzoli and Grasso, 2004). We compare the DNS of compressible decaying homogeneous and isotropic turbulence with a reference, well tested numerical solver for Direct Numerical Simulations of compressible turbulence by Pirozzoli and Grasso (2004); Bernardini and Pirozzoli (2009). For this comparison we fix the following initial parameters: p = R m = 1, γ m = 1.4, Pr = 1, Ma rms = 0.2, C 0 = 0, u 2 rms = 2K 0 = 0.056, k 0 = 4, λ T = 0.5, τ e 3.6596, µ = 5.885846 * 10 −4 , Re λ 116, Re k0 100. Thus a grid with N = 256 3 cells gives k max 127 and k max η K = 2π, big enough to have a DNS. The simulation has been performed on 1024 cores on the Fermi Blue Gene/Q infrastructure at Italian CINECA super-computing center (http://www.cineca.it), on which about 5 h are needed to complete the highest-resolution runs (256 3 of cells) up to time t/τ e = 5.465 (about 3500 time-steps). The average required total CPU time on 1024 Fermi cores is about 1-3 millions of cells per second, with the variability associated with the number of solid phases described by the model. This value is confirmed in all benchmark cases presented in this paper. Fig. 1 reports the parallel efficiency on both the Fermi and the PLX (a Linux cluster based on Intel Xeon esa-and quad-core processors at 2.4 GHz) machines at CINECA. The ASHEE code efficiency is very good (above 0.9) up to 512 cores (i.e., up to about 30000 cells per core), but it is overall satisfactory for 1024 cores, with efficiency larger than 0.8 on PLX and slightly lower (about 0.7) on Fermi probably due to the limited level of cache optimization and input/output scalability (Culpo, 2011). The code was run also on 2048 cores on Fermi with parallel efficiency of 0.45 (Dagna, 2013).
The so called Q-criterion (Garnier et al., 2009) allows indeed the identification of coherent vortices inside a three dimensional velocity field. In Fig. 3 we present a comparison of the energy spectrum E(k) obtained with the ASHEE model and the model by Bernardini and Pirozzoli (2009) after approximatively 1 eddy turnover time; the L 2 norm of the difference between the two spectra is 4.0 * 10 −4 . This validates the accuracy of our numerical code in the single-phase and shock-free case.    4 shows the evolution of several integral parameters describing the dynamics of the decaying homogeneous and isotropic turbulence. Fig. 4a displays the density fluctuations ρ rms = (ρ − ρ Ω ) 2 Ω , the density contrast ρ max /ρ min and the standard measure of compressibility C = |∇ · u| 2 Ω / |∇u| 2 Ω which takes value between 0 (incompressible flow) and 1 (potential flow) (Boffetta et al., 2007). All the quantities showed in Fig. 4a depend on the initial Mach number end compressibility. For the case shown, Ma rms = 0.2 and we obtain very similar result to those ported in Fig. 18 and 19 by Garnier et al. (1999). Fig. 4b shows the kinetic energy spectrum at t/τ e = 0, 1.093, 5.465. We notice that the energy spectrum widens from the initial condition until its tail reach k k max 127. Then system becomes to dissipate and the maximum of the energy spectrum decreases. The largest scales tend to lose energy slower than the other scales and the spectrum widens also in the larger scale direction. Fig. 4c presents the evolution of K (total turbulent kinetic energy), H (enstrophy), λ T (Taylor microscale). We notice that the total kinetic energy decreases monotonically and at t 5.5τ e just 15% of its initial value is conserved. On the other hand, enstrophy increases until it reaches a maximum at 1.5 < t/τ e < 2. It then starts to decrease monotonically. This behavior is related to the two different stages we have highlighted in the analysis of the energy spectrum evolution. In the first stage, viscous effects are negligible and enstrophy increases due to vortex stretching. During the second stage, viscous diffusion starts to have an important role and distorted dissipative structures are created (Garnier et al., 1999). Also the Taylor microscale reflects this behavior, reaching a minimum at the end of the first stage and increasing monotonically during the second stage of the evolution. It is a characteristic of the magnitude of the velocity gradients in the inertial range: by comparing it with δ we can have an idea of the broadness of the range of wave numbers where the flow is dissipative. In this DNS, we have λ T 10.2δ at t 5.5τ e .
In Fig. 4d we show the evolution of the Kolmogorov time scale τ K during the evolution of the decaying turbulence.
We finally compare in Fig. 5 the DNS described with simulations at lower resolution with N = 32 3 and N = 64 3 cells. In this case, it is expected that the spectra diverge from the DNS, unless an appropriate subgrid model is introduced to simulate the effects of the unresolved to the resolved scales. Several subgrid models have been tested (Cerminara, 2015b), both static and dynamic. Fig. 5 presents the resulting spectrum using the dynamic WALE model (Nicoud and Ducros, 1999;Lodato et al., 2009). In this figure, we notice how the dynamic WALE model works pretty well for both the 32 3 and 64 3 LES, avoiding the smallest scales to accumulate unphysical energy.

Two-phase isotropic turbulence
In this section we test the capability of our numerical code to correctly describe the decoupling between solid and gaseous phases when St j < 0.2 and to explore its behavior when the equilibrium-Eulerian hypothesis St j < 0.2 is not fulfilled so that a limiter to the relative velocity u g − u j is applied.
To this aim, we performed a numerical simulation of homogeneous and isotropic turbulence with a gas phase initialized with the same initial and geometric conditions described in Sec. 4.1. We added to that configuration 5 solid particle classes (j = 2 ÷ 6) chosen in such a way that St j ∈ [0.03,1],     dar, 2009). We set the material density of all the particles toρ j = 10 3 . In order to have a small contribution of the particle phases to the fluid dynamics -one way couplinghere we set the solid particles mass fraction to a small value, y j = 0.002, so that y g = 0.99.
In Fig. 6 we show a slice of the turbulent box at t/τ e 2.2. Panel a) displays the solid mass fraction, highlighting the preferential concentration and clustering of particles in (a) Mass fraction (b) Acceleration Fig. 6: Slice of the turbulent box at t/τ e 2.2. The two panels represent respectively a logaritmic color map of y 3 (St max = 0.5) and of |a g |.
response to the development of the acceleration field (panel b) associated with turbulent eddies. As described in Maxey (1987); Rani and Balachandar (2003), a good measure for the degree of preferential concentration in incompressible flows is the weighted average on the particle mass fraction of the quantity (|D| 2 − |S| 2 ), where S is the vorticity tensor, i.e. the skew symmetric part of the gas velocity gradient and D is its symmetrical part. For compressible flows, we choose to consider This is a good measure because (use integration by parts, Gauss theorem and Eq. (14) with w j = 0): Moreover, it is worth noting that P j vanishes in absence of preferential concentration. By dimensional analysis, preferential concentration is expected to behave as: because it must be proportional to τ j and have a dimension of [s −2 ]. As described by [Pope 2001], the typical time-scale corresponding to an eddy length-scale ξ in the inertial subrange, can be evaluated by means of the Kolmogorov's theory as:  . We obtain a good agreement between equilibrium-Eulerian LES/DNS and Lagrangian DNS simulations. The fit for the data by Rani and Balachandar (2003) is found in Eq. (66).
where the Taylor microscale λ T is defined by Eq. 47. Since the time based on the Taylor microscale is defined as we can evaluate the typical time at the smallest resolved LES scale ξ knowing the kinetic energy K(t) and λ T (t): In Fig. 7 we show the time-evolution of the degree of preferential concentration as a function of the Stokes number for both DNS with 256 3 cells and the LES with 32 3 cells. There, we multiply P j by τ ξ τ j in order to make it dimensionless and to plot on the same graph all the different particles at different times together.
At t = 0 the preferential concentration is zero for all Stokes number. Then, preferential concentration of each particle class increases up to a maximum value and then it decreases because of the decaying of the turbulent energy. The maximum degree of preferential concentration is reached by each particle class when τ K is minimum (at t/τ e 1.7, cf. Fig. 4d). Then, P j decreases and merges with the curve relative to the next particle class at the final simulation time, when τ K is about twice its minimum. Note that the expected behavior of Eq. (61) is reproduced for St j < 0.2 and in particular we find: Moreover, by comparing our results with the Eulerian-Lagrangian simulation described in Rani and Balachandar (2003), we note that our limiter for the preferential concentration when St > 0.2 is well behaving. For the sake of completeness, we found that the best fit in the range St < 2.5 for the data found by Rani and Balachandar (2003) is: with root mean square of residuals 8.5 * 10 −3 .
For what concerns the 32 3 LES simulation, Fig. 7 shows that the Stokes number of each particle class in the LES case is much smaller than its DNS counterpart. Accordingly with Balachandar and Eaton (2010), we have confirming that the equilibrium-Eulerian model widens its applicability under the LES approximation. We also notice that the presented LES is able to reproduce the expected degree of preferential concentration with a satisfactory level of accuracy when St < 0.2. In particular, the LES slightly overestimates preferential concentration and the time needed to reach the equilibrium and to "forget" the particle initial condition.

Turbulent forced plume
As a second benchmark, we discuss high-resolution, threedimensional numerical simulation of a forced gas plume, produced by the injection of a gas flow from a circular inlet into a stable atmospheric environment at lower temperature (and higher density). Such an experiment allows to test the numerical model behavior against some of the fundamental processes controlling volcanic plumes, namely density variations, non-isotropic turbulence, mixing, air entrainment, and thermal exchange. Our study is mainly aimed at assessing the capability of the numerical model to describe the timeaverage behavior of a turbulent plume and to reproduce the magnitude of large-scale fluctuations and large-eddy structures. We will mainly refer to laboratory experiments by George et al. (1977) and Shabbir and George (1994) and numerical simulations by Zhou et al. (2001) for a quantitative assessment of model results. Numerical simulations describe a vertical round forced plume with heated air as the injection fluid. The plume axis is aligned with the gravity vector and is subjected to a positive buoyancy force. The heat source diameter 2b 0 is 6.35 As discussed by Zhou et al. (2001) the development of the turbulent plume regime is quite sensitive to the inlet conditions: we therefore tested the model by adding a periodic perturbation and a non-homogeneous inlet profile to anticipate the symmetry breaking, and the transition from a laminar to a turbulent flow. The radial profile of vertical velocity has the form: where δ r is the thickness of the turbulent boundary layer at the plume inlet, that we have set at δ r = 0.1b 0 . A periodical forcing and a random perturbation of intensity 0.05U 0 (r) has been superimposed to mimic a turbulent inlet. The resulting average mass, momentum and buoyancy flux are The computational grid is composed of 360 × 180 × 180 uniformly spaced cells (deformed near the bottom plane to conform to the circular inlet) in a box of size 12.8×6.4×6.4 diameters. In particular, the inlet is discretized with 400 cells. The adaptive time step was set to keep the Courant number less than 0.2. Based on estimates by Plourde et al. (2008), the selected mesh refinement is coarser than the required grid to fully resolve turbulent scales in a DNS (which would require about 720×360×360 cells). Nonetheless, this mesh is resolved enough to avoid the use of a subgrid-scale model. This can be verified by analyzing the energy spectra of fluctuations on the plume axis and at the plume outer edges. In Fig. 8 we show the energy spectra of temperature and pressure as a function of the non-dimensional frequency: the Strouhal number Str = f * 2b 0 /u 0 (f is the frequency in [Hz]). We recover a result similar to Plourde et al. (2008), where the inertial-convective regime with the decay −5/3 and the inertial-diffusive regime with the steeper decay −3 are observable (List, 1982).
Model results describe the establishment of the turbulent plume through the development of fluid-dynamic instabilities near the vent (puffing is clearly recognized as a toroidal vortex in Fig. 9a). The breaking of large-eddies progressively leads to the onset of the developed turbulence regime, which is responsible of the mixing with the surrounding ambient air, radial spreading of the plume and decrease of the plume average temperature and velocity. Figure 9a displays the spatial distribution of gas temperature. Mixing becomes to be effective above a distance of about 4 diameters. Figure 9b displays the distribution of the vorticity, represented by values of the Q u invariant (Eq. 58). The figure clearly identifies the toroidal vortex associated to the first instability mode (puffing, dominant at such Reynolds numbers). We have observed the other instability modes (helical and meandering, Lesieur et al., 2005) only by increasing the forcing intensity (not shown).
Experimental observations by George et al. (1977) and Shabbir and George (1994) reveal that the behavior of forced  plumes far enough from the inlet can be well described by integral one-dimensional plume models (Morton et al., 1956;Morton, 1959) provided that an adequate empirical entrainment coefficient is used. In the buoyant plume regime at this Reynolds number George et al. (1977) obtained an entrainment coefficient of 0.153. To compare numerical result with experimental observations and one-dimensional average plume models, we have time-averaged the numerical results between 4 and 10 s (when the turbulent regime was fully developed) and computed the vertical mass Q(z), momentum M (z) and buoyancy F (z) fluxes as a function of the height. To perform this operation, we define the time averaging operation (·) and the radial domain where (x,y,z) = x are the spatial coordinates, y tracer is the mass fraction field of a tracer injected from the vent with initial mass fraction y tracer,0 and u z is the axial component of the velocity field. We use this definition for Ω(z) for coherence with integral plume models, where the mean velocity field is assumed to have the same direction of the plume axis (cf. Morton et al., 1956;Woods, 1988;Cerminara, 2015a,b). This hypothesis is tested in Fig. 10a, where it can be verified that the time-averaged streamlines inside the plume are parallel to the axis (Fig. 10b shows the instantaneous streamlines and velocity magnitude field).
The plume fluxes are evaluated as follows (cf. George et al., 1977;Shabbir and George, 1994;Kaminski et al., 2005): where ρ α = ρ α (z) is the atmospheric density. From these quantities it is possible to retrieve the main plume parameters where (·) is the derivative along the plume axis and T α is the atmospheric temperature profile. Figure 11 displays the average plume radius and velocity. As previously reported by Fanneløp and Webber (2003) and Plourde et al. (2008), the plume radius initially shrinks due to the sudden increase of velocity due to buoyancy (at z = 0.1 m). Above, turbulent mixing becomes to be effective and increases the plume radius while decreasing the average velocity. The upper inset in Fig. 11 represents the values of the vertical mass q = Q/Q 0 , momentum m = M/M 0 and buoyancy f = F/F 0 , normalized with the inlet values. All variables have the expected trends and, in particular, the buoyancy flux is constant (as expected for weak ambient stratification) whereas q and m monotonically increase and attain the theoretical asymptotic trends shown also in Fig. 12. Indeed, Fanneløp and Webber (2003) have shown that an integral plume model for non-Boussinesq regimes (i.e., large density contrasts) in the approximation of weak ambient stratification and adopting the Ricou and Spalding (1961) formulation for the entrainment coefficient, has a first integral such that q 2 is proportional to m 5/2 at all elevations. Figure 12 demonstrates that this relationship is well reproduced by our numerical simulations, as also observed in DNS by Plourde et al. (2008).
The lower inset in Fig. 11 shows the computed entrainment coefficient, which is very close to the value found in experiments (George et al., 1977;Shabbir and George, 1994) and numerical simulations (Zhou et al., 2001) of an analogous forced plume. We found a value around 0.14 in the buoyant plume region (6.4 < z/2b 0 < 16).
The analysis of radial profiles led to a similar conclusions: in Fig. 13, we show the evolution of the radial profiles for the mean vertical velocity field. In this figure, we also report the plume radius as evaluated from Gaussian fits of these profiles on horizontal slices: The slope of the function b fit (z) has been evaluated in the region 6.4 < z/2b 0 < 16, to obtain b fit /z = 0.142 ± 0.001 to be compared with the result of George et al. (1977): b fit /z = 0.135 ± 0.010. Finally, figure 14 reports the time-average values of the vertical velocity and temperature along the plume axis. As observed in laboratory experiments, velocity is slightly increasing and temperature is almost constant up to above 4 inlet diameters, before the full development of the turbulence. When the turbulent regime is established, the decay of the velocity and temperature follows the trends predicted by the one-dimensional theory and observed in experiments. The insets displays the average value of the vertical velocity and temperature fluctuations along the axis. Coherently with experimental results (George et al., 1977), velocity fluctuations reach their maximum value and a stationary trend (corresponding to about the 30% of the mean value) at a lower height (about 3 inlet diameters) with respect to temperature fluctuations (which reach a stationary value about the 40% above 4 inlet diameters).

Transonic and supersonic flows
Although not essential in the present application, the ability of solving transonic and supersonic regimes is also required for the full-scale simulation of volcanic processes. We here test the behavior of the ASHEE code in presence of shocks in the classical Sod's shock tube test case (Sod, 1978) describing the expansion of a compressible, single-phase gas having adiabatic index γ = 1.4. At t = 0 the domain of length 10 m is subdivided in two symmetric subsets. In the first subset (spatial coordinate x < 0) we set u = 0, p = 10 5 Pa, T = 348.432 K, so that ρ = 1. In the second subset (x > 0), we set u = 0, p = 10 4 Pa, T = 278.746 K, so that ρ = 0.125 kg/m 3 . We indicate with c = 374.348 m/s the speed of sound of the gas in the x < 0 part of the domain. We impose zero gradient boundary conditions (∂ x (·) = 0) for all the variables u, p, T . As described in Sod (1978), a reference analytic solution exists for this problem.
In Fig. 15 we show the density profile obtained with the ASHEE model after 0.007 s of simulation. We performed two simulations at different resolution. The first has 100 cells and it is compared with the OpenFOAM solver rhoCentralFoam with a second order semi-discrete, non staggered central scheme of Kurganov et al. (2001) for the fluxes, and a total variation diminishing limiter (Van Leer, 1997) for the interpolation. We refer to Greenshields et al. (2010) for a presentation of rhoCentralFoam and of the Sod's shock tube test case. The inset of Fig. 15 is the simulation with an higher resolution (1000 cells). In this figure, we notice that the code performs satisfactorily both at low and high resolution. It is capable to capture the shocks pretty well, with a diffusion that is comparable with that obtained with rhoCentralFoam, a solver conceived for simulating shocks.

3D simulation of a turbulent volcanic plume
Numerical simulations of volcanic plumes were conducted in the framework of the IAVCEI (International Association of Volcanology and Geochemistry of the Earth Interior) plume model intercomparison initiative (Costa et al., 2015), consisting in performing a set of simulations using a standard set of input parameters so that independent results could be meaningfully compared and evaluated, discuss different approaches, and identify crucial issues of state of the art models. We here discuss three-dimensional numerical simulation of a weak volcanic plume in a stratified, calm atmosphere,  whose input data were set assuming parameters and meteorological conditions similar to those of the 26 January 2011 Shinmoe-dake eruption (Suzuki and Koyaguchi, 2013). Initial conditions and injection parameters are reported in Table  2.  Table 2: Vent conditions for the weak volcanic plume simulation.
The particle size distribution is composed of two individual classes of pyroclasts in equal weight proportion representing, respectively, fine (diameter d = 0.0625 mm; den-sity ρ = 2700 kg/m 3 , volume fraction = 0.00086821) and coarse ash (diameter d = 1.0000 mm; density ρ = 2200 kg/m 3 , volume fraction = 0.00106553). With respect to the laboratory benchmark case of Section 4.3, volcanic plumes are characterized by non-Boussinesq regimes at the vent and buoyancy reversal (with the initial mixture density about 4 times larger than the atmospheric one) and by a stratified atmosphere (Fig. 16). However, the most relevant difference is due to the significant temperature contrast (900 K) and to the presence of a high particle content which may strongly affect the mixing properties of the plume.
The Stokes number of the solid particles is, in general, a complex function of time and space, since the turbulent flow is characterized by a wide spectrum of relevant time and length scales. Generally, the Stokes number is associated with the most energetic turbulent eddy scale which, for laboratory plumes, has a typical turnover time of the order of τ L ∼ Str Dv Uv ≈ 0.12 s, where D v and U v are the plume diameter and velocity at the vent, respectively, and Str is the Strouhal number, of the order Str = 0.3 (Zhou et al.,  2001). Based on this time scale, and computing the particle relaxation time from Eq. 3, the Stokes number for the two adopted particle classes is about St coarse ≈ 5 and St fine ≈ 0.2, so we expect to see non-equilibrium phenomena for both particles classes, with more evident effects on the coarsest phase. However, the Stokes number, as an average value in all the plume is not as high as calculated above. Indeed, by using Eq. (64) as reference time for the turbulent dynamics, we obtain St coarse ≈ 0.6 and St fine ≈ 0.03. It is worth recalling here that the equilibrium-Eulerian approach is accurate and advantageous for particles having St ≤ 0.2 and that, in our model, we numerically limit the acceleration field in order to keep the turbulent non-equilibrium within this limit, as explained in Sect. 3 and tested in Sect. 4.2 Fig. 7. The averaged value of this limit -measuring the importance of the decoupling limiter for this simulation -is approximately 0.6.
The computational domain is cylindrical and is extended 483b 0 × 765b 0 in the radial and vertical directions (b 0 being the vent radius). The numerical grid is non-uniform and non-orthogonal. The discretization of the vent is represented in Fig. (17a). For the highest resolution run, the cell size increases from a minimum grid size ∆r = 2b 0 /32 with no radial grading factor in the region where the plume is expected to develop (Fig. 17b), whose initial radius is equal to 2.5b 0 and increases linearly with an angle θ such that tanθ = 0.147, slightly larger than tanθ = 0.12 predicted by the Morton's plume theory with entrainment k = 0.1 (Ishimine, 2006). Outside this region, a radial grading factor of 1.0446 is applied. Along z, 2048 cells are utilized. The minimum vertical cell size is ∆z = 2b 0 /32, and a grading factor of 1.00187 is imposed. The azimuthal resolution is constant and equal to 1 32 π (5.625 degrees). The resulting total number of cells is 10,747,904. This numerical mesh guarantees accuracy of the results: the solution procedure utilizes 2 PISO and 2 PIMPLE loops to achieve an absolute residual PIMPLE = 10 −7 (see Sec. 3).
Simulation of 720 s of eruption required about 490,000 time steps (imposing a CFL constrain of 0.2, resulting in an average time-step dt ≈ 1.5 ms, with a maximum velocity at the vent of about 150 m/s) for a total run-time of about 25 days on 1024 cores on the Fermi architecture at CINECA (meaning about 2.25 millions of cells per second, consistently with estimates of Sec. 4). Figure 18 shows the development of the volcanic plume at t = 400 s. Because of the atmospheric stratification, the plume reaches a neutral buoyancy condition at about 10 km above the vent (i.e., 11.5 km above the sea level, still within the troposphere). Due to its inertia, the plume reaches its maximum plume height H max ≈ 12 km, higher than the neutral buoyancy level, before spreading radially to form the socalled volcanic umbrella. The two orthogonal sections highlight the different spatial distribution of the volumetric fraction of fine (right) and coarse (left) ash particles, due to the different coupling regime with the gas phase. Coarse particles has indeed a larger settling velocity w s = τ s g which causes a more intense proximal fallout from the plume margins and a reduced transport by the umbrella. This is highlighted by the plot of the streamlines of the mixture velocity along a vertical section (Fig. 19), showing that the plume updraft is surrounded by a shell of settling coarse particles, which also inhibit air entrainment while promoting particle re-entrainment into the plume.
Besides settling, the large inertia of the coarse ash is responsible for the kinematic decoupling, leading to preferential concentration and clustering of particles at the margins of turbulent eddies. To illustrate this phenomenon, in a nonhomogeneous flow, the instantaneous preferential concentration is computed as the (normalized) ratio between the jth particle concentration and the concentration of a tracer (in   Table 2). Isosurface and vertical sections of the fine (light white) and coarse (light sand) ash volume fractions. The two-dimensional plots represent the distribution of the volume concentration of coarse (left) and fine (right) particles across vertical orthogonal slices crossing the plume axis.
our case, water vapor), i.e., C j = y j y j,0 · y tracer,0 y tracer , where the 0 subscript corresponds to the value at the vent. Fig. 20 shows the distribution of C j for the coarsest particles at t = 400 s. The color scale is logarithmic and symmetric with respect to 1, which corresponds to the nil pref- Fig. 19: Vertical section of the instantaneous value of the mixture velocity modulus (in logarithmic scale) at t = 400 s and velocity streamlines. erential concentration. For C j < 1, the mixture is relatively depleted of particles (green to blue scale); for C j > 1, particles are clustered (green to red scale), with mass fraction up to 5 times larger and 20 times smaller than the value it would have in absence of preferential concentration. This behavior is expected to affect the mixing and entrainment process. It is also worth remarking that the more uniform red area beyond the plume margins corresponds to the region of settling particles below the umbrella region. On the other hand, the top of the plume is relatively depleted of coarse particles. The corresponding Figure 21 for fine particles confirms that these are tightly coupled to the gas phase and almost behave as tracers (value of C fine is everywhere around 1). These conclusions are coherent with the a-priori estimate of St j we gave at the beginning of this section, based on the Taylor microscale time (Eq. (64)).
Finally, we present the results obtained by averaging the volcanic plume flow field over time (in a time-window [300-720] s where the plume has reached statistically stationary conditions) and over the azimuthal angle, in order to allow comparison with one-dimensional integral models (e.g., Woods, 1988) and discuss the effect of numerical resolution. The averaging procedure is the generalization of that explained in Sect. 4 to the multiphase case (see Cerminara, 2015a). The form of the results presented are similar to those obtained in Fig. 11 for the laboratory plume test case. Figure 22 presents the results of the averaging procedure for three multiphase flow simulations at different resolution (panels a-c). In particular, panel a) has the highest resolution (minimum radial cell size ∆r = 2b 0 /32 with 2b 0 equal to the inlet diameter); panel b) has ∆r = 2b 0 /16; panel c) has ∆r = 2b 0 /8. In panel d) we present results at the lowestresolution obtained by imposing the full kinematic equilibrium between gas and particles, i.e., by adopting the dustygas model (Marble, 1970).
Results demonstrate that the numerical model is quite robust and accurate so that even low-resolution simulations are able to capture the main features of the volcanic plume development. However, the maximum plume height systematically decreases from 12100 m (a), to 11300 m (b) to 11000 m (c) when we decrease the resolution. Analogously, the Neutral Buoyancy Level (NBL) decreases from 7800 m (a) to 7200 m (b) to 7100 m (c). Although the lowest resolution run seems to underestimate the maximum plume height and the plume radius by about 10%, the average velocity profile (and the vertical profiles of q, m and f ) is consistent in the three runs, showing a jet-plume transition at about 2000 m above the vent, also corresponding to the transition to a super-buoyancy region (Woods, 2010). The computed entrainment coefficient is also consistent and relatively independent on the grid resolution and shows a different behavior with respect to the laboratory case, associated with the effect of the density contrast. In this case, a maximum value of about k ∼ 0.1 is obtained in the buoyant plume region be-  Interestingly, we find that in three-dimensional simulations the entrainment decreases near the NBL and it become negative above that level. This happens because the mass exit from the plume region defined in Eq. (69) moving from it to the umbrella cloud. In this way, the mass flow q of the plume decreases above the NBL and a stationary solution can be achieved. This is not the case in integral plume models with positive entrainment coefficient, where the maximum plume height is reached as a singularity point with divergent mass flow and plume radius (cf. Morton, 1959;Woods, 1988). We plan to study this behavior more thoroughly in future studies.
The dusty-gas model shows a significantly different behavior, with a larger plume radius, a slightly higher entrainment coefficient and a more marked jet-plume transition with no further acceleration (without a super buoyancy transition). The plume height is slightly lower than the non-equilibrium case at the same resolution having maximum plume height and neutral buoyancy level of 9900 m and 6100 m, respectively. Numerical simulations thus suggest that the effects of non-equilibrium gas-particle processes (preferential concentration and settling) on air entrainment and mixing are non-negligible. These effects are certainly overlooked in the volcanological literature and will be studied more thoroughly in future studies, by applying the present model to other realistic volcanological case studies.

Conclusions
We have developed a new, equilibrium-Eulerian model to numerically simulate compressible turbulent gas-particle flows. The model is suited to simulate relatively dilute mixtures (particle volume concentration 10 −3 ) and particles with Stokes number St 0.2. It is appropriate to describe the dynamics of volcanic ash plumes, with kinematic decoupling between the gas and the particles, assumed in thermal equilibrium.
We have tested the model against controlled experiments to assess the reliability of the physical and numerical formulation and the adequacy of the model to simulate the main controlling phenomena in volcanic turbulent plumes, and in particular: 1) multiphase turbulence (including preferential concentration and density effects); 2) buoyancy and compressibility effects; 3) stratification and density nonhomogeneity.
The model reproduces the main features of volcanic plumes, namely: 1) buoyancy reversal and jet-plume transition; 2) plume maximum height and spreading of the umbrella above the neutral buoyancy level; 3) turbulent mixing and air entrainment; 4) clustering of particles; 5) proximal fallout and re-entrainment of particles in the plume. Results demonstrate that the compressible equilibrium-Eulerian approach adopted in the ASHEE model is suited to simulate the three-dimensional dynamics of volcanic plumes, being able to correctly reproduce the non-equilibrium behavior of gas-particle mixtures with a limited computational cost.
Finally, the adopted open-source computational infrastructure, based on OpenFOAM, will make the model easily portable and usable and will ease the maintenance and implementation of new modules, making ASHEE suitable for collaborative research in different volcanological contexts.