Determining cooling rates from mica 40Ar/39Ar thermochronology data: Effect of cooling path shape

Tectonic models are commonly underpinned by metamorphic cooling rates derived from diffusive‐loss thermochronology data. Such cooling ages are usually linked to temperature via Dodson's closure temperature (TC) formulation, which specifies a 1/time‐shaped cooling path (Contributions to Mineralogy and Petrology, 1973, 40, 259). Geologists, however, commonly discuss cooling rates as a linear temperature/time shape. We present the results of a series of simple finite‐difference diffusion models for Ar diffusion in muscovite and biotite that show that the difference in recorded age between 1/t and linear cooling paths increases significantly with hotter starting temperatures, slower cooling rates and smaller grain sizes. Our results show that it is essential to constrain the cooling path shape in order to make meaningful interpretations of the measured data.


| INTRODUC TI ON
The ratio of parent to daughter radiogenic isotopes has been used for over a hundred years to constrain geological ages and time-scales (e.g. as reviewed by Condon & Schmitz, 2013). Minerals that host the radioactive parent element are commonly referred to as either "geochronometers", which record the timing of their crystallisation or "thermochronometers", which record the timing of cooling through an estimated temperature window at some point after their crystallisation (e.g. as reviewed by Reiners, 2005). The record of different time-temperature pairs in any one rock or tectonic region helps to constrain thermal history and thus provide clues about the mechanism(s) by which the rocks were exhumed to the surface.
Many thermochronometers are based on the premise that some of the daughter isotope concentration is lost via thermally activated diffusion at high temperatures, and that the resulting mineral age can be linked to temperature via the mathematics governing such diffusion. The temperature of a thermochronometer-bearing rock at the time the thermochronometer recorded its apparent (bulk, whole-grain average) cooling age is most commonly estimated using Dodson's closure temperature (T C ) formulation (Dodson, 1973), which, for thermally activated diffusion described by is given by: where D is the diffusion coefficient, D 0 is the diffusion preexponential factor, R is the gas constant, E A is the activation energy, a is the diffusion (or grain) radius, A is a grain-shape-related constant and τ relates the T C to cooling rate: This result of an analytical solution to the diffusion equation has had an enduring legacy due to its mathematical elegance and simplicity of application. However, the Dodson T C formulation is underpinned by several important assumptions and approximations: 1. that thermally activated volume diffusion was the only mechanism by which the daughter isotope was mobilised within the mineral; 2. that the mineral crystallised with no inherited daughter isotope; 3. that a daughter isotope concentration of zero was maintained at the mineral grain boundary throughout cooling; 4. that the starting temperature was high enough for diffusion of the daughter isotope to be efficient, and removal from the grain to be geologically instantaneous, and 5. that the cooling path from the time of crystallisation to the time of closure conformed to a 1/t-shape (where t is time).
These approximations have a major impact on the applicability of the formulation to any particular geological scenario. The further any scenario deviates from these assumptions, the greater the (commonly unquantified and unreported) interpretational uncertainties on the link between age and temperature. A refinement of the T C formulation to consider cases that did not conform to point (4) was proposed by Ganguly & Tirone, 1999, but has not been applied by the thermochronometer community to nearly the same extent that the original Dodson formulation has been. The Dodson closure temperature formulation is most commonly used to constrain cooling rates by linking the T C + time pair to a higher temperature + time pair linearly. However, T C has been derived explicitly for temperature histories that involve cooling proportional to 1/t ( Figure 1) as this creates a linear time dependence in the exponent in exp( −E A /RT ) and allows the analytic integration of the time dependence. To calculate a closure temperature using the Dodson T C formulation and then to use that result to calculate a linear cooling rate is therefore both circular (as also noted by e.g. Ganguly & Tirone, 2009) and ultimately incorrect.
Modern analytical equipment can now provide ever more precise isotope concentration (age) data, at ever-increasing spatial resolution. Furthermore, the diffusion equation can be solved numerically on any standard computer. Here we investigate the effects of linear and 1/t cooling path shapes on the bulk ages and core-rim age profiles of Ar in muscovite and biotite in grains of different radii that have cooled from different temperatures at different rates. The model results show that the ages recorded by muscovite and biotite that have cooled following these different simple end-member paths   differ significantly, especially at higher peak temperatures' increase and smaller grain sizes. The results also allow cooling rates to be determined directly if there is independent evidence for cooling path shape and so long as the time at which cooling started is known.

| THE D IFFARG P_ INVER S E CODE
The finite-difference code DiffArgP_inverse is a modified version of DiffArg (Wheeler, 1996) Further details of the DiffArg_Inverse code are presented in Data S1.

| ME THODS
The bulk (volume-integrated) 40 Ar/ 39 Ar ages of muscovite and biotite of different grain sizes were modelled for a variety of different starting temperatures and linear vs. inverse (1/t) cooling histories.
Muscovite and biotite were modelled with cylindrical geometry and grain radii of 1, 0.5 and 0.25 mm as these are the most typical grain sizes picked for metamorphic 40 Ar/ 39 Ar analyses. The diffusion parameters applied to each mineral are outlined in Table 1.
All minerals were modelled as "crystallising" then instantaneously cooling from starting temperatures of 700, 600, 500 and 450°C at a starting pressure of 1 GPa to represent a variety of metamorphic terranes exhuming from mid-crustal conditions (Tables 1,   2, 3, Table S2 and Table S4). A series of muscovite models was run at a starting pressure of 2 GPa to more closely match conditions found in subduction zones (c.f. ; Table 4, Table S3), and a further series of muscovite models was run with spherical geometry to allow comparison with the cylindrical geometry models (Table   S5 and Figure S6). Linear cooling rates of 5, 10, 25, 50 and 70°C/ Ma were run in order to compare results for typical rates of cooling in different tectonic terranes. 1/t cooling rate models were run for equivalent "time to reach 0°C" as the linear models, in order to compare results for different cooling path shapes. Model pressures were decreased to 0 GPa over the same time interval.
The grain-boundary conditions in all models were modelled as zero daughter element concentration, in order to investigate behaviour in an open system. Model ages were calculated for two-dimensional (cylindrical) diffusion geometry (Hames & Bowring, 1994) and the time integration was performed using the Crank-Nicholson solver, with a recommended time step that is 10 times larger than the value suggested for a stable fully explicit method (Table 1; Wheeler, 1996).  Table S7.

| RE SULTS
The model results are plotted in Figures 2 and 3 (muscovite modelled from pressures of 1 and 2 GPa respectively) and 4 (biotite from 1 GPa). Summary model results for the bulk (volume-averaged) ages are presented in Tables 2, 3 (muscovite at 1 and 2 GPa respectively) and 4 (biotite at 1 GPa). Full results including core-rim model age variations are presented in Table S2 (muscovite 1 GPa) Table S3 (muscovite 2 GPa), Table S4 (biotite) and Table S5 (muscovite 1 GPa with spherical geometry).
The graphs all show similar trends: 1. Faster cooling results in a smaller difference in time between the timing of maximum temperature attainment (cooling initiation) and the recorded cooling age (Δt).

3.
Smaller grain sizes result in larger Δt.

4.
Smaller Δt values are recorded for the 1/t models than for the linear models.
Results (1) and (3) provide consistency checks to show that the models are behaving as expected. Result (2) similarly matches the predictions of the modified formulation of Ganguly & Tirone, 1999. Result (4) clearly shows the importance of the cooling path shape on the resulting thermochronometer age-this will be discussed further below.
(e) Figure 2 shows that very little diffusive loss is expected in white mica grains that cool from relatively low peak temperatures of 450°C. The 40 Ar/ 39 Ar age of a 0.25 mm radius white mica grain cooling linearly at a rate of 5°C/Ma from 450°C and 1 GPa would be expected to be 2.4 Ma younger than the peak temperature age, whereas one cooling from 700°C would be expected to yield an age that is 52 Ma younger (Table 2). Similarly, the 40 Ar/ 39 Ar age of a 0.25 mm radius white mica grain cooling linearly at a rate of 5°C/ Ma from 600°C and 2 GPa would be expected to be ~25 Ma younger than the peak temperature age (Table 4). Similar-sized grains cooling to 0°C over the same time interval but following a 1/t path from 450 or 700°C at 1 GPa would only yield ages that were 0.6 or 16 Ma younger than the peak temperature age. A 1 mm radius grain cooling from 450°C, however, would be expected to record an age within uncertainty of the timing of peak metamorphism.
Models run using spherical diffusion geometry yield slightly younger ages (Δt of 54 Ma rather than 52 Ma for a 0.25 mm radius grain cooling from 700°C at 5°C/Ma, for example; Table S5; Figure S6). Figure 4 shows that biotite should yield significantly younger ages than muscovite for grains of the same radius, cooling from the same starting temperature and following the same cooling path. For example, a 0.25 mm radius grain cooling at 5°C/Ma from 450°C would be expected to be 26 Ma younger than the age of peak temperature metamorphism, whereas one cooling from 700°C at the same rate would be expected to yield an age that was 76 Ma younger (Table 4).

| D ISCUSS I ON
The results clearly show that the shape of the cooling path makes an increasingly important contribution to the recorded thermochronometer age as grain sizes and cooling rates decrease and peak temperatures increase. The uncertainty inherent in using the Dodson T C formulation to estimate (linear) cooling rates therefore also magnifies accordingly.
As predicted by examination of the Arrhenius rate [1], the model results are more sensitive to systematic uncertainties in the experimentally determined activation energy (E A ) than in the preexponential factor (D 0 ) for each mineral (Figures 2-4 and Table S7).
These figures show that uncertainties in the diffusion parameters have a significant, but systematic, effect on the recorded thermochronological ages. These uncertainties apply equally to both cooling history shapes discussed here.
The most recent diffusion parameters for muscovite (Harrison et al., 2009) were calculated for isotropic three-dimensional (spherical) diffusion geometry. Forster & Lister, 2017 suggested that modelling muscovite as a cylinder but using diffusion parameters calculated for spherical geometry is inappropriate. However, the overall difference in the diffusion coefficients obtained for spherical and cylindrical diffusion models is a factor of 2 in D 0 (Forster & Lister, 2017), which corresponds to an uncertainty in the activation energy of < 0.6 kcal/mol at 400 K. This is well below the

| APPLYING MODEL RE SULTS TO NATUR AL SYS TEMS
The results presented here can be used to constrain the cooling rates of natural systems if the following pieces of information are known or can be estimated:

A petrologically-based interpretation of the temperature at
which the dated grain(s) grew, and the portion of the metamorphic path along which the grain(s) grew (e.g. prograde peak or retrograde). This will inform and constrain the extent of diffusive opportunity that the grain could have experienced. For example, a grain growing during the prograde history will have longer residence at high temperatures, therefore allowing it more opportunity to lose argon.

2.
The peak temperature experienced by the grain(s), required for the ultimate determination of a cooling rate.

3.
The time at which the grain reached its peak temperature (constrained or estimated by independent geochronometers), required for the ultimate determination of a cooling rate. This is further discussed below.

4.
The thermochronometric ages of the grains of interest; different data collection methods are further discussed below.

5.
The grain size(s) of the dated grains.

(d) (h)
Note that only very simple cooling path shapes have been modelled here. Steady progress is being made in the development of modelling tools that can suggest a "best-fit" cooling path to U-Th-He, fission track and U-Pb rutile data, but currently none of these tools explicitly incorporate 40 Ar/ 39 Ar data: e.g. HeFTy (Ketcham, 2005), QTQt (Gallagher, 2012), UpBeat (Smye, Marsh, Vermeesch, Garber, & Stockli, 2018).

| Determining the timing of cooling initiation
Direct determination of a cooling rate (and cooling rate shape) from thermo-and geochronological data requires that at least two, and possibly three, T-t pairs are known. Timing of peak T in metamorphic rocks is commonly constrained by U-Pb ages of zircon, monazite, garnet, allanite and/or rutile, with secondary (higher-temperature cooling) T-t pairs provided by U-Pb rutile and/or titanite data. There are, of course, multiple uncertainties inherent in linking these ages to peak temperature because all of these minerals may crystallise at different stages of the metamorphic PT path. Careful petrochronological investigation is required to confirm that the ages yielded by any of these minerals relate to the timing of attainment of peak temperatures or higher-than-argon-closure cooling (e.g. Kohn, Engi, & Lanari, 2017). 40 Ar/ 39 Ar mica data can currently be collected in many different ways: by multiple-or single-grain step heating experiments (e.g. Turner, 1970), by single-grain fusion methods (e.g. Fleck & Carr, 1990) or by laser ablation (e.g. Kelley, Arnaud, & Turner, 1994). All methods have their advantages and disadvantages in terms of volume of material analysed, analytical precision and petrographic (location) control on age.

ME THODS
The model data presented here are compatible for assessment against the bulk (volume-averaged) ages-i.e. equivalent to singlegrain fusion 40 Ar/ 39 Ar data. We caution against using multiple-or single-grain step heating 40 Ar/ 39 Ar ages to compare against model results. Plateau ages imply no core-rim variation in Ar distribution (and thus an interpretation of rapid cooling), but a plateau result does not in itself guarantee that the calculated age is geologically meaningful, especially in high pressure metamorphic rocks (e.g. Sherlock & Arnaud, 1999). Non-plateau spectra can be produced by a variety of factors that complicate linking spectrum shapes to within-grain Ar distribution. Single-grain fusion populations can help provide an assessment of how homogeneous Ar is distributed across mica grains within individual samples (e.g. Uunk, Brouwer, ter Voorde, & Wijbrans, 2018).
In-situ, high-spatial precision 40 Ar/ 39 Ar data such as collected by laser ablation methods, and collected in grains large enough and cooled slowly enough from a high enough temperature to be able to detect such changes, can also be assessed against the core-rim model age predictions for simple linear and 1/t cooling histories presented in Tables S2, S3, S4.

| Comparing analytical data to model results
The time difference (Δt) between the timing of the thermal peak (or to be absolutely correct, the timing of cooling initiation) and the age recorded by the thermochronometer (Figures 2-4)  is not enough information to determine whether (a) the system was diffusively open (a fundamental requirement of any diffusive-based interpretative link between age, temperature and cooling rate is that there is effectively an infinite sink for the daughter element diffusing out of the mineral grain) and/or (b) whether the cooling path was overall linear or some other shape. Both of these can be resolved following a match between data and model predictions. At rapid cooling rates, the difference between the cooling ages predicted by a linear temperature decrease and a 1/t-shaped path would be indistinguishable within the typical uncertainties in analytical results and in the experimental diffusion parameters. At cooling rates < 10°C/ Ma, differences in the shapes of the cooling paths start to become important for distinguishing between exhumation mechanisms.
Small values of Δt e.g. <1 Ma are currently challenging to resolve analytically. The mica 40 Ar/ 39 Ar models for low starting temperatures confirm previous suggestions that rapidly cooled rocks that reached low peak temperatures (such as in subduction zones) will not yield ages that allow cooling rates to be determined.

| Other factors affecting daughter element distribution
Inheritance or loss of daughter product during recrystallisation and deformation during cooling can affect daughter element concentrations much more than diffusion (Allaz, Engi, Berger, & Villa, 2011;Villa, 1998;Villa, Bucher, Bousquet, Kleinhanns, & Schmid, 2014). It is also obvious that recrystallisation during exhumation means that the temperature that that particular grain cooled from may be lower than the peak temperature. In cases where thermochronometer minerals show signs of secondary recrystallisation or other chemical modification, the model results are almost certainly not applicable, and a link between temperature and age may be more difficult to constrain. The diffusion models are only applicable to rocks in which an open system can be assumed, and where both the timing and pressure-temperature conditions of the last episode of mineral crystallisation are known or can be estimated.
If the results presented here are used to estimate cooling rates or constrain cooling path shapes, each practitioner will need to estimate the geological uncertainty for their particular study, noting that this is almost certainly the largest current overall source of error in their interpretation. Our results are based on the assumption that cooling starts directly after the model grain has crystallised at peak temperatures. In reality, the minerals of interest may have grown along the prograde path and/or have resided at peak temperatures for a geologically significant period of time before cooling started. If temperatures were low enough for diffusion to be inefficient, some of that pre-cooling history may be recorded in the thermochronometer minerals. Thermochronologists should model the effect of pre-peak thermal history for their particular geological location to convince themselves whether or not the thermochronometer minerals in their study area may record this.

| CON CLUS IONS
The rates and time-scales over which rocks are buried, transformed, deformed and exhumed help constrain the tectonic mechanisms that act on them. 40 Ar/ 39 Ar data from mica have long been used to link time to temperature and thus constrain cooling rates. The Dodson closure temperature formulation (Dodson, 1973) provides an elegant analytical solution to the diffusion equation but its application for determining cooling rates is commonly based on assumptions that are a poor match to geological reality. Our results of a series of diffusion models that quantify the differences in age expected from a simple linear and 1/t-shaped cooling histories show that the cooling path shape exerts considerable influence on the resulting age at hotter starting temperatures, slower cooling rates and smaller grain sizes. If the cooling path shape and timing of cooling initiation are known, then our results also provide a simple way of estimating cooling rates and cooling rate shapes from the difference between the timing of cooling initiation at maximum temperature and the yielded thermochronometer age. Future incorporation of 40 Ar/ 39 Ar diffusion systematics into forward modelling packages that also consider other thermochronometer data provides the best future solution for constraining cooling rates and thus exhumation mechanisms. More precise diffusion data and more constraints on how real geological systems evolve over time will provide further benefit.

ACK N OWLED G EM ENTS
CSM was funded on an Open University Charter Studentship.
CJW acknowledges prior financial support from a NERC Advanced Fellowship (NE/H016279/1) and a NERC small grant (NE/ J013072/1). The authors thank Simon Kelley for many fruitful discussions about linking age to stage and open system behaviour, and Leo Ingvorsen for helping to run some of the early models.
CSM thanks Sarah Sherlock and Alison Halton for PhD supervision and guidance.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section at the end of the article.
Data S1. Instructions for operating DiffargP_inverse. Table S2. Full muscovite results for 1 GPa models. Table S3. Full muscovite results for 2 GPa models.   for the equivalent "time to surface" which is plotted on the x-axis.
The equivalent linear rate is plotted underneath the "time to surface" value on the linear model plots. The y-axis plots the difference between the time at which cooling starts and the recorded 40 Ar/ 39 Ar age: if this is, the grain size and the starting temperature are known for the analysed samples, then the cooling rate can be read off the graph directly. Note the differences in the y-axis scale between the linear and 1/t results.