Next Article in Journal
Mapping Potentially Acid Generating Material on Abandoned Mine Lands Using Remotely Piloted Aerial Systems
Next Article in Special Issue
Thermochronology of Alkali Feldspar and Muscovite at T > 150 °C Using the 40Ar/39Ar Method: A Review
Previous Article in Journal
Editorial for Special Issue “Mineralogy of Meteorites”
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Minerals
Volume 11
Issue 4
10.3390/min11040364
minerals-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Daniil V. Popov
1,2,* and
Richard A. Spikings
1
1
Department of Earth Sciences, University of Geneva, 13 Rue des Maraichers, CH-1205 Geneva, Switzerland
2
School of Environment, Earth and Ecosystem Sciences, The Open University, Walton Hall, Milton Keynes MK7 6AA, UK
*
Author to whom correspondence should be addressed.
Minerals 2021, 11(4), 364; https://doi.org/10.3390/min11040364
Submission received: 11 February 2021 / Revised: 26 March 2021 / Accepted: 28 March 2021 / Published: 31 March 2021
(This article belongs to the Special Issue Thermochronology at Temperatures Higher than 150 °C)
Browse Figures
Versions Notes

Abstract

:
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 infinite reservoir, (2) have a simple geometry and (3) are chemically homogeneous. Here we explore the significance 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 flux across their boundaries, and thus the thermal histories that are reconstructed assuming a fully open boundary may be significantly inaccurate, precluding a meaningful interpretation. We also find that when apatite boundaries are flux-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.

1. 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 geometry on intra-grain variations of U-Th-Pb dates in apatite inclusions and in apatite crystals whose boundaries are open to complete loss of Pb. Throughout this work we adhere to the assumption that the redistribution of radiogenic isotopes in minerals is dominated by volume diffusion, although we acknowledge that the role of fluid-induced dissolution-reprecipitation is frequently more important, including in case of Pb in apatite [15,16,17,18,19,20,21,22,23,24].

2. Previous Work

2.1. Behaviour of Radiogenic Isotopes at Crystal Boundaries

The assumption that radiogenic isotopes are completely lost to an infinite reservoir once they reach the boundaries of the dated crystal is the norm in thermochronological literature. This assumption is intrinsic to Dodson’s equations for the calculation of closure temperature (Tc Dodson) [1,2], as well as to their extensions that were developed by Ganguly and Tirone [3] and Guralnik et al. [4] and the related equations derived by Albarède [25] and Gardés and Montel [26]. This assumption is also intrinsic to software packages such as HeFTy [5] and QtQt [6], which are widely used to reconstruct the t-T paths of rocks by inversion modelling of isotopic dates (and other data). This assumption was made by the majority of previous thermochronological studies that utilised apatite U-Th-Pb dating [7,8,9,10,11], and other isotopic dating methods [12,13,14]. However, the validity of this assumption has not been systematically tested, and some studies propose different behaviour of radiogenic isotopes at crystal boundaries (see below).
A number of studies document anomalously old 40Ar/39Ar dates of micas and feldspars from metamorphic rocks and suggest that these reflect 40Ar uptake from grain boundary fluids or from other minerals during metamorphism [27,28,29,30,31,32,33,34]. Baxter proposed a mathematical model, in which the amount of 40Ar that can diffuse out of K-bearing minerals is limited by the 40Ar storage capacity of the surrounding minerals and the efficiency of 40Ar removal to an external sink [35]. This model assumes that on the scale of individual crystals every mineral in a rock remains in equilibrium with the grain boundary medium, and 40Ar is distributed between the minerals and grain boundary medium in accordance with equilibrium partition coefficients. On a larger scale 40Ar is assumed to diffuse along grain boundaries to some distant external sink (infinite reservoir) such as a shear zone.
Some studies report incomplete resetting of the Rb-Sr systematics of biotite during metamorphism, which they relate to the inability of biotite to re-equilibrate its Sr isotopic composition due to sluggish Sr diffusion through the surrounding minerals [36,37]. Giletti theorised that in any given polymineralic rock the mineral with the lowest Tc Dodson should be effectively closed with respect to diffusive re-equilibration of the Sr isotopic composition at the same temperature as the mineral with the second lowest Tc Dodson [38]. Jenkin et al. explored the idea of Giletti by numerical modelling of diffusive re-equilibration of Sr isotopic composition in a bi-mineralic rock, where the two constituent minerals were assumed to be in instantaneous chemical equilibrium with each other at the crystal boundaries [39]. They concluded that Giletti’s suggestion is only accurate when the modal abundance of the mineral with the lowest Tc Dodson approaches 100%. If the abundance of this mineral is lower, the temperature at which it becomes closed is lower than the Tc Dodson of the other mineral and tends to its own Tc Dodson as its modal abundance approaches 0%. The latter implies that Giletti’s suggestion is irrelevant for accessory minerals.
Several studies illustrate the potential inhibition of resetting U-Pb dates of a mineral due to it being surrounded by other minerals. Ault et al. showed that inclusions of radiation-damaged zircon in quartz are shielded from Pb loss by their host [40]. Antoine et al. showed that U-Pb dates of some apatite inclusions in zircon reflect their crystallisation age, while the U-Pb systematics of groundmass apatite crystals were completely reset by a later metamorphic event at temperatures that exceeded Tc Dodson of apatite [24]. Ibanez-Mejia et al. report U-Pb dates of apatite that are older than expected from independent t-T constraints [41]. They conclude that Pb diffusion in apatite was slower than suggested by experimental data, although they acknowledge the possibility that the removal of Pb from apatite was hampered by the surrounding minerals. Finally, Kohn et al. 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.

2.2. 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.

3. 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].
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).

4. Results and Discussion

4.1. Chemically Homogeneous Spheres

4.1.1. Forward Modelling: Do U-Th-Pb Dates of Apatite Inclusions Depend on Their Host Minerals?

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 (EaHost/EaApatite) and the preexponential factors (log(D0Host/D0Apatite)) 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 EaHost/EaApatite and high log(D0Host/D0Apatite) where Tc Inclusion = Tc Dodson, which we refer to hereafter as a field of open boundary behaviour. The second is a region of high EaHost/EaApatite and low log(D0Host/D0Apatite) 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(D0Host/D0Apatite) 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 EaHost/EaApatite 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. 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 100 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.

4.1.2. 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 bulk-grain 206Pb/238U 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 206Pb/238U 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 206Pb/238U 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 206Pb/238U 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 206Pb/238U 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 206Pb/238U 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.

4.2. 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 206Pb/238U 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.
The 206Pb/238U 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 206Pb/238U 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 206Pb/238U dates to apatite crystals in the first group and thus have lost similar fractions of Pb (Figure 5f). Yet, their 206Pb/238U 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 bulk-grain 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.

4.3. 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 206Pb/238U date maps and core to rim 206Pb/238U 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.
A comparison of the 206Pb/238U date maps and core to rim 206Pb/238U 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 206Pb/238U 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 206Pb/238U 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 206Pb/238U 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 206Pb/238U dates at their boundaries, while different inclusions in alkali feldspar yield different 206Pb/238U dates at their boundaries. Finally, it is notable that the youngest 206Pb/238U 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 206Pb/238U 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 206Pb/238U 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 206Pb/238U 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 206Pb/238U 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]).

5. 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 (DGB) exceed those for crystal volumes (DM) by factors of 104–108, 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 (MGB) and crystal volumes (MM) in regimes C and B and at the effective diffusivity (Deffective) 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 DGB/DM = 104 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 Deffective is only slightly greater that DM (Figure 8c). However, if DGB/DM = 108, 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 Deffective exceeds DM 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.

6. 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.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/min11040364/s1, Archive with MATLAB scripts that were used here, scripts.zip; archive with examples of HeFTy files and a supplementary figure illustrating our inversion modelling results, inversion.zip.

Author Contributions

Conceptualisation and MATLAB script development, data acquisition, analysis and presentation were carried out by D.V.P. with the advice and input from R.A.S. as a supervisor and project administrator. Funding for this work was acquired first by R.A.S. and later by D.V.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the Swiss National Science Foundation through research grant 200021_160052 awarded to R.S. followed by Early Postdoc.Mobility fellowship P2GEP2_191478 awarded to D.P.

Acknowledgments

We are grateful to Yu. Podladchikov for his useful comments and suggestions during early stages of writing the MATLAB scripts that were used here. We are also grateful to M. Odlum, the authors of three anonymous reviews and the editor whose criticisms helped to improve our manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Dodson, M.H. Closure Temperature in Cooling Geochronological and Petrological Systems. Contrib. Mineral. Petrol. 1973, 40, 259–274. [Google Scholar] [CrossRef]
  2. Dodson, M.H. Closure Profiles in Cooling Systems. Mater. Sci. Forum 1986, 7, 145–154. [Google Scholar] [CrossRef]
  3. Ganguly, J.; Tirone, M. Diffusion Closure Temperature and Age of a Mineral with Arbitrary Extent of Diffusion: Theoretical Formulation and Applications. Earth Planet. Sci. Lett. 1999, 170, 131–140. [Google Scholar] [CrossRef]
  4. Guralnik, B.; Jain, M.; Herman, F.; Paris, R.B.; Harrison, T.M.; Murray, A.S.; Valla, P.G.; Rhodes, E.J. Effective Closure Temperature in Leaky and/or Saturating Thermochronometers. Earth Planet. Sci. Lett. 2013, 384, 209–218. [Google Scholar] [CrossRef]
  5. Ketcham, R.A. Forward and Inverse Modeling of Low-Temperature Thermochronometry Data. Rev. Mineral. Geochem. 2005, 58, 275–314. [Google Scholar] [CrossRef]
  6. Gallagher, K. Transdimensional Inverse Thermal History Modeling for Quantitative Thermochronology. J. Geophys. Res. Solid Earth 2012, 117. [Google Scholar] [CrossRef] [Green Version]
  7. Chamberlain, K.R.; Bowring, S.A. Apatite–Feldspar U–Pb Thermochronometer: A Reliable, Mid-Range (450 °C), Diffusion-Controlled System. Chem. Geol. 2001, 172, 173–200. [Google Scholar] [CrossRef]
  8. Schoene, B.; Bowring, S.A. Determining Accurate Temperature–Time Paths from U–Pb Thermochronology: An Example from the Kaapvaal Craton, Southern Africa. Geochim. Cosmochim. Acta 2007, 71, 165–185. [Google Scholar] [CrossRef]
  9. Cochrane, R.; Spikings, R.A.; Chew, D.; Wotzlaw, J.-F.; Chiaradia, M.; Tyrrell, S.; Schaltegger, U.; Van der Lelij, R. High Temperature (>350 °C) Thermochronology and Mechanisms of Pb Loss in Apatite. Geochim. Cosmochim. Acta 2014, 127, 39–56. [Google Scholar] [CrossRef]
  10. Paul, A.N.; Spikings, R.A.; Ulianov, A.; Ovtcharova, M. High Temperature (>350 °C) Thermal Histories of the Long Lived (>500 Ma) Active Margin of Ecuador and Colombia: Apatite, Titanite and Rutile U-Pb Thermochronology. Geochim. Cosmochim. Acta 2018, 228, 275–300. [Google Scholar] [CrossRef]
  11. Paul, A.N.; Spikings, R.A.; Chew, D.; Daly, J.S. The Effect of Intra-Crystal Uranium Zonation on Apatite U-Pb Thermochronology: A Combined ID-TIMS and LA-MC-ICP-MS Study. Geochim. Cosmochim. Acta 2019, 251, 15–35. [Google Scholar] [CrossRef]
  12. Reiners, P.W.; Ehlers, T.A.; Zeitler, P.K. Past, Present, and Future of Thermochronology. Rev. Mineral. Geochem. 2005, 58, 1–18. [Google Scholar] [CrossRef] [Green Version]
  13. Smye, A.J.; Stockli, D.F. Rutile U–Pb Age Depth Profiling: A Continuous Record of Lithospheric Thermal Evolution. Earth Planet. Sci. Lett. 2014, 408, 171–182. [Google Scholar] [CrossRef]
  14. Smye, A.J.; Marsh, J.H.; Vermeesch, P.; Garber, J.M.; Stockli, D.F. Applications and Limitations of U-Pb Thermochronology to Middle and Lower Crustal Thermal Histories. Chem. Geol. 2018, 494, 1–18. [Google Scholar] [CrossRef] [Green Version]
  15. Villa, I.M. From Nanometer to Megameter: Isotopes, Atomic-Scale Processes, and Continent-Scale Tectonic Models. Lithos 2006, 87, 155–173. [Google Scholar] [CrossRef]
  16. Harrison, T.M.; Heizler, M.T.; McKeegan, K.D.; Schmitt, A.K. In Situ 40K–40Ca ‘Double-plus’ SIMS Dating Resolves Klokken Feldspar 40K–40Ar Paradox. Earth Planet. Sci. Lett. 2010, 299, 426–433. [Google Scholar] [CrossRef]
  17. Villa, I.M.; Hanchar, J.M. K-Feldspar Hygrochronology. Geochim. Cosmochim. Acta 2013, 101, 24–33. [Google Scholar] [CrossRef]
  18. Villa, I.M. Diffusion in Mineral Geochronometers: Present and Absent. Chem. Geol. 2016, 420, 1–10. [Google Scholar] [CrossRef]
  19. Holder, R.M.; Hacker, B.R. Fluid-Driven Resetting of Titanite Following Ultrahigh-Temperature Metamorphism in Southern Madagascar. Chem. Geol. 2019, 504, 38–52. [Google Scholar] [CrossRef]
  20. Popov, D.V.; Spikings, R.A. Diffusion vs. Fluid Alteration in Alkali Feldspar 40Ar/39Ar Thermochronology: Does Cross-Correlation of Log(r/R₀) and Age Spectra Validate Thermal Histories? Chem. Geol. 2020, 539, 119506. [Google Scholar] [CrossRef]
  21. Popov, D.V.; Spikings, R.A.; Scaillet, S.; O’Sullivan, G.; Chew, D.; Badenszki, E.; Daly, J.S.; Razakamanana, T.; Davies, J.H.F.L. Diffusion and Fluid Interaction in Itrongay Pegmatite (Madagascar): Evidence from in Situ 40Ar/39Ar Dating of Gem-Quality Alkali Feldspar and U-Pb Dating of Protogenetic Apatite Inclusions. Chem. Geol. 2020, 556, 119841. [Google Scholar] [CrossRef]
  22. Kirkland, C.L.; Yakymchuk, C.; Szilas, K.; Evans, N.; Hollis, J.; McDonald, B.; Gardiner, N.J. Apatite: A U-Pb Thermochronometer or Geochronometer? Lithos 2018, 318–319, 143–157. [Google Scholar] [CrossRef]
  23. Glorie, S.; Jepson, G.; Konopelko, D.; Mirkamalov, R.; Meeuws, F.; Gilbert, S.; Gillespie, J.; Collins, A.S.; Xiao, W.; Dewaele, S.; et al. Thermochronological and Geochemical Footprints of Post-Orogenic Fluid Alteration Recorded in Apatite: Implications for Mineralisation in the Uzbek Tian Shan. Gondwana Res. 2019, 71, 1–15. [Google Scholar] [CrossRef]
  24. Antoine, C.; Bruand, E.; Guitreau, M.; Devidal, J.-L. Understanding Preservation of Primary Signatures in Apatite by Comparing Matrix and Zircon-Hosted Crystals from the Eoarchean Acasta Gneiss Complex (Canada). Geochem. Geophys. Geosyst. 2020, 21. [Google Scholar] [CrossRef]
  25. Albarède, F. The Thermal History of Leaky Chronometers above Their Closure Temperature. Geophys. Res. Lett. 2003, 30, 15-1–15-4. [Google Scholar] [CrossRef]
  26. Gardés, E.; Montel, J.-M. Opening and Resetting Temperatures in Heating Geochronological Systems. Contrib. Mineral. Petrol. 2009, 158, 185–195. [Google Scholar] [CrossRef]
  27. Arnaud, N.O.; Kelley, S.P. Evidence for Excess Argon during High Pressure Metamorphism in the Dora Maira Massif (Western Alps, Italy), Using an Ultra-Violet Laser Ablation Microprobe 40Ar-39Ar Technique. Contrib. Mineral. Petrol. 1995, 121, 1–11. [Google Scholar] [CrossRef]
  28. Scaillet, S. Excess 40Ar Transport Scale and Mechanism in High-Pressure Phengites: A Case Study from an Eclogitized Metabasite of the Dora-Maira Nappe, Western Alps. Geochim. Cosmochim. Acta 1996, 60, 1075–1090. [Google Scholar] [CrossRef]
  29. Pickles, C.; Kelley, S.; Reddy, S.; Wheeler, J. Determination of High Spatial Resolution Argon Isotope Variations in Metamorphic Biotites. Geochim. Cosmochim. Acta 1997, 61, 3809–3833. [Google Scholar] [CrossRef]
  30. Kelley, S. Excess Argon in K–Ar and Ar–Ar Geochronology. Chem. Geol. 2002, 188, 1–22. [Google Scholar] [CrossRef]
  31. Camacho, A.; Lee, J.K.; Hensen, B.J.; Braun, J. Short-Lived Orogenic Cycles and the Eclogitization of Cold Crust by Spasmodic Hot Fluids. Nature 2005, 435, 1191–1196. [Google Scholar] [CrossRef]
  32. Warren, C.; Sherlock, S.; Kelley, S. Interpreting High-Pressure Phengite 40Ar/39Ar Laserprobe Ages: An Example from Saih Hatat, NE Oman. Contrib. Mineral. Petrol. 2011, 161, 991–1009. [Google Scholar] [CrossRef]
  33. McDonald, C.S.; Warren, C.J.; Mark, D.F.; Halton, A.M.; Kelley, S.P.; Sherlock, S.C. Argon Redistribution during a Metamorphic Cycle: Consequences for Determining Cooling Rates. Chem. Geol. 2016, 443, 182–197. [Google Scholar] [CrossRef] [Green Version]
  34. McDonald, C.S.; Regis, D.; Warren, C.J.; Kelley, S.P.; Sherlock, S.C. Recycling Argon through Metamorphic Reactions: The Record in Symplectites. Lithos 2018, 300, 200–211. [Google Scholar] [CrossRef]
  35. Baxter, E.F. Quantification of the Factors Controlling the Presence of Excess 40Ar or 4He. Earth Planet. Sci. Lett. 2003, 216, 619–634. [Google Scholar] [CrossRef]
  36. Kühn, A.; Glodny, J.; Iden, K.; Austrheim, H. Retention of Precambrian Rb/Sr Phlogopite Ages through Caledonian Eclogite Facies Metamorphism, Bergen Arc Complex, W-Norway. Lithos 2000, 51, 305–330. [Google Scholar] [CrossRef]
  37. Glodny, J.; Kühn, A.; Austrheim, H. akon Diffusion versus Recrystallization Processes in Rb–Sr Geochronology: Isotopic Relics in Eclogite Facies Rocks, Western Gneiss Region, Norway. Geochim. Cosmochim. Acta 2008, 72, 506–525. [Google Scholar] [CrossRef] [Green Version]
  38. Giletti, B.J. Rb and Sr Diffusion in Alkali Feldspars, with Implications for Cooling Histories of Rocks. Geochim. Cosmochim. Acta 1991, 55, 1331–1343. [Google Scholar] [CrossRef]
  39. Jenkin, G.R.; Rogers, G.; Fallick, A.E.; Farrow, C.M. Rb-Sr Closure Temperatures in Bi-Mineralic Rocks: A Mode Effect and Test for Different Diffusion Models. Chem. Geol. 1995, 122, 227–240. [Google Scholar] [CrossRef]
  40. Ault, A.K.; Flowers, R.M.; Mahan, K.H. Quartz Shielding of Sub-10 μm Zircons from Radiation Damage-Enhanced Pb Loss: An Example from a Metamorphosed Mafic Dike, Northwestern Wyoming Craton. Earth Planet. Sci. Lett. 2012, 339–340, 57–66. [Google Scholar] [CrossRef]
  41. Ibanez-Mejia, M.; Bloch, E.M.; Vervoort, J.D. Timescales of Collisional Metamorphism from Sm-Nd, Lu-Hf and U-Pb Thermochronology: A Case from the Proterozoic Putumayo Orogen of Amazonia. Geochim. Cosmochim. Acta 2018, 235, 103–126. [Google Scholar] [CrossRef]
  42. Kohn, M.J.; Penniston-Dorland, S.C.; Ferreira, J.C.S. Implications of Near-Rim Compositional Zoning in Rutile for Geothermometry, Geospeedometry, and Trace Element Equilibration. Contrib. Mineral. Petrol. 2016, 171. [Google Scholar] [CrossRef]
  43. Meesters, A.G.C.A.; Dunai, T.J. Solving the Production–Diffusion Equation for Finite Diffusion Domains of Various Shapes: Part I. Implications for low-temperature (U–Th)/He thermochronology. Chem. Geol. 2002, 186, 333–344. [Google Scholar] [CrossRef]
  44. Meesters, A.G.C.A.; Dunai, T.J. Solving the Production–Diffusion Equation for Finite Diffusion Domains of Various Shapes: Part II. Application to cases with a-ejection and nonhomogeneous distribution of the source. Chem. Geol. 2002, 186, 347–363. [Google Scholar] [CrossRef]
  45. Watson, E.B.; Wanser, K.H.; Farley, K.A. Anisotropic Diffusion in a Finite Cylinder, with Geochemical Applications. Geochim. Cosmochim. Acta 2010, 74, 614–633. [Google Scholar] [CrossRef]
  46. Gautheron, C.; Tassan-Got, L. A Monte Carlo Approach to Diffusion Applied to Noble Gas/Helium Thermochronology. Chem. Geol. 2010, 273, 212–224. [Google Scholar] [CrossRef]
  47. Farley, K.A.; Shuster, D.L.; Ketcham, R.A. U and Th Zonation in Apatite Observed by Laser Ablation ICPMS, and Implications for the (U–Th)/He System. Geochim. Cosmochim. Acta 2011, 75, 4515–4530. [Google Scholar] [CrossRef]
  48. Huber, C.; Cassata, W.S.; Renne, P.R. A Lattice Boltzmann Model for Noble Gas Diffusion in Solids: The Importance of Domain Shape and Diffusive Anisotropy and Implications for Thermochronometry. Geochim. Cosmochim. Acta 2011, 75, 2170–2186. [Google Scholar] [CrossRef]
  49. Fox, M.; McKeon, R.E.; Shuster, D.L. Incorporating 3-D Parent Nuclide Zonation for Apatite 4He/3He Thermochronometry: An Example from the Appalachian Mountains. Geochem. Geophys. Geosyst. 2014, 15, 4217–4229. [Google Scholar] [CrossRef] [Green Version]
  50. Crank, J. The Mathematics of Diffusion, 2nd ed.; Oxford University Press: Oxford, UK, 1975; p. 414. [Google Scholar]
  51. Cherniak, D.J. Diffusion in Accessory Minerals: Zircon, Titanite, Apatite, Monazite and Xenotime. Rev. Mineral. Geochem. 2010, 72, 827–869. [Google Scholar] [CrossRef]
  52. Cherniak, D.J. Cation Diffusion in Feldspars. Rev. Mineral. Geochem. 2010, 72, 691–733. [Google Scholar] [CrossRef]
  53. Cherniak, D.J.; Dimanov, A. Diffusion in Pyroxene, Mica and Amphibole. Rev. Mineral. Geochem. 2010, 72, 641–690. [Google Scholar] [CrossRef]
  54. Cherniak, D.J. Pb Diffusion in Rutile. Contrib. Mineral. Petrol. 2000, 139, 198–207. [Google Scholar] [CrossRef]
  55. Bindeman, I.N.; Davis, A.M.; Drake, M.J. Ion Microprobe Study of Plagioclase-Basalt Partition Experiments at Natural Concentration Levels of Trace Elements. Geochim. Cosmochim. Acta 1998, 62, 1175–1193. [Google Scholar] [CrossRef]
  56. Nielsen, R.L.; Ustunisik, G.; Weinsteiger, A.B.; Tepley, F.J.; Johnston, A.D.; Kent, A.J.R. Trace Element Partitioning between Plagioclase and Melt: An Investigation of the Impact of Experimental and Analytical Procedures. Geochem. Geophys. Geosyst. 2017, 18, 3359–3384. [Google Scholar] [CrossRef]
  57. Prowatke, S.; Klemme, S. Trace Element Partitioning between Apatite and Silicate Melts. Geochim. Cosmochim. Acta 2006, 70, 4513–4527. [Google Scholar] [CrossRef]
  58. Blundy, J. Mineral-Melt Partitioning of Uranium, Thorium and Their Daughters. Rev. Mineral. Geochem. 2003, 52, 59–123. [Google Scholar] [CrossRef] [Green Version]
  59. Fedele, L.; Lustrino, M.; Melluso, L.; Morra, V.; Zanetti, A.; Vannucci, R. Trace-Element Partitioning between Plagioclase, Alkali Feldspar, Ti-Magnetite, Biotite, Apatite, and Evolved Potassic Liquids from Campi Flegrei (Southern Italy). Am. Mineral. 2015, 100, 233–249. [Google Scholar] [CrossRef]
  60. Arzilli, F.; Fabbrizio, A.; Schmidt, M.W.; Petrelli, M.; Maimaiti, M.; Dingwell, D.B.; Paris, E.; Burton, M.; Carroll, M.R. The Effect of Diffusive Re-Equilibration Time on Trace Element Partitioning between Alkali Feldspar and Trachytic Melts. Chem. Geol. 2018, 495, 50–66. [Google Scholar] [CrossRef] [Green Version]
  61. Joesten, R. Grain-Boundary Diffusion Kinetics in Silicate and Oxide Minerals. In Diffusion, Atomic Ordering, and Mass Transport; Ganguly, J., Ed.; Advances in Physical Geochemistry; Springer: New York, NY, USA, 1991; Volume 8, pp. 345–395. ISBN 978-1-4613-9021-3. [Google Scholar]
  62. Dohmen, R.; Milke, R. Diffusion in Polycrystalline Materials: Grain Boundaries, Mathematical Models, and Experimental Data. Rev. Mineral. Geochem. 2010, 72, 921–970. [Google Scholar] [CrossRef]
  63. Kaur, I.; Mishin, Y.; Gust, W. Fundamentals of Grain and Interphase Boundary Diffusion, 3rd ed.; John Wiley: Chichester, UK; New York, NY, USA, 1995; ISBN 978-0-471-93819-4. [Google Scholar]
  64. Mishin, Y.; Herzig, C. Diffusion in Fine-Grained Materials: Theoretical Aspects and Experimental Possibilities. Nanostruct. Mater. 1995, 6, 859–862. [Google Scholar] [CrossRef]
  65. Mishin, Y.; Herzig, C. Grain Boundary Diffusion: Recent Progress and Future Research. Mater. Sci. Eng. A 1999, 260, 55–71. [Google Scholar] [CrossRef]
  66. Hiraga, T.; Kohlstedt, D.L. Equilibrium Interface Segregation in the Diopside–Forsterite System I: Analytical Techniques, Thermodynamics, and Segregation Characteristics. Geochim. Cosmochim. Acta 2007, 71, 1266–1280. [Google Scholar] [CrossRef]
  67. Farver, J.R.; Yund, R.A. Grain Boundary Diffusion of Oxygen, Potassium and Calcium in Natural and Hot-Pressed Feldspar Aggregates. Contrib. Mineral. Petrol. 1995, 118, 340–355. [Google Scholar] [CrossRef]
Minerals 11 00364 g001 550
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.
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.
Minerals 11 00364 g001
Minerals 11 00364 g002 550
Figure 2. 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.
Figure 2. 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.
Minerals 11 00364 g002
Minerals 11 00364 g003 550
Figure 3. A comparison of Tc Inclusion and Tc Dodson for chemically homogeneous spherical apatite inclusions in different host minerals. (a) The dependence of the difference between Tc Inclusion and Tc Dodson on diffusion properties of the host mineral for fixed r, Kd and CR. The green dotted line indicates the transition from open (Tc Inclusion = Tc Dodson) to flux-limited (Tc Inclusion > Tc Dodson) boundary behaviour. Data points for different minerals show the available experimental data on Pb diffusion parameters [51,52,53,54]. Marker fill colours are used to highlight separate determinations from the same mineral. (b) The dependence of the transition from open to flux-limited boundary behaviour on Kd for fixed r and CR. (c) The dependence of the transition from open to flux-limited boundary behaviour on r and CR for a fixed Kd. (d) The dependence of the difference between Tc Inclusion and Tc Dodson on EaHost/EaApatite for different r and CR and fixed Kd and log(D0Host/D0Apatite). Tc Inclusion–Tc Dodson values plateau towards higher EaHost/EaApatite values because Tc Inclusion reaches the modelled crystallisation temperature.
Figure 3. A comparison of Tc Inclusion and Tc Dodson for chemically homogeneous spherical apatite inclusions in different host minerals. (a) The dependence of the difference between Tc Inclusion and Tc Dodson on diffusion properties of the host mineral for fixed r, Kd and CR. The green dotted line indicates the transition from open (Tc Inclusion = Tc Dodson) to flux-limited (Tc Inclusion > Tc Dodson) boundary behaviour. Data points for different minerals show the available experimental data on Pb diffusion parameters [51,52,53,54]. Marker fill colours are used to highlight separate determinations from the same mineral. (b) The dependence of the transition from open to flux-limited boundary behaviour on Kd for fixed r and CR. (c) The dependence of the transition from open to flux-limited boundary behaviour on r and CR for a fixed Kd. (d) The dependence of the difference between Tc Inclusion and Tc Dodson on EaHost/EaApatite for different r and CR and fixed Kd and log(D0Host/D0Apatite). Tc Inclusion–Tc Dodson values plateau towards higher EaHost/EaApatite values because Tc Inclusion reaches the modelled crystallisation temperature.
Minerals 11 00364 g003
Minerals 11 00364 g004 550
Figure 4. A comparison of t-T paths that were obtained using HeFTy [5] for our synthetic bulk-grain 206Pb/238U dates of apatite inclusions in alkali feldspar (in red and orange) with the original thermal histories simulated in spherincl.m (black line). We originally assumed flux-limited boundary behaviour with Kd = 5, while inversion modelling assumed open boundary behaviour. (a) The case of monotonic cooling history. (b) The case of a thermal history with a reheating event. Good-fit reconstructed t-T paths have p > 0.5, while acceptable-fit ones have p > 0.05.
Figure 4. A comparison of t-T paths that were obtained using HeFTy [5] for our synthetic bulk-grain 206Pb/238U dates of apatite inclusions in alkali feldspar (in red and orange) with the original thermal histories simulated in spherincl.m (black line). We originally assumed flux-limited boundary behaviour with Kd = 5, while inversion modelling assumed open boundary behaviour. (a) The case of monotonic cooling history. (b) The case of a thermal history with a reheating event. Good-fit reconstructed t-T paths have p > 0.5, while acceptable-fit ones have p > 0.05.
Minerals 11 00364 g004
Minerals 11 00364 g005 550
Figure 5. A comparison of synthetic bulk-grain 206Pb/238U dates of chemically zoned apatite crystals with either open or flux limited boundaries. (a) Modelled core to rim U concentration profiles. (b,c) Modelled “cooler” and “hotter” thermal histories, respectively. (d) 206Pb/238U dates obtained from apatite crystals with open boundaries that experienced the “cooler” thermal history (the first group). (e,f) 206Pb/238U dates of alkali feldspar-hosted apatite inclusions that followed the “cooler” and “hotter” thermal histories, respectively (the second and third groups, respectively).
Figure 5. A comparison of synthetic bulk-grain 206Pb/238U dates of chemically zoned apatite crystals with either open or flux limited boundaries. (a) Modelled core to rim U concentration profiles. (b,c) Modelled “cooler” and “hotter” thermal histories, respectively. (d) 206Pb/238U dates obtained from apatite crystals with open boundaries that experienced the “cooler” thermal history (the first group). (e,f) 206Pb/238U dates of alkali feldspar-hosted apatite inclusions that followed the “cooler” and “hotter” thermal histories, respectively (the second and third groups, respectively).
Minerals 11 00364 g005
Minerals 11 00364 g006 550
Figure 6. A comparison of bulk-grain and intra-grain 206Pb/238U dates of 12 apatite crystals that have different geometries, U zonation and degree of boundary openness but the same t-T path. (a) The modelled t-T path. (b) U concentration maps for elongated and r-normalised spherical crystals. U concentration maps for S/V-normalised spherical crystals are obtained by multiplying the r and z coordinates in those for r-normalised spherical crystals by 1.4. (c) Assumptions about radiogenic Pb behaviour at the apatite boundaries. (d) 206Pb/238U date maps for chemically homogeneous elongated and r-normalised spherical crystals. (e) Core to rim 206Pb/238U date profiles and bulk-grain 206Pb/238U dates for all of the chemically homogeneous crystals. (f) 206Pb/238U date maps for chemically zoned elongated and r-normalised spherical crystals. (g) Core to rim 206Pb/238U date profiles and bulk-grain 206Pb/238U dates for all of the chemically zoned crystals.
Figure 6. A comparison of bulk-grain and intra-grain 206Pb/238U dates of 12 apatite crystals that have different geometries, U zonation and degree of boundary openness but the same t-T path. (a) The modelled t-T path. (b) U concentration maps for elongated and r-normalised spherical crystals. U concentration maps for S/V-normalised spherical crystals are obtained by multiplying the r and z coordinates in those for r-normalised spherical crystals by 1.4. (c) Assumptions about radiogenic Pb behaviour at the apatite boundaries. (d) 206Pb/238U date maps for chemically homogeneous elongated and r-normalised spherical crystals. (e) Core to rim 206Pb/238U date profiles and bulk-grain 206Pb/238U dates for all of the chemically homogeneous crystals. (f) 206Pb/238U date maps for chemically zoned elongated and r-normalised spherical crystals. (g) Core to rim 206Pb/238U date profiles and bulk-grain 206Pb/238U dates for all of the chemically zoned crystals.
Minerals 11 00364 g006
Minerals 11 00364 g007 550
Figure 7. A comparison of the modelled crystal geometries.
Figure 7. A comparison of the modelled crystal geometries.
Minerals 11 00364 g007
Minerals 11 00364 g008 550
Figure 8. Effects of the abundance and diffusion properties of grain boundaries on the transport of Pb in polycrystalline alkali feldspar. (a) The relationship between time, temperature and the kinetic regime of diffusion in polycrystalline alkali feldspar depending on the assumed diffusivity of Pb along grain boundaries. The blue line shows t’, the brown to orange lines show t’’ for different DGB/DM ratios, the lower and upper black dotted lines show time when Harrison’s condition for regime A becomes satisfied for crystal sizes of 50 μm and 1 mm, respectively (see text for details). The green region indicates the range of DGB/DM ratios for K and Ca in feldspars at a given temperature [52,67]. The grey lines indicate 1 ka, 1 Ma and 1 Ga of isothermal residence. (b) The contribution of grain boundaries to the total absorption of Pb (expressed in a form of MGB/MM ratios) depending on the spacing between grain boundaries for different DGB/DM ratios and kinetic regimes. The lower line for each DGB/DM ratio corresponds to regime C, while the upper line corresponds to regime B at the time of transition to regimes B4 or A. MGB/MM ratios gradually increase from the lower line to the upper line while within regime B. Green shading highlights crystal sizes ranging from 50 μm to 1 mm. The grey lines indicate the fraction of the total absorbed mass that was absorbed by grain boundaries. (c) The relationship between Deffective and the spacing between grain boundaries for different DGB/DM ratios. The green field shows crystal sizes that range from 50 μm to 1 mm.
Figure 8. Effects of the abundance and diffusion properties of grain boundaries on the transport of Pb in polycrystalline alkali feldspar. (a) The relationship between time, temperature and the kinetic regime of diffusion in polycrystalline alkali feldspar depending on the assumed diffusivity of Pb along grain boundaries. The blue line shows t’, the brown to orange lines show t’’ for different DGB/DM ratios, the lower and upper black dotted lines show time when Harrison’s condition for regime A becomes satisfied for crystal sizes of 50 μm and 1 mm, respectively (see text for details). The green region indicates the range of DGB/DM ratios for K and Ca in feldspars at a given temperature [52,67]. The grey lines indicate 1 ka, 1 Ma and 1 Ga of isothermal residence. (b) The contribution of grain boundaries to the total absorption of Pb (expressed in a form of MGB/MM ratios) depending on the spacing between grain boundaries for different DGB/DM ratios and kinetic regimes. The lower line for each DGB/DM ratio corresponds to regime C, while the upper line corresponds to regime B at the time of transition to regimes B4 or A. MGB/MM ratios gradually increase from the lower line to the upper line while within regime B. Green shading highlights crystal sizes ranging from 50 μm to 1 mm. The grey lines indicate the fraction of the total absorbed mass that was absorbed by grain boundaries. (c) The relationship between Deffective and the spacing between grain boundaries for different DGB/DM ratios. The green field shows crystal sizes that range from 50 μm to 1 mm.
Minerals 11 00364 g008
Table
Table 1. 206Pb/238U dates that were used for inversion modelling in HeFTy.
Table 1. 206Pb/238U dates that were used for inversion modelling in HeFTy.
Radius (μm)206Pb/238U Date ± 1σ (Ma) Monotonic Cooling206Pb/238U Date ± 1σ (Ma) History with Reheating
10120.0 ± 1.293.1 ± 0.9
20120.4 ± 1.2114.4 ± 1.1
30120.7 ± 1.2132.3 ± 1.3
50121.0 ± 1.2159.4 ± 1.6
100121.4 ± 1.2196.9 ± 2.0
150121.7 ± 1.2214.4 ± 2.1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Popov, D.V.; Spikings, R.A. Numerical Modelling of Radiogenic Ingrowth and Diffusion of Pb in Apatite Inclusions with Variable Shape and U-Th Zonation. Minerals 2021, 11, 364. https://doi.org/10.3390/min11040364

AMA Style

Popov DV, Spikings RA. Numerical Modelling of Radiogenic Ingrowth and Diffusion of Pb in Apatite Inclusions with Variable Shape and U-Th Zonation. Minerals. 2021; 11(4):364. https://doi.org/10.3390/min11040364

Chicago/Turabian Style

Popov, Daniil V., and Richard A. Spikings. 2021. "Numerical Modelling of Radiogenic Ingrowth and Diffusion of Pb in Apatite Inclusions with Variable Shape and U-Th Zonation" Minerals 11, no. 4: 364. https://doi.org/10.3390/min11040364

APA Style

Popov, D. V., & Spikings, R. A. (2021). Numerical Modelling of Radiogenic Ingrowth and Diffusion of Pb in Apatite Inclusions with Variable Shape and U-Th Zonation. Minerals, 11(4), 364. https://doi.org/10.3390/min11040364

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop