Resurgence initiation and subsolidus eruption of cold carapace of warm magma at Toba Caldera, Sumatra

Supervolcanoes like Toba Caldera, Sumatra, produce the largest eruptions on Earth. However, the magmatic conditions and processes during the period of recovery after catastrophic supereruptions, known as resurgence, are poorly understood. Here we use Bayesian statistical analysis and inverse thermal history modelling of feldspar argon-argon and zircon uranium-thorium/helium ages to investigate resurgence after the 74-thousand-year-old Youngest Toba Tuff eruption. We identify a discordance of up to around 13.6 thousand years between older feldspar and younger zircon ages. Our modelling suggests cold storage of feldspar antecrysts prior to eruption for a maximum duration of around 5 and 13 thousand years at between 280 °C and 500 °C. We propose that the solidified carapace of remnant magma after the Youngest Toba Tuff eruption erupted in a subsolidus state, without being thermally remobilized or rejuvenated. Our study indicates that resurgent uplift and volcanism initiated approximately 5 thousand years after the climactic caldera forming supereruption. Post-caldera eruptions at Toba Caldera, Sumatra, following the supereruption 74,000 years ago, sampled cooler margins of the warm magma reservoir and imply its heterogeneity, according to thermochronology analyses coupled with Bayesian statistics.

L arge calderas are collapse structures formed by the catastrophic evacuation of magma from the tops of upper crustal silicic magma reservoirs. Often called supervolcanoes, these calderas and their associated magmatic foundations provide valuable insight into outstanding questions about the nature of crustal magmatic systems, the construction of upper crustal batholiths, and the evolution of the continental crust. They are also the source of some of the most catastrophic geologic processes on Earth-supereruptions, and as such are considered among the most hazardous volcanoes on Earth. The period of recovery after a supereruption, known as resurgence, remains one of the most poorly understood volcanic phenomena, despite the fact that all large active calderas on Earth today like Yellowstone (USA), Campi Flegrei (Italy), and Toba (Indonesia) are resurgent (long-term activity) and restless (short-term activity), presenting evidence of continued unrest in the form of gas emissions, seismicity, structural adjustments, and eruptions 1,2 . As the volcanological community grapples with hazard assessment and monitoring of current activity with a view to forecasting future activity and mitigating hazards, a better understanding of the resurgent and restless history of large calderas is, therefore, an imperative.
Forecasting future activity and hazards from active calderas, and indeed from volcanoes in general, requires knowledge of the pre-eruptive magma storage. However, fueled partly by the lack of geophysical evidence for "eruptible magma" beneath active calderas and other volcanoes [3][4][5] , our discourse about the preeruptive state of magmatic systems is polarized. A view that espouses a "cold" non-eruptible, possibly subsolidus state with only a brief excursion to high temperatures required to mobilize and erupt [6][7][8] is opposed by those that find evidence for prolonged and fitful pre-eruptive histories at temperatures well above the solidus [9][10][11] . Like the famous Hindu proverb of the blind men and the elephant, these disparate views expose that our interpretations may be limited by the parts of the magmatic system that is sampled by the eruption and the tools being used 12,13 . Given that magmatic systems being compared range at least three or four orders of magnitude in size and are proportionally thermally and rheologically heterogeneous, the interpretations of crystal-scale petrochronology are likely to be non-unique [14][15][16] .
Herein we present a thermochronologic perspective on these issues from Toba Caldera on the island of Sumatra, the site of the Earth's largest recent supereruption. First, we resolve the ages of post-caldera eruptions that constrain the initiation of resurgence and recovery from the climactic supereruption of 74 ka. Then we show that these post-caldera eruptions sampled the cold "halo" of a long-lived warm magma system and erupted under subsolidus conditions without thermal remobilization. This work demonstrates that magmatic reservoirs are thermally heterogeneous, open, and incompletely sampled by eruptions, resulting in a spectrum of complementary interpretations about their preeruptive state. Different thermal histories from different parts of a magma reservoir should be expected and characterizing preeruptive histories as warm-or cold-stored is not useful 17 .
Toba Caldera. Toba Caldera is a composite caldera complex resulting from four caldera-forming eruptions that punctuate a 1.2 Myr cyclic history 18 . The most recent caldera cycle culminated with the eruption of the Youngest Toba Tuff (YTT) ca. 74 kyr ago based on sanidine 40 Ar/ 39 Ar (hereafter Ar) ages of 73.9 ± 0.6 19 and 75.0 ± 1.8 ka 20 (all uncertainties herein are reported at 2σ level). This eruption evacuated at least 2,800 km 3 of rhyolite magma from a "warm" reservoir several times that volume, where zircon was quasi-continuously saturated for several hundred thousand years maintained by fitful magma recharge 11 .
Since the YTT eruption, Toba Caldera has been in resurgence, i.e., the recovery phase after the supereruption where magmastatic, lithostatic, and isostatic equilibrium of the caldera superstructure is re-established after the catastrophe of the climactic caldera-forming eruption. At large calderas like Toba, this is commonly manifested by the uplift of the caldera floor as the magma rebounds (i.e., viscous rebound 21 ) and new magma is intruded. Eruptions often accompany structural uplift as extension creates pathways to the surface. Samosir Island, the structural uplift associated with the YTT stage of the Toba Caldera, has been found to have uplifted over a period of at least 36−42 kyr 22 and has been modelled to have been driven by the magmatic force of rebounding remnant magma 23 . The structural resurgence was accompanied by volcanic activity (Fig. 1a). On Samosir Island, lava domes erupted on the north and east coast, forming the rhyolitic Samosir lava domes [24][25][26] . These domes are generally geochemically indistinguishable from the YTT, although much more crystal-rich (YTT: 12-25 wt.%; domes: 35-50 wt.%), and have been proposed to be extrusions of remnant YTT magma 23,24,27 . At the southern end of the caldera, post-YTT lava domes are collectively known as the Pardepur lava domes. These are distinct from the Samosir domes as they are dacitic rather than rhyolitic, and do not have sanidine as a phenocryst phase. Other post-YTT eruptives include basaltic andesite to dacitic domes and lava flows in the north and along the west, forming Sipisopiso, Singgalang, and Pusuk Buhit, respectively 24,25 .
Recent work has shown that the Samosir domes are concordant in Ar age with the climactic eruption~74 ka, and this has been interpreted to indicate that the domes are pre-resurgent and may have been in place thousands of years before resurgence uplifted them to their current locations 26 . However, this interpretation is at odds with evidence from combined U−Th-disequilibrium and (U−Th)/He zircon (hereafter ZHe 28 ) geochronology in Mucek et al. 25 who resolved younger eruption ages for some of these post-caldera lava domes. This new chronology has important implications for the recent and current state of the Toba system, but the apparent discordance between the Ar and ZHe ages has not been addressed. This paper presents new Ar data, recalculated ZHe ages based on new U−Pb zircon data, and thermal history models to investigate this discordance and its implications for the resurgent history of Toba Caldera and pre-eruptive magma dynamics there.

Results and discussion
Recalculated ZHe eruption ages are presented in Table 1 and Figs. 1, 2. The ZHe eruption ages for two YTT samples recalculated based on newly acquired U−Pb zircon crystallization ages (Supplementary Data 1) are 75.7 ± 4.1 and 75.5 ± 4.7 ka. These are concordant with the published Ar ages of 73.9 ± 0.6 and 75.0 ± 1.8 ka on sanidine from YTT 19,20 , as well as with our new Ar ages (74.2 ± 0.4 and 75 ± 0.5 ka) reported in Supplementary Data 2. These new ZHe ages are also somewhat younger and more accurate than the ZHe ages of 79.5 ± 5.5 and 78.1 ± 5.8 ka reported for these samples by Mucek et al. 25 , which were corrected for disequilibrium based on U−Th-disequilibrium data only 25 . This suggests that the ZHe data reduction procedure employed in this study offers more accurate and more reliable results. Recalculated ZHe ages for Samosir (North), Samosir (Tuk Tuk), and Pardepur (South) post-YTT lava domes are 70.0 ± 3.2, 67.1 ± 2.8, and 61.8 ± 4.1 ka, respectively ( Fig. 1; Table 1; Supplementary Data 3).
Individual sanidine crystals from the two YTT samples returned concordant stacked plateau Ar ages of 74.2 ± 0.4 and 75.0 ± 0.5 ka (Supplementary Fig. 1a(i), (ii)). The corresponding stacked inverse isochron ages of 74.5 ± 0.4 and 75.2 ± 0.5 ka, respectively, concur with the stacked plateau ages (Supplementary Fig. 1c(i), (ii)). All our ages are concordant with published ages; our new YTT eruption ages are determined by sanidine Ar 19,20 and ZHe methods 25,29 . Bayesian inference results confirm that there is no credible difference between posterior mean distributions of our sanidine Ar and ZHe ages for the YTT ( Fig. 3; Table 1).
Thermochronologic constraints on magma dynamics. With respect to ZHe ages, our new Ar ages are concordant for the YTT samples, and older for the two Samosir and the Pardepur (South) lava domes (Fig. 3). The indistinguishable ZHe and Ar ages yielded by the two YTT samples imply rapid cooling through the temperature sensitivity ranges of the Ar system in feldspar (T c = 350°C; Ar partial retention zone (PRZ):~290-450°C 8,31 ) and He system in zircon (T c = 190°C; He PRZ:~160-260°C 32 ), which is confirmed by thermal history modelling results (Fig. 4a). This rapid cooling is consistent with catastrophic eruption. In contrast, the age discordance in the post-caldera domes implies a more complex thermal history.
The discordance with ZHe ages found in the lava dome samples insinuates that sanidine and plagioclase in the lava domes retained their YTT Ar signatures for several thousand years. Our hypothesis is that the post-caldera domes with discordant ages tapped the "cold halo" of the climactic 74 ka Toba magma system, where Ar retention was possible (Fig. 5).
The survival of YTT ages in feldspars from the post-caldera domes indicates that their pre-eruptive storage was under conditions where Ar was retained, indicating pre-eruptive storage temperatures of less than~350°C. Helium in zircon, in contrast, became immobile only after these lava domes erupted and cooled below~190°C. ZHe ages thus more reliably record the time of eruptive quenching than Ar ages. The interpretation of ZHe ages as cooling ages recording eruption ages is corroborated by (i) absence of evidence for any post-dome-eruption thermal events in the area and by (ii) a lack of correlation between ZHe ages and crystal size parameters ( Table 1) that would have caused partial rejuvenation of ZHe ages or indicated slow cooling. Based on the Bayesian analysis the ZHe eruption ages are significantly younger than feldspar Ar ages by 4.6 ± 2.4 to 13.6 ± 4.6 ka (Fig. 3). We interpret these Δt's as the duration of pre-eruptive storage between~350 and~190°C-the respective Ar and He closure temperatures.
To constrain time−temperature conditions of pre-eruptive crystal storage at higher resolution, thermal histories of the samples were modelled based on the measured feldspar Ar and ZHe thermochronology data using different modelling strategies designed to answer specific questions (see Supplementary Note 3 for more details). The thermal history modelling results suggest the following key points: I. Ar and ZHe data from dome samples require a retardation in cooling rate after the YTT event (Fig. 4a). In other words, these data cannot be explained by a simple, fast cooling history as in the case of YTT samples; II. Minimum and maximum temperatures the dome magmas could have experienced during post-YTT pre-eruptive storage and eventual eruption are in the range from >280°C (T min ) to <500°C (T max ) (Fig. 4b); III. Maximum duration of pre-eruptive magma storage, or, in other words, the maximum periods the dome magmas could have been stored at constant temperatures prior to eruption are constrained to 4.5 kyr and 7.4 kyr for Samosir domes, and to 12.7 kyr for Pardepur dome. Temperatures that allow the maximum duration of pre-eruptive storage period to be achieved, are in the range of~300−425°C (Fig. 4c).
The domes are eruptions of subsolidus remnant magma. At the calculated storage temperatures, the dome magma would have been subsolidus and uneruptible. Deeper, hotter magma might be invoked to rejuvenate "uneruptible magma" [33][34][35] and this would Mucek et al. 25 and single-crystal Ar data produced by incremental heating or total fusion from this study, colour-coded according to the Fig. 1. Coloured vertical bars correspond to 2 sigma uncertainties for individual analyses; grey analyses are not included in the weighted mean/plateau age calculation for the reasons given in Supplementary Data 2 and 3; translucent analyses are out of scale (full data are listed and shown in Supplementary Data 2). The thin grey horizontal lines through each population represent the weighted mean (for ZHe data) or plateau age (for Ar data), the outer thick grey horizontal lines mark the corresponding 2 sigma uncertainty. Labels: sample code, type of analysis (ZHe or Ar), age and 2 sigma uncertainty in ka, number of analysesconsidered for age calculation and total, respectively. Purple rectangle across the length of the figure represents the currently accepted YTT eruption ages based on sanidine Ar data of Mark et al. 20 and Storey et al. 19 . Note that (i) the weighted mean ZHe ages and Ar plateau ages for YTT samples overlap within uncertainty with the accepted YTT eruption age, which proves the accuracy of the techniques used. (ii) The weighted mean ZHe ages (interpreted to represent the time of eruption) for all post-YTT domes are younger (outside the uncertainties) than the corresponding Ar ages as well as the accepted YTT eruption age. This suggests that the post-YTT domes erupted significantly later than the closure of the Ar system and after the YTT event.
be consistent with the recharge-fueled thermal longevity of the Toba magma system 11 . However, our thermal modelling indicates that if this was the case at Toba, then the discordance in ages allows maximum temperatures of only <500°C (Fig. 4b, c and Supplementary Fig. 3). Under these constraints, we propose that the domes we have sampled were erupted in a subsolidus state and were not thermally rejuvenated. It is increasingly recognized that the largest catastrophic caldera-forming eruptions are driven by the foundering of the roof into the chamber [36][37][38][39][40] . Such eruptions arrest, once the foundering caldera block meets increasingly dense and viscous crystal-rich magma and either isostatic equilibrium, is reached between the subsiding block and the remnant magma or a maximum viscosity condition is reached 41,42 . Thus, we propose that after the YTT eruption the remnant magma was a coarsely crystalline mush that was barely mobile and non-eruptible. The upper margins of this magma intruded into fractures or eruptive conduits at the caldera margin and caldera floor where it would be rapidly chilled as dikes and plugs with interstitial quenched glass (Fig. 5a). This chilled zone would have enhanced the rheological barrier at the top of the chamber 43 . Within fourthousand years, the smallest Δt we have calculated (Fig. 4c), resurgent uplift of the caldera block was initiated by magmatic rebound and recharge. This would reopen preexisting weaknesses in the caldera floor and solidified "cold" remnant YTT magma was squeezed out as deeper mobile remnant YTT and hotter deeper recharge magma penetrated into the opening fractures and pushed the subsolidus magma out to form the effused domes. Reheating may have occurred, but was either not above 500°C or not pervasive or steady enough to reset the Ar clock in sanidine. Some partially or completely molten magma may have effused to endogenously grow the domes late in the process, but this would have been deep in the domes and those parts would not be sampled at the surface.
Effusion of solid dome-forming rock as "whalebacks" or spines lubricated or pushed by deeper hotter magma has been described at many recent and historic dome-forming eruptions [44][45][46][47][48]  In proposing this hypothesis for the Samosir and Pardepur (South) dome, we emphasize that the volume of individual dome masses on Samosir are tiny (<0.1 km 3 ) and the aggregated volume of magma represented by the Samosir domes is <1 km 3 ,~0.04% the volume of the 2,800 km 3 YTT eruption or 0.16% of the 600 km 3 of resurging magma calculated to be the motivating drive for the Samosir uplift 22 .
Timing of dome effusion and resurgent uplift at Toba. The eruption age of domes along faults associated with resurgent uplift provides a constraint on the initiation of resurgence. Our eruption ages of 70.0 ± 3.2 and 67.1 ± 2.8 ka for the Samosir (Tuk Tuk) and (North) domes, respectively, in context of the sedimentary and paleomagnetic evidence that resurgent uplift progressed for at least 40,000 years after the climactic eruption 22,23,25 and the evidence for a complex tectonic history until very recentlỹ 2 ka 23 , indicate that these domes are syn-resurgent and postdated the YTT eruption by 4.6 ± 2.4 kyr (Samosir (North)) and 8.0 ± 2.1 kyr (Samosir (Tuk Tuk)). These delays are also consistent with the recognition of lake sediment beneath the Samosir (Tuk Tuk) group of domes 25 .
The~13 kyr delay until the eruption of the Pardepur (South) dome is also consistent with the timeframe of resurgent uplift. Although Pardepur (South) is not part of the resurgent uplift of Samosir Island, it is located on a NW-SE fault that defines the western edge of Samosir on which other post-caldera effusions at Pusuk Buhit and a series of uplifted areas interpreted as cryptodomes are also located 24,49 . The clear spatial relationship between post-climactic eruptions and faults involved in the deformation of the resurgent Samosir Island further supports a magmatic drive for resurgent uplift within the caldera.
'Cold' versus 'warm' storage: different parts of the elephant? Our new data insinuates that the uppermost, cold halo of the (Note that these histograms are parameter values from the posterior distributions that are permissible or "credible" within the measured data. The histograms are not the measured data distributions.) Horizontal bars represent 95% high-density intervals (HDI); for exact values see Table 1. Since the 95% HDI of the difference of means falls above zero and 100% of the predicted credible values are >0, we can conclude that the means of Ar and ZHe datasets are credibly different (full data in Supplementary Fig. 2 magma reservoir was tapped and squeezed out along faults on the edges and within the uplifting resurgent dome to form the two Samosir domes, Pardepur (South), and maybe other post-caldera domes on Samosir (Fig. 5). Most of the mass of the domes would therefore have been in 'cold' storage and effused in the solid-state to form the Samosir and Pardepur (South) lava domes. These three domes, therefore, record the thermal history of the cold zones of the otherwise long-lived warm Toba magma 11 and represent a serendipitous archive of contrasting thermal histories in the Toba magma. Such contrasting thermal histories should be an expected condition of large magma systems with dominant molten magma being surrounded by a progressively cooler halo that would include subsolidus portions. Sampling of this coldstored halo has been described from Long Valley caldera where, on the basis of high precision single-crystal Ar dating, Andersen et al. 8 argued that sanidine antecrysts in the Bishop Tuff were derived from the cold margins of the warm pre-eruptive magma reservoir. In large eruptions like the Bishop Tuff and the Youngest Toba Tuff, crystals with "cold" and "warm" histories can be mixed in the same eruption, and the relatively small volume of "cold" crystals would only be identified by serendipity or through extensive sampling approaches coupled to thermochronology. In small eruptions like the Toba domes, these antecrysts escape syn-eruptive integration and only record the "cold- stored" part of the magma reservoir. They, therefore, provide a serendipitous but vital record that without the context of the YTT eruption, might lead to the conclusion that Toba was a cold-storage magma system. A comprehensive understanding of the complete magmatic history and well-contextualized sampling is therefore crucial if pre-eruptive storage is to be fully understood. Like the famous Hindu proverb alluded to above, sampling different parts of the elephant can lead to different conclusions. Thermochronology can provide a new and valuable direct thermal perspective to studies of pre-eruptive magma residence. The insights can complement and extend other detailed petrochronologic approaches that have been used to identify complex pre-eruptive magmatic histories at other young volcanic systems that might include rejuvenation (e.g., Taupo volcano 50 ).

Conclusion
Bayesian inference reveals that ZHe ages from three post-caldera domes at Toba Caldera are up to ca. 14 kyr younger than corresponding Ar ages from feldspars that therefore must be YTT antecrysts. Coupled diffusion models constrained by this discordance support a thermal history where post-YTT eruptions tapped the "cold halo" of the warm remnant post-climactic magma. Up to ca. 13 kyr of "cold" storage at temperatures between~300 and~425°C are permitted; minimum and maximum storage temperatures are 280 and 500°C, respectively. This thermal history is consistent with dome eruptions under subsolidus conditions motivated by magmatic pressure of the resurging remnant magma. No evidence for pervasive thermal remobilization or "rejuvenation" is recorded. Our work thus demonstrates a significant delay between the YTT eruption and eruption of these domes. We propose that eruptions of the domes signal the onset of resurgent uplift and associated opening of pathways to the surface through which remnant solidified conduit plugs and dikes were extruded to the surface by invading magma acting like the plunger in a syringe. These insights have only been possible by leveraging the unique ability of coupled thermochronology and the resolving power of Bayesian inference. Our study promotes cold and warm storage as part of a continuum of pre-eruptive thermal histories that should be expected in magmatic reservoirs and cautions that the answers to the questions we pose about magma storage will depend on what we sample and study. 40 Ar/ 39 Ar dating. Samples of three post-caldera lava domes previously dated by Mucek et al. 25 were dated by the 40 Ar/ 39 Ar (hereafter Ar) method applied to sanidine and plagioclase; two from the Samosir lava domes and one from the southern Pardepur lava dome (Fig. 1a). Additionally, two YTT samples, one from the north and the other from the south, were also processed to act as a procedural reference (Supplementary Data 2). Ar dates (Supplementary Fig. 1) were determined by incremental heating or total fusion of individual crystals on the ARGUS-VI multicollector mass spectrometer at Oregon State University Argon laboratory following the procedures detailed in Supplementary Note 1 and Supplementary Data 2.

Methods
U−Pb dating and (U-Th)/He data recalculation. The zircon (U−Th)/He (hereafter ZHe) dataset published in Mucek et al. 25 , which contains zircon crystallization ages based on U−Th-disequilibrium data, was complemented by additional U−Pb data obtained on crystals in secular equilibrium in order to better quantify the correction for U-series disequilibrium and pre-eruptive crystal residence time that needs to be applied to (U−Th)/He data 51 . Zircon from three samples was analyzed by U−Pb methods-LT12-001 and LT12-012-both from Youngest Toba Tuff and LT12-003 from Samosir lava dome. The U−Pb analysis was conducted on a CAMECA ims 1270 ion probe at UCLA and followed the procedures described in Supplementary Note 2. Previously published U−Th-disequilibrium and newly obtained U−Pb zircon crystallization data (Supplementary Data 1) were then used to correct (U−Th)/He data for disequilibrium using MCHeCalc software 52 and adopting D 230 and D 231 parameters used by Mucek et al. 25 . Disequilibrium corrected (U−Th)/He dates were then used to calculate the mean value, which is interpreted as the representative ZHe age of each sample. This was done by applying two different approaches. In the first, traditionally-used 'frequentist approach', the error-weighted mean with 95% confidence intervals were calculated using Isoplot v.4.15 Excel add-in 53 from the data after the rejection of statistical outliers identified based on the modified 2σ criterion function implemented in Isoplot 53 . In the second 'Bayesian approach' 54 , a Bayesian probability model (specified in Supplementary Note 3) was applied to predict posterior mean values and their 95% high-density intervals (HDI).
Bayesian inference for comparing two groups. The critical issue in this study is the discordance between the Ar and ZHe ages obtained on the same sample. To test whether there is a statistically significant difference between the Ar and ZHe ages obtained on the same sample and to quantify this difference, we developed a Bayesian analysis for comparing two groups which is described in Supplementary  Notes 3 and 4. The great potential of the Bayesian analysis to improve the resolution of geochronological datasets hampered by data precision has recently been demonstrated by Morgan et al. 55 , who, by employing Bayesian inference, were able to discriminate time intervals as small as~20 kyr in the eruptive history of the Fish Canyon Tuff (ca. 28 Ma) and associated units of the La Garita magmatic system 55 .
In brief, our Bayesian model is similar to the Bayesian t-test of Kruschke 30 , but was modified so that it considers uncertainties of the measured data and weights each measurement in proportion to its uncertainty following the approach of Gelman et al. 56 applied to the "eight schools" example. The model combines the data (i.e., Ar and disequilibrium corrected ZHe ages with 2σ uncertainties) and prior information to construct posterior distribution of the mean ages (and 95% credible intervals) for Ar and ZHe data. These are then used to calculate the posterior distribution of the difference between the mean Ar and ZHe ages, which can be used to answer questions such as, for example, what is the probability that Ar age is greater than ZHe age, and what is the plausible range of values of the difference? For computing the Bayesian analysis in this study, we report here a code written in R programming language 57 that uses JAGS MCMC library accessible from R via rjags package 58 . The R-code (rocks_ages_final.R) can be found in Supplementary Note 4. Given that we have no existing expectations for variance present, we assigned an uninformative prior expectation on the mean parameters to be a broad normal distribution with no negative values. Accordingly, the standard deviation of the prior on mean was set to 1,000 times the standard deviation of the pooled data. An MCMC method with 100,000 burn-in and 1,000,000 iterations was then used to compute posterior probabilities of means and difference of means for each datasets.
Thermal history modelling. Thermal histories for samples based on coupled Ar and He diffusion models were modelled using the HeFTy software 59 . The HeFTy software allows simultaneous reconstruction of time−temperature paths for multiple thermochronometric systems for which diffusion kinetics data are known. Details of the HeFTy modelling principles are described elsewhere 59 . In brief, in the so-called inverse modelling approach, used in this study, HeFTy generates a large number of random thermal histories that can be forced to pass through pre-defined time-temperature constraints. For each thermal trajectory, HeFTy calculates a thermochronological age that results from the competing processes of production, retention, and diffusive loss of daughter isotopes in modelled thermochronometric systems (in our case 4 He in zircon (U−Th)/He system and 39 Ar in feldspar 40 Ar/ 39 Ar system). The calculated thermochronological ages are then compared to the measured thermochronological ages. If the agreement between the modelled and the measured values is acceptable or good based on a p-value-based hypothesis test, the thermal histories are accepted as permissible by the measured data, in other case the thermal histories are rejected 59 .
Here the objective of the modelling was to reconstruct thermal histories that can reproduce the measured Ar and ZHe data and qualitatively constrain time −temperature conditions of feldspar and zircon crystal residence prior to the final eruption. Accordingly, four different models (Supplementary Note 5; Supplementary Data 4; Supplementary Figs. 3, 4), were employed in order to: (i) test for a difference in thermal histories between the YTT sample and post-YTT lava domes (modelling approaches termed "free run/unconstrained" and "Helium eruption age" in Supplementary Note 5); (ii) constrain minimum and maximum temperatures the dome samples could have experienced during their pre-eruptive storage ("resurgent reheating" model in Supplementary Note 5); (iii) constrain maximum duration of the pre-eruptive storage and corresponding temperatures ("storage duration" model in Supplementary Note 5).

Data availability
The authors declare that all data generated or analyzed during this study and supporting the findings of this study are included in this published article and its Supplementary Information. Supplementary Data 1-4 contain all data that was used to generate both the figures in the main text (Figs. [1][2][3][4][5] and in the supplementary files ( Supplementary  Figs. 1-4), following procedures as outlined in Supplementary Notes 1-5. Supplementary Data 1-4 can also be found on Dryad (datadryad.org), using the following DOI link https://doi.org/10.5061/dryad.9p8cz8wgs.

Code availability
The R-Code developed for computing the Bayesian analysis can be found in Supplementary Note 4. The R-Code can be run in any software that can run the R programming language. A suggested free and open-source integrated development environment for R programming language used by the authors is RStudio.