Research ArticleGEOCHEMISTRY

Transient rhyolite melt extraction to produce a shallow granitic pluton

See allHide authors and affiliations

Science Advances  19 May 2021:
Vol. 7, no. 21, eabf0604
DOI: 10.1126/sciadv.abf0604


Rhyolitic melt that fuels explosive eruptions often originates in the upper crust via extraction from crystal-rich sources, implying an evolutionary link between volcanism and residual plutonism. However, the time scales over which these systems evolve are mainly understood through erupted deposits, limiting confirmation of this connection. Exhumed plutons that preserve a record of high-silica melt segregation provide a critical subvolcanic perspective on rhyolite generation, permitting comparison between time scales of long-term assembly and transient melt extraction events. Here, U-Pb zircon petrochronology and 40Ar/39Ar thermochronology constrain silicic melt segregation and residual cumulate formation in a ~7 to 6 Ma, shallow (3 to 7 km depth) Andean pluton. Thermo-petrological simulations linked to a zircon saturation model map spatiotemporal melt flux distributions. Our findings suggest that ~50 km3 of rhyolitic melt was extracted in ~130 ka, transient pluton assembly that indicates the thermal viability of advanced magma differentiation in the upper crust.


Extraction of interstitial melt from partially crystallized, upper crustal magma reservoirs is considered to be a dominant process in generating high-silica rhyolite capable of erupting explosively (1, 2). Volcanism involving rhyolite poses a substantial societal hazard (3) and can alter global climate (4) when the volumes are sufficiently large. The segregation of rhyolitic melt from magma reservoirs is a process that results in chemical fractionation and should produce complementary crystalline residues that are observable in the upper crustal plutonic record (515). Evidence of these residues often appear in the crystal cargo that has been incorporated into erupted products as aggregates (16, 17). However, identifying plutons that are unequivocally the residuum of rhyolite extraction (i.e., silicic cumulates) can be difficult, perhaps because they are susceptible to subsequent partial remelting (18), remobilization (19), and/or preservation of trapped melt (10, 12), potentially masking diagnostic geochemical and textural signatures (7, 11). Accordingly, some models for rhyolite generation, mainly inferred from plutons, favor magma differentiation in the hotter, lower crust (15 to 45+ km deep) (2025). These models suggest an ephemeral nature of melt-rich bodies in the upper crust, whereby magmas sourcing large eruptions are viable briefly only during high-flux periods (25, 26). Alternative models suggest that magma reservoirs can be incubated in the upper crust and might spend 104 to 105 yr at high crystal fractions (2730) before melt may be removed or rapidly remobilized (31). Rhyolites can also be generated by other mechanisms, for example, via remelting of crustal material [e.g., (32)]. Our focus here is on rhyolite that may be generated via extraction from crystalline mush (1, 2).

Central to these models and the production of rhyolitic magma is an understanding of the time scales of magma accumulation, upper-crustal residence, and subsequent melt extraction. However, our knowledge of these time scales is dominated by erupted volcanic products, limiting cohesive models that seek connections between the volcanic and plutonic realms. Although silicic magma reservoirs may grow incrementally over ca. 104 to 106 year time scales, they can be remobilized before eruption much more rapidly (ca. 101 to 103 years) (33). Whereas typical long-term plutonic emplacement rates are generally 10−3 to 10−4 km3/year., the accumulation of large melt-rich magma bodies capable of erupting is widely thought to require magma fluxes one to two orders of magnitude greater (21). However, recently, some silicic plutonic systems have been shown to exhibit volcanic fluxes (34), providing important constraints on the growth rates and eruption potential of fossil magma reservoirs. Missing from these pioneering studies are constraints on the time scales of melt extraction within a plutonic system that preserves a rock record of the full sequence of differentiation building up to this evolved stage from initial mafic pulses to production of rhyolite and complementary silicic cumulate residues.

The late Miocene (7.2 to 6.2 Ma) Risco Bayo–Huemul (RBH) plutonic complex (~300 km3) in the Southern Andean Cordillera of Chile (Fig. 1A) solidified at depths of ~3 to 7 km and has been exhumed to expose its uppermost 1.5 km (9, 10, 35). Our recent work has used textural, petrologic, geochemical, and mineralogical data to hypothesize that domains of the younger, more silicic Huemul pluton represent near end-members of rhyolite melt segregation [now high-silica granite, 76 weight % (wt %) SiO2] and complementary residual silicic cumulate formation (now quartz monzonite, 62 wt % SiO2; Fig. 1A) (9, 10). The older RB pluton is an amalgamation of incrementally emplaced mafic-to-intermediate (51 to 66 wt % SiO2) magma batches that thermally primed the upper crust, thereby facilitating subsequent high-silica melt extraction in the Huemul pluton (10). Here, we test the temporal and thermal feasibility of the melt extraction model for high-silica rhyolite generation by combining U-Pb petrochronology (i.e., coupled age and compositional analysis) and 40Ar/39Ar thermochronology. We use magma emplacement and thermal evolution numerical models to track zircon saturation and growth to directly constrain the radioisotopic dates within thermally viable scenarios of melt evolution. Our results constrain upper crustal (as shallow as 3.5 km) magma differentiation via melt segregation from a mush, a process that may precede many rhyolitic eruptions (1, 2, 28, 30), over ~130 ka in a pluton.

Fig. 1 Geologic map and geochronology from the RBH plutonic complex.

(A) Geologic map of the RBH complex highlighting internal compositional domains [modified from the work of Schaen et al. (9)]. Tatara–San Pedro (TSP) and Pellado (VP) volcanoes. (B) Rank order plot of radioisotopic dates from hand samples within the RBH plutonic complex. Δt represents the time between oldest and youngest zircon date; ty indicates the youngest zircon date for each hand sample. All 206Pb/238U dates have been Th-corrected. Asterisks (*) indicate wallrock hornfels that preserve 17 to 9 Ma relic Ar ages from protolith within discordant portions of age spectra. Minimum laser ablation (LA) zircon age range of RB1503B from the work of Schaen et al. (9). CA-ID-TIMS, chemical abrasion–isotope dilution thermal ionization mass spectrometry.


Petrochronologic fingerprinting of zircon populations

Zircon is a common U-rich accessory mineral that partitions elements from its host melt, which permits determination of (i) the timing of crystallization (36) and (ii) the compositional evolution of the magma from which it grew (37, 38). Sixty-nine zircons were measured by U-Pb chemical abrasion–isotope dilution thermal ionization mass spectrometry (CA-ID-TIMS; Fig. 1B and table S1). These same zircons were previously characterized texturally via cathodoluminescence imaging before laser ablation split stream (LASS) analysis by Schaen et al. (9). This LASS study focused on rim data and determined pluton-scale age relationships and compositional evolution of zircon populations for the RBH complex (9). However, given the young ages, zircon LASS data are unable to resolve changes in zircon composition over time in this system (9). Inheritance in the form of scarce (n = 3; this study) 10 to 88 Ma xenocrystic cores was also revealed by LASS (9); in both cases, these data are easily identifiable and omitted from further discussion. The 206Pb/238U ID-TIMS dates from the RBH complex span ~1 Ma, from 7.2 to 6.2 Ma (Fig. 1B). Zircon crystallization within the northeast portion of RB spanned ~240 ka, between 7.193 ± 0.014 and 6.956 ± 0.053 Ma (2σ), with ≤20 ka between domains. Zircon dates support field observations of sharp cross-cutting relationships and confirm rapid, pulsed, incremental emplacement from RB gabbro to granodiorite (10). In contrast, zircon dates from the Huemul span ~190 ka from quartz monzonite to high-silica granite, between 6.384 ± 0.022 and 6.199 ± 0.022 Ma (n = 38; Fig. 1B). Huemul outcrops do not preserve evidence that permits relative age sequencing (10). Instead, contacts between domains are gradational with the high-silica granite structurally overlying the granite (10). Quartz monzonite U-Pb dates, which range from 6.384 ± 0.22 to 6.328 ± 0.052 Ma, are resolvable from a younger population of zircon (n = 12) in the granite and high-silica granite domains with dates from 6.261 ± 0.018 to 6.199 ± 0.022 Ma. In addition, the youngest zircon in the quartz monzonite is older than the youngest zircon from the granite and both high-silica granite by 100 to 160 ka. This implies a process that preserved a younger zircon population in the granite and in high-silica granite domains but not in the quartz monzonite, although there is a close spatial association among all three domains (Fig. 1A).

Trace element analyses (TEA) by solution inductively coupled plasma mass spectrometry (ICP-MS) on the wash products of ID column chemistry permit integration of diagnostic compositional information from the same aliquot as zircon dates (37, 39). The short crystallization durations of these plutons make zircon compositions critical to illuminate the time scales of magmatic processes and to isolate the zircon that Schaen et al. (9, 10) hypothesized might record melt extraction. With decreasing zircon age, TE concentrations, including all the rare earth elements (REE), apart from Eu, follow Y and decrease in RB samples while increasing in Huemul samples (Fig. 2A). These zircon TE trends, along with a narrow range of Yb/Dy (Fig. 2B), mirror similar fractional crystallization trends in Huemul rocks (9, 10), which indicate melt evolution from a common parental magma. Flat middle to heavy REE bulk rock trends (9) imply that Yb/Dy is being controlled by zircon crystallization and not accessory titanite. Thus, the larger range in Yb/Dy for RB samples indicates late zircon saturation at high crystallinities (80 to 90%) for these mafic to intermediate rocks. Distinct fields in the RB zircon TEA (Fig. 2A) reinforce the resolvable ages (Fig. 1B) and cross-cutting field relationships (9, 10), which together imply chemically distinct RB pulses that were rapidly (≤ 20 ka between pulses) and incrementally emplaced. Eu/Eu* trends in RBH rocks and minerals (Fig. 2, C and D) are a diagnostic indicator of feldspar fractionation and silicic cumulate formation (9, 10), which created the quartz monzonite domain. The quartz monzonite (bulk, 62 wt % SiO2) and granite (bulk, 68 to 70 wt % SiO2) zircons display a restricted range (~0.3 to 0.4 in the melt) in Eu/Eu* from two rocks with very different bulk compositions (9, 10), which warrants explanation.

Fig. 2 TEA of zircons analyzed by U-Pb CA-ID-TIMS-TEA.

(A) Y [parts per million (ppm)]. (B) Yb/Dy. (C) Eu/Eu* [EuN/(√SmN*GdN)] of equilibrium melt compositions through time, calculated from the zircon TEA data using the partition coefficients of Padilla and Gualda (40). (D) Detailed view of Huemul Eu/Eu* melt compositions through time within dotted box of (A). The purple band represents the range in Eu/Eu* of melt from the bulk rock fractional crystallization model of Schaen et al. (9) between 38 and 55% crystallinity. This range is insensitive to choice of zircon partition coefficients. Green bar in (A) represents the U-Pb zircon LA age range (not composition) of RB porphyritic diorite sample RB1503B [see Schaen et al. (9)]. All uncertainties are 2σ analytical.

We have cast the zircon TEA into estimates of equilibrium melt compositions by using zircon-melt partition coefficients (40). Estimates of the melt composition that the zircon grew from can be compared to a bulk rock fractional crystallization model for Huemul (purple band in Fig. 2D) (9). Huemul bulk rock Zr models indicate that zircon will begin to saturate at 38% crystallinity (9). Segregation of interstitial melt from crystalline magma mush (i.e., a continuous framework of crystals) is most efficient within a crystallinity window (~50 to 70% crystals) after convection has ceased but before rheologic crystal lockup (1, 2, 41). The quartz monzonite and granite TEA data reflect zircon crystallization upon saturation up to the lower bound of this optimal crystallinity window, recording parental mush conditions just before where melt segregation was viable (Fig. 2D). The compositional gap between these domains and the high-silica granite zircons is a clear indicator of a short-lived process like rhyolite melt segregation (Fig. 2D) (41, 42). Whereas the high-silica granite contains zircon inherited from this parental mush, with similar Eu/Eu* to the quartz monzonite and granite (“generation 1”; Fig. 2D), the zircon having the lowest Eu/Eu* records the highly fractionated rhyolitic melts formed via extraction (“generation 2”). We used TEA to isolate these particular high-silica granite zircons (n = 15) and performed statistical bootstrapping (random sampling with replacement) of this measured U-Pb age/uncertainty population. The difference between minimum and maximum bootstrapped ages by this method are then taken as modeled mean zircon crystallization durations. The mean of the mean durations of 15,000 trials from this bootstrapped population is 132 ka, with upper and lower bounds of 204 and 75 ka, respectively, at 95% confidence. To calculate a melt flux from these time constraints requires accurate estimates of rock volumes. The volume of the high-silica granite domain is between 34 and 112 km3, depending on the assumed thickness (1.2 to 4 km) and whether the pluton is scale invariant (43). This translates to melt fluxes between 0.00017 and 0.0015 km3/yr. These magma accumulation rates are minimum estimates, as they do not account for potential volume losses associated with dike propagation or eruption of the high-silica granite (35).

Thermal viability of melt extraction

A multiphase/multiscale pluton emplacement model was used to simulate thermally viable scenarios for Huemul melt reservoir evolution and extraction (Fig. 3A). These finite-volume simulations build upon previous multiphase approaches (41), magma-tectonic interaction and accommodation models (44), and stochastic methodology (45). In particular, these simulations extend the approach to three dimensions to account for the spatial arrangement of the mapped plutons (Fig. 4, A and B). A domain of 40 km by 40 km by 60 km (40 km in horizontal directions) is used to include the regional thermal evolution. An initial simplified layered crustal structure is used with a gabbro lower crust to an intermediate upper crust composed of metavolcanics and plutons (Fig. 4, A and B). Compositional proxies for each layer are based on regional outcrops (46) and seismic studies in the nearby Laguna del Maule volcanic field (47) and use the model of Abers and Hacker (48). Conservation of energy and mass is solved as described in Materials and Methods. These simulations are constrained by RBH rock compositions and thermodynamic calculations using rhyolite-MELTS (49). The size and three-dimensional (3D) structure of the computational domain necessitated conducting proxy momentum calculations as a subgrid model to compute the separation of melt from crystals (Fig. 4D). That is, melt extraction rates were computed using the local melt properties, melt fraction, and drag relationships between the melt and crystals (Fig. 4D). This one-way coupled approach does not calculate emergent melt/deformation feedback but does permit computing hundreds of 3D scenarios that are consistent with energy conservation and thermodynamics relationships. Uncertainty in compositional heterogeneities in the crust, thicknesses of the pluton units, crystal size distributions in the partially molten reservoir, and thermal structure in the surrounding crust motivated us to vary these properties in a stochastic fashion and to report on the range of outcomes (Fig. 3A). Each simulation begins with precursor incremental emplacement of RB-sized (26 to 88 km3) magma pulses over the time scales represented by their TIMS dates (Fig. 1B). Melt with the composition of the Huemul granite domain is then injected into the upper crust at a depth of 5 km, where crystallinity, residence time, and local chemistry of all melt are tracked throughout the subgrid model (Fig. 3A) (1, 2, 41). This procedure is repeated for a wide range of initial Huemul granite volumes and number of pulses spanning the range of quartz monzonite ages (Fig. 1B), producing the curve defined by circles in Fig. 3A. Across the 3D finite element space of each simulation, we have linked these calculations to a zircon saturation model (Fig. 4E) (9), which permits estimates of the spatiotemporal distribution of extracted melt that has achieved zircon saturation (square curve in Fig. 3A). We calculate zircon saturation as a proxy for time periods of crystal growth but do not explicitly calculate zircon growth rates as limited by diffusion. For any given simulation (identical colors in Fig. 3A), the residence time for the same volume of extractable melt is a factor of 2 to 3 shorter for only the melt that is zircon saturated (Fig. 3A). This emphasizes a critical point that the time scales represented by measured zircon dates are always inherently minimum reflections of the total magma accumulation (50). For example, the simulation that produced 48 km3 of extractable melt over ~275 ka would result in a zircon dataset spanning ~20 to 200 ka, with a mean of ~100 ka (highlighted simulation; Fig. 3A). All the melt (not just extractable) in the 48 km3 simulation is present for ~300 ka, although mostly at very small melt fraction.

Fig. 3 Melt volume and time scales derived from numerical simulations of pluton emplacement.

(A) Finite element numerical simulations of pluton emplacement that track the time scales and volumes of melt that is within the optimal magma crystallinity (~50 to 70%) to be extracted. The colored bar and associated symbols show the range of emplaced intrusion sizes that were simulated for the Huemul pluton. Circles represent total extractable melt for a given simulation, while squares represent only extractable melt that was zircon saturated for that same simulation (indicated by identical color). The gray horizontal bars are the 1σ SD interval of a zircon age distribution (squares are the mean). For both curves, symbols are sized by the duration of all the melt present (not just the extractable melt). (B) Comparison between zircon crystallization durations generated from numerical simulations (green) and the mean crystallization duration from a bootstrapped population of the measured TIMS dates (pink) using the subset of high-silica granite generation 2 zircons isolated in Fig. 2D. Green distribution comes from the 48-km3 melt simulation outlined in magenta in (A).

Fig. 4 Compilation of melt segregation observations in the Huemul pluton as related to the geology, zircon petrochronology, and the numerical modeling.

(A) Macroscale simulation of RBH intrusions (location of the intrusions shown as the blue, transparent surfaces) showing isosurfaces of temperatures (200°, 400°, 600°, and 800°C). Relative location indicated by the overlain digital elevation model (DEM). This plot occurs 50 ka after the start of the Huemul intrusion and corresponds to an intruded volume of 286 km3. Frontal-arc volcanoes TSP and VP correlate to Fig. 1A. (B) Geologic cross section of the Huemul pluton. A-A′ is indicated in Fig. 1A. (C) Schematic depiction of zircon chemical populations derived from melt extraction as isolated in Fig. 2D and loosely conforming to the cross-section B-B′ in (B). (D) Example depiction part of the domain of subgrid scale model calculation using information extracted from the magma reservoir simulations. Relative velocity is computed by determining the drag and buoyancy forces. Here, the melt fraction is shown with the background color overlain by isotherms. Vectors indicate relative separation velocity (the scale indicates the relative melt-crystal velocity produced by separation of 2 cm/yr). (E) The zircon saturation model of Schaen et al. (9) is solved via this method in each subgrid cell within the numerical models [(D) for example]. Zircon saturation experiments (36) are integrated with rhyolite-MELTS models to calculate Zr concentration required for zircon saturation (black curve). Major element content of rhyolite-MELTS model was also used to calculate the M parameter (blue curve). The Zr melt curve (green curve) assumes an initial Zr concentration equal to its bulk rock and uses the bulk partition coefficients and crystallization model of Schaen et al. (9). The Zr concentration required for zircon saturation in HC1301 (yellow star) yields a temperature of ~848°C compared to ~787°C using only the bulk saturation estimate.

To constrain the temporal viability of rhyolite segregation as the formation process for the Huemul pluton, the melt durations and cooling rates from the numerical simulations are compared to the measured U-Pb time scales (Fig. 3B) and 40Ar/39Ar thermochronology (Fig. 1B), respectively. The zircon ages from the 48-km3 extracted melt numerical simulation display a comparable distribution to the mean zircon crystallization durations from the bootstrapped age population (generation 2; Fig. 4C) that were isolated using U-Pb CA-TIMS-TEA (Fig. 3B). Whereas null hypothesis testing suggests that these two distributions are not identical, their weighted means are within the average measurable uncertainty of the TIMS dates (Fig. 3B). Given the analytical limitations for these very young low-U zircons, the good agreement between the melt extraction numerical simulations and our isotopic dates from the rock record is notable. This implies from a temporal perspective that extraction of ~48 km3 of melt from granitic mush is a viable process to produce the Huemul high-silica granite.


The predominant limitation of melt segregation in the cold upper crust is the inevitability of the magma cooling to rock. Whereas incubation of crystalline magmas in a state of cold storage (27, 28) can enhance mush longevity, segregation of melt (±crystals) is a precursor to some styles of rhyolitic eruptions (51). Accordingly, we examine the end-member scenario of pure conductive cooling without recharge (52) in the numerical simulations to compare maximum cooling rates and minimum thermal time scales to the measured isotopic constraints. Using the maximum Δt of U-Pb dates (Fig. 1B) as an estimated duration between zircon saturation (850°C) and the solidus (680°C) for Huemul rocks (9) results in magmatic cooling rates of about 650° to 1400°C/Ma. These high cooling rates are consistent with crystallization of the youngest zircons coevally to amphibole and biotite 40Ar/39Ar closure dates in the same hand samples (Fig. 1B). Resolvable biotite and orthoclase 40Ar/39Ar dates (Fig. 1B) also permit subsolidus cooling rates between biotite (~320°C) and orthoclase (~280°C) closure (53, 54) in the range of 80° to 170°C/Ma. The average magmatic cooling rate of the 48-km3 extracted melt simulation is about 600°C/Ma, in good agreement with those calculated using the geochronology. This suggests that while cooling was rapid in Huemul, segregation of ~50 km3 of melt could have occurred in ≤130 ka (Fig. 3B).

Magma accumulation rates calculated using the Huemul TIMS dates are lower than those from larger plutonic bodies such as the Golden Horn batholith (34), which contain zones (>424 km3) that exhibit the high fluxes (~0.0125 km3/yr) comparable to silicic supereruptions (50). Supereruption-sized rhyolitic systems preserve higher flux rates than the 0.00017 to 0.0015 km3/yr. rates estimated here that formed the Huemul high-silica granite. Eruption rates from Oruanui (New Zealand, 530 km3) average 0.85 km3/yr (17); young Yellowstone extrusives (Wyoming, USA, ~40 to 70 km3) are bracketed between 0.0067 and 0.07 km3/yr (55), and the Bishop Tuff (>600 km3, California, USA) preserves rates of 0.0075 km3/yr (56). The upper bounds of the fluxes that we calculate for Huemul are comparable to some volcanic systems of similar volume (30, 50, 57, 58). For example, ~20 km east of the RBH complex, the rear-arc Laguna del Maule volcanic field (>40 km3 of rhyolite erupted during the last 20 ka) preserves time-averaged magma accumulation rates of 0.0023 km3/yr (58) but has also experienced high-flux events driven by recharge rates up to 0.03 to 0.04 km3/yr (58). Although it remains uncertain whether Huemul produced an eruption, the accumulation rates here suggest that this scenario is unfavorable. Notwithstanding, mid-to-late Miocene ignimbrites in the surrounding area (59) are being studied to test directly whether they may represent erupted products related to this pluton.

Together, the petro/thermochronology and numerical simulations place temporal and thermal constraints on melt segregation and imply that ~50 km3 of rhyolite was extracted from an upper crustal mush to form the Huemul pluton within ~130 ka at depths as shallow as 3.5 km (10). Interstitial melt segregation is a widely invoked process for the generation of many styles of eruptible rhyolite from upper crustal magma reservoirs (51). Our findings imply that this process can occur to form plutons and operate over transient time scales while still preserving typical long-term emplacement rates. This confirms melt segregation as a viable process of upper crustal differentiation to produce both rhyolite and high-silica granite, highlighting shallow plutonic systems as important links to the volcanic realm whereby similar melt compositions are created via mush distillation before eruption.


U-Pb zircon geochronology

Zircons were separated from each rock sample using standard techniques and dated following Schoene et al. (60); only a brief summary is given here. Zircons were dated using CA-ID-TIMS at Princeton University. The methods for chemical abrasion are modified from those presented by Mattinson (61) to allow for single grain analysis (62). Grains were spiked with the EARTHTIME 205Pb-233U-235U isotopic tracer [i.e., ET535 (63, 64)], and U and Pb were purified from the dissolved sample using methods modified from Krogh (65). Column rinses were set aside for TEA (37), with methods described below.

Uranium and Pb isotopes were measured on an IsotopX Phoenix62 TIMS at Princeton University. Lead was run as a metal and measured by peak hopping on a Daly photomultiplier. Uranium was analyzed as UO2 and was measured statically on series Faraday cups. Measured ratios were corrected assuming an 18O/16O of 0.00205 ± 0.00004 (2σ), corresponding to the modern atmospheric value (66). Corrections for mass-dependent fractionation of Pb were done using measurements of the 202Pb/205Pb ratio on other samples measured with the ET2535 tracer solution. Uranium mass fractionation was monitored in real time using the 233U/235U in the ET535 tracer and assuming a 238U/235U of the zircon of 137.818 ± 0.045 (2σ), which represents the mean value of 238U/235U measured in natural zircon (67).

The zircons were measured over two time periods, which requires slight modifications in several aspects of data reduction. First, the dead time determination for the Daly ion counter, which drifts by ~1 ns as a function of time over months to years, was adjusted to match the time when specific zircons were analyzed. Second, the Pb blank composition used for common Pb subtraction also was monitored over time, and the data were reduced with the corresponding blank composition. This shift corresponds to the use of side filament heating for organic interference reduction beginning in 2017; see methods in the work of Schoene et al. (60).

A correction for initial secular disequilibrium in the 238U-206Pb system due to the exclusion of Th during zircon crystallization (68) was made for each analysis using a ratio of zircon/melt partition coefficients (fThU) that is specific to the composition of the magma. While these partition coefficients are not perfectly known, this choice has little impact on the results from this study. For gabbroic to dioritic compositions, an fThU value of 0.4007 is used (69). For granodioritic compositions, an fThU value of 0.2455 is used (69). For quartz monzonite compositions, an fThU value of 0.2862 l is used (69). For granitic compositions, an fThU of 0.1379 is used (53).

All data reduction was done with the Tripoli and ET_Redux software packages (70) using the algorithms presented by McLean et al. (71). The U decay constants are from the work of Jaffey et al. (72). All uncertainties in this study are presented as ±2σ and include internal uncertainties only. Uncertainties in the tracer composition and decay constants are negligible in this particular case for examining cooling rates through comparison with the 40Ar/39Ar thermochronology.

Zircon TEA

The TE composition of the unknown zircons was analyzed using the TIMS-TEA method of Schoene et al. (37). It uses the wash solutions obtained during anion exchange column chemistry during U-Pb separations. This solution contains all of the TE from the dissolved volume of zircon, except U and Pb, and can be analyzed by solution ICP-MS. To accomplish this, the TE aliquots were dried down and then redissolved in 3% HNO3 + 0.2% HF + In (1 part per billion). The resulting solution was analyzed on a Thermo Fisher Scientific iCAP quadrupole ICP-MS using a Teledyne Cetac ASX-100 autosampler. Uptake time was 20 s, and 60 s was used to wash the line with 3% HNO3 + 0.2% HF between each analysis. Measured elements include Y, Zr, Hf, REE, Th, and In (used as an internal standard). A dilution series of a synthetic zircon solution was used to generate a concentration-intensity calibration curve over the range of concentrations observed in most zircon TIMS-TEA analyses. Reproducibility was assessed using a homogeneous solution of Plešovice zircon (73) in addition to a solution with a known Zr/Hf ratio of 50. Measurements of procedural blanks were used to monitor for laboratory TE contamination. Sets of four unknowns were bracketed by measurements of the Zr/Hf and Plešovice standard, and a new calibration curve was made for every 20 unknowns. Measurement of these unknowns was done concurrently with the dataset presented by Chambers et al. (74), and the total number of standards run during these sessions is reported here. Solution measurements were converted to zircon concentrations by assuming that all of the measured TE substitute for Zr4+, such that Σ = Zr + Hf + Sc + Y + Nb + Ta + REE = 497,646 parts per million (ppm). The results from the analyses are presented in table S2. Procedural blanks show no significant laboratory TE contamination for the elements of interest. Repeat analyses of the Zr/Hf ratio in our Zr/Hf standard solution give a weighted mean value of 50.53 ± 0.29 (2σ; mean square weighted deviation, 0.30), within ~1% of the reference value of 50. Repeat runs of the Plešovice zircon solution are highly reproducible for REE heavier than Nd (fig. S1), with one outlier, and provide confidence in our reproducibility of natural zircon TE concentrations.

40Ar/39Ar thermochronology

Twenty-nine 40Ar/39Ar incremental heating experiments were conducted on amphibole, biotite, and orthoclase from the RBH plutonic complex along with eight hornfels wallrock samples (groundmass and plagioclase) in the WiscAr Laboratory at the University of Wisconsin-Madison. Mineral and groundmass separates were isolated by crushing, sieving from 180 to 500 μm, magnetic sorting, density separation using methylene iodide, and then hand-picking under a binocular microscope to remove altered fractions. The purified separates were wrapped in aluminum foil, placed in 2.5-cm aluminum disks, and irradiated along with the 28.201-Ma Fish Canyon sanidine standard (75) at the Oregon State University TRIGA Reactor in the Cadmium-Lined in-Core Irradiation Tube (CLICIT). Approximately 5 to 20 mg of each separate was placed in either a 3 mm by 20 mm trough or a 2 mm by 2 mm well within a copper planchette and incrementally heated with a 25-W CO2 laser following the procedures by Jicha and Brown (76). Gas was analyzed using a MAP 215-50 mass spectrometer. 40Ar/39Ar dates (table S3) are calculated using the decay constants of Min et al. (77), and analytical uncertainties, including J contributions, are reported at 2σ. Corrections for radioactive decay of 39Ar and 37Ar were made using the decay constants of Stoenner et al. (78). The atmospheric 40Ar/36Ar ratio of 298.56 ± 0.31 is by Lee et al. (79), and interfering isotope production ratios are by Jicha and Brown (76) and Renne et al. (80). Each 40Ar/39Ar experiment yields an isochron with an atmospheric intercept; thus, we discuss and interpret the plateau dates (table S3). Uncertainties in the decay constants are negligible in this particular case for examining cooling rates through comparison with the U-Pb data. Data reduction was performed using ArArCALC software (81). The bulk of the 40Ar/39Ar incremental heating experiments on minerals within RBH rocks exhibit well-defined concordant age spectra comparable to volcanic experiments as a result of rapid cooling. In some instances, initial low-temperature steps display younger apparent ages typical of argon loss and/or slightly older apparent ages due to either excess argon or 39Ar recoil effects (45). In both cases, these steps are excluded from concordant plateaus within age spectra that display atmospheric intercepts in isochron plots. Some hornfels samples yield saddle-shaped spectra with older apparent ages at the high- and low-temperature steps with no clear plateau age (44), while others form concordant plateaus. The majority of the hornfels experiments display older steps (~17 to 9 Ma) toward the high-temperature end of experiments (table S3). These older tails in apparent age approach the inferred age of the Eocene to Miocene volcanic protoliths (59). Hence, the lower bounds of these “saddle” spectra are plotted in Fig. 1B as an estimate of the timing of protolith resetting, but they are not considered plateau dates.

Numerical modeling

We use a 3D finite-volume thermal model to examine the residence time of melt in the RBH system. This model builds upon previous models that have examined response of the crust to the flux of basaltic magma over long time scales (82, 83), extraction of melt from a crystalline residue (41), and the interaction of tectonics and melt flux in modifying the thermal state of the upper crust (84). The previous approaches used 2D geometries, but here, we use a 3D geometry to describe the melt residence time, compositional evolution, and extraction for prescribed intrusion scenarios. We use a multiscale approach and, over the full domain, compute the enthalpy balance, composition, and displacement of intrusions to conserve mass in response to successive intrusions. We use a multiscale modeling approach to make the computation tractable. In the macroscale model, we specify that the sum of melt and crystal volume fractions is equal to 1 (Eq. 1) and solve the conservation of mass and energy (Eqs. 2 and 3), respectivelykϕk=1(1)t(ϕkρk)+xi(ϕkρku¯i)=Rk(2)ϕkρkck[T¯t+u¯iT¯xi]=q¯xi+ϕkRkL(3)

Here, ϕ denotes volume fraction, ρ denotes density, u denotes velocity terms, q denotes heat flux, R are mass exchange terms associated with phase change, and L is latent heat. The k subscripts denote the distinct phases, and the i subscripts denote orthogonal directions (i = 1, 2, and 3 in these 3D simulations). The overbar denotes spatially averaged quantities over a simulated grid cell. We note that, in this context, u¯i, the average velocity, is solely implemented to displace the crust and intrusions to accommodate successive intrusions (84). There is some uncertainty in the accommodation mechanism, and for the purposes of these calculations, we assume that the intrusions are accommodated by crustal shortening resulting in vertical translation (35).

To determine the sensible and latent enthalpy balance during melting and crystallization and mass partitioning during phase change (R terms in Eqs. 2 and 3), we use the rhyolite-MELTS software (49) and have a proxy composition for the surrounding crust and all intrusions based on specified major oxide composition. In this way, we can associate composition and phase information to each location in the entire computational volume similar to the approach used by Dufek and Bachmann (41). During crystallization, the Rk value for the magma phase would be positive and have an equal but opposite complementary term in the crystal phase equation. While this approach admits the possibility of examining kinetics in crystallization, here, we use rhyolite-MELTS calculations to determine the appropriate mass exchange and set the rates such that, over each time step, we maintain equilibrium. Likewise, the latent heat terms are derived from the MELTS calculations for self-consistency. We use a temperature-dependent thermal conductivity of average continental crust as determined by Whittington et al. (85).

The large overall computational volume in these calculations (40 km by 40 km by 60 km) makes the detailed examination of mixing of crustal and intrusion melts prohibitive, and hence, we have neglected this mixing here. For these conditions, post facto examination of the simulations found minimal upper crustal melting (<5% total volume of all melts in the simulations) at the margins of the intrusions, in line with isotopic evidence that suggests that the absence of radiogenic crustal melting contributed to the composition evolution of these melts (10). This enabled an optimized scheme for tracking the compositional proxies in the simulation domain.

As the large, 3D computational volume currently precludes resolving crystal scale momentum conservation, we have adopted a one-way coupled approach to examine a proxy for melt extraction in subgrid scale domains. Here, we use the melt fraction–dependent drag relationships (86, 87) used by Dufek and Bachmann (41) primarily developed in multiphase flow engineering applications to compute the time-dependent flux of melt and enthalpy. Over each time step, a subgrid model is called for each location in the spatial domain to compute melt extraction rates. The drag between magma and crystals, D, is dependent on the volume fraction of crystals. We assume that, above the close packing volume fraction of crystals, the single particle drag factor is based on particle Reynolds number using the Schiller and Naumann correlation, and the total drag factor of a collection of particles is given by the correlation of Di Felice to account for hindered settling at higher particle volume fraction (86). The close packing volume fraction of crystals is unlikely to be a single value over the range of shapes, grain sizes, and stress conditions experienced, and in these calculations, we use a value of 0.64 ± 0.2, where the values are selected randomly in this range for each subgrid calculation. For the range of volume fractions of crystals that we considered here, we examine only the interplay between melt-crystal drag and buoyancy. Here, drag, D, is given byDi=ftc(um,iuc,i)(4)

Here, f is the total drag correlationf=f0ϕmβ(5)

The single particle drag correlation is given by (87)f0=1+.15 Rep.687(6)which is appropriate for the relatively limited particle Reynolds numbers (Rep) encountered in magmatic flow. The exponent, β, is given asβ=3.70.65exp[(1.5log(Rep)22.0](7)

The crystal’s response time, tc, istc=ρcdc218 μm(8)

Here, dc refers to the crystal diameter, and μm is the melt viscosity (computed on the basis of temperature and composition using rhyolite-MELTS). When the volume fraction of crystals exceeds the close packing limit, the drag force is given asDi=μmK(um,iuc,i)(9)where the permeability, K, is given by the Blake-Karmen-Kozeny relationshipK=ϕm3dc28M(1ϕm2)(10)

Here, we use a constant M value of 50.

These relationships produce the greatest relative melt-crystal velocities after crystal lockup (41), with little relative motion at very low melt fraction (where permeability is low) and at high melt fraction due to the low Stokes number of the crystals in these viscous magmas. Physical properties of the melts (such as viscosity) are computed using rhyolite-MELTS (49) for the temperature and composition in each grid domain. As no spatial information of the phases is resolved below the grid scale, we assume that, internally, the melt and solid are homogeneous for each subgrid scale calculation. This is obviously a simplification, and future work should focus on evaluating a range of subgrid arrangements. The current model does not examine deformation associated with tectonic stresses, nor does it resolve crystal clustering or deformation at the small scale or other features at the crystal scale such as nonlocal granular interactions (88). Hence, the extraction time scales here should be viewed as end-members permitted by the enthalpy constraints of the problem.

The entire computational volume extends 60 km below the surface (the base of the crust was taken as 50 km in these calculations) and 40 km by 40 km in horizontal extent (movies S1 to S5). The large-scale model resolution is 100 m in each direction, and simulations are conducted for 1-Ma durations. The initial crust is composed of a layering of metavolcanics and plutons, subvolcanic plutons, and basement plutons, respectively, and the inferred density structure is consistent with seismic profiles from the nearby Laguna del Maule volcanic field, Chile (47). While there is uncertainty in the initial crustal structure, we examined varying this structure in both composition and extent, and it had little effect on the results presented here, as the upper crustal intrusions generated little crustal melting. The initial geotherm considered here has a surface heat flux of 70 mW/m2 but rises substantially in response to the intrusion scenarios. The initial geotherm assumes a representative radiogenic heat production of ~8.0 × 10−10 W/kg near the surface, and we assume thatw the concentration of heat-producing elements decreases exponentially in the crust with a characteristic length scale of 16 km.

The RB intrusions have a fixed geometry based on the mapped extent of these plutons with depths ranging from 6.4 to 3 km using the composition of the gabbro (sample RB1604). The Huemul intrusion(s) were varied in a series of stochastic simulations with a horizontal footprint of 16 km by 5 km. The proxy composition for these intrusions was based on Huemul granite (9). Because of uncertainty in the thickness of the pluton and to evaluate potential multi-intrusion formation scenarios, we examined various volumes by varying the total thickness of the intrusion package and evaluated single intrusion scenarios and protracted intrusion scenarios lasting up to 100 ka. As discussed earlier, we also examine a range of close packing values and crystal sizes ranging from ~0.001 to 0.005 m. Each of these is randomly sampled every time step/location. We conducted a set of 540 separate simulations to evaluate the range of melt volumes and extraction rates through time. Using the major oxide composition and calculated temperatures, we also evaluated where in each cell zircon saturation is first achieved (assuming that Zr is only incorporated in zircon) (36) on the basis of the zircon saturation model (table S4) of Schaen et al. (9). Video files capturing temperature and melt evolution over the course of selected simulations are provided in the Supplementary Materials (movies S1 to S5).


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: A.J.S. appreciates laboratory assistance from A. Pamukcu and K. Samperton during visits to Princeton University for U-Pb TIMS work, along with the support and companionship of N. Garibaldi during several remote field seasons in Chile. This work benefitted from insightful reviews by C. Wilson and L. Caricchi along with editorial handling and comments by C.-T. Lee. Funding: This research was supported by NSF grants EAR-1411779, EAR-1650232, EAR-16502265, EAR-1650156, and EAR-1940994. Author contributions: A.J.S. and B.S.S. conceived the project, conducted field mapping, and collected rock samples. A.J.S. developed the interpretations and wrote the manuscript. A.J.S. generated U-Pb TIMS data under the supervision of B.S. B.S. also generated U-Pb TIMS data. J.D. designed and performed the numerical simulations. M.P.E. generated the zircon trace element data on U-Pb TIMS solutions. A.J.S. generated the 40Ar/39Ar data. B.R.J also generated some 40Ar/39Ar data. J.M.C. edited the manuscript. All authors contributed to editing the manuscript and interpreting the results. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article