An optimization tool for identifying multiple-diffusion domain model parameters

The multiple-diffusion domain (MDD) model empirically describes the diffusive behavior of noble gases in some terrestrial materials and has been commonly used to interpret ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ stepwise degassing observations in K-feldspar. When applied in this manner, the MDD model can be used to test crustal exhumation scenarios by identifying the permissible thermal paths a rock sample could have undergone over geologic time, assuming the diffusive properties of Ar within the mineral are accurately understood. More generally, the MDD model provides a framework for quantifying the temperature-dependent diffusivity of noble gases in minerals. However, constraining MDD parameters that successfully predict the results of step-heating diffusion experiments is a complex task, and the assumptions made by existing numerical methods used to quantify model parameters can bias the absolute temperatures permitted by thermal modeling. For example, the most commonly used method assumes that no domains lose more than 60 % of their gas during early heating steps (Lovera et al., 1997). This assumption is unverifiable, and we show that the Lovera et al. (1997) procedure may bias predicted temperatures towards lower values when it is violated. To address this potential bias and to provide greater accessibility to the MDD model, we present a new open-source method for constraining MDD parameters from stepwise degassing experimental results, called the “MDD Tool Kit” (https://github.com/dgorin1/mddtoolkit, last access: 11 October 2024). This software optimizes all MDD parameters simultaneously and removes any need for user-defined E_a or regression fitting choices used by other tools. In doing so, this new method eliminates assumptions about the domain size distribution. To test the validity of our thermal predictions, we then use the MDD Tool Kit (https://github.com/dgorin1/mddtoolkit) to interpret ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ results from the Grayback Fault, AZ, USA. Although the resulting thermal histories are consistently ∼ 60–75 °C higher than those found in previous studies, they agree with independent observations from apatite fission track, zircon fission track, and $(\mathrm{U}-\mathrm{Th})/\mathrm{He}$.

1 Introduction

${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ thermochronology is a valuable tool for studying Earth's crustal exhumation because it constrains a mineral's continuous thermal history through geologic time (McDougall and Harrison, 1999). While ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ geochronology was initially developed to quantify crystallization timing of rapidly cooled igneous rocks (Turner, 1968), its application to thermochronology did not begin in earnest until Dodson's (1973) quantitative description of closure theory, which allowed one to determine when a mineral cooled below a mineral-specific “closure temperature” for radiogenic ⁴⁰Ar retention (Reiners, 2005). During this early period of thermochronology research, minerals and isotopic systems with disparate closure temperatures were combined to constrain a given rock's permissible thermal histories (e.g., Wagner, 1977). While this approach yielded valuable insight into Earth's surface processes (McDougall and Harrison, 1999), it requires multiple mineral phases – which are not always present – to coexist in a sample. Further, the thermal histories inferred from this technique are discontinuous, with one age–temperature measurement for each mineral phase, thus limiting the ability to quantify cooling rates at specific times (McDougall and Harrison, 1999) and potentially introducing biases. Additionally, this approach does not fully accommodate the complexity of the diffusive properties observed in some minerals (Zeitler, 1987; McDougall and Harrison, 1999).

Since the late 1960s, it had been understood that complex retention properties pose challenges for interpreting the thermal histories of some samples. More specifically, the results of some stepwise degassing experiments were not consistent with volume diffusion from a single domain in some minerals, requiring a new interpretative framework (Zeitler, 1987; Lovera et al., 1989). Progress came first with Zeitler's (1987) observations in K-feldspar that the anomalous behavior seen in ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age spectra and associated Arrhenius plots could be explained by outward diffusion of ³⁹Ar from several “domains” (distinct, non-interacting geometries in which volume diffusion takes place) of varying sizes simultaneously. A formalism for this model, called the multiple-diffusion domain (MDD) model, was then described by Lovera et al. (1989). This model enabled more nuanced interpretations of sample-specific diffusion kinetics – and thereby thermal histories – of minerals.

While this model has received criticism (e.g., Villa, 1994; Parsons et al., 1999; Popov and Spikings, 2020; Popov et al., 2020a, b; Spikings and Popov, 2021), it has largely been adopted by the thermochronology community since the 2000s due to its ability to constrain time–temperature (t–T) histories that are consistent with both an observed ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age spectrum and calculated diffusion kinetics of a given sample (Reiners, 2005; Harrison and Lovera, 2014). However, we show that the method used by Lovera et al. (1997) to tune the MDD model parameters to match the results of a step-heating diffusion experiment potentially introduces a systematic bias towards colder temperatures through geologic time. This was recognized, but not addressed, at the time of their publication, likely due to insufficient computing power to resolve the issue (Lovera et al., 1997).

2 The MDD model

The MDD model is best understood through analogy to volume diffusion through a single domain. In this simple case, the temperature-dependent diffusivity is described by

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{D\left(T\right)}{a^{\mathrm{2}}}}={\displaystyle \frac{D_{\mathrm{0}}}{a^{\mathrm{2}}}}\cdot e^{\frac{-E_{\mathrm{a}}}{RT}},\end{array}$$

where T is the absolute temperature (K), D is diffusivity at T (cm² s⁻¹), D₀ is diffusivity at infinite T, E_a is activation energy (kJ mol⁻¹), R is the gas constant (kJ mol⁻¹ K⁻¹), and a is the radius of the diffusion domain (cm) (see Table 1 for variable definitions).

Lovera et al. (1997)Lovera et al. (1997)Lovera et al. (1997)(Lovera et al., 1997)

Table 1 Table of variables.

This relationship is typically determined empirically for noble gas diffusion in minerals through a stepwise degassing diffusion experiment. In such experiments, the diffusant is first produced in situ by a proton or neutron irradiation to ensure a homogeneous initial distribution. The sample is then placed under static vacuum where it is repeatedly heated to a known temperature, and the quantity of the diffusant is measured (Fechtig and Kalbitzer, 1966; McDougall and Harrison, 1999). By assuming a geometry for the diffusion domain, the fractional loss at each step can be used to calculate a corresponding $\frac{D}{a^{\mathrm{2}}}$ with the equations of Fechtig and Kalbitzer (1966) or Crank (1975). Ginster and Reiners (2018) summarized and propagated measurement uncertainties through these equations, and we use their forms in this work (Table 2).

Table 2 Equations used by the MDD Tool Kit (https://github.com/dgorin1/mddtoolkit) to calculate f, the fractional loss from each domain during each heating step in a stepwise degassing experiment (Crank, 1975; Ginster and Reiners, 2018).

When consistent with volume diffusion through a single diffusion domain, a step-heating experiment will produce a linear Arrhenius relationship between calculated values of $\log \left(\frac{D}{a^{\mathrm{2}}}\right)$ and $\frac{\mathrm{1}}{T}$. In these cases, a linear regression can then be fit to the results to determine E_a and $\frac{D_{\mathrm{0}}}{a_{\mathrm{0}}^{\mathrm{2}}}$ (Lovera et al., 1997) for a sample. Several mineral–diffusant pairs such as ³He in apatite and ³He in olivine exhibit such behavior (Shuster et al., 2004).

However, not all minerals exhibit diffusive behavior consistent with volume diffusion from a single diffusion domain. The diffusive behavior of some minerals appears to be more consistent with diffusion from several non-interacting domains of varying sizes, within a given mineral, diffusing simultaneously. This behavior was later identified and roughly quantified by Gillespie et al. (1982) and Zeitler (1987) in K-feldspar, orthoclase, and microcline and then formalized by Lovera et al. (1989).

In contrast to the single-diffusion domain model, the MDD model can be imagined as a series of non-interacting diffusion domains of varying diffusive length scale (e.g., different radii, a) that all diffuse simultaneously. The choice of geometry for the domains is somewhat arbitrary, as it remains unclear what these diffusion domains physically represent in a mineral sample (Harrison and Lovera, 2014; Parsons et al., 1999). Some authors use a plane sheet, while others prefer a spherical or cylindrical geometry. Regardless of this choice, it is assumed that each diffusion domain has the same geometry, E_a, and D₀. Therefore, a diffusion domain can be described completely by three parameters: (i) E_a (common to all domains), (ii) its pre-exponential term $(\frac{D_{\mathrm{0}}}{a^{\mathrm{2}}}{)}_i$, and (iii) the proportion of the total gas it contains (ϕ). Thus, a multiple-diffusion domain model with n domains has 2n−1 free parameters because ${\sum }_{i=\mathrm{1}}^n{\mathit{\varphi }}_n=\mathrm{1}$.

An inherent challenge in applying this model is that each of these 2n−1 parameters need to be optimized to accurately predict the results of a given diffusion experiment. Because authors using the MDD model fit as many as 10 domains, this optimization is regularly performed on 19-dimensional vectors or larger, making the exercise nontrivial.

2.1 Critiques of the MDD model

Validating the MDD model's assumptions is beyond this paper's scope; here we briefly review the assumptions and their prior critiques. The foundational assumption of the MDD model is that the transport of ⁴⁰Ar within minerals over geologic time occurs primarily by volume diffusion. This assumption predicts that low ⁴⁰Ar concentrations should exist near the outer edge of a crystal (i.e., for a boundary condition of nearly zero ⁴⁰Ar concentration external to the crystal) and that higher concentrations should exist towards the mineral interior, translating to apparently younger and older ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages in those locations, respectively. Studies have documented such spatial correlations (Flude et al., 2014). However, others have shown that polyphasic samples can violate this expectation (Popov and Spikings, 2020); such minerals are not suitable candidates for MDD modeling.

It is also assumed that diffusive behavior observed during laboratory step-heating experiments is the same as under natural conditions over geologic time. The MDD model proposes that the diffusive behavior of Ar within some minerals can be described by numerous non-interacting, infinite sheets simultaneously diffusing within the same mineral. As evidence, Lovera et al. (2002) argued that a correlation between the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age spectra and log($r/r_{\mathrm{0}}$) plots (see Lovera et al., 2002, for a description of log($r/r_{\mathrm{0}}$) plots) validates, or is at least consistent with, these assumptions. They find, however, that only ∼ 40 % of K-feldspars demonstrate sufficient correlation for MDD modeling and suggest that the remaining samples have been affected by recrystallization or other mineral inclusions (Lovera et al., 1997). Other authors have gone further, asserting that recrystallization within minerals is nearly always responsible for this behavior (Popov and Spikings, 2020; Popov et al., 2020a, b; Spikings and Popov, 2021). While complex mineralogy is likely responsible in some cases, detailed petrologic examinations and a strong correlation between an ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age spectrum and a log($r/r_{\mathrm{0}}$) plot can increase confidence that the MDD model framework is applicable to an individual sample. Indeed, such petrologic investigations are commonplace in MDD studies of K-feldspars (e.g., Wong et al., 2023).

Others have questioned the primary assumption that Ar diffusion under laboratory conditions functions the same way as natural diffusion over geologic time by suggesting that mineral structures may be altered during diffusion experiments. Popov et al. (2020a) observed the development of cracks within crystals during some experiments. However, Lovera et al. (1993) performed a double-irradiation experiment where they measured the diffusion kinetics of a grain up to ∼ 850 °C, re-irradiated, it and then successfully reproduced the diffusion kinetics of the first experiment, indicating that cracking or other annealing effects were not a primary cause of MDD behavior in that sample. The Lovera et al. (1993) study also indicated that the diffusion kinetics of the mineral were not altered by heating during the diffusion experiment.

Some further MDD model assumptions are made which require knowledge of the physical structure of the domains. For example, it is assumed that Ar is instantaneously removed from the mineral after reaching a domain boundary, that diffusion in each domain is isotropic, that the domains are all formed at the same time, and that the parent isotope, ⁴⁰K, is uniformly distributed within each domain (Parsons et al., 1999; Lovera et al., 1989). Because it is unclear what the diffusion domains physically represent within a mineral, these assumptions are currently unverifiable. Although it is unsatisfying that a physical representation of the diffusion domains is not clearly identifiable within K-feldspar and other minerals, the MDD model need not be rejected as a tool for thermochronology. Indeed, ample evidence suggests that the MDD model reliably predicts thermal histories supported by independent thermochronological data (e.g., Warnock and Zeitler, 1998; Reiners and Farley, 1999; Axen et al., 2000; Spell et al., 2000; Kirby et al., 2002; Reiners et al., 2004; Shirvell et al., 2009). Though empirical, the MDD model remains a useful tool for constraining such histories through geologic time.

3 Existing optimization methods

The Lovera et al. (1997) method for identifying MDD model parameters is the only published algorithm for this purpose that we are aware of and has been used in many studies (e.g., Harrison et al., 1995; Quidelleur et al., 1997; Grove et al., 2003; Weirich et al., 2012; Wong et al., 2023). The routine begins by defining a reference domain which is used to calculate an E_a by fitting an uncertainty-weighted linear regression to the ${\log }_{\mathrm{10}}\left(\frac{D}{a^{\mathrm{2}}}\right)$ values resulting from a subset of the low-temperature heating steps in an Arrhenius plot corrected for excess Ar released from fluid inclusions (Harrison et al., 1994); the resultant E_a is assumed to be applicable to all domains. This is done because – even in samples exhibiting MDD-like behavior – these diffusivities tend to approximate a line (Fig. 1a).

Figure 1 (a) Arrhenius plot resulting from synthetic Experiment A. Gray circles show the expected $\frac{D}{a^{\mathrm{2}}}$ values calculated from the pre-determined diffusion kinetics (Table 4), and transparent gray circles show the heating steps excluded from the fitting exercise. Black circles represent the values predicted by the MDD Tool Kit (https://github.com/dgorin1/mddtoolkit) method. (b) Same as (a) except showing the expected $\frac{D}{a^{\mathrm{2}}}$ values calculated from the pre-determined diffusion kinetics for Experiment B. (c) Value of q⋅N as the number of points included in the uncertainty-weighted linear regression increases for Experiment A. The Lovera et al. (1997) algorithm for determining E_a maximizes the value of this function to determine the appropriate number of points to include. The black dot indicates this value. (d) The predicted E_a as one increases the number of points included in the unweighted linear regression. The black point indicates the E_a selected by the Lovera et al. (1997) algorithm for Experiment A. In this case, the algorithm selects the correct E_a. (e) Same as panel (c) but showing the results for synthetic Experiment B. (f) Same as panel (d) but for synthetic Experiment B. In this case, the Lovera et al. (1997) algorithm underestimates the correct E_a.

To determine the number of points to include in the uncertainty-weighted linear regression, the lowest-temperature heating steps are sequentially added to the regression. After adding each step, a chi-square misfit static is calculated and is then used to calculate a goodness-of-fit probability (q) (Lovera et al., 1997). These relationships are described as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\mathit{\chi }}_{\mathrm{lovera}}^{\mathrm{2}}=\sum _{i=\mathrm{1}}^N{\left({\displaystyle \frac{{\log }_{\mathrm{10}}(\frac{D}{a^{\mathrm{2}}}{)}_i-{\log }_{\mathrm{10}}(\frac{\widehat{D}}{a^{\mathrm{2}}}{)}_i}{{\mathit{\sigma }}_i}}\right)}^{\mathrm{2}},\text{(3)}& {\displaystyle }q=\mathrm{gammq}\left({\displaystyle \frac{N-\mathrm{2}}{\mathrm{2}}},{\displaystyle \frac{{\mathit{\chi }}^{\mathrm{2}}}{\mathrm{2}}}\right),\end{array}$$

where gammq is the incomplete gamma function, $(\frac{D}{a^{\mathrm{2}}}{)}_i$ is the observed pre-exponential term, $(\frac{\widehat{D}}{a^{\mathrm{2}}}{)}_i$ is the modeled pre-exponential term, σ is the uncertainty in the observed pre-exponential term, and N is the number of steps included (Lovera et al., 1997; Press, 2007). This calculation is repeated until the value of q⋅N is maximized. Once the E_a is determined from the above process, this value is held fixed while a Levenberg–Marquardt method (Press, 2007) is used to adjust the gas fraction (ϕ) and ${\log }_{\mathrm{10}}\left(\frac{D_{\mathrm{0}}}{a^{\mathrm{2}}}\right)$ for each domain to minimize the above chi-square quantity. The reader is directed to Appendix B of Lovera et al. (1997) and the referenced sections of Press (2007) for additional information.

While this routine is generally robust and capable of accurately quantifying E_a in synthetic data experiments where the kinetics are defined, it fails to do so when the smallest domains lose greater than 60 % of their gas in early heating steps (Fig. 1). Thus, use of this routine makes the implicit assumption that no domains diffuse most of their gas during the initial heating steps. When this assumption is violated, the regression will overestimate E_a and ${\log }_{\mathrm{10}}\left(\frac{D_{\mathrm{0}}}{a_{\mathrm{0}}^{\mathrm{2}}}\right)$ (Fig. 1). Because there is currently no way to know the domain size distribution a priori, it is not possible to know whether this assumption is valid for any given sample (Lovera et al., 1997).

Furthermore, Lovera et al. (1997) found evidence that the variability of predicted E_a values from different aliquots of the same mineral decreased when they were able to use a higher percentage of the total gas released in their linear regressions. Simply put, the more gas included in the calculation, the less variation in predicted E_a. This observation suggests that a routine which maximizes the fraction of ³⁹Ar used in the fitting exercise might lead to more precise estimates of diffusion kinetics.

4 The MDD Tool Kit approach

With the increased computational power since the publication of the original routine for optimizing these models, we apply SciPy's implementation of the differential evolution algorithm to fit all the diffusion kinetics parameters simultaneously (Storn and Price, 1997; Lampinen, 2002; Qiang and Mitchell, 2014; Virtanen et al., 2020). This algorithm is a population-based genetic algorithm for optimizing over continuous search spaces in high dimensions with nonlinear and non-differentiable misfit functions (Storn and Price, 1997).

4.1 Differential evolution

Our implementation of SciPy's differential evolution is a method for iteratively improving a randomly selected “population” of guesses until the best-fitting diffusion kinetics are found (Fig. 2). The algorithm improves these guesses by combining elements of the most successful vectors with those remaining. In this case, success refers to the value of a misfit statistic calculated between the guess's forward-modeled gas releases and the observed results. Through this process, each successive generation achieves a misfit equal to or lower than the previous one.

Figure 2 (a) Results of an MDD Tool Kit (https://github.com/dgorin1/mddtoolkit) optimization applied to Wong et al. (2023) sample GR-27 demonstrating the iterative nature of the improvement using the %_frac misfit statistic shown in Arrhenius space. (b) Same as (a), but shown in the space where the optimization is performed – the fractional release at each heating step. Red circles show the optimization halted after 8 iterations (%_frac=15.06), blue circles show the same optimization halted after 28 iterations (%_frac=8.33), and black circles show the same optimization results after complete convergence after 18 419 iterations (%_frac=1.82). The number of iterations is independently determined by the differential evolution algorithm using the convergence criteria outlined by Virtanen et al. (2020). The gray circles show the observed results.

We apply two such misfit statistics to the data analyzed in this study: one error-weighted, χ², and another where all points are weighted equally (%_frac). Our χ² misfit accounts for the uncertainty in the ³⁹Ar measurements by performing its calculations in units of moles so that the measurement uncertainty can readily be included. This is accomplished by multiplying the forward-modeled gas fractions (F_i) by the total number of moles released during the experiment. Because this value is not directly measured, we add it as a parameter to our model (${\widehat{M}}_{\mathrm{tot}}$) and allow our optimization algorithm to solve for it when using this misfit statistic. We allow the model to choose any ${\widehat{M}}_{\mathrm{tot}}$ value within 3σ of the observed value. Our χ² misfit statistic is thereby defined as follows:

$$\begin{array}{}\text{(4)}& {\mathit{\chi }}^{\mathrm{2}}=\sum _{i=\mathrm{1}}^N{\left({\displaystyle \frac{M_i-\widehat{F_i}{\widehat{M}}_{\mathrm{tot}}}{{\mathit{\sigma }}_i}}\right)}^{\mathrm{2}}.\end{array}$$

By contrast, %_frac weights all heating steps equally, regardless of their associated measurement uncertainty. This misfit reports a percent difference between the measured and modeled gas fraction released at each step and is defined as follows:

$$\begin{array}{}\text{(5)}& {\mathit{\%}}_{\mathrm{Frac}}=\sum _{i=\mathrm{1}}^N{\displaystyle \frac{F_i-\widehat{F_i}}{F_i}},\end{array}$$

where F_i is the measured gas fraction at a given heating step, and ${\widehat{F}}_i$ is the modeled gas fraction at a given heating step.

Regardless of which misfit statistic is used, differential evolution optimizes the parameters in the same way and begins by generating its initial population. Given a vector of parameters to be tuned for an n-domain model, $\mathit{X}({\widehat{M}}_{\mathrm{tot}},E_{\mathrm{a}},{\frac{D_{\mathrm{0}}}{a^{\mathrm{2}}}}_{(i,i+\mathrm{1},\mathrm{\dots },n)},{\mathit{\varphi }}_{(i,i+\mathrm{1},\mathrm{\dots },n-\mathrm{1})})$, an initial population of candidate vectors is quasi-randomly generated using the Latin hypercube method (Iman et al., 1981). These vectors ideally capture the full range of the sample space. To avoid user bias, we prescribe search ranges much larger than we imagine to be realistic for each variable based on prior work (Table 3; Lovera et al., 1997).

Table 3 Parameter search ranges used in this publication. These values can be adjusted in the open-access code.

Table 4 Setup and results of synthetic diffusion experiment optimizations. Experiments A and B are the same except that a domain with $\ln (\frac{D_{\mathrm{0}}}{a^{\mathrm{2}}}{)}_{\mathrm{1}}$ has been added to Experiment B. To redistribute the gas, 1 % was removed from all other domains such that ϕ₁=0.01.

From here, the improvement process is iterative. To begin, a target vector, X_i, which is to be potentially improved, is selected. Next, a multi-step process attempts to replace X_i with an improved offspring vector, ${\mathit{X}}_i^{\prime }$, as defined by the misfit function, g. First, the best-fit vector, X_best (i.e., lowest value of g(X_i)), in the current population is copied and modified by adding a scaled value of the difference between two other randomly selected vectors in the population (${\mathit{X}}_{r_{\mathrm{1}}}$ and ${\mathit{X}}_{r_{\mathrm{2}}}$) as follows:

$$\begin{array}{}\text{(6)}& {\mathit{U}}_i={\mathit{X}}_{\mathrm{best}}+\mathit{\beta }({\mathit{X}}_{r_{\mathrm{1}}}-{\mathit{X}}_{r_{\mathrm{2}}}).\end{array}$$

Here, β is a uniformly random value between 0.5 and 1.0. Next, U_i is combined with X_best to produce the offspring vector ${\mathit{X}}_i^{\prime }$. This combination is performed by, for each element of U_i, sampling a uniform distribution on [0,1) and replacing the element of X_i with the corresponding element of U_i if the value is less than 0.7, the default recombination constant (see Storn and Price, 1997; Virtanen et al., 2020). In this roundabout way, the trial vector's generation is informed by the existing population.

At this point, X_i has been generated, and if $g\left({\mathit{X}}_i^{\prime }\right)<g\left({\mathit{X}}_i\right)$, X_i is replaced with ${\mathit{X}}_i^{\prime }$ in the population. Once the above improvement steps have been performed for every vector in the population, we advance a generation, and population metrics are calculated. Based on the results of these metrics, either another generation is calculated or the algorithm returns the best-fit vector from the population.

5 Synthetic data experiments

Synthetic step-heating diffusion experiments allow for the validation of optimization methods like those described above. For a given set of diffusion kinetics parameters, one can calculate the expected number of moles released at each step, and these release fractions can then be used as an input to the optimization algorithm to search for the known diffusion kinetics parameters. If the correct parameters are returned, the optimization algorithm is validated. To evaluate the ability of our differential evolution routine (Storn and Price, 1997) in comparison to existing methods (Lovera, 1992; Lovera et al., 1997), we present the results of two such synthetic data experiments (Fig. 1).

5.1 Synthetic data methods

In calculating the number of moles released for each step of a given heating schedule, two assumptions were made. First, E_a was assumed to be common to all domains (Lovera et al., 1997). Second, we required the $\ln (\frac{D_{\mathrm{0}}}{a^{\mathrm{2}}}{)}_i$ values of all domains to differ by at least 0.25 natural log units. For the purposes of this experiment, any two domains with $\ln (\frac{D_{\mathrm{0}}}{a^{\mathrm{2}}}{)}_i$ values differing by less than 0.25 were considered to be well represented by a model with one fewer domain.

To guide our choices for the synthetic heating schedule and prescribed diffusion kinetics parameters, we used existing literature on K-feldspar diffusion experiments (e.g., Lovera et al., 1997). We begin by defining Experiment A where the heating schedule and M_tot were selected from the Lovera et al. (1997) N13 K-feldspar experiment because this experiment has been published many times (e.g., Lovera et al., 1997; Harrison and Lovera, 2014; Reiners et al., 2017) and because it represents a typical K-feldspar diffusion experiment. E_a was prescribed by randomly sampling from a Gaussian distribution with a mean of 192.5 kJ mol⁻¹ and σ of 25, mirroring the Lovera et al. (1997) database of K-feldspar E_a values. The $\ln (\frac{D_{\mathrm{0}}}{a^{\mathrm{2}}}{)}_i$ values were selected uniformly randomly from a range of 5–25 natural log units. Finally, the ϕ_i values were prescribed randomly, requiring only that the sum of the gas fractions equal 1. The largest of these ϕ_i values were intentionally placed in the largest domains, reflecting common K-feldspar behavior.

Experiment B was then designed to intentionally violate the assumption made by the Lovera et al. (1997) fitting algorithm that no domain should significantly degas during early heating steps. This was done by taking the kinetics from Experiment A, removing 1 % of the total gas from each domain, and then placing this gas in a new highly diffusive domain with $\ln (\frac{D}{a^{\mathrm{2}}}{)}_{\mathrm{1}}=\mathrm{23.4}$, and ϕ₁=0.01.

To calculate the fractional releases and number of moles released after each heating step, we began with the equations for plane sheet geometry outlined in Ginster and Reiners (2018), which relate each heating step's duration and temperature as well as the domain's $\ln (\frac{D_{\mathrm{0}}}{a^{\mathrm{2}}}{)}_i$ to a fractional release from that domain ((F_dom)_i). To determine the total gas fraction released for a sample from each heating step (F_i), and not just for a specific domain ((F_dom)_i), we summed the contributions from each domain as follows: $F_i={\sum }_{i=\mathrm{1}}^n(F_{\mathrm{dom}}{)}_i{\mathit{\varphi }}_i$. In finding the number of moles released at each step, we calculated $M_i=F_i\cdot M_{\mathrm{tot}}$, where M_tot is the total number of atoms measured in the Lovera et al. (1997) N13 K-feldspar experiment. To approximate uncertainties, we simply multiplied the percent error in the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age from each step of the N13 experiment by M_i.

Using these synthetic degassing datasets, we attempted to solve for the prescribed diffusion kinetics of both experiments using our MDD Tool Kit (https://github.com/dgorin1/mddtoolkit) optimization method (Lovera et al., 1997). Because K-feldspar is known to melt above 1100 °C (Luo et al., 2014), we excluded all calculated heating steps above 1100 °C in our misfit calculations.

5.2 Synthetic data results and discussion

The MDD Tool Kit (https://github.com/dgorin1/mddtoolkit) method successfully quantified the diffusion kinetics of Experiments A and B, returning the correct E_a to within 0.02 kJ mol⁻¹ (Table 4, Fig. 1). While the MDD Tool Kit routine did not perform as well in capturing ϕ₆₋₈ this is unsurprising since a higher percentage of the gas in those domains was released during the high-temperature steps excluded from the fitting exercise.

The Lovera et al. (1997) algorithm correctly quantified the E_a of Experiment A, but it underestimated that of Experiment B by 9% (Table 4). While we did not implement a Levenberg–Marquardt method (Press, 2007) to solve for each $\ln (\frac{D_{\mathrm{0}}}{a^{\mathrm{2}}}{)}_i$ and ϕ_i value, these parameters could not be correct given the incorrectly predicted E_a.

Our synthetic experiment results clearly demonstrate that the Lovera et al. (1997) algorithm can underestimate E_a (Fig. 1). Any sample that contains at least one domain which loses a significant portion of its gas during initial heating steps is prone to this bias. Importantly, this error is not bidirectional; the Lovera algorithm will systematically underestimate E_a. In contrast, because the optimized E_a does not depend on a linear regression but instead fits E_a as a free parameter, the MDD Tool Kit method appears to avoid this bias. This finding suggests that the use of the MDD Tool Kit to reanalyze any MDD model thermochronology study based on the Lovera approach would either produce similar results or systematically higher temperatures through geologic time.

6 Case study: Wong et al. (2023) field validation of the MDD model

To assess the accuracy of the MDD model, Wong et al. (2023) conducted a field validation study of K-feldspar ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ thermochronology at the Grayback normal fault block, AZ, USA. This field site is well studied, and several independent thermochronometers have been used there to measure geothermal gradients for application to models of continental extension (Howard and Foster, 1996; Wong et al., 2015).

The field validation primarily relies on four samples from various stratigraphic positions within the Tea Cup pluton, which intruded into the overlying ∼ 1.4 Ga Ruin Granite at ∼ 70 Ma (Banks et al., 1972). Samples originally from ∼ 12 km below the paleo-surface are accessible since these formations were tilted ∼ 90° to the east during mid-Oligocene extension of the Grayback Fault (Banks et al., 1972; Wong et al., 2023). The straightforward stratigraphy and availability of independent thermochronometry data make this location optimal for the validation study.

Wong et al. (2023) found good agreement between their best-fitting ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ t–T results and existing estimates of the paleo-geothermal gradient prior to mid-Oligocene extension (Howard and Foster, 1996; Wong et al., 2015, 2023). A secondary assessment was provided by fission track and $(\mathrm{U}-\mathrm{Th})/\mathrm{He}$ ages in zircon and apatite. Although the observed ages generally do not directly overlap in time with the K-feldspar MDD histories, the results do not preclude one another. In total, their study appears to validate the MDD model. To determine whether the MDD Tool Kit (https://github.com/dgorin1/mddtoolkit) method of fitting diffusion kinetics is consistent with these independent observations, we reanalyze their ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ stepwise degassing data using our new method.

6.1 Methods

6.1.1 Diffusion kinetics

In determining the diffusion kinetics of each sample, Wong et al. (2023) used a modified version of the Lovera et al. (1997) fitting algorithm. Their primary modification was to fit an unweighted linear regression to the beginning heating steps from each sample to determine E_a and $\ln \frac{D_{\mathrm{0}}}{a_{\mathrm{0}}^{\mathrm{2}}}$. Instead of using the goodness-of-fit metric of Lovera et al. (1997), Wong et al. (2023) varied how many heating steps they included in the regression and explored the resultant effect of E_a on the constrained t–T path. They then chose the resultant E_a which most closely agreed with independent thermochronological data, noting that the choice of E_a mainly affected the absolute temperature predictions and not the form of a t–T history. Although this was a routine choice given the tools available at the time of their publication, such user-defined choices may introduce bias, a concern that the MDD Tool Kit (https://github.com/dgorin1/mddtoolkit) inherently circumvents.

In our application of the MDD Tool Kit (https://github.com/dgorin1/mddtoolkit) to Wong et al. (2023) diffusion experiment results, we apply both the χ² and %_frac misfit statistics defined above. As previously described, heating steps greater than 1100°C were excluded from our misfit calculations.

6.1.2 Thermal paths

When applying K-feldspar MDD modeling, a number of thermal histories are generated using a Monte Carlo approach to predict an ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age spectrum for each prescribed t–T path and using the sample's apparent diffusion kinetics. A misfit is then calculated between the modeled and measured ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age spectra to determine the fitness of a particular t–T path. In our study, we generated 30 000 thermal paths per sample and calculated radiogenic ⁴⁰Ar production and diffusive loss for each t–T path using a Crank–Nicolson discretization of the diffusion equation for an infinite sheet geometry (Crank, 1975). We then focus our analysis on the 100 best-fitting t–T paths for each sample based on the following misfit statistic:

$$\begin{array}{}\text{(7)}& {\mathit{\chi }}_{\mathrm{Age}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N}}\cdot \sum _{i=\mathrm{1}}^N{\left({\displaystyle \frac{({\mathrm{Age}}_{\mathrm{Measured}}{)}_i-({\mathrm{Age}}_{\mathrm{Modeled}}{)}_i}{({\mathit{\sigma }}_{{\mathrm{Age}}_{\mathrm{Measured}}}{)}_i}}\right)}^{\mathrm{2}},\end{array}$$

where N is the number of steps included in the thermal modeling, and σ is the reported uncertainty in the final ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age. As is common practice, Wong et al. (2023) ignored measurements thought to be contaminated with excess ⁴⁰Ar (⁴⁰Ar_E; Lovera et al., 2002); for comparison with their results, we exclude the same data.

6.2 Results

Our predicted E_a values (Tables 5–7, Fig. 3), regardless of misfit statistic, are systematically higher than those found by Wong et al. (2023) and are at the upper range of those published in the Lovera et al. (1997) database of K-feldspar diffusion experiments (Fig. 1). Given the bias inherent to the regression fitting method (Fig. 1), it is not surprising that our predicted E_a values are systematically higher, since the values of E_a were optimized using the entire dataset, rather than prescribed by linear regression to a subset of user-defined heating steps tuned to agree with independent thermochronological data. And, although we defined a large search space for E_a, we find a smaller inter-sample range in predicted E_a than published by Wong et al. (2023), consistent with the Lovera et al. (1997) observation that the variation of K-feldspar E_a values decreases when more gas is included in the E_a calculation (Tables 5–7; Fig. 3).

Wong et al. (2023)

Table 5 Results from sample GR-1.

Wong et al. (2023)

Table 6 Results from sample GR-2.

Wong et al. (2023)

Table 7 Results from sample GR-27.

Figure 3 Arrhenius plots showing our reanalysis of Wong et al. (2023) samples with both the χ² and %_frac misfit statistics. Red lines represent the diffusion kinetics of each individual domain, and their thicknesses are proportional to the domain's ϕ value. Gray circles show the experimental results, and black dots show the MDD model predictions. All E_a values are given in kJ mol⁻¹.

Despite this finding, the choice of misfit statistic appears to influence the calculated diffusion kinetics of a given sample. In our analysis, choosing a χ² misfit instead of the %_frac misfit statistic led to intra-sample differences in the predicted E_a between 4.4–22.4 kJ mol⁻¹ (Tables 8 and 5, respectively). For example, there is a 22.4 kJ mol⁻¹ difference between the best-fit E_a value for sample GR-2 when using the χ² misfit compared to the %_frac misfit (Table 6). This disparity may provide an estimate of the uncertainty in the E_a, with the true E_a lying anywhere between those values. We recommend that investigators using the MDD Tool Kit consider model fits generated by both the χ² and the %_frac misfit statistics as equally plausible, given the lack of justification for choosing one over the other.

Wong et al. (2023)

Table 8 Results from sample GR-8.

An additional noteworthy model behavior emerged during our analysis. In three cases, the MDD Tool Kit (https://github.com/dgorin1/mddtoolkit) identified at least one retentive domain that exhibited no gas diffusion throughout the simulated experiment, as indicated in Tables 4–7. The total ϕ value(s) of these domains typically equaled the gas quantity contained in the heating steps above 1100 °C. While the result that no gas diffused from the most retentive domain (or domains) during the simulated experiment may appear to suggest a poor model choice, this is not necessarily the case. Since K-feldspar begins to melt above 1100 °C, the MDD model is simply not applicable above this temperature. Furthermore, any domain retaining all its gas during simulation will yield identical release fractions, thereby making the model's fit insensitive to retentivities above a specific threshold (assuming all other parameters remain constant). Given our deliberate inclusion of a wide search space for each diffusion parameter and the stochastic nature of our optimization algorithm, such behavior is to be expected.

6.3 Reinterpretation of field calibration

Although our reinterpretation of these data largely finds t–T paths of similar form to those of Wong et al. (2023), our absolute temperature predictions are about ∼ 60–75 °C higher (Fig. 4). To demonstrate that our approach yields results that are equally permissible, we compare them to the independent t–T constraints that exist at the Grayback Fault.

Figure 4 Resulting best-fit thermal pathways and ⁴⁰Ar-³⁹Ar age spectra from our reanalysis of the Wong et al. (2023) field validation of the MDD method in K-feldspar. Blue shaded regions in the thermal history plots correspond to the χ² misfit statistic. The orange shaded regions correspond to the %_frac misfit statistic. More specifically, the shaded regions represent a 1σ deviation from the median of the top 100 best-fitting thermal paths for a given sample (plotted in bold). The dotted line and gray regions represent Wong et al. (2023) predictions for the same samples. The ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age spectrum plots show the predicted ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages for the median χ² (blue) and %_frac (orange) thermal paths. Gray boxes represent the measured values from Wong et al. (2023). Heating steps excluded from the modeling exercise are made transparent.

Wong et al. (2023) proposed three criteria to assess the validity of an MDD model: (i) the thermal histories should be consistent with the stratigraphic heights of the samples; (ii) the form of the predictions should match prior work, including the timing, rate, and duration of cooling events; and (iii) the absolute temperatures should agree with those predicted by previous estimates of the pre-extensional geothermal gradient (Howard and Foster, 1996; Wong et al., 2015). The constrained t–T paths shown in Wong et al. (2023) and our reanalysis of the same dataset using the MDD Tool Kit (https://github.com/dgorin1/mddtoolkit) both meet these criteria (Fig. 4). The predicted t–T paths for all four samples generally align with their respective stratigraphic positions; samples from higher stratigraphic heights reflect the oldest portions of the thermal histories (Fig. 4). Our findings reveal that the sample situated at the highest stratigraphic level (GR-1) cooled below its Ar closure temperature around ∼ 55 Ma, indicating relatively rapid cooling between ∼ 70–55 Ma. This trend is consistent with the Howard and Foster (1996) interpretation that the most shallow depths of the pluton experienced rapid cooling as it equilibrated with ambient temperatures.

MDD results for the next sample (GR-2), slightly deeper than GR-1, suggest a gradual cooling rate of approximately ∼ 5 °C Ma⁻¹ between ∼ 55–30 Ma, consistent with the prediction of ∼ 4–6 °C Ma⁻¹ by Howard and Foster (1996).

The two samples at the lowest stratigraphic levels, GR-27 and GR-8, exhibit similar temperature histories. The median pathway for GR-27 calculated with the χ² misfit indicates slightly higher temperatures compared to GR-8, despite the latter occupying a slightly lower stratigraphic depth. However, given that their thermal pathways are within 1σ of each other, it suggests that this technique may not be capable of resolving subtle differences between these samples. Moreover, the presence of significant Ar_E (Lovera et al., 2002) in sample GR-8 introduces additional uncertainty into its Paleogene thermal history. Despite the presence of Ar_E (Lovera et al., 2002), both of these deep stratigraphic samples demonstrate rapid cooling commencing at around ∼ 27–28 Ma, consistent with the predicted timing of this cooling based on apatite fission track (AFT) ages (Howard and Foster, 1996).

Finally, all four of our newly calculated MDD models predict similar absolute temperatures to those estimated for the paleo-geothermal gradient along the Tea Cup pluton prior to onset of extension along the Grayback Fault. The paleo-gradient has been estimated from temperature predictions made at 27 Ma at a variety of paleodepths using AFT (Howard and Foster, 1996). While our estimates are consistently warmer than those calculated by Wong et al. (2023), they are still within the uncertainty bounds of the predicted thermal gradient (Fig. 5). The overall agreement between the rates and timing of cooling, as well as the relationship between relative temperatures over time and the samples' stratigraphic order, indicates that the MDD Tool Kit (https://github.com/dgorin1/mddtoolkit) generates t–T results that are not excluded by the existing geologic and independent thermochronology data.

Figure 5 Predicted paleotemperatures at 27 Ma from samples GR-2, GR-27, and GR-8. The gray symbols represent Wong et al. (2023) results, the orange symbols represent our results from the %_frac misfit statistic, and the blue symbols represent our results from the χ² misfit statistic.

7 Conclusions

The multiple-diffusion domain model for ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ thermochronology is valuable because it allows one to constrain a mineral's continuous thermal history through geologic time (McDougall and Harrison, 1999). However, the methods previously used to empirically fit stepwise degassing data required unverifiable assumptions about the grain size distribution of a sample. Further, the commonly used method of quantifying E_a by fitting an uncertainty-weighted linear regression to the lowest-temperature degassing steps of an experiment does not reliably return the correct values but instead inadvertently introduces a user bias. When this method fails to predict the correct value, it underestimates the true value. To address these limitations, we present a new numerical routine that does not require fitting a linear regression to a user-defined subset of heating steps to quantify a sample's E_a. Our new method utilizes a differential evolution algorithm to robustly search the MDD parameter space and solve for all parameters simultaneously. The code, entitled MDD Tool Kit (https://github.com/dgorin1/mddtoolkit), is open-source, pip-installable, and available on GitHub (https://github.com/dgorin1/mddtoolkit). To evaluate its validity, we apply this new method to reinterpret the dataset published in the Wong et al. (2023) field validation of the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ K-feldspar MDD system at the Grayback Fault. The diffusion kinetics fit by the MDD Tool Kit (https://github.com/dgorin1/mddtoolkit) predict t–T paths that are consistent with independent observations and geologic constraints and, on average, 60–75 °C warmer than those previously published. We attribute this temperature difference to biases potentially introduced by the previously used modeling strategies.

Appendix A: Supplemental tables

For completeness, we include the optimization inputs used to reanalyze the experiments performed by Wong et al. (2023) (Tables A1–A4).

Table A1 Sample GR-1 input (Wong et al., 2023).

Table A2 Sample GR-2 input (Wong et al., 2023).

Table A3 Sample GR-8 input (Wong et al., 2023).

Table A4 Sample GR-27 input (Wong et al., 2023).