Part III: Climate & Oceans | Chapter 4

Paleoclimatology

Ice cores, oxygen isotope proxies, Milankovitch cycles, and deep-time climate events

3.1 Oxygen Isotope Paleothermometry

The ratio of $^{18}$O to $^{16}$O in geological archives (ice cores, marine carbonates, speleothems) records past temperatures. The isotopic composition is expressed in delta notation:

$$\delta^{18}\text{O} = \left(\frac{(^{18}\text{O}/^{16}\text{O})_{\text{sample}}}{(^{18}\text{O}/^{16}\text{O})_{\text{standard}}} - 1\right) \times 1000 \text{ \textperthousand}$$

Derivation 1: Temperature Dependence of Fractionation

Harold Urey (1947) showed that the equilibrium fractionation between calcite and water is temperature-dependent. The fractionation factor $\alpha$ relates to the equilibrium constant of the isotope exchange reaction:

$$\alpha = \frac{(^{18}\text{O}/^{16}\text{O})_{\text{CaCO}_3}}{(^{18}\text{O}/^{16}\text{O})_{\text{H}_2\text{O}}} \approx 1 + \frac{A}{T^2}$$

The paleotemperature equation of Epstein et al. (1953) gives:

$$T(°\text{C}) = 16.5 - 4.3(\delta^{18}\text{O}_c - \delta^{18}\text{O}_w) + 0.14(\delta^{18}\text{O}_c - \delta^{18}\text{O}_w)^2$$

where $\delta^{18}\text{O}_c$ is the carbonate value (PDB standard) and$\delta^{18}\text{O}_w$ is the seawater value (SMOW standard). A 1 per mille increase in $\delta^{18}\text{O}_c$ corresponds to approximately 4$°$C of cooling.

Cesare Emiliani (1955) applied this technique to deep-sea foraminifera, discovering the cyclic pattern of glacial-interglacial oscillations recorded in marine sediments. Nicholas Shackleton later showed that much of the signal reflects ice volume (ice sheets preferentially store $^{16}$O, enriching seawater in $^{18}$O) rather than temperature alone.

3.2 Milankovitch Cycles

Derivation 2: Insolation from Orbital Parameters

Milutin Milankovitch (1941) calculated that variations in Earth's orbital parameters modulate the seasonal and latitudinal distribution of solar radiation, pacing the ice ages. The three parameters are:

Eccentricity ($e$): Shape of orbit, periods ~100 kyr and ~400 kyr, range 0.005-0.058
Obliquity ($\varepsilon$): Axial tilt, period ~41 kyr, range 22.1$°$-24.5$°$
Precession: $e\sin\varpi$ where $\varpi$ is the longitude of perihelion, periods ~19 and ~23 kyr

The mean daily insolation at the top of the atmosphere at latitude $\phi$ on a given day depends on the solar declination $\delta$ and the hour angle of sunset $H_0$:

$$Q = \frac{S_0}{\pi}\left(\frac{a}{r}\right)^2 (H_0 \sin\phi \sin\delta + \cos\phi \cos\delta \sin H_0)$$

where $r$ is the Earth-Sun distance and $H_0 = \arccos(-\tan\phi\tan\delta)$. The critical quantity is summer insolation at 65$°$N: when it is low, winter snow survives through summer and ice sheets grow.

Historical Context

James Croll (1864) first proposed that orbital variations drive ice ages. Milankovitch (1920-1941) computed detailed insolation curves. The theory was confirmed by Hays, Imbrie, and Shackleton (1976), who found the 100 kyr, 41 kyr, and 23 kyr periodicities in deep-sea sediment cores, matching the orbital periods. This paper, often called the "Pacemaker of the Ice Ages," is a landmark in paleoclimatology.

3.3 Ice Core Records

Derivation 3: Ice Core Depth-Age Relationship

The depth-age relationship in an ice core depends on the accumulation rate and ice flow dynamics. For a simple model with constant accumulation rate $b$ and uniform vertical strain rate (Nye model), the depth-age relation is:

$$z(t) = H\left(1 - e^{-bt/H}\right)$$

where $H$ is the ice sheet thickness and $t$ is age. The inverse gives:

$$t(z) = -\frac{H}{b}\ln\left(1 - \frac{z}{H}\right)$$

This means annual layers thin exponentially with depth: near the surface, layers are close to $b$ thick, but at great depth they become vanishingly thin, ultimately limiting temporal resolution. The EPICA Dome C core in Antarctica reaches 800,000 years at 3,270 m depth.

Ice cores preserve not only $\delta^{18}$O (temperature) but also trapped air bubbles recording past atmospheric CO$_2$ and CH$_4$ concentrations, dust, volcanic aerosols, and cosmogenic isotopes. The EPICA record shows CO$_2$ varied between ~180 ppm (glacials) and ~280 ppm (interglacials) over the past 800 kyr, closely tracking Antarctic temperature.

3.4 Snowball Earth and the PETM

Derivation 4: CO$_2$ Threshold for Snowball Escape

During the Neoproterozoic (717-635 Ma), geological evidence suggests Earth was completely or nearly completely ice-covered. Escape from a snowball state requires massive CO$_2$buildup from volcanic outgassing (with no silicate weathering sink on ice-covered continents).

The required CO$_2$ for deglaciation can be estimated from the energy balance. With high albedo ($\alpha \approx 0.6$), the absorbed solar flux is:

$$\frac{(1-0.6) \times 1361}{4} = 136 \text{ W/m}^2$$

Using the logarithmic radiative forcing law, the CO$_2$ needed to raise temperature above the melting threshold (~30 K of greenhouse warming above effective temperature):

$$\Delta F \approx 5.35 \times \ln\left(\frac{C}{280}\right) \approx 50 \text{ W/m}^2$$

This requires $C \approx 30,000$-100,000 ppm (~10-30%), achievable after ~10 Myr of volcanic outgassing at typical rates. Joseph Kirschvink (1992) proposed the snowball Earth hypothesis; Paul Hoffman et al. (1998) provided geological evidence including cap carbonates deposited immediately after glaciation.

The Paleocene-Eocene Thermal Maximum (PETM)

At 56 Ma, a massive carbon release (~3,000-10,000 Gt C) caused 5-8$°$C of global warming in ~20,000 years. The carbon isotope excursion ($\delta^{13}$C dropped by ~3 per mille) indicates a light-carbon source, possibly methane hydrate dissociation, volcanic CO$_2$, or thermogenic methane. Deep ocean temperatures rose to ~20$°$C, and the event lasted ~170,000 years before carbon cycle feedbacks (silicate weathering) restored equilibrium.

The Great Oxidation Event

At ~2.4 Ga, atmospheric oxygen rose from trace levels ($< 10^{-5}$ PAL) to perhaps 1-10% of present levels. Evidence includes the disappearance of detrital uraninite and pyrite (oxidized at the surface), the appearance of red beds (oxidized iron), and the loss of mass-independent sulfur isotope fractionation (which requires UV photolysis of SO$_2$ in an ozone-free atmosphere). The rise of oxygen was likely driven by cyanobacterial photosynthesis, but the ~300 Myr delay between the oldest cyanobacterial fossils (~2.7 Ga) and the GOE suggests that reducing sinks (volcanic gases, reduced iron) had to be overcome before oxygen could accumulate.

Cenozoic Cooling Trend

The Cenozoic Era (66 Ma to present) records a long-term cooling trend from greenhouse conditions (no ice caps, tropical forests at high latitudes) to the current icehouse. Key transitions include: the Eocene-Oligocene boundary (~34 Ma) when Antarctica glaciated, coinciding with the opening of Drake Passage and thermal isolation of Antarctica; the middle Miocene (~14 Ma) expansion of the East Antarctic Ice Sheet; and the intensification of Northern Hemisphere glaciation at ~2.7 Ma. The long-term cooling is attributed primarily to declining atmospheric CO$_2$ from enhanced silicate weathering driven by Himalayan uplift and changes in ocean circulation.

3.5 Marine Isotope Stages and Climate Proxies

Derivation 5: Spectral Analysis of Climate Records

To identify periodicities in climate records, we use spectral analysis. The power spectrum of a discrete time series $x_n$ sampled at interval $\Delta t$ is:

$$P(f_k) = \frac{\Delta t}{N}\left|\sum_{n=0}^{N-1} x_n e^{-2\pi i f_k n \Delta t}\right|^2$$

where $f_k = k/(N\Delta t)$. Applied to benthic $\delta^{18}$O records, three dominant peaks emerge at periods of ~100, ~41, and ~23 kyr, corresponding to eccentricity, obliquity, and precession.

The marine isotope stage (MIS) system numbers glacial-interglacial cycles back through time. Odd numbers are interglacials (MIS 1 = Holocene, MIS 5 = last interglacial at ~125 ka), even numbers are glacials (MIS 2 = Last Glacial Maximum at ~21 ka). The LR04 benthic stack (Lisiecki and Raymo, 2005) provides the standard reference curve extending to 5.3 Ma.

A major puzzle is the Mid-Pleistocene Transition (~1 Ma): ice age cycles switched from 41 kyr (obliquity-dominated) to ~100 kyr periodicity, even though eccentricity forcing is the weakest of the three. Proposed explanations include changes in CO$_2$, erosion of regolith exposing bedrock, and nonlinear ice sheet dynamics.

3.6 Additional Climate Proxies

Beyond oxygen isotopes, paleoclimatologists use a wide array of proxies to reconstruct past environments. Each proxy has strengths, limitations, and applicable timescales:

Deuterium in Ice Cores

The deuterium-hydrogen ratio ($\delta$D) in ice cores provides an independent temperature proxy. The relationship is approximately linear:

$$\delta\text{D} \approx 8 \times \delta^{18}\text{O} + 10$$

This is the Global Meteoric Water Line (Craig, 1961). The deuterium excess$d = \delta\text{D} - 8\delta^{18}\text{O}$ provides information about evaporation conditions at the moisture source region. Antarctic ice cores show temperature variations of ~8-10$°$C between glacial and interglacial periods.

Mg/Ca Paleothermometry

The Mg/Ca ratio in foraminiferal calcite is temperature-dependent because Mg substitution into the CaCO$_3$ lattice is endothermic. The relationship is exponential:

$$\text{Mg/Ca} = B \exp(A \cdot T)$$

where $A \approx 0.09$ per $°$C and $B$ depends on species. This proxy isolates the temperature signal from the ice volume signal that complicates$\delta^{18}$O interpretation. Elderfield and Ganssen (2000) used Mg/Ca to separate temperature and ice volume contributions to the $\delta^{18}$O record.

Alkenone Thermometry

Marine phytoplankton (coccolithophores) produce long-chain unsaturated ketones (alkenones) whose degree of unsaturation depends on growth temperature. The U$^{K'}_{37}$ index, defined as:

$$U^{K'}_{37} = \frac{C_{37:2}}{C_{37:2} + C_{37:3}}$$

correlates linearly with sea surface temperature: $T(°C) = (U^{K'}_{37} - 0.044)/0.033$(Prahl et al., 1988). Alkenones are preserved in sediments for tens of millions of years, providing SST records through the Cenozoic.

Tree Ring Dendroclimatology

Tree ring widths and densities record annual climate variations over the past ~12,000 years (oldest living trees: bristlecone pines, ~5000 years; subfossil wood extends further). Ring width depends on temperature, moisture, and growing season length. Standardization removes age-related growth trends, and cross-dating ensures accurate chronology. The Medieval Climate Anomaly (~950-1250 CE) and Little Ice Age (~1450-1850 CE) are well documented in tree ring records.

Speleothem Records

Cave stalagmites provide high-resolution paleoclimate records through their$\delta^{18}$O and $\delta^{13}$C composition, growth rate, and trace element content. U-Th dating gives precise chronologies back to ~500,000 years. The$\delta^{18}$O of speleothem calcite records both cave temperature and the isotopic composition of drip water (which reflects rainfall amount and source). In the Asian monsoon region, speleothem records from Hulu, Dongge, and Sanbao caves in China have provided the most detailed records of monsoon variability, revealing strong orbital (precession) pacing and abrupt millennial-scale events (Dansgaard-Oeschger oscillations) that correlate precisely with Greenland ice core records.

Pollen and Biomarker Proxies

Fossil pollen preserved in lake sediments and peat bogs records vegetation change and, by inference, climate. Quantitative transfer functions relate modern pollen assemblages to climate parameters (temperature, precipitation), then apply these relationships to fossil assemblages. The European Pollen Database contains records from thousands of sites, documenting the northward migration of forests after the last glaciation at rates of ~200-500 m/yr.

Biomarkers (molecular fossils) provide additional constraints. The TEX$_{86}$proxy uses the distribution of glycerol dialkyl glycerol tetraethers (GDGTs) produced by marine archaea to reconstruct sea surface temperature. The BIT index distinguishes terrestrial from marine organic matter input to sediments. These organic proxies complement inorganic methods and extend temperature reconstructions through the Cenozoic and beyond.

Computational Lab: Paleoclimatology

Milankovitch Cycles, Isotope Records, and Ice Core Analysis

Python

script.py334 lines

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

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

# ============================================================
# 1. Milankovitch Orbital Parameters
# ============================================================
print()
print("1. Milankovitch Orbital Parameters Over Time")
print("-" * 65)

# Simplified orbital parameter oscillations
t_kyr = np.linspace(0, 800, 8000)  # 0 to 800 ka

# Eccentricity (100 kyr and 400 kyr periods)
ecc = 0.028 + 0.02 * np.sin(2*np.pi*t_kyr/100) + 0.012 * np.sin(2*np.pi*t_kyr/413)

# Obliquity (41 kyr period)
obliq = 23.3 + 1.2 * np.sin(2*np.pi*t_kyr/41)

# Precession (23 kyr and 19 kyr periods)
prec = 0.02 * np.sin(2*np.pi*t_kyr/23) + 0.015 * np.sin(2*np.pi*t_kyr/19)

print(f"  {'Time (ka)':12} {'Eccentricity':14} {'Obliquity (deg)':18} {'Precession index'}")
print("  " + "-" * 55)

for idx in range(0, len(t_kyr), 1000):
    t = t_kyr[idx]
    print(f"  {t:10.0f} {ecc[idx]:12.4f} {obliq[idx]:16.2f} {prec[idx]:12.4f}")

# ============================================================
# 2. Insolation Calculation
# ============================================================
print()
print()
print("2. Summer Insolation at 65N")
print("-" * 55)

S0 = 1361.0
lat = 65.0  # degrees

# Simplified: Q_summer ~ S0/pi * (1+e*cos(w)) * f(obliquity, latitude)
# The critical parameter is a combination of precession and eccentricity

print(f"  {'Time (ka)':12} {'Q_65N (W/m2)':16} {'Anomaly':14} {'Notes'}")
print("  " + "-" * 55)

# Compute insolation (simplified model)
Q_65N = np.zeros_like(t_kyr)
for i in range(len(t_kyr)):
    eps = np.radians(obliq[i])
    # Summer solstice declination = obliquity
    delta = eps
    phi = np.radians(lat)
    # Hour angle at sunset
    cos_H0 = -np.tan(phi) * np.tan(delta)
    cos_H0 = np.clip(cos_H0, -1, 1)
    H0 = np.arccos(cos_H0)
    # Daily insolation
    Q_65N[i] = S0/np.pi * (1 + ecc[i]) * (H0*np.sin(phi)*np.sin(delta) +
               np.cos(phi)*np.cos(delta)*np.sin(H0))

Q_mean = np.mean(Q_65N)
for idx in range(0, len(t_kyr), 800):
    t = t_kyr[idx]
    anomaly = Q_65N[idx] - Q_mean
    note = ""
    if anomaly < -15:
        note = "glaciation favorable"
    elif anomaly > 15:
        note = "deglaciation favorable"
    print(f"  {t:10.0f} {Q_65N[idx]:14.1f} {anomaly:12.1f} {note}")

# ============================================================
# 3. Synthetic d18O Record
# ============================================================
print()
print()
print("3. Synthetic Benthic d18O Record (simplified LR04-like)")
print("-" * 60)

# Combine orbital signals with nonlinear ice sheet response
# d18O roughly follows insolation with slow ice growth and fast decay
d18O = np.zeros_like(t_kyr)
ice_volume = 0.0
dt = t_kyr[1] - t_kyr[0]

for i in range(1, len(t_kyr)):
    # Insolation anomaly drives ice changes
    Q_anom = (Q_65N[i] - Q_mean) / Q_mean
    if Q_anom < 0:
        # Slow ice growth
        dice = -0.3 * Q_anom * dt  # ice grows when insolation low
    else:
        # Fast ice decay (asymmetry -> sawtooth)
        dice = -0.8 * Q_anom * dt
    ice_volume += dice
    ice_volume = max(0, ice_volume)
    d18O[i] = 3.0 + 0.5 * ice_volume  # baseline + ice volume effect

# Smooth
kernel_size = 50
kernel = np.ones(kernel_size) / kernel_size
d18O_smooth = np.convolve(d18O, kernel, mode='same')

print(f"  {'Time (ka)':12} {'d18O (permil)':16} {'Ice vol (arb)':14} {'MIS stage'}")
print("  " + "-" * 50)

mis_map = {0: "1 (Holocene)", 21: "2 (LGM)", 125: "5e (Eemian)",
           200: "7", 330: "9", 410: "11", 500: "13"}

for t_target in [0, 10, 21, 50, 80, 125, 200, 330, 410, 500, 630, 780]:
    idx = np.argmin(abs(t_kyr - t_target))
    mis = mis_map.get(t_target, "")
    print(f"  {t_target:10} {d18O_smooth[idx]:14.2f} {ice_volume:12.2f} {mis}")

# ============================================================
# 4. Ice Core Depth-Age Model
# ============================================================
print()
print()
print("4. Ice Core Depth-Age Relationship (Nye Model)")
print("-" * 55)

H = 3200.0   # ice thickness (m)
b = 0.03     # accumulation rate (m/yr ice equivalent)

print(f"  Ice thickness: {H} m, Accumulation: {b*1000:.0f} mm/yr")
print()
print(f"  {'Depth (m)':14} {'Age (kyr)':14} {'Layer thick (mm)':18} {'Resolution'}")
print("  " + "-" * 55)

for z in [0, 100, 500, 1000, 1500, 2000, 2500, 2800, 3000, 3100, 3150, 3190]:
    if z < H:
        age = -H/b * np.log(1 - z/H) / 1000  # kyr
        layer = b * (1 - z/H) * 1000  # mm/yr
        if layer > 10:
            res = "annual"
        elif layer > 1:
            res = "decadal"
        elif layer > 0.1:
            res = "centennial"
        else:
            res = "millennial"
        print(f"  {z:12} {age:12.1f} {layer:16.2f} {res}")

# ============================================================
# 5. Paleotemperature Equation
# ============================================================
print()
print()
print("5. Paleotemperature from Carbonate d18O")
print("-" * 60)

print(f"  {'d18O_c (permil)':18} {'d18O_w assumed':16} {'T (C)':10} {'Environment'}")
print("  " + "-" * 55)

environments = [
    (-1.0, 0.0, "Warm tropical surface"),
    (0.0, 0.0, "Subtropical"),
    (1.0, 0.0, "Temperate"),
    (2.0, 0.0, "Cool/deep"),
    (3.0, 0.0, "Cold deep water"),
    (4.0, 0.0, "Near freezing"),
    (0.0, -1.0, "Glacial tropical (ice effect)"),
    (3.0, 1.0, "Glacial deep (ice effect)"),
]

for d18Oc, d18Ow, env in environments:
    diff = d18Oc - d18Ow
    T = 16.5 - 4.3 * diff + 0.14 * diff**2
    print(f"  {d18Oc:16.1f} {d18Ow:14.1f} {T:8.1f} {env}")

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

# Plot 1: Orbital parameters
ax1 = axes[0, 0]
ax1b = ax1.twinx()
ax1.plot(t_kyr, ecc, 'b-', linewidth=1, alpha=0.8, label='Eccentricity')
ax1b.plot(t_kyr, obliq, 'r-', linewidth=1, alpha=0.8, label='Obliquity (deg)')
ax1.set_xlabel('Time (ka)')
ax1.set_ylabel('Eccentricity', color='blue')
ax1b.set_ylabel('Obliquity (deg)', color='red')
ax1.set_title('Orbital Parameters (last 800 kyr)')
ax1.legend(loc='upper left', fontsize=8)
ax1b.legend(loc='upper right', fontsize=8)
ax1.invert_xaxis()
ax1.grid(True, alpha=0.3)

# Plot 2: Insolation at 65N
ax = axes[0, 1]
ax.plot(t_kyr, Q_65N, 'darkorange', linewidth=1)
ax.fill_between(t_kyr, Q_mean, Q_65N, where=Q_65N > Q_mean,
                alpha=0.3, color='red', label='Above mean')
ax.fill_between(t_kyr, Q_mean, Q_65N, where=Q_65N < Q_mean,
                alpha=0.3, color='blue', label='Below mean')
ax.axhline(y=Q_mean, color='gray', linestyle='--', linewidth=1)
ax.set_xlabel('Time (ka)')
ax.set_ylabel('Insolation (W/m^2)')
ax.set_title('Summer Insolation at 65N')
ax.legend(fontsize=9)
ax.invert_xaxis()
ax.grid(True, alpha=0.3)

# Plot 3: Synthetic d18O record
ax = axes[1, 0]
ax.plot(t_kyr, d18O_smooth, 'b-', linewidth=1)
ax.fill_between(t_kyr, 3.0, d18O_smooth, where=d18O_smooth > 3.5,
                alpha=0.3, color='lightblue', label='Glacial')
ax.fill_between(t_kyr, 3.0, d18O_smooth, where=d18O_smooth < 3.5,
                alpha=0.3, color='lightyellow', label='Interglacial')
ax.set_xlabel('Time (ka)')
ax.set_ylabel('d18O (permil)')
ax.set_title('Synthetic Benthic d18O Record')
ax.invert_yaxis()  # Convention: more positive = colder/more ice
ax.invert_xaxis()
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# Plot 4: Power spectrum
ax = axes[1, 1]
# Compute FFT of d18O record
N = len(d18O_smooth)
dt_spec = t_kyr[1] - t_kyr[0]  # kyr
freqs = np.fft.rfftfreq(N, d=dt_spec)
fft_vals = np.fft.rfft(d18O_smooth - np.mean(d18O_smooth))
power = np.abs(fft_vals)**2

# Convert frequency to period
mask = freqs > 0
periods = 1.0 / freqs[mask]
power_plot = power[mask]

# Only plot reasonable period range
range_mask = (periods > 10) & (periods < 500)
ax.semilogy(periods[range_mask], power_plot[range_mask], 'b-', linewidth=1)

# Mark Milankovitch periods
for period, label, color in [(100, '100 kyr', 'red'), (41, '41 kyr', 'green'),
                               (23, '23 kyr', 'orange'), (19, '19 kyr', 'purple')]:
    ax.axvline(x=period, color=color, linestyle='--', alpha=0.7, label=label)

ax.set_xlabel('Period (kyr)')
ax.set_ylabel('Spectral Power')
ax.set_title('Power Spectrum of d18O Record')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3, which='both')
ax.set_xlim(10, 500)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
plt.close()
print()
print()
print("6. Mg/Ca Paleothermometry")
print("-" * 55)
print(f"  {'Mg/Ca (mmol/mol)':20} {'T (C)':12} {'Environment'}")
print("  " + "-" * 45)

A_mgca = 0.09
B_mgca = 0.38

for mgca in [0.5, 0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0]:
    T_mgca = np.log(mgca / B_mgca) / A_mgca
    if T_mgca < 5:
        env = "Deep/polar"
    elif T_mgca < 15:
        env = "Temperate"
    elif T_mgca < 25:
        env = "Subtropical"
    else:
        env = "Tropical"
    print(f"  {mgca:18.1f} {T_mgca:10.1f} {env}")

# ============================================================
# 7. Alkenone SST Proxy
# ============================================================
print()
print()
print("7. Alkenone U37' Sea Surface Temperature")
print("-" * 55)
print(f"  {'UK37':12} {'SST (C)':12} {'Region'}")
print("  " + "-" * 35)

for uk in [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]:
    T_alk = (uk - 0.044) / 0.033
    if T_alk < 10:
        region = "Subpolar"
    elif T_alk < 20:
        region = "Temperate"
    elif T_alk < 28:
        region = "Subtropical"
    else:
        region = "Tropical"
    print(f"  {uk:10.2f} {T_alk:10.1f} {region}")

# ============================================================
# 8. CO2 and Temperature Correlation
# ============================================================
print()
print()
print("8. Ice Core CO2-Temperature Relationship")
print("-" * 60)

# Synthetic CO2 correlated with d18O
CO2_glacial = 180
CO2_interglacial = 280
CO2_record = CO2_interglacial - (d18O_smooth - 3.0) / (np.max(d18O_smooth) - 3.0) * (CO2_interglacial - CO2_glacial)
CO2_record = np.clip(CO2_record, 170, 300)

print(f"  {'Time (ka)':12} {'CO2 (ppm)':14} {'dT (C, approx)':16} {'Regime'}")
print("  " + "-" * 55)

for t_target in [0, 10, 21, 50, 80, 125, 200, 330, 410, 500, 630, 780]:
    idx = np.argmin(abs(t_kyr - t_target))
    co2_val = CO2_record[idx]
    dT_approx = (co2_val - 230) * 0.08  # rough scaling
    regime = "Interglacial" if co2_val > 250 else ("Glacial" if co2_val < 200 else "Transition")
    print(f"  {t_target:10} {co2_val:12.0f} {dT_approx:14.1f} {regime}")

print()
print()
print("[Chart saved: orbital parameters, insolation, d18O record, power spectrum]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Paleotemperature Proxies and Climate History

Python

script.py80 lines

import numpy as np

print("PALEOTEMPERATURE PROXY CALCULATIONS")
print("=" * 70)

print()
print("1. Multi-Proxy Temperature Estimates")
print("-" * 65)

print(f"  {'Proxy':20} {'Range':20} {'Resolution':16} {'Accuracy'}")
print("  " + "-" * 60)

proxies = [
    ("d18O (foram)", "0 - 200 Ma", "100-1000 yr", "+/- 1-2 C"),
    ("Mg/Ca (foram)", "0 - 100 Ma", "100-1000 yr", "+/- 1 C"),
    ("Alkenone UK37", "0 - 50 Ma", "100-1000 yr", "+/- 1.5 C"),
    ("TEX86 (GDGT)", "0 - 200 Ma", "1000-10000 yr", "+/- 2 C"),
    ("d18O (ice core)", "0 - 0.8 Ma", "1-100 yr", "+/- 0.5 C"),
    ("Tree rings", "0 - 12 ka", "1 yr", "+/- 0.5 C"),
    ("Speleothems", "0 - 0.5 Ma", "1-100 yr", "+/- 1 C"),
    ("Clumped isotopes", "0 - 500 Ma", "10000 yr", "+/- 2 C"),
]

for name, rng, res, acc in proxies:
    print(f"  {name:18} {rng:18} {res:14} {acc}")

# Cenozoic temperature history
print()
print()
print("2. Key Cenozoic Climate Events")
print("-" * 65)

events = [
    (0, 15, 288, "Present day (pre-industrial)"),
    (0.02, -6, 282, "Last Glacial Maximum"),
    (0.125, 17, 289, "Last Interglacial (Eemian)"),
    (3, 17, 290, "Mid-Pliocene Warm Period"),
    (14, 20, 293, "Mid-Miocene Optimum"),
    (34, 23, 296, "Pre-Eocene-Oligocene glaciation"),
    (40, 25, 298, "Middle Eocene"),
    (50, 28, 301, "Early Eocene Climatic Optimum"),
    (56, 30, 303, "PETM peak"),
    (65, 22, 295, "End-Cretaceous"),
]

print(f"  {'Age (Ma)':12} {'Global T (C)':14} {'CO2 (ppm est.)':18} {'Event'}")
print("  " + "-" * 70)

for age, T, T_K, event in events:
    # Very rough CO2 estimate from temperature
    if age < 1:
        co2 = 280
    elif age < 5:
        co2 = 400
    elif age < 20:
        co2 = 500
    elif age < 40:
        co2 = 800
    elif age < 55:
        co2 = 1200
    elif age < 57:
        co2 = 2000
    else:
        co2 = 600
    print(f"  {age:10} {T:12} {co2:16} {event}")

# Ice volume estimation
print()
print()
print("3. d18O to Ice Volume and Temperature Decomposition")
print("-" * 60)
print(f"  Total benthic d18O change (glacial-interglacial): ~1.8 permil")
print(f"  Temperature component: ~0.8 permil (~3 C deep ocean cooling)")
print(f"  Ice volume component: ~1.0 permil")
print(f"  Sea level equivalent: ~120 m (from coral reef evidence)")
print(f"  -> 1 permil d18O_sw ~ 120 m sea level change")

print()
print("Complete: paleoclimate proxy analysis finished.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Ocean Circulation Next: Geochemistry →