Thermochron

ApatiteFT(T::Type{<:AbstractFloat} = Float64; 
    age::Number = NaN,              # [Ma] fission track age
    age_sigma::Number = NaN,        # [Ma] fission track age uncertainty
    offset::Number = zero(T),       # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),       # [m] height relative to other samples / surface (negative for depth)
    dpar::Number = NaN,             # [um] diameter parallel to track
    F::Number = NaN,                # [APFU] F concentration, in atoms per formula unit
    Cl::Number = NaN,               # [APFU] Cl concentration, in atoms per formula unit
    OH::Number = NaN,               # [APFU] OH concentration, in atoms per formula unit
    rmr0::Number = NaN,             # [unitless] annealing parameter
    name::String = "",              # Sample or grain name
    notes::String = "",             # Sample notes
    agesteps::AbstracVector | tsteps::AbstracVector, # Temporal discretization
)

Construct an ApatiteFT chronometer representing a apatite with a fission track age of age ± age_sigma [Ma] and a relative annealing resistance specified by rmr0, and optionally a constant temperature offset (relative to other samples) of offset [C].

If not provided directly, rmr0 will be calculated, in order of preference:

1. from F, Cl, and OH together, via the apatite_rmr0model function
2. from Cl alone, via the apatite_rmr0fromcl function
3. from dpar, via the apatite_rmr0fromdpar functions
4. using a default fallback value of 0.83, if none of the above are provided.

Temporal discretization follows the age steps specified by agesteps (age before present) and/or tsteps (forward time since crystallization), in Ma, where tsteps must be sorted in increasing order.

ApatiteHe(T=Float64;
    age::Number = T(NaN),                   # [Ma] raw helium age
    age_sigma::Number = T(NaN),             # [Ma] raw helium age uncertainty (one-sigma)
    offset::Number = zero(T),               # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),               # [m] height relative to other samples / surface (negative for depth)
    r::Number,                              # [um] spherical radius 
    dr::Number = one(T),                    # [um] radial step size
    U238::Number,                           # [ppm] apatite U-238 concentration
    Th232::Number,                          # [ppm] apatite Th-232 concentration
    Sm147::Number = zero(T),                # [ppm] apatite Sm-147 concentration
    U238_matrix::Number = zero(T),          # [ppm] matrix U-238 concentration
    Th232_matrix::Number = zero(T),         # [ppm] matrix Th-232 concentration
    Sm147_matrix::Number = zero(T),         # [ppm] matrix Sm-147 concentration
    grainsize_matrix::Number = one(T),      # [mm] average grain size of matrix rock
    volumeweighting::Symbol=:cylindrical,   # (:spherical, :cylindrical, or :planar) relative volume proportions of each radial model shell, for averaging purposes
    name::String = "",                      # Sample or grain name
    notes::String = "",                     # Sample notes
    agesteps::AbstractVector | tsteps::AbstractVector, # Temporal discretization
)

Construct an ApatiteHe chronometer representing an apatite with a raw helium age of age ± age_sigma [Ma], a spherical 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 halfwidth step of dr [μm], while temporal discretization follows the age steps specified by agesteps (age before present) and/or tsteps (forward time since crystallization), in Ma, where tsteps must be sorted in increasing order.

ApatiteTrackLength(T::Type{<:AbstractFloat}=Float64; 
    length::Number = NaN,                   # [um] fission track length
    offset::Number = zero(T),               # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),               # [m] height relative to other samples / surface (negative for depth)
    l0::Number = NaN,                       # [um] Initial track length
    l0_sigma::Number = NaN,                 # [um] Initial track length unertainty
    dpar::Number = NaN,                     # [um] diameter parallel to track
    F::Number = NaN,                        # [APFU] F concentration, in atoms per formula unit
    Cl::Number = NaN,                       # [APFU] Cl concentration, in atoms per formula unit
    OH::Number = NaN,                       # [APFU] OH concentration, in atoms per formula unit
    rmr0::Number = NaN,                     # [unitless] annealing parameter
    name::String = "",                      # Sample or grain name
    notes::String = "",                     # Sample notes
    agesteps::AbstracVector | tsteps::AbstracVector, # Temporal discretization
)

Construct an ApatiteTrackLengthOriented chronometer representing a single apatite fission track length um long with a relative annealing resistance specified by rmr0, optionally at a constant temperature offset (relative to other samples) of offset [C].

If not provided directly, rmr0 will be calculated, in order of preference:

1. from F, Cl, and OH together, via the apatite_rmr0model function
2. from Cl alone, via the apatite_rmr0fromcl function
3. from dpar, via the apatite_rmr0fromdpar functions
4. using a default fallback value of 0.83, if none of the above are provided.

If not provided directly, l0 and l0_sigma will be estimated using the apatite_l0fromdpar function, where dpar in turn is estimated using the apatite_dparfromrmr0 function if not provided directly.

Temporal discretization follows the age steps specified by agesteps (age before present) and/or tsteps (forward time since crystallization), in Ma, where tsteps must be sorted in increasing order.

ApatiteTrackLengthOriented(T::Type{<:AbstractFloat}=Float64; 
    length::Number = NaN,                   # [um] fission track length
    angle::Number = NaN,                    # [degrees] track angle from the c-axis
    lcmod::Number = lcmod(length, angle),   # [um] model length of an equivalent c-axis parallel rack
    offset::Number = zero(T),               # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),               # [m] height relative to other samples / surface (negative for depth)
    l0::Number = NaN,                       # [um] Initial track length
    l0_sigma::Number = NaN,                 # [um] Initial track length unertainty
    dpar::Number = NaN,                     # [um] diameter parallel to track
    F::Number = NaN,                        # [APFU] F concentration, in atoms per formula unit
    Cl::Number = NaN,                       # [APFU] Cl concentration, in atoms per formula unit
    OH::Number = NaN,                       # [APFU] OH concentration, in atoms per formula unit
    rmr0::Number = NaN,                     # [unitless] annealing parameter
    name::String = "",                      # Sample or grain name
    notes::String = "",                     # Sample notes
    agesteps::AbstracVector | tsteps::AbstracVector, # Temporal discretization
)

Construct an ApatiteTrackLengthOriented chronometer representing a single apatite fission track length um long, oriented at angle degrees to the c-axis, with a relative annealing resistance specified by rmr0, optionally at a constant temperature offset (relative to other samples) of offset [C].

If not provided directly, rmr0 will be calculated, in order of preference:

1. from F, Cl, and OH together, via the apatite_rmr0model function
2. from Cl alone, via the apatite_rmr0fromcl function
3. from dpar, via the apatite_rmr0fromdpar functions
4. using a default fallback value of 0.83, if none of the above are provided.

If not provided directly, l0 and l0_sigma will be estimated using the apatite_l0fromdpar function, where dpar in turn is estimated using the apatite_dparfromrmr0 function if not provided directly.

Temporal discretization follows the age steps specified by agesteps (age before present) and/or tsteps (forward time since crystallization), in Ma, where tsteps must be sorted in increasing order.

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
)

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
)

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
)

C-axis projected 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
)

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

MDiffusivity(
    D0::NTuple{N,T}                 # [cm^2/sec] Maximum diffusivity
    D0_logsigma::NTuple{N,T}        # [unitless] log uncertainty
    Ea::NTuple{N,T}                 # [kJ/mol] Activation energy
    Ea_logsigma::NTuple{N,T}        # [unitless] log uncertainty
    volume_fraction::NTuple{N,T}    # [unitless] Volume fraction of each domain
)

Multiple diffusivities for multiple domains

MSDiffusivity(
    model::DiffusivityModel{T}      # Underlying (wrapped) diffusivity model
    scale::NTuple{N,T}              # [unitless] relative domain size 
    scale_logsigma::NTuple{N,T}     # [unitless] log uncertainty 
    volume_fraction::NTuple{N,T}    # [unitless] Volume fraction of each domain
)

One diffusivity, rescaled multiple times

MonaziteFT(T::Type{<:AbstractFloat} = Float64; 
    age::Number = NaN,              # [Ma] fission track age
    age_sigma::Number = NaN,        # [Ma] fission track age uncertainty
    offset::Number = zero(T),       # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),       # [m] height relative to other samples / surface (negative for depth)
    name::String = "",              # Sample or grain name
    notes::String = "",             # Sample notes
    agesteps::AbstracVector | tsteps::AbstracVector, # Temporal discretization
)

Construct a MonaziteFT chronometer representing a monazite with a fission track age of age ± age_sigma [Ma], optionally at a constant temperature offset (relative to other samples) of offset [C].

Temporal discretization follows the age steps specified by agesteps (age before present) and/or tsteps (forward time since crystallization), in Ma, where tsteps must be sorted in increasing order.

MonaziteTrackLength(T::Type{<:AbstractFloat} = Float64; 
    length::Number = NaN,                   # [um] fission track length
    offset::Number = zero(T),               # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),               # [m] height relative to other samples / surface (negative for depth)
    l0::Number = 10.60,                     # [um] Initial track length
    l0_sigma::Number = 0.19,                # [um] Initial track length unertainty
    name::String = "",                      # Sample or grain name
    notes::String = "",                     # Sample notes
    agesteps::AbstracVector | tsteps::AbstracVector, # Temporal discretization
)

Construct a MonaziteTrackLength chronometer representing a single monazite fission track length um long, optionally at a constant temperature offset (relative to other samples) of offset [C].

Temporal discretization follows the age steps specified by agesteps (age before present) and/or tsteps (forward time since crystallization), in Ma, where tsteps must be sorted in increasing order.

MultipleDomain(T=Float64, C=PlanarAr;
    step_age::AbstractVector,                       # [Ma] measured Ar-40/Ar-39 ages at each degassing step
    step_age_sigma::AbstractVector,                 # [Ma] measured Ar-40/Ar-39 age uncertainties (one-sigma) at each degassing step
    fraction_experimental::AbstractVector,          # [unitless] cumulative fraction of total Ar-39 released each degassing step
    fraction_experimental_sigma=fill(T(0.005), size(fraction_experimental)),     # [unitless] uncertainty in degassing fraction
    tsteps_experimental::AbstractVector,            # [s] time steps of experimental heating schedule
    Tsteps_experimental::AbstractVector,            # [C] temperature steps of experimental heating schedule
    fit::AbstractVector,                            # [Bool] Whether or not each degassing step should be used in inversion
    offset::Number = zero(T),                       # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),                       # [m] height relative to other samples / surface (negative for depth)
    fuse::Bool = true,                              # [Bool] Treat the grain as having fused (released all remaining Ar)
    volume_fraction::AbstractVector,                # [unitless] fraction of total volume represented by each domain
    r::Number = 100,                                # [um] model domain radius (spherical) or half-width (planar)
    dr::Number = one(T),                            # [um] model domain radius step
    name::String = "",                              # Sample or grain name
    notes::String = "",                             # Sample notes
    agesteps::AbstractVector | tsteps::AbstractVector, # Temporal discretization
    kwargs...                                       # Any further keyword arguments are forwarded to the constructor for the domain `C`
)

Construct a MultipleDomain diffusion chronometer given an observed release spectrum and degassing schedule, where each domain is in turn represented by a Chronometer object (e.g. PlanarAr, SphericalAr).

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

See also: MDiffusivity, PlanarAr, SphericalAr, degas!

PlanarAr(T=Float64;
    age::Number = T(NaN),                   # [Ma] Ar-40/Ar-39 age
    age_sigma::Number = T(NaN),             # [Ma] Ar-40/Ar-39 age uncertainty (one-sigma)
    offset::Number = zero(T),               # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),               # [m] height relative to other samples / surface (negative for depth)
    r::Number,                              # [um] planar half-width
    dr::Number = one(T),                    # [um] radial step size
    K40::Number = 16.34,                    # [ppm] mineral K-40 concentration
    K40_matrix::Number = zero(T)            # [ppm] matrix K-40 concentration
    grainsize_matrix::Number = one(T),      # [mm] average grain size of matrix rock
    name::String = "",                      # Sample or grain name
    notes::String = "",                     # Sample notes
    agesteps::AbstractVector | tsteps::AbstractVector, # Temporal discretization
)

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], while temporal discretization follows the age steps specified by agesteps (age before present) and/or tsteps (forward time since crystallization), in Ma, where tsteps must be sorted in increasing order.

PlanarHe(T=Float64;
    age::Number = T(NaN),                   # [Ma] raw helium age
    age_sigma::Number = T(NaN),             # [Ma] raw helium age uncertainty (one-sigma)
    offset::Number = zero(T),               # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),               # [m] height relative to other samples / surface (negative for depth)
    stoppingpower::Number = T(1.189),       # [unitless] alpha stopping power relative to apatite
    r::Number,                              # [um] planar half-width
    dr::Number = one(T),                    # [um] radial step size
    U238::Number,                           # [ppm] mineral U-238 concentration
    Th232::Number,                          # [ppm] mineral Th-232 concentration
    Sm147::Number = zero(T),                # [ppm] mineral Sm-147 concentration
    U238_matrix::Number = zero(T),          # [ppm] matrix U-238 concentration
    Th232_matrix::Number = zero(T),         # [ppm] matrix Th-232 concentration
    Sm147_matrix::Number = zero(T),         # [ppm] matrix Sm-147 concentration
    grainsize_matrix::Number = one(T),      # [mm] average grain size of matrix rock
    name::String = "",                      # Sample or grain name
    notes::String = "",                     # Sample notes
    agesteps::AbstractVector | tsteps::AbstractVector, # Temporal discretization
)

Construct an PlanarHe chronometer representing a mineral with a raw helium age of age ± age_sigma [Ma], uniform diffusivity, a planar half-width 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 halfwidth step of dr [μm], while temporal discretization follows the age steps specified by agesteps (age before present) and/or tsteps (forward time since crystallization), in Ma, where tsteps must be sorted in increasing order.

RDAAM(
    D0_L::T=0.6071               # [cm^2/sec] Maximum diffusivity
    D0_L_logsigma::T=log(2)/2    # [unitless] log uncertainty (default = log(2)/2 = a factor of 2 two-sigma)
    Ea_L::T=122.3                # [kJ/mol] Activation energy
    Ea_L_sigma::T=12.23          # [kJ/mol] uncertainty
    Ea_trap::T=34.0              # [kJ/mol] Activation energy
    Ea_trap_sigma::T=3.4         # [kJ/mol] uncertainty
    psi::T=1e-13                # empirical polynomial coefficient
    omega::T=1e-22              # empirical polynomial coefficient
    etaq::T=0.91                # Durango ηq
    rhoap::T=3.19               # [g/cm3] Density of apatite 
    L::T=0.000815               # [cm] Etchable fission track half-length
    lambdaf::T=8.46e-17         # [1/a] U-238 spontaneous fission decay constant
    lambdaD::T=1.55125e-10      # [1/a] U-238 total decay constant
    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)

SDiffusivity(
    model::DiffusivityModel{T}      # Underlying (wrapped) diffusivity model
    scale::T                        # [unitless] relative domain size (default = 1.0)
    scale_logsigma::T               # [unitless] log uncertainty (default = 1.0 = a factor of ℯ, one-sigma)
)

One diffusivity, rescaled

SingleDomain(T=Float64, C=ApatiteHe;
    step_age::AbstractVector,                       # [Ma] measured Ar-40/Ar-39 ages at each degassing step
    step_age_sigma::AbstractVector,                 # [Ma] measured Ar-40/Ar-39 age uncertainties (one-sigma) at each degassing step
    fraction_experimental::AbstractVector,          # [unitless] cumulative fraction of total Ar-39 released each degassing step
    fraction_experimental_sigma=fill(T(0.005), size(fraction_experimental)),     # [unitless] uncertainty in degassing fraction
    tsteps_experimental::AbstractVector,            # [s] time steps of experimental heating schedule
    Tsteps_experimental::AbstractVector,            # [C] temperature steps of experimental heating schedule
    fit::AbstractVector,                            # [Bool] Whether or not each degassing step should be used in inversion
    offset::Number = zero(T),                       # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),                       # [m] height relative to other samples / surface (negative for depth)
    fuse::Bool = true,                              # [Bool] Treat the grain as having fused (released all remaining Ar)
    name::String = "",                              # Sample or grain name
    notes::String = "",                             # Sample notes
    agesteps::AbstractVector | tsteps::AbstractVector, # Temporal discretization
    kwargs...                                       # Any further keyword arguments are forwarded to the constructor for the domain `C`
)

Construct a SingleDomain diffusion chronometer given an observed experimental release spectrum and degassing schedule, where the single modelled domain is represented by a Chronometer object (e.g. ApatiteHe for He-4/He-3).

See also: degas!

SphericalAr(T=Float64;
    age::Number = T(NaN),                   # [Ma] Ar-40/Ar-39 age
    age_sigma::Number = T(NaN),             # [Ma] Ar-40/Ar-39 age uncertainty (one-sigma)
    offset::Number = zero(T),               # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),               # [m] height relative to other samples / surface (negative for depth)
    r::Number,                              # [um] equivalent spherical radius 
    dr::Number = one(T),                    # [um] radial step size
    K40::Number = 16.34,                    # [ppm] mineral K-40 concentration
    K40_matrix::Number = zero(T)            # [ppm] matrix K-40 concentration
    grainsize_matrix::Number = one(T),      # [mm] average grain size of matrix rock
    name::String = "",                      # Sample or grain name
    notes::String = "",                     # Sample notes
    agesteps::AbstractVector | tsteps::AbstractVector, # Temporal discretization
)

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

Spatial discretization follows a halfwidth step of dr [μm], while temporal discretization follows the age steps specified by agesteps (age before present) and/or tsteps (forward time since crystallization), in Ma, where tsteps must be sorted in increasing order.

SphericalHe(T=Float64;
    age::Number = T(NaN),                   # [Ma] raw helium age
    age_sigma::Number = T(NaN),             # [Ma] raw helium age uncertainty (one-sigma)
    offset::Number = zero(T),               # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),               # [m] height relative to other samples / surface (negative for depth)
    stoppingpower::Number = T(1.189),       # [unitless] alpha stopping power relative to apatite
    r::Number,                              # [um] spherical radius 
    dr::Number = one(T),                    # [um] radial step size
    U238::Number,                           # [ppm] mineral U-238 concentration
    Th232::Number,                          # [ppm] mineral Th-232 concentration
    Sm147::Number = zero(T),                # [ppm] mineral Sm-147 concentration
    U238_matrix::Number = zero(T),          # [ppm] matrix U-238 concentration
    Th232_matrix::Number = zero(T),         # [ppm] matrix Th-232 concentration
    Sm147_matrix::Number = zero(T),         # [ppm] matrix Sm-147 concentration
    grainsize_matrix::Number = one(T),      # [mm] average grain size of matrix rock
    name::String = "",                      # Sample or grain name
    notes::String = "",                     # Sample notes
    agesteps::AbstractVector | tsteps::AbstractVector, # Temporal discretization
)

Construct a SphericalHe chronometer representing a mineral with a raw helium age of age ± age_sigma [Ma], uniform diffusivity, a spherical 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 halfwidth step of dr [μm], while temporal discretization follows the age steps specified by agesteps (age before present) and/or tsteps (forward time since crystallization), in Ma, where tsteps must be sorted in increasing order.

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
)

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

ZRDAAM(
    D0_z::T = 193188.0          # [cm^2/sec] Maximum diffusivity, crystalline endmember
    D0_z_logsigma::T=log(2)/2   # [unitless] log uncertainty (default = log(2)/2 = a factor of 2 two-sigma)
    Ea_z::T=165.0               # [kJ/mol] Activation energy, crystalline endmember
    Ea_z_sigma::T=16.5          # [kJ/mol] uncertainty 
    D0_N17::T = 6.367E-3        # [cm^2/sec] Maximum diffusivity, amorphous endmember
    D0_N17_logsigma::T=log(2)/2 # [unitless] log uncertainty (default = log(2)/2 = a factor of 2 two-sigma)
    Ea_N17::T=71.0              # [kJ/mol] Activation energy, amorphous endmember
    Ea_N17_sigma::T=7.1         # [kJ/mol] uncertainty 
    lint0::T=45920.0            # [nm]
    SV::T=1.669                 # [1/nm]
    Ba::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)

ZirconFT(T::Type{<:AbstractFloat} = Float64; 
    age::Number = NaN,              # [Ma] fission track age
    age_sigma::Number = NaN,        # [Ma] fission track age uncertainty
    offset::Number = zero(T),       # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),       # [m] height relative to other samples / surface (negative for depth)
    name::String = "",              # Sample or grain name
    notes::String = "",             # Sample notes
    agesteps::AbstracVector | tsteps::AbstracVector, # Temporal discretization
)

Construct a ZirconFT chronometer representing a zircon with a fission track age of age ± age_sigma [Ma], optionally at a constant temperature offset (relative to other samples) of offset [C].

Temporal discretization follows the age steps specified by agesteps (age before present) and/or tsteps (forward time since crystallization), in Ma, where tsteps must be sorted in increasing order.

ZirconHe(T=Float64;
    age::Number = T(NaN),                   # [Ma] raw helium age
    age_sigma::Number = T(NaN),             # [Ma] raw helium age uncertainty (one-sigma)
    offset::Number = zero(T),               # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),               # [m] height relative to other samples / surface (negative for depth)
    r::Number,                              # [um] spherical radius 
    dr::Number = one(T),                    # [um] radial step size
    U238::Number,                           # [ppm] zircon U-238 concentration
    Th232::Number,                          # [ppm] zircon Th-232 concentration
    Sm147::Number = zero(T),                # [ppm] zircon Sm-147 concentration
    U238_matrix::Number = zero(T),          # [ppm] matrix U-238 concentration
    Th232_matrix::Number = zero(T),         # [ppm] matrix Th-232 concentration
    Sm147_matrix::Number = zero(T),         # [ppm] matrix Sm-147 concentration
    grainsize_matrix::Number = one(T),      # [mm] average grain size of matrix rock
    volumeweighting::Symbol=:cylindrical,   # (:spherical, :cylindrical, or :planar) relative volume proportions of each radial model shell, for averaging purposes
    name::String = "",                      # Sample or grain name
    notes::String = "",                     # Sample notes
    agesteps::AbstractVector | tsteps::AbstractVector, # Temporal discretization
)

Construct a ZirconHe chronometer representing a zircon with a raw helium age of age ± age_sigma [Ma], a spherical 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 halfwidth step of dr [μm], while temporal discretization follows the age steps specified by agesteps (age before present) and/or tsteps (forward time since crystallization), in Ma, where tsteps must be sorted in increasing order.

ZirconTrackLength(T::Type{<:AbstractFloat}=Float64; 
    length::Number = NaN,                   # [um] fission track length
    offset::Number = zero(T),               # [C] temperature offset relative to other samples / surface (positive is warmer)
    height::Number = zero(T),               # [m] height relative to other samples / surface (negative for depth)
    l0::Number = 11.17,                     # [um] Initial track length
    l0_sigma::Number = 0.051,               # [um] Initial track length unertainty
    name::String = "",                      # Sample or grain name
    notes::String = "",                     # Sample notes  
    agesteps::AbstracVector | tsteps::AbstracVector, # Temporal discretization
)

Construct a ZirconTrackLength chronometer representing a single zircon fission track length um long, optionally at a constant temperature offset (relative to other samples) of offset [C].

Temporal discretization follows the age steps specified by agesteps (age before present) and/or tsteps (forward time since crystallization), in Ma, where tsteps must be sorted in increasing order.

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::AbstractVector{<:Chronometer}, params::NamedTuple, npoints::Int, path.agepoints::Vector, path.Tpoints::Vector, constraint::Constraint, boundary::Boundary, [detail::DetailInterval];
    rp::RegionalParameters = RegionalParameters(),
    rng::AbstractRNG = Random.default_rng(),
    liveplot::Bool = false, 
    maxplots::Int = 512, 
    maxplotsburnin::Int = maxplots÷2, 
    maxplotscollection::Int = maxplots÷2
)

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 params, 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, params::NamedTuple, constraint::Constraint, boundary::Boundary, [detail::DetailInterval])

MCMC_varkinetics(chrons::AbstractVector{<:Chronometer}, params::NamedTuple, npoints::Int, path.agepoints::Vector, path.Tpoints::Vector, constraint::Constraint, boundary::Boundary, [detail::DetailInterval];
    rp::RegionalParameters{T} = RegionalParameters{T}(),
    rng::AbstractRNG = Random.default_rng(),
    liveplot::Bool = false, 
    maxplots::Int = 512, 
    maxplotsburnin::Int = maxplots÷2, 
    maxplotscollection::Int = maxplots÷2
)

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 params, 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, params::NamedTuple, constraint::Constraint, boundary::Boundary, [detail::DetailInterval])

alphacorrectionslab(alpharadii, rparent, redges, parentₒ, redgesₒ)

Effective parent concentrations accounting for alphas per decay at secular equilibrium as well as alpha ejection (given internal radial rparent concentration vector and spatial discretization redges) and injection (given uniform external parentₒ concentration and spatial discretization redgesₒ) for an infinite slab mineral grain.

alphacorrectionspherical(alpharadii, rparent, redges, relvolumes, parentₒ, redgesₒ, relvolumesₒ)

Effective parent concentrations accounting for alphas per decay at secular equilibrium as well as alpha ejection (given internal radial rparent concentration vector and spatial discretization redges, relvolumes) and injection (given uniform external parentₒ concentration and spatial discretization redgesₒ, relvolumesₒ) for an infinite slab mineral grain.

alphastoppingpower(mineral)

Alpha particle stopping power relative to that of apatite, as calculated from the mean alpha stopping distances of Ketcham et al. 2011 (doi: 10.1016/j.gca.2011.10.011)

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

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

anneal!(data::AbstractVector{<: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

apatite_dparfromrmr0(rmr₀)

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

dpar = (log(1 - rmr₀) + 1.834)/0.647 + 1.75

l0 = apatite_l0fromdpar(dpar)

Calculate l0 as a function of dpar for "multikinetic" apatite fission track following the relation (equation 1) of Carlson et al. 1999 (doi: 10.2138/am-1999-0901), that is

$l_{0,m} = 15.63 + 0.283*D_{par}$

The results may be used along with a constant uncertainty of 0.1367 based on the scatter around the appropriate line in Figure 1 of Carlson et al. 1999.

l0 = apatite_l0modfromdpar(dpar)

Calculate c-axis-projected l0 as a function of dpar for "multikinetic" apatite fission track following the relation (equation 2) of Carlson et al. 1999 (doi: 10.2138/am-1999-0901), that is

$l_{0,c,mod} = 16.10 + 0.205*D_{par}$

The results may be used along with a constant uncertainty of 0.1311 based on the scatter around the appropriate line in Figure 1 of Carlson et al. 1999.

apatite_rmr0fromcl(Cl)

Calculate rmr₀ as a function of chlorine content Cl [APFU] for "multikinetic" apatite fission track following the relation (Fig. 7a) of Ketcham et al. 1999 (doi: 10.2138/am-1999-0903)

rmr₀ = 1 - exp(2.107(1 - abs(Cl - 1)) - 1.834)

apatite_rmr0fromdpar(dpar)

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

rmr₀ = 1 - exp(0.647(dpar-1.75) - 1.834)

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

Calculate rmr₀ 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)

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

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

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

function empiricaluncertainty!(σcalc, chrons, C::Type{<:HeliumSample};
    fraction::Number = 1/sqrt(2),
    sigma_eU::Number = (C<:ZirconHe) ? 100.0 : (C<:ApatiteHe) ? 10.0 : 25.0,
    sigma_offset::Number = 10.0,
)

Given a vector of chronometers chrons, update the uncertainties in σcalc to reflect the uncertainty implied by the observed (empirical) scatter in the helium ages for samples with similar eU and temperature offset (i.e., elevation).

By default, half the empirical variance (1/sqrt(2) of the standard deviation) is interpreted as external uncertainty which should be reflected in σcalc.

Similarity in eU and offset is assessed on the basis of Gaussian kernels with bandwidth equal to sigma_eU for eU and sigma_offset for temperature offset.

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)
modelage(mineral::SphericalAr, Tsteps, dm::Diffusivity)
modelage(mineral::PlanarAr, Tsteps, dm::Diffusivity)

Calculate the predicted bulk age of a noble gas chronometer that has experienced a given t-T path (specified by tsteps_geol(mineral) for time and Tsteps for temperature), at a time resolution determined by tsteps_geol(mineral) 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::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 determined by mineral.tsteps 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::ApatiteTrackLengthOriented, 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 determined by mineral.tsteps 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, TK, 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.

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.