Bull Volcanol (2013) 75:738

DOI 10.1007/s00445-013-0738-x

COLLECTION: MONOGENETIC VOLCANISM

Modeling thickness variability in tephra deposition

Emily Kawabata · Mark S. Bebbington ·
Shane J. Cronin · Ting Wang

Received: 5 November 2012 / Accepted: 9 June 2013 / Published online: 2 August 2013
© Springer-Verlag Berlin Heidelberg 2013

Abstract The attenuation of tephra fall thickness is most eruption. Of the distributions considered to describe vari-
commonly estimated after contouring isolated and often ability in these measurements, the lognormal performed
irregular field measurements into smooth isopachs, with poorly, due to its tendency to predict a small number of
varying degrees of subjectivity introduced in the process. greatly over-thickened deposits. The Weibull and gamma
Here, we present an explicit description of the variability distributions fitted the data to a very similar degree and
introduced into a semiempirical tephra attenuation relation. produced very similar estimates for the “effective volume,”
This opens the way to fitting models to actual tephra obser- mean wind direction, and mass/thickness attenuation rate.
vations through maximum likelihood estimation, rather than The latter can be inverted to obtain an estimate of the mean
using weighted least-squares estimation on the isopachs. column height. The estimated wind direction, and the col-
The method is illustrated for small-scale basaltic explosive umn height derived from the estimated thickness attenuation
eruptions using a simple, but typical, data set of the actual parameter, agreed very well with the direct observations
tephra thickness data published from the 1973 Heimaey made during the eruption. Augmented by a mixture frame-
work allowing for the incorporation of multiple lobes and/or
vents, the model was able to identify the source and direc-
Editorial responsibility: I.E.M. Smith, Guest Editor tion of tephra deposition for the 1977 Ukinrek Maars
This paper constitutes part of a topical collection: eruptions from only the tephra thickness data.

Smith IEM, Nemeth K, and Ross P-S (eds) Monogenetic volcanism Keywords Tephra · Ukinrek Maars · Aleatory
and its relevance to the evolution of volcanic fields. uncertainty · Statistical method
E. Kawabata
Institute of Fundamental Sciences–Statistics,
Massey University, Introduction
Private Bag 11222, Palmerston North 4442, New Zealand
e-mail: e.kawabata@massey.ac.nz
One of the major hazards associated with explosive erup-
M. S. Bebbington () · S. J. Cronin tions is the dispersal of tephra. Apart from the hazard to
Volcanic Risk Solutions, aviation from fine particles suspended in the upper atmo-
Massey University,
Private Bag 11222, Palmerston North 4442, New Zealand
sphere (Miller and Casadevall 2000), tephra hazards are
e-mail: m.bebbington@massey.ac.nz associated with its deposited depth, loading, grain size,
S. J. Cronin
and electromagnetic and chemical properties. Tephra fall
e-mail: s.j.cronin@massey.ac.nz may cause respiratory illness, damage to buildings and
storm-water infrastructure, render roads impassable, disable
T. Wang electrical distribution networks, contaminate water supplies,
Department of Mathematics and Statistics,
University of Otago,
destroy crops, kill livestock, and drastically change land-
P.O. Box 56, Dunedin 9054, New Zealand scape stability and flood vulnerability (Baxter et al. 1981;
e-mail: Ting.Wang@otago.ac.nz Heiken et al. 1995; Cronin et al. 1998; Stewart et al. 2006).
738, Page 2 of 14 Bull Volcanol (2013) 75:738

Tephra falls may blanket many tens to thousands of square models and the numerical simulations is provided by the
kilometers with deposit geometry typically established by class of “semiempirical models.” These include the model
spot analysis of thickness (or mass) and other properties of Gonzalez-Mellado and De la Cruz-Reyna (2010), which
at individual sites. Rapid local redeposition of tephra via treats the eruptive volume as a parameter, with a wind-based
fluvial and aeolian processes is common and, along with radial dependence term. However, the suggested fitting pro-
compaction, means that considerable measurement error cedure was via estimated isopachs. Rhoades et al. (2002)
can be introduced, increasing with time from deposition. examined a data set of multiple eruptions from Taupo Vol-
The spot tephra thicknesses are used to construct isopach cano, thus incorporating volume as a variable and develop-
maps, with interpolation typically carried out by hand draw- ing a linear relationship between the logarithms of thickness
ing or increasingly by GIS interpolation methods. Isopachs and volume, modulated by an elliptical term for wind. The
(contours of ash thickness) are, in turn, the basis, along with slope and intercept, and the wind-based terms were all esti-
grain-size distribution, for the most widely applied empiri- mated along with their variances. The result was used by
cal methods used for estimation of volcanic eruption proper- Bebbington et al. (2008) to produce a probabilistic tephra
ties and hazard (i.e., total erupted volume, eruption column hazard model for Mt Taranaki.
height, and mass ejection rate; Carey and Sparks 1986; The advantage of the model of Rhoades et al. (2002)
Sparks 1986; Pyle 1989; Fierstein and Nathenson 1992; is that it is simple to invert, both to calculate a volume–
Bonadonna et al. 1998; Sparks et al. 1997; Pyle 2000; frequency curve (Bebbington et al. 2008), and to estimate
Sulpizio 2005). These all involve application of a mathe- the relative likelihood of a tephra being sourced from dif-
matical function representing tephra attenuation along its ferent vents (Bebbington and Cronin 2011). The potential
longest axis, indexed by a small number of parameters, disadvantage of this approach is that fitting is done by
which is fitted to the data using a least-squares approach. using linear regression on the logarithm of the thickness.
Tephra dispersion and attenuation are also estimated Hence, the variability in the deposited thickness at a given
using numerical models based on advection–diffusion point is assumed to be lognormally distributed. It is not
equations, such as HAZMAP (Macedonio et al. 1988; known whether this is the most accurate representation of
Barberi et al. 1990), ASHFALL (Hurst and Turner 1999), variability in tephra deposition. Obviously, the distribution
Tephra 2 (Connor et al. 2001; Bonadonna et al. 2005), must be limited to nonnegative values, but this leaves open
or FALL3D (Costa et al. 2006). While these can be used many questions regarding its skewness and tail character-
for hazard forecasting as part of a Monte Carlo proce- istics. For example, the lognormal distribution has a thick
dure along with probabilistic models of eruption size and tail, implying that very large amounts of overthickening
meteorological conditions (Hurst and Smith 2004; can occur at large distances from the volcano. Other pos-
Bonadonna et al. 2005; Costa et al. 2009), the inverse prob- sible alternatives to describe tephra fall attenuation are the
lem is less amenable due to difficulty in finding the optimal Weibull (Bonadonna and Costa 2012) and gamma distribu-
fit (Connor and Connor 2006). Approaches to inverting tions, which are both capable of having a mode at some
the observed tephra dispersal to estimate the eruptive positive thickness.
parameters usually use a specified wind profile (e.g., Scollo In this paper, we will use the semiempirical model(s) of
et al. 2007, 2008; Kratzmann et al. 2010), neglect wind Gonzalez-Mellado and De la Cruz-Reyna (2010) to estimate
(e.g., Volentik et al. 2010), or apply an exogenously speci- the mean tephra thickness to be expected at a given dis-
fied average (e.g., Pfeiffer et al. 2005; Johnston et al. 2012), tance and azimuth from an eruption vent, combined with a
and often fix other parameters as well. The parameters variety of distributions for the actual thickness given this
being inverted for are either specified in an experimental mean. We will show how the parameters in the semiem-
design, or optimized in a “one-at-a-time” (Johnston et al. pirical model and the variance in this “error distribution”
2012) or downhill simplex scheme (Connor and Connor can be simultaneously estimated using maximum likelihood
2006). There is no objective measure of “best fit”; typically, methods. In a significant difference from common practice,
some form of weighted least-squares error (Costa et al. all our estimation is based on the actual individual tephra
2009) is minimized. The results are that the solutions are thickness measurements, rather than on the derived
possibly nonoptimum, not least due to leaving wind out isopachs. For this purpose, the 1973 Heimaey eruption
of the design, and can be nonunique (Pfeiffer et al. 2005; (Thorarinsson et al. 1973) is used as an example of explosive
Scollo et al. 2007; Kratzmann et al. 2010; Bonasia et al. basaltic (Strombolian) eruption, because it was observed
2010; Volentik et al. 2010; Johnston et al. 2012) particularly throughout and its published isopach and thickness data are
for relatively sparse deposit data, due to the dependencies very typical (Self et al. 1974), in being both sparse and not
among the parameters involved. entirely consistent with drawn isopachs.
For modeling tephra fall attenuation from single erup- The remainder of the paper is structured as follows: we
tions, a compromise between the simple tephra attenuation will next review the case study data set, that of the 1973
Bull Volcanol (2013) 75:738 Page 3 of 14, 738

Heimaey eruption (Thorarinsson et al. 1973), and construct Meteorological Service at the Storhofdi Lighthouse. While
a nonparametric estimate of the variability in the deposited the former wind direction resulted in measurable deposition
thickness. The model(s) are developed in the following sec- on land, an estimated 65 % of the tephra deposited under the
tion, with the more mathematical details relegated to the influence of the westerly winds fell into the sea and were
Appendix. The results of fitting the models to the Heimaey not measured. The resulting single-lobe pattern is ideal for
data are shown in the “Results” section, along with an this study. The data consist of observed thicknesses at 36
assessment of how well the various models explain the locations, ranging between 3 and 450 cm, as shown in Fig. 1.
data and its variation. A sensitivity analysis is presented in The actual thickness observed at any given point differs
section “Sensitivity analyses.” The model is extended in sec- from an “ideal” thickness due to small random effects of
tion “Multiple lobes and/or vents” to tackle the problem of wind strength and direction, but primarily due to local redis-
identifying multiple lobes from possibly multiple sources, tribution of fallen tephra by wind, rain, and slope processes
illustrated on data from the 1977 Ukinrek Maars eruption. (especially in highly irregular topography). This difference
This is followed by discussion. can be termed the “sampling error,” which describes the
inherent variability of the observations. The absolute error
(or residual) at a point i is Ai = Oi − Ei , where Oi is
Data the observed thickness at location i, and Ei is the expected
thickness at location i. However, a multiplicative error struc-
In order to determine a suitable distribution for the vari- ture is assumed here—such that the size of the error is
ability in deposited tephra thickness for basaltic explosive proportional to the expected thickness. In other words, the
eruptions, a data set is required without the complications of focus is the “relative error” Ai /Ei ≥ −1. For reasons that
multiple lobes or vents, thus the selection of the eruption on will become obvious, it is shifted to become nonnegative,
24 January 1973 eruption of Eldfell on Heimaey, Iceland. and the (shifted) relative error is considered Ri = Ai /Ei +
The tephra dispersal data (Self et al. 1974) were collected 1 = Oi /Ei . Thus, an observation that aligns exactly with
from 24 January to 1 February 1973. Winds were blow- the expectation has a relative error of one. Values less than
ing predominantly to the northwest and the northeast during or greater than one indicate under- or over-thickening of the
the eruption, according to measurements from the Icelandic tephra, respectively.

Fig. 1 Heimaey: Observed

thickness (in centimeter) and N
locations. The vent is indicated 63°27’N
by the triangle
11

24

45
33 44
63
28 60 150
20 40
50 98
Vestmannaeyjar 250
110 118
5 68
Vestmannaeyjabaer 230
450
52 400
6 13 81 Eldfell
63°26’N

250
7 18

32
200
5
Vestmannaey Jar

7 33

7 8

0 0.4km

20°17’W 20°15’W
738, Page 4 of 14 Bull Volcanol (2013) 75:738

The relative errors are all defined in terms of an expected a

value that is to be calculated from a parametric model. As 1.4
lognormal distribution
the error is thus dependent on the model, and selection Weibull distribution
1.2 gamma distribution
amongst multiple models is required, a model indepen-
dent (nonparametric) of the error distribution is initially
1
estimated in order to discover the candidate parametric
distributions that might be suitable. 0.8

Density
A leave-one-out cross-validation approach to estimating
the expected thickness was applied. The observation at loca- 0.6
tion i is omitted, and a surface is fitted to the remaining 35
data points using a triangulated C1 -continuous interpolating 0.4

surface. The value of this surface at location i is taken to be

Ei . This is repeated for all locations inside the convex hull 0.2

of the data to obtain relative errors, as shown by the his-

0
togram in Fig. 2a. It is clear that the relative sampling error 0 0.5 1 1.5 2
is nonnegative and right-skewed, with a mode less than one. Relative error (= observed thickness/expected thickness)

If we denote our relative error by R, then the lognormal

b 1
 
1 (log R − μ)2
f (R) = √ exp − , (1) 0.9

R 2πσ 2 2σ 2 0.8
Cumulative distribution function

Weibull 0.7

  η  0.6
ηR η−1 R
f (R) = exp − , (2) 0.5
νη ν
0.4

and gamma 0.3

  0.2 Observed errors

R κ−1 −R lognormal distribution
f (R) = κ exp (3) 0.1 Weibull distribution
ω (κ) ω gamma distribution
0
0 0.5 1 1.5 2

distributions are good exemplars of distributions that can Relative error (= observed thickness/expected thickness)

be used to model such errors. All have a location param- Fig. 2 Nonparametric relative sampling error. a Density. b
eter (μ, ν, ω) and a shape parameter (σ, η, κ). As can be Distribution
seen from (1) to (3), they differ mainly in the rate at which
the tail decays. The gamma decays much more rapidly
than the lognormal, and the Weibull may decay faster or
slower than the gamma, depending on the value of the shape Methodology
parameter η. Fitting these distributions to the nonparametric
relative errors via maximum likelihood estimation produces The nonparametric relative errors have been shown to be
the densities superimposed on the data histogram in Fig. 2a. consistent with the Weibull and gamma, and perhaps log-
The lognormal density is less faithful to the data than normal, distributions. Hence, these error distributions were
the other two candidates, due primarily to the righthand tail embedded in the framework of a tephra dispersal model.
of the data not being thick enough, and this is quantified For the reasons outlined in the “intro” section, an empirical
more clearly in Fig. 2b. The maximum distance between the model is required that has the facility to include possi-
data (the step function) and the curve is the Kolmogorov– ble wind effect. Since the Rhoades et al. (2002) model
Smirnov statistic for the distance between distributions, assumes a lognormal error distribution, which may not be
with values of 0.156, 0.122, and 0.124 for the lognormal, the best description, the model of Gonzalez-Mellado and
Weibull, and gamma distributions, respectively. While none De la Cruz-Reyna (2010) is preferred here.
of these are large enough (> 0.264) to reject the distribu- The basic form of the Gonzalez-Mellado and De la Cruz-
tion at the 5 % significance level, they do indicate that the Reyna (2010) model is that the expected tephra fall deposit
Weibull and gamma distributions both fit the data well, and thickness is the product of a power-law decay with dis-
the lognormal slightly less well. tance (Bonadonna and Houghton 2005) and a noncircular
Bull Volcanol (2013) 75:738 Page 5 of 14, 738

term based on wind direction. A power-law decay with dis- depends on the estimated uncertainties (Costa et al. 2009)
tance was preferred to an exponential primarily because and can strongly affect the results. In effect, the MLE
it fits well in both the near and far field. It has also estimation treats the uncertainties explicitly, rather than
been shown that a simple exponential decay may not well approximately.
describe well-preserved tephra deposits, due to distal ash The error distributions are incorporated by treating the
settling differently (Sparks et al. 1992; Rose 1993). Using tephra thickness formula (5) or (6) as a link function giv-
the three exponential segments necessary to model accu- ing the mean of the thickness distribution at the location
rately the thinning of well-preserved deposits (Bonadonna (r, ξ ). The shape parameter of the distribution (σ for the
and Houghton 2005) is also undesirable, because it would lognormal, η for the Weibull, κ for the gamma) becomes an
add another four parameters to the estimation problem. additional parameter to be estimated. The likelihood formu-
The attenuation relation by Gonzalez-Mellado and De la lae are derived in the Appendix. In each case, we see that
Cruz-Reyna (2010) is the coefficient of variation (standard deviation divided by
mean) is a function solely of the shape parameter. In other
words, the coefficient of variation is a constant in all direc-
T (r, U, θ) = γ exp [−βU r(1 − cos θ)] r −α , (4)
tions and at all distances, which is exactly the multiplicative
where T is the tephra thickness (in centimeter) at a distance error structure we require.
from the vent r (in kilometer) in direction θ relative to the Two baseline models (with and without wind) and three
wind direction. If the supposed wind direction is given by φ, error distributions makes for a total of six models to be fit-
and that the direction from the vent to the deposit location is ted. This leaves the question of which model best describes
ξ (both measured in degrees anticlockwise from East), then the data. As all the error distributions have the same num-
θ = ξ - φ, and (4) becomes ber of parameters, this can be decided on the basis of the
likelihood. Model (6) is nested within the model (5), and
T (r, ξ ) = γ exp{−βU r[1 − cos(ξ − φ)]}r −α (5) so the model with more parameters can be justified using
−1
Note that the wind direction φ and speed U (km h ) are a likelihood ratio test. A shorthand for this is the Akaike
considered as the mean (with respect to the eruption rate) Information Criterion (Akaike 1977):
values of the predominant wind, and as such will be esti-
AIC = 2p − 2 log L, (7)
mated from the data (the left-hand side of the equation). The
other parameters to be estimated are α, β and γ . The latter where p is the number of parameters, and log L the log like-
is the expected thickness at 1 km from the vent along the lihood. Smaller Akaike Information Criterions (AICs) indi-
dispersal axis, which is a proxy for the eruption size. The cate better models, with the effect of additional parameters
dimensionless attenuation parameter α is nonnegative, and being compensated for in order to avoid over-fitting.
β is inversely related to the diffusion coefficient, and thus
can be regarded as a proxy for the grain-size distribution.
However, there is an identifiability issue, in that βU cannot Results
be separately estimated, and so βU is considered as a single
variable, reducing the number of parameters to be estimated Table 1 shows the estimated parameters from fitting each
to four. Moreover, βU (along with φ) is also a nuisance of the six possible models to the Heimaey tephra thick-
parameter, i.e., one that may or may not be present in the ness data. To allow the shape parameters of the three
model. If there is no significant wind, then U = 0 and (5) error distributions to be meaningfully compared, they have
reduces to been converted into a coefficient of variation (CV, the
standard deviation divided by the mean), using Eqs. 12,
T (r, ξ ) = γ r −α , (6)
14, and 16 given in the Appendix. We see that the esti-
with only two parameters to be estimated. mated parameters are consistent for all distributions. In
Gonzalez-Mellado and De la Cruz-Reyna (2010) fitted particular, the estimated wind direction φ is approxi-
their model to isopach data, thus avoiding the question of mately NNW, in agreement with the meteorological data
sampling error. However, the actual observed tephra thick- (Self et al. 1974).
nesses and locations are used in our formulation, rather The Weibull distribution has the smallest AIC, indicating
than isopachs, and thus sampling error must be included that it is the best of the three distributions, but the gamma
in our fitting. By incorporating an error distribution, the distribution is not significantly worse. However, the differ-
model can be fitted using standard statistical maximum like- ence in AICs is sufficient to reject the lognormal distribution
lihood estimation (MLE) methods. This avoids the necessity as a description of the sampling error. Moreover, the lognor-
of choosing a weighting scheme in the least-squares mini- mal distribution systematically overestimates the volume of
mization procedure. The choice of such a weighting scheme the eruption (γ ) with respect to the other two distributions.
738, Page 6 of 14 Bull Volcanol (2013) 75:738

Table 1 Estimated parameters of the models improvement from the model (6) without wind, compared
to the model (5) incorporating wind. The chi-squared (2
Baseline model Parameter Error distribution
degrees of freedom) statistic of twice the difference in the
lognormal Weibull gamma
loglikelihood ratios is significant for all three error distribu-
(6) γ 54.0 48.4 48.7 tions, with P values less than 10−8 indicating that the model
α 2.00 1.97 1.97 fit is significantly improved by including wind effects.
CV 1.18 0.70 0.78 The AIC and likelihood ratio test tell us which model
log L −179.4 −175.5 −176.3 is better, but the question of whether the model is a good
AIC 364.8 357.0 358.5 description of the data can only be answered by resid-
(5) γ 83.8 71.0 76.0 ual analysis, i.e., by examining the discrepancy between
βU 1.42 1.48 1.43 the data and the model. For a first visual inspection,
φ 114.1 125.3 119.2 Fig. 3 shows the residual error (observation divided by the
α 1.61 1.88 1.75 expected mean from the model) at each observation loca-
CV 0.59 0.45 0.49 tion. The means and standard deviations of the residual
log L −160.1 −157.0 −157.9 errors are calculated in Table 2. We see that the residuals
AIC 330.2 324.0 325.7 are much closer to one (smaller standard deviation) under
the model (5). The lognormal overestimates thicknesses (the
mean residual, or ratio of observed to estimated, is less than
one), for reasons discussed above, while the Weibull and the
It also has the highest CV, indicating a poorer absolute fit. gamma error distributions appear to be unbiased.
Note that the CV decreases for all the error distributions Although some of the residuals from our models are
with the introduction of wind effects. There is definitely an large, the aim is to explicitly quantify this variation.

2 2 2

Relative Errors
North displacement (km)

North displacement (km)

North displacement (km)

1 1 1
< 0.4
0.4 − 0.75
0.75 − 1
1 − 1.2
0 0 0 1.2 − 1.6
> 1.6

−1 −1 −1

−2 −1 0 −2 −1 0 −2 −1 0
East displacement (km) East displacement (km) East displacement (km)

2 2 2
North displacement (km)

North displacement (km)

North displacement (km)

1 1 1

0 0 0

−1 −1 −1

−2 −1 0 −2 −1 0 −2 −1 0
East displacement (km) East displacement (km) East displacement (km)

Fig. 3 Relative errors. The top row is the model (6) without wind; the bottom row, the model (5) with wind. The columns are for the lognormal,
Weibull, and gamma error distributions, reading left to right. Displacements are relative to the vent location, indicated by the triangle
Bull Volcanol (2013) 75:738 Page 7 of 14, 738

Table 2 Residual error statistics error distributions exhibit behavior consistent with the data,
although it should be noted that the thick tail (producing
Baseline model Error distribution Mean SD
occasional very large relative residuals) in the lognormal
(6) lognormal 0.899 0.599 distribution is not being examined by this procedure.
Weibull 1.004 0.669
gamma 1.000 0.666
(5) lognormal 0.978 0.461 Sensitivity analyses
Weibull 1.001 0.458
gamma 1.000 0.463 The fitted model may be sensitive to either parameters or
the data. We will investigate the latter by Monte Carlo
simulation and refitting. The Heimaey data consist of 36
Whether the residuals are incompatible with the model can tephra thicknesses in the fall direction. We randomly delete
be evaluated by obtaining confidence bounds via Monte N = 3, 6, 12, or 24 of these measurements and refit the
Carlo simulation. This uses repeated random sampling from best model (5) with Weibull error distribution. The results
the model with the estimated parameters to simulate pos- are shown in Table 3. We see that the estimates are con-
sible thicknesses consistent with the model at each loca- sistent in that the means vary little for N ≤ 12, and that
tion. Figure 4 shows the ratio between the simulated and the variability increases with decreasing amounts of data (N
observed thicknesses. A perfect fit is given by the horizon- increasing). Deleting two thirds of the data (N = 24) is
tal line, which can be compared with the 90 % point-wise too much for the model, inducing a tendency to find more
coverage band. For the model (6) without wind, there is eccentric (higher βU ) lobes in variable directions. It appears
a systematic trend with the ratio decreasing with observed that a number of the simulations retain locations suggest-
thickness. This indicates that small thicknesses are being ing different fall directions, or possibly more than one such
overfitted, and large thicknesses are underfitted. On the direction, with consequent instability in the other (coupled)
other hand, the baseline model (5) does not exhibit this parameters. The error distributions in the parameters are rea-
behavior; it also shows generally tighter bounds. All three sonably symmetric, apart from α, as would be expected. We

1 1 1
10 10 10
Fitted/Observed

Fitted/Observed
Fitted/Observed

0 0 0
10 10 10

−1 −1 −1
10 10 10

1 2 1 2 1 2
10 10 10 10 10 10
Observed (cm) Observed (cm) Observed (cm)

1 1 1
10 10 10
Fitted/Observed

Fitted/Observed

Fitted/Observed

0 0 0
10 10 10

−1 −1 −1
10 10 10

1 2 1 2 1 2
10 10 10 10 10 10
Observed (cm) Observed (cm) Observed (cm)

Fig. 4 Monte Carlo residual bounds for tephra thickness. The top row is the model (6) without wind; bottom row, the model (5) with wind. The
columns are for the lognormal, Weibull, and gamma error distributions reading left to right. Medians are shown by the dashed lines, 90 % bounds
by the dotted lines
738, Page 8 of 14 Bull Volcanol (2013) 75:738

Table 3 Parameter estimates

for the sensitivity analysis on No. of data Parameter
Heimaey data. Errors are 1σ
deleted, N γ βU φ α

3 71.3 ± 2.5 1.48 ± 0.07 125.0 ± 4.3 1.88 ± 0.08

6 71.8 ± 4.6 1.49 ± 0.12 124.8 ± 6.9 1.87 ± 0.14
12 72.2 ± 6.0 1.51 ± 0.19 125.9 ± 9.4 1.88 ± 0.19
24 81.4 ± 22.1 1.84 ± 1.04 124.0 ± 19.2 1.85 ± 0.44

can conclude that the model is robust to the amount of data, lobes, vents, and eruptions through a mixture framework. In
provided that the data represent the fall pattern correctly; as order to demonstrate this, we will briefly analyze the tephra
ever, the garbage-in-garbage-out principle applies. dispersal from the Ukinrek Maars eruption of March and
The sensitivity to the individual parameters is investi- April 1977.
gated following an idea analogous to that of Scollo et al. We will suppose that we have m vents, and that the ith
(2008). Starting at the optimum solution, we vary each vent has ni components (lobes). The model (4) to prescribe
parameter, one at a time, in the best model of (5) with the thickness at a given location then becomes
Weibull error distribution. Figure 5 shows the resulting
likelihoods, demonstrating that the parameter values are T (r1 , . . . , rm , U, θ1, . . . , θm )
robustly estimated.  m  ni
  −α
=γ Pi,j exp −{βU }i,j ri (1 − cos θi ) ri i,j , (8)
i=1 j =1
Multiple lobes and/or vents
where (ri , θi ) is the distance and azimuth (relative to the
The statistical model readily lends itself to investigation prevailing wind direction) to the location from the ith vent,
of prehistorical eruptions, allowing us to identify multiple and each component has its own α and βU parameters.

Fig. 5 Sensitivity analysis on

the estimated parameters. The −156 −156
maximum likelihood estimates −158 −158
are shown by the dashed lines
−160 −160
−162 −162
logL

logL

−164 −164
−166 −166
−168 −168
−170 −170
−172 −172
60 65 70 75 80 1 1.2 1.4 1.6 1.8
γ βU

−156 −156
−158 −158
−160 −160
−162 −162
logL

logL

−164 −164
−166 −166
−168 −168
−170 −170
−172 −172
110 120 130 140 1.6 1.8 2 2.2
φ α
Bull Volcanol (2013) 75:738 Page 9 of 14, 738

There is only one γ as the differing sizes of the com- (Kienle et al. 1980; Self et al. 1980). So, in mixture model
ponents are identified through the presence of the mixing terms, the source and direction of the observed components
distribution Pi,j , where i j Pi,j = 1. This is equiva- are
lent to summing over components with γi,j = γ Pi,j . The
West Maar
models can then be fitted to the data using maximum like-
lihood, in a minor elaboration of the technique described W1: SE (fall)
above. W2: SW-WSW (surge plus fall)
The tephra thickness data shown in Fig. 6 are taken from
Self et al. (1980). Note that we have deleted one thick- East Maar
ness measurement of 37 cm approximately 1 km southeast E1: NNE (fall 2 April)
of the East Maar, as it is the only measurement from that E2: NW-NNW (fall 1 April, 5 April, and 7–9 April)
described and partly sketched southerly tephra lobe in Self E3: WNW (surges throughout eruption)
et al. (1980). It is not possible to fit a component with E4: E (7–9 April)
four parameters to a single observation. High-level winds E5: SSE (7–9 April)
distributing fine tephra to distal areas were apparently not
related to the proximal tephra distribution which was con- As we have no useable data between SW and ESE, we
trolled only by the low-level winds under 2 km of altitude expect that we will be unable to detect components W1, W2,
(Kienle et al. 1980). The West Maar erupted first (30–31 and E5. Hence, we will limit ourselves to checking for, at
March 1977), with two tephra fall units directed to the SE most, four components.
and SW, overlain by surge deposits, which are thickest on Recall that the number of parameters in the distribution
the WSW edge of the crater (Kienle et al. 1980). The map is 4NC , where NC is the number of components. This is
of the deposits is in Self et al. (1980), which includes a too many to estimate from the 24 data available (Fig. 6).
slightly different chronology to that in Kienle et al. (1980). Hence, we will adopt the earlier idea of fitting a surface to
The eruption of the East Maar probably began on 1 April the observed data using a triangulated C1 -continuous inter-
(Self et al. 1980), with tephra fallout to the NNE and NNW, polating surface, overlaying a grid (at a spacing of 0.25 km)
with further eruptions on 5 April dispersing tephra to the and using the interpolated thicknesses at the grid points
NW and N. Later, Strombolian phases of eruptions from 7– (see Fig. 6) as our fitting data. Note that we restrict the
9 April produced tephra fall lobes to the NW, E, and SSE grid within the convex hull of the observed data to avoid

Fig. 6 Ukinrek Maars:

Observed thickness (in N
centimeter) and locations are in
large type. Interpolated
thicknesses (see text) are in
small type. The vents are
indicated
2
5
4 5 5 5 5 4 4 4 4
2.5 5 Becharof Lake
3 3 4 4 5 5 5 5 5 6 6 6 5 3
57°51’N
4 4 3 3 4 4 5 5 5 5 6 6 7 8 8 7 5 3
5
4 4 3 4 5 6 6 6 6 6 6 8 10 11 9 7 4 2
6 10
4 4 3
3 3 6 8 9 9 9 8 6 11 13 14 11 9 6 4
4 8 6 17
4 4 5 7 9 12 14 15 14 14 14 17 18 17 14 11 8 5 3
10 20
4 4 6 8 13 18 21 23 23 19 24 25 22 21 17 13 9 6 4 2
20 22 2
4 5 6 9 17 25 20 15 10 6 4 3

3 4 5 15 31 24 17 11 6 3
3
7
57°50’N 36 26 17 10 4

West East 37 24 14 5
Maar Maar 16 6

4
0 0.6km

156°33’W 156°30’W
738, Page 10 of 14 Bull Volcanol (2013) 75:738

extrapolation. This results in c. 150 data, enough to fit the Discussion

required models.
For the purposes of illustration, we will only present the We have considered only the tephra thickness, as this is the
Weibull error distribution results. The lognormal has already most common directly measured quantity. In many respects,
been shown to fit poorly, and while the gamma distribution particularly the risk to built structures, the tephra loading
is a viable alternative for one or two components, it does can be a more important measure. Often, thickness is con-
increasingly poorly relative to the Weibull as the number of verted into loading by assuming an average particle bulk
components increases. This is understandable, as the thin- density and packing density; measured values are often in
ner tail of the Weibull better describes the reduced sampling the range of 1,000 to 1,400 kg/m3 (Cronin et al. 1998). It
variation as excess thicknesses are ascribed to multiple com- is possible that a method similar to that above could be
ponents. Hence, the preferred model, as identified by AIC, developed to model tephra loading directly, provided that an
will give the number of components (on each Maar). As the attenuation model can be formulated.
East Maar is approximately 16 times the volume of the West Gonzalez-Mellado and De la Cruz-Reyna (2010) derived
Maar (Kienle et al. 1980), we have required at least one a number of empirical relationships based on fitting their
component on the East Maar. model (4) to 14 well-documented eruptions. In particular,
The results are shown in Table 4. The best model is they found that
clearly that with one component on West Maar (S), and
α = 2.535 − 0.051H, (9)
three on East Maar (E, NE and NNW), the last three of
which correspond quite closely to the sought-after E4, E1, giving an estimate of H , the height (in kilometer) of the
and E2, respectively. The West Maar component appears eruption column. It should be borne in mind that both α
to have (correctly) identified the aggregation of W1 and and H have uncertainties associated with them besides those
W2 from the scattering opposite to the dispersal axis. This allowed for in the regression analysis, and so this relation is
will have been overlaid with E3. The remaining parameter an approximation at best. Of course, the attenuation parame-
estimates are given in Table 5, although we will restrict fur- ter α is also a function of total grain-size distribution, which
ther comment to the preferred 1 West/3 East components is linked to column height. A second relationship between
model, which is illustrated in Fig. 7. All the models assign β and the column height,
a minimum weighting of 65 % to the East Maar, while the
1 114.407 − 4.189H H < 15.5,
preferred model assigns 80 % weight to the East Maar. This = (10)
2β −770.17 + 52.822H H ≥ 15.5,
volume apportioning is consistent with the relative sizes of
the two maars and with the fact that the bulk of activity which might otherwise allow us to separate the wind speed
was observed from the latterly erupting East maar (Kienle U from βU , appears to be an artifact of the inclusion of
et al. 1980). The decay with distance (α) and the eccen- the Pinatubo-2 eruption, with its column height of 43 km.
tricity (βU ) indicate a wide proximal fall from the West Without this point, a quadratic regression has an adjusted
maar, while the East maar contributes two narrower falls R 2 of zero, indicating no relationship.
(E4 and E1) and a smaller wide deposit slowly thinning For Heimaey, calculating H from the estimated α via (9),
with distance (E2). Again, these are in excellent agreement using the example of the Weibull error distribution and base-
with observations (Kienle et al. 1980, Self et al. 1980). The line model (5) which had the best AIC, yields an estimated
residual plot for the best fitting model is shown in Fig. 8; column height of H = 2.44 km, in line with the observed 2–
there appear to be no systematic patterns. 3 km (Wilson et al. 1978). We note that, in this instance, the

Table 4 Ukinrek models,

AICs, and fall directions Model (no. of components) AIC γ Directions (deg. anticlockwise from E)

West Maar East Maar West Maar East Maar

0 2 570.0 39.5 216, 39

1 1 538.9 44.8 90 30
0 3 514.5 78.7 27, 158, 96
1 2 498.9 67.0 94 12, 52
2 1 516.2 66.0 89, 42 18
1 3 468.9 74.0 264 12, 52, 113
2 2 477.0 78.3 297, 98 12, 52
0 4 484.9 78.7 15, 157, 99, 55
Bull Volcanol (2013) 75:738 Page 11 of 14, 738

Table 5 Ukinrek models, estimated parameters

Model (no. of components) P βU α

West Maar East Maar West Maar East Maar West Maar East Maar West Maar East Maar

0 2 0.64, 0.36 0.1, 3.7 1.9, 2.3

1 1 0.28 0.72 0.3 1.8 1.5 2.5
0 3 0.46, 0.46, 0.08 1.9, 4.6, 0.9 2.6, 2.9, 0.4
1 2 0.21 0.56, 0.23 0.4 4.8, 5.2 1.6 3.1, 1.7
2 1 0.20, 0.15 0.65 0.4, 10.5 4.0 1.5, 1.0 3.2
1 3 0.20 0.51, 0.21, 0.07 0.2 5.0, 6.1, 0.6 3.2 3.1, 1.7, 0.4
2 2 0.23, 0.06 0.49, 0.21 0.3, 0.6 5.3, 5.6 2.9, 0.3 3.2, 1.7
0 4 0.42, 0.35, 0.13, 0.10 3.4, 4.5, 0.7, 9.9 2.7, 2.7, 1.1, 1.3

particle density is near-uniform across the deposit (mean ± height. The observed column heights (Kienle et al. 1980) of
standard deviation of 2.08 ± 0.29 g/cm3 from 99 measure- 6.5 km (West Maar) and 3.5 km (East Maar) would require
ments). Hence, the isopachs are a constant function of the α’s of 2.2 and 2.4, respectively. We note that the less com-
isopleths which are usually used to estimate column height. plex models in Table 5 have α’s much closer to these ideal
As the Heimaey eruption was observed, we have the aver- values. Thus, the empirical relation (9) appears more likely
age wind velocity to the NW, which seems to have been to be valid where the deposit is a single lobe correspond-
on the order of 50 km/h (Self et al. 1974). From the esti- ing to the maximum column height, which are the condi-
mated βU = 1.48, this gives a value of β = 0.03 and, thus, tions under which it was derived by Gonzalez-Mellado and
an “effective diffusion coefficient” (Gonzalez-Mellado and De la Cruz-Reyna (2010). The higher βU values observed
De la Cruz-Reyna 2010) D = 17 km2 /h = 4,700 m2 /s, for some of the components in Table 5 may reflect the gen-
which appears reasonable. erally higher wind velocity environment at Ukinrek, wind
The relation (9) appears to break down for the Ukinrek shear causing bent-over plumes, or directional control in the
model with α’s greater than 2.535 (negative column height), individual eruptions.
or small enough (α = 0.3) to correspond to a 44-km column The diffusion coefficient is also an empirical para-
meter in such “semianalytical” models such as TEPHRA2
N WM1 and HAZMAP, “describing complex plume and atmos-
Becharof Lake
EM1 pheric processes not captured in the physical model”
EM2
EM3 (Volentik et al. 2010). Hence, the difference in “physical
reality” between physical and statistical models is one of
degree rather than a hard boundary. The advantage of the
2
5
5 N
57°51’N 2.5
4 5
6 10
4 3 8 6
17
10 0.48
20 1.23
20 22
2 Becharof Lake
15 57°51’N 0.70 1.13
7 3 1.39
3 East 0.99
57°50’N
West Maar 37 0.78 0.99
Maar 0.58 1.10
1.00 0.62 1.27
4 1.19 1.05
0.86 0.77 0.86
1.13 0.82
0.58 0.61
57°50’N 0.75
0 0.6km West East
Maar Maar
0.69
156°33’W 156°30’W
0 0.6km

Fig. 7 Ukinrek Maars: The 3-cm isocontours for each component 156°33’W 156°30’W
are shown. WM1 is the West Maar component, while the East Maar
components are denoted EM1–EM3, as detailed in Tables 4 and 5. Fig. 8 Ukinrek: Relative errors at measurement locations for the 1
Observed thickness (in centimeter) and locations are also shown. The West/3 East components model. The text is scaled by size of error. The
vents are indicated vents are indicated
738, Page 12 of 14 Bull Volcanol (2013) 75:738

statistical model is that there is an objective, consistent Applied to the 1973 Heimaey eruption, we were able to
method of parameter estimation through maximum likeli- decisively reject the model without wind effects. The wind
hood. The advantage of the physical model is that “feasible,” direction obtained with the other model corresponds well
but perhaps not consistent, estimates of column height and with the mean wind direction from the onshore wind that
mass eruption rate can be obtained, provided that grain-size deposited the measured tephra thicknesses. The estimated
information and/or eruption duration is also available. attenuation parameter indicated a mean column height equal
Heimaey and Ukinrek are typical examples of basaltic to that observed for the eruption (Wilson et al. 1978).
systems, scaled appropriately for the typical eruption Elaborating the model further in a mixture framework
types occurring in monogenetic volcanic fields. Hence, the enables the model to fit multiple lobes and/or vents. Applied
lessons from this study can be carried over to reconstruc- to the 1977 Ukinrek Maars eruption, it was able to identify
tion of parameters of tephra falls in prehistoric eruptions lobes in the correct direction from each vent. This has obvi-
in such fields, after allowance is made for tephra set- ous utility in studying unobserved eruptions with multiple
tling and compaction. The successful inferences made from vents.
the very sparse (and partially incomplete) thickness data
for Ukinrek is particularly encouraging, showing that this Acknowledgments The first author (E.K.) is supported by the New
approach is applicable to the typical quality of data obtain- Zealand Earthquake Commission and GNS Science. M.B. and S.C. are
able through paleo-volcanology studies. In particular, the supported by the New Zealand Natural Hazards Research Platform,
statistical model for tephra variation provides a platform project “Living with Volcanic Risk.” We thank Christina Magill and
an anonymous reviewer for improvements suggested to the original
for maximum likelihood estimation and, hence, the use manuscript.
of penalized fitting criteria such as the AIC, which per-
mits consistent estimation of the number of lobes. While
multiple tephra lobes can be inferred in complex multi- Appendix
event eruption episodes using this method, the order of
these cannot be extracted without additional stratigraphic Here, we derive the likelihood formulae for the various com-
information. On the other hand, we can estimate erup- binations of model and error distribution. Let T denote the
tive parameters including volume and identify statisti- tephra thickness obtained from (5) or (6).
cally the likely sources among multiple eruptive centers
(cf. Bebbington and Cronin 2011) and, thus, incorporate this Lognormal distribution
into hazard estimates.
In estimating the actual volume, the fact that the power Assume T ∼ lognormal (μN , σN ) where μN is the location
law produces an infinite thickness at r = 0 is an issue. parameter and σN is the scale parameter of the conjugate
To get around this, we could possibly follow the lead of normal distribution. Thus, μLN = exp μN + σN2 /2 and
Rhoades et al. (2002) and replace r −α by (r + δγ )−α ,  
which is dimensionally correct, and gives a finite thickness σLN = exp σN2 − 1 exp 2μN + σN2 are the lognor-
at r = 0, at the price of adding one more parameter to the mal parameters by definition. Then, μLNi = Ti (ri , ξi ) for
estimation problem. However, it is easier to simply exclude the ith observation, and so μNi = log Ti − σN2 /2. Note that
the area of the crater from calculations following, in effect, the location parameter is no longer a constant, as it varies
the suggestion of Bonadonna and Houghton (2005). Vol- for each observation.
umes can most easily be calculated by numerical integration The loglikelihood function for the complete sample i =
of the Eq. 4, using the estimated parameters. 1, 2, . . . , n is then obtained from (1) as follows:

n    n
log L = − log 2πσN2 − log Ti
Conclusion 2
i=1
n 2
i=1 log Ti − μNi
We have shown that a semiempirical model of tephra depo- − , (11)
2σN2
sition can be combined with an error distribution to produce
a statistical model capable of being fitted to actual measure- and the coefficient of variation is
ments rather than isopachs constructed from these measure- σLN  
ments. This provides a ready-made inversion formula for CV = = exp σN2 − 1 exp 2μNi + σN2 /
μLN
the volume of the eruption and the average dispersal axis.  
In addition, it can be used to forecast the distribution of exp μNi + σN2 /2
tephra at locations without data, and it provides an objective
measure of goodness of fit to the data. = exp σN2 − 1. (12)
Bull Volcanol (2013) 75:738 Page 13 of 14, 738

Weibull distribution Baxter PJ, Ing R, Falk H, French J, Stein GF, Bernstein RS, Merchant
JA, Allard J (1981) Mount St. Helens eruptions, May 18 to June
12, 1980: an overview of the acute health impact. J Am Med Assoc
Assume T ∼ Weibull (η, ν) where η is the shape parameter
246:2585–2589
and ν is  the scale parameter. Then μW = ν (1 + 1/η) and Bebbington M, Cronin SJ (2011) Spatio-temporal hazard estimation in
σW = ν  (1 + 2/η) −  2 (1 + 1/η) by definition. Thus, the Auckland Volcanic field, New Zealand, with a new event-order
μWi = Ti (ri , ξi ) for the ith observation, and νi = Ti / (1+ model. Bull Volcanol 73:55–72
1/η). Again, the scale parameter is no longer a constant, as Bebbington M, Cronin S, Chapman I, Turner M (2008) Quantifying
volcanic ash fall hazard to electricity infrastructure. J Volcanol
it varies for each observation. Geotherm Res 177:1055–1062
The loglikelihood function for the complete sample i = Bonadonna C, Houghton BF (2005) Total grain-size and volume of
1, 2, . . . , n is then obtained from (2) as follows: tephra-fall deposits. Bull Volcanol 67:441–456
Bonadonna C, Costa A (2012) Estimating the volume of tephra
deposits: a new simple strategy. Geology 40:415–418
n
 n
 Bonadonna C, Ernst GGJ, Sparks RSJ (1998) Thickness variations and
log L = n log η + (η − 1) log Ti − η log νi volume estimates of tephra fall deposits: the importance of particle
i=1 i=1 Reynolds number. J Volcanol Geotherm Res 81:173–187
n  η
 Bonadonna C, Connor CB, Houghton BF, Connor L, Byrne M, Laing
Ti
− , (13) A, Hincks T (2005) Probabilistic modeling of tephra dispersion:
νi hazard assessment of a multi-phase eruption at Tarawera, New
i=1
Zealand. J Geophys Res 110:B03203
and the coefficient of variation is Bonasia R, Macedonio G, Costa A, Mele D, Sulpizio R (2010) Numer-
σW ical inversion and analysis of tephra fallout deposits from the 472
CV = = νi  (1 + 2/η) −  2 (1 + 1/η)/ AD sub-Plinian eruption at Vesuvius (Italy) through a new best-fit
μW procedure. J Volcanol Geotherm Res 189:238–246
Carey S, Sparks RSJ (1986) Quantitative models of the fall and dis-
(νi (1 + 1/η))
persal of tephra from volcanic eruption columns. Bull Volcanol
 48:109–125
(1 + 2/η) Connor CB, Hill BE, Winfrey B, Franklin NM, LaFemina PC (2001)
= − 1. (14)
 2 (1 + 1/η) Estimation of volcanic hazards from tephra fallout. Nat Haz Rev
2:33–42
Connor LJ, Connor CB (2006) Inversion is the key to disper-
Gamma distribution sion: understanding eruption dynamics buy inverting tephra fall-
out. In: Mader H et al. (ed) Statistics in volcanology, vol 1:
Assume T ∼ gamma (κ, ω) where κ is the shape param- special publications of IAVCEI. Geological Society, London,
eter and ω is the scale parameter. Then μG = κω and pp 231–242
√ Costa A, Macedonio G, Folch A (2006) A three-dimensional Eulerian
σG = κω by definition. So, μGi = Ti (ri , ξi ) for the ith model for transport and deposition of volcanic ashes. Earth Planet
observation, and ωi = Ti /κ. Again, the scale parameter is Sci Lett 241:634–647
no longer a constant. Costa A, Dell’Erba F, Di Vito MA, Isaia R, Macedonio G, Orsi G,
The loglikelihood function for the complete sample i = Pfeiffer T (2009) Tephra fallout hazard assessment at the Campi
Flegrei caldera (Italy). Bull Volcanol 71:259–273
1, 2, . . . , n is then obtained from (3) as follows: Cronin SJ, Hedley MJ, Neall VE, Smith G (1998) Agronomic impact
n
 n
 of tephra fallout from 1995 and 1996 Ruapehu volcano eruptions,
log L = (κ − 1) log Ti − κ log ωi New Zealand. Environ Geol 34:21–30
i=1 i=1 Fierstein J, Nathenson M (1992) Another look at the calculation of
n fallout tephra volumes. Bull Volcanol 54:156–167
 Ti Gonzalez-Mellado AO, De la Cruz-Reyna S (2010) A simple semi-
−n log (κ) − , (15)
ωi empirical approach to model thickness of ash-deposits for differ-
i=1 ent eruption scenarios. Nat Haz Earth Syst Sci 10:2241–2257
and the coefficient of variation is Heiken G, Murphy M, Hackett W, Scott W (1995) Volcanic haz-
√ ards on energy infrastructure of the United States. United States
σG κωi 1 Department of Energy, Washington DC, LA-UR 95-1087
CV = = =√ . (16)
μG κωi κ Hurst AW, Turner R (1999) Performance of the program ASHFALL
for forecasting ashfall during the 1995 and 1996 eruptions of
Ruapehu volcano. NZ J Geol Geophys 42:615–622
Hurst T, Smith W (2004) A Monte Carlo methodology for modelling
References ashfall hazards. J Volcanol Geotherm Res 138:393–403
Johnston EN, Phillips JC, Bonadonna C, Watson IM (2012) Recon-
Akaike H (1977) On entropy maximization principle. In: structing the tephra dispersal pattern from the bronze age eruption
Krishnaiah PR (ed) Applications of statistics. North-Holland, of Santorini using an advection-diffusion model. Bull Volcanol
Amsterdam, pp 27–41 74:1485–1507
Barberi F, Macedonio G, Pareschi MT, Santacroce R (1990) Mapping Kienle J, Kyle PR, Self S, Motyka RJ, Lorenz V (1980) Ukinrek
the tephra fallout risk: an example from Vesuvius, Italy. Nature Maars, Alaska, I. April 1977 eruption sequence, petrology and
344:142–144 tectonic setting. J Volcanol Geotherm Res 7:11–37
738, Page 14 of 14 Bull Volcanol (2013) 75:738

Kratzmann DJ, Carey SN, Fero J, Scasso RA, Naranjo J-A (2010) Self S, Sparks RSJ, Booth B, Walker GPL (1974) The 1973 Heimaey
Simulations of tephra dispersal from the 1991 explosive eruptions Strombolian scoria deposit, Iceland. Geol Mag 111:539–548
of Hudson volcano, Chile. J Volcanol Geotherm Res 190:337– Self S, Kienle J, Huot J-P (1980) Ukinrek Maars, Alaska, II. Deposits
353 and formation of the 1977 craters. J Volcanol Geotherm Res 7:39–
Macedonio G, Pareschi MT, Santacroce R (1988) A numerical simu- 65
lation of the Plinian fall phase of 79 AD eruption of Vesuvius. J Sparks RSJ (1986) The dimension and dynamics of volcanic eruption
Geophys Res 93:14817–14827 columns. J Volcanol Geotherm Res 48:13–15
Miller TP, Casadevall TJ (2000) Volcanic ash hazards to aviation. Sparks RSJ, Bursik MI, Ablay GJ, Thomas RME, Carey SN (1992)
In: Sigurdsson H (ed) Encyclopedia of volcanoes. Academic, San Sedimentation of tephra by volcanic plumes. 2: controls on thick-
Diego, pp 915–930 ness and grain-size variations of tephra fall deposits. Bull Volcanol
Pfeiffer T, Costa A, Macedonio G (2005) A model for the numeri- 54:685–695
cal simulation of tephra fall deposits. J Volcanol Geotherm Res Sparks RSJ, Bursik M, Carey SN, Gilbert JS, Glaze LS, Sigurds-
140:273–294 son H, Woods AW (1997) Volcanic plumes. Wiley, Chichester,
Pyle DM (1989) The thickness, volume and grain size of tephra fall pp 574
deposits. Bull Volcanol 51:1–15 Stewart C, Johnston DM, Leonard G, Horwell C, Thordarsson T,
Pyle DM (2000) Sizes of volcanic eruptions. In: Sigurdsson H, Cronin SJ (2006) Contamination of water supplies by volcanic
et al. (eds) Encyclopedia of volcanoes. Academic, San Diego, ashfall: a literature review and simple impact modelling. J Vol-
pp 263–269 canol Geotherm Res 158:296–306
Rhoades DA, Dowrick DJ, Wilson CJN (2002) Volcanic haz- Sulpizio R (2005) Three empirical methods for the calculation of distal
ard in New Zealand: scaling and attenuation relations for volume of tephra-fall deposits. J Volcanol Geotherm Res 145(3–
tephra fall deposits from Taupo volcano. Nat Hazards 26:147– 4):315–33
174 Thorarinsson S, Steinthorsson S, Einarsson Th, Kristmannsdottir H,
Rose WI (1993) Comment on ‘another look at the calculation of fallout Oskarsson N (1973) The eruption on Heimaey, Iceland. Nature
tephra volumes’ by Judy Fierstein and Manuel Nathenson. Bull 241:372–375
Volcanol 55:372–374 Volentik ACM, Bonadonna C, Connor CB, Connor LJ, Rosi M (2010)
Scollo S, Del Carlo P, Coltelli M (2007) Tephra fallout of 2001 Etna Modeling tephra dispersal in absence of wind: insights from the
flank eruption: analysis of the deposit and plume dispersion. J climatic phase of the 2450 BP Plinian eruption of Pululagua
Volcanol Geotherm Res 160:147–164 volcano (Ecuador). J Volcanol Geotherm Res 193:117–136
Scollo S, Tarantola S, Bonadonna C, Coltelli M, Saltelli A (2008) Wilson L, Sparks RSJ, Huang TC, Watkins ND (1978) The control of
Sensitivity analysis and uncertainty estimation for tephra dispersal volcanic column heights by eruption energetics and dynamics. J
models. J Geophys Res 113:B06202 Geophys Res 83:1829–1836