- Research
- Open Access

# AshCalc–a new tool for the comparison of the exponential, power-law and Weibull models of tephra deposition

- Matthew L Daggitt
^{1}Email author, - Tamsin A Mather
^{2}Email author, - David M Pyle
^{2}and - Stephen Page
^{3}

**3**:7

https://doi.org/10.1186/2191-5040-3-7

© Daggitt et al.; licensee Springer. 2014

**Received:**3 October 2013**Accepted:**14 April 2014**Published:**9 May 2014

## Abstract

This paper presents a new tool, AshCalc, for the comparison of the three most commonly used models for the calculation of the bulk volume of volcanic tephra fall deposits: the exponential model, the power law model and the Weibull model. AshCalc provides a simple and intuitive tool to speed up the analysis of tephra deposits and compare and contrast the fits for each model. Two improvements in terms of computational performance are implemented in AshCalc for the estimation of the parameters for the Weibull model. The first is an analytic method for reducing the number of free parameters, whilst the second exaggerates the minima in parameter space, leading to a more robust solution. We show that AshCalc provides volume estimates in line with other previously published estimates and hence can be used with a high degree of confidence. We include the open source python code for Ashcalc with the intention that it can be used both as a stand-alone program and integrated into other python projects.

## Keywords

- AshCalc
- Tephra volumes
- Exponential
- Power-law
- Weibull

## Introduction

Various models have been proposed for thickness as a function of isopach area, T(x). The most widely used to date are the multi-segment exponential model (Pyle, 1989;1990;1995), the power-law model (Bonadonna et al. 1998; Bonadonna and Houghton, 2005) and more recently a Weibull distribution has been proposed as an alternative model (Bonadonna & Costa, 2012;2013). Each model has its own advantages and drawbacks and their suitability may depend on the characteristics of the deposit. Despite the multiple models, there is, as yet, no routine way to facilitate comparisons between the different approaches. For this reason, we have developed AshCalc, which is designed to do just this.

An alternative approach to volume estimation has recently been proposed by Burden et al. (2013), which uses a purely statistical method to determine the volume of fall deposits without relying on the construction of isopach maps. This is certainly a viable and complementary approach, but it is not one that we shall explore further. As Burden et al. (2013) showed, the results of their approach are consistent with those derived from isopach-based models.

### Introduction to ash volume estimate methodologies

#### Exponential model

where c represents the theoretical maximum thickness at the vent and m characterises the rate of decrease in tephra thickness. Taking the logarithm of both sides linearises the equation and the parameters, c and m, are found by applying least squares regression.

By fitting multiple exponential segments to the data, the exponential law is capable of modelling deposits which have a variable rate of thinning with distance away from the vent–a feature which is common to many well constrained deposits (e.g. Pyle, 1989;1990; Fierstein and Nathanson, 1992; Bonadonna and Houghton 2005; Watt et al. 2009). This approach lends itself to a variety of simple approximations, for example to determine the minimum volume of a deposit from sparse information (e.g. Pyle, 1995;1999; Legros, 2000; Sulpizio, 2005).

#### Power law model

Where c is a linear scaling factor and m characterises the rate of decrease in tephra thickness. Taking the log of both sides linearises the equation and the parameters, c and m, are found by applying least squares regression.

The main disadvantage of the power law model is that T(x) is not integrable between 0 and ∞ and therefore proximal and distal limits of integration have to be chosen (e.g., Bonadonna and Costa 2012).

#### Weibull model

The additional parameter, λ, allows the Weibull model to capture the variation in thinning rate that necessitates multiple segments in the exponential model.

## Methods

### General inputs and outputs

As input, AshCalc takes a list of thickness(m), area(km^{2}) pairs for each isopach as well as additional model specific parameters (see Sections 2.2-2.4 below). For all models, AshCalc displays a plot of thickness against $\phantom{\rule{0.25em}{0ex}}\sqrt{\mathrm{Area}}$, the total estimated tephra deposit, the estimates for the model parameters, the mean relative squared error in the fit of the model to the data and the equation for T(x). In addition to generating optimal parameters it allows the user to experiment with their own values for the parameters, updating the graphs and volume and error estimates accordingly. The user manual, containing installation instructions, is provided as Additional file 3. Specific details relating to running AshCalc for each model are highlighted below.

### Exponential model

AshCalc uses the scientific python library, SciPy, to perform multiple-segment least squares regression to determine the segment parameters. When more than a single segment is used, AshCalc chooses the segment boundaries as to ensure the continuity of the thickness function where possible.

Model specific inputs: The number of exponential segments

Model specific outputs: A graph of the least squares regression fit, as well as individual segment volumes, equations and estimates.

### Power law model

AshCalc performs simple least squares regression to determine the power law parameters.

Model specific inputs: Proximal and distal limits of integration

Model specific outputs: Graph of the least squares regression fit, graph of the error surface for parameters c and m, an estimate of the proximal bound as described in Bonadonna and Houghton (2005).

### Weibull model

The Weibull model is a three parameter model and hence, unlike the exponential and power law distributions, optimal values for λ, k and θ cannot be determined analytically. Instead an optimisation algorithm must explore the parameter space and come up with a suitable estimation of their true values. Optimisation algorithms work by trying to minimise a function of the parameters. The function should be chosen such that the minimal parameters result in the distribution with the best fit to the data.

Decreasing the number of parameters in an optimisation problem greatly reduces the size of the parameter space needing to be searched. This significantly increases the performance of the optimisation algorithm when finding solutions. The magnitude of the derivatives of the function (gradient) around the minima also affects the performance of algorithm in finding minima.

The RSE is greater than or equal to zero for all parameter values, with zero representing a fit that goes through every data point perfectly.

- 1)
Eliminating theta by deriving an optimal value for it in terms of λ and k and the given isopach data. This turns the three parameter function RSE(λ,k,θ) into the two parameter function RSE(λ,k), resulting in increased performance of the optimisation algorithm.

- 2)
Minimising a new quantity $\mathrm{Error}\left(\mathrm{\lambda},\mathrm{k}\right)=\mathrm{RSE}\left(\mathrm{\lambda},\mathrm{k}\right)+\mathrm{ln}\left(\mathrm{RSE}\left(\mathrm{\lambda},\mathrm{k}\right)\right)$ as alternative to

*RSE*(λ,k). The new function has much larger gradients around the minimum of*RSE*(λ,k), improving the performance of the optimisation algorithms.

In this instance the performance of an algorithm is measured by how many computation steps it has to perform to get to within some pre-defined range around the optimal parameters. As it is measured in computation steps rather than computation time this measure is independent of the hardware the algorithm is run on.

#### Eliminating theta

Eq 1 is a quadratic in θ and hence has a single analytically derivable minimum with respect to θ.

The full derivation of Eq 2 can be found in Additional file 4.

The elimination of θ as a free parameter from the model reduces it to a two parameter model. As well as increasing the performance of the optimisation algorithm this also allows the parameter space to be plotted in three dimensions. This allows a visual inspection of the parameter space by users, aiding in the process of finding appropriate bounds to be set for the parameters in the optimisation algorithm.

#### Emphasising the minima

The function RSE has the drawback that the denominator of each term, the observed thickness ${\mathrm{T}}_{\mathrm{observed}}\left(\mathrm{x}\right)$, maps the squared errors smaller than the observed thickness into the range [0,1]. This mapping results in large flat basins around the minimum of the RSE function, decreasing the performance of optimisation algorithms.

This means that if $\left({\mathrm{\lambda}}_{1},{\mathrm{k}}_{1}\right)$ is a minimum for $\mathrm{ln}\left(\mathrm{RSE}\left({\mathrm{\lambda}}_{1},{\mathrm{k}}_{2}\right)\right)$ then it is also a minimum for $\mathrm{RSE}\left({\mathrm{\lambda}}_{1},{\mathrm{k}}_{1}\right)$. Hence minimising $\mathrm{ln}\left(\mathrm{RSE}\left({\mathrm{\lambda}}_{1},{\mathrm{k}}_{2}\right)\right)$ is equivalent to minimising $\mathrm{RSE}\left({\mathrm{\lambda}}_{1},{\mathrm{k}}_{1}\right)$.

An undesirable property of $\mathrm{ln}\left(\mathrm{RSE}\left({\mathrm{\lambda}}_{1},{\mathrm{k}}_{2}\right)\right)$ is that it flattens the error surface in the range [1, ∞] (Figure 2d). Ideally, for $\mathrm{RSE}\left({\mathrm{\lambda}}_{1},{\mathrm{k}}_{1}\right)\ge 1$ the properties of $\mathrm{RSE}\left({\mathrm{\lambda}}_{1},{\mathrm{k}}_{1}\right)$ are desirable, whilst for $\mathrm{RSE}\left({\mathrm{\lambda}}_{1},{\mathrm{k}}_{1}\right)\le 1$ the properties of $\mathrm{ln}\left(\mathrm{RSE}\left({\mathrm{\lambda}}_{1},{\mathrm{k}}_{2}\right)\right)$ are desired.

This function preserves the ordering of all points on the original error surface as the logarithm is a monotonic function. Therefore minimising Eq 3 is equivalent to minimising Eq 1.

When RSE is in the range [0,1], the logarithmic term dominates Eq 3 mapping it to the range [‒ ∞, 1] (Figure 2e) and when the RSE is in the range [1, ∞] the non-logarithmic term dominates, preserving the greater magnitude of the derivatives (Figure 2f). This transformation reduces the size and flatness of any basins, hence increasing the probability that the optimisation algorithm will find better values for λ and k for a given number of iterations.

#### AshCalc and the Weibull model

AshCalc takes advantage of the two improvements above and uses a variant of random-restart stochastic hill-climbing for its optimisation algorithm. Within each run, the maximum amount each parameter can change decreases after each iteration.

Model specific inputs: Number of hill-climbing runs, number of iterations per hill-climbing run and upper and lower bounds for λ and k. Guidelines on what to choose for these inputs can be found in the user manual in Additional file 2.

Model specific outputs: Graph of the error surface for parameters λ and k.

## Results

In order to assess the accuracy of AshCalc, the program’s outputs were compared to published values for four eruptions. The four were chosen with the aim of being spread over a range of eruptive magnitudes and deposit shapes and, in order of eruptive magnitude, they are the 1992 Cerro Negro eruption (Bonadonna & Costa, 2012), the 2008 Chaitén eruption (Watt et al. 2009), the Rungwe pumice eruption (Fontijn et al. 2011) and the Minoan eruption (Pyle, 1990).

**A comparison between values produced by AshCalc and previously published values for three eruptions of varying magnitude**

Eruption | Source | Exponential (single segment) eruption volume (km | Exponential (two segments) eruption volume (km | Power law eruption volume (km | Weibull eruption volume |
---|---|---|---|---|---|

Chaitén 2008–Whole deposit (isopachs closed) | Watt et al. 2009 | - | 0.17 | - | - |

AshCalc | - | 0.1713 | 0.2158 | 0.2033 | |

Chaitén 2008–Whole deposit (mapped area) | Watt et al. 2009 | - | 0.15 | - | - |

AshCalc | - | 0.1518 | 0.1548 | 0.1566 | |

Chaitén 2008–3rd of May lobe | Watt et al. 2009 | - | 0.05 | - | - |

AshCalc | - | 0.0478 | 0.0506 | 0.0535 | |

Chaitén 2008–6th of May lobe | Watt et al. 2009 | - | 0.03 | 0.04 | - |

AshCalc | - | 0.0347 | 0.0420 | 0.0361 | |

Rungwe pumice eruption–Malawi (1 cm) | Fontijn et al. 2011 | 2.2 | 3.2 | 3.3 | - |

AshCalc | 2.145 | 3.236 | 3.251 | - | |

Rungwe pumice eruption–Malawi (3 cm) | Fontijn et al. 2011 | - | - | 5.8 | - |

AshCalc | - | 4.301 | 5.818 | - | |

Rungwe pumice eruption–Malawi (5 cm) | Fontijn et al. 2011 | - | 5.6 | 8.5 | - |

AshCalc | - | 5.559 | 8.493 | - | |

Minoan eruption | ,1990 | - | 38.5 | - | - |

AshCalc | - | 38.83 | 39.01 | 42.34 |

AshCalc provides values equal, within the number of significant figures provided, to all values calculated using the exponential model by Watt et al. (2009) for Chaitén. Additionally AshCalc’s Weibull and power law fits provide reasonable volumes compared to the exponential estimates by Watt et al. (2009).

For the Rungwe pumice eruption, AshCalc again provides volumes equal to the published values of all fits, within the number of published significant figures. The true value of the single segment exponential fit to 2 s.f. has been acknowledged by K. Fontijn to be 2.1 km^{3} rather than the 2.2 km^{3} published (pers. comm.).

For the Minoan eruption, AshCalc calculates a two segment exponential volume within 1% of the original value.

**A comparison between values produced by AshCalc and previously published values for the 1992 Cerro Negro eruption (Bonadonna & Costa,**
2012
**)**

Resulting Weibull fit | Bonadonna & Costa., 2008 | AshCalc |
---|---|---|

Volume (km | 0.045 | 0.04556 |

Λ | 4.2 | 4.131 |

k | 0.78 | 0.7667 |

Θ | 1.0004 | 1.02325 |

Relative mean squared error | 0.000541037 | 0.00050026 |

In this paper, we have focussed on the analysis of deposit thickness as a tool for determining erupted volume as deposit thicknesses are most readily and commonly measured in the field. An exactly similar approach applies for the analysis of deposit mass per unit area as a way of determining total erupted mass, and the programme can readily be adapted for this purpose by multiplying the final calculated volume by the estimated mean density.

## Conclusions

AshCalc has the potential to be a powerful tool for volcanologists to calculate and compare the various models (exponential, power law and Weibull) for tephra deposits as matches for their data and calculation of volumes. It has been shown that AshCalc provides volume estimates in line with other previously published estimates and hence can be used with a high degree of confidence. We implement two improvements in terms of computational performance for the estimation of the parameters for the Weibull model in AshCalc. The first is an analytic method for reducing the number of free parameters and the second exaggerates the minima in parameter space, leading to a more robust solution. As both a stand-alone program and a series of python modules that can be integrated into other python projects, it is hoped that AshCalc will become an important tool in the volcanology community.

## Declarations

### Acknowledgements

This work was carried out during the tenure of National Environment Research Council (NERC) research experience placements (to MLD and SP). Mather and Pyle acknowledge support from NERC/ESRC grants NE/J020001/1 and NE/J020052/1 (STREVA) and NERC grant NE/I013210/1. Mather further acknowledges funding from the Leverhulme Trust. The authors would like to thank Karen Fontijn and Harriet Rawson for the valuable user feedback during the creation of AshCalc. Karen Fontijn is also acknowledged for helpful discussion on the Rungwe pumice eruption datasets.

## Authors’ Affiliations

## References

- Alfano F, Bonadonna C, Volentik ACM, Connor CB, Watt SFL, Pyle DM, Connor LJ: Tephra stratigraphy and eruptive volume of the May, 2008, Chaitén eruption, Chile.
*Bull Volcanol*2011, 73: 613–630. 10.1007/s00445-010-0428-xView ArticleGoogle Scholar - Bonadonna C, Costa A: Estimating the volume of tephra deposits: a new simple strategy.
*Geology*2012, 40: 415–418. 10.1130/G32769.1View ArticleGoogle Scholar - Bonadonna C, Costa A: Plume height, volume, and classification of explosive volcanic eruptions based on the Weibull function.
*Bull Volcanol*2013, 75: 742.View ArticleGoogle Scholar - Bonadonna C, Houghton BF: Total grain-size distribution and volume of tephra-fall deposits.
*Bull Volcanol*2005, 67: 441–456. 10.1007/s00445-004-0386-2View ArticleGoogle Scholar - Bonadonna C, Ernst G, Sparks R: Thickness variations and volume estimates of tephra fall deposits: the importance of particle Reynolds number.
*J Volcanol Geother Res*1998, 81: 173–187. 10.1016/S0377-0273(98)00007-9View ArticleGoogle Scholar - Burden RE, Chen L, Phillips CJ: A statistical method for determining the volume of volcanic fall deposits.
*Bull Volcanol*2013, 75: 707.View ArticleGoogle Scholar - Fierstein JN: Another look at the calculation of fallout tephra volume.
*Bull Volcanol*1992, 54: 156–167. 10.1007/BF00278005View ArticleGoogle Scholar - Fontijn K, Ernst GG, Bonadonna C, Elburg MA, Mbede E, Jacobs P: The ~4-ka Rungwe Pumice (South-Western Tanzania): a wind-still Plinian eruption.
*Bull Volcanol*2011, 73: 1353–1368. 10.1007/s00445-011-0486-8View ArticleGoogle Scholar - Legros F: Minimum volume of a tephra fallout deposit estimated from a single isopach.
*J Volcanol Geotherm Res*2000, 96: 25–32. 10.1016/S0377-0273(99)00135-3View ArticleGoogle Scholar - Pyle D: Quaternary tephra in Africa.
*Glob Planet Chang*1999, 21: 95–112. 10.1016/S0921-8181(99)00009-0View ArticleGoogle Scholar - Pyle DM: The thickness, volume and grainsize of tephra fall deposits.
*Bull Volcanol*1989, 51: 1–15.View ArticleGoogle Scholar - Pyle DM: New volume estimates for the Minoan eruption.
*Thera Aegean World III*1990, 2: 113–121.Google Scholar - Pyle DM: Assessment of the minimum volume of tephra fall deposits.
*J Volcanol Geotherm Res*1995, 69: 379–382. 10.1016/0377-0273(95)00038-0View ArticleGoogle Scholar - Sulpizio R: Three empirical methods for the calculation of distal volume of tephra-fall deposits.
*J Volcanol Geotherm Res*2005, 145: 315–336. 10.1016/j.jvolgeores.2005.03.001View ArticleGoogle Scholar - Watt SFL, Pyle DM, Mather TA, Martin RS, Matthews NE: Fallout and distribution of volcanic ash over Argentina following the May 2008 explosive eruption of Chaitén, Chile.
*J Geophys Res*2009, 114: B04207.Google Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited.