Detrital zircon resolve longevity and evolution of silicic magmatism in extinct volcanic centers: A case study from the East Fjords of Iceland

Silicic magmatism is a relatively common feature of volcanic provinces dominated volumetrically by basaltic lavas. The bimodal association of abundant basalt with varying amounts of silicic magmatism is observed in a range of tectonomagmatic settings, including intracontinental hotspot tracks (e.g., Steens–Columbia River–Yellowstone, USA); continental rifts and large igneous provinces (e.g., Paraná-Etendeka, Namibia-Brazil; Africa-Arabian volcanic province, Ethiopia-Yemen; Karoo, southern Africa); and oceanic islands (e.g., Azores) and plateaus (e.g., Kerguelen) (e.g., Christiansen, 1984; Peate et al., 1992; Wallace et al., 2002; Snyder et al., 2004; I.U. Peate et al., 2005; Pankhurst et al., 2011; Coble and Mahood, 2012). Characterizing the silicic magmas and elucidating their genesis can offer insights into processes within large igneous provinces and the formation of juvenile continental crust (e.g., Pankhurst et al., 2011).

Iceland, located at the intersection of a mantle hotspot and the Mid-Atlantic Ridge, offers a unique example of a basaltic province that produces significant volumes of silicic magma. Rhyolite is estimated to compose ∼10% of the erupted products of the island’s volcanoes (e.g., Walker, 1966; Gunnarsson et al., 1998; Jónasson, 2007). The longevity of the dominantly mafic volcanic centers is poorly constrained because dating of basalts by radiometric methods is notoriously difficult due the general absence of minerals which can be precisely dated using U-Pb systematics. Zircon, a mineral that is very rare in mafic magmas but commonly present in rhyolitic melts, is especially valuable because it can provide both precise and accurate radiometric ages as well as elemental and isotopic tracers that provide petrogenetic information (e.g., I.U. Peate et al., 2005; Kemp et al., 2007; Claiborne et al., 2010a; Barboni and Schoene, 2014).

In this paper we evaluate the longevity and evolution of silicic magmatism at two large, closely spaced, dissected volcanic centers in northeastern Iceland using stream-sampled detrital zircon. We compare this life history to that inferred for central volcanoes in Iceland and, more generally, to silicic magmas within basaltic provinces globally. This approach, i.e., dating and chemical analysis of stream-sampled zircon, is an efficient way to evaluate the history of dissected volcanic centers.

Of the erupted material from Iceland’s volcanoes, ∼10% is rhyolite (e.g., Walker, 1966; Gunnarsson et al., 1998; Jónasson, 2007). Much of this material is concentrated in large central volcanoes, where the relative abundance of silicic material increases to ∼20%–30% of the exposed rocks (Walker, 1963; Gunnarsson et al., 1998). These major systems have been estimated to have a maximum lifespan of ∼0.5–1.5 m.y. (e.g., Saemundsson, 1979; Jakobsson, 1979; Thordarson and Larsen, 2007). There is ambiguity in the literature as to whether this estimate of longevity solely reflects active volcanism, or whether it includes the magmatic underpinnings of the system. Furthermore, these estimates are mostly based on K-Ar ages and magnetostratigraphy; zircon dates are notably lacking.

The most extensively investigated volcanic systems in Iceland are those that have recently erupted. Among these active volcanic systems, Torfajökull in southern Iceland (Fig. 1) has the largest areal exposure of silicic rock of any center (active or extinct) in Iceland, and it is the only one that is dominantly silicic (e.g., Gunnarsson et al., 1998). With an eruptive record spanning ∼400 k.y., it also has the longest well-documented history of any historically active Icelandic volcano (McGarvie et al., 2006). Kerlingarfjöll (∼300 k.y.; Flude et al., 2010), and Öræfajökull (Saemundsson, 1979) may span similar periods; Krafla has a somewhat shorter recognized history of ∼100 k.y. (Jónasson, 1994). The relatively young Ljósufjöll, on Snæfellsnes, is the only volcanic system that has a documented lifespan approaching 1 m.y., from ca. 1.1 Ma to after 100 ka (Flude et al., 2008). A recent study of the Austurhorn intrusive complex in southeastern Iceland (Padilla et al., 2016) revealed a duration of ∼300 k.y. for the zircon-saturated history of the system, suggesting comparable longevity of Icelandic intrusive and volcanic systems.

Volcanic systems whose silicic products erupted over a substantial timespan (>∼100 k.y.) include Torfajökull, Kerlingarfjöll, Thingmuli, and Ljósufjöll (Fig. 1). All of these volcanoes reveal substantial diversity in the compositions of silicic rocks, and distinct compositional series have been identified at Torfajökull (peralkaline versus subalkaline, metaluminous; McGarvie et al., 2006), Kerlingarfjöll (high-Nb versus low-Nb, all subalkaline, metaluminous; Flude et al., 2010), and Thingmuli (high-Fe versus low-Fe, all subalkaline, metaluminous; Charreteur et al., 2013). None of these volcanoes, however, exhibits a monotonic evolution through time from less to more silicic or vice versa (e.g., McGarvie et al., 2006; Flude et al., 2010; Charreteur et al., 2013). At Torfajökull, volcanism was initially peralkaline but shifted to subalkaline and metaluminous after ca. 70 ka. At Kerlingarfjöll, early low-Nb lavas gave way to high-Nb lavas ca. 250 ka. High- and low-Fe magmas erupted at Thingmuli with no discernible correlation between composition and time.

Prior to the ongoing magmatic activity at Torfajökull, the Borgarfjörđur eystri region in the northernmost East Fjords was the site of the most voluminous outpourings of silicic magma recognized in Iceland’s history (Gustafsson et al., 1989; Johannesson and Saemundsson, 2009; Berg, 2016; Vogler 2014). Four closely spaced volcanic centers, active ca. 13–12 Ma (Martin et al., 2011), comprise this terrane: Dyrfjoll, Breiðuvík, Kækjuskörð, and Herfell (Fig. 2). The landscape is dominated by pale-colored mountains, reflecting the silicic substrate that constitutes >20% of the region (Berg, 2016; Vogler, 2014; Figs. 3 and 4).

The concept of a volcanic system as proposed by Walker (e.g., 1958, 1963) for other Icelandic silicic centers is perhaps not strictly applicable in the Borgarfjörđur eystri region. Both the size of the silicic region and the abundance of differentiated rocks greatly exceed other localities, except for Torfajökull. Furthermore, the absence of a recognized fissure swarm, which typically characterizes volcanic systems elsewhere (cf. Thordarson and Larsen, 2007), suggests a unique magmatic-tectonic setting for Borgarfjörđur eystri. The volcanic systems described by Walker and others farther south in the East Fjords, and elsewhere in Iceland, are thought to have been principally formed in an extensional tectonic environment (Walker, 1958, 1963, 1964, 1966; Gibson and Walker, 1963; Carmichael, 1964). However, the Breiðuvík area may be related to a fracture zone and propagating rift segment of the modern Northern Volcanic Zone (NVZ) (Martin et al., 2011). For these reasons, we use the term volcanic center for each of the volcanic (silicic) sites within the Borgarfjörður eystri region (center < complex < system).

This study focuses on the volcanic centers Breiðuvík and Kækjuskörð. The surface geology in the field area is composed exclusively of extrusive rocks (e.g., Gustafsson et al., 1989; Vogler, 2014); intrusive bodies are not observed, other than the prevalent north-trending regional basalt dikes and high-level dikes crisscrossing the eroded silicic centers (e.g., cone sheets; cf. Dyrfjöll to the northwest; Burchardt et al., 2011). Breiðuvík, the largest of the silicic centers in the Borgarfjörður eystri region, is exposed over ∼30 km² and is a downsag caldera ∼10 km in diameter with estimated subsidence >600 m (Fig. 3; Vogler, 2014). Within this area, silicic rocks constitute as much as 35% of the bedrock; intermediate rocks compose ∼15%, and the remaining 50% is basalt (Vogler, 2014). A majority of the silicic rocks are lavas, but the most prominent unit is the large Hvítserkur ignimbrite, which likely filled the Breiðuvík caldera before extensive glacial erosion. Immediately to the west, the smaller Kækjuskörð center also exposes abundant silicic lavas and smaller ignimbrites. These units are overlain by another large (≥5 km³) ignimbrite (Gustafsson et al., 1989) that likely originated from the nearby Mount Herfell, southwest of Kækjuskörð.

Detrital zircon from modern sediments provides information about eroded terranes that may be otherwise difficult or impossible to investigate (e.g., Fedo et al., 2003; Moecher and Samson, 2006; Barth et al., 2013). We posit that strategically selected sediment samples from specific, confined, localities, like the flanks of a volcanic edifice, or within the bounds of a caldera, can provide an overview of the magmatic history of a specific volcanic center. A detrital study of this nature lacks the detail that comes with close observations of field relationships and targeted sampling of specific igneous units (e.g., Berg, 2016; Vogler, 2014; Sharman et al., 2015). However, detrital sampling has the benefit of minimizing human-introduced biases inherent to sampling in situ rock units (e.g., issues of desirability related to accessibility or alteration and decisions about selecting samples that are unique versus representative). The detrital zircon spectra could be strongly biased toward one or more zircon-rich units in the drainage basin. This overrepresentation is an unavoidable but acceptable risk in a reconnaissance study such as this.

We use trace elements, U-Pb dates, and O and Hf isotopes from detrital zircon sampled in streams that drain large exposures of Breiðuvík and Kækjuskörð to undertake a reconnaissance study of the two volcanic centers. The objectives of this detrital investigation are to: (1) determine if the closely spaced Breiðuvík and Kækjuskörð centers have distinguishable silicic histories; (2) estimate the total lifespan of silicic magmatism at large, extinct, Icelandic volcanoes; (3) investigate the character and evolution of silicic magmas at these volcanoes through time; and (4) evaluate the effectiveness and efficiency of stream sampling and analysis of zircon for resolving the evolution of eroded volcanic centers.

Rocks in the Borgarfjörður eystri region have undergone extensive metasomatic alteration by hydrothermal fluids (Vogler, 2014), limiting petrogenetic insight using whole-rock geochemical data. Therefore, we focus our analyses on zircon, a mineral that is generally resistant to hydrothermal alteration (Schmitt and Hulen, 2008; Watts et al., 2011; Bindeman et al., 2012; Milicich et al., 2013; Klemetti et al., 2014).

We collected sand samples (∼10 kg each) from two small rivers with drainage basins within the exposures of the Breiðuvík and neighboring Kækjuskörð volcano (Figs. 2 and 3). Sample ISS comes from the Stóraá River, a braided stream with a catchment entirely within the margins of the Breiðuvík caldera. The drainage area upstream of our sampling location encompasses the southern, western, and northwestern sectors of the caldera (>30% of the mapped caldera area). Approximately 50% of the catchment region is composed of silicic material. Sample ISKK was collected at the confluence of the Krossá and Kækjudalsá Rivers. The catchment area includes portions of prominent peaks associated with Breiðuvík (Hvítserkur and Hvítuhnjúkar) and minor exposures of the Herfell ignimbrite, but it is dominated by eruptive units from Kækjuskörð. Sample ISS therefore represents the Breiðuvík area while sample ISKK is more representative of the broader Borgarfjörður eystri region, specifically providing insight into the silicic material exposed at Kækjuskörð.

The catchment areas of the Stóraá and Krossá-Kækjudalsá Rivers drain only Breiðuvík and Kækjuskörð volcanic centers. Glacially derived sediment and modern river sediments in this confined region are locally derived, and assumed to originate only from the two volcanic centers (plus possible contributions from the local Herfell ignimbrite). Therefore, all zircon grains should relate directly to the magmatic history of Borgarfjörður eystri. No younger strata are present to complicate the detrital record. Any contributions from ancient (i.e., pre-Icelandic crustal contributions, e.g., Paquette et al., 2006, 2007; Torsvik et al., 2015) should be detectible by anomalously old zircon ages, high δ¹⁸O, or low ε_Hf values. Explosive silicic eruptions with extensive tephra distribution are known in Iceland, but we hope that the relative paucity of zircon in a typical Icelandic eruptive unit (e.g., Carley et al., 2011), coupled with aerial fractionation of dense zircon crystals, will minimize zircon contamination of the field area. Any zircon that did not originate in the region of study should be recognized by distinctive ages, ε_Hf or δ¹⁸O values.

In addition to river sands, we collected a representative pumice sample from the Hvítserkur ignimbrite within the Breiðuvík caldera. The Hvítserkur ignimbrite is the most prominent eruption product in this area, and we consider this pumice to be a representative whole-rock sample of Breiðuvík magma.

Zircon crystals were separated from river sands (sieved to <500 µm) using density (water and heavy liquid), magnetic, and hand-picking techniques. When hand-picking zircon grains, great care was taken to pick grains indiscriminately, to ensure that the diversity of zircon crystals present in the samples was preserved, without preference given to crystals of a certain size (large or small), shape (euhedral) or condition (e.g., broken, clear or colored). Grains were mounted in epoxy, ground and polished to expose crystal centers, and imaged by cathodoluminescence (CL) using a Tescan Vega 3 LM variable pressure scanning electron microscope at Vanderbilt University (Nashville, Tennessee). These CL images were used to guide ion microprobe and laser ablation analysis spot locations.

A portion (∼20 g) of an individual piece of pumice from the Hvítserkur ignimbrite was powdered in an agate mortar and pestle for whole-rock isotope analysis. Another small portion (∼20 g) was coarsely crushed using an agate mortar and pestle, and glass was removed by rinsing in concentrated hydrofluoric acid for 30 s at room temperature. Cleaned material was put through the standard density, magnetic, and hand-picking zircon separation steps. Two zircon crystals were found in the Hvítserkur pumice sample. These crystals were mounted, polished, imaged, and analyzed with the detrital zircon grains.

U-Pb ages were determined for 247 detrital zircon grains (83 from sample ISS and 164 from sample ISKK) using the sensitive high-resolution ion microprobe with reverse geometry (SHRIMP-RG) at Stanford University (Stanford, California). Zircon crystals separated from the Hvítserkur pumice were too small for isotopic analysis using the SHRIMP-RG. Analyses were performed using an O₂^– primary ion beam accelerated at 10 kV, with an intensity varying from 5.0 to 6.8 nA. The primary ion beam spot had a diameter between 25 and 28 µm and a depth of ∼2–3 µm. This high spatial resolution allowed us to target specific zones within zircon grains, guided by CL images. We preferentially targeted the interiors and rims of grains to capture the oldest and youngest zircon age populations. To efficiently do this assessment of age range, the majority of ages in this detrital study were collected using three scans (peak hopping cycles from mass 188 through 254) per analysis. A subset of ages was collected with five or six scans through the acquisition routine. The uncertainty of the calculated ages measured using the short scans is 2–3 times larger than for long scans. Measurement of U-Pb age standards (R33: 419 ± 1 Ma; Black et al., 2004) were done every 5^th analysis for 5 scan data, and every 6^th or 7^th for the 3 scan data. Data reduction for geochronology follows the methods described by Williams (1997) and Ireland and Williams (2003) and uses the Microsoft Excel add-in programs Squid 2.51 and Isoplot3.76 of Ludwig (2009, 2012). The data are not corrected for ²³⁰Th disequilibrium because adjustments to ages would be minor (<0.1 m.y.; within analytical uncertainty). Furthermore, in this detrital context, the necessary Th/U value of the parent melt is unconstrained and would only be an estimate. (Further information can be found in Supplemental Table S1¹.)

Elemental compositions were measured on 69 zircon crystals (32 from sample ISS, 35 from sample ISKK, and 2 from the Hvítserkur pumice) using the SHRIMP-RG. The primary beam energy ranged from ∼1.5 to 2.5 nA O₂^– and the analytical spots had a diameter of ∼15 µm and sputter depth of ∼1 µm. Trace element abundances were calculated relative to those on in-house standard MAD (Madagascar Green zircon; Barth and Wooden, 2010). We used methods following those described in Claiborne et al. (2006, 2010b), and Barth and Wooden (2010). Elemental compositions for Borgarfjörður eystri zircon (ISS and ISKK) were presented by Carley et al. (2014) as part of a larger data set. Because these compositions are very relevant to this study, they are also available in Supplemental Table S2².

Oxygen isotope analyses were collected from 134 zircon crystals (50 from sample ISS, 83 from sample ISKK, and 1 from the Hvítserkur pumice). Oxygen isotope ratios for Breiðuvík area zircon (ISS and ISKK) were previously included as part of a larger data compilation by Carley et al. (2014). Because these data are highly relevant to this study, they are included in Table 1. Additional information (e.g., analytical session, standards) can be found in Supplemental Table S3³.

Hf isotope compositions were measured for 71 detrital zircon crystals (30 for sample ISS and 41 for sample ISKK) using a Thermo-Scientific Neptune multicollector–inductively coupled plasma–mass spectrometer (MC-ICP-MS) interfaced to a GeoLas 193 laser ablation (LA) MC-ICP-MS, fired for 600 shots at 10 Hz and 5 J/cm² during a single analytical session at the Memorial University of Newfoundland (St. John’s). Zircon crystals separated from the Hvítserkur pumice were too small for Hf isotope analysis by LA-MC-ICP-MS. The laser ablation pits had a diameter of ∼50 µm, with a typical depth of ∼45 µm. The Hf isotopic composition of zircon crystals is the last measurement to be collected in the analytical sequence because the ∼50 µm spot and the depth of the ablation pit typically accounts for most of the analyzable space on Icelandic zircon crystals, and typically little remains after analysis is complete.

The analytical methods closely followed those of Fisher et al. (2011) and Souders et al. (2012), the only modification being the addition of 4 mL/min of nitrogen gas to improve the sensitivity of the MC-ICP-MS. Reference zircons (cycled between Plešovice, R33, MUNZirc Batch 142 and 144, and FC1) were analyzed from an in-house standard mount after every eight Icelandic zircon unknowns. Results are presented as ε_Hf values, which were calculated using the present-day CHUR (chondrite uniform reservoir) ¹⁷⁶Hf/¹⁷⁷Hf ratio of 0.282785 (Bouvier et al., 2008). Further information can be found in Supplemental Table S4⁴.

The bulk-rock Hf isotopic composition of a single piece of pumice from the Hvítserkur ignimbrite was determined at the Radiogenic Isotope and Geochronology Laboratory (RIGL) at Washington State University (WSU, Pullman), using the RIGL-WSU Thermo Finnigan MC-ICP-MS and following methods described in detail by McDowell et al. (2016). Analytical results were corrected for mass fractionation using ¹⁷⁹Hf/¹⁷⁷Hf = 0.7325 and normalized using Hf standard JMC-475 (¹⁷⁶Hf/¹⁷⁷Hf = 0.282160). Further information can be found in Supplemental Table S5⁵.

The 83 dated zircon grains from sample ISS (Breiðuvík only) have a median age of 12.9 Ma, and a weighted mean age of 13.0 ± 0.1 Ma (mean square of weighted deviates, MSWD = 3.4; Fig. 5). The youngest and oldest individual ages measured in this zircon population are 11.2 ± 0.8 Ma and 14.5 ± 0.4 Ma, respectively (individual age errors are 1σ). The 5 youngest ISS zircon crystals have a weighted mean age of 11.8 ± 0.4 Ma (MSWD = 0.28), while the oldest 5 ISS crystals had a weighted mean of 14.3 ± 0.3 Ma (MSWD = 0.17). The calculated errors for weighted means do not account for MSWD values; the MSWD values are interpreted to reflect scatter in the age of the zircon population and not analytical error.

The 164 dated zircon crystals from sample ISKK (primarily Kækjuskörð drainage, partially Breiðuvík) have a median age of 12.9 Ma and a weighted mean age of 12.9 Ma ± 0.1 (MSWD = 3.1). The youngest and oldest individual zircon crystals are 11.2 ± 0.7 Ma and 15.0 Ma ± 0.9, respectively, while the weighted means of the youngest and oldest 5 grains are 11.3 ± 0.7 Ma (MSWD = 0.04) and 14.7 ± 0.6 Ma (0.23), respectively.

Elemental compositions of detrital zircon grains from Breiðuvík-Kækjuskörð are highly variable, but define a coherent population. The most obvious characteristics generally reflect those of magmatic zircon. They are highly enriched in heavy rare earth elements and Y relative to light rare earth elements; have relatively large positive Ce anomalies and relatively moderate negative Eu anomalies (Fig. 6); and they have very high Hf, U, and Th concentrations (e.g., Hoskin and Schaltegger, 2003; Claiborne et al., 2006, 2010b; Grimes et al., 2007, 2015; Carley et al., 2014; Fig. 7; Table S2 [see ^{footnote 2}]).

Rare earth element concentrations vary by a factor of 10–100 (La 0.1–10 × chondrite, Lu 10³-10⁴ × chondrite), but patterns are very similar between the two samples (Fig. 6). Likewise, U and Th are high and quite variable (∼20–1000 ppm and 10–1100 ppm, respectively; Fig. 7), with U concentrations covarying with Th (Th/U = 0.25–1.5) and Nb (∼4–400 ppm).

The concentrations of Hf in Breiðuvík-Kækjuskörð zircon crystals range from 7000 to 14,000 ppm, and correlate negatively with Ti (Fig. 8), which has concentrations of 4–36 ppm. Hf and Ti concentrations define two populations: ∼60% of the analyses have Hf >10,000 ppm and Ti <10 ppm, while the rest have Hf <10,000 ppm and Ti >10 ppm. The wide range of Ti suggests crystallization at variable temperatures (Watson and Harrison, 2005; Ferry and Watson, 2007). The median (7.5 ppm) and mean (10.8 ppm) of Ti are relatively high by global standards but somewhat lower than average for Iceland (Fu et al., 2008; Carley, 2014).

The ISS zircon population (n = 50) has δ¹⁸O values that range from 0.7‰ ± 0.6‰ to 3.7‰ ± 0.4‰ (errors are 1σ), with a median value of 3.1‰ and a weighted mean of 3.1‰ ± 0.1‰ (MSWD = 0.77; Fig. 9). The 5 ISS analyses with the lowest δ¹⁸O values have a weighted mean of 2.4‰ ± 0.02‰ (MSWD = 2.8); the 5 highest analyses have a weighted mean of 3.5‰ ± 0.4‰ (MSWD = 0.03). The δ¹⁸O of ISKK zircon (n = 83) ranges from 1.3‰ ± 0.6‰ to 4.2 ± 0.4‰, with a median value of 3.3‰ and a weighted mean of 3.5‰ ± 0.1‰ (MSWD = 2.1). The weighted mean for the 5 lowest values is 2.1‰ ± 0.8‰ (MSWD = 2.4); for the 5 highest values, the weighted mean is 3.9‰ ± 0.2‰ (MSWD = 0.15).

These oxygen values are consistent with the range of 1.1‰–4.7‰ reported by Berg (2016) for Borgarfjörður Eystri; they are significantly lower than whole-rock analyses of devitrified ignimbrite from the Borgarfjörður eystri region (Martin and Sigmarsson, 2010; Berg, 2016). The heavier oxygen signatures in the bulk rock samples are probably a consequence of post-eruptive, low-temperature, metasomatic alteration of the devitrified ignimbrite. The discrepancy in δ¹⁸O demonstrates the inherent complications in O isotope analyses of altered whole-rock samples, and emphasizes the robustness of zircon in recording the primary magmatic conditions (Valley, 2003).

The range of ε_Hf for the ISS zircon population (n = 31) is +13.2 ± 0.5 to +16.8 ± 1.0 (1 σ standard error). The 31 analyses have a median ε_Hf of 14.7 and weighted mean of 14.8 ± 0.3 (MSWD = 2.8; Fig. 10). The 5 analyses with the lowest ε_Hf have a weighted mean of 13.9 ± 0.4 (MSWD = 0.68), and 16.2 ± 0.5 (MSWD = 0.47) for the 5 highest analyses.

The 42 Hf analyses done on sample ISKK zircon grains have a minimum ε_Hf value of 14.7 ± 1.1, a maximum of 16.3 ± 1.6, and a median of 14.7. The zircon grains analyzed have a weighted mean of 14.6 ± 0.2 ε_Hf (MSWD = 1.9). The 5 analyses with the lowest ε_Hf have a weighted mean of 13.6 ± 0.4 (MSWD = 0.9) and 15.8 ± 0.5 (MSWD = 0.39) for the 5 highest analyses.

A high-precision solution analysis of the Hvítserkur pumice yielded a bulk-rock ε_Hf value of 15.1 ± 0.1 (1σ).

Similarities in zircon populations from samples ISS (Breiðuvík) and ISKK (mostly Kækjuskörð, contributions from Breiðuvík) suggest that zircon from these two closely spaced volcanic centers shared similar magmatic histories in terms of source rock, magmatic processes, and timing of magmatic activity and zircon crystallization. Both samples are zircon rich compared to other Icelandic volcanic and detrital samples (e.g., Carley et al., 2011, 2014; Carley, 2014). For ISS and ISKK, the weighted means of age (13.0 ± 0.1 Ma versus 12.9 ± 0.1 Ma) and Hf isotopes (14.8 ± 0.3 ε_Hf versus 14.6 ± 0.2 ε_Hf; Fig. 11) are within error of one another, and they are also very similar in weighted mean of δ¹⁸O (3.1‰ ± 0.1‰ versus 3.5‰ ± 0.1‰).

Similarities between ages, O isotopes, and Hf isotopes in zircon populations from samples ISS and ISKK were quantified using Kolmogorov-Smirnov (K-S) testing (Fletcher et al., 2007). In a K-S test, populations are considered to be subsamples from the same distribution as long as the probability values (P) are >0.05. In the case of Breiduvik-Kækjuskörð, the K-S tests yielded high probabilities of similarity between samples ISS and ISKK for age (P = 1.00), O isotopes (P = 0.31) and Hf isotopes (P = 0.95). Other statistical tests (e.g., cross-correlation coefficient, similarity coefficient, and likeness coefficient; Saylor and Sundell, 2016) lead essentially to the same conclusion. Although K-S testing reveals striking similarities between zircon populations from these two drainages, there is no discernible correlation between Hf or O and age (Fig. 12).

Despite the fact that the detrital zircon samples represent mixed populations of volcanic phenocrysts, which potentially include antecrysts and xenocrysts, they define a coherent range of ε_Hf and δ¹⁸O compositions, notable when compared to those for other Icelandic systems (Fig. 13). However, the statistical overlap in ages and isotopic compositions between ISS and ISKK suggests that zircon from these two volcanic centers were sourced from a system or systems that generated silicic magmas from the same or similar source rocks, influenced by similar modification processes in the crust, during the same prolonged period (discussed in the Geochemical Character, Petrogenesis, and Evolution of Silicic Magmas at Breiðuvík-Kækjuskörð section). Furthermore, the weighted means for ε_Hf for ISS (14.8 ± 0.3) and ISKK zircon (14.6 ± 0.2) are within 2σ error of the bulk-rock ε_Hf of Hvítserkur pumice (15.1 ± 0.1), which supports the interpretation that detrital zircon crystals provide a robust representation of the system that produced these volcanoes.

Published ages for the Borgarfjörður eystri region range from 13.5 to 12.1 Ma (Berg, 2016; Martin et al., 2011; cf. Paquette et al., 2006, 2007). The difference between these ages is within the 0.5–1.5 m.y. age range that is typically associated with the lifespan of a volcanic system in Iceland. However, this detrital study reveals U-Pb zircon ages from Breiðuvík-Kækjuskörð that range from ca. 15.0 to 11.2 Ma, suggesting a substantially longer period of magmatism during which zircon crystallized (Figs. 5, 11, 12, 14, and 15; Tables 1 and 2).

To assess the implications of the spread of dates obtained in our data set, we first evaluate whether the wide range of measured ages reflects error-based scatter. The MSWD value for ²⁰⁶Pb/²³⁸U ages is 3.4 for ISS (83 analyses) and 3.1 for ISKK (164 analyses). These high MSWD values preclude the possibility of a single age population of zircon (Mahon, 1996). Thus, we infer that the range of zircon ages within each sample reflects a prolonged magmatic history during which zircon growth occurred. The continuous spread in ages (Figs. 5 and 15) could reflect (1) nearly continuous zircon growth within an evolving magmatic system; or (2) multiple episodes of punctuated, short-lived zircon growth (e.g., Bindeman and Melnik, 2016) that are shorter than the analytical precision of individual analyses. The analytical limitations of the current data set prevent us from confidently ruling out either alternative.

To quantify the longevity of the zircon-saturated history of this magmatic system (whether punctuated or continuous), we first combine ages for samples ISS and ISKK to achieve a greater n value, and thus a more statistically robust data set. This tactic is justified by isotopic similarities in age and other isotopes (Fig. 11; Table 2). Similar to the individual samples, the combined zircon population has a weighted mean age of 12.9 ± 0.1 Ma (MSWD = 3.2, n = 247) and a probability density curve that closely resembles a normal distribution (Fig. 15).

The longevity suggested by this combined zircon population may be >3 m.y., or more conservatively ≥2.8 m.y., depending on the approach used to evaluate the observed spread of ages. The absolute difference between the minimum (11.2 ± 0.7 Ma) and maximum (15.0 ± 0.9 Ma) ages for individual zircon crystals yields a lifespan of 3.8 ± 1.1 m.y. Limiting consideration to measured ages with analytical precision <0.5 m.y., the difference between the high-precision minimum (11.5 ± 0.3 Ma) and maximum (14.9 ± 0.4 Ma) is 3.4 ± 0.5 m.y. The difference between weighted mean ages for the 5 oldest and 5 youngest ages in our combined data set supports an age span of 3.4 ± 0.9 m.y. The time captured by two standard deviations (1.4 m.y.) about the mean (12.9 Ma) of the combined ISS and ISKK data set is ≥2.8 m.y., the most conservative assessment of magmatic longevity in this study.

Whether the detrital zircon ages from eroding volcanic terranes at Breiðuvík-Kækjuskörð represent a record of continuous or punctuated zircon saturation, the ≥2.8 m.y. history convincingly documents greater longevity than has been previously recognized in Iceland. The previously proposed maximum lifetime of an Icelandic volcanic system (≥1.5 m.y.) has essentially been doubled (e.g., Saemundsson, 1979; Jakobsson, 1979; Thordarson and Larsen, 2007).

To evaluate the confidence that we can place in this estimated duration of silicic magmatic activity at Breiðuvík-Kækjuskörð, we performed Monte Carlo modeling in MATLAB (https://www.mathworks.com/). In our modeling, we theoretically replicate our detrital zircon sampling and dating by randomly selecting ages and errors from within their respective distributions (in this case, normal distribution for ages, gamma distribution for errors). We tested sample sizes of 83 and 247 based on the number of ages we measured for the individual sample ISS and for the combined ISS and ISKK samples, respectively. These simulations were repeated many times (i.e., 100× , 1000×). Modeling results are summarized in Table 2, and details of our procedures are described in Supplemental Item S1⁶.

Using Monte Carlo modeling to emulate our sampling and dating of sample ISS (n = 83), we determine that if we were to repeat this study with a similar sample size there is a 97% likelihood that we would find a zircon that exceeds a 3 m.y. time span about the mean age (13.0 Ma). In revisiting the sampling of all dated Breiðuvík-Kækjuskörð zircon (ISS + ISKK, n = 247), we find that there is a 98% likelihood that we would find a zircon that exceeds a 3 m.y. time span about the mean age (13.0 Ma). Therefore, the Monte Carlo model supports the assessment that the age span captured by zircon at Breiðuvík-Kækjuskörð is ≥2.8 m.y.

Our results suggest that Breiðuvík-Kækjuskörð detrital zircon grains document a longer magmatic history than has been reported for any other Icelandic volcanic system. Zircon-saturated silicic activity likely spanned an interval of ≥2.8 m.y., from ca. 14.5 to 11.5 Ma. While crystallization was active throughout this prolonged interval of time, the majority of zircon growth likely occurred during a shorter window of time that better aligns with previous estimates of Icelandic volcano longevity: an interval of 1.4 m.y. (2σ standard deviations) centered around the mean age of 13 Ma.

We contend that our data provide a useful estimate of longevity of the silicic portion of the Breiðuvík-Kækjuskörð system, but this does not necessarily imply that we have defined the lifetime of the system as a whole. Two possible complications are evident: (1) basaltic magmatism may have been more protracted than silicic magmatism, initiating earlier and/or terminating later, but left no recognizable zircon record; or (2) the eruptive portion of the Breiðuvík-Kækjuskörð history may have been shorter than suggested by the zircon record if a prolonged period of silicic intrusion preceded the first eruptions.

We consider the first possibility (prolonged basaltic magmatism) to be likely. In Iceland and other basalt-dominated provinces where volcanic systems also produce silicic magma, dominantly or entirely mafic eruptions tend to give way late in a center’s history to mixed or dominantly silicic eruptions (see following discussion, “Breiðuvík-Kækjuskörð in a Global Context”). This is consistent with the inference that silicic magmatism in these settings is either a consequence of mafic magma undergoing extreme fractionation, or mafic magma prompting anatexis of basalt-heated crust.

The second possibility (zircon captures a long intrusive but short eruptive history) is also plausible. It is common for antecrystic zircon to crystallize early in plutonic portions of systems, before entrainment in later eruptive events (e.g., Charlier et al., 2005; Claiborne et al., 2010a; Schmitt et al., 2011). However, we doubt that this could account for a substantial portion of the ≥2.8 m.y. history of zircon ages observed at Breiðuvík-Kækjuskörð. Documented examples of erupted antecrysts are rarely more than a few hundred thousand years older than the eruption in which they reached the surface. Furthermore, the zircons in this study do not have characteristics that suggest growth and storage in a plutonic environment and subsequent entrainment. They are typically small and acicular (<∼150 µm long) with simple zoning, and few or none show evidence of resorption (Supplemental Item S2⁷).

Geochemical Character, Petrogenesis, and Evolution of Silicic Magmas at Breiðuvík-Kækjuskörð

Breiðuvík-Kækjuskörð zircon compositions are variable but within the established elemental and isotopic range of Icelandic zircon (Figs. 13, 16, and 17; Carley et al., 2014). Similarities reveal that Breiðuvík-Kækjuskörð magmas broadly share petrogenetic origins with those generated throughout Iceland’s history. Elemental and isotopic compositions of Breiðuvík-Kækjuskörð zircon show no discernible change through time (Fig. 12; Supplemental Item S3⁸), suggesting negligible changes in petrogenesis of silicic magmas over the long history of the centers, in terms of either sources or magma evolution.

The wide range in zircon trace element concentrations, including an abundance of relatively low Ti, high Hf analyses that cluster near the most evolved end of the Iceland range, suggests considerable fractionation of the zircon-saturated host magmas (Fig. 16; cf. Carley et al., 2014; Grimes et al., 2015; Claiborne et al., 2006, 2010b). Time independence of compositions, however, suggests repeating cycles of fractionation rather than a continuing monotonic process (see Supplemental Item 3 [^{footnote 8}]).

Narrow ranges of ε_Hf and δ¹⁸O indicate little variability in source materials that contributed to the Breiðuvík-Kækjuskörð magmas. At ∼+15, ε_Hf for Breiðuvík-Kækjuskörð zircon is above average among Icelandic volcanic zircon (Fig. 13; Carley, 2014). This ε_Hf is, however, typical of rift zone basalts and matches closely values for basalts of similar age from this region (Eastern Iceland rift zone; Kitagawa et al., 2008; Peate et al., 2010). These high ε_Hf values indicate that the Breiðuvík-Kækjuskörð magmas were derived from a depleted mantle source that was more similar to the sources of mid-oceanic ridge basalt and Iceland rift basalts than to the relatively enriched mantle proposed for Icelandic off-rift-axis basalts and rhyolites. The tight clustering of δ¹⁸O values near the Icelandic zircon median of +3.2‰ (Carley et al., 2014) suggests that either the source material was mildly, yet uniformly, hydrothermally altered crust, or that parental magmas were uniformly contaminated by more highly altered crustal material (cf. Carley et al., 2014; Padilla et al., 2016).

The compositional similarities of zircon from both Breiðuvík and Kækjuskörð volcanic centers, and their identical age distributions, strongly suggest a common magmatic and volcanic origin. It is likely that Breiðuvík and Kækjuskörð together are part of a single volcanic complex (center < complex < system).

Silicic Magmas in the Breiðuvík-Kækjuskörð Complex Compared to those in Active Icelandic Volcanic Systems

Zircon compositions from three active systems, Krafla, Öræfajökull, and Torfajökull, are compared to the extinct Breiðuvík-Kækjuskörð complex (Figs. 1, 13, 16, and 17). Krafla is an active volcanic system situated on the actively rifting Northern Volcanic Zone (Fig. 1). Öræfajökull had the most voluminous explosive eruption of notably homogeneous silicic magma in Iceland’s recorded history, in A.D. 1362 (Thorarinsson, 1958; Larsen et al., 1999; Selbekk and Tronnes, 2007; Sharma et al., 2008). Its explosivity and homogeneity prompt a comparison with Breiðuvík-Kækjuskörð. Torfajökull is included in the comparison as the only system in Iceland that has erupted as much or more silicic material than the Borgarfjörður eystri region, which includes Breiðuvík and Kækjuskörð (Gustafsson et al., 1989; Johannesson and Saemundsson, 2009; Berg, 2016; Vogler, 2014).

The δ¹⁸O values of Breiðuvík-Kækjuskörð zircons are within a narrow range (∼3.5‰–4‰) relative to those from Krafla (∼0.5‰–2‰), Öræfajökull (∼2.5‰–5‰), and Torfajökull (∼0.5‰–3‰; Fig. 13). This coherency indicates little variability in the contributions of altered crustal source materials to the Breiðuvík-Kækjuskörð magmas throughout the long (likely ≥2.8 m.y.) history of the complex relative to the three modern systems. The lighter oxygen signature at both Krafla and Torfajökull (Fig. 13) suggests a greater influence of altered crustal material at these active rift zone volcanoes than at Breiðuvík-Kækjuskörð.

The ε_Hf of Breiðuvík-Kækjuskörð zircon grains (ε_Hf ∼13.5–16.5) are among the highest Hf isotope values yet discovered for Icelandic zircon (Fig. 13), similar only to zircon from Krafla (ε_Hf ∼13.5–17). These two zircon populations, which share a connection with the Northern Volcanic Zone, are markedly more radiogenic than those of Torfajökull (ε_Hf ∼10–12) and Öræfajökull (ε_Hf ∼10.5–15; Figs. 1 and 13). Torfajökull is situated on the propagating Eastern Volcanic Zone and Öræfajökull is on the off-rift Öræfi Volcanic Belt (Fig. 1).

Elemental compositions reveal that the crystallization conditions of zircon at Breiðuvík-Kækjuskörð, while not unusual for Icelandic systems, differ in detail from those at the active systems in this comparison. For similar Nb concentrations, Breiðuvík-Kækjuskörð zircons tend to have lower values of Th/U than those from Torfajökull, Krafla, or Öræfajökull (Fig. 16). Breiðuvík-Kækjuskörð zircon compositions also differ from the three active systems in extending to markedly higher Hf (suggestive of more highly fractionated magmas) and lower Ti (indicating lower crystallization temperatures) (Figs. 16 and 17; Carley et al., 2011, 2014).

The estimated lifespan of silicic magmatism at Breiðuvík-Kækjuskörð is consistent with eruptive durations reported for silicic magmas in systems dominated by basalt elsewhere on Earth. Well-constrained estimates for individual systems range from ∼1 to 3 m.y. (Steens–Columbia River–Yellowstone, USA, large rhyolite caldera systems; e.g., Christiansen 1984; Coble and Mahood 2012, 2016; Benson and Mahood, 2016; Karoo rhyolites; e.g., Riley et al., 2004; and Yemen rhyolites; I.U. Peate et al., 2005). Silicic volcanic rocks in these systems are commonly interbedded with basalt, but increase in abundance toward the upper part of the volcanic sections (e.g., Pankhurst et al., 2011). Basalt-rhyolite sequences in Iceland follow a similar pattern, with silicic volcanism appearing, and often dominating, the late eruptive sequence at systems such as Krafla (Jónasson, 1994), Hekla (Sigmarsson et al., 1992; Sverrisdottir, 2007), Kerlingarfjöll (Flude et al., 2010), and Thingmuli (Charreteur et al., 2013). Despite this similar temporal relationship between basalt and rhyolite in mafic-dominant provinces, this study of Breiðuvík-Kækjuskörð is the first to suggest a longevity for Icelandic magma systems that is comparable to those that are typical globally.

A focused detrital zircon study, as presented here for Breiðuvík-Kækjuskörð, is a powerful approach for illuminating a comprehensive history of silicic magmatism at an individual volcanic center or in a focused eruptive region. This is particularly true when the nuanced details provided by individual eruptive units and their relationships are not critical to the regional-scale history. This detrital approach combats sampling bias (e.g., limiting numbers of samples for analytical practicalities, preferentially selecting samples that are easily accessible) that might otherwise obscure important chapters of a system’s history and general geochemical characteristics. Furthermore, this sampling technique is essentially nondestructive to the natural landscape, and it is effective and efficient for understanding the magmatic underpinnings of the geology of a drainage basin (particularly when measured results are further corroborated by modeling, such as with the Monte Carlo approach to assessing age heterogeneity exercised in this study). A focused detrital accessory mineral approach to investigating silicic systems has widespread applicability. It can be applied globally in locations where rocks have undergone pervasive alteration, where accessibility is a barrier (particularly for systems that are remote, rugged, or obscured by glacial ice), or where silicic histories are overwhelmed in dominantly basaltic terranes (e.g., ocean islands, flood basalt provinces).

Constraining the longevity of Icelandic volcanoes has important implications for better understanding rhyolitic centers that are active, and potentially hazardous. Before now, the approximate lifetime of Icelandic volcanic systems was estimated to be 0.5–1.5 m.y. (Saemundsson, 1979), based on K-Ar dating and correlation with the paleomagnetic time scale. Zircon geochronology adds precise, robust age data to Iceland’s geologic record that was not achievable by earlier methods. Extinct silicic complexes such as Breiðuvík-Kækjuskörð provide a lifecycle perspective that is unachievable with volcanic systems that are still active today. With a zircon record that captures ≥2.8 m.y. of silicic activity, we demonstrate that Breiðuvík-Kækjuskörð underwent magmatic activity longer than has ever been previously reported for an Icelandic volcano. It also establishes that Icelandic volcanoes are capable of achieving magmatic longevity comparable to other long-lived silicic systems in otherwise mafic settings, such as Karoo (Pankhurst et al., 2011) and Yemen (I.U. Peate et al., 2005).

If Breiðuvík-Kækjuskörð is representative of silicic systems in Iceland, it is important that we recalibrate our assessment of magmatic lifespans in Iceland. If Breiðuvík is unique in Iceland, it begs the question, why is Breiðuvík-Kækjuskörð different from typical silicic systems in Iceland, and does its unique nature (including its longevity) explain the great concentration of silicic magma found in the region?

This work was supported by funding from National Science Foundation (NSF) Graduate Research fellowship DGE-0909667 to Carley, NSF Research Grant EAR-1220523 to Miller, National Geographic Society Grant 9030–12 to Miller, Rannís (Rannsóknamiðstöð Íslands) Iceland Science Fund Grant 130308–052 to Sigmarsson, and a Natural Sciences and Engineering Research Council of Canada Discovery Grant to Hanchar. The ion microprobe facility at University of California, Los Angeles is partly supported by a grant from the Instrumentation and Facilities Program, Division of Earth Sciences, NSF. We thank Erla Vogler for her guidance and support in the field, John Wilson (Lafayette College) for assistance with figures, David Furbish (Vanderbilt University) for his help in developing the Monte Carlo code, Oscar Lovera for implementation of the Kolmogorov–Smirnov code, and Joy Zhou (Lafayette College) for help executing code for additional statistical analyses. We greatly appreciate the analytical assistance provided by Aaron Covey (Vanderbilt University), Joe Wooden and Brad Ito (Stanford University), Rebecca Lam (Memorial University of Newfoundland), Jeff Vervoort, Diane Wilford, and Charles Knaack (Washington State University). We also thank editor Shanaka deSilva as well as Erik Klemetti and an anonymous reviewer for their thoughtful suggestions and improvements to this manuscript.