Thermochron

ApatiteHe(T=Float64;
    age::Number = T(NaN),
    age_sigma::Number = T(NaN),
    offset::Number = zero(T),
    r::Number, 
    dr::Number = one(T), 
    U238::Number, 
    Th232::Number, 
    Sm147::Number = zero(T), 
    U238_matrix::Number = zero(T), 
    Th232_matrix::Number = zero(T), 
    Sm147_matrix::Number = zero(T), 
    volumeweighting::Symbol = :spherical,
    agesteps::AbstractRange,
)

Construct an ApatiteHe chronometer representing an apatite with a raw helium age of age ± age_sigma [Ma], a radius of r [μm], and uniform U, Th and Sm concentrations specified by U238, Th232, and Sm147 [PPMW]. A present day U-235/U-238 ratio of 1/137.818 is assumed.

Spatial discretization follows a radius step of dr [μm], and temporal discretization follows the age steps specified by the agesteps range, in Ma.

Diffusivity(
    D0::T = 59.98               # [cm^2/sec] Maximum diffusion coefficient
    D0_logsigma::T = log(2)/2   # [unitless] log uncertainty (default = log(2)/2 = a factor of 2 two-sigma)
    Ea::T = 205.94              # [kJ/mol] Activation energy
    Ea_logsigma::T = log(2)/2   # [unitless] log uncertainty (default = log(2)/2 = a factor of 2 two-sigma)
)

A generic diffusivity model, with user-specified D0 and Ea. Default values are appropriate for argon in k-feldspar.

Guenthner2013FC(
    C0 = 6.24534         # Guenthner et al. 2013 re-fit of Yamada et al. 2007 zircon
    C1 = -0.11977        # Guenthner et al. 2013 re-fit of Yamada et al. 2007 zircon
    C2 = -314.93688      # Guenthner et al. 2013 re-fit of Yamada et al. 2007 zircon
    C3 = -14.2868        # Guenthner et al. 2013 re-fit of Yamada et al. 2007 zircon
    alpha = -0.057206897 # Guenthner et al. 2013 re-fit of Yamada et al. 2007 zircon
    l0 = 11.17           # [um] Initial track length
    l0_sigma = 0.051     # [um] Initial track length uncertainty
)

Fanning Curvilinear zircon annealing model with simplified Box-Cox transform from Guenthner et al. 2013 (doi: 10.2475/03.2013.01)

Jones2021FA(
    C0 = 1.374           # Annealing parameter
    C1 = -4.192812e-5    # Annealing parameter
    C2 = -22.70885029    # Annealing parameter
    C3 = 0.0             # Annealing parameter
    l0 = 10.60           # [um] Initial track length
    l0_sigma = 0.19      # [um] Initial track length uncertainty
)

Parallel Arrhenius monazite annealing model modified from Jones et al., 2021 (doi: 10.5194/gchron-3-89-2021)

Ketcham1999FC(
    C0::T = -19.844     # "Simultaneous fit" from Ketcham et al. 1999 apatite
    C1::T = 0.38951     # "Simultaneous fit" from Ketcham et al. 1999 apatite
    C2::T = -51.253     # "Simultaneous fit" from Ketcham et al. 1999 apatite
    C3::T = -7.6423     # "Simultaneous fit" from Ketcham et al. 1999 apatite
    alpha::T = -0.12327 # Box-Cox transform parameter
    beta::T = -11.988   # Box-Cox transform parameter
    l0::T = 16.38       # [um] Initial track length
    l0_sigma::T = 0.09  # [um] Initial track length unertainty
)

Fanning Curvilinear apatite annealing model from Ketcham, 1999 (doi: 10.2138/am-1999-0903)

Ketcham2007FC(
    C0 = 0.39528     # "Simultaneous fit" from Ketcham et al. 2007 apatite
    C1 = 0.01073     # "Simultaneous fit" from Ketcham et al. 2007 apatite
    C2 = -65.12969   # "Simultaneous fit" from Ketcham et al. 2007 apatite
    C3 = -7.91715    # "Simultaneous fit" from Ketcham et al. 2007 apatite
    alpha = 0.04672  # "Simultaneous fit" from Ketcham et al. 2007 apatite
    l0 = 16.38       # [um] Initial track length
    l0_sigma = 0.09  # [um] Initial track length unertainty
)

Fanning Curvilinear apatite annealing model with simplified Box-Cox transform from Ketcham, 2007 (doi: 10.2138/am.2007.2281)

MDDiffusivity(
    D0::NTuple{N,T}             # [cm^2/sec] Maximum diffusivity
    D0_logsigma::NTuple{N,T}    # [unitless] log uncertainty (default = 1/2 = a factor of ℯ two-sigma)
    Ea::T                       # [kJ/mol] Activation energy
    Ea_logsigma::T              # [unitless] log uncertainty (default = 1/2 = a factor of ℯ two-sigma)
)

Multiple diffusivities for multiple domains

MultipleDomain(T=Float64, C=PlanarAr;
    age::AbstractVector,
    age_sigma::AbstractVector,
    fraction_experimental::AbstractVector,
    fraction_experimental_sigma::Number=T(0.01),
    tsteps_experimental::AbstractVector,
    Tsteps_experimental::AbstractVector,
    fit::AbstractVector,
    offset::Number = zero(T),
    fuse::Bool = true,
    volume_fraction::AbstractVector,
    r::Number = 100,
    dr::Number = one(T),
    K40::Number = 16.34, 
    agesteps::AbstractRange,
    tsteps::AbstractRange=reverse(agesteps),
)

Construct a MultipleDomain diffusion chronometer given an observed argon release spectrum, degassing schedule, where domain is represented by a PlanarAr or SphericalAr chronometer.

Domain diffusivity and volume parameters must be supplied as vectors Ea [kJ/mol], lnD0a2 [log(1/s)], and volume_fraction [unitless] obtained by separately fitting the release spectrum (the former two as an MDDiffusivity object).

See also: MDDiffusivity, PlanarAr, SphericalAr, degas!

PlanarAr(T=Float64;
    age::Number = T(NaN),
    age_sigma::Number = T(NaN),
    offset::Number = zero(T),
    r::Number, 
    dr::Number = one(T), 
    K40::Number=16.34, 
    agesteps::AbstractRange,
)

Construct an PlanarAr chronometer representing a mineral with a raw argon age of age ± age_sigma [Ma], a uniform diffusivity, a radius of r [μm], and uniform K-40 concentrations specified by K40 [PPM].

Spatial discretization follows a halfwidth step of dr [μm], and temporal discretization follows the age steps specified by the agesteps range, in Ma.

PlanarHe(T=Float64;
    age::Number = T(NaN),
    age_sigma::Number = T(NaN),
    offset::Number = zero(T),
    stoppingpower::Number = T(1.189),
    r::Number, 
    dr::Number = one(T), 
    U238::Number, 
    Th232::Number, 
    Sm147::Number = zero(T), 
    U238_matrix::Number = zero(T), 
    Th232_matrix::Number = zero(T), 
    Sm147_matrix::Number = zero(T), 
    agesteps::AbstractRange,
)

Construct an PlanarHe chronometer representing a mineral with a raw helium age of age ± age_sigma [Ma], uniform diffusivity, a halfwidth of r [μm], and uniform U, Th and Sm concentrations specified by U238, Th232, and Sm147 [PPM]. (A present day U-235/U-238 ratio of 1/137.818 is assumed)

Spatial discretization follows a halfwidth step of dr [μm], and temporal discretization follows the age steps specified by the agesteps range, in Ma.

RDAAM(
    D0L::T=0.6071               # [cm^2/sec] Maximum diffusivity
    D0L_logsigma::T=log(2)/2    # [unitless] log uncertainty (default = log(2)/2 = a factor of 2 two-sigma)
    EaL::T=122.3                # [kJ/mol] Activation energy
    EaL_logsigma::T=log(2)/4    # [unitless] log uncertainty (default = log(2)/4 = a factor of sqrt(2) two-sigma)
    EaTrap::T=34.0              # [kJ/mol] Activation energy
    EaTrap_logsigma::T=log(2)/4 # [unitless] log uncertainty (default = log(2)/4 = a factor of sqrt(2) two-sigma)
    psi::T=1e-13                # empirical polynomial coefficient
    omega::T=1e-22              # empirical polynomial coefficient
    etaq::T=0.91                # Durango ηq
    rhoap::T=3.19               # Density of apatite [g/cm3]
    L::T=0.000815               # Etchable fission track half-length, cm
    lambdaf::T=8.46e-17         # 
    lambdaD::T=1.55125e-10      # 
    beta::T=0.04672             # Apatite annealing parameter. Also caled alpha, but equivalent to beta in ZRDAAM
    C0::T=0.39528               # Apatite annealing parameter
    C1::T=0.01073               # Apatite annealing parameter
    C2::T=-65.12969 - LOG_SEC_MYR # Apatite annealing parameter. Includes conversion factor from seconds to Myr for dt, in addition to traditional C2 value
    C3::T=-7.91715              # Apatite annealing parameter
    rmr0::T=0.83                # Damage conversion parameter
    rmr0_sigma::T=0.15          # Damage conversion parameter uncertainty
    kappa::T=1.04-0.83          # Damage conversion parameter
    kappa_rmr0::T=1.04          # Damage conversion parameter (the sum of kappa and rmr0)
)

Apatite Radiation Damage Accumulation and Annealing Model (RDAAM) of Flowers et al. 2009 (doi: 10.1016/j.gca.2009.01.015)

SphericalAr(T=Float64;
    age::Number = T(NaN),
    age_sigma::Number = T(NaN),
    offset::Number = zero(T),
    r::Number, 
    dr::Number = one(T), 
    K40::Number=16.34, 
    agesteps::AbstractRange,
)

Construct an SphericalAr chronometer representing a mineral with a raw argon age of age ± age_sigma [Ma], a uniform diffusivity, a radius of r [μm], and uniform K-40 concentrations specified by K40 [PPM].

Spatial discretization follows a radius step of dr [μm], and temporal discretization follows the age steps specified by the agesteps range, in Ma.

SphericalHe(T=Float64;
    age::Number = T(NaN),
    age_sigma::Number = T(NaN),
    offset::Number = zero(T),
    stoppingpower::Number = T(1.189),
    r::Number, 
    dr::Number = one(T), 
    U238::Number, 
    Th232::Number, 
    Sm147::Number = zero(T), 
    U238_matrix::Number = zero(T), 
    Th232_matrix::Number = zero(T), 
    Sm147_matrix::Number = zero(T), 
    agesteps::AbstractRange,
)

Construct a SphericalHe chronometer representing a mineral with a raw helium age of age ± age_sigma [Ma], uniform diffusivity, a radius of r [μm], and uniform U, Th and Sm concentrations specified by U238, Th232, and Sm147 [PPM]. (A present day U-235/U-238 ratio of 1/137.818 is assumed)

Spatial discretization follows a radius step of dr [μm], and temporal discretization follows the age steps specified by the agesteps range, in Ma.

Yamada2007PC(
    c0p = -63.37     # Yamada et al. 2007 zircon
    c1p = 0.212      # Yamada et al. 2007 zircon
    bp = 43.00       # Yamada et al. 2007 zircon
    l0 = 11.17       # [um] effective initial track length (μmax)
    l0_sigma = 0.051 # [um] effective initial track length uncertainty (σ)
)

Parallel Curvilinear zircon annealing model of Yamada, 2007 (doi: 10.1016/j.chemgeo.2006.09.002)

ZRDAAM(
    DzD0::T = 193188.0          # [cm^2/sec] Maximum diffusivity, crystalline endmember
    DzD0_logsigma::T=log(2)/2   # [unitless] log uncertainty (default = log(2)/2 = a factor of 2 two-sigma)
    DzEa::T=165.0               # [kJ/mol] Activation energy, crystalline endmember
    DzEa_logsigma::T=log(2)/4   # [unitless] log uncertainty (default = log(2)/4 = a factor of sqrt(2) two-sigma)
    DN17D0::T = 6.367E-3        # [cm^2/sec] Maximum diffusivity, amorphous endmember
    DN17D0_logsigma::T=log(2)/2 # [unitless] log uncertainty (default = log(2)/2 = a factor of 2 two-sigma)
    DN17Ea::T=71.0              # [kJ/mol] Activation energy, amorphous endmember
    DN17Ea_logsigma::T=log(2)/4 # [unitless] log uncertainty (default = log(2)/4 = a factor of sqrt(2) two-sigma)
    lint0::T=45920.0            # [nm]
    SV::T=1.669                 # [1/nm]
    Bα::T=5.48E-19              # Amorphous material produced per alpha decay [g/alpha]
    Phi::T=3.0                  # unitless
    beta::T=-0.05721            # Zircon anealing parameter
    C0::T=6.24534               # Zircon anealing parameter
    C1::T=-0.11977              # Zircon anealing parameter
    C2::T=-314.937 - LOG_SEC_MYR # Zircon anealing parameter. Includes conversion factor from seconds to Myr for dt (for performance), in addition to traditional C2 value
    C3::T=-14.2868              # Zircon anealing parameter
    rmin::T=0.2                 # Damage conversion parameter
    rmin_sigma::T=0.15          # Damage conversion parameter uncertainty
)

Zircon Radiation Damage Accumulation and Annealing Model (ZRDAAM) of Guenthner et al. 2013 (doi: 10.2475/03.2013.01)

ZirconHe(T=Float64;
    age::Number = T(NaN),
    age_sigma::Number = T(NaN),
    T(offset),
    r::Number = one(T),
    dr::Number, 
    U238::Number,
    Th232::Number,
    Sm147::Number = zero(T),
    U238_matrix::Number = zero(T), 
    Th232_matrix::Number = zero(T), 
    Sm147_matrix::Number = zero(T), 
    volumeweighting::Symbol=:cylindrical,
    agesteps::AbstractRange,
)

Construct a ZirconHe chronometer representing a zircon with a raw helium age of age ± age_sigma [Ma], a radius of r [μm], and uniform U, Th and Sm concentrations specified by U238, Th232, and Sm147 [PPMw]. A present day U-235/U-238 ratio of 1/137.818 is assumed.

Spatial discretization follows a radius step of dr μm, and temporal discretization follows the age steps specified by the agesteps range, in Ma.

image_from_paths(tT::TTResult; 
    xresolution::Int=1800, 
    yresolution::Int=1200, 
    xrange = nanextrema(tT.tpointdist), 
    yrange = nanextrema(tT.Tpointdist), 
)

Produce a 2d image (histogram) of path densities given a TTResult

julia> imgcounts, xq, yq = image_from_paths(tT; xrange=boundary.agepoints, yrange=boundary.T₀)

MCMC(chrons::Vector{<:Chronometer}, model::NamedTuple, npoints::Int, path.agepoints::Vector, path.Tpoints::Vector, constraint::Constraint, boundary::Boundary, [detail::DetailInterval])

Markov chain Monte Carlo time-Temperature inversion of the thermochronometric chrons specified as a vector chrons of Chronometer objects (ZirconHe, ApatiteHe, ApatiteFT, etc.) and model parameters specified by the named tuple model, with variable diffusion kinetics.

Returns a TTResult object containing posterior time-temperature paths,

See also MCMC_varkinetics for a variant with variable diffusion kinetics.

Examples

tT = MCMC(chrons::NamedTuple, model::NamedTuple, constraint::Constraint, boundary::Boundary, [detail::DetailInterval])

MCMC_varkinetics(chrons::Vector{<:Chronometer}, model::NamedTuple, npoints::Int, path.agepoints::Vector, path.Tpoints::Vector, constraint::Constraint, boundary::Boundary, [detail::DetailInterval])

Markov chain Monte Carlo time-Temperature inversion of the thermochronometric chrons specified as a vector chrons of Chronometer objects (ZirconHe, ApatiteHe, ApatiteFT, etc.) and model parameters specified by the named tuple model, with variable diffusion kinetics.

Returns a TTResult object containing posterior time-temperature paths, and a KineticResult object containing the posterior kinetic parameters.

See also MCMC for a variant with constant diffusion kinetics.

Examples

tT, kinetics = MCMC_varkinetics(chrons::NamedTuple, model::NamedTuple, constraint::Constraint, boundary::Boundary, [detail::DetailInterval])

ρᵣ = anneal(dt::Number, tsteps::Vector, Tsteps::Matrix, [model::DiffusivityModel=ZRDAAM()])

ZirconHe damage annealing model as in Guenthner et al. 2013 (AJS)

anneal!(data::Vector{<:Chronometer}, ::Type{<:HeliumSample}, tsteps, Tsteps, dm::DiffusivityModel)
anneal!(mineral::ZirconHe, Tsteps, dm::ZRDAAM)
anneal!(mineral::ApatiteHe, Tsteps, dm::RDAAM)
anneal!(ρᵣ::Matrix, dt::Number, tsteps, Tsteps, [dm::DiffusivityModel=ZRDAAM()])

In-place version of anneal

chronometers([T=Float64], data, model)

Construct a vector of Chronometer objects given a dataset data and model parameters model

lcmod(l, θ)

Calculate the model c-axis equivalent length ("lc,mod") given a measured "confined" fission track length l [microns] and angle from the c-axis θ [degrees] following the approach of Donelick et al. 1999 (doi: 10.2138/am-1999-0902)

modelage(mineral::ZirconHe, Tsteps, [ρᵣ], dm::ZRDAAM)
modelage(mineral::ApatiteHe, Tsteps, [ρᵣ], dm::RDAAM)
modelage(mineral::SphericalHe, Tsteps, dm::Diffusivity)
modelage(mineral::PlanarHe, Tsteps, dm::Diffusivity)

Calculate the predicted bulk U-Th/He age of a zircon, apatite, or other mineral that has experienced a given t-T path (specified by mineral.tsteps for time and Tsteps for temperature, at a time resolution of step(mineral.tsteps)) using a Crank-Nicolson diffusion solution for a spherical (or planar slab) grain of radius (or halfwidth) mineral.r at spatial resolution mineral.dr.

Spherical implementation based on the the Crank-Nicolson solution for diffusion out of a spherical mineral crystal in Ketcham, 2005 (doi: 10.2138/rmg.2005.58.11).

modelage(mineral::SphericalAr, Tsteps, dm::Diffusivity)
modelage(mineral::PlanarAr, Tsteps, dm::Diffusivity)

Calculate the precdicted bulk K/Ar age of a mineral that has experienced a given t-T path (specified by mineral.tsteps for time and Tsteps for temperature, at a time resolution of step(mineral.tsteps)) using a Crank-Nicholson diffusion solution for a spherical (or planar slab) grain of radius (or halfwidth ) mineral.r at spatial resolution mineral.dr.

Spherical implementation based on the the Crank-Nicolson solution for diffusion out of a spherical mineral crystal in Ketcham, 2005 (doi: 10.2138/rmg.2005.58.11).

modelage(mineral::ZirconFT, Tsteps, am::ZirconAnnealingModel)
modelage(mineral::MonaziteFT, Tsteps, am::MonaziteAnnealingModel)
modelage(mineral::ApatiteFT, Tsteps, am::ApatiteAnnealingModel)

Calculate the precdicted fission track age of an apatite that has experienced a given t-T path (specified by mineral.tsteps for time and Tsteps for temperature, at a time resolution of step(mineral.tsteps)) and given annealing model parameters am.

Possible annealing model types and the references for the equations which they respetively implement include Ketcham1999FC Fanning Curvilinear apatite model of Ketcham et al. 1999 (doi: 10.2138/am-1999-0903) Ketcham2007FC Fanning Curvilinear apatite model of Ketcham et al. 2007 (doi: 10.2138/am.2007.2281) Yamada2007PC Parallel Curvilinear zircon model of Yamada et al. 2007 (doi: 10.1016/j.chemgeo.2006.09.002) Guenthner2013FC Fanning Curvilinear zircon model of Guenthner et al. 2013 (doi: 10.2475/03.2013.01) Jones2021FA Fanning Arrhenius (Fanning Linear) model adapted from Jones et al. 2021 (doi: 10.5194/gchron-3-89-2021)

modellength(track::ApatiteTrackLength, Tsteps, am::ApatiteAnnealingModel)

Calculate the predicted mean and standard deviation of the distribution of fission track lengths of an apatite that has experienced a given t-T path (specified by track.tsteps for time and Tsteps for temperature, at a time resolution of step(mineral.tsteps)) and given annealing model parameters am.

Possible annealing model types and the references for the equations which they respetively implement include Ketcham1999FC Fanning Curvilinear apatite model of Ketcham et al. 1999 (doi: 10.2138/am-1999-0903) Ketcham2007FC Fanning Curvilinear apatite model of Ketcham et al. 2007 (doi: 10.2138/am.2007.2281)

reltrackdensity(t, T, am::AnnealingModel)

Calculate the relative track density ρ corresponding to a given relative track length r

Follows the relations of Ketcham et al. (2000), equations 7a and 7b for apatite and Tagami et al. (1999) for zircon

See also: reltracklength.

rmr0fromdpar(dpar)

Calculate rmr0 as a function of dpar for "multikinetic" apatite fission track following the relation (Fig. 7) of Ketcham et al. 1999 (doi: 10.2138/am-1999-0903)

rmr0 = 1 - exp(0.647(dpar-1.75) - 1.834)

rmr0model(F, Cl, OH, Mn=0, Fe=0, others=0)

Calculate rmr0 as a function of composition (specified in terms of atoms per fomula unit, or APFU) for "multikinetic" apatite fission track thermochronology.

Implements equation 11 from Ketcham et al. 2007 (doi: 10.2138/am.2007.2281)

rmr0 = (-0.0495 -0.0348F +0.3528|Cl - 1| +0.0701|OH - 1| 
        -0.8592Mn -1.2252Fe -0.1721Others)^0.1433

simannealT(n::Integer, Tₐ::Number, λₐ::Number)

To avoid getting stuck in local optima, increase the probability of accepting new proposals at higher annealing "temperature"

slabsphereintersectiondensity(redges::Vector, ralpha, d)

Calculate the fractional intersection density of an alpha stopping sphere of radius ralpha with each concentric slab-shell (with edges redges and relative volumes rvolumes) of a planar slab crystal where the two are separated by distance d

slabsphereintersectionfraction(rₚ, rₛ, d)

Calculate the fraction of the surface area of a sphere s with radius rₛ that intersects the interior of a planar slab p of halfwidth rₚ if the two are separated by distance d.

sphereintersectiondensity(redges::Vector, rvolumes::Vector, ralpha, d)

Calculate the volume-nomalized fractional intersection density of an alpha stopping sphere of radius ralpha with each concentric shell (with shell edges redges and relative volumes rvolumes) of a spherical crystal where the two are separated by distance d

sphereintersectionfraction(r₁, r₂, d)

Calculate the fraction of the surface area of a sphere s₂ with radius r₂ that intersects the interior of a sphere s₁ of radius r₁ if the two are separated by distance d.