gdtchron.he

Module for forward modeling of (U-Th)/He ages.

This code follows the workflow from Ketcham (2005) and includes the alpha correction from Ketcham et al. (2011)

Functions

alpha_correction(stop_distance, radius)	Calculate alpha ejection correction factor, after Ketcham et al. (2011).
calc_age(he_molg, u238_molg, u235_molg, th_molg)	Calculate (U-Th)/He age from U, Th, and He concentrations.
calc_beta(diffusivity, node_spacing, ...)	Calculate beta, a substitution term from Ketcham (2005).
calc_diffusivity(temperature, system)	Calculate diffusivity from temperature and system diffusion parameters.
calc_he_production_rate(u238_molg, ...)	Calculate instantaneous He production rate as a function of U and Th.
calc_node_positions(node_spacing, radius)	Calculate node positions given spacing and radius.
forward_model_he(temps, tsteps, system, u, ...)	Forward model a (U-Th)/He age for a particular time-temperature path.
he_profile(avg_temps, tsteps, system, ...[, ...])	Solve for the helium profile of the grain at the present day.
model_alpha_ejection(node_positions, ...)	Model retained fraction of He after alpha ejection.
profile_to_age(x, node_positions, radius, ...)	Calculate age based on final helium profile.
sum_he_shells(x, node_positions, radius)	Sum He produced within all nodes of the modeled crystal.
tridiag_banded(a, b, c, diag_length[, dtype])	Set up tridiagonal matrix in banded form from values for the 3 diagonals.
u_th_ppm_to_molg(u_ppm, th_ppm)	Convert concentrations of U and Th from ppm to mol / g.

gdtchron.he.alpha_correction(stop_distance, radius)

Calculate alpha ejection correction factor, after Ketcham et al. (2011).

Uses Equation 2 of Ketcham et al. (2011)

Parameters:

• stop_distance (float or NumPy array of floats) – Stopping distance(s) for particular isotopic system(s) (micrometers).

• radius (float) – Radius of the grain (micrometers)

Returns:

tau – Alpha correction factor(s) (F_T in Ketcham et al. (2011))

Return type:

float or NumPy array of floats

gdtchron.he.calc_age(he_molg, u238_molg, u235_molg, th_molg)

Calculate (U-Th)/He age from U, Th, and He concentrations.

Uses Equation 15 from Ketcham (2005). Note that no alpha correction is applied here. Instead, the alpha correction is applied to the amounts of each parent isotope fed into this function, following Ketcham et al. (2011).

Parameters:

• he_molg (float) – Amount of He (mol / g)

• u238_molg (float) – Amount of U238 (mol / g)

• u235_molg (float) – Amount of U235 (mol / g)

• th_molg (float) – Amount of Th232 (mol / g)

Returns:

age – Calculated (U-Th)/He age (Ma)

Return type:

float

gdtchron.he.calc_beta(diffusivity, node_spacing, time_interval)

Calculate beta, a substitution term from Ketcham (2005).

The equation uses diffusivity, the spacing of nodes within the modeled grain, and the timestep duration. It comes from the in-text equation between equations 20 and 21 in Ketcham (2005).

Parameters:

• diffusivity (float) – Diffusivity (micrometers^2 / Myr)

• node_spacing (float) – Spacing of nodes in the modeled crystal (micrometers)

• time_interval (float) – Timestep in the thermal model (Myr)

Returns:

beta – Beta (unitless), after Ketcham (2005).

Return type:

float

gdtchron.he.calc_diffusivity(temperature, system)

Calculate diffusivity from temperature and system diffusion parameters.

After Reiners and Brandon (2006), with PV_a term assumed to be 0.

Parameters:

• temperature (float) – Temperature (K)

• system (string) – String indicating whether to use parameters for the apatite system (‘AHe’) or the zircon system (‘ZHe’)

Returns:

kappa – Diffusivity (micrometers^2 / Myr)

Return type:

float

gdtchron.he.calc_he_production_rate(u238_molg, u235_molg, th_molg)

Calculate instantaneous He production rate as a function of U and Th.

Parameters:

• u238_molg (float) – U238 (mol / g)

• u235_molg (float) – U235 (mol / g)

• th_molg (float) – Th232 (mol / g)

Returns:

he_production – He production (mol * g^-1 * Myr^-1)

Return type:

float

gdtchron.he.calc_node_positions(node_spacing, radius)

Calculate node positions given spacing and radius.

Follows Ketcham (2005); see Figure 8 for an example. The first node is a half-spacing away from the center and the last node is a full spacing away from the edge of the grain.

Parameters:

• node_spacing (float) – Distance between nodes in the crystal (micrometers)

• radius (float) – Radius of the grain (micrometers)

Returns:

node_positions – Radial positions of each modeled node (micrometers)

Return type:

NumPy array of floats

gdtchron.he.forward_model_he(temps, tsteps, system, u, th, radius, nodes=513, initial_x=None, return_all=False)

Forward model a (U-Th)/He age for a particular time-temperature path.

Uses finite difference method for diffusion within a sphere as described in Ketcham (2005). Applies alpha ejection correction after Ketcham et al. (2011). Returns corrected age and optionally additional relevant values.

Parameters:

• temps (NumPy array of length n, where n > 1) – List of temperatures (K) for the time-temperature path

• tsteps (NumPy array of floats with length n, where n > 1) – Array of times (Myrs BP) in chronological (descending) order. First time is start of first timestep, last time is end of last timestep. Each pair of adjacent times composes a timestep. Each time corresponds to a temperature in temps.

• system (string) – Isotopic system. Current options are ‘AHe’ and ‘ZHe’.

• u (float) – U concentration (ppm)

• th (float) – Th concentration (ppm)

• radius (float) – Radius of the grain (micrometers)

• nodes (float, optional) – Number of nodes to model within the crystal. (default: 513)

• initial_x (NumPy array of floats or None, optional) – Array containing the initial values for x (mol * micrometers / g). If all values in the array are np.nan or initial_x is set to None, this function assumes that the initial values for x are 0 for all nodes. Note that x is referred to as u in Ketcham (2005). (default: None)

• return_all (bool) – Boolean indicating whether to return additional values calculated by the functions as a part of the forward model. If False, only the corrected age is returned. (default: False)

Returns:

• age_corrected (float) – Corrected (U-Th) / He age

• age_uncorrected (float (optional)) – (U-Th) / He age without alpha correction. Returned if return_all is True.

• he_nmolg (float (optional)) – Total amount (nmol / g) of He within the modeled crystal. Returned if return_all is True.

• position_normalized (NumPy array of floats (optional)) – Radial positions of each modeled node, normalized so that the radius has a radial position of 1. Returned if return_all is True.

• v_normalized (NumPy array of floats (optional)) – Radial profile of He, normalized so that the position with the highest concentration has a value of 1. Returned if return_all is True.

• x (NumPy array of floats (optional)) – Matrix x solved for using Equation 21 Ketcham (2005). Equivalent to the concentration times the node position. Returned if return_all is True.

gdtchron.he.he_profile(avg_temps, tsteps, system, radius, uth_molg, node_information, initial_x=None)

Solve for the helium profile of the grain at the present day.

Solves for final x in Equation 21 of Ketcham (2005). Note that in Ketcham (2005), this variable is referred to as u, but x is used here to distinguish this variable from the uranium-related variables.

Parameters:

• avg_temps (NumPy array of floats with length n) – List of average temperatures (K) for each timestep for the time-temperature path

• tsteps (NumPy array of floats with length n + 1) – Array of times (Myrs BP) in chronological (descending) order. First time is start of first timestep, last time is end of last timestep. Each pair of adjacent times composes a timestep.

• system (string) – Isotopic system. Current options are ‘AHe’ and ‘ZHe’.

• radius (float) – Radius of the grain (micrometers)

• uth_molg (tuple) –

  u238_molgfloat

  The concentration of U238 (mol / g),

  u235_molgfloat

  The concentration of U235 (mol / g)

  th_molgfloat

  The concentration of Th (mol / g)

• node_information (tuple) –

  nodesfloat

  Number of nodes to model within the crystal. The default is 513.

  node_spacingfloat

  Spacing of nodes in the modeled crystal (micrometers)

  node_positionsNumPy array of floats

  Radial positions of each modeled node (micrometers)

• initial_x (NumPy array of floats or None, optional) – Array containing the initial values for x (mol * micrometers / g). If all values in the array are np.nan or initial_x is set to None, this function assumes that the initial values for x are 0 for all nodes (default: None)

Returns:

x – Matrix x solved for using Equation 21 Ketcham (2005). Equivalent to the concentration times the node position.

Return type:

NumPy array of floats

gdtchron.he.model_alpha_ejection(node_positions, stop_distance, radius)

Model retained fraction of He after alpha ejection.

Calculations from in-text equations in Ketcham (2005).

Parameters:

• node_positions (NumPy array of floats) – Radial positions of each modeled node (micrometers)

• stop_distance (float) – Stopping distance for particular isotopic system (micrometers).

• radius (float) – Radius of the grain (micrometers)

Returns:

retained_fraction_edge – Fraction of He retained after alpha ejection for each node position

Return type:

NumPy array of floats

gdtchron.he.profile_to_age(x, node_positions, radius, uth_molg, stop_distances)

Calculate age based on final helium profile.

Calculate corrected and uncorrected age based on final He profile. Also calculates normalized He profile and node positions.

Parameters:

• x (NumPy array of floats) – Matrix x (mol * micrometers / g) solved for using Equation 21 from Ketcham (2005). Equivalent to the concentration times the node position. In Ketcham (2005), this variable is referred to as u, but x is used here to distinguish this variable from the uranium-related variables.

• node_positions (NumPy array of floats) – Radial positions of each modeled node (micrometers)

• radius (float) – Radius of the grain (micrometers)

• uth_molg (tuple) –

  u238_molgfloat

  The concentration of U238 (mol / g),

  u235_molgfloat

  The concentration of U235 (mol / g)

  th_molgfloat

  The concentration of Th (mol / g)

• stop_distances (NumPy array of floats) – Stopping distances (micrometers) for (in order) U238, U235, and Th232

Returns:

• age_corrected (float) – Corrected (U-Th) / He age

• age_uncorrected (Float) – (U-Th) / He age without alpha correction

• he_molg (float) – Total amount (mol / g) of He within the modeled crystal.

• position_normalized (NumPy array of floats) – Radial positions of each modeled node, normalized so that the radius has a radial position of 1

• v_normalized (NumPy array of floats) – Radial profile of He, normalized so that the position with the highest concentration has a value of 1

gdtchron.he.sum_he_shells(x, node_positions, radius)

Sum He produced within all nodes of the modeled crystal.

Uses substition for He concentration after Ketcham (2005). Converts radial profile of He to system of shells, so He is weighted by volume of shell.

Parameters:

• x (NumPy array of floats) – Matrix x (mol * micrometers / g) solved for using Equation 21 from Ketcham (2005). Equivalent to the concentration times the node position. In Ketcham (2005), this variable is referred to as u, but x is used here to distinguish this variable from the uranium-related variables.

• node_positions (NumPy array of floats) – Radial positions of each modeled node (micrometers)

• radius (float) – Radius of the grain (micrometers)

Returns:

• he_molg (float) – Total amount (mol / g) of He within the modeled crystal.

• v (NumPy array of floats) – Radial profile of He (mol / g)

gdtchron.he.tridiag_banded(a, b, c, diag_length, dtype=<class 'numpy.float32'>)

Set up tridiagonal matrix in banded form from values for the 3 diagonals.

For example, a tridiagonal matrix in banded form with values a, b, c for its diagonals and with a principal diagonal of length 6 appears as follows: [0, a, a, a, a, a] [b, b, b, b, b, b] [c, c, c, c, c, 0]

Parameters:

• a (float) – First diagonal value

• b (float) – Second diagonal value

• c (float) – Third diagonal value

• diag_length (float) – Length of principal diagonal

• dtype (type) – Type of numbers in the matrix (default: np.float32). 32-bit floats are preferred to save memory when running large numbers of forward models.

Returns:

tridiag_matrix – Tridiagonal matrix

Return type:

NumPy ndarray

gdtchron.he.u_th_ppm_to_molg(u_ppm, th_ppm)

Convert concentrations of U and Th from ppm to mol / g.

Parameters:

• u_ppm (float) – U concentration (ppm)

• th_ppm (float) – Th concentration (ppm)

Returns:

• u238_molg (float) – U-238 (mol / g)

• u235_molg (float) – U-235 (mol / g)

• th_molg (float) – Th-232 (mol / g)