Zircon survival in shallow asthenosphere and deep lithosphere

Starting materials

The zircon crystals used in this study are from the ~730 Ma Mud Tank carbonatites (Crohn and Moore 1984; Gain et al. 2019) provided as crystals ~1 cm in size by Sebastien Meffre (UTAS, Hobart, Australia). These crystals were first prepared as doubly polished ~1000 mm thick sections, then made into smaller double-polished chips by slight pressing of the section. The mid-ocean ridge basalt (MORB) glass used in our experiments is a typical moderately differentiated (8.2 wt% of MgO, M = 3.3; G = 2.5, where M and G are factors controlling solubility of zircon in silicate melts, e.g., Harrison and Watson 1983; Gervasoni et al. 2016 and references therein, Supplementary Material^[1]) glassy tholeiitic basalt (number 3786/3) from Knipovich ridge of the Mid-Atlantic Ridge sampled by dredging during 38th Research Vessel Adademic Mstislav Keldysh expedition (Sushchevskaya et al. 2000). The MORB glass (that contains no traces of olivine and/or chromite phenocrysts) has been crushed to powder (<100 mm glass size).

The haplobasatic glass, corresponding the anorthite-diopside eutectic in the system (CaO-MgO-Al₂O₃-SiO₂) was synthetized using mixtures of reagent grade oxides and carbonates (SiO₂, Al₂O₃, MgO, CaCO₃) following the method of Toplis et al. (1994).

Experimental techniques

The reaction of natural Mud Tank zircon with basaltic liquids has been investigated at pressures from 0.1 MPa to 0.7 GPa corresponding to the depths of crustal basaltic chamber (up to 20 km), typical of mafic magma transport conditions in oceanic settings. Seven experimental runs have been completed in this pressure range, at temperatures in the range 1250–1300 °C and low ratios of zircon to basalt (0.04 to 0.16) (Table 1, Fig. 1).

Figure 1

(a) Backscattered electron image of Z1 sample, demonstrating congruent dissolution of zircon (Zrn) in natural basaltic melt at 0.5 GPa. (b and c) Backscattered electron image of Z2 sample illustrating incongruent zircon dissolution in natural basaltic melt. A reaction with basaltic melt at 0.1 MPa causes zircon replacement by baddeleyite and SiO₂-rich melt. (d and e) Backscattered electron image of Z3 sample demonstrating incongruent zircon dissolution in natural basaltic melt. A reaction with basaltic melt at 0.1 MPa causes zircon replacement by baddeleyite (Bdy) and liberation of SiO₂-rich melt (L_SiO2). (f) Z4 sample: the bulk zircon dissolution in natural basaltic melt at 0.7 GPa according to the congruent mechanism. (g) Backscattered electron image of Z6 sample illustrating incongruent zircon dissolution in haplobasaltic melt. A reaction with the melt at 0.7 GPa causes zircon replacement by baddeleyite and liberation of SiO₂-rich melt. (h) Backscattered electron image of Z7 sample illustrating incongruent zircon dissolution in haplobasaltic melt. A reaction with the melt at 0.7 GPa results in zircon replacement by baddeleyite and SiO₂-rich melt (L_SiO2). (i) Backscattered electron image of Z10 sample, demonstrating congruent dissolution of zircon (Zrn) in natural basaltic melt at 0.2 GPa.

Experiments were performed in Au₈₀Pd₂₀ capsules into which a double-polished zircon chip was placed at the bottom, the upper part of the capsule being filled with basaltic glass powder. Starting materials were weighed carefully to control zircon/melt ratio. Distilled water^MQ was added in one experiment (Z7) to study the effect of hydrous basalt (Table 1). This water was introduced with a micro-syringe to the basaltic powder when the sample was already loaded into the capsule. Because of the fast kinetics of zircon dissolution at high temperatures and to avoid the complete dissolution of zircon chip, run durations had to be shorter than the time necessary to reach oxygen fugacity equilibrium with piston-cylinder double-capsule techniques using mineral buffers (48 h: Matjuschkin et al. 2015). The redox conditions in our kinetic experiments are thus considered to be controlled by the initial Fe²⁺/Fe³⁺ ratios of the natural samples used as starting materials (MORB glass), possibly affected by partial exchange of Al₂O₃ with pressure media (in piston-cylinder experiments). The redox conditions established were estimated to range from FMQ(+1.5) to FMQ(+3.6) (where values indicate logarithmic units compared to the FMQ, fayalite-magnetite-quartz buffer) (Borisova et al. 2020). These redox conditions were estimated from the olivine-chromite assemblages using the equations of Ballhaus et al. (1991) and from the mole fraction of ferric iron Fe³⁺ obtained in the basaltic glass by XANES (X‑ray absorption near edge structure).

Two piston-cylinder systems were used in our experiments. The experiment (Z1) used the Max Voggenreiter end-loaded Boyd-England piston-cylinder apparatus at the Bavarian Research Institute of Experimental Geochemistry and Geophysics (BGI), Bayreuth, Germany. Talc cells 3/4 inch in diameter with Pyrex sleeves were used. A tapered graphite furnace was inserted in each cell (Fig. 2a). Alumina (Al₂O₃) spacers were used as pressure-transmitting medium. An Au₈₀Pd₂₀ capsule loaded with starting materials was set in the central part of the assembly. A 20% pressure correction was applied for the friction between the talc cell and pressure vessel. A molybdenum disulfide (MoS₂) lubricant was introduced to minimize the friction. The temperature in the upper part of the capsules was controlled by a EUROTHERM (2404) controller via either W₃Re₉₇/W₂₅Re₇₅ (type D) or Pt₆Rh₉₄/Pt₃₀Rh₇₀ (type B) thermocouple accurate to ±0.5 °C. The uncertainty on pressure was estimated to be ±0.04 GPa. The sample was compressed to 0.5 GPa during a period of 20 min and then heated up to the run temperature (1300 °C) at a rate of 100 °C/min. The samples were maintained at run conditions for the desired duration, then quenched by switching off the power supply. Decompression then lasted 20 min to 2 h. The quench rate to room temperature was up to 80 °C/min.

Figure 2

(a) Principal schema of assemblage of piston cylinder applied at BGI (Germany). (b) Principal schema of assemblage of piston cylinder applied at KIEM (Russia).

Experiments (Z4, Z6, Z7) used the end-loaded Boyd-England piston-cylinder apparatus at the Korzhinskii Institute of Experimental Mineralogy, Chernogolovka, Russia. Standard talc-Pyrex cells 3/4 inch in diameter, equipped with tube graphite heaters and inserts made of MgO ceramics were used as pressure-transmitting medium (Fig. 2b). The pressure at elevated temperatures was calibrated against two reactions of brucite = periclase + H₂O and albite = jadeite + quartz equilibria. A pressure correction (12%) was introduced for the friction between the cell and hard-alloy vessel. To minimize friction, a Pb foil and molybdenum disulfide (MoS₂) lubricant were used. Au₈₀Pd₂₀ capsules with starting mixtures were mounted in the central parts of the cells. The temperature in the upper part of the capsules was controlled to ±1 °C using a MINITHERM controller via a W₉₅Re₅/W₈₀Re₂₀ thermocouple insulated by mullite and Al₂O₃ without pressure correction. For the Pyrex-bearing assemblies, the sample was heated to 550–600 °C at low confining pressure (0.15–0.2 GPa) for a few minutes to soften the Pyrex glass, then both temperature and pressure were increased simultaneously up to the desired run conditions. The uncertainty on pressure was estimated to be ±0.04 GPa. Samples were maintained at run conditions for desired durations then quenched by switching off the power supply. The quench rate was 100–200 °C/min.

Experiments Z2–Z3 used a vertical gas-mixing furnace at one atmosphere pressure, 1300 °C and controlled oxygen fugacity at IRAP, Toulouse, France. Redox conditions in the furnace were controlled by mixtures of CO and CO₂ corresponding to the fayalite-magnetite-quartz (FMQ) buffer. The closed, but unsealed, Au₈₀Pd₂₀ capsules loaded with starting materials were held on a Pt wire (0.3 mm in diameter) and introduced into the hot spot of the gas-mixing furnace. Samples were rapidly quenched (~1000 °C/s) by dropping the sample into the cold zone of the furnace, in the CO-CO₂ atmosphere, by melting a thin Pt suspension wire with an electric shortcut. After evacuation of the CO, the samples were recovered and mounted in epoxy and polished with SiC papers.

For the experiment Z10, we used an internally heated gas pressure vessel at the Korzhinskii Institute of Experimental Mineralogy, Chernogolovka, Russia. The pressure in the system was created by pure Ar gas. The system was heated by a furnace with two windings (to minimize the thermal gradient). The temperature was controlled by a TRM-101 OVEN controller through two S-type (Pt₉₀Rh₁₀ vs. Pt₁₀₀) thermocouples. The thermocouples were mounted at the top and close to the bottom of the run hot spot to monitor the temperature gradient. The duration of experiment was 5 h and the experiment was quenched by switching off the furnace. The pressure during the quench was maintained constant down to 550 °C, and then slowly released. The cooling rate from 1250 to 1000 °C was 167 °C/min, and then 90 °C/min down to 550 °C. Runs Z5, Z8, and Z9 were unsuccessful because of technical issues. After the runs, all successful run capsules were cut in two parts using a diamond saw, mounted in epoxy, and then polished using SiC papers and diamond pastes.

Microanalytical methods

To study the trace element profile of Zr in the reaction zones we used electron probe microanalysis (CAMECA SX-Five). Major element analyses of minerals and glasses (Fig. 3, Supplementary Material^[1]) and sample imaging were made at the Géosciences Environnement Toulouse (GET, Toulouse, France) laboratory and at the Centre de Microcaractérisation Raimond Castaing (Toulouse, France). The main experimental phases (baddeleyite, zircon, and glasses) in the samples were identified by the EDS technique using a scanning electron microscope (SEM) JEOL JSM-6360 LV with energy-dispersive X‑ray spectroscopy (EDS) in GET, Toulouse, France.

Figure 3

(a) The measured Zr contents vs. distance (in micrometers) from the zircon crystal in the Z1 sample (0.5 GPa, 1300 °C, 5 h, Table 1). (b) The corrected Zr profile in the Z1 sample as the function of erf^-1 (1 – C_x/C_o) vs. distance (in micrometers) from the zircon crystal, suggesting zircon dissolution to be controlled by Zr diffusion in the basaltic melt. The uncertainty of the Zr concentrations (in parts per million) are in the limit of 5% relative standard deviation (RSD).

Major and minor (Zr) elements in the crystals and glasses were analyzed using a CAMECA SX-Five microprobe in Centre Castaing, Toulouse, France. For the electron microprobe, operating conditions were: an accelerating voltage 15 kV, currents of 20 nA for the major elements depending on the resistance of the material under the electron beam, and 100 nA for the minor Zr. Focused (~2 mm) beam conditions were used for the major and trace (Cr, Zr) element analyses of the glass phases. The following synthetic and natural standards were used for calibration: natural albite (for Na), natural corundum (Al), natural wollastonite (Si, Ca), natural sanidine (K), synthetic pyrophanite (Mn, Ti), natural hematite (Fe), natural periclase (Mg), synthetic Cr₂O₃ (Cr), and a reference natural zircon (Zr) provided by the Micro-Analysis Consultants Ltd. Element and background counting times were 5 s for Na and K and 10 s for other major elements, whereas peak counting times were 120 s for Cr and 240 s for Zr. The detection limit was 70 ppm for Cr and Zr, this value being derived from the Cameca SX Five software based on the MPI-DING standard glasses. The synthetic reference MPI-DING glasses of mafic and ultramafic composition (KL2-G, ML3B-G, GOR132-G, GOR128-G, Jochum et al. 2006) obtained from natural rock powders were analyzed as unknown samples to monitor the accuracy of the major and trace element analyses. The accuracy estimated on the reference glasses ranges from 0.5 to 5% (1σ RSD = relative standard deviation), depending on the element contents in the reference glasses. The starting MORB glass containing 0.56 wt% H₂O was analyzed according to the method of Bindeman et al. (2012).

Secondary fluorescence effect and kinetic modeling

Given that our overall aim is to model the behavior of zircon in natural systems, we need to extract relevant parameters from our experimental data that can be used in generalized numerical models of dissolution. In this respect, there are two fundamental parameters that control mineral dissolution: (1) the composition of liquid that is saturated in the phase of interest (that can be simplified to the concentration of a limiting element, in our case Zr, at the crystal-liquid interface), and (2) the diffusion coefficient of that limiting element in the liquid. For both of these issues, a potential complication is secondary fluorescence during microprobe analyses that can introduce a spatially variable signature of the element of interest around host crystals (Borisova et al. 2018). Particular attention must thus be paid to eliminate this analytical artifact.

When possible, the interface melt content (C_o) was obtained from all experiments. The linear interface between the zircon grain and glass allows subtraction of the secondary fluorescence effect and to get the correct Zr concentrations in the interface glasses (Table 1; Fig. 1, Supplementary Materials^[1]). An estimate of the Zr diffusion coefficient was extracted from one diffusion profile of the sample Z1 following the method of Harrison and Watson (1983). For that, Zr (parts per million) profiles perpendicular to the zircon-glass interface in the Z1 sample were measured, correcting for secondary fluorescence effects measured between zircon and the same (MORB) basaltic glass (Borisova et al. 2018). The corrected profile taking into consideration of the effect was obtained by linear interpolation of the secondary fluorescence effect with the function “numpy.interp” of Python language (ver. 1.17). (e.g., Oliphant 2006) and by simple subtraction from the measured Zr concentration. The mathematic transformation of the corrected Zr contents in the profiles to the error function (erf^-1) correlates linearly with the distance from the interface with a slope of 0.0022 (Fig. 3, Supplementary Material^[1]), suggesting diffusional control of the Zr distribution due to zircon dissolution in the basaltic melt according to the criterion given by Harrison and Watson (1983). Unfortunately, there is no way to extract the diffusion coefficient values from experiments where zircon dissolves incongruently, compared to congruent dissolution in natural basaltic system at ≥0.2 GPa.

Kinetic modeling of the bulk dissolution of a spherical zircon was completed using the MATLAB software of Bindeman and Melnik (2016) using spherical coordinates with the extracted diffusion coefficient D_Zr = 2.87E-08 cm²/s at 1300 °C and 0.5 GPa (Supplementary Materials^[1]). The parameters used in the calculations were the interface concentration of Zr in the basaltic melt of 20 840 ± 104 ppm and the initial Zr concentration of 94.3 ppm in the starting basaltic melt (Borisova et al. 2020).

Zircon vs. baddeleyite stability during dissolution

In general terms the relative stability of zircon and baddeleyite in terrestrial and extraterrestrial systems can be rationalized in terms of the silica activity of the silicate melt associated with the reaction:

where Zrn is zircon, Bdy is baddeleyite, and L is silicate melt. Equation 1 implies that melts of higher silica activity will favor zircon stability, consistent with the fact that zircon prefers to crystallize in felsic, intermediate, and silica-saturated (high aSiO2,rather than mafic and ultramafic systems (low aSiO2. Indeed, Equation 1 predicts baddeleyite stability in low-a_SiO2 systems, such as ultramafic (e.g., carbonatite, kimberlite) and mafic (haplobasaltic) systems. For example, Gervasoni et al. (2017) investigated the system of zircon in low-aSiO2carbonatite melt at 0.7 GPa and these authors described the stability of baddeleyite and incongruent dissolution of zircon crystals with formation of corona of baddeleyite. In contrast, no baddeleyite crystallization was observed by Shao et al. (2018) in a synthetic basalt system at pressures 0.5–1.5 GPa, implying baddeleyite stability only at lower pressures (<0.5 GPa) in their basaltic systems. The observed effect of pressure on the relative stability of zircon and baddeleyite is consistent with the fact that increasing pressure tends to increase silica activity of mafic silicate melts, as illustrated by the displacement of the olivine-orthopyroxene cotectic. In any case, decompression will favor incongruent rather than congruent dissolution. We also note that comparison of Z2 and Z3 samples from the same time series experiments implies an approximately constant rate of baddeleyite rim formation of 0.9 ± 0.2 mm/min at constant temperature, pressure, and water contents. This rate may be used as a kinetic indicator to calculate approximate timescales of the zircon megacryst transport to the surface.

Comparison of runs Z6 and Z7 suggests that hydrous conditions (Table 1) favor congruent dissolution. Taken at face value, this would suggest that degassing leads to a tendency for incongruent dissolution, but we note that addition of water is generally accepted to decrease silica activity, a discrepancy that requires more study.

Zr diffusivity and zircon solubility in silicate melts

Most existing models considering zircon solubility and dissolution are based on felsic and intermediate and synthetic systems (zircon in felsic and intermediate silicate liquids, e. g., Harrison and Watson 1983; Boehnke et al. 2013; Zhang and Xu 2016; and references therein). In those cases, zircon is apparently less soluble and stable compared to the natural mafic and ultramafic systems. The main factors controlling solubility of zircon in silicate melts are temperature and the melt composition expressed in terms of ratio of nonbridging O atoms per tetrahedrally coordinated cations (NBO/T), which is a simple measure of silicate melt structure based on the ratio between network formers and network modifiers. For example, the factors such as M {[(Na + K + 2∙Ca)/(Al∙Si)] (all in cation ratio)}, G (where numerator (3∙Al₂O₃ + SiO₂) represents the network formers, while the denominator (Na₂O + K₂O + CaO + MgO + FeO) represent the network modifiers of the melt) or B expressed in oxide mole fractions in the melts 0.14XTiO2/XSiO2+1.3XCaO/XSiO2+1.5XNa2O/XSiO2−4.5XK2O/XSiO2−2.7XAl2O3/XSiO22+XMgO/XSiO22−3.7XCaO/XSiO22+75XK2O/XSiO22as well as the melt water contents have been proposed as factors controlling solubility of zircon in silicate melts (e.g., Harrison and Watson 1983; Gervasoni et al. 2016, 2017; Borisov and Aranovich 2019 and references therein). It is expected that the synthetic (Fe- and Ti-absent) and natural Fe- and Ti-containing basaltic melts would not be similar in respect to zircon solubility. Our natural basaltic liquid (8.2 wt% of MgO, M = 3.3; G = 2.5, B = 0.15) provides a theoretical zircon solubility of 23 078 ppm (or 2.3 wt% of Zr) according to model of Borisov and Aranovich (2019), which is in a very good accordance with the experimentally measured solubility of 20 840 ± 104 ppm. The theoretical solubility of zircon in MORB melts is thus comparable to experimental ones (2.2–3.6 wt% of Zr) obtained by Boehnke et al. (2013) for basaltic liquids at 1225 °C and 1 GPa. The factors such as M, G, and B are indicators of zircon solubility in silicate melts but may be also used for prediction of Zr diffusivity, although the most efficient factor controlling the elemental diffusivity in silicate melts besides temperature remains the melt viscosity controlled primarily by the melt SiO₂ content (e.g., Mungall 2002).

To predict zircon survival in natural basaltic systems, the diffusion coefficient of zirconium has been estimated based on the Z1 experiment, as described above (Fig. 3). Using the zirconium diffusion coefficient derived in this way (2.87E-08 cm²/s) and the kinetic model of Bindeman and Melnik (2016), it is predicted that a sphere of zircon 100 mm in diameter will survive for ~10 h, a sphere of 50 mm diameter for ~3 h and a grain 10 mm in diameter for only ~10 min. At the other extreme, the dissolution of a 1 cm crystal such as Mud Tank zircon, would require ~96 000 h (~11 yr).

However, before reaching firm conclusions concerning the timescales of zircon dissolution, we note that literature data on Zr diffusion in synthetic melts (Fig. 5) vary over more than 5 orders of magnitude. Despite this diversity, Zhang and Xu (2016) have proposed a predictive model of Zr diffusivity as a function of the melt composition and water contents (but not taking into account halogen and CO₂ contents). Our experimentally measured Zr diffusivity has been compared to Equation 9 of Zhang and Xu (2016) for our studied melt composition and experimental conditions (Supplementary Material^[1]) as well as to the measured and calculated values from the other experimental works (Fig. 6). We note that our measured Zr diffusivity in natural tholeiitic basalt is comparable to those of synthetic haplobasalt (LaTourrette et al. 1996) and synthetic basalt (Holycross and Watson 2016), close to the line 1:1 on Figure 6, suggesting a good fit between the experiments for dry basaltic systems and the theoretical model for the basaltic systems. Additionally, other experimental data (except for those of Koepke and Behrens 2001) broadly spread along the line 1:1. We thus conclude that the model of Zhang and Xu (2016) may be used to predict Zr diffusivity in volatile-poor natural basaltic melts such as those studied here.

Figure 5

All available and new (this work) experimental data on Zr diffusivity (D) in dry and hydrous silicate melts [in log(D), D in cm²/s] vs. temperature (10 000/T in K^-1). The (haplo-) basaltic systems are represented by data of 1996 LaT = LaTourrette et al. 1996, 2016HW = Holycross and Watson 2016, and this work. The other systems, such as 1983HW representing data of Harrison and Watson 1983, 1998NK = Nakamura and Kushiro 1998, 1999M = Mungall et al. 1999, 2001KB = Koepke and Behrens 2001, 2002B = Baker et al. 2002, 2016ZX = Zhang and Xu 2016, 2018HW = Holycross and Watson 2018 are dry to hydrous intermediate to felsic systems.

Figure 6

All available and new data on the measured Zr diffusivity [log(D), D in cm²/s] in dry and hydrous silicate melts vs. calculated Zr diffusivity [log(D), D in cm²/s] after Zhang and Xu (2016) (calculated on hydrous basis). The (haplo-) basaltic systems are represented by data of 1996 LaT = LaTourrette et al. 1996, 2016HW = Holycross and Watson 2016, and this work. The other systems, such as 1983HW representing data of Harrison and Watson 1983, 1998NK = Nakamura and Kushiro 1998, 1999M = Mungall et al. 1999, 2001KB = Koepke and Behrens 2001, 2002B = Baker et al. 2002, 2016ZX = Zhang and Xu 2016, 2018HW = Holycross and Watson 2018 are dry to hydrous intermediate to felsic systems.