Magma reservoir dynamics at Toba caldera, Indonesia, recorded by oxygen isotope zoning in quartz

Quartz is a common phase in high-silica igneous rocks and is resistant to post-eruptive alteration, thus offering a reliable record of magmatic processes in silicic magma systems. Here we employ the 75 ka Toba super-eruption as a case study to show that quartz can resolve late-stage temporal changes in magmatic δ18O values. Overall, Toba quartz crystals exhibit comparatively high δ18O values, up to 10.2‰, due to magma residence within, and assimilation of, local granite basement. However, some 40% of the analysed quartz crystals display a decrease in δ18O values in outermost growth zones compared to their cores, with values as low as 6.7‰ (maximum ∆core−rim = 1.8‰). These lower values are consistent with the limited zircon record available for Toba, and the crystallisation history of Toba quartz traces an influx of a low-δ18O component into the magma reservoir just prior to eruption. Here we argue that this late-stage low-δ18O component is derived from hydrothermally-altered roof material. Our study demonstrates that quartz isotope stratigraphy can resolve magmatic events that may remain undetected by whole-rock or zircon isotope studies, and that assimilation of altered roof material may represent a viable eruption trigger in large Toba-style magmatic systems.

prior to eruption (cf. ref. 12). Indeed, Matthews et al. (ref. 12) previously applied Ti-diffusion chronometry to quartz crystals from the Young Toba Tuff (YTT) and suggested that short, episodic magma recharge events on a ca. 100-year timescale were important in priming the large-volume YTT magma reservoir for eruption. In the current paper, we present detailed SIMS oxygen isotope analysis of quartz crystals from the YTT, including several of the quartz crystals analysed by ref. 12 for diffusion chronometry to i) test for magma recharge, mixing, and pre-eruptive crustal recycling in the YTT reservoir as recorded in the crystal-scale oxygen isotope record, and ii) assess the utility of oxygen isotope stratigraphy in quartz as a tracer of pre-eruptive magmatic processes in large-volume rhyolitic systems.

Geological setting
The Toba caldera, north Sumatra ( Fig. 1), is thought to originate from subduction of the Investigator Ridge transform fault sited on the Indo-Australian Plate beneath the Sunda Trench, leading to formation of a deep crustal hot zone 13,14 that feeds the Toba volcano. The Toba system has produced four major explosive eruptions over the last 1.2 million years 15,16 , the oldest being the 35 km 3 Harrangoal Dacite Tuff (HDT). This event was followed by the quartz-rich 1,000 km 3 Old Toba Tuff (OTT) ~0.79 Ma 17 , and the 60 km 3 rhyolitic Middle Toba Tuff (MTT) ~0.5 Ma 18 . Finally, the cataclysmic 2,800 km 3 rhyolitic Youngest Toba Tuff (YTT) erupted ca. 75 ka 19 . The YTT deposits extend over 20,000 km 2 on Sumatra and are recorded as far as 2,000 km away (e.g. in India) 20 . The Toba super-eruption is the largest and best-known of the late Quaternary, and was speculated to have brought mankind to near extinction 21 .
High-silica caldera systems like Toba (Indonesia), Yellowstone (USA) and Taupo (New Zealand) produce some of the largest volume eruptions on Earth 22 . Recent work on magma petrogenesis of the Yellowstone hot spot track, for example, has revealed oxygen isotope diversity in zircon that includes both high-and low-δ 18 O values relative to 'normal' rhyolites 9 . For instance, the higher δ 18 O values in Yellowstone zircon are thought to represent assimilation of unaltered crust by ascending magmas, whereas the lower end of the δ 18 O zircon spectrum has been attributed to recycling of hydrothermally-altered (low-δ 18 O) rhyolitic material into the Yellowstone magma system. Similar intra-crystalline isotopic diversity is observed in zircon from other high-silica magmatic systems worldwide, including Toba, where comparable processes are hypothesised to have been in operation ( Fig. 2 and ref. 9). To fully resolve the complexity of the Toba crystal record, however, isotopic heterogeneities need to be resolved in more detail. We therefore use Toba as a case-study to assess the utility of quartz as a tracer of pre-eruptive magmatic processes, especially since Toba eruptive rocks are particularly quartz rich 23,24 . Furthermore, widespread low-δ 18 O crust is not expected at Toba's latitude 25 , meaning that we can also test for the role of late-stage additions of volcanically produced high-temperature hydrothermally-altered crust in driving large silicic eruptions.

Petrography and geochemistry
The YTT samples used in this study are from stratigraphically constrained tuff units exposed along the caldera margin and from outcrops surrounding the Toba caldera. Bulk tuff and pumice samples exhibit no signs The red ellipse designates the crustal swell of the Batak Tumor 15 , upon which the Toba caldera sits (in blue). The topographic swell is believed to result from upward pressure from the buoyant Toba magma system 15  To provide comparative data for our SIMS approach, we also performed oxygen isotope analyses on YTT whole-rock tuff and pumice clasts by conventional fluorination and on separated crystal phases by bulk crystal laser fluorination (LF) (see Methods for details of analytical procedures). The whole-rock data have δ 18 O values between 8.2 and 9.9‰ (x = 9.0 ± 0.15‰, 1σ , n = 11; Supplementary Table 2), whereas the older HDT (1.2 My; 61 wt.% SiO 2 ) has a δ 18 O whole-rock value of 8.6‰ (n = 1). The YTT bulk crystal analyses by LF give δ 18 O values of 8.2 to 9.3‰ for quartz (x = 8.9 ± 0.15‰, 1σ , n = 13), 7.6 to 8.1‰ for feldspar (x = 7.8 ± 0.15‰, 1σ , n = 6), 6.8 to 8.9‰ for biotite (x = 7.2 ± 0.15‰, 1σ , n = 8), and 5.5 to 6.4‰ for amphibole (x = 5.9 ± 0.15‰, 1σ , n = 4). In comparison, the δ 18 O quartz values from basement granitoids near Toba range from 6.4 to 9.5‰ (± 0.15‰, 1σ , n = 5).
Quartz crystals in the YTT are typically fragmented due to the high explosivity of the eruption 26 , but a record of multiple dissolution and reprecipitation events can nevertheless be discerned (e.g., ref. 24). To resolve spatial δ 18 O zonation within individual quartz crystals, we coupled cathodoluminescence (CL) and SIMS oxygen isotope analyses of 15 quartz crystals from pumice (reflecting pre-eruptive magma batches) and from bulk tuff (encapsulating mixed material on eruption). We present a total of 191 individual SIMS spot analyses that define a range of δ 18 O quartz from 6.7 to 10.2‰ (x = 8.8 ± 0.4‰, 1σ ; Supplementary Table 3). In order to assess the meaningfulness of our data in a wider Toba context, we furthermore analysed two clinopyroxene samples from the 1.2 Ma HDT for their 3 He/ 4 He ratios. We note that the YTT is virtually devoid of pyroxene (see ref. 23) and hence it was unfeasible to extract sufficient crystals for helium isotope analysis. The 3 He/ 4 He values obtained (0.67 R A and 1.76 R A , where R A = air 3 He/ 4 He) are lower than typical mantle values (8 ± 1 R A ) or contemporary geothermal fluid values at Toba (6.6 R A ; ref. 27), and instead approach crustal values (< 0.5 R A ). Moreover, one whole-rock YTT sample was also analysed for Sr-Nd isotopes and yielded an 87 Sr/ 86 Sr value of 0.714033 (± 13) and a 143 Nd/ 144 Nd value of 0.512113 (± 13). These radiogenic isotope data, although few in number, are consistent with previous studies that used radiogenic isotopes and which suggested that substantial crustal material was indeed involved in YTT petrogenesis 28,29 .

Resolving crustal differentiation processes at Toba
As is often the case in high-resolution analytical approaches 2,30 , our new SIMS data extend the established δ 18 O range for the YTT crystal data towards both higher and lower values, with an overall range of 3.5‰ for the full set of YTT quartz crystals (Fig. 3). Employing a crystal-melt fractionation factor (Δ 18 O quartz-rhyolite ) of 0.5‰ (e.g., ref. 31), the δ 18 O value of the Toba magma(s) can be calculated from the derived crystal data. The average calculated δ 18 O magma value is 8.3 ± 0.4‰ (n = 191, 1σ ) for SIMS data and 8.4 ± 0.2‰ (n = 13, 1σ ) for the LF data. The datasets show significant overlap, thus ruling out standardisation bias. Moreover, they overlap with calculated δ 18 O magma values determined using YTT zircon data from the literature 9 (average δ 18 O magma value = 8.6 ± 0.2‰; n = 50, 2σ ) ( Fig. 3). We note that the bulk crystal conventional fluorination quartz data (ref. 32; data obtained pre-1990) produce a slightly higher average δ 18 O magma value of 9.1 ± 0.2‰ (n = 12, 1σ ) compared to our bulk quartz data (Fig. 3). This is likely due to our use of the more modern laser fluorination technique that provides more accurate results. Nevertheless, it is important to note that all YTT δ 18 O quartz and corresponding magma values are elevated relative to primitive Sunda arc basalt, which is expected to have δ 18 O values between 5.5 and 6.1‰ (refs 33-35), and ~6.7‰ for its closed-system rhyolitic differentiates (ref. 36). Therefore, we suggest that values ≥ 7‰ in the YTT quartz crystals indicate addition of an external high δ 18 O component.
Turning our attention to the single crystal data available prior to this study, there is evidence for intra-crystal isotopic heterogeneity from a reconnaissance investigation of zircon from the YTT (beyond 1σ ; maximum ∆ 18 O core−rim = 0.72‰; Fig. 2 and ref. 9) but more data are required to fully resolve the issue of intra-and inter-crystal variability at Toba. Here, we take the opportunity to better resolve isotopic diversity within and between single crystal grains of our extensive sample suite, with the added advantage that our quartz crystals are larger than the zircon grains of the reconnaissance study 9 , thus permitting more detailed δ 18 O profiling through the crystals. Further, by adopting the SIMS approach, we overcome the blending effect that is usually associated with analysis of bulk rock or bulk crystals (cf. ref. 32). Indeed, individual YTT quartz crystals demonstrate core to rim variation of up to 1.8‰, and, importantly, ~40% of the analysed quartz crystals show a significant decrease (beyond 1σ ) in their δ 18 O values from their cores to their outer growth zones (Fig. 4). The remaining ~60% of the analysed quartz crystals display limited resolvable core to rim variation in δ 18 O values along the analysis profiles. The possibility to conduct a large number of intra-crystal analyses in single quartz grains therefore allows us to record isotopic changes over small length-and hence time-scales, and so counters the averaging effect of whole-grain crystal data or spatially-limited intra-crystal traverses, as is often carried out for co-magmatic zircon (cf. refs 37 and 38) (Fig. 2).
The δ 18 O magma values of the YTT previously recorded in zircon, and now also in quartz, are consistent with an open and dynamic magma reservoir that received variable crustal and magmatic inputs during its evolution (Fig. 5a). To explore the early quartz crystallisation history at Toba, we modeled how much assimilation of granitic crust with an average δ 18 O value of 10‰ would be required to achieve our average magma value of 8.3‰, equivalent to the average measured quartz core value of 8.8‰. We employ the following mass-balance mixing model: is the value of the starting magma, assumed here to have a δ 18 O magma value of 7‰, δ O c 18 is the oxygen isotope ratio of the new component that enters the system, set here as 10‰ (i.e. a whole-rock value in equilibrium with our measured crustal quartz of 9.5‰), and x is the mass fraction of this component. Assuming a rhyolite magma with a δ 18 O starting value of 7‰ and a crustal contaminant with 10‰, then approximately 43% assimilation of high δ 18 O granitic crustal material would be required to produce the magma δ 18 O value in equilibrium with our average quartz core value. However, the exact amount of assimilation is difficult to constrain due to the wide variation in crustal δ 18 O values, and because it is likely that the highest quartz core values may represent xenocrystals of crustal origin, found at Toba and elsewhere in the Sunda arc (e.g. refs 24 and 39). Our new Sr-Nd-He isotopic data corroborate the oxygen isotope data and provide further evidence for addition of significant amounts of continental crust to the Toba magmatic sys-   to intra-reservoir convection and remobilisation processes 42 . The quartz crystals in the YTT have thus inherited a range of δ 18 O values, which when averaged produce a similar value overall. While we cannot always be sure that we have sampled the very outer rims of a crystal, the fact that intra-crystal changes tend to progress to lower values in some 40% of crystals analysed would imply that this incoming magma was recorded only in certain parts of the gigantic YTT magma reservoir, and thus effectively during the closing stages of pre-eruptive YTT chamber evolution.
Such a change towards relatively low δ 18 O values is also consistent with the concept of late-stage compositional diversification proposed previously for Toba 12,38,43 , and similar high-silica systems, such as Yellowstone or the Snake River Plain 7,9 . It is important to note that such a decrease in δ 18 O values cannot result from normal closed-system crystal fractionation. In addition, due to the preserved primary textures and intra-crystal compositional variations, we can exclude any post-crystallisation fluid overprint, which would result in irregular domains that reflect fluid fronts. Therefore, either basaltic replenishment with mantle-like δ 18 O values could have occurred 12 or low-δ 18 O hydrothermally altered materials present in the shallow crust were added to the crystallising magma during the final stages of magmatic evolution, perhaps in a fashion similar to what was recently suggested for Yellowstone and for the Snake River Plain rhyolites 7,9,44 . In order to quantify the addition of a low-δ 18   mass fraction of this component. To provide a first order assessment, we used a relatively normal mantle-type δ 18 O value of 5.5‰ for basalt replenishment 35 , and we consider hydrothermally-altered volcanic material from the caldera roof to have δ 18 O = 0‰, as has been suggested for comparable high silica systems elsewhere 7,44 .
Our modelling suggests that replenishment by basaltic magma with δ 18 O = 5.5‰ would be required in a proportion of ~75% basalt to 25% rhyolite (at δ 18 O = 8.3‰) to shift host magma to a δ 18 O magma composition identical to values calculated from the lowest quartz rim values of our study (lowest quartz rim δ 18 O magma = 6.2‰; Supplementary Figure 1). However, this proportion of mafic magma is unrealistic because quartz would no longer crystallise from a mixed magma of such an intermediate composition since the SiO 2 content would be far too low to form quartz 45 . Indeed, a melt SiO 2 content exceeding 70 wt.% is required to crystallise quartz at the pressure of the YTT magma chamber 23,45 . Furthermore, mafic replenishment would be expected to produce low-δ 18 O rims grown directly onto resorption surfaces due to contact with mafic magma 12 and this is not observed. More generally, the role of mafic replenishment in driving super-eruptions is debated and indeed they may not affect super-eruption type magma reservoirs significantly because incremental mafic replenishments would not usually over-pressurise or pervasively alter the composition of such a large reservoir (cf. refs 46 and 47). However, mafic replenishments can act as a heat source if there is underplating or injection into the lowermost parts of a silicic magma reservoir 42 . Crystal zones that formed following underplating could reflect a magma hotter than recorded by previous rims. In other words, the parental magma batches would have been heated to temperatures above zircon and quartz saturation at a given δ 18 O value, such that this hotter magma would have promoted assimilation of roof material. Subsequent cooling below zircon and quartz saturation would crystallise these minerals with diverse δ 18 O values, which would explain, at least in part, the degree of the isotopic diversity observed in YTT zircon and quartz. High-temperature hydrothermally-altered roof and wall material surrounding the magma reservoir is an therefore alternative low-δ 18 O component that could have been added, especially given the observation that widespread low-δ 18 O crustal country rocks are probably not abundant beyond larger volcanic centres in this part of the globe 36 . Hydrothermally-altered roof rocks can have δ 18 O values of ~0‰ or lower due to prolonged high-temperature interaction with meteoric waters 7,9,44,48 . In the case of Toba, the reservoir roof and wall rocks are either composed of silica-rich magmatic rocks or pre-existing granitoids and are probably heavily fractured due to previous caldera collapses 15 . These rock types are fusible over short timescales, especially when hydrated and hydrothermally overprinted 43,[48][49][50] . Our model suggests that admixing ~25% of a rhyolitic melt derived from a high-temperature altered roof rock with a δ 18 O value of 0‰ can reproduce the lowest δ 18 O rim values of the analysed YTT quartz (Fig. 5b). This amount of assimilation is energetically feasible and in line with estimates from other large caldera systems (e.g., ref. 44). For these reasons, we favour pre-eruptive addition of low-δ 18 O hydrothermally altered crustal melts to explain our new SIMS oxygen data.
A corollary of the roof assimilation model is the initiation of a positive feedback loop wherein explosive shattering of H 2 O-rich silicic rock on contact with magma could accelerate assimilation 51 . Moreover, after several cycles of caldera-forming activity at Toba, the reservoir roof (i.e. the caldera floor rocks) would be highly fractured and liable to mechanical break-up by piecemeal intra-caldera structures on a variety of scales (e.g., refs 15 and 52). Such a faulted and fractured roof would provide pathways for infiltrating hydrothermal fluids, and would provide an easily mobilised reservoir of low-δ 18 O material spatially adjacent to the magma reservoir. Indeed, mechanical roof failure and roof subsidence into the magma reservoir is a highly probable eruption trigger at large caldera systems such as Toba 52,53 , especially as this caldera sits upon the Batak Tumor, an enormous swell or tumescence area around Lake Toba, which is related to emplacement of the large Toba reservoir and probably also to buoyant magma rise 15 (Fig. 1). Associated roof uplift from emplacement and rise of silicic magma within this swell will have imparted an extensional regime onto the reservoir roof (e.g., refs 15, 46, 47 and 52), which would thus have been intensely fractured and weakened already prior to the YTT event and hence made susceptible to partial collapse and rapid recycling during caldera unrest.

Linking assimilation with timescales
Seven of the 15 YTT quartz crystals examined here for their δ 18 O values were previously analysed for Ti-in-quartz diffusion to derive residence timescales 12 . Ti-in-quartz showed that the outer YTT quartz rims on the studied crystals had an average residence time on the order of 10 2 years, which most likely reflects the period of final crystal growth prior to eruption. On the other hand, chemical diversification already commenced some 35 ka before the final YTT event 38 , and likely involved mafic underplating and small to medium-sized replenishments (cf. ref. 54). However, our new δ 18 O data lead us to suggest that chemical diversification of the YTT magma may have been influenced by crustal assimilation and progressively more pronounced self-assimilation that recycled part of the volcanic superstructure due to increasing roof instability in the run-up to the cataclysmic YTT event. Numerical experiments have already highlighted that hydrated volcanic material can melt over very short timescales even in large silicic caldera settings and at rates of many metres per year 43 . Such timescales are in agreement with those established for the YTT system 12,38 .
Late-stage addition of high-silica, hydrated, and hence fertile crustal material from the reservoir roof, as supported by our δ 18 O quartz data, would likely supply extra volatiles to the resident magma and so could promote build-up of sufficient volatile over-pressure to initiate an eruption (e.g., refs 55-57; Supplementary Figure 2). The source of these extra volatiles could be (i) water-filled vesicles in assimilated roof material (e.g. ref. 58), or (ii) assimilation of hydrothermally altered clay-rich intra-caldera materials with high H 2 O contents (e.g. > 6 wt.%, ref. 59). Furthermore, stoping may have been an important mechanism to swiftly introduce hydrous roof material, leading to initial degassing from clays and volatile-filled vesicles of the stoped material, followed by partial hydrous melting, and, if left for sufficient time at high temperature, bulk melting of the solid residue of the stoped materials. Due to its high silica nature, and upper crustal storage level, we envisage that the Toba system would have rapidly achieved volatile oversaturation upon assimilation of hydrous roof material, which may explain why only 40% of our quartz crystals record a late-stage drop in their δ 18 O values (e.g., Fig. 6). A large portion of the The wider implications of our results are that the large size of quartz crystals and its high oxygen diffusion rates at magmatic temperatures (10 4 times faster than zircon 10,11 ) now allow us to resolve events on the order of hundreds of years prior to eruption and to assess the degree of late-stage crustal assimilation. As illustrated in Fig. 6, the spatial resolution of SIMS analysis of quartz can provide a level of detail that complements whole-rock and zircon studies (cf. refs 9 and 38), and, in the case of Toba, records an otherwise undetected change in the δ 18 O values of the outer zones of a portion of the studied quartz crystals.

Methods Summary
The YTT whole rock and pumice samples were crushed, sieved, cleaned and picked for pristine quartz, feldspar, biotite and amphibole. For CL characterisation, quartz grains were mounted in epoxy resin and analysed with a CL detector attached to a JEOL microprobe at the Department of Earth Sciences at Bristol University. Whole rock powders were analysed for oxygen isotopes using a conventional line and mineral grains by using a laser fluorination line at the University of Cape Town. For SIMS analysis, quartz crystals were fixed in epoxy mounts and coated with gold. In-situ O-isotope measurements were made using a CAMECA IMS 1280 SIMS instrument at the Swedish Museum of Natural History, Stockholm (NordSIM). A 20 keV Cs + primary beam of ca. 2.5 nA was used in critically-focussed mode together with a 5 μ m raster to sputter a ca. 10 μ m sample area. The runs comprised a 90 second pre-sputter period with a raster of 20 μ m, and field aperture centering using the 16 O signal followed by 64 seconds of data acquisition using two Faraday detectors in the multicollector system operating a common mass resolution of ca. 2500. Corrections for instrumental mass fractionation were determined using international reference material NBS-28 (silica sand), with further details given in refs 60 and 61. Full descriptions of all methods employed, including for Sr-, Nd-, and He-isotope analyses, are provided in the Supplementary Information.