No evidence for tephra in Greenland from the historic eruption of Vesuvius in 79 CE: implications for geochronology and paleoclimatology

Volcanic fallout in polar ice sheets provides important opportunities to date and correlate ice-core records as well as to investigate the environmental impacts of eruptions. Only the geochemical characterization of volcanic ash (tephra) embedded in the ice strata can confirm the source of the eruption, however, and is a requisite if historical eruption ages are to be used as valid chronological checks on annual ice layer counting. Here we report the investigation of ash particles in a Greenland ice core that are associated with a volcanic sulfuric acid layer previously attributed to the 79 CE eruption of Vesuvius. Major and trace element composition of the particles indicates that the tephra does not derive from Vesuvius but most likely originates from an unidentified eruption in the Aleutian arc. Using ash dispersal modeling, we find that only an eruption large enough to include stratospheric injection is likely to account for the sizable (24–85 µm) ash particles observed in the Greenland ice at this time. Despite its likely explosivity, this event does not appear to have triggered significant climate perturbations, unlike some other large extratropical eruptions. In light of a recent re-evaluation of the Greenland ice-core chronologies, our findings further challenge the previous assignation of this volcanic event to 79 CE. We highlight the need for the revised Common Era ice-core chronology to be formally accepted by the wider ice-core and climate modeling communities in order to ensure robust age linkages to precisely dated historical and paleoclimate proxy records.

1 Introduction

Volcanism is now widely accepted as the most significant natural driver of high-frequency (i.e., interannual to centennial timescales) climate variability (Crowley, 2000; Schurer et al., 2013). Polar ice cores provide long, continuous records of volcanism through their capture of volcanic aerosols, halogens and particles (Hammer et al., 1980; Zielinski et al., 1994). The precise dating of the volcanic signals proffered by the ice cores enables the climate and societal impacts of volcanic eruptions to be evaluated with respect to historical records and has contributed to the discernment of strong chronological correlations between past volcanism, climate perturbations and societal upheaval (Sigl et al., 2015; Büntgen et al., 2016; Toohey et al., 2016).

Volcanic sulfate (or sulfur) remains the primary proxy in the ice cores used for reconstructing past volcanism, as continuous records can be constructed using a continuous ice-core analytical system (McConnell et al., 2002, 2017) in which longitudinal samples are melted sequentially and the meltwater is analyzed for a broad range of elements and chemical species, including sulfur, in real time. Volcanic sulfates cannot, however, pinpoint the source of the eruption, a critical factor needed to determine the relationship between sulfate flux in the ice cores and atmospheric burden, thus also determining the climatic effectiveness of an eruption. In particular, stratospheric injection height of material is a key variable controlling the radiative forcing potential of an eruption (Marshall et al., 2019; Toohey et al., 2019). In the absence of documented evidence, eruption column height can be estimated from isopleth analysis of proximal tephra fallout (Carey and Sparks, 1986; Burden et al., 2011). Volcanic ash (fine-grained tephra, often referred to as cryptotephra when shard concentrations are too low to be visible to the naked eye) in ice cores can also enable the source volcano to be identified, thus improving the scope for modeling aerosol distribution and climate impact. The geochemical composition of tephra is a product of the magma composition from which it derives, which is determined by the local geology. As the magma source is subject to changes in supply, temperature and pressure, tephra chemistry may be unique to each eruption or phase of the eruption, but it frequently carries a geochemical signature in its major or trace element composition that enables it to be attributed to a specific source (Westgate and Gorton, 1981).

Figure 1 Key sites and ice-core records referred to in this paper: (a) location of Greenland ice-core drill sites and volcanic sources; (b) non-sea-salt sulfur (nssS) from four Greenland ice cores showing a prominent peak at 88 CE; (c) detail of nssS peak in the NEEM-2011-S1 ice core and position of the samples examined in this study (inset shows tephra shard found in sample QUB-1832).

Here we examine a tephra deposit associated with a large volcanic sulfate signal in Greenland ice previously attributed to the 79 CE eruption of Vesuvius (Fig. 1; Hammer et al., 2003; Vinther et al., 2006; Barbante et al., 2013). The event has no coinciding Antarctic signal, which would support the theory of a Northern Hemisphere extratropical source. The sulfate signal in Greenland is quite strong: the global radiative forcing from the eruption was estimated to rank 23rd of all eruptions of the past 2500 years (Sigl et al., 2015). Under the assumption of an extratropical source, the stratospheric sulfur injection of the eruption is estimated to be 10.6 Tg S (Toohey and Sigl, 2017), which would be the fifth largest injection from an extratropical NH eruption over the same time span, comparable to the upper range of estimates for the 1912 Katmai eruption (Stothers, 1996) and greater in sulfur injection than the 1991 Mt. Pinatubo eruption. The attribution of the event to Vesuvius was initially based upon the acid layer's apparent temporal coincidence with the 79 CE eruption (Zielinski, 1995) and later appeared to be confirmed by the major element analysis of diminutive tephra particles (2–5 µm particle size) in the GRIP ice core whose geochemistry was found to be consistent with Vesuvius tephra (Barbante et al., 2013).

The Greenland Ice Core Chronology 2005 (GICC05) provides a synchronized dating for three Greenland ice cores (NGRIP, GRIP and Dye-3) in which three volcanic marker events were believed to correspond to historic volcanic eruptions (i.e., Vesuvius 79 CE, Hekla 1104 CE and Öraefajökull 1362 CE). The Vesuvius and Öraefajökull horizons were both supported by unpublished tephra evidence when GICC05 was released (Vinther et al., 2006) and subsequently published (Coulter et al., 2012; Barbante et al., 2013), although attempts to isolate Hekla 1104 CE tephra have not, so far, been successful (Coulter et al., 2012; Plunkett, unpubl. data, 2013). These horizons were subsequently assigned an age uncertainty of ± 0 years (Vinther et al., 2006). This ice-core chronology has since been revised, however, following the recognition of distinct solar proton events at 774 and 993 CE that highlighted a temporal offset in GICC05, prompting the construction of the revised ice-core chronologies NS1-2011 (Sigl et al., 2015) and DRI_NGRIP2 (McConnell et al., 2018). The offset is further supported by an independent evaluation of radionuclide anomalies in Greenland ice cores and tree-ring records (Adolphi and Muscheler, 2016). As a result, the date of the sulfate layer formerly assigned to Vesuvius 79 CE is now placed at 88 CE, with a smaller volcanic sulfate concentration peak at 80 CE in NGRIP2 and GISP2 (Fig. 1b). We identify in the northern Greenland NEEM-2011-S1 ice core a distinct tephra associated with the 88 CE event. We consider the potential source of the tephra and examine the consequences of misattributing ice-core volcanic signals to specific events in terms of geochronology and paleoclimate reconstructions.

2 Methods

Three annual ice samples bracketing a sulfate peak at 410.56 m in the NEEM-2011-S1 core were cut for tephra analysis between 410.20 and 410.85 m using a bandsaw in the ice-core laboratory at the Desert Research Institute (DRI), Nevada (Fig. 1c). Samples were transferred to precleaned Nalgene bottles and sent to Queen's University Belfast (QUB) for tephra extraction where they were assigned unique identifier codes (QUB-1832: 410.85–410.65 m, 87–88 CE; QUB-1833: 410.65–410.40 m, 88–89 CE; QUB-1834: 410.40–410.20 m, 89–90 CE). The meltwater was centrifuged to concentrate particulates, which were then pipetted onto pre-ground glass slides that were placed on a hotplate in a laminar flow cupboard. The samples were covered in Buehler epoxy resin. Tephra shard content was quantified using ×100–400 magnification on a light microscope with the aid of cross-polarization. Samples containing tephra were ground using 12 µm alumina powder and polished using 6, 3 and 1 µm diamond paste until the surfaces of the shards were exposed for analysis.

To determine the source of the tephra, major element geochemical analysis was performed on a JEOL FEGSEM 6500F electron microscope at QUB, and trace element composition was analyzed by laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) on a Thermo SCIENTIFIC iCAP Q coupled to a Photon Machines Analyte Excite Excimer UV laser at the iCRAG Lab at Trinity College Dublin. Methods and operating conditions are presented in Dataset 1 in the Supplement. Major element compositions have been normalized to 100 % to enable comparison with published material (raw and normalized data are presented in Dataset 1 in the Supplement).

To evaluate possible source eruptions for the recorded tephra, a list of all known eruptions within the period 100 BCE–300 CE with an estimated volcanic explosivity index (VEI) of ≥ 4 was collated from Siebert et al. (2011), supplemented by data from the Smithsonian's Global Volcanism Program (GVP; http://volcano.si.edu/ (last access: 21 February 2018); Dataset 2 in the Supplement). Stevenson et al. (2013) have demonstrated that fine ash (> 32 µm) can be dispersed > 1000 km from low plumes (<5 km). Plunkett and Pilcher (2018) concluded that only larger (≥ VEI 5) eruptions from more distant volcanic regions are presently represented in the ultra-distal cryptotephra records in northwestern Europe. Consequently, while acknowledging that small to moderate eruptions in Icelandic or the Jan Mayen islands are within range of dispersing tephra to Greenland, we consider VEI 4 a sufficient threshold for exploring potential sources at distances of 4000 km or more. The wide time span allows for chronological uncertainty of the events, of which only Vesuvius 79 CE is historically dated. In view of Iceland's proximity to Greenland, we also considered smaller-magnitude eruptions from this region, excluding volcanic systems that typically produce homogenous basaltic glasses, such as Grímsvötn and Bárdarbunga. Where possible, reference glass geochemical data for the events were collated from published records for comparison with the tephra identified in the NEEM-2011-S1 ice core and were added to the QUB tephra database. For possible correlatives, we analyzed the geochemical composition of reference glass samples under the same operating conditions as the NEEM-2011-S1 tephra (see the Supplement for sample details).

Figure 2 Selected geochemical biplots and a spider plot showing compositional similarities of glass shards from samples QUB-1832 and QUB–1833 from NEEM-2011-S1. Plate setting differentiations are based on Tomlinson et al. (2015), and the MORB-normalized spider diagram follows Pearce (1983).

On the basis of three potential source volcanoes suggested by the tephra major element geochemistry (namely Aniakchak, Chikurachki and Popocatépetl; see Fig. 1), we modeled ash distribution using the advection–dispersion–sedimentation software Ash3d (Schwaiger et al., 2012) to evaluate the necessary eruption and meteorological parameters needed to transport ash as far as the NEEM coring site (77^∘27^′ N, 51^∘3.6^′ W; 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). Ash3d simulates the source plume by inserting ash in a column of nodes above the vent. This source ash is then advected as it sediments out of the atmosphere. No accelerated fallout mechanisms were considered in these simulations such as aggregation of ash or removal due to interactions with hydrometeors. Ash3d calculates tephra deposition as mass per unit area at sites, even if the fallout is vanishingly small. Hence, even exceedingly small fallout values contributing very low concentrations of tephra can be identified as nonzero values (Dunbar et al., 2017). Model parameters were selected on the basis of each volcanic source under consideration (see Sect. 3.4). As the meteorological conditions at the time of the eruption are not known, we used data drawn from the 2.5^∘ NCEP-NCAR Reanalysis II dataset (https://climatedataguide.ucar.edu/climate-data/ncep-reanalysis-r2, last access: 1 March 2021) that captures 3D time-varying atmospheric structure (including temperature and horizontal and vertical winds) up to a pressure level of 10 mb. This pressure level corresponds to a geopotential height of approximately 31 km. Start dates and times for the events were selected randomly from the period between 1950 and 2010 but limited to just the winter months (November–February), reflecting the likely season of the eruption based on sulfate and tephra deposition in the ice core (Sigl et al., 2015; McConnell et al., 2020). For each of the three sources, 350 simulations were run as an initial sensitivity test on the basis of approximated eruption characteristics. Additionally, to distinguish meteorological differences from differences in the eruption source parameters and allow a broader consideration of potential source parameters, we randomly selected 1000 start times from 1950 to 2010 and simulated events using identical source parameters at the three locations of interest.

3 Results

3.1 The NEEM-2011-S1 tephra

Cryptotephra particles were identified in samples QUB-1832 (n = 18) and QUB-1833 (n = 5), demonstrating that ash deposition preceded the maximum of sulfuric aerosol depositions at this location in Greenland. Shards were pale brown, ranged between 24 and 85 µm in longest axis, and displayed a varied morphology, including blocky, bubble-wall and platy shards, some with elongated or occasionally spherical vesicles (the Supplement). Several shards contained minerals, and one was microlite-rich, but most shards were aphyric. Major element analyses of the glass were obtained for six and three shards in these samples, respectively, and trace element data for five and three of the shards. The data reveal a homogenous andesitic glass that evidently derives from a single source, and the sample is hereafter named QUB-1832/3 (representing ice spanning 410.85–410.40 m, or ∼ 87–89 CE but very likely occurring in the winter of 87–88 CE; Sigl et al., 2015). Major and trace element compositions suggest a source along an active arc (see Tomlinson et al., 2015; Fig. 2).

The GVP database lists five VEI 6, three VEI 5 and 31 VEI 4 eruptions dating to between 100 BCE and 300 CE. The VEI 6 eruptions include Churchill (Alaska), Okmok (Alaska) and Ksudach (Kamchatka), whose chronologies have all been recently revised and no longer encompass the first century CE (Ponomareva et al., 2017; McConnell et al., 2020; Reuther et al., 2020). The magnitude of Ambrym, Vanuatu – previously thought to have entailed a cataclysmic collapse of the caldera – has also been questioned and was likely a less pronounced, gradual event (Németh et al., 2009). Of the three VEI 5 eruptions, only Masaya, Nicaragua, is located on an active arc. We have been unable to find major element glass geochemistry for many of the VEI 4 eruptions, but we draw upon trace element data (for glass or whole rock), tectonic setting and information from other eruptives to evaluate their correspondence with QUB-1832/3 (Dataset 2 in the Supplement).

Figure 3 Comparison of QUB-1832/3 glass geochemistry with known VEI 4–6 eruptions in the approximate timeframe. (a) Total alkali–silica plot (following Le Bas et al., 1986) and (b) FeO_total CaO biplot, showing major element geochemical differentiation of QUB1832/3 from known VEI 4 to 6 eruptions in the period 100 to 300 BCE and the GRIP 492.2 m tephra recorded by Barbante et al. (2013). 1 – Churchill WRAn (Preece et al., 2014); 2 – Ksudach KS₁ (Plunkett et al., 2015); 3 – Taupo (Newnham et al., 1995); 4 – Okmok (McConnell et al., 2020); 5 – Apoyeque Chiltepe (Kutterolf et al., 2008); 6 – Masaya Tuff (Kutterolf et al., 2008); 7 – Furnas C (Wastegård et al., 2020); 8 – Vesuvius (Santacroce et al., 2008); 9 – Shiveluch units 18 and 21 (Ponomareva et al., 2015); 10 – Glacier Peak A (Foit et al., 2004); 11 – Tacaná Mixcun (ICP analysis; Macías et al., 2000); 12 – Pelée P3 (Martel and Poussineau, 2007; 13 – Cotopaxi KA1 (Garrison et al., 2011); 14 – Merapi (Andreastuti, 1999); 15 – El Misti (Legros 2001) [note: Fe value shown is Fe₂O₃]; 16 – Katla SILK-Yn (Dugmore et al., 2000); 17 – Snaefellsjökull Sn-1 (Larsen et al., 2001); 18 – Torfajökull Grakolla (Oládóttir et al., 2011); 19 – Askja A-2000/Glen Garry (Barber et al., 2008).

3.2 Major element geochemistry

The QUB-1832/3 glass geochemistry is clearly distinct from Vesuvius glass and from the data presented for the GRIP ash particles by Barbante et al. (2013). On the basis of major element compositions, we can also firmly reject correlations with the VEI 6 eruptions from Ksudach, Okmok, Churchill (WRAn tephra), Apoyeque (Chiltepe tephra; Nicaragua) and Taupo (New Zealand), as well as the VEI 5 eruptions from Furnas (Azores) and Masaya (Fig. 3a and b). We lack a full suite of data for Ambrym, but glasses from this eruption do not include an andesitic component (Robin et al., 1993). We find no geochemical matches with any of the VEI 4 eruptions for which we have data. We also exclude Iceland as a source, as the geochemistry matches neither the known eruptions of this time period nor any other Icelandic tephras for which we have data (e.g., Gudmundsdóttir et al., 2012; 2016; Plunkett and Pilcher, 2018). Similarly, we discount Jan Mayen from consideration as only basaltic tephras have been reported from this source (Gjerløw et al., 2016), and no eruption is known in the timeframe under consideration. The geochemistry does not correlate with any known Kamchatkan source (Kyle et al., 2011; Ponomareva et al., 2017; V. Ponomareva, personal communication, 2017).

Figure 4 Individual glass shard major element composition of QUB-19832/3 compared with close correlatives, including new point analyses (symbols) on OSC1–5 and Popocatépetl (including samples from the Pink and Lorenzo Pumices, as well as an unspecified “Popocatépetl” ash sample from the Smithsonian Institute). Compositional fields are based on published data: Aniakchak (Kaufman et al., 2012; Davies et al., 2016); OSC1–5 (Hasegawa et al., 2011); Popocatépetl matrix glass (Siebe et al., 1996). Mount Spurr points represent published mean values (Riehle et al., 1985; Béget et al., 1994; Child et al., 1998). QUB-1859 data are from Sigl et al. (2015).

Extending our search beyond the age range of the QUB-1832/3 tephra, we find geochemical correlations with two shards (shards 5 and 7, for which there are two and four point analyses, respectively) within the mixed population of the previously published QUB-1859 from NEEM-2011-S1, dated to 536 CE (NS1-2011 chronology; Sigl et al., 2015; Fig. 4). A comparison with published data reveals a close similarity to the andesitic component of Aniakchak (Kaufman et al., 2012; Davies et al., 2016) and to material from Mount Spurr (Riehle, 1985; Béget et al., 1994; Child et al., 1998). QUB-1832/3 has a slightly lower TiO₂ content than the Aniakchak tephra, however. Although no large Aniakchak eruptions are known from the timeframe of interest, the GVP reports an eruption with no assigned VEI at 255 ± 200 CE. Microbeam data are not available for the Mount Spurr tephras, but the average FeO_total content also tends to be lower than the NEEM-2011-S1 tephra (Fig. 4). No activity at this volcano has been reported in the timeframe of the 88 CE event.

Other close major element correlations include OSC1–5 (a series of five eruptions dating to between 1650 BCE and 400 CE) from Chikurachki, Kuril Islands (Hasegawa et al., 2011), and matrix glasses from Popocatépetl, Mexico (Siebe et al., 1996). The Chikurachki eruptions are not listed in the GVP database, but large Popocatépetl eruptions are dated to 200 ± 300 BCE (VEI 5: Lorenzo Pumice) and 822–823 CE (VEI 4: Pink Pumice). The latter date is based on an assumed linkage with a sulfate peak in the Greenland ice core GISP2 (Siebe et al., 1996, 2017) at this time (817 CE on the NS-2011 timescale), but the event has also been reported as dating to ∼ 750 CE (Siebe and Macías, 2004). We reanalyzed OSC1–5 glass and reference samples from Popocatépetl, including material from the Lorenzo and Pink Pumices, and from reference ash samples provided by the Smithsonian Institution on the JEOL FEGSEM 6500F (Fig. 4). Our results show that QUB-1832/3 differentiates from OSC1–5 in FeO_total and CaO content, rendering Chikurachki an unlikely source of the NEEM-2011-S1 tephra. Small but critical divergences between our data and published glass values for Chikurachki (Hasegawa et al., 2011) likely reflect differences in data quality. Analyses from the Pink and Lorenzo Pumices from Popocatépetl demonstrate homogenous dacitic to rhyodacitic glass compositions that do not overlap with the andesitic data points of QUB-1832/3. The Smithsonian samples – from an unspecified event or events – show greater compositional variability along a trend from rhyolitic to basaltic trachyandesitic. The majority of analyses obtained from these samples have a rhyolitic to dacitic composition, seemingly on a compositional trend towards the QUB-1832/3 population except in FeO_total concentrations. Some rhyolitic shards of QUB-1859 (shards 1, 10, 11 and 2: Sigl et al., 2015) attributed to the 536 CE fallout closely match the high-silica Popocatépetl analyses, but divergence in TiO₂ (shards 1, 10, 11) or K₂O (shard 2) is evident.

Figure 5 Trace element composition of QUB-1832/3 compared with tephras from known large eruptions between 100 BCE and 300 CE. WR indicates data from whole rock analysis; for averaged glass analyses and compositional ranges, the number of data points is given. (a–d) Normalized primitive mantle composition (following Sun and McDonough, 1989) of QUB-1832/3 compared with potential source regions and/or volcanoes: (a) Lesser Antilles field (Toothill et al., 2007; Dufrane et al., 2009; Cassidy et al., 2012; Labanieh et al., 2012; Howe et al., 2014, 2015), with average values for Liamuiga (Toothill et al., 2007) and Pelée (Labanieh et al., 2012); (b) whole rock composition for the first century CE eruptions of Tacaná (Mixcun Flow Deposit; Mora et al., 2013) and Turrialba (TU37/1; DiPiazza et al., 2015) with the glass trace compositional field of more recent Turrialba eruptives (DeVitre et al., 2019); (c) averaged whole rock values for the first to second century CE eruptions of Guagua Pichincha (Samiengo et al., 2010) and Cotopaxi (Garrison et al., 2011), as well as the compositional field of Calbuco, including the eighth century CE Ca13 eruption (Watt et al., 2011); (d) Japanese volcanic zones after Albert et al. (2019); (e, f) selected biplots showing variance in selected high-field-strength elements.

3.3 Trace element geochemistry

The trace element signature of QUB-1832/3 is clearly distinct from candidate eruptions in Latin America (Fig. 5a–c). Lesser Antilles (Caribbean) sources have a Rb–U and Ta > Na signature that distinguishes them from the cryptotephra. Although the QUB-1832/3 glass is less evolved than the Japanese tephras reported by Albert et al. (2019), its trace element composition is inconsistent with Japanese sources (Fig. 5d).

Figure 6 Trace element composition of QUB-1832/3 against close correlatives in major element glass composition. WR indicates data from whole rock analysis. Aniakchak field based on individual glass shard analyses in Kaufman et al. (2012). Mount Spurr whole rock (WR) data based on Nye et al. (1995) and George et al. (2003). (a) Normalized primitive mantle composition (following Sun and McDonough, 1989); (b–d) selected biplots of high-field-strength elements (HFSEs).

Our best matches using major element geochemistry are with two andesitic shards in QUB-1859, various Holocene andesitic eruptives of Aniakchak and two recent eruptives from Mount Spurr; our analysis of Popocatépetl is inconclusive. Trace element data were collected from the pertinent shards in QUB-1859 and confirm that these shards derive from the same source as QUB-1832/3 (Fig. 6). Attempts to analyze other shards within QUB-1859, specifically those lying close to the rhyolitic end of the Popocatépetl compositional field, were unsuccessful due to the entrapment of minerals.

We compared our results with published trace element data for Aniakchak glass (Kaufman et al., 2012) and with newly obtained data from the four Smithsonian reference samples from Popocatépetl. Our analyses show that the ice-core tephras are lower in trace element concentrations than Aniakchak but share consistent high-field-strength element (HFSE) ratios, and a link with this source or source region cannot be excluded (Fig. 6b–d). In contrast, clearly distinctive HFSE ratios (e.g., Th: Nb; Fig. 6d) differentiate the tephra from Popocatépetl and allow this source to be eliminated. We were unable to locate trace element data for Mount Spurr glass but have compared our results with whole rock data from the 1953 and 1992 eruptions.

Figure 7 Ash3d simulations of ash fallout from Aniakchak (left), Chikurachki (center) and Popocatépetl (right) following randomized start dates between 1950 and 2010: (a) 16 km high column, eruption duration 9 h, erupted ash volume 0.3 km³; (b) 20 km high column, eruption duration 12 h, erupted ash volume 0.41 km³; (c) 30 km high column, eruption duration 8 h, erupted ash volume 1.75 km³. The spectrum denotes deposit thickness (red being thickest).

3.4 Modeling ash dispersal

Our geochemical analysis of the QUB-1832/3 tephra demonstrates that it derives from an eruption other than a known event of VEI 4 or higher listed in the GVP database. The geochemical composition also enables us to eliminate from consideration well-characterized volcanoes in Iceland, Kamchatka, and, on the basis of trace element composition, the Lesser Antilles and Japan. Other likely source regions include Alaska, the Kuril arc and Mexico. We therefore explore the eruption parameters needed for tephra to be transported from these sources to Greenland by modeling ash dispersal from the three candidates (Aniakchak, Chikurachki and Popocatépetl) highlighted by major element geochemical similarities.

Despite a lack of clear correspondence in the major and trace element composition, we cannot firmly reject Aniakchak or a neighboring system as a possible source of the QUB-1832/3 tephra. Aniakchak is not known to have erupted in the timeframe of interest, implying that any such eruption is unlikely to have been very large or deposits would have been preserved locally. Since the caldera-forming eruption (CFE) at ∼ 1640 BCE (GICC05 chronology) whose tephra has been securely identified in Greenland ice (Coulter et al., 2012; Plunkett et al., 2017), at least two VEI ≥ 4 eruptions of Aniakchak have occurred. To simulate ash dispersal, we therefore used parameters for a high VEI 4 event erupting a volume of 0.5 km³ dense rock equivalent (DRE) with a 12 km plume height over a duration of 15 h with a grain size distribution of 500 µm (70 %), 80 µm (20 %) and 30 µm (10 %). This distribution approximates the fine component observed in other andesitic Aleutian volcanos, while the coarse component is consolidated into one grain size bin that is rapidly deposited. The particle fall model used in the Ash3d simulations is that of Wilson and Huang (1979) whereby a shape factor, F (the ratio of the average of the minor axes to the major axis), determines the drag coefficient. We used F = 0.44 for all simulations, which is the average value of volcanic ash from the Wilson and Huang study. Particle size for this fall model corresponds to the arithmetic average of the ellipsoid axes. In ∼ 10 % of the simulations, fine tephra (∼ 30 µm) reached the NEEM ice-core drilling site, but not larger particles. The model suggests that a VEI 4 eruption could, under suitable meteorological conditions, disperse ash to Greenland, but probably not a smaller event (Fig. S1 in the Supplement).

Our major element data suggest that Chikurachki is an unlikely source of QUB-1832/3, but large (VEI 4) eruptions from this and other Kuril eruptions are capable of injecting sulfur into the stratosphere, as demonstrated by the 1986 Chikurachki eruption (0.75 Tg of SO₂ emitted to an elevation of 11 km; Global Volcanism Program, 2013) and the 2019 Raikoke (13–17 km plume height erupting 1.5 ± 0.2 Tg SO₂; de Leeuw et al. 2021). The volcanic history of the Kuril Islands is poorly constrained prior to the settlement of the region in the 18th century, but at least four other VEI 4 eruptions from Chikurachki are known to have occurred during the middle to late Holocene (Siebert et al., 2011). We therefore examined the potential ash dispersal of a comparable event – the 1986 eruption (0.25 km² DRE, 11 km high plume, eruption duration 7 h, with grain size distribution of 500 µm (70 %), 80 µm (20 %) and 30 µm (10 %) as described by Belousova and Belousov, 2002) – to transport tephra as far as northern Greenland. In none of the 350 Ash3d simulations did tephra reach Greenland (Fig. S2 in the Supplement).

Thirdly, we examine the possible dispersal of Popocatépetl ash. Recent eruptions have been of a modest nature here (VEI 3 or less), injecting small volumes of ash and sulfur into the troposphere. Larger eruptions dating to 200 ± 300 BCE (VEI 5) and ∼ 750 CE (VEI 4) are inferred from tephra fall deposits (known respectively as the Lorenzo and Pink Pumices; Siebe and Macías, 2004). For the Ash3d simulations, we used the grain size parameters outlined for an older VEI 6 event, the ∼ 3700 BCE Ochre Pumice (Arana-Salinas et al., 2010), in combination with estimates for the ∼ 200 BCE VEI 5 Lorenzo Pumice (Panfil et al., 1999): a plume height of 30 km, eruption duration of 20 h, an eruptive volume of 1 km³ DRE, with a grain size distribution of 250 µm (15 %), 125 µm (20 %), 63 µm (40 %) and 31 µm (25 %). For the occasional meteorological conditions in which the maximal height of the airborne ash was higher than data provided by the 10 mb pressure level, wind vectors were held constant at the values of the 10 mb level and temperatures were calculated from the US Standard Atmosphere. Such instances seldom occurred and affected the top several hundred meters of the simulation domain. In 32 % of the Ash3d simulations, ash reached Greenland, including larger grains (63 µm) in eight (2 %) of the simulations. Ash also dispersed over the Africa, Europe, central Asia and the North American circumpolar region, demonstrating the potential reach of ejecta from large Mexican eruptions (Fig. S3 in the Supplement).

These simulations used eruption source parameters (plume height, eruptive volume, duration, grain size distribution) that were drawn from past events for Chikurachki and Popocatépetl (see Sect. 3.4). To address the question of how big or how small an eruption would need to be or could be for ash from each volcano to fall out in Greenland, we randomly selected 1000 start times and applied the same eruption source parameters to each of the three volcanoes, with eruptions commencing at randomly selected times between 1950 and 2010. Eruptive volume was selected from a lognormal distribution between 0.01 and 2.0 km³ (dense rock equivalent). Plume heights were calculated using an empirical best-fit relationship between eruptive volume and plume height (H = 25.9 + 6.64 ⋅ log 10(V), where V is erupted volume in DRE, and H is height above the vent in km) with a random adjustment (Gaussian with μ = 0 km and σ = 2.9 km; Mastin et al., 2020, Eq. 3). Mass eruption rate was calculated from a best-fit relationship with plume height (Mastin et al., 2020, Eq. 2) and duration using the mass eruption rate and eruptive volume (see Mastin et al., 2020). For all of these 1000 cases, the same grain size distribution (GSD) was used (15 % 0.250 mm, 20 % 0.125 mm, 40 % 0.063 mm and 25 % 0.031 mm). Using a GSD with this greater fine component and considering plume heights in a broader range (7–30.5 km), hypothetical eruptions at each of the three volcanoes resulted in ash deposited in Greenland. The outputs of all simulations are archived at Plunkett et al. (2021).

Figure 8 Summary of Ash3d modeling of tephra fallout in Greenland: (a) distribution of total runs (black line) for 1000 randomized eruption scenarios along with the number of simulations resulting in 31 µm ash (colored lines) and 63 µm ash (dots) fallout at NEEM for each hypothetical eruption; (b) the percentage of cases with fallout at NEEM for all months; (c) distribution of plume heights for randomized eruption scenarios during November to February, along with the number of corresponding cases with 63 µm fallout at NEEM for each hypothetical eruption (dots); (d) the percentage of simulated cases with fallout at NEEM for the winter months (November to February).

Generally, for both Aniakchak and Chikurachki, the greater the plume height of the hypothetical eruption, the more likely ash will fall out in Greenland (Fig. 8). With Popocatépetl, the trend is not as pronounced. In all cases, no 31 µm ash was deposited in Greenland for plume heights less than 10 km for Aniakchak or 12 km for Chikurachki and Popocatépetl. In less than 1 % of cases for both Aniakchak and Popocatépetl were 63 µm particles deposited, and exceptionally large plume heights were required to achieve this (27 and 19 km, respectively). In none of the 1000 simulations of Chikurachki did 63 µm ash reach Greenland.

For these hypothetical cases, 43 % of the simulations of Aniakchak resulted in fine ash deposited in Greenland as opposed to 24 % and 13 % from Chikurachki and Popocatépetl, respectively. There is a weak sensitivity to season, with a slight increase in the odds of ash reaching Greenland for each of these three volcanoes during wintertime (50 %, 26 % and 20 % for Aniakchak, Chikurachki and Popocatépetl; Fig. 8d). Again, these likelihoods are only with consideration of the source location with respect to the meteorology. Additionally, there are other mechanisms that could play a role that were not considered in these simulations for both accelerated removal mechanisms, such as ash aggregation, and mechanisms for prolonging the transport, such as vertical mixing. Nevertheless, when considering the geologic history, Chikurachki becomes much less likely due to its historically lower plume heights and coarser grain sizes.

4 Discussion

4.1 Deconstructing the ice-core “Vesuvius” marker horizon

We identify in the NEEM-2011-S1 ice core a significant concentration of large cryptotephra shards in association with the sulfate layer previously attributed to the Vesuvius 79 CE eruption. Our geochemical characterization of the tephra conclusively demonstrates that it is not derived from Vesuvius. Coupled with the revised NS1-2011 ice-core chronology that now places this event in winter 87–88 CE (Sigl et al., 2015), our findings call into question the robustness of the geochemical data from glass shards in the GRIP ice core that were reported by Barbante et al. (2013). Instrumental inaccuracy arising from the use of energy-dispersive spectrometry (EDS) on unpolished shards in the GRIP sample cannot explain the significant geochemical differences between the GRIP and NEEM tephras, which are clearly not from the same source. The Barbante et al. (2013) data lie close to the Vesuvius compositional field but differ in Mg and Ti content, elements that should not be greatly altered by EDS analysis. In contrast, alkalis, which are mobile and prone to loss under a focused electron beam, map perfectly onto the Vesuvius range after correction for the expected loss of Na in a phonolitic tephra such as Vesuvius (Barbante et al., 2013). It is therefore possible that the geochemical attribution of the GRIP tephra to Vesuvius is simply erroneous, arising from the analysis of unpolished shards and correction of Na in line with expected values. Alternatively, if the GRIP tephra is indeed from Vesuvius, the validity of very small ash particles (< 5 µm) as time-synchronous markers in polar ice cores is questionable. In our simulations, ash was deposited almost entirely within days to weeks following the event, although the model is not sufficiently tuned to consider extremely fine particles approaching dust size. We consider it nevertheless unlikely that Vesuvius ash was being deposited in Greenland up to 9 years after the eruption and propose that the GRIP misattribution is better explained by issues stemming from the methodological approach used to analyze the GRIP tephra. Under the revised NS1-2011 and DRI_NGRIP2 chronologies (Sigl et al., 2015; McConnell et al., 2018), a small volcanic sulfate peak in NGRIP2 and GISP2 at 80 CE is observed, which may derive from the VEI 5 eruption of Vesuvius in 79 CE or equally an unknown eruption of similar age (Fig. 1b). Drilling of the NEEM-2011-S1 core was stopped at an age of 86 CE, so we could not sample the corresponding depth range for tephra. We suggest that additional screening for potential tephra should be extended by the ice-core community to the corresponding core sections of archived ice cores (e.g., NGRIP, NEEM, Dye-3 and others).

Although we have been unable to determine the source of our 88 CE tephra, we rule out on the basis of major and trace element glass geochemistry the likelihood that it was from Iceland, Greenland's closest volcanic source, or from any known event of VEI ≥ 4 presently dated to within ∼ 200 years of the event, emphasizing that eruption inventories are incomplete and often imprecisely dated (Brown et al., 2014, Rougier et al., 2016). Our ash models suggest that tephra from smaller-magnitude eruptions in Mexico or the North Pacific region is unlikely to reach Greenland and that plumes from the northwestern Pacific must attain elevations of 10 km or more to have the potential to disperse ash as far as the North Atlantic. Stevenson et al. (2015) similarly found that ash from 10 km high plumes in Iceland could rapidly travel 850 km, but our findings extend the range of moderately high ash clouds to considerably greater distances. Presently, ultra-distal records of cryptotephra in northwestern Europe appear limited to VEI ≥ 5 eruptions (Plunkett and Pilcher, 2018), but there is clearly potential for tephras from moderate eruptions to disperse to Greenland. For ash particles as large as those in QUB-1832/3 and QUB-1859 to disperse such distances, however, much higher (> 19 km) plumes are likely needed. The glass geochemistry of the 88 CE tephra most closely resembles products from the eastern Aleutian range, although no candidate eruption can be identified despite the likelihood that it was a large-magnitude event. This serves to emphasize the incompleteness of the volcanic record, a salutary reminder that large, unknown eruptions may be candidates for prominent volcanic events in Greenland ice cores even when there are close temporal associations with documented eruptions.

4.2 Geochronological implications

The misattribution of the 88 CE signal can be traced back to early ice-core research in which Vesuvius was posited as a possible source because the ice layer counts were within 10 years of the historical eruption (Zielinski, 1995). Subsequent adherence to this tie point within the Greenland ice-core chronology, seemingly supported by the presence of tephra attributed to Vesuvius (Barbante et al., 2013), impeded the recognition of coeval historical and proxy signals of past environmental change despite calls for a reconsideration of the ice-core dating on the basis of similar intervals between events in the ice-core and tree-ring records (Baillie, 1996, 2008; Baillie and McAneney, 2015; Adolphi and Muscheler, 2016). Nonetheless, the marker event played a pivotal role in constraining the Common Era section of the GICC05 chronology (Vinther et al., 2006), which is widely considered to be the reference chronology for Northern Hemisphere paleoclimate (Rasmussen et al., 2014). All major deep Greenland ice cores are aligned on GICC05 (Rasmussen et al., 2013; Seierstad et al., 2014; Vinther et al., 2006, 2008), including the most recently drilled EGRIP ice core (Mojtabavi et al., 2020), similar to the way in which Antarctic ice cores are increasingly aligned to the WD2014 chronology (Baggenstos et al., 2018; Buizert et al., 2018; Cole-Dai et al., 2021; Sigl et al., 2016; Winski et al., 2019).

Figure 9 Volcanic forcing and temperature response in the mid-sixth century CE: annual mean stratospheric aerosol optical depth between 30 and 90^∘ N (grey; Toohey and Sigl, 2017), JJA temperature anomalies from a network of extratropical tree-ring sites in Eurasia (red; Büntgen et al., 2020), and mean annual temperature anomalies from three stacked annual-layer-dated Greenland ice cores (blue; Dye-3, NGRIP1 and GRIP; Vinther et al., 2006; PAGES 2k Consortium, 2013) using (a) the revised NS1-2011 chronology (Sigl et al., 2015) and (b) the original GICC05 chronology (Vinther et al., 2006).

The development of the independent NS1-2011 chronology, unconstrained by this tie point (and another previously assumed tie point at 1105 CE ascribed to the Hekla 1104 CE eruption), has since enabled a more precise dating, confirmed by tephra linkages, of past volcanic eruptions (e.g., Ilopango 431 ± 2 CE – Smith et al., 2020; Katla 822 CE – Plunkett et al., 2020; Churchill 853 CE – Jensen et al., 2014; Eldgjá 939–940 CE – Oppenheimer et al., 2018; Changbaishan 946 CE – Sun et al., 2014; Veiðivötn 1477 CE – Abbott et al., 2021). The chronology has also enabled a more accurate appraisal of volcanic effects on climate and society, including on Northern Hemisphere summer temperatures (Büntgen et al., 2020; Helama et al., 2021; Zhu et al., 2020), Nile summer flood variability and socioeconomic effects in ancient Egypt (Manning et al., 2017), and Chinese dynastic collapses (Gao et al., 2021). It is largely consistent with the DRI_NGRI2 chronology, an independent annual-layer-counted chronology from the NGRIP2 ice core (McConnell et al., 2018, 2019) used to constrain the date and climate effects of the caldera-forming Okmok II eruption in winter 44–43 BCE (McConnell et al., 2020).

Over the last decade, numerous continental- to global-scale temperature reconstructions have been made based on multi-proxy records that include ice cores (McKay and Kaufman, 2014; PAGES 2k Consortium, 2013, 2017, 2019; Neukom et al., 2019; Kaufman et al., 2020; Konecky et al., 2020) and data assimilation employing these databases (Steiger et al., 2018; Tardif et al., 2019). These reconstructions are increasingly used as a benchmark in assessments of anthropogenic global warming (Steig et al., 2013; Abram et al., 2016; Neukom et al., 2019b) or to explore the role of natural forcing on past climate (Tejedor et al., 2021). In these reconstructions, the relative contribution of ice-core records to the overall proxies is generally increasing with increasing age. In the most recent PAGES 2k reconstruction (PAGES 2k Consortium, 2017), 48 ice cores (8 %) are included for 1800–1850 CE and 22 (24 %) for 500–550 CE. Of the latter 22 records, 10 have been aligned to the GICC05 chronology using volcanic synchronization and only two have been aligned to the new WD2014 chronology (Sigl et al., 2016). Because of the reduced fraction of proxies with annual dating precision (i.e., tree rings) in the first millennium CE and in combination with coarse resolution, dating uncertainties and biases in other proxies, multi-proxy reconstructions consequently depict smaller amplitudes of natural climate variability in the first millennium CE than in the second millennium CE under comparable external climate forcings (Büntgen et al., 2020). Accurate alignment of the available annual resolution data – tree rings and polar ice cores – is critical to redress this issue and hinges on the integration of data based on a robust ice-core chronology (see, for example, Fig. 9). Given our firm refutation of the Vesuvius 79 CE tie point in the GICC05 chronology, we urge that corresponding databases (e.g., PAGES 2k, PAGES Iso2k, Last Millennium Reanalysis, Paleo Hydrodynamics Data Assimilation) be updated to accurately capture the timing and full amplitude of paleoclimatic variability, pending the forthcoming revision of GICC05 (Sune Rasmussen, personal communication, 2021).

4.3 Sulfate deposition for 88 CE and implications for climate forcing

Unsupported linkages between volcanic signals and their source have consequences for the understanding of volcanic impacts. For pre-satellite-era eruptions, determination of atmospheric aerosol loading – the leading mechanism by which volcanoes impact climate – hinges upon models of aerosol dispersal and deposition based on suspected source and ice-core sulfur concentrations (Gao et al., 2007). Injection height may also determine the duration of aerosol persistence in the stratosphere and, consequently, its climate effectiveness (Toohey et al., 2019; Marshall et al., 2019; Aubry et al., 2020). Some prior knowledge of the source and eruption parameters is therefore needed to refine the modeling and interpretation of climate response.

Here, our ash dispersal models enable us to constrain the minimum injection height needed for eruptions in three volcanic regions to extend tephra to Greenland. Our results demonstrate that stratospheric injection was very likely required to disperse large ash particles to northern Greenland from either Aniakchak in the high latitudes or Popocatépetl in the tropics. These two scenarios have implications for the determination of stratospheric sulfate loading, with a considerably larger sulfur emission expected from the tropical source to account for the aerosol deposition observed in Greenland. Consequently, were this event to be attributed to Popocatépetl rather than high-latitude (e.g., Aniakchak) or midlatitude (e.g., Vesuvius) sources, we would anticipate a larger global mean forcing similar to the Huaynaputina eruption of 1600 CE. Instead, an extratropical location implies an event of smaller stratospheric sulfur injection and therefore lesser global climate significance, roughly equivalent to some estimates of the sulfur injection from the 1912 Katmai eruption (Stothers, 1996). Annual-resolution climate records for this period are rare, but an analysis of limited tree-ring records from northern Eurasia (Büntgen et al., 2020) and ice-core records from Greenland (PAGES 2k Consortium, 2013) suggests limited and spatially variable cooling signals associated with the 88 CE eruption that are indistinguishable from the usual variability.

Although inconclusive, the tephra evidence points to the eastern Aleutians as the most likely origin of the 88 CE eruption. Whatever the source, the same volcano contributed to the volcanic signal in Greenland ice cores in 536 CE, following which there was a major climate perturbation with societal consequences (Sigl et al., 2015; Büntgen et al., 2016). Other tephras in the NEEM-2011-S1 ice core at this time illustrate that the 536 CE sulfur signal is potentially a composite of three simultaneous eruptions of uncertain origin (Sigl et al., 2015), confounding the estimation of atmospheric aerosol loading. Like QUB-1832/3, the large (> 50 µm) shard sizes in the 536 CE ice imply stratospheric injection. The Okmok II eruption of 43 BCE has also been tied to climate impacts (McConnell et al., 2020), clearly demonstrating the significant forcing potential of high-latitude eruptions.

Table 1 Volcanic sulfate deposition and duration of above-background sulfate levels for three volcanic events in Greenland ice cores (Sigl et al., 2013, 2015; McConnell et al., 2018, 2020; Clausen et al., 1997; Zielinski, 1995; Toohey and Sigl, 2017).

The 43 BCE, 88 CE and 536 CE events seem to share a similar volcanic source region and magnitude (VEI 4+ or greater), but their aerosol representation in Greenland ice cores is spatially divergent (Table 1). Ice-core records of cumulative volcanic sulfate deposition in Greenland for the 88 CE event are available from NEEM-2011-S1 (Sigl et al., 2015), NGRIP2 (McConnell et al., 2018), GISP2 (Zielinski, 1995), GRIP and Dye-3 (Clausen et al., 1997). These records demonstrate a strong spatial gradient in volcanic sulfate deposition following this eruption: high sulfate mass deposition is observed in northern Greenland (NEEM: 83 $\mathrm{kg}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), but only moderate deposition is recorded in central and southern Greenland (NGRIP2: 32 $\mathrm{kg}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$; GRIP: 26 $\mathrm{kg}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$; GISP2 29 $\mathrm{kg}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$; Dye 3: 31 $\mathrm{kg}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). The mean deposition (40 $\mathrm{kg}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) over the five ice-core sites is much lower than the deposition following the 536 CE (75 $\mathrm{kg}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) or Okmok II (117 $\mathrm{kg}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) events. Moreover, the duration of volcanic sulfate deposition is markedly shorter (0.9 years) for the 88 CE eruption compared to that of the 536 CE (1.7 years) and Okmok II (2.3 years) events, suggesting more rapid atmospheric aerosol removal. The relative importance of other variables – duration of eruption, contributing effects of other eruptions, prior state of the climate system – cannot presently be determined.

5 Conclusions

The NS1-2011 revised ice core and closely related DRI_NGRIP2 ice-core chronologies that differentiate the conspicuous first century CE volcanic signal in Greenland ice from the historical age of Vesuvius are now supported by tephra evidence that demonstrates the sulfate was from an eruption other than Vesuvius. Although we have been unable to pinpoint the precise source eruption, the geochemical signature of its glass shards most strongly points to the Aleutian arc as the source region. Shard characteristics, combined with ash dispersal models, implicate a large-magnitude eruption with stratospheric injection. Constraining the location and minimum eruption height of the event promises to improve the estimation of stratospheric aerosol loading and radiative forcing when combined with understanding of the dependence of aerosol lifetime on these eruption parameters (Marshall et al., 2019; Toohey et al., 2019). We note, however, that unlike notable events such as the 43 BCE Okmok II and 536 CE unidentified eruptions, the 88 CE eruption is not presently associated with any known major climatic perturbations.

Increasingly, the recognition of solar proton events in tree-ring and ice-core archives is demonstrating offsets in the ice-core chronology, and it is imperative that ice-core chronologies are now updated in recognition of the known dating bias in GICC05. Multi-proxy compilations (PAGES 2k, LMR, PHYDA and others) in which ice-core-derived δ¹⁸O and δD records from Greenland (and also some from Antarctica) are still constrained by the erroneous 79 CE tie point need to be updated to accurately capture natural climate variability. Previous work has demonstrated that only under the revised NS1-2011 chronology is the magnitude of post-volcanic cooling inferred from tree-ring reconstructions over the course of the last 2 millennia consistent (Sigl et al., 2015; Büntgen et al., 2020). Unfounded attributions of volcanic events in ice cores, irrespective of their close temporal coincidence with known eruptions, should be avoided to eliminate unconscious bias in the construction of ice-core chronologies and the interpretation of volcanic impacts.