Diffusion of helium in zircon and apatite

Chemical Geology 268 (2009) 155–166

Contents lists available at ScienceDirect

Chemical Geology j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / c h e m g e o

Diffusion of helium in zircon and apatite D.J. Cherniak ⁎, E.B. Watson, J.B. Thomas Department of Earth and Environmental Sciences, Rensselaer Polytechnic Institute, Troy, NY 12180, USA

a r t i c l e

i n f o

Article history: Received 10 June 2009 Received in revised form 20 August 2009 Accepted 28 August 2009 Editor: R.L. Rudnick Keywords: Zircon Apatite Helium Diffusion Nuclear reaction analysis Thermochronology

a b s t r a c t Diffusion of helium has been characterized in natural zircon and apatite. Polished slabs of zircon and apatite, oriented either normal or parallel to c were implanted with 100 keV 3He at a dose of 5 × 1015 3 He/cm2. Diffusion experiments on implanted zircon and apatite were run in Pt capsules in 1-atm furnaces. 3He distributions following experiments were measured with Nuclear Reaction Analysis using the reaction 3He(d,p)4He. For diffusion in zircon we obtain the following Arrhenius relations: −7

D⊥c = 2:3 × 10

−5

D∥c = 1:7 × 10

expð−146 F 11 kJ mol expð−148 F 17 kJ mol

−1

−1

2 −1

= RTÞm s

2 −1

= RTÞm s

:

:

Although activation energies for diffusion normal and parallel to c are comparable, there is marked diffusional anisotropy, with diffusion parallel to c nearly 2 orders of magnitude faster than transport normal to c. These diffusivities bracket the range of values determined for He diffusion in zircon in bulk-release experiments, although the role of anisotropy could not be directly evaluated in those measurements. In apatite, the following Arrhenius relation was obtained over the temperature range of 148–449 °C for diffusion normal to c: −6

D = 2:10 × 10

expð−117 F 6 kJ mol

−1

2 −1

= RTÞm s

:

In contrast to zircon, apatite shows little evidence of anisotropy. He diffusivities obtained in this study fall about an order of magnitude lower than diffusivities measured through bulk release of He through step-heating, and within an order of magnitude of determinations where ion implantation was used to introduce helium and He distributions measured with elastic recoil detection. Since the diffusion of He in zircon exhibits such pronounced anisotropy, helium diffusional loss and closure cannot be modeled with simple spherical geometries and the assumption of isotropic diffusion. A finite-element code (CYLMOD) has recently been created to simulate diffusion in cylindrical geometry with differing radial and axial diffusion coefficients. We present some applications of the code in evaluating helium lost from zircon grains as a function of grain size and length to diameter ratios, and consider the effects of “shape anisotropy”, where diffusion is isotropic (as in the case of apatite) but shapes of crystal grains or fragments may depart significantly from spherical geometry. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Radiogenic helium has long been employed as a low-temperature chronometer. Given the relatively rapid diffusion rates of He and consequent low closure temperatures, ages determined from measurement of helium can provide a range of information about the

⁎ Corresponding author. Department of Earth and Environmental Sciences, Science Center 1W19, Rensselaer Polytechnic Institute, 110 8th St., Troy, NY 12180, USA. Tel.: +1 518 276 8827; fax: +1 518 276 2012. E-mail address: [email protected] (D.J. Cherniak). 0009-2541/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.chemgeo.2009.08.011

time–temperature evolution of rocks inaccessible by other means. Measured (U–Th)/He ages for apatite, zircon and titanite, when coupled with a quantitative understanding of the effects of temperature on measured ages due to diffusional loss and conditions for He retention in minerals, can be used to map thermal structure to estimate the timing, rate, and structural details of the process of exhumation and various near-surface tectonic and geomorphic processes (e.g., Zeitler et al., 1987; Wolf et al., 1996a,b; House et al., 1997, 1998; Wolf et al., 1998; Reiners and Farley, 1999; Reiners et al., 2000; Reiners and Farley, 2001). With advancement of analytical techniques these methods have found increasing use (e.g., Wolf et al., 1998; Farley, 2000, 2002; Reiners et al, 2002, 2004, 2005a; Shuster

156

D.J. Cherniak et al. / Chemical Geology 268 (2009) 155–166

et al., 2006; Shuster and Farley, 2009). Further understanding of He diffusion kinetics is clearly essential for improved characterization of conditions under which He is lost or retained in these mineral phases and in refining interpretations of time–temperature histories. Much of the present knowledge of diffusion behavior in the minerals apatite and zircon has been obtained from outgassing studies, where samples are heated in a stepwise manner and the He released is measured by mass spectrometry (e.g., Farley et al., 1996; Wolf et al., 1996a,b; House et al., 1997; Wolf et al., 1998; Reiners et al., 2000). While the application of these methods has provided invaluable information about the He diffusion systematics of apatite, zircon, titanite, REE phosphates and other minerals, several limitations exist. These “bulk-release” patterns do not provide information about initial distributions of helium. Helium may be nonuniformly distributed, depending on the distribution of U and Th in mineral grains, alpha ejection from mineral grains during radioactive decay, and differential He losses due to past thermal histories. In addition, microstructures within mineral grains may create diffusion fast paths, which cannot be readily distinguished from lattice diffusion because profiles are not measured directly in these analyses. Refinements of step-release techniques (e.g., Shuster et al., 2003; Shuster and Farley, 2004, 2005a,b) which involve the generation of a uniform 3He distribution by proton irradiation, coupled with measurement of the ratio of simultaneously released 4He/3He, have helped to surmount some of these difficulties and provide more detailed information about the 4He distributions within mineral grains. However, non-Arrhenian and other poorly understood diffusion behaviors have been observed for helium in some of these stepheating studies. These include departure from linear trends for apatite at temperatures above ~280 °C (e.g., Farley, 2000, 2002) and the appearance of anomalously high He diffusion in early stages of stepheating of zircon (Reiners et al., 2004). It has been suggested (e.g., Shuster et al., 2006; Shuster and Farley, 2009) that radiation-damaged apatite will be more retentive of He because of trapping of He in traps created through lattice disruption resulting from damage produced by alpha-recoil particles, fission fragments, or other species produced during radioactive decay. Further, anisotropy of diffusion may be pronounced in some mineral phases. Since the existence and extent of anisotropy may be difficult to evaluate through bulk-release techniques, direct profiling of He distributions to evaluate diffusion coefficients may provide additional understanding of diffusion behavior and processes. Direct profiling to characterize the diffusion of implanted 3He and 4 He in apatite and related ceramics, zirconia and other oxides, phosphates and silicates has been undertaken using the techniques of nuclear reaction analysis (NRA) and Elastic Recoil Detection (ERD) (e.g., Ouchani et al., 1998; Costantini et al., 2002; Gosset et al., 2002; Trocellier et al., 2003a,b; Miro et al., 2006, 2007). These investigations complement the findings from bulk-release experiments and provide additional constraints on He diffusion. A similar approach to measuring He diffusion in apatite and zircon is taken here, with ion implantation used to introduce 3He, and 3He concentrations measured by nuclear reaction analysis. In addition to applications in low-temperature geochronometry, the diffusion findings reported here may have significance for radioactive waste storage. Zircon, apatite and related ceramics have been considered as candidates for storage of various grades of nuclear waste. For example, apatite has been proposed as a primary matrix for actinide elements separated from high-level radioactive waste, or for use in engineered barriers for waste disposal in deep geological repositories (e.g., Weber et al., 1984), and zircon has been suggested as a material for storage of weapons-grade plutonium (e.g., Ewing et al., 1995; Weber et al., 1997). Since helium is produced in these materials through radioactive decay of actinides and lanthanides, it is important to characterize helium transport rates in order to better assess the long-term stability of these phases.

2. Experimental procedure 2.1. Sample preparation The apatite used in this study is a natural fluorapatite from Durango, Mexico. The chemical composition of the Durango fluorapatite, as reported by Young et al. (1969), is Ca9.83Na0.08Sr0.01RE0.09 (PO4)5.87(SO4)0.05(SiO4)0.06(AsO4)0.01(CO3)0.01F1.90Cl0.12(OH)0.01. This apatite has average concentrations of Th and U of 202 and 11 ppm, respectively, and has an age of 31.44 ± 0.18 Ma (McDowell et al., 2005). Among the rare earth elements, La, Ce, and Nd are each present in concentrations of a few thousand ppm. Two natural zircons are used in this study: one specimen from the Mud Tank carbonatite, and an optically clear, natural megacrystic zircon from an unspecified locality in Australia. The latter is referred to in this paper as DR Zircon. The Mud Tank zircon has 11 and 5 ppm U and Th, respectively (e.g., Watson and Cherniak, 1997), and has an age of 732 ± 5 Ma (Black and Gulson, 1978). The zircons were cut parallel or normal to the c-axis into slabs about 0.5 mm thick and polished to 0.3 µm alumina, followed by a chemical polish with colloidal silica. Specimens of apatite were also cut either normal or parallel to c, and polished using a procedure similar to that for the zircon. Following polishing, samples were cleaned ultrasonically in distilled water and ethanol. Polished samples were not heat-treated prior to implantation and heating in diffusion experiments. 2.2. Ion implantation and diffusion experiments Polished samples were mounted for ion implantation on an aluminum plate using carbon paint. The mounted samples were implanted with 100 keV 3He ions produced in the Extrion ion implanter at the Ion Beam Laboratory at the University at Albany. Implant doses were 5 × 1015 3He/cm2. The implanted zircon and apatite slabs were then cut into sections ~2 to 3 mm on a side, and heated in Pt capsules placed in Kanthalwound 1-atm furnaces for times ranging from several minutes to several weeks at temperatures from 148 to 604 °C. Temperatures in furnaces were monitored with chromel–alumel (type K) thermocouples, with temperature uncertainties ~± 2 °C. 2.3. NRA analysis Nuclear reaction analysis was carried out at the 4 MeV Dynamitron accelerator at the University at Albany. Samples were analyzed using the 3He(d,p)4He reaction (e.g., Pronko and Pronko, 1974; Dieumegard et al., 1979; Payne et al., 1989; Paszti, 1992). This reaction has been used to measure 3He in a variety of minerals and ceramic materials, including apatite (e.g., Miro et al., 2006), zirconia (Gosset et al., 2002; Costantini et al., 2003; Trocellier et al., 2003a,b) and britholite (Gosset and Trocellier, 2005; Gosset et al., 2002; Costantini et al., 2002; Trocellier et al., 2003a,b), and garnet (Roselieb et al., 2006). The protons produced in the reaction, along with backscattered deuterons and products of various (d,p) and (d,α) reactions induced with light elements comprising the sample, are detected with a solid state surface barrier detector with 1500 µm depletion depth positioned at 167.5° with respect to the incident beam. Because the protons produced in the 3He(d,p)4He reaction are so energetic, they stand apart from other contributions to the spectrum, with very little background (Fig. 1). The cross-section of the reaction has a maximum around 430 keV, but the peak is broad so it is not possible to obtain high depth resolution for 3He profiling using typical approaches for resonant or non-resonant NRA. For these analyses, we adapt the techniques used in other work (e.g., Costantini et al., 2002; Miro et al., 2006) by performing analyses over a range of energies to better define the profile, and comparing the proton yield from the annealed sample (i.e., a sample from a diffusion experiment) to an implanted, unannealed

D.J. Cherniak et al. / Chemical Geology 268 (2009) 155–166

Fig. 1. Example spectrum from a 3He-implanted apatite showing yields from the 3He(d, p)4He reaction. The high-energy (high channel number) signal is from protons produced in the reaction. The lower-energy signal is due to backscattered deuterium and products of various (d,p) and (d,α) reactions induced with light elements in the sample. Incident deuterium beam energy is 0.5 MeV.

sample at each energy step. For each sample, analyses at 4 or 5 different energies were performed for energies ranging from 0.5 to 1.0 MeV, using an energy step of 0.1 to 0.2 MeV. The ratio of total protons detected from the 3He(d,p)4He reaction for annealed samples to those measured for an implanted unannealed sample of the same material was calculated for each energy. Fig. 2 shows relative count yields from implanted zircons for various incident deuteron energies, comparing those from an implanted, unannealed sample to those from samples annealed at 450 °C for 4 and 14 h. These measured ratios were converted into diffusivities by first evaluating the fractional loss of diffusant from an implanted profile as a function of Dt. For a semi-infinite medium with the concentration of diffusant equal to zero at x = 0, the distribution of the implanted species can be described as a function of depth x and time t as (Ryssel and Ruge, 1986): Nm = 2 Nðx; tÞ = rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   2Dt 1 + ΔR 2 " + exp −

ð

0 pffiffiffiffiffiffi pffiffi 13 # 2 Rpffiffi4Dt 2ΔR + xpffiffiffiffiffiffi ðx−RÞ2 2ΔR 4Dt A5 4 @ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p × 1 + erf exp − 2ΔR2 + 4Dt 2ΔR2 + 4Dt "

0 pffiffiffiffiffiffi pffiffi # 2 Rpffiffi4Dt 2ΔR − xpffiffiffiffiffiffi ðx + RÞ2 2ΔR 4Dt 41 + erf @pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi × 2 2 2ΔR + 4Dt 2ΔR + 4Dt

Þ

31 5A

ð1Þ

157

where D is the diffusion coefficient, Nm is the maximum concentration of the implanted species (in an unannealed sample), R is the range (depth in the material) of implanted species, and ΔR is the range straggle (FWHM width of the initial implanted distribution). We obtain values of R and ΔR for 3He in zircon and fluorapatite from the Monte-Carlo simulation program SRIM 2006 (Ziegler and Biersack, 2006); these are 3800 and 830 Å, and 4800 and 920 Å, respectively; values for Nm and for zircon and apatite are ~ 0.30 and ~0.35 at.%. This expression describes the 3He distribution, but the proton yield measured will be a function not only of the 3He concentration but also of the depth in the material over which the incident deuteron beam can induce the nuclear reaction, the cross-section of the reaction as a function of energy (and depth in the material), the number of deuterons impinging on the target (Nd), and the solid angle subtended by the detector (Ω). The number of detected protons for a profile analyzed with a given beam energy E0 can be determined from the expression

∞

Np ðE0 Þ = Nd ðE0 ÞΩ∫x = 0

dσðEd ðxÞÞ ρðxÞdx dΩ

ð2Þ

where dσ(Ed(x))/dΩ is the differential cross-section for the 3He(d,p) 4 He reaction at energy Ed, which describes the probability of the reaction occurring at this deuteron energy. The deuteron energy Ed(x) will be reached at some particular depth x in the material, which will depend on the incident deuteron energy (E0) and the rate of energy loss for the deuterons with depth within the material. In the integral above, ρ(x) represents the depth distribution of the diffusant. The above integral is evaluated numerically for each incident deuteron energy by calculating the variation with depth of the reaction crosssection, using curves of cross-section as a function of energy derived from the values of Moller and Besenbacher (1980) and Mayer et al. (2005) corrected for beam-detector angle, with the depth scales determined by calculating energy loss with depth for incident deuterons in apatite and zircon targets using stopping powers obtained from the software SRIM-2006 (Ziegler and Biersack, 2006). The crosssection curves as a function of energy are then used to determine the proton yield from the 3He concentrations ρ(x) across the profile as a function of depth; these are summed to determine the proton yield for the entire profile for the value of incident proton energy E0. The asimplanted yields for unannealed samples are calculated and then can be ratioed to the proton yields calculated for profiles for annealed implanted samples with given values of Dt (derived from the equations above), resulting in a direct relationship at each incident beam energy between the ratio of proton yield for implanted untreated and annealed samples and a value of Dt. Uncertainties in diffusivities are determined from counting statistics from the detected proton signals and the variance of calculated diffusivities among incident beam energies for each sample. 3. Results

Fig. 2. Count yields of protons produced in the 3He(d,p)4He reaction for various incident deuteron energies, comparing spectra from a 3He implanted unannealed zircon to those from implanted zircons annealed at 450 °C for 4 and 14 h. The ratio of the count yields of the annealed and unannealed samples at each energy will be a function of the diffusion coefficient and heating time. See text for additional details.

3 He diffusion coefficients for apatite and zircon are presented in Tables 1 and 2 and plotted in Figs. 3 and 4. Over the temperature range of 148–449 °C, the data for apatite appear to define linear Arrhenius behavior. From a least-squares fit to the apatite data, we obtain an activation energy of 117 ± 6 kJ/mol and pre-exponential factor of 2.10 × 10− 6 m2/s (log Do = −5.678 ± 0.599) for diffusion normal to c. Diffusion parallel to c does not differ significantly; for these data we obtain an activation energy of 129 ± 12 kJ/mol and preexponential factor of 4.21 × 10− 5 m2/s (log Do = − 4.375 ± 1.136). For zircon, a fit to the data for the DR zircon for diffusion normal to c over the temperature range of 300–600 °C yields an activation energy of 154±13 kJ/mol and pre-exponential factor of 6.25×10− 7 m2/s

158

D.J. Cherniak et al. / Chemical Geology 268 (2009) 155–166

Table 1 3 He diffusion in apatite. Time (s)

D (m2/s)

log D

+/−

Diffusion normal to c He3Ap-6 449 He3Ap-5 401 He3Ap-4 349 He3Ap-2 300 He3Ap-13 302 He3Ap-14 302 He3Ap-3 250 He3Ap-7 200 He3Ap-8a 148

4.20 × 102 1.20 × 103 1.80 × 103 1.20 × 103 3.60 × 103 1.08 × 104 1.80 × 103 2.59 × 105 1.83 × 106

4.43 × 10− 15 2.53 × 10− 15 4.11 × 10− 16 4.56 × 10− 17 4.60 × 10− 17 3.46 × 10− 17 1.61 × 10− 17 1.54 × 10− 19 9.26 × 10− 21

− 14.35 − 14.60 − 15.39 − 16.34 − 16.34 − 16.46 − 16.79 − 18.81 − 20.03

0.24 0.32 0.26 0.30 0.17 0.36 0.35 0.29 0.47

Diffusion parallel to c He3Ap-10 402 He3Ap-9 354 He3Ap-12 301 He3Ap-8b 148

1.20 × 103 1.80 × 103 3.60 × 103 1.83 × 106

1.00 × 10− 14 2.22 × 10− 15 4.89 × 10− 17 1.07 × 10− 20

− 14.00 − 14.65 − 16.31 − 19.97

0.41 0.35 0.18 0.47

T (°C)

(log Do =−6.204±0.992). A fit to the DR and Mud Tank diffusion data for this orientation results in an activation energy of 146±11 kJ/mol and pre-exponential factor of 2.27 ×10− 7 m2/s (log Do =−6.645 ± 0.799). Diffusion parallel to c appears considerably faster (by nearly two orders of magnitude), but with similar activation energy for diffusion. From a fit to these data, we obtain an activation energy for diffusion 148± 17 kJ/mol and a pre-exponential factor of 1.67×10− 5 m2/s (log Do = −4.778±1.324). The time series (Fig. 5) for both apatite and zircon indicates similar diffusivities over times ranging over about an order of magnitude at 300 and 450 °C, respectively, suggesting that the dominant process being measured is volume diffusion, and that diffusion behavior is consistent with the assumptions made in the diffusion model.

Table 2 3 He diffusion in zircon. T (°C) Mud Tank zircon Diffusion normal to c He3MTZ-10 551 He3MTZ-7 499 He3MTZ-9 501 He3MTZ-5 449 He3MTZ-8 451 DR zircon Diffusion normal to c He3DRZ-7 604 He3DRZ-3 551 He3DRZ-2 501 He3DRZ-1 451 He3DRZ-9 450 He3DRZ-10 450 He3DRZ-8 398 He3DRZ-4 390 He3DRZ-5 350 He3DRZ-6 317 Diffusion parallel to c He3DRZ-12 549 He3DRZ-14 502 He3DRZ-21 500 He3DRZ-15 450 He3DRZ-18 449 He3DRZ-16 402 He3DRZ-17 401 He3DRZ-11 347 He3DRZ-19 348 He3DRZ-20 298

Time (s)

D (m2/s)

log D

+/−

1.80 × 103 9.00 × 102 2.70 × 103 9.00 × 102 4.50 × 103

5.53 × 10− 17 8.47 × 10− 17 2.19 × 10− 17 2.58 × 10− 17 1.60 × 10− 17

− 16.26 − 16.07 − 16.66 − 16.59 − 16.79

0.24 0.31 0.16 0.30 0.23

1.20 × 103 1.80 × 103 2.70 × 103 4.50 × 103 1.44 × 104 5.04 × 104 3.17 × 105 1.71 × 105 4.14 × 105 3.63 × 106

6.17 × 10− 16 1.56 × 10− 16 1.73 × 10− 17 1.05 × 10− 17 3.46 × 10− 18 3.61 × 10− 18 8.55 × 10− 19 2.81 × 10− 19 1.39 × 10− 19 1.18 × 10− 20

− 15.21 − 15.81 − 16.76 − 16.98 − 17.46 − 17.44 − 18.07 − 18.55 − 18.86 − 19.93

0.44 0.36 0.25 0.25 0.36 0.28 0.50 0.27 0.27 0.43

2.10 × 103 2.70 × 103 6.00 × 102 1.44 × 104 1.20 × 103 1.55 × 105 3.60 × 103 4.32 × 105 6.66 × 104 3.17 × 105

6.33 × 10− 15 1.44 × 10− 15 2.27 × 10− 15 2.31 × 10− 16 2.84 × 10− 16 4.25 × 10− 17 6.86 × 10− 17 2.28 × 10− 18 6.88 × 10− 18 3.95 × 10− 19

− 14.20 − 14.84 − 14.64 − 15.64 − 15.55 − 16.37 − 16.16 − 17.64 − 17.16 − 18.40

0.45 0.42 0.40 0.36 0.33 0.45 0.22 0.50 0.25 0.20

Fig. 3. Arrhenius plot for helium diffusion in apatite. From a fit to the data for diffusion normal to c (white circles) an activation energy of 117 ± 6 kJ/mol and pre-exponential factor of 2.10 × 10− 6 m2/s are obtained. Diffusion parallel to c (black circles) appears similar; for these data we obtain an activation energy of 129 ± 12 kJ/mol and preexponential factor of 4.21 × 10− 5 m2/s.

4. Discussion 4.1. Comparison with other diffusion data Farley (2000) measured 4He diffusion in Durango fluorapatite by stepped heating of apatite aliquots, obtaining an activation energy of 33 kcal/mol (138 kJ/mol) and pre-exponential factor of 1.5× 10− 4 m2/s for temperatures below ~300 °C. He observed a departure from linear Arrhenius behavior for temperatures above 265 °C, as had been reported previously (e.g., Zeitler et al., 1987; Wolf et al., 1996b). Shuster et al. (2003) have measured both 3He (produced through spallation by bombardment of materials with 150 MeV protons) and 4He diffusion in Durango fluorapatite, and obtained activation energies of 148 kJ/mol for both 3He and 4He, and Do/a2 values of 8.9 × 106 s− 1 and 7.3 × 106 s− 1 for 3 He and 4He, respectively. Apatite shards analyzed in this work were

Fig. 4. Arrhenius plot for helium diffusion in zircon. For diffusion normal to c we obtain an activation energy of 146 ± 11 kJ/mol and pre-exponential factor 2.27 × 10− 7 m2/s from a fit to the data for the DR (white circles) and Mud Tank (grey circles) zircons. Diffusion of He in zircon appears anisotropic, with diffusion parallel to c (black circles) considerably faster than diffusion normal to c, although activation energies for diffusion are similar. For diffusion parallel to c, the activation energy and pre-exponential factor determined from fits to the data are 148 ± 17 kJ/mol and 1.67 × 10− 5 m2/s, respectively.

D.J. Cherniak et al. / Chemical Geology 268 (2009) 155–166

Fig. 5. Time series for He diffusion in apatite (a) and zircon (b). He diffusivities are similar for times differing by up to a factor of 9 for apatite at 300 °C, and by times ranging over an order of magnitude for zircon at 450 °C, indicating that volume diffusion of He is the dominant process being measured.

sieved between 160 and 180 µm, so the values for the pre-exponential factor Do (taking the effective diffusion radius as 80 µm) would be 5.7 × 10− 2 m2/s and 4.7 × 10− 2 m2/s for 3He and 4He, respectively; given the stated uncertainties these are statistically indistinguishable. Shuster et al. (2003) speculate that this similarity of diffusivities of He isotopes, rather than a difference proportional to the square root of the ratios of isotope masses that might be expected from kinetic theory, results because He movement is limited by diffusion of crystallographic defects rather than intrinsic properties of the diffusing He atoms. Such observations may have important implications for mechanisms of noble gas migration in general and the relative behavior of He vs. other noble gases. Our diffusivities fall about an order of magnitude lower than data from both Farley (2000) and Shuster et al. (2003) (Fig. 6). Further, we see little evidence of a departure from linear Arrhenius behavior, as a fit to the data 300 °C and below yields Arrhenius parameters similar to those of the full data set up to 450 °C. Our results are consistent with the lack of anisotropy observed in these earlier studies. We also show data from the recent study of Shuster and Farley (2009) where apatite underwent neutron irradiation in order to evaluate the effects of radiation damage on helium retention. Arrhenius relations reported in this paper (using a value of a of 85 µm to calculate the pre-exponential factors Do from values for ln(Do/a2) given the reported apatite cross-section of 170 µm) for Durango fluorapatites with no neutron irradiation and those exposed to a 90 hour irradiation are plotted. These values are 139 and 161 kJ/mol for activation energies and 5.5 × 10− 3 and 1.66 m2/s for the pre-exponential factors, respectively. The Arrhenius trend for the unirradiated material plots near previous values from Farley (2000) and Shuster et al. (2003). The trend for the irradiated material plots slightly above this at higher temperatures (above ~200 °C). Also shown are Arrhenius relations for synthetic apatite from Shuster and Farley (2009) (with Do calculated from ln(Do/a2) values using the reported synthetic apatite crosssection of 300 µm). The synthetic apatites investigated by Shuster and Farley (2009) exhibit significantly slower helium diffusivities

159

Fig. 6. Summary of measurements of diffusion of He in apatite, comparing the findings from the present study with the ERD measurements of Ouchani et al. (1998) on ionimplanted fluorapatite, NRA measurements of Miro et al. (2006, 2007) on ionimplanted fluorapatite and Ca–Nd silicate apatite ceramics, and the step-heating experiments on natural and synthetic fluorapatite from Farley (2000), Shuster et al. (2003), and Shuster and Farley (2009). Our diffusivities are similar to those found by Ouchani et al. (1998), and approach those of Farley (2000) and Shuster et al. (2003) at low temperatures but deviate at higher temperatures because of differences in activation energy. The results for He diffusion in Durango fluorapatite from Shuster and Farley (2009) fall near the trends measured by Farley (2000) and Shuster et al. (2003), with the Arrhenius line for neutron-irradiated Durango fluorapatite from Shuster and Farley (2009) plotting slightly above these. Diffusivities measured in synthetic fluorapatite by Shuster and Farley (2009) are significantly slower than those measured for Durango fluorapatite, but these fall near the measurements of He diffusion in Durango fluorapatite by Ouchani et al. (1998) and somewhat below the findings reported in the present work. See text for additional discussion.

than Durango fluorapatite, whether or not they have been neutronirradiated. These diffusivities for synthetic apatite approach values for He diffusivities in other synthetic and natural apatites obtained by profiling diffusional broadening of depth distributions of ion-implanted helium using accelerator-based ion beam techniques. Ouchani et al. (1998) measured 4He diffusion in ion-implanted single-crystal Durango fluorapatite and synthetic Nd and Na bearing fluorapatites with 0.1% and 2% Nd (formulae Ca9.99Na0.01Nd0.01(PO4)6F2 and Ca9.8Na0.2 Nd0.2(PO4)6F2, respectively). Helium profiles were measured with elastic recoil detection analysis using an incident C beam. From the data for Durango fluorapatite, an activation energy of 120 kJ/mol and pre-exponential factor of 1.5 × 10− 6 m2/s are obtained. These values are very similar to our findings (Fig. 6). Diffusion of He in the synthetic (with Nd and Na) apatites appears similar to Ouchani's findings for Durango apatite, suggesting that He diffusion in apatite is not strongly affected by this range of compositional variation. This is also supported by the findings of Shuster et al. (2006) from a study of natural apatites from numerous localities. As in the other studies of apatite, little anisotropy for diffusion was noted. Miro et al. (2006) measured 3He diffusion in sintered fluorapatite ceramics. The 3He profiles from implanted annealed samples were measured by nuclear reaction analysis using the 3He(d,p)4He reaction. For polycrystalline fluorapatite, they obtained an activation energy of 122 kJ/mol (1.27 eV) and pre-exponential factor of 7.0 × 10− 7 m2/s (which yields a “normalized” value of 6.1 × 10− 7 m2/s when multiplied

160

D.J. Cherniak et al. / Chemical Geology 268 (2009) 155–166

by (4/3)− 1/2, their correction for isotopic effects to compare with data for 4He diffusion). However it should be noted that given the findings of Shuster et al. (2003) of the similarity of diffusion of 3He and 4He, the form of this mass correction, which is commonly applied for gases, may not be appropriate for diffusion in these materials. This is also supported by recent findings on isotopic fractionation of a range of elements in condensed phases (Richter et al., 2008, 2009). Arrhenius parameters reported for the sintered Nd silicate fluorapatite (formula Ca4Nd6(SiO4)6F2) are an activation energy of 86 kJ/mol (0.89 eV) and pre-exponential factor 6.7 × 10− 10 m2/s (resulting in a value of 5.8 × 10− 10 m2/s with the normalization factor of (4/3)− 1/2). While these reported diffusivities do not differ greatly from results for singlecrystal material, it is not clear how the relative contributions from grain boundary and lattice diffusion affect transport and measurements of diffusion in these polycrystalline materials. Grain sizes were reported to be ~ 20 µm, and beam spot sizes for analyses were about 0.5 mm, so spectra sampled multiple grains and grain boundaries in the material. He diffusion in zircon has been evaluated by Reiners et al. (2004) by step-heating of natural zircon grains from various sources. They obtained activation energies ranging from 163 to 174 kJ/mol and preexponential factors from 9.0 × 10− 6 to 1.5 × 10− 4 m2/s. Lines generated from these parameters yield diffusivities that bracket our diffusion data (Fig. 7); the variances in their data appear to reflect the anisotropy of diffusion and the differences in aspect ratio and other properties of the whole zircon grains used in the step-release studies. Reich et al. (2007) have simulated the diffusion of He in perfectly crystalline zircon using molecular dynamics simulations and empirical force field and quantum mechanical calculations. They find strong anisotropy, with diffusion in the [001] direction the most favorable pathway with a low activation energy for diffusion (13.4 kJ/mol). Their calculated diffusivities for this orientation, as well as the [100] direction, are considerably faster than those measured in the present study (Fig. 8). However, the authors point out that these simulations represent “tracer” diffusivities, i.e., a single atom diffusing through a perfectly crystalline lattice free of defects, which may differ significantly from volume diffusion coefficients representing an aggregate movement of atoms under a chemical potential gradient; such simulations therefore may not be predictive of bulk lattice diffusion rates. A more recent study by Saadoune et al. (2009) obtains similar results

Fig. 8. Comparison of Arrhenius relations from our experimental measurements of He diffusion in zircon with diffusivities calculated by Reich et al. (2007) from molecular dynamics simulations. Their calculated diffusivities for [001] and [100] are considerably faster than experimental determinations. However, these simulations represent “tracer” diffusivities, i.e., a single atom diffusing through a perfectly crystalline lattice free of defects, which may differ significantly from volume diffusion coefficients that represent an aggregate movement of atoms under a chemical potential gradient.

for “perfect” zircons, also predicting that diffusion along [001] will be most favorable for He transport, although there are differences among the results from these studies for migration along (100)/(010) directions. The calculations of Saadoune et al. (2009) suggest that the presence of point defects in zircon may affect He mobility, and note the possibility of defect trapping playing a role in He transport. However, these observations have not yet been quantitatively related to He diffusivities, and the activation energies calculated (~75–80 kJ/mol) are considerably lower than those measured in this study and that of Reiners et al. (2004). 5. Geological implications 5.1. He diffusion in accessory minerals

Fig. 7. Comparison of measurements of He diffusion in zircon parallel and normal to c (black and white circles, respectively) from the present study with those of Reiners et al. (2004) (dashed lines). The lines are calculated using the range in Arrhenius parameters reported by Reiners et al. (2004), i.e., activation energies 163–174 kJ/mol and preexponential factors 9.0 × 10− 6–1.5 × 10− 4 m2/s. These lines bracket our data and appear to reflect the effects of contributions from diffusion normal and parallel to c in the bulkrelease experiments of Reiners et al. (2004).

Fig. 9 shows a plot of He diffusion data for apatite and zircon (obtained from the present study) along with sphene (Shuster et al., 2003), monazite (Boyce et al., 2005), xenotime, and monazite and zircon-structure REE phosphates (Farley, 2007). Helium diffusion is faster in apatite than in zircon, while both sphene and CePO4 fall between the Arrhenius lines for diffusion parallel and normal to c in zircon. He diffusion in GdPO4 (like CePO4, a monazite-structure REE orthophosphate) is slower than in zircon, while He diffusion in zirconstructure REE phosphates and xenotime appears extremely fast (~5–7 orders of magnitude faster than He diffusion in apatite, and 7–10 orders of magnitude faster than in zircon). It is not clear why these differences exist. Although zircon and the medium to heavy rare earth phosphates are isostructural, the substitution of Zr and Si for REE and P, respectively, may influence available sites and diffusional pathways in such a manner to significantly influence He diffusion. Farley (2007) notes that natural materials tend to be more He-retentive than these synthetic REE phosphates and speculates that this may be a consequence of lattice disruption due to radiation damage in the natural specimens, or differences in defect densities between natural and the flux-grown synthetic minerals. However, it should be noted that in

D.J. Cherniak et al. / Chemical Geology 268 (2009) 155–166

Fig. 9. Comparison of He diffusion in apatite, zircon, sphene and REE orthophosphates. Sources for data: REE phosphates—Farley (2007); monazite—Boyce et al. (2005); sphene—Shuster et al. (2003); apatite and zircon—this study.

apatite the opposite has been observed, where the results of Shuster and Farley (2009) indicate that He retentivity is considerably higher in synthetic apatite than in natural specimens. The presence of flux inclusions or lattice-incorporated Pb from the Pb pyrophosphate flux used to grow the REE phosphates used in Farley's (2007) study may influence diffusional behavior as well. 5.2. Potential effects of radiation damage on He diffusion Helium diffusion may be influenced by radiation damage produced through decay of U, Th and their daughter products. Since both apatite and zircon incorporate these elements (which is obviously a key reason why these minerals are useful in He thermochronometry), they may experience these effects. In zircon, helium diffusion may be enhanced by alpha-recoil damage, when the damaged regions form a connected network of “fast paths” that permit more rapid transport (Nasdala et al., 2004; Reiners et al., 2005b). However, evidence suggests that radiation damage effects may not have a significant effect on He diffusion in zircon up to doses of ~4 × 1018 alpha decay events per gram (Reiners et al., 2004). Numerous studies indicate that radiation damage is annealed less readily and at higher temperatures in zircon than in apatite and other phosphates (e.g. Weber et al., 1997; Meldrum et al., 1997, 1999; Palenik and Ewing, 2002; Ewing et al., 2003; Palenik et al., 2003). It has been proposed (Shuster et al., 2006; Shuster and Farley, 2009; Flowers et al., 2009) that radiation damage (both natural and induced) in apatite may slow diffusion rates. Shuster et al. (2006) suggest that this is the result of defect trapping, with the authors drawing an analogy with this mechanism to the behavior of a gas bubble or fluid inclusion within a solid; models have been developed to attempt to quantify and better understand this behavior (Shuster et al., 2006; Flowers et al., 2009), incorporating the effects of accumulation and annealing of radiation damage. Although these models appear to be successful in reconciling discrepancies in interpretation of apatite (U–Th)/He thermochronometry from several localities, several outstanding questions remain. While point or other localized defects may provide traps for atoms that could limit diffusion rates, the extended defects typically created through radiation damage (predominantly alpha recoil) would seem

161

more likely to enhance, rather than inhibit, diffusion. It seems plausible in the case of apatite, however, that significant annealing of extended defects may occur at modest temperatures. However, it is unclear that there occurs any aggregation of He in extended defects, as would be needed to support the analogy used by Shuster et al. (2006) in suggesting that trapped He atoms behave like a gas in diffusing through the apatite structure. In addition, the significantly lower helium diffusivities measured by Shuster and Farley (2009) in synthetic apatite (with no radiation damage) compared with those measured for Durango and other natural fluorapatites with natural radiation damage seem difficult to reconcile with the simple conclusion that radiation damage generally tends to decrease diffusivities. It is important as well to distinguish between various types of radiation damage, both natural and induced, which may have significantly different effects on materials. In the present study, in which we have implanted light ions at relatively low energies, radiation damage is largely restricted to ionization rather than significant displacement of lattice atoms so radiation damage effects that might enhance diffusion should be minimal (e.g., Tesmer and Nastasi, 1995; Ziegler and Biersack, 2006). Although these are natural materials, the Mud Tank zircon (and likely the DR zircon, from a similar locality) has very low U and Th so has little natural radiation damage. The apatite has higher Th, but, as noted above, radiation damage in apatite is likely to be repaired at low temperature. The good agreement of our diffusivities with those from earlier work also supports the contention that radiation damage has little effect on measured He diffusivities in these experiments. In addition, as has been noted by Ewing and Wang (2002), the ionization by which alpha particles dissipate most of their energy causes enhanced annealing of alpha-recoil damage. This process, which is the dominant means of recrystallization recovery, explains the preservation of the crystalline state of natural apatites that contain actinide elements. This effect has also been demonstrated by Chaumont et al. (2002), who investigated the recovery from radiation-induced defects in fluorapatite irradiated with He and Pb, with the latter simulating the effects of alpha-recoil damage. 5.3. He closure temperatures From our diffusion data, we can calculate mean closure temperatures using the formalism of Dodson (1973, 1986). In Fig. 10 we plot mean closure temperatures for apatite and zircon as a function of effective diffusion radius for a cooling rate of 10 °C/Ma. For apatite, we plot curves using spherical geometry and Arrhenius parameters for diffusion both parallel and normal to c. For zircon, given the large anisotropy of diffusion, closure temperature curves are plotted using Arrhenius parameters for the “slow” direction (diffusion normal to c), with spherical geometry, and with planar geometry (the case of 1dimensional diffusion) for the “fast” direction. For comparison, we plot closure temperature estimates for Durango fluorapatite for given grain radii from Farley (2000) and Ouchani et al. (1998), for synthetic fluorapatite (Shuster and Farley, 2009), and for zircon (Reiners et al., 2004). There is reasonably good agreement for apatite (although the value for synthetic apatite plots higher than the others given its comparatively slow He diffusion), and for zircon there is overlap with the range of closure temperatures estimated by Reiners et al. (2004). The range of Tc values for zircon appears a consequence of contributions to diffusion from both axial and radial diffusion in the zircon, with contributions varying depending on the aspect ratio of the grains relative to diffusion rates. We will further consider implications of anisotropic diffusion in a later section. In Fig. 11, we plot curves for various values of closure temperature against cooling rates and grain radius for apatite, again using the equation of Dodson (1973) and diffusion parameters for He diffusion normal to c. For grain sizes and cooling rates typical in nature, closure temperatures for apatite will be in the range of 55–75 °C.

162

D.J. Cherniak et al. / Chemical Geology 268 (2009) 155–166

ages will be essentially zero. In the regions between, there will be partial retention of He, so He ages will vary with both temperature (depth) and time, and over time a steady-state condition will be achieved where radiogenic production and diffusive loss are balanced at each depth. Fig. 12 shows curves of apatite He ages in Ma for isothermal holding times of 50 Ma at the indicated temperatures for crystal radii of 150, 100 and 50 µm. The curves are generated using the expression of Wolf et al. (1998), assuming an initial He age of zero and the diffusion parameters measured in the present study. The helium PRZ is defined by Wolf et al. (1998) as conditions under which He ages fall between 10 and 90% of the holding time. For the parameters above, these values would be ~30–72 °C, for 100 μm grain radii, corresponding to depths between 1 and 3.1 km, assuming a geothermal gradient of 20 °C/km and surface temperature of 10 °C; values for 50 and 150 µm radii grains would be 21–61 °C and 35–80 °C, respectively. We can also perform similar calculations using the Arrhenius parameters of Farley (2000). Since diffusivities at low temperatures are similar for these two datasets, conditions for the PRZ are also similar, with a temperature range of 33–69 °C for 100 μm grain radii. 5.5. Implications of diffusive anisotropy of diffusion in zircon—Modeling with the CYLMOD code Fig. 10. Mean closure temperatures for apatite and zircon as a function of effective diffusion radius for a cooling rate of 10 °C/Ma, using the expression of Dodson (1973). For apatite, we plot curves using spherical geometry and Arrhenius parameters for diffusion both parallel and normal to c. For zircon, given the large anisotropy, closure temperature curves are plotted using Arrhenius parameters for diffusion normal to c with spherical geometry, and with planar geometry (the case of 1-dimensional diffusion) for diffusion parallel to c. For comparison, closure temperature estimates for Durango fluorapatite from Farley (2000) and Ouchani et al. (1998), for synthetic fluorapatite (Shuster and Farley, 2009), and for zircon (Reiners et al., 2004) for given grain radii are shown.

Another way to examine the expected diffusion behavior of helium in natural settings is to explore the variation of He ages with depth (and thus temperature). As temperature increases with depth, He ages will decrease, defining a helium partial retention zone (HePRZ). At the earth's surface and at low temperatures helium is retained and measured He ages will correlate with calendar time. At temperatures above ~100 °C, He is lost as quickly as it is produced, so measured He

Because of the anisotropy of He diffusion in zircon, simple calculations assuming spherical geometry and isotropic diffusion cannot adequately describe helium's diffusion behavior. In order to obtain a more thorough understanding, we can approach the problem by considering diffusion in a cylinder of finite length with differing diffusivities in the axial and radial directions. This problem has been solved analytically (Watson et al., submitted for publication), and a finite-difference code (CYLMOD) has been created to simulate anisotropic diffusion in this geometry. We present a few simple applications of the CYLMOD code, with isothermal heating, to evaluate the retentivity of He in zircons of various diameters and aspect ratios (here defined as the ratio of length of the zircon to its diameter). The simulations are performed for varying times at temperatures of 350 and 200 °C, with zircon radii of 10, 50 and 100 µm (diameters 20, 100, and 200 µm, respectively), and aspect ratios of 1:1, 2:1 and 5:1 for each zircon radius. In all cases, the caxis is the axial direction. Curves in Fig. 13 are shown for the case of 50 µm zircons, where the fraction of He remaining in the zircon is

Fig. 11. Curves for various values of closure temperature (in °C, noted on curves) as a function of cooling rate and grain radius for apatite, again using the equation of Dodson (1973) and diffusion parameters for He diffusion normal to c. For grain sizes and cooling rates typical in nature, closure temperatures for apatite will be in the range of 55–75 °C.

Fig. 12. Curves of apatite He ages (in Ma) for holding times of 50 Ma at the indicated temperatures for crystal radii of 50, 100 and 150 µm. The curves are generated using the expression of Wolf et al. (1998), assuming an initial He age of zero, and the diffusion parameters measured in the present study. The helium partial retention zone (PRZ) is defined by Wolf et al. (1998) as the conditions under which He ages fall between 10 and 90% of the holding time. For 100 µm radii grains, this zone is over the temperature range of 30–72 °C, corresponding to depths between 1 and 3.1 km, assuming a geothermal gradient of 20 °C/km and surface temperature of 10 °C. See text for further details.

5.4. Helium partial retention zones

D.J. Cherniak et al. / Chemical Geology 268 (2009) 155–166

Fig. 13. Fraction of helium retained as a function of time for various aspect ratios (ratio of grain length to diameter) modeled using the CLYMOD code for anisotropic diffusion in a cylinder. Calculations were performed for a temperature of 200 °C using the diffusion parameters for He diffusion in zircon normal and parallel to c.

plotted against the duration of isothermal heating at 200 °C. These curves confirm the intuitive conclusion that zircons will retain a higher fraction of He for a given time as the length along the faster-diffusing direction is increased. For example, 100 µm long zircons (1:1 aspect ratio) will lose 46% of their He in 10,000 years, while 200 and 500 µm zircons of the same radius would lose 27% and 15%, respectively. The curves for fraction of He remaining can also be normalized with respect to zircon radius. For a given aspect ratio, regardless of radius, the curves for the He fraction remaining will coincide if plotted as a function of time divided by the zircon radius squared. In addition, because diffusivities parallel and normal to c have parallel trends on an Arrhenius plot, they will always differ by the same relative amount regardless of temperature. Hence, we can also apply these findings to other temperatures by relating the time (t2) for a given amount of He retained at temperature T2 (with diffusivity D2) by simply multiplying the times (t1) at temperature T1 (with diffusivity D1) by the ratios of the diffusivities, viz. (D1/D2) ⁎ t1 = t2. Curves for the fraction of helium retained can also be calculated from analytical expressions presented in Watson et al. (submitted for publication). For the case of diffusion in cylindrical geometry with isothermal heating, we can define a “spherical equivalent” Ds/R2, where Ds is the effective diffusion coefficient for spherical geometry, with R the effective radius of the sphere. " rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi# Ds 4 D11 D33 D11 D33 = + + 2 9 a2 R2 h2 a2 h2

163

solution exists (Crank, 1975) relating the fraction of diffusant remaining to the dimensionless parameter Dt/a2. In order to compare the CYLMOD-determined values to those for a sphere, we “normalize” them by multiplying t/a2 for the 1:1 aspect ratio curve by a value of D to create a dimensionless parameter. In this case, we generate curves normalized by both the “fast” D (diffusion parallel to c) and the “slow” D (diffusion normal to c) at the temperature (200 °C) at which the curves (and therefore diffusion times for a given amount of He loss) were generated. These curves are plotted in Fig. 14. As can be seen in the figure, the two curves for the cylindrical model bracket that for the sphere, again a result that makes sense intuitively since the 1:1 cylinder approximates the sphere dimensionally but has contributions to total He loss from both the fast and slow diffusion directions. We can also generate an effective Ds using the expression above (Eq. (4)); the curve “normalized” to this value of D falls close to the trend for spherical geometry. Another instructive example is to calculate criteria for “center retention” for zircon grains of varying aspect ratios. In this case we use the criterion that preservation of the initial composition of grain centers corresponds to the concentration in the center node of the modeled zircon grain being at 99.9% of the original value; this is similar to the “center retention” criterion for a sphere. In the case of a sphere (isotropic diffusion) the value of Dt/a2 for preservation of initial concentration at its center is 0.03; in the case of a long cylinder (radial diffusion only) this value is 0.037 for preservation of initial concentration along the center axis. Using the CYLMOD code, we calculate values of t/a2 for center retention for various cylinder aspect ratios for the case of isothermal heating at 350 °C. The results of these simulations are plotted in Fig. 15. They are also compared with the limiting cases of diffusion in both a sphere and infinite cylinder, which are plotted using the criteria above (Dt/a2 equal to 0.03 and 0.037, respectively), calculating t/a2 for each geometry for both the “fast” (diffusion parallel to c) and “slow” (diffusion normal to c) diffusivities for He in zircon. The curve for the cylinder with anisotropic diffusion approaches the “fast diffusion” value of t/a2 for a 1:1 aspect ratio (when the aspect ratio for the cylinder is close to that of a sphere, and diffusion will be dominated by the faster transport in the axial direction); for large aspect ratios it approaches “slow diffusion” t/a2

ð3Þ

a and h are the radius and height of the cylinder, respectively, and D11 and D33 are the radial and axial diffusion coefficients (at a given temperature). Values of Ds/R2 for given aspect ratio, temperature and diffusion parameters can be determined from this expression. For the case of R = a = h/2 (cylinder of 1:1 aspect ratio and radius equal to that of the sphere), this expression becomes Ds =

  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 D D11 + 33 + D11 D33 9 4

ð4Þ

If we consider the case where axial diffusion is two orders of magnitude faster than radial diffusion (i.e., D33 = 100D11, similar to the difference for zircon) the effective diffusivity for a sphere of radius R = a is 16D11. Such an “averaged” diffusivity would be the value estimated for the diffusion coefficient in a material with anisotropic diffusion if spherical geometry and isotropic diffusion are assumed. We can directly compare He retention curves for a cylinder having a 1:1 aspect ratio with those for a sphere, for which an analytical

Fig. 14. Fraction of He retained as a function of the dimensionless parameter Dt/a2 for a sphere and for a cylinder with 1:1 aspect ratio (length = diameter). The dashed curve is calculated using analytical expressions for loss of diffusant from a sphere (Crank, 1975). The curves above and below were generated by determining times for loss of He at 200 °C using the CYLMOD code. The times were then transformed to a dimensionless parameter by multiplying by the diffusivities of He in zircon at 200 °C, using diffusion parallel to c for the “fast” curve and diffusion normal to c for the “slow” curve, and dividing through in both cases by the square of the grain radius. The dotted curve is calculated using the “effective diffusivity” Ds, derived from Eq. (4), to normalize. The “fast” and “slow” curves bracket the curve for spherical geometry, indicating an effective diffusivity comprised of contributions from both diffusion directions.

164

D.J. Cherniak et al. / Chemical Geology 268 (2009) 155–166

Fig. 15. Conditions for preservation of composition at center of cylindrical zircon grains of various aspect ratios, derived from modeling with the CYLMOD code. Times to meet this condition were determined at 350 °C for both anisotropic (circles) and isotropic (triangles) diffusion; in the latter case using both the “fast” (axial, parallel to c) (triangles down) and “slow” (radial, normal to c) (triangles up) diffusivities for He in zircon. Also plotted for comparison are t/a2 conditions for center retention for a sphere and infinite cylinder (radial diffusion only) at 350 °C for both “fast” and “slow” diffusion. As might be anticipated, conditions for center retention in a cylindrical zircon with equant dimensions approach the criterion for center retention in a sphere, but as aspect ratios increase axial diffusion becomes less dominant and conditions approach those for radial diffusion only. When diffusivities do not differ with orientation, the variation of center retention criteria with aspect ratio is significantly less pronounced.

value for an infinite cylinder. This is intuitively reasonable since diffusion in the axial direction will have less influence on diffusional loss as the axial dimension increases, so diffusion will be dominated by radial transport as cylinder length tends to infinity. Also shown for comparison are CYLMOD calculations for cylinders with isotropic diffusion, using both the “fast” and “slow” diffusivities at 350 °C; these bracket the curve for anisotropic diffusion. While this example is plotted for diffusion at 350 °C, we can also determine center retention times at other temperatures as described above by multiplying the times for diffusion at 350 °C by the ratio of diffusion coefficients at 350 °C and those at the temperature of interest. These are only a few examples of application of this code to simple modeling. More complex thermal histories can be simulated, and the straightforward extension of the relationships to different grain sizes and temperatures in the case of the parallel Arrhenius relations for He in zircon facilitates evaluation of a broad range of parameters. 5.6. Dependence of He release on aspect ratio for cases of isotropic diffusion It is sometimes the case that although diffusion in a particular mineral is isotropic, crystals of the mineral phase commonly have large length to diameter ratios. Apatite is a good example, since it appears that He diffusion is essentially isotropic, but its crystals are often elongate along the c-axis. This “shape anisotropy” has the potential to affect interpretation of He diffusional loss patterns if fragments or grains investigated depart significantly from equant dimensions, as the percentage of He released for given time–temperature conditions may differ from that for a simple spherical model. We can use the CYLMOD code to model He loss from apatite of various aspect ratios, and compare total fraction of He lost from these geometries with that lost from a sphere of equal radius to evaluate these effects. The calculations plotted in Fig. 16 are performed for isothermal heating at 300 °C, using the diffusion parameters for He in apatite reported in this paper, and an apatite radius (a) of 50 µm. Times (t) used in calculations are determined from

Fig. 16. Effect of “shape anisotropy” (differences in aspect ratio) where He diffusion is isotropic, using the example of diffusion in apatite. The curves correspond to different fractional amounts of He loss, plotting the differences between He fractional losses for spheres and those for cylinders with various aspect ratios (length/grain diameter). As aspect ratios increase, a greater fraction of He is retained in the cylindrical grains compared with the sphere because of the increasing diffusion distance in the axial direction. It can also be observed that for small fractions of He release, the relative differences between released fraction for the cylinder and sphere are considerably larger than in cases of greater total fraction of release. See text for further details.

specific values of Dt/a2 for 1%, 10%, 50%, 75% and 90% loss of diffusant from a sphere of radius a, derived from the expression presented by Crank (1975, Eq. 6.20) for the total amount of diffusing substance entering or leaving a sphere. The CYLMOD code is run using these times, and aspect ratios for the cylindrical grains range from 1:1 to 8:1. From the fraction of He released for each simulation, we determine a “release difference”, defined as (Hers − Herc)/Hers where Hers is the fraction of He released from the sphere (e.g., 0.1 for 10% loss of diffusant) heated for a specific time, and Herc is the fraction released for a cylinder of specific aspect ratio heated for the same time as the sphere. For aspect ratios near one, “the release difference” is small, as might be expected as the cylinder approaches the dimensions of a sphere. As aspect ratios increase, a greater fraction of He is retained in the cylindrical grains compared with the sphere because of the increasing diffusion distance in the axial direction. Another important observation is that although the total amounts are low for small fractions of He release, the relative differences between released fraction for the cylinder and sphere are considerably larger than for greater total fractions of release. For example, for 1% release the “release difference” for a grain of 8:1 aspect ratio is ~33%, while for 90% release the “release difference” for the same dimension is ~ 14%. While this represents a specific example, the results can be generalized to describe the effects of shape anisotropy at other temperatures (and diffusivities) and grain radii, provided the dimensionless parameters Dt/a2 have equal values. For the case of “shape anisotropy” with cylinder having isotropic diffusion, one can also determine an “effective radius” for a sphere which would have equivalent diffusional loss of He. It has been demonstrated analytically (e.g., Watson et al., submitted for publication) that equivalent losses will occur for spheres and cylinders of equal surface area/volume ratios. Using the standard expressions for surface areas and volumes of spheres and cylinders, and setting these ratios equal, we have the expression for the effective radius, Re, where h and a are the height and radius of the cylinder, respectively: Re =

3ah 2ðh + aÞ

ð5Þ

It may also be noted that this expression can be derived from Eq. (3) above by setting Ds = D11 = D33. From Eq. (5), it can be seen that the effective radius will be equal to a when h = 2a (equant

D.J. Cherniak et al. / Chemical Geology 268 (2009) 155–166

cylinder). For the case of a cylinder of length 2 times its diameter (i.e., h = 4a) the effective radius will be 1.2 times the cylinder radius; as the cylinder length approaches infinity the effective radius approaches 1.5a. Diffusional loss of He (determined from the expressions of Crank, 1975) in spheres of radii Re calculated using Eq. (5) for various cylinder aspect ratios can be compared with the values for He loss from cylinders calculated with the CYLMOD code as described earlier in this section. He losses from the cylinders (with aspect ratios ranging from 2 to 16) and spheres of effective radii Re agree within less than 1% for He losses up to 50%. Acknowledgements We thank Wayne Skala for assistance with ion implantation. We appreciate the constructive comments of Ken Farley and an anonymous reviewer, which helped in improving the manuscript. This work was supported by grant no. EAR-0440228 from the National Science Foundation (to E.B. Watson). References Black, L.P., Gulson, B.L., 1978. The age of the Mud Tank carbonite, Strangways Range, Northern Territory. BMR J. Aust. Geol. Geophys. 3, 227–232. Boyce, J.W., Hodges, K.V., Olszewski, W.J., Jercinovic, M.J., 2005. He diffusion in monazite: implications for (U–Th)/He thermochronometry. Geochem. Geophys. Geosyst. 6, Q12004. Chaumont, J., Soulet, S., Krupa, J.C., Carpena, J., 2002. Competition between disorder creation and annealing in fluorapatite nuclear waste forms. J. Nucl. Mater. 301, 122–128. Costantini, J.-M., Trocellier, P., Haussy, J., Grob, J.-J., 2002. Nuclear reaction analysis of helium diffusion in britholite. Nucl. Instrum. Methods B 195, 400–407. Costantini, J.-M., Grob, J.-J., Haussy, J., Trocellier, P., Trouslard, P., 2003. Nuclear reaction analysis of helium migration in zirconia. J. Nucl. Mater. 321, 281–287. Crank, J., 1975. The Mathematics of Diffusion (2nd ed.) Clarendon Press, Oxford, 414 pp. Dieumegard, D., Dubreuil, D., Amsel, G., 1979. Analysis and depth profiling of deuterium with the D(3He, p)4He reaction by detecting the protons at backward angles. Nucl. Instrum. Methods 166, 431–445. Dodson, M.H., 1973. Closure temperatures in cooling geological and petrological systems. Contrib. Mineral. Petrol. 40, 259–274. Dodson, M.H., 1986. Closure profiles in cooling systems. Mat. Sci. Forum 7, 145–154. Ewing, R.C., Wang, L.M., 2002. Phosphates as nuclear waste forms. In: Kohn, M.J., Ravokan, J., Hughes, J.M. (Eds.), Phosphates—Geochemical, Geobiological, and Materials Importance: Reviews in Mineralogy and Geochemistry, vol. 48, pp. 673–699. Ewing, R.C., Lutze, W., Weber, W.J., 1995. Zircon: a host-phase for the disposal of weapons plutonium. J. Mater. Res. 22, 243–246. Ewing, R.C., Meldrum, A., Wang, L.M., Weber, W.J., Corrales, L.R., 2003. Radiation effects in zircon. In: Hanchar, J.M., Hoskin, P.W.O. (Eds.), Zircon: Reviews in Mineralogy and Geochemistry, vol. 53, pp. 387–425. Farley, K.A., 2000. Helium diffusion from apatite: general behavior as illustrated by Durango fluorapatite. J. Geophys. Res. 105, 2903–2914. Farley, K.A., 2002. (U–Th)/He dating: techniques, calibrations and applications. In: Porcelli, D., Ballentine, C.J., Wieler, R. (Eds.), Noble Gases in Geochemistry and Cosmochemistry: Reviews in Mineralogy and Geochemistry, vol. 47, pp. 819–844. Farley, K.A., 2007. He diffusion systematics in minerals; evidence from synthetic monazite and zircon structure phosphates. Geochim. Cosmochim. Acta 71, 4015–4024. Farley, K.A., Wolf, R.A., Silver, L.T., 1996. The effects of long-stopping distances on (U–Th)/ He ages. Geochim. Cosmochim. Acta 60, 4223–4229. Flowers, R.M., Ketcham, R.A., Shuster, D.L., Farley, K.A., 2009. Apatite (U–Th)/He thermochronometry using a radiation damage accumulation and annealing model. Geochim. Cosmochim. Acta 73, 2347–2365. Gosset, D., Trocellier, P., 2005. Determination of the helium thermal diffusion coefficient in britholite using a NRA method: new results. J. Nucl. Mater. 336, 140–144. Gosset, D., Trocellier, P., Serruys, Y., 2002. Determination of the helium diffusion coefficient in nuclear waste storage ceramics by a nuclear reaction analysis method. J. Nucl. Mater. 303, 115–124. House, M.A., Wernicke, B.P., Farley, K.A., Dumitru, T.A., 1997. Cenozoic thermal evolution of the central Sierra Nevada, California, from (U–Th)/He thermochronometry. Earth Planet. Sci. Lett. 151, 167–179. House, M.A., Wernicke, B.P., Farley, K.A., 1998. Dating topography of the Sierra Nevada, California, using apatite (U–Th)/He ages. Nature 396, 66–69. Mayer, M., Alimov, V.K., Roth, J., 2005. Differential cross-section of the D(3He, p)4He nuclear reaction and depth profiling of deuterium up to large depths. Nucl. Instrum. Methods B234, 169–175. McDowell, F.W., McIntosh, W.C., Farley, K.A., 2005. A precise 40Ar–39Ar reference age for the Durango apatite (U–Th)/He and fission-track dating standard. Chem. Geol. 214, 249–263. Meldrum, A., Boatner, L.A., Ewing, R.C., 1997. Displacive radiation effects in the monaziteand zircon-structure orthophosphates. Phys. Rev. B56, 13805–13814. Meldrum, A., Zinkle, S.J., Boatner, L.A., Ewing, R.C., 1999. Heavy-ion irradiation effects in the ABO4 orthosilicates: decomposition, amorphization, and recrystallization. Phys. Rev. B 59, 3981–3992.

165

Miro, S., Studer, F., Costantini, J.-M., Haussy, J., Trouslard, P., Grob, J.-J., 2006. Effect of composition on helium diffusion in fluoroapatites investigated with nuclear reaction analysis. J. Nucl. Mater. 355, 1–9. Miro, S., Studer, F., Costantini, J.-M., Berger, P., Haussy, J., Trouslard, P., Grob, J.-J., 2007. Effect of gold ion irradiation on helium migration in fluoroapatites investigated with nuclear reaction analysis. J. Nucl. Mater. 362, 445–450. Moller, W., Besenbacher, F., 1980. A note on the 3He + D nuclear-reaction cross section. Nucl. Instrum. Methods 168, 111–114. Nasdala, L., Reiners, P.W., Garver, J.I., Kennedy, A.K., Stern, R.A., Balan, E., Wirth, R., 2004. Incomplete retention of radiation damage in zircon from Sri Lanka. Am. Mineral. 89, 219–231. Ouchani, S., Dran, J.-C., Chaumont, J., 1998. Exfoliation and diffusion following helium ion implantation in fluorapatite: implications for radiochronology and radioactive waste disposal. Appl. Geochem. 13, 707–714. Palenik, C.S., Ewing, R.C., 2002. Microanalysis of radiation damage across a zoned zircon crystal. Scientific Basis for Nuclear Waste Management XXV. Symposium, Materials Research Society Symposium Proceedings, vol. 713, pp. 521–526. Palenik, C.S., Nasdala, L., Ewing, R.C., 2003. Radiation damage in zircon. Am. Mineral. 88, 770–781. Paszti, F., 1992. Microanalysis of He using charged particle accelerators. Nucl. Instrum. Methods B66, 83–106. Payne, R.S., Clough, A.S., Murphy, P., Mills, P.J., 1989. Use of the D(3He, p)4He reaction to study polymer diffusion in polymer melts. Nucl. Instrum. Methods B42, 130–134. Pronko, P.P., Pronko, J.G., 1974. Depth profiling of 3He and 2H in solids using the 3He(d, p) 4He resonance. Phys. Rev. B 9, 2670–2678. Reich, M., Ewing, R.C., Ehlers, T.A., Becker, U., 2007. Low-temperature anisotropic diffusion of helium in zircon: implications for zircon (U–Th)/He thermochronometry. Geochim. Cosmochim. Acta 71, 3119–3130. Reiners, P.W., Farley, K.A., 1999. Helium diffusion and (U–Th)/He thermochronometry of titanite. Geochim. Cosmochim. Acta 63, 3845–3859. Reiners, P.W., Farley, K.A., 2001. Influence of crystal size on apatite (U–Th)/He thermochronology: an example from the Bighorn Mountains, Wyoming. Earth Planet. Sci. Lett. 188, 413–420. Reiners, P.W., Brady, R., Farley, K.A., Fryxell, J.E., Wernicke, B., Lux, D., 2000. Helium and argon thermochronometry of the Gold Butte Block, South Virgin Mountains, Nevada. Earth Planet. Sci. Lett. 178, 315–326. Reiners, P.W., Farley, K.A., Hickes, H.J., 2002. He diffusion and (U–Th)/He thermochronometry of zircon; initial results from Fish Canyon Tuff and Gold Butte. Tectonophysics 349, 297–308. Reiners, P.W., Spell, T.L., Nicolescu, S., Zanetti, K.A., 2004. Zircon (U–Th)/He thermochronometry; He diffusion and comparisons with 40Ar/ 39Ar dating. Geochim. Cosmochim. Acta 68, 1857–1887. Reiners, P.W., Ehlers, T.A., Zeitler, P.K., 2005a. Past, present, and future of thermochronology. In: Reiners, P.W., Ehlers, T.A. (Eds.), Low-temperature Thermochronology: Techniques, Interpretations, and Applications. Reviews in Mineralogy and Geochemistry, vol. 58. Mineralogical Society of America and Geochemical Society, Washington, DC, pp. 1–18. Reiners, P.W., Campbell, I.H., Nicolescu, S., Allen, C.M., Hourigan, J.K., Garver, J.I., Mattinson, J.M., Cowan, D.S., 2005b. (U–Th)/(He–Pb) double dating of detrital zircons. Am. J. Sci. 305, 259–311. Richter, F.M., Watson, E.B., Mendybaev, R.A., Teng, F.-Z., Janney, P.E., 2008. Magnesium isotope fractionation in silicate melts by chemical and thermal diffusion. Geochim. Cosmochim. Acta 72, 206–220. Richter, F.M., Watson, E.B., Mendybaev, R.A., Dauphas, N., Georg, B., Watkins, J., Valley, J., 2009. Isotopic fractionation of the major elements of molten basalt by chemical and thermal diffusion. Geochim. Cosmochim. Acta 73, 4250–4263. Roselieb, K., Rauch, F., Dersch, O., Buttner, H., 2006. Diffusivity and solubility of He in garnet: an exploratory study using nuclear reaction analysis. Nucl. Instrum. Methods B 244, 412–418. Ryssel, H., Ruge, I., 1986. Ion Implantation. J. Wiley, New York. Saadoune, I., Purton, J.A., de Leeuw, N.H., 2009. He incorporation and diffusion pathways in pure and defective zircon ZrSiO4: a density functional theory study. Chem. Geol. 258, 182–196. Shuster, D.L., Farley, K.A., 2004. 4He/3He thermochronometry. Earth Planet. Sci. Lett. 217, 1–17. Shuster, D.L., Farley, K.A., 2005a. 4He/3He thermochronometry; theory, practice, and potential complications. In: Reiners, P.W., Ehlers, T.A. (Eds.), Low-temperature Thermochronology: Techniques, Interpretations, and Applications. Reviews in Mineralogy and Geochemistry, vol. 58. Mineralogical Society of America and Geochemical Society, Washington, DC, pp. 181–203. Shuster, D.L., Farley, K.A., 2005b. Diffusion kinetics of proton-induced Ne-21, He-3, and He-4 in quartz. Geochim. Cosmochim. Acta 69, 2349–2359. Shuster, D.L., Farley, K.A., 2009. The influence of artificial radiation damage and thermal annealing on helium diffusion kinetics in apatite. Geochim. Cosmochim. Acta 73, 183–196. Shuster, D.L., Farley, K.A., Sisterson, J.M., Burnett, D.S., 2003. Quantifying the diffusion kinetics and spatial distribution of radiogenic 4He in minerals containing protoninduced 3He. Earth Planet. Sci. Lett. 217, 19–32. Shuster, D.L., Flowers, R.M., Farley, K.A., 2006. The influence of natural radiation damage on helium diffusion kinetics in apatite. Earth Planet. Sci. Lett. 249, 148–161. Tesmer, J.R., Nastasi, M., 1995. Handbook of Modern Ion Beam Materials Analysis. Materials Research Society, Pittsburgh, PA. Trocellier, P., Gosset, D., Simeone, D., Costantini, J.M., Deschanels, X., Roudil, D., Serruys, Y., Grynszpan, R., Saudé, S., Beauvy, M., 2003a. Application of nuclear reaction geometry for 3He depth profiling in nuclear ceramics. Nucl. Instrum. Methods B206, 1077–1082.

166

D.J. Cherniak et al. / Chemical Geology 268 (2009) 155–166

Trocellier, P., Gosset, D., Simeone, D., Costantini, J.-M., Deschanels, X., Roudil, D., Serruys, Y., Grynszpan, R., Saudé, S., Beauvy, M., 2003b. 3He thermal diffusion coefficient measurement in crystalline ceramics by NRA depth profiling. Nucl. Instrum. Methods B210, 507–512. Watson, E.B., Cherniak, D.J., 1997. Oxygen diffusion in zircon. Earth Planet. Sci. Lett. 148, 527–544. Watson, E.B., Wanser, K.H., Farley, K.A., submitted for publication. Anisotropic diffusion in a finite cylinder, with geochemical applications. Submitted to Geochimica et Cosmochimica Acta. Weber, W.J., Pederson, L.R., Gray, W.J., McVay, G.L., 1984. Radiation effects on nuclear waste storage materials. Nucl. Instrum. Methods B 229, 527–533. Weber, W.J., Ewing, R.C., Meldrum, A., 1997. Kinetics of alpha-decay-induced amorphization in zircon and apatite containing weapons-grade plutonium or other actinides. J. Nucl. Mater 250, 147–155.

Wolf, R.W., Farley, K.A., Silver, L.T., 1996a. Helium diffusion and low temperature thermochronometry of apatite. Geochim. Cosmochim. Acta 60, 4231–4240. Wolf, R.W., Farley, K.A., Silver, L.T., 1996b. Assessment of (U–Th)He thermochronometry: the low-temperature history of the San Jacinto Mountains, California. Geology 25, 65–68. Wolf, R.W., Farley, K.A., Kass, D.M., 1998. Modeling of the temperature sensitivity of the apatite (U–Th)/He thermochronometer. Chem. Geol. 148, 105–114. Young, E., Myers, A., Munson, E., Conklin, N., 1969. Mineralogy and geochemistry of fluorapatite from Cerro de Mercado, Durango, Mexico. U.S. Geol. Surv. Prof. Pap., vol. 650-D, pp. D84–D93. Zeitler, P.K., Herczig, A.L., McDougall, I., Honda, M., 1987. U–Th–He dating of apatite: a potential thermochronometer. Geochim. Cosmochim. Acta 51, 2865–2868. Ziegler, J.F., Biersack, J.P., 2006. The stopping and range of ions in matter. Computer code SRIM 2006, http://www.srim.org.