Hekla Revisited: Fractionation of a Magma Body at Historical Timescales

Hekla is an elongate volcano that lies at the intersection of the South Iceland Seismic Zone and the Eastern Volcanic Zone. We report major and trace element, oxygen isotopic, and H 2 O analyses on rocks, glass, melt inclusions, and minerals from almost all of the historical lavas and tephra deposits. This new dataset conﬁrms the remarkable observation that not only are many eruptions compositionally zoned from felsic to maﬁc, but the extent of zoning relates directly to the length of repose since the previous eruption. Compositional data are consistent with the origin of the basaltic andesites and andesites by fractional crystallization, with no measurable crustal interaction once basaltic andesite has been produced. Although the 1104 CE Plinian rhyolite and 1158 CE dacite are also created by fractional crystallization, uranium–thorium isotopic disequilibria measured by others require that they evolved in a separate body, where magma is stored in a molten state for > 10 4 years. Consistent trace element trends and ratios, as well as oxygenisotopic data, preclude signiﬁcant crustal input into the evolving magma. The phenocryst assemblages are dominated by crystals that formed from their host melt; an exception is the 1158 CE dacite, which contains abundant crystals that formed from the 1104 CE rhyolite melt. A suite of thermobarometers indicates that most crystals formed in the lower crust at temperatures ranging from (cid:2) 1010 to 850 (cid:3) C. Hekla’s unique and systematic petrological time series and geophysical activity are attributed to the unusual geometry of the magma body, which we propose to be a tabular, vertically elongate macro-dike, extending from the lower to the upper crust. The vertical body is recharged with basaltic andesite magma at the end of each eruption, which then undergoes cooling and crystallization until the subsequent eruption. The entire system is supplied by a lower-crustal body of basaltic andesite, which is produced by fractional crystallization of basaltic magma in a reservoir that is thermo-chemically buffered to (cid:2) 1010 (cid:3) C. Cooling and crystallization of recharged basaltic andesite magma in a background geothermal gradient from the lower to the shallow crust accounts for the systematic relationship between repose and composition.


INTRODUCTION
Hekla volcano in Iceland is one of Earth's best examples of the link between petrological and volcanic processes (Fig. 1). Not only do Hekla eruptions initiate with an explosive phase and culminate in effusion of lava, but most of the 18 historically documented eruptions are normally zoned, with more evolved magmas erupting first, followed by increasingly mafic products (Thorarinsson, 1950(Thorarinsson, , 1967. Moreover, and most importantly, the extent of magmatic differentiation is a function of the repose time between eruptions ( Fig. 2; Thorarinsson, 1967), providing a unique opportunity to explore rates of magmatic differentiation over human timescales.
Owing to its classic volcanological and petrological characteristics, Hekla volcano has been subjected to nearly every modern geophysical and analytical method, including seismic and geodetic studies, analyses of actinide-series isotopes (Sigmarsson et al., 1992a;Chekol et al., 2011), melt inclusion measurements (Moune et al., 2007;Portnyagin et al., 2012;Lucic et al., 2016), stable isotopes (Schuessler et al., 2009;Savage et al., 2011;Prytulak et al., 2017aPrytulak et al., , 2017bDeng et al., 2019;Inglis et al., 2019), compositional analyses of individual zircon grains (Carley et al., 2011;Bindeman et al., 2012;Carley et al., 2020), and radiogenic isotopes (Sigmarsson et al., 1991(Sigmarsson et al., , 1992aChekol et al., 2011;Carley et al., 2020). What have been missing are modern, state-of-the-art major and trace element measurements and crystal-scale studies of all the historical eruptive products, which are the focus of this study. Although much of the previous work on Hekla cited above includes elemental analyses, most studies include only a few analyses per eruption or only a few eruptive units, which misses two of the most important features of this volcano and its historical eruptive record: the zoned eruptions and the correlation between differentiation and repose (Fig. 2). Moreover, the datasets are from different laboratories, using a variety of methods and standards. In addition to analyzing several samples from most of the historical eruptions in the same laboratories using consistent methods and standards, we also have performed detailed petrography on the rocks and measured the compositions of many of the phenocryst phases. We also include data on melt inclusions, including volatile contents, and oxygen isotopic ratios of olivine crystals.

Geological setting
Hekla volcano is located near the intersection of the South Iceland Seismic Zone and the western margin of the Eastern Volcanic Zone (Fig. 1). In southern Iceland, spreading partitions into the Eastern and Western Volcanic Zones, which are connected by the South Iceland Seismic Zone, a series of north-south-trending right-lateral 'bookshelf faults' that accommodate leftlateral transfer of plate spreading motion between the Reykjanes Ridge and the Eastern Volcanic Zone (Einarsson, 1991). Hekla diagonally cuts the inside corner of the ridge-transform intersection (Fig. 1). Northeast of Hekla, the Eastern Volcanic Zone is an active rift zone characterized by NE-trending normal faults and Holocene eruptive fissures that have been the source of large volume (2-15 km 3 ) eruptions of tholeiitic basalt (Jó nsson et al., 1997;Einarsson, 2008;Pedersen et al., 2019), as well as numerous NE-trending   Thorarinsson (1967;dashed line), showing that the extent of differentiation and compositional zoning within individual eruptions is dictated by the repose between eruptions. hyaloclastite ridges that formed as subglacial fissure eruptions during the Pleistocene. Torfajö kull volcano, which has a large caldera and is the largest silicic volcano in Iceland, lies $25 km east of Hekla, on the western flank of the Eastern Volcanic Zone. Pyroclastic rocks erupted from Torfajö kull crop out around the base of Hekla.
The Hekla edifice is an $10 km long ENE-trending volcanic ridge bisected by a 5-6 km long fissure extending the length of the ridge (Fig. 1). The ridge formed within the last 3000 years through the eruption and accumulation of lavas erupted from the summit fissure (Thorarinsson, 1950;Sverrisdottir, 2007). The modern edifice is constructed on top of a basaltic hyaloclastite ridge, which is exposed along the north, west, and east sides of the volcano. Basalt has erupted from fissure systems a few kilometers east and NE of Hekla as recently as 1878 and 1913, and to the SW of Hekla in 1554 (Thorarinsson, 1967). The basalts erupted on the peripheral fissures are distinct from the intermediate and silicic rocks produced from the Hekla edifice itself, but because of their transitional alkalic character and geochemical similarities to Hekla products, the compositions of the peripheral deposits are considered by most investigators as parental to the Hekla magmatic system (Jakobsson, 1979;Sigmarsson et al., 1992a). Rhyolite lavas with glassy, brecciated margins related to Torfajö kull volcanic center are exposed at the surface 10 km east of Hekla. Xenoliths of rhyolite ignimbrite ejected during Hekla's 1970 eruption (Thorarinsson & Sigvaldason, 1972;Sigmarsson et al., 1992a) and clasts of rhyolitic hyaloclastite ejected from fissures on the NE shoulder of Hekla (our observations) indicate that rocks related to the Torfajö kull system, or an older equivalent, are also present at depth beneath Hekla.
In addition to its common intermediate effusive eruptions, Hekla has experienced five silicic Plinian eruptions, the oldest of which occurred $7000 years ago. The tephra from these eruptions have volumes of 0Á5-2Á2 km 3 dense-rock equivalent (DRE) and have been deposited at intervals of 1000-2500 years (Larsen & Thorarinsson, 1977). Several of the prehistoric Plinian tephra deposits exhibit compositional zoning from rhyolite at the base and middle of the deposits to dacite and andesite in the upper parts of the deposits (Larsen & Thorarinsson, 1977;Sverrisdottir, 2007;Larsen et al., 2020). The most recent of the explosive silicic eruptions occurred in 1104 CE, with the eruption of 0Á9-1Á0 km 3 of tephra (Thorarinsson, 1967;Larsen & Thorarinsson, 1977;Janebo et al., 2016aJanebo et al., , 2016b. In addition to the silicic Plinian events, more than 10 explosive eruptions, which are zoned from andesite or dacite to basaltic andesite, erupted between 1050 BCE and 1104 CE (Larsen et al., 2020). In this study, we restrict all discussion to the historical eruptions, and therefore consider only the 1104 CE rhyolite eruption among the rhyolitic Plinian events in Hekla's history.
In a seminal work combining observations from the tephra record, the geology of the volcano, and historical records, Thorarinsson (1967) documented 13 eruptions of Hekla in the period following the 1104 CE rhyolite eruption to the 1947-1948 eruption. These eruptions occurred after repose periods ranging from 15 to 120 years, with an average repose of $64 years. Thorarinsson (1967) noted that most eruptions had a consistent pattern: (1) they begin with a highly explosive, presumably volatile-rich, initial phase of activity, which changes to effusive activity that may last weeks to months; (2) eruptions consistently exhibit compositional zoning in which the earliest erupted tephra is the most evolved and subsequent tephra and lavas are progressively less evolved. Thorarinsson (1967) also observed that the maximum SiO 2 content of the initial tephra is positively correlated to the length of the preceding repose period, and most eruptions terminated with a baseline composition of basaltic andesite (Fig. 2).

Geophysical constraints on Hekla's magmatic plumbing
Hekla has been monitored by modern geophysical methods since the mid-20th century, providing important constraints on the transport and storage of magma. Geirsson et al. (2012) used continuous and campaign GPS data between July 2000 (after the most recent eruption, which had a volume of 0Á19 km 3 ) and January 2010 to model the reinflation of the magma chamber. The geodetic data are best modeled by a sill-like chamber at 24 6 4 km depth. Subsidence on the flanks of the volcano further suggests a vertical finger-like chamber extending between 17-28 km deep at its base and 9-13 km at its top. Pre-eruptive inflation measured by InSAR is modeled as a 14-20 km deep magma chamber, inflating at 0Á003-0Á02 km 3 a -1 (Ofeigsson et al., 2011). Post-eruptive deflation in 2000 had a source centered at 14-18 km depth, and there is evidence of excess magma stored in a dike between 5Á8 km and the surface. One study that took place around the 2000 eruption and combined observations from InSAR, GPS, tilt, and a borehole strainmeter suggested that the depth of the pressure source is only about 10 km (Sturkell et al., 2013). These results also indicate that the dike supplying the eruptive fissure extends to <500 m depth, but that it does not connect along its entire strike to the deep magma reservoir.
Seismic studies reveal that Hekla magma does not slowly fracture its way to the surface preceding eruptions. Instead, Hekla displays a striking absence of precursory seismicity and strain. For example, seismicity began only 79 min before the initial breakout of magma in 2000, and a nearby strainmeter registered contraction, thought to be caused by ascent of the eruptionfeeding dike, 34 min earlier (Soosalu & Einarsson, 2005;Sturkell et al., 2013). Interestingly, the earthquakes' hypocenters appear to have propagated downward from 1-3 km to 12-15 km in the hour preceding eruption, which is not the signal expected from a dike rising from a deep magma reservoir. A similar phenomenon has been described during the 2010 CE eruption of Eyjafjallajö kull, where downward-propagating seismicity is attributed to a decompression wave caused by the initial tapping of a shallow magma body (Tarasewicz et al., 2012). In 1991, seismic activity began only 32 min prior to the eruption, with tremor detected only after the eruption commenced (Soosalu & Einarsson, 2002;Soosalu et al., 2003). Foci clustered in two zones, at 13-8 km and <4 km. The total seismic energy released over the course of the eruption only amounted to the equivalent of a magnitude 3Á4 earthquake. Soosalu & Einarsson (2004) used S-wave attenuation to conclude that no sizable body of magma exists between 4 and 14 km, indicating that any major reservoir must lie in the lower crust.

METHODS
In this study, we sampled all the historical tephras and lavas that we could confidently identify to construct a complete and internally consistent dataset with modern analytical tools. We sampled rocks from all the exposed historical eruptive units produced by Hekla, guided by geological maps (Thorarinsson, 1967;Hö skuldsson et al., 2007;Pedersen et al., 2019) and Thorarinsson's (1950Thorarinsson's ( , 1967 extraordinarily detailed stratigraphic columns from dozens of excavated pits and natural outcrops (Supplementary Data Electronic Appendix 1; supplementary data are available for downloading at http://www.petrology.oxfordjournals.org). To sample the widest compositional range for each historical eruption, care was taken to identify and collect tephra deposits along dispersal axes. Eruptive order was constrained by centimeter-scale stratigraphic documentation, extraction and analysis of single pumice fragments, and eruption narratives. We also report analyses of historically erupted basalts that surround the Hekla edifice, but we do not discuss their role in the petrogenesis [already thoroughly discussed by Sigmarsson et al. (1992a)].
All major elements are reported and plotted in units of per cent oxide by weight; they are normalized to 100 % but the analytical totals are reported. Trace elements are reported as parts per million by weight. Tephra clasts for major element analysis by X-ray fluorescence (XRF) were hand-crushed in an alumina mortar, and lava samples were chipped in a tungstencarbide jaw crusher before pulverization in a tungstencarbide shatterbox. Single pumice clasts were analyzed where possible (>2Á5 g). Major element analyses by XRF were undertaken at the Peter Hooper GeoAnalytical Laboratory at Washington State University (Johnson et al., 1999;Supplementary Data Electronic Appendix 2). Unground glass handpicked from the exterior rinds of individual tephra clasts and ground lava samples was dissolved for trace element analysis by inductively coupled plasma-mass spectrometry (ICP-MS) on a Varian 820MS quadrupole ICP-MS system at Colgate University, using methods described by Harpp et al. (2003) (Supplementary Data Electronic Appendix 2). Reproducabilities for XRF and ICP-MS are estimated from separate analyses of multiple tephra clasts from the same horizon in tephra fall and are reported in Supplementary Data Electronic Appendix 2 and shown in the variation diagrams.
Major element analyses of glass, including glassy melt inclusions, were determined on Cameca Camebax and JOEL 8500F microprobes at Washington State University (Supplementary Data Electronic Appendix 3). Glass analyses were performed with a beam voltage of 15 kV, a 5-10 mm spot size, and a beam current of 10 nA. Sodium loss was minimized by measuring its abundance immediately after exposing the sample to the beam and correcting for time-dependent volatilization. Precisions based on replications of analyses are reported in Supplementary Data Electronic Appendix 3.
Mineral modes were determined by counting $1000 points on an automated stage mounted on a petrographic microscope. Modes are reported as the proportion of crystals >0Á2 mm in major dimension without interpretation of genetic origin. Most Hekla rocks have <6 % crystals, thus mineral concentrates were mounted in epoxy disks for electron micropobe analysis. Crushed rocks were sieved to 0Á35-1Á0 mm and, after separation of the magnetic fraction, treated with concentrated HBF 4 to remove adhering groundmass; optical and back-scatter electron imaging of grain boundaries and zoning indicates that the phenocrysts were not corroded by this treatment, only glass and fine-grained material. The grains were then mounted in epoxy and analyzed on three microprobes: a Cameca Camebax (Washington State University; WSU) and 2 JOEL 8500 instruments (WSU and University of Hawaii). Beam voltage was set to 15-20 kV, current to 15 nA, and the beam focused to $1 mm for mineral analyses.
Oxygen isotopic analyses were carried out at the University of Oregon by laser fluorination using BrF 5 as a reagent and a MAT253 mass spectrometer according to the methods of Bindeman et al. (2008b;Supplementary Data Electronic Appendix 4). Owing to their small grain size, many analyses were performed on 2-14 olivine grains. In one basalt sample and the 1104 CE rhyolite tephra, olivine is large enough to permit analysis of single grains. Many olivine grains contain magnetite inclusions, which were culled with a magnet for selecting samples for oxygen isotopic analysis. Data were normalized to SCO ¼ 5Á25 & and UWG garnet ¼ 5Á80 &. Precision is estimated to be 60Á1 & (one standard deviation) on the basis of concurrently run standards.
Naturally quenched melt inclusions (MIs) in olivine were doubly intersected for analysis by Fourier transform infrared spectroscopy (FTIR) at the University of Oregon (Supplementary Data Electronic Appendix 5). Only olivine from tephra samples was used, with the exception of sample H-07-1, which is a glassy spatter bomb, to ensure rapid quenching and to minimize the effect of post-eruptive hydrogen loss that occurs in slowly cooled lava. Infrared spectra were collected between 6000 and 1000 cm À1 , and the main peaks of interest for this study were the total OH stretching vibration at 3570 cm À1 and the carbonate doublet at 1435 and 1515 cm À1 , although carbonate was detected in only two melt inclusions. Two to four replicate spectra were acquired for each MI. Thickness was measured using both a digital micrometer (62 lm) and the reflectance interference fringe method described by Wysoczanski & Tani (2006). Absorbances are converted to H 2 O concentrations using the Beer-Lambert Law. We used an absorption value of 63 6 3 l mol -1 cm -1 for H 2 O (Dixon et al., 1988). An absorption coefficient of $350 l mol -1 cm -1 was used for the CO 2 doublet, calculated using the glass composition (Dixon & Pan, 1995). Peak heights were calculated using a straight-line background correction , and glass densities, including the effect of H 2 O, were calculated from the major element compositions as described by Luhr (2001). Combined uncertainties in density, absorption coefficient, and thickness cause precision to be $10 % (Dixon et al., 1988), where most of this error is due to the thickness measurement. Consequently, the average standard deviations (1r) for H 2 O are 60Á25 wt% (one standard deviation).
Apatite grains in a basaltic andesite and two andesites from the 1947 eruption were analyzed by laser ablation (LA-)ICP-MS in the Geoanalytical Laboratory at Washington State University, using a New Wave 213 nm laser and Element2 mass spectrometer (Supplementary Data Electronic Appendix 10). Ablation troughs were $10 lm wide and tens of micrometers long with an irregular geometry. Counts per second are internally normalized to stoichiometric CaO contents and measured 43 Ca counts per second, with BCR-2g and NIST610 glasses as external standards. Relative precision based on replicate analyses is estimated to be 10 % for Th and 11 % for U.

Major elements
Our bulk-rock major element data from historical Hekla tephra and lavas (Supplementary Data Electronic Appendix 2) are consistent with those reported in other studies (Baldridge et al., 1973;Sigmarsson et al., 1992a;Chekol et al., 2011). The volcano is constructed of a transitional subalkaline suite, with whole-rock compositions ranging from basaltic andesite through rhyolite (54-70 % SiO 2 ; all oxides are reported as weight per cent). Rocks with <54 % SiO 2 have not erupted from the Hekla edifice in historical time, but basaltic material has been produced by nearby fissures as recently as 1878 CE and1913 CE (Thorarinsson, 1967). In our study, any rock referred to as a Hekla rock erupted from vents at the summit or along the flanks of Hekla volcano; we exclude rocks erupted from surrounding fissures and monogenetic vents that are not directly on the edifice.
Hekla lava and tephra exhibit remarkably uniform major element trends at both the inter-eruption and intra-eruption level (Fig. 3). Variations in MgO, TiO 2 , FeO, CaO, and P 2 O 5 all define negative correlations with SiO 2 , whereas K 2 O and Na 2 O correlate positively with SiO 2 . There is a strong inflection in Al 2 O 3 at $60 % SiO 2 (Fig. 3). Many of the elements, especially TiO 2 , Al 2 O 3 , and P 2 O 5 , do not define linear trends. Silica contents of glasses tend to be slightly more evolved than the bulkrocks (Fig. 4), especially the 1104 CE rhyolite. Otherwise, in comparison with whole-rock compositions, glasses tend to be poorer in elements that are concentrated in the major phenocryst phases or microlites (Fig. 4).
The 1104 CE and 1158 CE rocks are exceptional as the sole historical eruptions of Hekla to produce only silicic material, with no known late phase of basaltic andesite lava (Thorarinsson, 1967;Sigmarsson et al., 1992a;Larsen et al., 1999;Fig. 3). The 1104 tephra is relatively homogeneous at the two localities we sampled; bulk-rock analyses of individual pumice clasts vary from 70Á0 % SiO 2 in the bottom third of the deposit to 68Á9 % SiO 2 in the top third of the deposits (on a normalized, volatile-free basis). Glasses from the 1104 deposits range from 71Á0 to 72Á5 % SiO 2 (Fig. 4).
Pumice clasts from the 1158 CE eruption have uniform dacitic ($65 % SiO 2 ) compositions despite a range of color and are similar to analyses reported by Larsen et al. (1999). The Haahraun lava flow on the south side of Hekla, now exposed as a kipuka in the 1991 CE lava flow field (Fig. 1), is the only known silicic lava on Hekla and is thought to have also erupted in 1158 (Sigmarsson et al., 1992a).
The most recent eruptions with substantial compositional zoning are the 1845 CE and the 1947-1948 CE events, which occurred after repose periods of 77 and 102 years. Well-preserved deposits of the 1845 and 1947 tephra layers are characterized by a brown to black color transition from the bottom to the top of the tephra deposits, corresponding to a compositional change from $63 to 59 % SiO 2 (Fig. 5). Lavas emitted during the waning phase of the eruptions have a minimum of $54 % SiO 2 . The 1947-1948 eruption produced materials with a compositional range nearly as large as that of all compositions erupted since the 1158 CE dacite (Fig. 3).
Following the 1947-1948 CE eruption, Hekla erupted again after a repose of only 22 years. The products from the 1970 CE eruption have a much narrower compositional range, $1 % SiO 2 difference between tephra phase and last-erupted lava. Chekol et al. (2011) and Thorarinsson & Sigvaldason (1972) reported analyses of silicic pumice from the 1970 deposits, and Sigmarsson et al. (1992) reported a silicic ignimbrite xenolith. We observed white frothy inclusions in 1970 tephra, but, as discussed below, we interpret these to be fused rocks that were incorporated into the 1970 magma very near the surface and are volumetrically unimportant. Since 1970, Hekla has shifted to more regular eruptions, with repose periods of only 9-10 years, although the volcano has been quiet for more than 20 years since the 2000 CE eruption (Gronvold et al., 1983;Hö skuldsson et al., 2007;Gudnason et al., 2018). The products of these most recent eruptions exhibit minimal compositional variability compared with the previous historical eruptions (Fig. 3). The most mafic lavas analyzed in this study are from the last phases of the 1947, 1970, 1980, 1991, and 2000 eruptions, which have SiO 2 contents approaching 54 %.

Phenocrysts
Hekla lava and tephra are notably crystal-poor and have an anhydrous assemblage consisting of plagioclase þ olivine þ clinopyroxene þ magnetite þ apatite (irrespective of bulk-rock composition) 6 orthopyroxene 6 ilmenite 6 zircon. The basaltic andesites have a maximum crystallinity of 9 % (by volume), and crystallinity decreases with increasing SiO 2 content to the andesites (Supplementary Data Electronic Appendix 3; Fig. 6). In contrast, the silicic andesites, dacites, and rhyolites are richer in crystals. The relationship between crystallinity and composition is also reflected in zoned eruptions.
For the 1947-1948 event, the first-erupted, most evolved andesite contains $3 modal% crystals, which decreases to nearly 0 % in the upper part of the tephra deposit (Sigvaldason, 1974). This finding is also consistent with crystallinity data from the 1104 CE rhyolite tephra and prehistoric H3 and H4 tephra, which also exhibit a correlation between crystal content and composition (Sigvaldason, 1974;Sverrisdottir, 2007). The size of mafic phenocrysts also correlates with rock composition. Olivine phenocrysts in basaltic andesites vary from 0Á3 to 0Á5 mm; pyroxene phenocrysts are from 0Á5 to 1 mm. Olivine and pyroxene crystals in the andesite, dacite, and rhyolite samples are twice these sizes. Basaltic andesites and early erupted andesites have euhedral phenocrysts characterized by sharp, well-defined grain boundaries, whereas phenocrysts in latererupted, crystal-poor andesite tephra and lava have subhedral to anhedral habits. Most pyroxene and plagioclase in the 1104 CE rhyolite is subhedral, although most faces are faceted (Supplementary Data Electronic Appendices 7 and 8). Many grains have concentric growth zones. Some grains enclose inclusions of apatite, zircon, magnetite, and pyrrhotite. The higher-Mg pyroxene in the 1158 CE dacite tephra is mostly euhedral, although many grains have inclusion-riddled and corroded lower-Mg cores. Some of these grains have zircon, magnetite, and pyrrhotite inclusions. Plagioclase is also faceted in the dacite and much is intergrown with anhedral olivine and pyroxene in glomerocrysts.
Plagioclase is the dominant phenocryst phase in all Hekla rocks, typically constituting between 50 and 80 % of the phenocryst assemblage. Most rocks contain clinopyroxene, and orthopyroxene is present in some andesite and basaltic andesite. Orthopyroxene is relatively abundant in the upper parts of the 1845 CE and 1947 CE andesite tephra deposits, where it is more abundant than clinopyroxene. Sparse orthopyroxene Fig. 5. Compositions of rocks, matrix glasses, melt inclusions, and mineral phases change consistently over the course of the 1947 eruption. Eruptive order is lowest tephra, middle tephra, upper tephra, main lava phase, waning lava phase. Analytical precision is indicated by horizontal bars on each of the plots where it is greater than the size of the symbols. also occurs in the 1389 CE and 1970 CE tephra, as well as late-stage spatter from the 1947-1948 eruption.
Magnetite phenocrysts are ubiquitous in Hekla rocks and occur as solitary phenocrysts, in glomerocrysts, and abundant inclusions in olivine and pyroxene phenocrysts. Olivine-magnetite intergrowths are especially abundant in basaltic andesites. Ilmenite occurs as sparse $100 mm euhedral inclusions enclosed within, or closely associated with, magnetite, clinopyroxene, and olivine phenocrysts. All rocks contain apatite, typically as 50-100 mm euhedral microphenocrysts along the margins of and as inclusions in magnetite, olivine, and clinopyroxene phenocrysts and, to a lesser degree, plagioclase. Apatite inclusions are larger and constitute a greater proportion of the assemblage in the basaltic andesites than in more evolved compositions. Zircon is a trace phase in the 1104 CE rhyolite tephra and 1158 CE dacite tephra but is otherwise absent. It occurs as 20-100 mm prismatic microphenocrysts on the edges of and as inclusions within clinopyroxene, magnetite, and olivine. Pyrrhotite is present as a trace phase throughout the Hekla suite. It occurs as 5-20 mm anhedral inclusions in magnetite, olivine, and pyroxene phenocrysts, as well as some plagioclase crystals from the 1947 CE andesite.
Olivine compositions range in composition from Fo 11 to Fo 59 and correlate with rock compositions and form a trend that is sub-parallel to the equilibrium olivine field predicted by standard olivine-liquid Fe-Mg partitioning models, but at lower Fo contents in the more silicic rocks (Rhodes et al., 1979;Putirka, 2008;Fig. 7a). Experimental data on iron-rich compositions indicate that the K D Fe/Mg increases to 0Á41 for olivine of composition Fo 27 (Toplis & Carroll, 1995;compare Fo 11 in the 1104 CE rhyolite), which fits the Hekla data better than the standard value of 0Á30.
One 1766 CE andesite contains sparse pl-px-ol-mt glomerocrysts with anomalously low-Fo olivine (Fo 29-31 ), compositions similar to some olivine from the 1158 CE dacite. Euhedral olivine phenocrysts within the host lava have Fo 37-39 (Fig. 7a). The 1104 CE rhyolite and 1158 CE dacite tephra have a restricted range of olivine compositions, and the crystals are unzoned (Fig.  7b). Most olivine grains in the first-erupted andesite tephra from the 1845 and 1947 zoned eruptions are compositionally uniform, whereas olivine grains from the later-erupted tephra are reversely zoned, with spatially broad zoning (Fig. 7c). Basaltic andesite lava and spatter from the final phase of the 1947-1948 eruption is unzoned or has olivine with narrow, normally zoned rims (Fig. 7b). Basaltic andesite tephras from the last Fig. 6. The crystallinity of Hekla lavas is strongly related to rock composition, with a crystal content minimum observed in the silicic andesite. Crystal content of the 1104 CE rhyolite is from Sigvaldason (1974). Dacites are the most crystal-rich rocks and contain a large proportion of xenocrysts. four Hekla eruptions contain unzoned olivine similar to the olivine from late in the 1947-1948 eruption (Fig. 7b).
Pyroxene compositions also correlate with rock composition but are more complexly zoned than olivine (Fig. 8). Pyroxenes in the basaltic andesites exhibit no regular zoning pattern and have sectors that vary by <2 % in Mg# (Fig. 8c). Pyroxene from low-SiO 2 andesites from the upper parts of the 1947 and 1845 tephra have $50 mm reversely zoned rims that trend to compositions of augite from the basaltic andesites. One 1766 CE andesite lava sample contains sparse clinopyroxene in pl-px-ol-mt glomerocrysts, which have orthopyroxene reaction rims and are similar in composition to pyroxenes in dacites ( Fig. 7; H05-08). Pyroxene in the 1158 CE dacite is both normally and reversely zoned, but all measured grains have similar rim compositions (Mg# $43; Fig. 8b). The most Fe-rich cores in dacite pyroxene have the same composition as the rims and cores of pyroxene in the 1104 CE rhyolite. These evolved cores contain abundant glassy inclusions. The 1104 CE rhyolite contains pyroxene that has either oscillatory zoning ($5 % variation in Mg#) or normal zoning across the entire grain, but rim compositions converge to Mg# $27 (Fig. 8b).
The compositions of plagioclase phenocrysts also correlate strongly with the composition of their host rocks (Fig. 9a). Most of the plagioclase grains in the basaltic andesites are normally zoned, with cores of An 54-60 and rims extending to An 40 . Plagioclase in the 1158 CE dacitic tephra has oscillatory zoning, ranging in molar anorthite by $4 % (Fig. 9b). Plagioclase in the 1104 CE rhyolite is broadly normally zoned, with rims up to 300 mm thick (Fig. 9b). Unlike the reversely zoned pyroxene, plagioclase in the early erupted andesitic 1947 tephra lacks zoned rims, but the cores are oscillatory zoned from An 51 to An 56 .

Melt inclusions
We analyzed glassy olivine-hosted melt inclusions in basaltic andesite tephra from the 1947, 1970, and 2000 CE eruptions, the zoned eruptions from 1510, 1845, and 1947 CE, and the 1104 CE rhyolite. Major element compositions of the olivine-hosted melt inclusions are similar to those of the whole-rock compositions, the same finding as in earlier studies (Moune et al., 2007;Portnyagin et al., 2012;Lucic et al., 2016), although melt inclusions in olivine of similar composition have a small range in Mg#, suggesting post-entrapment crystallization. For most rocks, the range in glass Mg# is <5 and can be accounted for by <10 % crystallization of olivine. Because of uncertainty in the Fe-Mg exchange coefficient for fayalitic olivine and the small amount of postentrapment crystallization, we do not correct the melt inclusion compositions (Figs 4 and 5). In contrast to the Hekla suite, melt inclusions in olivine from primitive Icelandic basalts tend to be compositionally heterogeneous, differing from the carrier melt, which has been interpreted as entrainment of a cumulate mush that was precipitated from compositionally distinct melts (e.g. Maclennan, 2008). In general, melt inclusions in olivine from the basaltic andesites tend to be slightly poorer in SiO 2 than the carrier melts ( Fig. 4a), but the absence of any basaltic melt inclusions is notable. The melt inclusions in the 1104 CE rhyolite olivine are more silicic than the whole-rock but are similar in composition to the carrier melt ( Fig. 4a-d).
H 2 O concentrations in melt inclusions increase with SiO 2 and incompatible element concentrations, especially the maximum concentrations from any given deposit (Fig. 10). This follows the pattern described by The sulfur concentrations in melt inclusions from the basaltic andesites agree well with those measured by Moune et al. (2007) and match the sulfur concentrations of pyrrhotite-saturated Hekla melts determined experimentally (Moune et al., 2009). The sulfur concentrations also agree well with the predicted sulfide solubility for models with a FeS ¼ 1 (Smythe et al., 2017). Although the sulfur concentrations in the melt inclusions in the andesites are 100-150 ppm richer in sulfur than the thermodynamic model predicts, they are $100 ppm poorer in sulfur than experimental glasses with 9 % FeO* that are pyrrhotite-saturated (Moune et al., 2009).

Trace elements
Rb is a highly incompatible element in the Hekla system; other elements that are typically incompatible in mafic systems are affected by the observed phenocryst phases, including apatite, ilmenite, and ferroaugite, as well as zircon in the most evolved rocks. Barium, U, Pb, and Th correlate almost perfectly with Rb (Fig. 11); the slope and intercept on these plots indicate that Ba, U, Pb, and Th are slightly less incompatible than Rb, if they are related by crystallization. Yttrium (and the heavy  rare earth elements; HREE) forms three distinct linear trends ( Fig. 11): tephra from the four late 20th century eruptions have distinctly higher Y at a given Rb concentration in comparison with the pre-1970 CE eruptions. Also, the 1158 CE dacite and 1104 CE rhyolite are offset to low Y concentrations. Strontium, much like Al 2 O 3 , correlates with Rb then inflects to define a negative relationship in the andesites. The tephra are systematically richer in Sr than the lavas for a given Rb concentration (Fig. 11). Zirconium and Nb correlate positively with Rb, except for the 1104 CE rhyodacitic tephra, which have anomalously low Nb and Zr concentrations; in contrast, Fig. 11. Trace element concentrations determined by ICP-MS plotted against Rb, the most incompatible element in the Hekla suite. It should be noted that linear trends of Zr, Nb, La, and Ba intersect the y-axis, indicating that element is less incompatible than Rb. Lava and tephra samples are separated by a dashed line: for many of the elements, tephra clasts have distinct trace element compositions compared with the lava. Zirconium variations in 1104 CE lavas indicate zircon saturation and removal from the system. the 1158 CE dacite is unusually rich in these elements. The trend of the 1947 zoned eruption is offset to slightly low Zr and Nb concentrations, especially at intermediate compositions. Uranium and Th have a strongly coherent linear trend, with an r 2 of 0Á99 (Fig. 12). Because the linear array has a non-zero intercept, the U/Th changes with the extent of evolution.
Owing to the large number of analyses, we apply ratios of REE rather than the standard REE pattern (Fig.  13). The ratio Dy/Yb is a measure of the fractionation across the middle REE (MREE) to HREE, and Dy* (Davidson et al., 2013) quantifies the concavity of the REE pattern; Dy* values less than unity indicate patterns that are concave-up. In Hekla rocks, Dy/Yb and Dy* are strongly correlated. The patterns are increasingly concave-up with evolution, whereas Dy/Yb are greater than one and decrease with evolution, indicating decreasing MREE-HREE fractionation. Hekla rocks also exhibit a correlation between La/Sm and La (Fig.  13b), and La/Sm increases steadily with evolution of the Hekla suite. Curiously, the strongly zoned 1947 CE and 1845 CE tephra and lava are systematically shifted to higher La/Sm at a given La concentration.
Compatible trace elements have complex patterns, many of which define two trends; samples from 1970 to 2000 CE exhibit a different slope from the rest of the historical suite (Fig. 11). For example, Sc variations have an overall negative slope, consistent with a system crystallizing ferroaugite. In contrast, rocks from 1970 to 2000 have a positive slope, indicating incompatible behavior by Sc. In fact, this correlation is largely due to tephra samples systematically having greater Sc concentrations for a given Rb concentration. Likewise, V decreases exponentially with respect to Rb, but the tephra samples are systematically offset to higher Rb and V than the lava samples.

Oxygen isotopes
As is observed in many igneous rocks from Iceland (e.g. Bindeman et al., 2008b), olivine from Hekla rocks has 18 O/ 16 O below typical mantle values (Fig. 14). Moreover, individual olivine grains display within-sample heterogeneity. These olivine separates compare with wholerock values at Hekla of 4Á7-5Á1 & ( Fig. 14; Sigmarsson et al., 1992a). Single crystals of olivine from the 1104 CE rhyolite range from 3Á4 to 4Á2 6 0Á1 &, but there is no systematic change over the course of the eruption. The 1158 CE dacite has similar d 18 O values to the andesites. Olivine from the strongly zoned 1947 CE eruption ranges from 3Á5 to 4Á6 6 0Á1 & as the eruption proceeded from andesite to basaltic andesite, but, again, there is no systematic change over the course of the  We also analyzed olivine from several historical and Holocene basalts erupted in the vicinity of Hekla, which range from 3Á7 to 4Á7 6 0Á1 & (Supplementary Data Electronic Appendix 4; Fig. 14). The most notable observation is that olivine in the most primitive basalts has the highest d 18 O, and d 18 O decreases strongly with differentiation while still in the basaltic realm, a trend similar to that for Laki and other large-volume basaltic systems in Iceland (Bindeman et al., 2008b).

Zoned eruptions
As noted by Thorarinsson (1967), many Hekla eruptions are zoned, and the extent of the zoning depends on the repose between eruptions (Fig. 2). We focus on the 1947 CE eruption as a type-example (Fig. 5), although the 1845 eruption exhibits similar compositional diversity. The first erupted tephra is highest in SiO 2 , with 62Á4 % SiO 2 , corresponding to a glass composition of 63Á6 % SiO 2 . Thirteen melt inclusions trapped in 10 different olivine crystals have between 61 and 64 % SiO 2 , with an average composition of 62Á5 % SiO 2 , virtually the same as the rock. Over the course of the eruption, whole-rocks and glasses steadily decrease to 54Á3 % SiO 2 ; the glassy melt inclusions have practically the same composition as the whole-rocks and glasses. Melt inclusions from the upper tephra, which is also poor in crystals, are more compositionally variable.

DISCUSSION
A wealth of geochemical data have been published on Hekla, which have been interpreted primarily in two ways. One line of hypotheses is that crystal fractionation of basaltic parent magma creates a range of intermediate to silicic magma compositions (Baldridge et al., 1973). The second class of models contends that the silicic magmas contain a large component of crustal rocks. In the latter models, intermediate compositions are attributed to magma mixing of end-member compositions, a combination of crystal fractionation and magma mixing, or a combination of assimilation and fractional crystallization (Sigvaldason, 1974;Sigmarsson et al., 1991Sigmarsson et al., , 1992aSverrisdottir, 2007;Martin & Sigmarsson, 2010;Chekol et al., 2011;Lucic et al., 2016). The comprehensive dataset presented here permits evaluation of these hypotheses, particularly given the constraints of timing and eruptive order.

Conditions of magma storage Magma storage temperatures
Equilibration temperatures and oxygen fugacities from Fe-Ti exchange in magnetite-ilmenite pairs using the calibration of Ghiorso & Evans (2008) are calculated after verifying equilibrium using the Mg-Mn partitioning test of Bacon & Hirschmann (1988). Temperatures range from 1011 6 5 C for basaltic andesite (single oxide pair from 1970 tephra) and 953 6 14 C for andesite (three oxide pairs from 1947 tephra). Analyses of an oxide pair included within a clinopyroxene phenocryst from the 1510 CE andesite tephra yield a temperature of 897 C. This phenocryst also contains zircon inclusions, which otherwise are present only in dacite and rhyolite at Hekla. Calculated oxygen fugacities for all magnetite-ilmenite pairs fall within 0Á25 log units of the fayalite-magnetite-quartz (FMQ) buffer, except for the slightly more oxidized [log(fO 2 ) ¼ FMQ þ 0Á5] results for the 1970 basaltic andesite oxide pair. These temperatures and oxygen fugacities compare with an average of 870 6 26 C (fO 2 ¼ 0Á4 log units more reduced than FMQ) for the prehistoric H3 tephra (Portnyagin et al., 2012).
We apply the two-pyroxene geothermometer of Putirka [2008;equation (36)] to all pyroxene pairs from each rock for which we have data. Pairs that are not in Fe-Mg equilibrium are discarded and the remaining pairs are averaged. The only pyroxenes found to pass this equilibrium test are in the basaltic andesites and andesites, and they yield temperatures of 1003 6 8 C and 978 6 14 C, in close agreement with the oxide thermometer (Fig. 15). Forty-three clinopyroxene analyses from basaltic andesites are calculated to be in Fe-Mg equilibrium with their host glass or whole-rock (K D ¼ 0Á27 6 0Á3; Fig. 8), and application of the clinopyroxene-melt geothermometer [Putirka, 2008, equation (33)] yields T ¼ 1029 6 15 C. Forty-seven pyroxene analyses from andesite have FeO and MgO concentrations that meet the equilibrium test and yield T ¼ 952 6 20 C. Application of a more stringent filter involving four separate exchange equilibria removes 23 pyroxene analyses from the basaltic andesites and all but two analyses Oxygen isotope ratios in olivine crystals separated from Hekla rocks and basalts erupted in the vicinity of Hekla compared with bulk-rock data of Sigmarsson et al. (1992a).
Large symbols indicate composites of 2-14 grains; smaller symbols are single-grain analyses. Illustrated uncertainty is 2r. Arrow shows predicted trend for olivine owing to fractional crystallization: although melt d 18 O increases owing to Fe-Ti oxide fractionation, the fractionation factor between olivine and silicic melt is larger in magnitude. from the andesites. Calculated temperatures of this filtered subset are 1030 6 14 C and 952 to 984 C for the basaltic andesites and andesites.
The apatite saturation thermometer (Harrison & Watson, 1984) applies to all Hekla rocks, owing to the ubiquitous presence of apatite phenocrysts. Calculated temperatures for basaltic andesite, andesite, dacite, and rhyolite are 1043, 997, 933, and 877 C (Fig. 15). Only the 1104 CE rhyolites contain zircon in contact with glass, and it yields a saturation temperature of 878 6 6 C (Boehnke et al., 2013), almost precisely the same as the other geothermometers (Fig. 15). Carley et al. (2011) applied the zircon saturation geothermometer to 1104 CE rhyolite tephra, which yields 920 C (we note that their interpretation was that their sample is an 1158 CE rock, on the basis of the glass composition, whereas we believe it is a 1104 CE rhyolite, which is confirmed by Calvin Miller, personal communication, 2020). The Ti-inzircon geothermometer in the same rocks results in temperatures ranging from 750 to 930 C.
In summary, a series of geothermometers yield consistent results, with the most mafic basaltic andesites erupting between 1040 and 1000 C, the andesites ($60 % SiO 2 ) between 1000 and 950 C, the 1158 CE dacite from 950 to 900 C, and the 1104 CE rhyolite from 890 to 870 C (Fig. 15).

Magma storage pressures
Our interpretation of the melt inclusion H 2 O concentrations is that the maximum values in any given sample represent the magmatic concentrations at the site where the crystals grew prior to eruption, owing to the excellent correlation of the maximum H 2 O concentrations with incompatible element contents (Fig. 10). The same interpretation has been made by Moune et al. (2007) and Portnyagin et al. (2012). Lower values are interpreted to be either the result of entrapment during ascent and degassing of the magma or diffusive hydrogen loss from the melt inclusions as the host melt dehydrates. If the magmas were undersaturated in H 2 O-rich vapor as the crystals were forming, calculated vapor saturation pressure is a minimum estimate of the trapping pressure for the inclusions.  Papale et al., 2006). Another solubility model (Ghiorso & Gualda, 2015) does not converge for the mafic and felsic melts, but the calculated saturation pressure for the andesite melt is 212 MPa. The H 2 O solubility model of Moore et al. (1998) yields pressures of 215, 140, and 65 MPa for the rhyolite, andesite, and basaltic andesite (Fig. 15). Given these very shallow depths relative to the geophysically determined depth to the magma reservoir and the mineral-melt geobarometric results reported below, it is clear that the Hekla magma cools and evolves in the absence of a water-rich fluid. Major vesiculation probably does not occur until the magma is at very shallow depths, <2 km for the basaltic andesite.
We screened 143 clinopyroxene-melt pairs in the basaltic andesites and andesites for Fe-Mg, CaTs, EnFs, and DiHd equilibrium (bound set at 620 %) and found 90 analyses to be in Fe-Mg equilibrium with their host glass (or whole-rock for samples lacking glass), and 22 analyses to be in equilibrium for all four of these exchange components. H 2 O concentrations of melts are estimated on the basis of constant K 2 O/H 2 O, and we apply Putirka's (2008) equation (31). Kernel density estimation (with a Gaussian filter) yields a maximum density of 570 6 130 MPa for the basaltic andesites and 520 and 540 MPa for the andesites, when the more stringent filter is applied (interquartile range; total range is 490-710 MPa). When only the Fe-Mg filter is applied, estimated pressures are 510 6 110 MPa for the basaltic andesites and 400 6 130 MPa for the andesites, consistent with the depth to the main magma reservoir as determined by surface deformation since 2000 CE (Geirsson et al., 2012). Putirka's (2008) equation (30) yields pressures $100 MPa lower. Thus, the clinopyroxene compositions are consistent with the basaltic andesite magmas partially crystallizing at depths of at least 12 km and probably >17 km, in the lower Icelandic crust. The large number of crystals that are not in equilibrium according to these four exchange reactions might be due to rapid crystal growth relative to diffusion rates, or alternatively the crystals grew from melts that no longer host them and they were picked up by the carrier melts as magma ascended over the course of ascent and eruption. Moreover, the exchange reactions are calibrated mostly on basalts and may be imperfectly calibrated for the evolved Hekla compositions.
H 2 O concentrations of the melts can also be estimated with the plagioclase hygrometer (Waters & Lange, 2015). Average plagioclase and glass (where available) or whole-rock compositions are applied, along with temperatures from geothermometers and a pressure of 0Á5 GPa. The hygrometer is only weakly dependent on temperature and pressure (Waters & Lange, 2015). The plagioclase compositions agree with the melt-inclusion data and indicate a steady increase in H 2 O with magmatic evolution (Fig. 10) and an excellent correlation with K 2 O, but at values that are 0Á5-1Á0 % H 2 O lower than the maximum melt inclusion concentrations. The cause of this systematic difference is unclear. Waters & Lange (2015) reported an uncertainty in their calibration of 60Á35 % H 2 O. Uncertainty in temperature and pressure estimates results in increasing the uncertainty by an additional 0Á2 %. The congruence of the plagioclase-melt equilibrium calculations and the melt inclusion measurements confirms that melt H 2 O correlates strongly with incompatible elements in the Hekla suite, however.
Melting experiments on H3 tephra rhyodacite, which is a prehistoric Plinian deposit (about 1000 BCE) similar in composition to the 1104 CE tephra studied here, show the melt to be saturated with olivine, plagioclase, clinopyroxene, and vapor at approximately 850 6 15 C and P H2O ¼ 130-175 MPa (Weber & Castro, 2017). This implies storage in the shallow crust.

Magmatic evolution at Hekla
We approach the problem of the petrogenesis of Hekla's diverse suite by breaking the problem into three parts: the relationship between the basaltic andesites and andesites (including the zoned eruptions), the significance of low-d 18 O values, and the origin of the dacite and rhyolite.

Basaltic andesites to andesites (54-62 % SiO 2 )
The strikingly coherent relationships among the major elements across Hekla eruptions suggest that crystallization is probably the dominant process relating the basaltic andesite to the andesites (Fig. 3). Basaltic andesites at Hekla contain the phenocryst assemblage pl þ aug þ ol þ mt þ ilm þ ap. We apply a least-squares mass-balance method to determine whether these phases can account for the chemical variation of Hekla magmas with 54-62 % SiO 2 , Our results indicate that the eruptive products of the 1947 eruption, which range from 54Á7 % (H07-01) to 62Á5 % (H05-18) SiO 2 can be explained by fractionation of 4Á2 % olivine, 12Á8 % plagioclase, 6Á1 % augite, 4Á0 % magnetite, 2Á0 % apatite, and 1Á1 % ilmenite, with a sum of squares of residuals of 0Á02; mineral compositions are averages from the parental basaltic andesite (Supplementary Data Electronic Appendices 6-9). This amounts to a total of 30 % crystallization, and according to the geothermometers, it is driven by cooling of the melt of 50 6 20 C.
As noted by others (Chekol et al., 2011), the basaltic andesite to andesite compositional changes can be modeled by the MELTS algorithm. The best-fit models are those at lower-crustal pressure conditions (400-500 MPa) and H 2 O concentrations of 0Á4 % in the parental basalt. We also modeled the liquid line of descent of the 1947 basaltic andesite, 1947 andesite, and 1158 dacite using H 2 O concentrations defined by the melt inclusion K 2 O-maximum H 2 O correlation. Models run at pressures of 200-600 MPa using rhyolite-MELTS (Gualda et al., 2012) closely match the major element trends, with the exception of Al 2 O 3 and P 2 O 5 . The highest Al 2 O 3 concentrations at Hekla are 15Á4 %, in andesites with $62 % SiO 2 , and these rocks contain plagioclase phenocrysts. Rhyolite-MELTS predicts that the inflection in Al 2 O 3 is at 16Á3 % at 200 MPa or 16Á8 % at 600 MPa, alumina concentrations that exceed those in any Hekla rocks, most of which contain plagioclase phenocrysts. Rhyolite-MELTS also does not predict saturation with apatite in Hekla andesite and dacite, despite the rocks containing euhedral apatite phenocrysts and having the compositional fingerprint of having crystallized apatite. Villiger et al. (2004) and Balta & McSween (2013) noted that MELTS is successful at predicting experimental liquid lines of descent to only within 1 wt% for most oxide components, and with that criterion, the Hekla suite is well modeled by rhyolite-MELTS with the caveat that the model underestimates the temperature of plagioclase saturation at these higher pressures and H 2 O concentrations.
The coherent linear trends between incompatible elements are also consistent with crystallization processes to explain the relationship between basaltic andesite and andesite (Fig. 11). The array of data between compatible and incompatible elements (e.g. Rb versus V, Sc, and Sr) are curved, also consistent with crystallization. To model all of the trace elements quantitatively, we adopt the partition coefficients in table 4 of Nielsen (1992), supplemented by more recent data from the Geochemical Earth Reference Model (https://earth ref.org/GERM/) (Supplementary Data Electronic Appendix 7). Mineral and melt fractions determined by the least-squares major element model are applied to calculate the bulk partition coefficient, and the Rayleigh equation is used. The concentrations predicted by this method match the incompatible element contents in Hekla lavas closely (Fig. 16), certainly within the uncertainty in the partition coefficients. The poorest fits are for U and Th (see below) and the MREE to HREE, both of which we attribute to poorly determined partition coefficients for apatite, which are probably higher in the Hekla suite than those adopted here (see Prowatke & Klemme, 2006). Also, the distribution coefficients for HREE are known to be substantially higher for iron-rich pyroxene (Olin & Wolff, 2010) than for the compositions we used.
Mineral compositions in Hekla lava and tephra support the fractional crystallization model. With few exceptions, olivine, pyroxene, and plagioclase solid-solution compositions correlate strongly with their host rock compositions (Figs 7-9). This would not be the case if the intermediate rocks resulted from magma mixing (assuming some crystals survived hybrization), or if the andesites were formed by crustal melting or wall rock assimilation (assuming that some crystals would have derived from the protolith or wall rock). With few exceptions, the cores of most crystals have only minor oscillatory zoning, fluctuating by only small variations in Mg# or An (Figs 8, 10, and 12), and the most common style of zoning is narrow, normally zoned rims.
Melt inclusions in olivine in the basaltic andesites tend to have the same composition as the matrix glass or to be slightly more mafic (<2 % SiO 2 ), which is consistent with the olivine crystallizing from a slightly hotter ancestor of the host melt, whose composition shifted owing to crystallization of the phenocrysts themselves. If the evolution from basaltic andesite to andesite is caused by 60 C cooling (Fig. 15), then the <2 % change in SiO 2 between the melt inclusions and host glass could be driven by <30 C of cooling. A few of the melt inclusions in basaltic andesites have anomalously low MgO concentrations, no doubt owing to small amounts of post-entrapment crystallization (Fig. 4), but the vast majority of them have compositions matching those of the host glass, indicating a near absence of exotic melts that would be indicative of hybridization.
The sulfur concentrations in olivine-hosted melt inclusions in both the basaltic andesites and the andesites are close to those predicted at sulfide saturation (Smythe et al., 2017) and produced experimentally (Moune et al., 2009). Moreover, most Hekla rocks contain pyrrhotite, generally as inclusions in phenocrysts. The measured sulfur concentrations in melt inclusions in the andesites are consistent with fractionation of 0Á10-0Á15 % by mass pyrrhotite from a basaltic andesite melt.

Origin of low 18 O/ 16 O
Because the oxygen isotope ratios of nearby basalts overlap with those of the evolved rocks erupted from Hekla (Fig. 14) (Sigmarsson et al., 1991(Sigmarsson et al., , 1992a. Models of crystallization that include changing isotopic fractionation factors with temperature and mineral and melt composition predict an $0Á2 & increase in melt d 18 O as the SiO 2 concentration increases from 55  Nielsen's (1992) partition coefficients to the mineral assemblage and melt fraction calculated from least-squares massbalance of the major elements. Compositions are normalized to a primitive basaltic andesite. It should be noted that the vertical axis is linear, not logarithmic. Propagated analytical uncertainty is $3 % RSD for most of these elements, thus smaller than the symbols, with the exception of V. (b) Same as above but showing fractionation model for dacite from andesite and normalized to andesite. A noteworthy feature is the positive excursions of U, Th, and Hf, which are attributable to the xenocrystic zirdon in the dacite. (c) Fractionation model for dacite to rhyolite, applying the partition coefficients of Bacon & Druitt (1988). Zircon abundance is estimated by mass-balance of zirconium.
to 70 % (Bindeman et al., 2004). The model predictions are in close agreement with empirical data from suites of magmas that have evolved to dacite and rhyolite by fractional crystallization at the Galá pagos Spreading Center ($0Á4 & difference between basaltic andesite and rhyolite; Muehlenbachs & Byerly, 1982) and the Galá pagos Islands (Geist et al., 1995). Superimposed on the melt trend is an increase in D 18 O(melt-olivine) by 1 &, because of compositional and temperature effects on olivine-melt fractionation: olivine is more Fe-rich and the melt more Si-rich as the magma cools. Therefore, d 18 O in olivine is expected to decrease by about 0Á6 & from basalt to rhyolite, or 0Á2 & as basaltic andesite magma cools from 1010 C to 875 C and the melt evolves from 55 to 70 % SiO 2 , precisely what is observed in the Hekla basaltic andesite to rhyolite suite ( Fig. 14; Bindeman et al., 2008a). The moderate scatter in d 18 O in olivine from the same eruption probably involves inheritance from isotopically heterogeneous basaltic melts, comparable with those in the large-volume fissure system of Iceland (e.g. Bindeman et al. 2008b). In both cases recycled olivine will preserve heterogeneous isotopic compositions of melts even if they have equilibrated with respect to Fe-Mg, owing to slow diffusion of oxygen relative to Fe-Mg diffusion (Ryerson et al., 1989).

Origin of dacite and rhyolite
The origin of the dacites is particularly important, because they have been proposed to be the result of crustal melting or assimilation coupled with crystallization, on the basis of evidence of U/Th isotopic disequilibrium (Sigmarsson et al., 1992a;Chekol et al., 2011). Dacite magma is posited to be the parental magma to Hekla's Plinian rhyolites, which formed by fractional crystallization of the dacite. Our findings, in contrast, indicate that the dacites and rhyolites adhere nearly perfectly to a fractional crystallization model, whereby the dacites derive from crystallization of silicic andesite. The composition of sample H05-63, from the 1158 CE dacite (whose composition is the same as that of 1158 tephra, H09-01), can be created by removal of 2Á3 % olivine, 6Á8 % plagioclase, 1Á1 % augite, 1Á0 % magnetite, 0Á4 % apatite, and 0Á4 % ilmenite from an andesite melt (H05-18) with 62Á5 % SiO 2 (mineral compositions are averages of compositions from the parental andesite; Supplementary Data Electronic Appendices 6-9). This leaves 88 % dacitic melt, with a sum of the squares of the residuals of 0Á86. We note, however, that there is strong petrographic evidence that the petrogenesis of the dacite is more complicated, because it contains xenocrystic crystals, including zircon, that we interpret as having formed from the 1104 CE rhyolite and being incorporated in a dacitic melt. Trace element systematics also supports the crystallization model, although the partition coefficients for Sc, V, and Sr appear to be underestimated (Fig. 16b), and high U, Th, and Hf concentrations no doubt reflect the xenocrystic zircon observed in these rocks. As with the basaltic andesite to andesite model, U and Th are predicted to be richer in the evolved rocks than is observed, which we attribute to underestimation of the partition coefficients of U and Th in apatite.
The major element mass-balance model relating the 1158 CE dacite to the 1104 CE rhyolite (H07-25) is an almost perfect fit, with an r 2 < 0Á01, and previous work by Sigmarsson et al. (1992a) has also proposed that Hekla's rhyolites are generated by fractional crystallization of dacite melt. The model predicts crystallization of 3 % olivine, 10 % plagioclase, 1Á4 % ferroaugite, 1Á3 % magnetite, and 0Á3 % apatite (mineral compositions are averages of analyses of crystals in the parental dacite; Supplementary Data Electronic Appendix 6). The low concentrations of sulfur in the melt inclusions in olivine from the rhyolite (96 6 50 ppm) also require that pyrrhotite, which is observed in the rocks, also fractionates, but we do not include it in the trace element model. We adopt the partition coefficients of Bacon & Druitt (1988) for dacite and interpolate for elements not included in their dataset (Supplementary Data Electronic Appendix 10). The amount of zircon fractionation was calculated by mass-balance of Zr, and partition coefficients calculated from average Hekla zircon from Carley et al. (2011) were applied. Again, the fit is excellent except for the compatible elements (which are particularly sensitive to the uncertain values of the partition coefficients) as well as U and Th. Thus, our proposal for the formation of the 1158 CE dacite requires a combination of fractional crystallization and incorporation of the cumulate complement to the 1104 CE melt-rich magma that erupted 54 years previously.
Pyroxene from the 1158 CE dacite contains reversely, normally, and un-zoned grains, indicating that many of the grains did not simply crystallize from their host melts. Some plagioclase and pyroxene in the 1104 CE rhyolite is broadly normally zoned, consistent with its growth from a cooling melt, whereas olivine and some pyroxene and plagioclase is unzoned (Figs 7-9). Melt inclusions in olivine in the 1104 CE rhyolite are slightly more evolved (by about 4 % SiO 2 ) than the host rock, but the compositions are the same as the host glass (Fig. 4), indicating that these crystals formed from their carrier liquids.

Sparse evidence for an open system
In our experience, the Hekla rocks form one of the most compositionally coherent igneous suites identified to date. With few exceptions, major and trace elements form strikingly consistent trends, and most phenocryst compositions correlate strongly with those of the carrier melts. There is a complete absence of exotic melt inclusions (Fig. 4), and crystals that might have been inherited from mafic and silicic endmembers are absent. If the intermediate rocks are hybrids between basaltic andesite and a silicic magma, nearly all the crystals would have had to have formed after hybridization, and we know of no locality where hybrids do not contain inherited crystals. In the thousands of rocks we examined in the field, we observed only one xenolith locality (silicic hyaloclastites in a prehistoric fissure on the NE shoulder of the volcano) and only a single silicic rock with mafic enclaves.
That said, there are several specific examples of crystals that are distinctly out of equilibrium with their carrier melts, as follows.
1. The uppermost tephra of the 1947 eruption, which is andesite, contains reversely zoned olivine, with rims that are more magnesium rich than the cores, by about 5 % Fo. These reversely zoned olivine have the same composition as the olivine in the first-erupted andesitic tephra (Figs 5d and 7), and they contain evolved melt inclusions that have the composition of the first-erupted rocks (Figs 4 and 5b). The late tephra also contain reversely zoned pyroxene whose cores are the composition of pyroxene in the early erupted andesite, but with rims that are $7 mol% richer in diopside component than the cores (Fig. 8).
Our interpretation is that the crystals initially grew from the evolved andesitic melt, and then were immersed subsequently in a more Mg-rich andesite melt shortly before eruption; there must have been convective stirring of the magma to some degree. Application of Fe-Mg diffusion coefficients of Chakraborty (2010) at 1000 C yields residence timescales of 70 days for the olivine to account for the thickness of the boundary layer (assuming a wellmixed open boundary). The late-erupted, lessevolved tephra also has anomalous trace element compositions, with Nb and Zr concentrations that are $6 and 15 ppm lower than expected from their Rb concentrations (Fig. 11). This offset can be explained by mixing with a small fraction (<10 %) of cooler melt that had fractionated ilmenite and zircon, a composition similar to the 1104 CE rhyolite, or assimilation of a rock of similar composition. 2. The 1158 CE dacitic lava contains both normally and reversely zoned ferroaugite, whose cores are about 6 mol% richer and poorer in diopside component than the rims. This suggests that some pyroxene cores grew from slightly less evolved dacite and rhyolite melts, followed by growth from the carrier dacitic melt. The boundary layers separating these zones are 20-30 mm thick. As noted above, some of the phenocrysts in both the dacite lava and tephra contain zircon, even though the melt was not zirconsaturated (the zircon saturation temperature for this rock is calculated to be 619 C, >200 C below its eruption temperature). Our interpretation is that the 1158 CE dacite contains phenocrysts that originally crystallized from the 1104 CE rhyolite, and the 1158 magma is a hybrid between dacitic melt and crystal mush residual from the 1104 rhyolite. This makes sense given the close temporal relationship between the two eruptions.
3. A single plagioclase from the 2000 tephra has a core that is An 47 , then is reversely zoned with a wide rim close to equilibrium with the carrier melt (An 57-60 ). This plagioclase was probably inherited from a dacite (Chekol et al., 2011), but it resided in the basaltic andesite magma long enough for a 500 mm rim to grow. This plagioclase crystal is only one of dozens we analyzed from this eruption, the others having close to equilibrium compositions. 4. The only mafic enclaves we observed and sampled in the historical rocks of Hekla are in the 1766 CE andesitic lava, which contains the most forsteritic olivines of any sample, up to Fo 86 . In all likelihood, these olivines crystallized from a basaltic melt that has not been sampled otherwise at Hekla. As other researchers have pointed out (Sigmarsson et al., 1992a, and others), the likely parent to Hekla basaltic andesite is a composition similar to that of the nearby 1913 Lambafitarhraun basaltic fissure eruption. 5. A pyroxene in the 1510 CE andesite contains magnetite, ilmenite, and zircon inclusions and clearly formed from a dacitic magma. This is the only example of such a crystal in the rocks we inspected. Its origin may be from mixing with silicic magma that had survived for 352 years or partial assimilation of solidified 1158 CE dacite.
Subtle differences in the trace element concentrations of the 1947 CE eruptive products suggest that its magma might have assimilated a small amount of silicic rocks or mixed with silicic magma, comparable in composition with the 1104 CE rhyolite. Specifically, the intermediate 1947 rocks are characterized by anomalously low Nb/Rb and Zr/Rb, which could be accounted for by addition of <5 % rhyolitic melt (or granite) with the composition of the 1104 CE tephra. On the basis of U-Th isotopic disequilibrium and a correlation with SiO 2 , Chekol et al. (2011) also found that the 1947 intermediate magmas also are affected by assimilation of young rocks, although they invoked interaction with rocks that we attribute to Torfajokull.
In summary, phenocrystic and melt-inclusion evidence for magma mixing, assimilation, or reentrainment of cumulate crystals in Hekla's silicic rocks is sparse. Most crystals in Hekla magmas formed from the melt that hosts them at the time of eruption. No cumulate xenoliths and exceedingly few 'antecrysts' have been sampled from the volcano, and no material that might be attributed to a cumulate mush has been observed. Consequently, crystal-liquid segregation must be very efficient at all stages of magmatic evolution at Hekla.

The case against crustal melting
Previous studies have invoked several pieces of evidence to argue that Hekla dacites and rhyolites result from crustal melting or extensive crustal assimilation and are related to the basaltic andesites and andesites only via the transfer of heat.

Hekla's silicic rocks have low 18 O/ 16 O compared to
the Earth's mantle (Sigmarsson et al., 1991;Bindeman et al., 2012). The origin of the low d 18 O in olivine is discussed above: all evidence points to it being inherited from low d 18 O basalts, followed by isotopic fractionation during crystallization. 2. Dacites erupted in 2000 have distinct Sr, Nd, and Th isotopic ratios (Chekol et al., 2011). 3. Rocks from the large-volume silicic eruptions have fractionated Th/U relative to the mafic rocks (Sigmarsson et al., 1992a). U-Th isotopic disequilibrium indicates that this fractionation must have occurred tens of thousands of years ago, In addition to hypotheses that propose that Hekla's silicic rocks are produced by anatexis, Lucic et al. (2016) proposed that the basaltic andesites are due to melting of gabbro. They based their hypothesis on the very low abundance of CO 2 and variable H 2 O concentrations in melt inclusions in Hekla rocks.
The exotic Sr and Nd isotopic compositions of dacites from the 2000 CE eruption noted by Chekol et al. (2011) prove that these rocks could not have originated by fractional crystallization of the co-eruptive basaltic andesites. We question, however, the ancestry of these dacites and suggest that they might be very local and low-volume, remelted pieces of silicic rocks from around the vent in the shallow subsurface, as also found in the 1970 deposits (Sigmarsson et al., 1992a). Thus, their origin bears little on the petrogenesis of the voluminous 1158 CE dacite or the Plinian rhyolites. None of the five lava and tephra samples from the 2000 eruption that we analyzed and no other published analyses (Hö skuldsson et al., 2007;Moune et al., 2007;Portnyagin et al., 2012) show any sign of a silicic component. The dacites reported by Chekol et al. (2011) have the isotopic composition of rhyolites from Torfajö kull (Martin & Sigmarsson, 2010), and young Torfajö kull tephra are exposed all around Hekla. All other analyses of 2000 rocks indicate that the eruption was strongly homogeneous basaltic andesite (Fig. 3), as is predicted by the short repose between 1991 and 2000. Although this contamination is local and not of great petrogenetic significance, the compositions of the rare contaminated rocks in the 2000 CE deposits indicate that crustal contributions are readily detectable at Hekla. Sigmarsson et al. (1992a) made the strong case that recent fractionation of U/Th between the andesite and dacite requires that the dacite be derived by anatexis of amphibolite from the Icelandic crust. The covariation of U and Th of the basaltic andesite to andesite sequence is extraordinarily linear and nearly the same as the andesite to dacite to rhyolite sequence, especially considering that the rhyolites have crystallized a small amount of zircon, and the dacites have accumulated it. Previous investigators have agreed that the transition from basaltic andesite to andesites is caused by fractional crystallization (e.g. Baldridge et al., 1973;Lucic et al., 2016) or in combination with lesser extents of assimilation (AFC; Sigmarsson et al., 1992a;Chekol et al., 2011): the interpretation that is disagreed upon is the origin of the dacites and rhyolites. The linear trend of U and Th in our suite could be explained by extraction of a cumulate with low but constant U and Th concentrations. It should be noted that linear regression of U versus Th results in a line that intersects the U axis (Fig. 12), meaning that if the rocks are related by fractional crystallization, U is slightly less incompatible than Th. Because of the non-zero intercept, the U/Th changes as both U and Th abundances increase with differentiation, a phenomenon observed and emphasized by Sigmarsson et al. (1992a). Our LA-ICP-MS data on apatite are consistent with this hypothesis, although we recognize that the uncertainty is large (Fig. 12). Chekol et al. (2011) noted that literature data indicate that D U apatite > D Th apatite , which well explains the Hekla U-Th trend, although experimental data indicate that the relative partitioning is not well known [Prowatke & Klemme (2006) found values of D U apatite ranging from 0Á2 to 1Á4 and D Th apatite ranging from 0Á7 to 1Á3 between apatite and basaltic andesite and andesite melt].
The hypothesis that the dacites are the result of anatexis (Sigmarsson et al., 1992a) is untenable for other reasons. The compositions of Hekla dacites are vastly different from those of melts produced by melting of natural amphibolite at pressures of between 0Á3 and 1Á0 GPa, which are typically tonalitic and peraluminous. The experimental melts have Al 2 O 3 ranging from 14 to 19 % in dacitic compositions and K 2 O concentrations ranging from 0Á9 to 2Á8 % (e.g. Beard & Lofgren, 1991;Rushmer, 1991;Zhang et al., 2013).
The experiments also indicate that amphibole is in the residue up to 20-40 % melting, at temperatures >1000 C, and trace element concentrations also disprove an origin of the dacites by melting of amphibolite. For example, Brophy (2008) has shown that partial melting of amphibolite yields distinctive REE patterns. Specifically, partial melting of amphibolite yields La concentrations that are very slightly enriched over the concentrations in the source and Yb concentrations that are depleted. In contrast, Hekla dacites and rhyolites have La concentrations that are up to 150 % enriched relative to the most mafic rocks at Hekla and even more so compared to regional basalts (Figs 3 and 16). Yb is also higher in the silicic rocks than the mafic ones. Thus, the REE are inconsistent with melting an amphibolite source to generate the silicic magmas at Hekla.
The basement beneath Hekla is unlikely to have the same trace element and radiogenic isotopic composition as Hekla rocks. The oldest crust around Hekla formed at the Western Volcanic Zone, $60 km west of the volcano. At a half-spreading rate of 18 km Ma -1 (Einarsson, 2008), the oldest rocks beneath Hekla should have ages >3 Ma. Many cycles of volcanism over that time span will result in an isotopically heterogenous crust [e.g. Hemond et al. (1993) found a range in 87 Sr/ 86 Sr from 0Á7029 to 0Á7032 in the Western Volcanic Zone]. In contrast, the Sr and Nd isotopes of Hekla's dacites and rhyolites are the same as those of the mafic rocks (Sigmarsson et al., 1992a), and it would be fortuitous if the rocks that melted to form the dacitic magma had the precise Sr-isotopic composition of Hekla basaltic andesite. Moreover, the near-perfect correlation between incompatible elements across the Hekla suite (Fig. 3) further argues against older crustal sources for the silicic Hekla material. One possibility is that crustal sources of the silicic magmas are older Hekla magmas, thus having isotopic and elemental compositions similar to the material erupted historically (Sigmarsson et al., 1992a). Most potential protoliths in the buried crustal section at Hekla will have undergone some combination of weathering (Patino et al., 2003), hydrothermal alteration, and zeolite-, greenschist-and amphibolite-facies metamorphism, each of which would have been accompanied first by the introduction of fluids at temperatures up to 400 C, followed by production and loss of fluids by dehydration reactions beginning at the greenschist facies (e.g. Oskarsson et al., 1982). For example, a series of Icelandic rocks that have undergone simple burial metamorphism (to 90 C) in the upper few kilometers of the crust (Neuhoff et al., 1999) have undergone large amounts of mineralogical and major element transformation, and mobile trace elements have been altered by more than 10-fold: their composition is unlike that of the unmetamorphosed basalt. Amphibolite-facies metamorphism would then alter the compositions of the original basalts and gabbros even more. If amphibolites melted, they would have left a residue comprising amphibole, pyroxene, and plagioclase up to 20-40 % melting and 1000 C (Rushmer, 1991;Zhang et al., 2013;Jagoutz & Klein, 2018). It is highly unlikely that such metamorphic and anatectic processes could have occurred without strongly fractionating Pb/Rb, Zr/Rb, K/Rb, Th/Rb, U/Rb, Dy*, and Ba/Rb (Figs 14 and 15). It is challenging to speculate on what the sulfur concentration of a basaltic andesite melt produced by melting of gabbro would be (Lucic et al., 2016), and likewise for melting of amphibolite to produce the dacite (Sigmarsson et al., 1992). If the crustal rocks contained residual sulfide, then the observed concentrations would be produced, but we know of no way of testing that hypothesis.
Hekla has been investigated for the stable isotopic composition of Li, Si, Ti, V, Fe, Mo, Zr, and Tl, owing to the large compositional range of the eruptive rocks (Schuessler et al., 2009;Savage et al., 2011;Prytulak et al., 2017aPrytulak et al., , 2017bYang et al., 2015;Deng et al., 2019;Inglis et al., 2019). The results fall into two broad categories: nuclides that show isotopic fractionation (Si, Ti, V, Zr, and Fe) and those that do not (Li, Mo, and Tl). A crustal melting model for the origin of the dacites would require that hydration and metamorphism of the mafic crust be without any stable isotope fractionation, then partial melting of this hydrated metamorphic rock without isotopic fractionation of lithium, molybdenum, or thallium. Moreover, thallium and lithium concentrations correlate with other incompatible elements in the Hekla suite (Schuessler et al., 2009;Prytulak et al., 2017a), which would be exceedingly unlikely if hydrated crustal rocks partially melted. In contrast, iron, titanium, zirconium, and silicon isotopes change in the basalt to rhyolite sequence by an amount predicted by fractional crystallization models (Schuessler et al., 2009;Savage et al., 2011;Deng et al., 2019;Inglis et al., 2019). Vanadium isotopes also correspond to crystallization models for the basaltic andesite to andesite sequence (Prytulak et al., 2017b), but unfortunately the silicic rocks are too poor in V to measure the isotopic composition. A crustal melting model for the origin of the dacites would require that hydration and metamorphism of the mafic crust, followed by partial melting of this hydrated metamorphic rock, should fractionate the stable isotopes of Fe, Ti, Zr, and Si in exactly the same magnitude as predicted for crystallization.
The case against crustal melting is further supported through examination of a nearby center whose silicic rocks have been shown convincingly to be due to anatexis. In contrast to the tight, coherent trend of incompatible elements at Hekla, compositions of Torfajö kull rocks, especially the silicic ones, are strongly scattered ( Fig. 17; Gunnarsson et al., 1998;Zellmer et al., 2008), particularly when variations in relatively immobile incompatible elements such as Zr are compared with those of mobile ones, such as Rb. Moreover, amphibole-bearing xenoliths have been described at Torfajö kull (Gunnarsson et al., 1998) but have not been discovered at Hekla. Although the two volcanoes are close, they lie in different tectonic settings, and Torfajö kull has a more prolonged hydrothermal and volcanic history. Nevertheless, it provides an excellent example of how crustal melting and extensive assimilation can be identified via petrological and geochemical tools.

Constraints from volatiles
The correlation of the maximum H 2 O concentrations in melt inclusions and the plagioclase hygrometer with K 2 O (Fig. 10) confirms that H 2 O is an incompatible element in the Hekla system (Moune et al., 2007;Portnyagin et al., 2012;Lucic et al., 2016). The calculated saturation pressures indicate the minimum pressure at which the inclusion was trapped by the host grain (Figs 10 and  15). The inclusions with lower H 2 O contents are likely to have lost hydrogen via diffusion through the host as the surrounding magma degasses during ascent. Although the diffusivity of H þ is poorly known for fayalitic olivine, existing data (Bucholz et al., 2013;Barth et al., 2019) suggest that the H 2 O contents in Hekla melt inclusions could be reduced by $50 % during ascent in $2 days (for 1020 C basaltic andesite) to $12 days (for 870 C rhyolite; assuming a 40 mm inclusion that is 460 mm from the grain boundary along the a-axis of olivine; A. Barth, written communication, 2020). In fact, most of the melt inclusions in this study are <100 mm from the grain boundaries, thus the diffusive timescales could be far shorter (Supplementary Data Electronic Appendix 5). The very narrow range of plagioclase compositions in rocks that have a wide range of melt inclusion H 2 O concentrations is evidence against the measured range of the melt inclusions being due to different magmatic H 2 O concentrations ( Fig. 10; Waters & Lange, 2015). The melt inclusions measured by Lucic et al. (2016) appear to be more strongly affected by degassing than those in other melt inclusion studies and our sample set, although this is not their interpretation.
The concentration of H 2 O in melt inclusions in the silicic rocks and that calculated with the plagioclase hygrometer agree with H 2 O concentrations predicted by the crystallization model but could not have been generated by the melting of crustal materials (Portnyagin et al., 2012). A typical amphibolite with a basaltic or gabbroic protolith has $50 % amphibole and thus contains $1Á0 % H 2 O. It is theoretically possible to generate a dacite with 5 % H 2 O by melting the protolith by 20 % if H 2 O is perfectly incompatible, but because amphibole is in the residue, melting would more probably require melting $5 % to create a melt with 5 % H 2 O (Jagoutz & Klein, 2018). Melting of either a Hekla-related rock or one from the Western Volcanic Zone <5 % would create a melt that is enriched in incompatible trace elements many times that of Hekla dacite, and trace element ratios indicative of residual amphibole would be strongly fractionated. Similar arguments can be made against the hypothesis of melting of gabbro to create basaltic andesite magma (Lucic et al., 2016): gabbro contains too little H 2 O to produce basaltic andesite magma with 2-3 % H 2 O, especially if it contains only the anhydrous minerals known to have crystallized from Hekla magmas and regional basalts.
The very low CO 2 concentrations in the melt inclusion glasses are perplexing: if one assumes a constant CO 2 /Nb of 400 in Icelandic magmas [within the range estimated by Maclennan (2017)], then one would predict CO 2 concentrations of 2Á6 % in a typical Hekla basaltic andesite, far greater than the solubility of CO 2 at the conditions of magma storage and olivine growth at Hekla (calculated to be $2500 ppm at 0Á45 GPa; Shishkina et al., 2010). In other words, at any depth in the crust and well into the upper mantle, Hekla's magmas are probably saturated with a CO 2 -rich fluid. One way to explain the low CO 2 concentrations is that the basalts that are parental to Hekla basaltic andesite are created by remelting of strongly hydrated rocks in the lower crust, as hypothesized by Bindeman et al. (2008b). This hypothesis is supported by the oxygen isotope data in both the basalts and the Hekla rocks (Fig. 14), which demonstrate that individual olivine grains crystallized from melts with variable oxygen isotope compositions. For the Fa-rich olivine in the rhyolite, melt inclusions could have trapped melt that had evolved by large extents of vapor-saturated fractional crystallization from less differentiated magma at >200 MPa (Wallace, 2005), and this would have produced rhyolitic magma with little dissolved CO 2 . The differentiation to form rhyolitic melt could also have occurred at higher pressures than this, as long as the melts ascended to $200 MPa and crystallized their phenocrysts at that depth, consistent with experimental phase equilibrium evidence. The lack of CO 2 in the more mafic melt inclusions could potentially be the result of post-entrapment bubble formation (e.g. Wallace et al., 2015), as bubbles are common in the melt inclusions we analyzed, but we note that Lucic et al. (2016) analyzed only melt inclusions without vapor bubbles and found similarly low to below detection values. A final possibility is decrepitation of the melt inclusions during ascent (e.g. Maclennan, 2017). For example, it is estimated that fewer than 10 % of the melt inclusions in olivine at Laki record magmatic CO 2 concentrations. Further study of Hekla melt inclusions by microanalysis and Raman measurements of vapor compositions in included bubbles are called for to solve this CO 2 paradox.
The low CO 2 concentrations might also be explained by circulation of magma to shallower depths, degassing, and return to the lower crust, as has been proposed at Kilauea (Wallace & Anderson, 1998;Barsanti et al., 2009), although this is unlikely to occur in the more evolved rocks, which are buoyant relative to the basaltic andesite irrespective of volatile content. de/georoc/). We attribute the scatter in the Torfajö kull data to contributions of crustal rocks that have been altered by hydrothermal activity and metamorphism, then fractionated by partial melting in the presence of amphibole and other metamorphic phases.
Thorium isotopic disequilibrium and origin of the 1104 CE rhyolite A feature of the Hekla magmatic suite is yet to be explained and is inconsistent with a crystallization origin for the silicic rocks: model ages from Th isotope analyses of the dacites and rhyolites, and their constituent zircon, are much greater than the repose times between eruptions (Thorarinsson, 1950(Thorarinsson, , 1967our Fig. 2). Model ages based on whole-rock Th isotopic analyses of the dacites and rhyolites are as high as 76 000 6 14 000 years (Sigmarsson et al., 1992;Chekol et al., 2011). Model ages of zircon that crystallized from a rhyolite magma range from modern to 43 000 6 7000 years, and they have heterogeneous oxygen isotope compositions (Carley et al., 2011;Bindeman et al., 2012). These model ages are in stark contrast to the observation that the rhyolite erupted no more than 480-720 years after the previous explosive eruption from Hekla, which is zoned from dacite to basaltic andesite (Larsen et al., 2020). In addition to the incontrovertible evidence of zircon recycling, the relatively low ( 238 U/ 232 Th) in the rhyolites almost certainly involves removal of zircon, but the old model ages of both wholerocks and individual zircon crystals require that this occurred tens of thousands of years before the eruptions. The effect of the incorporation of zircon into the 1158 CE dacite and the control on the model ages are unknown and would require detailed single-crystal study of those rocks. Simple mass-balance (Fig. 12) suggests that only a small fraction of the U and Th in the dacite derives from xenocrystic zircon.
As noted by Sigmarsson et al. (1992a), there is little difference between thorium isotopic compositions of the basaltic andesites [( 230 Th/ 232 Th) ¼ 1Á014 6 0Á020] and the andesites [( 230 Th/ 232 Th) ¼ 1Á003 6 0Á020], meaning that the isotopes are consistent with the andesites originating by fractional crystallization of basaltic andesite over the course of a century. The U-Th isotopic disequilibria of the dacites and rhyolites and the model ages derived from them preclude that they could be generated by recent (<10 000 years) fractionation of andesite melts. Given that the strong correlations of trace elements, including U and Th (Fig. 12), preclude an origin of the dacites and rhyolites by extensive assimilation or anatexis, we propose that the 1104 CE rhyolite was produced in a separate, crustal magma body by extensive crystal fractionation of the more commonly erupted mafic magmas from Hekla. One observation that may support the hypothesis of a separate magma body is that there is no known late-erupted basaltic andesite lava associated with the 1104 CE and 1158 CE eruptions, although we recognize that it is possible that they could be buried by younger rocks.
We propose that once produced by extensive fractionational crystallization, either the silicic magma or mostly solidified granodiorite to granite is stored over the course of >10 4 years, conditions that are respectively referred to as 'hot' and 'cold' storage. The cold storage hypothesis is motivated by the commonly observed discrepancy between 238 U-230 Th model ages and other geochronometers (Cooper & Kent, 2014), precisely the observation at Hekla, although in this case the discrepancy is between U-Th isotopic model ages of rocks (and zircon) and eruption interval. Although we recognize that the U-Th isotopic disequilibrium can be explained in other ways, we contend that the Hekla rhyolite must be stored in a largely molten state for at least tens of thousands of years, under hot, largely molten conditions. Beyond zircon with old model ages, there is no evidence for remobilized mush (or rock) in the 1104 CE rhyolite. Most systems that we are aware of that have been attributed to cold storage (e.g. Mt Hood, Okataina, Kneeling Nun) are characterized by rocks with an abundance of phenocrysts that bear evidence of being reheated and mobilized, including reverse zoning and multiple crystal populations. Moreover, the relatively H 2 O-rich compositions of Hekla's silicic rocks are nearly impossible to have been derived by remelting of solid rock. Storage of molten magma for long periods is facilitated by the high temperatures (600-800 C) of the Iceland lower crust, as well as recharge by the steadystate basaltic andesite magma Fig. 18).

Difference between tephra and lava
The systematic difference in compositions between tephra and lava in all eruptions (Fig. 11) remains enigmatic, and we do not have a satisfactory explanation for it. We believe that it is not an analytical artifact, because we analyzed the tephra and lava samples in random  Turcotte & Schubert (2002) for conductive cooling of a semi-infinite half-space is applied, with the modification that a term is added to the heat capacity to account for the latent heat of crystallization. Model curves are for different dike half-thicknesses and far-field wall-rock temperatures ranging from 400 to 800 C. Horizontal lines indicate temperature estimates for basaltic andesite and andesite. Vertical line indicates the typical hiatus observed at Hekla preceding eruptions that first erupt andesite; the intersection of this line with the cooling curves indicates that andesite can be made from a cooling basaltic andesite dike that is 100 m wide in 100 years, if the background geothermal gradient is 500-600 C over the height of the magma body. order and by three different methods (XRF, ICP-MS, and electron microprobe). Contamination owing to different grinding protocols is unlikely, as the tephra is systematically richer in Sr, Y, V, and Sc, but it never came in contact with metal during sample preparation. Analyses were completed on scoria and pumice clasts 0Á5 to >2 cm in diameter, so this is not a manifestation of crystal-glass fractionation in the eruption plume, although it is possible that the rinds of clasts, which we preferentially sampled, have fewer crystals than the interiors. The primary compositional differences are that the tephra is richer in both compatible (Sc, V, and Sr) and incompatible (Rb, Nb, La) elements than the lavas that erupted shortly after the tephra. In terms of major elements, the major differences are that tephra has less TiO 2 , MgO, and P 2 O 5 . One possibility is that part of the magma body assimilated a rock with these characteristics, but we are not familiar with rocks that are rich in TiO 2 , MgO, P 2 O 5 , Sr, Sc, V, and Rb. Alternatively, the earlier-erupted rocks may have crystallized and fractionated olivine, magnetite, and apatite (but not feldspar or clinopyroxene), or the later-erupted (lava) magma accumulated these crystals. Least-squares mass-balance on the most extreme lava-tephra pair from the 1991 eruption indicates that the major element differences are explained by 2Á8 % fractionation of olivine and apatite (about 2:1 proportions), thus whatever process was responsible for the systematic differences between the lava and tephra it had only a minor effect. The greater abundance of Sc and Rb in the tephra, however, requires up to 27 % accumulation of olivine and apatite in the lava. The sample used for the mass-balance calculations (H05-02) has only 4 % crystals dominated by plagioclase and clinopyroxene. The most credible explanation is that the tephra derived from a magma body that crystallized the assemblage plagioclase þ clinopyroxene þ olivine þ magnetite þ apatite 6 ilmenite, and the olivine and apatite (6 magnetite) were removed from this magma, perhaps by density sorting.

Hekla's magmatic plumbing
A number of distinctive aspects of Hekla require consideration in any hypothesis that addresses its magmatic evolution. First, a multitude of geophysical and petrological evidence indicates that inter-eruptive cooling and crystallization over the past 50 years occurs in the middle and lower crust, centered at about 20 km depth (pressures $0Á6 GPa). Second, deformation and seismicity related to an eruption progresses extraordinarily rapidly, in the matter of 30-60 min (e.g. Soosalu & Einarsson, 2002;Soosalu et al., 2003;Soosalu & Einarsson, 2004;Geirsson et al., 2012). Third, chemical evidence is most consistent with fractional crystallization, and the eruptive rocks are crystal-poor (Fig. 6), indicating that crystal-liquid differentiation is efficient. Fourth, the chemical processing of the magma occurs on a timescale of decades to centuries.
Because most eruptions return to a basaltic andesite composition during the waning stages of the eruptions, and basaltic andesite is the most voluminous eruptive product, there must be a large supply of this $1020 6 20 C magma in the lower crust at Hekla. The basaltic andesites are almost certainly formed by cooling and crystallization of basaltic magma (Sigmarsson et al., 1992a). The constancy of the baseline composition and deformation data suggests that basalt intrudes into the base of the crustal magmatic system. The heat input from the basaltic recharge magma must balance (at least over the course of recorded history) that lost to the surroundings to keep the magma body in a thermochemical steady state, buffered at about 1020 6 20 C, and the melt composition evolving to basaltic andesite. We emphasize that the Hekla plumbing system is special, in that it leads to extremely efficient crystal-liquid segregation: there is no evidence in the erupted basaltic andesites for the existence of the basaltic parents or the cumulate complement, including crystals with primitive compositions, mafic melt inclusions, and mafic enclaves (one exception in an andesite lava).
Previous studies (e.g. Thorarinsson & Sigvaldason, 1972;Sigmarsson et al., 1992a;Gudmundsson et al., 1992) have proposed a distinctive geometry for the Hekla magma reservoir, a tall inverted-funnel-shaped cupola that extends $7 km vertically and $5 km wide, located above a lower-crustal sill. At least one geodetic study (Geirsson et al., 2012) suggested a different geometry: an extended vertical reservoir, located from 21 to 10 km depth. We hypothesize that a vertically elongate tabular shape is responsible for the unique characteristics of the volcano: time-dependent, rapid ($1 century) zonation of the magma between eruptions (Fig. 19). In turn, the tall, thin geometry of the magma body is a manifestation of Hekla's tectonic location at the inside corner of the intersection of the Eastern Volcanic Zone and the South Iceland Fault Zone. An elongate, upright tabular magma body is also consistent with the seismic observations (Soosalu & Einarsson 2002;Soosalu et al., 2003;Soosalu & Einarsson, 2004). If molten magma persists between eruptions from a depth of 20 km to the shallow crust, pre-eruptive seismicity would be minimized, because the magma is stored and ascends through a long-established, trans-crustal molten conduit.
Heat flow measurements indicate that the background geothermal gradient in the shallow crust near Hekla is $100 C km -1 (Arnó rsson et al., 2008). The gradient must be significantly lower in the deep crust, as it approaches the adiabat of the shallow mantle beneath Iceland. For example, the steady-state geotherm in the lower crust of 5 Ma oceanic lithosphere is $45 C km -1 (e.g. Turcotte & Schubert, 2002). Thus, the temperature difference between basaltic andesite melt and the farfield wall-rock is likely to be between 300 and 700 C greater at the top of Hekla's magma reservoir than it is at the bottom (recognizing that the geothermal gradient around a volcano is exceedingly complex owing to transfer of heat by precursory magma and active hydrothermal circulation).
We propose that more cooling and crystallization takes place in the upper part of the tall magma storage reservoir than in the lower section, owing to different far-field temperatures with depth, resulting in variable extents of differentiation. To test this hypothesis, we apply a simple thermal model (Fig. 18), based on equation (4-113) of Turcotte & Schubert (2002), which describes conductive cooling of a semi-infinite halfspace. To account for latent heat of crystallization, we add an energy term: the latent heat divided by the amount of crystallization (0Á3) over a specified temperature range (60 C). This simplified model assumes the magma is stagnant and develops a thermal boundary layer. If convection in the magma body was vigorous, the Nusselt number would be two (Marsh, 1989), and the resulting heat loss would be at twice the rate calculated below.
Our thermal model parameters are based on the 1947 CE eruption. Geothermometry indicates that the most silicic andesites of the 1947 eruption, which erupted explosively and first, were $50 C cooler than the later-erupted and parental, effusive basaltic andesites. Preceding this eruption was the well-documented 1845 CE eruption, which was similarly zoned (up to 62 % SiO 2 ) after a 79 year hiatus. Because the 1845 eruption ended with the eruption of basaltic andesite, we assume that the conduit extending from the surface to the lower-crustal magma reservoir was completely filled with $1010 C basaltic andesite magma at the end of Fig. 19. Conceptual cross-section of the Hekla magmatic plumbing system, drawn perpendicular to the NE-SW strike of the volcano's fissure system and using the same scale vertically and horizontally. Left panel shows the thermal and compositional evolution after 20 years of repose. Right panel shows the thermal and compositional evolution after a century of repose. The geothermal gradient is such that the far-field temperature is 400 C near the top of the system and 800 C near the base. The thin body extending from the lower-crustal magma reservoir is only 100 m wide, and the inset shows it expanded 20Â in the horizontal. This body barely cools in 20 years, but over a century its upper segment cools 65 C, and the melt evolves from basaltic andesite (red) to andesite (purple shades). The rhyolite from the 1104 CE and dacite from the 1158 CE eruptions derived by fractional crystallization of basaltic andesite, but in an independent magma body that is molten for >10 000 years. Its location in the crust is poorly constrained but >6 km depth. that event. Geochemical mass-balance indicates that 50 of cooling corresponds to $30 % crystallization. For the thermal model (Fig. 18), we apply standard values of thermal conductivity of 1Á4 J m -1 s -1 K -1 , heat capacity of 1 J K -1 , a density of 2440 kg m -3 , and heat of crystallization of 400 J g -1 m -1 .
This thermal model suggests that the interior of a vertical tabular body slightly wider than 100 m (50 m half-width) would cool to the temperature of the erupted andesite with 600 C wall rocks after $100 years, whereas the magma in contact with 800 C wall rock would have cooled only about 25 C during that interval (Fig. 18). If the igneous body flared at depth to a half-width >100 m, as envisioned by Sigmarsson et al. (1992a), cooling over the same period would be insignificant (Fig. 18). The volume of the 1947 eruption has recently been recalculated by Pedersen et al. (2018) to be 711 6 122 million m 3 dense rock equivalent. A tabular magma body that is 100 m wide and 7 km tall would have to be only 1 km long along-strike to accommodate this volume, which is realistic given that the vents of the 1947 eruption are >10 km apart along-strike ( Fig. 1). This thickness calculation does not account for the volume of crystals that forms as the result of cooling, which is $30 % for the differentiation of basaltic andesite to silicic andesite. The strike length therefore would have to be 1300 m, again a small proportion of the distance between the vents. If the tall magma body were more cylindrical, as proposed by Geirsson et al. (2012), the cooling times would be slightly longer.
The 1947 CE eruption is the best-documented eruption that followed a long hiatus and consequently is strongly zoned. In contrast, the more recent eruptions in the 20th century followed much shorter repose, between 9 and 23 years. Application of the dimensions of a tall and thin magma body cited above (1000 m Â 100 m Â 7000 m) to the cooling model indicates that heat loss over 10 years would be about 1 C (Fig. 18), which would lead to inconsequential extents of fractionation, as has been observed. Observations of the repetitive eruption sequences, each eruption returning to a baseline composition of basaltic andesite, suggest that recharge after each eruption resets the magma in this vertical body to $1010 C.
The increasing buoyancy of the melt owing to crystallization is a critical factor that controls the zonation of the tall and thin Hekla magma reservoir, because the more evolved melts are progressively less dense, even though they are cooler, owing to crystallization of irontitanium oxides and removal of iron from the melt. Basaltic andesite melt at 1010 C has a density of 2450 kg m -3 , whereas the evolved andesite that erupted early in 1947 CE has a density of 2350 kg m -3 at its eruption temperature of 950 C (Lange & Carmichael, 1987). Thus, if the evolved melts are made by greater extents of cooling in the upper, cooler part of the tall magma body, the reservoir will be stably stratified (Fig. 19). Consequently, there will be no driving force for compositional convection, preserving the heterogeneity from top to bottom of the elongate reservoir.
Differentiation of the magma by cooling in a narrow body that extends for >10 km vertically also might explain the variable H 2 O concentrations and very low CO 2 concentrations. Replenishment into this body is always of the steady-state basaltic andesite composition, and if some of this magma were stored at depths of only a few kilometers, the melt would be poor in both CO 2 and H 2 O (Fig. 19). Basaltic andesite magma that is progressively richer in volatiles would extend down to the steady-state reservoir in the lower crust. The magma stored in the vertical magma body supplies the magma that is variably degassed and erupts in the early explosive phase of Hekla eruptions. We note that all of our melt inclusion analyses are from the early explosive phase, a bias that is introduced because it is well known that tephra preserves melt inclusion H 2 O concentrations better than lava. One way to test this hypothesis is that it predicts that melt inclusions in the more voluminous later-erupted lavas should have higher CO 2 concentrations. Evidence against this hypothesis is provided by the mineral compositions: early erupted materials have pyroxene and plagioclase that are no different in composition from the late-erupted minerals, and crystallization at systematically different pressures and in the presence of different volatile concentrations would produce heterogenous populations.
The greatest conundrum in the evolution of Hekla's magmas is the extraordinary efficiency of the crystalmelt segregation and the absence of evidence for a cumulate mush. Crystal-rich clots, phenocrysts out of equilibrium with the carrier melt, exotic melt inclusion compositions, and cumulate xenoliths are exceedingly rare. We speculate that Hekla's suite of compositions might be due to the unusual magma chamber geometry. As pointed out by a number of researchers (e.g. Marsh, 1989;Maclennan et al., 2001;Gudmundsson, 2012;Geist et al., 2014), most magma reservoirs emplaced in the oceanic crust are sill-like, especially those at shallow levels, and that geometry may be more conducive to remobilization of crystal-rich mush. Some have argued that the advance of a solidification front is more rapid than the gravitational separation of residual liquids in a sill (Marsh, 2006(Marsh, , 2013, although this principle is controversial (e.g. Latypov et al., 2015;Bachmann & Huber, 2018). Owing to the dominance of horizontally extensive bodies (i.e. sill-like reservoirs), most modern work on crystal-liquid segregation considers mushy and crystal-settling systems dominated by gravitational forces (e.g. Bachmann & Huber, 2018). In contrast, there has been little research into crystallization on the sides of tall magma bodies since the study by McBirney et al. (1985). In fact, we envision the differentiation in the Hekla magma body to be much like that invoked by McBirney et al. (1985), although we attribute the accumulation of silicic melt at the top of the system to result from crystallization in the presence of a geothermal gradient as opposed to buoyant boundary layers of melt rising along the sidewall.

CONCLUSIONS
On the basis of our detailed analysis of Hekla materials, the preponderance of evidence is that the sequence of magmas ranging from 54 to 70 % SiO 2 at Hekla volcano results from fractional crystallization, with at most minor contributions from crustal melting or assimilation of crustal materials. On the basis of evidence from zoned eruptions, most of the range of magmatic compositions must exist contemporaneously in the plumbing system at any given time. Intra-eruption mixing is rare, suggesting that the eruptive sequence taps a stably stratified plumbing system. The systematic relationship between repose time and the extent of evolution is remarkable and is explained by the marginal cooling of a tall, thin magma body that is surrounded by wall rock that is hundreds of degrees cooler at its top than its base.
The strong correlations between trace elements of different mobilities during hydrothermal and metamorphic processes in the Hekla suite, along with the incompatible behavior of H 2 O, argue strongly against significant crustal interaction in the basaltic andesite to rhyolite compositional spectrum. The basaltic andesite to andesite sequence develops over the course of a century in the tabular conduit that extends from the lower to middle crust, perhaps more than 10 km high. Actinide-series disequilibria suggest that the rhyolite developed in a separate magma body, one that also evolved by crystallization, but for a protracted time (>50 000 years). The depth at which the silicic magma is stored is poorly constrained, beyond that H 2 O-rich vapor is saturated in the melt at 7 km depth, providing a minimum estimate. The low CO 2 concentrations in melt inclusions suggest, however, that the silicic magma has undergone substantial degassing in the middle to upper crust. The dacites that followed the eruption of the rhyolites by 54 years contain abundant crystals inherited from the crystal mush complement to the rhyolites and are one of the very few candidates for hybridization within the Hekla system.