 Methodology
 Open Access
 Published:
A stepbystep evaluation of empirical methods to quantify eruption source parameters from tephrafall deposits
Journal of Applied Volcanology volume 8, Article number: 1 (2019)
Abstract
This paper describes the stepbystep process of characterizing tephrafall deposits based on isopach, isomass and isopleth maps as well as thickness transects at different distances from their source. It covers the most frequently used empirical methods of integration (i.e., exponential, power–law and Weibull) and provides a description of the key physical parameters that can be retrieved from tephrafall deposits. To streamline this process, a Matlab function called TephraFits is proposed, which is highly customizable and also guides the interpretation of the results. The function calculates parameters such as the deposit volume/mass, the VEI/magnitude, and the rates of thickness–decay away from the source and assists in eruption classification using deposit–based schemes. The function also contains a stochastic mode that can be used to propagate the uncertainty from field data to the quantification of eruption source parameters. The use of this function is demonstrated using the the 1180 ±80 years B.P. andesitic subplinian/Plinian tephra deposit Layer 5 of Cotopaxi volcano, Ecuador. In addition, we constrain the often delicate choice of the distal integration limit of the power–law method from synthetic deposits produced with the advection–diffusion model Tephra2.
Introduction
The geometry of tephra deposits reflects a combination of eruptive style, intensity and environmental conditions. Quantifying this geometry supplies parameters that constrain physical processes of the generation, transportation and sedimentation of the ejecta. The geometry is typically described by field mapping of the spatial distribution of thickness or mass/area, maximum size of pyroclasts and grainsize distribution, which are then interpolated into isopach, isomass and isopleth maps and fitted based on various strategies in order to characterize the deposit in 3 dimensions. Parameters retrieved through this process are the basis of various quantitative, field–based characterizations of tephra deposits (e.g. Thorarinsson 1954, Walker 1980, 1973, Wilson and Walker 1987, Newhall and Self 1982, Sparks 1986, Carey and Sparks 1986, Cas and Wright 1988). To this day, no publication compiles a detailed step–by–step guideline to the calculation of physical parameters derived from field–based studies of tephrafall deposits, which has induced some degree of confusion on the use of the range of available techniques (e.g. calculation of the tephra volume using the exponential technique of Pyle (1989) versus Fierstein and Nathenson (1992); Nathenson (2017)).
The scope of this paper is, therefore, to provide practical guidelines from field observations to physical constraints. First, it describes the critical steps of the characterization of tephrafall deposits using the most common methods found in the literature, namely exponential, power–law and Weibull approaches. This first step attempts to clarify the entire process and provide the enduser with both practical examples and references to the original seminal papers. Second, it presents a new Matlab function named TephraFits. This function streamlines the calculations behind the characterization of tephrafall deposits and provides guidance to critically interpret the results. It is a commandbased implementation of existing techniques and allows the use of a wide range of input types to quantify physical and empirical parameters required for the reconstruction of eruption source parameters.
Third, it attempts to quantify the distal integration limit of the power–law integration method. This parameter can have a first–order importance when estimating the tephra volume using the power–law method but has never been systematically constrained. To overcome this limitation and support the integration of a power–law fit, we used the advection–diffusion model Tephra2 (Bonadonna et al. 2005) to produce synthetic deposits for various eruption conditions, which were used to quantify the ratio of the mass contained within a given distal isopach to the total mass simulated.
To sum up, this manuscript provides practical guidelines to the characterization of tephrafall deposits using empirical techniques. We illustrate the main methodology by applying TephraFits to a well–studied subplinian–Plinian tephrafall deposit (i.e. Layer 5 of Cotopaxi volcano, Ecuador; see “Case–study” section). This application is meant only to show the use of TephraFits and is not meant to explore the limitations associated with a range of deposit exposure, which will follow the same limitations associated with the individual strategies (see Bonadonna et al. 2015, for a review).
For clarity, three points should be noted. First, this manuscript does not present new methods per se and assumes that all methods implemented in TephraFits are validated in the original papers cited throughout the text. Second, only empirical methods based on the subjective interpretation of isolines are considered. The use of valuable alternative methods (e.g. Connor and Connor 2006, Burden et al. 2011, 2013, Engwell et al. 2015, Yang and Bursik 2016) is outside the scope of this paper. Third, for specific details of individual methods the reader is referred to the original papers; a list of references covering the theories and applications of these methods is provided as an Appendix.
Background
The method of Pyle (1989) is based on the observation that “various parameters of tephra–fall deposits, notably thickness, maximum clast size and median diameter, decrease in a linear manner when the logarithm of the parameter is plotted against distance” (Thorarinsson 1954). Pyle (1989) also introduced the use of the square root of the area A\(^{\frac {1}{2}}\) of contour maps (e.g. isopach, isopleth) as a proxy for distance from the source, which eliminates complications due to complex contour geometries, mostly due to wind effects.
TephraFits adds to a collection of codes dedicated to the characterization of tephrafall deposits such as AshCalc (Daggitt et al. 2014) and TError (Biass et al. 2014). TephraFits is command–linebased, highly customizable and combines updates on computational aspects of AshCalc and the stochastic approach to uncertainty assessment of TError, but extends these capabilities to work not only with isopach data but also on isomass maps, isopleth maps and thickness transects. It is maintained on GitHub at https://github.com/e5k/TephraFits and all new updates are presented at https://e5k.github.io .
Function design
Table 1 presents a comprehensive list of all input arguments controlling the behaviour of the function.
Input datasets
The TephraFit function is designed to work with four different types of data:

Isopach thickness (cm) as a function of the square–root of isopach area (km) for calculating parameters such as the tephra volume (km^{3}), b_{t} (Pyle 1989) or λ_{TH} (Bonadonna and Costa 2013);

Isomass load (kg m ^{−2}) as a function of the square–root of isomass area (km) for calculating the tephra mass (kg);

Isopleth diameter (cm) as a function of the square–root of isopleth area (km) for calculating such parameters as b_{c} (Pyle 1989) or λ_{ML} (Bonadonna and Costa 2013);

Outcrop thickness (cm) as a function of the distance from the vent (km) for calculating thinning rates along single transects (Houghton et al. 2004).
Fitting strategies
Three widely used fitting methods are implemented in TephraFit, which can be used with any of the four types of dataset described above. These are:

The exponential method (Fierstein and Nathenson 1992; Bonadonna and Houghton 2005);

The power–law method (Bonadonna and Houghton 2005);

The Weibull method (Bonadonna and Costa 2012; 2013).
Case–study
We use the 1180 ±80 years B.P. Layer 5 of Cotopaxi volcano (Ecuador) as a case study to illustrate the function’s capabilities (Fig. 1; Table 2). Layer 5 is a scoriaceous lapilli fall with a silica content of 58 wt.% (Barberi et al. 1995). Figure 1 shows the deposit as mapped by Biass and Bonadonna (2011). Isopleth maps were compiled by averaging three axes of the five largest lithics from a 0.5 m^{2} area and results in a plume height of 25–30 km when the method of Carey and Sparks (1986) is applied. Biass and Bonadonna (2011) also suggest an eruption duration of 1–2 h. The total grain–size distribution presented by Tsunematsu and Bonadonna (2015) result in median and sorting coefficients of 2.1 and 1.8 ϕ. The application of TephraFit to a range of deposit exposures is outside the scope of the paper, as the limitations associated with deposit exposure are intrinsic to individual models (see Bonadonna et al. 2015 for a review).
Function usage & method
To use the function, navigate to the location of tephraFits.m in the Matlab command. The general usage is:
where xdata and ydata are two vectors of the same size containing the main field–derived values, which vary for each type of calculation. To calculate the volume of Layer 5 from isopach data, xdata and ydata should be entered as:
Additional arguments are considered either required or optional (Table 1). Optional arguments permit the behaviour of the function to be customized and are described in Table 1. Required arguments are those specific to either the type of input data or the fitting approach and are briefly described below.
Fit–specific parameters
The third input argument of the function (fitType) defines the fitting method (or methods) applied to the dataset. Each fitting approach requires different arguments described below. The following examples illustrate how to use one or more fitting methods:
Exponential integration strategy
The exponential approach describes the trend of tephrafall deposits using one or multiple exponential segments (Fierstein and Nathenson 1992; Bonadonna and Houghton 2005). If multiple exponential segments are identified, a break–in–slope can be specified by using the argument BIS followed by the index of xData after which the break–in–slope occurs (Table 1). For the isopach data for Layer 5, a break–in–slope occurs around an A\(^{\frac {1}{2}}\) value of 9.6 km, which is after the second value in the area vector and results in an index value of 2 (Table 2). For the sake of the example, let us assume that a second break–in–slope occurs at an A\(^{\frac {1}{2}}\) value of 18 km, which corresponds to an index value of 4. If only one segment is identified, BIS is not specified.
Alternatively, it is possible to request TephraFits to find the best segments to fit by either minimizing the root–mean square error or by maximising the r^{2}. To do this, the segment argument should be specified and followed by the range of segments to test.
Regardless of some confusion in the literature regarding the definition of thickness half–distance and, therefore, the estimation of the volume of tephra deposits using the exponential approach, we found that equation 12 of Fierstein and Nathenson (1992) and equation 3 of Pyle (1989) provide very similar results. Nonetheless, equation 12 of Fierstein and Nathenson 1992 is used in TephraFit because it is mathemetically more rigourous and therefore applicable to a more general case. For multiple segments, equation 3 of Bonadonna and Houghton (2005) is used. TephraFits also calculates the thickness half–distance of Fierstein and Nathenson (1992) and Nathenson (2017) used in the volume calculation as \(b_{tA}=\frac {\log (2)}{k}\), and the thickness (b_{t}) and maximum clast (b_{c}) half–distance of Pyle (1989) used for eruption classification as \(\frac {\log (2)}{k\sqrt {\pi }}\). For distance transects, a generic distance b at which the y value halves is calculated as \(b=\frac {\log (2)}{k}\), which represents a true distance from the vent.
Power–law integration strategy
Unlike an exponential fit, a power–law function cannot be integrated across the interval [0, inf] and proximal and distal integration limits must be defined. Here, the proximal integration limit is computed from the maximum thickness at the source, which can be derived from the ordinate (i.e. T_{0}) of the most proximal exponential segment (equation 7 of Bonadonna and Houghton (2005)). If only a power–law fit is requested, the argument T_{0} must be specified independently. If a power–law fit is requested along with an exponential fit, TephraFits automatically retrieves the T_{0} value from the most proximal exponential segment. The distal integration limit (km) must be specified using the argument C.
If the power–law fit is requested on isopach or isomass data, TephraFits automatically outputs a sensitivity analysis showing the variability of the volume (or mass) to the distal integration limit. If the power–law fit results in an exponent m<2, the integration of the total deposit volume/mass is highly sensitive to the distal integration limit and TephraFits will display a warning.
Choosing the distal integration limit
Amongst all uncertainties described above, the discrepancy introduced by the choice of the distal integration limit of the power–law method has never been systematically quantified. As already mentioned in Bonadonna and Houghton (2005) and Bonadonna et al. (2015), the volume is particularly sensitive to the proximal or to the distal integration limit when the power–law exponent is > 2 (typically for small deposits) or < 2 (typically for large or poorly exposed deposits), respectively. The proximal integration limit is typically easier to constrain based on the first exponential segment (e.g. Bonadonna and Houghton (2005)). However, the distal integration limit is more difficult to determine when no observations on the distal deposits are available. As a guideline to the choice of C, we used synthetic deposits generated using Tephra2 (Bonadonna et al. 2005) to assess the area covered by the 1 kg/m^{2} isomass (i.e. thickness of 1 mm if a bulk and uniform deposit density of 1000 kg/m^{2} is assumed) as a function of critical eruption source parameters such as plume height, erupted mass, eruption duration and total grainsize distribution (TGSD). Given the complex relations amongst TGSDs with eruption intensity and magnitude (e.g. Costa et al. (2016); Rust and Cashman (2011)), we identified ranges of published eruption source parameters in order to provide general guidance to the choice of the distal integration limit for the spectrum of explosive styles considered in this paper. In particular, sixty–six simulations were performed using plume heights of 5–30 km and eruption durations of 1, 6 and 12 h. The mass of each simulation was calculated as the product of mass eruption rate and duration, where the mass eruption rate is calculated from the plume height using the model of Degruyter and Bonadonna (2012). The effect of the total grainsize distribution (TGSD) was also assessed by varying the mode of a unimodal Gaussian distribution. Available TGSDs typically range between 2.5—2.5 ϕ. These include the 1875 Askja D deposit, Iceland (2.3 ϕ, Plinian; Sparks et al. 1981), the 1996 eruption of Ruapehu, New Zealand (0.8 ϕ, Vulcanian/subplinian; Bonadonna and Houghton, 2005), the 1974 eruption of Fuego, Guatemala (0.05—0.58 ϕ, Vulcanian/subplinian; Rose et al. 2007), the 2001 eruption of Etna, Italy (2.0 ϕ, Stombolian paroxysm/subplinian; Scollo et al. 2007) and the 2008 eruption of Chaiten, Chile (2.6 ϕ, subplinian; Alfano et al. 2016), where ϕ=− log2(d) and d is the bin diameter in mm. To provide an illustration of the effect of coarse and fine TGSD, we modelled two distributions with modes of −1 and 1 ϕ. All simulations used a standard wind profile (Bonadonna and Phillips 2003) with a wind speed of 20 m/s at the tropopause.
Figure 2a shows the square–root of the area of the 1 kg/m^{2} isomass as a function of plume height, eruption duration and TGSD. Figure 2 shows the relative amount of mass contained within the 1 kg/m^{2} isomass and should be viewed as a first–order guide to the choice of the distal integration limit. As an example, a C value of 150 km was defined for Layer 5 of Cotopaxi volcano.
Weibull integration strategy
Following Daggitt et al. (2014), TephraFits fits a 2–parameter Weibull function by minimizing RSE(λ,n)+ln(RSE(λ,n)) and expressing θ as a function of λ and n (Bonadonna and Costa 2012; 2013). In order to constrain the fitting algorithm, it is necessary to define initial search ranges of λ and n. If the Weibull fit is used alongside any other method on either isopach or isomass data, TephraFits sets ranges of λ and n based on the VEI Table 2 of Bonadonna and Costa (2013). If isomass data are used, the VEI is estimated by converting the erupted mass to a tephra volume using a bulk deposit density of 1000 kg m^{−2}. Alternatively, it is possible to specify custom ranges by using the arguments lambdaRange and nRange, which are two 1×2 vectors containing the minimum and maximum intervals of each parameter. Additionally, when the Weibull method is used on isopleth data, TephraFits also estimates the plume height (km above mean sampling elevation) using equation 7 of Bonadonna and Costa (2013).
Note that as identified by Bonadonna and Costa (2012, 2013), not all deposits are equally sensitive to the initial ranges λ and n. Identifying the appropriate lambdaRange and nRange is an iterative process requiring a critical interpretation of the results in order to thoroughly assess the sensitivity of the Weibull fit. Additionally, Weibull fits are highly sensitive to the influence of individual data points, particularly for distal values. A sensitivity analysis of the effect of the most distal value point on the volume calculation is presented in Table 3 of Bonadonna and Costa (2013).
Accessing results
The raw output values computed in TephraFits can be saved by specifying one output argument:
The variable out is a Matlab structure (type doc struct in the Matlab command for detailed documentation). A structure format contains fields, which can be accessed by typing out.fieldname. The output structure from TephraFits contains the field fitProps, which contains fit–specific input parameters, along with one additional field for each fit. All output fields generated by TephraFits are summarized in Table 3.
Examples
The following section provides series of examples that illustrate the characterization of Layer 5. Each example is ready to be used and can be copied and pasted in the Matlab command.
Isopach map and volume calculation
To work with isopach data, specify ’dataType’, ’isopach’ as an input argument and enter the square–root of isopach area (km) as the first argument and the isopach thickness (cm) as the second argument. Alternatively, the default behaviour of TephraFits is to assume isopach data if dataType is not specified. The following example calculates the tephra volume using all fitting approaches. In the first example, neither lambdaRange nor nRange are specified to the Weibull fit and initial search ranges are estimated from the mean VEI resulting from the exponential and power–law methods (i.e. 2–200 and 5–100, respectively; Bonadonna & Costa 2013), which results in an unsatisfactory solution. The second example uses refined ranges resulting in a satisfactory fit.
Isomass map and mass calculation
To work with isomass data, specify ’dataType’, ’isomass’ as an input argument and enter the square–root of isomass area (km) as the first argument and the isomass value (kg m ^{−2}) as the second argument.
Isopleth map
To work with isopleth data, specify ’dataType’, ’isopleth’ as an input argument and enter the square–root of isopleth area (km) as the first argument and the isopleth diameter (cm) as the second argument. Note that in the case of isopleth data, the area under the fitted curve is not integrated; it is therefore not necessary to specify the argument C when used with a power–law fit.
Note that an implementation of the method of Carey and Sparks (1986) to calculate the height of strong plumes using isopleth maps available at https://github.com/e5k/CareySparks86_Matlab.
Downwind and crosswind transects
TephraFits also allows plotting user–defined transects on isopach, isomass or isopleth maps to visualize the decay of a chosen parameter (e.g. thickness, load or maximum clast diameter) as a function of the distance from the vent. This functionality allows comparison of 2D transects of the same deposit, which was used first by Houghton et al. (2004) to compare the thinning trends of the circular 1912 deposit of Novarupta and subsequently by May et al. (2015) and Biass et al. (2018) to compare the thinning of discrete lobes in multilobate deposits of K\(\bar {i}\)lauea. Although this functionality is mostly designed for a visual comparison of transects, the parameters of the exponential fits (e.g. decay trend and half–decay) can quantify differences. To work with thickness transects, specify ’dataType’, ’transect’ as an input argument and enter the outcrop distance from the source (km) as the first argument and the outcrop thickness (cm) as the second argument. As for isopleth data, the C argument of the power–law does not need to be specified. For the sake of the example, we compared here the exponential thinning along the downwind (red points in Fig. 1) and crosswind axes (white points in Fig. 1). Two exponential segments can be recognized on the crosswind transect, and only one in the downwind transect.
Uncertainty assessment
If the argument runMode is set to probabilistic, TephraFits uses the stochastic method of Biass et al. (2014) to characterize the uncertainty of all computed parameters. This requires the definition of an uncertainty around each parameter, which are used by TephraFits to performs nbRuns simulations while randomly sampling input parameters in each run (Table 1). This results in distributions of final values characterized by a median value and a confidence interval defined by the argument errorBounds (the 5^{th}–95^{th} interval by default). Note that error and uncertainty are here considered synonymous. As a guidance, Table 4 provides an overview of the typical range of uncertainties reported in the literature.
Error values on all parameters are specified in TephraFits following the Probabilistic section of Table 1. The input error represents an uncertainty range around each input parameter reference value and is given in percent. The shape of the input error envelope can be controlled by the argument errorType. If set to normal, the input error envelope is a Normal distribution with a 3σ value defined by the user–defined input error. If set to uniform, the error envelope is uniform and bounded by one user–defined input error on each side of the central value.
The error on the xData and yData input data can be specified in two ways. If xError and yError are input as single values, the same error is uniformly applied to all separate values of xData and yData. Alternatively, if xError and yError are input as vectors of the same size as xData and yData, they then define specific error on each point. In the context of isopach data, this is useful since both the thickness and area of distal isopachs are typically affected by larger uncertainties than those of proximal isopachs. The following examples illustrate the use of the uncertainty assessment on the calculation of the volume of Layer 5 with 100 runs of Monte Carlo simulation.
Outputs of the probabilistic runs contain duplicates of any field names appended by the letter P (e.g. volume_km3P) that contain the values of each simulations. Output distributions can easily be explored by using Matlab’s plotting function histogram. The following example illustrates how to display the distribution of the final m exponent of the power–law fit:
Output
Figure 3 is the output of TephraFits for the example shown in the “Uncertainty assessment” section. Figure 3a contains the fits. If runMode is single, only the solid lines are output. If runMode is probabilistic, dashed lines represent the confidence interval defined by errorBounds (Table 1) over all separate fits of a given method. Note that plotting all fits with their confidence intervals can result in illegible plots, in which case the argument fits2plot becomes useful (Table 1). Figure 3b shows the relationship between observed and computed values resulting from any fit, which is valuable to estimate the quality of fit. The 1:1 line represents a theoretical solution where the fit would exactly reproduce the observed values. Figure 3c and d are only displayed if isopach or isomass data are used, in which case they represent volume (or mass; Fig. 3c) and VEI (or magnitude; Fig. 3d). When ran probabilistically, Fig. 3c displays the volume/mass as box–and–whiskers plots, where the black horizontal line is the median, the box represents the interquartile range (i.e. 25–75^{th} percentiles) and the whiskers the interval specified in errorBounds.
If the power–law method is used, a sensitivity analysis of the effect of the distal integration limit C is also output (Fig. 4). This plot is useful to assess the sensitivity of the final value of volume/mass as a function of C and expresses both the variability of absolute values of volume/mass (left y–axis) and a discrepancy relative to the user–input value of C (right y–axis, expressed in % where a negative value represents an underestimation of the user–input value).
Eruption classification
Classifying volcanic eruptions is a difficult task for reasons summarized by Bonadonna et al. (2016). Since the 70’s, classification schemes based on tephrafall deposits have become prominent for dominantly explosive eruptions (e.g. Self and Sparks, 1978, Cas and Wright, 1988, Wright et al., 1980, Walker, 1973, 1980, Newhall and Self, 1982, Pyle, 1989, 2000). Most commonly, tephra–based classifications describe eruptions in terms of magnitude (i.e. ejected volume/mass) and intensity (i.e. volume/mass of ejecta per unit time, often scaled to the plume height; Mastin et al. (2009)). In addition to the VEI and magnitude scales (Newhall and Self 1982; Pyle 2000; Houghton et al. 2013), TephraFits can calculate magnitude–related (i.e. b_{t} and λ_{TH}) and intensity–related (i.e. b_{c} and λ_{ML}) parameters used in the classification schemes of Pyle (1989) and Bonadonna and Costa (2013) (Fig. 5).
If both isopach and isopleth data are available and fitted with the exponential and/or Weibull method, it is possible to use TephraFits to plot the deposit on the classification schemes of Pyle (1989) and Bonadonna and Costa (2013). For this, enter the isopach and isopleth structures obtained from previous TephraFits runs as the first two arguments:
Regarding the classification scheme of Pyle (1989), note that his Figure 6 relies on the b_{tr} and b_{cr} values computed from the most distal exponential segment of both isopach and isopleth data (his Fig. 3 and 10, respectively). TephraFits plots all combinations of available segments, and the user is trusted to critically interpret the results. According to Pyle (1989), most tephra fall deposits show b_{cr}/ b_{tr} ratios of ∼1–2. The classification of Bonadonna and Costa (2013) includes intermediate regions between endmember eruption styles based on an uncertainty of 20% on the plume height. TephraFits plots the ±20% interval around plume heights of 10 km (small–moderate), 14 km (subplinian), 24 km (Plinian) and 41 km (ultraplinian). For instance, should an eruption plot between the subplinian (i.e. Ht ≥14+.2×14=16.8 km) and Plinian regimes (i.e. Ht ≥24−.2×24=21.2 km), the classification of Bonadonna and Costa (2013) would result in an eruption style at the interface between subplinian and Plinian eruptions. In addition, TephraFits also plots a ±30% uncertainty on λ_{TH} and a ±50% uncertainty on λ_{ML}/λ_{TH} (Bonadonna and Costa 2013).
Results and discussion
Interpretation of layer 5
The volume of Layer 5 as calculated by TephraFits varies between 0.23 km^{3} (Weibull), 0.29 km^{3} (2 exponential segments) and 0.38 km^{3} (power–law, C=150 km), all resulting in a VEI of 4. Using a probabilistic approach, the 90% confidence intervals (i.e. 5−95^{th} percentiles interval) are 0.23–0.43 km^{3} (2 exponential segments), 0.31–0.48 km^{3} (power–law) and 0.19–3.53 km^{3} (Weibull). These observations show that volumes derived with the exponential and the Weibull methods are more consistent and lower (although in the same order of magnitude) than the power–law method. This result agrees with Fig. 3a, which shows how the power–law overestimates the thickness in both the extrapolated proximal and distal parts. However, the probabilistic approach reveals that the volume calculated using a Weibull fit is more sensitive to uncertainties on isopach thickness and square–root of isopach areas, with the 95^{th} percentile of the distribution resulting in a VEI of 5, whereas the 5−95^{th} interval of both the exponential and the power–law methods result in a consistent VEI of 4. The power–law fit of Layer 5 has an exponent m<2, suggesting a high sensitivity to the choice of the distal integration limit C. Varying C between 100–1000 km result in volumes of 0.35–0.50 km^{3}, which correspond to –8%–+31% of the value computed using C=150 km, respectively (Fig. 4).
Thinning trends of Layer 5 result in b_{t} values of 1.1 and 2.9 km and b_{tA} values of 1.9 and 5.3 km for proximal and distal exponential segments, respectively. The trend of decay of clast size with distance from the vent results in a b_{c} value of 2.9 km. The λ_{TH} and λ_{ML} values are 13.6 and 16.7 km, respectively, which equate to a plume height of ∼24 km above the mean sampling elevation (Eq. 7 of Bonadonna and Costa (2013)). Converting this plume height to a height above the vent using a mean sampling altitude of 3000 m asl, we applied equation 3 of Bonadonna and Phillips (2003) to calculate a plume corner located 4.5 km away from the vent (e.g. x0=0.2×H_{t}, where x0 is the downwind position of the plume corner and H_{t} is the height of the top of the plume relative to the vent height).
When plotted on the classification scheme of Pyle (1989), combinations of b_{c} and b_{t} values for both isopach segments of Layer 5 suggest a subplinian eruption. The classification of Bonadonna and Costa (2013) suggests a Plinian eruption.
The thickness transect along the crosswind axis reveals a break in slope at a distance of 5.2 km (Fig. 6). Considering a plume corner located 4.5 km away from the source, the proximal segment of the crosswind transect is likely to record a component of material sedimented from the convective plume, whereas the distal segment is dominated by material sedimented from the umbrella cloud. Along the downwind transect, the most proximal outcrop is located at a distance of ∼10 km from the source. Only a single exponential segment can be identified, which probably does not accurately represent the proximal thinning.
The distal integration limit of the power–law integration method
Figure 2 is designed to serve as a first–order guide to the choice of the distal integration limit as a function of plume height, eruption duration and total grain–size distribution. Considering the 1 kg/m^{2} isomass only, it shows how a fine total–grain–size distribution covers a wider area than a coarse one. The 1 kg/m^{2} isomass (i.e. 1 mm isopach contour considering a deposit density of 1000 kg/m^{3}), although not representing the true edge of a deposit, was chosen to permit computations on a domain with a reasonable size; it represents >90% of the total erupted mass for events characterized by coarse TGSDs and 80–95% for fine TGSDs (Fig. 2b).
More specifically, Fig. 2 suggests how a coarsegrained tephrafall deposit produced by a plume height of 25–30 km lasting 1–6 h covers a square–root of the area of ∼100–250 km, with the 1 kg/m^{2} isomass representing >98% of the total erupted mass. In the case of a deposit associated with an eruptive plume of 15 km, the corresponding C would vary between 50 and 170 km depending on TGSD and duration. As a result, the mass or volume could be assessed using the power–law integration providing a range (i.e. mass or volume calculated using C between 50 and 170 km). If more information on TGSD and eruption duration is available, the integration limit could be assessed with more certainty.
For Layer 5 of Cotopaxi volcano, we defined a C value of 150 km. This is justified by an eruption duration of 1–2 h and a median value of the total grain–size distribution of 2.1 ϕ (Biass and Bonadonna 2011; Tsunematsu and Bonadonna 2015), which is coarser than the range of grain–size distributions modelled here.
Limitations of datasets and empirical methods
The characterization of tephrafall deposits as illustrated throughout this manuscript is typically applied on (usually old) deposits that have been subjected to erosion, reworking and burial. This process is often the only method to quantify physical eruption source parameters for the majority of past eruptions. Field mapping documents and quantifies discrete observations of a deposit, from which subjective contouring is the basis for a reconstruction of the entire deposit in three dimensions Klawonn et al. (2014a). This is used to constrain the processes that generated, transported and sedimented the tephra using largely empirical models. Therefore, numerous sources of uncertainties interact and propagate to the final results, and one must interpret and discuss the results critically in the view of the limitations of both the source dataset and the empirical methods used.
In particular, the Weibull method is more sensitive to the distribution of observed data than the exponential and the powerlaw methods as it has more free parameters. As a result, the associated integration is more complex and can introduce purely computational issues that contribute to discrepancies between results. For instance, the optimization algorithms used to fit a Weibull function adopt different approaches and assumptions when implemented in Excel, Python (i.e. AshCalc) or Matlab (i.e. TError, TephraFits). Bonadonna and Costa (2012) already pointed out that not all deposits were equally sensitive to the choice of initial optimization parameters, and whereas some result in a unique well–constrained numerical solution that reflects a geological reality, others can be described by multiple solutions that drastically change the interpretation of the deposit.
Warnings are your friends!
In TephraFits, numerous warnings were implemented to help the user to interpret the data and to inform when the calculations are performed in unconstrained conditions or at the limit of the assumptions of the methods. The most common warnings concern:

1.
The power–law method, when the exponent m<2, making volume estimates highly sensitive to the distal integration limit;

2.
The Weibull method, when solutions for λ or n consistently converge towards the edges of lambdaRange or nRange, implying that the initial ranges should be expanded.
We encourage the users to review the various warnings to critically interpret the results of TephraFits.
Conclusion
This paper aims at providing practical guidelines for the characterization of tephrafall deposits, from the processing of field observations to the quantification of key eruption source parameters (i.e. erupted volume and mass) and eruption classification. The proposed function, TephraFits, is designed to facilitate this process and provides keys to the interpretation of tephrafall deposits and critical warnings when conditions are at the border line of the application of individual models. It also facilitates the interpretation of tephrafall deposits by combining various models for their characterization in one platform with the following functions:

Fit isopach, isomass and isopleth data with exponential, power–law and Weibull functions;

Calculate the eruption mass/volume (Bonadonna and Houghton 2005; Bonadonna and Costa 2012; Fierstein and Nathenson 1992);

Classify eruptions (Bonadonna and Costa 2013; Pyle 1989);

Analyse downwind and crosswind transects for a better analysis of tephra geometry and eruptive style (Houghton et al. 2004);

Propagate the uncertainty of field–based parameters to the calculation of the eruption source parameter (Biass et al. 2014).
We encourage users to contribute to the development of the function and report suggestions and problems to the authors using the dedicated GitHub page.
Appendix
Table 5 contains a non–exhaustive list of publications that covers the background material used in TephraFits.
References
Alfano, F, Bonadonna C, Watt S, Connor C, Volentik A, Pyle DM. Reconstruction of total grain size distribution of the climactic phase of a longlasting eruption: the example of the 2008–2013 Chaitén eruption. Bull Volcanol. 2016; 78(7):46. ISSN 14320819. https://doi.org/10.1007/s0044501610405.
Barberi, F, Coltelli M, Frullani A, Rosi M, Almeida E. Chronology and dispersal characteristics of recently (last 5000 years) erupted tephra of Cotopaxi (Ecuador): implications for longterm eruptive forecasting. J Volcanol Geotherm Res. 1995; 69(3–4):217–39.
Biass, S, Bagheri G, Aeberhard W, Bonadonna C. TError: towards a better quantification of the uncertainty propagated during the characterization of tephra deposits. Stat Volcanol. 2014; 1(2):1–27. https://doi.org/10.5038/2163338X.1.2.
Biass, S, Bonadonna C. A quantitative uncertainty assessment of eruptive parameters derived from tephra deposits: the example of two large eruptions of Cotopaxi volcano, Ecuador. Bull Volcanol. 2011; 73(1):73–90.
Biass, S, Swanson DA, Houghton BF. New perspective on the nineteenthcentury golden pumice deposit of Klauea Volcano In: Poland, MP, Garcia MO, Camp VE, Grunder A, editors. Field Volcanology: A Tribute to the Distinguished Career of Don Swanson. Boulder: Geological Society of America: 2018. p. 227–46. ISBN 9780813725383 https://doi.org/10.1130/2018.2538(10).
Bonadonna, C, Biass S, Costa A. Physical characterization of explosive volcanic eruptions based on tephra deposits: Propagation of uncertainties and sensitivity analysis. J Volcanol Geotherm Res. 2015; 296:80–100. ISSN 03770273. https://doi.org/10.1016/j.jvolgeores.2015.03.009.
Bonadonna, C, Cioni R, Costa A, Druitt T, Phillips J, Pioli L, Andronico D, Harris A, Scollo S, Bachmann O, Bagheri G, Biass S, Brogi F, Cashman K, Dominguez L, Dürig T, Galland O, Giordano G, Gudmundsson M, Hort M, Höskuldsson A, Houghton B, Komorowski JC, Küppers U, Lacanna G, Le Pennec JL, Macedonio G, Manga M, Manzella I, Vitturi MdM, Neri A, Pistolesi M, Polacci M, Ripepe M, Rossi E, Scheu B, Sulpizio R, Tripoli B, Valade S, Valentine G, Vidal C, Wallenstein N. MeMoVolc report on classification and dynamics of volcanic explosive eruptions. Bull Volcanol. 2016; 78(11):84. ISSN 14320819. https://doi.org/10.1007/s004450161071y.
Bonadonna, C, Cioni R, Pistolesi M, Connor C, Scollo S, Pioli L, Rosi M. Determination of the largest clast sizes of tephra deposits for the characterization of explosive eruptions: a study of the IAVCEI commission on tephra hazard modelling. Bull Volcanol. 2013; 75(1):1–15. ISSN 02588900. https://doi.org/10.1007/s0044501206803.
Bonadonna, C, Connor CB, Houghton BF, Connor L, Byrne M, Laing A, Hincks TK. Probabilistic modeling of tephra dispersal: Hazard assessment of a multiphase rhyolitic eruption at Tarawera, New Zealand. J Geophys Res. 2005; 110(B3):2156–202.
Bonadonna, C, Costa A. Estimating the volume of tephra deposits: A new simple strategy. Geology. 2012; 40(5):415–8.
Bonadonna, C, Costa A. Plume height, volume, and classification of explosive volcanic eruptions based on the Weibull function. Bull Volcanol. 2013; 75(8):1–19.
Bonadonna, C, Houghton B. Total grainsize distribution and volume of tephrafall deposits. Bull Volcanol. 2005; 67(5):441–56.
Bonadonna, C, Phillips JC. Sedimentation from strong volcanic plumes. J Geophys Res Solid Earth. 2003; 108(B7):2340. ISSN 21562202. https://doi.org/10.1029/2002JB002034.
Burden, R, Phillips J, Hincks T. Estimating volcanic plume heights from depositional clast size. J Geophys Res. 2011; 116(B11):B11206.
Burden, RE, Chen L, Phillips JC. A statistical method for determining the volume of volcanic fall deposits. Bull Volcanol. 2013; 75(6):1–10. ISSN 02588900. https://doi.org/10.1007/s0044501307074.
Carey, S, Sparks R. Quantitative models of the fallout and dispersal of tephra from volcanic eruption columns. Bull Volcanol. 1986; 48(2–3):109–25.
Cas, R, Wright J. Volcanic successions, modern and ancient: a geological approach to processes, products, and successions. London: Allen & Unwin; 1988.
Cioni, R, Bertagnini A, Andronico D, Cole PD, Mundula F. The 512 AD eruption of Vesuvius: complex dynamics of a small scale subplinian event. Bull Volcanol. 2011; 73(7):789–810. ISSN 02588900. https://doi.org/10.1007/s0044501104543.
Connor, LJ, Connor CB. Inversion is the key to dispersion: understanding eruption dynamics by inverting tephra fallout In: Mader, HM, Coles SG, Connor CB, Connor LJ, editors. Statistics in volcanology. London: Geological Society of London: 2006. p. 231–42.
Costa, A, Pioli L, Bonadonna C. Assessing tephra total grainsize distribution: Insights from field data analysis. Earth Planet Sci Lett. 2016; 443:90–107. ISSN 0012821X. https://doi.org/10.1016/j.epsl.2016.02.040.
Daggitt, ML, Mather TA, Pyle DM, Page S. AshCalc–a new tool for the comparison of the exponential, powerlaw and Weibull models of tephra deposition. J Appl Volcanol. 2014; 3(1):7. ISSN 21915040. https://doi.org/10.1186/2191504037.
Degruyter, W, Bonadonna C. Improving on mass flow rate estimates of volcanic eruptions. Geophys Res Lett. 2012; 39:16. https://doi.org/10.1029/2012GL052566.
Engwell, SL, Aspinall WP, Sparks RSJ. An objective method for the production of isopach maps and implications for the estimation of tephra deposit volumes and their uncertainties. Bull Volcanol. 2015; 77(7):1–18. ISSN 02588900. https://doi.org/10.1007/s004450150942y.
Engwell, SL, Sparks RSJ, Aspinall WP. Quantifying uncertainties in the measurement of tephra fall thickness. J Appl Volcanol. 2013; 2(1):1–12. https://doi.org/10.1186/2191504025.
Fierstein, J, Nathenson M. Another look at the calculation of fallout tephra volumes. Bull Volcanol. 1992; 54(2):156–67.
Houghton, B, Carey RJ. Pyroclastic fall deposits In: Sigurdsson, H, Houghton B, McNutt S, Rymer H, Stix J, editors. Encyclopedia of Volcanoes, 2nd Edition. London: Academic Press: 2015. p. 599–615. ISBN 9780123859389.
Houghton, BF, Swanson DA, Rausch J, Carey RJ, Fagents SA, Orr TR. Pushing the volcanic explosivity index to its limit and beyond: Constraints from exceptionally weak explosive eruptions at Kilauea in 2008. Geology. 2013; 41(6):627–30. ISSN 00917613. https://doi.org/10.1130/G34146.1.
Houghton, BF, Wilson CJN, Fierstein J, Hildreth W. Complex proximal deposition during the Plinian eruptions of 1912 at Novarupta, Alaska. Bull Volcanol. 2004; 66(2):95–133. ISSN 14320819. https://doi.org/10.1007/s0044500302977.
Klawonn, M, Houghton B, Swanson D, Fagents S, Wessel P, Wolfe C. Constraining explosive volcanism: subjective choices during estimates of eruption magnitude. Bull Volcanol. 2014a; 76(2):1–6. ISSN 02588900. https://doi.org/10.1007/s0044501307933.
Klawonn, M, Houghton B, Swanson D, Fagents S, Wessel P, Wolfe C. From field data to volumes: constraining uncertainties in pyroclastic eruption parameters. Bull Volcanol. 2014b; 76(7):1–16. ISSN 02588900. https://doi.org/10.1007/s0044501408391.
Le Pennec, JL, Ruiz GA, Ramón P, Palacios E, Mothes P, Yepes H. Impact of tephra falls on Andean communities: The influences of eruption size and weather conditions during the 1999–2001 activity of Tungurahua volcano, Ecuador. J Volcanol Geotherm Res. 2012; 217–218:91–103. ISSN 03770273. https://doi.org/10.1016/j.jvolgeores.2011.06.011.
Legros, F. Minimum volume of a tephra fallout deposit estimated from a single isopach. J Volcanol Geotherm Res. 2000; 96(1–2):25–32.
Mastin, L, Guffanti M, Servranckx R, Webley P, Barsotti S, Dean K, Durant A, Ewert J, Neri A, Rose WI, Schneider D, Siebert L, Stunder B, Swanson G, Tupper A, Volentik A, Waythomas C. A multidisciplinary effort to assign realistic source parameters to models of volcanic ashcloud transport and dispersion during eruptions. J Volcanol Geotherm Res. 2009; 186(1–2):10–21.
May, M, Carey RJ, Swanson DA, Houghton BF. ReticuliteProducing Fountains From Ring Fractures in Klauea Caldera ca. 1500 CE In: Carey, R, Cayol V, Poland M, Weis D, editors. Hawaiian Volcanoes. Hoboken: Wiley: 2015. p. 351–67. agu geophy edition. ISBN 9781118872079. https://doi.org/10.1002/9781118872079.ch16 https://doi.org/10.1002/9781118872079.ch16.
Nathenson, M. Revised tephra volumes for Cascade Range volcanoes. J Volcanol Geotherm Res. 2017; 341(Supplement C):42–52. ISSN 03770273. https://doi.org/10.1016/j.jvolgeores.2017.04.021.
Newhall, CG, Self S. The volcanic explosivity index (VEI) An estimate of explosive magnitude for historical volcanism. J Geophys Res. 1982; 87(C2):1231–8.
Pyle, D. The thickness, volume and grainsize of tephra fall deposits. Bull Volcanol. 1989; 51(1):1–15.
Pyle, DM. Sizes of volcanic eruptions. San Diego: Academic Press; 2000, pp. 263–9.
Rose, WI, Self S, Murrow PJ, Bonadonna C, Durant AJ, Ernst GGJ. Nature and significance of small volume fall deposits at composite volcanoes: Insights from the October 14, 1974 Fuego eruption, Guatemala. Bull Volcanol. 2007; 70(9):1043–67. ISSN 02588900. https://doi.org/10.1007/s0044500701875.
Rust, A, Cashman K. Permeability controls on expansion and size distributions of pyroclasts. J Geophys Res. 2011; 116(B11):B11202.
Scollo, S, Del Carlo P, Coltelli M. Tephra fallout of 2001 Etna flank eruption: Analysis of the deposit and plume dispersion. J Volcanol Geotherm Res. 2007; 160(1–2):147–64. ISSN 03770273. https://doi.org/10.1016/j.jvolgeores.2006.09.007.
Self, S, Sparks RS. Characteristics of widespread pyroclastic deposits formed by the interaction of silicic magma and water. Bull Volcanol. 1978; 41(3):196–212. ISSN 0366483X. https://doi.org/10.1007/BF02597223.
Sparks, R. The dimensions and dynamics of volcanic eruption columns. Bull Volcanol. 1986; 48(1):3–15.
Sparks, RSJ, Wilson L, Sigurdsson H. The pyroclastic deposits of the 1875 eruption of Askja, Iceland. Phil Trans R Soc London Ser A Math Phys Sci. 1981; 299:241–73.
Sulpizio, R. J Volcanol Geotherm Res. 2005; 145(3–4):315–36.
Thorarinsson, S. The eruption of Hekla, 194748, 3, The tephrafall from Hekla, March 29th, 1947. Visindafélag ĺslendinga. 1954; 1:3.
Tsunematsu, K, Bonadonna C. Grainsize features of two large eruptions from Cotopaxi volcano (Ecuador) and implications for the calculation of the total grainsize distribution. Bull Volcanol. 2015; 77(7):1–12. ISSN 02588900. https://doi.org/10.1007/s0044501509494.
Walker, GPL. Explosive volcanic eruptions — a new classification scheme. Geol Rundsch. 1973; 62(2):431–46. ISSN 14321149. https://doi.org/10.1007/BF01840108.
Walker, GPL. The Taupo pumice: product of the most powerful known (ultraplinian) eruption. J Volcanol Geotherm Res. 1980; 8(1):69–94.
Wilson, L. Explosive Volcanic EruptionsII The Atmospheric Trajectories of Pyroclasts. Geophys J Int. 1972; 30(4):381–92. https://doi.org/10.1111/j.1365246X.1972.tb05822.x.
Wilson, L, Walker G. Explosive volcanic eruptions  VI. Ejecta dispersal in plinian eruptions: the control of eruption conditions and atmospheric properties. Geophys J R Astr Soc. 1987; 89(2):657–79.
Wright, JV, Smith AL, Self S. A working terminology of pyroclastic deposits. J Volcanol Geotherm Res. 1980; 8:315–36. https://doi.org/10.1016/03770273(80)901110.
Yang, Q, Bursik M. A new interpolation method to model thickness, isopachs, extent, and volume of tephra fall deposits. Bull Volcanol. 2016; 78(10):68. ISSN 14320819. https://doi.org/10.1007/s0044501610610.
Acknowledgements
This research was supported by NSF EAR1427357 (SB and BFH). This work comprises Earth Observatory of Singapore contribution no. 196. This research is supported by the National Research Foundation Singapore and the Singapore Ministry of Education under the Research Centres of Excellence initiative. We are grateful to Thomas Wilson and Jan Lindsay for their excellent editorial work, Kristi Wallace and one anonymous reviewer whose comments and suggestions greatly improved the manuscript. TephraFits relies on dependencies available on Matlab File Exchange website. We are particularly grateful to J. Landsey for the function bplot, J. D’Errico for the function fminsearchbnd, J. Wells for the function RMSE and A. Horchler for the function waittext, who retain the full copyright for their work.
Availability of data and materials
TephraFits is maintained on GitHub at https://github.com/e5k/TephraFits .
Author information
Affiliations
Contributions
SB wrote the function, which was imagined by SB, CB and BFH. SB wrote the manuscript with inputs from CB and BFH. All authors read and approved the final manuscript
Corresponding author
Correspondence to Sébastien Biass.
Ethics declarations
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Biass, S., Bonadonna, C. & Houghton, B.F. A stepbystep evaluation of empirical methods to quantify eruption source parameters from tephrafall deposits. J Appl. Volcanol. 8, 1 (2019) doi:10.1186/s1361701800811
Received
Accepted
Published
DOI
Keywords
 Tephrafall deposits
 Volume
 Mass
 Isopach
 Isomass
 Isopleth
 Eruption classification