Skip to main content

A step-by-step evaluation of empirical methods to quantify eruption source parameters from tephra-fall deposits


This paper describes the step-by-step process of characterizing tephra-fall 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 tephra-fall 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.


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 grain-size 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 tephra-fall 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 tephra-fall 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 end-user 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 tephra-fall deposits and provides guidance to critically interpret the results. It is a command-based 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 tephra-fall deposits using empirical techniques. We illustrate the main methodology by applying TephraFits to a well–studied subplinian–Plinian tephra-fall 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.


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 tephra-fall deposits such as AshCalc (Daggitt et al. 2014) and TError (Biass et al. 2014). TephraFits is command–line-based, 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 and all new updates are presented at .

Function design

Table 1 presents a comprehensive list of all input arguments controlling the behaviour of the function.

Table 1 Summary of all input arguments available to customize the behaviour of the TephraFits 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 (km3), bt (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 bc (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:


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 m2 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).

Fig. 1
figure 1

Deposit of Layer 5 of Cotopaxi volcano showing a isopach (cm) and b wall-rock lithic isopleth (cm) maps (Biass and Bonadonna 2011). Dashed lines suggest a poor degree of constraint. White dots are outcrops used for crosswind transects and red dots for downwind transects

Table 2 Deposit characteristics of Layer 5 of Cotopaxi volcano as described by Biass and Bonadonna (2011)

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 tephra-fall 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 r2. 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 (bt) and maximum clast (bc) 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. T0) of the most proximal exponential segment (equation 7 of Bonadonna and Houghton (2005)). If only a power–law fit is requested, the argument T0 must be specified independently. If a power–law fit is requested along with an exponential fit, TephraFits automatically retrieves the T0 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/m2 isomass (i.e. thickness of 1 mm if a bulk and uniform deposit density of 1000 kg/m2 is assumed) as a function of critical eruption source parameters such as plume height, erupted mass, eruption duration and total grain-size 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 grain-size 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/m2 isomass as a function of plume height, eruption duration and TGSD. Figure 2 shows the relative amount of mass contained within the 1 kg/m2 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.

Fig. 2
figure 2

a Square–root of the 1 kg/m2 isomass as a function of the plume height for eruption durations of 1, 6 and 12 h and modes of the Gaussian total grain–size distribution of -1 (open symbol) and 1 (solid symbols) and b corresponding percentage of the total erupted mass comprised within the 1 kg/m2 isomass. Small dots represent values for which the 1 kg/m2 isomass extends outside of the computation domain

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.

Table 3 Output variables of each fit type. nbSegments is the number of exponential segments


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

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 5th–95th 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.

Table 4 Ranges of uncertainties on parameters used for the compilation of isopach and isopleth maps as quantified 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:


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–75th percentiles) and the whiskers the interval specified in errorBounds.

Fig. 3
figure 3

Output figure of the probabilistic run on isopach data of Layer 5. a Fits showing the central value (solid line) and the confidence interval specified by errorBounds (dashed lines). The error on each point is shown by grey bars. b Goodness of fit between observed and computed values. c Volume/mass plot shown as single values when ’runMode’, ’single’ or as a box–and–whiskers plot when ’runMode’, ’probabilistic’. d VEI/magnitude plot

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).

Fig. 4
figure 4

Sensitivity analysis of the distal integration limit of the power–law method on the final volume/mass. The left y–axis shows the absolute values of volume/mass as a function of C, whereas the right y–axis shows the discrepancy (%) relative to the user–specified value of C

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 tephra-fall 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. bt and λTH) and intensity–related (i.e. bc and λML) parameters used in the classification schemes of Pyle (1989) and Bonadonna and Costa (2013) (Fig. 5).

Fig. 5
figure 5

Classification of Layer 5 of Cotopaxi volcano using the schemes of aPyle (1989) and bBonadonna and Costa (2013). When using Pyle (1989), one point is plotted for each isopach segment. For Layer 5, the blue and red dots are the proximal and distal segments, respectively. When using Bonadonna and Costa (2013), TephraFits plots a ±30% uncertainty on λTH and a ±50% uncertainty on λML/λTH

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 btr and bcr 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 bcr/ btr ratios of 1–2. The classification of Bonadonna and Costa (2013) includes intermediate regions between end-member 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 km3 (Weibull), 0.29 km3 (2 exponential segments) and 0.38 km3 (power–law, C=150 km), all resulting in a VEI of 4. Using a probabilistic approach, the 90% confidence intervals (i.e. 5−95th percentiles interval) are 0.23–0.43 km3 (2 exponential segments), 0.31–0.48 km3 (power–law) and 0.19–3.53 km3 (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 95th percentile of the distribution resulting in a VEI of 5, whereas the 5−95th 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 km3, which correspond to –8%–+31% of the value computed using C=150 km, respectively (Fig. 4).

Thinning trends of Layer 5 result in bt values of 1.1 and 2.9 km and btA 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 bc 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×Ht, where x0 is the downwind position of the plume corner and Ht 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 bc and bt 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.

Fig. 6
figure 6

Exponential fits for the downwind (blue; 1 segment) and crosswind (red; two segments) transects for Layer 5. The location of the transects is on Fig. 1. The location of the plume corner is calculated following Bonadonna and Phillips (2003)

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/m2 isomass only, it shows how a fine total–grain–size distribution covers a wider area than a coarse one. The 1 kg/m2 isomass (i.e. 1 mm isopach contour considering a deposit density of 1000 kg/m3), 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 coarse-grained tephra-fall 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/m2 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 tephra-fall 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 power-law 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. 1.

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

  2. 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.


This paper aims at providing practical guidelines for the characterization of tephra-fall 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 tephra-fall deposits and critical warnings when conditions are at the border line of the application of individual models. It also facilitates the interpretation of tephra-fall deposits by combining various models for their characterization in one platform with the following functions:

We encourage users to contribute to the development of the function and report suggestions and problems to the authors using the dedicated GitHub page.


Table 5 contains a non–exhaustive list of publications that covers the background material used in TephraFits.

Table 5 Selection of publications covering the characterization of tephra–fall deposits


  • Alfano, F, Bonadonna C, Watt S, Connor C, Volentik A, Pyle DM. Reconstruction of total grain size distribution of the climactic phase of a long-lasting eruption: the example of the 2008–2013 Chaitén eruption. Bull Volcanol. 2016; 78(7):46. ISSN 1432-0819.

    Article  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 long-term eruptive forecasting. J Volcanol Geotherm Res. 1995; 69(3–4):217–39.

    Article  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.

    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.

    Article  Google Scholar 

  • Biass, S, Swanson DA, Houghton BF. New perspective on the nineteenth-century 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

    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.

    Article  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 1432-0819.

    Article  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 0258-8900.

    Article  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.

    Article  Google Scholar 

  • Bonadonna, C, Costa A. Estimating the volume of tephra deposits: A new simple strategy. Geology. 2012; 40(5):415–8.

    Article  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.

    Article  Google Scholar 

  • Bonadonna, C, Houghton B. Total grain-size distribution and volume of tephra-fall deposits. Bull Volcanol. 2005; 67(5):441–56.

    Article  Google Scholar 

  • Bonadonna, C, Phillips JC. Sedimentation from strong volcanic plumes. J Geophys Res Solid Earth. 2003; 108(B7):2340. ISSN 2156-2202.

    Article  Google Scholar 

  • Burden, R, Phillips J, Hincks T. Estimating volcanic plume heights from depositional clast size. J Geophys Res. 2011; 116(B11):B11206.

    Article  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 0258-8900.

    Article  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.

    Article  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 0258-8900.

    Article  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 grain-size distribution: Insights from field data analysis. Earth Planet Sci Lett. 2016; 443:90–107. ISSN 0012821X.

    Article  Google Scholar 

  • Daggitt, ML, Mather TA, Pyle DM, Page S. AshCalc–a new tool for the comparison of the exponential, power-law and Weibull models of tephra deposition. J Appl Volcanol. 2014; 3(1):7. ISSN 2191-5040.

    Article  Google Scholar 

  • Degruyter, W, Bonadonna C. Improving on mass flow rate estimates of volcanic eruptions. Geophys Res Lett. 2012; 39:16.

    Article  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 0258-8900.

    Article  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.

    Article  Google Scholar 

  • Fierstein, J, Nathenson M. Another look at the calculation of fallout tephra volumes. Bull Volcanol. 1992; 54(2):156–67.

    Article  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.

    Article  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 1432-0819.

    Article  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 0258-8900.

    Article  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 0258-8900.

    Article  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.

    Article  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.

    Article  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 ash-cloud transport and dispersion during eruptions. J Volcanol Geotherm Res. 2009; 186(1–2):10–21.

    Article  Google Scholar 

  • May, M, Carey RJ, Swanson DA, Houghton BF. Reticulite-Producing 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.

    Google Scholar 

  • Nathenson, M. Revised tephra volumes for Cascade Range volcanoes. J Volcanol Geotherm Res. 2017; 341(Supplement C):42–52. ISSN 0377-0273.

    Article  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.

    Article  Google Scholar 

  • Pyle, D. The thickness, volume and grainsize of tephra fall deposits. Bull Volcanol. 1989; 51(1):1–15.

    Article  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 0258-8900.

    Article  Google Scholar 

  • Rust, A, Cashman K. Permeability controls on expansion and size distributions of pyroclasts. J Geophys Res. 2011; 116(B11):B11202.

    Article  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.

    Article  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.

    Article  Google Scholar 

  • Sparks, R. The dimensions and dynamics of volcanic eruption columns. Bull Volcanol. 1986; 48(1):3–15.

    Article  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.

    Article  Google Scholar 

  • Sulpizio, R. J Volcanol Geotherm Res. 2005; 145(3–4):315–36.

    Article  Google Scholar 

  • Thorarinsson, S. The eruption of Hekla, 1947-48, 3, The tephra-fall from Hekla, March 29th, 1947. Visindafélag ĺslendinga. 1954; 1:3.

    Google Scholar 

  • Tsunematsu, K, Bonadonna C. Grain-size features of two large eruptions from Cotopaxi volcano (Ecuador) and implications for the calculation of the total grain-size distribution. Bull Volcanol. 2015; 77(7):1–12. ISSN 0258-8900.

    Article  Google Scholar 

  • Walker, GPL. Explosive volcanic eruptions — a new classification scheme. Geol Rundsch. 1973; 62(2):431–46. ISSN 1432-1149.

    Article  Google Scholar 

  • Walker, GPL. The Taupo pumice: product of the most powerful known (ultraplinian) eruption. J Volcanol Geotherm Res. 1980; 8(1):69–94.

    Article  Google Scholar 

  • Wilson, L. Explosive Volcanic Eruptions-II The Atmospheric Trajectories of Pyroclasts. Geophys J Int. 1972; 30(4):381–92.

    Article  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.

    Article  Google Scholar 

  • Wright, JV, Smith AL, Self S. A working terminology of pyroclastic deposits. J Volcanol Geotherm Res. 1980; 8:315–36.

    Article  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 1432-0819.

    Article  Google Scholar 

Download references


This research was supported by NSF EAR-1427357 (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 .

Author information

Authors and Affiliations



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(, 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.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Biass, S., Bonadonna, C. & Houghton, B.F. A step-by-step evaluation of empirical methods to quantify eruption source parameters from tephra-fall deposits. J Appl. Volcanol. 8, 1 (2019).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: