 Methodology
 Open Access
A stepbystep evaluation of empirical methods to quantify eruption source parameters from tephrafall deposits
 Sébastien Biass^{1, 2}Email authorView ORCID ID profile,
 Costanza Bonadonna^{3}View ORCID ID profile and
 Bruce F. Houghton^{1}
https://doi.org/10.1186/s1361701800811
© The Author(s) 2019
 Received: 10 May 2018
 Accepted: 30 November 2018
 Published: 15 January 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.
Keywords
 Tephrafall deposits
 Volume
 Mass
 Isopach
 Isomass
 Isopleth
 Eruption classification
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
Summary of all input arguments available to customize the behaviour of the TephraFits function
Variable name  Description  Variable type  

Required  
fitType  Fitting strategie(s), required to be entered as the third input argument. Accepts:  String or cell array  
∙’exponential’: Fit using a single or multiple exponential segments as  of strings  
developed by Pyle (1989) and Fierstein and Nathenson (1992)  
∙’powerlaw’: Power–law fit as developed by Bonadonna and Houghton (2005)  
∙’Weibull’: Weibull fit as developed by Bonadonna and Costa (2012)  
dataType  Type of input data. Accepts:  String  
∙’isopach’ (default)  
−xData: Square root of isopach area (km)  
−yData: Isopach thickness (cm)  
∙’isomass’  
−xData: Square root of isomass area (km)  
−yData: Isomass load (kg m ^{−2})  
∙’isopleth’  
−xData: Square root of isopleth area (km)  
−yData: Clast diameter (cm)  
∙’transect’  
−xData: Distance from source (km)  
−yData: Deposit thickness (cm)  
Fitspecific  Arguments required for specific fitting methods. Entered as the fifth argument onwards as a namevalue pair.  
Exponential  BIS  Location of the break(s)inslope for multiple exponential segments specified as a numeric value for 2 segments or a n−1 vector for n segments. The location specifies the value of xData after which the breakinslope occurs. If BIS is not specified, only one exponential segment is used.  Double or vector of double 
segments  The segments option enables the automatic fitting of exponential segment by minimising residuals. The value should be a 1×2 vector containing the minimum and maximum number of segments to fit.  Double or vector of double  
optimize  Defines what parameter to optimize when the segments option is used (optional). Either ’rms’ to minimise the root–mean square error (default) or ’r2’ to maximise the r^{2}.  String  
Power–law  C  Distal integration limit (km). Used only if fitType contains ’powerlaw’ and if dataType is either ’isopach’ or ’isomass’  Double 
T _{0}  Ordinate of the most proximal segment of the exponential fits. The ordinate should be expressed as a linear value of yData (i.e. not on a logarithmic scale). It is not necessary to define T_{0} when the power–law fit is parsed along the exponential.  Double  
Weibull  lambdaRange  Range of λ values entered as a 1×2 vector containing [min,max] used during the optimisation of the Weibull parameters. lambdaRangemust be specified along with nRange. In the specific case of isopachs where the Weibull method is requested along with any other fit type, the function uses the ranges of λ and n defined by Bonadonna and Costa (2013) as a function of the VEI obtained from the other fits.  1×2 vector of double 
nRange  Same as lambdaRange  
Optional  Optional arguments controlling the behaviour of the function. Entered as the fifth argument onward as a namevalue pair.  
Probabilistic  runMode  Defines if the probabilistic mode is enabled for the characterization of uncertainties. Accepts:  String 
∙’single’: A single fit is performed (default)  
∙’probabilistic’: Multiple runs are performed using Monte Carlo simulations  
The following arguments are only used if ’runMode’ is set to ’probabilistic’  
nbRuns  Number of runs of the probabilistic mode  Double  
xError  Error (in %) on xData. xError can be specified either as a single value, which assumes an equal error for all xData, or as a vector of the same size as xData containing errors on individual points  Double or vector of size of xData  
yError  Same as xError for yData  Double or vector of size of yData  
CError  Error (in %) on the distal integration limit C. Used only if fitType contains ’powerlaw’ and if dataType is either ’isopach’ or ’isomass’  Double  
errorType  Probability density function of the error around the central value used for Monte Carlo simulations. Accepts:  String  
∙’normal’: Gaussian distribution of errors using userdefined error as 3 sigma of the  
distribution (default)  
∙’uniform’: Uniform distribution of errors using userdefined error as extreme  
values  
errorBounds  Percentiles used to express the spread of the final values. Should be specified as a 1×2 vector containing [min,max]. By default, using the 5^{th} and 95^{th} percentiles  1×2 vector of double  
Plotting  Scale  Scale of the yaxis for plotting. Accepts:  String 
∙’log10’: Log 10 logarithm (default)  
∙’ln’: Natural logarithm  
∙’linear’: Linear  
maxDistance  Maximum extent of curve extrapolation in distal part for plotting. 1 means 100%, i.e. the distance to the most distal point is doubled (default)  Double  
fits2plot  Defines which fits to plot. Parsed as a 1 ×length(fitTypes) boolean vector. For example, if fitType = ’exponential’, ’powerlaw and ’fits2plot’ is [1,0], only the exponential fit will be plotted  Boolean vector of size of fitType  
plotType  Plot type. Accepts:  String  
∙’subplot’: Multiple plots in one figure (default)  
∙’separate’: Individual figure for each plot  
∙’none’: No plot 
Input datasets

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

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
Deposit characteristics of Layer 5 of Cotopaxi volcano as described by Biass and Bonadonna (2011)
Isopach  Isopleth  Isomass  Transect DW  Transect XW  

Area\(^{\frac {1}{2}}\)  Thickness  Area\(^{\frac {1}{2}}\)  Diameter  Area\(^{\frac {1}{2}}\)  Load  Distance  Thickness  Distance  Thickness 
km  cm  km  cm  km  kg m ^{−2}  km  cm  km  cm 
7.0  100  23.7  0.8  7.0  1000  13.1  35  4.4  115 
8.9  50  22.9  1  8.9  500  13.4  30  5.0  42 
12.3  30  20.4  1.6  12.3  300  14.2  32  5.3  33 
17.4  20  18.3  2  17.4  200  16.0  18  8.3  23 
21.4  10  14.9  3  21.4  100  16.3  23  8.9  13 
25.5  5  14.1  3.2  25.5  50  23.9  8  
12.2  4  24.8  6 
Function usage & method
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.
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:
Output variables of each fit type. nbSegments is the number of exponential segments
Variable  Description  Size  

Global  
X, Y  X and Y values of the fit used for plotting  n×100  
Ym  Computed Y values at observed X locations  n× number of X values  
Depositspecific  
Isomass  mass_kg  Mass (kg)  1×n 
magnitude  Magnitude (Pyle 2000)  1×n  
Isopach  volume_km3  Volume (km^{3})  1×n 
VEI  Volcanic explosivity index (Newhall and Self 1982; Houghton et al. 2013)  1×n  
Fitspecific  
Exponential  T0  Ordinate of the exponential segments  nbSegments×n 
k  Slope of the exponential segments  nbSegments×n  
I  Sqrt(Area) of intersection of the exponential segments  nbSegments×n  
b  Distance at which the y value halves. See text for details  nbSegments×n  
r2  R^{2} value of the exponential fits  nbSegments×n  
Power–law  m  Power–law coefficient  1×n 
TPL  Power–law exponent  1×n  
r2  R^{2} value of the power–law fit  1×n  
Weibull  theta  θ parameter of the Weibull fit  1×n 
lambda  λ parameter of the Weibull fit (km)  1×n  
n  n parameter of the Weibull fit  1×n  
r2  R^{2} of the Weibull fit  1×n  
H  When dataType is isopleth, plume height (km asl) calculated with equation 7 of Bonadonna and Costa (2013)  1×n  
Other  
range  Range of volume (km^{3}; isopach) or mass (kg; isomass) corresponding to the interval of percentiles specified by errorBound (Table 1)  1×2 
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
Ranges of uncertainties on parameters used for the compilation of isopach and isopleth maps as quantified in the literature
Isopach  Natural variance  30%  
Observational error  9%  
4% (Proximal)  
8% (Medial)  
21% (Distal)  
Data contouring  7%  
1540% (Proximal)  
<10% (Medial)  
2025% (Distal)  
Isopleth  Clast characterization  10%  
Averaging technique  Up to 100% 
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
Eruption classification
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 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!
 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

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.
Selection of publications covering the characterization of tephra–fall deposits
Characterization of tephrafall deposits  
Volume calculation  
Uncertainty assessment  
Declarations
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 .
Authors’ 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
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.
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.
Authors’ Affiliations
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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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).Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 Bonadonna, C, Costa A. Estimating the volume of tephra deposits: A new simple strategy. Geology. 2012; 40(5):415–8.Google Scholar
 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.Google Scholar
 Bonadonna, C, Houghton B. Total grainsize distribution and volume of tephrafall deposits. Bull Volcanol. 2005; 67(5):441–56.Google Scholar
 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.Google Scholar
 Burden, R, Phillips J, Hincks T. Estimating volcanic plume heights from depositional clast size. J Geophys Res. 2011; 116(B11):B11206.Google Scholar
 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.Google Scholar
 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.Google Scholar
 Cas, R, Wright J. Volcanic successions, modern and ancient: a geological approach to processes, products, and successions. London: Allen & Unwin; 1988.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 Fierstein, J, Nathenson M. Another look at the calculation of fallout tephra volumes. Bull Volcanol. 1992; 54(2):156–67.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 Legros, F. Minimum volume of a tephra fallout deposit estimated from a single isopach. J Volcanol Geotherm Res. 2000; 96(1–2):25–32.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 Pyle, D. The thickness, volume and grainsize of tephra fall deposits. Bull Volcanol. 1989; 51(1):1–15.Google Scholar
 Pyle, DM. Sizes of volcanic eruptions. San Diego: Academic Press; 2000, pp. 263–9.Google Scholar
 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.Google Scholar
 Rust, A, Cashman K. Permeability controls on expansion and size distributions of pyroclasts. J Geophys Res. 2011; 116(B11):B11202.Google Scholar
 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.Google Scholar
 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.Google Scholar
 Sparks, R. The dimensions and dynamics of volcanic eruption columns. Bull Volcanol. 1986; 48(1):3–15.Google Scholar
 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.Google Scholar
 Sulpizio, R. J Volcanol Geotherm Res. 2005; 145(3–4):315–36.Google Scholar
 Thorarinsson, S. The eruption of Hekla, 194748, 3, The tephrafall from Hekla, March 29th, 1947. Visindafélag ĺslendinga. 1954; 1:3.Google Scholar
 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.Google Scholar
 Walker, GPL. Explosive volcanic eruptions — a new classification scheme. Geol Rundsch. 1973; 62(2):431–46. ISSN 14321149. https://doi.org/10.1007/BF01840108.Google Scholar
 Walker, GPL. The Taupo pumice: product of the most powerful known (ultraplinian) eruption. J Volcanol Geotherm Res. 1980; 8(1):69–94.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar
 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.Google Scholar