Helium diffusion in natural zircon: Radiation damage, anisotropy, and the interpretation of zircon (U-Th)/He thermochronology

Positive Date-eU Correlations

Due to similarities between our positive correlations and those observed with the apatite (U-Th)/He thermochronometer (for example, Flowers and others, 2007; Flowers and others, 2009; Flowers and Kelley, 2011), the effects of radiation damage on He diffusion in apatite provide context for interpreting zircon He positive correlations. Shuster and others (2006) showed that He diffusivity in apatite decreased with increasing damage and subsequent studies (for example, Flowers and others, 2007) supported this conclusion with observations of positive date-eU correlations. We suggest that similar behavior occurs in zircon and leads to positive date-eU correlations. For example, in the context of a thermal history like that depicted in figure 9, grains with different amounts of radiation damage, as well as an initial span of uniform dates, are reset to varying degrees during a reheating event (fig. 9A). If eU is a proxy for the total accumulated radiation damage (that is, all the zircons in a sample have experienced the same t-T history), then those grains with low eU lose a larger fraction of their He, resulting in a younger date, than grains with high eU, leading to a positive correlation. Samples that have undergone slow, monotonic cooling may also exhibit positive correlations. If a sample spends a significant amount of time at temperatures low enough for damage accumulation without annealing, but high enough to be in the PRZ, then arrays of zircon He dates may form a positive correlation (fig. 9A). Both scenarios demonstrate that this correlation results from a sample with grains that span a range in eU, and have resided in the PRZ after disparate amounts of damage have accumulated in those grains.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 9.

Schematic of possible time-temperature paths that produce date-eU relationships. The arrows on each plot's axes indicate the direction in which the given value increases. Varying eU contents lead to differential accumulations of radiation damage, which in turn lead to differential He diffusivities. (A) If accumulation is relatively low, and a sample experiences a thermal pulse, or cools slowly through the PRZ so that damage in-growth and diffusion happen simultaneously, then a positive correlation results. (B) Conversely, if damage accumulation is relatively high (due to older zircons), then a thermal pulse results in a negative correlation, or significant He loss may occur at low temperatures.

A key difference between our interpretation of the damage-diffusivity relationship in zircon and the interpretation for the apatite system, though, is the mechanism that we suggest causes diffusivity to decrease with increasing damage. In the apatite He system, the decrease in He diffusivity is hypothesized to result from accumulation of crystal defects caused by alpha recoil damage that act as He traps. These traps sequester He (governed by an equilibrium partition coefficient) and prevent or slow its diffusive migration out of the crystal (Farley, 2000). In contrast to He diffusion in apatite, various authors (Farley, 2007; Reich and others, 2007; Saadoune and others, 2009) have demonstrated that He diffusion in an ideal or defect-free zircon should occur almost solely along c-axis parallel channels [0 0 1]. Any disruption of these pathways would force He through c-axis orthogonal openings ([1 0 0], [0 1 0], and [1 0 1]), which are much less energetically favorable (Reich and others, 2007). We argue that decreases in He diffusivity in zircon are largely due to the increasing disruption of diffusion fast-paths (c-axis parallel channels) by radiation damage, an effect similar to road blocks being placed on a major highway. As these barriers are erected inside the zircon, a He atom's path becomes more tortuous and the effective diffusivity of the grain decreases. Evidence for this increasing disruption of c-axis parallel channels comes from figure 5. The D₀ values for orthogonal oriented slabs (diffusion predominantly in the c-axis parallel direction) decrease across the damage spectrum whereas the D₀ values for the parallel oriented slabs (diffusion predominantly in the c-axis orthogonal direction) remain the same, with both sets of D₀ values becoming similar at high damage. In other words, diffusion kinetics in the orthogonal direction begin to more closely resemble diffusion kinetics in the parallel direction with increasing damage, and this increasing similarity can be explained by tortuosity. Tortuosity may also contribute to lowering He diffusivities in apatite, but we envision this phenomenon is more important for He diffusion in zircon due to its strongly anisotropic behavior in specimens with little or no accumulated damage. We do not rule out that damage zones in zircon may also trap some amount of He, but we suggest that this effect is secondary compared to the closing of preferred diffusion directions, which are probably not present in apatite.

Negative Date-eU Correlations

Like the positive date-eU correlations, previous research provides some context for interpreting our negative date-eU correlations. Various authors (Holland, 1954; Hurley, 1954; and Nasdala and others, 2004a) have suggested that He diffusivities increase abruptly once zircon reaches a certain threshold of radiation damage. As damage increases, Reiners (2005) proposed that zircon He dates scale with the remaining crystalline fraction of the zircon as determined by the double-overlapping cascade model (Gibbons, 1972). This further suggests that above a threshold, interconnected damage zones form through-going channels in the zircon lattice and create fast diffusion pathways for He. In order to reach this interconnection or percolation threshold, zircons must sustain long-term damage accumulation at temperatures low enough to prevent annealing. Zircons may achieve a heavily damaged state in which either significant He loss and resetting occurs at surface temperatures, or some brief, low-temperature reheating event may cause resetting (fig. 9B). In detail, some amount of less diffusive material must remain in heavily damaged zircons as most negative correlations are gradual (that is Minnesota River Valley, Bighorn suite) and not so abrupt. However, in general, both scenarios could result in negative correlations and are plausible explanations for the datasets plotted in figure 2. Two of these samples, the Minnesota River Valley and Bighorn suites, come from Archean rocks that have likely been within a few kilometers of the surface for 10⁸ to 10⁹ years and have consequently accumulated large amounts of radiation damage. The Sri Lankan zircons are not as old but some do have high eU, which in some cases resulted in the complete breakdown of crystal structure (for example N17, Nasdala and others, 2004a).

The Sri Lankan dataset also shows a potential percolation threshold effect whereby dates are fairly reproducible up to a critical eU concentration, in this case about 2000 ppm eU, above which they decrease with increasing eU. The Cooma date-eU trend reinforces the notion of a transition in diffusion behavior as it displays both a positive and negative correlation with an abrupt transition between each. Damage in-growth since granulite facies metamorphism at ∼433 Ma (Williams, 2001) and a large disparity in eU concentration (∼100-1300 ppm) produced a suite of zircons in which two different processes affected He diffusion in different parts of the eU spectrum. In this particular dataset, the transition between the two types of He diffusion occurs at an eU of ∼1000 ppm. Similarly, nearly all zircons in our positive correlations contain eU concentrations below 1000 ppm. The Cooma dataset serves as an important demonstration of how radiation damage may have contrasting effects on He diffusivity in a suite of grains from a single sample, and shows the approximate concentration of eU over which a transition from one process to the other may occur. Because both types of diffusion mechanism may operate on the same sample, these samples underscore the necessity for understanding the damage-diffusivity relationship across the entire damage spectrum.

Arrhenius Trends

There are a number of possible explanations for the non-linear portions of the Arrhenius plots in figure 4. Some of these plots show concave-up trends in the initial prograde steps, in which less than ∼1 to 2 percent of the He is released. In contrast, subsequent temperature steps show much less variation and are more nearly linear. These concave-up trends can either plot above or below the linear Arrhenius trends. This feature, or ones similar to it, has been observed previously in zircon (Reiners and others, 2002; Reiners and others, 2004) and several other minerals, including titanite (Reiners and Farley, 1999), goethite (Shuster and others, 2005), magnetite (Blackburn and others, 2007), and apatite (Farley, 2000). It has been attributed to inhomogeneous He distributions due to zonation, alpha ejection, or initial diffusion rounding; surface roughness; multiple diffusion domains; and radiation damage (for example, Reiners, 2005).

Despite the myriad possible causes, samples with concave-up curves that lie above the linear trends (RB140, BR231, and M127) are still difficult to explain. Anisotropy could produce non-Arrhenius behavior of this type in step-heating results if activation energies of the contrasting diffusion directions are different (Reich and others, 2007; Watson and others, 2010; Bengston and others, 2012). However, our data, as well as all other experimental He diffusion data on zircon or zircon-structure phases (Farley, 2007; Cherniak and others, 2009), indicate that anisotropy is manifest as differences in frequency factor, not activation energies, which would not lead to significant departure from linearity on an Arrhenius plot. Furthermore, if anisotropy in zircon was caused by differences in activation energy, then large changes in slope should always be evident at low temperatures regardless of how many prograde and retrograde cycles have been performed. For the most part, we do not observe this. We therefore rule out anisotropy as the source of non-Arrhenius behavior in our initial prograde temperature steps.

Reiners and others (2002) suggested that this type of non-Arrhenius behavior could be due to the interaction between radiation damage zones and the zircon surface, with the damage zones acting as grain boundaries or fast diffusion pathways. Recent evidence of two diffusion pathways for Ar in quartz (Clay and others, 2010), as well as investigations of fast path diffusion in other minerals (Yund and others, 1981; Yund and others, 1989; Yurimoto and others, 1989; Worden and others, 1990; Hacker and Christie, 1991) are consistent with this. For example, because lattice diffusion generally has a higher E_a than grain boundary diffusion (Chakraborty, 2008), lattice diffusion can be faster at high temperatures, so diffusion from grain-boundary-like domains could dominate gas release in early, prograde steps, yielding initially high diffusivities in Arrhenius plots. With higher temperature steps increasing fractions of gas would derive from lattice diffusion as a result of the difference in values of E_a, but also because surficial grain-boundary-like sites would be rapidly depleted in the initial few tenths to one percent of gas released (fig. 4). If grain-boundary-like sites were reoccupied with He during high temperature steps, this could also explain persistent release of gas via grain boundary diffusion in later steps. Arrhenius trends for M127, Mud Tank, RB140, and BR231 appear to exhibit slight curvature at the lowest temperatures of the post-high temperature heating steps (after initial 500 °C step), consistent with this explanation.

Concave-up curves that lie below the linear Arrhenius trends (Mud Tank, G3, and N17) probably result from a combination of the behavior described above and initially rounded He profiles. N17, with a closure temperature below 0 °C, likely possessed a rounded concentration profile prior to being step-heated, and, despite our selection of interior parts of the Mud Tank and G3 zircons, some portion from the rounded diffusion profile of these two samples was included as well. Although the initial He released from these zircons does not appear to conform to ideal expectations of simple Arrhenius behavior, its effect after the first few percent of gas release (except for N17) is negligible.

Functional Form for Damage-Diffusivity Relationship

For the remainder of our discussion, we directly relate alpha dose to structural damage, and derive a mathematical parameterization of the damage-diffusivity relationship that fits the data in figure 7. This requires a number of assumptions. As previously stated, we assume that alpha doses calculated from He ages sufficiently reflect the total accumulation of alpha-decay events since structural damage was last annealed. Because most of our diffusion samples have thermal histories involving a phase of relatively rapid cooling, and because we have direct and independent measurements of a proxy for structural damage (FWHM) for a subset of them, this is a reasonable assumption.

We must also assume that He diffusion in zircon scales with both alpha dose and structural damage in the same way. This is an important consideration as alpha dose is a calculated damage proxy and not a direct measure of damage (we detail below the reasons for using alpha dose). We again rely on observations from the apatite He system to support this assumption. Shuster and Farley (2009) showed that the amount of ionizing kinetic energy released into a sample (kerma) and the fission-track density had similar effects on He diffusivity. They demonstrated that He diffusivity in apatite changes systematically with both increasing and decreasing kerma (through annealing) caused by either artificial irradiation or natural, self-irradiation as monitored by fission track density. Flowers and others (2009) expanded on this and derived a term, effective spontaneous fission track density, that directly related alpha dose, accumulation and annealing of fission tracks, and He diffusivity. Unfortunately, experimental observations analogous to those of Shuster and Farley (2009) are currently lacking for zircon. But both of these studies suggest that, at least to first order, fission track density (which is itself a measure of one type of structural radiation damage) scales with “effective alpha dose,” and can therefore be related to bulk He diffusivity.

Finally, we are assuming that the negative correlation between alpha dose and diffusivity at low doses is not due to the effects of He (or Pb) concentration on He diffusivity. Shuster and Farley (2009) clearly showed that radiation damage, not He concentration, controls diffusivity in apatite. But as stated above, the types of experiments performed by Shuster and Farley (2009) do not yet exist for zircon, so we cannot completely rule out a concentration-dependency for He diffusion. However, the data of Shuster and Farley (2009) provide some confidence that similar processes are occurring in zircon, and this could be confirmed in the future by applying their methodology to zircon samples.

Although these assumptions are required, we suggest that linking alpha doses (preferentially those calculated using FT or He ages) to radiation damage is an appropriate choice for our purposes as these assumed “effective alpha doses” provide a straight-forward method for calculating damage accumulation through time. Furthermore, as we detail below, it provides a crucial link between equations describing He diffusivity and those describing damage annealing. The alpha dose is therefore the best measurement for integrating He diffusion, damage accumulation, and damage annealing over geologic timescales into an easily accessible thermochronometric modeling tool, which is our primary objective in this section. With this in mind, we first derive a parameterization for the damage-diffusivity relationship and then integrate this parameterization into a He diffusion and damage annealing numerical model.

A mathematical description of the damage-diffusivity relationship must account for both the initial decrease and ultimate increase in diffusivity across the spectrum of damage from at least Mud Tank to N17 (∼10¹⁶-10¹⁹ α/g). These two samples represent the damage range encountered for the vast majority of zircon crystals sampled at or near the Earth's surface. We desire a functional form consistent with the hypothesis that decreasing diffusivity at low damage is caused by accumulation of isolated damage zones that block crystallographically preferred He transport pathways and increase the tortuosity of He migration. At high alpha doses, increasing diffusivity would be due to decreasing effective domain size of undamaged zircon volumes, which are increasingly separated by interconnected fast-diffusing damage zones.

In order to explain decreasing diffusivity at low damage, we introduce an effective diffusivity (D_e) that represents He diffusion in a damage-free lattice modified by increasing tortuosity. Increasing radiation damage blocks or constricts easy He migration paths (for example, c-axis parallel channels) forcing He to take a more tortuous path out of the zircon. This is analogous to diffusion in a porous medium where effective diffusivity is expressed as (Cussler, 1984): D_z is the diffusion coefficient within the pores, or in our case, diffusion along c-axis pipes in a pristine zircon, and τ is the tortuosity. We represent D_z with the diffusion kinetics from a minimally damaged zircon. For this zircon's D₀, we fit the ORT_C slabs in figure 5 with a power-law relationship that yields an equation of y = 134.89*x^−1.578, where y is D₀ and x is dose. Projecting this relationship down to 1 × 10¹⁴ α/g gives a D₀ of 193188 cm²/s (all constants and their values described in the remaining text are listed in table 5). For the E_a, we average the activation energies from the samples with minimal amorphous fractions (all samples excluding G3 and N17) and the previously published results. Although these parameters are both extrapolations, the following equations could be easily modified to account for future diffusion data from less damaged or more appropriate zircon specimens.

Tortuosity τ could be represented in a variety of ways, and the atomic-scale processes by which radiation damage may affect migration pathway dimensions and diffusivity are complex. Damage may displace atoms into open channels that, depending on which species is displaced, could cause variable decreases in porosity. C-axis channels of initially high ionic porosity could be completely blocked and a larger fraction of the He migration path would be forced to occur orthogonal to the c-axis. Although the exact geometry or porosity of damage zones is hard to constrain, we predict that most zones should act as He barriers and τ should increase as the chance increases for diffusing He atoms to encounter these barriers. To model this behavior, we use a metric introduced by Ketcham and others (2013) for characterizing the undamaged portion of the lattice, mean intercept length, l_int, which is the average distance a particle can travel in a single direction without encountering a damage zone. The expression for τ relates the calculated l_int in a given zircon to the mean intercept length in our extrapolated, minimally damaged zircon, which also displays high diffusivity (l_{int 0}): The right-hand side of equation (2) is squared in part to improve the fit to our diffusion data (see below), but tortuousity is often mathematically expressed as the square of pore spacing or geometry (for example, Epstein, 1989). Ketcham and others (2013) have derived an empirical relation for l_int by modeling the accumulation and percolation of chains of connected, capsule-shaped alpha recoil tracks. They express l_int as a function of fraction amorphous (f_a): where SV corresponds to the surface to volume ratio of the capsules (1.669 nm⁻¹). In turn, f_a is described using the direct impact model (Gibbons, 1972): where B_a is the mass of amorphous material produced per alpha decay (5.48 × 10⁻¹⁹ g/α-event), and α is the alpha dose.

An explanation for the increase in diffusivity at high damage requires a derivation that accounts for interconnected amorphous zones. We hypothesize that at sufficiently high self-irradiation levels, amorphous zones caused by damage become interconnected and connected to the grain's surface. This represents an important shift from damage zones acting as barriers to damage zones acting as fast paths. The processes by which this might occur are not entirely apparent. Various macroscopic and long-range order properties change at high damage and track with amorphous fraction (for example, Ewing and others, 2003) and we expect that diffusivity should as well. However, several different atomistic processes may cause these changes. Devanathan and others (2006) showed that damage zones are characterized by an amorphous core surrounded by a high interstitial density rind. These amorphous cores could become connected at high damage and may form fast paths while the rinds cause increasing tortuosity at lower damage. Interconnected fission tracks, which only reach a percolation threshold at high damage (Ketcham and others, 2013) are another candidate for creating these fast paths (see below). Regardless of the exact mechanism, the net effect of this process is creation of an increasing number of progressively smaller, undamaged zones that are increasingly isolated from one another by increasingly interconnected and progressively larger damage zones with much higher diffusivity. These two effects can be accounted for by 1) modeling a decrease in the size of the diffusion domain of the undamaged portion of the grain, and 2) describing the bulk diffusivity as a harmonic average (as appropriate for an average of rates) in the undamaged and damaged portions of the grain. In detail, if the grain had a heterogeneous spatial distribution of U and Th, increasing damage could conceivably produce an apparent spectrum of diffusion domain sizes that might manifest itself in step-heating data as decreasing diffusivity or concave-up Arrhenius trends.

We first parameterize our effective diffusivity using an harmonic average and to do this, we recast equation (1) as: where D_N17 corresponds to the diffusivity of amorphous N17, and f′_c represents fraction crystalline and is equal to 1 − f′_a. We again use the direct impact model to describe f′_a and f′_c, but to improve the fit to our diffusion data we include an additional term (φ) within the exponential: A value greater than 1 for φ causes f′_a to increase more rapidly at lower alpha doses. As we demonstrate below, a value of 3 for φ is necessary for a proper fit to the data, however, we currently have no physical explanation for why this term should be required. It is possible that, while the direct impact model may adequately describe the buildup of amorphous zones in zircon, it does not fully account for the interconnection of certain parts of the zones (for example high-vacancy damage core vs. damage rim), or for the contribution of an additional interacting effect such as accumulation of fission tracks (Ketcham and others, 2013).

In order to account for the reduction of the domain size of the undamaged portion of the grain, we modify equation (5) by dividing each D by domain size (a): In this set-up, a is equal to the initial grain size in an undamaged zircon, which effectively decreases as f′_a increases following equation (6).

The φ term in equation (6) and the scaling of domain size with f′_a in equation (7) are admittedly somewhat empirical and heuristic, respectively, and their relationship to actual physical phenomena are tenuous. Recent work by Ketcham and others (2013) may offer some additional insight into this issue. These authors suggest that fission track interconnection may play an important, but until recently, underappreciated role in creating the macroscopic and long-range order qualities typically attributed to f_a. Furthermore, they derive an equation for a term that seems to more accurately reflect the effects of radiation damage on decreasing effective domain size: mean distance to nearest fission track d_n. This number decreases with increasing fission track percolation and replacing (a*f′_c)² with d_n in equation (7) results in a similar functional shape. Unfortunately, doing so does not provide a better fit to the real data. The current calculations do not account for whether the nearest fission track is part of a network connected to the outside of the grain, and improving the model in this respect may result in a better fit. For our present discussion, though, we proceed with equation (7) as it is currently derived.

Our combined equation for effective diffusivity is: where l_{int 0} is equal to 45920 nm (the value of l_int calculated from equation (3) at an alpha dose of 1 × 10¹⁴ α/g). Figure 10A shows equation (8) with our step-heating results. Our parameterization adequately captures the decrease in diffusivity by ∼3 orders of magnitude at low damage and the subsequent increase in diffusivity by ∼11 orders of magnitude at high damage. We also show a comparison between the T_c data from figure 8 and T_c values calculated using equation (8) and an initial grain radius of 60 microns (fig. 10B). To obtain effective E_a and D₀ values for this curve, we use a method that relies on pseudo-Arrhenius trends calculated at discrete doses with equation (8). These trends provide the kinetic parameters necessary to then calculate the T_c at the corresponding dose.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 10.

(A) Diffusivity as a function of alpha dose calculated at temperatures that roughly bracket the nominal zircon PRZ. Data points are the diffusivities of the various samples included in figure 7 (calculated with a domain size of 100 microns) and curves represent effective diffusivity as defined by equation (8). The various parameters used for this plot and their values are listed in table 5. (B) Closure temperature as a function of alpha dose calculated in a manner similar to figure 8. Gray points in both plots represent single sample that does not fall along our observed trend.

Although diffusional anisotropy may manifest itself in our model through the importance of increasing tortuosity at low damage, our model does not directly account for anisotropic diffusion. As figure 6 suggests, the effect of anisotropy on He diffusivity is minimal compared to the effect of radiation damage and we have not focused on it in this section. However, our derivation makes several predictions for a relationship between anisotropy and radiation damage. With increasing damage, c-axis parallel channels might be expected to become increasingly blocked and, correspondingly, anisotropy should decrease (Farley, 2007). Specifically, Mud Tank—our least damaged, oriented zircon—should be more anisotropic than all other samples. Figure 5 shows a large disparity between the D₀ values for the two Mud Tank slabs, and the two trends (constant for PAR_C slabs, decreasing for ORT_C slabs) seem to converge at sample M127. This suggests that, at least in terms of D₀, anisotropic differences decrease with increasing damage. However, the results in figure 6 complicate this observation. At high temperatures (500 °C), the difference in diffusivity between the Mud Tank slabs is ∼2 orders of magnitude, but at lower temperatures the degree of anisotropy between these two slabs is comparable to the anisotropy between slabs in the “moderate to high damage” group (∼1 order of magnitude difference, note that the Mud Tank lines have different slopes in figure 6). This suggests that anisotropy is relatively invariant across the damage spectrum from Mud Tank to M127, which seems to contradict our hypothesis that damage decreases the degree of anisotropy. Perhaps the anisotropy in Mud Tank is similar to the anisotropy in M127 because Mud Tank has a large number of other defects that also contribute to decreasing anisotropy. But why then would these defects not also lower diffusivity? We do not have a satisfactory answer for this yet, but point defects may have a relatively small effect on diffusivity compared to radiation damage, due to damage's more chaotic nature. Actual characterization of this difference requires more work. For now we focus our remaining discussion on using equation (8) to forward model date-eU correlations.

Implementation of Damage-Diffusivity Parameterization

As a demonstration of the utility of our parameterization, we link equation (8) to another equation describing damage annealing in zircon and forward model date-eU correlations in a manner similar to the RDAAM (Flowers and others, 2009). Model inputs consist of an effective diffusion domain lengthscale (assumed to be radius of a sphere with an equivalent surface-area to volume ratio as the grain), and U and Th concentrations for each grain, and a discretized t-T history for the entire dataset. With these, our model calculates the total He production, damage accumulation, damage annealing, He diffusivity, and He loss at each time step.

Damage accumulation and annealing are quantified with a series of equations similar to those described by Flowers and others (2009), and rely upon the kinetics of fission track annealing in zircon. Although we treat alpha damage as the primary factor in creating tortuosity and interconnections, using a fission track annealing model means we must assume that alpha damage anneals in a fashion similar to fission tracks. This is a potential problem in the apatite system as well, but the RDAAM's ability to predict date-eU correlations from reasonable t-T histories (for example, Flowers and others, 2007; Flowers and Kelley, 2011) shows that, to first order, modeling damage annealing with fission track kinetics is valid. What works well in apatite, though, may not be suitable for zircon. Especially problematic is the well-documented observation that alpha damage annealing in zircon occurs via two disparate processes—epitaxial recrystallization and ZrO₂ nano-crystal formation—at different temperatures and at different initial damage concentrations (Meldrum and others, 1998; Capitani and others, 2000; Zhang and others, 2000; Nasdala and others, 2002; Zhang and others, 2010). Furthermore, Garver and others (2005) suggested that zircon fission track annealing processes are most likely damage-dependent as well. A possible solution would be to use an alpha damage annealing model that accounts for multiple annealing processes, but despite an extensive literature on alpha damage accumulation and annealing in zircon, to the best of our knowledge, no kinetic model has been parameterized to describe alpha damage annealing at time scales beyond hours or days. If such a model is developed, equation (8) could be easily coupled to it due to our modular model design. More importantly, our main objective here is not to contrast and compare various annealing models. Rather, we simply desire a quantitative approximation of damage annealing kinetics in zircon that seems reasonable for geologic time scales. Given the currently available annealing models, a fission track model best satisfies this requirement.

In order to combine the total alpha damage produced during a series of discrete time steps with a fission track annealing model, we introduce α_e or equivalent alpha dose (α/g): where the initial factor serves to convert from nmol to decays, α is the number of alpha decays (in nmol) per gram produced in each time step, and ρ_r is reduced (normalized) fission-track density of fission tracks that formed during that time step. The alpha dose from each time step (t₂, t₁) is calculated as: where [²³⁸U], et cetera are in nmol/g. Our derivation of ρ_r starts with length reduction r, for which we use a simplified version of the fanning curvilinear fit (for example, Ketcham and others, 2007) to the ZFT annealing data of Yamada and others (2007): where β = −0.05721, C₀ = 6.24534, C₁ = −0.11977, C₂ = −314.937, and C₃ = −14.2868. To convert from reduced length to reduced density (ρ_r), we use the relation based on data reported by Tagami and others (1990), which begins at one and is truncated at a ρ/ρ₀ value of 0.36, below which there are no data:

Once calculated, values of α_e are linked to a diffusion model via equation (8). With an estimate of the diffusion coefficient from equation (8), we then solve the diffusion equation numerically for a spherical geometry with the Crank-Nicholson finite difference scheme used in the thermal modeling software package HeFTy (Ketcham, 2005). Although this is not the best representation of He diffusion in zircon (it implies isotropic diffusion), a spherical model captures the first-order features of the damage-diffusivity relationship. Future versions could incorporate the cylindrical finite element scheme of Watson and others (2010), which accounts for the effects of anisotropy. Alternatively, a spherical calculation scheme can be employed using the “active radius” method introduced by Gautheron and Tassan-Got (2010), who described how to determine the radius of an isotropic sphere that replicates diffusive loss from a prism with a given aspect ratio and diffusive anisotropy.

Our model demonstration consists of 6 different thermal histories, each designed to capture aspects of the date-eU correlations in figures 1 and 2. As inputs, we use eU ranging from ∼50 to ∼5000 ppm, grain radii of 60 μm, and t-T paths that begin at 600 Ma and end at the present. For comparison, we also model the single zircon He date that results from using the kinetics of Reiners and others (2004) with each of our thermal histories (black diamonds in fig. 11B). These t-T paths encapsulate 6 representative scenarios: 1) slow, monotonic cooling from 600 Ma to the present, 2) early cooling followed by a pulse of early-stage reheating, 3) early cooling followed by a pulse of late-stage reheating, 4) early cooling followed by prolonged time spent in the PRZ and then subsequent late cooling, 5) long term early heating and subsequent late-stage cooling, and 6) early cooling and prolonged exposure to low temperatures (fig. 11). Model outputs that result from scenarios 5 and 6 show mostly flat date-eU correlation, which are typical of most zircon He datasets. These two scenarios also result in He dates that are almost identical to dates obtained using the kinetics of Reiners and others (2004). For thermal histories with a single, rapid pulse of cooling from high temperatures, our new model does not alter interpretations of zircon He datasets made with previously published kinetics as this style of cooling will most likely not result in a date-eU correlation. Interestingly though, scenario 6 shows a steep negative correlation at high damage, despite having never been reheated above 20 °C post-500 Ma. These zircons have entered the PRZ without changing temperature and we note that a similar process may have occurred in the Sri Lankan dataset to produce the young dates for zircons K6 and N17.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 11.

Representative thermal histories (A) with corresponding forward modeled date-eU correlations (B). Each numbered t-T path in (A) corresponds to the similarly numbered correlation in (B). Thermal history number 6 is partially obscured by 2 and 3 because it follows the same path but lacks a post-500 Ma reheating event. For reference, we also plot the zircon He dates that result from each thermal history in (A) if the kinetics of Reiners and others (2004) are used. These are represented by the numbered black diamonds on the y-axis of (B).

Scenarios 1, 2, 3, and 4 demonstrate the various thermal histories that may produce positive, negative, or both types of correlations in the same sample. The thermal history for scenario 1 is characterized by slow cooling through the PRZ and results in a broad and relatively confined positive date-eU correlation (dates increase from 161-210 Ma). Damage in-growth and He diffusion occur simultaneously over a prolonged time span, which causes the damage amounts for the zircons in this scenario to be relatively similar while He is diffusing. Given high enough eU concentrations, the t-T path for scenario 1 also produces negative correlations in the same sample. Although this particular negative correlation is subtle, an older initial age for the start of slow cooling could produce a more pronounced negative correlation, as we further demonstrate with our Minnesota dataset below. Scenario 4 is somewhat similar to scenario 1 as both samples spend prolonged time periods in the PRZ. However, for most of their history, the zircons in scenario 4 are held at a lower temperature (180 °C) than those in scenario 1. At this temperature, damage in-growth outpaces annealing and the thermal history produces a negative correlation at relatively low eU concentrations once the sample is finally cooled.

In contrast, relatively short-lived pulses of reheating and cooling result in more marked positive or negative correlations. The thermal histories for both scenarios 2 and 3 contain reheating events that take place after the zircons have been held at low temperatures (20 °C) for at least 100 my. During this time span, the zircons in each scenario have acquired large differences in damage. In scenario 2, this results in a positive correlation that spans a date range from 428 to 538 Ma, but also a negative correlation that drops to zero at the highest eU. Like scenario 6, the PRZ for zircons with eU in excess of ∼3500 ppm is at very low temperatures and they no longer retain He. In scenario 3, we have chosen both a later start time (100 Ma) and a lower maximum temperature for the reheating event than scenario 2 (130 °C as opposed to 180 °C). The resulting plateau of dates at ∼550 Ma followed by a steep negative correlation is somewhat similar to scenario 6, except the steep drop off in dates occurs over much lower eU concentrations (from ∼1500 to 2000 ppm). The late-stage reheating event also produces a small plateau of ∼50 Ma dates at eU concentrations greater than ∼2000 ppm. These dates correspond to the initiation of cooling and illustrate that, given a large span in radiation damage, multiple pulses of reheating and subsequent cooling could be recorded in zircons from the same hand sample.

As a final demonstration, we forward model the thermal history of two real datasets shown in figures 1 and 2. We use one of the Apennines datasets (AP54B) as a representative positive date-eU correlation, and the Minnesota dataset as a representative negative date-eU correlation. For each dataset, we present a realistic thermal history that could have produced the observed data and show the resulting model output in figure 12. In the Apennines example, we rely on the thermochronometric data from previous publications to guide our t-T path construction. Bernet and others (2001) obtained ZFT dates on the same grains shown in figures 1 and 12A, which yield a peak date of 20.7±3.6 Ma (these author's P1 population). This translates to a lag time of 6.8 my as the depositional age of this particular unit in the Marnoso-arenacea Formation is 13.9 Ma. Zattin and others (2002) conducted a detailed regional study of this Formation and found that fully and partially reset AFT dates from their hinterland samples (most internal to the thrust belt) suggest maximum burial temperatures of 120 to 125 °C and post-depositional exhumation between 4 and 6 Ma. If we consider these multiple t-T constraints as model inputs, then we can produce the black line in figure 12A, a positive date-eU correlation that is in good agreement with the real dataset. In figure 12A, we also show the results from a simpler thermal history (dotted lines in t-T history and date-eU correlation) in order to demonstrate that our preferred thermal history (black line) produces a distinctive date-eU trend that best fits the data. Our new model therefore tightly constrains the thermal history of a given sample and discriminates between potential t-T paths.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 12.

Forward modeled date-eU correlations matched to datasets shown in figures 1 and 2. (A) Positive date-eU correlation (detrital Apennines sample AP54B) with forward modeled correlation resulting from thermal history shown in the left-hand panel. This thermal history is characterized by a period of rapid exhumation in the source terrain (as constrained by ZFT dates from Bernet and others, 2001), deposition and burial in the Apennine foreland basin, and final exhumation to the surface in the latest Miocene. We used a grain radius of 48 microns for each model zircon, which is the average for this dataset. Dotted light gray lines in both panels represent an alternative thermal history and the corresponding date-eU correlation. (B) Negative date-eU correlation (igneous Minnesota River Valley sample) with forward model correlation resulting from the thermal history shown in the left-hand panel. This thermal history begins at the approximate end of tectonism associated with the Penokean Orogeny in this area and proceeds at a cooling rate of .06 °C/my until 1100 Ma. This time corresponds with the initiation of the Keweenaw Rift System (the Minnesota River Valley was situated near the rift shoulder) and we model this as an increase in the long-term cooling rate to ∼.17 °C/my. Again, we used the average grain radius from this dataset (52 microns) as an input for each model zircon. Dashed and dotted lines in both panels represent two alternative thermal histories and the corresponding date-eU correlations.

Our Minnesota dataset has fewer t-T constraints and our modeled t-T path is slightly more speculative. A negative date-eU correlation and an oldest date of ∼925 Ma though suggests that these zircons have experienced very slow cooling rates since the Proterozoic. With this in mind, we explored several possible slow cooling paths since the Penokean orogeny at 1870 to 1820 Ma, which is the most recent episode of regional metamorphism to affect this area (Holm and others, 1998). Our best fit to the data consists of an initial cooling event beginning at 1850 Ma and 250 °C that proceeds at a rate of .06 °C/my. In order to reproduce both the pseudo-plateau of dates at low eU concentrations and the steep negative correlation at high eU concentrations, we accelerate our cooling rate from .06 °C/my to ∼.17 °C/my at 1100 Ma, which matches the age of opening for the failed Keweenawan Rift system. Again, we include the results from a couple of simpler histories (dotted and dashed lines in fig. 12B) to show that our choice for the Minnesota thermal history is not arbitrary. Slow cooling rates with a single value yield date-eU correlations that match either the high eU or low eU trends in the real data, but not both. We find a good fit to the data only by changing the cooling rate at a specific time (in this case, 1100 Ma). Despite being somewhat speculative, our model t-T constraints for the Minnesota dataset agree with the regional geologic history, produce a reasonable fit to the data, and demonstrate that the Minnesota sample has experienced very slow cooling at low temperatures for the past 1.8 by.

Impact of eU Zonation on Zircon Date-eU Correlations

In the above models, we assumed that both our real and model zircons are homogenous in their U and Th (hence eU) concentrations. Typical zircons, though, usually possess some degree of U and Th zonation. Although previous studies have discussed the effects of parent zonation on apatite (U-Th)/He dates in depth (for example, Farley and others, 1996; Meesters and Dunai, 2002; Hourigan and others, 2005; Farley and others, 2011; Ault and Flowers, 2012; Gautheron and others, 2012), the effects of strong parent zonation in zircon may be significantly different because of the reversal in damage-diffusivity relationship with progressive damage accumulation. This means that, for some thermal histories, different parts of the same zircon grain may have extremely different behavior, both of which may be very different from that expected from a grain with equivalent bulk eU homogeneously distributed.

Farley and others (2011) detailed three ways in which heterogeneous eU in zoned apatite (lacking the damage-diffusivity reversal potential) can affect measured ages and interpreted t-T histories. If an homogenous alpha ejection correction factor (F_TH) is used for a zoned grain, then the resulting He date will be either too young or too old depending on where the majority of eU is concentrated (rim or core, respectively) (Farley and others, 1996; Hourigan and others, 2005). This issue can be dealt with by using a zoned alpha ejection correction (F_TZ) that accounts for alpha particle redistribution, which is an option in HeFTy. Zonation also affects He diffusivity by altering the He concentration profile, and, because He diffusion is damage dependent, by creating distinct domains with different diffusion kinetics inside the crystal. In the apatite He system, these three factors can contribute to He date scatter (Flowers and Kelley, 2011). Ault and Flowers (2012) suggested, however, that for typical apatites eU concentrations between zones do not vary by more than a factor of ∼2, and for typical thermal histories the resulting relative date difference between zoned and homogenous apatites with the same bulk eU is no more than ∼10 percent.

We expect these zonation issues to result in greater fractional date differences for the zircon He system. Order of magnitude differences in eU zonation are not uncommon in zircons (for example, fig. 13 in Reiners and others, 2004), and these can cause large discrepancies in both damage and He concentration. Furthermore, unlike apatite, He diffusivity in zircon will either decrease or increase depending on the eU concentration of a given zone and the specific t-T path of the zircon. The interplay amongst the alpha ejection correction factors, He concentration profiles, and damage profiles in zoned zircons with different thermal histories is therefore complicated and a detailed examination of real He datasets with zoned grains is beyond the scope of our current study. However, we discuss below the results from several HeFTy models to show how zonation affects differential damage accumulation within a zircon and sample date-eU correlations.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 13.

Time-temperature histories and corresponding date-eU correlations for unzoned zircons, zircons with high eU rims, and zircons with high eU cores. The numbering for each t-T path is similar to figure 11, except we have omitted number 6 and replaced number 5 with a different thermal history. Core eU concentrations are either enriched or depleted by a factor of 7 relative to the bulk concentration of the whole grain. For example, if the bulk concentration is 500 ppm, then the core concentration in the high eU core zircon is 3500 ppm and it is ∼71 ppm in the high eU rim zircon. Both the concentration in the rim and the radial position of the rim are determined by maximizing the zircon zonation factor (see text for details). All zircons have a radius of 60 microns and the core is composed of either the inner 20 or 40 microns of the grain. Large bold symbols connected by horizontal curves represent F_TH corrected dates, small bold symbols represent F_TZ corrected dates, and small transparent symbols represent uncorrected dates. Lightly colored vertical bars have been added to aid in connecting corresponding F_TH and F_TZ corrected dates. See text for details on uncorrected, F_TH, and F_TZ corrected zircons. Note that the scale for each y-axis is different.

These model simulations consist of five representative thermal histories and the resulting date-eU correlations for a suite of unzoned zircons and zircons with simple concentric zonation of high eU cores or high eU rims (fig. 13). Both the zoned and unzoned zircons in each plot have bulk eU concentrations ranging from 250 to 1250 ppm with radii of 60 microns. The zonation profiles for the zircons in our models have two main inputs, the core eU concentration and the core's radial position (rim concentration is set by the bulk eU). In order to provide some uniformity to our choice of variables, we use a zonation impact index (mass of eU for the whole grain divided by the mass difference between the core and rim), which provides a rough estimate of magnitude of the effect that eU zonation has on a zircon's He concentration and damage profile. For a given bulk and core eU concentration, this value reaches a maximum at a certain core radial position and rim eU concentration. In turn, this particular radial position and rim eU maximizes the difference in diffusion kinetics between rim and core that result from He concentration and radiation damage disparities. Because we want to show a worst-case scenario for each date-eU correlation, we have chosen the rim eU concentration and core radial position that correspond to the maximum zonation impact number. For core eU concentration, we pick values that differ from the bulk eU concentration by a factor of seven (a bulk eU of 500 ppm leads to a core eU of either 71.4 or 3500 ppm), which results in a core radial position that is 20 microns from the center (total radius of 60 microns) for the enriched core zircons, and a core radial position that is 40 microns from the center for the enriched rim zircons. In zircons with enriched rims, the ratio between eU concentrations in the rim and core is 9.52:1, while the same ratio is 1:9.1 for zircons with enriched cores. A final consideration is the choice of a F_TH or a F_TZ alpha ejection correction. For unzoned zircons, we model both the uncorrected date at a given eU and the F_TH corrected date (F_TH = F_TZ for the unzoned case). For each zoned zircon, we model the uncorrected date, the F_TZ corrected date using the “redistribution” option in HeFTy, and the F_TH corrected date. The F_TH for each zoned zircon is calculated by taking the ratio between the unzoned F_TH corrected date at equivalent eU and the unzoned uncorrected date at equivalent eU. For zoned grains, the F_TH date is equivalent to measuring a raw grain date on a zoned zircon and applying a naive alpha ejection correction assuming no parent zonation.

The model results for all corrected and uncorrected zircons are shown in figure 13. Thermal histories 1 through 4 are the same as in figure 11. We have omitted t-T paths 5 and 6 from figure 11 as both of these yield nearly flat date-eU correlations for zoned and unzoned zircons. Instead, we have added a new t-T path 5 to figure 13 that could be appropriate for zircons from Laramide uplifts of the US Rocky Mountains (for example, our Bighorn sample). In all t-T scenarios, the F_TH corrected zoned dates (large bold symbols connected by curves in fig. 13) for zircons with high eU rims are younger than their unzoned counterparts. Only some of this discrepancy is due to an improper alpha ejection correction, as the F_TZ corrected zoned dates (small bold symbols) are also younger than unzoned zircons with the same bulk eU. The other causes for younger dates—which in scenarios 3, 4 and 5, are the predominant ones—are the combined effects of a heterogeneous He concentration profile and radiation damage. Rims that are enriched in eU relative to the core cause an increase in effective bulk diffusivity (and younger dates) because more He is located near the grain boundary, and, at high eU concentrations, the rims become heavily damaged. This damage effect is particularly apparent in scenario 3, where the grains with the highest bulk eU concentrations have accumulated enough damage such that the rim acts as a diffusion fast path. If we model the high-eU-rim zircon with a bulk eU of 1250 ppm using the kinetics of Reiners and others (2004), then this scenario yields a date of 544 Ma (as opposed to 256 Ma using the kinetics presented here), which further suggests that radiation damage is the primary cause of these younger dates.

Despite being systematically younger, most of the high-eU-rim dates have a similar style of date-eU correlation as the unzoned dates. One exception is scenario 5, where a date-eU correlation that is monotonically negative in the unzoned case is positive at low bulk eU concentrations. Although the rims of the high-eU-rim zircons have accumulated high degrees of damage prior to cooling, the cores range from low (∼5.5 × 10¹⁶ α/g at 50 Ma) to moderate (∼1.1 × 10¹⁷ α/g at 50 Ma) amounts of damage. This difference in damage is enough to cause a more retentive core and lower bulk diffusivity in the 500 ppm bulk eU grain than the 250 ppm bulk eU grain, which in turn results in a positive date-eU correlation over this eU range.

To explain the date-eU correlations for the high-eU-core zircons, we must similarly consider the degree of damage in both the core and the rim and how those two damage domains combine to affect bulk diffusivity. In scenario 1, dates are either older or younger than the unzoned dates at the same eU, and the style of date-eU correlation changes from a positive correlation in the unzoned grains to a positive-negative-plateau correlation in the high-eU-core grains. The onset of a negative correlation has shifted as the zircon core accumulates a high degree of damage (and therefore has a high diffusivity) at relatively low bulk eU concentrations. Interestingly though, a plateau at the highest bulk eU concentrations suggests that the high diffusivity of the core is somewhat mitigated by the degree of damage in the zircon's rim. A similar phenomenon occurs in scenario 2, and is especially apparent in scenario 3 if we compare the high-eU-rim grains with the high-eU-core grains. In scenarios 1 and 2, the rims of the high-eU-core zircons at high bulk eU concentrations have accumulated a moderate amount of damage (roughly 4 × 10¹⁷ α/g prior to any thermal event) such that the rim decreases diffusivity and acts to retard He diffusing out of the crystal, similar to the effects of RDAAM on zoned apatites with eU enriched rims (Farley and others, 2011; Ault and Flowers, 2012). For the rims of the high-eU-core grains in scenario 3, the damage accumulation is more substantial, but the final reheating event occurs at a relatively low temperature. This temperature is low enough to cause these zircon rims to be more retentive than the corresponding rims in the high-eU-rim zircons (and thus yield older dates).

In contrast to scenarios 1, 2, and 3, the thermal histories in 4 and 5 result in damage higher than ∼5 × 10¹⁷ α/g, and relatively high diffusivity in both the cores and rims for most bulk eU concentrations. In terms of damage, this places almost all of the domains in these high-eU-core zircons to the high-damage side of the point of lowest diffusivity in figure 10A. The rims of the high-eU-core zircons with 250 ppm bulk eU concentrations are the only domains with damage lower than 5 × 10¹⁷ α/g (both are ∼3 × 10¹⁷ α/g prior to final cooling at 50 Ma) in these two scenarios. The net result in both scenarios is a date-eU correlation for the high-eU-core grains that is similar in style to the unzoned zircons, but with systematically younger dates. An exception to this is the positive correlation for bulk eU concentrations of 250 ppm and 500 ppm in scenario 5. Among the high-eU-core grains, the rim for the 250 ppm bulk eU grain has a damage amount of ∼2.9 × 10¹⁷ α/g at 50 Ma (the age prior to final cooling), while the rim for the 500 ppm bulk eU grain has a damage amount of ∼5.8 × 10¹⁷ α/g at 50 Ma. This difference in damage causes the 500 ppm bulk eU grain to have a more retentive rim and thus lower bulk diffusivity relative to the 250 ppm bulk eU grain, which in turn produces a slight positive correlation.

These five t-T scenarios illustrate that date-eU correlations for zoned zircons may be complex because of reasons already discussed in previous studies, but also because of the reversal in He diffusion behavior with progressive damage accumulation. In the absence of a priori knowledge of parent zonation patterns, the effects of parent zonation on He diffusivity may complicate thermochronologic interpretations of date-eU correlations. We again note that we have purposely chosen worst-case zoning scenarios in order to convey the full scope of this issue. Real zircons may possess a less extreme degree or pattern of zonation, which will mitigate some of the date dispersion observed in figure 13. Consistent styles but variable degrees of zonation (for example, enriched rim but with varying rim thickness or enrichment factor relative to the core) may be expected among zircons from some rock types. This would produce dispersion in a band of date-eU correlations between the unzoned case (black curves) and an extreme zoned case (blue or red curves, depending on zonation style) in figure 13. However, suites of zircons from some samples may not possess systematic zonation patterns like those in figure 13. Especially in detrital settings, one may expect to date zircons with a range of zonation patterns: some may have high eU cores, some may have high eU rims, and others might be unzoned. In this context, the curves in figure 13 should be interpreted as bounding constraints for the total range of zonation variability. For a given t-T path, and no consistent zonation style, real date-eU correlations could potentially plot anywhere within the black, blue, or red curves.

Despite the apparent severity of the zonation problem, however, parent zonation may be characterized from laser ablation depth profiles or other techniques, prior to bulk grain dating. These observations also point to the potential for exploiting parent zonation to provide date-eU trends within individual grains, for example with in-situ laser ablation dating (Vermeesch and others, 2012). Under certain conditions, zircon grains with zoned parent concentrations may behave similarly to crystals with multiple He diffusion domains. In-situ measurements and/or step-heating analyses could potentially be used to interrogate the distribution of He and dates among these domains, providing powerful constraints on thermal histories from single grains.