Inverting sediment bedforms for evaluating the hazard of dilute pyroclastic density currents in the field

Pyroclastic density currents are ground hugging gas-particle flows associated to explosive volcanic eruptions and moving down a volcano's slope, causing devastation and deaths. Because of the hostile nature they cannot be analyzed directly and most of their fluid dynamic behavior is reconstructed by the deposits left in the geological record, which frequently show peculiar structures such as ripples and dune bedforms. Here, a set of equations is simplified to link flow behavior to particle motion and deposition. This allows to construct a phase diagram by which impact parameters of dilute pyroclastic density currents, representing important factors of hazard, can be calculated by inverting bedforms wavelength and grain size, without the need of more complex models that require extensive work in the laboratory.

www.nature.com/scientificreports/ is gas-particle mixture density, contrasts the resistance of buildings to the flow 13 . Values lower than 1 kPa are not expected to cause damage to buildings 20 . The current starts breaking openings between 1 and 5 kPa, and at values higher than 5 kPa breaking of walls starts to occur. Particle volumetric concentration is a factor of hazard because ash particles in suspension in air are very harmful to breath 21,22 , even at temperatures lower than 200 °C and at concentration as low as 0.001 (that is typical of dilute PDCs), if exposure time is longer than a couple of minutes 21,23 . The sedimentation rate helps in defining hazard because it can be used for approximating flow www.nature.com/scientificreports/ duration (t) of dilute pyroclastic density currents 23,24 (see the method section), which lasts until sedimentation is fed from turbulent suspension. In this paper we use data from bedforms of dilute pyroclastic density currents for reconstructing the impact parameters and contribute to hazard assessment.
Bedforms of the types of dunes and ripples have been widely recognized in the deposits of dilute PDCs since the pioneering observations of Richards 25 , Moore 26 and Fisher and Waters 27 . Further details on other types of bedforms are nowadays emerging from observation of recent eruptions 28,29 . Differently from what has been done for fluvial and turbiditic currents, only very few attempts have been made to construct phase diagrams defining the stability fields of bedforms as a function of PDCs flow regimes. Only very recently Smith et al. 12 have proposed a phase diagram for highly concentrated volcanic granular currents. Douillet et al. 8 carried wind tunnel experiments to relate the grain size of different types of volcanic particles to their threshold for motion at various slopes. Dellino et al. 30 proposed a phase diagram ( Fig. 1) in which volcanic deposits are classified according to the sedimentation rate and bedload transportation rate (Q b ), which are modelled after large-scale experiments 13,31 . This agrees with the approach used in the field of sedimentary currents, where it is widely recognized that the proportion of bedload to suspended load and the sediment size are the major controlling factors on bedforms formation 5 . The experimentally validated model equations of the sedimentation and bedload transportation rates, on which the phase diagram of Fig. 1 is constructed, are defined and discussed in the method section. The lower right portion of the diagram represents the field of massive deposits, which are not considered further in this study. The upper portion of the diagram represents the field of stratified deposits with ripple and dune bedforms, which are related to highly expanded, fast-moving, dilute and turbulent PDCs, which are the focus of this paper. In the original version of the diagram, a link between the wavelength of bedforms and the ratio between the sedimentation and bedload transportation rates emerged, with ripples having a ratio larger than 0.05 and dunes smaller than 0.05. However, a well-defined correlation could not be found because of the paucity of data.
In this paper, we further populate the diagram in the portion of dilute PDCs by adding 88 points relative to various eruptions of Vesuvius, Campi Flegrei and Vulcano in Italy. With this addition, the new dataset consists of 98 deposits (Fig. 1) and covers a wide span of the sedimentation vs bedload transportation rates space, allowing a more detailed analysis of bedforms characteristics in terms of the flow parameters.
In the enlarged dataset the bedform wavelength ranges from ripples ( Fig. 2a), starting at 10 cm, to dunes (Fig. 2b), up to 250 cm.
We never found antidunes, in fact their interpretation has always been questioned in volcanic deposits 7,29,32 . The software PYFLOW 2.0 by Dioguardi and Mele 33 has been used to plot data in Fig. 1. It was implemented here so to obtain both the impact parameters of the current and also the sedimentation and bedload transportation rates, as defined by (11), (12) and (13) of the method section. The software employs sediment data that result from time-consuming laboratory analyses, which involve technologies and calculation resources not available to all scientists. The aim of this paper is to rearrange and simplify the dataset in order to construct a  30 , which are also included. The legend of volcanoes from which deposits were analysed is inserted. Modified after Dellino et al. 30  www.nature.com/scientificreports/ phase diagram where by means of deposits wavelength and grain size the impact parameters of dilute PDCs can be reconstructed directly in the field. Four fitting laws (Fig. 3), showing good correlation coefficients, were obtained by means of a regression analysis of the variables appearing in the equations of the sedimentation and bedload transportation rates. They represent relationships between four main fluid dynamic parameters of the current (namely shear velocity, u * , particle volumetric concentration, sedimentation rate and bedload transportation rate) with bedform median grain size and wavelength. In the fitting laws many extra terms of the original equations do not appear anymore, thus allowing a reduction of the complexity of the formulas, and still evidencing the main characteristics of deposit formation process in terms of the current fluid dynamics.
A relationship between the bedload transportation rate and the product u * 2 C is evidenced (Fig. 3a). Given that C is directly proportional to the density of the gas-particle mixture (ρ mix ) and ρ mix u * 2 is the turbulent shear stress of the current 34 (τ), it implies that the bedload transportation rate is proportional to the shear stress, which confirms findings of sedimentary currents 29 . A relationship between D 0.5u * and u * 2 C (Fig. 3b) implies that  www.nature.com/scientificreports/ shear stress is proportional to bedforms median grain size, again in agreement with sedimentary deposits 30 . A relationship can be evidenced between the product of u * 0.4 C 0.624 and the sedimentation rate (Fig. 3c). Since the exponents of C and u * are both lower than 1, while in the fitting with the bedload transportation rate they are 1 and 2, respectively ( Fig. 3a), it means that with an increase of C and u * the difference between the sedimentation rate and bedload transportation rate increases, and their ratio decreases. This justifies the decrease of the ratio as the bedform wavelength increases, as shown by the fitting of Fig. 3d. It gives a quantitative definition of the relationship between wavelength and the ratio between sedimentation and bedload transportation rates, which could not be precisely defined in the original diagram of Dellino et al. 30 . This new fitting has been obtained by plotting the wavelength, as measured in the field, of 32 deposits characterized by well exposed bedforms, against the value of the ratio of sedimentation and bedload transportation rates of the same deposits, as calculated by PYFLOW2. The continuity in the trend of Fig. 3d is at odd since ripples and dunes are not supposed to represent a continuum, being them separate by a hydrodynamic discontinuity 34 . Indeed, small ripples do not interfere with the upper current surface, while dunes, being larger, interfere with it. This discontinuity does not appear in the diagram of Fig. 3d, likely because a true interface between the current and the surrounding atmosphere does not exist in PDCs, at least in well-developed dilute currents, which are instead characterized by a very gradual passage between the two 6 (see the method section for our model of dilute PDCs). By means of the fitting laws, the following system of equations is formed: It can be solved numerically, once grain size (D) and wavelength (λ) are obtained in the field, and the current shear velocity ( u * ), particle concentration (C), sedimentation rate (S rw ) and bedload transportation rate (Q b ), are obtained. The dynamic pressure, which is also relevant for the hazard of dilute PDCs, is subsequently calculated according to (1) and (2) by means of the particle volumetric concentration and shear velocity. The former is used for the calculation of flow velocity (V) as defined by the law of the wall of a turbulent boundary layer 34 which is the physical model of PDCs that we employ (see the method section), where V(y) is the velocity profile of the stratified flow 34 , and k s is the substrate roughness.
When comparing results obtained by (3) with those resulting from PYFLOW 2, the average absolute error of shear velocity is 28% and that of particle volumetric concentration is 30%, meaning that a good approximation of the two parameters can be achieved by means of the simplified formulas. In the case of sedimentation rate, the error is larger, about 45%, but still it is important to discuss its role as it allows to approximate the calculation of flow duration, t 23 , which is a main factor of hazard (see the method section).
As to show the relationship of bedforms' median size and wavelength with current fluid dynamics and impact parameters, the system (3) was used to construct the phase diagrams of Fig. 4, where focus is put on bedforms of wavelength between 10 and 300 cm, although bedforms with larger wavelength can be found in the geologic record of volcanic deposits. The latter scenario is out of the range of applicability of our model and only bedforms that develop on an almost flat surface are considered here. Much larger bedforms typically develop as an interplay between the current's flow dynamics and the ground morphology elements 29 . The range of median size (D) of Fig. 4 was set to typical grain-size values for laminated pyroclastic bedforms (between 4 and − 2 phi i.e. 0.0064 mm and 4 mm respectively) since larger sizes (coarse lapilli and bombs), which typically form lenticular beds, are thought to represent highly concentrated traction-carpets at the base of dilute PDCs 29,35,36 , to which our model does not apply.
The current velocity, V, of Fig. 4a, represents the average value of the lowest 1000 cm of the current and was obtained by integrating (4) over flow height, with k s = 10 cm. We chose this depth-averaging height of the current because, in dilute PDCs, the portion responsible for the dynamic impact is the lowermost one (the shear flow) and 1000 cm represent a reasonable estimate of an average building height 37 . V increases at increasing λ but with different trends and rates depending on particle size D. In the diagram, velocity ranges from about 10 m/s to about 130 m/s. The trends are significantly different for the finer grain sizes (1 to 4 phi) at the smaller wavelengths (up to 50 cm), which can be interpreted as the smaller the wavelength, the smaller are the sedimentation and bedload transportation rates and shear stress, therefore the higher the velocity required to develop bedforms with the finer grain sizes. This is because for fine ash, due to the very low Reynolds number of shear (R e* ), the initiation of motion at the bedload occurs at a very high critical Shield's number (θ t ) 24,38 . Naturally, larger clasts have a higher threshold for motion, and thus involve a steeper shear stress gradient in the basal flow. This configuration increases the turbulence level, and thus the sediment discharge, allowing bedforms of longer wavelength to develop. The volumetric concentration, (C), ranges from less than 0.001 to about 0.017 (Fig. 4b). It decreases as wavelength increases, and it does so for all grain sizes, although at a rate that decreases at decreasing grain size D, because a higher concentration favors a higher sedimentation rate (see (11) of the method section) and a larger ratio between sedimentation and bedload transportation rates, hence a smaller wavelength. The change in trend www.nature.com/scientificreports/ with decreasing grain size can be explained by the fact that the finer are is particles, the lower the concentration required to develop bedforms with small wavelengths. Current density, which was calculated by means of (2) and fixing ρ s = 2000 kg/m 3 and ρ f = 0.9 kg/m 3 (which is reasonable if the fluid, made up of volcanic gas plus entrained cold atmosphere, is at about 200 °C 23 ) (Fig. 4c) follows the trend of concentration, and varies from less than 2 kg/ m 3 to about 35 kg/m 3 . The trend of dynamic pressure (P dyn ) (Fig. 4d) varies from less than 1 kPa with smaller wavelengths and finer grain sizes, which is a value that does not cause severe damages to buildings 10,20 , to almost 30 kPa with larger wavelengths and coarser grain sizes, a value that can destroy even the more resistant, modern buildings of reinforced concrete 10,20 . The sedimentation rate (Fig. 4e) increases as grain size coarsens, meaning that with finer sizes flow duration is longer because smaller particles result in a smaller settling velocity. As far as the wavelength is concerned, for the finest sizes, the sedimentation rate increases at increasing wavelength, meaning a decrease of flow duration with longer bedforms. With the coarsest sizes, instead, the sedimentation rate decreases as wavelength increases, meaning a longer flow duration with longer bedforms. The ranges of wavelength and grain size of Fig. 4 replicate the ranges of our dataset of dilute PDCs deposits, and result in impact parameters that span from currents that do not impact severely on structures, to values of devastating effects. Such a range well represents the situation of large-scale dilute PDCs whose strength decreases along runout 23 , and change from totally destructive flows around the volcano to residual currents that in the distal outreach do not possess a high strength but can still be rich in ash, as is the case of the dilute PDCs of the AD79 Vesuvius eruption at Pompeii 23 . Such fine material can be dangerous to breath even at concentrations lower than 0.001 39 , if flow duration lasts more than a couple of minutes.
As an operative help, to complement the diagrams of Fig. 4, the system (3) was solved at discrete intervals of grain size and wavelength, to construct a table where the stability fields of dynamic pressure, particle concentration and sedimentation rate, are represented inside a grid ( Table 2). The values are averaged among the four neighboring grid points and the uncertainty is expressed in terms of one standard deviation. Dynamic pressure  . 4) where the values of u * , C and S r are set at half phi intervals of D in relation to λ. By means of these data, and specifying in (1) and (4) the value of ks, ρ s , and H at which to integrate V, more precise data of the impact parameters can be obtained.
With our diagrams and tables at hand it is thus possible to invert bedforms of past eruptions, and reconstruct the impact parameters of dilute pyroclastic density current, contributing to hazard assessment. Various approximations are introduced in our model as, for example, a fixed average density of particles, which in fact can be highly variable among the components of the grain-size mixture, or the median value as representing a single characteristic value of the grain-size distribution, which in some cases can be quite complex. Due to these approximations, we recommend to carefully consider the uncertainty ranges of Table 2 (and of the supplementary tables) when judging the values of impact parameters. Furthermore, it is true that bedforms are not always well exposed in their complete longitudinal profile, because of truncations due to erosion. Sometimes they are also difficult to measure precisely, because a direct access to the deposit is hard. However, in the case that these can be accessed, well-preserved bedforms are a common signature of dilute PDC deposition. In such case our model can be used for quantifying the impact parameters of dilute pyroclastic density currents and contribute to the hazard assessment of active volcanoes.

Method
The reconstruction of the impact parameters of PDCs is based on a flow dynamics model that starts with the assumption that the turbulent current is velocity and density stratified 19,37 . In the stratified multiphase gas-particle current, the basal part is a shear flow that moves attached to the ground and has a density higher than atmosphere (Fig. 5). The upper part is buoyant, because particle concentration decreases with height down to a value that, combined with the effect of gas temperature, makes the mixture density lower than the surrounding atmosphere.
In a PDC, particles are mainly transported by turbulent suspension and sedimentation is controlled by a balance between flow shear velocity u * , which is controlled by fluid turbulence and favors suspension, and particle settling velocity, w t = (4gD(ρ s -ρ mix )/3C d ρ mix ) 0.5 , which favors sedimentation, where C d is drag coefficient. During sedimentation, it is assumed that particles of different composition, i.e. crystals and glass, settle at the same aerodynamic conditions, e.g., with the same terminal fall velocity 37 . Therefore, by equating the settling velocity of the glass and crystal components in the deposit, and assuming that sedimentation starts when P n = 2.5 40 , hence when w t = u * , flow shear velocity and density ρ sf of the shear flow can be calculated after D, ρ s and C d are measured Table 2. Phase diagram in which the stability fields of the impact parameters P dyn , C and S rw , are expressed as a function of λ and D of bedforms. The values inside the grid represent the average between the four neighboring grid points and the uncertainty is expressed as the standard deviation. ks = 10 cm, ρ s = 2000 kg/m 3 . www.nature.com/scientificreports/ in the laboratory 37,41 . These are the main input data in the PYFLOW_2.0 code 33 , which allows reconstructing the current parameters. The code is based on a model that assumes PDCs behave as turbulent boundary layer shear flows moving over a rough surface 37 , which velocity profile is given by (9). The model has been validated by experiments 13 and already applied to other eruptions 42 . Here it is summarized as adapted from Dellino et al. 23 .

Wavelength (m) Grain size (ϕ)
The maximum volumetric concentration of particles that can be transported in turbulent suspension, i.e. the maximum current capacity, is a function of the Rouse number of the particulate mixture taken in suspension. The profile of volumetric concentration over current height is regulated by the Rouse model 40 .
where C 0 is the particle volumetric concentration at a reference height y 0 and H is the total current thickness. Assuming steady sedimentation, H is obtained by the ratio H dep /C sf where H dep is deposit thickness and C sf is the depth-averaged concentration in the basal shear flow, which can be calculated by ρ sf = ρ s C sf + ρ f (1-C sf ), when ρ sf and ρ f are known.
The shear-flow height and density are obtained by solving the system of (6) and (7), which is valid for a turbulent current where τ is the shear-driving stress of the flow moving down an inclined slope of angle α.
The density profile, which is a function of concentration, particle density and gas density, is: The gas density and Rouse number are obtained by solving numerically the following system: Equation (9) states that atmospheric density, ρ a , is reached at the top of the shear flow, H sf , and Eq. (10) states that the average density of the shear flow, ρ sf refers to the part of the flow that goes from the reference level, y 0 , to the shear flow top height, H sf .
By combining the velocity and density profiles, the dynamic pressure profile is finally obtained. The profiles of the flow parameters are expressed in terms of a probability density function that depends on the variance of particle characteristics. www.nature.com/scientificreports/ In dilute PDCs, sedimentation occurs by continuous aggradation during the passage of the current, and the total time of aggradation is a proxy of flow duration, t 23 , which is equal to deposit thickness, H dep , divided by the aggradation rate A r . The latter is equal to S rw divided by one meter, which is the reference width of the sedimentation rate per unit width, (see Dellino et al. 30 ). Therefore, flow duration, which approximates the time in which harmful concentrations of ash are suspended in the current to which a human being can be exposed, can be calculated by means of S rw .
The software has been here implemented, therefore, as to obtain also the sedimentation rate (S r ) as defined by the experimental model of Dellino et al. 24 , as: with the subscript i referring to the ith particle-size class and n being the number of size classes of the grain-size distribution of the sediment, where P ni = w ti /ku * is the Rouse number of the i th size fraction of the solid material suspended in the current, with k the Von Karman constant (0.4) and w ti the terminal velocity of the i th size fraction. P n * = P navg /P nsusp is the normalized Rouse number of the current, i.e. the ratio between the average Rouse number of the solid material in the current and the Rouse number at maximum suspension capacity. φ i , ρ s i and P ni are the weight fraction, the density and the Rouse number of the i th grain-size fraction, respectively. The sedimentation rate was transformed in the sedimentation rate per unit width, S rw in order to make it comparable with Q b dimension 30 .
The software has been also implemented with the calculation of the bedload transportation rate (Q b ) by the use of the experimental formulation of Dellino et al. 30 : where q bi is the volumetric bedload transport rate of the i th size fraction per unit width of the flow, and ξ = τ/τ ri is the normalized shear stress, where τ ri is the minimum shear stress needed to move the i th size fraction at bedload.