{'fraction_non_alphanumeric': 0.030299094552471425, 'fraction_numerical': 0.01705135817129286, 'mean_word_length': 4.793507470708374, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reconstructing eruptive source parameters from tephra deposit: a numerical approach for medium-sized explosive eruptions

Since the seventies, several reconstruction techniques have been proposed, and are currently used, to extrapolate and quantify eruptive parameters from sampled deposit datasets. Discrete numbers of tephra ground loadings or stratigraphic records are usually processed to estimate source eruptive values. Reconstruction techniques like Pyle, Power law and Weibull are adopted as standard to quantify the erupted mass (or volume) whereas Voronoi for reconstructing the granulometry. Reconstructed values can be affected by large uncertainty due to complexities occurring within the atmospheric dispersion and deposition of volcanic particles. Here we want to quantify the sensitivity of reconstruction techniques, and to quantify how much estimated values of mass and grain size differ from emitted and deposited ones. We adopted a numerical approach simulating with a dispersal code a mild explosive event occurring at Mt. Etna, with eruptive parameters similar to those estimated for eruptions occurred in the last decade. Then we created a synthetic deposit by integrating the mass on the ground computed by the model over the computational domain (>50000 km2). Multiple samplings of the simulated deposit are used for generating a large dataset of sampling tests afterwards processed with standard reconstruction techniques. Results are then compared and evaluated through a statistical analysis, based on 2000 sampling tests of 100 samplings points. On average, all the used techniques underestimate deposited and emitted mass. A similar analysis, carried on Voronoi results, shows that information on the total grain size distribution is strongly deteriorated. Here we present a new method allowing an estimate of the deficiency in deposited mass for each simulated class. Finally a sensitivity study on eruptive parameters is presented in order to generalize our results to a wider range of eruptive conditions.


Introduction
During explosive volcanic eruptions, solid particles of various dimensions, ranging from several cm to a few microns (tephra) are injected into the atmosphere. Particle are subsequently settled toward the ground under the action of atmospheric dynamics, aggregative processes, and gravity. Atmospheric transport process, depending on particle characteristics (size, density, sphericity) and altitude of release, can last from a few hours up to weeks. Area affected by ash fallout can extend thousands of square kilometers around the vent (Sparks et al., 1997). A primary objective of most volcanologists is to quantify and to measure the magnitude and the intensity of an eruption. Erupted mass is fundamental in order to deepen our comprehension of the mechanisms triggering a volcanic eruption, the column evolution, and the ash cloud dynamics. During a volcanic crisis several direct measurements can provided information about the ongoing activity. Remote sensing instruments are currently adopted by volcanologists to give estimates of column height (Rose, Kostinski, and Kelley, 1995; Petersen and Arason, 2010), duration of the event (Johnson, Aster, and Kyle, 2004), gas and particles exit velocity (Dubosclard et al., 2004), and plume composition (Rose, Bluth, and Ernst, 2000; Spinetti et al., 2013). Some of these methods provide estimates with good resolution and details on the measured quantity. Nevertheless, after an eruption, the deposit is still one of the main products to analyze in order to attempt a cumulative estimation of erupted mass and grain-size distribution. Historically, the integration of discrete samplings of deposit thickness values (or mass loading values) has been used to quantify the volume (or mass) of solid erupted material and, in this way, to assign a size and style to an eruption (Fisher, 1964;Walker, 1973). In past decades, several techniques have been proposed to optimize the integration of field data and to obtain an estimate of the emitted material (see for example Pyle (1989), and Fierstein and Nathenson (1992)). More recently, Bonadonna and Houghton (2005) and Bonadonna and Costa (2012) introduced other methodologies for estimating volume and total grain size distribution, currently adopted by volcanologists during their field studies. Furthermore, Burden and Chen (2013); Engwell, Sparks, and Aspinall (2013) used statistical methods to study the uncertainty associated with tephra thickness measurements and volume estimation . Indeed, when estimating the amount of erupted material and grain-size distribution from discrete sampling the information that we can obtain directly from the tephra deposit is just a part of the total emitted by the source. Very few attempts have been made to interpret and revise the ground information in order to better constrain the initial eruptive conditions: Connor and Connor (2006) used an inversion techniques and recently Gudmundsson et al. (2012); Stevenson et al. (2015) performed an integration of ground measurements and satellite observations. The importance of estimating eruptive source parameters (GSD and erupted mass) with a good accuracy is also due to the growing use of dispersal codes for ash hazards assessment purposes (Sparks et al., 1997;Textor et al., 2005;Folch, 2012; Fagents, Gregg, and Lopes, 2013). Indeed, numerical models for ash dispersion are generally initialized with source conditions inferred from field observations. Numerical codes can also help in constraining eruptive source conditions by comparing numerical results with measured data. In Burden, Phillips, and Hincks (2011) and , the estimation of column height as a function of depositional clast size has been inferred by the use of a numerical code. Klawonn, Wolfe, Frazer, and Houghton (2012), also inverted deposit data to better constrain eruptive parameters with a dispersal model. Quite recently,Woodhouse et al. (2013); Bursik et al. (2012) investigated wind effect on column dynamics and height revealing the importance of its consideration to avoid strong under-estimation of mass flow-rate. Indeed, a numerical simulation has the advantage to keep all main initial conditions under control, and to follow the evolution of the main macroscopic features of investigated process. Here we use a computational approach to explore the quality of information contained in a tephra ground deposit and to test the efficacy of reconstruction techniques in quantifying eruptive parameters, used as code input. We pursue this objective by simulating a fictitious eruptive event at Mt. Etna for which initial source conditions are defined a priori and calculating its deposit as described in the Fig. A4. Deposited masses for each grain size (DM i ) are computed and compared with the known corresponding masses emitted at the vent (EM i ). Furthermore, we randomly sample the deposit and apply several reconstruction techniques to obtain an estimate (DM * i ) of the deposited mass. This is repeated multiple times, in order to perform a statistical analysis on the results.  Figure 1: Schematic representation of the computational strategy adopted.
2 Adopted method

Numerical model
To create a synthetic deposit we adopt the dispersal code VOL-CALPUFF (Barsotti, A, and J, 2008) simulating volcanic ash transport and deposition in a transient and 3D atmosphere. The model couples an Eulerian description of the initial plume rise phase, where plume theory equations are solved (Bursik, 2001), with a Lagrangian description of the transport of material leaving the eruptive column. At each step particles are released from the column as a series of Gaussian packets (puffs) which are transported and diffused by the wind during their fall toward the ground due to gravity. Tracking puff movements within the 3D domain, the code computes at each time step the amounts of mass advected out of the domain, still suspended in the atmosphere and deposited on the ground. VOL-CALPUFF model has been tested and adopted to simulate volcanic ash dispersal and deposition at several volcanoes worldwide (Barsotti and Neri, 2008; Barsotti et al., 2011;Spinetti et al., 2013), confirming how this model is a well-established tool for simulating eruptive events of variable sizes and occurring over a large range of spatial scales. For the application shown in this paper, the VOL-CALPUFF model has been initialized with meteorological input data produced by the non-hydrostatic code LAMI (Doms and Schattler, 2002) with 7km horizontal resolution and 23 vertical levels and refined down to one km horizontally and one hour in time by the processor CALMET (Scire et al., 2000).

Recontruction techniques 2.2.1 Erupted mass
The determination of the amount of erupted material is historically made starting from field measurements collected over the domain at discrete points (Walker, 1973;Pyle, 1989;Fierstein and Nathenson, 1992). These data can be either measurement of deposit thickness for historical eruptions or of mass loading for more recent ones. Several methods have been proposed to estimate erupted volume from deposit data (Rose et al., 1973;Pyle, 1989;Fierstein and Nathenson, 1992;Sulpizio, 2005;Bonadonna and Costa, 2012;Burden and Chen, 2013). Estimation process consists on several steps, as suggested by the IAVCEI Commission on Tephra Hazard Modelling . First step is drawing a discrete number (N) of isopach (or isomass) contours, usually using an interpolation techniques. Second step is calculating the squared Table 1: Reconstruction methods presented with explicit formulas for mass loading (t) as a function of corresponding isomass areas (x), and residual cost functions. N is the number of isomass. For Weibull methods (referring to Weibull 1, 2 and 3), parameters have been chose, according with (Bonadonna and Costa, 2012), in the best initial range λ ∈ [.1, 1000]km, θ ∈ [.1, 5000]cm and n ∈ [0. 2,2].
roots of the obtained isopach areas (x i ) corresponding to t i thickness. The main idea is to express t i as a function of the squared roots x i of its corresponding isopachs area. This is done by a fitting procedure, where the function t(x) is chosen from a fixed family, described by one or more parameter, by minimizing a residual function. Finally the eruptive volume is calculated by integrating the function t(x) over the domain. Recently an alternative to these steps, avoiding drawing isomass maps, has been proposed in Burden and Chen (2013) using a statistical method on the sampling points measurements.
Within the first category of methods, in this work we test the following families of fitting functions (see also Table 1): 1) Exponential or Pyle's method (Pyle, 1989): This method is based on the assumption that the thickness, calculated on elliptical isopachs, follows an exponential decay with distance from the vent. In Fierstein and Nathenson (1992) and Bonadonna, Ernst, and Sparks R (1998) A modification of the method is proposed, mainly to account for the break-in-slope of some tephra deposits, by using various exponential segments. This method is now widely used even though it is very sensitive to the number and the choice of straight segments. Here we utilize three segments and an automatic procedure is defined to compute their bounds in order to minimize the residual function, thus reducing the arbitrariness of the choice.
2) Power law method (Bonadonna and Houghton, 2005): In their paper the authors suggest a power-law best fit of field data to integrate the total erupted volume. Associated volume can be calculated as: V = 2t pl x m+1 dx , but this integral, when evaluated over the interval [0, +∞] , is infinite. To avoid this problem such an interval is replaced by the arbitrary one [x 0 , x 1 ]. As already stated by Bonadonna and Costa (2013) Power Law results are strongly influenced by this choice. For this reason, in our statistical analysis average values of these distances are fixed at 1km and 300 km (accordingly with the eruption size) with a range of variability equal to ±10% .
3) Weibull method (Bonadonna and Costa, 2012; Bonadonna and Costa, 2013): Recently a new method based on Weibull functions has been proposed by Bonadonna and Costa. As main assumption the model describes the function xT (x) with a Weibull distribution and three main parameters (λ, θ, n). λ represents the characteristic decay length scale of deposit thinning, θ represents a thickness scale and n is a shape parameter (dimensionless). In this work, model outcomes are expressed as ground loadings, so we apply all the previous methods for the mass estimation by considering isomass contours instead of isopachs as frequently done in case of fresh deposit (Scollo, Del Carlo, and Coltelli, 2007;Andronico et al., 2009).

Granulometry
Total grain-size distribution (TGSD) of tephra-fall deposits is a crucial eruptive parameter for tephra-dispersal modelling and hazard mitigation plans Folch, 2012) and, several techniques have been developed to reconstruct it on the basis of ground information (Murrow, Rose, and Self, 1980; Bonadonna and Houghton, 2005). Usually, grain size are expressed in φ scale where φ = −log 2 (d/d 0 ) , d is the particle diameter and d 0 is a reference diameter of 1 mm. After Bonadonna and Houghton (2005), Voronoi technique has been adopted by the volcanological community as a benchmark for reconstructing TGSD. For this reason, it has been chosen in our work for a comparison with the exact emitted and deposited granulometry. Voronoi, or nearest-neighbor, technique is a constant piecewise interpolation method, built upon the Voronoi tessellation. Given a set of points called seeds, the plane, over which they lay, is partitioned into convex cells, each one consisting of all those points closer to that seed than to any others. Then, as proposed in Bonadonna and Houghton (2005), the grain size and mass load of each sample point is assigned to the corresponding cell and total grain-size distribution is obtained from the mass-weighted (where the mass is defined as the product between the cell area and the cell ground loading) average of all the sampled values over the whole deposit. The main requirement to apply Voronoi method is to fix the deposit extent, adding to the original data set also deposit zero values. This step is essential to prevent the external Voronoi cells from having unlimited area. Dependence of the reconstructed granulometry on the location of the zeros has been previously recognized and discussed in (Bonadonna and Houghton, 2005).

Application to Mt. Etna
Mt. Etna is the most active volcano in Europe, and each year it shows a wide variability in its eruptive activity. It has been extensively studied, and a wide literature exists regarding past activities, eruptive styles and products (Branca and P, 2005 (Andronico et al., 2009). At that time a quite persistent wind was blowing toward south-east with intensities up to 10 m/s at 3 km above the vent. Due to these meteorological conditions, the investigated domain, here fixed to 272x277 km 2 , is displaced south-east of the vent. All the main volcanological parameters, (e.g. mass eruption rate) are fixed a priori in order to simulate a mild eruption producing a 3-km high tilted column above the vent lasting for 6 hours. To this aim, the mass flow-rate has been fixed to 1 · 10 5 kg/s, with total emitted mass of 2.5 · 10 9 kg. Simulated total grain size distribution is in φ range [-5, 6] ( 32 mm to 16µm). Resulting column height and ground loading values are consistent with observation (Andronico, Cristaldi, and Scollo, 2008) and more general with those observed for recent eruptive activities at  Fig. 2A shows, for different size classes, the vertical distribution of cumulative mass released along the column. Mass lost along the column is calculated for each class as a function of settling velocities. Column height m agl Settling velocity is computed as a function of particle Reynolds number and depends on particle characteristics (dimension, density, sphericity) as well as atmospheric properties (air density and viscosity), according with the modification of   A sensitivity analysis has been carried out assuming a normal TGSD with parameters µ and σ ranging in the intervals [−3, 5] and [0.5, 6] respectively. Initial distribution has been then partitioned in 24 classes equally spaced in the φ scale, while the other input parameters have been kept fixed as described before. Fig. 3B presents the results of the sensitivity analysis, showing that column height is slightly sensitive to the size distribution. For this reason, without any loss in result generality, we assumed an initial granulometry with a uniform distribution. This choice does not affect results also because, assuming particles do not interact during the transport, for each class the percentage of deposited mass with respect to the emitted one is independent from the assumed initial granulometry (the independent deposits for each class are presented in the Auxiliary Material). As a consequence our results can be easily extended to different distributions of particles leaving the column.   Dispersal model produces a numerical result on 250×255 computational cells. At each cell center is associated a ground loading value (mass per unit area), obtained from mass integration of Gaussian packets deposited on the ground over the cell. This integration is performed by using a Taylor formula with a maximum numerical error fixed a priori (smaller than machine precision). In the following we will refer to this integration of Gaussian functions with the term "exact integration". In this way we calculate, for each simulated particle size, the cumulative amount of mass deposited over the computational domain. Fig. 3A shows the simulated tephra deposit calculated after 24 hours. This time frame is a reasonable time to allow the ash cloud leaving the computational domain and in real events to prevent further modification of the deposit by wind action. Only a small part of the deposit interests the inland while a large amount falls into the sea. Deposit cross-wind extension mainly depends on wind temporal stability and direction during the eruptive event. Due to wind variability on direction with altitude and with time the simulated eruption produces a wide deposit, which enlarges downwind the vent (Fig. 3A). In our case the column injects material up to 3 km above the vent and the particles falling close to the vent produce ground loading values larger than 10 2 kg/m 2 . Far away ground loading decreases but still shows values larger than 1 kg/m 2 at about 30 km downwind of the volcano (close to the main city of Catania and the International Airport "Fontanarossa").
To test reconstruction techniques the whole set of 250×255 integral loading values is sampled to create 2000 smaller data sets ("sampling tests") of 100 points over which different techniques are applied. For each sampling test, points are generated using a cylindrical coordinate system and the sampling procedure is performed using a uniform probability distribution along the supposed major dispersal inferred main direction axis and a normal distribution for the angular coordinate. Major axis is calculated based on output atmospheric concentration and wind direction as usually done in sampling field work for recent eruptions (Bonadonna and Houghton, 2005). Due to wind direction variability during the eruptive event, and therefore to the difficulty to find a major dispersal axis, a small random variability is considered on the angle in each sampling test. The area proximal to the vent is excluded by the sampling, but at least a point within 4 kilometers from the vent is considered in each test. We also sampled only inside the 10 −2 kg/m 2 isomass contour (corresponding to 0.01mm) represented by black points. This value is chosen as a threshold in agreement to studies on eruptions of similar size ( Two examples of sampling tests are shown in Fig. 3B and 3C where the Natural Neighbor (NN) interpolation method (Sibson, 1981) is applied to the sampled points to reconstruct the deposit all over the domain and 20 logarithmically spaced isomass contours are drawn.

RESULTS
After an eruption, tephra deposit is often the only product available to volcanologists to estimate physical quantities in order to characterize the explosive event. The eruption "size" and style is estimated by generally inferring the mass and the particle size distribution, emitted and deposited. In this section mass (or volume) and granulometry values, calculated using "exact integration" and classical reconstruction techniques are compared with the known values used as Scenario input. Several well-established reconstruction techniques have been tested over 2000 sampling tests.

Erupted and deposited mass
Here, as already stated, we first perform an "exact integration" of the ground loading produced by the numerical simulation on each cell of the computational domain. Avoiding intermediate steps as interpolation techniques and fitting functions we reduce all the subjective choices. In this way, the simulated deposit is analyzed to quantify the exact amount of emitted material falling in the considered domain of 250×255 square cells after 24 hours. In Fig. 4 cumulative mass and cumulative number of cells vs. ground loading are plotted (dots).Total emitted mass, equal to 2.5 · 10 9 kg, is reported on the left y-axis, and this reveals how the exact calculation of deposited mass is lower than the emitted one.
In particular the deposited mass corresponds to about 82% of the emitted one even though the calculation has been obtained considering ground loading values down to 10 −4 kg/m 2 . In these figures we also plotted, with black dashed lines, the 10%, 50% and 90% of the total deposited mass and the corresponding ground loading which is, respectively of 3.0 · 10 2 , 30 and 4 · 10 −2 kg/m 2 .  The dashed lines denotes the 10%, 50% and 90% of the total deposited mass and the corresponding ground loading. The black star on the y-axis corresponds to the emitted mass equal to 2.5 · 10 9 kg. In gray is expressed the cumulative pixels number calculated for a given mass loading.

Comparison with reconstruction techniques
In the bottom panel of Fig. 5 are plotted the corresponding values of the mass loading for each of the contours versus the square root of the area enclosed by the isomass for all the 2000 sampling tests (colored dots). Here the results of a statistical analysis are plotted. As a reference also the two selected tests of Figure 3 are reported with green and yellow circles. Dark gray line represents the median of the isomass squared area. Whereas the gray and light gray regions are respectively the 25 th − 75 th and the 5 th − 95 th percentile.
Reconstructed isomass area values (Fig. 5C) are underestimated on average for ground loadings larger than 10 −1 kg/m 2 , which is where most of the deposited mass In addition to all these techniques also a numerical integration (Num Int) of the reconstructed deposit using the NN is presented (orange). The 95 th −5 th intervals are also reported. As a reference the dashed line represents deposited mass whereas continuous black line the emitted one.
is concentrated (Fig. 4). On the contrary the sampling tests, on average, overestimate the areas for ground loading values smaller than 10 −1 kg/m 2 .
Finally, using the different functions presented in Table 1, we fit ground loading values as function of squared area values and the best fitting function is integrated to calculate the deposited mass, as shown in Fig. 5A-B. Black dashed line represents the mass loading vs squared area plot, obtained with the exact integration over the whole domain, as done for the plots in Fig. 4. In Fig. 5D the estimate total mass is presented in a synthetic way, expressed for each of the techniques by the average values (averaged over the 2000 tests) with the 5 th and 95 th percentiles and the standard deviations. For an immediate comparison the values of the emitted mass (solid line) and exact deposited mass (dashed line) are also reported. All the techniques except Weibull3, on average, underestimate both the deposited and emitted mass. Weibull3, on average, overestimates total mass up to 1 order of magnitude, and in general using the Weibull functions a large variability in the results is produced. Due to this high variability, Bonadonna and Costa (2013)

Deposited Granulometry
In analogy with the previous section, here we first compare the emitted granulometry EM i with the deposited (DM i ) obtained from an exact integration of the simulated deposit over the 24 grain size classes. Afterwards, results are compared with the TGSD obtained using the commonly adopted Voronoi tessellation method.
Results are shown in Fig. 6A, where deposited mass values obtained with the exact integration over four spatial domains (represented by different color bars) are reported in a bar plot as percentages of the emitted mass. Blue bars correspond to the deposit integrated over the whole computational domain. Other bars represent smaller domains, obtained excluding two squares centered on the vent with side length of 4 (Subdo-main_1) and 2 km (Subdomain_2) and they can be seen as representative of samplings performed excluding very-proximal regions around the vent. Deposited mass (DM i ) (Fig. 6A), down to φ ≤ 3, reaches 100% of the emitted one ( EM i ). For φ = 4 this percentage reduces to about 90% and abruptly down to 5% for φ = 5. Finest simulated class is almost absent in the computed deposit. Results obtained excluding the material deposited inside a square centered on the vent with a side smaller than 2 km (yellow bars) do not significantly differ from those obtained considering the entire computational domain (Andronico et al., 2014). On the other hand, on Subdomain_1 (cyan bars) all the classes larger than φ = 0 are strongly underestimated in the deposit, and DM tot is only 68% of EM tot . Largest simulated class loses about 80% of the initial mass. As expected, the deposited mass estimation strongly depends also on the minimum ground loading (hereafter referred as threshold) considered for the integration: the higher the threshold, the larger the underestimation of deposited mass. Fig. 6B presents for each simulated grain-size the deposited mass DM i when different thresholds are considered. Classes in the φ range interval [1,2] are underestimated of about 5-10% considering a 0.1 kg/m 2 threshold (cyan bars); these percentages rise to 70-90% when 10 kg/m 2 thresholds are considered (orange bars in Fig. 6A). We can notice that for our scenario a large amount of the emitted mass for classes φ > 4 is advected out the domain. In order to better highlight the depositional patterns for particles of different size, Fig. 7 shows the computed ground loading contours [g/m 2 ] at the end of the simulation for different classes: A) φ = 2 and B) φ = 4. For each class, the black isomass line encloses the area with mass loading larger than the 50% of the maximum ground loading over the domain (max cells (DM i (cell))/2). In this way we quantify for each class both the spreading of the deposit and the amount of mass close to the peak location. This information is reported in Fig. 7C, where markers size corresponds to the area enclosed by the two isomass curves. In addition, the mean distance from the vent is calculated for each class averaging the cells with ground loading larger than the 50% maximum (square) and as a reference inside the 10 −2 g/m 2 isomass (circle). If this last one is open the last closed isomass is shown with a dashed circle. Colors represent the cells corresponding to ground loading [ g/m 2 ] and the dimension of the symbols the corresponding area [ km 2 ]. Pink area on the y-axis of panel C denotes the distance at which ash particles reach the domain border. Smaller particles, due to the low settling velocity, are carried farther and diffuse more by the action of the wind. These two effects are respectively visible on the increase on the mean peak distance from the vent and on the increase of the corresponding cells number with the consequent reduction of the max ground loading for the finer classes.   As we have already shown the effect of the domain restriction is drastically reducing the amount of deposited mass for finer classes. Here this is well visible for classes smaller than φ = 4 where the peak is leaving the domain border. As consequence only the finer material falling off the margin of the column is found. The cutoff effect due to a limited domain is not easily overcoming, in fact the area interested by the fallout of particle finer then φ = 4 in the case of the described Scenario is larger than 10 6 km 2 . However the spatial distribution of small class is more homogeneous due to the larger diffusion, this is visible from the reduced distance between the last closed isomass (circles) and the peak (square). As direct consequence, for a smaller class, the error committed using the interpolation techniques for reconstruct DM i is small. Since column's height depends minimally on the emitted GSD, we can conclude that for each class the peak distance from the vent is minimally depending on the emitted percentage. Similar results have been found using analytical model for mono-disperse plume released from height altitude (Tirabassi and al., 2009). Because the peak distance, for fixed column height is not depending on the emitted mass but only on the atmospheric condition, allow us to easily estimate it. This would help to find the best region to sample, in order to better reconstruct the eruptive scenario.

Comparison with the Voronoi technique
We test here Voronoi technique described in section 2.2.2 for reconstructing TGSD. The main requirement to apply Voronoi method is to fix the deposit extent, adding to the original data set also deposit zero values. We remark that fixing zero value is a strong assumption because, as shown in Fig. 7, zeros are not the same for each class. For our purposes, here we fix the zeros of the simulated deposit to correspond to the 10 −2 kg/m 2 isomass contours (about 0.01mm). As shown in Fig. 6, finer particles estimation strongly depends on the thresholds defining the external boundary of the deposit, and therefore by the zero position.   For TestA_1, the TGSD presents a multi-modal distribution, which is neither in the emitted granulometry nor in those obtained with the exact integration over the whole deposit (Fig. 6), representing thus an artificial effect of the sampling, and/or the reconstruction technique. One possible explanation for this is related with the fact that different classes, having different settling velocities, will generate a peak in the deposit at different distances form the vent. Indeed, if a sampling point is close to a peak of a specific class, this maximum value of the deposit will be extended over the entire Voronoi cell. This could have the final effect of overestimating the reconstructed class mass and, consequently, to cause a peak in the GSD. For particle with φ ≥ 3 the peaks is located beyond the domain as we can see in Fig 7. Due to the diffusion smaller classes present a more uniform distribution on the ground and a smaller maximum load. For this reason also a linear method is well approximating the deposit mass.
Peak, associated with the samples and the reconstruction technique, is in the relative amount of mass and not in the absolute value, and thus could be also due to an underestimation of the other classes.
For the classes corresponding to the peak a larger gap is present between the 5 th and the 95 th percentiles (light-gray area). A large variability in the results obtained for the different sampling tests is also indicated by the large standard deviation. Conversely, for finer particles, a better accordance between values obtained with exact integration and values reconstructed with Voronoi technique is found. For these classes, also the gap between the 5 th percentile and the 95 th percentile and the standard deviation are small, and thus the reconstruction is more stable with the respect to the choice of sampling points. To quantify the error associated with the reconstructions obtained for each sampling test, a mean relative error averaged over the classes is defined as The average of the relative errors is larger than 0.8 for both scenarios and 95% of the sampling tests have a relative error larger than 0.35, denoting that for the most of the test a large error is associated with the TGSD reconstruction.

DISCUSSION
Analysis presented in the previous sections confirms that inferring quantitative data from the deposits can be very difficult and entails misleading results. Estimates values are far from representing a full picture of the initial eruptive condition at the vent. On one hand, the reconstruction of the total mass and the granulometry with an exact integration over the considered computational domain (an area larger than 50000 km 2 for the study presented here) shows that a large amount of the emitted mass is not found on the ground. On the other hand, the information present in the deposit is further degraded when measured by using a finite number of sampling values and then processed by reconstruction methods. In particular, informations on coarser particles are partially lost when very proximal area is excluded (Fig. 6A). In addition, due to the finite area over which field samplings are performed, here simulated by a finite computational domain, and the larger area over which the dispersal phenomenon takes place, the information on the finest material is only partially reconstructed (Fig. 6B). A combination of these two effects, here obtained considering a ground loading larger than 0.1 kg/m 2 (roughly corresponding to 0.1 mm) and excluding an area within about 2 km from the vent, causes an underestimation of the deposited mass up to about 50%. This reflects in an almost uni-modal TGSD, reconstructed from the deposit, where the tails are both depleted with respect to the initial ones. We remark that the results and the underestimations discussed so far are not an effect of the procedure, since the integration is exact, but only of the choice of the integration domain, although quite extended in this study. In this case, the two constraints on the integration domain adopted to figure out the TGSD, i.e. ground loading > 0.1 kg/m 2 and distance from the vent > 2 km, can reflect a reasonable minimum value sampled in the field and the objective difficulties in reaching areas very close to the vent. Despite the strong assumption of a uniform distribution at the vent, the deposited GSD looks similar to those reconstructed on the base of field studies for similar real eruptions at Mt. Etna ( , (Scollo, Del Carlo, and Coltelli, 2007). Results are obtained for a φ range spanning from -5 to 6, but the same investigations, considering finer particles, would produce similar outcomes. Deposited mass, obtained by an exact integration, always results in an underestimation of the emitted mass EM. This highlight the need for some correction factors to fill the gap between the estimated and the emitted eruptive parameters. In order to generalize our results to other scenarios we conduct a sensitivity study varying the exit velocity and the radius of the column over 2400 combinations. Fig. 9A shows the average column height calculated over six hours for different scenarios. Each scenario is the result of a run made using different initial velocities and radii for the eruptive column. Initial velocity ranges between 10 m/s and 160 m/s, and the radius between 5 m and 55 m, resulting in a mass flow rate between 6.03 · 10 7 and 1.15 · 10 11 kg/h. Yellow square represents the initial condition of the scenario described in the previous section, with v=30 m/s and r=18m. Color plot is obtained interpolating the results over 2400 simulations (black dots). Similarly, in Fig. 9B, a color plot shows the percentage of deposited mass over the emitted (DM i /EM i ) for φ = 4. If we look at the reference scenario discussed in the paper (yellow square) on Fig 9B, we see that the deposited mass over the considered domain for φ = 4 corresponds only to the 90% ± 5% of the emitted one (see also Fig.6). In Fig 9C the deposited mass (DM i ) for different class sizes in the interval [2,6] in the φ scale is expressed as a percentage of the emitted (EM i ). This plot allows defining a correction factors for the other classes, and more general for different eruptive scenarios. In this way we can relate the deposited TGSD with the emitted one. Before applying the method we need to be careful considering the uncertainty related to the reconstruction techniques applied to samplings dataset. This is highlighted in Figs. 5D and 8B where the erupted mass (EM tot ) and the GSD are obtained respectively with Pyle, Power law and Weibull methods and with the Voronoi technique. In particular, despite the large number of samples considered for each sampling test (100 points could be considered a very well sampled deposit), the erupted mass values obtained show a large variability and, for a significant percentage of tests, an overestimation (Fig. 5D). In general, the Pyle and Power law fitting methods seems to produce more stable results than the Weibull function, providing a smaller standard deviation within the 2000 sampling tests performed, and making more reliable the application of a correction factor. Also when Voronoi technique is used, TGSD presents fictitious peaks in the resultant distribution (Fig. 8B) and a large variability in the results are produced depending on the choice of the sampling points with regard to locations. As already discussed, the grain size distribution is expressed as fractions of the total mass and not as absolute values of deposited mass for each class. In this way the error of a single class affects all the others as individual mass fractions depend on total mass. In this work measurements are considered not affected by errors. Ultimately, measurements errors could have the effect of a further deterioration of the information present in the deposit, decreasing the margin of confidence in the results.
Depositional models tend to underestimate ground loading in the proximal region and consequently overestimate distal deposited mass. This is mostly due to processes neglected by the models, for example aggregation  and gravity currents derived by instability (Manzella et al., 2015). In both cases finer particles are settling down with a higher settling velocity and consequently the amount of finer particles deposited is larger. On the other side smaller particles aggregated on a larger one modify the drag coefficients of the aggregate, reducing its settling velocity. As a consequence larger particles involved in the aggregation can travel longer.

CONCLUSION
In order to better understand a deposit footprint and its thoroughness in quantifying eruptive source parameters like mass and GSD, a synthetic tephra deposit has been generated by using a dispersal code, and then analyzed. We performed this analysis by estimating the emitted and deposited values of mass and GSD by using an exact integration over the computational domain as well as using standard reconstruction techniques. These values have been then compared with the emitted mass fixed a priori as input for the numerical simulations. Exact integration showed an underestimation of the initial quantities mainly due to the finite extension of the domain, but still consistent with the physical processes occurring during the transport and the sedimentation, for example the longer distance and time-travel associated with fine particles.
Quantitative reconstruction of eruptive conditions based on deposit sampling solely can be affected by large uncertainty, and can address result in misleading interpretation. Even by using a large and well-distributed dataset of sampled points (100 points) over the modelled domain the results can be affected by large errors.
In this case, estimated mass was not representative of emitted one. Reconstruction techniques, performed over 2000 different sampling tests, showed a gap in the mass values up to an order of magnitude between the 5 th percentile and the 95 th . Large standard deviation values and large relative errors affected reconstructed grain size distribution obtained with Voronoi. The technique has been applied to 2000 sampling tests of 100 points, revealing the strong sensitivity of the results on the considered samples and often producing fictitious peaks in the distribution. In this case, an exact integration performed considering a threshold on the ground loading and excluding a proximal area (reflecting a reasonable minimum value sampled in the field and the objective difficulties in reaching areas very close to the vent), produced an almost unimodal GSD where the tails are both depleted.
In particular the emitted fraction of finer particles was underestimated, with consequently a bias on the associate assessment of risk, especially for airborne.
We point out that obtained outcomes are not the image of a single eruption. Results have been extended using a sensitivity study on eruptive parameters, generalizing the analysis conducted so far depending on the eruption scale. Finally, this work showed how, in support of field-based studies, numerical models can represent a useful tool to better understand the complexity of the deposition process and, consequently, to manage the associated uncertainty.
On one hand, exact integration performed over the simulated deposit, which resulted in an underestimation of emitted mass has been used to define some correction factors for estimation of source eruptive parameters. On the other hand, results obtained with the exact integration provided reference values to compare the different integration techniques and assess their reliability.
In order to extend and generalize presented results other tests have been performed and presented on the auxiliary material

Recommendations
In the following part we summarized article result in a best practice list. In order to limit arbitrary choices on estimation process, our recommendations are: • Collecting values along the major dispersion axis with a log spaced distance from the vent.
• Collecting values along minors dispersion axis.
• Collecting zeros values to constrain deposit limits.
• Using an interpolation technique as Natural Neighbor for each class (no need of assigning zeros values as Voronoi).
• Integrating class by class.
• Using correction factors to account for domain restrictions.  This auxiliary material contains:

ACKNOWLEDGMENTS
• A sensitivity study of eruptive initial parameters.
• An analysis of the uncertainty for the volume estimation techniques based on a Montecarlo method.
• A sensitivity study of a simple aggregation model

A.1 Sensitivity analysis
In order to generalize the results presented in the article, a sensitivity study of VOL-CALPUFF has been performed. Tab A1 resumes main volcanological input parameters used to initialize the different runs. In the following we will refer to the main example presented on the article with Scenario_A. Fig. A1 shows a detailed deposit map for Scenario_A, where deposit (DM i ) is plotted for different particle classes: A) φ = −5 ... E) φ = 6. Isomass contours are log-spaced from 10 −5 to 10 [kg/m 2 ]. Although smaller particles have a larger diffusion coefficient they are carried farther by the wind due to the low settling velocity. Therefore most of the finer particles (φ >4) are leaving the domain and only a few percent of the emitted mass (EM i ) is deposited on the ground over the considered domain. Finest simulated class (φ = 6), as shown in Fig A1n, is almost absent in the computed deposit. Moreover, we can notice that different classes present different main axis directions on the deposit. This segregation process is mostly due to a shear on the wind direction with the altitude. Smaller particles are more influenced by higher wind directions and consequently also their deposit trace. This suggests to be careful when inferring the main dispersal axis and consequently when choosing sampling locations. Wind measurements and reliable weather data are fundamental if we want to proper characterize the footprint left on the ground by an eruption.
In order to generalize the results presented in the main paper, we investigate here results dependency on the emitted grain size distribution (GSD). location. Results are shown for three different classes φ = 2 (A), φ = 4 (B) and φ = 6 (C). For each of these classes the variability of DM i /EM i is really small (3−4% respect to the mean value), and for φ > 4 this ratio is even smaller. Therefore the percentage of deposited mass with respect to the emitted one is almost independent from the emitted GSD, allowing to generalize the results obtained for Scenario_A also to GSDs different from the original one presented in the manuscript.
In addition to the analysis on the effect of the initial particles distribution, here we want to investigate the role of particle shape factor ψ and density ρ on the transport and deposition process. Several experimental studies have highlighted how an irregular shape of falling particles could affect their velocity toward the ground (Walker, 1971;Wilson, 1972;Riley, Rose, and Bluth, 2003;Dellino et al., 2005). As introduced in Riley, Rose, and Bluth (2003) the sphericity parameter ψ is defined as the ratio between the projected area A p and the square of the projected perimeter P p : In Fig A3 deposited mass DM i for different class sizes (φ = 2 (A), φ = 4 (B), φ = 6 (C)) is expressed as a percentage of the emitted mass EM i . Each blue dot represents a run made using as emitted GSD a uniform distribution but with different factor shape ψ ∈ [0.5, 1] and density ρ ∈ [1500, 2000] kg/m 3 . In Fig A3F is presented a schematic example of different particles with different sphericity values ψ. For spherical particles ψ = 1. Colored contour plot is obtained interpolating the results over 121 simulations. In Fig. A3D and A3E) are plotted respectively as input response surfaces the emitted mass EM tot and averaged column height. We can notice that for a fixed φ class the ratio DM i /EM i increases with increasing density and sphericity. In fact more irregular and lighter particles are transported farther. Besides, relative changes in the DM i /EM i increases with the decreasing size: 0.3% for φ = 2, 7% φ = 4 and 8% for φ = 6. As expected EM i (Fig A3D) increases with particles density and is not influenced by the sphericity, whereas denser particles produce lower column heights (see Fig A3E). Finally, less spherical are the particles, less rapid is their vertical displacement allowing a longer staying in the atmosphere.

A.2 Uncertainty analysis
In order to quantify the interval of confidence of the fitting parameters and consequently of the estimated volume, different techniques have been recently introduced ( Figure A4: Schematic representation of the bootstrap method. Here in a similar way of (Biass et al., 2013) we used a bootstrap Monte-Carlo method (Efron and Tibshirani, 1993) as presented in Fig A4. In order to create a large dataset, the first step is to randomly perturb the calculated isomass area. If errors in the isomass squared area calculation (y points) are available, then the random variation in each point is determined from its given error. Otherwise, the random variation is determined from average standard deviation of y points. A large number of random data sets are then produced and fitted, and the variance of fit results is used as the error for the fit parameters. Even if estimations of observational errors have been produced for mass loading measurement (Bernard, 2013), it not trivial to predetermine the error distribution on corresponding isomass area, especially for sparse sampling points dataset.
Bootstrap method has been applied to the volume reconstruction techniques from two selected tests (TestA_1, TestA_2) using 1000 "perturbed" y-data sets. Analysis results are presented in table A2. We can notice the following main facts: • a small residual norm (i.e. a good fitting) does not guarantee a small interval of confidence.
• estimate error and consequently the estimate confidence interval are not representative of the real error committed.
Moreover, comparing the three different Weibull methods, we observe that results strongly depend on the choice of the residual function, and generally the volume mean is smaller than the volume errors. Unfortunately, it is difficult to find a non-subjective criterion for choosing the appropriate residual function. Even if we can find a satisfactory fitting of the mass loading vs squared area, this does not ensure a good volume estimation when errors are also associated with area values. However the estimation is based on the assumption that the fitting model is valid, and does not account for error in the interpolation of the isomass contours from sparse data. Finally we have to stress that, to a lesser extent, the mass estimation methods proposed are also sensitive to the values and the number of isomasses.

A.3 Aggregation
As stated in several papers (Carey and Sigurdsson, 1982;; Taddeucci et al., 2011; Brown, Bonadonna, and Durant, 2012), particle aggregation can significantly affect the dispersal process and the tephra depositional pattern. For this reason we extended our results considering a simple aggregation process. As proposed by Textor et al. (2006) aggregation is supposed occurring within the margin of the column. To this aim, we adopt a model based on (Cornell, Carey, and Sigurdsson, 1983;Sulpizio et al., 2012), where 10 − 50% of emitted mass for φ > 2 particles aggregates as φ = 2 (250 µm ) particles. As a first approximation particle density is not modified and is equals to 2000 kg/m 3 . As stated by Brown, Bonadonna, and Durant (2012)), aggregates, due to the collision on the ground and during measurements are supposed to be found as completely break apart on the deposit.
In Fig.A5A cumulative mass is plotted as a function of ground loading. Total emitted mass (EM tot ), equal to 2.5 · 10 9 kg, is reported on the left y-axis. This reveals how DM tot corresponds to about 82% of EM tot , even though the calculation has been obtained considering ground loading values down to 10 −4 kg/m 2 over an extended domain. In these figures, for the simulations without aggregation, we also plotted, with black dashed lines, the 10%, 50% and 90% of the total deposited mass ( DM tot ) and the corresponding ground loading which is, respectively of 3.0·10 2 , 30 and 4·10 −2 kg/m 2 . Results are shown considering no aggregation (dots) or with 50% degree of aggregation (crosses). When 50% of aggregation is introduced in the model, a larger amount of mass is deposited within the domain, resulting in an increase of 10% for Scenario A . Fig. A5B shows cumulative pixels number calculated as a function of mass loading for simulation with and without aggregation. Three intervals can be identify in ground loading values, where the crosses (simulations with aggregation) are respectively below (interval I) (at distance), above (interval II) (intermediate) or coincident (interval III) (proximal) with the dots (simulations without aggregation). These intervals depend on the way aggregation is modelled (in this study on the percentages assumed to aggregate) and, above all, on the dimension of the aggregates and their fall velocities. Most of the particles involved in the aggregation process do not fall in proximal areas. Consequently, close to the vent, where the cumulative ground loading is larger than 10 kg/m 2 aggregation has a negligible effect (interval III). Below these values aggregation reveals its main effects. In particular, in the range [10 −1 , 10] kg/m 2 (interval II) a large portion of particles falls over a smaller area with higher loadings. At the same time, the area corresponding to low loading values, here smaller than 10 −1 kg/m 2 (interval I), is reduced. This is also shown in Fig. A6A and A6C, where the TGSD is plotted for different aggregation percentages (0% -black, 10% -blue and 50% -yellow).
In Fig. A6B it is also plotted, for each aggregating i-class (finer than φ =3), the mass deposited as aggregate (ADM i ) expressed as percentage of the deposited one (DM i ). As expected the amount of deposited fine particles increases proportionally with aggregation. Nevertheless, even when 10% and 50% of aggregation is considered, for φ = 5 and φ = 6 classes deposit only as aggregates ( ADM i /DM i =100% ) in the considered domain, This plot reveals how using these ratio as indicative of the degree of aggregation of sampled particles can be misleading. Indeed, the observation of a particular class in the deposit only as aggregate could be erroneously interpreted as evidence that 100% of the emitted mass for the i-class aggregated during the transport.
vii A A B Figure A5: A) Cumulative mass versus mass loading. Dashed lines denote the 10% , 50% and 90% of the total deposited mass and the corresponding ground loading. Black segment corresponds to the EM tot equal to 2.5 · 10 9 kg. B) Figures shows results without (dots) or with 50% degree of aggregation (crosses). We can identify three intervals in ground loading values where the crosses (simulations with aggregation) are respectively below (interval I) (at distance), above (interval II) (intermediate) or coincident (interval III) (proximal).