Modulation of zircon solubility by crystal–melt dynamics

Zircon dating is commonly used to quantify timescales of magmatic processes, but our ap-preciation of the consequences of internal magma body dynamics lags behind ever-increasing analytical capabilities. In particular, it has been shown that crystal accumulation and melting of cumulates by recharge-delivered heat may affect melt chemistry within magma bodies. We considered the effect of such processes on zircon solubility in highly evolved silicate melts of diverse chemical affinities. Our modeling shows that in most cases cumulate melting perpetuates the zircon saturation behavior of the first melts emplaced at shallow storage levels. Once cumulate melting is established, the ease of saturating in zircon is controlled by cumulate mineralogy, with a particular effect of the amount of cumulate zircon and its availability for resorption. The fidelity of zircon as a recorder of magma system history thus depends on both the system’s chemical affinity and mineralogy, and


INTRODUCTION
Knowledge of zircon stability is key to the modern understanding of magmatic systems. Endowed with both the ability to incorporate U and Th and an exceptional resistance to breakdown, zircon is a robust geochronometer and a powerful tracer of magmatic processes. Indeed, its properties have facilitated major advances in our understanding of the absolute timescales and rates of magma-body growth, maturation, and eruption through U-Th-Pb geochronology of volcanic (Reid et al., 1997;Wotzlaw et al., 2013) and plutonic rocks (Coleman et al., 2004;Samperton et al., 2015). However, magmatictimescale studies produce results that are notoriously difficult to interpret and show first-order disparities, e.g., between the long, 10 5 -10 6 yr apparent periods of zircon crystallization in water-rich, subduction-related volcanic systems (Claiborne et al., 2010;Tierney et al., 2019) and shorter, 10 3 -10 4 yr apparent timescales, mostly from drier melts (Rivera et al., 2014;Wotzlaw et al., 2014). The existing framework for interpreting zircon saturation, growth, and dissolution (built on pioneering work by Harrison and Watson [1983] and Watson and Harrison [1983]) successfully explains most common zircon crystallization scenarios; however, the internal complexity of magma bodies evident from detailed petrological studies suggests that further prog-ress in interpreting complex geochronological data might require an updated conceptual framework for zircon stability in magmas.
A number of observations such as the cooccurrence of crystal-rich and crystal-poor silicic rocks, compositional and mineralogical zoning of ignimbrites, and compositional gaps in cogenetic rock suites have been used to propose the "mush model" (Bachmann and Bergantz, 2004;Hildreth, 2004;cf. Glazner et al., 2015). In this model, much of the evolved melt's compositional identity is decided through the dynamics of solids (crystals), silicate liquids (melts), and fluids (exsolved volatiles) occurring during storage in upper-crustal magma reservoirs (Sparks et al., 2019). One important consequence is the physical separation of crystals from melts resulting in the formation of a high-crystallinity cumulate mush. Such crystal-rich regions are not only passive crystal graveyards but they may be reactivated and melted due to heating by invading recharge magma, thus generating compositionally distinct melts that may contribute to magma-chamber zoning (Wolff et al., 2015). While it is possible to predict how certain chemical parameters would behave as a result of this likely ubiquitous process, such predictions are rarely incorporated into interpretations of zircon data. In this paper, we address this issue by exploring the influence of such crystal-melt dynamics on the composition of evolved silicate melts and their zircon saturation behavior.

ZIRCON SATURATION IN EVOLVED MELTS
Experimental studies of zircon saturation in silicate melts have identified its dependence on two main parameters: temperature and melt composition, with minor water-content and pressure effects Boehnke et al., 2013;Gervasoni et al., 2016).
One key melt-compositional variable is the concentration of Zr, a major constituent of zircon. The other major variable is the degree of silicate melt polymerization, or the balance between network-modifying cations (e.g., Na, K, Ca, Mg, Fe 2+ ) and silicate melt network formers (Si, Al), traditionally represented with parameters such as 1/Si or A/CNK [Al/(Ca + Na + K)]. To simplify the description of melt polymerization, zircon saturation models have expressed these relations with a single compound parameter such as M  or G (Gervasoni et al., 2016). For the purpose of the present argument, we will use the most commonly applied model of Watson and Harrison (1983) and their parameter M, a cation ratio of (Na + K + 2Ca) / (Al × Si), approximately corresponding to (1/Si)/(A/CNK). As M effectively describes melt depolymerization, high values of M correspond to enhanced zircon solubility (Fig. 1A).
The form of the zircon saturation isotherms ( Fig. 1A) points to two complementary processes that can lead to zircon saturation developing from an initially undersaturated melt: (1) decreasing solubility due to cooling; and (2) evolving melt composition, generally toward lower M and higher Zr, by closed-or open-system magmatic differentiation. Magmatic evolution across tectonic settings produces fractionated melts converging on a limiting trajectory (Fig. 1B) linking two end points, the granite and the phonolite minima in the quartz-nepheline-kalsilite system (Wolff, 2017;Schmidt and Weidendorfer, 2018). Whenever that limit is reached, fractional crystallization redirects the residual melts toward one of the two minima and away from a thermal barrier, the alkali-feldspar thermal divide (see Tuttle and Bowen [1958] and Hamilton and MacKenzie [1965] for the shape of the feldspar saturation surface) through crystallization of a feldspardominated mineral assemblage (Wolff, 2017). Closed-system magmatic processes provide no means to cross that thermal barrier; it is only rarely overstepped in cases of significant sialic contamination of magma (e.g., Riishuus et al., 2008). Most common magmatic evolutionary trends generate melts that approach the appropriate zircon saturation curves at M between ∼1.3 and 2.5, and their most differentiated, zirconundersaturated derivatives generally trend away from the composition of their dominant mineral phase, alkali feldspar (M ∼1.7, Zr = 0). We will show that the ability of these melts to reach and maintain zircon saturation is defined largely by their placement in M-Zr space during crustal storage, and the mineralogy of their cumulates.

CRYSTAL-MELT DYNAMICS AND ZIRCON SOLUBILITY
Most evolved magmas capable of saturating in zircon are emplaced in the middle to upper crust, where they undergo storage and further differentiation culminating either in eruption or in cooling to form a pluton (Putirka, 2017). During crustal storage, a host of processes of crystal-liquid separation generate a typical geometry of supernatant melt overlying a crystalline mush (Bachmann and Huber, 2018;Holness, 2018). Thermal energy necessary to maintain shallow magma reservoirs is delivered through repeated recharge of less-evolved, hotter melts generated deeper within the magmatic column (Cashman et al., 2017). In this framework, melt chemistry largely depends on crystal-melt dynamics, i.e., the relative importance of crystallization, physical separation of melts and crystals, melting by added heat, and input of mafic components.
The melt-compositional space of these interactions is defined by (1) the chemical lineage of the magmas ( Fig. 1) and (2) the composition of the cumulates. The mineralogy of the latter varies at different stages of magmatic evolution, but is predictable for the most evolved melts in the zircon saturation region. Most common cumulate mineral phases incorporate small to negligible amounts of Zr (Table 1), effectively making zircon the sole mineral carrier of that element. M values of any cumulate are largely controlled by its most abundant phase, feldspar, whose M varies systematically for plagioclases but is nearly invariable for alkali feldspars (Fig. S1 in the Supplemental Material 1 ). Melts in the vicinity of either the granite or the phonolite minimum typically crystallize one or two feldspars ranging from high-K sanidine to low-Ca plagioclase (anorthite content [An] < 25); all such feldspars have essentially identical M (1.7 ± 0.1). In contrast, M of ferromagnesian phases varies widely from values close to that of feldspar (e.g., biotite) up to well over 100 for Ca-rich pyroxenes; clinopyroxenes in particular have the potential of leveraging cumulate compositions to very high M values.
Knowledge of cumulate phase proportions and mineral chemistry allows us to predict their position in M-Zr space (Fig. 2). Cumulates (brown-gray in Fig. 2) form the fraction-ated solid counterpart to the supernatant, highly evolved, "extracted" melts (red in Fig. 2). Overprinting this, recharge-induced melting of cumulates may generate a family of mixed-composition melts (blue in Fig. 2) that would occupy a broad trajectory between pure cumulate solids and pure extracted melts; their detailed position will depend on the original mineralogy of the cumulate, the degree of partial cumulate melting, and the amount of mixing with the resident melts. Chemical differences reveal an important divide between alkaline and subalkaline melts (sensu Irvine and Baragar, 1971) in their ability to reach and maintain zircon saturation: • In evolved alkaline systems ( Fig. 2A), cumulate mineralogy is heavily dominated by alkali feldspars (Wolff, 2017) and zircon free. Any amount of such low-Zr cumulate melt will tend to dilute Zr in the resulting melt; together with the associated increase in temperature (and allowing for a potential mafic addition), alkaline systems should be expected to only move deeper into zircon undersaturation once cumulate melting is established. • Subalkaline melts evolve toward an M-Zr region where zircon saturation has low Zr concentration requirements (Fig. 2B). Here, all but the earliest cumulates likely contain zircon, therefore most cumulate melts in these systems will be zircon saturated irrespective of the associated temperature increase. Given the high Zr contents of these cumulates, such magmas should be characterized by large-scale inheritance of old zircon crystal domains. In our model, the only ways to reach significant zircon undersaturation of granitic melts are through substantial additions of either (1) mafic to intermediate melts, or (2) melts of cumulates that are zircon poor (or where zircon dissolution is kinetically inhibited) or have particularly high M (e.g., melts rich in dissolved pyroxene component).

CUMULATE MELTING SCENARIOS
We tested our predictions on three chemically distinct cases of volcanic deposits containing 1 Supplemental Material. Analytical data, descriptions of units in Figure 3, details of MELTS modeling, and additional diagrams. Please visit https:// doi .org/10.1130/GEOL.26213S.12221927 to access the supplemental material, and contact editing@ geosociety.org with any questions.  Note: M = (Na + K + 2Ca) / (Al × Si). An-anorthite content. Estimates are based on GEOROC (http:// georoc.mpch-mainz.gwdg.de/georoc/) compilations of mineral chemistry data. In extreme cases, some silicates may reach up to a few thousand parts per million Zr; these cases are limited to zirconundersaturated melts of alkaline affinity. evidence for cumulate melting, where we compared observed melt compositions with the results of thermodynamic models simulating closed-system cumulate melting (Fig. 3). In each case, the dominance of alkali feldspar in the phase assemblage results in strong Ba depletions in extracted melts and corresponding strong enrichments in compositions representing cumulate melts (Wolff et al., 2020).
The zoned, phonolitic-trachytic Campanian Ignimbrite (Campi Flegrei, Italy; Fig. 3A) contains sanidine and low-Ca plagioclase as the dominant cumulate phases (Forni et al., 2016). As a result, the cumulate melting path is dominated by feldspar resorption with the associated Zr dilution. The absence of zircon in the cumulate implies that any amount of cumulate melting in this system should act against zircon saturation.
Carpenter Ridge Tuff (Colorado, USA; Fig. 3B; Bachmann et al., 2014) represents a calc-alkaline system, which saturates in zircon easily, resulting in early zircon removal into the cumulate phase. Cumulate melting therefore occurs under conditions of zircon saturation, with only minimal changes to melt M values as the near-eutectic phase assemblage is progressively digested. This is in agreement with strong Zr enrichments observed in samples interpreted as cumulate melts (Bachmann et al., 2014), the magnitude of which suggests high-degree melting at ∼900 °C.
The calc-alkaline rhyolites of Lipari (Italy; Fig. 3C) are also zircon saturated, but they have been shown to contain clinopyroxene-rich, Zrpoor cumulate fragments (Forni et al., 2015). The low Zr contents of the cumulates lead in our model to the cumulate melts becoming zircon undersaturated following resorption of feldspar at a relatively modest temperature increase (800 °C and 50% melt). This model illustrates that undersaturation may be reached easily (at low degrees of melting) even for zircon-bearing cumulates in cases where Zr supply to the cumulate melt is limited.

CONSEQUENCES FOR ZIRCON IN UPPER-CRUSTAL MAGMAS
Our modeling of cumulates and melts generated through cumulate melting allows predictions to be made about zircon saturation behavior of magmas of diverse chemical affinities. Specifically, we conclude that (1) Cumulate melting perpetuates the zircon saturation behavior of most melts. Alkaline magmas evolve almost exclusively in the zircon undersaturated field ( Figs. 2A and 3A), but the restricted compositional range of their cumulates (M ∼1.7-2.2, Zr ∼0), the likely temperature increase accompanying cumulate melting, and the undersaturated nature of added mafic components may drive such systems into a vicious circle of zircon undersaturation. While on occasion phonolitic melts do saturate zircon (e.g., early erupted portions of Laacher See Tephra, Eifel, Germany; Wörner and Schmincke, 1984), we propose that within a single active magmatic system, the likelihood of zircon crystallization decreases with time as the cumulate melting feedback is established.
In contrast, metaluminous and peraluminous rhyolites have very low requirements for zircon saturation. In such systems, early-crystallized zircon is likely to become part of any cumulate, making any prospective cumulate-derived melt, or its mixtures with the resident melt, zircon saturated-irrespective of the associated temperature increase (Fig. 3B). In these conditions, some cumulate zircon should always be stable through cumulate melting, which increases the potential for long-term preservation of zircon crystals (i.e., long U-Pb crystallization timescales, survival of incorporated country-rock zircon). Melting initial, zircon-free cumulates or substantial mafic additions may bring about transient zircon undersaturation, but for low-M melts, the ease of reaching saturation would quickly deposit zircon-bearing cumulates, effectively starting the cumulate melting feedback and so buffering the system in the zircon-saturated field.
(2) A particular kind of behavior can be expected from systems whose cumulates are zircon poor (e.g., Lipari, Fig. 3C), cases where cumulate zircons are shielded from melting by other mineral phases (e.g., Reid and Vazquez, 2017), or where their dissolution is kinetically inhibited (e.g., Bryan et al., 2008). Such conditions would result in the melt becoming zircon undersaturated as a result of fairly modest temperature increases, even for nominally easily saturated, rhyolitic magmas. This behavior   cumulates: fsp ± qtz, ox + (hbl, bt, cpx, opx) ± zircon

Alkaline melts
cumulates: fsp ± fsptd, bt, hbl, cpx, ox no zircon A B may be exhibited on a variety of scales from local mush environments to entire magma reservoirs, and has the potential of erasing large portions of the zircon crystallization record, thus shortening the apparent pre-eruptive timescales recovered by geochronology. This combination of cumulate mineralogy and repeated cumulate melting may explain the conundrum of supereruptive volumes of high-silica rhyolites apparently generated in geologically short times, such as the Yellowstone-Snake River Plain rhyolites (western United States; Rivera et al., 2014;Wotzlaw et al., 2014).
(3) It follows that the fidelity of zircon as a recorder of magmatic history is system dependent, and may need to be evaluated on a case-bycase basis. A long "pre-history" of a system may be wiped from the zircon record by a specific recharge history operating on a system of particular magmatic chemical affinity and mineralogy. ] is controlled by major mineral phase stability, and Zr by its partitioning into minerals and by temperature-dependent zircon solubility in melt. Gray lines represent solubility at magmatic storage conditions (A, 750 °C in B), at apparent cumulate melting temperature (900 °C in B), and at a modeled transition between zircon saturation and undersaturation (C). fsp-feldspars (san-sanidine, pl-plagioclase); qtz-quartz; bt-biotite; hbl-hornblende; cpxclinopyroxene; opxorthopyroxene; ox-Fe-Ti oxides; ol-olivine; ap-apatite.