The closure temperature(s) of zircon Raman dating

Zircon Raman dating based on irradiation damage is a debated concept but not an established geo-/thermochronological method. One issue is the temperature range of radiation-damage annealing over geological timescales. We conducted isochronal and isothermal annealing experiments on radiation-damaged zircons between 500 and 1000 ^∘C for durations between 10 min and 5 d to describe the annealing kinetics. We measured the widths (Γ) and positions (ω) of the ν₁(SiO₄), ν₂(SiO₄), and ν₃(SiO₄) internal Raman bands, and the external rotation Raman band at ∼974, 438, 1008, and 356 cm⁻¹ after each annealing step. We fitted a Johnson–Mehl–Avrami–Kolmogorov and a distributed activation energy model to the fractional annealing data, calculated from the widths of the ν₂(SiO₄), ν₃(SiO₄), and external rotation bands. From the kinetic models, we determined closure temperatures T_c for damage accumulation for each Raman band. T_c ranges from 330 to 370 ^∘C for the internal ν₂(SiO₄) and ν₃(SiO₄) bands; the external rotation band is more sensitive to thermal annealing (T_c∼260 to 310 ^∘C). Our estimates are in general agreement with previous ones, but more geological evidence is needed to validate the results. The T_c difference for the different Raman bands offers the prospect of a multi-closure-temperature zircon Raman thermochronometer.

1 Introduction

Zircon (ZrSiO₄) is used with several geochronometers because of the substitution of U and Th for Zr in its lattice. Its occurrence in various types of rocks and high chemical and mechanical resistance make it useful for geochronological applications. The α disintegration of U and Th creates lattice disorder by the impact of α particles and the recoil of daughter nuclei. The zircon Raman spectrum is sensitive to lattice damage: the downshift and broadening of the Raman bands provide a quantitative measure for the radiation damage (Nasdala et al., 1995, 1998, 2001; Palenik et al., 2003; Váczi and Nasdala, 2017).

Zircon dating by measuring the accumulated radiation-damage density and the U and Th content was first suggested by Holland and Kulp (1950) and Hurley and Fairbairn (1953), who proposed the quantification of the damage density from the gravimetric density or X-ray diffraction (XRD) analyses. Deliens et al. (1977) calculated “infrared apparent ages” for a suite of Precambrian zircons from the radiation-damage-induced changes in the infrared (IR) spectrum. Pidgeon et al. (1998) analogously suggested zircon dating based on Raman measurements of the radiation damage. The advantages of the determination of radiation damage with Raman over XRD or IR spectroscopy are the simple sample preparation and standardization, short measurement time, high spatial resolution, and the sensitive damage response of the Raman bands over a broad range of damage densities (Nasdala et al., 1995; Pidgeon et al., 1998). Zircon Raman dating has been discussed in several studies (Nasdala et al., 2001; Pidgeon, 2014; Pidgeon et al., 2016, 2018; Hueck et al., 2018; Jonckheere et al., 2019), but it is not an established geochronological method. One of the key issues for the geological application of a zircon Raman geo-/thermochronometer is the quantification of radiation-damage annealing over geological timescales (Nasdala et al., 2001, 2002; Geisler and Pidgeon, 2002).

Radiation damage is annealed at elevated temperatures with the exact temperature depending on the annealing duration (Zhang et al., 2000a; Geisler et al., 2001; Nasdala et al., 2001; Pidgeon et al., 2016). Figure 1 shows the Raman spectra of a zircon subjected to progressive annealing. The Raman bands shift to higher wavenumbers, towards the band positions of well-ordered zircon, and become narrower and more intense. This loss of damage with temperature and time is a problem for the interpretation of zircon Raman dates as crystallization ages (Nasdala et al., 2001, 2002) but unlocks the potential for determining cooling ages and analyzing the thermal histories of natural zircon samples (Resentini et al., 2020). Annealing of radiation damage also affects He diffusion and is thus a process that needs to be taken into account in the interpretation of (U–Th) $/$ He dates of zircon (Ginster et al., 2019; Anderson et al., 2020). Since the Raman analysis of zircon is non-destructive, there is a high potential for including zircon Raman dating into a multi-method geochronological workflow on zircon.

Figure 1 Raman spectrum of an unannealed, radiation-damaged zircon (black, damage density $\approx \mathrm{80}\times {\mathrm{10}}^{\mathrm{16}}\mathit{\alpha }$/g, calculated after Váczi and Nasdala, 2017) compared with the spectra after its cumulative annealing for 1 h at 600, 1 h at 700, and 1 h at 800 ^∘C (red), and additional heating for two 1 h steps at 900 and 1000 ^∘C (orange). The intensities of the latter two spectra are reduced by a factor of 5 for comparison.

A thermochronometer is characterized by its closure temperature T_c, the temperature of the dated sample at the time of its apparent age (Dodson, 1973, 1979). T_c estimates for α-damage annealing range from ∼130 ^∘C for natural samples at isothermal conditions in the German Continental Deep Drilling Program (KTB) borehole (Jonckheere et al., 2019) to ∼650 ^∘C for the re-crystallization of metamict zircon based on the retention of Pb in zircons that were heated to these temperatures (Mezger and Krogstad, 1997). Pidgeon (2014) placed T_c between 230 and 320 ^∘C based on the comparison of Raman ages with other thermochronological data of the same geological units.

Previous laboratory annealing experiments distinguished several annealing stages based on the changes of lattice constants measured by XRD (Weber, 1993; Colombo and Chrosch, 1998a, b) and on changes in the relationship of Raman shift (ω₃) to bandwidth (Γ₃, as full width at half maximum; FWHM) of the ν₃(SiO₄) Raman band (Geisler et al., 2001; Geisler, 2002; Ginster et al., 2019; Dias et al., 2020). Figure 2 plots ω₃ against Γ₃ for the experiments of Geisler (2002) and Ginster et al. (2019). The offset and difference in slope between the damage accumulation and annealing trends are evident. Breaks in slope of the annealing trend mark transitions between the annealing stages. A sharp break separates the steep stage (I) and the flat stage (II), but a more gradual transition occurs between stage II and a stage III assumed by Geisler (2002). Stage I is dominated by elimination of point defects in the basal plane of the crystal structure by tilting and twisting of the ZrO₈ polyhedra (Ríos et al., 2000); stage II is ascribed to crystallization of amorphous domains (Colombo and Chrosch, 1998a; Capitani et al., 2000; Geisler et al., 2001; Geisler, 2002; Ginster et al., 2019); stage III is related to diffusion of point defects that form due to extension in the c direction. These point defects are assumed to be more stable due to the strong linkage of the SiO₄ and ZrO₈ polyhedra in the c direction, and their elimination includes a re-alignment of tilted SiO₄ tetrahedra (Ríos et al., 2000; Geisler, 2002).

Figure 2 Position-bandwidth (ω–Γ) plot of the zircon ν₃(SiO₄) band, showing the stage I–III annealing data of Geisler (2002) and Ginster et al. (2019) compared to the radiation-damage accumulation trend (dashed gray arrow), based on the unannealed starting material of our experiment. The black annealing trajectory shows the difference in slope between the stages leading up to the values of synthetic zircon.

We aim to investigate the change of the major Raman bands due to annealing and to estimate T_c. We track the changes of ω and Γ of the ν₁(SiO₄), ν₂(SiO₄), and ν₃(SiO₄) internal Raman bands, and the external rotation band at ∼974, 438, 1008, and 356 cm⁻¹ (Kolesov et al., 2001), for isochronal annealing runs at different temperatures. The Raman data for these bands and three other external bands that we do not use for further analysis are listed in the Supplement (Tables S1 and S2). We fit two kinetic models to the widths of the three most intense Raman bands for isothermal annealing for different time intervals and temperatures and consider their extrapolation to geological timescales. We discuss the closure temperatures calculated from the models in comparison to previous T_c estimates.

2 Methods and materials

2.1 Zircon samples

We separated zircon from a sample of the Oederan Forest subtype of the Late Carboniferous Schweddey ignimbrite from the Flöha Basin in Saxony, Germany. Löcse et al. (2019) report U–Pb zircon ages of 309.0 ± 1.8 and 309.4 ± 2.6 Ma for the samples (MfNC-2014-01 and MfNC-2014-02) taken from this lithology. We selected these samples because we assume that the zircon has retained the radiation damage accumulated since their crystallization. Zircon separation was carried out as described in Sperner et al. (2014). The zircon grains for the annealing experiments were hand-picked under a binocular microscope. The grains are prismatic with dominant {100} prisms and {101} pyramids. Their lengths range from ∼150 to 300 µm. They are transparent to translucent and colorless to brownish red. Löcse et al. (2019) described oscillatory and sector zoning in zircons from this ignimbrite. Some zircons show inclusions of apatite or quartz. The damage densities calculated from the measured Γ₃ using the calibration of Váczi and Nasdala (2017) range from ∼22 to $\sim \mathrm{200}\times {\mathrm{10}}^{\mathrm{16}}\mathit{\alpha }$/g; they can thus be classified as slightly to highly damaged. We discarded grains with cracks that might fall apart upon heating or cooling. For comparison, we measured a zircon synthesized as pure ZrSiO₄ by Guillong et al. (2015).

2.2 Raman spectrometry

We measured the Raman spectra using a TriVista spectrometer (Princeton Instruments) in single mode connected to an Olympus BX51WI microscope with a 50× (numerical aperture of 0.75) objective. The spot size on the sample was ∼2 µm; the depth likely exceeds ∼5 µm. The power of the 488 nm incident laser light on the sample was ∼12 mW. We chose this laser wavelength to increase the Raman intensity compared to the more traditionally used 514, 532, or 633 nm. Replicate measurements on in-house standards show that the laser power does not affect the lattice damage. The wavenumber calibration used the 219.2, 520.7, and 1001.4 cm⁻¹ bands of sulfur, silicon, and polystyrene. The spectral resolution is ∼0.8 cm⁻¹ and the pixel resolution on the detector is ∼0.2 cm⁻¹. We measured the zircon spectra as a single measurement per spot, acquiring 10 spectra for 20 s each in step-and-glue mode. In step-and-glue mode, several spectra are recorded at different detector positions to cover a larger range of wavenumbers at high pixel resolution. We measured three steps spanning 170 to 1100 cm⁻¹, adding up to a total acquisition time of 10 min. We cut the spectra into three Raman band groups and fitted the bands with Lorentz functions using a third-order polynomial for background subtraction. We measured one spot in most grains, but in some grains we measured one spot at the center and one at the rim. We also excluded spectra showing Raman bands of inclusions or asymmetric bands resulting from low- and high-damage domains (Nasdala et al., 2005). We corrected the bandwidth for the instrumental function following Tanabe and Hiraishi (1980).

2.3 Annealing experiments

We performed isothermal and isochronal annealing runs in a Linn LM111.06 and a Nabertherm LT3/11 muffle oven. For each run, we prepared a set of 6 to 12 zircon grains covering a broad range of initial α damage to be annealed together. The zircon grains were individually wrapped in Monel 400 foil (a nickel/copper alloy) and inserted in the pre-heated oven in a ceramic crucible. The temperatures of the isothermal annealing runs ranged from 500 to 1000 ^∘C for cumulative annealing times of 10, 30, 90 min, 5 h, 1, and 5 d. The experiments followed the approach of Geisler et al. (2001) with each zircon being annealed in consecutive steps at a constant temperature and cooled for the measurement between the steps.

In the isochronal experiments, we annealed the zircons for 1 h runs at 600 to 1000 ^∘C with a 100 ^∘C interval. The Raman spectrum was measured at room temperature after each run. The locations of the measurement spots on the zircon grains were recorded before each annealing step to assure measurements at the same locations. Grains that disintegrated during annealing were discarded. We rejected measurements when we had reason to doubt that they were in the same spot as the previous measurements on the same grain.

3 Results and discussion

3.1 Changes in band position and width

The changes in band position, bandwidth, and intensity are different for each Raman band (Fig. 1) as reported in earlier studies (Zhang et al., 2000a; Geisler, 2002; Dias et al., 2020). Figure 3 shows plots of ω vs. Γ for the isochronal annealing runs (Supplement Table S1).

Figure 3 ω–Γ plots of the ν₁(SiO₄) (a), ν₂(SiO₄) (b), ν₃(SiO₄) (c), and external rotation (ER) (d) Raman bands for the isochronal annealing runs for 1 h between 600 and 1000 ^∘C. Gray circles represent the unannealed zircons from which the annealed samples were selected. The shaded arrows serve as guides, encompassing the range of subparallel annealing trajectories for the annealed zircons from different starting points. The dashed lines approximately mark the lowest- and highest-damage trajectories.

The different bands exhibit a common trend of decreasing Γ and increasing ω with increasing temperature but follow distinct trajectories through the ω–Γ space. Overall, those for the ν₁(SiO₄), ν₃(SiO₄), and external rotation bands resemble that of ν₃(SiO₄) in Fig. 2 with a steep segment at the beginning, followed by a change in the slope towards a flatter trend. We observed the greatest changes for ω₃ and Γ_ER. We interpret the change in the slope at ∼700 ^∘C as the transition from stage I to stage II annealing (Geisler et al., 2001). The slopes of stage I and stage II segments are different for each band. A striking difference exists between the ν₂(SiO₄) band near 438 cm⁻¹ and the other bands. During stage I, ω₂ decreases, whereas the other three bands shift to higher wavenumbers. The decrease of ω₂ reverses at higher temperatures at the onset of stage II of the other bands. Γ₂ values decrease throughout the annealing process like the other bandwidths. The scatter of the ω–Γ data around each common trend is limited, producing a well-defined trend for all Raman bands for stage II, irrespective of the different radiation-damage densities in the unannealed zircons. Our data do not show the gradual steepening between stages II and III as in Fig. 2.

Figure 4 traces the ω–Γ trends for two zircon grains with similar initial radiation damage through isothermal annealing at 600 and 1000 ^∘C (Supplement Table S2). As expected, annealing proceeds faster at 1000 ^∘C. The ω–Γ trends follow the same trajectories as in Fig. 3. The stage I sections of the 1000 ^∘C trajectories must be assumed because the first annealing step already reached stage II.

Figure 4 ω–Γ plots of the ν₁(SiO₄) (a), ν₂(SiO₄) (b), ν₃(SiO₄) (c), and ER (d) Raman bands for the isothermal annealing runs of two zircons at 600 and 1000 ^∘C. The initial damage density was $\sim \mathrm{70}\times {\mathrm{10}}^{\mathrm{16}}\mathit{\alpha }$/g for both zircons. Annealing durations were 10, 30, and 90 min, 5 h, 1, and 5 d for the 600 ^∘C experiment and 10, 30, and 90 min, and 1 d for the 1000 ^∘C. Gray circles represent the unannealed zircons from which the annealed samples were selected. The arrows are orientation guides for the annealing trajectories. The dashed line represents an expected trajectory for the stage I annealing at 1000 ^∘C, for which the zircon already reached stage II during the first annealing run.

Figure 5 compares the change in band positions of the isochronal runs with the results of Zhang et al. (2000a) and in the ω₃ position of Geisler (2002) and Ginster et al. (2019; Fig. 5c). The ω₁, ω₃, and ω_ER data define a rising trajectory up to 800 ^∘C and a flat annealing trend at higher temperatures. This reflects the trends for stages I and II (Figs. 2, 3, 4). For ω₃, most of our data are slightly lower than those of Ginster et al. (2019), consistent with the values of Geisler (2002) but are more strongly annealed than those of Zhang et al. (2000a). We also observed more annealing of ω compared with Zhang et al. (2000a) for the other three bands. The small difference between our isochronal runs for 1 h and the 90 min runs of Ginster et al. (2019) can be attributed to the difference in annealing time.

Figure 5 ω–temperature plots of the ν₁(SiO₄) (a), ν₂(SiO₄) (b), ν₃(SiO₄) (c), and ER (d) Raman bands for the isochronal experiments compared with the 1 h isochronal annealing data of Zhang et al. (2000a) and the 90 min annealing experiments of Geisler (2002) and Ginster et al. (2019).

The main difference between our results and those of Zhang et al. (2000a) and Geisler (2000) relates to ω₂, for which Geisler (2002) reported a sharp increase only in stage III. The data of Zhang et al. (2000a) show a slight decrease of ω₂ during the first two stages that reverses in stage III. Our data show an initial drop of ω₂, followed by a shift to higher wavenumbers at the onset of stage II.

The differences between the annealing trends for the different Raman bands can be interpreted as a result of the different Raman modes. We present a hypothesis for the downshift and reversal of ω₂ during annealing in Appendix B.

Figure 6 plots ω₃ vs. Γ₃ for our isochronal and isothermal runs, superimposed on the data of Geisler et al. (2001), Geisler (2002), and Ginster et al. (2019). From stage II on, the annealing data define a well-defined common trend, even for zircon samples with different initial damage densities. The samples follow subparallel trajectories through stage I. The convergence towards a common stage II is also apparent for the other Raman bands (Figs. 3 and 4). Stage II is interpreted as representing a state of the zircon lattice that is independent of the damage accumulation history. We assume, based on the interpretation of the successive annealing stages of Colombo and Chrosch (1998a), Ríos et al. (2000), and Geisler et al. (2001), that stage II describes zircons in which the lattice has lost most of its point defects and is predominantly strained by the amorphous domains caused by α recoils. In this case, the position of a zircon along the stage II trend represents the remaining amorphous fraction.

Figure 6 Composite ω–Γ plot of the ν₃(SiO₄) annealing data from this study compared with published data (Geisler, 2002; Ginster et al., 2019). The data are subdivided into isochronal and isothermal runs and samples larger and smaller than an initial Γ₃ of 12 cm⁻¹ ($\sim \mathrm{70}\times {\mathrm{10}}^{\mathrm{16}}\mathit{\alpha }$/g). Gray arrow: radiation-damage accumulation trajectory based on the unannealed zircons from this study; dashed black lines: stage I annealing; bold black line: stage II annealing.

3.2 Kinetic modeling and closure temperature

For estimating the temperatures at which annealing takes place on geological timescales, we fitted kinetic models to the Raman bandwidth data for the isothermal annealing runs. We fitted Γ₂, Γ₃, and Γ_ER, but not Γ₁ which shows lower bandwidths than the other bands, implying lower sensitivity to radiation damage. We quantified the fractional lattice repair $\mathrm{\Phi }(t,T),$ following isothermal annealing for a time t and a temperature T equivalent to the parameter α of Geisler et al. (2001):

$$\begin{array}{}\text{(1)}& \mathrm{\Phi }\left(t,T\right)={\displaystyle \frac{{\mathrm{\Gamma }}_i-\mathrm{\Gamma }\left(t,T\right)}{{\mathrm{\Gamma }}_i-{\mathrm{\Gamma }}_{\mathrm{0}}}}.\end{array}$$

Γ_i is the bandwidth of the unannealed sample, and Γ(t,T) that after annealing for a (cumulative) time t at temperature T. Γ₀ is the bandwidth of undamaged zircon; we assumed 5.0, 1.9, and 3.6 cm⁻¹ for the ν₂(SiO₄), ν₃(SiO₄), and external rotation bands, based on the values of the synthetic zircon we measured. Φ=0 indicates no annealing; Φ=1 is complete annealing. We listed values for Φ of the isothermal annealing runs in Supplement Table S2. A Pearson correlation test with Bonferroni correction for multiple testing (Abdi, 2007) showed that none of the co-annealed zircons exhibited a significant dependence of Φ on the initial damage Γ_i. We calculated the arithmetic mean Φ(t,T) values of each experimental condition (t, T) to approximate equal weighting of the different isothermal annealing runs. Mean values for Φ are listed in Supplement Table S3.

Figure 7 plots Φ against annealing time for the three Raman bands; Φ ranges from ∼0.1 to ∼0.7 for Γ₂ and Γ₃ and from ∼0.2 to ∼0.8 for Γ_ER. The trends are approximately linear with logarithmic time and roughly parallel to each other. As expected, Φ increases with time and temperature. The values for ν₃(SiO₄) are consistent with those of Ginster et al. (2019). Ginster et al. (2019) worked with zircons from Sri Lankan gem gravels which have a crystallization age of ∼550–570 Ma (Holland and Gottfried, 1955; Nasdala et al., 2004) but have only retained ∼55 % of all radiation damage accumulated since their formation (Holland and Gottfried, 1955; Murakami et al., 1991; Nasdala et al., 2004). The time of their partial annealing and the geological processes involved are debated (Nasdala et al., 2004). Our volcanic zircons are analogous in age and tectonic setting to the unannealed Saxonian rhyolite zircons of Nasdala et al. (1998) and thus different in provenance, age, and thermal history. Based on the consistency of ν₃(SiO₄) annealing data, we assume that the annealing kinetics of our samples and those of Ginster et al. (2019) are applicable to a broad range of zircons.

Figure 7 Plots of the annealed fraction Φ against ln annealing time for the ν₂(SiO₄) (a), ν₃(SiO₄) (b), and ER (c) Raman bands. The arithmetic means are connected with lines for visual guidance.

We fitted two models to the mean fractions of annealing from the isothermal runs: a Johnson–Mehl–Avrami–Kolmogorov (JMAK) model (Kolmogorov, 1937; Avrami, 1939; Johnson and Mehl, 1939) and a distributed activation energy (DAE) model (Lakshmanan et al., 1991; Lakshmanan and White, 1994). The JMAK model is described by

$$\begin{array}{}\text{(2)}& \mathrm{\Phi }\left(t,T\right)=\mathrm{1}-\exp \left[-{\left(kt\right)}^n\right],\end{array}$$

where n is the Avrami exponent and k is a temperature-dependent rate factor that follows an Arrhenius law:

$$\begin{array}{}\text{(3)}& k=k_{\mathrm{0}}\exp \left({\displaystyle \frac{-E_{\mathrm{A}}}{\mathit{\kappa }T}}\right).\end{array}$$

k₀ is a frequency factor, E_A an activation energy, and κ the Boltzmann constant. JMAK models are used for describing crystallization processes (Avrami, 1939; Johnson and Mehl, 1939). Since crystallization of amorphous domains takes place during radiation-damage annealing, Geisler et al. (2001) and Geisler (2002) used this model for estimating the activation energies of the radiation-damage annealing stages II and III of zircon. The DAE model assumes that the annealing process draws from a distribution of activation energies. It is applied to processes involving sub-reactions with different activation energies and has been used for describing hydrocarbon decomposition and fission-track annealing (Lakshmanan et al., 1991; Lakshmanan and White, 1994). The fractional repair is expressed as follows:

$$\begin{array}{}\text{(4)}& \mathrm{\Phi }\left(t,T\right)=\mathrm{1}-\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}G\left(E\right)\exp \left[-t{k}_{\mathrm{0}}\exp \left({\displaystyle \frac{-E}{\mathit{\kappa }T}}\right)\right]\mathrm{d}E.\end{array}$$

k₀ is a frequency factor and G(E) a Gaussian distribution of activation energies E with mean E₀ and standard deviation σ. We fitted both models by minimizing the sum of squared Φ residuals (SSRs). Table 1 lists and Fig. 8 shows the results. The (mean) activation energies are between 2.7 and 3.0 eV for the three bands and both models. In contrast, the k₀ values span 3 orders of magnitude. The Avrami exponent is similar for Γ₂ and Γ₃ (n=0.11) and lower for Γ_ER (n=0.08). The standard deviations of G(E) are ∼1 eV for the three Raman bands. The best-fit SSRs are comparable for all models, with the lowest values for Γ_ER. The overall agreement of predicted and measured Φ values is close to 1:1 for all experimental conditions (Fig. 8). The SSR surfaces plotted against log k₀ and mean E_A show a distinct trough of low SSR that includes the best-fit parameters.

Table 1 Parameter estimates for the Johnson–Mehl–Avrami–Kolmogorov (JMAK) and distributed activation energy (DAE) annealing models and calculated closure temperatures T_c for cooling rates of 1, 10, 30, and 100 K Myr⁻¹.

Figure 8 Overview of the JMAK (a–c) and DAE (d–f) modeling results for Γ₂, Γ₃, and Γ_ER. The left panels compare the predicted and measured Φ values. The right panels show the sum of squared residual (SSR) surfaces as a function of log k₀ and E_A (JMAK) or E₀ (DAE) for optimal values of n and σ (Table 1).

We estimate the closure temperatures T_c with the approach of Dodson (1979) for fission tracks:

$$\begin{array}{}\text{(5)}& t_{\mathrm{50}}={\displaystyle \frac{-\mathit{\kappa }{T_{\mathrm{c}}}^{\mathrm{2}}}{E\left(\mathrm{d}T/\mathrm{d}t\right)}}.\end{array}$$

The equation considers cooling through the closure temperature following a linear increase of $\mathrm{1}/T$ with time; t₅₀ is the time at which half the damage is retained. E is the activation energy (E_A; JMAK) or its distribution (G(E); DAE). Equation (5) is equated to the model Eqs. (2), (3), and (4), rearranging and substituting 0.5 for Φ(t,T) (Appendix A):

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{E_{\mathrm{A}}}{\mathit{\kappa }{T}_{\mathrm{c}}}}=ln\left[{\displaystyle \frac{-\mathit{\kappa }{T}_{\mathrm{c}}^{\mathrm{2}}{k}_{\mathrm{0}}}{E_{\mathrm{A}}\left(\mathrm{d}T/\mathrm{d}t\right)\sqrt[n]{\ln \mathrm{2}}}}\right]\end{array}$$

for the JMAK model and

$$\begin{array}{}\text{(7)}& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}G\left(E\right)\exp \left[{\displaystyle \frac{\mathit{\kappa }{T_{\mathrm{c}}}^{\mathrm{2}}{k}_{\mathrm{0}}}{E\left(\mathrm{d}T/\mathrm{d}t\right)}}\exp \left({\displaystyle \frac{-E}{\mathit{\kappa }{T}_{\mathrm{c}}}}\right)\right]\mathrm{d}E=\mathrm{0.5}\end{array}$$

for the DAE model.

Equations (6) and (7) are solved iteratively for T_c. Table 1 lists the results for cooling rates of 1, 10, 30, and 100 K Myr⁻¹. As expected, T_c is higher for the faster cooling; the difference is ∼20–30 ^∘C for a 10-fold increase in cooling rate. Values for T_c at 10 K Myr⁻¹ cooling rate range from 260 to 370 ^∘C. The T_c values for all Raman bands are higher for the JMAK than for the DAE models. For both, T_c is highest for Γ₂ and lowest for Γ_ER; T_c for Γ₃ is slightly lower than for Γ₂. The more sensitive response of Γ_ER to annealing compared to Γ₂ is in agreement with the low temperature annealing results of Dias et al. (2020). We interpret the difference in sensitivity between the external rotation band and the internal bands to the stronger Si–O bonds within the SiO₄ tetrahedra and the weaker Zr–O bonds between the tetrahedra (Dawson et al., 1971).

The linear troughs in Fig. 8 reflect a trade-off between E_A (E₀ for DAE models) and k₀. Different parameter pairs fit the data equally well due to the limited range of laboratory annealing times and temperatures (Mialhe et al., 1988; Lakshmanan et al., 1991). The trade-off is a problem for the extrapolation of the experimental data to geological timescales, since T_c varies along the trough.

The kinetic parameters of our JMAK model for Γ₃ can be compared to the JMAK-model results of the Geisler et al. (2001) and Geisler (2002) for the stage II and III annealing. Their Avrami exponent (n=0.11) agrees with ours for Γ₃. Their activation energy is ∼3.8 eV for stage II and ranges from 6.4 to 6.9 eV in stage III; their log k₀ is 9.3 for stage II and >15 in stage III. The differences result mainly from the trade-off between E_A (E₀) and log k₀ and do not necessarily reflect different kinetics. Moreover, most of our data are from annealing stages I and II, whereas those of Geisler et al. (2001) and Geisler (2002) are from stage II and stage III. Their results and those of Ginster et al. (2019) suggest that stage I annealing requires a lower activation energy than stages II and III, which could also in part account for the lower activation energies obtained from our models. The variation of T_c along the SSR troughs in Fig. 8 is also the probable reason for the different T_c estimates for the JMAK and DAE models. There is no independent physical evidence for either model, and both models fit our experimental data equally well (Table 1). Therefore, we assume the DAE value as the lower limit and the JMAK value as the upper limit of the T_c range for each Raman band.

Figure 9a compares our T_c values with previous estimates from geological and experimental evidence. The wide range of T_c (160 to 650 ^∘C) is in part due to the different approaches. That of Deliens et al. (1977) resulted from comparing radiation-damage ages of Precambrian zircons, calculated from an internal bending IR band, with the ages determined with established geochronometers. The IR ages tended to be lower than the corresponding titanite U–Pb ages (T_c≳650 ^∘C; Stearns et al., 2015) and whole-rock Rb–Sr ages but were mostly higher than mica and feldspar Rb–Sr ages (T_c∼320 to 575 ^∘C; Harrison and McDougall, 1980; Giletti, 1991). The zircon radiation-damage T_c estimate of Mezger and Krogstad (1997) is based on the observation that zircons that remained below 600 to 650 ^∘C during parts of their geological history experienced Pb loss by Pb leaching from metamict zones.

Figure 9 (a) Zircon radiation-damage closure temperatures T_c based on annealing experiments (blue) and geological data (green). (b) Partial annealing zones for the JMAK models for Γ₂, Γ₃, and Γ_ER for residence times between 1 Ma and 1 Ga. The shaded colors correspond to the extent of the partial annealing zones for the three Raman bands.

Jonckheere et al. (2019) measured Γ₃ for isothermal holding for ∼80 Myr at increasing temperatures in the KTB borehole and interpreted its downhole decrease as due to stage I annealing. Hueck et al. (2018) and Pidgeon (2014) dated zircons with Raman based on the Γ₃ vs. radiation-dose calibration of Palenik et al. (2003). Hueck et al. (2018) compared their results with corresponding (U–Th) $/$ He zircon ages (T_c≈170 to 190 ^∘C; Reiners et al., 2004) and age–eU modeling results, finding that their Raman ages were consistently higher than the (U–Th) $/$ He ages. Pidgeon (2014) dated zircons from various Australian Precambrian rocks, whose Raman dates were consistent with the biotite Rb–Sr cooling ages (T_c∼320 ^∘C, Harrison and McDougall, 1980) for the same units. Pidgeon (2014) placed the onset of stage I zircon radiation-damage annealing at ∼230 ^∘C. We calculated a closure temperature for the three Γ₃ fanning-linear Arrhenius models of Ginster et al. (2019) at 50 % damage retention. This gives T_c values between 420 and 480 ^∘C. Our model estimates a range from 330 to 370 ^∘C for Γ₂ and Γ₃, reflecting the model-dependent extrapolation of the experimental data to geological timescales.

The lower T_c (260 to 310 ^∘C) for Γ_ER suggest that geological zircon radiation-damage annealing cannot be described by a single T_c. Instead, different Raman bands record different parts of the thermal history of a zircon. The dearth of independent experimental and geological data for Γ₂ and Γ_ER makes it difficult to be certain that their closure temperatures are different from that of Γ₃, as the annealing data suggest. For the best-studied Raman parameter Γ₃, our experimental data favor a closure temperature between 330 and 360 ^∘C in the region between the estimates of Pidgeon (2014) and Deliens et al. (1977).

Figure 9b shows the partial annealing zones for the JMAK models for Γ₂, Γ₃, and Γ_ER. The partial annealing zone temperatures are highest for Γ₂ and lowest for Γ_ER. Under isothermal holding for >1 Ma, partial annealing occurs at temperatures as low as 200 ^∘C, and full annealing requires temperatures above 450 ^∘C. The low-temperature boundary is in agreement with the stage I annealing temperature of Pidgeon (2014) but higher than that of Jonckheere et al. (2019); the upper boundary is consistent with full annealing at 600–650 ^∘C assumed by Mezger and Krogstad (1997).

4 Conclusions

The results of our isochronal and isothermal annealing experiments indicate that the ν₁(SiO₄), ν₂(SiO₄), ν₃(SiO₄), and external rotation Raman bands at 974, 438, 1008, and 356 cm⁻¹ of radiation-damaged zircon anneal differently with respect to the bandwidth (Γ) and band position (ω). Γ decreases for all Raman bands during annealing while ω increases, but ω₂ drops to lower wavenumbers during the first annealing stage, increasing again from the second annealing stage onward. The different annealing trajectories can help to detect partial annealing in natural zircons.

Our ν₃(SiO₄) annealing data on volcanic zircons are consistent with those of Ginster et al. (2019), obtained on samples with different provenance, age, and thermal history. This suggests that the results of the annealing experiments are representative of a wide range of natural zircons.

The JMAK and distributed activation energy models yield closure temperatures (T_c) between 260 and 370 ^∘C. Overall, this temperature range overlaps with most earlier estimates. The different T_c values for the two kinetic models show that model selection and the trade-offs between the model parameters play a significant role for the extrapolation of laboratory annealing data to geological timescales. This uncertainty in extrapolation emphasizes the need for geological data to constrain T_c. Independent of the model, the calculated T_c is comparable for Γ₂ and Γ₃ (330 to 370 ^∘C) but lower for Γ_ER (260 to 310 ^∘C). This difference offers the prospect of multi-T_c zircon Raman dating using several Raman bands.

Appendix A: Closure temperature equations

Dodson (1973) defined the closure temperature T_c as the temperature of a parent–daughter system at the apparent age of a rock. For fission tracks, Dodson (1979) equates the time for 50 % annealing at the closure temperature to the cooling time constant, defined for a cooling following a linear increase of $\mathrm{1}/T$:

$$\begin{array}{}\text{(A1)}& t_{\mathrm{50}}={\displaystyle \frac{-\mathit{\kappa }{T_{\mathrm{c}}}^{\mathrm{2}}}{E_{\mathrm{A}}\left(\mathrm{d}T/\mathrm{d}t\right)}}.\end{array}$$

κ is the Boltzmann constant and E_A is the activation energy of the annealing process at 50 % annealing. For the JMAK model, the fraction of annealing Φ(t,T) is given by

$$\begin{array}{}\text{(A2)}& \mathrm{\Phi }\left(t,T\right)=\mathrm{1}-\exp \left[-{\left(t{k}_{\mathrm{0}}\exp \left({\displaystyle \frac{-E_{\mathrm{A}}}{\mathit{\kappa }T}}\right)\right)}^n\right].\end{array}$$

k₀ is a frequency factor, E_A the activation energy, and n the Avrami exponent.

Rearranging Eq. (A2) yields

$$\begin{array}{}\text{(A3)}& {\displaystyle }-{\left(t{k}_{\mathrm{0}}\exp \left({\displaystyle \frac{-E_{\mathrm{A}}}{\mathit{\kappa }T}}\right)\right)}^n=\ln \left[\mathrm{1}-\mathrm{\Phi }\left(t,T\right)\right],\text{(A4)}& {\displaystyle \frac{E_{\mathrm{A}}}{\mathit{\kappa }T}}=\ln \left[{\displaystyle \frac{t{k}_{\mathrm{0}}}{\sqrt[n]{\ln \left(\frac{\mathrm{1}}{\mathrm{1}-\mathrm{\Phi }\left(t,T\right)}\right)}}}\right].\end{array}$$

Substituting Eq. (A1) for t and 0.5 for Φ(t,T) yields

$$\begin{array}{}\text{(A5)}& {\displaystyle \frac{E_{\mathrm{A}}}{\mathit{\kappa }{T}_{\mathrm{c}}}}=\ln \left[{\displaystyle \frac{-\mathit{\kappa }{T_{\mathrm{c}}}^{\mathrm{2}}{k}_{\mathrm{0}}}{E_{\mathrm{A}}\left(\mathrm{d}T/\mathrm{d}t\right)\sqrt[n]{\ln \mathrm{2}}}}\right].\end{array}$$

Equation (A5) is solved iteratively for T_c, given the model parameters E_A, k₀, and n in Table 1 and assuming a cooling rate $\mathrm{d}T/\mathrm{d}t$.

For the distributed activation energy model, Φ(t,T) is given by

$$\begin{array}{}\text{(A6)}& \mathrm{\Phi }\left(t,T\right)=\mathrm{1}-\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}G\left(E_{\mathrm{0}},\mathit{\sigma }\right)\exp \left[-t{k}_{\mathrm{0}}\exp \left({\displaystyle \frac{-E}{\mathit{\kappa }T}}\right)\right]\mathrm{d}E.\end{array}$$

G(E₀,σ) is the Gaussian distribution of activation energies with mean E₀ and standard deviation σ; k₀ is a frequency factor. Substituting Eq. (A1) for t and 0.5 for Φ(t,T) gives

$$\begin{array}{}\text{(A7)}& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}G\left(E_{\mathrm{0}},\mathit{\sigma }\right)\exp \left[{\displaystyle \frac{\mathit{\kappa }{T_{\mathrm{c}}}^{\mathrm{2}}{k}_{\mathrm{0}}}{E\left(\mathrm{d}T/\mathrm{d}t\right)}}\exp \left({\displaystyle \frac{-E}{\mathit{\kappa }{T}_{\mathrm{c}}}}\right)\right]\mathrm{d}E=\mathrm{0.5}.\end{array}$$

Equation (A7) is solved iteratively for T_c, given the model parameters E₀, σ, and k₀ in Table 1 and assuming a cooling rate $\mathrm{d}T/\mathrm{d}t$.

Appendix B: Hypothesis for the downshift of ω₂

Other than the ν₁(SiO₄) and ν₃(SiO₄) stretching modes, ν₂(SiO₄) is a bending mode whose frequency depends on the angle between the Si–O bonds in the SiO₄ tetrahedron (Geisler, 2002). The O–Si–O angle is related to the ratio of the unit cell parameters a and c (Tokuda et al., 2019).

Figure B1 plots c vs. a for the XRD data for the annealing experiments of Colombo and Chrosch (1998a) and Geisler (2002). During damage accumulation, c increases more than a (Salje et al., 1999; Zhang et al., 2000b). The increased ratio $c/a$ reduces the O–Si–O angle between the oxygen atoms shared by Si and Zr (Tokuda et al., 2019), shifting the ν₂(SiO₄) Raman band to lower wavenumbers. During stage I annealing, the unit cell shrinks anisotropically, reducing a more than c, causing a further increase of $c/a$ and lowering of the O–Si–O angle. The anisotropic shrinkage is thought to be due to the preferential diffusion of point defects in the basal plane of zircon during recovery (Ríos et al., 2000; Colombo and Chrosch, 1998a). We interpret the further decrease of the O–Si–O angle to cause the decrease of ω₂ during stage I. The annealing trend in the c vs. a plot changes to a preferential reduction of c during stage II until the values of well-ordered zircon are reached. The decrease of the $c/a$ values is accompanied by the opening of the O–Si–O angle, which we associate with the reversal of ω₂ during stage II. The scatter in the unit cell data in Fig. B1 partly masks the trends so that the downshift of ω₂ due to the increase in $c/a$ in stage I is not definite but needs future investigation by combined Raman and XRD analyses of annealed zircon.

Figure B1 Zircon unit cell measurements of Colombo and Chrosch (1998a) and Geisler (2002). The colored lines are the annealing trajectories, starting from different damage densities. The dashed arrow is an accumulation trend based on the unannealed zircons. The continuous black and gray lines indicate constant $c/a$ ratios.