Part III: Climate & Oceans | Chapter 1

Climate Science

Energy balance, the greenhouse effect, climate sensitivity, and feedback mechanisms

3.1 Earth's Energy Balance

Earth's climate is fundamentally governed by the balance between absorbed solar radiation and emitted thermal (longwave) radiation. The solar constant is $S_0 \approx 1361$ W/m$^2$, but Earth intercepts this over a disk of area $\pi R^2$ and distributes it over a sphere of area $4\pi R^2$.

Derivation 1: Effective Temperature of Earth

The absorbed solar power equals:

$$P_{\text{abs}} = (1 - \alpha) \frac{S_0}{4} \cdot 4\pi R^2 = (1 - \alpha) S_0 \pi R^2$$

where $\alpha \approx 0.30$ is Earth's planetary albedo. The emitted thermal radiation follows the Stefan-Boltzmann law:

$$P_{\text{emit}} = \sigma T_e^4 \cdot 4\pi R^2$$

Setting $P_{\text{abs}} = P_{\text{emit}}$ and solving for the effective temperature:

$$T_e = \left[\frac{(1 - \alpha) S_0}{4\sigma}\right]^{1/4}$$

With $\alpha = 0.30$ and $\sigma = 5.67 \times 10^{-8}$ W/(m$^2$$\cdot$K$^4$):

$$T_e = \left[\frac{0.70 \times 1361}{4 \times 5.67 \times 10^{-8}}\right]^{1/4} = 255 \text{ K} = -18°\text{C}$$

Earth's actual mean surface temperature is about 288 K (15$°$C). The 33 K difference is the greenhouse effect, first quantified by Joseph Fourier (1824) and explained physically by John Tyndall (1859) and Svante Arrhenius (1896).

3.2 The Greenhouse Effect

Derivation 2: Single-Layer Atmospheric Model

The simplest greenhouse model treats the atmosphere as a single layer that is transparent to shortwave solar radiation but absorbs a fraction $\epsilon$ of the longwave radiation emitted by the surface. Energy balance at the surface:

$$\frac{(1-\alpha)S_0}{4} + \epsilon \sigma T_a^4 = \sigma T_s^4$$

Energy balance for the atmospheric layer (absorbs $\epsilon$ from surface, emits $\epsilon$both up and down):

$$\epsilon \sigma T_s^4 = 2\epsilon \sigma T_a^4$$

From the atmospheric balance: $T_a^4 = T_s^4/2$, giving $T_a = T_s/2^{1/4}$. Substituting back into the surface balance:

$$T_s = T_e \left(\frac{2}{2 - \epsilon}\right)^{1/4}$$

For a perfectly absorbing atmosphere ($\epsilon = 1$): $T_s = 2^{1/4} T_e = 1.19 T_e \approx 303$ K. The real Earth has $T_s \approx 288$ K, corresponding to an effective emissivity $\epsilon \approx 0.77$.

Multi-Layer Model

Extending to $N$ perfectly absorbing layers, the surface temperature is:

$$T_s = (N + 1)^{1/4} T_e$$

Venus, with its thick CO$_2$ atmosphere, effectively has many absorbing layers, yielding a surface temperature of ~735 K despite receiving less solar radiation per unit area than Earth at its cloud tops.

3.3 Climate Sensitivity and Forcing

Derivation 3: Radiative Forcing and Climate Sensitivity

A perturbation to the energy balance (a radiative forcing $\Delta F$ in W/m$^2$) changes the surface temperature. Linearizing the Stefan-Boltzmann law around the current state:

$$\Delta F = \frac{dF}{dT}\Delta T = 4\sigma T_e^3 \Delta T$$

The no-feedback (Planck) climate sensitivity parameter is:

$$\lambda_0 = \frac{\Delta T}{\Delta F} = \frac{1}{4\sigma T_e^3} \approx 0.27 \text{ K/(W/m}^2\text{)}$$

Doubling CO$_2$ produces a forcing of $\Delta F_{2\times} \approx 3.7$ W/m$^2$(from the logarithmic absorption law: $\Delta F = 5.35 \ln(C/C_0)$ W/m$^2$). The Planck-only response would give $\Delta T \approx 1.0$ K. Feedbacks amplify this.

Derivation 4: Feedback Amplification

With feedbacks, the total response is amplified by the feedback factor $f$:

$$\Delta T = \frac{\lambda_0 \Delta F}{1 - f} = \frac{\Delta F}{\lambda_0^{-1} - \sum_i \lambda_i^{-1}}$$

where $f = \sum_i f_i$ is the sum of individual feedback parameters. The major feedbacks are:

Water vapor feedback ($f_{\text{WV}} \approx +0.3$): Clausius-Clapeyron relation dictates ~7% more moisture per K warming
Ice-albedo feedback ($f_{\text{ice}} \approx +0.1$): Melting ice reduces albedo, absorbing more solar radiation
Lapse rate feedback ($f_{\text{LR}} \approx -0.1$): Tropical upper troposphere warms faster, increasing emission
Cloud feedback ($f_{\text{cloud}} \approx +0.1$ to $+0.3$): Most uncertain; depends on cloud type and altitude changes

The total feedback factor $f \approx 0.4$-0.7 gives an equilibrium climate sensitivity (ECS) of $\Delta T_{2\times} \approx 2.5$-4.0 K for doubled CO$_2$. The IPCC AR6 (2021) assessed ECS as likely between 2.5 and 4.0 K, with a best estimate of 3.0 K.

3.4 The Ice-Albedo Feedback

Derivation 5: Stability Analysis of Ice-Albedo Feedback

The ice-albedo feedback can be modeled by making albedo a function of temperature. The energy balance becomes:

$$C \frac{dT}{dt} = \frac{(1 - \alpha(T))S_0}{4} - \sigma T^4$$

where $C$ is the heat capacity and $\alpha(T)$ transitions from high albedo ($\alpha_i \approx 0.6$) at cold temperatures to low albedo ($\alpha_w \approx 0.3$) at warm temperatures:

$$\alpha(T) = \begin{cases} \alpha_i & T < T_i \\ \alpha_i + (\alpha_w - \alpha_i)\frac{T - T_i}{T_w - T_i} & T_i \leq T \leq T_w \\ \alpha_w & T > T_w \end{cases}$$

The equilibrium states are found by setting $dT/dt = 0$. This yields three equilibria: a warm state, a cold (snowball) state, and an unstable intermediate state. The system exhibits hysteresis: once the planet enters a snowball state, a much larger forcing is required to escape than to enter it. Budyko (1969) and Sellers (1969) independently developed this analysis.

The stability of each equilibrium depends on the sign of the effective feedback parameter. At the unstable equilibrium, a small perturbation will drive the climate either toward the warm or snowball state, demonstrating the potential for rapid climate transitions.

3.5 Transient Climate Response

The ocean's large heat capacity delays the climate response to forcing. The transient climate response (TCR) measures the warming at the time of CO$_2$ doubling under a 1%/yr increase scenario. A simple energy balance model:

$$C_{\text{eff}} \frac{dT}{dt} = \Delta F - \frac{T}{\lambda}$$

where $C_{\text{eff}}$ is the effective heat capacity of the climate system (~17 W$\cdot$yr/(m$^2\cdot$K) for the mixed layer) and $\lambda$ is the total climate sensitivity parameter. The TCR is generally 50-80% of the ECS. IPCC AR6 assesses TCR as likely 1.4-2.2 K, with best estimate 1.8 K.

The CO$_2$ concentration has risen from ~280 ppm (pre-industrial) to ~425 ppm (2024), a forcing of $\Delta F \approx 5.35 \times \ln(425/280) \approx 2.2$ W/m$^2$. Including other greenhouse gases (CH$_4$, N$_2$O, halocarbons) and aerosol effects, the net anthropogenic forcing is approximately 2.7 W/m$^2$, consistent with the observed warming of ~1.2 K since pre-industrial times.

Emission Scenarios and Warming Projections

The IPCC Shared Socioeconomic Pathways (SSPs) project future warming depending on emissions trajectories. Under SSP1-2.6 (strong mitigation), warming is limited to ~1.8 K (range 1.3-2.4 K) by 2100. Under SSP5-8.5 (fossil-fuel intensive), warming reaches ~4.4 K (range 3.3-5.7 K). The remaining carbon budget for 1.5$°$C warming (from 2020) is approximately 400-500 Gt CO$_2$, which at current emission rates (~40 Gt CO$_2$/yr) would be exhausted by ~2030. Even under ambitious mitigation, net-negative emissions (carbon removal) are likely needed in the second half of the century to stabilize temperatures.

The committed warming from emissions already in the atmosphere depends on Earth system inertia. The ocean has absorbed about 90% of the excess heat, creating a thermal imbalance of ~0.8 W/m$^2$ at the top of the atmosphere. Even if emissions stopped today, temperatures would remain near current levels for decades to centuries, and sea level would continue rising for millennia due to thermal expansion and ice sheet dynamics.

3.6 The Global Carbon Cycle

Understanding climate change requires understanding the carbon cycle. The major reservoirs and their approximate carbon content are: atmosphere (~870 Gt C as of 2023), ocean (~38,000 Gt C), terrestrial biosphere (~2000 Gt C in vegetation, ~1500 Gt C in soils), and fossil fuels (~10,000 Gt C remaining). The geological reservoir (carbonates + organic carbon in rocks) contains ~100 million Gt C.

The Revelle Buffer Factor

The ocean's capacity to absorb CO$_2$ is governed by carbonate chemistry. The Revelle (buffer) factor describes the fractional change in dissolved CO$_2$relative to the fractional change in total dissolved inorganic carbon (DIC):

$$\mathcal{R} = \frac{\Delta[\text{CO}_2]/[\text{CO}_2]}{\Delta\text{DIC}/\text{DIC}} \approx 10\text{--}15$$

The large Revelle factor means that only ~1/10 of added CO$_2$ dissolves into the ocean at equilibrium. As atmospheric CO$_2$ rises, the Revelle factor increases, making the ocean an increasingly less efficient carbon sink. Roger Revelle and Hans Suess (1957) famously noted that humanity was conducting a "large-scale geophysical experiment" by burning fossil fuels.

Silicate Weathering Thermostat

On geological timescales ($> 10^5$ years), atmospheric CO$_2$ is regulated by the silicate weathering feedback. Chemical weathering of silicate rocks consumes CO$_2$:

$$\text{CaSiO}_3 + \text{CO}_2 \rightarrow \text{CaCO}_3 + \text{SiO}_2$$

Weathering rates increase with temperature (through the Arrhenius dependence of reaction kinetics and increased rainfall), creating a negative feedback: higher CO$_2$ $\rightarrow$ warmer climate $\rightarrow$ faster weathering $\rightarrow$CO$_2$ drawdown. This thermostat has maintained habitable conditions on Earth for ~4 billion years despite a ~30% increase in solar luminosity (the faint young Sun paradox). Walker, Hays, and Kasting (1981) formalized this feedback mechanism.

Computational Lab: Climate Science

Energy Balance, Greenhouse Effect, and Climate Sensitivity

Python

script.py379 lines

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

# ============================================================
# CLIMATE SCIENCE SIMULATION
# ============================================================
print("CLIMATE SCIENCE SIMULATION")
print("=" * 70)

sigma = 5.67e-8  # Stefan-Boltzmann constant
S0 = 1361.0      # solar constant (W/m^2)

# ============================================================
# 1. Effective Temperature of Planets
# ============================================================
print()
print("1. Effective and Surface Temperatures of Planets")
print("-" * 70)

planets = [
    ("Mercury", 9116, 0.12, 0.0, 440),
    ("Venus", 2613, 0.77, 0.99, 735),
    ("Earth", 1361, 0.30, 0.77, 288),
    ("Mars", 589, 0.25, 0.01, 218),
    ("Jupiter", 50.3, 0.50, 0.0, 124),
]

print(f"  {'Planet':12} {'S (W/m2)':12} {'Albedo':10} {'T_eff (K)':12} {'T_surf (K)':12} {'GH effect (K)'}")
print("  " + "-" * 65)

for name, S, alpha, eps, T_obs in planets:
    T_eff = ((1 - alpha) * S / (4 * sigma))**0.25
    if eps > 0:
        T_s_model = T_eff * (2 / (2 - eps))**0.25
    else:
        T_s_model = T_eff
    GH = T_obs - T_eff
    print(f"  {name:10} {S:10.1f} {alpha:8.2f} {T_eff:10.1f} {T_obs:10.0f} {GH:10.1f}")

# ============================================================
# 2. Multi-Layer Greenhouse Model
# ============================================================
print()
print()
print("2. Multi-Layer Greenhouse Model")
print("-" * 55)

T_eff_earth = ((1 - 0.30) * 1361 / (4 * sigma))**0.25
print(f"  Earth T_eff = {T_eff_earth:.1f} K")
print()
print(f"  {'N layers':12} {'T_surface (K)':16} {'T_surface (C)':16} {'GH warming (K)'}")
print("  " + "-" * 55)

for N in range(0, 11):
    T_s = (N + 1)**0.25 * T_eff_earth
    print(f"  {N:10} {T_s:14.1f} {T_s - 273.15:14.1f} {T_s - T_eff_earth:10.1f}")

# ============================================================
# 3. Radiative Forcing from CO2
# ============================================================
print()
print()
print("3. CO2 Radiative Forcing and Temperature Response")
print("-" * 65)

lambda_0 = 1.0 / (4 * sigma * 255**3)  # K/(W/m^2)
print(f"  Planck sensitivity parameter: {lambda_0:.3f} K/(W/m^2)")
print()

C0 = 280.0  # pre-industrial CO2 (ppm)
print(f"  {'CO2 (ppm)':14} {'Ratio':10} {'dF (W/m2)':14} {'dT_Planck (K)':16} {'dT_ECS (K)'}")
print("  " + "-" * 65)

ECS_lambda = 0.8  # K/(W/m^2) for ECS ~ 3K

for C in [280, 350, 400, 450, 500, 560, 700, 840, 1000, 1120]:
    ratio = C / C0
    dF = 5.35 * np.log(ratio)
    dT_planck = lambda_0 * dF
    dT_ecs = ECS_lambda * dF
    print(f"  {C:12} {ratio:8.2f} {dF:12.2f} {dT_planck:14.2f} {dT_ecs:10.2f}")

# ============================================================
# 4. Feedback Analysis
# ============================================================
print()
print()
print("4. Climate Feedback Analysis")
print("-" * 55)

feedbacks = [
    ("Planck (reference)", -3.2),
    ("Water vapor", 1.8),
    ("Lapse rate", -0.5),
    ("WV + LR combined", 1.3),
    ("Ice-albedo", 0.35),
    ("Cloud (SW)", 0.45),
    ("Cloud (LW)", 0.25),
    ("Total cloud", 0.70),
]

print(f"  {'Feedback':24} {'lambda_i (W/m2/K)':20} {'Contribution to ECS'}")
print("  " + "-" * 55)

total = 0
for name, lam in feedbacks:
    total += lam
    if "Total" not in name and "Planck" not in name and "combined" not in name:
        continue
    ecs_contrib = ""
    print(f"  {name:22} {lam:18.2f} {ecs_contrib}")

# Compute ECS from total feedback
total_net = -3.2 + 1.3 + 0.35 + 0.70
ECS = 3.7 / abs(total_net)
print(f"  Net feedback parameter: {total_net:.2f} W/(m^2.K)")
print(f"  ECS for 2xCO2: {ECS:.1f} K")

# ============================================================
# 5. Ice-Albedo Stability Analysis
# ============================================================
print()
print()
print("5. Ice-Albedo Feedback: Equilibrium States")
print("-" * 60)

alpha_ice = 0.6
alpha_warm = 0.3
T_i = 260.0  # ice threshold
T_w = 290.0  # warm threshold

print(f"  {'S/S0':10} {'T_warm (K)':14} {'T_cold (K)':14} {'T_unstable (K)':16} {'States'}")
print("  " + "-" * 60)

for S_ratio in [0.90, 0.92, 0.94, 0.96, 0.98, 1.00, 1.02, 1.04, 1.06, 1.08]:
    S = S_ratio * S0
    # Find equilibria numerically
    T_range = np.linspace(200, 330, 1000)
    equilibria = []
    for j in range(len(T_range) - 1):
        T = T_range[j]
        if T < T_i:
            alpha_T = alpha_ice
        elif T > T_w:
            alpha_T = alpha_warm
        else:
            alpha_T = alpha_ice + (alpha_warm - alpha_ice) * (T - T_i) / (T_w - T_i)
        absorbed = (1 - alpha_T) * S / 4
        emitted = sigma * T**4
        diff1 = absorbed - emitted

T2 = T_range[j+1]
        if T2 < T_i:
            alpha_T2 = alpha_ice
        elif T2 > T_w:
            alpha_T2 = alpha_warm
        else:
            alpha_T2 = alpha_ice + (alpha_warm - alpha_ice) * (T2 - T_i) / (T_w - T_i)
        diff2 = (1 - alpha_T2) * S / 4 - sigma * T2**4

if diff1 * diff2 < 0:
            equilibria.append((T + T2) / 2)

n = len(equilibria)
    if n >= 3:
        print(f"  {S_ratio:8.2f} {equilibria[2]:12.1f} {equilibria[0]:12.1f} {equilibria[1]:14.1f} 3 (bistable)")
    elif n == 2:
        print(f"  {S_ratio:8.2f} {equilibria[1]:12.1f} {equilibria[0]:12.1f} {'---':14} 2 (transition)")
    elif n == 1:
        state = "warm" if equilibria[0] > 280 else "snowball"
        print(f"  {S_ratio:8.2f} {equilibria[0]:12.1f} {'---':12} {'---':14} 1 ({state})")

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

# Plot 1: Energy Balance Diagram
ax = axes[0, 0]
T_range = np.linspace(200, 330, 500)

# Fixed albedo
absorbed_fixed = np.ones_like(T_range) * (1 - 0.3) * S0 / 4
emitted = sigma * T_range**4
ax.plot(T_range, absorbed_fixed, 'r-', linewidth=2, label='Absorbed (fixed albedo)')
ax.plot(T_range, emitted, 'b-', linewidth=2, label='Emitted (Stefan-Boltzmann)')

# With ice-albedo feedback
alpha_T_arr = np.where(T_range < T_i, alpha_ice,
              np.where(T_range > T_w, alpha_warm,
              alpha_ice + (alpha_warm - alpha_ice) * (T_range - T_i) / (T_w - T_i)))
absorbed_fb = (1 - alpha_T_arr) * S0 / 4
ax.plot(T_range, absorbed_fb, 'g-', linewidth=2, label='Absorbed (ice-albedo FB)')

# Mark equilibria
for T in T_range:
    if T < T_i:
        a = alpha_ice
    elif T > T_w:
        a = alpha_warm
    else:
        a = alpha_ice + (alpha_warm - alpha_ice) * (T - T_i) / (T_w - T_i)

ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Energy flux (W/m^2)')
ax.set_title('Energy Balance: Absorbed vs Emitted')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# Plot 2: CO2 forcing and temperature response
ax = axes[0, 1]
CO2_range = np.linspace(200, 1200, 200)
dF_range = 5.35 * np.log(CO2_range / 280)

for ecs, label in [(1.5, 'ECS=1.5K'), (3.0, 'ECS=3.0K'), (4.5, 'ECS=4.5K')]:
    lam_ecs = ecs / 3.7
    dT = lam_ecs * dF_range
    ax.plot(CO2_range, dT, linewidth=2, label=label)

ax.axvline(x=280, color='gray', linestyle=':', alpha=0.5, label='Pre-industrial')
ax.axvline(x=420, color='red', linestyle=':', alpha=0.5, label='Current (~420 ppm)')
ax.axvline(x=560, color='orange', linestyle=':', alpha=0.5, label='2x CO2')
ax.set_xlabel('CO2 concentration (ppm)')
ax.set_ylabel('Temperature change (K)')
ax.set_title('Global Warming vs CO2 Level')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# Plot 3: Transient response
ax = axes[1, 0]
years = np.arange(0, 300)
C_eff = 17.0  # W*yr/(m^2*K) effective heat capacity
lam_total = 1.0 / 0.85  # W/(m^2*K)

# 1% per year CO2 increase
CO2_t = 280 * 1.01**years
dF_t = 5.35 * np.log(CO2_t / 280)

# Numerical integration of dT/dt = (dF - dT/lam) / C_eff
dt = 1.0  # year
T_response = np.zeros(len(years))
for i in range(1, len(years)):
    dTdt = (dF_t[i-1] - T_response[i-1] * lam_total) / C_eff
    T_response[i] = T_response[i-1] + dTdt * dt

ax.plot(years, T_response, 'r-', linewidth=2, label='Temperature response')
ax.axhline(y=3.0, color='blue', linestyle='--', alpha=0.5, label='ECS = 3.0 K')
ax.axvline(x=70, color='green', linestyle=':', alpha=0.5, label='2x CO2 (year 70)')
ax.fill_between(years, T_response * 0.8, T_response * 1.2, alpha=0.15, color='red')

ax.set_xlabel('Years')
ax.set_ylabel('Temperature change (K)')
ax.set_title('Transient Climate Response (1%/yr CO2)')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 300)

# Plot 4: Ice-albedo hysteresis
ax = axes[1, 1]
S_range = np.linspace(0.85, 1.15, 300)
warm_branch = []
cold_branch = []
unstable_branch = []

for Sr in S_range:
    S = Sr * S0
    T_scan = np.linspace(200, 330, 2000)
    equil = []
    for j in range(len(T_scan) - 1):
        T = T_scan[j]
        if T < T_i:
            a = alpha_ice
        elif T > T_w:
            a = alpha_warm
        else:
            a = alpha_ice + (alpha_warm - alpha_ice) * (T - T_i) / (T_w - T_i)
        d1 = (1-a)*S/4 - sigma*T**4

T2 = T_scan[j+1]
        if T2 < T_i:
            a2 = alpha_ice
        elif T2 > T_w:
            a2 = alpha_warm
        else:
            a2 = alpha_ice + (alpha_warm - alpha_ice) * (T2 - T_i) / (T_w - T_i)
        d2 = (1-a2)*S/4 - sigma*T2**4

if d1 * d2 < 0:
            equil.append((T + T2) / 2)

if len(equil) >= 3:
        cold_branch.append((Sr, equil[0]))
        unstable_branch.append((Sr, equil[1]))
        warm_branch.append((Sr, equil[2]))
    elif len(equil) >= 1:
        if equil[-1] > 270:
            warm_branch.append((Sr, equil[-1]))
        else:
            cold_branch.append((Sr, equil[0]))

if warm_branch:
    wb = np.array(warm_branch)
    ax.plot(wb[:, 0], wb[:, 1], 'r-', linewidth=2, label='Warm state (stable)')
if cold_branch:
    cb = np.array(cold_branch)
    ax.plot(cb[:, 0], cb[:, 1], 'b-', linewidth=2, label='Snowball state (stable)')
if unstable_branch:
    ub = np.array(unstable_branch)
    ax.plot(ub[:, 0], ub[:, 1], 'k--', linewidth=1.5, label='Unstable equilibrium')

ax.axvline(x=1.0, color='green', linestyle=':', alpha=0.5, label='Current S/S0')
ax.set_xlabel('Solar constant ratio (S/S0)')
ax.set_ylabel('Equilibrium Temperature (K)')
ax.set_title('Ice-Albedo Hysteresis (Budyko-Sellers)')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
plt.close()
print()
print()
print("6. Carbon Cycle Reservoirs and Fluxes")
print("-" * 65)

reservoirs = [
    ("Atmosphere", 870, "CO2, CH4"),
    ("Ocean surface", 900, "DIC"),
    ("Ocean deep", 37100, "DIC"),
    ("Terrestrial biosphere", 2000, "Vegetation"),
    ("Soils", 1500, "Organic C"),
    ("Fossil fuels (remaining)", 10000, "Coal, oil, gas"),
    ("Carbonate rocks", 6e7, "CaCO3"),
    ("Organic C rocks", 1.5e7, "Kerogen"),
]

print(f"  {'Reservoir':30} {'C (Gt)':16} {'Type'}")
print("  " + "-" * 55)
for name, C, ctype in reservoirs:
    print(f"  {name:28} {C:14.0f} {ctype}")

print()
print("  Annual fluxes:")
fluxes = [
    ("Fossil fuel emissions", 10.0),
    ("Land use change", 1.1),
    ("Ocean uptake", -2.5),
    ("Land uptake", -3.1),
    ("Atmospheric growth", 5.5),
    ("Photosynthesis (gross)", -120),
    ("Respiration + fire", 120),
    ("Ocean-atmosphere exchange", 80),
    ("Volcanic outgassing", 0.15),
    ("Weathering sink", -0.3),
]

print(f"  {'Flux':30} {'Gt C/yr':12} {'Direction'}")
print("  " + "-" * 50)
for name, flux in fluxes:
    direction = "source" if flux > 0 else "sink"
    if abs(flux) < 0.01:
        direction = "~balanced"
    print(f"  {name:28} {flux:10.2f} {direction}")

print()
print("  Revelle buffer factor calculation:")
print("  " + "-" * 40)
for pCO2 in [280, 350, 400, 450, 500, 560, 700, 1000]:
    # Approximate: R increases with CO2
    R = 9.0 + 0.01 * (pCO2 - 280)
    ocean_fraction = 1.0 / R
    print(f"    pCO2 = {pCO2} ppm: R = {R:.1f}, ocean absorbs {ocean_fraction*100:.1f}% of added CO2")

print()
print()
print("[Chart saved: energy balance, CO2 forcing, transient response, ice-albedo hysteresis]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Carbon Budget and Emission Scenario Analysis

Python

script.py61 lines

import numpy as np

print("CARBON BUDGET AND WARMING SCENARIOS")
print("=" * 70)

print()
print("1. Remaining Carbon Budget for Temperature Targets")
print("-" * 60)

# TCRE: Transient Climate Response to Cumulative Emissions
# ~1.65 K per 1000 Gt CO2 (median IPCC AR6)
TCRE = 1.65 / 1000  # K per Gt CO2

current_warming = 1.2  # K above pre-industrial
emission_rate = 40  # Gt CO2/yr (2023)

print(f"  TCRE = {TCRE*1000:.2f} K per 1000 Gt CO2")
print(f"  Current warming = {current_warming} K, Emission rate = {emission_rate} Gt/yr")
print()
print(f"  {'Target (K)':14} {'Remaining dT':14} {'Budget (Gt CO2)':18} {'Years at current rate'}")
print("  " + "-" * 60)

for target in [1.5, 1.7, 2.0, 2.5, 3.0, 4.0]:
    remaining_dT = target - current_warming
    budget = remaining_dT / TCRE
    years = budget / emission_rate
    print(f"  {target:12.1f} {remaining_dT:12.1f} {budget:16.0f} {years:10.0f}")

print()
print()
print("2. CO2 Pathway for Different Scenarios")
print("-" * 60)

print(f"  {'Year':10}", end="")
scenarios = [
    ("SSP1-2.6", -0.03),
    ("SSP2-4.5", 0.005),
    ("SSP3-7.0", 0.015),
    ("SSP5-8.5", 0.025),
]
for name, _ in scenarios:
    print(f"{name:14}", end="")
print()
print("  " + "-" * 60)

for year in range(2025, 2105, 10):
    print(f"  {year:8}", end="")
    for name, rate in scenarios:
        dt = year - 2025
        if name == "SSP1-2.6":
            co2 = 420 * np.exp(rate * dt) if dt < 30 else 420 * np.exp(rate*30) * np.exp(-0.01*(dt-30))
        else:
            co2 = 420 * np.exp(rate * dt)
        print(f"  {co2:12.0f}", end="")
    print(" ppm")

print()
print("  Note: SSP1-2.6 assumes net-zero by ~2070, decline after")
print()
print("Complete: carbon budget analysis finished.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Geodynamics Next: Atmospheric Physics →