Non-traditional stable isotope behaviors in immiscible silica-melts in a mafic magma chamber

Non-traditional stable isotopes have increasingly been applied to studies of igneous processes including planetary differentiation. Equilibrium isotope fractionation of these elements in silicates is expected to be negligible at magmatic temperatures (δ57Fe difference often less than 0.2 per mil). However, an increasing number of data has revealed a puzzling observation, e.g., the δ57Fe for silicic magmas ranges from 0‰ up to 0.6‰, with the most positive δ57Fe almost exclusively found in A-type granitoids. Several interpretations have been proposed by different research groups, but these have so far failed to explain some aspects of the observations. Here we propose a dynamic, diffusion-induced isotope fractionation model that assumes Si-melts are growing and ascending immiscibly in a Fe-rich bulk magma chamber. Our model offers predictions on the behavior of non-traditional stable isotope such as Fe, Mg, Si, and Li that are consistent with observations from many A-type granitoids, especially those associated with layered intrusions. Diffusion-induced isotope fractionation may be more commonly preserved in magmatic rocks than was originally predicted.


Model
Silica liquid immiscibility was suggested as a method of magmatic differentiation a century ago 35 . However, it is only recently, thanks to progress in experimental and petrographic studies 36 , that A-type felsic rocks at the top of many layered intrusions and in many tholeiitic provinces have been regarded as immiscible silica melts [18][19][20][21] . The process of silica liquid immiscibility is complex and not well understood and is beyond the scope of the present study. Here we only address the isotopic behavior of immiscible  57 FeO or 54 FeO line represents a hypothetical diffusional profile for the corresponding pure FeO isotope endmember. The thickness of the interfacial melt is exaggerated for illustration purpose. See text for details.
Scientific RepoRts | 5:17561 | DOI: 10.1038/srep17561 Si-melts in a mafic magma. Once nucleated, Si-melts, will grow at the expense of the bulk melt. The less dense Si-melts rise toward the top of the magma chamber. The growth rate will depend on the position of the bulk magma within the Si-liquid and Fe-liquid field. A schematic illustration of the FeO and SiO 2 compositional profile during a convective growth of ascending Si-rich melts is presented in Fig. 1.
According to the spinodal diagram of tholeiitic basalt, the interfacial melt of a growing Si-melt equilibrates chemically with a Fe-melt. The growth rate of immiscible Si-melt is determined by the compositional difference between its interfacial melt (Fe-melt) and the bulk melt 37 (Fig. 1). For example, if we assume that SiO 2 is the principal equilibrium-determining component for the Si-melt 38 , the growth rate of a Si-melt is determined by ∆(SiO 2 bulk melt -SiO 2 Fe-melt ), where SiO 2 bulk and SiO 2 Fe-melt are silica contents of the bulk melt and the interfacial Fe-melt, respectively.
Under this growth scenario, a local thermodynamic equilibrium between the growing Si-melt and its interfacial melt (Fe-melt) holds 39 . Therefore, the isotopic composition of a growing Si-melt is equal to its interfacial melt (the Fe-melt) if we ignore the small equilibrium fractionation effect 40 . According to Fick's first law, the diffusive flux of an element for a growing Si-melt is which can be rewritten as where J is the diffusive flux, D is the diffusion coefficient, C is the concentration, x is distance, ∂C/∂x is the concentration gradient, C ∞ is the concentration in the far-field (magma in the chamber), C 0 is the interfacial melt concentration, and BL is thickness of the compositional boundary layer. Due to the difference in diffusion coefficient among isotopes of an element, the compositional difference in isotopes at the interface (∆C) can be expressed as where D 1 and D 2 are diffusion coefficients for heavy and light isotope, respectively. For example, the diffusivity (D 1 ) of the heavy isotope 57 Fe is less than that (D 2 ) of the lighter 54 Fe, therefore the Si-melt becomes enriched in 57 Fe (Fig. 1). Ignoring the small equilibrium fractionation factor for isotopes at high temperatures, we can obtain a solution for a steady-state isotope composition of a crystal growing in an infinite medium: where δ is the difference between isotope ratio in a crystal or a liquid and of the bulk magma (in ‰) following Watson and Miller 40 . Since D 1 /D 2 = (m 2 /m 1 )β 41 , equation 4 can be rewritten as where m 1 and m 2 are masses of heavy and light isotopes, respectively. β is an empirical parameter obtained from experiments 29 and theoretical computation 42 .
The isotopic compositions of a major element of a growing Si-melt can be expressed as where C i,∞ is the concentration of element i in the bulk melt, C i, Fe-melt is the concentration of element i in the interfacial melt, i.e., Fe-melt.
Diffusion of trace elements is often complicated by processes such as: (1) the fact that interface-melt concentration is not fixed by thermodynamic equilibrium, and (2) uphill diffusion 37 (diffusion of a component against its concentration gradient caused by decoupling of concentration and chemical activity of an element 43 ). Therefore, the isotope behavior of a trace element cannot be represented simply by Eq. (6). Here we ignore multicomponent effects and treat trace element diffusion as binary. We obtain the interfacial concentration for a growing Si-melt according to Zhang 37  where C i,∞ is the concentration of trace element i in the far-field melt (magma in the chamber); C i,Fe-melt is the interfacial melt concentration of trace element i, erfc is the complimentary error function, k i is simple partition coefficient of element i between crystal and melt; γ = α(D/D i ) 1/2 , where D and D i are diffusion coefficients of the major component and trace element i, respectively; α is a parameter related to composition or growth rate and is determined by the major component, for example, SiO 2 in a Si-melt, and can be solved by The above solution shows that the isotope fractionation for a trace element is independent of its own concentration and is controlled by the growth rate of a crystal or immiscible liquid (equivalent to α), its partition coefficient, and its diffusivity.
Model predictions. Now we have quantitative models for isotope behavior of both major and trace elements in immiscibly growing Si-melts. We now examine two major elements Fe and Si and two trace elements Mg and Li.
Applying Eq. (6) we have calculated the Fe isotope effect using published experimental data by Charlier and Grove 39 , a β value of 0.015 41 , an initial δ 57 Fe value of δ 57 Fe MORBs , and a C Fe,∞ value taken as (C Fe-liquid + C Si-liquid )/2 and (2/3× C Fe-liquid + 1/3 × C Si-liquid ), respectively (Table 1). Because the concentration of interface Fe is always higher than that of the far-field magma during the growth of a Si-melt, i.e., C Fe, Fe-melt > C Fe , ∞ , the first two terms on the right hand side of Eq. (6) are both positive. Thus, qualitatively, the Si-melt will always be enriched in heavy Fe isotopes (Fig. 1). The calculated result shows that the  Table 1. Calculated results of diffusional fractionation of Fe isotopes in immiscible Si-melts. Note: the first three columns in Table 1 to Table 6 are the same, which are experimental data from Charlier and Grove 39  δ (Si-liquid -bulk magma) increases with increasing SiO 2 content in the Si-melt (Fig. 2a). Comparing with the Fe isotope data of igneous systems published in recent years (Fig. 2a), the unusually heavy Fe isotope enrichment in some of the Si-rich A-type granitoids is consistent with the model predictions.
Examining the data further, it is apparent that the calculated δ 57 Fe, assuming a far-field Fe concentra- Fig. 2a), is higher than those observed in A-type granitoids. However, when the C Fe,∞ value is taken as "2/3× C Fe-liquid + 1/3 × C Si-liquid " (green dash line in Fig. 2a), the observed data match with our modeled ones very well. Additional factors may play a role for values lower than the model-predicted δ 57 Fe values for SiO 2 % higher than 70% (Fig. 2a). For example, there may be some degree of isotope re-equilibrium between the Si-liquid and the bulk melt. In other words, if a Si-liquid cannot be separated from the rest of magma effectively, it will be isotopically homogenized, such as is likely for the Skaergaard intrusion 44 . Also, some A-type granitoids were not formed by the process envisioned by our model but rather by fractional crystallization and re-melting of tholeiitic material 45 , these two process are unlikely to produce the observed Fe isotope fractionation in A-type granitoids, as suggested 16,23 .
Si is also a major element in melts. According to Eq. (6), qualitatively, the δ should be opposite to that of Fe and the Si-melts should be enriched in light Si isotopes because C Si, Fe-melt < C Si,∞ in the right-hand term of the equation. In fact, the ratio of C Si, Fe-melt /C Si,∞ in Si-melts produced by experiments can range from 0.67 to 0.92, if we use the same published experimental data 39 and bulk C Si,∞ = (C Si, Fe-melt + C Si, Si-melt )/2, the calculated δ 30 Si ranges from − 0.52 to − 1.49‰ (Table 2) if we use a Si β factor of 0.047 42 . However, experiments have yielded a near-zero β value for Si isotopes during chemical diffusion 29 , which is very different from 0.047, a value obtained from classical molecular dynamics calculations of a simple SiO 2 -MgO system 42 . One possible explanation for the observed near-zero β value is that the diffusing species of Si is a network former and diffuses as a large species, e.g., as [SiO 4 ] n . If β is near-zero, the  15 and references therein (compiled in Supplementary  Information). The δ 57 Fe value of terrestrial basalts (MORBs) is from Teng, et al. 10 . The diffusional fractionation trend 1 and 2 (orange and green dash lines) are calculated using Eq. (6) based on data reported in Table 1. (b) Mg isotopes from Telus, et al. 16 and references therein (Supplementary Information). The δ 26 Mg value of terrestrial basalts (MORBs) is from Teng, et al. 3 . The diffusional fractionation trend 1 and 2 (orange and green dash lines) are calculated using Eq. (9) based on data reported in Table 2. (c) Plot of δ 30 Si vs. SiO 2 displaying an "igneous array" (blue line) for Si isotopes from Savage, et al. 46 , Savage, et al. 59 , and Zambardi, et al. 23 (Supplementary Information). The diffusional fractionation trend 1 and 2 (orange and green dash lines) are calculated using Eq. (9) based on data reported in Table 3. (d) Li isotope data from Li, et al. 11 and Teng, et al. 60 (Supplementary Information). The diffusional fractionation trend 1 and 2 (orange and green dash lines) are calculated using Eq. (9) based on data reported in Table 4.
δ will be close to zero regardless the value of C Si, Fe-melt /C Si , ∞ in the second term of Eq. (6). Therefore, diffusional enrichment of lighter Si isotopes in Si-melts should be negligible.
Indeed, the observed pattern is different between Si isotopes and Fe isotopes. The δ 30 Si -[SiO 2 ]% plot displays a positive correlation that is shared by A-, I-types of granitoids and basalts (Fig. 2c), and by samples of different locations with distinctly different mineral assemblages, such as Hekla 46 and Cedar Butte volcano 23 . This can be explained by the equilibrium silicate melt structure being an overwhelming control on Si isotope composition, see Fig. 4 in Zambardi, et al. 23 . Although the equilibrium Si fractionation factor between two conjugate immiscible silicate melts has not been calculated or measured, qualitative evidence indicates that heavier Si isotopes are enriched in the more polymerized melts, i.e. heavier Si isotopes increase as the ratio of NBO (non-bridging oxygen) to T (tetrahedron) decreases (Fig. 2c). It is possible that, at equilibrium, bonding with a BO (bridging oxygen) prefers slightly heavier Si isotopes than bonding with NBO in silicate melts. This feature is consistent with the fact that 18 O is also preferred in the immiscible Si-melts, the more polymerized structure melt, a phenomenon observed in experiments 47,48 .
The S-type granitoids are slightly enriched in light Si isotopes with respect to I-and A-type granitoids (Fig. 2c) because the main source of S-type granitoids is sediments 49 which are commonly enriched in light Si isotopes relative to igneous rocks 50 . Overall, diffusion does not seem to play any significant role in Si isotope distribution during magmatic processes.
In A-type granitoids, Mg can be treated as a trace element whose chemical properties are similar to Fe during silicate melt unmixing 39 . Thus Mg's isotope behavior should be similar to Fe's. Indeed, our calculation using Eq. (9) under the same magmatic conditions shows that δ 26 Mg increases with increasing SiO 2 content in ascending Si-melts (Table 3), a prediction in close agreement with the observed trend (Fig. 2b). The reasons for the calculated values being higher than the observed ones are similar to the reasons given for Fe isotopes.
Li is a trace element known to have a high diffusivity in melts. According to Eq. (9), at a high diffusion rate, γ (γ = α(D Si /D Li ) 1/2 ) approaches zero because D Li ≫ D Si , which leads the second term on the right to approach zero as well, resulting in a near-zero δ value (Table 4). So far, observed data do not display any correlation between δ 7 Li and [SiO 2 ]% or among A-type, I-type, and S-type granitoids (Fig. 2d). This is consistent with our model prediction for any element with a high diffusivity. The observed large spread of Li isotope composition of A-type granitoids must be due to other processes.
It is worth noting that our model treats the bulk magma as an infinite reservoir for Fe and Mg. The rationale for this is: (1) Experiments have shown that the evolved silicate melt produced by fractional crystallization from a MORB basalt just prior to silicate melt unmixing is less than 30% of the total volume of the MORB basalt 51        For example, A-type granites can also form via extreme fractional crystallization of a basaltic magma or partial melt of a basaltic parent rock 45 . In this case the Fe and Mg isotope compositions will be controlled by equilibrium isotope effects which generate smaller degrees of isotope fractionation than diffusion induced isotope effects associated with Si-melts formed through immiscibility. Apart for the aforementioned isotope systems, our model is consistent with Zn and Mo isotope behaviors in Hekla rhyolitic melts 52,53 . Although the β parameters for Zn and Mo have not been determined by experiments, the isotope behaviors of Zn in Hekla rhyolitic melts should be similar to those of Fe and Mg considering the association of Zn with Fe-melts 54 and its similar atomic weight with Fe. According to Richter, et al. 29 , the atomic weight of element Mo is too large to have a sizeable β value. Therefore Mo isotope fractionation in Hekla rhyolitic melts should be absent, as has been observed 53 . Our calculated DIE for Zn and Mo with these assumptions can be found in Fig. S1 in Supplementary Information. Our model also gives a testable prediction on isotope behaviors for other systems in A-type granitoids. For example, we predict a large DIE for Ti during silicate melt unmixing due to the larger 50 Ti/ 46 Ti mass ratio and an expected larger β value for Ti 29 . We also predict that highly compatible elements in immiscible Si-melts, e.g. K 54 , should have an isotopic pattern opposite to those of Fe and Mg. Nevertheless, similar to Li, the high diffusion rate of K in a basalt 55 may result in little to no apparent isotope fractionation. In addition, experimental results indicated that diffusion-induced Ca isotope fractionation depends on silicate liquid's composition 30 . This composition-dependent DIE is not fully understood at a molecular level. Our model, combined with the large variations in chemical compositions of the A-type granitoids, may shed light on the puzzling Ca isotope behavior in melts during silicate melt unmixing.
While our predictions await testing by new isotope measurements, we would like to point out one broader implications of the immiscibility-based model. Other processes, such as bubble growth in melts or liquids 56 and carbonatite genesis 57 are controlled by immiscibility. Diffusional isotope effects in erupted volcanic gases or in carbonatites could bear information on the dynamics of igneous processes, as has already been speculated 58 .