Earth and Planetary Science Letters The role of phosphates for the Lu–Hf chronology of meteorites

The Hf isotopic system is widely used for dating and tracing cosmochemical and geological processes, but still suffers from two uncertainties. First, Lu–Hf isochrons for some early Solar System materials have excess slope of unknown origin that should not be expected for meteorites with ages precisely determined with other isotopic chronometers. This observation translates to an apparent Lu decay constant higher than the one calculated by comparing ages obtained with various dating methods on terrestrial samples. Second, unlike the well constrained Sm/Nd value (to within 2%) for the chondritic uniform reservoir (CHUR), the Lu/Hf ratios in chondrites vary up to 18% when considering all chondrites, adding uncertainty to the Lu/Hf CHUR value. In order to better understand the Lu–Hf systematics of chondrites, we analyzed mineral fractions from the Richardton H5 chondrite to construct an internal Lu–Hf isochron, and set up a numerical model to investigate the effect of preferential diffusion of Lu compared to Hf from phosphate, the phase with the highest Lu–Hf ratio in chondrites, to other minerals. The isochron yields an age of 4647 ± 210 million years (Myr) using the accepted 176 Lu decay constant of 1 . 867 ± 0 . 008 × 10 − 11 yr − 1 . Combining this study with the phosphate fractions measured in a previous study yields a slope of 0 . 08855 ± 0 . 00072, translating to a 176 Lu decay constant of 1 . 862 ± 0 . 016 × 10 − 11 yr − 1 using the Pb–Pb age previously obtained, in agreement with the accepted value. The large variation of the Lu/Hf phosphates combined with observations in the present study identify phosphates as the key in perturbing Lu–Hf dating and generating the isochron slope discrepancy. This is critical as apatite has substantially higher diffusion rates of rare earth elements than most silicate minerals that comprise stony meteorites. Results of numerical modeling depending of temperature peak, size of the grains and duration of the metamorphic event, show that diffusion processes in phosphate can produce an apparently older Lu–Hf isochron, while this effect will remain negligible in perturbing the Sm–Nd chronology. Our results suggest that only type 3 chondrites with the lowest metamorphic grade and large minerals with minimal diffusive effects are suitable for determination of the Lu–Hf CHUR values and the Lu decay constant respectively. © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

The 176 Lu-176 Hf isotopic system is widely used for dating and tracing cosmochemical and geological processes, but still suffers from two uncertainties. First, Lu-Hf isochrons for some early Solar System materials have excess slope of unknown origin that should not be expected for meteorites with ages precisely determined with other isotopic chronometers. This observation translates to an apparent Lu decay constant higher than the one calculated by comparing ages obtained with various dating methods on terrestrial samples. Second, unlike the well constrained Sm/Nd value (to within 2%) for the chondritic uniform reservoir (CHUR), the Lu/Hf ratios in chondrites vary up to 18% when considering all chondrites, adding uncertainty to the Lu/Hf CHUR value. In order to better understand the Lu-Hf systematics of chondrites, we analyzed mineral fractions from the Richardton H5 chondrite to construct an internal Lu-Hf isochron, and set up a numerical model to investigate the effect of preferential diffusion of Lu compared to Hf from phosphate, the phase with the highest Lu-Hf ratio in chondrites, to other minerals.
The isochron yields an age of 4647 ± 210 million years (Myr) using the accepted 176 Lu decay constant of 1.867 ± 0.008 × 10 −11 yr −1 . Combining this study with the phosphate fractions measured in a previous study yields a slope of 0.08855 ± 0.00072, translating to a 176 Lu decay constant of 1.862 ± 0.016 × 10 −11 yr −1 using the Pb-Pb age previously obtained, in agreement with the accepted value. The large variation of the Lu/Hf phosphates combined with observations in the present study identify phosphates as the key in perturbing Lu-Hf dating and generating the isochron slope discrepancy. This is critical as apatite has substantially higher diffusion rates of rare earth elements than most silicate minerals that comprise stony meteorites. Results of numerical modeling depending of temperature peak, size of the grains and duration of the metamorphic event, show that diffusion processes in phosphate can produce

Introduction
The 176 Lu-176 Hf isotopic system is extensively used for dating cosmochemical and geological processes, and for studying planetary evolution. However, two uncertainties in the Lu-Hf systematic still need to be explained. First, Lu-Hf isochrons for many meteorites have steeper slopes than can be expected using the well-established decay rate of 176 Lu in terrestrial rocks. Isochron comparisons performed on terrestrial and extraterrestrial geological objects give similar values for λ 176 Lu ("terrestrial" average of Mesh indicates the size of the sieve in micron, and + or − indicate if the fraction was larger or smaller than the sieve; IF is the current (in ampere) used for magnetic separation. See analytical method for details. discrepancy, as even the most recent studies (Grinyer et al., 2003;Luo and Kong, 2006;Nir-El and Haquin, 2003;Nir-El and Lavi, 1998) show a considerable dispersion and give an average value of 1.848 ± 0.212 × 10 −11 yr −1 with a relatively large uncertainty that encompasses the two groups. Several processes have been proposed to explain this discrepancy between terrestrial and extraterrestrial objects, such as accelerated decay of 176 Hf that may have occurred at the beginning of the solar system under the influence of irradiation by γ -rays and cosmic rays (Albarede et al., 2006;Thrane et al., 2010) and branched decay of 176 Lu (Amelin and Davis, 2005;Norman et al., 2004). Accelerated decay of 176 Lu is ruled out by the absence of variations in Lu isotopic composition between the meteorites of different ages that formed before and after the proposed irradiation event (Wimpenny et al., 2015).
Branched decay is ruled out by the absence of Yb isotopic variations correlated with Lu/Yb ratios in ancient minerals (Amelin and Davis, 2005), and by γ -spectrometry measurements (Amelin and Davis, 2005;Norman et al., 2004). Hence there is a need to find another explanation for the observed excess slopes of meteorite Lu-Hf isochrons (e.g., Bizzarro et al., 2003;Thrane et al., 2010). It has recently been suggested that terrestrial contamination could cause excess radiogenic Hf and steeper Lu-Hf isochrons (Bast et al., 2017). However, it is not clear if this effect would be observable in all meteorites besides hot desert ones. Finally, because different minerals have different closure temperatures, they will provide different ages depending on the cooling rate (e.g., Göpel et al., 1994). In that case, the whole-rock isochron will be shifted towards lower slopes and thus younger ages (i.e. an average value between the oldest and the youngest minerals), in contrast to a steeper slope that translates to higher λ 176 Lu as investigated here. Second, unlike the Sm/Nd ratios that show little variation in type 1-6 bulk chondrite samples and allow the Sm/Nd ratio in the chondritic uniform reservoir (CHUR) to be constrained within ∼2% (Bouvier et al., 2008;Jacobsen and Wasserburg, 1984;Patchett et al., 2004), the Lu/Hf ratios in chondrites vary up to 28% when all petrologic types (1-6) are considered (Bizzarro et al., 2003;Bouvier et al., 2008;Patchett et al., 2004). A difference in the variation range of Lu-Hf was also observed between un-metamorphosed (type 1-3) ( Lu/Hf ∼4%) and metamorphosed chondrites ( Lu/Hf ∼28%) (Bouvier et al., 2008), as well as a subsequent shift in the slope of the isochron (Bouvier et al., 2015). The initial 176 Hf/ 177 Hf in the early solar system has been constrained by direct analysis of meteorite zircons (Iizuka et al., 2015), but the uncertainty of the CHUR value for the Lu/Hf ratio remains. The discrepancy is also observed between the average chondritic 143 Nd/ 144 Nd vs. 176 Hf/ 177 Hf ratio (see compilation in Martin et al., 2013). This problem has been addressed by Bouvier et al. (2008) who suggested that only type 1-3 chondrites with the lowest metamorphic grade should be used to determine the Lu decay constant and the CHUR values, because they provide the smallest Lu-Hf variation amongst chondrites of 4%. It should be noted that the case is not totally closed, as even some type 3 chondrites were not considered in proposing this well-defined average (Bouvier et al., 2008). Bouvier et al. (2008) also proposed a new Lu decay constant of 1.884 ± 0.060 × 10 −11 yr −1 , which is intermediate between the "terrestrial" and "meteoritic" values but has an uncertainty too large for usage of the 176 Lu-176 Hf system in geochemistry. Finally, Martin et al. (2013) observed a large redistribution of rare earth elements (REE) among different phases during metamorphism, while Hf is retained in its initial host phases. This raises the possibility that, under certain circumstances, preferential loss of Lu from phosphates, compared to Hf, may affect Lu/Hf isochron ages, even though this seems to be at odds with the study of Amelin (2005) who obtained a precise "terrestrial" Lu decay constant exclusively on phosphate grains from the Richardton H5 chondrite. We articulate in this paper that mineral grain size plays a significant role in determining whether or not a meaningful isochron can be obtained.
In order to better understand the Lu-Hf systematics of chondrites, various silicate fractions from the Richardton H5 chondrite have been analyzed to construct an internal Lu-Hf isochron. This meteorite has been well characterized, with an accurate U-Pb age on the pyroxene fraction of chondrules of 4.5627 ± 0.0017 billion years (Gyr), and one on phosphate of 4.5507 ± 0.0026 Gyr . Besides, a Lu-Hf isochron has also been obtained on pure phosphate fractions, giving, by comparison with the U-Pb age, a "terrestrial" like (e.g., Scherer et al., 2001) Lu decay constant of 1.8640 ± 0.0146 × 10 −11 yr −1 (Amelin, 2005). The data obtained here on the silicate fractions are first compared to these previous data. In a second step, diffusion modeling has been performed to evaluate the possibility that metamorphism can lead to an apparent variation in the Lu decay constant.

Analytical methods
The different fractions have been obtained by gently crushing the sample, followed by sieving and magnetic separation, as described in Table 1. Lutetium and Hf have been chemically purified at the Université Libre de Bruxelles (ULB) using the procedure described in Debaille et al. (2007), without any leaching. In brief, after dissolution, a 5% aliquot was separated to be spiked with a mixed 179 Hf-176 Lu spike. Hf and REE were separated on a cationic resin in HCl. Rare earth elements were subsequently purified on HDEHP resin in HCl in order to collect the Yb-Lu cut . Finally, the Hf cut was first purified from Fe on an anionic column in 6N HCl, and the remaining matrix (notably Ti) was eliminated on Ln.Spec resin with 6N HCl plus traces of H 2 O 2 . Hafnium was collected in 4N HF. Hafnium cuts and spiked Lu and Hf cuts have been measured at ULB in 0.05N HNO 3 (+0.05 N HF for Hf cuts) on Nu-Plasma MC-ICP-MS equipped with a DSN-100 desolvating nebulizer. The repeated measurements of JMC stan- in the numerical model. The acceptance or loss of REE and Hf by silicate minerals is driven by the difference between the initial concentration and the chemical equilibrium for each mineral for an element. While Lu (in grey) can be accepted in pyroxene and olivine, Hf (in black) can be accepted in plagioclase (see Table 3). Lu and Hf diffusing out of minerals can also be transferred to aqueous fluids during metamorphism, and will likely be incorporated in another newly-formed phosphate, but outside from the small subsample considered here. (b) Schematic evolution of the εHf of the simplified system with time (not at scale). T 0 is the beginning of the solar system; T 1 is the accretion of the chondrite parent body, and should represent the age of the chondrite if no metamorphic event occurred; T 2 is the closure time of the chronometer, i.e. the end of the metamorphic event. New concentrations and 176 Hf/ 177 Hf ratios are modeled as instantaneous diffusion at T 2 without considering radiogenic ingrowth during the metamorphic event. Note that only the εHf of plagioclase is modified after the metamorphic event because it is the only mineral accepting Hf (not seen on the cartoon). dard led to an average 176 Hf/ 177 Hf value of 0.282150 ± 10 (2σ , n = 15) and measured values were re-normalized to the JMC-475 accepted value of 0.282163 . Internal reproducibility is better than the total external reproducibility of 36 ppm. Blank values during the analytical campaigns were better than 20 pg for Hf and 17 pg for Lu. No correction was applied as those blank values are negligible. No neutron capture correction was applied on the data, as no resolvable anomalies in 178 Hf/ 177 Hf or 180 Hf/ 177 Hf from a terrestrial standard were found, in agreement with the study of Sprung et al. (2010). Spike stripping has been performed by iteration as in Debaille et al. (2007). Results are presented in Table 1.

Diffusion model set-up
The model set up here is based on the equations developed by Watson and Cherniak (2013). This model considers for a metamorphic event a parabolic curve for the temperature vs. time profile (Watson and Cherniak, 2013). This model allows us to consider an asymmetric temperature profile, as is likely the case for Richardton, with a relatively slow cooling rate of 26 • C/Myr from 4.5627 Gyr to 4.5507 Gyr . The final diffusive loss has been shown to be identical for symmetric and asymmetric parabolic curves, provided that the peak temperature and total heating + cooling timescales are the same in both cases (Watson and Cherniak, 2013). Richardton meteorite is an H5 type that has been well characterized, with a metamorphic peak of 1055 ± 106 K , in agreement with the temperature range of 1000-1140 K estimated for type 5 chondrites generally (Monnereau et al., 2013). The cooling rates and peak temperatures for Richardton are based on two U-Pb isochrons, one on chondrules and pyroxene fractions, and one on phosphate fractions .
A simple system is used, with only four minerals (phosphate, plagioclase, olivine and pyroxene) that lose or gain elements by diffusion in response to a metamorphic event (Fig. 1a). No distinction is made between the different phosphates (apatite, merrillite) or pyroxenes (low-Ca or high-Ca). The mineral abundances and initial concentrations used in the model are from Martin et al. (2013), based on the values obtained on the CK4 NWA 5798 only for Lu and Hf concentrations, and from different samples for Sm and Nd when data were missing (Table 2). Values from a lower metamorphic grade have not been considered, as plagioclase and phosphate are nonexistent to extremely small (<15 μm for phosphates), e.g., Rubin and Grossman (1985) in type 3 chondrites. The concentrations considered in the diffusion model are for illustrative purposes only, and will be different for different compositional groups and petrologic types of chondrites. However, comparing the different types of chondrites is beyond the scope of this paper. The amount of Lu and Hf leaving phosphate by diffusion is calculated following the equation from Watson and Cherniak (2013): where log ζ represents (indirectly) the fraction (F ) of element lost due to diffusion; D 0 is the pre-exponential factor as determined by diffusion experiments; E a , also determined by diffusion experiments, is the activation energy in the Arrhenius expression for diffusion of the element of interest in the mineral of interest; t is the duration of the metamorphic event; a is the radius of the mineral which is approximated as a sphere; T max is the maximum temperature reached in the metamorphic event; and R is the ideal gas constant. The translation from log ζ to F is not linear and has to be calculated over small intervals (Watson and Cherniak, 2013). It is important to clarify the meaning of F in the present context. Because Lu is compatible in phosphate, the factor F does not apply to the total amount of Lu that is held in this mineral, but instead to the fraction of Lu that can be lost by diffusion in order to reach chemical equilibrium between phosphate and surrounding phases. The equilibrium concentrations can be calculated by using the partition coefficients for each mineral in the system (Table 3). One significant uncertainty here is the role of metamorphic/hydrothermal fluids. While carbonaceous chondrites are known to contain a few percent of water, ordinary chondrites have been considered relatively dry (Hutchison, 2004). However, hydrous and/or halogen-rich metamorphic fluids may have played an important role before and during metamorphism, even in ordinary chondrites (Jones et al., 2014(Jones et al., , 2016. In order to calculate the equilibrium concentration, partitioning of Lu between phosphate and metamorphic fluids also has to be taken into account. Estimating the amount of metamorphic fluids that could have circulated through a chondrite is relatively complicated, so it is assumed here that this amount is equal to the amount of hydrated mineral (i.e. (chloro-)apatite) in OC, ∼0.1% (Jones et al., 2014). This value is a minimum, as the water content in type 3 ordinary chondrites could be as high as 1 to 2.5% (Rubin and Grossman, 1985), and up to 10% in carbonaceous chondrites (Hutchison, 2004). However, the exact value will have little influence on the model presented here. From the partition coefficient between minerals and hydrothermal fluids (Table 3), it is possible to calculate the ideal concentration between each mineral at equilibrium by equations (2) and (3): where C is the concentration, X is the proportion of the mineral (or fluid) in the system (see Table 2, recalculated at 100% with 0.1% of hydrothermal fluids), and K mineral/mineral D the partition coefficient between two minerals (Table 3). Table 4 shows the concentrations of Lu, Hf, Sm and Nd at equilibrium for the different minerals. The main observation between Table 2 (measured concentrations) and Table 4 (equilibrium con- The mineral starting composition are presented in Table 2, where the fluid starting composition is assumed to be equal to 0, and its proportion 0.1%. Concentration at chemical equilibrium of one element reflects the ideal composition of a mineral based solely on the partition coefficient of one element within this mineral. centration) is that in the presence of hydrothermal fluids, phosphate can lose up to 86% of Lu, while plagioclase is susceptible to accept Hf, and olivine and pyroxene to accept very small amounts of REE, including Lu. The factor F representing the amount of element lost by diffusion as in eq. (1) concerns only the difference between the measured and the equilibrium concentrations. A factor F of 1 represents the most severe loss of Lu from the phosphate, i.e. 86% from the starting composition in the present case. In the simplified system considered here, phosphate is the mineral with the highest diffusion rates, with log D REE of ∼ −23 at 750 • C (Cherniak, 2000), compared to olivine (∼ −24) (Chakraborty, 2010;Cherniak and Dimanov, 2010), plagioclase (∼ −25) (Cherniak, 2010a(Cherniak, , 2003, and pyroxene (∼ −24/−26) (Bloch et al., 2017;Van Orman et al., 2001) at the same temperature. Part of the Lu, Hf, Sm and Nd fluxes going out from minerals will diffuse within other minerals able to accept them according to their own diffusion rate and chemical equilibrium, while the rest is considered to be taken away by metamorphic fluids, and will likely be incorporated in newly-formed phosphates (Jones et al., 2014) -but outside the small system considered here (Fig. 1a). As stated earlier, the role of metamorphic fluids can be very important. Classically, a distinction is made between open vs. closed systems. In the first case, metamorphic fluids leave the studied subsample, while in a closed system elements are transferred among phases but remain within the considered subsample. In the second case, the bulk composition of the subsample does not change. Of course, whether open or closed system behavior occurred is mainly dependent on the size of the subsample considered. In the case of the simplified system built here any new phosphate will grow outside of the considered subsample, i.e. in an open system, because of the small size of the subsample. While growth of phosphates over increasing metamorphic conditions is well documented (e.g., Jones et al., 2014;Van Schmus and Ribbe, 1969), there is a complex relationship between loss by diffusion and incorporation of Lu brought by aqueous fluids during the simultaneous growth of phosphates. This complexity is not considered here, as it is complicated to assess. In addition, no other phase has been considered, such as the matrix or sulfide where REE content can be high (e.g., Martin et al., 2013).
Concerning the diffusion of Hf, no data exist for the three minerals considered here (plagioclase, olivine and phosphate). Diffusion of tetravalent cations in most minerals where experimental diffusion data exist is known to be slower than diffusion of trivalent cations such as REE. This is illustrated by the diffusion of Hf in diopside (Bloch and Ganguly, 2014) and Ti in olivine and diopside Liang, 2012, 2014). Concerning phosphate, the most important mineral considered here, the diffusion rate of U can be taken as a proxy, as both U and Hf are tetravalent cations. Among tetravalent cations, smaller cations such as Ti will diffuse faster than larger ones such as Hf or U (Cherniak et al., 1997;Cherniak and Liang, 2012). This is observed in diopside where at 750 • C, log D U is ∼ −27 (Van Orman et al., 1998) and log D Hf is This run has been performed for a peak temperature of 750 • C, for a metamorphic event of 10 Myr, and phosphate size of 50 μm. For the purpose of modeling the flux of Lu and Hf from other minerals, it has been considered that olivine size = 200 μm; pyroxene and plagioclase size = 100 μm, and the diffusion data are, for plagioclase: Lu (Yb, Cherniak, 2003), Hf (from pyroxene (Bloch and Ganguly, 2014)); pyroxene: Lu (Bloch et al., 2017), Hf (U, (Van Orman et al., 1998)); phosphate: Lu (Yb, Cherniak, 2000), Hf (U, Cherniak, 2005); olivine: Lu (Yb, Cherniak, 2010b), Hf (Ti, Cherniak and Liang, 2014). SSF: solar system formation, starting values are the CHUR parameters (Bouvier et al., 2008). A distinction is made at T 2 before and after the metamorphic event, as it is modeled that chemical re-equilibration is instantaneous at T 2 , for the sake of simplicity. a The bulk amount of Lu and Hf is not the same before and after metamorphism in the small subsample considered here because metamorphic fluids took away the elements. However, at larger scale, the sample will likely preserve its bulk chondritic composition. The diffusion of Lu and Hf from and to olivine, pyroxene and plagioclase has been taken into account but is barely seen for Lu. Acceptance of Lu and Hf is depending of the diffusion coefficient until the chemical equilibrium in the considered mineral would be reached (this is never the case in the present model). As plagioclase is the only mineral able to accept Hf, the resulting Hf isotopic composition is mass-balanced between the 176 Hf/ 177 Hf in plagioclase before the metamorphic event, and the average isotopic composition of Hf from olivine, phosphate and pyroxene. The resulting λ 176 Lu from this simulation is 1.8907 × 10 −11 yr −1 .
∼ −24 (Bloch and Ganguly, 2014). In the absence of available data for Hf diffusion in phosphate, the diffusion of U has been considered with log D U of ∼ −26 at 750 • C (Cherniak, 2005), i.e. 3 orders of magnitude slower than Lu. The diffusion of Ti in olivine has been taken as a proxy for Hf (Cherniak and Liang, 2014) and the diffusion of Hf in plagioclase has been taken as in diopside (Bloch and Ganguly, 2014). The diffusion history of this system is divided in the following steps (Fig. 1b): T 0 = 4.5685 Gyr (Bouvier et al., 2007) T 1 = 4.5664 Gyr; accretion of the chondrite parent body. This time has to be delayed relative to T 0 because of the presence of 26 Al that could melt the asteroid parent body if accreted earlier. It is arbitrarily considered here that the accretion time corresponds to 3 half-lives of 26 Al after T 0 as this value has no real importance on the model. Each mineral then evolves in 176 Hf/ 177 Hf at its own rate, depending on its Lu/Hf ratio.
T 2 = 4.5564 Gyr (Amelin, 2005); closure time of the system for diffusion, at the end of the metamorphic event. It is considered here that the closure time is the same for all minerals. It is difficult to evaluate simultaneously the Hf radiogenic ingrowth occurring during Lu diffusion. For simplicity, the new Lu/Hf and 176 Hf/ 177 Hf ratio are thus calculated at T 2 , when the temperature is too low to allow diffusion anymore. This simplification is reasonable as the duration of the metamorphic event (0-10 Myr) is short compared to the half-life of 176 Lu (35.6 Gyr). The metamorphic event starts sometime between T 1 and T 2 , and can last maximum 10 Myr, i.e. the difference between T 1 and T 2 . Finally, after T 2 , the isotopic composition of each mineral evolves until the present time (Fig. 1b), when an isochron can be calculated. In this simplified system, as explained previously, mainly phosphate will be impacted by the metamorphic event.
Knowing T 1 , the actual age of the chondrite without any metamorphic event, an apparent Lu decay constant can be calculated by comparing the slope of the isochron after the metamorphic event.
The model is influenced by the three main parameters of eq. (1): maximum temperature reached, duration of the metamorphic event and grain size of the phosphate. Each simulation has been run with one parameter fixed and two parameters varying. Temperature is varying from 400 to 950 • C, duration from 0 to 10 Myr and size of phosphate from 10 to 100 μm. When those parameters are fixed in each simulation, these are 700 • C, 5 Myr and 50 μm respectively. As an illustration, a single run has been performed at a peak temperature of 750 • C, for a duration of 10 Myr and a phosphate size of 50 μm, and is presented in Table 5.

Isotope analyses
The Lu-Hf regression obtained on silicate fractions of Richardton (Table 1) yields an age of 4647 ± 210 Myr (MSWD = 9.4) using the λ 176 Lu of 1.867 × 10 −11 yr −1 (Fig. 2a). Low precision is caused by a relatively small spread of the 176 Lu/ 177 Hf ratios from 0.02 to 0.05, and by scatter of the data points exceeding analytical uncertainty.
The old age of 4647 ± 210 Myr agrees with the currently accepted age of the solar system only within its large uncertainty limits. Interestingly, the slope of our regression is also consistent with the slope of the Lu-Hf isochron defined by multiple fractions of phosphate minerals (apatite and merrillite) from the Richardton meteorite (Amelin, 2005). Combining all the fractions from the two studies yields a slope of 0.08855 ± 0.00072, similar to but more precise than the values obtained with each data set separately (Fig. 2b). Richardton phosphates show an extreme variation of the 176 Lu/ 177 Hf ratios from 0.8 to 143 (Amelin, 2005). This slope is consistent with the U-Pb ages obtained previously  and using the terrestrial Lu decay constant of 1.867 ± 0.008 × 10 −11 yr −1 (Scherer et al., 2001). Fig. 3 shows the results for different simulations obtained from the numerical model detailed above and showing the effect of temperature, phosphate grain size and duration of the metamorphic event on the apparent Lu decay constant. The idea is to observe the change of the isochron slope after the metamorphic event due to the preferential diffusion of Lu over Hf from phosphate. Then, this slope change is translated into an apparent Lu decay constant change that would be required to obtain the correct age of  (Ludwig, 2003), and (b) by combining the silicate fractions of this study with the phosphate fractions (circle) from Amelin (2005). Error bars are 2σ . the sample. The apparent Lu decay constant can vary over a large range, systematically higher than the expected 176 Lu decay constant 1.867 ± 0.008 × 10 −11 yr −1 (Scherer et al., 2001), depending on the metamorphic conditions and size of phosphates. All the curves reach a plateau where no more variation of 176 λ is observed. This is related to the fact that the fraction F of element that can be lost by diffusion from the phosphate is finite, and the plateau indicates that F has reached 1. In that case, the remaining amount of that element is the amount predicted by chemical equilibrium, and no more variation can be recorded. Considering the simulation example shown in Table 5 and Fig. 4, it should be noted that for the metamorphic event considered here, F = 0.98 (out) and 0.05 (out) in phosphate for Lu and Hf respectively. On the other hand, for plagioclase, F = 0.02 (out) and 0.19 (in; this last value is poorly constrained as taken from diffusion in diopside) for Lu and Hf respectively; for olivine, F = 0.14 (in) and . After the metamorphic episode, phosphate will be perturbed by the preferential loss of Lu over Hf. The resulting slope of the new isochron is expressed in term of apparent λ 176 Lu using the theoretical age that should be obtained if no metamorphic event occurred (T 1 ). Note that at higher temperature, when Hf is finally diffusing out from phosphate, the λ 176 Lu will slightly drop down, but still remains higher than 1.867 ± 0.008 × 10 −11 yr −1 (Scherer et al., 2001).  Table 5. In order to process the artificial isochrons in Isoplot (Ludwig, 2003), the associated error has been set up to 0.5% (2RSD) on the 176 Lu/ 177 Hf ratio, and 35 ppm (2RSD) on the 176 Hf/ 177 Hf. Data points for the artificial "normal" isochron are not shown, and white circles corresponds to the isochron after metamorphism. The "normal" isochron is characterized by a MSWD of 0 as this isochron has no scatter. On the other hand, the metamorphic event resulted in both an apparent older age, an initial 176 Hf/ 177 Hf associated with a larger uncertainty and a larger MSWD. From the new slope of the isochron after metamorphism, the apparent λ 176 Lu is 1.8907 × 10 −11 yr −1 . 0.19 (out); and for pyroxene, F = 0.04 (in) and 0.19 (out). The direction of the flux (in or out) is depending to the difference between the initial and the chemical equilibrium concentrations. Since diffusion of Lu into olivine and pyroxene and out of plagioclase is so slow, no variation is recorded. The resulting apparent 176 Lu decay for the parameters entered in the numerical model is 1.891 ± 0.002 × 10 −11 yr −1 .

The role of phosphates in disturbing the Lu-Hf chronology
The large variation in Lu/Hf ratio compared to the range obtained in the present study focusing on silicate fractions identifies the key role of phosphate in the Lu-Hf isotope systematics of chondrites and other meteorites, as the major host of Lu. Calcium phosphate minerals are common accessory minerals in both chondrites (e.g., Van Schmus and Ribbe, 1969) and achondrites (e.g., Hsu and Crozaz, 1996) that concentrate parent elements of several isotopic chronometers: Lu, Sm, U and Th. An important difference between the role of phosphates in the Sm-Nd systematics, on one hand, and Lu-Hf and U-Th-Pb, on the other, is that the phosphates effectively exclude Hf and Pb during crystal growth, and are similar in these sense to the "classical" geochronometer mineral systems such as U-Pb in zircon, and Rb-Sr and K-Ar in mica. Sm and Nd in phosphates are, however, only moderately fractionated from the silicate minerals (Amelin and Rotenberg, 2004;Brannon et al., 1988;Jacobsen and Wasserburg, 1984;Martin et al., 2013).
Extreme fractionation between parent and daughter elements by phosphates in the Lu-Hf and U-Pb systems can be either bene-ficial or detrimental for isotopic dating, depending on whether the phosphates remained closed to diffusion. By increasing the spread of the Lu/Hf ratios, this fractionation helps to obtain precise ages using the isochrons that contain phosphate minerals with very high Lu/Hf ratios, silicate minerals with moderate to low Lu/Hf ratios, and possibly oxide minerals (ilmenite and spinels) with very low Lu/Hf ratios. But it has also been observed that the presence of phosphates can induce poor reproducibility between duplicates in Lu-Hf measurements (Bizzarro et al., 2003;Bouvier et al., 2008;Patchett et al., 2004), which critically depends on bulk rock sampling size.
Apatite has substantially higher diffusion rates of elements such as REE (Cherniak, 2000) and Pb (Cherniak et al., 1991) than most silicate minerals that comprise meteorites. It can, therefore, become open to diffusion under metamorphic conditions and lose both parent and daughter elements if conditions are met. As a result, the phosphate Lu-Hf age would be perturbed, depending on the respective loss of Lu and Hf. At first thought, isotopic systems of the minerals with low Lu and Hf concentrations, such as plagioclase and olivine, might be anticipated to be more susceptible to the gain of radiogenic Hf than relatively Hf-rich pyroxene and oxide minerals. This would be consistent with the interpretation of Sm-Nd isotopic exchange between phosphates and silicate minerals in achondrites (Prinzhofer et al., 1992) that indicated that the 143 Nd/ 144 Nd ratios in plagioclase can be elevated due to such exchange, while the mafic minerals are unaffected because they are already rich in Nd. However, the final flux for one element will depend on the mineral assemblage considered in the system, both for the diffusion rate of this element in each mineral, and for its concentration in each mineral at chemical equilibrium, determined by the partition coefficient. The conditions of metamorphism (such as dry vs. wet, and temperature) will also influence the outcome of the perturbation. Martin et al. (2013) observed that while REE are borne first by phosphate and then silicate in type 3 CK chondrites, the proportion of REE borne by phosphates decreases with metamorphism, despite the compatibility of these elements in phosphate. On the other hand, the Hf content remains uniform over the whole metamorphic range. These authors also observed that Sm and Nd behave coherently. The behavior of an element under metamorphic conditions is an interplay between diffusion rate and partitioning between the minerals, the first being independent of the mineralogical assemblage, while the second strongly depends of the different possibilities for an element to be incorporated in mineral phases at equilibrium concentration.
The magnitude of the Hf isotopic disturbance in the phosphates is thus directly related to the loss of slowly diffusing but incompatible radiogenic Hf compared to the loss of fast diffusing but compatible Lu. As stated previously, the main parameters will be the duration of the metamorphic event, the temperature reached and the size of the phosphate crystals. Analysis of coarse phosphate crystals with diameter >50-80 μm that are typically handpicked from meteorites for geochronological studies (e.g.,  are relatively little affected by diffusion-related loss of Lu, and yield more reliable ages than would be obtained by extraction of tiny phosphate grains closely intergrown with silicate minerals by selective acid dissolution. On the other hand, small grain size of phosphates in chondrites of low metamorphic grade is beneficial for determinations of the composition of CHUR, because it reduces the spread of Lu/Hf between the specimens caused by uneven occurrence of phosphate grains (Amelin and Davis, 2005).

The role of diffusion in disturbing the Lu-Hf chronology
The shift towards larger apparent 176 λ values after the metamorphic event as seen in Fig. 3 corresponds to what is observed in several types of meteorites in previous studies (Bizzarro et al., 2003(Bizzarro et al., , 2012Blichert-Toft et al., 2002;Patchett and Tatsumoto, 1980;Thrane et al., 2010). The system presented here is very simple compared to the reality, as (i) addition of Lu to phosphate by aqueous fluids is not modeled, and (ii) no other component has been considered, although a component such as the matrix could be important. The model also suffers from missing information about element partitioning between minerals and metamorphic fluids. However, it can address the essence of diffusion-related Lu-Hf isochron modification. The main observation is that the preferential diffusion of Lu over Hf from phosphate, the mineral in which Lu is most concentrated in chondrites, can increase the isochron slope and thereby generate a spuriously large apparent Lu decay constant. Such a conclusion has also been reached by Bloch and Ganguly (2015) and Bloch et al. (2017) with a numerical model using garnet and pyroxene respectively. Even though garnet is not a relevant phase in chondrites, it behaves similarly to phosphate regarding the preferential loss of Lu over Hf by diffusion and can lead to erroneous older ages (Bloch and Ganguly, 2015). Concerning pyroxene, the same process certainly occurred in chondrites and achondrites, but we suggest that because (i) pyroxenes are generally much larger compared to phosphates, hence being less affected by diffusion; (ii) Lu diffuses faster of one order of magnitude in phosphate compared to pyroxene (Bloch et al., 2017;Cherniak, 2000); and (iii) pyroxene tends to accept Lu instead of losing it according to chemical equilibrium concentration (Table 4), phosphates are the main mineral causing Lu-Hf isotopic disturbance in chondrites and also possibly in achondrites. With the model set up here, the apparent λ 176 Lu can go as high as 1.895 × 10 −11 yr −1 , resolvable above the accepted value of 1.867 ± 0.008 × 10 −11 yr −1 (Scherer et al., 2001). On Fig. 3a, it can be clearly seen that beyond a certain size of phosphate (∼60 μm), the induced variation on the apparent Lu decay constant is small. This corroborates the robustness of a previous study that focused only on large phosphate grains, where no shift was observed for λ 176 Lu (Amelin, 2005). Fig. 3b shows that below 600 • C, for any grain size of phosphate, no λ 176 Lu shift is observed, which reinforces the conclusion from Bouvier et al. (2008) that only type 1-3 chondrites should be used to evaluate Lu-Hf systematics. Finally, it appears that the maximum temperature reached is a very important factor in disturbing the Lu-Hf systematics (Fig. 3c), even though this effect could be counterbalanced by the fact that phosphate crystals may grow as the metamorphic temperature increases, and can incorporate Lu transported by metamorphic fluids (Jones et al., 2014). At higher temperatures (>750 • C), when Hf is finally starting to diffuse out from phosphate, the shift in the apparent λ 176 Lu is to slightly lower values (Fig. 3b). Beside the effect on the apparent Lu decay constant, it should be noted that the scatter around the isochron after the metamorphic event also increases, due to the decrease in the Lu/Hf spread between mineral separates, and uneven loss of Lu and Hf from phosphate (Table 5). This is shown in Fig. 4 where two isochrons, before and after metamorphism (Table 5), are plotted, with an increasing scatter after the metamorphic event. This is also observed in the present study with the Lu-Hf obtained on the mineral separates of Richardton, showing a shifted isochron (i.e. larger λ 176 Lu) associated with a larger error. Even if the present study disagrees on the origin of the disturbance, it agrees with the conclusions of Bast et al. (2017) that using perturbed fractions can produce erroneous Lu-Hf isochrons.

The case of Sm-Nd and U-Pb systematics
When the numerical modeling is applied to the Sm-Nd systematics (Fig. 5), a shift between an apparent λ 147 Sm and the actual λ 147 Sm can be produced, but is very small, less than 0.02% in relative variation. This is much lower than the error on λ 147 Sm of 0.7% (Lugmair and Marti, 1978). In comparison, relative variation observed here for λ 176 Lu can be up to 1.5%, larger than the error in λ 176 Lu of 0.4% (Amelin, 2005) obtained from age comparison studies. The goal here is not to propose a best value for λ 176 Lu, but to demonstrate that the excess slope of meteorite Lu-Hf isochrons and hence higher apparent λ 176 Lu could possibly be explained by the conjunction of thermal metamorphism and small-sized subsamples.
Despite its simplicity, this model clearly demonstrates that the preferential diffusion of Lu over Hf from phosphate can affect the Lu-Hf chronology, while the effect will remain minor for the Sm-Nd chronology, because (1) Sm and Nd have similar diffusion coefficients, and (2) the difference of Sm/Nd ratios between phosphate and silicate minerals is much smaller than the difference of Lu/Hf ratios between these minerals.
Considering the U-Pb chronometer, Lu-Hf and U-Pb radiogenic couples should behave the same during diffusion within phosphate because in both cases, the parent element is compatible in phosphate while the daughter element is not. However, Pb is divalent and will diffuse faster than U, compared to Hf which is tetravalent and diffuses slower than Lu. As such, the U-Pb dating will be easier to reset during metamorphic events towards younger ages, in contrast to the Lu-Hf dating, which could give older ages as explained in the present study. This has been observed by Barfod et al. (2003Barfod et al. ( , 2005, who measured Lu-Hf ages systematically older than U-Pb ages in metamorphic terrestrial apatites. They concluded that the closure temperature of Lu-Hf is higher (675 ± 25 • C) than that for U-Pb (480-525 • C) (Cherniak et al., 1991). Note that these closure temperatures cannot be extrapolated to the present study, as the terrestrial apatites are larger than 0.5 cm, much larger than the phosphates existing in meteorites and modeled in the present study (∼50 μm).
Diffusion is a process that may occur on a cm-scale over the duration of a metamorphic event. Aqueous fluids will entrain REE that diffuses out of phosphate, and these REE are likely to be incorporated into newly-formed phosphates located further away. There is some debate over the open-vs. closed-system processes during metamorphism in chondrites (Jones et al., 2014), but over a cm-scale, open-system behavior can be expected, and at a smaller scale it is almost certain. The size of chips and samples used for meteorite chronology is often smaller than one centimeter. This model implies that these chips are thus not equilibrated with regard to Lu-Hf chronology. The situation is different for type 1-3 chondrites that show low metamorphic disturbance. It should be noted that type 1-2 chondrites are also metamorphosed, but low temperature aqueous alteration seems to create less disturbance in the Lu-Hf system and consequently show less disturbance in their Lu-Hf systematics (Bouvier et al., 2008). This is illustrated by the conclusion of Bouvier et al. (2008), where a more robust average is obtained for the Lu-Hf system when using only type 1-3 chondrites although even some type 3 chondrites show some perturbation. This conclusion could also affect terrestrial chronology for metamorphic rocks. Garnet shows a similar behavior to phosphate in the sense that Lu is strongly compatible in garnet while Hf is not. Smit et al. (2013) observed a systematic difference between Sm-Nd and Lu-Hf ages in garnet, and interpreted it as reflecting the preferential loss of radiogenic Hf over Lu from the garnet (while Sm and Nd would be lost relatively equally). They proposed the net loss of Hf instead of Lu from the garnet because Lu is compatible in garnet. However, Baxter and Scherer (2013) indicated that the preferential loss of Lu (fast diffusion but compatible) over Hf (slow diffusion but incompatible) is possible depending of the equilibrium partitioning between the mineral and its surroundings. Finally, Bloch and Ganguly (2015) confirmed these conclusions by modeling older ages due to the preferential loss of Lu over Hf in garnet. We thus reach a similar conclusion that the diffusion rate of an element leaving a mineral can be a parameter more important than the compatibility of this element in this mineral for disturbing the Lu-Hf chronology of small-scale subsample. We also show that the size of the subsample might be critical for considering open vs. closed-system with regard to metamorphic fluids.

Conclusions
The isochron obtained for the silicate fractions of Richardton leads to an age consistent with the age of the solar system only within its large uncertainty. However, when those fractions are combined with the phosphate fractions previously obtained on Richardton, the slope is consistent with a "terrestrial" Lu decay constant for Richardton. We identify phosphate as a key to the perturbation of the Lu-Hf system. A diffusion model was set up to investigate this hypothesis. Despite its simplicity, the diffusion model developed here provides a robust explanation for an apparent larger meteoritic Lu decay constant observed in some meteorites, because of the preferential diffusion of Lu over Hf from phosphates grains. In order to obtain precise and robust isochrons, mineralogically pure handpicked phosphates should be used as they remain preserved from diffusive loss because of their large size (over ∼80 μm). On the other hand, selective acid leaching of meteorites containing dispersed micron sized and smaller phosphates could lead to preferential dissolution of disturbed phosphate grains. Low grade metamorphic events, below 600-650 • C, should prevent diffusion and preserve the Lu-Hf systematics. This conclusion reinforces the conclusion from Bouvier et al. (2008) that only type 3 meteorites should be used for constraining the average Lu-Hf CHUR value. On the other hand, the Sm-Nd systematics will remain robust during a metamorphic event, because diffusive exchange of Sm and Nd will be similar.
Two additional pieces of information are required for more accurate constraints on the role of phosphates in Lu-Hf systematics of meteorites: 1) Experimental data for diffusion of Hf in phosphate. There are also no experimental diffusion data for merrillite or chloroapatite. It is commonly assumed that the diffusion rates of various elements in these minerals are similar to those in apatite, but this may not be the case when considering compositional and structural differences. 2) Detailed mineral inventory of Lu and Hf, including fine-grained interstitial minerals if present, has to be determined in every chronological study. Despite the wide use of Lu-Hf chronology, caution should be taken, and a double-check using a second dating system should always be performed.