Part IV: Applied Earth Science | Chapter 2

Isotope Geology

Radioactive decay, geochronology, and the isochron method

4.1 Radioactive Decay Law

Derivation 1: The Decay Equation

Radioactive decay is a first-order process where the number of parent atoms $N$decreases exponentially with time:

$$\frac{dN}{dt} = -\lambda N$$

where $\lambda$ is the decay constant (probability of decay per unit time). Integrating:

$$N(t) = N_0 e^{-\lambda t}$$

The half-life is the time for half the parent atoms to decay:

$$t_{1/2} = \frac{\ln 2}{\lambda} = \frac{0.6931}{\lambda}$$

The number of daughter atoms produced is:

$$D^* = N_0 - N = N(e^{\lambda t} - 1)$$

If the system also contained initial daughter atoms $D_0$, the total daughter is:

$$D = D_0 + N(e^{\lambda t} - 1)$$

This is the fundamental geochronology equation. To determine age $t$, we need to measure $D$ and $N$ in the sample and either know or eliminate $D_0$.

Historical Context

Ernest Rutherford and Frederick Soddy (1902) formulated the decay law. Arthur Holmes (1911) was the first to use radioactive decay for geological dating. Clair Patterson (1956) determined the age of Earth as 4.55 Ga using lead isotopes from the Canyon Diablo meteorite, a landmark measurement that also led to his campaign against lead pollution.

4.2 The Isochron Method

Derivation 2: The Rb-Sr Isochron

The Rb-Sr system exploits the decay of $^{87}$Rb to $^{87}$Sr ($\lambda = 1.42 \times 10^{-11}$ yr$^{-1}$, $t_{1/2} = 48.8$ Gyr). Normalizing to the stable isotope $^{86}$Sr:

$$\frac{^{87}\text{Sr}}{^{86}\text{Sr}} = \left(\frac{^{87}\text{Sr}}{^{86}\text{Sr}}\right)_0 + \frac{^{87}\text{Rb}}{^{86}\text{Sr}}(e^{\lambda t} - 1)$$

This is a linear equation in $y = {^{87}\text{Sr}}/{^{86}\text{Sr}}$ vs$x = {^{87}\text{Rb}}/{^{86}\text{Sr}}$. Minerals from a single rock that crystallized at the same time from a homogeneous source will plot on a straight line (the isochron) with:

Slope = $(e^{\lambda t} - 1) \Rightarrow t = \frac{1}{\lambda}\ln(\text{slope} + 1)$
Y-intercept = $(^{87}\text{Sr}/^{86}\text{Sr})_0$ (initial ratio)

The beauty of the isochron method is that it simultaneously determines both the age and the initial isotope ratio, without requiring knowledge of $D_0$ a priori. It also provides a test of the closed-system assumption: if samples have remained closed since crystallization, they will be collinear.

4.3 The U-Pb System and Concordia

Derivation 3: The Concordia Diagram

The U-Pb system is the gold standard of geochronology because uranium has two long-lived isotopes that decay to different lead isotopes:

$$^{238}\text{U} \rightarrow ^{206}\text{Pb} + 8\alpha + 6\beta^- \quad (\lambda_{238} = 1.551 \times 10^{-10} \text{ yr}^{-1})$$

$$^{235}\text{U} \rightarrow ^{207}\text{Pb} + 7\alpha + 4\beta^- \quad (\lambda_{235} = 9.849 \times 10^{-10} \text{ yr}^{-1})$$

The two daughter/parent ratios define a parametric curve in time (the concordia):

$$x = \frac{^{206}\text{Pb}^*}{^{238}\text{U}} = e^{\lambda_{238}t} - 1$$

$$y = \frac{^{207}\text{Pb}^*}{^{235}\text{U}} = e^{\lambda_{235}t} - 1$$

A closed system plots on the concordia curve at the point corresponding to its age. If Pb loss occurred at time $t_1$ from a mineral crystallized at $t_0$, the data point moves off concordia along a discordia line. The upper intercept gives $t_0$ (crystallization) and the lower intercept gives $t_1$(Pb loss event). This dual-chronometer makes U-Pb uniquely powerful for detecting and correcting for open-system behavior.

4.4 K-Ar, Ar-Ar, and Radiocarbon Dating

Derivation 4: The K-Ar Age Equation

Potassium-40 decays by two pathways: $\beta^-$ emission to $^{40}$Ca (89.5%) and electron capture to $^{40}$Ar (10.5%). The total decay constant is$\lambda = \lambda_\beta + \lambda_e = 5.543 \times 10^{-10}$ yr$^{-1}$. The $^{40}$Ar produced is:

$$^{40}\text{Ar}^* = \frac{\lambda_e}{\lambda}\,^{40}\text{K}(e^{\lambda t} - 1)$$

Solving for age:

$$t = \frac{1}{\lambda}\ln\left(\frac{\lambda}{\lambda_e}\frac{^{40}\text{Ar}^*}{^{40}\text{K}} + 1\right)$$

The $^{40}$Ar/$^{39}$Ar method improves on K-Ar by irradiating the sample in a reactor to convert $^{39}$K to $^{39}$Ar, allowing both K and Ar to be measured as gases on the same instrument. Step-heating experiments reveal complex thermal histories.

Derivation 5: Radiocarbon Dating

Radiocarbon ($^{14}$C) is produced in the upper atmosphere by cosmic ray neutrons: $^{14}$N(n,p)$^{14}$C. It decays with $t_{1/2} = 5730$ years ($\lambda = 1.21 \times 10^{-4}$ yr$^{-1}$). The age of organic material is:

$$t = -\frac{1}{\lambda}\ln\left(\frac{A}{A_0}\right) = -8267 \ln\left(\frac{A}{A_0}\right)$$

where $A$ and $A_0$ are the current and initial $^{14}$C activities. Willard Libby (1949) developed radiocarbon dating, for which he received the Nobel Prize in Chemistry (1960). The method is applicable to ~50,000 years. Calibration against tree rings and corals corrects for past variations in atmospheric$^{14}$C production.

Accelerator Mass Spectrometry

Traditional radiocarbon dating measured the decay rate of $^{14}$C (requiring large samples and long counting times). Accelerator mass spectrometry (AMS), developed in the 1970s, directly counts $^{14}$C atoms, requiring only milligram samples and enabling dating of individual seeds, pollen grains, or amino acids. AMS has also enabled high-precision measurements of other cosmogenic nuclides ($^{10}$Be,$^{26}$Al, $^{36}$Cl) at natural abundance levels of $\sim 10^{-15}$.

The radiocarbon calibration curve (IntCal20) reveals significant departures from the simple decay model due to variations in cosmic ray flux (modulated by solar activity and geomagnetic field strength) and carbon cycle changes. The "radiocarbon plateau" during the Younger Dryas (~12,800-11,500 cal BP) complicates dating of this critical climate transition. The "bomb spike" from nuclear weapons testing (peaking at ~2x natural $^{14}$C in 1963) provides a precise chronostratigraphic marker and is used in forensic, environmental, and biomedical applications.

4.5 The Age of Earth and Oldest Materials

Clair Patterson (1956) used Pb isotopes from iron meteorites (representing primordial Pb) and oceanic sediments (representing modern Pb) to construct a Pb-Pb isochron giving an age of $4.55 \pm 0.07$ Ga. The Pb-Pb age equation:

$$\frac{(^{207}\text{Pb}/^{204}\text{Pb}) - (^{207}\text{Pb}/^{204}\text{Pb})_0}{(^{206}\text{Pb}/^{204}\text{Pb}) - (^{206}\text{Pb}/^{204}\text{Pb})_0} = \frac{1}{137.88}\frac{e^{\lambda_{235}t} - 1}{e^{\lambda_{238}t} - 1}$$

where 137.88 is the present-day $^{238}$U/$^{235}$U ratio. The oldest terrestrial minerals are 4.374 Ga zircons from Jack Hills, Australia (Wilde et al., 2001). The oldest whole-rock ages are ~4.03 Ga from the Acasta Gneiss, Canada. These ancient zircons preserve a record of the earliest crust and, through oxygen isotope studies, suggest liquid water existed on Earth's surface by 4.3 Ga.

4.6 Stable Isotope Geochemistry

Beyond radiogenic isotopes for dating, stable isotope ratios (H, C, N, O, S) provide information about temperatures, fluid sources, and biogeochemical processes. The fractionation factor between two phases A and B is:

$$\alpha_{A-B} = \frac{R_A}{R_B} = \frac{1000 + \delta_A}{1000 + \delta_B}$$

where $R$ is the heavy-to-light isotope ratio. The fractionation$\Delta_{A-B} = \delta_A - \delta_B \approx 1000\ln\alpha_{A-B}$ is temperature-dependent for equilibrium reactions:

$$1000\ln\alpha \approx \frac{A}{T^2} + \frac{B}{T} + C$$

The $\delta^{13}$C system is particularly useful in paleoenvironmental studies. Photosynthesis preferentially incorporates $^{12}$C, giving organic matter$\delta^{13}$C of about -25 per mille (C3 plants) or -13 per mille (C4 plants). Marine carbonates have $\delta^{13}$C near 0 per mille. Excursions in the$\delta^{13}$C record indicate perturbations to the global carbon cycle, such as the negative excursion at the PETM.

Sulfur Isotopes and Redox History

Bacterial sulfate reduction fractionates sulfur isotopes by 20-45 per mille, producing $^{34}$S-depleted sulfides. The $\delta^{34}$S record of seawater sulfate (from evaporites and carbonate-associated sulfate) tracks the balance between oxidative weathering and burial of reduced sulfur, providing constraints on atmospheric oxygen levels through Earth history. Mass-independent fractionation (MIF) of sulfur isotopes ($\Delta^{33}$S anomalies) in Archean rocks indicates the absence of an ozone layer before the Great Oxidation Event (~2.4 Ga).

The Sm-Nd System

The $^{147}$Sm-$^{143}$Nd system ($t_{1/2} = 106$ Gyr) is valuable because Sm and Nd are both rare earth elements with similar chemical behavior, making the system resistant to metamorphic disturbance. The $\varepsilon_{\text{Nd}}$ notation expresses the deviation from the chondritic evolution curve:

$$\varepsilon_{\text{Nd}} = \left(\frac{(^{143}\text{Nd}/^{144}\text{Nd})_{\text{sample}}}{(^{143}\text{Nd}/^{144}\text{Nd})_{\text{CHUR}}} - 1\right) \times 10^4$$

Positive $\varepsilon_{\text{Nd}}$ indicates derivation from a depleted mantle source (MORB typically +8 to +12), while negative values indicate crustal recycling or enriched sources. The model age $T_{\text{DM}}$ (depleted mantle age) gives the time when the crustal protolith separated from the depleted mantle.

Computational Lab: Isotope Geology

Radioactive Decay, Isochrons, and U-Pb Concordia

Python

script.py343 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================================================
# ISOTOPE GEOLOGY SIMULATION
# ============================================================
print("ISOTOPE GEOLOGY SIMULATION")
print("=" * 70)

# ============================================================
# 1. Radioactive Decay of Major Systems
# ============================================================
print()
print("1. Major Geochronological Systems")
print("-" * 70)

systems = [
    ("238U->206Pb", 1.551e-10, "4.468 Ga"),
    ("235U->207Pb", 9.849e-10, "0.704 Ga"),
    ("232Th->208Pb", 4.948e-11, "14.01 Ga"),
    ("87Rb->87Sr", 1.42e-11, "48.8 Ga"),
    ("147Sm->143Nd", 6.54e-12, "106 Ga"),
    ("40K->40Ar", 5.543e-10, "1.25 Ga"),
    ("176Lu->176Hf", 1.867e-11, "37.1 Ga"),
    ("187Re->187Os", 1.666e-11, "41.6 Ga"),
    ("14C->14N", 1.21e-4, "5730 yr"),
]

print(f"  {'System':20} {'lambda (yr^-1)':18} {'Half-life':14} {'Useful range'}")
print("  " + "-" * 65)

for name, lam, t_half in systems:
    if lam > 1e-5:
        useful = "0-50 kyr"
    elif lam > 1e-9:
        useful = "10 kyr - 4.6 Ga"
    else:
        useful = "100 Ma - 4.6 Ga"
    print(f"  {name:18} {lam:16.3e} {t_half:12} {useful}")

# ============================================================
# 2. Radioactive Decay Curves
# ============================================================
print()
print()
print("2. Decay and Ingrowth Curves")
print("-" * 55)

lam_238 = 1.551e-10
lam_235 = 9.849e-10

print(f"  {'Time (Ga)':12} {'N/N0 (238U)':14} {'N/N0 (235U)':14} {'235U/238U'}")
print("  " + "-" * 50)

for t_Ga in [0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 4.567]:
    t = t_Ga * 1e9
    N238 = np.exp(-lam_238 * t)
    N235 = np.exp(-lam_235 * t)
    ratio = N235 / N238 * (1/137.88)  # present day 235/238 is 1/137.88
    # At time t ago, the ratio was higher
    ratio_past = (1/137.88) * np.exp((lam_235 - lam_238) * t)
    print(f"  {t_Ga:10.3f} {N238:12.4f} {N235:12.4f} {ratio_past:10.4f}")

# ============================================================
# 3. Rb-Sr Isochron
# ============================================================
print()
print()
print("3. Rb-Sr Isochron Example")
print("-" * 65)

# Synthetic isochron: age = 2.0 Ga, initial 87Sr/86Sr = 0.7020
age_Ga = 2.0
lam_Rb = 1.42e-11
t = age_Ga * 1e9
slope = np.exp(lam_Rb * t) - 1
initial_ratio = 0.7020

print(f"  True age: {age_Ga} Ga")
print(f"  Isochron slope: {slope:.6f}")
print(f"  Initial 87Sr/86Sr: {initial_ratio}")
print()

# Synthetic mineral data
minerals = [
    ("Plagioclase", 0.05),
    ("Whole rock", 0.8),
    ("Biotite", 3.0),
    ("K-feldspar", 5.0),
    ("Muscovite", 12.0),
]

print(f"  {'Mineral':15} {'87Rb/86Sr':14} {'87Sr/86Sr':14} {'Calculated age (Ga)'}")
print("  " + "-" * 55)

for name, Rb_Sr in minerals:
    Sr87_86 = initial_ratio + Rb_Sr * slope
    # Calculate age from single measurement
    calc_age = np.log(1 + (Sr87_86 - initial_ratio) / Rb_Sr) / lam_Rb / 1e9
    print(f"  {name:13} {Rb_Sr:12.3f} {Sr87_86:12.6f} {calc_age:10.3f}")

# ============================================================
# 4. U-Pb Concordia
# ============================================================
print()
print()
print("4. U-Pb Concordia Points")
print("-" * 60)

print(f"  {'Age (Ga)':12} {'206Pb*/238U':14} {'207Pb*/235U':14}")
print("  " + "-" * 40)

for t_Ga in [0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5]:
    t = t_Ga * 1e9
    x = np.exp(lam_238 * t) - 1
    y = np.exp(lam_235 * t) - 1
    print(f"  {t_Ga:10.1f} {x:12.4f} {y:12.4f}")

# Discordant zircons
print()
print("  Discordant zircons (crystallized 3.0 Ga, Pb loss at 0.5 Ga):")
t0 = 3.0e9
t1 = 0.5e9

x_conc = np.exp(lam_238 * t0) - 1
y_conc = np.exp(lam_235 * t0) - 1
x_lower = np.exp(lam_238 * t1) - 1
y_lower = np.exp(lam_235 * t1) - 1

print(f"  Upper intercept: ({x_conc:.4f}, {y_conc:.4f}) -> {t0/1e9:.1f} Ga")
print(f"  Lower intercept: ({x_lower:.4f}, {y_lower:.4f}) -> {t1/1e9:.1f} Ga")

for frac in [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]:
    x_disc = x_conc * (1-frac) + x_lower * frac
    y_disc = y_conc * (1-frac) + y_lower * frac
    print(f"    Pb loss {frac*100:.0f}%: ({x_disc:.4f}, {y_disc:.4f})")

# ============================================================
# 5. Radiocarbon Dating
# ============================================================
print()
print()
print("5. Radiocarbon Dating Examples")
print("-" * 60)

lam_C14 = 1.21e-4  # yr^-1
t_half_C14 = 5730

print(f"  {'Sample':25} {'14C activity (%)':18} {'Age (yr BP)':14} {'Calendar ~'}")
print("  " + "-" * 65)

samples_c14 = [
    ("Modern wood", 100.0, "present"),
    ("Medieval artifact", 88.0, "~1000 AD"),
    ("Roman artifact", 78.0, "~200 BC"),
    ("Neolithic grain", 50.0, "~4000 BC"),
    ("Cave painting", 15.0, "~13500 BC"),
    ("Ice Age bone", 5.0, "~22000 BC"),
    ("Mammoth tusk", 1.0, "~35000 BC"),
    ("Near detection limit", 0.1, "~53000 BC"),
]

for name, activity_pct, cal in samples_c14:
    if activity_pct > 0:
        A_ratio = activity_pct / 100.0
        age = -np.log(A_ratio) / lam_C14
        print(f"  {name:23} {activity_pct:16.1f} {age:12.0f} {cal}")

# ============================================================
# PLOTS
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 11))

# Plot 1: Decay curves
ax = axes[0, 0]
t_range = np.linspace(0, 5, 500)  # Ga

for name, lam, color in [("238U", lam_238, 'blue'), ("235U", lam_235, 'red'),
                           ("87Rb", 1.42e-11, 'green'), ("40K", 5.543e-10, 'orange')]:
    N_frac = np.exp(-lam * t_range * 1e9)
    ax.plot(t_range, N_frac, color=color, linewidth=2, label=name)

ax.axvline(x=4.567, color='gray', linestyle=':', alpha=0.5, label='Age of Earth')
ax.set_xlabel('Time (Ga)')
ax.set_ylabel('N/N0 (fraction remaining)')
ax.set_title('Radioactive Decay Curves')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 5)

# Plot 2: Rb-Sr Isochron
ax = axes[0, 1]
x_iso = np.linspace(0, 15, 100)
y_iso = initial_ratio + slope * x_iso

# Plot isochron line
ax.plot(x_iso, y_iso, 'b-', linewidth=2, label=f'Isochron: t = {age_Ga} Ga')

# Plot mineral points
for name, Rb_Sr in minerals:
    Sr87_86 = initial_ratio + Rb_Sr * slope
    ax.plot(Rb_Sr, Sr87_86, 'ro', markersize=10)
    ax.annotate(name, (Rb_Sr, Sr87_86), textcoords="offset points",
                xytext=(5, 5), fontsize=8)

# Initial ratio
ax.axhline(y=initial_ratio, color='gray', linestyle='--', alpha=0.5,
           label=f'Initial ratio = {initial_ratio}')

ax.set_xlabel('87Rb/86Sr')
ax.set_ylabel('87Sr/86Sr')
ax.set_title('Rb-Sr Isochron Diagram')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# Plot 3: U-Pb Concordia
ax = axes[1, 0]
t_conc = np.linspace(0.05, 4.6, 500)
x_conc_arr = np.exp(lam_238 * t_conc * 1e9) - 1
y_conc_arr = np.exp(lam_235 * t_conc * 1e9) - 1

ax.plot(x_conc_arr, y_conc_arr, 'b-', linewidth=2, label='Concordia')

# Mark ages
for t_mark in [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5]:
    x_m = np.exp(lam_238 * t_mark * 1e9) - 1
    y_m = np.exp(lam_235 * t_mark * 1e9) - 1
    ax.plot(x_m, y_m, 'bo', markersize=6)
    ax.annotate(f'{t_mark:.1f} Ga', (x_m, y_m), textcoords="offset points",
                xytext=(5, -10), fontsize=7)

# Discordia line
x_upper = np.exp(lam_238 * 3.0e9) - 1
y_upper = np.exp(lam_235 * 3.0e9) - 1
x_low = np.exp(lam_238 * 0.5e9) - 1
y_low = np.exp(lam_235 * 0.5e9) - 1

fracs = np.linspace(0, 1, 50)
x_disc_arr = x_upper * (1-fracs) + x_low * fracs
y_disc_arr = y_upper * (1-fracs) + y_low * fracs
ax.plot(x_disc_arr, y_disc_arr, 'r--', linewidth=2, label='Discordia (3.0-0.5 Ga)')

# Discordant points
for f_pt in [0, 0.2, 0.4, 0.6, 0.8, 1.0]:
    xp = x_upper * (1-f_pt) + x_low * f_pt
    yp = y_upper * (1-f_pt) + y_low * f_pt
    ax.plot(xp, yp, 'rs', markersize=8)

ax.set_xlabel('206Pb*/238U')
ax.set_ylabel('207Pb*/235U')
ax.set_title('U-Pb Concordia Diagram')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# Plot 4: Radiocarbon decay
ax = axes[1, 1]
t_c14 = np.linspace(0, 60000, 500)
activity = 100 * np.exp(-lam_C14 * t_c14)

ax.plot(t_c14/1000, activity, 'b-', linewidth=2)
ax.axhline(y=50, color='red', linestyle=':', alpha=0.5, label=f't_1/2 = {t_half_C14} yr')
ax.axhline(y=1, color='gray', linestyle='--', alpha=0.5, label='Detection limit (~1%)')

# Mark events
events = [(0, 'Present'), (5730, '1 half-life'), (11460, '2 half-lives'),
          (21000, 'LGM'), (40000, 'Practical limit')]
for t_ev, label in events:
    a_ev = 100 * np.exp(-lam_C14 * t_ev)
    ax.plot(t_ev/1000, a_ev, 'ro', markersize=6)
    ax.annotate(label, (t_ev/1000, a_ev), textcoords="offset points",
                xytext=(5, 5), fontsize=8)

ax.set_xlabel('Time (kyr BP)')
ax.set_ylabel('14C Activity (% modern)')
ax.set_title('Radiocarbon Decay Curve')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 60)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
plt.close()
print()
print()
print("6. Sm-Nd Model Ages and Epsilon Notation")
print("-" * 60)

lam_Sm = 6.54e-12  # yr^-1
Nd_CHUR_0 = 0.512638  # present CHUR
Sm_Nd_CHUR = 0.1967   # Sm/Nd ratio of CHUR

print(f"  {'Sample':18} {'147Sm/144Nd':16} {'143Nd/144Nd':16} {'epsilon_Nd':14} {'T_DM (Ga)'}")
print("  " + "-" * 70)

samples_nd = [
    ("MORB", 0.220, 0.51315),
    ("OIB (Hawaii)", 0.185, 0.51290),
    ("Archean gneiss", 0.110, 0.51080),
    ("Continental sed.", 0.115, 0.51190),
    ("Avg upper crust", 0.118, 0.51200),
]

# Depleted mantle present-day values
Nd_DM = 0.51330
Sm_Nd_DM = 0.2137

for name, SmNd, NdNd in samples_nd:
    eps = ((NdNd / Nd_CHUR_0) - 1) * 10000
    # Model age: when did sample's Nd evolution diverge from DM?
    # (143Nd/144Nd)_sample - (143Nd/144Nd)_DM at time T_DM = 0
    # Iterate to find T_DM
    T_DM = 0
    for t_try_Ga in np.linspace(0, 4.5, 1000):
        t_try = t_try_Ga * 1e9
        Nd_sample_t = NdNd - SmNd * (np.exp(lam_Sm * t_try) - 1)
        Nd_DM_t = Nd_DM - Sm_Nd_DM * (np.exp(lam_Sm * t_try) - 1)
        if Nd_sample_t <= Nd_DM_t:
            T_DM = t_try_Ga
            break
    print(f"  {name:16} {SmNd:14.3f} {NdNd:14.5f} {eps:12.1f} {T_DM:8.2f}")

# ============================================================
# 7. Stable Isotope Fractionation
# ============================================================
print()
print()
print("7. Oxygen Isotope Fractionation: Quartz-Water")
print("-" * 55)
print(f"  {'T (C)':10} {'T (K)':10} {'1000*ln(alpha)':16} {'delta_qtz - delta_w'}")
print("  " + "-" * 50)

for T_C in [100, 200, 300, 400, 500, 600, 700, 800, 900, 1000]:
    T_K = T_C + 273.15
    # Clayton et al. (1972): 1000*ln(alpha) = 3.38*(10^6/T^2) - 3.40
    frac = 3.38 * 1e6 / T_K**2 - 3.40
    print(f"  {T_C:8} {T_K:8.1f} {frac:14.2f} {frac:10.2f} permil")

print()
print()
print("[Chart saved: decay curves, Rb-Sr isochron, U-Pb concordia, radiocarbon]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

4.7 Cosmogenic Nuclide Dating

Cosmic rays produce radioactive nuclides ($^{10}$Be, $^{26}$Al, $^{36}$Cl,$^{14}$C, $^{3}$He) in rocks exposed at Earth's surface. The concentration of these nuclides increases with exposure time:

$$N(t) = \frac{P}{\lambda + \rho\varepsilon/\Lambda}(1 - e^{-(\lambda + \rho\varepsilon/\Lambda)t})$$

where $P$ is the production rate (atoms/g/yr), $\lambda$ is the decay constant,$\varepsilon$ is erosion rate, $\rho$ is rock density, and $\Lambda$ is the cosmic ray attenuation length (~160 g/cm$^2$). For negligible erosion:$N \approx Pt$ for short exposures, saturating at $P/\lambda$ after several half-lives.

Applications include dating glacial moraines, lava flows, and fault scarps; measuring erosion rates from $^{10}$Be/$^{26}$Al ratios in river sediments; and determining burial ages from the decay of cosmogenic nuclides after sediment is shielded from cosmic rays. The technique has been pivotal in establishing the chronology of Quaternary glaciations with precision of ~5-10%.

The $^{10}$Be/$^{26}$Al ratio provides an additional constraint because these nuclides have different half-lives ($t_{1/2}$ of 1.39 Myr and 0.72 Myr respectively). A surface at production steady-state has a fixed ratio of ~6.75, while burial reduces this ratio as $^{26}$Al decays faster than $^{10}$Be. This two-nuclide method constrains both exposure and burial history.

Cosmogenic Nuclide Dating and Exposure Age Calculations

Python

script.py76 lines

import numpy as np

print("COSMOGENIC NUCLIDE DATING")
print("=" * 70)

print()
print("1. 10Be Exposure Age Calculation")
print("-" * 55)

# 10Be production rate and decay
P_Be = 4.0  # atoms/g/yr at sea level high latitude
lambda_Be = np.log(2) / (1.39e6)  # decay constant (yr^-1)
rho_rock = 2.7  # g/cm^3
Lambda = 160  # g/cm^2 attenuation length

print(f"  Production rate: {P_Be} atoms/g/yr (SLHL)")
print(f"  Half-life: 1.39 Myr, lambda = {lambda_Be:.4e} yr^-1")
print()
print(f"  {'Exposure (kyr)':18} {'[10Be] (at/g)':16} {'Saturation %'}")
print("  " + "-" * 45)

N_sat = P_Be / lambda_Be

for t_kyr in [1, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000]:
    t = t_kyr * 1000
    N = P_Be / lambda_Be * (1 - np.exp(-lambda_Be * t))
    sat_pct = N / N_sat * 100
    print(f"  {t_kyr:16} {N:14.2e} {sat_pct:8.1f}")

# Erosion effect
print()
print()
print("2. Effect of Erosion on Apparent Ages")
print("-" * 60)

print(f"  {'True age (kyr)':18} {'e=0 mm/kyr':14} {'e=1 mm/kyr':14} {'e=10 mm/kyr'}")
print("  " + "-" * 55)

for t_kyr in [10, 50, 100, 200, 500, 1000]:
    t = t_kyr * 1000
    concs = []
    for e_mm in [0, 1, 10]:
        e = e_mm * 1e-4  # mm/kyr -> cm/yr
        erosion_term = rho_rock * e / Lambda
        N = P_Be / (lambda_Be + erosion_term) * (1 - np.exp(-(lambda_Be + erosion_term) * t))
        # Apparent age if erosion ignored
        if N < N_sat:
            t_app = -np.log(1 - N * lambda_Be / P_Be) / lambda_Be / 1000
        else:
            t_app = float('inf')
        concs.append(t_app)
    print(f"  {t_kyr:16} {concs[0]:12.0f} {concs[1]:12.0f} {concs[2]:12.0f} kyr apparent")

# 10Be/26Al burial dating
print()
print()
print("3. 10Be/26Al Burial Dating")
print("-" * 55)

lambda_Al = np.log(2) / (0.72e6)
ratio_prod = 6.75  # production ratio 26Al/10Be

print(f"  Production ratio: 26Al/10Be = {ratio_prod}")
print()
print(f"  {'Burial time (Ma)':20} {'26Al/10Be ratio':18} {'Ratio/production'}")
print("  " + "-" * 50)

for t_Ma in [0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0]:
    t = t_Ma * 1e6
    ratio = ratio_prod * np.exp(-(lambda_Al - lambda_Be) * t)
    frac = ratio / ratio_prod
    print(f"  {t_Ma:18.1f} {ratio:16.2f} {frac:8.3f}")

print()
print("Complete: cosmogenic nuclide calculations finished.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Geochemistry Next: Hydrogeology →