New developments in the analysis of volcanic pyroclastic density currents through numerical simulations of multiphase flows

Introduction Conclusions References


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 the magma fragmentation process, which occurs in the volcanic conduit during magma ascent (Sparks, 1976).If the density of the column remains greater than the atmospheric one, the column collapses in a fountain from which radially spreading PDCs form (Woods, 1995).The fountains are characterized by significant unsteady interactions between the jet leaving the vent and the collapsing column, which produces recycling of eruptive material into the jet and oscillations in their height (Valentine et al., 1991;Neri and Dobran, 1994).The complexity of the processes is due to the multiphase and multi-component nature of the eruptive columns, as well as to the transient and multidimensional features of the fountains.
During the PDCs propagation away from the vent, coignimbritic plumes form above the flow, while particle sedimentation and sorting occur inside (Sparks, 1976;Druitt, 1998).Their unsteady behaviour is caused by ash dragged inward the top by eruption induced winds, and by water vapour buoyancy effects occurring over it (Valentine, 1998).Both processes are characterized by transient and multidimensional dynamics.
Since the mid-1970s, theoretical works have been combined with geological studies in the attempt to describe the PDCs dynamics on a more physical basis, as shown by Neri et al. (2003).Important processes such as fluidization and sedimentation of particles in the flow, as well as entrainment of air, were recognized as crucial.
The development of numerical multiphase codes helped in attempting the analysis of explosive volcanism.A two-dimensional, transient, two-phase flow model, initially developed to numerically simulate a caldera-forming 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 interpenetrating continua with specific constitutive equations, multiphase flow models became particularly suitable for describing strongly transient and multidimensional non-equilibrium processes.The description of the particles segregation, air entrainment, and air elutriation, came out from the set of constitutive equations, instead of defining new parameters (Valentine, 1998;Macedonio and Neri, 2000;Burgisser and Bergantz, 2002;Valentine et al., 2002).Contemporaneously, PDCs were treated as granular flows, which are defined as a collection of discrete solid particles dispersed in a moving interstitial fluid.These flows own all the common 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 (Syamlal, 1987;Besnard and Harlow, 1998;Brey et al., 1999;Kashiwa and Heyden, 2000;Lakehal, 2002).To account for the whole spectrum of granular rheologies, the multiphase computer code GMFIX (Geophysical Multiphase Flow with Interphase Exchanges) has been employed, which can successfully simulate a large span of 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 on a large scale was developed, in which particles were considered 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 PDCs depositional processes and boundary layer reported in the literature (Dartevelle, 2004;Dartevelle et al., 2004).PDCs are considered as granular flows 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 possible new features in the PDCs.The longitudinal transformations will be examined through the plots of the pyroclasts concentration isolines as a function of the distance and height from the vent, while the vertical transformations through the plots of the concentration vs the PDC height.Moreover, by the study of the sedimentation rate as a function of time, the features of the basal part of the PDC will be used to characterise the deposit originating under the PDC (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).Finally, merging together the studies of the concentration trends, sedimentation rates as a function of time, and geological observations, the numerical simulations will be related with some phases of strombolian to plinian eruptions occurred in the Phlegraean and Vesuvius areas (South Italy) as the Neapolitan Yellow Tuff (Cole and Scarpati, 1993), the Capodimonte Tuff and the Trentaremi Tuff (Cole et al., 1994), and the 79 AD Vesuvius Plinian (Luongo et al., 2003).

Numerical Technique
Granular flows are made up of a large number of particles that inelastically interact among them.Therefore, being impractical to solve singularly the motion of each individual particle, the Implicit Multifield Formalism (IMF), which treats all phases in the flow as interpenetrating continua, has been employed.Each point variable (mass, velocity, temperature, pressure, etc.) is volume-averaged over a region greater than the particles 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 volume-averaged variable for each phase.The three equations state that: 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 to 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, in the averaging process, some information that may involve the bulk flow behaviour are lost, and therefore it is necessary to supply the following constitutive equations for interfacial drag, viscous stress tensor, heat capacity, exchange heat, and heat conduction (Syamlal et al., 1993;Crowe et al., 1996).
( ) In Eqs.(1), s denotes the solid phase, g the gas phase, α the concentration, ρ the density, v the velocity, d s the diameter of the solid particles, τ the viscous stress tensor, σ the stress, φ the angle of internal friction (Srivastava and Sundaresan, 2003), C v the heat capacity, β a numerical coefficient, T the temperature, Q the exchange heat between phases, k the thermal conductivity, N u the Nusselt number (Gunn, 1978), µ the viscosity, and q the heat conduction.
In the 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. ( 2), Φ is the dissipation function, P the pressure, ℑ the density variation in time, and Γ the viscous dissipation involved in slipping, collisions, and dragging.It has been recognized experimentally, numerically, and theoretically, that three granular behaviours can be discerned, as shown by Dartevelle (2004).In the dilute part of the flow, particles fluctuate and translate (1) (2) randomly: this motion produces a viscous dissipation called kinetic.At higher volumetric concentration (1 % < α s < 50 %), the collisions among particles become predominant, and then the collisional dissipation takes the place of the kinetic.The collisions among particles and between particles and walls are characterized by dimensionless restitution coefficients.At very high volumetric concentration (> 50 %), 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).Granular flows are also turbulent flows.Turbulence is a very complex phenomenon and, in general, modelling it in details is difficult without some empiricism.It strongly depends on mean flow properties, and so it is time and space dependent.To analyze in depth the features, a division in cells of the computational domain (Fig. 1) is needed, which entails the creation of a grid.Turbulence makes necessary to add transient terms to the conservation equations, which become non-linear nor analytically resolvable.Thus, numerical solution procedures must be adopted, in which the advance in time and space is acquired by means of iterations on each quantity.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 simulates the large Kolmogorov scale effects of turbulence fields, composed by eddies having diameter from 10 m to ~1 km.The model describes eddies dynamics within the subgrid scale, which goes from cm to a few meters.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 turbulence is given by the following equation: 3), ξ is the energy cascade rate and k the conductivity.From the energy spectrum, simple scaling laws can be employed to deduce, in an appropriate way, 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. ( 4), r is the geometric mean 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 twenty years, as shown by Dartevelle (2004).
The RANS model of turbulence averages out the whole set of unsteadiness, considering all the instabilities as part of the turbulence.That is reasonable, since at high Reynolds numbers there is a cascade of energy from the largest to the smallest scales, with a dissipation of the transferred energy.In this model, every quantity Ω is written as the sum of a time averaged value and a fluctuation about that value: In Eq. 5, the averaging interval must be greater than the typical time scale of the fluctuations and the effects of turbulence are assumed to be represented as an increased viscosity.The application of this factorization in the Navier-Stokes momentum equation leads to the following new one, where the turbulent kinetic energy has been introduced: ( (5) In Eq. ( 6), ℜ is the Reynolds constant and ε the dissipation.The turbulent kinetic energy and the dissipation are connected through the following equation (Ferziger and Períc, 2002): In Eq. ( 7), 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 been used often in numerical simulations: that is the reason why this model has been chosen here to analyze PDCs dynamics.
All the above described equations are solved numerically by means of the recently validated GMFIX (Geophysical Multiphase Flow with Interphase eXchanges) software (Dartevelle and Valentine, 2007).To that, it is necessary to define 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 an extensively validated (Boyle et al., 1998) FORTRAN 90 general purpose code (Syamlal, 1998), assembled in the Los Alamos National Laboratory from the K-FIX (Kachina with Fully Implicit eXchange) code (Rivard and Torrey, 1977).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: it manages, in fact, the discretization of the conservation equations in space and time (Rivard and Torrey, 1977;Lakehal, 2002).For the space discretization, GMFIX uses a finite volume method, where the physical domain is divided into discrete three-dimensional volumes, over which the conservation equations are integrated.This integration procedure ensures global conservation of mass, momentum, and energy, independently of the grid size.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, because an excessively coarse grid may neglect particle settling and deposit building (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 as respect to the second-order schemes (Dartevelle et al., 2004).The products obtained through numerical simulations performed with GMFIX have been processed by using the software MATLAB, to generate isolines contour plots (for example, of concentration).
(6) (7) 3 Results To perform the numerical simulations here presented, the turbulence is described through the RANS model: the boundary and the initial conditions are detailed in Table 1.
Temperature, solid particles frictional concentration, turbulent kinetic energy and dissipation are parameters which vary between the two extreme values reported in Table 1.A detailed description of the results obtained from one simulation is reported as an example.The results from the analysis of the other simulations will be synthesized at the end of this paragraph.
A plot of the solid particles concentration as a function of height and distance from the vent at different time is shown in Fig. 2. At 20 s, the pyroclastic fountain begins to collapse and expands to 900 m.The flow sectors at low concentration (≤ 5 × 10 -4 ) ascend in the atmosphere up to 500 m.At 40 s, the PDC generated from the pyroclastic fountain propagates to 1900 m from the vent.Its basal part, reported as red-coloured, propagates to 1800 m, with an average velocity of 45 m/s.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 flow sectors with a low concentration ascend in the atmosphere up to 950 m.At increasing time, from 60 s to 100 s, the PDC continues to propagate to 4 km from the vent, and its basal part to 3900 m.The average velocity decreases, reaching 43 m/s.The maximum height of the tail decreases at increasing time to 160 m.Similarly, the height of the body and head reduces to 90 m.The flow sectors at low concentration continue to ascend in the atmosphere up to 2500 m.At 100 s, the PDC propagates to 4100 m from the vent.The basal part propagates to 4 km, with an average velocity of 40 m/s.The tail, which reaches a maximum height of 150 m, is still the thickest part.The body and the head have become 85 m high.Several thermal plumes are forming along the tail and the body, produced by the loss of the momentum within the PDC and by the dilution due to air ingestion coming from the head.The flow sectors at low concentration ascend in the atmosphere to a height little more than 2500 m.From 120 s to 160 s, the PDC continues to propagate to 6 km from the vent, and its basal part to 5900 m.The average velocity keeps decreasing, reaching 37 m/s.The maximum height of the tail diminishes to 160 m.The height of the body and head reduces to 80 m as well.The thermal plumes above the tail and body keep increasing in number and prolong their ascending in the atmosphere up to 300 m.The flow sectors at low concentration continue to ascend in the atmosphere, well beyond the height of 2500 m.At 180 s, the PDC propagates to 8 km from the vent, while its basal part propagates to 6.5 km with an average velocity of 36 m/s.Its tail, which reaches a maximum height of 100 m, is now the thinnest part.The thermal plumes present above it, produced by the loss of momentum, 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, high 310 m at maximum, is now the thickest part.This abrupt increase in the height is caused by a large accretion of air ingestion in the last twenty seconds.The flow sectors at low concentration ascend in the atmosphere abundantly over the height of 2500 m.
The changes in the behaviour of PDC dynamics are well described by the plot in Fig. 3 of the sedimentation rate (SR) vs time.The SR is defined as the product of the average velocity of solid particles and their concentration within the basal part, at a certain time and distance from the vent.The chosen value is 500 m, since, at that distance, the PDC body is present in all the snapshots of Fig. 2. Between 0 s and 20 s, the SR shows a steep slope caused by a continuous increase of the pyroclastic fountain collapse.Then, between 20 s and 40 s, the SR decreases because the generation of the PDC from the pyroclastic fountain reduces its collapsing rate.Up to 80 s, the SR keeps decreasing because the behaviour of the PDC has not changed markedly.Between 80 s and 160 s, the SR displays a plateau, meaning that the formation of thermal plumes does not affect the settlement of the solid particles in the basal part of the PDC.At last, between 160 s and 180 s, the SR again increases steeply because of a quick intensification in the dilution rate of the upper part of the PDC and in the supply rate of solid particles into its basal part.
The trends of solid particles 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. 4. The concentration at the base (about 15 m high) of the PDC increases quickly between 20 s and 40 s, then slowly until 180 s, reaching a maximum value at about 3 %.This part of the PDC, therefore, can be considered as a boundary layer (BL) with low solid particles 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, instead, 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 K to 900 K. Besides, the BL concentration increases from 3 % to 3.5 % as the solid particles 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 numerical simulations described show a vertical stratification of the solid particles concentration and a quick development of a thin BL under the more dilute and turbulent PDC portion.In addition, the SR does not exceed 2.0 m/s 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 traction BL within a turbulent PDC suggests a connection with the mechanisms leading to the formation of stratified facies.They can be cross or diffuse in relation 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).Although turbulent eddies are unable to penetrate into the BL, they exert fluctuating shear stresses on the BL.The deposition will be unsteady (stepwise aggradation) and produce thin, nonpersistent and variably stratified diffuse bedding.On the contrary, when the turbulent eddies and the fluctuating shear stresses they produce are not transmitted down to the BL at all, their action is dampened by a relatively high-concentration fluid layer dominated by gas escape.Deposit stepwise aggradation is thus fairly steady despite overriding perturbations, and massive lapilli tuff is formed (Branney and Kokelaar, 2002).
The formation of stratified facies is never correlated with the magnitude of the eruption, but only with the development of appropriate depositional conditions within the PDCs.Stratified facies have been observed inside the PDCs succession from small to very large eruptions 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 (Figs.6a and 6b), 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, a large, 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 deposition from PDCs has been recorded in the succession (Fig. 6d) occurred during the 79 AD Plinian Vesuvius eruption.From the study of the distribution of the damages provoked by the impact of the PDCs in Pompeii it comes out that the intensity of the destruction is elevated when an open space is present between two next walls.Therefore, while the concentration increased in the basal part of the PDC, the overhanging part became increasingly dilute, and the solid particles settled in traction BL from which stratified deposits formed in open spaces (Luongo et al., 2003).
Another characteristic, repeatedly found in the execution of numerical simulations, is the development, above the PDC tail and body, of convective thermal plumes, which influence the concentration trend.These low concentration plumes (less than 10 -4 ) bring small dimensions solid particles up in the atmosphere (to 1 km).A fall deposition of pyroclastic fragments occurs from them, causing the formation of vesiculated ash layers.Indeed, these layers form when a pause in the flow of the underlying PDC allows the deposition from the overhanging ash cloud.Ash thin levels formed exclusively by small particles of volcanic glass, and of coignimbritic ash strata between a depositional unit and the following one, classified as vesiculated ash layers, have been seen in the successions of Neapolitan Yellow Tuff (Cole and Scarpati, 1993).Moreover, in the deposits of the Plinian Vesuvius eruption, plane-parallel ash layers and well-stratified ash strata, also classified as vesiculated ash layers, are present (Luongo et al., 2003).
Finally, it must be noted that the formation of thermal plumes above the body and tail does not affect the BL behaviour.In other words, there is a net dissociation between the PDC boundary layer, where a kinetic-collisional regime develops, and the overlying diluted remnant part, in which a pure kinetic regime holds.Indeed, the model developed and applied to the study of the PDCs dynamics shows that the denser and basal underflow within the BL never outruns the upper more diluted suspension.This has been noticed in various eruptions, as the Neapolitan Yellow Tuff (Cole and Scarpati, 1993) and the Vesuvius Plinian (Luongo et al., 2003).
The simulations features here described are rather different from those obtained by Dartevelle et al. (2004).In the latter case, the head of the PDC is very high since the beginning, while the simulations reported here show that only at 180 s the head exceeds the rest of the PDC.In addition, at 80 s, the PDC head concentration has decreased by a 10 3 factor, while in the present simulations it remains almost constant.At 100 s, the basal underflow has outrun the head of the PDC and has a concentration of ∼ 30 %, while in our case the underflow never outruns the head and the basal concentration is ∼ 2.75 %.Finally, at 180 s, the underflow concentration is ∼ 55 % and so the BL shows plastic-frictional behaviour, while the underflow concentration we get is ∼ 3 ÷ 4 % and then the BL displays kinetic-collisional behaviour.Besides, the maximum height reached by the thermal plumes is rather greater than that attained in the present case.The differences found between the simulations here reported and those performed by Dartevelle et al. (2004) are mainly 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 cutoff 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 the small-scale is related to the large-scale.This difference has strong consequences on the overall simulation, and very dissimilar results are obtained, indeed.From the analysis of the solid particles concentration within the BL at 180s in the two cases, it derives that the forming deposit will show stratified facies in our simulations, and, on the contrary, 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 the Neapolitan Yellow Tuff (Cole and Scarpati, 1993), the Capodimonte and the Trentaremi eruptions (Cole et al., 1994), and the 79 AD Plinian Vesuvius eruption (Luongo et al., 2003).On the contrary, the simulations performed by using the LES model can be best applied to the study of the dynamics of Montserrat, Mount Pinatubo, and Lascar eruptions (Druitt, 1998;Calder et al., 2000).

Conclusions
Two-dimensional numerical simulations of the dynamics of fountains and associated PDCs have been performed by using a granular multiphase transient model, in order to validate and compare results here obtained with some phases of historical eruptions.PDCs 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 spatial and temporal description of the granular flow macroscopic dynamics, but also the depiction of sedimentation into the BL and deposition from it.On a large-scale, low concentration (≤ 5 × 10 -4 ) sectors of the flow lie in the upper part of the granular flow, above the fountain, and above the pyroclastic current tail and body as thermal plumes.The high concentration sectors, on the contrary, constitute the fountain and remain along the ground of the granular flow.Granular flows are therefore formed by a high-concentrated BL underlying a low-concentrated suspension.This structure of the flow is present from the proximal to the distal regions.On a small-scale, the dynamics in the BL is strongly affected by interactions between solid particles, whereas in the overlying dilute suspension it is controlled by the dragging of the gas phase that produces particles dispersion.Significant thermal disequilibrium effects are evidenced 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 BL implies that its dynamics is tractive, and that suggests a connection with the mechanisms leading to the formation of stratified facies, which can be cross or diffuse in relation to the values assumed by the concentration.A stratified facies is present when the concentration in the BL is sufficiently high to prevent turbulent sorting of fine ash.The results from numerical simulations appear to be in qualitative, and in some respect also quantitative, agreement with field survey observations collected from the deposits of the Neapolitan Yellow Tuff, the Capodimonte and Trentaremi eruptions, and the 79 AD Plinian Vesuvius eruption.Further understanding of the processes investigated should be attained by the development of more accurate multiphase transient models, as well as by the applications of these models to other well-documented eruptions.