the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# An optimization tool for identifying multiple-diffusion domain model parameters

### Andrew L. Gorin

### Joshua M. Gorin

### Marie Bergelin

### David L. Shuster

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}$.

- Article
(2920 KB) - Full-text XML
- BibTeX
- EndNote

${}^{\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 ^{40}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 ^{39}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).

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

where *T* is the absolute temperature (K), *D* is diffusivity at *T* (cm^{2} s^{−1}), *D*_{0} is diffusivity at infinite *T*, *E*_{a} is activation energy (kJ mol^{−1}), *R* is the gas constant (kJ mol^{−1} K^{−1}), and *a* is the radius of the diffusion domain (cm) (see Table 1 for variable definitions).

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

When consistent with volume diffusion through a single diffusion domain, a step-heating experiment will produce a linear Arrhenius relationship between calculated values of $\mathrm{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 ^{3}He in apatite and ^{3}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*_{0}. 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 2*n*−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 2*n*−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 ^{40}Ar within minerals over geologic time occurs primarily by volume diffusion. This assumption predicts that low ^{40}Ar concentrations should exist near the outer edge of a crystal (i.e., for a boundary condition of nearly zero ^{40}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, ^{40}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.

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 ${\mathrm{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).

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:

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 ${\mathrm{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 ${\mathrm{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 ^{39}Ar used in the fitting exercise might lead to more precise estimates of diffusion kinetics.

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.

We apply two such misfit statistics to the data analyzed in this study: one error-weighted, *χ*^{2}, and another where all points are weighted equally (*%*_{frac}). Our *χ*^{2} misfit accounts for the uncertainty in the ^{39}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 *χ*^{2} misfit statistic is thereby defined as follows:

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:

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

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:

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.

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 $\mathrm{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 $\mathrm{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^{−1} and *σ* of 25, mirroring the Lovera et al. (1997) database of K-feldspar *E*_{a} values. The $\mathrm{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 $\mathrm{ln}(\frac{D}{{a}^{\mathrm{2}}}{)}_{\mathrm{1}}=\mathrm{23.4}$, and *ϕ*_{1}=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 $\mathrm{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^{−1} (Table 4, Fig. 1). While the MDD Tool Kit routine did not perform as well in capturing *ϕ*_{6−8} 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 $\mathrm{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.

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 $\mathrm{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 *χ*^{2} 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 ^{40}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:

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 ^{40}Ar (^{40}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).

Despite this finding, the choice of misfit statistic appears to influence the calculated diffusion kinetics of a given sample. In our analysis, choosing a *χ*^{2} misfit instead of the *%*_{frac} misfit statistic led to intra-sample differences in the predicted *E*_{a} between 4.4–22.4 kJ mol^{−1} (Tables 8 and 5, respectively). For example, there is a 22.4 kJ mol^{−1} difference between the best-fit *E*_{a} value for sample GR-2 when using the *χ*^{2} 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 *χ*^{2} and the *%*_{frac} misfit statistics as equally plausible, given the lack of justification for choosing one over the other.

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.

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^{−1} between ∼ 55–30 Ma, consistent with the prediction of ∼ 4–6 *°*C Ma^{−1} 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 *χ*^{2} 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.

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.

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

The code released in this publication is available at our GitHub link https://github.com/dgorin1/mddtoolkit (last access: 14 October 2024) and Zenodo https://doi.org/10.5281/zenodo.13921573 Gorin and Gorin, 2024. All other data in this paper are either provided in tabular format or are referenced.

ALG and DLS conceived of the study. ALG and JMG designed the software. ALG, DLS, and MB contributed to the theory behind the software. ALG and DLS led manuscript writing. ALG constructed figures. All authors reviewed and edited the manuscript.

The contact author has declared that none of the authors has any competing interests.

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

We acknowledge Oscar Lovera for insightful discussions about the theory of the MDD model. We also thank Martin Wong for providing us with data for calculating measured ^{39}Ar uncertainties for his dataset.

This research has been supported by the Ann and Gordon Getty foundation.

This paper was edited by Klaus Mezger and reviewed by Simon P. Kelley and one anonymous referee.

Axen, G. J., Grove, M., Stockli, D., Lovera, O. M., Rothstein, D. A., Fletcher, J. M., Farley, K., and Abbott, P. L.: Thermal evolution of Monte Blanco dome: Low‐angle normal faulting during Gulf of California rifting and late Eocene denudation of the eastern Peninsular Ranges, Tectonics, 19, 197–212, https://doi.org/10.1029/1999TC001123, 2000. a

Banks, N. G., Cornwall, H. R., Silberman, M. L., Creasey, S. C., and Marvin, R. F.: Chronology of Intrusion and Ore Deposition at Ray, Arizona; Part I, K-Ar Ages, Econ. Geol., 67, 864–878, https://doi.org/10.2113/gsecongeo.67.7.864, 1972. a, b

Crank, J.: The Mathematics of Diffusion, Oxford University Press,2nd Edn., ISBN 0198533446, 1975. a, b, c

Dodson, M. H.: Closure temperature in cooling geochronological and petrological systems, Contrib. Miner. Petrol., 40, 259–274, https://doi.org/10.1007/BF00373790, 1973. a

Fechtig, H. and Kalbitzer, S.: The Diffusion of Argon in Potassium-Bearing Solids, in: Potassium Argon Dating, pp. 68–107, Springer Berlin Heidelberg, Berlin, Heidelberg, ISBN 978-3-642-87895-4, https://doi.org/10.1007/978-3-642-87895-4_4, 1966. a, b

Flude, S., Halton, A. M., Kelley, S. P., Sherlock, S. C., Schwanethal, J., and Wilkinson, C. M.: Observation of centimetre-scale argon diffusion in alkali feldspars: implications for ^{40}Ar/^{39} Ar thermochronology, Geol. Soc. Lond. Spec. Publ., 378, 265–275, https://doi.org/10.1144/SP378.25, 2014. a

Gillespie, A. R., Huneke, J. C., and Wasserburg, G. J.: An assessment of ^{40}Ar–^{39}Ar dating of incompletely degassed xenoliths, J. Geophys. Res.-Sol. Ea., 87, 9247–9257, https://doi.org/10.1029/JB087iB11p09247, 1982. a

Ginster, U. and Reiners, P. W.: Error Propagation in the Derivation of Noble Gas Diffusion Parameters for Minerals From Step Heating Experiments, Geochem. Geophy. Geosy., 19, 3706–3720, https://doi.org/10.1029/2018GC007531, 2018. a, b, c

Gorin, J. and Gorin, D.: dgorin1/mddtoolkit: MDD Tool Kit (0.1.0), Zenodo [code], https://doi.org/10.5281/zenodo.13921573, 2024. a

Grove, M., Lovera, O., and Harrison, M.: Late Cretaceous cooling of the east-central Peninsular Ranges Batholith (33 degrees N); relationship to La Posta Pluton emplacement, Laramide shallow subduction, and forearc sedimentation, in: Tectonic evolution of northwestern Mexico and the Southwestern USA, Geol. Soc. Am., ISBN 978-0-8137-2374-7, https://doi.org/10.1130/0-8137-2374-4.355, 2003. a

Harrison, T. M. and Lovera, O. M.: The multi-diffusion domain model: past, present and future, Geol. Soc. Lond. Spec. Publ., 378, 91–106, https://doi.org/10.1144/SP378.9, 2014. a, b, c

Harrison, T. M., Heizler, M. T., Lovera, O. M., Chen Wenji, and Grove, M.: A chlorine disinfectant for excess argon released from K-felsspar during step heating, Earth Planet. Sc. Lett., 123, 95–104, https://doi.org/10.1016/0012-821X(94)90260-7, 1994. a

Harrison, T. M., Copeland, P., Kidd, W. S. F., and Lovera, O. M.: Activation of the Nyainqentanghla Shear Zone: Implications for uplift of the southern Tibetan Plateau, Tectonics, 14, 658–676, https://doi.org/10.1029/95TC00608, 1995. a

Howard, K. A. and Foster, D. A.: Thermal and unroofing history of a thick, tilted Basin‐and‐Range crustal section in the Tortilla Mountains, Arizona, J. Geophys. Res.-Sol. Ea., 101, 511–522, https://doi.org/10.1029/95JB02909, 1996. a, b, c, d, e, f, g

Iman, R. L., Helton, J. C., and Campbell, J. E.: An Approach to Sensitivity Analysis of Computer Models: Part I – Introduction, Input Variable Selection and Preliminary Variable Assessment, J. Qual. Technol., 13, 174–183, https://doi.org/10.1080/00224065.1981.11978748, 1981. a

Kirby, E., Reiners, P. W., Krol, M. A., Whipple, K. X., Hodges, K. V., Farley, K. A., Tang, W., and Chen, Z.: Late Cenozoic evolution of the eastern margin of the Tibetan Plateau: Inferences from ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}$Ar and (U‐Th)/He thermochronology, Tectonics, 21, 1-1–1-20, https://doi.org/10.1029/2000TC001246, 2002. a

Lampinen, J.: A constraint handling approach for the differential evolution algorithm, in: Proceedings of the 2002 Congress on Evolutionary Computation, CEC'02 (Cat. No.02TH8600), 2, 1468–1473, IEEE, Honolulu, HI, USA, ISBN 978-0-7803-7282-5, https://doi.org/10.1109/CEC.2002.1004459, 2002. a

Lovera, O. M.: Computer programs to model ^{40}Ar/^{39}Ar diffusion data from multidomain samples, Comput. Geosci., 18, 789–813, https://doi.org/10.1016/0098-3004(92)90025-M, 1992. a

Lovera, O. M., Richter, F. M., and Harrison, T. M.: The ^{40}Ar/^{39} Ar thermochronometry for slowly cooled samples having a distribution of diffusion domain sizes, J. Geophys. Res.-Sol. Ea., 94, 17917–17935, https://doi.org/10.1029/JB094iB12p17917, 1989. a, b, c, d

Lovera, O. M., Heizler, M. T., and Harrison, T. M.: Argon diffusion domains in K-feldspar II: kinetic properties of MH-10, Contrib. Miner. Petrol., 113, 381–393, https://doi.org/10.1007/BF00286929, 1993. a, b

Lovera, O. M., Grove, M., Mark Harrison, T., and Mahon, K.: Systematic analysis of K-feldspar step heating results: I. Significance of activation energy determinations, Geochim. Cosmochim. Ac., 61, 3171–3192, https://doi.org/10.1016/S0016-7037(97)00147-6, 1997. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac, ad, ae, af, ag, ah, ai

Lovera, O. M., Grove, M., and Harrison, T.: Systematic analysis of K-feldspar 40Ar/39Ar step heating results II: relevance of laboratory argon diffusion properties to nature, Geochim. Cosmochim. Ac., 66, 1237–1255, https://doi.org/10.1016/S0016-7037(01)00846-8, 2002. a, b, c, d, e

Luo, G., Nie, X., Wu, S., and Wang, Y.: Influence of compound silicate gangue containing potassium in baiyunebo iron ore on solid-phase reactions during sintering process, in: Recent Advances in Structural Integrity Analysis – Proceedings of the International Congress (APCF/SIF-2014), 473–476, Elsevier, ISBN 978-0-08-100203-2, https://doi.org/10.1533/9780081002254.473, 2014. a

McDougall, I. and Harrison, T. M.: Geochronology and Thermochronology by the ^{40}Ar/^{39}Ar Method, Oxford University Press, ISBN 9780195109207, 1999. a, b, c, d, e, f

Parsons, I., Brown, W. L., and Smith, J. V.: 40Ar/39Ar thermochronology using alkali feldspars: real thermal history or mathematical mirage of microtexture?, Contrib. Miner. Petrol., 136, 92–110, https://doi.org/10.1007/s004100050526, 1999. a, b, c

Popov, D. V. and Spikings, R. A.: Diffusion vs. fluid alteration in alkali feldspar ^{40}Ar/^{39}Ar thermochronology: Does cross-correlation of log(r/r_{0}) and age spectra validate thermal histories?, Chem. Geol., 539, 119506, https://doi.org/10.1016/j.chemgeo.2020.119506, 2020. a, b, c

Popov, D. V., Spikings, R. A., and Kouzmanov, K.: Pathways for 39Ar loss during step-heating of alkali feldspar megacrysts from the Shap granite (UK): Combined evidence from diffusion experiments and characterisation of heating-induced texture modifications, Chem. Geol., 547, 119677, https://doi.org/10.1016/j.chemgeo.2020.119677, 2020a. a, b, c

Popov, D. V., Spikings, R. A., Scaillet, S., O'Sullivan, G., Chew, D., Badenszki, E., Daly, J. S., Razakamanana, T., and Davies, J. H.: Diffusion and fluid interaction in Itrongay pegmatite (Madagascar): Evidence from in situ ^{40}Ar/^{39}Ar dating of gem-quality alkali feldspar and U Pb dating of protogenetic apatite inclusions, Chem. Geol., 556, 119841, https://doi.org/10.1016/j.chemgeo.2020.119841, 2020b. a, b

Press, W. H. (Ed.): Numerical recipes: the art of scientific computing, Cambridge University Press, Cambridge, UK, New York, 3rd Edn., ISBN 978-0-521-88068-8, 978-0-521-88407-5, 978-0-521-70685-8, 2007. a, b, c, d

Qiang, J. and Mitchell, C.: A Unified Differential Evolution Algorithm for Global Optimization, https://www.osti.gov/servlets/purl/1163659 (last access: 11 October 2024), 2014. a

Quidelleur, X., Grove, M., Lovera, O. M., Harrison, T. M., Yin, A., and Ryerson, F. J.: Thermal evolution and slip history of the Renbu Zedong Thrust, southeastern Tibet, J. Geophys. Res.-Sol. Ea., 102, 2659–2679, https://doi.org/10.1029/96JB02483, 1997. a

Reiners, P. W.: Past, Present, and Future of Thermochronology, Rev. Miner. Geochem., 58, 1–18, https://doi.org/10.2138/rmg.2005.58.1, 2005. a, b

Reiners, P. W. and Farley, K. A.: Helium diffusion and (U–Th)/He thermochronometry of titanite, Geochim. Cosmochim. Ac., 63, 3845–3859, https://doi.org/10.1016/S0016-7037(99)00170-2, 1999. a

Reiners, P. W., Spell, T. L., Nicolescu, S., and Zanetti, K. A.: Zircon (U-Th)/He thermochronometry: He diffusion and comparisons with ^{40}Ar/^{39}Ar dating, Geochim. Cosmochim. Ac., 68, 1857–1887, https://doi.org/10.1016/j.gca.2003.10.021, 2004. a

Reiners, P. W., Carlson, R. W., Renne, P. R., Cooper, K. M., Granger, D. E., McLean, N. M., and Schoene, B.: Geochronology and thermochronology, John Wiley & Sons, https://doi.org/10.1002/9781118455876, 2017. a

Shirvell, C. R., Stockli, D. F., Axen, G. J., and Grove, M.: Miocene‐Pliocene exhumation along the west Salton detachment fault, southern California, from (U‐Th)/He thermochronometry of apatite and zircon, Tectonics, 28, 2007TC002172, https://doi.org/10.1029/2007TC002172, 2009. a

Shuster, D. L., Farley, K. A., Sisterson, J. M., and Burnett, D. S.: Quantifying the diffusion kinetics and spatial distributions of radiogenic 4He in minerals containing proton-induced 3He, Earth Planet. Sc. Lett.s, 217, 19–32, https://doi.org/10.1016/S0012-821X(03)00594-6, 2004. a

Spell, T. L., McDougall, I., and Tulloch, A. J.: Thermochronologic constraints on the breakup of the Pacific Gondwana margin: The Paparoa metamorphic core complex, South Island, New Zealand, Tectonics, 19, 433–451, https://doi.org/10.1029/1999TC900046, 2000. a

Spikings, R. A. and Popov, D. V.: Thermochronology of Alkali Feldspar and Muscovite at *T*>150 °C Using the ^{40}Ar/^{39}Ar Method: A Review, Minerals, 11, 1025, https://doi.org/10.3390/min11091025, 2021. a, b

Storn, R. and Price, K.: Differential Evolution – A simple and Efficient Heuristic for Global Optimization over Continuous Spaces, J. Global Optim., 11, 341–359, https://doi.org/10.1023/A:1008202821328, 1997. a, b, c, d

Turner, G.: The Distribution of Potassium and Argon in Chondrites, in: Origin and Distribution of the Elements, Elsevier, 387–398, ISBN 978-0-08-012835-1, https://doi.org/10.1016/B978-0-08-012835-1.50039-8, 1968. a

Villa, I. M.: Multipath Ar transport in K-feldspar deduced from isothermal heating experiments, Earth Planet. Sc. Lett., 122, 393–401, https://doi.org/10.1016/0012-821X(94)90010-8, 1994. a

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., Van Der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, I., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., Van Mulbregt, P., SciPy 1.0 Contributors, Vijaykumar, A., Bardelli, A. P., Rothberg, A., Hilboll, A., Kloeckner, A., Scopatz, A., Lee, A., Rokem, A., Woods, C. N., Fulton, C., Masson, C., Häggström, C., Fitzgerald, C., Nicholson, D. A., Hagen, D. R., Pasechnik, D. V., Olivetti, E., Martin, E., Wieser, E., Silva, F., Lenders, F., Wilhelm, F., Young, G., Price, G. A., Ingold, G.-L., Allen, G. E., Lee, G. R., Audren, H., Probst, I., Dietrich, J. P., Silterra, J., Webber, J. T., Slavič, J., Nothman, J., Buchner, J., Kulick, J., Schönberger, J. L., De Miranda Cardoso, J. V., Reimer, J., Harrington, J., Rodríguez, J. L. C., Nunez-Iglesias, J., Kuczynski, J., Tritz, K., Thoma, M., Newville, M., Kümmerer, M., Bolingbroke, M., Tartre, M., Pak, M., Smith, N. J., Nowaczyk, N., Shebanov, N., Pavlyk, O., Brodtkorb, P. A., Lee, P., McGibbon, R. T., Feldbauer, R., Lewis, S., Tygier, S., Sievert, S., Vigna, S., Peterson, S., More, S., Pudlik, T., Oshima, T., Pingel, T. J., Robitaille, T. P., Spura, T., Jones, T. R., Cera, T., Leslie, T., Zito, T., Krauss, T., Upadhyay, U., Halchenko, Y. O., and Vázquez-Baeza, Y.: SciPy 1.0: fundamental algorithms for scientific computing in Python, Nat. Meth., 17, 261–272, https://doi.org/10.1038/s41592-019-0686-2, 2020. a, b, c

Wagner, A.: Cooling ages derived by apatite fission-track, mica Rb-Sr and K-Ar dating; the uplift and cooling history of the Central Alps, Mem. Inst. Geol. Mineral., Univ. of Padova, 30, 1–27, https://cir.nii.ac.jp/crid/1570854175293805312 (last access: 11 October 2024), 1977. a

Warnock, A. C. and Zeitler, P. K.: ^{40}Ar/^{39}Ar thermochronometry of K-feldspar from the KTB borehole, Germany, Earth Planet. Sc. Lett., 158, 67–79, https://doi.org/10.1016/S0012-821X(98)00044-2, 1998. a

Weirich, J., Isachsen, C., Johnson, J., and Swindle, T.: Variability of diffusion of argon in albite, pyroxene, and olivine in shocked and unshocked samples, Geochim. Cosmochim. Ac., 77, 546–560, https://doi.org/10.1016/j.gca.2011.10.040, 2012. a

Wong, M. S., Gleason, D. M., O’Brien, H. P., and Idleman, B. D.: Confirmation of a low pre-extensional geothermal gradient in the Grayback normal fault block, Arizona: Structural and AHe thermochronologic evidence, Geol. Soc. Am. Bull., 127, 200–210, https://doi.org/10.1130/B31033.1, 2015. a, b, c

Wong, M. S., Gans, P. B., and Roesler, D.: Field calibration of ^{40}Ar/^{39}Ar K-feldspar multiple diffusion domain (MDD) thermal histories at the Grayback normal fault block, Arizona, USA, Geology, 51, 885–889, https://doi.org/10.1130/G51319.1, 2023. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac, ad, ae, af, ag

Zeitler, P. K.: Argon diffusion in partially outgassed alkali feldspars: Insights from analysis, Chem. Geol., 65, 167–181, https://doi.org/10.1016/0168-9622(87)90071-6, 1987. a, b, c, d

- Abstract
- Introduction
- The MDD model
- Existing optimization methods
- The MDD Tool Kit approach
- Synthetic data experiments
- Case study: Wong et al. (2023) field validation of the MDD model
- Conclusions
- Appendix A: Supplemental tables
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References

^{40}Ar/

^{39}Ar dataset (Wong, 2023) to showcase its efficacy.

- Abstract
- Introduction
- The MDD model
- Existing optimization methods
- The MDD Tool Kit approach
- Synthetic data experiments
- Case study: Wong et al. (2023) field validation of the MDD model
- Conclusions
- Appendix A: Supplemental tables
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References