New developments in the analysis of column-collapse pyroclastic density currents through numerical simulations of multiphase flows

A granular multiphase model has been used to evaluate the action of differently sized particles on the dynamics of fountains and associated pyroclastic density currents. The model takes into account the overall disequilibrium conditions between a gas phase and several solid phases, each characterized by its own physical properties. The dynamics of the granular flows (fountains and pyroclastic density currents) has been simulated by adopting a Reynolds-averaged Navier-Stokes model for describing the turbulence effects. Numerical simulations have been carried out by using different values for the eruptive column temperature at the vent, solid particle frictional concentration, turbulent kinetic energy, and dissipation. The results obtained provide evidence of the multiphase nature of the model and describe several disequilibrium effects. The low concentration (≤5× 10−4) zones lie in the upper part of the granular flow, above the fountain, and above the tail and body of pyroclastic density current as thermal plumes. The high concentration zones, on the contrary, lie in the fountain and at the base of the current. Hence, pyroclastic density currents are assimilated to granular flows constituted by a low concentration suspension flowing above a high concentration basal layer (boundary layer), from the proximal regions to the distal ones. Interactions among the solid particles in the boundary layer of the granular flow are controlled by collisions between particles, whereas the dispersal of particles in the suspension is determined by the dragging of the gas phase. The simulations describe well the dynamics of a tractive boundary layer leading to the formation of stratified facies during Strombolian to Plinian eruptions.


Introduction
Pyroclastic density currents (PDCs) are among the most complex processes occurring during explosive volcanic eruptions (Branney and Kokelaar, 2002).They originate from eruptive columns formed by magma fragmentation processes, which arise in the volcanic conduit when tensile inner stress overcomes the magma breaking strength (Zhang, 1999).If the density of the eruptive column remains greater than the atmospheric one, the column collapses in a fountain, from which PDCs extend radially (Woods, 1995).The fountains are characterized by considerable unsteady interactions between the jet and the collapsing part of the columns, which produces recycling of eruptive material into the jet and oscillations in their heights (Valentine et al., 1991;Neri and Dobran, 1994).The complexity of the recycling processes is due to the multiphase nature of the eruptive columns, as well as to the transient and multidimensional properties of the fountains.
During the propagation of PDCs away from the vent, thermal plumes rise above the flow, while solid particles settle in its basal part (Druitt, 1998).The unsteady behaviour of the PDCs is caused by water vapour buoyancy effects occurring over them and by ash dragged toward the top by eruptioninduced winds (Valentine, 1998).Both processes are characterized by transient and multidimensional dynamics.
In the effort to reach a quantitative understanding of the dynamics of column-collapse PDCs, since the mid-1970s theoretical studies have been put together with geological analyses, as shown by Neri et al. (2003).Critical importance was attributed to processes such as fluidization and sedimentation of particles in the flow, as well as entrainment of air.
Published by Copernicus Publications on behalf of the European Geosciences Union.

S. Lepore and C. Scarpati: A new analysis of pyroclastic density currents
The development of numerical multiphase codes allowed the implementation of new procedures in the analysis of explosive volcanism.A two-dimensional and two-phase flow model, initially developed to numerically simulate a calderaforming eruption (Wohletz et al., 1984), was then improved to reproduce PDCs by Valentine and Wohletz (1989).Next, a two-component description of the gas phase and a kinetic depiction for the dense gas-particle regime (Dobran et al., 1993) and non-equilibrium effects between particles of two different sizes (Neri and Macedonio, 1996) were included.Thereafter, treating the gas and the solid phases as permeable continua, described by constitutive equations, multiphase flow models became particularly suited for describing transient and multidimensional non-equilibrium processes.The description of solid particle sedimentation, as well as of air entrainment and elutriation, was carried out directly from the set of constitutive equations without the need to define new parameters (Valentine, 1998;Macedonio and Neri, 2000;Burgisser and Bergantz, 2002;Valentine et al., 2002).At the same time, PDCs were treated as granular flows, defined as moving interstitial fluids in which an assemblage of discrete solid particles is dispersed (Dartevelle, 2004).These flows hold all the ordinary properties of multiphase flows and show a wide variety of behaviours and features.Depending on the loading conditions, the flows are highly dissipative because of frictions, inelastic collisions, and multiphase turbulence.Finally, they display a wide range of grain concentrations, as well as complex (non-linear, non-uniform, and unsteady) rheologies (Dartevelle, 2004, and references cited therein).To account for the whole spectrum of rheologies, the multiphase computer code GMFIX (Geophysical Multiphase Flow with Interphase Exchanges) has been devised, which can successfully simulate several pyroclastic phenomena and related eruptive processes (Dartevelle, 2004;Dartevelle et al., 2004).Multiphase flow models have been extensively tested through laboratory experiments and numerical simulations (Gidaspow, 1994;Boyle et al., 1998;Crowe et al, 1998;Dartevelle and Valentine, 2007).Meanwhile, a transient three-dimensional flow model of pyroclastic dispersion has been proposed, in which solid particles are considered to be in dynamic equilibrium with the gas phase (Oberhuber et al., 1998;Fadlun et al., 2000;Suzuki et al., 2005;Ongaro et al., 2007).
The aim of the present work is to improve the description of the PDC depositional processes and boundary layer, as inferred from multiphase numerical models reported in the literature (Dartevelle, 2004;Dartevelle et al., 2004).PDCs are considered as granular flows (Neri et al., 2003), and the GMFIX code is used to develop two-dimensional multiphase numerical simulations.Employing the RANS (Reynolds Average Navier-Stokes) model to describe the turbulence (Ferziger and Períc, 2002;Liu and Chow, 2002), numerical simulations are useful for bringing to light new information about the transport and deposition mechanisms of PDCs.The longitudinal transformations will be examined 1.00 × 10 −3 Heat capacity of solids, C s (J K −1 kg −1 ) (1.00 ÷ 1.30) × 10 3 Heat capacity of gas, C g (J K −1 kg −1 ) (3.30 ÷ 3.60) × 10 3 Maximum viscosity, µ max (kg m −1 s −1 ) 1.00 × 10 3 Viscosity of gas, µ g (kg m −1 s −1 ) 1.00 × 10 −5 Viscosity of solids, µ s (kg m −1 s −1 ) 5.00 1.00 ÷ 10.0 through the isolines of pyroclasts concentration, as a function of the distance and height from the vent, while the vertical transformations will be investigated through the concentration vs. the PDC height.Moreover, the study of the sedimentation rate in the basal part of the PDC as a function of time will be used to define the features of the aggrading deposit (Giordano et al., 2008).The results obtained will be compared with those of Dartevelle et al. (2004), to bring out the differences due to the use of the LES (Large Eddy Simulations) model (Ferziger and Períc, 2002;Liu and Chow, 2002).
Our simulations are performed using parameters to develop a Plinian-type eruption (Table 1) with the collapse of the eruptive column and the propagation of a PDC.The rheological conditions in the depositional system at the base of the PDC produce a tractive regime (see below).Finally, merging together the studies of the concentration trends, sedimentation rates as a function of time, and geological observations, we discuss how these numerical simulations can be useful to explain the presence of stratified facies in Strombolian to Plinian deposits.

Numerical technique
Granular flows are made up of a large number of particles that inelastically interact among each other.Therefore, being unfeasible to solve singularly the dynamics of each particle, the Implicit Multi-Field formalism (IMF), which handles all phases in the flow as permeable continua, has been employed.Each point variable (mass, velocity, temperature, pressure, etc.) is volume-averaged over a region greater than the particle dimension, but much smaller than the area of the whole flow domain (Syamlal et al., 1993).Thus, the detailed small-scale fluctuations within the flow are not analytically solved (they are somewhat smoothed out), and all the point variables are replaced by local average variables.
As for the averaged part, the physical description is made through the conservation equations of mass, momentum, and energy, which are formulated in terms of the local volumeaveraged variable for each phase.The complete list of symbols, units, and constants is provided in Table 1.
The Eqs. (1a), (1b), and (1c) state the following: the density change with time is equal to the momentum gradient; the momentum change in time and space equals the sum of drag force (friction between solid and gas), pressure gradient, viscous forces, and gravity force; the energy change is equivalent to the sum of the heat exchange between phases, the heat conduction of each phase, the work done by the drag force due to the frictional contacts within the flow, and the work associated with the volume change of the gas phase (Valentine and Wohletz, 1989;Dobran et al., 1993;Neri and Macedonio, 1996;Neri et al., 2003;Dartevelle, 2004).Unfortunately, during the averaging process, some information involving the bulk flow behaviour could be lost, and therefore the following constitutive equations are needed for interfacial drag D, viscous stress tensor τ , heat capacity C, exchange heat Q, and heat conduction q (Syamlal et al., 1993;Crowe et al., 1996): In Eq. ( 2), the subscript "t" denotes the mathematical operation of transpose, φ the angle of internal friction that describes frictional contacts between solid particles (Srivastava and Sundaresan, 2003;Kelfoun, 2011; Table 1), I the unit tensor, and N u the Nusselt number (Gunn, 1978), and β 0 , β 1 , and β 2 are numerical constants.
In granular flows, a viscous dissipation within the solid phase is also present.It is related to the particles volumetric concentration and is described through variations of the granular energy E γ in time and space, as follows: In Eq. ( 3), is the dissipation function, the density variation in time, and the viscous dissipation involved in slipping, collisions, and dragging.From experimental, numerical, and theoretical studies, it was inferred that three granular behaviours can be discerned (Dartevelle, 2004).In the low concentration part of the flow, particles fluctuate and translate randomly, thus producing a viscous dissipation called kinetic.At higher volumetric concentration (1 % < α s < 50 %, Table 1), the collisions among particles become predominant; hence, the viscous dissipation becomes collisional-like.The collisions among particles and between particles and walls are characterized by dimensionless restitution coefficients.At very high volumetric concentration (50 % < α s < 65 %, Table 1), particles endure sliding and rubbing contacts that are the source of frictional dissipation.The frictional contacts between particles and walls are described through the angle of wall friction θ (Srivastava and Sundaresan, 2003; Kelfoun, 2011; Table 1).
Granular flows are also turbulent flows showing time and space dependent properties.With the aim of analysing accurately the turbulence features, the computational domain, 10 km long (L x ) and 2.5 km high (L y ), is divided into 1000 cells, each of 10 m (R x )× 2.5 m (R y ) (boundary conditions in Table 1).Moreover, it is necessary to add transient terms to the conservation equations, which become neither non-linear nor analytically resolvable.Thus, a numerical iterative solution procedure must be adopted to describe the fluctuations of each quantity in time and space due to the turbulence.To facilitate the procedure, under-relaxation parameters in the iterations could be used (Ferziger and Períc, 2002).To describe the effects on the motion of solid particles within the granular flow, two fundamental models are reported in the literature: namely the LES and the RANS (Moeng, 1984;Smagorinsky, 1993;Leith, 1993;Dartevelle et al., 2004).
The LES model depicts the large Kolmogorov scale effects of turbulence fields, composed of eddies having diameter from 10 m to ∼1 km.The model represents separately the eddy dynamics within the subgrid scale, which goes from cm to a few metres.Moreover, it assumes that the turbulent kinetic energy cascades from the largest to the smallest eddies, until the granular viscosity is able to dissipate the transferred kinetic energy (Moeng, 1984).The energy spectrum of the www.solid-earth.net/3/161/2012/Solid Earth, 3, 161-173, 2012 turbulence is given by the following equation: In Eq. ( 4), ξ is the energy cascade rate and k the conductivity.From the energy spectrum, simple scaling laws allow deducing appropriately the eddy viscosity and the eddy thermal conductivity.These two quantities are used to define the turbulent subgrid shear stress and the turbulent heat flux (Smagorinsky, 1993;Leith, 1993).The turbulent viscosity is related to the Smagorinsky constant ℘ as follows: In Eq. ( 5), r is the geometric average of the grid size and ||ψ|| the Euclidian norm of the rate-of-strain tensor (Nieuwstadt et al., 1991).Several authors have used the LES model in the last 20 years (Dartevelle, 2004, and references cited therein).
The RANS model of turbulence averages out the whole set of unsteadiness, considering all the instabilities as part of the turbulence.The averaging process is reasonable at high Reynolds numbers: in that case, in fact, there is a cascade of energy from the largest to the smallest scales, where viscous effects prevail and the energy is dissipated (McWilliams, 1990).Therefore, every quantity can be written as the sum of a time-averaged value and a fluctuation around that value: In Eq. ( 6), the averaging interval must be greater than the typical time scale of the fluctuations, and the effects of turbulence are represented as an increased viscosity.The application of this factorization in the Navier-Stokes momentum equation leads to the following one, where the turbulent kinetic energy K tur (Table 1) has been introduced: In Eq. ( 7), is the Reynolds constant and ε the dissipation (Table 1).The turbulent kinetic energy and the dissipation are connected through the following equation (Ferziger and Períc, 2002): In Eq. ( 8), L is a mixing length scale (Hoffman and Chiang, 2000) introduced to relate the turbulent viscosity with the mean velocity of the flow (Odier et al., 2009).Unlike the LES, the RANS model has not often been used in numerical simulations of geophysical multiphase granular flows.However, from a mathematical point of view, the RANS describes more satisfactorily the interactions between particles and walls (Benhamadouche and Laurence, 2003) and represents more simply the flow transformations due to the turbulence (Cantero et al., 2008).For these reasons, it has been preferred to LES to analyze PDCs dynamics.
All the above described equations are solved numerically by means of the recently validated GMFIX software (Dartevelle and Valentine, 2007).That requires the definition of a computational domain, the boundary conditions and the initial ones.The domain used in numerical simulations is reported in Fig. 1, where cylindrical coordinates are used: the symmetry axis corresponds to the main axis of the eruptive pyroclastic column, and the boundary conditions are defined as reflecting at left, outflow at right and top, and no-slip at bottom.
The GMFIX software derives from the computer code MFIX (Multiphase Flow with Interphase eXchanges).MFIX is a FORTRAN 90 general purpose code (Boyle et al., 1998;Syamlal, 1998;Guenther and Syamlal, 2001), assembled in the Los Alamos National Laboratory from the K-FIX (Kachina with Fully Implicit eXchange) code (Dartevelle et al., 2004).MFIX has been updated into the Geophysical version GMFIX to deal with typical geophysical applications, maintaining all the previous capabilities and adding new ones, such as the work associated with volumetric variations of the gas phase, the standard atmospheric profiles, the LES and RANS turbulence models, and the subgrid turbulent heat flux (Dartevelle, 2004;Dartevelle et al., 2004).The "FIX" family codes have been used successfully in volcanological analyses, starting from Valentine and Wohletz (1989).
The IMF formalism adopted by the "FIX" codes allows handling the main features of multiphase flows, as it manages the discretization of the conservation equations in space Solid Earth, 3, 161-173, 2012 and time (Lakehal, 2002).For the space discretization, GM-FIX uses a finite volume method, which entails that the physical domain is divided into discrete three-dimensional cells, over which the conservation equations are integrated.This integration procedure ensures conservation of mass, momentum, and energy, on the entire domain.Scalar quantities, such as mass and temperature, are computed at the cell centre, whereas velocity components are calculated along the cell boundaries.The cell dimension is critical, as solid particles settling and deposition could be neglected by an excessively (larger than 100 m 2 ) coarse grid (Patankar, 1980;Dobran et al., 1993;Neri et al., 2003;Dartevelle et al., 2004).As for the time discretization, GMFIX uses an implicit backward Euler method, and includes various first-order (e.g.FOU) and second-order (e.g.Superbee, Smart, and Minmod) accurate schemes for discretizing the convection terms (Syamlal, 1998).The FOU discretization scheme was favoured for its stability, better convergence, and because significant differences were not seen with respect to the secondorder schemes (Dartevelle et al., 2004).The products obtained through numerical simulations performed with GM-FIX have been processed by using the software MATLAB to generate isolines contour plots.

Results
To perform the numerical simulations presented here, the turbulence is described through the RANS model; the boundary and the initial conditions are detailed in Table 1  Temperature, heat capacity, frictional concentration of solid particles (defined as the volume of solids divided by the total volume of the flow), turbulent kinetic energy and dissipation are parameters that vary between two extreme values (Table 1), while the others are assumed to be constant.The values for velocity and vent dimension are determined from field evaluations of mass flux (Dartevelle et al., 2004).The density of solid particles is the average between the values for lithics and pumices, while the frictional and the maximum concentrations of solid particles are assigned according to Srivastava and Sundaresan (2003).The values for turbulent kinetic energy and dissipation are defined in agreement with Ferziger and Períc (2002).A detailed description of the results obtained from one simulation is reported as an exam-ple.The results from the analysis of the other simulations will be synthesized at the end of this section.
A plot of the solid particle concentration (reported in decimal logarithmic units) as a function of height and distance from the vent at different times is shown in Fig. 2. At 20 s, a pyroclastic fountain begins to collapse and a PDC generated from it expands laterally to 900 m.Above the fountain, a low concentration (≤5 × 10 −4 ) thermal plume ascends in the atmosphere up to 500 m.At 40 s, the PDC propagates to 1900 m from the vent.Its high concentration basal part, reported as red-coloured, propagates to 1800 m, with an average velocity of 45 m s −1 .The tail of the PDC, which reaches a maximum height of 200 m, is the thickest part.The body and the head, instead, are 120 m high.The plume above the fountain ascends in the atmosphere up to The maximum height of the tail diminishes to 110 m.The height of the body and head reduces to 80 m as well.The thermal plumes above the tail and body increase continuously in number and rise in the atmosphere up to 300 m.The plume above the fountain continues to ascend in the atmosphere beyond the height of 2500 m.At 180 s, the PDC propagates to 8 km from the vent while its high concentration basal part propagates to 6.5 km with an average velocity of 36 m s −1 .Its tail, which reaches a maximum height of 100 m, is now the thinnest part.The thermal plumes present above it reach a maximum height of 450 m.The PDC body, instead, has become 150 m high, and the thermal plumes above the tail and body attain a maximum height of 850 m.The PDC head, 310 m high, is now the thickest part.This abrupt increase in the height is caused by a large accretion of air ingestion in the last 20 seconds.The plume above the fountain ascends in the atmosphere over the height of 2500 m.
The changes in the features of fountain and PDC dynamics are well described by the plot in Fig. 3 of the sedimentation rate (SR) vs. time.SR is defined as the product of the average horizontal velocity of solid particles and their concentration within the basal part, at a certain time and distance from the vent.The chosen distance is 500 m, since, at that distance, the PDC body is present in all the snapshots of Fig. 2.During the first 40 s the SR increases rapidly, showing an unsteady behaviour.Then this parameter fluctuates only slightly, denoting a quasi-steady behaviour between 40 s and 160 s.Unsteadiness appears again during the final 20 s of this simulation.Waxing unsteadiness is related to a strong increase of pyroclastic concentration possibly due to vertical segregation of clasts during the initial phase of an overloaded current (Hiscott, 1994).
The trends of solid particle concentration as a function of the PDC height, evaluated at the distance of 500 m from the vent and at different times, are analyzed in Fig. 4a.The concentration at the base (about 12 m high) of the PDC (Fig. 4b) increases quickly between 20 s and 40 s, then slowly until 180 s, reaching a maximum value at about 3 %.In the basal part of the PDC, therefore, solid particles settle continuously reaching a boundary layer (BL) with a low particle concentration.Above the base, the concentration of the solid particles decreases progressively passing from 20 s to 80 s. Between 100 s and 180 s, a sudden jump in the concentration, characterized by some oscillations due to the formation of the thermal plumes, is present in the overall gradual decreasing trend.
Finally, the analysis of the results obtained from the other simulations shows that the BL concentration increases from 3 % as maximum value, to 4 % as the temperature decreases from 1200 to 900 K. Besides, the BL concentration increases from 3 to 3.5 % as the solid particle frictional concentration grows from 0.55 up to 0.65, while the BL concentration remains almost unaltered around 3 % as the turbulent kinetic energy increases from 0.01 to 0.2 m 2 s −2 and the dissipation from 1.0 to 10.0 m 2 s −3 .

Discussion
The initial and boundary conditions reported in Table 1, used in all our numerical simulations, pertain to the activity of collapsing eruptive columns and associated PDCs occurring during Plinian eruptions, as stated by the value of the mass discharge rate, which is 2.3 × 10 7 kg s −1 .The set of the results obtained from all the simulations shows a vertical stratification of the solid particle concentration and a quick development of a thin BL under the more dilute and turbulent PDC portion.In addition, the SR shows low values and the concentration in the BL increases at increasing time, first quickly, and then slowly, up to a maximum value in the 3/4 % range.
According to the literature, the maximum concentration value reached implies that the BL dynamics is tractive (Lowe, 1982).The existence of a tractive BL within a turbulent PDC suggests that the dynamics in its basal  part can be assimilated to those producing stratified facies.They can be crossed or diffused according to the values assumed by the volumetric concentration of solids in the BL (Chough and Sohn, 1990;Branney and Kokelaar, 2002).The formation of a stratified deposit is suggested by the low maximum value of the sedimentation rate (Giordano et al., 2008).A stratified facies is present when the concentration in the BL is sufficiently high to prevent turbulent sorting of fine ash (Cole and Scarpati, 1993).Turbulent eddies are unable to penetrate into the BL; however, they exert fluctuating shear stresses on it.The deposition (stepwise aggradation) will be unsteady and produce thin, nonpersistent and variably stratified deposits (Branney and Kokelaar, 2002).The formation of stratified facies is not related to the magnitude of the eruption, but only to the development of appropriate depositional conditions within the PDCs.Stratified facies have been observed inside the PDC successions of Strombolian to Plinian eruptions.Here, we briefly report some examples in the Campanian volcanic area.The different areas covered by selected pyroclastic deposits are reported in Fig. 5. Stratified facies are well developed in Strombolian deposits covering an area of less than 1 km 2 as well as Plinian deposits up to 500 km 2 .Trentaremi and Capodimonte tuff cones (Cole et al., 1994) are small monogenetic vents of phreatomagmatic style.Their deposits show a limited dispersion.The flanks of both cones are formed by stratified tuff with several sand-wave structures (Fig. 6a and b), which derive from turbulent low-concentration PDCs.Different stratified facies (Fig. 6c) have been recognized in the deposits formed by PDCs during Neapolitan Yellow Tuff, which was a large and phreatoplinian eruption.These facies have been interpreted as the product of deposition from turbulent low-concentration PDCs, due to traction sedimentation (Cole and Scarpati, 1993).Finally, a stratified facies has been recorded in the succession (Fig. 6d) that occurred during the 79 AD Plinian Vesuvius eruption (Luongo et al., 2003).
As previously shown, convective thermal plumes develop above the PDC tail and body.These low concentration plumes (less than 10 −4 ) bring small-dimensional solid particles up in the atmosphere (to 1 km height above the sea level).A fall deposition of pyroclastic fragments occurs from them (Sparks and Walker, 1977), causing the formation of vesiculated ash layers (Cole and Scarpati, 1993) at the end of the eruptive phase or when a gap in the emplacement of PDC occurs.The formation of thermal plumes above the body and tail does not affect the BL behaviour.Namely, a net dissociation exists between the PDC boundary layer, where a kinetic-collisional regime develops, and the low concentration suspension that flows above, in which a pure kinetic regime holds.Indeed, the model employed and applied to the study of the PDCs dynamics shows that the high concentration boundary layer never outruns the upper more diluted suspension.This has been noticed in various eruptions, such as the Neapolitan Yellow Tuff (Cole and Scarpati, 1993) and the 79 AD Vesuvius Plinian (Luongo et al., 2003).
The simulation features described here are rather different from those obtained by Dartevelle et al. (2004).In our case, results show the following: the head thickness exceeds the rest of the PDC only at 180 s; at 100 s BL never outruns the head, and its concentration is ∼2.75 %; at 180 s, the BL concentration we get is ∼3/4 % and then it displays kinetic-collisional behaviour.In simulations of Dartevelle et al. (2004), conversely the following was found: the head of the PDC is very high from the beginning; at 100 s, the high concentration BL has outrun the head of the PDC and has a concentration of ∼30 %; at 180 s, the BL concentration is ∼55 % and so it shows plastic-frictional behaviour.Presumably, these differences are due to the different models used to describe the turbulence.In the RANS model, the main interest is to capture the bulk properties of the PDC, namely the macroscopic features (large-scale) of the granular flow.The turbulent fluctuations (small-scale) are interpreted as deviations from the averaged values.Instead, in the LES model the large-scale properties of the granular flow are described independently from the small-scale features.The cut-off between the large scales and the subgrid scales is supposed to take place in the stationary subgrid range (Dartevelle, 2005).The main difference between the two models, therefore, is the way in which subgrid and macroscopic scales are correlated.This difference has strong consequences on the overall simulation, and different results are obtained.From the analysis of the solid particle concentration within the BL at 180 s in the two cases, it follows that the associated deposit will show stratified facies in our simulations and massive facies in the Dartevelle et al. (2004) simulation.Consequently, the simulations performed by using the RANS model can be best applied to the study of the dynamics of eruptive phases in which associated deposits display stratified facies.On the contrary, the simulations performed by using the LES model can be best applied to the analysis of the dynamics of eruptive phases where associated deposits show massive facies.

Conclusions
Two-dimensional numerical simulations of the dynamics of fountains and associated PDCs have been performed by using a granular multiphase model, with the aim of reproducing the dynamics of stratified lithofacies recurrent in pyroclastic density current deposits of both historical and prehistorical eruptions.Pyroclastic density currents are considered as granular flows, and the GMFIX code, together with the RANS model for describing turbulence, is used to carry out numerical simulations.The analysis of the results allowed not only the description of the granular flow macroscopic dynamics, but also the depiction of sedimentation into the boundary layer and deposition from it.On a large scale, the rising mixture above the collapsing fountain, the upper part of the granular flow and thermal plumes above the pyroclastic current tail and body are characterized by a low concentration (≤ 5 × 10 −4 ) of solid particles.The collapsing fountain and the basal part of the density stratified pyroclastic current, instead, show a higher concentration of solid particles.From our simulations, it appears that granular flows are therefore formed by a high concentration boundary layer underlying a low concentration suspension.This structure of the flow is present from the proximal to the distal regions.On a small scale, the dynamics in the boundary layer is strongly affected by interactions between solid particles, whereas in the overlying low concentration suspension, the dispersion of the solid particles is controlled by the dragging of the gas phase.Significant thermal disequilibrium effects are observed between gas and solid particles in the granular flow especially in the regions of effective air entrainment, as in the head of the flow.The analysis of the maximal values of solid particles volumetric concentration reached in the boundary layer implies that its dynamics is tractive and suggests a connection with the mechanisms leading to the formation of stratified facies.The results from our numerical simulations reproduce some of the features that lead to the stratified facies observed in the field and associated with a wide range of physical conditions (e.g.magmatic and phreatomagmatic styles of Strombolian to Plinian events).A deeper knowledge of the processes investigated should be attained by the development of more detailed granular multiphase models, combining field-based investigation on deposits of past eruptions with modern computational fluid dynamics techniques.

Fig. 2a .
Fig. 2a.Snapshots of the solid particle concentration (logarithmic scale) as a function of distance and height at 20 s, 40 s, and 60 s.

Fig. 2b .
Fig. 2b.Snapshots of the solid particle concentration (logarithmic scale) as a function of distance and height at 80 s, 100 s, and 120 s.

Fig. 2c .
Fig. 2c.Snapshots of the solid particle concentration (logarithmic scale) as a function of distance and height at 140 s, 160 s, and 180 s.

Fig. 3 .
Fig. 3. Sedimentation rate in the PDC as a function of time at 500 m from the vent.

Fig. 4a .
Fig. 4a.Temporal curves of solid particles logarithmic concentration as a function of height at 500 m from the vent.

Fig. 4b .
Fig. 4b.Zoom of Fig. 4 illustrating the solid particle concentration at the base of the PDC.

Fig. 5 .
Fig. 5. Areas covered by some of the stratified pyroclastic current deposits of Phlegraean Fields and Vesuvius.

Fig. 6 .
Fig. 6.Stratified deposits: (A) Trentaremi Tuff showing abundant cross-stratification and sand-wave structures; (B) panoramic view of the Capodimonte Tuff, with undulating thin ash and fine lapilli layers; (C) section through the Neapolitan Yellow Tuff.Note the complex multilayered nature of this deposit; (D) sedimentary structures in the AD 79 deposit in the ancient town of Pompeii.Fragments of tiles and walls are present in the deposit.

Table 1 .
Boundary and initial conditions used for numerical simulations.