Skip to content


  • Methodology
  • Open Access

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

Journal of Applied VolcanologySociety and Volcanoes20198:1

  • Received: 10 May 2018
  • Accepted: 30 November 2018
  • Published:


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.


  • Tephra-fall deposits
  • Volume
  • Mass
  • Isopach
  • Isomass
  • Isopleth
  • Eruption classification


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


Variable name


Variable type




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)



Type of input data. Accepts:



 ∙’isopach’ (default)


  −xData: Square root of isopach area (km)


  −yData: Isopach thickness (cm)




  −xData: Square root of isomass area (km)


  −yData: Isomass load (kg m −2)




  −xData: Square root of isopleth area (km)


  −yData: Clast diameter (cm)




  −xData: Distance from source (km)


  −yData: Deposit thickness (cm)



Arguments required for specific fitting methods. Entered as the fifth argument onwards as a name-value pair.



Location of the break(s)-in-slope 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 break-in-slope occurs. If BIS is not specified, only one exponential segment is used.

Double or vector of double



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



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




Distal integration limit (km). Used only if fitType contains ’powerlaw’ and if dataType is either ’isopach’ or ’isomass’



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 T0 when the power–law fit is parsed along the exponential.




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



Same as lambdaRange



Optional arguments controlling the behaviour of the function. Entered as the fifth argument onward as a name-value pair.



Defines if the probabilistic mode is enabled for the characterization of uncertainties. Accepts:



 ∙’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’



Number of runs of the probabilistic mode




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



Same as xError for yData

Double or vector of size of yData



Error (in %) on the distal integration limit C. Used only if fitType contains ’powerlaw’ and if dataType is either ’isopach’ or ’isomass’




Probability density function of the error around the central value used for Monte Carlo simulations. Accepts:



 ∙’normal’: Gaussian distribution of errors using user-defined error as 3 sigma of the


distribution (default)


 ∙’uniform’: Uniform distribution of errors using user-defined error as extreme





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 5th and 95th percentiles

1×2 vector of double



Scale of the y-axis for plotting. Accepts:



 ∙’log10’: Log 10 logarithm (default)


 ∙’ln’: Natural logarithm


 ∙’linear’: Linear



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)




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



Plot type. Accepts:



 ∙’subplot’: Multiple plots in one figure (default)


 ∙’separate’: Individual figure for each plot


 ∙’none’: No plot


Note that when specified, default denotes the behavior of the function if a given argument is not defined

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




Transect DW

Transect XW

Area\(^{\frac {1}{2}}\)


Area\(^{\frac {1}{2}}\)


Area\(^{\frac {1}{2}}\)











kg m −2






































































Note that Biass and Bonadonna (2011) did not characterize isomass data in the field. For this example, the isomass data are derived from isopach data using a bulk deposit density of 1000 kg m−3. DW and XW refer to thickness transects along downwind and crosswind dispersal axes, respectively (Fig. 1). Note the lack of proximal data along the downwind axis

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







X, Y

X and Y values of the fit used for plotting




Computed Y values at observed X locations

n× number of X values





Mass (kg)




Magnitude (Pyle 2000)




Volume (km3)




Volcanic explosivity index (Newhall and Self 1982; Houghton et al. 2013)






Ordinate of the exponential segments




Slope of the exponential segments




Sqrt(Area) of intersection of the exponential segments




Distance at which the y value halves. See text for details




R2 value of the exponential fits




Power–law coefficient




Power–law exponent




R2 value of the power–law fit




θ parameter of the Weibull fit




λ parameter of the Weibull fit (km)




n parameter of the Weibull fit




R2 of the Weibull fit




When dataType is isopleth, plume height (km asl) calculated with equation 7 of Bonadonna and Costa (2013)





Range of volume (km3; isopach) or mass (kg; isomass) corresponding to the interval of percentiles specified by errorBound (Table 1)


If the probabilistic mode is disabled, n=1. If the probabilistic mode is activated, TephraFits also outputs a second set of variables containing the outputs of all single Monte Carlo simulations, in which case variable names are appended by a P (e.g. VEIP) and n is the number of runs (nbRuns; Table 1)


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


Natural variance


Engwell et al. (2013)


Observational error


Engwell et al. (2013)


4% (Proximal)

Le Pennec et al. (2012)


8% (Medial)

Le Pennec et al. (2012)


21% (Distal)

Le Pennec et al. (2012)


Data contouring


Engwell et al. (2013)


15-40% (Proximal)

Klawonn et al. (2014a, b)


<10% (Medial)

Klawonn et al. (2014a, b)


20-25% (Distal)

Klawonn et al. (2014a, b)


Clast characterization


Bonadonna et al. (2013)


Averaging technique

Up to 100%

Bonadonna et al. (2013)

Natural variance relates to the variability of a deposit around an outcrop as a consequence of both primary (e.g. heterogeneous sedimentation) and secondary (e.g. reworking, erosion) processes. Observational error relates to the variability of thickness measurement by various operators. Data contouring relates to the variability of subjective choices when contouring isopach maps. Clast characterization relates to the variability in the identification of three characteristic length of a clast. Averaging technique relates to the variability of various techniques use to estimate the maximum clast (e.g. the mean value of a variable number of clasts)

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



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 .

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

Department of Geology and Geophysics, University of Hawaii at Manoa, Honolulu, HI, USA
Earth Observatory of Singapore/Asian School of the Environment, Nanyang Technological University, Singapore, Singapore
Department of Earth Sciences, University of Geneva, Geneva, Switzerland


  1. 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. Scholar
  2. 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.Google Scholar
  3. 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. Scholar
  4. 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
  5. 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 Scholar
  6. 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. Scholar
  7. 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. Scholar
  8. 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. Scholar
  9. 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
  10. Bonadonna, C, Costa A. Estimating the volume of tephra deposits: A new simple strategy. Geology. 2012; 40(5):415–8.Google Scholar
  11. 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
  12. Bonadonna, C, Houghton B. Total grain-size distribution and volume of tephra-fall deposits. Bull Volcanol. 2005; 67(5):441–56.Google Scholar
  13. Bonadonna, C, Phillips JC. Sedimentation from strong volcanic plumes. J Geophys Res Solid Earth. 2003; 108(B7):2340. ISSN 2156-2202. Scholar
  14. Burden, R, Phillips J, Hincks T. Estimating volcanic plume heights from depositional clast size. J Geophys Res. 2011; 116(B11):B11206.Google Scholar
  15. 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. Scholar
  16. 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
  17. Cas, R, Wright J. Volcanic successions, modern and ancient: a geological approach to processes, products, and successions. London: Allen & Unwin; 1988.Google Scholar
  18. 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. Scholar
  19. 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
  20. 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. Scholar
  21. 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. Scholar
  22. Degruyter, W, Bonadonna C. Improving on mass flow rate estimates of volcanic eruptions. Geophys Res Lett. 2012; 39:16. Scholar
  23. 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. Scholar
  24. Engwell, SL, Sparks RSJ, Aspinall WP. Quantifying uncertainties in the measurement of tephra fall thickness. J Appl Volcanol. 2013; 2(1):1–12. Scholar
  25. Fierstein, J, Nathenson M. Another look at the calculation of fallout tephra volumes. Bull Volcanol. 1992; 54(2):156–67.Google Scholar
  26. 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
  27. 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. Scholar
  28. 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. Scholar
  29. 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. Scholar
  30. 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. Scholar
  31. 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. Scholar
  32. 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
  33. 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.Google Scholar
  34. 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. Scholar
  35. Nathenson, M. Revised tephra volumes for Cascade Range volcanoes. J Volcanol Geotherm Res. 2017; 341(Supplement C):42–52. ISSN 0377-0273. Scholar
  36. 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
  37. Pyle, D. The thickness, volume and grainsize of tephra fall deposits. Bull Volcanol. 1989; 51(1):1–15.Google Scholar
  38. Pyle, DM. Sizes of volcanic eruptions. San Diego: Academic Press; 2000, pp. 263–9.Google Scholar
  39. 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. Scholar
  40. Rust, A, Cashman K. Permeability controls on expansion and size distributions of pyroclasts. J Geophys Res. 2011; 116(B11):B11202.Google Scholar
  41. 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. Scholar
  42. 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. Scholar
  43. Sparks, R. The dimensions and dynamics of volcanic eruption columns. Bull Volcanol. 1986; 48(1):3–15.Google Scholar
  44. 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
  45. Sulpizio, R. J Volcanol Geotherm Res. 2005; 145(3–4):315–36.Google Scholar
  46. 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
  47. 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. Scholar
  48. Walker, GPL. Explosive volcanic eruptions — a new classification scheme. Geol Rundsch. 1973; 62(2):431–46. ISSN 1432-1149. Scholar
  49. Walker, GPL. The Taupo pumice: product of the most powerful known (ultraplinian) eruption. J Volcanol Geotherm Res. 1980; 8(1):69–94.Google Scholar
  50. Wilson, L. Explosive Volcanic Eruptions-II The Atmospheric Trajectories of Pyroclasts. Geophys J Int. 1972; 30(4):381–92. Scholar
  51. 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
  52. Wright, JV, Smith AL, Self S. A working terminology of pyroclastic deposits. J Volcanol Geotherm Res. 1980; 8:315–36. Scholar
  53. 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. Scholar


© The Author(s) 2019