Reproducing pyroclastic density currents deposits of the AD 79 eruption of the Somma-Vesuvius volcano using the box-model approach

In this study we use PyBox, a new numerical implementation of the box-model approach, to reproduce pyroclastic density current (PDC) deposits from the Somma-Vesuvius volcano (Italy). Our simplified model assumes inertial flow front dynamics and mass deposition equations, and axisymmetric conditions inside circular sectors. Tephra volume and density, 15 and Total Grain Size Distribution of EU3pf and EU4b/c, two well-studied PDC units from different phases of the AD 79 Pompeii eruption of Somma-Vesuvius (Italy) are used as input parameters. Such units correspond to the deposits from variably dilute, turbulent PDCs. We perform a quantitative comparison and uncertainty quantification of numerical model outputs with respect to the observed data of unit thickness, inundation areas, and grain size distribution as a function of the radial distance to the source. The simulations that we performed with PyBox were done considering: i) polydisperse 20 conditions, given by the total grain size distribution of the deposit, or monodisperse conditions, given by the mean Sauter diameter of the deposit; ii) round-angle axisymmetrical collapses or divided in two circular sectors. We obtain a range of plausible initial volume concentrations of solid particles from 2.5% to 6%, depending on the unit and the circular sector. Optimal modelling results of flow extent and deposit thickness are reached on the EU4b/c unit in a polydisperse and sectorialized situation, indicating that using total grain size distribution and particle densities as close as possible to the real 25 conditions significantly improve the performance of the PyBox code. The study findings suggest that the box model simplified approaches adopted have promising applications in constraining the plausible range of the input parameters of more computationally expensive models.

Although the 1D kinetic approaches cannot capture the multidimensional features of dynamics, they represent an important tool for several purposes.Firstly, it is practical to rely on simplified and fast numerical codes, which can be run 10 4 -10 6 times without an excessive computational expense, in order to produce statistically robust probabilistic hazard maps (Neri et al., 2015;Bevilacqua et al., 2017;Aravena et al., 2020).Furthermore, since 2D or 3D multiphase models require high computational times, often on the order of days or weeks for a single simulation, it is convenient to use simplified approaches, such as the box-model, in order to constrain the input space (Ogburn and Calder, 2017;Bevilacqua et al., 2019a).Finally, extensively testing the numerical models in a statistical framework, and evaluating the difference between model outputs and actual observations, also allows estimating the effect of the various modeling assumptions under uncertain input conditions (e.g., Patra et al., 2018Patra et al., , 2020;;Bevilacqua et al., 2019b).Model uncertainty is probably the most difficult class of epistemic uncertainty to evaluate robustly, but it is indeed a potentially large component of the total uncertainty affecting PDC inundation forecasts.
In this paper, we test the suitability of box-model approach (described in section 2.1 and in Appendix A), as implemented numerically in the PyBox code (Biagioli et al., 2019), by quantifying its performance when reproducing some key features of the remarkably well-known PDC deposits from one of the best studied and documented volcanic events, the AD 79 eruption of the Somma-Vesuvius (SV) volcano (detailed in section 2.2).It is indeed accepted that the box model is able to describe the main features of large-volume (VEI 6 to 8;Newhall and Self, 1982), low aspect ratio ignimbrites, whose dynamics is dominantly inertial (Dade and Huppert, 1996).However, the model has never been tested again PDC generated by VEI 5 Plinian eruptions.The procedure involves the calculation of the difference between model output and field data in terms of i) thickness profile, ii) areal invasion overlapping and iii) grain sizes (GS) volume fractions at various distances from the source.Similar approaches have been adopted in literature (Dade and Huppert, 1996;Kelfoun, 2011;Charbonnier et al., 2015).Tierz et al. (2016a,b) and Sandri et al. (2018) proposed a quantification of the uncertainty derived from the energy cone approach that relies on the comparison between invaded area and maximum runout of model output and field data.Our approach aims at the more detailed comparison of physical parameters (especially thickness and grain sizes, section 5), which allows a further investigation of the strengths and limitations of the PyBox model when used to simulate different PDC types (section 6).The PyBox code, is a numerical implementation of the box-model integral formulation for axisymmetric gravity-driven particle currents based on the pioneering work of Huppert and Simpson (1980).The theory is detailed in Bonnecaze et al. (1995) and Hallworth et al. (1998).The volume extent of gravity currents is approximated by an ideal geometric element, called "box", which preserves its volume and geometric shape class, and only changes its height/base ratio through time (see Figure 1).The box does not rotate or shear, but only stretches out as the flow progresses.In this study the geometric shape of the box is assumed to be a cylinder, or cylindrical sector, i.e. we assume axisymmetric conditions.The model describes the propagation of a turbulent particle-laden gravity current, i.e. a homogeneous fluid with suspended particles.Inertial effects are assumed to have a leading role with respect to viscous forces and particle-particle interactions.
Thus, the particle sedimentation is modelled and modifies the current inertia during propagation.In this study we assume the classical dam break configuration, in which a column of fluid instantaneously collapses and propagates, under gravity, in a surrounding atmosphere with uniform density ρ atm .Other authors (Bonnecaze et al., 1995;Dade andHuppert, 1995a, b, 1996) have instead considered gravity currents produced by the constant flux release of dense suspension from a source.Our approach differs from considering constant stress acting on the basal area (Dade and Huppert, 1998).Constant stress dynamics have been further explored in literature, and they can lead to different equations if the basal area grows linearly or with the square of the radius (Kelfoun et al., 2009;Kelfoun, 2011;Ogburn and Calder, 2017;Aspinall et al., 2019).PDCs are driven by their density excess with respect to the surrounding air: the density of the current ρ c is defined as the sum of the density of an interstitial gas, ρ g , and the bulk densities of the pyroclasts carried by the flow, (ρ s i ) i=1,…,N .In this study we assume ρ atm ≠ ρ g , i.e. the interstitial gas is hotter than surrounding atmosphere, differently from Neri et al. (2015) and Bevilacqua et al. (2017).The code allows ρ atm > ρ g , but thermal properties remain constant for the duration of the flow, and in this study we assumed ρ atm = ρ g .The thermodynamics of cooling effects are explored in Bursik and Woods (1996).A proper way to express the density contrast between the current and the ambient fluid is given by the reduced gravity  ′ , that can be rewritten in terms of the densities and the volume fractions described above (see Biagioli et al., 2019).That said, we make some additional simplifying hypotheses.First of all, we assume that the mixture flow regime is incompressible and inviscid, since we assume that the dynamics of the current is dominated by the balance between inertial and buoyancy forces.
The assumption of incompressibility implies that the initial volume V 0 remains constant.Moreover, we assume that, within the current, the vertical mixing, due to turbulence, produces a vertically uniform distribution of particles.The particles are assumed to sediment out of the current at a rate proportional to their constant terminal (or settling) velocity (w s i ) i=1,…,N .Once deposited, they cannot be re-entrained by the flow; the converse is explored in (see Fauria et al., 2016).Finally, surface effects of the ambient fluid are neglected.
Under these hypotheses, the box model for particle-laden gravity currents states that the velocity of the current front (u) is related to the average depth of the current (h) by the von Kármán equation for density currents  = ( ′ ℎ) 1 2 , where Fr is the Froude number, a dimensionless number expressing the ratio between inertial and buoyancy forces (Benjamin, 1968;Huppert and Simpson, 1980) and  ′ is the reduced gravity.In addition, we assume that particles can settle to the ground and this process changes the solid particle fractions (ε i ) i=1,…,N .
The box model for axisymmetric currents thus reads: Eq. ( 1) By solving these equations, we computed the amount of mass loss by sedimentation, per unit area, per time step, for each particle class.So, the thickness profile of the i th particle class is the ratio of the i th deposited mass to the i th solid density multiplied by the packing fraction α measured in the deposit.More details on the numerical solver are provided in Appendix A.
In the calculation of the region invaded by a PDC, first we calculate the maximum flow runout over flat, i.e. the distance at which ρ c = ρ atm .The flow stops propagating when the solid fraction becomes lower than a critical value, and the remaining mixture of gas and particles lifts off, possibly generating a phoenix cloud if hot gas is assumed.In the case of monodisperse systems there are analytical solutions for the maximum flow runout (Bonnecaze et al., 1995;Esposti Ongaro et al., 2016;Bevilacqua, 2019).Then, once set a vent location, we assess the capability of topographic reliefs to block the current.In particular, the invasion areas are obtained by using the so-called energy-conoid model, based on the assumption of nonlinear, monotonic decay of flow kinetic energy with distance (Neri et al., 2015;Bevilacqua, 2016;Esposti Ongaro et al., 2016;Bevilacqua et al., 2017;Aspinall et al., 2019;Aravena et al., 2020).In more detail, first we determine the maximum height ℎ  of an obstacle the flow can overcome.Then we compare the kinetic energy of the current front and the potential energy associated to the obstacle top.In this approach we are neglecting returning waves.When investigating the current flow on complex topographies, we finally consider that the flow may start from positive elevation or encounter upward slopes after downward slopes.In this case, we compare h max at a given distance from the vent and the difference in level experienced by the current between the previous and the present sampled positions.
In the PyBox code, the main input parameters are summarized by: a) the total collapsing volume (expressed in terms of the dimension of the initial cylinder/rectangle with height=ℎ 0 and radius/base= 0 ); b) the initial concentration of solid particles, subdivided (for polydisperse simulations) into single particle volumetric fractions (ε 0 ), with respect to the gas; c) the density of single particles ρ s ; d) ambient air density (ρ atm =1.12 kg/m 3 ) and gravity current temperature; e) Froude number of the flow, experimentally measured by Esposti Ongaro et al. (2016) as Fr=1.18;g) gravity acceleration (g=9.81 m/s 2 ).With respect to points b) and c), more details are provided in section 3.2.
The PDC deposits of the EU3pf unit record the phase of total column collapse closing the Plinian phase of the eruption.
These are ca. 1 m thick on average, radially dispersed (up to 10 km from vent area) and moderately controlled by local topography, so resulting in a complex vertical and lateral facies variability (Gurioli, 1999;Gurioli et al., 1999) possibly related to local variation in turbulence, concentration and stratification of the current.Median clast size tends to diminish gradually from proximal to distal locations, and the coarsest deposits (generally present as breccia lenses in the EU3pf sequence) are located within paleodepressions.Gurioli et al. (1999) pointed out that: i) in the southern part of the SV area the relatively smooth paleo-topography controlled only locally the overall deposition of this PDC; ii) in the eastern sector of SV, the interaction of the current with the ridge representing the remnants of the old Mount Somma caldera (Fig. 2a) possibly triggered a general increase of the current turbulence and velocity and a more efficient air ingestion which resulted in the local deposition of a thinly stratified sequence; iii) in the western sector of SV, the presence of a breach in the caldera wall and of an important break in slope in the area of Piano delle Ginestre (Fig. 2a), possibly increased deposition from the PDC, producing a large, several meters thick depositional fan toward the sea-facing sectors (like in Herculaneum -Fig.2a); iv) in The AD 79 EU4 marks the reappraisal of the eruption after the end of the Plinian phase and was related by Cioni et al. (1999) to the onset of the caldera collapse.This complex unit has been subdivided into three distinct layers (Cioni et al., 1992): a thin basal fallout layer ("EU4a"), a PDC deposit derived from the collapse of the short-lived column that emplaced the EU4a layer ("EU4b"), and the products of the co-ignimbritic plume mainly derived by ash elutriation from the current depositing EU4b ("EU4c").Gurioli (1999) illustrates how the EU4 unit is furthermore complicated, since it actually presents a second fallout bed ("EU4a2") interlayered within the level "EU4b".This fallout bed can be clearly recognized only in distal sections of the southern sector, while in the north and in the west it is represented by a discontinuous level of ballistic ejecta.Level "a2" divides level "b" in two parts, which are approximately 2/3 (the lower one) and 1/3 (the upper one) of the total thickness of level "b" (Gurioli, 1999).Runout of the EU4b PDC is one of the largest runouts observed for the SV PDCs; it was maximum toward south (up to about 20 km from vent area; Gurioli et al., 2010).This unit has been extensively studied by Gurioli (1999) who highlighted that the high shear rate exerted by the EU4b is clearly evidenced by the formation of "traction carpets" and local erosion of the pumice-bearing layer of the underlying EU4a.The EU4b deposit can be interpreted as derived from a short-lived sustained, unsteady, density stratified current.From a sedimentological point of view, EU4b shows clear vertical grain size and textural variations, from cross bedded, fine lapilli to coarse ash laminae at the base up to a massive, fine ash-bearing, poorly sorted, matrix-supported bed at top (Gurioli, 1999).During deposition of EU4b, ash elutriated from the current formed a convective plume dispersed from the prevailing winds in a south-eastern direction, which deposited EU4c mainly by fallout.The clear field association of these two deposits (indicated as EU4b/c) gives here the uncommon possibility to evaluate with a larger accuracy two of the most important PDC source parameters: erupted volume and TGSD.

Model input parameters and field data for comparison
The main properties of the EU3pf and EU4b/c units, i.e. thicknesses/total volume, maximum runout and total grain size distribution -TGSD, have been calculated in Cioni et al. (2020) and partially processed to fit with PyBox input requirements (see sections 3.1.1and 3.1.3).Densities of single grain sizes (Barberi et al., 1989; section 3.1.2)and emplacement temperatures of PDCs (T=600 K for both EU3pf and EU4; Cioni et al., 2004) are derived from the literature.
In summary, the total volume, the TGSD, the densities and the temperature obtained from the field data are used as the main inputs of PyBox.Thus, the model produces several outputs: (i) mean unit thickness as a function of the radial distance from the source, (ii) inundated area, (iii) grain size distribution as a function of radial distance from the source.All these outputs https://doi.org/10.5194/se-2020-138Preprint.Discussion started: 9 October 2020 c Author(s) 2020.CC BY 4.0 License.
are finally compared to the corresponding field data.The initial volumetric fraction ε 0 of the solid particles over the gas is the main tuning parameter that is explored to fit the outputs with the field data.This procedure is repeated under monodisperse and polydisperse conditions, and by performing round-angle axisymmetrical collapses or sectorialized collapses, i.e. divided in two circular sectors with different input parameters.

Thickness, maximum runout and volumes
Cioni et al. ( 2020) recently revised and elaborated a large amount of field data from EU3pf and EU4b/c (106 and 102 stratigraphic sections, respectively), tracing detailed isopach maps and defining the maximum runout distance (the ideal 0 m isopach) and the related uncertainty.Given the objective difficulty to trace the exact position of a 0 m isopach for the deposit of a past eruption, Cioni et al. (2020) proposed to define three different outlines of PDC maximum runouts, namely "5 th percentile", "50 th percentile" and "95 th percentile" (called Maximum Runout Lines, MRLs), basing on the uncertainty associated to each segment of the proposed 0 m isopach.The MRLs of EU3pf and EU4b are shown in Figure 3c and 3d respectively.
Cioni et al. ( 2020) also calculated the volumes of both EU3pf and EU4b/c, using these maps to derive a digital elevation model of the deposits with the triangular irregular network (TIN) method (Lee and Schachter, 1980).In this study, we considered volume estimations (Table 1) related to the MRL 50 , i.e. the 50 th percentile of the maximum runout distance.
Given the asymmetric shape of unit EU4b/c and, partially, of unit EU3pf, we have calculated also the volumes dividing each unit in two circular sectors: N and S for EU3pf, NW and SE for EU4b/c.These subdivisions have been also used to calculate the related TGSDs (see unit 3.1.3)and to perform sectorialized simulations (see section 4). Figure 3c,d

Density data
In order to provide density values for each GS, we used the mass fractions of the different components (juveniles, lithics and crystalssee Table S1 from the Supporting Information) calculated by Gurioli (1999).Such values were associated to the averaged density measurements for these three components presented in Barberi et al. (1989), through which we extrapolated the weighted mean (with respect to mass fraction) mean density of each grain size class for both EU3pf and EU4b/c units (Table 2).Table 2. Calculated mean densities for each grain size for both the EU3pf and EU4b/c units.

Grain size data: total grain size distribution (TGSD) and mean Sauter diameter (MSD)
The  In the simulations under monodisperse conditions, we used the value of mean Sauter diameter (MSD) of the volumetric TGSD (e.g., Neri et al., 2015).According to Fan and Zhu (1998), the Sauter diameter of each particle class size is also called d 32 (see also Breard et al., 2018), and it is the diameter of a sphere having the same ratio of external surface to volume as the particle, which is given by: Eq. ( 4) where V is the particle volume, S is the particle surface, d v is the diameter of a sphere having the same volume as the particle and d s is the diameter of a sphere having the same external surface as the particle.In order to obtain a value for the MSD instead, given a deposit sample divided in N grain size classes, we have initially calculated the number of particles of each grain size i=1,…,N, that is: Eq. ( 5) where   is the cumulative volume of the i-th grain size class, and   is the radius of the i-th grain size.The mean MSD is finally derived as: Eq. ( 6) where   and   are the diameters of, respectively, the i-th and j-th grain sizes.

Comparison between field data and simulation outputs
Since the PyBox code assumes axisymmetric conditions, the thickness outputs are equal along all the radial directions of the collapse, and only vary as a function of the distance to the source.These output data were compared with the mean radial profiles of unit thickness (for both EU3pf and EU4b/c) as derived from the digital models of deposit in Cioni et al. (2020).
For building the radial profiles, the average thickness was estimated over concentric circles drawn with a 100-m step of distance.The radial thickness profiles were drawn starting from a distance of 3 km from the vent, as no thickness data are available for sites closer than 3 km.Due to the lack of reliable data, we have moreover excluded from our analyses the portions of the circles located in marine areas.In order to describe the variation range of the thicknesses of the deposits, we are providing minimum and maximum thicknesses along each circle in Appendix B (Fig. A1).
Concerning the inundation area, the methodology adopted is similar to the one used by Tierz et al. (2016b) and relies on the The False Negative case (no deposit nor simulation) has not been obviously calculated.In statistical literature, the True Positive value is also called Jaccard Index of similarity (Tierz et al., 2016b;Patra et al., 2020).While the TP/TN/FP approach, and in general the Jaccard Index, focus on the areal overlapping, other metrics can specifically focus on the distance between the boundaries of the inundated areas, i.e. the Hausdorff distance, detecting and comparing channelized features in the deposit (Aravena et al., 2020).However, PyBox is not specifically aimed at the replication of such features, and we focus on the areal overlapping properties.The scarcity of stratigraphic sections in the N sector (for the EU3pf unit) and the NW sector (for the EU4b/c unit) negatively affects the availability of comparisons with respect to volume fractions, which are forcedly limited to sections at 4 km of distance from the hypothetical vent area, most of which have been collected at the bottom of paleovalleys.Moreover, for the EU3pf unit, even in the S sector the available samples are mostly concentrated in the area of Herculaneum (5 samples).

Results
The results of 6 simulations (4 for the EU3pf unit and 2 for the EU4b/c unit) are discussed here (see Table 4 for the main input parameters).These simulations are the result of an extensive investigation in which a wide range of different values of ε 0 have been tested, following a trial-and-error procedure aimed at reproducing more closely the thickness profile of the deposit.In particular, we performed several simulations varying ε 0 between 0.5% and 6% (for EU3pf) and between 0.1% and 5% (for EU4b/c).The values in Table 4 represent the optimal combinations.simulations performed using the unmodified DEM did not produce major differences.

Unit Simulation
In the EU3pf case study, we performed both axisymmetrical simulations over a round-angle (given the quasi-circular shape of the deposit) and also axisymmetrical-sectorialized simulations to investigate possible sheltering effects of the Mount Somma scarp (Fig. 2a).In particular, two distinct column collapses, one to the N and the other to the S, each of which has a collapsed volume corresponding to the actual deposit volume in that sector.In the EU4b/c case study, we performed only axisymmetrical-sectorialized simulations, to reproduce more closely the dynamics of the related collapse, as indicated by the different dispersal in the NW and SE sectors of the PDC deposit.In particular, two distinct collapses for the same simulation, one to the NW and the other to the SE) In summary we provide: a) the thickness comparison between deposit and modelled results (Figure 6) and between simulations done with different initial volumetric fraction of solid particles (ε 0 -Figure 7); b) the inundation areas, including the quantitative matching of simulations and actual deposit (Figure 8 and Table 5); c) the grain size distribution comparison, between deposit and modelled values, i.e. the volume fractions of ash vs lapilli (Figure 9) and of all the grain size classes (Figure 10).

General considerations
Testing PyBox with respect to field data is aimed at two main objectives: i) quantifying the degree of reproduction of the real PDC deposit of Plinian eruptions in terms of thickness, inundation area and grain size and ii) evaluating the reliability of the code when considering different assumptions, i.e. polydisperse vs. monodisperse situations, and 360° axisymmetric conditions vs dividing circular sectors.Before commenting our results, two main general considerations, in common for both EU3pf and EU4b/c, deserve a special discussion.

Runout truncation and non-deposited material
We remind that PyBox code produces the map of the inundated area (Neri et al., 2015;Bevilacqua, 2016), by truncating the runout wherever the kinetic energy of the flow is lower than the potential energy associated to a topographic obstacle (Section 2.1 and Appendix A).In this way, however, the material that lies beyond the truncation is neither redistributed nor considered any more.However, depending on the topography in our case study, this amount of material is not extremely high.For instance, EU4_poly_AS (Table 4), in its SE part, has several truncations due to the intersection of the decay function of kinetic energy with several topographic barriers, i.e. the Apennines to the ENE and the Sorrentina Peninsula to the SE (Figs. 2 and 7).For the whole SE part of the deposit, the topographic barriers are located between 11.85 km and 19.25 km from vent area, with a mean vale of 15 km.If we truncate PyBox deposit in correspondence of these three limits, the non-deposited volume is between 3.46x10 6 m 3 (cut at 19.25 km) and 2.3x10 7 m 3 (cut at 11.85 km), with a mean value of 1.27x10 7 m 3 (cut at 15 km).Considering that the volume collapsed to the SE is 1.5x10 8 m 3 , the non-deposited volume corresponds therefore to a value between 2% and 15%, with a mean of 8%.The amount of volume effectively "lost" is relatively small, also considering that the total volume of the collapsing mixture is inclusive of the EU4c unit (co-ignimbritic part).However, further development of the code might consider a strategy to redistribute this non-deposited material (e.g., Aravena et al., 2020).

Initial volumetric fraction of solid particles
The value of the initial volumetric fraction of solid particles (ε 0 ) in the PDC represents one of the most uncertain parameters, for which few constraints exist.Recently, Valentine (2020) performed several multiphase simulations using mono-or bidisperse distributions to investigate the initiation of PDCs from collapsing mixtures, and to derive criteria to determine when either a depth-averaged model or a box-model are best suited to be employed for hazard modelling purposes.The author concluded that, among other factors (e.g.impact speed or relative proportion of fine to coarse particles), a volumetric concentration of particles of around 1% (slightly lower than those used in this paper) is generally capable of producing a dense underflow and a dilute, faster overriding flow.For such cases, and considering an impacting mixture consisting of at least ca.50% coarse particles (> 1 cm diameter) relative to fines (< 1 cm diameter), Valentine (2020) suggests that a depthaveraged granular flow model well approximates such PDCs, and could be reasonably used for hazard assessment purposes.For the units studied here, the sedimentological features show that there are clear evidences of the formation of a dense underflow in, respectively, the N part of Somma-Vesuvius volcano (EU3pf unit; Gurioli et al., 1999) and in correspondence of the urban settlements of Herculaneum and Pompeii (EU4 unit; Cioni et al., 1999;Gurioli et al., 2002).We however think the employment of a box-model is justified for at least the unit EU4b/c, which can be considered intermediate between a dilute, turbulent and a granular concentrated current, in the sense of Branney and Kokelaar (2002), but closer to the dilute end-member type.In this view, the box-model can be effectively employed to describe the overriding dilute part units similar to the EU4, following a two-layer approach (Kelfoun, 2017;Valentine, 2020).
For the box-model used here, it should be kept in mind the variation of the ε 0 value might have an important effect on the simulated deposit thicknesses, as seen in Figure 7.In both units, in fact, the model results for thickness at the beginning of the simulated area (i.e. 3 km from vent area), vary from ca. 1 to ca. 2 m (for EU3pf) or from ca. 1.2 m to ca. 3.6 m (for EU4b/c) if ε 0 is varied, respectively, from 1.5% to 6% and from 0.5% to 5%.

Thickness comparison
The first parameter that we compare between the deposit and the modelled results is the thickness variation with the distance to the source, an approach already adopted, for instance, by Dade and Huppert (1996).
Our comparison focuses on the average thickness calculated over concentric circles drawn with a 100-m step of distance.
However, the thickness variation of the deposit in different radial directions describes two different situations for the EU3pf and EU4b/c units and deserves a brief discussion, detailed in Appendix B.
The average thickness of the EU3pf deposit mean profile initially shows an increasing trend (between 3 to 4 km to the N and between 3 to 6 km to the S -Fig.6a) followed by a slow, constant decrease.This situation could highlight a lower capability of the current to deposit in more proximal areas, allowing the mass to be redistributed toward more distal sections.This could also be motivated with a spatial variation of the PDC flux regime, which was more turbulent in proximal areas than in distal ones, as also testified by the abundance of lithofacies typical of dilute and turbulent PDCs (//LT to xsLT; see Fig. 2b and Gurioli et al., 1999).Instead, the spatial homogeneity of lithofacies for the EU4b/c unit (Cioni et al., 1992) suggests a higher uniformity of its parent PDC.Moreover, the trend of the mean deposit thickness profile has a steep and rapid decrease of thickness up to 5-6 km, followed (after a break in slope) by a "tail" with an increasing gentler decrease of thickness.This peculiar trend is in agreement with the lithofacies association in the unit EU4b/c (Cioni et al., 1992), which indicate a progressive dilution of the current through time and a progressive aggradation of the deposit.
That said, the degree of matching between the modelled and the real thickness of the EU3pf unit is less accurate than in the EU4b/c case study.However, the mean thickness profile of the actual deposit is roughly parallel with the model, in some parts.Under polydisperse conditions, PyBox does not improve its performance in replicating the thickness profile of EU3pf.
The difficulties of PyBox in reproducing the thickness average profile testifies the strongest interaction of the EU3pf unit with the non-homogeneous topography (see also Gurioli, 1999;Cioni et al., 2020) and the likely dominant role of the density stratification and granular transport in the deposition process.To the North there was in fact an extremely rough topography, https://doi.org/10.5194/se-2020-138Preprint.Discussion started: 9 October 2020 c Author(s) 2020.CC BY 4.0 License.
similar to the present one, where the interaction of the PDC with the surface produced largely variable lithofacies.To the South instead there was a gentler topography, with a topographic high on which the town of Pompeii (see Fig. 2a) was built.This latter aspect is also evident from Vogel and Märker (2010), who reconstructed the pre-AD79 paleotopography of the plain to the SE of the SV edifice.From this work, it is possible to appreciate how the modelled depth of the pre-AD79 surface is 0-1 m lower with respect to the present surface in correspondence of the present town of Pompeii and the ancient Pompeii excavations (due to the presence of piles of tephra fallout deposits up to 2 m thick), while it is up to 6-7 m deeper to the NW of these sites.
The thickness comparison of the EU4b/c unit, on the contrary, suggests that this unit was likely deposited under inertial flow conditions, dominated by turbulent transport.The SE "tail" part of the deposit is particularly very well reproduced by the polydisperse simulations, where the simulated profile is almost coincident with the deposit profile (Fig. 6bright).
Conversely, to the NW the modelled thickness in the initial part overestimates a bit the real deposit (Fig. 6b).The polydisperse simulations (blue dashed lines in Fig. 6b) are much closer to the measured trend of the mean thickness profile than under the monodisperse conditions (i.e.MSD), demonstrating the key role of the grain-size distribution in gas-particle turbulent transport.

Comparison of inundated areas
The areal overlapping between the model output area and the actual deposit (True Positive -TP) is discussed together with the quantification of model overestimation (False Positive -FP) and underestimation (True Negative -TN).In Table 5 we also provided the TP/FP/TN estimates also for the 5 th and 95 th percentiles of the maximum runout lines (MRLs), i.e. a measure of the spatial uncertainty affecting the actual deposit.We remark that the TN instances could be interesting from a hazard point of view, because they actually represent the underestimation of the model: a conservative approach is therefore to use the lowest value of the TN instances as a threshold to evaluate the reliability of a model.
As said above, the polydisperse simulations of the EU3pf unit poorly fit with the deposit thickness, and the inundated area is significantly larger than the deposit area.Thus, they are not included in the quantitative estimation of area match/mismatch.For instance, while the maximum runouts of the deposit are on the order of 8-10 km, the maximum runout given by the model (in absence of topography) is ca.13-15 km.The monodisperse simulations perform better, in this sense, and maximum runouts are slightly different (ca.7-10 km) from the real ones: for this reason, only the monodisperse simulations for the EU3pf case have been considered in Figure 8 and Table 5.More precisely, the axisymmetrical (EU3pf_mono_AX) and the sectorialized (EU3pf_mono_AS) share a similar degree of TP instances (between 63% and 75% -Table 5), but have opposite properties for what concerns overestimation/underestimation.EU3pf_mono_AX has in fact a higher tendency to underestimate (FP < TN -Table 5) while EU3pf_mono_AS tends to overestimate the actual deposit (FP > TN -Table 5).
For what concerns the EU4b/c simulations (Fig. 8), we report the quantitative matching of both the simulations under polydisperse and monodisperse conditions.The most striking feature that could be seen from Figure 8 is that, while to the SE https://doi.org/10.5194/se-2020-138Preprint.Discussion started: 9 October 2020 c Author(s) 2020.CC BY 4.0 License.a good matching is obtained, to the NW the polydisperse simulation tends to sensibly overestimate the inundation area.
Conversely, the monodisperse simulation is more equilibrate between NW and SE.This could be motivated, for the SE part, with the surrounding morphology of the Sorrentina Peninsula and the Apennines, which acts as a natural barrier, and, for the NW sector, to the absence of morphological constraints especially to the N.The results presented in Fig. 8 and Table 5 show that the TP values for the simulation EU4_poly_AS are in the interval 61%-73%, while TN values range from 0.7% and 2%, and FP values range from 24% and 38%.Thus, while the degree of overlapping between model and deposit is at an acceptable value and the percentage of model underestimation is below 2%, the model tends to appreciably overestimate the median outline of the deposit.On the opposite, the simulation EU4_mono_AS shows the highest TP values (73%-80%) and the lowest FP (3%-8%).Despite these better performances, it should be always reminded how the thickness profile is less accurate under monodisperse conditions.We remark that, beyond 14 km (ca.2-3 km beyond the deposit MRL 95 ) the thickness provided by the model under polydisperse conditions is < 1 mm (see Figure 8c).Shallow deposits might be possibly affected by erosion, and the actual deposit in the NW sector might in fact resemble the PyBox results.We also remark that the MRLs defined by Cioni et al. (2020) have been defined up the 95 th percentile, meaning that there is still a 5% chance that the actual MRL could be placed further away from the source.This is very significant in the NW part of the EU4b/c deposit, where no or very few outcrops can be found beyond 5-6 km from vent area.

Grain size comparison
Finally, we consider the volume fraction of the grain sizes of the actual deposits versus those derived from PyBox.We present the results in two different ways.Firstly, we provide a general overview of what are the relative proportions of ash/lapilli with distance to the source (Fig. 9), and then we provide more complete volumetric grain size comparisons for each Φ unit (Fig. 10).We note that this comparison is one of the most uncertain because of some inherent epistemic uncertainties in the data, that are: i) the complete lack of ultra-proximal sites possibly enriched in coarse grained particles that influenced the calculated TGSD; ii) the fact that the sections used for TGSD calculation and data comparison are (for both units) located mainly along the aprons of the volcano, in many cases in correspondence of the lower parts of valleys or paleovalleys.This could have led to have an under-representation of the finer-grained deposits located in high or paleo-high morphological locations.
The data presented in Fig. 9 confirms the differences between EU3pf and EU4b/c.EU3pf (Fig. 9a) shows that the simulated and real volumetric contents of ash/lapilli are similar only up to 4 km (both to the N and to the S).Then, the relative proportions of ash/lapilli in the simulations indicate that, after 6 km, the simulated grain sizes are made almost entirely (>90%) by ash, with a sensitive difference with respect to field data (only to the S, as to the N there are no available measurements).The most extreme situation could be seen at 9 km, where the modelled grain sizes are composed for > 80% in volume by the two finest ones (4Φ-5Φ), while deposit data indicate a more equal distribution of grain sizes.In Fig. 10a we https://doi.org/10.5194/se-2020-138Preprint.Discussion started: 9 October 2020 c Author(s) 2020.CC BY 4.0 License.observe that at 4 km (both N and S) the grain size distributions are similar between the actual deposit and the model, although there is a shift of ca. 2 Φ toward the finer grain sizes in the modelled data.
For the EU4b/c unit, we observe that the general proportions between ash and lapilli (Fig. 9b) are more similar between the model and the deposit (especially at 4 km from the vent area to the N).However, in Fig. 10b we see that at 4 km to the N, the situation is opposite to EU3pf, since the modelled grain size is richer in coarse particles than the actual deposit.Such difference might be motivated with the above-mentioned roughness of the topography, which might favour the deposition of coarser particles at locations < 4 km.In the SE sector the differences between modelled and observed grain sizes are lower at 6 km and 9 km distance to the source, while are greater at 14 km and 20 km, where the 2 finest modelled grain sizes account for > 80% of the volume.

Conclusions
We have evaluated the suitability of the box-model approach implemented in the PyBox code to reproduce the deposits of EU3pf and EU4b/c, two well-studied PDC units from different phases of the AD 79 Pompeii eruption of Somma-Vesuvius (Italy).The total volume, the TGSD, the grain densities, and the temperature obtained from the field data are used as the main inputs of PyBox.The model produces several outputs that can be directly compared with the inundation areas and radially-averaged PDC deposit features, namely the unit thickness profile and the grain size distribution as a function of the radial distance to the source.We have performed simulations either under polydisperse or monodisperse conditions, given by, respectively, the total grain size distribution and the mean Sauter diameter of the deposit.We have tested axisymmetrical collapses either round-angle or divided in two circular sectors.
The initial volumetric fraction ε 0 of the solid particles over the gas is the main tuning parameter (given its uncertainty) that is explored to fit the outputs with the field data.In this study, we obtained the best fit of deposit data with a plausible initial volume concentration of solid particles from 3% to 6% for EU3pf (depending on the circular sector) and of 2.5% for EU4b/c.These concentrations optimize the reproduction of the thickness profile of the actual deposits.
Concerning the EU3pf unit: 1) the average thickness of the EU3pf deposit initially shows an increasing trend, from 3 to 4 km to the N and from 3 to 6 km to the S, followed by a slow, constant decrease.The simulated thickness poorly resembles the actual deposit, although the maximum values are comparable and the two profiles are roughly parallel, in some parts.Under polydisperse conditions, PyBox does not improve its performance in reproducing the thickness profile of EU3pf; 2) in the monodisperse simulations of EU3pf the maximum runouts are slightly different from the real ones, but overall consistent.
The round-angle and sectorialized simulations share a similar degree of TP instances (between 63% and 75%), but have opposite properties for what concerns overestimation/underestimation.The round-angle axisymmetric simulation underestimates the actual deposit (FP < TN) while the sectorialized simulation overestimates the actual deposit (FP > TN); 3) the simulated and real volumetric contents of ash/lapilli in EU3pf are similar only up to 4 km.Then, the relative proportions https://doi.org/10.5194/se-2020-138Preprint.Discussion started: 9 October 2020 c Author(s) 2020.CC BY 4.0 License. of ash/lapilli in the simulations indicate that the simulated grain sizes are made almost entirely (>90%) by ash, with a sensitive difference with respect to field data after 6 km.We observe that at 4 km the grain size distributions are similar between the actual deposit and the model, although there is a shift of ca. 2 Φ toward the finer grain sizes in the modelled data.
Concerning instead the EU4b/c unit: 1) this unit has a steep and rapid decrease of thickness up to 5-6 km, followed, after a break in slope, by a "tail" with a gentler decrease of thickness.The polydisperse box-model simulations are much closer to the measured trend of the mean thickness profile than under the monodisperse conditions.The SE thickness profile of the polydisperse simulation is almost coincident (within the uncertainty range) with the corresponding part of the deposit (specifically after 6 km, and with a ca.0.5 m overestimation between 3.5-6 km), while to the NW the modelled thickness slightly overestimates the real deposit in the initial part (up to ca. 6 km); 2) in the simulations of EU4b/c, a good matching of inundated area towards the South-East is obtained.Towards the North-West the polydisperse simulation sensibly overestimates the inundation area.On the opposite, the simulation under monodisperse conditions shows the highest TP values (73%-80%) and the lowest FP (3%-8%).However, the thickness profile is less accurate under monodisperse conditions.Moreover, shallow deposits in the NW sector might be possibly affected by erosion, and the actual deposit in the NW sector might in fact resemble the PyBox results obtained under polydisperse conditions; 3) the general proportions between ash and lapilli in EU4b/c are similar between the model and the deposit.However, at 4 km to the N, the situation is opposite to EU3pf, since the modelled grain size is richer in coarse particles than the actual deposit.In the SE sector the differences between modelled and observed grain sizes are lower at 6 km and 9 km distance to the source, while are greater at 14 km and 20 km, where the 2 finest modelled grain sizes account for > 80% of the volume; 4) in the SE sector, because of model runout truncation, we evaluated an average non-deposited volume of 1.2710 7 m 3 (cut at 15 km).Considering that the volume collapsed to the SE is 1.510 8 m 3 , the average non-deposited volume corresponds therefore to a value of 8%.
Thus, the amount of volume effectively "lost" with the PyBox approach is relatively small, also considering that the total volume of the collapsing mixture is inclusive of the co-ignimbritic part.
Pyroclastic density currents generated by Plinian eruptions span over a wide range of characters, and can display very different behaviour and interaction with the topography.During the AD 79 eruption of Somma-Vesuvius, two PDC units, despite both emplaced after column collapses, display significantly different sedimentological features and should likely be better described by different models.The study findings indicates that the box-model, which is suited to describe turbulent particle-laden inertial gravity currents, well describes the EU4b/c PDC unit but is not able to accurately catch some of the main features of the EU3pf unit.This is probably due to its strongly stratified features, which make the interaction with the topography of the basal concentrated part of the flow a controlling factor in the deposition process.In particular, results highlight again the key role of the grain-size distribution in the description of inertial PDCs: while the final runout is mostly controlled by the finest portion of the distribution, the total grain size distribution strongly affects the thickness profile (e.g., Fig. 6b) and it is an essential ingredient for proper modelling of the PDC dynamics.

Figure 1 .
Figure 1.a) Schematic diagram of an inertial gravity current with a depth h c , flow front velocity u c and density ρ c in an ambient fluid of density ρ 0 (modified from Roche et al., 2013); b) Evolution of channelized currents through a series of equal-area rectangles, according to the model (hence the name "box model").
https://doi.org/10.5194/se-2020-138Preprint.Discussion started: 9 October 2020 c Author(s) 2020.CC BY 4.0 License.Our model consists of a set of ordinary differential equations, that provide the time evolution of flow front distance from the source, l(t), together with the current height h(t) and the solid particle volume fraction (ε i ) i=1,…,N , N being the number of particle classes considered.The volume fractions refer to a constant volume of the mixture flow, not reduced by the deposition.

Figure 2 .
Figure 2. a) Location of the Somma-Vesuvius volcano.Coordinates are expressed in the UTM WGS84-33N system; b) the https://doi.org/10.5194/se-2020-138Preprint.Discussion started: 9 October 2020 c Author(s) 2020.CC BY 4.0 License. the northern sector of SV, the deeply eroded paleotopography (with many valleys cut on steep slopes) favoured the development within the whole current of a fast-moving, dense basal underflow able to segregate the coarse, lithic material and to deposit thick lobes in the main valleys, and of a slower and more dilute portion travelling and depositing thin, stratified beds also on morphological highs.

Figure 3 .
Figure 3. Thicknesses and isopach lines for the a) EU3pf and b) EU4b/c units; MRLs of the c) EU3pf and d) EU4b/c units.Inferred position of AD 79 vent (red triangle) and SV caldera outline (dark orange dashed line) after Tadini et al. (2017).Light green dashed lines delimit the sectors (N-S for EU3pf and NW-SE for EU4) of the different column collapses.Background DEM fromTarquini et al. (2007).220 TGSD estimations are necessary to do simulations under polydisperse conditions.The present version of PyBox takes as input the volumetric TGSD (i.e. in terms of volumetric percentages), while TGSD data from Cioni et al. (2020) are in weight percentages.These latter values have been therefore converted into volumetric percentages by considering the abovementioned densities (section 3.1.3). Figure 4 displays the volumetric TGSDs employed for EU3pf (Total, N and S) and the EU4b/c (Total, NW and SE).

Figure 4 .
Figure 4. Volumetric total grain size distributions for the EU3pf and EU4b/c units.
approach described byFawcett (2006)  and implemented byCepeda et al. (2010) for landslide deposit back-analysis.This method is based on the quantification of the areal overlapping between the measured deposit (true classes) and the modelled deposit (hypothesized classes) (Figure5).In particular, we quantify: a) the areal percentage of model intersecting the actual deposit (true positive -TP); b) the areal percentage of model overestimating the actual deposit (false positive -FP); c) the percentage of model underestimating the actual deposit (true negative -TN).More precisely:

Figure 5 .
Figure 5. Sketch representing the three areas used for the validation procedure (the model output outline is drawn in dashed black line).
https://doi.org/10.5194/se-2020-138Preprint.Discussion started: 9 October 2020 c Author(s) 2020.CC BY 4.0 License.We adopted a simplified version of the paleotopography prior to the AD 79 eruption starting from the 10-m resolution Digital Elevation Model of Tarquini et al. (2007) and from the reconstruction given in Cioni et al. (1999) and Santacroce et al. (2003) (Fig. 8).Particularly, the present Gran Cono edifice and part of the caldera morphology have been replaced with a flat area, and a simplified reconstruction of the southern part of the Mount Somma scarp has been inserted.However, https://doi.org/10.5194/se-2020-138Preprint.Discussion started: 9 October 2020 c Author(s) 2020.CC BY 4.0 License.

Figure 6 .
Figure 6.Mean thickness comparison between the simulations (dashed lines) and the actual deposit (solid line) of a) EU3pf and b) EU4b/c units.Different boxes concern different circular sectors.330

Figure 7 .
Figure 7.Comparison between simulations (dashed lines) assuming different initial volumetric fractions of solid particles (ε 0 ), and the actual deposit (solid line), of the a) EU3pf unit S and b) EU4b/c unit SE.In (b), the inset is a magnification of the thicknesses further than 9 km from the vent.

Figure 8 .
Figure 8. Inundation area of the simulations of the EU3pf (a-b) and EU4b/c (c-d) units.The dashed lines represent the theoretical isopachs (in m) of the simulated deposit.Vent location (red triangle), vent uncertainty area (red line) and SV caldera (orange dashed line) as in Tadini et al. (2017).MRLs as in Figure 3.The DEM used in the simulations and as a background derives from Tarquini et al. (2007) according to the modifications explained in section 4. 340

Figure 9 .
Figure 9. Volumetric content of ash/lapilli of model/deposit with distance to the source, of the units a) EU3pf N/S (left and 345 right respectively) and b) EU4b/c NW/SE (left and right respectively).

Figure 10 .
Figure 10.Comparison of volumetric grain sizes of the a) EU3pf and b) EU4b/c units.Different boxes concern different distances to the source.

Table 1 .
Volume of the EU3pf and EU4b/c units.

Table 3
summarizes the calculated MSDs for the studied units (in Φ), along with the corresponding density values (obtained interpolating those in Table2).

Table 3 .
MSD values and related densities for the different units studied.

Table 5 .
True Positive (TP), False Positive (FP) and True Negative (TN) instances of the simulations in Figure 8.