Compressibility of single-crystal zircon up to 19 GPa: implications for the partitioning coefficient of trace elements

Zircon is a widely studied accessory mineral that helps us understand the evolution of the Earth. The unique nature of zircon is due to its chemical and physical strength and the characteristic crystallographic zirconium site, which is compatible with rare Earth elements and actinoids. According to the lattice strain model, trace element partitioning coefficients can be predicted by zircon's compressibility of the cation site. However, the crystal structure of zircon has been precisely determined at less than 5 GPa by single-crystal X-ray diffraction measurements. Here, we have precisely determined the crystal structure of zircon under high-pressure conditions of up to 19 GPa, covering the whole stability field of zircon using a diamond anvil cell. We have also theoretically calculated the unit cell parameters of zircon at high pressure by molecular dynamics simulations. The unit cell parameters and the bond length showed abnormal trends above the high-pressure stability limit of zircon. The previous discrepancy in the compressibility is likely due to the non-hydrostatic condition and instability of the sample. We determined the compressibility of the Zr–O bond and ZrO₈ polyhedra based on the determined crystal structure. The averaged length of the Zr–O bond showed a consistent modulus, with the partitioning coefficients predicted by the lattice strain model. On the contrary, the model with a Poisson solid assumption is inconsistent with the partitioning coefficients. The averaged lengths of the cation–anion bonds under high pressure can help in understanding the partitioning coefficient between minerals and melt.

1 Introduction

ZrSiO₄, zircon, is the most important accessory mineral for understanding the ancient history of the Earth and planets since it possesses essential information on the genesis of rocks older than billions of years. Indeed, the oldest fragments on the Earth found so far are dated to 4.3–4.4 Ga by U–Pb decay geochronology using detrital zircon that likely formed in the ancient magma chamber (Froude et al., 1983; Wilde et al., 2001). Moreover, zircon can be used to measure the oxygen fugacity of ancient Earth by means of oxybarometry using Ce and Eu anomalies in the rare Earth element (REE) pattern as old as the Hadean (Burnham and Berry, 2012; Yang et al., 2014). In addition to characterizing igneous rocks, zircon can host radioactive elements such as U and Th, which are major heat sources of the Earth. Since zircon can preserve actinoids for a geological timescale, it is also known to be a promising host material for nuclear waste storage (Ewing, 1999). Indeed, zircon containing > 10 wt % of U was found in the reactant of nuclear fuel and zirconium cladding, likely formed during the accident at the Chornobyl nuclear power plant in 1986 (Anderson et al., 1993; Geisler et al., 2005). On the other hand, the stability of the zircon is used to estimate the shock history of meteorites based on the high-pressure and high-temperature phase transition in ZrSiO₄ (Reid and Ringwood, 1969; Timms et al., 2017). In the field of high-pressure experiments, zircon crystal is often used as a pressure calibrant in hydro-thermal diamond anvil cells thanks to its chemical and mechanical strength (Schmidt et al., 2013).

The unique crystal chemistry of the zircon is due to the characteristic valence state and crystallographic site of the zirconium (Zr) in the zircon crystal. The eight-fold ZrO₈ distorted dodecahedral sites are composed of four shorter Zr–O bonds (2.13 Å at 1 atm) and four longer Zr^′–O bonds (2.27 Å at 1 atm), and the polyhedra are share edges with each other. The compressibility of the zircon has been determined by pioneering single-crystal X-ray diffraction (XRD) studies up to 4.8 GPa using a diamond anvil cell (Hazen and Finger, 1979). Although the most recent study reported unit cell parameters of zircon by single-crystal XRD up to 8.5 GPa, the fine structures, such as atomic coordinates, were not fully determined (Ehlers et al., 2022). While the compressibility of zircon has previously been determined by powder XRD data under pressure lower than 27 GPa (Marqués et al., 2006; Ono et al., 2004a; Van Westrenen et al., 2004) and by theoretical data based on ab initio calculations up to 30 GPa (Stangarone et al., 2019), the compressibility results are not consistent with each other. In general, the single-crystal XRD measurements can resolve the crystal structure on a finer scale than the powder XRD measurements can, even at high pressure (Boffa Ballaran et al., 2013).

Previous experimental studies have shown that zircon undergoes a phase transition to a high-pressure polymorph with CaWO₄-type (scheelite) structure above 10–15 GPa and 1000–2000 K (Ono et al., 2004b; Reid and Ringwood, 1969). The high-pressure polymorph of zircon was found in the Eocene impact ejecta layer of an ocean drilling sample and named reidite (Glass and Liu, 2001). Reidite is believed to be decomposed into ZrO₂ and SiO₂ at pressure higher than 20.5 GPa and temperature higher than 1773 K (Tange and Takahashi, 2004). A diamond anvil cell study showed that zircon transformed into reidite at 23 GPa by cold compression (Knittle and Williams, 1993). Single-crystal XRD study has so far been limited to half of the stable pressure condition of zircon (Hazen and Finger, 1979).

To reconstruct the rock-forming process, laboratory experiments have determined the partitioning coefficient of trace elements between zircon and melt at 1 atm (Burnham and Berry, 2012) and in high-pressure experiments up to 2 GPa–1573 K (Luo and Ayers, 2009; Rubatto and Hermann, 2007; Taylor et al., 2015). Although the partitioning coefficient might be altered if the zircon is formed under a high-pressure condition, partitioning experiments of the trace elements have never been determined above 2 GPa because of technical difficulty. It has been argued that the effect of pressure on the REE pattern is not clear in zircon (Rubatto and Hermann, 2007), while a moderate effect is reported in amphibole and clinopyroxene (Adam and Green, 2003). The ambiguous effect of pressure is partly because of the laboratory experiments' relatively short timescale, which may cause disequilibrium partitioning coefficients (Luo and Ayers, 2009). A crystallographic approach, such as using the lattice strain model, may provide alternative information on the partitioning coefficient of trace elements.

Here, we have determined the crystal structure of zircon up to 19 GPa, covering the whole stability field of zircon by single-crystal X-ray diffraction (XRD) and a diamond anvil cell. We found that the compressibility of the averaged bond length can help in understanding the partitioning coefficients of trace elements between zircon and melt.

2 Experimental procedures

We used gem-quality zircon (annealed at 900 °C for 16 h) from Sri Lanka. The zircon sample was crushed under pestle and mortar, and three grains of zircon were randomly picked and polished for chemical analysis. The chemical composition of the zircon was determined by wavelength-dispersive X-ray spectroscopy (WDX) analysis using a field-emission-type electron microprobe analyzer (JEOL, JXA-8500F) and analyzer crystals of LIFH, PETH, TAPH, LDE1, and LDE6H. Compositional images were taken to check the chemical homogeneity. Qualitative WDX analysis was conducted at 15 kV and 100 nA with an electron beam of 10 µm in diameter. The analysis covered a wide wavelength of 1.2–102 Å. The elements with detectable amounts were selected for further quantitative WDX analysis, which was conducted at 15 kV and 50 nA with an electron beam of 10 µm in diameter.

The zircon crystal was further crushed under pestle and mortar, and a crystal was carefully picked using a polarization optical microscope so that the sample solely contained a single domain. We conducted single-crystal XRD measurements of the samples of about $\mathrm{50}\times \mathrm{50}\times \mathrm{10}$ µm³. Two independent high-pressure experiments were conducted using a Boehler-Almax-type diamond anvil cell with a wide opening angle. The cell is capable of X-ray diffraction measurements under high-pressure conditions. The culets of the diamond anvils were 300 µm in diameter. The samples were put into a 100 µm diameter hole made in the rhenium gasket. The opening angle of the diamond anvil cell was ± 35°. We used a methanol : ethanol (4:1) pressure medium for runs 1 and 2 and a nitrogen pressure medium for run 3, which preserves hydrostatic conditions in the diamond anvil cell up to slightly above 10 GPa (Chen et al., 2021; Klotz et al., 2009).

We used micro-focused X-ray diffractometers (D8 VENTURE, Bruker) at Meiji University with a target of Mo, set to 50 kV and 1.4 mA for single-crystal XRD measurements (Maitani et al., 2024). The software package APEX4 (Bruker) was used for data processing after the measurement, and SHELXL was used for crystal structure refinement (Sheldrick, 2008). After refinement, we calculated the volume, distortion index, and effective coordination number of ZrO₈ polyhedra by VESTA software (Baur, 1974; Momma and Izumi, 2011).

The pressure was determined by ruby fluorescence and Raman spectroscopy measurements of the zircon (Mao et al., 1986; Schmidt et al., 2013). We used a confocal laser Raman spectrometer system (NRS-4500, JASCO) composed of a red-enhanced CCD detector (DR324B, Andor) and Czerny–Turner-type spectrometer (f=200 mm) with 1800 G mm⁻¹ grating for the optical measurements. High spatial resolution along the laser axis was achieved (∼ 1.5 µm) by the confocal microscope system used in the NRS-4500. The laser was a TEM₀₀ single-mode 532 nm laser with an output of 50 mW. The exposure time was less than 1 s for the ruby fluorescence and 30–180 s for the zircon. The laser spot was focused on the sample surface to about 1 µm in diameter. The Si and polypropylene standards were used for the calibration before and after experiments. Moreover, we measured zircon and ruby in 1 atm conditions after the pressure determination to double-check the calibration.

We have theoretically calculated the crystal structure of the zircon at high pressure by molecular dynamics (MD) simulation using MXDORTO (Kawamura, 1996). The interatomic potential models employed in the program offer good explanations about and knowledge of the target material (Kumagai et al., 1994). The velocity Verlet algorithm (Verlet, 1967) and Ewald method (Ewald, 1921) were used to calculate the atomic motions and electrostatic interactions, respectively. For the MD calculations, a fundamental tetragonal cell of zircon in the crystalline state consisting of 375 ZrSiO₄ molecules was used as the initial structure. To equilibrate the fundamental cell, the MD code was run using an NTP (constant temperature and pressure) ensemble at 10 K. All MD calculations were performed with a time step of 0.5 fs.

The calculations were performed using an interatomic potential model (Ikeda-Fukazawa, 2016). The model has two-body interactions for Si–O and a three-body force for the Si–O–Si and system. The two-body interactions consist of Coulomb forces, short-range repulsion, and covalent bonds. The interatomic potential energy of the two-body interaction between atoms i and j, u_ij, is given by

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{u}_{ij}\left(r_{ij}\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }{\mathit{\epsilon }}_{\mathrm{0}}}}{\displaystyle \frac{z_i{z}_j{e}^{\mathrm{2}}}{r_{ij}}}+f_{\mathrm{0}}\left(b_i+b_j\right)\exp \left({\displaystyle \frac{a_i+a_j-r_{ij}}{b_i+b_j}}\right)\\ & +D_{\mathrm{1}ij}\exp \left(-{\mathit{\beta }}_{\mathrm{1}ij}{r}_{ij}\right)+D_{\mathrm{2}ij}\exp \left(-{\mathit{\beta }}_{\mathrm{2}ij}{r}_{ij}\right),\end{array}\end{array}$$

where r_ij is the distance between atoms i and j; ε₀ is the vacuum permittivity; e is the elementary charge; z_i, a_i, and b_i are the partial charge, repulsion diameter, and softness coefficient of atom i, respectively; and f₀ is a constant for the unit adaptation between these terms (41.865 kJ nm⁻¹ mol⁻¹). D_1ij, β_1ij, D_2ij, and β_2ij are parameters for adjusting the radial covalent bonds.

The potential energy of the three-body part, u_jij, is given by

$$\begin{array}{}\text{(2)}& u_{jij}\left({\mathit{\theta }}_{jij},r_{ij}\right)=-f_k\left[\cos \left\{\mathrm{2}\left({\mathit{\theta }}_{jij}-{\mathit{\theta }}_{\mathrm{0}}\right)\right\}-\mathrm{1}\right]\sqrt{k_{\mathrm{1}}{k}_{\mathrm{2}}},\end{array}$$

where

$$\begin{array}{}\text{(3)}& k_i={\displaystyle \frac{\mathrm{1}}{\exp \left\{g_r\left(r_{ij}-r_m\right)+\mathrm{1}\right\}}}.\end{array}$$

θ_jij is the angle between atoms i and j, and θ₀, f_k, r_m, and g_r are parameters for adjusting the angular part of the covalent bond. The potential parameters were determined by constraining the model to reproduce the experimental results of density, unit cell parameters, and compressibility for zircon. Although the parameters were previously proposed by several studies using different potential models (Devanathan et al., 2004; Du et al., 2006), no reported parameter reproduces the unit cell volume and unit cell parameters at high pressure simultaneously. Therefore, we constructed the parameters based on the compression data by single-crystal XRD measurements, referring to the potential model from a previous study (Kubo et al., 2020). We have determined the Zr–O bond length by fitting the pair correlation functions of the equilibrium state.

3 Results

The electron microprobe analysis showed that the zircon was almost pure ZrSiO₄. The chemical composition of the zircon was determined to be Zr = 49.0(2) wt %, Si = 15.9(1) wt %, and O = 32.0(1) wt % by quantitative analysis using WDX. The trace amount of Hf was detected to be 0.7 wt %. The total concentration is slightly less than 100 %, likely because of the surface conditions of the sample. No elements other than Zr, Si, O, or Hf were detected by qualitative WDX analysis conducted across a wide wavelength. The compositional images showed excellent chemical homogeneity of the sample.

We compare the pressure values determined by zircon and the ruby calibrant in run 1. The pressure values were closely consistent with each other under lower-pressure conditions at < 8 GPa (Table 1), as proposed by a previous study (Schmidt et al., 2013). Above the applicable pressure of the zircon calibrant (> 10 GPa), the intensities of the peaks of zircon in the Raman spectra were weak, and the calculated pressure was significantly lower than that of the ruby calibrant (Table 1). Therefore, we preferentially used the ruby scale as a measure for pressure. The pressure values higher than 8 GPa according to the zircon scale (two data points in run 2) are not used for any further analysis because of their lower reliability.

Table 1 Unit cell parameters of zircon at high pressure.

^a Pressure determined by Raman spectra of zircon. ^b Pressure determined by ruby fluorescence.

Table 2 Zr–O bond length of zircon at high pressure.

We have determined the crystal structure of the zircon at high pressure up to 19 GPa by single-crystal XRD using a diamond anvil cell. The obtained number of reflections was typically around 600, covering 6–25°, with a range of −4 < h < 4, −7 < k < 7, and −6 < l < 6. We obtained crystal structure with a tetragonal crystal system and space group I4₁/amd in all measurements, despite not using any crystal structure known a priori during the analysis. The R values were typically less than ∼ 0.05 in any conditions. The obtained unit cell parameters are summarized in Figs. 1, 2, and 3 and Table 1. The Zr–O bond lengths and ZrO₈ polyhedra are shown in Figs. 4 and 5 and Table 2. The unit cell volume (V=260.069(18) Å) at 1 atm was slightly smaller than that reported in previous single-crystal XRD studies: V=260.803 Å (Hazen and Finger, 1979) and V=261.088 Å (Ehlers et al., 2022). This may be due to metamictization by radiational damage, which induces the expansion of the unit cell volume (Woodhead et al., 1991). To clear the metamictization, we annealed the zircon before the experiments. Our unit cell volume is close to that with minimum damage (Woodhead et al., 1991). Regardless of the metamictization, the compressional behavior of the unit cell parameters, bond length, and ZrO₈ polyhedra is in excellent agreement with the results of previous studies using single-crystal XRD (Hazen and Finger, 1979) and later powder XRD (Marqués et al., 2006) (Figs. 3–5). The unit cell parameters agree well with those of a recent study (Ehlers et al., 2022). In contrast, the results of earlier powder XRD measurements are significantly different from our data (Ono et al., 2004a; Van Westrenen et al., 2004) (Fig. 3).

Figure 1 The relationships between the pressure and unit cell volume of zircon. Filled black symbols show the data that are likely taken in the zircon stability field, and open symbols show those in metastable pressure conditions. Gray symbols represent the data calculated by MD simulations. Triangles are run 1, squares are run 2, and circles are run 3 in this study.

Figure 2 The relationships between pressure and unit cell parameter for the (a) a axis and (b) c axis of zircon. The vertical axis is normalized using the values in ambient conditions. Filled black symbols show the data that are likely taken in the zircon stability field, and open symbols show those in metastable pressure conditions. Gray symbols represent the data calculated by MD simulations. Triangles are run 1, squares are run 2, and circles are run 3 in this study.

We have fitted the obtained unit cell volume to the Birch–Murnaghan equation of state with a fixed value of $K^{\prime }=\mathrm{4}$. When we used all data up to 19 GPa, the unit cell parameters showed clear deviation below and above around 8 GPa (Figs. 1, 2, and 3). If we fitted K^′ as a free parameter, we obtained an abnormally high K^′ value, reaching 30. This may be because the experimental pressure exceeds the high-pressure stability limit of the zircon (Ono et al., 2004b; Tange and Takahashi, 2004). Therefore, we have fitted the obtained unit cell volume to the Birch–Murnaghan equation of state using data below 8 GPa (Fig. 3). The bulk modulus K was 234.80(4) GPa with a fixed value of $K^{\prime }=\mathrm{4}$, which agrees well with the 227(2) GPa reported by previous experiments up to 4.8 GPa (Hazen and Finger, 1979).

Based on the refined crystal structure data, we have analyzed the Zr–O bond length in zircon (Fig. 4), distortion index and effective coordination number (Fig. S1), and volume of ZrO₈ polyhedra (Fig. 5) under high-pressure conditions using VESTA software (Momma and Izumi, 2011). The bulk modulus of the ZrO₈ polyhedral volume is 171 GPa with a fixed value of $K^{\prime }=\mathrm{4}$ (Fig. 5). Although it has been reported in the literature that the bulk modulus of ZrO₈ polyhedra is ∼ 280 GPa (Hazen and Finger, 1979), it has been argued that the value should be re-examined because it is anomalously high compared to the value (225 GPa) of SiO₄ polyhedra (Marqués et al., 2006). Our new fitting of previous crystallographic data (Hazen and Finger, 1979) showed that the bulk modulus of ZrO₈ polyhedra is 202 GPa and is comparable with that of SiO₄. The value is fairly consistent with our data.

Figure 3 The fitting results of the Birch–Murnaghan equation of state in the unit cell volume of the zircon. Filled symbols represent our results in the stability field of the zircon, and open symbols represent data in the metastable pressure conditions. The red line is the results obtained from the data below 8 GPa, including a previous single-crystal XRD study (K=234.80(4) GPa with a fixed value of $K^{\prime }=\mathrm{4}$) (Hazen and Finger, 1979). The solid black line is the result using all data. The dash-dotted line, dotted line, and dashed line represent previous studies by powder XRD measurements (Marqués et al., 2006; Ono et al., 2004a; Van Westrenen et al., 2004).

Figure 4 The relationships between pressure and the (a) Zr^′–O bond and (b) Zr–O bond length of the zircon. Filled symbols represent our results in the stability field of the zircon, and open symbols represent data in the metastable pressure conditions. Lines show the results of MD calculations. The vertical axis is normalized using the values in ambient conditions.

The bond length and volume of ZrO₈ polyhedra decreased monotonically with pressure up to ∼ 8 GPa and showed an abnormal trend above the stability field (Figs. 4 and 5). We also determined unit cell parameters and the volume of zircon at high pressure by MD calculations. Both unit cell parameters and the volume are well consistent with the data taken by the single-crystal XRD measurement (Figs. 1 and 2). The parameters used in the present MD calculations are listed in Table 3. This is the first parameter data set for Zr, Si, and O, which is fully consistent with the compressibility of zircon covering all stable pressure conditions.

Table 3 Potential parameters for MD calculations.

4 Discussions

The obtained compressibility of the zircon was not consistent with some previous studies (Fig. 3). The discrepancy may be due to the hydrostatic condition and instability of the sample. The methanol–ethanol(–water) pressure medium was used in two previous studies (Hazen and Finger, 1979; Marqués et al., 2006) and this study. On the other hand, neon (Van Westrenen et al., 2004) and NaCl (Ono et al., 2004a) were used in other studies. While methanol–ethanol and neon preserve homogeneous pressure conditions in the sample chamber up to around 10 GPa, the pressure effect on the line widths of the ruby R1 fluorescence is different (Klotz et al., 2009). This may indicate that the local stress environment is not equal between methanol–ethanol and neon. Since the experimental data are closely consistent with the results of theoretical calculations in the ideal hydrostatic conditions (Figs. 1, 2, and 3), our experimental samples were in hydrostatic conditions below approximately 10 GPa. In the current study, the compressibility of zircon was determined based on the data taken under relatively lower pressure (< 8 GPa). Below this pressure, it is most likely that zircon is stable and under hydrostatic conditions. Alternatively, the pressure calibrant may also cause the discrepancy. Ono et al. (2004a) used equations of state of gold, and others used the ruby pressure scale (Hazen and Finger, 1979; Marqués et al., 2006; Van Westrenen et al., 2004). The metamictization may also have an effect on compressibility because the metamictization induces the disordered structure in the crystal (Ehlers et al., 2022; Özkan et al., 1974; Özkan and Jamieson, 1978). If the number of defects increases due to metamictization, the compressibility may be lower than in the results of this study.

Figure 5 The fitting results of the Birch–Murnaghan equation of state in the ZrO₈ polyhedral volume of the zircon.

We determined the compression behavior of Zr–O bonds (2.13 Å at 1 atm) and Zr^′–O bonds (2.27 Å at 1 atm) of zircon. Both experiments and MD calculations showed that the Zr^′–O bonds are less compressible than Zr–O bonds (Fig. 4). This suggests that the ZrO₈ polyhedra are more distorted under higher-pressure conditions. Further distortion likely causes the instability of the zircon crystal structure toward the phase transition to reidite. The anisotropic compressibility of the a and c axes of the unit cell is due to the strength of the Zr–O and Zr^′–O bonding (Fig. 2), as Zr–O bonds are along the a axis and Zr^′–O bonds are along the c axis of the unit cell. In addition, as Zr–O bonds share corners with SiO₄ polyhedra, the Zr–O bond may become shorter because of the adjacent rigid Si–O bonds.

The partitioning coefficient of trace elements has been theoretically modeled as a function of ionic radii based on strain in the lattice (Brice, 1975; Onuma et al., 1968). Such a relationship between the partitioning coefficients and the ionic radius is known as the Onuma diagram (Onuma et al., 1968). The lattice strain model considers the substitution of trace elements to cause an expansion or contraction of the crystallographic site, which induces an excess or deficit of free energy. According to the model (Blundy and Wood, 1994), partitioning coefficients can be expressed by the equation below.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{D}_i(P,T,X)& =D_{\mathrm{0}}(P,T,X)\\ & \times \exp \left[{\displaystyle \frac{-\mathrm{4}\mathit{\pi }E{N}_{\mathrm{A}}\left[\frac{r_{\mathrm{0}}}{\mathrm{2}}{\left(r_i-r_{\mathrm{0}}\right)}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{3}}{\left(r_i-r_{\mathrm{0}}\right)}^{\mathrm{3}}\right]}{RT}}\right],\end{array}\end{array}$$

where D_i denotes the partitioning coefficients of trace elements with an ionic radius of r_i under given conditions of pressure (P), temperature (T), and composition (X). D_i takes a maximum value (D₀) with the most compatible (fictive) element with an ionic radius of r₀. E is Young's modulus, N_A is the Avogadro constant, and R is the gas constant.

In earlier work, Young's modulus E was estimated from the relationships between E, site bulk moduli K_site, and the Poisson ratio ν: $E=\mathrm{3}\times K_{\mathrm{site}}\times$ (1 – 2×ν), with an assumption of a Poisson solid (ν=0.25) (Blundy and Wood, 1994). Our results showed that the site bulk modulus of ZrO₈ in the zircon was K_site=167(25) GPa, resulting in E=251(38) GPa, according to the Poisson solid assumption. To test the validity of the Poisson solid assumption, we compared the partitioning coefficients between the lattice strain model and experimental data taken at 1 atm and 1573 K (Burnham and Berry, 2012) (Fig. 6). For modeling the eight-fold-coordinated M³⁺ ion, we used D₀=3.32 and r₀=0.956 Å, as proposed by previous studies. The value of E=251(38) GPa clearly does not match the value based on the partitioning coefficients (Fig. 6). The assumption of the Poisson solid may not be appropriate for modeling the lattice strain of the ion site.

Figure 6 The relationship between the ionic radius and partitioning coefficients (Onuma diagram) between zircon and melt at 1 atm and 1573 K. Filled symbols represent the experimental results of the M³⁺ cation, and the open symbols represent those of the M⁴⁺ cation (Burnham and Berry, 2012). Lines are the partitioning coefficient for an eight-fold cation based on the lattice strain model with various Young's moduli (see text).

Alternatively, as a first step, we preliminarily assumed that the ZrO₈ polyhedra are in a similar condition to the incompressible fluid, in which the mechanical pressure is equal to the average of the normal stress – P= $({\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}}+{\mathit{\sigma }}_{\mathrm{33}})/\mathrm{3}$ – because the strains ε of the Zr–O bond are small. The pressure medium preserves hydrostatic conditions in this study, and thus, $P={\mathit{\sigma }}_{\mathrm{11}}={\mathit{\sigma }}_{\mathrm{22}}={\mathit{\sigma }}_{\mathrm{33}}$, and the shear stress is zero. In this preliminary model, the stress applied to the Zr–O bond is equal to the experimentally determined pressure. Young's modulus of the Zr–O bond can be determined by a slope of the strain–stress curve via Hooke's law, $\mathit{\sigma }=E\times \mathit{\epsilon }$ (Fig. S2). The linear regression of strain in the average bond length of Zr–O by this study and previous data (Hazen and Finger, 1979) at < 8 GPa showed Young's modulus E=606(40) GPa. The bond length of Zr^′–O (2.27 Å at 1 atm) showed a high value of E=789(91) GPa, and Zr–O (2.13 Å at 1 atm) showed a lower value of E=412(41) GPa. We summarize Young's moduli estimated by this study in Table S1 in the Supplement. The linear relationship between stress and strain indicates that the partitioning coefficient does not strongly depend on the pressure according to the lattice strain model. We tested various Young's moduli by comparison with the experimentally determined partitioning coefficients (Fig. 6). The obtained Young's modulus of the average bond length (E=606 GPa) agrees well with the experimentally determined partitioning coefficient of the 14 M³⁺ ions determined at 1 atm and 1573 K (Burnham and Berry, 2012) (Fig. 6). In addition to the M³⁺ ion, we fitted D_i for four M⁴⁺ ions with a fixed E value of 606 GPa. The fitting results showed D₀=44(8) and r₀=0.905(9) Å for an eight-fold M⁴⁺ ion. The results are fairly consistent with experiments except for the results of Ti⁴⁺, which may be more compatible with the four-fold-coordinated SiO₄ site. Although the model is based on a preliminary assumption, the compressibility of the average bond length agrees well with the partitioning coefficient of the trace elements. The model can predict the partitioning coefficient based on the compression behavior of the bond length, which is relatively straightforward to measure by single-crystal XRD measurements.

Figure 7 The Onuma diagram between zircon and melt at (a) 0.7 GPa and 1223 K (Taylor et al., 2015) and (b) 2 GPa and 1223 K (Rubatto and Hermann, 2007). Lines are partitioning coefficients by the lattice strain model with various Young's moduli. Filled symbols represent the experimental results of the M³⁺ cation, and the open symbols represent those of the M⁴⁺ cation.

We further compared the model with two available experimental data sets under high-pressure conditions (Rubatto and Hermann, 2007; Taylor et al., 2015) (Fig. 7). We modeled the partitioning coefficient between zircon and melt at 1223 K and high pressure based on the elastic parameters of the zircon that were obtained by this study. For modeling the eight-fold-coordinated M³⁺ ion, we used D₀ and r₀ as proposed by previous studies (Rubatto and Hermann, 2007; Taylor et al., 2015). We found the D₀ and r₀ of eight-fold-coordinated M⁴⁺ ion to be D₀=1329(200) and r₀=0.86(2) Å by fitting the partitioning coefficients of Hf, Zr, U, and Th, since this was not determined in the previous study of Rubatto and Hermann (2007). The partitioning coefficient did not match the data when it was modeled by Young's modulus calculated from the site compressibility with an assumption of a Poisson solid (E=251 GPa). In contrast to the result at 1 atm (Fig. 6), the model using a modulus of averaged Zr–O bond length considerably deviates from the experimentally determined partitioning coefficient (Fig. 7). The value is in between the moduli of a shorter Zr–O bond length (2.13 Å at 1 atm) and averaged bond length. The reason for this shift is not clear. The elements may preferentially push/pull the shortest Zr–O bond during the substitution at high pressure. Alternatively, soft Zr–O bonds are preferentially deformed during the substitution of the cation. To test this hypothesis, the compressibility of the bond length and partitioning coefficient in various minerals at high pressure must be compared. Moreover, recently, Karato (2016) revised the lattice strain model by including the bulk modulus of the trace elements, which is a conceptual parameter that infers the stiffness of the trace elements. The bulk modulus of the trace elements may be largely different under high-pressure conditions, although it is far more difficult to determine than the conventional bulk modulus of the crystal. This may explain the deviation of the model under high-pressure conditions, while the bulk modulus of the trace elements has not been well studied so far.

5 Concluding remarks

We have precisely determined the crystal structure of zircon under high-pressure conditions up to 19 GPa, covering the whole stability field of zircon using a diamond anvil cell and single-crystal XRD measurements. We also theoretically calculated the crystal structure of zircon at high pressure using molecular dynamics simulations. We determined the compressibility of the Zr–O bond and ZrO₈ polyhedra based on the determined crystal structure under pressure conditions lower than 8 GPa, in which zircon is certainly stable. We estimated Young's moduli of the cation site from the bond length to predict the partitioning coefficients of the trace elements via the lattice strain model. The moduli determined based on the averaged Zr–O bond length match the experimentally determined partitioning coefficients well.