Part IV: Applied Earth Science | Chapter 1

Geochemistry

Thermodynamics of mineral equilibria, element classification, and rare earth elements

4.1 Thermodynamics of Mineral Equilibria

Geochemistry applies chemical thermodynamics to understand which minerals are stable under given conditions of temperature, pressure, and composition. The fundamental criterion for equilibrium is minimization of the Gibbs free energy.

Derivation 1: Gibbs Free Energy and Mineral Stability

The Gibbs free energy of a reaction is:

$$\Delta G = \Delta H - T\Delta S$$

where $\Delta H$ is enthalpy change and $\Delta S$ is entropy change. At constant temperature and pressure, a reaction proceeds spontaneously when $\Delta G < 0$. The pressure and temperature dependence of $\Delta G$ is:

$$\Delta G(T, P) = \Delta G^0 + \int_{T_0}^{T} \Delta C_P \, dT' - T\int_{T_0}^{T} \frac{\Delta C_P}{T'} dT' + \int_{P_0}^{P} \Delta V \, dP'$$

For solids at moderate pressures, the volume term $\Delta V \Delta P$ is often small. The equilibrium boundary in P-T space is found by setting $\Delta G = 0$. The Clausius-Clapeyron equation gives the slope:

$$\frac{dP}{dT} = \frac{\Delta S}{\Delta V}$$

For example, the diamond-graphite transition has a positive Clapeyron slope because diamond has smaller volume ($\Delta V < 0$ for graphite$\rightarrow$diamond) and lower entropy ($\Delta S < 0$), so $dP/dT > 0$. Diamond is stable at high P (above ~4 GPa at 1000$°$C).

Historical Context

J. Willard Gibbs (1876-1878) established the thermodynamic framework. Norman Bowen (1928) applied it to igneous petrology with his reaction series. H.J. Greenwood, J.B. Thompson, and others developed the thermodynamic databases used in modern phase equilibrium software (THERMOCALC, Perple_X).

4.2 Chemical Potential and Activity

Derivation 2: The Equilibrium Constant from Chemical Potentials

For a reaction involving species with chemical potentials $\mu_i$, equilibrium requires:

$$\sum_i \nu_i \mu_i = 0$$

where $\nu_i$ are stoichiometric coefficients (positive for products, negative for reactants). The chemical potential of component $i$ in a solution is:

$$\mu_i = \mu_i^0(T, P) + RT \ln a_i$$

where $a_i$ is the thermodynamic activity. Substituting into the equilibrium condition:

$$\Delta G^0 = -RT \ln K, \quad K = \prod_i a_i^{\nu_i}$$

This is the fundamental relationship between the standard Gibbs energy change and the equilibrium constant $K$. The van't Hoff equation gives the temperature dependence:

$$\frac{d\ln K}{dT} = \frac{\Delta H^0}{RT^2}$$

For endothermic reactions ($\Delta H^0 > 0$), $K$ increases with temperature, favoring products at higher T.

4.3 Goldschmidt Classification

Derivation 3: Partition Coefficients and Element Distribution

Victor Goldschmidt (1937) classified elements into four groups based on their geochemical affinity during planetary differentiation:

Lithophile (rock-loving): Prefer silicate/oxide phases. Si, Al, Ca, Na, K, Mg, Ti, Mn, Cr, U, Th, REE
Siderophile (iron-loving): Prefer metallic iron. Fe, Co, Ni, Pt, Pd, Au, Re, Os, Ir, W, Mo
Chalcophile (sulfide-loving): Prefer sulfide phases. Cu, Zn, Pb, Ag, As, Sb, Bi, Cd, Hg
Atmophile (atmosphere-loving): Concentrated in gas phase. N, noble gases, C (partially), H

The distribution of an element between two phases is quantified by the partition (distribution) coefficient:

$$D_i^{\text{crystal/melt}} = \frac{C_i^{\text{crystal}}}{C_i^{\text{melt}}}$$

If $D \gg 1$, the element is compatible (concentrated in the solid); if $D \ll 1$, it is incompatible (concentrated in the melt). The bulk partition coefficient for a multi-phase assemblage is the weighted average:

$$D_0 = \sum_j X_j D_j$$

where $X_j$ is the weight fraction and $D_j$ is the partition coefficient for mineral phase $j$. During fractional crystallization, the concentration in the residual melt follows:

$$C_L = C_0 F^{(D-1)}$$

where $F$ is the fraction of melt remaining. This is the Rayleigh fractionation equation, which explains why incompatible elements (e.g., K, Rb, U, Th) become enriched in evolved melts and ultimately in the continental crust.

4.4 Rare Earth Element Geochemistry

Derivation 4: REE Patterns from Partial Melting

The rare earth elements (La through Lu) have smoothly varying ionic radii and uniformly $3+$ charge (except Ce$^{4+}$ and Eu$^{2+}$). Their concentrations normalized to chondritic meteorite values reveal diagnostic patterns. For batch (equilibrium) partial melting:

$$\frac{C_L}{C_0} = \frac{1}{D_0 + F(1 - D_0)}$$

where $D_0$ is the bulk partition coefficient and $F$ is the melt fraction. Since light REE (LREE: La, Ce, Nd) are more incompatible than heavy REE (HREE: Er, Yb, Lu), partial melts of the mantle are LREE-enriched relative to chondrites.

The Eu anomaly ($\text{Eu}/\text{Eu}^* = \text{Eu}_N / \sqrt{\text{Sm}_N \times \text{Gd}_N}$) is diagnostic: negative anomalies indicate plagioclase fractionation (Eu$^{2+}$substitutes for Ca$^{2+}$ in plagioclase), characteristic of evolved continental crust.

Derivation 5: The Lanthanide Contraction and Ionic Radii

The systematic decrease in ionic radius across the lanthanide series arises from imperfect shielding of $4f$ electrons. As nuclear charge increases from La ($Z = 57$) to Lu ($Z = 71$), each added $4f$ electron provides incomplete shielding of the increasing nuclear charge, resulting in a contraction:

$$r_{\text{ionic}}(\text{REE}^{3+}) \approx r_{\text{La}^{3+}} - k(Z - 57)$$

with $r_{\text{La}^{3+}} = 1.16$ angstrom and $r_{\text{Lu}^{3+}} = 0.977$ angstrom (Shannon, 1976). This smooth variation means that minerals preferentially incorporate either LREE (large sites: apatite, monazite) or HREE (small sites: garnet, zircon), producing the characteristic REE patterns used in petrogenetic modeling.

4.5 Applications: Phase Diagrams and Geothermobarometry

Geothermobarometry uses the temperature and pressure dependence of mineral equilibria to determine the P-T conditions of rock formation. The garnet-biotite Fe-Mg exchange thermometer is based on:

$$K_D = \frac{(Fe/Mg)^{\text{Grt}}}{(Fe/Mg)^{\text{Bt}}}$$

which is temperature-dependent: $\ln K_D = A/T + B$. Similarly, the Al content of orthopyroxene coexisting with garnet is pressure-sensitive, providing a barometer. Together, thermobarometry constrains the P-T paths of metamorphic rocks, revealing burial, exhumation, and tectonic processes.

Modern Computational Petrology

Software packages like THERMOCALC (Powell and Holland), Perple_X (Connolly), and MELTS (Ghiorso and Sack) compute phase equilibria from internally consistent thermodynamic databases. Pseudosections (P-T phase diagrams calculated for a specific bulk composition) show the stable mineral assemblage at each point in P-T space, enabling direct comparison with observed mineral assemblages.

The concept of chemical potential diagrams extends phase diagrams to open systems. For example, T-$f$O$_2$ diagrams show the stability of iron-bearing minerals as a function of temperature and oxygen fugacity, critical for understanding ore formation and metamorphic redox reactions. Activity-activity diagrams for aqueous systems predict mineral stability in weathering and hydrothermal environments.

The trace element and isotopic composition of zircon (ZrSiO$_4$) has revolutionized geochemistry. Zircon's chemical and physical robustness preserves records of crystallization age (U-Pb), temperature (Ti-in-zircon thermometer: $T(°C) = -4800/[\log(\text{Ti/ppm}) - 5.711] - 273$), source composition ($\varepsilon_{\text{Hf}}$), and oxygen isotopes ($\delta^{18}$O). Detrital zircon studies have reconstructed the growth history of the continental crust back to 4.37 Ga.

4.6 Multi-Element Diagrams and Mantle Reservoirs

Extended trace element diagrams (spider diagrams) plot mantle-normalized concentrations of many elements ordered by increasing incompatibility. Different tectonic settings produce characteristic patterns:

N-MORB: Flat to LREE-depleted, reflecting the depleted mantle source
E-MORB and OIB: LREE-enriched with smooth patterns, from enriched mantle plumes
Island arc basalt: Enriched in LILE (K, Rb, Ba, Sr) with Nb-Ta negative anomalies, due to slab fluid addition
Continental flood basalt: Variable, often with crustal contamination signatures

The mantle is geochemically heterogeneous, comprising several end-member reservoirs identified by radiogenic isotopes: DMM (Depleted MORB Mantle), EM1 and EM2 (Enriched Mantle, from recycled sediments or lower continental lithosphere), HIMU (High-$\mu$, from recycled oceanic crust), and FOZO (focal zone of mantle plumes). These reservoirs reflect 4.5 Gyr of plate tectonic recycling and mantle mixing.

The Continental Crust Paradox

The continental crust has an average composition close to andesite ($\sim$60% SiO$_2$), yet it is ultimately derived from basaltic melts of the mantle. Rudnick and Gao (2003) proposed that the crust differentiates through intracrustal melting and delamination of dense mafic residues. The heat-producing elements (U, Th, K) are strongly enriched in the upper crust, with $\sim$40% of Earth's total heat production concentrated in just 0.6% of its mass.

Gibbs Phase Rule in Petrology

The phase rule constrains the number of degrees of freedom $F$ in a system:

$$F = C - P + 2$$

where $C$ is the number of components and $P$ is the number of phases. In a one-component system (e.g., SiO$_2$), two phases coexist along a univariant line (F=1), and three phases at an invariant point (F=0). The Al$_2$SiO$_5$ system (kyanite-sillimanite-andalusite) has a triple point at approximately 500$°$C and 3.8 kbar, used as a key constraint in metamorphic petrology.

Computational Lab: Geochemistry

Thermodynamics, REE Patterns, and Fractional Crystallization

Python

script.py345 lines

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

# ============================================================
# GEOCHEMISTRY SIMULATION
# ============================================================
print("GEOCHEMISTRY SIMULATION")
print("=" * 70)

R = 8.314  # gas constant

# ============================================================
# 1. Gibbs Free Energy and Equilibrium
# ============================================================
print()
print("1. Mineral Reaction Equilibria")
print("-" * 60)

reactions = [
    ("Calcite -> Aragonite", 1.0, -4.1, -12.7),
    ("Graphite -> Diamond", 2.9, 3.3, -1.88),
    ("Albite -> Jadeite + Quartz", 17.3, -35.6, -17.0),
    ("Kyanite -> Sillimanite", 7.1, 7.0, -0.53),
    ("Muscovite dehydration", 100.0, 150.0, -7.0),
]

print(f"  {'Reaction':35} {'dH (kJ)':12} {'dS (J/K)':12} {'dV (cm3)':12} {'T_eq (C)'}")
print("  " + "-" * 75)

for name, dH_kJ, dS_JK, dV_cm3 in reactions:
    dH = dH_kJ * 1000
    # T_eq at 1 atm: dG = dH - T*dS = 0 -> T = dH/dS
    if abs(dS_JK) > 0.1:
        T_eq = dH / dS_JK
        T_eq_C = T_eq - 273.15
    else:
        T_eq_C = float('inf')
    print(f"  {name:33} {dH_kJ:10.1f} {dS_JK:10.1f} {dV_cm3:10.2f} {T_eq_C:10.0f}")

# Clapeyron slopes
print()
print("  Clapeyron slopes (dP/dT):")
for name, dH_kJ, dS_JK, dV_cm3 in reactions:
    if abs(dV_cm3) > 0.01:
        dV_m3 = dV_cm3 * 1e-6  # cm^3 -> m^3
        dPdT = dS_JK / dV_m3  # Pa/K
        dPdT_bar = dPdT / 1e5  # bar/K
        print(f"    {name}: {dPdT_bar:.1f} bar/K = {dPdT_bar/10:.2f} kbar/K")

# ============================================================
# 2. Partition Coefficients and Fractional Crystallization
# ============================================================
print()
print()
print("2. Rayleigh Fractional Crystallization")
print("-" * 65)

elements = [
    ("Ni", 10.0, 2000),   # compatible
    ("Cr", 7.0, 400),
    ("Sr", 2.0, 200),
    ("Ba", 0.1, 50),      # incompatible
    ("Rb", 0.01, 10),
    ("U", 0.005, 0.5),
    ("La", 0.02, 5.0),
    ("Yb", 0.5, 3.0),
]

print(f"  {'Element':10} {'D':8}", end="")
for F in [1.0, 0.8, 0.6, 0.4, 0.2, 0.1, 0.05]:
    print(f"{'F=' + str(F):10}", end="")
print()
print("  " + "-" * 80)

for name, D, C0 in elements:
    print(f"  {name:8} {D:6.3f}", end="")
    for F in [1.0, 0.8, 0.6, 0.4, 0.2, 0.1, 0.05]:
        CL = C0 * F**(D - 1)
        print(f"  {CL:8.1f}", end="")
    print()

# ============================================================
# 3. REE Patterns
# ============================================================
print()
print()
print("3. REE Chondrite-Normalized Patterns")
print("-" * 65)

# Chondrite values (CI, ppm) - McDonough & Sun 1995
ree_names = ["La", "Ce", "Pr", "Nd", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu"]
chondrite = [0.237, 0.613, 0.0928, 0.457, 0.148, 0.0563, 0.199, 0.0361, 0.246, 0.0546, 0.160, 0.0247, 0.161, 0.0246]

# Sample compositions (ppm)
samples = {
    "MORB": [3.2, 10.0, 1.7, 10.0, 3.5, 1.3, 5.0, 0.9, 6.0, 1.3, 3.8, 0.55, 3.5, 0.5],
    "OIB": [25.0, 55.0, 7.0, 30.0, 7.0, 2.3, 6.5, 1.0, 5.5, 1.1, 2.8, 0.4, 2.5, 0.36],
    "Granite": [35.0, 70.0, 8.0, 28.0, 5.0, 0.8, 4.0, 0.6, 3.5, 0.7, 2.0, 0.3, 2.0, 0.3],
    "UCC": [30.0, 64.0, 7.1, 26.0, 4.5, 0.88, 3.8, 0.64, 3.5, 0.80, 2.3, 0.33, 2.2, 0.32],
}

print(f"  {'REE':6} {'Chondrite':12}", end="")
for sname in samples:
    print(f"{sname + '_N':12}", end="")
print()
print("  " + "-" * 60)

for i, ree in enumerate(ree_names):
    print(f"  {ree:4} {chondrite[i]:10.4f}", end="")
    for sname, vals in samples.items():
        norm = vals[i] / chondrite[i]
        print(f"  {norm:10.1f}", end="")
    print()

# Eu anomaly
print()
for sname, vals in samples.items():
    Sm_N = vals[5-1] / chondrite[5-1]  # Sm index
    Eu_N = vals[5] / chondrite[5]      # Eu
    Gd_N = vals[6] / chondrite[6]      # Gd (index 6)
    Eu_star = np.sqrt(Sm_N * Gd_N)
    Eu_anom = Eu_N / Eu_star
    print(f"  {sname}: Eu/Eu* = {Eu_anom:.2f}")

# ============================================================
# 4. Goldschmidt Classification Table
# ============================================================
print()
print()
print("4. Goldschmidt Classification: Key Elements")
print("-" * 60)

gold_class = [
    ("Fe", "Siderophile", 26, 1811, "Core: 85% by mass"),
    ("Ni", "Siderophile", 28, 1728, "Core: 5% by mass"),
    ("Au", "Siderophile", 79, 1337, "Highly depleted in mantle"),
    ("Si", "Lithophile", 14, 1687, "Mantle: olivine, pyroxene"),
    ("Al", "Lithophile", 13, 933, "Crust enriched"),
    ("Ca", "Lithophile", 20, 1115, "Plagioclase, pyroxene"),
    ("K", "Lithophile", 19, 337, "Incompatible, crust-enriched"),
    ("Cu", "Chalcophile", 29, 1358, "Porphyry deposits"),
    ("Pb", "Chalcophile", 82, 601, "Sulfide ores"),
    ("Zn", "Chalcophile", 30, 693, "VMS deposits"),
]

print(f"  {'Element':10} {'Class':15} {'Z':6} {'Notes'}")
print("  " + "-" * 50)
for name, cls, Z, mp, notes in gold_class:
    print(f"  {name:8} {cls:13} {Z:4} {notes}")

# ============================================================
# 5. Geothermobarometry
# ============================================================
print()
print()
print("5. Garnet-Biotite Thermometry")
print("-" * 55)

# ln(Kd) = A/T + B  (Ferry & Spear 1978)
A_FS = -2109.0
B_FS = 0.782

print(f"  {'KD (Fe/Mg)':14} {'T (C)':12} {'Metamorphic grade'}")
print("  " + "-" * 45)

for KD in [1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15.0]:
    ln_KD = np.log(KD)
    T_K = A_FS / (ln_KD - B_FS)
    T_C = T_K - 273.15
    if T_C < 400:
        grade = "Greenschist"
    elif T_C < 550:
        grade = "Amphibolite"
    elif T_C < 700:
        grade = "Upper amphibolite"
    elif T_C < 850:
        grade = "Granulite"
    else:
        grade = "Ultra-high T"
    if T_C > 0:
        print(f"  {KD:12.1f} {T_C:10.0f} {grade}")

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

# Plot 1: REE patterns
ax = axes[0, 0]
x_ree = np.arange(len(ree_names))
colors_ree = {'MORB': 'blue', 'OIB': 'red', 'Granite': 'green', 'UCC': 'orange'}

for sname, vals in samples.items():
    norm_vals = [vals[i]/chondrite[i] for i in range(len(ree_names))]
    ax.semilogy(x_ree, norm_vals, 'o-', linewidth=2, markersize=5,
                color=colors_ree[sname], label=sname)

ax.set_xticks(x_ree)
ax.set_xticklabels(ree_names, fontsize=8)
ax.set_ylabel('Sample / Chondrite')
ax.set_title('Chondrite-Normalized REE Patterns')
ax.legend()
ax.grid(True, alpha=0.3, which='both')
ax.set_ylim(1, 500)

# Plot 2: Rayleigh fractionation
ax = axes[0, 1]
F_range = np.linspace(0.01, 1.0, 200)

for name, D, C0 in [("Ni (D=10)", 10, 1), ("Cr (D=7)", 7, 1),
                      ("Sr (D=2)", 2, 1), ("Ba (D=0.1)", 0.1, 1),
                      ("Rb (D=0.01)", 0.01, 1)]:
    CL = F_range**(D - 1)
    ax.semilogy(F_range, CL, linewidth=2, label=name)

ax.set_xlabel('Fraction of melt remaining (F)')
ax.set_ylabel('C_L / C_0')
ax.set_title('Rayleigh Fractional Crystallization')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3, which='both')
ax.set_xlim(0, 1)

# Plot 3: Phase diagram (simplified)
ax = axes[1, 0]
# Al2SiO5 phase diagram (simplified)
T_range_PT = np.linspace(200, 900, 200)

# Kyanite-Sillimanite boundary (positive slope)
P_ky_sill = 0.003 * T_range_PT + 1.5  # kbar, approximate
# Kyanite-Andalusite boundary
P_ky_and = 0.005 * (T_range_PT - 400) + 4.0
# Andalusite-Sillimanite boundary
T_and_sill = np.linspace(500, 900, 100)
P_and_sill = 0.002 * (T_and_sill - 500) + 3.5

# Simplified boundaries
T_trip = 500  # triple point T (C)
P_trip = 3.8  # triple point P (kbar)

T1 = np.linspace(200, T_trip, 50)
P1 = P_trip + 0.025 * (T_trip - T1)  # Ky-And

T2 = np.linspace(T_trip, 900, 50)
P2 = P_trip + 0.004 * (T2 - T_trip)  # Ky-Sill

T3 = np.linspace(T_trip, 900, 50)
P3 = P_trip - 0.002 * (T3 - T_trip)  # And-Sill

ax.plot(T1, P1, 'b-', linewidth=2)
ax.plot(T2, P2, 'r-', linewidth=2)
ax.plot(T3, P3, 'g-', linewidth=2)
ax.plot(T_trip, P_trip, 'ko', markersize=10, label='Triple point')

ax.text(350, 6, 'Kyanite', fontsize=14, fontweight='bold', color='blue')
ax.text(650, 6, 'Sillimanite', fontsize=14, fontweight='bold', color='red')
ax.text(600, 2, 'Andalusite', fontsize=14, fontweight='bold', color='green')

ax.set_xlabel('Temperature (C)')
ax.set_ylabel('Pressure (kbar)')
ax.set_title('Al2SiO5 Phase Diagram')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(200, 900)
ax.set_ylim(0, 10)

# Plot 4: Batch vs fractional melting
ax = axes[1, 1]
F_melt = np.linspace(0.001, 0.5, 200)

for D0, label in [(0.01, 'D=0.01 (very incompat.)'),
                    (0.1, 'D=0.1 (incompat.)'),
                    (1.0, 'D=1.0 (neutral)'),
                    (5.0, 'D=5.0 (compat.)')]:
    # Batch melting
    CL_batch = 1.0 / (D0 + F_melt * (1 - D0))
    # Fractional melting (aggregated)
    CL_frac = (1.0 / F_melt) * (1 - (1 - F_melt)**(1.0/D0)) if D0 != 0 else np.ones_like(F_melt)

ax.plot(F_melt*100, CL_batch, linewidth=2, label=f'Batch {label}')

ax.set_xlabel('Melt fraction F (%)')
ax.set_ylabel('C_L / C_0')
ax.set_title('Batch Partial Melting')
ax.legend(fontsize=7)
ax.grid(True, alpha=0.3)
ax.set_ylim(0.01, 100)
ax.set_yscale('log')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
plt.close()
print()
print()
print("6. Spider Diagram: Mantle-Normalized Trace Elements")
print("-" * 65)

# Primitive mantle values (McDonough & Sun 1995)
spider_elements = ["Rb", "Ba", "Th", "U", "Nb", "K", "La", "Ce", "Sr", "Nd", "Sm", "Zr", "Ti", "Y", "Yb"]
PM_values = [0.6, 6.6, 0.08, 0.02, 0.66, 240, 0.65, 1.68, 19.9, 1.25, 0.41, 10.5, 1205, 4.3, 0.44]

# MORB and IAB compositions
morb_vals = [0.5, 10, 0.1, 0.04, 2.0, 500, 3.0, 10, 100, 10, 3.5, 80, 8000, 30, 3.5]
iab_vals = [15, 200, 1.5, 0.5, 1.0, 8000, 8, 20, 400, 12, 3.0, 60, 5000, 20, 2.0]

print(f"  {'Element':10} {'PM':10} {'MORB_N':12} {'IAB_N':12} {'IAB anomaly'}")
print("  " + "-" * 50)

for i, el in enumerate(spider_elements):
    morb_n = morb_vals[i] / PM_values[i]
    iab_n = iab_vals[i] / PM_values[i]
    anomaly = ""
    if el in ["Nb", "Ti"] and iab_n < morb_n * 0.5:
        anomaly = "NEGATIVE"
    elif el in ["Rb", "Ba", "K", "Sr"] and iab_n > morb_n * 3:
        anomaly = "POSITIVE (fluid)"
    print(f"  {el:8} {PM_values[i]:8.2f} {morb_n:10.1f} {iab_n:10.1f} {anomaly}")

# Phase rule examples
print()
print()
print("7. Gibbs Phase Rule Examples")
print("-" * 55)

print(f"  {'System':25} {'C':6} {'P':6} {'F':6} {'Meaning'}")
print("  " + "-" * 50)

phase_examples = [
    ("H2O at triple point", 1, 3, 0, "Fixed T and P"),
    ("H2O ice+water", 1, 2, 1, "Along melting curve"),
    ("H2O liquid only", 1, 1, 2, "T and P free"),
    ("Al2SiO5 triple point", 1, 3, 0, "Ky-And-Sill point"),
    ("Granite solidus", 3, 4, 1, "Univariant curve"),
    ("Basalt (4 phases)", 4, 4, 2, "T and P free"),
]

for name, C, P, F in phase_examples:
    meaning = phase_examples[[x[0] for x in phase_examples].index(name)][4]
    print(f"  {name:23} {C:4} {P:4} {F:4} {meaning}")

print()
print()
print("[Chart saved: REE patterns, fractionation, Al2SiO5 diagram, partial melting]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Geothermobarometry and Mineral Stability Calculations

Python

script.py69 lines

import numpy as np

# ============================================================
# GEOTHERMOBAROMETRY AND MINERAL STABILITY
# ============================================================
print("GEOTHERMOBAROMETRY CALCULATIONS")
print("=" * 70)

R = 8.314

# ============================================================
# 1. Diamond-Graphite Phase Boundary
# ============================================================
print()
print("1. Diamond-Graphite Phase Boundary (P vs T)")
print("-" * 55)

# dG = dH - T*dS + P*dV = 0 at equilibrium
dH_dg = 1900.0   # J/mol (graphite -> diamond)
dS_dg = -3.36    # J/(mol.K)
dV_dg = -1.88e-6 # m^3/mol

print(f"  dH = {dH_dg} J/mol, dS = {dS_dg} J/(mol.K)")
print(f"  dV = {dV_dg*1e6:.2f} cm^3/mol")
print(f"  dP/dT = dS/dV = {dS_dg/dV_dg/1e9:.2f} GPa/K")
print()
print(f"  {'T (C)':12} {'P_eq (GPa)':14} {'P_eq (kbar)':14} {'Diamond stable above'}")
print("  " + "-" * 50)

for T_C in [500, 700, 900, 1000, 1200, 1400, 1600, 1800, 2000]:
    T_K = T_C + 273.15
    # P = (T*dS - dH) / dV  (from dG = 0)
    P_eq = (T_K * dS_dg - dH_dg) / dV_dg
    P_GPa = P_eq / 1e9
    P_kbar = P_GPa * 10
    depth_km = P_GPa / 0.033  # approx 33 km per GPa
    print(f"  {T_C:10} {P_GPa:12.2f} {P_kbar:12.1f} {depth_km:.0f} km depth")

# ============================================================
# 2. Garnet-Biotite Fe-Mg Exchange Thermometer
# ============================================================
print()
print()
print("2. Garnet-Biotite Fe-Mg Thermometer (Ferry & Spear 1978)")
print("-" * 60)

# ln(KD) = -2109/T + 0.782
A_FS = -2109.0
B_FS = 0.782

print(f"  {'Fe/Mg_Grt':14} {'Fe/Mg_Bt':14} {'KD':10} {'ln(KD)':10} {'T (C)':10}")
print("  " + "-" * 55)

sample_pairs = [
    (2.5, 0.8), (3.0, 1.0), (4.0, 1.5), (5.0, 2.0),
    (6.0, 2.5), (8.0, 3.0), (10.0, 3.5), (15.0, 5.0),
]

for FeMg_grt, FeMg_bt in sample_pairs:
    KD = FeMg_grt / FeMg_bt
    ln_KD = np.log(KD)
    T_K = A_FS / (ln_KD - B_FS)
    T_C = T_K - 273.15
    if T_C > 0 and T_C < 1200:
        print(f"  {FeMg_grt:12.1f} {FeMg_bt:12.1f} {KD:8.2f} {ln_KD:8.3f} {T_C:8.0f}")

print()
print("Complete: geothermobarometry calculations finished.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Paleoclimatology Next: Isotope Geology →