A GIS-based application for volume estimation and spatial distribution analysis of tephra fallout: a case study of the 122 BC Etna eruption

In this study, we demonstrate the use of a geo-spatial information system for volume estimation of fallout deposits and for identification of syneruptive and post-eruptive depositional mechanisms. For the first time, we present thickness distribution maps, isopachs maps, and fallout deposit volumes for the single stratigraphic units (A-F) of the 122 BC Plinian eruption of Etna, the most powerful eruption of this volcano in historical times. Thickness data collected during the field survey were organized into a geo-referenced database, and several interpolation algorithms were used to calculate the volumes of the six fallout layers of eruption (units A-F). The results are compared with those obtained using the Pyle method, which bases volume calculations on the exponential thickness-decay law of the deposits. Differences in the two methods are analyzed through applying two-dimensional (2D) and 3D geo-statistical analysis to thickness data, with an 'ideal' fallout deposit used as reference. Our approach allowed both identification of stratigraphic sections where the deposits were affected by secondary erosional or accumulation phenomena, and assessment of whether the secondary processes were caused by local morphologic conditions or by variations in eruptive dynamics (e.g., rotation of dispersal axis direction).


Introduction
To define mass eruption rates and volcanic hazards of pyroclastic fallout, it is important to reconstruct the distribution of tephra fallout deposits and calculate their volumes. Assessment of tephra volumes is complex and it depends on several variables. First, field data are scarce and of poor accuracy due to erosion or the inaccessibility of the deposit (mainly for prehistoric eruptions). Secondly, the methods used to calculate volumes introduce different assumptions, which can lead to 50% or greater errors in the estimated volumes [Fierstein and Nathenson 1992].
Several methods have been proposed to attenuate these errors [e.g., Rose et al. 1973, Walker 1980, Froggatt 1982, Pyle 1989, Fierstein and Nathenson 1992, Legros 2000, Sulpizio 2005, Bonadonna and Costa 2012], but each has its limitations. Methods based on the exponential decay of de-posit thickness with distance from the vent [e.g., Pyle 1989, Fierstein andNathenson 1992] have been widely used, even though they sometimes underestimate the total volume of tephra, and mostly in the case of eruptions that produce a significant amount of fine particles. This is because many of the best preserved historical fall deposits do not show simple exponential thinning, but rather two or more straightline segments on a semi-log plot of thickness versus area 1/2 [Bonadonna and Houghton 2005]. The method introduced by Walker [1980] takes into account the concentration of loose crystals in the deposit, and it is dependent on an accurate knowledge of the phenocryst/ glass ratio in both the pumice and the deposits [Houghton et al. 2000]; this method can also only be applied to deposits that consist of porphyritic juvenile fragments. All volume estimation methods introduce some errors, which suggests that a precision of two or three significant digits is generally meaningful [Sulpizio 2005].
GIS-based methodologies have been used in volcanology mainly for mapping areas at risk of dangerous syn-eruptive and post-eruptive phenomena, such as lava flow, tephra fallout, pyroclastic currents, lahars and floods [Pareschi et al. 2000, Bisson et al. 2007, Toyos et al. 2007, Bisson et al. 2010, Vicari et al. 2011a, Vicari et al. 2011b, Ganci et al. 2012]. These techniques have been little used for volume estimation of volcanic deposits on the basis of field data. In rare cases, GIS methods have been used for volume calculations of pyroclastic flow, lava flow and lahars only [e.g., Isaia et al. 2004, Behncke et al. 2005, Neri et al. 2008, Mũnoz-Salinas et al. 2009]. To date, no studies have used GIS to investigate the spatial distributions of measured thicknesses of pyroclastic fallout deposits, nor to calculate their volumes. The only studies that have dealt with this issue are based on numerical dispersion models. In particular, the study of Connor and Connor [2006] proposed inversion modeling to search for an optimal set of parameters that can explain some physical observations; e.g., the mass of erupted material.
The present study proposes that for each of the fallout units of the eruptive sequence of the 122 BC Plinian eruption of Etna, comparisons are made between volume calculations based on spatial interpolations assisted by GIS, and volumes obtained through the Pyle method [Pyle 1989]. The differences in the volume calculations obtained using these two methods are discussed by means of geo-statistical analyses, which allow the identification and better interpretation of syn-eruptive and post-eruptive depositional mechanisms for each fallout unit. Both spatial interpolations and geo-statistical analysis tools are available in the GIS software platform used in this study (ESRI®, ArcGis 9.2 software).

Background on the 122 BC Etna eruption
The 122 BC Plinian eruption was the most powerful eruption of Etna in historical times, and it led to the formation of a thick pyroclastic sequence that is well preserved on the south-eastern flank of the volcano. Historical chronicles report that the fallout of lapilli over the ancient town of Catania ( Figure 1) caused fires and the collapse of roofs, and hid the sun behind a cloud of ash for some days. The damage was so severe that the inhabitants were exempted from paying taxes to Rome for 10 years (Coltelli et al. [1998], and references therein).
The 122 BC eruption sequence included seven eruptive units ( Figure 2) that consisted of pyroclastic fall (A-F) and flow deposits (G); their lithological features have been described previously [Coltelli et al. 1998, Houghton et al. 2004, Sable et al. 2006]. In brief, unit A is a thin layer of well-sorted, coarse, black ash fallout. Units C and E consist of well-sorted scoria and subordinate wall-rock lithic lapilli, and these are interpreted as Plinian fallout deposits. Units B, D and F are bedded to massive ash-fall tuffs that were formed during the phreatomagmatic phases of activity. Lastly, unit G is a massive ash flow a few meters in thickness; however, as it crops out over a small area, we were not able to calculate its volume.  The only volume data for the 122 BC eruption were published by Coltelli et al. [1998]. According to these authors, the volume of units C and E is 0.285 km 3 , which they calculated as the cumulative thickness of the two layers from the ln(thickness) versus the square root of the isopach area plot [Pyle 1989]. In contrast, the estimated volume of units D and F is 0.1 ±0.02 km 3 . The total volume of the deposits was estimated to be about 0.4 km 3 . No volume estimation of units A and B have been reported previously.

Isopach maps and volume estimation using the Pyle method
For each volcanic unit (A-F), the spatial distribution of the measured thickness and the hand-drawn isopach maps are shown in Figure 3. Hand-drawn isopach maps are built by joining data points with the same, or very similar, thickness. In this way we draw curves that move away from the vent that represent the decreasing trend of the fallout thickness. On the basis of these new isopach maps, we estimated the volume of each fall-out unit (A-F) using the method based on the exponential decay of deposit thickness with distance from the vent [Pyle 1989]. These new maps were obtained by introducing new data derived from a recent field study ( Figure  3, red points), with respect to the previous study [Coltelli et al. 1998]. In particular, we discovered deposits that belong to unit C in some cores from the Augusta area ( Figure 3c); this indicated that the unit was 3 cm to 4 cm thick at 70 km from the vent [De Martini et al. 2010, Smedile et al. 2011. As far as the other units are concerned, we found new stratigraphic sections distributed on the volcano flanks and located at a maximum distance of 20 km from the vent (Figure 3d, e).

Volume estimation performed in the GIS environment
We first organized the field data using ArcGIS 9.2 software (ESRI® platform). The data were geo-coded according to the World Geodetic System 84 (WGS84) Universal Transverse Mercator (UTM) Zone 33 reference coordinate system. The thickness distribution of each unit was represented by a vector point layer that localized the exact position (x, y) of the analyzed stratigraphic sections, and defined the thickness value of each unit. There were 91 sections and a total of 252 sampled thickness points (40, 28, 91, 39, 34 and 20, for units A, B, C, D, E and F, respectively; see Figure 3). The number of sections did not always coincide with the number of measured thicknesses for each unit, because in some sections, one or more units might have been eroded or not deposited. For each studied section, we therefore checked whether the volcanic unit was overlain by the subsequent unit; if so, the section was labeled as complete, and if not, it was labeled as incomplete.
Starting from the distribution of the thickness deposit data points and from the hand drawn isopachs, we applied several types of spatial interpolation techniques to create 3D surfaces that reproduced the variations in thickness of each volcanic unit. In particular, we used a deterministic approach that applied four different algorithms that are available in the Arc tool box of ArcGIS 9.2: TIN (triangular irregular network), NNI (natural neighbor interpolation), IDW (inverse distance weighted), and SPLINE. We first applied a TIN interpolation. A TIN is a vector surface model that consists of a network of triangles that are formed by connecting nodes (input data) according to the Delaunay criterion [Lee and Schachter 1980]. In this model, each point on the edge of or inside the triangle has a specific value that is calculated through linear interpolation.
We then applied NNI, which differs from TIN interpolation because it builds a raster surface model. To calculate the thickness values in unmeasured points, this local interpolation finds the closest subset of input data points surrounding an unmeasured point and applies weights to these that are proportional to the overlapping areas between the Voronoi polygons of the data subset and the Voronoi polygon created around the unmeasured point [Okabe et al. 2000]. The surface passes through the input samples and is smoothed everywhere, except at the location of the input samples. The NNI is also known as the Sibson or 'area-stealing' interpolation [Sibson 1981].
The last two interpolations, the IDW and SPLINE interpolations, are quite similar because they assign values to unmeasured points according to the measured surrounding values and to specified mathematical formulae that determine the smoothness of the resulting raster surface. The IDW determines the values using a linearly weighted combination of a set of sample points, where the weight is a function of the inverse distance: where Vp j is the estimated thickness at location J, V i is the measured sample value in each point i, and d i is the distance between Vp j and V i . The SPLINE interpolation calculates the values with a mathematical function that minimizes the overall surface curvature, which produces a smooth surface that passes exactly through the input points described by the following equation: where N is the number of used points, m j are coefficients that are determined through the solution of a system of linear equations, r j is the distance from the point (x,y) to j, and T(x,y) and R(r) are defined according to the selected mode.
These last three interpolations (NNI, IDW, SPLINE) produce a raster model that is represented by a grid (a matrix of rows and columns) in which the thickness data are attributed to each cell, instead of the single points in the TIN model. The spatial resolution of these raster models is 150 m. The four spatial interpolations (SInt) were applied in three different modes, using: (i) the available thickness points and the most distal isopach; (ii) isopachs alone; and (iii) the available thickness points and isopachs. The isopachs were clearly the same as those used for the volume calculation according to the Pyle method. Due to the lack of measured deposits in the most proximal areas, an additional value was introduced at 1 km downwind from the vent. This value was extrapolated from the mathematical curve that best fit the data points in the plot representing deposit thickness along the dispersal axis versus distance from the vent. Figure 4 reports the best fit of the points, which in some cases is represented by a power law curve and in others by a straight line (Table 1, see R 2 ).
The volumes estimated using the four spatial interpolations in the three modes described above are reported in Table 2, where they are compared with the values obtained using the Pyle method. The volumes of each unit calculated using TIN, NNI, IDW and SPLINE interpolations have the same order of magnitude as, and are always less than, those determined using the Pyle method. The differences be-tween the Pyle and SInt volumes, expressed as * 100, where vi is the volume of each volcanic unit and SInt j is the spatial interpolation considered), are reported in Table 3. In general, the lowest D values for units B and E were obtained when processing only the thickness points (mode (i)), whereas the lowest D for units A, C, D and F were obtained when processing isopachs alone (mode (ii)). Note that GIS-BASED TEPHRA VOLUME ESTIMATION  To explain the reasons for this significant discrepancy, we combined volcanological and stratigraphical information on the fallout deposits studied with a detailed geo-exploration of the data, as a GIS-based approach that analyzed the spatial distribution of geo-referenced information through 2D-3D maps and interactive plots. The findings were validated through analysis of a hypothetical set of data with an 'ideal' fallout deposit used as a reference, and they are presented in the following section.

Preliminary considerations of the differences in the volume calculations
In this section, we examine some factors with the aim of providing possible explanations about the differences that resulted from the volume estimations obtained by the various interpolation methods used: 1) Stratigraphic considerations: Units B, D and F crop out in a large number of complete sections (<90%), whereas units A, C and E are incomplete in ca. 50% of the sections.
This lack of information is probably due to syn-depositional and post-depositional erosion processes that might be the source of the D volume, noting that such erosion is more likely to occur in the unconsolidated ash and scoriaceous lapilli layers (e.g., units A, C and E) than in the consolidated tuff beds (e.g., units B, D and F). Spatial algorithms reconstructed a thickness surface by interpolation of all of the available measurement points. This surface represents the deposit on the ground and might contain errors due to under/ overestimation of the original thickness data that resulted from erosion, whereas the Pyle method takes into account the overall exponential thinning of the deposits, and this tends to smooth the effects of isolated under/ overestimated data.
2) Input data density for volume computation: Another issue to consider is how much the sampling density of the individual units affected the volume calculation. Below a threshold density, it is clearly impossible to use any of these methods. In our data, the lowest density value was 0.05 point/km 2 (unit F), which allowed us to obtain consistent results using different methodologies. Tables 2 and 3 show that when the highest density of input data was used (mode (iii), points and isopachs), there were no relevant differences in the volumes of all of the units. This shows that the major differences between the Pyle and SInt volumes are independent of the densities.    3) Distal thickness: A further point to investigate is the possible relationship between the measured thickness of the distal deposit and the volume. Volcanic units characterized by a distal thickness of <2 cm generally show TIN volumes in reasonable agreement with volumes estimated using the Pyle method (the distal isopachs of units B, D and F are 0.5, 1 and 2 cm, respectively). The only exception was unit A: although a 1-cm-isopach was reconstructed, the D ranged between 33% and 40%. For the volcanic units characterized by a distal thickness >2 cm (units C and E), there was no agreement between the TIN and the Pyle volumes. On the basis of these data, the tephra volume beyond the distal isopach, which is considered in the Pyle method but not when using the SInt, was not always able to explain D.
4) The size of the eruptive event: Discrepancies in volumes are apparently not correlated to the size of the explosive events that produced each unit. For example, units A and C derive from explosions of the smallest and the largest magnitude, respectively, and these show the largest differences between the Pyle and SInt volumes.

Geo-exploration of data and reconstruction of an 'ideal' fallout dispersion map
As the previous considerations did not fully explain the issue of the different volume calculations, we adopted a new approach for investigating the data. This was based on GIS tools that allowed the analysis of the spatial distributions of the information (in our case, the measured thickness) through 2D maps and interactive plots, such as frequency histograms, correlation plots and 3D perspective graphs, built using the Geo-statistical module of ArcGIS 9.2 software.
Frequency histograms offer an insight into the variability of the data through the analysis of three elements: location, spread and shape. Correlation diagrams (QQplots) showed the quantiles of the input dataset against those of the standard normal distribution, where points distant from the straight line are inconsistent with a normal distribution (outliers). The 3D perspective graphs visualize the data by projecting the thickness of each unit onto the xz and yz planes. These projections can be fitted with polynomial curves, which allows the identification of possible global directional trends.
As mentioned in section 3.2, these analyses were applied to an 'ideal' fallout dispersion to obtain the reference plots for comparison with those derived from the real fallout deposit.
These comparisons might suggest explanations for the discrepancies in the volumes between the Pyle and SInt methods.
We used a grid of points (Figure 5a, one point per km) to create a map of the fallout from an explosive eruption in which the thickness and dispersion of the deposits followed an exponential thinning law that is typical of deposits derived from volcanic plumes affected by winds of constant direction. The grid was built using the 122 BC eruptive vent located at the summit of Etna (3250 m a.s.l.) as a reference source, and an on-land dispersion area of 308 km 2 characterized by a southern dispersal axis and by a fallout thickness ranging from 1.8 m (T max ) to 1 cm (T min ). The volumes calculated by using the Pyle and SInt methods are reported in Table 4, and show that for the units studied, the TIN volumes (mode (ii)) are most similar to those from the Pyle method, as already mentioned in section 3.2. Figure 5b-e shows the geo-statistical plots for this 'ideal' dataset, and it clearly demonstrates that the thickness of the deposits does not follow a normal distribution. The frequency histogram (Figure 5b) shows a strongly positively skewed curve with a very thin tail. The QQplot (Figure 5c) shows a family of curves of the type y = e x −1. The trend of these functions has a very flat portion along the x-axis and evident outliers, especially in the right upper portion of the plot, above the straight line. These outliers correspond to the thick strata that occurred in proximal areas. The 3D perspective graphs (Figure 5d, e) show a bell-shaped trend in the XZ plane and a NNW-SSE exponential trend in the YZ plane. The bell-shaped trend describes a more rapid decrease in thickness from the vent and it indicates the deposits located along the semi-minor dispersal axis. The second trend describes an exponential or linear decay of thickness with distance from the vent, and it represents the deposits along the main dispersal axis (Figure 5d).
These data indicate that the spatial trends of an 'ideal' fallout dispersion followed a regular behavior that was characterized by well-defined statistical distributions. Deviations from these trends might help us to assess the syn-depositional and post-depositional processes that modified the original thicknesses of the primary fallout deposits.

Inferences from geo-spatial analysis
Two-dimensional maps (Figure 3) that were obtained by overlaying the point thickness layer of each unit, the munic-  ipality boundaries, and the contour lines on the digital elevation model of Mount Etna [Tarquini et al. 2007] revealed that the areas most affected by the 122 BC tephra deposits were the municipalities of Nicolosi, Pedara, Trecastagni and the northern sector of Belpasso and Ragalna. Five of the six studied deposits (units A, B, D, E, F) were dispersed on land over an area above an altitude of 500 m; only unit C was partly dispersed offshore (Figure 3c), which reached a thickness of 4 cm south of Priolo, which was located about 70 km SE of the vent. These maps also highlighted the areas where no field data are available (e.g., proximal areas, the Valle del Bove depression, and the sector between Nicolosi and Catania; Figure 3). This lack of data is due to the post-122 BC volcanic deposits that cover the 122 BC deposits (e.g., the 1669 eruption on the south-eastern flank of Etna Volcano). When the geo-statistical methods used to analyze the 'ideal' dispersion were applied to our dataset ( Figure 6), we observed a non-normal distribution of the data. The curves in the frequency histograms and QQplots appeared to distinguish two different groups of data: the first group (units B, D and F) shows trends similar to those of the 'ideal' dispersion, while the second group (units A, C and E) shows trends that differ considerably.
The histograms for the first group reproduced the positively skewed curve of the 'ideal' dispersion despite the lack of thickness intervals. The correlation plots reveal a general (more or less marked) exponential trend. The 3D diagrams highlight a bell-shaped and an exponential trend in the XZ and YZ planes, respectively. In particular, the volcanic unit that best fits the behavior of the 'ideal' eruption was unit D, which had the lowest D, of 0.01% (mode (ii); see Table 3).
In the second group, units A and E show more or less marked polymodal histograms and QQplots that are characterized by an evident clustering of the data (unit A) or by nearly linear correlations (unit E). In both cases the 3D perspective graphs show trends that are not always consistent with those of the 'ideal' dispersion. For instance, the projections of thickness points define a linear trend in the XZ plane for unit A, which indicates a likely rotation of the dispersal axis from SW to SE. Such behavior has been observed previously in the case of, for instance, the 2002-2003 eruption of Etna, during which a change in the wind direction caused the dispersal axis to rotate within a few days, from S to SE [Andronico et al. 2009].
These findings indicate that the thicknesses are under-estimated and over-estimated, which suggested a possible complex depositional mechanism in which erosion processes or rotation of the dispersal axis due to changes in wind direction during the eruption might have had important roles. Note that D TIN (the Pyle-TIN volume difference) was about 40% for both units, and that this value did not change significantly when the volumes were estimated using the IDW algorithm, which is the spatial interpolation that best approximates the expo-nential thinning law (on which the Pyle method is based).
As units A and E do not follow the 'ideal' thickness distribution, this indicates that either the Pyle method gave incorrect volume estimates or the missing distal information represented an important portion of these eruptive units. Both of these hypotheses must be taken into account to explain the discrepancies between the Pyle and TIN volumes. Considering that the major distal thickness belonged to unit E (5 cm), whereas the highest data clustering belonged to unit A, we can hypothesize that the D is due to the relevant lack of distal deposits in the first case, and erosional and reworking processes in the second.
For unit C, the geo-statistical diagrams describe trends that are only partially similar to those of the 'ideal' dispersion. The histogram has a positively skewed curve that is characterized by a very thick tail that lacks some intervals of thickness. The QQplot shows an intermediate pattern between the exponential curve of the ideal dispersion and the linear trend that is typical of a normal distribution; there were several outliers along the upper and lower portions of the straight line. If the outliers above the straight line are not proximal thickness data (as in the QQplot of the 'ideal' dispersion), they indicate overestimated values that were due to accumulation phenomena; vice versa, if outliers below the straight line are proximal data, they indicate eroded deposits.
Finally, the 3D perspective graph highlights an incoherent trend in the depositional mechanisms that occurred along the main dispersal axis: the exponential or linear decay of the thickness with distance from the vent was not clearly reproduced in the YZ plane, whereas the bell-shaped trend that described a more rapid decrease of the thickness along the semi-dispersal axis was well defined in the XZ plane.
All of the above observations suggest that the depositional mechanism for unit C was probably affected by erosion processes or changes in wind direction during sedimentation. However, D for unit C (about 26%) was significantly smaller than that of units A and E (see Table 3). This suggests that the Pyle-SInt differences in volume calculations are caused by strong contributions from the distal deposits that were not considered by the SInt method, because the minimum thickness interpolated was 4 cm, which was recovered at about 70 km from the vent. Thus the distal deposit with thickness <4 cm was lost, and this might represent a significant fraction of the volume of the erupted products.

Conclusions
This study presents new thickness distribution maps and isopach maps that were processed using a GIS-based approach as an alternative to the traditional Pyle method to calculate the volume of fallout deposits (units A-F) of the 122 BC Etna Plinian eruptive sequence. These data indicate that among all the interpolation algorithms used, the TIN interpolation gives the volume estimations of fallout deposits that are the closest to the Pyle values.
Furthermore, a new analysis based on the interactive 2D-3D geo-statistical diagrams was applied to the spatial distribution thickness. Two groups of deposits were identified: (i) The first group (units B, D and F) is characterized by very similar Pyle and SInt volumes (D <13%), and this defines the products of explosive eruptions in which the dispersal axis of the plume did not change direction significantly during sedimentation. This is confirmed by the geo-statistical diagrams that show trends that were very similar to those of the 'ideal' dispersion, and it is supported by the well-preserved deposits, as shown by the field survey (the high percentage of complete sections is related to the lithological characteristics of the tuff layers). Under these conditions, D was very low, due to the negligible amount of missing distal data (the minimum distal thickness recovered was always <2 cm).
(ii) The second group (units A, C and E) is characterized by significantly different Pyle and SInt volume estimations (D >25%), and consists of deposits formed through more complex depositional processes in which global and/or local mechanisms of erosion/accumulation have had significant roles. This is highlighted by the geo-statistical diagrams that show trends that are very different from those of the 'ideal' dispersion. Marked polymodal histograms and clustering of the data in the QQplots indicate possible reworking, overthickening and erosional processes. Three-dimensional perspective graphs that show linear patterns in the XZ plane indicate a likely rotation of the dispersal axis. Outliers in the QQplot identify deposit accumulation or erosion, which was confirmed by the occurrence of lithological facies, as loose scoriaceous lapilli and ash beds, which were more liable to syn-depositional and post-depositional erosion. Under such conditions, the Pyle and SInt volumes differed significantly (D was always >25%). This discrepancy can also be explained to a certain extent by a lack of distal data, as the minimum distal thickness recovered was always >2 cm.
In conclusion, our analyses provide a preliminary case study of the volumes and dispersal areas of fallout deposits, with the aim to better understanding the syn-eruptive and/or post-eruptive processes that were caused by local (e.g., erosion, over-thickening, reworking) and global (e.g., rotation of the dispersal axis due to changes in wind direction) processes that modified the regular/ ideal distributions of the pyroclastic deposits.