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

Abstract

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.

Introduction

The geometry of tephra deposits reflects a combination of eruptive style, intensity and environmental conditions. Quantifying this geometry supplies parameters that constrain physical processes of the generation, transportation and sedimentation of the ejecta. The geometry is typically described by field mapping of the spatial distribution of thickness or mass/area, maximum size of pyroclasts and 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.

Background

The method of Pyle (1989) is based on the observation that “various parameters of tephra–fall deposits, notably thickness, maximum clast size and median diameter, decrease in a linear manner when the logarithm of the parameter is plotted against distance” (Thorarinsson 1954). Pyle (1989) also introduced the use of the square root of the area A$^{\frac {1}{2}}$ of contour maps (e.g. isopach, isopleth) as a proxy for distance from the source, which eliminates complications due to complex contour geometries, mostly due to wind effects.

TephraFits adds to a collection of codes dedicated to the characterization of 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 https://github.com/e5k/TephraFits and all new updates are presented at https://e5k.github.io .

Function design

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

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 (km³), b_t (Pyle 1989) or λ_TH (Bonadonna and Costa 2013);

• Isomass load (kg m ⁻²) as a function of the square–root of isomass area (km) for calculating the tephra mass (kg);

• Isopleth diameter (cm) as a function of the square–root of isopleth area (km) for calculating such parameters as b_c (Pyle 1989) or λ_ML (Bonadonna and Costa 2013);

• Outcrop thickness (cm) as a function of the distance from the vent (km) for calculating thinning rates along single transects (Houghton et al. 2004).

Fitting strategies

Three widely used fitting methods are implemented in TephraFit, which can be used with any of the four types of dataset described above. These are:

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

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

• The Weibull method (Bonadonna and Costa 2012; 2013).

Case–study

We use the 1180 ±80 years B.P. Layer 5 of Cotopaxi volcano (Ecuador) as a case study to illustrate the function’s capabilities (Fig. 1; Table 2). Layer 5 is a scoriaceous lapilli fall with a silica content of 58 wt.% (Barberi et al. 1995). Figure 1 shows the deposit as mapped by Biass and Bonadonna (2011). Isopleth maps were compiled by averaging three axes of the five largest lithics from a 0.5 m² 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

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 r². To do this, the segment argument should be specified and followed by the range of segments to test.

Regardless of some confusion in the literature regarding the definition of thickness half–distance and, therefore, the estimation of the volume of tephra deposits using the exponential approach, we found that equation 12 of Fierstein and Nathenson (1992) and equation 3 of Pyle (1989) provide very similar results. Nonetheless, equation 12 of Fierstein and Nathenson 1992 is used in TephraFit because it is mathemetically more rigourous and therefore applicable to a more general case. For multiple segments, equation 3 of Bonadonna and Houghton (2005) is used. TephraFits also calculates the thickness half–distance of Fierstein and Nathenson (1992) and Nathenson (2017) used in the volume calculation as $b_{tA}=\frac {\log (2)}{k}$, and the thickness (b_t) and maximum clast (b_c) half–distance of Pyle (1989) used for eruption classification as $\frac {\log (2)}{k\sqrt {\pi }}$. For distance transects, a generic distance b at which the y value halves is calculated as $b=\frac {\log (2)}{k}$, which represents a true distance from the vent.

Power–law integration strategy

Unlike an exponential fit, a power–law function cannot be integrated across the interval [0, inf] and proximal and distal integration limits must be defined. Here, the proximal integration limit is computed from the maximum thickness at the source, which can be derived from the ordinate (i.e. T₀) of the most proximal exponential segment (equation 7 of Bonadonna and Houghton (2005)). If only a power–law fit is requested, the argument T₀ must be specified independently. If a power–law fit is requested along with an exponential fit, TephraFits automatically retrieves the T₀ value from the most proximal exponential segment. The distal integration limit (km) must be specified using the argument C.

If the power–law fit is requested on isopach or isomass data, TephraFits automatically outputs a sensitivity analysis showing the variability of the volume (or mass) to the distal integration limit. If the power–law fit results in an exponent m<2, the integration of the total deposit volume/mass is highly sensitive to the distal integration limit and TephraFits will display a warning.

Choosing the distal integration limit

Amongst all uncertainties described above, the discrepancy introduced by the choice of the distal integration limit of the power–law method has never been systematically quantified. As already mentioned in Bonadonna and Houghton (2005) and Bonadonna et al. (2015), the volume is particularly sensitive to the proximal or to the distal integration limit when the power–law exponent is > 2 (typically for small deposits) or < 2 (typically for large or poorly exposed deposits), respectively. The proximal integration limit is typically easier to constrain based on the first exponential segment (e.g. Bonadonna and Houghton (2005)). However, the distal integration limit is more difficult to determine when no observations on the distal deposits are available. As a guideline to the choice of C, we used synthetic deposits generated using Tephra2 (Bonadonna et al. 2005) to assess the area covered by the 1 kg/m² isomass (i.e. thickness of 1 mm if a bulk and uniform deposit density of 1000 kg/m² 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/m² isomass as a function of plume height, eruption duration and TGSD. Figure 2 shows the relative amount of mass contained within the 1 kg/m² 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

a Square–root of the 1 kg/m² 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/m² isomass. Small dots represent values for which the 1 kg/m² 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⁻². 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

Examples

The following section provides series of examples that illustrate the characterization of Layer 5. Each example is ready to be used and can be copied and pasted in the Matlab command.

Isopach map and volume calculation

To work with isopach data, specify ’dataType’, ’isopach’ as an input argument and enter the square–root of isopach area (km) as the first argument and the isopach thickness (cm) as the second argument. Alternatively, the default behaviour of TephraFits is to assume isopach data if dataType is not specified. The following example calculates the tephra volume using all fitting approaches. In the first example, neither lambdaRange nor nRange are specified to the Weibull fit and initial search ranges are estimated from the mean VEI resulting from the exponential and power–law methods (i.e. 2–200 and 5–100, respectively; Bonadonna & Costa 2013), which results in an unsatisfactory solution. The second example uses refined ranges resulting in a satisfactory fit.

Isomass map and mass calculation

To work with isomass data, specify ’dataType’, ’isomass’ as an input argument and enter the square–root of isomass area (km) as the first argument and the isomass value (kg m ⁻²) as the second argument.

Isopleth map

To work with isopleth data, specify ’dataType’, ’isopleth’ as an input argument and enter the square–root of isopleth area (km) as the first argument and the isopleth diameter (cm) as the second argument. Note that in the case of isopleth data, the area under the fitted curve is not integrated; it is therefore not necessary to specify the argument C when used with a power–law fit.

Note that an implementation of the method of Carey and Sparks (1986) to calculate the height of strong plumes using isopleth maps available at https://github.com/e5k/CareySparks86_Matlab.

Downwind and crosswind transects

TephraFits also allows plotting user–defined transects on isopach, isomass or isopleth maps to visualize the decay of a chosen parameter (e.g. thickness, load or maximum clast diameter) as a function of the distance from the vent. This functionality allows comparison of 2D transects of the same deposit, which was used first by Houghton et al. (2004) to compare the thinning trends of the circular 1912 deposit of Novarupta and subsequently by May et al. (2015) and Biass et al. (2018) to compare the thinning of discrete lobes in multilobate deposits of K$\bar {i}$lauea. Although this functionality is mostly designed for a visual comparison of transects, the parameters of the exponential fits (e.g. decay trend and half–decay) can quantify differences. To work with thickness transects, specify ’dataType’, ’transect’ as an input argument and enter the outcrop distance from the source (km) as the first argument and the outcrop thickness (cm) as the second argument. As for isopleth data, the C argument of the power–law does not need to be specified. For the sake of the example, we compared here the exponential thinning along the downwind (red points in Fig. 1) and crosswind axes (white points in Fig. 1). Two exponential segments can be recognized on the crosswind transect, and only one in the downwind transect.

Uncertainty assessment

If the argument runMode is set to probabilistic, TephraFits uses the stochastic method of Biass et al. (2014) to characterize the uncertainty of all computed parameters. This requires the definition of an uncertainty around each parameter, which are used by TephraFits to performs nbRuns simulations while randomly sampling input parameters in each run (Table 1). This results in distributions of final values characterized by a median value and a confidence interval defined by the argument errorBounds (the 5^th–95^th interval by default). Note that error and uncertainty are here considered synonymous. As a guidance, Table 4 provides an overview of the typical range of uncertainties reported in the literature.

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:

Output

Figure 3 is the output of TephraFits for the example shown in the “Uncertainty assessment” section. Figure 3a contains the fits. If runMode is single, only the solid lines are output. If runMode is probabilistic, dashed lines represent the confidence interval defined by errorBounds (Table 1) over all separate fits of a given method. Note that plotting all fits with their confidence intervals can result in illegible plots, in which case the argument fits2plot becomes useful (Table 1). Figure 3b shows the relationship between observed and computed values resulting from any fit, which is valuable to estimate the quality of fit. The 1:1 line represents a theoretical solution where the fit would exactly reproduce the observed values. Figure 3c and d are only displayed if isopach or isomass data are used, in which case they represent volume (or mass; Fig. 3c) and VEI (or magnitude; Fig. 3d). When ran probabilistically, Fig. 3c displays the volume/mass as box–and–whiskers plots, where the black horizontal line is the median, the box represents the interquartile range (i.e. 25–75^th percentiles) and the whiskers the interval specified in errorBounds.

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

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. b_t and λ_TH) and intensity–related (i.e. b_c and λ_ML) parameters used in the classification schemes of Pyle (1989) and Bonadonna and Costa (2013) (Fig. 5).

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 b_tr and b_cr values computed from the most distal exponential segment of both isopach and isopleth data (his Fig. 3 and 10, respectively). TephraFits plots all combinations of available segments, and the user is trusted to critically interpret the results. According to Pyle (1989), most tephra fall deposits show b_cr/ b_tr ratios of ∼1–2. The classification of Bonadonna and Costa (2013) includes intermediate regions between 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 km³ (Weibull), 0.29 km³ (2 exponential segments) and 0.38 km³ (power–law, C=150 km), all resulting in a VEI of 4. Using a probabilistic approach, the 90% confidence intervals (i.e. 5−95^th percentiles interval) are 0.23–0.43 km³ (2 exponential segments), 0.31–0.48 km³ (power–law) and 0.19–3.53 km³ (Weibull). These observations show that volumes derived with the exponential and the Weibull methods are more consistent and lower (although in the same order of magnitude) than the power–law method. This result agrees with Fig. 3a, which shows how the power–law overestimates the thickness in both the extrapolated proximal and distal parts. However, the probabilistic approach reveals that the volume calculated using a Weibull fit is more sensitive to uncertainties on isopach thickness and square–root of isopach areas, with the 95^th percentile of the distribution resulting in a VEI of 5, whereas the 5−95^th interval of both the exponential and the power–law methods result in a consistent VEI of 4. The power–law fit of Layer 5 has an exponent m<2, suggesting a high sensitivity to the choice of the distal integration limit C. Varying C between 100–1000 km result in volumes of 0.35–0.50 km³, which correspond to –8%–+31% of the value computed using C=150 km, respectively (Fig. 4).

Thinning trends of Layer 5 result in b_t values of 1.1 and 2.9 km and b_tA values of 1.9 and 5.3 km for proximal and distal exponential segments, respectively. The trend of decay of clast size with distance from the vent results in a b_c value of 2.9 km. The λ_TH and λ_ML values are 13.6 and 16.7 km, respectively, which equate to a plume height of ∼24 km above the mean sampling elevation (Eq. 7 of Bonadonna and Costa (2013)). Converting this plume height to a height above the vent using a mean sampling altitude of 3000 m asl, we applied equation 3 of Bonadonna and Phillips (2003) to calculate a plume corner located 4.5 km away from the vent (e.g. x0=0.2×H_t, where x0 is the downwind position of the plume corner and H_t is the height of the top of the plume relative to the vent height).

When plotted on the classification scheme of Pyle (1989), combinations of b_c and b_t values for both isopach segments of Layer 5 suggest a subplinian eruption. The classification of Bonadonna and Costa (2013) suggests a Plinian eruption.

The thickness transect along the crosswind axis reveals a break in slope at a distance of 5.2 km (Fig. 6). Considering a plume corner located 4.5 km away from the source, the proximal segment of the crosswind transect is likely to record a component of material sedimented from the convective plume, whereas the distal segment is dominated by material sedimented from the umbrella cloud. Along the downwind transect, the most proximal outcrop is located at a distance of ∼10 km from the source. Only a single exponential segment can be identified, which probably does not accurately represent the proximal thinning.

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/m² isomass only, it shows how a fine total–grain–size distribution covers a wider area than a coarse one. The 1 kg/m² isomass (i.e. 1 mm isopach contour considering a deposit density of 1000 kg/m³), 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/m² 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. The power–law method, when the exponent m<2, making volume estimates highly sensitive to the distal integration limit;

• 2. The Weibull method, when solutions for λ or n consistently converge towards the edges of lambdaRange or nRange, implying that the initial ranges should be expanded.

We encourage the users to review the various warnings to critically interpret the results of TephraFits.

Conclusion

This paper aims at providing practical guidelines for the characterization of 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:

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

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

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

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

• Propagate the uncertainty of field–based parameters to the calculation of the eruption source parameter (Biass et al. 2014).

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

Appendix

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

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