Introduction

Large, caldera-forming eruptions pose a potentially devastating hazard to human populations and the environment. For example, the tephra generated by a magnitude (M) 7 or greater eruption can be dispersed across millions of square kilometres, posing immediate and long-term hazards (e.g. Tambora, 1815: Gertisser and Self 2015). In this paper, the term ‘tephra’ is used to refer only to fallout; deposits associated with pyroclastic density currents (ignimbrites) are discussed separately. Deducing the nature and scale of the hazards from large eruptions is challenging as their long repose period means eruptions of this scale have not been observed in modern times. Therefore, we are reliant on tephra deposits from prehistoric eruptions to inform our interpretations. Past eruptions are typically studied by mapping the tephra deposit in terms of thickness and grain size and then fitting the data using a range of models to approximate deposit extent, erupted volume, plume height and eruption magnitude (e.g. Carey and Sparks 1986; Pyle 1989, 2000; Bonadonna and Costa 2013). However, these methods often have limitations for large prehistoric eruptions that have dispersed fine ash over vast areas (e.g. Biass and Bonadonna 2011; Burden et al. 2013; Engwell et al. 2013).

Estimating the volume of prehistoric eruptions is key to assessing the magnitude-frequency relationships which may provide insight into the driving forces of volcanism and inform future risk mitigation strategies (Mason et al. 2004; Self 2006; Newhall et al. 2018; Rougier et al. 2018). However, eruptive volume estimates of unobserved eruptions require detailed interpretations of tephra deposits (Froggatt 1982; Pyle 1989; Fierstein and Nathenson 1992). This is challenging for large eruptions where the distal tephra is often deposited offshore (e.g. Tambora, 1815: Kandlbauer and Sparks 2014) limiting our interpretive power unless an extensive offshore record is available to supplement the terrestrial record (e.g. the Campanian tephra; Engwell et al. 2014). Owing to its widespread, on-land distal tephra deposit (Fig. 1), the ~ 7.7 ka eruption of Mount Mazama in Oregon, USA, is an ideal, large (M > 7; Crosweller et al. 2012; the LaMEVE database, VOLGRIPA; www.bgs.ac.uk/vogripa/view/controller.cfc?method=lameve) caldera forming eruption to study from this perspective (Williams 1942; Bacon 1983; Young 1990). Pacific deep-sea cores only record the Mazama tephra in turbidite deposits, not as a primary layer, supporting predominantly on-land tephra deposition (Fig. 2; Adams 1990). However, whilst the terrestrial Mazama tephra provides a unique opportunity to study a large prehistoric eruption, interpretation of the deposit is complicated by post-eruptive processes such as tephra remobilisation.

Remobilisation processes occur in the distal tephra deposits of large eruptions because of their wide spatial coverage and varied depositional environments. Tephra remobilisation is a secondary hazard acting over longer timescales than the primary hazards of the eruption. The remobilisation of tephra deposits is observed at a number of scales and can be initiated by fluvial, aeolian and gravitational processes. Documenting and studying sites where these processes are observed is key to enhancing our understanding of the processes and resulting hazards (Manville et al. 2000; Hadley et al. 2004; Lowe 2011; Gatti et al. 2013; Shapley and Finney 2015). For example, in the Younger Toba Tuff (YTT, M ~ 9.1) tephra, there is evidence of fluvial reworking and slumping at a number of distal sites in Malaysia and India > 350 km from source which altered the thickness and grain size (Gatti et al. 2011, 2013). Here we combine new field data with the extensive published data available for the distal Mazama tephra to better constrain the deposit extent and erupted volume. We identify and account for post-eruptive processes before producing new isopach maps and eruptive volume estimates.

Background

The Mazama tephra

Mount Mazama was a ~ 3700-m-high stratovolcano built primarily of andesite and dacite lavas from ~ 400 to 40 ka (Bacon and Lanphere 2006). During the Holocene, the climactic Mazama eruption formed modern day Crater Lake in Oregon, USA (Bacon 1983; Bacon and Druitt 1988; Bacon and Lanphere 2006). The eruption began with an increasingly violent Plinian phase before the onset of caldera collapse that produced large volume pyroclastic density currents (Williams 1942; Bacon 1983; Young 1990). The tephra from the climactic eruption is found as a visible layer across the western USA and Canada and as cryptotephra (a non-visible tephra layer) in Greenland (Hammer et al. 1980; Zdanowicz et al. 1999), Newfoundland (Pyne-O’Donnell et al. 2012), and the Great Lakes (Spano et al. 2017; Fig. 2c). Despite the extensive record of Mazama tephra and the importance of the tephra layers as an isochron, considerable uncertainty remains about the thickness distribution of the deposit.

Isopach maps are the most common way to represent the thickness of a tephra deposit and calculate erupted volume (Pyle 1989; Fierstein and Nathenson 1992; Bonadonna et al. 1998). However, constructing isopachs is challenging when there is insufficient and unreliable field data (Engwell et al. 2013). Therefore, for prehistoric eruptions such as the Mazama, the isopach maps contain a high degree of often unquantified uncertainty (Engwell et al. 2015), the effect of which is amplified when isopachs are hand drawn (Klawonn et al. 2014). The hand-drawn isopachs of the distal Mazama tephra in Fig. 1 demonstrate considerable variability reflecting the subjective interpretation by each of the authors. Additionally, all the isopach maps (Fig. 1) use a relatively small number of ash thickness measurements (less than 55) to reconstruct the deposit over an area that is > 1 million square kilometres, equivalent to ~ 1 measurement per 20,000 km2 (Lidstrom 1971; Matz 1987; Young 1990).

Fig. 1
figure 1

Compilation of Mazama isopach maps. a Location map. Red rectangle outlines the extent of maps (b–e). b Isopach map from Lidstrom (1971), tephra thicknesses shown in centimetres. c Isopach map from Matz (1987) with zone of ‘Distal Thickening’ in black. d Isopach map by Young (1990). For the published isopachs (b, c) some thicknesses are not clearly stated in the relevant thesis and have been marked as ‘ND’. e Minimum extent maps of the Mazama tephra from (1) Mullineaux (1974), (2) Westgate and Gorton (1981) and (3) Egan et al. (2015)

Full size image

The challenges of interpreting the distal Mazama tephra have long been recognised and as a result, the better-preserved deposits closer to source have been favoured for mapping and compositional studies (Williams 1942; Young 1990). In this study, we define distal as > 100 km from source and proximal as ≤ 100 km from the volcano. Mapping of the Plinian fall deposit by Williams (1942), Lidstrom (1971), and Young (1990) produced a high density of proximal tephra localities. Two distinct units were identified within the fall deposit, which is then capped by a fine ash layer, interpreted to be of co-ignimbrite origin (Fig. 3a; Young 1990). Here we integrate the proximal data collected by Young (1990) with previously reported and new distal localities to enhance our understanding of the primary tephra transport and deposition. Assimilating data sets from multiple disciplines, such as archaeology, palaeoclimatology and volcanology, requires an assessment of data quality and compatibility (see Lowe 2011). For example, a palaeoclimate study may only require that the tephra layer be identified for dating purposes and will therefore record limited physical information such as thickness (Heusser 1974). However, the observation that distal tephra exists at a given location, even without further details, is useful information.

The volume of the Mazama eruption

The challenge of quantifying the volume of the climactic Mazama eruption has been approached in a number of ways and has produced a range of volume estimates (Table 1). Some of the first estimates of dense-rock equivalent (DRE) volume were based on the volume of the caldera and reconstructed edifice. For example, Williams (1942) used the caldera depth to calculate an approximate caldera volume of ~ 50 km3. Adding this to the assumed ~ 30 km3 edifice, Williams estimated that > 70 km3 of material must have been ejected or removed by the climactic eruption. Williams and Goles (1968) and later Bacon (1983) revised those estimates to ~ 60 km3. However, more recently, Bacon and Lanphere (2006) reconstructed the Mount Mazama edifice and estimated a volume of ~ 112 km3. They attribute the discrepancy between the edifice and erupted volume to an incomplete knowledge of the collapse deposits in the caldera and the additional volume erupted by pre-climactic eruptions such as the Llao Rock and Cleetwood eruptions (Bacon and Lanphere 2006).

Table 1 Erupted volume estimates from the literature
Full size table

To estimate the bulk volume of the Mazama tephra deposit, multiple authors have used tephra thickness data (Williams 1942; Williams and Goles 1968; Lidstrom 1971; Young 1990; Table 1). Williams (1942) first estimate of ~ 15 km3 was based on the thickness of the proximal fall deposit (< 100 km from source). The estimate was revised by Williams and Goles (1968) who included three distal thickness measurements which increased the estimate by ~ 30 km3 (Table 1). The first estimate that used distal isopachs was ~ 120 km3 by Lidstrom (1971); Fig. 1b). Using different isopachs and incorporating the thick proximal fallout, Young (1990); Fig. 1d) also calculated a volume of ~ 120 km3. However, we would expect a discrepancy between these estimates because Lidstrom (1971) did not include the proximal fallout. This highlights the inconsistencies in volume estimates made from different isopach data (Fig. 1).

Converting the bulk erupted volume to a DRE volume requires assumptions of deposit density, which in turn depends on packing density and the vesicularity and composition of the clasts. The conversion to DRE volumes for deposits such as the Mazama tephra are complicated further by variations in the deposits. The Mazama tephra is composed of pumice, lithics and crystals, which have unique densities, morphologies and vesicularites, and the proportion of each component changes throughout the eruption (Lidstrom 1971; Bacon and Druitt 1988; Young 1990). The degree of deposit compaction also changes with grain size, depositional environment and time since deposition (Lidstrom 1971). These caveats for DRE calculations may explain some of the variability in DRE estimates found in Table 1.

Another source of variability in the bulk deposit volume estimates relates to the inclusion or exclusion of the ignimbrite volume (Table 1). Early work by Williams (1942) mapped the extent of the ignimbrite (‘glowing avalanches’) and estimated a bulk ignimbrite volume ~ 29 km3, which he included in his bulk erupted volume (total ~ 44 km3). More recent estimates of total erupted volume, by Young (1990), consider only the Plinian fall and distal tephra, or they are unclear as to whether the ignimbrite volume has been included (Bacon and Lanphere 2006; Johnston et al. 2014). It is likely that a significant proportion (> 10%) of the total erupted volume is contained in the ignimbrite, as flows reached over 70 km from source and are > 80 m thick around the caldera (Williams 1942). For this reason, the inclusion of this volume is important when calculating eruption magnitude.

Tephra remobilisation

A major challenge of using tephra layers to understand prehistoric eruptions is recognising modifications that occurred after the initial deposition. Processes that can remobilise primary tephra range from small scale bioturbation (Griggs et al. 2015), local slumping and grain flow down slopes (Boygle 1999) to large scale debris flows (Manville et al. 2000) and resuspension by surface winds (Wilson et al. 2011). Even if a tephra layer is not remobilised, varied depositional settings lead to differential preservation and compaction that result in variations in the tephra thickness preserved (Blong et al. 2017). The effectiveness of transferring a tephra deposit to the stratigraphic record will depend on the background sedimentation, vegetation cover and slope angle (Cutler et al. 2018; Dugmore et al. 2018). For this reason, peat bogs and lakes are often favoured depositional settings for studying tephra stratigraphy (Watson et al. 2016; McNamara et al. 2019). The Mazama tephra is recorded not only in peat bogs (Harward and Youngberg 1969) and lakes (Long et al. 2014) but also in dry land sections (Young 1990; this study) and aeolian sediments (Sweeney et al. 2005). This means that each tephra locality will have experienced different preservation mechanisms and compaction, as well as being exposed to a variety of remobilisation processes. Furthermore, the effectiveness of preservation and remobilisation processes vary with the physical tephra properties (thickness, grain size, density; Cutler et al. 2018; Dugmore et al. 2018).

Methods

We implemented a variety of methods to constrain the distribution of the Mazama tephra and estimate the total bulk erupted volume. We compiled a database of distal tephra localities containing data from multidisciplinary sources and new field sites. Each tephra locality was assessed for evidence of tephra mobilisation based on information in the original publication, field observation or what is implied for the location in relation to surrounding topography. Finally, we used thickness data from the newly optimised locality database to generate a reproducible isopach map using a cubic B-spline method. The new isopachs were then used for erupted volume calculations.

Compiling Mazama tephra localities

We have amalgamated and extended several existing databases of Mazama tephra localities to produce a single database that records > 250 occurrences of Mazama tephra (Lidstrom 1971; Matz 1987; Young 1990; Hallett et al. 1997; Egan et al. 2015; Jensen et al. 2019; J A Westgate, pers. commun.). The database records physical information about the tephra layer at each locality, which complements that of Egan et al. (2015), whose focus was radiocarbon dating. The information used was gathered from scientific literature, Masters and PhD theses, geological maps and field guides, as well as new field sites outlined in this study. The coordinates, locality name, reference(s), sampling method and research discipline of the original paper are provided for each locality. Where available, the thickness, grain size and the degree of remobilisation are reported with the accompanying metadata. If separate studies report different thicknesses at the same locality, an average thickness is used. Similarly, if an author ‘corrected’ the thickness without explaining the rationale or method, the original thickness is used rather than the corrected value (Matz 1987). Key field observations are included for terrestrial sites and for lake core localities, the depth of the tephra in the core is reported. If multiple cores were taken from the same site, only one locality is entered into the database, but the number of cores at each locality and any variability in thickness and grain size is documented in the metadata.

New field localities

This study includes new observations from ten field localities (all > 100 km from source; Supplementary S1). The Mount Bachelor locality (site 46; Fig. 2b) is approximately equivalent to a site mapped by Young (1990). The site required a tephra pit 1.5 m deep to expose the Mazama tephra stratigraphy (Fig. 3a). The total thickness and thickness of different horizons were measured, and samples were taken in 5 cm intervals from the top down. The stratigraphy was split according to the units in Young (1990) with a lower and upper pumice unit and final fine co-ignimbrite ash.

Fig. 2
figure 2

Distal Mazama tephra sites compiled for this study. a Map of the western USA and Canada with sites of Mazama tephra > 100 km from source. b Inset map with field sites from this study highlighted in white symbols with the site number adjacent. c Map of the USA, Canada and Greenland with sites of Mazama cryptotephra (Hammer et al. 1980; Zdanowicz et al. 1999; Pyne-O’Donnell et al. 2012; Spano et al. 2017)

Full size image
Fig. 3
figure 3

Field photographs of selected sites described in this study. a Tephra pit at Mount Bachelor (site 46). Star symbols reflect sampling sites. The tephra sequence records the climatic pumice fall and an upper fine-grained unit which has been interpreted as fine fallout associated with the ignimbrite phase (co-ig ash). The top mixed layer was dominantly soil but contained occasional pumice clasts. b Mazama tephra reworked at the surface near Prineville reservoir (site 50). c A reversely graded, discontinuous layer of Mazama tephra near Mitchell, Oregon (site 55). d A 70-cm-thick layer of Mazama tephra at Juniper Canyon (site 80) on top of alluvial fan deposits. e Dune-like discontinuous deposit of Mazama tephra at Spring Gulch (site 81). f 3 m lens of remobilised Mazama tephra at Pole Bridge on top of alluvial deposits (site 69). g Zoom on remobilised Mazama at site 69 showing absence of sedimentary structures

Full size image

At the Prineville Reservoir locality (site 50; Fig. 2b), the Mazama tephra had been reworked at the surface (Fig. 3b). This area was targeted based on previous studies (Harward and Youngberg 1969; Lidstrom 1971) that recorded Mazama ash in the area but provided no measurements of thickness or grain size. The Mazama tephra at the Mitchell locality (site 55; Fig. 2b) outcrops as a discontinuous layer which is interpreted as evidence of remobilisation (Fig. 3c). However, the presence of reverse grading in places may reflect the original stratigraphy.

We also examined seven sites > 450 km from Crater Lake (Fig. 2b). Some localities were previously documented in field guides (Carson and Pogue 1996) and Quaternary studies of the region (Waitt 1980). The sample sites were typically from riverbank, valley cut and roadcut exposures that required little excavation (Fig. 3d–g). No stratigraphic horizons could be identified in these deposits; therefore, only the total thickness was measured and bulk samples taken.

The grain size distributions (GSDs) of the tephra from the new field sites were measured by dry sieving from – 3 φ to 3 φ (8 mm to 125 μm) in half φ intervals. The fine material was then measured using the Mastersizer 3000™ by Malvern Instruments Ltd. in the School of Geographical Sciences at the University of Bristol. We recombined the 3 φ sieve fraction with the fine material (< 125 μm) prior to laser diffraction so the two methods of measuring grain size can be combined using the overlap in the 3 φ fraction (Eychenne et al. 2012; McNamara et al. 2018). For sites > 450 km from source, only laser diffraction was required due to the fine grain size. The grain size statistics—mode, median (Md), sorting (σ), skewness (Sk) and kurtosis (K)—were calculated using GRADISTAT (Folk and Ward 1957; Blott and Pye 2001).

Interpolated isopachs

We used a method of spline interpolation to interpret the thickness data collated in the database. Following the method outlined in Engwell et al. (2015), a cubic B-spline under tension was used to fit a model surface to the log thickness data (Inoue 1986; Supplementary S3). The FORTRAN code, developed by Inoue (1986) and modified by Engwell et al. (2015), outputs a gridded dataset across a specified x-y domain. This dataset was then processed in R using the ‘raster’, ‘ggplot2’ and ‘rgal’ packages to produce contoured plots of thickness (isopach maps).

Erupted volume calculations

From the spline fitted isopachs, the volume of distal Mazama tephra deposit was determined by fitting an exponential thickness decay to a plot of thickness versus square root isopach area (Pyle 1989). The same method was applied to previously published isopachs for comparison (Fig. 1).

The resulting log thickness versus square root area data was fit using ‘AshCalc’, a python tool for calculating erupted volumes (Daggitt et al. 2014). The exponential model assumes that tephra thickness can be described as a simple exponential decay away from the point of maximum thickness (Pyle 1989). For large and/or complex deposits, this exponential relationship breaks down and typically log thickness versus square root area plots are better described by multi-segment exponential decays (Pyle 1989; Fierstein and Nathenson 1992). In this work, we combine the distal isopachs with the spline isopachs of the proximal data (Young 1990) and fit multi-segment exponential decays using ‘AshCalc’ (Daggitt et al. 2014).

Results

The information amassed for each Mazama tephra locality and new field site is available in the database of tephra localities (Supplementary S1). The map in Fig. 2 shows all distal (> 100 km) localities in the database. Each thickness measurement has been classed as either ‘primary’, ‘secondary’ or ‘mixed’ based on the information available in the original reference, field observations and drainage analysis. Where there was no description or insufficient data, the locality was left as ‘unclassified’.

Identifying remobilisation

Thicknesses were classified as ‘primary’ where the measured thickness likely records the initial tephra thickness immediately after deposition. For example, site 46 (Fig. 3a) has been confidently classified as primary due to the well-preserved internal stratigraphy and topsoil development. Another primary site visited, site 73, preserves 30 cm of fine Mazama tephra with a sharp basal contact and well-developed topsoil. The locality is situated on a major drainage divide between the Snake River to the east and the John Day and Columbia Rivers to the west; therefore, we infer that the deposit at this site provides a minimum estimate of the primary fall deposit thickness, as no tephra has been remobilised from upslope (see Supplementary S2).

Thicknesses are considered ‘secondary’ if there is evidence for remobilisation such as slumping or bioturbation. An example from the literature is the 30 cm of Mazama tephra at site 174, Rockslide Lake in British Columbia, which contains multiple layers of Mazama tephra and has surrounding radiocarbon dates considerably younger than the eruptive date (Foit et al. 2004). A number of the sites are classed as secondary based on their field characteristics (Fig. 3b–g). Observations at site 80 (Fig. 3d) suggest that the 70 cm of the Mazama tephra overlying alluvial deposits has been remobilised (Sweeney et al. 2005; this study). Similarly, site 69 (Fig. 3f) records > 3 m of overthickened Mazama tephra (Carson and Pogue 1996; this study). Analysis of the upstream drainage (Supplementary S2) shows that, unlike site 73, there is substantial upstream drainage area for sites 69 and 80 which likely explains the overthickening.

‘Mixed’ thickness localities include both primary and secondary ash deposits. For example, the Mazama tephra at site 230 along the South Saskatchewan River has a sharp basal contact but is overlain by remobilised tephra that grades into the overlying sediment, which was included in the thickness reported (David 1970). Some sites where no thickness has been recorded were still classified based on the description of the stratigraphy.

Plotting log thickness against distance from source (Fig. 4) highlights the anomalous thickness values that have been classed as secondary and mixed. Typically, non-primary deposits tend to be thicker than the primary deposits at the same distance from source. The scatter in the thickness measurements with distance also reflects the scarcity of measurements that are on-axis, particularly > 500 km from source.

Fig. 4
figure 4

Log tephra thickness (cm) versus true distance from source (km). The solid line is a two-segment exponential fit to all the primary data and the dashed line is a two-segment fit to primary data < 200 km off-axis

Full size image

The grain size distributions (GSDs) of the ash samples collected for this study (Fig. 5) show that all the ash analysed from sites > 450 km from source is finer than 1 mm. All the samples display a unimodal distribution with a mode between 38 and 54 μm (4.8–4.3 φ) and median (Md) grain size of 33–49 μm (4.9–4.4 φ). The sorting (σ) of all samples is ~ 1.5 and the distributions are fine(negatively) skewed. The coarsest sample, with a median grain size of 49 μm, is from site 101 where the tephra had been excavated by a badger, which likely caused the fines to be lofted away explaining the fines depletion. The GSD from site 73, a primary tephra locality, overlaps with all other (remobilised) deposits.

Fig. 5
figure 5

Grain size distribution of distal samples. Site 73 (red) is the only primary (not remobilised) sample

Full size image

Cubic B-spline isopachs

The cubic B-spline method was applied to the primary thickness data (n = 45, where n is the number of measurements; Fig. 6a). Visually, the isopachs generated for the modelled thickness surface are similar to the published isopachs in terms of the dominant dispersal direction (Figs. 1 and 6), as forced by the inferred upwind thicknesses (Supplementary S3). However, even with the (inferred) upwind thicknesses, the spline isopachs do not match the tightness of the published upwind isopachs. Based on the lack of real upwind thickness data (SE of Crater Lake), we cannot assess whether the spline or published isopachs best represent the decay of thickness with distance in the upwind direction. Another difference between the spline and published isopachs is the extent and value of the thickest isopach. Specifically, the published isopachs show 30 cm of tephra reaching ~ 500 km downwind (Fig. 1b–d), but this is not reflected in the cubic B-spline isopachs (Fig. 6a). The spline isopach that best matches the published 30 cm isopachs is 20 cm thick. Whilst there are some important differences between the published and spline isopachs, the overall shape agrees with Lidstrom (1971) and Young (1990). In particular, along the main dispersal axis the extent of the 5 cm isopach is in close agreement.

Fig. 6
figure 6

Cubic B-spline isopachs. The thickness measurements used to generate each isopach are shown as filled black circles, the inferred thicknesses are plotted as open circles and the isopach thickness is in cm. a Primary thickness data and spline isopachs. b All thickness data and (mixed) spline isopachs. c Spline isopachs of proximal Young (1990) data. In a and b, the domain over which the spline interpolation was performed is shown by the dashed box and the proximal domain is a filled white square

Full size image

The cubic B-spline method was also applied to all the thickness data (n = 138) recorded in the database (Fig. 6b), including thickness measurements that have been classed as secondary, mixed and unclassified. The resulting isopachs record a more complex and convolute distribution of tephra thickness. In particular, the spline isopachs highlight an area of overthickening in north-eastern Washington, which broadly overlaps with the area designated as ‘distal thickening’ by Matz (1987). Another key difference between the published and spline isopachs is the 1 cm isopach. All of the published isopach maps include a 1-cm limit. However, the spline interpolation method had insufficient data to draw a closed contour of 1 cm.

The cubic B-spline isopachs of the proximal data from Young (1990) are shown in Fig. 6c. The overall distribution of thickness is close to his hand-drawn isopachs, although, as with the distal spline isopachs, the upwind decay of thickness is much tighter in the hand-drawn isopachs. A south-east lobe in the cubic spline isopachs (Fig. 6c) corresponds to a mapped SE lobe in the lower pumice unit of the Plinian fall that was observed and interpreted by Williams (1942), Lidstrom (1971), and Young (1990) to correspond to a different wind direction and or plume height at the onset of the eruption.

The goodness of fit of the spline surface to the thickness data is output as a root mean square (RMS) residual (see Table 2). Although the RMS value quantifies the fit of the whole spline surface, we are interested in the difference between the modelled and measured thicknesses at individual localities. To assess this difference, we calculate a percentage error between the measured thickness and modelled thickness at each locality:

$$ \%\mathrm{error}=\frac{\left|\mathrm{modelled}\ \mathrm{thickness}-\mathrm{measured}\ \mathrm{thickness}\right|}{\mathrm{modelled}\ \mathrm{thickness}}\times 100\kern4.5em $$
(1)
Table 2 Calculated volumes from spline and published isopachs
Full size table

Importantly, in Eq. 1, we divide the absolute difference by the modelled thickness because we want to highlight localities where the measured value deviates strongly from the spline predicted value. We expect the spline surface to poorly match thick proximal sites due to the majority of the data fitted with the cubic B-splines being thinner distal data. For example, at a proximal locality (site 26), the measured tephra thickness is 233 cm, but the spline surface predicts 38 cm which we calculate as a 486% error. Whilst the fit appears poor to these few proximal sites, they were necessary to include in the fitting of the cubic B-spline model to locate the area of maximum thickness (Fig. 7a). For the mixed spline isopachs (Fig. 7b), the RMS residual and number of points with > 400% difference between the spline and measured thickness is greater than in Fig. 7a because the spline fitting has smoothed highly irregular data. Comparing the measured thicknesses of secondary and mixed sites to the primary spline surface (Fig. 7c) clearly highlights areas of overthickening (secondary and mixed deposits), where the percentage difference is typically > 400%. Comparing the unclassified thicknesses to the primary spline surface (Fig. 7d) shows that 43 unclassified thicknesses are within 100% of the value predicted by the spline surface.

Fig. 7
figure 7

Percentage difference between the actual thickness measurements and the spline interpolated thicknesses. Isopach thicknesses are 20, 10 and 5 cm. a Primary thickness only. b Mixed spline thicknesses. c Difference between mixed thicknesses and primary spline surface. d Difference between primary and unclassified thicknesses and primary spline surface

Full size image

Volume calculations

The volume of the distal Mazama tephra calculated from the individual isopach maps in Figs. 1 and 6 are shown in Table 2. For Fig. 6b, the total area enclosed by the thickest isopach is the summation of the areas enclosed by the two closed 20 cm contours. The volumes calculated from the primary and mixed spline isopachs are 134 km3 and 161 km3 respectively. The mixed spline isopachs thus suggest a ~ 30-km3 increase in volume compared to the isopachs drawn using only the primary data. The largest calculated volume is 242 km3, which includes the isopachs that show a zone of distal thickening and a large 1-cm isopach (Fig. 1c; Matz 1987). The smallest volume (131 km3) was calculated from the distal Young (1990) isopachs. Importantly, all of the volumes calculated using the distal data alone are minima, as they do not include the volume in the thick proximal fallout.

The volumes in Table 2 calculated using both the proximal (Young 1990) and distal data were fit using one-, two- and three-segment exponential fits, as shown in Fig. 8. Using two exponential segments produces a better fit to the proximal and medial data for all of the isopach combinations compared to the single segment exponential decay which underestimates the thickness close to the vent and overestimates intermediate thicknesses. We also tested three-segment exponential fits, but conclude that they ‘overfit’ the data, as in most cases there are only two points controlling each exponential segment.

Fig. 8
figure 8

Square root area versus log thickness plots for proximal and distal Mazama isopachs. a–e Dashed line shows single segment fit, thick line shows two-segment fit and thin line shows three-segment fit. Corresponding volume estimates can be found in Table 2. All proximal data corresponds to cubic B-spline fit of Young (1990) data (Fig. 6c). a Primary spline isopachs. b Mixed spline isopachs. c Lidstrom (1971). d Young (1990). e Matz (1987). f Two-segment exponential fits to log thickness versus true distance data from Fig. 4 overlain on two-segment fits from a to e

Full size image

Discussion

Sources of variability in tephra thickness measurements

The amalgamation of > 250 Mazama tephra occurrences highlights the wide range of disciplines that make use of the tephra layer from volcanology (Westgate and Dreimanis 1967; Mullineaux 1974) to archaeology (Kittleman 1973; McFarland 1989). The multidisciplinary use of the Mazama tephra has resulted in different methods of data collection and variable data quality. However, this alone cannot explain the extreme variations in tephra thickness data (Figs. 2 and 4). The large variability of tephra thickness across the deposit also reflects the extensive remobilisation of the distal Mazama tephra, as the material has been eroded and redeposited over the past ~ 7700 years. Despite the paucity of primary distal Mazama tephra localities, detailed reporting and publication of all tephra sites is important as remobilised localities provide a record of tephra deposition close to that location.

Direct field observations of remobilised deposits identify common features of non-primary Mazama tephra, such as discontinuous beds, lens and fan-shaped outcrops and the co-occurrence with alluvial deposits (Fig. 3). Where we observed remobilised Mazama tephra, however, there were no obvious sedimentary structures; this differs from the remobilised YTT deposits which contain fine scale cross-beds (Gatti et al. 2013). The remobilised YTT also exhibits different GSDs compared to the primary YTT (Gatti et al. 2011, 2013); in contrast, the GSD of the remobilised and primary Mazama tephra are nearly identical, at least ~ 150–500 km from source (Fig. 5). Therefore, without broad deposit context, distinguishing the remobilised Mazama tephra deposits from primary Mazama tephra is difficult. The lack of sedimentary structures, and seemingly unaltered GSD, suggests that the remobilised distal samples have not experienced significant contamination from background sediments, or hydrodynamic sorting, even when found on top of alluvial deposits, implicating water in the remobilisation process. This is in contrast to the prominent sedimentary structures in remobilised YTT that are clearly the result of transport by water (Gatti et al. 2013). We hypothesise that downslope remobilisation of Mazama tephra occurred en-masse with minimal water entrainment and shortly after the primary deposition, such that little foreign material was incorporated and no sedimentary structures were developed.

During the mid-Holocene (~ 8–5 ka) North America was slightly warmer and drier than present (Dean et al. 1996; Thompson et al. 2016) and reconstructions of the vegetation cover (Fig. 9; Adams et al. 1997) suggest that most of the Mazama tephra was deposited in a semi-desert to tundra dominated environment. This supports our interpretation that the downslope processes that led to the substantially overthickened deposits (Fig. 3) could have been relatively dry. However, pollen suggest that the Mazama eruption occurred in northern hemisphere autumn (Mehringer et al. 1977). This could have resulted in snow melt the following spring remobilising the deposits in some areas (Manville et al. 2000) and the unimodal grain size distribution of the distal tephra (Fig. 5) may be the reason for the absence of obvious sedimentary structures. As such, the processes that lead to extreme overthickening (3 m from ~ 30 cm primary fall at site 69) require further study.

Fig. 9
figure 9

Mazama tephra distribution in relation to palaeovegetation cover. Vegetation zone limits are from Adams et al. (1997)

Full size image

To avoid the challenges of redeposited terrestrial tephra deposits, numerous studies have used lake sediments to reconstruct volcanic histories as lake cores can provide excellent records of tephra stratigraphy (e.g. Lowe et al. 1980; Hopkins et al. 2015; Fontijn et al. 2016; McNamara et al. 2019). However, we find that lake cores, as with land sections, do not always provide reliable records of tephra thickness. Remobilised tephra in lake cores can be identified by diffuse upper or lower contacts with the background sediments, altered GSDs and anomalous radiocarbon dates surrounding the layer (Hopkins et al. 2015). For example, Swiftcurrent Lake (site 150) in Montana contains 48 cm of Mazama tephra but with up to 25% clastic contamination (MacGregor et al. 2011), suggesting the layer measured was a mixture of primary and reworked material. Unfortunately, many of the lake core localities lack detailed observations of the tephra layer making it difficult to confidently classify the tephra as remobilised or primary.

A number of factors can affect the thickness of tephra recorded in lake cores, including the vegetation and steepness of the surrounding hillslopes (Dugmore et al. 2018), catchment size (Supplementary S2), lake depth, aeolian transport, the position of major inlets relative to the coring sites and the background sedimentation (McNamara et al. 2018). Factors independent of the lake setting may also influence the preservation of tephra in lake cores. For example, the rate of tephra deposition and the quantity of accumulated tephra will influence how the tephra is transported from the lake surface to the lake bottom (Engwell et al. 2014; McNamara et al. 2019). The density and grain size of the tephra could be another external factor affecting tephra preservation in lakes. The distal Mazama tephra is fines dominated (~ 32–63 μm; Fig. 5) which may make it more susceptible to syn-deposition remobilisation as the slower settling velocities will mean it has a longer residence time in the water column. Furthermore, once it has reached the lake bottom, fine-grained tephra experience additional remobilisation by lake bottom currents and bioturbation (Dashtgard et al. 2008).

Tephra thicknesses recorded in both lake and terrestrial settings may also have been influenced by wind remobilisation. The modal grain size of the distal Mazama tephra (Fig. 5) falls within the range measured for resuspended modern ash in Iceland and Argentina (Liu et al. 2014; Panebianco et al. 2017) and is also similar to (wind-transported) loess (Sun et al. 2004). Semi-desert and tundra environments, which we infer a substantial proportion of the tephra was deposited in (Fig. 9), provide the ideal dry and windy conditions for tephra resuspension. Wind remobilisation could thin the deposit in some places but could also explain overthickening in other sites. For example, at Marys Pond, Montana (site 102; Foit et al. 1993), the small glacial lake must have had an input of Mazama tephra other than purely upstream to explain the accumulation of 9–90 cm of tephra > 700 km from source. Other windblown particulates, such as loess, are known to accumulate in lakes (Pye 1995); therefore, windblown tephra may help explain some overthickened lake deposits.

Understanding mechanisms that remobilise tephra after a large eruption are crucial for anticipating post-eruptive hazards. Ash resuspension, for example, has closed airports (Hadley et al. 2004) and could pose hazards to human health (Horwell and Baxter 2006; Thorsteinsson et al. 2012). Resuspension events are a persistent hazard in Alaska over 100 years after the magnitude 6.5 Novarupta (Katmai) eruption in 1912 (Hadley et al. 2004) and it is likely that the Mazama tephra (M > 7) posed resuspension hazards over similar time scales or longer. Archaeological records provide evidence for long term ash remobilisation, revealing that communities abandoned parts of the north-western Plains (modern Wyoming, Montana, North and South Dakota, southern Alberta Saskatchewan and Manitoba) following the Mazama eruption and only returned permanently after 500 years (Oetelaar and Beaudoin 2016).

Validity of isopachs for large eruptions

The published isopachs for the distal Mazama tephra (Fig. 1) show the challenges in creating isopachs for large and older tephra deposits. Our study shows that the isopach variability derives from overthickening, sparse coverage and different environments of deposition (Figs. 2, 3 and 4). This raises the question of whether isopach maps are an appropriate way to present thickness data for large deposits, or if the uncertainties and limitations invalidate this approach. The disadvantages of hand-drawn isopachs have been discussed (Klawonn et al. 2014; Engwell et al. 2015), whilst acknowledging the value of intrinsic knowledge of the authors, which may not be possible to quantify or transfer by means other than hand-drawn isopachs. However, it is still important to clearly document the data used to construct the isopachs in a best effort to transfer this knowledge. As Fig. 1 shows, often the data are poorly documented and are therefore unusable except to determine the areal extent of tephra coverage.

Alternatively, studies of large and poorly documented deposits may report thickness against true distance from source (e.g. Campanian, Engwell et al. 2014; Toba, Gatti and Oppenheimer 2012; Taupo, Matthews et al. 2012). The equivalent plot of the Mazama thickness data in Fig. 4 shows the value of this approach. Here the large scatter of thickness data at the same distance from source reflects both the distance from the main dispersal axis and the inclusion of non-primary thickness measurements. Using the data in this form for volume calculations comes with caveats (see below and Supplementary S4); however, qualitatively comparing the decay of thickness with distance between large tephra deposits provides a means of contrasting the magnitudes of different eruptions without the need to draw isopachs. We demonstrate this in Fig. 10, where we compare exponential fits to the distal thickness data for the Campanian (Engwell et al. 2014), Rangitawa (Taupo; Matthews et al. 2012), YTT (Matthews et al. 2012) and Mazama eruptions. Whilst the data exhibit substantial scatter due to the reasons previously discussed, the simple exponential fits highlight the different thinning trends of ‘super-eruptions’ (M > 8) and magnitude ~ 7 eruptions such as the Mazama and Campanian.

Fig. 10
figure 10

Comparing the exponential decay of distal tephra thickness with distance from large magnitude eruptions. Younger Toba Tuff (YTT), Campanian and Mazama magnitude values (M) are taken from the LaMEVE database, (VOLGRIPA; www.bgs.ac.uk/vogripa/view/controller.cfc?method=lameve) and Rangitawa M is assuming DRE > 1500 km3 (Matthews et al. 2012)

Full size image

Spline fitted isopachs

To overcome the disadvantages of hand-drawn isopachs, we applied a spline interpolation method (Engwell et al. 2015) to the distal and proximal thickness data (Fig. 6). This method generates reproducible isopachs of the distal Mazama tephra with reduced subjectivity. However, the spline fitting still requires user input, therefore necessitating prior knowledge of the deposit. For example, fitting the cubic B-spline to the raw primary thickness data required inferred upwind thicknesses to characterise the asymmetry of the deposit (Supplementary S3). This user input can only be avoided entirely when the thickness dataset is spatially dense, with zero values delimiting the edge of the deposit (e.g. Mount Saint Helens, 1980; Sarna-Wojcicki et al. 1981; Engwell et al. 2015). Fitting the cubic B-spline model also requires a choice of the fitting parameters: roughness, tension, area divisions and weightings (Supplementary S3). For this study, we followed the guidelines of Engwell et al. (2015) and found the values that best visually recreated the expected dispersal; however, this introduces another area of subjectivity.

Irrespective of the limitations of the spline method, there are benefits to using models to generate isopachs. Notably, the same isopachs can now be reproduced given the same thickness dataset and model parameters regardless of user bias or experience which cannot be said for hand-drawn isopachs (Klawonn et al. 2014). The spline fit does not recreate the exact fallout of the Mazama tephra, as it involves no physical processes relating to ash dispersal and transport, but it generates a ‘model’ tephra surface from which anomalous thicknesses and unusual deposition can be highlighted (Fig. 7). Figure 7c shows the correlation between the percentage difference and the thickness classification. Additionally, Fig. 7d shows that this approach could be used to make best estimates of otherwise ‘unclassified’ localities which could inform future field studies of primary tephra localities. Another benefit of using the spline fitting approach is the capacity to include non-primary thicknesses to explore the sensitivity of isopach construction to these measurements (Fig. 6b).

Estimating the erupted volume from large tephra deposits

The bulk erupted volumes of the distal Mazama tephra calculated from published and spline isopachs using an exponential fit range from 131 to 242 km3 (Table 2). The range in estimates arises from the inclusion of non-primary thicknesses (Fig. 6b) and convolute published isopachs (Fig. 1c; Matz 1987). Including non-primary, overthickened thickness measurements in isopach construction is useful for highlighting areas of tephra remobilisation. However, using isopachs that include overthickened data to estimate erupted volumes is problematic as this negates the assumption of deposit thinning with increasing distance from source, which is fundamental to the methods used to estimate volumes from isopach maps (Pyle 1989; Fierstein and Nathenson 1992). For this reason, direct integration of the volume below the spline surface may provide a better way to estimate volume from non-primary isopachs (Supplementary S4). However, the spline fit underestimates the thickness at source and the integration range is limited by the thinnest isopach; this method provides only volume minima. For these reasons, we do not use the convolute (Matz 1987; Fig. 1c) and mixed spline isopachs (Fig. 6b) to estimate the total erupted volume of the Mazama tephra. Whilst the range in estimates remains significant (distal only 131–153 km3; Table 2), large uncertainties on erupted volumes are expected for large eruptions.

Combining the primary distal spline isopachs (Fig. 6a) and the proximal spline isopachs (Fig. 6c) gives the best estimate of 148 km3 bulk erupted volume. When added to the ~ 29 km3 assumed to be in the ignimbrite (Williams 1942), we estimate a total erupted volume of ~ 176 km3 which is ~ 20 km3 greater than the estimates in Table 1. To convert this to a dense-rock equivalent (DRE) volume, we assign average densities to the three portions of the deposit and use 2200 kg m−3 as the bulk density of the magma (Lidstrom 1971; Young 1990). For the proximal deposit, we use an average deposit density of 500 kg m−3, taking into account average pumice densities (Young 1990) and packing density. We measured a bulk density of 700 kg m−3 for the distal ash, which is in agreement with the value used by Lidstrom (1971), Bacon (1983) and Young (1990), with the low density likely reflecting the micro-vesicularity of the fine ash. For the ignimbrite we use 1200 kg m−3, the higher density reflecting the large lithic component in the ignimbrite flows (Bacon 1983). Using this density distribution and the relative proportion of volume in the proximal (19 km3) and distal (129 km3) segments, we calculate a DRE volume of ~ 61 km3 and calculate an eruption magnitude of 7.1 (M = log10(DRE × 2200) − 7; Pyle 2000). This value is in agreement with DRE estimates by Lidstrom (1971) and Bacon (1983) and is less than the caldera and edifice volume estimate by Bacon and Lanphere (2006) allowing for the volume contained in collapse deposits and pre-caldera forming eruptions (Table 1).

Conclusions

Estimating the eruptive volume and magnitude of prehistoric large magnitude eruptions is vital to interpreting the hazards posed by these eruptions. However, the tephra deposits required to estimate parameters, such as eruptive volume, are often poorly preserved and exhibit low sampling density and diverse data quality, especially at large distances from source (> 100 km). These limitations are encountered when interpreting the distal Mazama tephra. We combined an extensive literature search, field observations and new methods of interpolation and data manipulation, to develop an extensive record of Mazama tephra occurrences (Supplementary S1). Using the compilation of primary tephra sites, we constructed new reproducible isopachs of the deposit and revised the bulk erupted volume estimate to ~ 176 km3, including the ignimbrite volume. The new isopach map (Fig. 9) provides a modelled thickness distribution which broadly agrees with past hand-drawn isopachs and are a resource that can be used to inform future field studies. Specifically, further work is needed to constrain the upwind deposit and the limits of the distal tephra, for example the limit of the 1 cm isopach which has only been approximated (Fig. 9). This will improve our volume estimates as well as our understanding of the spatial footprint of hazards from the primary tephra fallout during large caldera forming eruptions.

We use our compilation of Mazama tephra occurrences and isopach maps to highlight areas where the Mazama tephra has been remobilised and overthickened and show that incorporating non-primary thicknesses dramatically influences the isopachs and the bulk erupted volume estimates. However, we emphasise that the reporting of all tephra sites is important as they are still a record of deposition at a specific location, even if the thickness cannot be used in isopach construction. Furthermore, understanding the mechanisms and conditions that encourage tephra remobilisation is crucial for interpreting the post-eruptive hazards posed by a Mazama-like eruption as they have the potential to have long-lasting effects (100+years) over large areas (millions of km2).