Numerical Modelling of Radiogenic Ingrowth and Diffusion of Pb in Apatite Inclusions with Variable Shape and U-Th Zonation

: The fundamental premise of apatite U-Th-Pb thermochronology is that radiogenic Pb is redistributed by volume diffusion. In practice, it is often additionally assumed that crystals (1) lose radiogenic Pb to an inﬁnite reservoir, (2) have a simple geometry and (3) are chemically homogeneous. Here we explore the signiﬁcance of the latter three assumptions by numerical modelling of Pb radiogenic ingrowth and diffusion in apatite inclusions within other minerals. Our results indicate that the host minerals are likely to hamper diffusive Pb loss from the apatite inclusions by limiting the Pb ﬂux across their boundaries, and thus the thermal histories that are reconstructed assuming a fully open boundary may be signiﬁcantly inaccurate, precluding a meaningful interpretation. We also ﬁnd that when apatite boundaries are ﬂux-limited, heterogeneities in U and Th concertation within apatite have subordinate effect on bulk-grain U-Th-Pb dates and can cause intra-grain U-Th-Pb dates to increase towards the boundaries. Finally, we show that it is important to correctly account for crystal geometry when modelling intra-grain U-Th-Pb dates. We suggest that the effect of surrounding minerals on diffusive Pb loss from apatite (and loss of other radiogenic isotopes from other minerals) should be examined more closely in future research. A heterogeneous distribution of U and Th in such apatite inclusions has a relatively minor to insigniﬁcant effect on their bulk-grain U-Th-Pb dates and can cause an increase in intra-grain U-Th-Pb dates towards the rim. We show that thermochronologists should make appropriate assumptions about the openness of apatite boundaries with respect to diffusive Pb loss when reconstructing t-T paths of rocks by inversion modelling of apatite U-Th-Pb data. We further show that it is important to account for apatite crystal geometry and U and Th zonation. The latter is challenging when dealing with intra-grain U-Th-Pb date variations, because these cannot be accurately predicted by modelling radiogenic ingrowth and diffusion of Pb in equivalent spherical crystals. Finally, apatite crystals that are not completely encapsulated in crystals of other minerals may also have ﬂux-limited boundaries, although there is insufﬁcient experimental data to make a fully constrained assessment. We conclude that the likelihood of ﬂux-limited boundary behaviour should be carefully considered in future research on apatite U-Th-Pb thermochronology, as well as thermochronology using other methods of isotopic dating.


Introduction
Many thermochronological methods use isotopic dates acquired by mass spectrometry, such as U-Th-Pb dates of apatite. These methods rely on a number of assumptions, violation of which may undermine the accuracy of the inferred time-temperature (t-T) paths. The only fundamental assumption is that volume diffusion dominates the isotopic transport in the analysed material, such that apparent resetting of the obtained isotopic dates can be quantitatively interpreted using some solution to the production-diffusion equation. In practice, however, additional assumptions are needed to solve the production-diffusion equation, and it is common to assume that [1][2][3][4][5][6][7][8][9][10][11][12][13][14]: (1) radiogenic isotopes are completely lost form the system once they reach crystal boundaries, (2) parent isotopes are homogeneously distributed within a crystal, (3) the crystal geometry can be approximated by a simple shape, such as a sphere.
Here we explore the validity of these additional assumptions using the example of apatite U-Th-Pb thermochronometer. We approach this problem by numerical modelling of radiogenic ingrowth and diffusion of Pb in apatite crystals that are completely encapsulated within crystals of other minerals, i.e., in inclusions of apatite within host minerals. First, we explore whether diffusive Pb loss from apatite inclusions can be sufficiently inhibited by their host minerals to affect the accuracy of t-T paths inferred by inversion modelling. Second, we assess the influence of a heterogeneous intra-grain distribution of U and Th on bulk-grain U-Th-Pb dates of apatite inclusions. Third, we examine the effect of crystal document disparate concentrations of trace elements at the rims of different rutile crystals within the same rock and argue that the disparity was caused by neighbouring minerals that limited the fluxes of trace elements proximal to the rutile crystals [42].
The cited examples highlight that there is yet much to be understood about the behaviour of radiogenic isotopes at crystal boundaries, and that assuming that they are completely lost to an infinite reservoir immediately after reaching crystal boundaries is not always valid. This assumption seems to be particularly questionable when the dated crystal is completely encapsulated within another crystal, such that radiogenic isotopes cannot be removed to an infinite reservoir through a network of grain boundaries. Here we provide MATLAB scripts to simulate radiogenic ingrowth and diffusion of Pb in such a system, albeit with certain limitations, and use these scripts to explore the effect of surrounding minerals on diffusive Pb loss from apatite inclusions.

Complex Geometry and Parent Isotope Zonation
The importance of accounting for crystal geometry and parent isotope zonation when modelling diffusive loss of radiogenic isotopes has been discussed in many studies [43][44][45][46][47][48][49], including some that used apatite U-Pb dating [8,11]. Following these studies, a chemically homogeneous spherical crystal of apatite is expected to yield a younger bulk-grain U-Th-Pb date than a chemically homogeneous spherocylindrical crystal of apatite with the same radius and thermal history. Its bulk-grain U-Th-Pb date is also expected to be younger than that of a spherical apatite crystal with identical dimensions and thermal history but increased U and Th concentrations in the core relative to the rim. The most accurate way to account for the observed crystal geometry and parent isotope zonation during thermochronological modelling is to perform three-dimensional simulations of production and diffusion of a radiogenic isotope. However, this is time-consuming and has not been widely used in thermochronological inversion modelling. Alternatively, a chemically zoned crystal with complex geometry can be represented using a 1-dimensional equivalent spherical model, in which case bulk-grain isotopic dates can be calculated with reasonable accuracy [43][44][45][46][47][48]. However, it remains underexplored whether analogous simplifications can be made to model intra-grain variations of isotopic dates. It is also unclear how crystal geometry and parent isotope zonation might affect the diffusive loss of a radiogenic isotope if this process is limited by an adjacent mineral. Our MATLAB scripts simulate radiogenic ingrowth and diffusion of Pb in spheres and finite cylinders with heterogeneous distributions of U and Th, and we use these scripts to investigate if one-dimensional spherical models can accurately predict the U-Th-Pb dates of chemically zoned non-spherical apatite inclusions.

Method
We have developed two MATLAB scripts for numerical modelling of radiogenic ingrowth and diffusion of Pb in U-Th-bearing inclusions of one mineral within another mineral. The first script is called spherincl.m and performs 1-dimensional simulations for a symmetric sphere. The second script is called cylincl.m and performs 2-dimensonal simulations for an axisymmetric finite cylinder. Both scripts utilise the Crank-Nicolson method on a regularly spaced grid and were programmed following Crank's book on diffusion [50]. Both scripts are provided in the supplementary archive scripts.zip, which also includes detailed notes on the equations used in the scripts and example input files for cylincl.m (spherincl.m does not need input files).
The key feature of our scripts is that they simulate radiogenic ingrowth and diffusion of Pb in a pair of juxtaposed minerals as shown in Figure 1. The mineral at the centre of the model, the inclusion, is completely surrounded by the other mineral, the host. At the interface these minerals are assumed to be in equilibrium according to the partition coefficient Kd, which remains constant through time. Any Pb that leaves the inclusion is assumed to diffuse through the host towards the outer boundary of the model without accumulating at the interface. Thus, we assume that there is no significant segregation of Pb into the grain boundary between the inclusion and the host, implying that one or both of these minerals do not have a strong tendency to expel Pb from their structure. The outer boundary of the model, which is the surface of the modelled composite sphere or finite cylinder, is assumed to completely lose any Pb that reaches it by diffusion and accumulate any Pb that is produced at it by decay. Thus, we assume that the volume of the inclusion is considerably smaller than the volume of the host, such that all of the Pb that escapes the inclusion can be absorbed by the host without significantly changing its average Pb concentration. As required by symmetry, we assume that there is no Pb flux across the inner boundary of the model, which is the centre of the modelled composite sphere or the central axis of the modelled composite finite cylinder. The diffusivity of Pb within each mineral is assumed to be independent of any intensive parameters other than temperature. U and Th are assumed to be immobile, since U and Th diffusivities are significantly lower than Pb diffusivity in apatite and in other minerals, for which there is experimental data [51].
is considerably smaller than the volume of the host, such that all of the Pb that escapes th inclusion can be absorbed by the host without significantly changing its average Pb con centration. As required by symmetry, we assume that there is no Pb flux across the inne boundary of the model, which is the centre of the modelled composite sphere or the cen tral axis of the modelled composite finite cylinder. The diffusivity of Pb within each min eral is assumed to be independent of any intensive parameters other than temperature. U and Th are assumed to be immobile, since U and Th diffusivities are significantly lowe than Pb diffusivity in apatite and in other minerals, for which there is experimental dat [51].
To evaluate the performance of our scripts we calculated the closure temperatur profiles for chemically homogeneous spherical apatite inclusions that are open to com plete loss of radiogenic Pb at the crystal boundary and compared these with the predic tions of Dodson's equation [2]. The diffusion properties of apatite were taken from [51 while the diffusion properties of the host mineral were arbitrarily chosen to ensure tha there is no accumulation of Pb next to the boundary with apatite (i.e., to mimic the absenc of the host mineral). The inclusions had radii of 150 μm and cooled at rates of 100 °C·Ma − 10 °C·Ma −1 and 1 °C·Ma −1 following Dodson's definition. The closure temperature profile that are calculated using our scripts are within 3 % of those calculated using Dodson' equation ( Figure 2).

Figure 1.
A graphical explanation for the chosen approach to simulate diffusive loss of radiogenic Pb from apatite inclusions hosted by other minerals. The left side shows a composite sphere modelled by spherincl.m and a composite finite cylinder modelled by cylincl.m. The right side provides a simplified description of how Pb diffusion across the apatite-host mineral interface was modelled. Black squares indicate the final three nodes in apatite, while circles show the first three nodes in the host mineral. C: concentration; dx: distance between the nodes; dt: duration of a time step; D: diffusion coefficient; J: diffusive flux; P: production by radioactive decay; Kd: partition coefficient; subscripts Ap and H designate apatite and the host mineral, respectively. Black squares indicate the final three nodes in apatite, while circles show the first three nodes in the host mineral. C: concentration; dx: distance between the nodes; dt: duration of a time step; D: diffusion coefficient; J: diffusive flux; P: production by radioactive decay; Kd: partition coefficient; subscripts Ap and H designate apatite and the host mineral, respectively.
To evaluate the performance of our scripts we calculated the closure temperature profiles for chemically homogeneous spherical apatite inclusions that are open to complete loss of radiogenic Pb at the crystal boundary and compared these with the predictions of Dodson's equation [2]. The diffusion properties of apatite were taken from [51], while the diffusion properties of the host mineral were arbitrarily chosen to ensure that there is no accumulation of Pb next to the boundary with apatite (i.e., to mimic the absence of the host mineral). The inclusions had radii of 150 µm and cooled at rates of 100 • C·Ma −1 , 10 • C·Ma −1 and 1 • C·Ma −1 following Dodson's definition. The closure temperature profiles that are calculated using our scripts are within 3 % of those calculated using Dodson's equation (Figure 2). Here we investigate whether or not the U-Th-Pb dates of apatite inclusions can be affected by their host minerals. We do this by comparing the temperatures at which chemically homogeneous spherical apatite crystals stop leaking Pb by diffusion when surrounded by an infinite reservoir (Tc Dodson) and when encapsulated within crystals of other minerals (Tc Inclusion). Tc Dodson were calculated using Dodson's equation [1]. Tc Inclusion were determined using sperincl.m by calculating U-Th-Pb dates of apatite inclusions in different minerals that cool along the t-T paths assumed for them by Dodson's equation and then identifying the temperatures that correspond to these dates. Pb diffusion parameters for apatite were taken from [51], while Pb diffusion parameters for the host phase were defined as a fraction or multiple of those for apatite. Calculations were run for hundreds of hypothetical host minerals to characterise the difference between Tc Inclusion and Tc Dodson as a function of ratios between the activation energies (Ea Host /Ea Apatite ) and the preexponential factors (log(D0 Host /D0 Apatite )) for Pb diffusion in apatite and the host mineral ( Figure 3a). The dependence of Tc Inclusion-Tc Dodson on diffusion properties of the host mineral was characterised for different inclusion radii (r = 150 and 50 μm), different partition coefficients (Kd = 0.01, 0.2, 5 and 100) and different cooling rates (CR = 50 °C Ma −1 and 0.5 °C Ma −1 ). Figure 3a shows the typical relationship between Tc Inclusion-Tc Dodson and diffusion properties of the host mineral. Two regions can be defined in this figure. The first is a region of low Ea Host /Ea Apatite and high log(D0 Host /D0 Apatite ) where Tc Inclusion = Tc Dodson, which we refer to hereafter as a field of open boundary behaviour. The second is a region of high Ea Host /Ea Apatite and low log(D0 Host /D0 Apatite ) where Tc Inclusion > Tc Dodson, which is referred to hereafter as a field of flux-limited boundary behaviour. The boundary between these regions is linear. The intercept of the boundary line with the log(D0 Host /D0 Apatite ) axis strongly depends on Kd, while its slope is nearly independent of Kd ( Figure 3b). Neither r nor CR have a significant effect on the slope and the intercept of the boundary line ( Figure 3c). The rate of increase of Tc Inclusion-Tc Dodson with increasing Ea Host /Ea Apatite is likewise nearly independent of r and CR (Figure 3d).
Our modelling indicates that the open boundary behaviour required for Tc Dodson to be accurate for apatite inclusions can only be attained if Pb diffusion in the host mineral is sufficiently rapid and Kd is sufficiently high. To see if this is likely to occur in natural setting we have plotted in Figure 3a,b the available experimental Pb diffusivity data for different minerals [51][52][53][54]. All of these minerals fall into the field of flux-limited boundary behaviour when Kd = 5, and nearly all of them remain in the same field when Kd = 100. A comparison of the closure temperature profiles calculated for spherical crystals with radii of 150 µm using our MATLAB scripts and Dodson's equation [2]. (a,b) Absolute values for closure temperatures obtained using spherincl.m and cylincl.m, respectively. (c,d) The ratios between the theoretical values for closure temperatures calculated using Dodson's equation [2] and those obtained using our MATLAB scripts. CR: cooling rate. Here we investigate whether or not the U-Th-Pb dates of apatite inclusions can be affected by their host minerals. We do this by comparing the temperatures at which chemically homogeneous spherical apatite crystals stop leaking Pb by diffusion when surrounded by an infinite reservoir (Tc Dodson ) and when encapsulated within crystals of other minerals (Tc Inclusion ). Tc Dodson were calculated using Dodson's equation [1]. Tc Inclusion were determined using sperincl.m by calculating U-Th-Pb dates of apatite inclusions in different minerals that cool along the t-T paths assumed for them by Dodson's equation and then identifying the temperatures that correspond to these dates. Pb diffusion parameters for apatite were taken from [51], while Pb diffusion parameters for the host phase were defined as a fraction or multiple of those for apatite. Calculations were run for hundreds of hypothetical host minerals to characterise the difference between Tc Inclusion and Tc Dodson as a function of ratios between the activation energies (E a Host /E a Apatite ) and the preexponential factors (log(D 0 Host /D 0 Apatite )) for Pb diffusion in apatite and the host mineral ( Figure 3a). The dependence of Tc Inclusion -Tc Dodson on diffusion properties of the host mineral was characterised for different inclusion radii (r = 150 and 50 µm), different partition coefficients (Kd = 0.01, 0.2, 5 and 100) and different cooling rates (CR = 50 • C Ma −1 and 0.5 • C Ma −1 ). Figure 3a shows the typical relationship between Tc Inclusion -Tc Dodson and diffusion properties of the host mineral. Two regions can be defined in this figure. The first is a region of low E a Host /E a Apatite and high log(D 0 Host /D 0 Apatite ) where Tc Inclusion = Tc Dodson , which we refer to hereafter as a field of open boundary behaviour. The second is a region of high E a Host /E a Apatite and low log(D 0 Host /D 0 Apatite ) where Tc Inclusion > Tc Dodson , which is referred to hereafter as a field of flux-limited boundary behaviour. The boundary between these regions is linear. The intercept of the boundary line with the log(D 0 Host /D 0 Apatite ) axis strongly depends on Kd, while its slope is nearly independent of Kd ( Figure 3b). Neither r nor CR have a significant effect on the slope and the intercept of the boundary line ( Figure 3c). The rate of increase of Tc Inclusion -Tc Dodson with increasing E a Host /E a Apatite is likewise nearly independent of r and CR ( Figure 3d).

Results and Discussion
would provide direct estimates for Kd at relevant physicochemical conditions. However, similar partition coefficients in the order of 10 −1 to 10 0 have been obtained for apatite-melt and feldspar-melt pairs [55][56][57][58][59][60], and thus it seems reasonable to assume that the upper limit for Kd is close to 5. If this is accurate, then feldspar-hosted apatite inclusions lie well inside the field of flux-limited boundary behaviour. Consequently, we conclude that host minerals are likely to inhibit diffusive Pb loss from apatite inclusions and thus significantly influence their U-Th-Pb dates. Our modelling indicates that the open boundary behaviour required for Tc Dodson to be accurate for apatite inclusions can only be attained if Pb diffusion in the host mineral is sufficiently rapid and Kd is sufficiently high. To see if this is likely to occur in natural setting we have plotted in Figure 3a,b the available experimental Pb diffusivity data for different minerals [51][52][53][54]. All of these minerals fall into the field of flux-limited boundary behaviour when Kd = 5, and nearly all of them remain in the same field when Kd = 100. Some minerals, such as zircon, fall very far from the field of open boundary behaviour even when Kd is clearly unrealistically high, and the U-Th-Pb dates of apatite inclusions in them reflect the crystallisation ages for the modelled t-T paths. This parallels the observation of Antoine et al., who present U-Pb dates from zircon-hosted apatite inclusions that record their age of crystallisation, despite experiencing a much younger metamorphic event at temperatures (575-745 • C) that exceed Tc Dodson for apatite (350-550 • C) [24].
Other minerals, such as feldspars, titanite and rutile, fall close to the field of open system behaviour when Kd = 100, and it is not so clear whether this Kd value is unrealistically high for them. Out of these the most probable host minerals for apatite are feldspars, including both plagioclase and alkali feldspar (others are typically accessory phases). We are not aware of any systematic studies of Pb partitioning of between apatite and feldspars that would provide direct estimates for Kd at relevant physicochemical conditions. However, similar partition coefficients in the order of 10 −1 to 10 0 have been obtained for apatite-melt and feldspar-melt pairs [55][56][57][58][59][60], and thus it seems reasonable to assume that the upper limit for Kd is close to 5. If this is accurate, then feldspar-hosted apatite inclusions lie well inside the field of flux-limited boundary behaviour. Consequently, we conclude that host minerals are likely to inhibit diffusive Pb loss from apatite inclusions and thus significantly influence their U-Th-Pb dates.

Inversion Modelling: Is the Knowledge of Boundary Behaviour Needed to Infer t-T Paths?
Here we explore the effect of making incorrect assumptions about the boundary behaviour of apatite inclusions during thermochronological modelling of their U-Th-Pb dates and whether this effect depends on the topology of the true thermal history. We first used forward modelling in spherincl.m to generate two sets of synthetic bulkgrain 206 Pb/ 238 U dates for chemically homogeneous spherical apatite inclusions in alkali feldspar. The first set was produced for an arbitrarily chosen monotonic cooling history (Figure 4a), while the second set was created for a more complex thermal history that includes a reheating event (Figure 4b), which was also arbitrarily chosen. Each of the two sets included 206 Pb/ 238 U dates of 6 inclusions with radii of 10, 20, 30, 50, 100 and 150 µm. Kd was kept at a value of 5 following the same rationale based on the mineral-melt partition coefficients from [57][58][59][60] as above. The diffusion parameters for Pb in apatite and alkali feldspar were taken from [51,52]. The obtained 206 Pb/ 238 U dates were assumed to have 1σ uncertainties of 1% (Table 1). We subsequently attempted to reconstruct the original t-T path for each set of 206 Pb/ 238 U dates by inversion modelling using HeFTy version 1.9.3.74 [5], which incorrectly assumed that the apatite boundaries are open to complete Pb loss (the only option available when using HeFTy). During the inversion modelling we used a priori t-T constraints to reduce the time needed to find 100 good-fit t-T paths with a goodness-of-fit parameter p > 0.5. Two constraints were introduced when inverting the apatite 206 Pb/ 238 U dates obtained for the monotonic cooling history, which are crystallisation at 1000 • C and 125 Ma and a broad t-T field between 0 and 400 • C and 0 and 125 Ma. The inversion gave 100 good-fit t-T paths after 3565 trials (Figure 4a). The 206 Pb/ 238 U dates obtained for the thermal history with reheating were inverted using six constraints, which are crystallisation at 770 • C and 255 Ma and five t-T boxes that were obtained from a series of preliminary calculations with inflated uncertainties. The latter five constraints allowed to find~100 good-fit t-T paths after~8,350,000 trials using two computers that ran simultaneously (Figure 4b). Examples of HeFTy files with input data and results of inversion modelling are provided in the supplementary archive inversion.zip.
The good-fit t-T paths inferred for the apatite inclusions that cooled monotonically have similar topology to the original thermal history, although they suggest that cooling to <350 • C occurred up to several Ma earlier than in the original thermal history (Figure 4a). In contrast, the good-fit t-T paths obtained for the apatite inclusions that experienced reheating are strikingly different from the original thermal history and do not reveal the simulated reheating event (Figure 4b). This is not a result of using tight a priori constraints because t-T solutions with reheating to >400 • C at~95 Ma are also not found when larger uncertainties and broader t-T boxes are used (a figure illustrating this is provided in the supplementary archive inversion.zip). Clearly, an incorrect assumption about the openness of the apatite boundaries to diffusive Pb loss can lead to profound inaccuracies in thermochronological reconstructions that were obtained by U-Th-Pb dating of apatite inclusions, particularly when dealing with complex thermal histories that feature reheating to temperatures comparable with Tc Dodson of apatite. The reconstructed t-T paths can have different topology compared to the original thermal histories, and thus even the amount of cooling or reheating that they suggest between two arbitrary points in time can be significantly inaccurate. Hence, we conclude that it is critically important to make appropriate assumptions about the type of boundary behaviour during thermochronological modelling.

Chemically Zoned Spheres: What Is the Effect of U and Th Zonation on Bulk-Grain U-Th-Pb Dates?
Here we assess how bulk-grain U-Th-Pb dates of apatite crystals with flux-limited boundaries are affected by heterogeneous intra-grain distributions of U and Th. Previous work showed that intra-grain variations in U and Th concentration should be accounted for when apatite boundaries are fully open with respect to diffusive Pb loss [11]. However, is this also true when apatite boundaries are flux-limited? We address this question by comparing synthetic bulk-grain 206 Pb/ 238 U dates within three groups of spherical apatite crystals with variable U zoning. Each group contained six chemically homogeneous crystals, six crystals with elevated U concentration at the rims and six crystals with U-enriched cores (Figure 5a). The crystals had radii of 10, 20, 30, 50, 100 and 150 μm. The first group of crystals was assumed to be completely open to diffusive Pb loss at the boundary (Kd = ∞). Crystals in the second and third groups were assumed to have alkali feldspar at the boundary, and Kd was fixed at 5 following the rationale based on the mineral-melt partition coefficients from [57][58][59][60] that is described above. The first and second groups followed a "cooler" t-T path (Figure 5b), while the third group followed a "hotter" t-T path (Figure 5c). The 'cooler' t-T path was arbitrarily chosen, while the "hotter" t-T path was created by adding 53 °C to each point in the 'cooler' one to ensure that similar fractional Pb loss is achieved by apatite crystals in the first and third groups. The diffusion parameters for Pb in apatite and alkali feldspar were taken from [51,52], and calculations were made using spherincl.m.  Here we assess how bulk-grain U-Th-Pb dates of apatite crystals with flux-limited boundaries are affected by heterogeneous intra-grain distributions of U and Th. Previous work showed that intra-grain variations in U and Th concentration should be accounted for when apatite boundaries are fully open with respect to diffusive Pb loss [11]. However, is this also true when apatite boundaries are flux-limited? We address this question by comparing synthetic bulk-grain 206 Pb/ 238 U dates within three groups of spherical apatite crystals with variable U zoning. Each group contained six chemically homogeneous crystals, six crystals with elevated U concentration at the rims and six crystals with Uenriched cores (Figure 5a). The crystals had radii of 10, 20, 30, 50, 100 and 150 µm. The first group of crystals was assumed to be completely open to diffusive Pb loss at the boundary (Kd = ∞). Crystals in the second and third groups were assumed to have alkali feldspar at the boundary, and Kd was fixed at 5 following the rationale based on the mineral-melt partition coefficients from [57][58][59][60] that is described above. The first and second groups followed a "cooler" t-T path (Figure 5b), while the third group followed a "hotter" t-T path (Figure 5c). The 'cooler' t-T path was arbitrarily chosen, while the "hotter" t-T path was created by adding 53 • C to each point in the 'cooler' one to ensure that similar fractional Pb loss is achieved by apatite crystals in the first and third groups. The diffusion parameters for Pb in apatite and alkali feldspar were taken from [51,52], and calculations were made using spherincl.m. (e,f) 206 Pb/ 238 U dates of alkali feldspar-hosted apatite inclusions that followed the "cooler" and "hotter" thermal histories, respectively (the second and third groups, respectively).

Spheres and Spherocylinders: What Is the Effect of Crystal Geometry on Intra-Grain U-Th-Pb Dates?
Here we investigate the effect of crystal geometry on intra-grain U-Th-Pb date variations in apatite and whether it depends on the U and Th zoning pattern and the type of boundary behaviour. Previous work has shown that crystal geometry influences the bulkgrain isotopic dates of minerals [43][44][45][46][47][48][49], but its effect on intra-grain variations of isotopic dates has not been investigated. Therefore, it remains unclear how significant this effect is and whether it depends on the intra-grain distribution of parent isotopes and the degree of boundary openness with respect to diffusive Pb loss. We address this issue by comparing synthetic intra-grain 206 Pb/ 238 U date maps and core to rim 206 Pb/ 238 U date profiles in 4 groups of apatite crystals, where each group contains three crystals that have different geometries but identical boundary behaviour and character of intra-grain U distribution. The groups differ from each other by boundary behaviour and character of intra-grain U distribution of crystals within them.
All of the crystals followed the same arbitrarily chosen thermal history that is shown in Figure 6a. The crystals within each group had one of the following geometries ( Figure  7): an elongated shape that is similar to a spherocylinder, a sphere whose radius is equal to the half-width of the elongated shape (referred to hereafter as r-normalised), or a sphere whose radius is 1.4 times longer than the half-width of the elongated shape, such that it has the same surface area (S) to volume (V) ratio as the elongated shape (referred to hereafter as S/V-normalised). The use of r-normalised spherical crystals is consistent with how previous studies modelled intra-grain U-Pb date variations in apatite [9,11]. The use of S/V-normalised spherical crystals copies how previous studies modelled intra-grain U-Pb date variations in rutile [13,14] and also follows the methodology of Meesters and Dunai to account for crystal geometry when modelling bulk-grain isotopic dates by using an equivalent spherical crystal [43,44]. Two of the four groups of crystals were assumed to be chemically homogeneous, while the remaining two groups had heterogeneous distri- (e,f) 206 Pb/ 238 U dates of alkali feldspar-hosted apatite inclusions that followed the "cooler" and "hotter" thermal histories, respectively (the second and third groups, respectively).
The 206 Pb/ 238 U dates of apatite crystals from the first group (Kd = ∞) are significantly scattered about the trend defined by chemically homogeneous crystals, and crystals with U-enriched rims appear to be younger by up to~50 Ma than crystals with identical radii and elevated U concentrations in the cores (Figure 5d). This observation is consistent with previous studies that examined the effect of parent isotope zonation on isotopic dates of apatite and other minerals [11,44]. The 206 Pb/ 238 U dates of apatite inclusions from the second group (Kd = 5) are less scattered, and the maximum difference between crystals of the same size is 17 Ma (Figure 5e). These inclusions are considerably older than those from the first group, so it may seem that the reduction in scatter is solely related to the reduction in fractional Pb loss. However, if this was the only factor at play, a greater scatter would be observed for the third group of apatite inclusions (Kd = 5), which yield similar 206 Pb/ 238 U dates to apatite crystals in the first group and thus have lost similar fractions of Pb (Figure 5f). Yet, their 206 Pb/ 238 U dates do not scatter at all about the trend defined by chemically homogeneous inclusions. This suggests that there is another factor that contributed to the reduction in scatter, which is the presence of alkali feldspar at the boundaries of the apatite inclusions.
Overall, our results can be explained by invoking two competing processes to describe the loss of radiogenic Pb that is produced in apatite near its boundary, which are the transfer across the boundary and the diffusion towards the core. In case of open boundary behaviour, the transfer across the boundary is infinitely faster than the diffusion towards the core, so the former process dominates. Consequently, apatite crystals with U-Th-enriched rims lose greater fractions of Pb than apatite crystals with U-Th-enriched cores, which gives rise to scattered bulk-grain U-Th-Pb dates. In contrast, in case of flux-limited boundary behaviour, the transfer across the boundary becomes slower than the diffusion towards the core. Therefore, the latter process dominates, particularly at temperatures exceeding Tc Dodson of apatite, which are necessary for achieving large fractional Pb loss. Consequently, apatite inclusions with variable U and Th zoning partially or completely homogenise with respect to radiogenic Pb before they lose it, which reduces or eliminates the scatter in bulkgrain U-Th-Pb dates. We therefore conclude that the effect of heterogeneous distributions of U and Th in apatite crystals on bulk-grain U-Th-Pb dates is relatively minor to insignificant when their boundaries are flux-limited.
Our results suggest that it may be possible to identify the type of boundary behaviour by observing whether the bulk-grain U-Th-Pb dates of apatite crystals with variable U and Th zoning are scattered. However, this will only work when the rock contains apatite and one other mineral. Apatite crystals in polymineralic rocks may share their boundaries with varying mineral phases, each with different diffusion properties and ability to incorporate Pb into structure, and hence with different potential to inhibit diffusive loss of Pb. This will likely result in scattered bulk-grain U-Th-Pb dates for any given crystal radius. Therefore, any scatter in bulk-grain U-Th-Pb dates of chemically zoned apatite crystals form a polymineralic rock, such as that observed in some natural samples from the Northern Andes [11], does not automatically imply that their boundaries are completely open to Pb loss.

Spheres and Spherocylinders: What Is the Effect of Crystal Geometry on Intra-Grain U-Th-Pb Dates?
Here we investigate the effect of crystal geometry on intra-grain U-Th-Pb date variations in apatite and whether it depends on the U and Th zoning pattern and the type of boundary behaviour. Previous work has shown that crystal geometry influences the bulk-grain isotopic dates of minerals [43][44][45][46][47][48][49], but its effect on intra-grain variations of isotopic dates has not been investigated. Therefore, it remains unclear how significant this effect is and whether it depends on the intra-grain distribution of parent isotopes and the degree of boundary openness with respect to diffusive Pb loss. We address this issue by comparing synthetic intra-grain 206 Pb/ 238 U date maps and core to rim 206 Pb/ 238 U date profiles in 4 groups of apatite crystals, where each group contains three crystals that have different geometries but identical boundary behaviour and character of intra-grain U distribution. The groups differ from each other by boundary behaviour and character of intra-grain U distribution of crystals within them.
All of the crystals followed the same arbitrarily chosen thermal history that is shown in Figure 6a. The crystals within each group had one of the following geometries ( Figure 7): an elongated shape that is similar to a spherocylinder, a sphere whose radius is equal to the half-width of the elongated shape (referred to hereafter as r-normalised), or a sphere whose radius is 1.4 times longer than the half-width of the elongated shape, such that it has the same surface area (S) to volume (V) ratio as the elongated shape (referred to hereafter as S/V-normalised). The use of r-normalised spherical crystals is consistent with how previous studies modelled intra-grain U-Pb date variations in apatite [9,11]. The use of S/V-normalised spherical crystals copies how previous studies modelled intra-grain U-Pb date variations in rutile [13,14] and also follows the methodology of Meesters and Dunai to account for crystal geometry when modelling bulk-grain isotopic dates by using an equivalent spherical crystal [43,44]. Two of the four groups of crystals were assumed to be chemically homogeneous, while the remaining two groups had heterogeneous distributions of U. The U concentration maps for elongated and r-normalised spherical crystals are shown in Figure 6b, while the U concentration maps for S/V-normalised spherical crystals are obtained by multiplying by 1.4 the coordinates in those for r-normalised spherical crystals. Finally, one group of the chemically homogeneous crystals and one group of the chemically zoned crystals were assumed to have open boundaries (Kd = ∞), while the crystals in the remaining two groups were assumed to be surrounded by alkali feldspar with Kd = 5 following [57][58][59][60], as explained above. The diffusion properties for Pb in apatite and alkali feldspar were taken from [51,52]. The radiogenic ingrowth and diffusion of Pb in the elongated crystals was modelled using cylincl.m, while spherincl.m was used to model the spherical crystals. rounded by the same mineral. Furthermore, the youngest U-Th-Pb date in apatite crystals with open boundaries will always be at their outermost rim, while apatite crystals with flux-limited boundaries may yield their youngest U-Th-Pb dates in more interior zones. This latter finding is relevant to the discussion of Villa [18] on how to distinguish between diffusive and fluid-mediated loss of radiogenic isotopes. According to him, any increase of isotopic dates towards the rim of a dated crystal invalidates the assumption that radiogenic isotopes were predominantly lost by volume diffusion. On the contrary, our modelling indicates that such an increase of isotopic dates may be caused by a heterogeneous intra-grain distribution of parent isotopes in a crystal with flux-limited boundaries (an example of apatite where this could be the case was reported in [21]).

Further Considerations
We have demonstrated that the U-Th-Pb dates of apatite inclusions can be significantly affected by their host minerals. However, it is clear that not every apatite crystal in a rock is entirely encapsulated within a crystal of another mineral. In our experience, most of them are located between several crystals of different minerals and connect to extensive networks of grain boundaries. Such networks are traditionally assumed, explicitly or implicitly, to provide pathways for the efficient removal of radiogenic isotopes to an infinite reservoir or to sink minerals [38,39]. This assumption is based on the fact that diffusion A comparison of the 206 Pb/ 238 U date maps and core to rim 206 Pb/ 238 U date profiles of all four groups (Figure 6d-g) reveals that the two groups of apatite inclusions in alkali feldspar (Kd = 5) have significantly older 206 Pb/ 238 U dates than the two groups of apatite crystals that have open boundaries (Kd = ∞), which is consistent with our previous observations. The difference between the intra-grain 206 Pb/ 238 U dates of any two crystals with the same geometry and U zonation but different boundary behaviour reaches hundreds of Ma (hundreds of %). Notably, some intra-grain 206 Pb/ 238 U dates in the zoned crystals exceed the assumed age of crystallisation by up to hundreds of Ma (tens of %), particularly in those with flux-limited boundaries. A more subtle observation is that all of the crystals with open boundaries have identical 206 Pb/ 238 U dates at their boundaries, while different inclusions in alkali feldspar yield different 206 Pb/ 238 U dates at their boundaries. Finally, it is notable that the youngest 206 Pb/ 238 U date in the crystals with open boundaries always occurs at their boundary, while this is not always true for the apatite inclusions within alkali feldspar. A more detailed look at each individual group shows that the bulk-grain and most of the intra-grain 206 Pb/ 238 U dates of the r-normalised spherical crystals are noticeably younger than those of the elongated crystals. Both types of dates may differ by up to tens of Ma (tens of %) between crystals that have the same core to rim U concentration profiles and degree of boundary openness. This effect is less pronounced in the chemically homogeneous crystals than in the crystals with heterogeneous distributions of U, and it is also less pronounced in the alkali feldspar-hosted inclusions than in the crystals with open boundaries. The bulk-grain 206 Pb/ 238 U dates of S/V-normalised spherical crystals are only slightly older than those yielded by the elongated crystals with identical boundary behaviour and core to rim U concentration profile. The difference between these stayed below 10 Ma (10%) in our calculations. However, intra-grain 206 Pb/ 238 U dates of the same pair of spherical and elongated crystals at the same normalised distance from the core may differ by tens of Ma (tens of %) in either direction.
A comparison of the differences in bulk-grain and intra-grain 206 Pb/ 238 U dates caused by variations in crystal geometry within individual groups shows that the effect of crystal geometry on diffusive Pb loss from apatite is significant. This effect becomes particularly pronounced in crystals with open boundaries and in crystals that have a heterogeneous intra-grain distribution of U and Th. Neither bulk-grain nor intra-grain U-Th-Pb dates of non-spherical apatite crystals can be reliably modelled using r-normalised spherical crystals. The use of S/V-normalised spherical crystals following the approach of Meesters and Dunai [43,44] is justified when modelling bulk-grain U-Th-Pb dates of apatite crystals that have non-spherical shapes. More accurate results could probably be achieved by using more sophisticated methods of accounting for crystal geometry, which were proposed by Farley et al. [47] and Huber et al. [48]. However, we show that intra-grain variations of U-Th-Pb dates in non-spherical crystals cannot be reliably modelled using S/V-normalised sphere equivalents, particularly when intra-grain distributions of U and Th are heterogeneous. Ideally, three-dimensional calculations are required to simulate intra-grain variations of U-Th-Pb dates in these cases, although two-dimensional calculations may be sufficient for axisymmetric crystals. Both options are computationally expensive and would considerably increase the machine time needed for thermochronological inversion modelling. This increase could be mitigated by first integrating spatial U-Th-Pb data to obtain bulk-grain U-Th-Pb dates and then modelling radiogenic ingrowth and diffusion of Pb in equivalent spheres to reproduce these dates. The resulting t-T solutions could subsequently be validated by forward modelling using more appropriate multi-dimensional calculations.
Our modelling also suggests that it is possible to distinguish between open and flux-limited boundary behaviour by examining intra-grain variations of U-Th-Pb dates in several apatite crystals from the same rock. Apatite crystals with open boundaries will always have identical U-Th-Pb dates at their boundaries, while apatite crystals with flux-limited boundaries can have different U-Th-Pb dates at their boundaries, even if they are all surrounded by the same mineral. Furthermore, the youngest U-Th-Pb date in apatite crystals with open boundaries will always be at their outermost rim, while apatite crystals with flux-limited boundaries may yield their youngest U-Th-Pb dates in more interior zones. This latter finding is relevant to the discussion of Villa [18] on how to distinguish between diffusive and fluid-mediated loss of radiogenic isotopes. According to him, any increase of isotopic dates towards the rim of a dated crystal invalidates the assumption that radiogenic isotopes were predominantly lost by volume diffusion. On the contrary, our modelling indicates that such an increase of isotopic dates may be caused by a heterogeneous intra-grain distribution of parent isotopes in a crystal with flux-limited boundaries (an example of apatite where this could be the case was reported in [21]).

Further Considerations
We have demonstrated that the U-Th-Pb dates of apatite inclusions can be significantly affected by their host minerals. However, it is clear that not every apatite crystal in a rock is entirely encapsulated within a crystal of another mineral. In our experience, most of them are located between several crystals of different minerals and connect to extensive networks of grain boundaries. Such networks are traditionally assumed, explicitly or implicitly, to provide pathways for the efficient removal of radiogenic isotopes to an infinite reservoir or to sink minerals [38,39]. This assumption is based on the fact that diffusion along grain boundaries is generally several orders of magnitude faster than diffusion within crystals [61,62]. However, grain boundary diffusion does not occur in isolation from intra-crystal diffusion and thus has limited capacity to influence the chemical transport in rocks, the implications of which for geomaterials have been reviewed by Joesten [61] and Dohmen and Milke [62]. Therefore, it is important to have an understanding of whether grain boundary diffusion can sufficiently enhance the removal of radiogenic Pb from apatite crystals that are not entirely encapsulated in other crystals before assuming open boundary behaviour for them.
As described in [61][62][63][64][65], diffusion of an isotope in a polycrystalline medium involves transport through both grain boundaries and crystal volumes. The relative importance of grain boundaries depends on many factors, including isotope partitioning into grain boundaries, diffusion properties of crystals and grain boundaries, crystal size, grain boundary width and timescales. Over short timescales, the isotope can travel along grain boundaries while diffusion through the surrounding crystals remains insignificant. This is kinetic regime C. As time progresses, significant quantities of the isotope start to diffuse from grain boundaries into the surrounding crystals, reducing the effectiveness of transport along grain boundaries and the length of diffusion profiles that can develop within them. Nevertheless, the isotope can still travel much further away from its source along grain boundaries than through crystal volumes. This is kinetic regime B. Eventually, the diffusion fronts in grain boundaries and within crystals catch up and start to advance at the same speed. This may be kinetic regime B4 or A. When the abundance of grain boundaries is low, such that the medium can be classified as coarse-grained [63,64], the total mass transfer in all kinetic regimes is dominated by diffusion through crystal volumes. Grain boundaries start to dominate the total mass transfer as their abundance increases, although the bulk diffusion coefficient for the medium in regime A remains a fraction of that of grain boundaries until the medium can be classified as special ultrafine-grained [63,64].
Unfortunately, most of the parameters required to characterise Pb diffusion though mineral matrices that can potentially surround apatite crystals remain unconstrained. For example, there are no published data regarding the partitioning of Pb into grain boundaries and Pb diffusivity along them for alkali feldspar, which we assumed to be the host mineral of apatite inclusions in most of our calculations. Grain boundary-mineral partition coefficients, also known as segregation factors, are in many cases similar to the inverse of mineral-melt partition coefficients [66]. Assuming that this relationship is valid for Pb in alkali feldspar, we can estimate the Pb segregation factor to be close to unity based on alkali feldspar-melt partitioning data [58][59][60]. A comparison of the available data on grain boundary [67] and volume [52] diffusion in feldspars shows that K and Ca diffusivities for grain boundaries (D GB ) exceed those for crystal volumes (D M ) by factors of 10 4 -10 8 , and perhaps it is reasonable to assume a similar range for Pb. Grain boundary widths are usually assumed to be 1 nm in the geoscience literature [61,62], and for simplicity we will assume that grain boundaries are parallel and regularly spaced. With these assumptions we can use the equations for characteristic times t' (C to B transition) and t" (B to B4 transition) and Harrison's condition for regime A (B/B4 to A transition) [63] to infer the kinetic regime that will likely be attained in polycrystalline alkali feldspar around a crystal of apatite. The impact of grain boundaries on the total mass transfer can be assessed by estimating the ratio of masses absorbed by grain boundaries (M GB ) and crystal volumes (M M ) in regimes C and B and at the effective diffusivity (D effective ) in regimes B4 and A, all of which depend on the size of the alkali feldspar crystals [63].
The obtained estimates for t', t" and Harrison's condition for regime A suggest that the transition from regime C to regime B should occur over any geologically relevant timescale, while the transition from regime B to regimes B4 or A may never occur, especially at 350 • C when apatite with open boundaries closes to diffusive Pb loss (Figure 8a). If D GB /D M = 10 4 and the spacing between grain boundaries lies between 50 µm and 1 mm, then grain boundary diffusion accounts for less than 10% of the total mass transfer ( Figure 8b) and D effective is only slightly greater that D M (Figure 8c). However, if D GB /D M = 10 8 , then grain boundaries can account for significantly more than 10% of the total mass transfer in alkali feldspar with similarly spaced grain boundaries (Figure 8b), while its D effective exceeds D M by two to three orders of magnitude (Figure 8c). A comparison of these estimates with our results from Section 4.1.1 suggests that the enhancement of bulk diffusion properties of polycrystalline alkali feldspar by grain boundaries may not be sufficient to ensure that apatite crystals within it become fully open with respect to diffusive Pb loss. We therefore conclude that the U-Th-Pb dates of apatite crystals that are not fully encapsulated in other crystals may be significantly affected by the surrounding mineral matrix. This potential effect needs to be carefully addressed in future work by experimentally constraining grain boundary diffusivities and segregation factors for Pb in various minerals and by developing software that is suitable for modelling diffusive Pb loss into a polycrystalline matrix that operates in regime B (our scripts are not). Furthermore, it is also important to examine how different mineral matrices can affect diffusive loss of radiogenic isotopes from other thermochronometers.
While a detailed investigation of Pb diffusion in polycrystalline matrices may reveal that it is appropriate to assume open boundary behaviour for apatite crystals connected to grain boundary networks, it is still important to recognise that apatite inclusions in other minerals are likely to have flux-limited boundaries. Apatite inclusions can potentially help with circumventing the increasingly recognised problem that fluid-induced dissolution-reprecipitation frequently dominates the redistribution of radiogenic isotopes in minerals [15][16][17][18][19][20][21], including Pb in apatite [22][23][24]. As highlighted by Antoine et al. for apatite [24] and Ault et al. for zircon [40], inclusions can be shielded from interaction with fluids by their host minerals and thereby preserve their U-Th-Pb systematics. Such inclusions, particularly if zoned with respect to U and Th, represent an ideal target for thermochronological studies, and our work provides a basis for interpreting their U-Th-Pb dates. Minerals 2021, 11, x FOR PEER REVIEW 15 of 18

Conclusions
Our numerical modelling shows that diffusive loss of radiogenic Pb from apatite inclusions can be significantly inhibited by their host minerals. The assumption that radiogenic Pb is completely and instantaneously removed to an infinite reservoir after reaching the boundaries of apatite inclusions is only valid for a limited array of potential host minerals that have high Pb diffusivities and host mineral-apatite Pb partition coefficients. The available Pb diffusivity and mineral-melt partitioning data for various minerals imply that most do not fall into this array and will limit the flux of radiogenic Pb across the boundaries of apatite inclusions. Diffusive Pb loss from apatite with flux-limited boundaries proceeds in a very different manner compared to the case of open boundary behaviour. A heterogeneous distribution of U and Th in such apatite inclusions has a relatively minor to insignificant effect on their bulk-grain U-Th-Pb dates and can cause an increase in intragrain U-Th-Pb dates towards the rim. We show that thermochronologists should make appropriate assumptions about the openness of apatite boundaries with respect to diffusive Pb loss when reconstructing t-T paths of rocks by inversion modelling of apatite U-Th-Pb data. We further show that it is important to account for apatite crystal geometry and U and Th zonation. The latter is challenging when dealing with intra-grain U-Th-Pb date variations, because these cannot be accurately predicted by modelling radiogenic ingrowth and diffusion of Pb in equivalent spherical crystals. Finally, apatite crystals that are not completely encapsulated in crystals of other minerals may also have flux-limited boundaries, although there is insufficient experimental data to make a fully constrained assessment. We conclude that the likelihood of flux-limited boundary behaviour should be carefully considered in future research on apatite U-Th-Pb thermochronology, as well as thermochronology using other methods of isotopic dating.

Conclusions
Our numerical modelling shows that diffusive loss of radiogenic Pb from apatite inclusions can be significantly inhibited by their host minerals. The assumption that radiogenic Pb is completely and instantaneously removed to an infinite reservoir after reaching the boundaries of apatite inclusions is only valid for a limited array of potential host minerals that have high Pb diffusivities and host mineral-apatite Pb partition coefficients. The available Pb diffusivity and mineral-melt partitioning data for various minerals imply that most do not fall into this array and will limit the flux of radiogenic Pb across the boundaries of apatite inclusions. Diffusive Pb loss from apatite with flux-limited boundaries proceeds in a very different manner compared to the case of open boundary behaviour. A heterogeneous distribution of U and Th in such apatite inclusions has a relatively minor to insignificant effect on their bulk-grain U-Th-Pb dates and can cause an increase in intra-grain U-Th-Pb dates towards the rim. We show that thermochronologists should make appropriate assumptions about the openness of apatite boundaries with respect to diffusive Pb loss when reconstructing t-T paths of rocks by inversion modelling of apatite U-Th-Pb data. We further show that it is important to account for apatite crystal geometry and U and Th zonation. The latter is challenging when dealing with intra-grain U-Th-Pb date variations, because these cannot be accurately predicted by modelling radiogenic ingrowth and diffusion of Pb in equivalent spherical crystals. Finally, apatite crystals that are not completely encapsulated in crystals of other minerals may also have flux-limited boundaries, although there is insufficient experimental data to make a fully constrained assessment. We conclude that the likelihood of flux-limited boundary behaviour should be carefully considered in future research on apatite U-Th-Pb thermochronology, as well as thermochronology using other methods of isotopic dating.