Module 6

Carbon Cycle & Climate Feedbacks

Global carbon budget, ocean and land sinks, feedback mechanisms, and tipping points

The carbon cycle is the master regulator of Earth’s climate on timescales from years to millennia. Each year, humanity emits roughly 10 Gt C from fossil fuels, yet only about 44% remains in the atmosphere — the rest is absorbed by oceans and terrestrial ecosystems. Understanding where this carbon goes, how those sinks will respond to warming, and which feedback mechanisms could amplify or dampen climate change is essential for projecting future temperatures and designing effective mitigation strategies.

6.1 The Global Carbon Budget

The fundamental equation of the global carbon budget describes the rate of change of atmospheric carbon as a balance between sources and sinks:

$\frac{dC_{\text{atm}}}{dt} = E_{\text{fossil}} + E_{\text{LUC}} - S_{\text{ocean}} - S_{\text{land}}$

where $E_{\text{fossil}}$ is fossil fuel and cement emissions (~9.5 Gt C/yr in 2023),$E_{\text{LUC}}$ is land-use change emissions (~1.1 Gt C/yr), $S_{\text{ocean}}$ is the ocean carbon sink (~2.8 Gt C/yr), and $S_{\text{land}}$ is the terrestrial carbon sink (~3.1 Gt C/yr). The budget imbalance is ~0.3 Gt C/yr.

Airborne Fraction

The airborne fraction (AF) is defined as the fraction of total emissions that remains in the atmosphere:

$\text{AF} = \frac{dC_{\text{atm}}/dt}{E_{\text{fossil}} + E_{\text{LUC}}} = 1 - \frac{S_{\text{ocean}} + S_{\text{land}}}{E_{\text{total}}} \approx 0.44$

Remarkably, the airborne fraction has remained approximately constant at 44% over the past six decades despite emissions roughly tripling. This constancy arises because both ocean and land sinks have scaled roughly proportionally with emissions. The ocean sink responds to the growing disequilibrium between atmospheric pCO$_2$ and dissolved CO$_2$, while land photosynthesis benefits from CO$_2$ fertilization. However, this proportionality is not guaranteed to continue — both sinks show signs of saturation.

We can derive the constancy by writing sink strengths as proportional to emissions with a time lag. If $S_{\text{ocean}} \approx \alpha_o \cdot E_{\text{total}}$ and$S_{\text{land}} \approx \alpha_l \cdot E_{\text{total}}$, then AF $= 1 - \alpha_o - \alpha_l$. The stability of AF implies$\alpha_o + \alpha_l \approx 0.56$ has been roughly constant. This is a coincidence of competing effects: the ocean sink efficiency is declining (Revelle factor increases) while the land sink has been increasing.

Carbon Reservoirs

Atmosphere: ~870 Gt C (as of 2023, ~420 ppm CO$_2$)
Ocean surface: ~900 Gt C (dissolved inorganic carbon)
Deep ocean: ~37,100 Gt C (vast but slow exchange)
Vegetation: ~450 Gt C (living biomass)
Soils: ~1,700 Gt C (top 1 m), ~3,000 Gt C including permafrost
Fossil fuels: ~4,000 Gt C remaining (coal, oil, gas)

6.2 The Ocean Carbon Sink

The ocean absorbs roughly 25% of anthropogenic CO$_2$ emissions despite containing 50 times more carbon than the atmosphere. Why only 25%? The answer lies in ocean chemistry, specifically the Revelle factor.

Ocean Carbon Chemistry

When CO$_2$ dissolves in seawater, it undergoes a series of equilibrium reactions:

$\text{CO}_2(\text{g}) \rightleftharpoons \text{CO}_2(\text{aq}) + \text{H}_2\text{O} \rightleftharpoons \text{H}_2\text{CO}_3 \rightleftharpoons \text{H}^+ + \text{HCO}_3^- \rightleftharpoons 2\text{H}^+ + \text{CO}_3^{2-}$

The total dissolved inorganic carbon (DIC) is:

$\text{DIC} = [\text{CO}_2(\text{aq})] + [\text{HCO}_3^-] + [\text{CO}_3^{2-}]$

At the pH of seawater (~8.1), about 90% of DIC is in the form of bicarbonate (HCO$_3^-$), ~9% as carbonate (CO$_3^{2-}$), and only ~1% as dissolved CO$_2$. This buffer chemistry means the ocean can absorb large amounts of CO$_2$ but with diminishing efficiency.

The Revelle Factor

The Revelle factor (or buffer factor) quantifies how much the ocean resists additional CO$_2$ uptake:

$R = \frac{\Delta p\text{CO}_2 / p\text{CO}_2}{\Delta \text{DIC} / \text{DIC}} \approx 10$

A Revelle factor of ~10 means that a 10% increase in atmospheric CO$_2$ produces only a ~1% increase in ocean DIC. Despite the ocean’s enormous DIC reservoir, the buffer chemistry limits the fraction of emissions that the ocean can absorb. As atmospheric CO$_2$ rises, the Revelle factor increases (the ocean becomes less effective as a sink), a positive feedback.

We can derive why the ocean absorbs only ~25%: if the atmosphere has ~870 Gt C and the surface ocean ~900 Gt C in DIC, a naive partitioning would give ~50% ocean uptake. But the Revelle factor reduces this by a factor of $\sim 1/R$:

$f_{\text{ocean}} \approx \frac{C_{\text{ocean,surface}}/R}{C_{\text{atm}} + C_{\text{ocean,surface}}/R} = \frac{900/10}{870 + 900/10} = \frac{90}{960} \approx 0.094$

When we account for the thermocline and deep ocean mixing on decadal timescales, the effective ocean volume increases, bringing the fraction closer to the observed ~25%. The biological pump (sinking organic matter and carbonate shells) also transfers carbon from surface to deep ocean, enhancing the sink.

Ocean Acidification

As the ocean absorbs CO$_2$, pH decreases. Since pre-industrial times, ocean pH has dropped from ~8.2 to ~8.1, a 26% increase in hydrogen ion concentration (pH is logarithmic). This ocean acidification threatens calcifying organisms (corals, shellfish, pteropods) and has cascading effects through marine food webs. See the Climate & Biodiversity course for ecological impacts.

6.3 The Terrestrial Carbon Sink

Terrestrial ecosystems absorb roughly 30% of anthropogenic CO$_2$ emissions. The net flux is determined by the balance between photosynthetic uptake and respiratory release:

$\text{NEP} = \text{GPP} - R_{\text{auto}} - R_{\text{hetero}} \approx 3 \text{ Gt C/yr}$

where GPP (gross primary production) is total photosynthesis (~120 Gt C/yr), $R_{\text{auto}}$ is autotrophic respiration by plants (~60 Gt C/yr), and $R_{\text{hetero}}$ is heterotrophic respiration by decomposers (~57 Gt C/yr). The small residual NEP (net ecosystem production) ~3 Gt C/yr is the current land carbon sink.

CO$_2$ Fertilization

Rising CO$_2$ enhances photosynthesis, particularly in C3 plants (which include most trees and crops). The relationship between photosynthesis rate A and CO$_2$ concentration follows a saturating curve:

$A = A_{\max} \cdot \frac{c_i - \Gamma^*}{c_i + K_c(1 + O/K_o)}$

This is the Farquhar–von Caemmerer–Berry model where $c_i$ is the intercellular CO$_2$ concentration, $\Gamma^*$ is the CO$_2$ compensation point, and $K_c$, $K_o$ are Michaelis constants for carboxylation and oxygenation.

FACE Experiments

Free-Air CO$_2$ Enrichment (FACE) experiments expose ecosystems to elevated CO$_2$ (~550 ppm) and measure responses over years. Key findings:

Initial GPP increase of 20–30% in forests, declining over 3–5 years
Nitrogen limitation constrains sustained fertilization in most ecosystems
Phosphorus limitation dominant in tropical forests and old soils
Increased water use efficiency (stomata partially close at high CO$_2$)
Soil carbon response highly variable, sometimes negative

The implication is that the land sink may not continue to grow proportionally with emissions. Nutrient limitation, drought stress, increased fire frequency, and deforestation all threaten to weaken or reverse the terrestrial sink this century. For more on soil carbon dynamics, see Climate & Biodiversity Module 8 on adaptation strategies.

Global Carbon Cycle

This diagram shows the major carbon reservoirs (boxes) and annual fluxes (arrows) in the global carbon cycle. Anthropogenic perturbations are shown in red.

Figure 6.1: Global carbon cycle showing major reservoirs (Gt C) and annual fluxes (Gt C/yr). Red arrows indicate anthropogenic perturbations. Values from Global Carbon Project (2023).

Carbon Budget Projection

This simulation projects atmospheric CO$_2$ concentration from 1959 to 2100 using the carbon budget equation. We model emission scenarios and compute the resulting airborne fraction and Revelle factor evolution.

Global Carbon Budget: Historical Fit & Future Projections

Python

script.py131 lines

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

# ======== Historical emissions data (simplified) ========
years_hist = np.arange(1959, 2024)
# Fossil emissions: exponential growth ~2%/yr until 2000, slower after
e_fossil_hist = 2.5 * np.exp(0.025 * (years_hist - 1959))
e_fossil_hist[years_hist > 2000] = e_fossil_hist[years_hist == 2000] * (1 + 0.015 * (years_hist[years_hist > 2000] - 2000))
e_fossil_hist = np.clip(e_fossil_hist, 0, 10.5)

# Land-use change: declining from 1.5 to 1.0
e_luc_hist = 1.5 - 0.008 * (years_hist - 1959)
e_luc_hist = np.clip(e_luc_hist, 0.9, 1.5)

# ======== Future scenarios (2024-2100) ========
years_future = np.arange(2024, 2101)
scenarios = {
    'SSP1-2.6': {'peak': 2025, 'decline_rate': 0.06, 'color': '#4ade80', 'luc_end': 0.3},
    'SSP2-4.5': {'peak': 2030, 'decline_rate': 0.025, 'color': '#fbbf24', 'luc_end': 0.5},
    'SSP5-8.5': {'peak': 2080, 'decline_rate': -0.015, 'color': '#f87171', 'luc_end': 0.8},
}

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# ======== Panel 1: Emissions ========
ax1 = axes[0, 0]
ax1.plot(years_hist, e_fossil_hist + e_luc_hist, 'w-', linewidth=2, label='Historical')

for name, params in scenarios.items():
    e_fossil_fut = np.copy(e_fossil_hist[-1] * np.ones(len(years_future)))
    for i, yr in enumerate(years_future):
        if yr <= params['peak']:
            e_fossil_fut[i] = e_fossil_hist[-1] * (1 + 0.01 * (yr - 2024))
        else:
            e_fossil_fut[i] = e_fossil_hist[-1] * (1 + 0.01 * (params['peak'] - 2024)) * \
                              np.exp(-params['decline_rate'] * (yr - params['peak']))
    e_fossil_fut = np.clip(e_fossil_fut, 0, 25)
    e_luc_fut = np.linspace(e_luc_hist[-1], params['luc_end'], len(years_future))
    total = e_fossil_fut + e_luc_fut
    ax1.plot(years_future, total, color=params['color'], linewidth=2, label=name)

ax1.set_xlabel('Year', color='white')
ax1.set_ylabel('Emissions (Gt C/yr)', color='white')
ax1.set_title('Total Carbon Emissions', color='white', fontsize=14)
ax1.legend(fontsize=9)
ax1.set_facecolor('#0a0a1a')
ax1.tick_params(colors='white')
ax1.grid(True, alpha=0.2)

# ======== Panel 2: Atmospheric CO2 ========
ax2 = axes[0, 1]
# Simple model: CO2(t+1) = CO2(t) + AF * E_total * (2.12 ppm per Gt C)
ppm_per_gtc = 2.12
co2_hist = [315.0]  # 1959 Mauna Loa
for i in range(1, len(years_hist)):
    af = 0.44
    co2_hist.append(co2_hist[-1] + af * (e_fossil_hist[i] + e_luc_hist[i]) * ppm_per_gtc)
co2_hist = np.array(co2_hist)

ax2.plot(years_hist, co2_hist, 'w-', linewidth=2, label='Historical model')

for name, params in scenarios.items():
    co2_fut = [co2_hist[-1]]
    e_fossil_fut = np.copy(e_fossil_hist[-1] * np.ones(len(years_future)))
    for i, yr in enumerate(years_future):
        if yr <= params['peak']:
            e_fossil_fut[i] = e_fossil_hist[-1] * (1 + 0.01 * (yr - 2024))
        else:
            e_fossil_fut[i] = e_fossil_hist[-1] * (1 + 0.01 * (params['peak'] - 2024)) * \
                              np.exp(-params['decline_rate'] * (yr - params['peak']))
    e_fossil_fut = np.clip(e_fossil_fut, 0, 25)
    e_luc_fut = np.linspace(e_luc_hist[-1], params['luc_end'], len(years_future))
    for i in range(1, len(years_future)):
        # Revelle-adjusted airborne fraction: increases with CO2
        revelle = 10 + 5 * (co2_fut[-1] - 280) / 280
        af_adj = 0.44 * (revelle / 10)
        af_adj = np.clip(af_adj, 0.35, 0.70)
        co2_fut.append(co2_fut[-1] + af_adj * (e_fossil_fut[i] + e_luc_fut[i]) * ppm_per_gtc)
    ax2.plot(years_future, co2_fut, color=params['color'], linewidth=2, label=name)

ax2.set_xlabel('Year', color='white')
ax2.set_ylabel('CO$_2$ (ppm)', color='white')
ax2.set_title('Atmospheric CO$_2$ Concentration', color='white', fontsize=14)
ax2.legend(fontsize=9)
ax2.set_facecolor('#0a0a1a')
ax2.tick_params(colors='white')
ax2.grid(True, alpha=0.2)

# ======== Panel 3: Airborne Fraction ========
ax3 = axes[1, 0]
# Historical AF
af_hist = np.gradient(co2_hist) / ppm_per_gtc / (e_fossil_hist + e_luc_hist)
af_hist = np.clip(af_hist, 0.2, 0.7)
ax3.plot(years_hist, af_hist, 'w-', linewidth=1.5, alpha=0.7)
ax3.axhline(y=0.44, color='#818cf8', linestyle='--', linewidth=1.5, label='Mean AF = 0.44')
ax3.fill_between(years_hist, 0.35, 0.55, alpha=0.15, color='#818cf8')
ax3.set_xlabel('Year', color='white')
ax3.set_ylabel('Airborne Fraction', color='white')
ax3.set_title('Airborne Fraction (Historical)', color='white', fontsize=14)
ax3.legend(fontsize=9)
ax3.set_facecolor('#0a0a1a')
ax3.tick_params(colors='white')
ax3.grid(True, alpha=0.2)
ax3.set_ylim(0.2, 0.7)

# ======== Panel 4: Revelle Factor ========
ax4 = axes[1, 1]
co2_range = np.linspace(280, 1000, 200)
revelle = 10 + 5 * (co2_range - 280) / 280
ax4.plot(co2_range, revelle, color='#60a5fa', linewidth=2.5)
ax4.axvline(x=280, color='#4ade80', linestyle=':', linewidth=1.5, label='Pre-industrial (280 ppm)')
ax4.axvline(x=420, color='#fbbf24', linestyle=':', linewidth=1.5, label='Present (420 ppm)')
ax4.fill_between(co2_range, revelle, alpha=0.1, color='#60a5fa')
ax4.set_xlabel('Atmospheric CO$_2$ (ppm)', color='white')
ax4.set_ylabel('Revelle Factor', color='white')
ax4.set_title('Revelle Factor vs CO$_2$', color='white', fontsize=14)
ax4.legend(fontsize=9)
ax4.set_facecolor('#0a0a1a')
ax4.tick_params(colors='white')
ax4.grid(True, alpha=0.2)

fig.patch.set_facecolor('#0a0a1a')
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')
plt.close()
print("Carbon budget projection complete.")
print(f"Historical CO2 range: {co2_hist[0]:.0f} - {co2_hist[-1]:.0f} ppm")
print(f"Revelle factor at 420 ppm: {10 + 5*(420-280)/280:.1f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6.4 Climate Feedback Mechanisms

Climate feedbacks amplify or dampen the initial warming from greenhouse gas forcing. The total temperature response to a radiative forcing $\Delta F$ is:

$\Delta T = \frac{\Delta F}{\lambda_{\text{Planck}}} \cdot f = \Delta T_0 \cdot f$

where $\lambda_{\text{Planck}} \approx 3.2$ W/m$^2$/K is the Planck response (the no-feedback sensitivity),$\Delta T_0 = \Delta F / \lambda_{\text{Planck}} \approx 1.2$ K for CO$_2$ doubling ($\Delta F = 3.7$ W/m$^2$), and $f$ is the feedback factor.

Feedback Factor Derivation

Each feedback process $i$ contributes a feedback parameter $\lambda_i$ (W/m$^2$/K) that modifies the radiative balance. The total feedback parameter is $\lambda = \sum_i \lambda_i$. The feedback factor is:

$f = \frac{1}{1 - \sum_i \lambda_i / \lambda_{\text{Planck}}} = \frac{1}{1 - g}$

where $g = \sum_i \lambda_i / \lambda_{\text{Planck}}$ is the total gain. If $g < 1$, the system is stable; if $g \to 1$, the system approaches a runaway. Current estimates give $g \approx 0.5\text{--}0.7$, so $f \approx 2\text{--}3$, amplifying the no-feedback warming by a factor of 2–3.

Major Feedback Mechanisms

Feedback

$\lambda_i$ (W/m$^2$/K)

Sign

Mechanism

Water vapor

+1.8

Positive

Clausius–Clapeyron: 7%/K more moisture

Ice–albedo

+0.3

Positive

Melting ice exposes dark ocean/land

Lapse rate

−0.6

Negative

Tropics warm more aloft than surface

Cloud

0 to +0.7

Uncertain

Low cloud reduction vs. high cloud increase

The water vapor feedback is the strongest single feedback. The Clausius–Clapeyron equation dictates that saturation vapor pressure increases exponentially with temperature (~7%/K). Since water vapor is a greenhouse gas, this creates a powerful positive feedback that roughly doubles the warming from CO$_2$ alone.

The lapse rate feedback is the main negative feedback: in the tropics, convection distributes warming more to the upper troposphere, which radiates more efficiently to space. This partially (but not fully) cancels the water vapor feedback. The combined water vapor + lapse rate feedback is $\approx +1.2$ W/m$^2$/K.

The cloud feedback remains the largest source of uncertainty in climate sensitivity. Low marine clouds may thin as temperatures rise (reducing their reflective cooling effect), while high clouds may rise (trapping more longwave radiation). Recent observations and high-resolution simulations suggest the net cloud feedback is likely positive, narrowing the long-standing uncertainty.

Feedback Analysis Simulation

This simulation computes equilibrium climate sensitivity (ECS) as a function of individual feedback strengths and visualizes the nonlinear amplification of the feedback factor.

Climate Feedback Analysis: ECS & Feedback Factor

Python

script.py133 lines

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

# ======== Feedback parameters ========
lambda_planck = 3.2  # W/m^2/K (Planck response)
delta_F_2xCO2 = 3.7  # W/m^2 (forcing from CO2 doubling)

# Individual feedbacks (W/m^2/K)
feedbacks = {
    'Water vapor': 1.8,
    'Ice-albedo': 0.3,
    'Lapse rate': -0.6,
    'Cloud (best est.)': 0.4,
}

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# ======== Panel 1: Feedback contributions ========
ax1 = axes[0, 0]
names = list(feedbacks.keys())
values = list(feedbacks.values())
colors = ['#60a5fa' if v > 0 else '#f87171' for v in values]
bars = ax1.barh(names, values, color=colors, edgecolor='white', linewidth=0.5)
ax1.axvline(x=0, color='white', linewidth=1)
ax1.set_xlabel('Feedback parameter (W/m$^2$/K)', color='white')
ax1.set_title('Individual Feedback Strengths', color='white', fontsize=14)
for bar, val in zip(bars, values):
    ax1.text(val + 0.05 if val > 0 else val - 0.15, bar.get_y() + bar.get_height()/2,
             f'{val:+.1f}', va='center', color='white', fontsize=11)
ax1.set_facecolor('#0a0a1a')
ax1.tick_params(colors='white')

# ======== Panel 2: Feedback factor vs gain ========
ax2 = axes[0, 1]
g_range = np.linspace(-0.5, 0.95, 300)
f_factor = 1 / (1 - g_range)

ax2.plot(g_range, f_factor, color='#818cf8', linewidth=2.5)
ax2.axhline(y=1, color='white', linestyle=':', alpha=0.5)
ax2.axvline(x=0, color='white', linestyle=':', alpha=0.5)

# Mark current best estimate
g_current = sum(feedbacks.values()) / lambda_planck
f_current = 1 / (1 - g_current)
ax2.plot(g_current, f_current, 'o', color='#fbbf24', markersize=10, zorder=5)
ax2.annotate(f'Current: g={g_current:.2f}, f={f_current:.1f}',
             xy=(g_current, f_current), xytext=(g_current - 0.3, f_current + 3),
             color='#fbbf24', fontsize=11,
             arrowprops=dict(arrowstyle='->', color='#fbbf24'))

ax2.set_xlabel('Total gain g = $\\Sigma\\lambda_i / \\lambda_{Planck}$', color='white')
ax2.set_ylabel('Feedback factor f = 1/(1-g)', color='white')
ax2.set_title('Nonlinear Feedback Amplification', color='white', fontsize=14)
ax2.set_ylim(0, 15)
ax2.set_facecolor('#0a0a1a')
ax2.tick_params(colors='white')
ax2.grid(True, alpha=0.2)

# ======== Panel 3: ECS distribution ========
ax3 = axes[1, 0]
# Monte Carlo sampling of ECS
np.random.seed(42)
n_samples = 50000
wv = np.random.normal(1.8, 0.3, n_samples)
ia = np.random.normal(0.3, 0.1, n_samples)
lr = np.random.normal(-0.6, 0.15, n_samples)
cl = np.random.normal(0.4, 0.35, n_samples)  # cloud most uncertain

total_lambda = wv + ia + lr + cl
gain = total_lambda / lambda_planck
gain = np.clip(gain, -1, 0.95)  # prevent runaway
f_samples = 1 / (1 - gain)
ecs = delta_F_2xCO2 / lambda_planck * f_samples

# Remove extreme outliers
ecs = ecs[(ecs > 0) & (ecs < 10)]

ax3.hist(ecs, bins=80, color='#818cf8', alpha=0.7, edgecolor='#0a0a1a', density=True)
ax3.axvline(x=np.median(ecs), color='#fbbf24', linewidth=2, label=f'Median: {np.median(ecs):.1f} K')
ax3.axvline(x=np.percentile(ecs, 5), color='#4ade80', linewidth=1.5, linestyle='--',
            label=f'5th pct: {np.percentile(ecs, 5):.1f} K')
ax3.axvline(x=np.percentile(ecs, 95), color='#f87171', linewidth=1.5, linestyle='--',
            label=f'95th pct: {np.percentile(ecs, 95):.1f} K')
ax3.set_xlabel('ECS (K)', color='white')
ax3.set_ylabel('Probability density', color='white')
ax3.set_title('Equilibrium Climate Sensitivity Distribution', color='white', fontsize=14)
ax3.legend(fontsize=9)
ax3.set_facecolor('#0a0a1a')
ax3.tick_params(colors='white')
ax3.grid(True, alpha=0.2)

# ======== Panel 4: Cumulative feedback effect ========
ax4 = axes[1, 1]
delta_T0 = delta_F_2xCO2 / lambda_planck  # no-feedback warming

feedback_names = ['No feedback', '+ Water vapor', '+ Ice-albedo', '+ Lapse rate', '+ Cloud']
cumulative = [0]
running = 0
for v in feedbacks.values():
    running += v
    cumulative.append(running)

delta_Ts = []
for c in cumulative:
    g = c / lambda_planck
    if g >= 1:
        g = 0.95
    delta_Ts.append(delta_T0 / (1 - g))

bar_colors = ['#818cf8', '#60a5fa', '#4ade80', '#f87171', '#fbbf24']
bars = ax4.bar(feedback_names, delta_Ts, color=bar_colors, edgecolor='white', linewidth=0.5)
for bar, val in zip(bars, delta_Ts):
    ax4.text(bar.get_x() + bar.get_width()/2, val + 0.1, f'{val:.1f} K',
             ha='center', color='white', fontsize=11)

ax4.set_ylabel('$\\Delta T$ (K) for 2xCO$_2$', color='white')
ax4.set_title('Cumulative Effect of Feedbacks', color='white', fontsize=14)
ax4.set_facecolor('#0a0a1a')
ax4.tick_params(colors='white', rotation=25)
ax4.grid(True, alpha=0.2, axis='y')

fig.patch.set_facecolor('#0a0a1a')
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')
plt.close()
print("Feedback analysis complete.")
print(f"No-feedback warming (2xCO2): {delta_T0:.2f} K")
print(f"With all feedbacks: {delta_Ts[-1]:.2f} K")
print(f"Feedback factor: {f_current:.2f}")
print(f"ECS median: {np.median(ecs):.1f} K, 5-95% range: {np.percentile(ecs, 5):.1f}-{np.percentile(ecs, 95):.1f} K")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6.5 Tipping Points & Early Warning Signals

Climate tipping points are thresholds beyond which self-reinforcing feedbacks drive a subsystem into a qualitatively different state. The concept is formalized through bifurcation theory: a system parameter crosses a critical value, and the system transitions abruptly to a new attractor.

Major Tipping Elements

Amazon rainforest dieback: Warming + deforestation reduces rainfall recycling. At ~3–4°C or 40% deforestation, savannification may become self-sustaining. Loss of ~50 Gt C from vegetation.
AMOC shutdown: Freshwater from Greenland melt disrupts the Atlantic Meridional Overturning Circulation. Could cool northern Europe by 5–10°C while warming the Southern Hemisphere. Recent AMOC weakening detected.
Permafrost carbon release: ~1,500 Gt C in permafrost soils. Thawing releases CO$_2$ and CH$_4$, creating a positive feedback. Could add 0.3°C by 2100 in high-emission scenarios.
West Antarctic Ice Sheet: Marine ice sheet instability (see Module 7) could cause irreversible collapse, raising sea level by ~3.3 m over centuries.
Monsoon disruption: Nonlinear shifts in Indian and West African monsoons could affect food and water security for billions.

Critical Slowing Down

Near a tipping point, the system loses resilience — perturbations take longer to decay. This is called critical slowing down. Consider a simple system near a fold bifurcation:

$\frac{dx}{dt} = -\frac{dV}{dx} = r + x^2 \quad \Rightarrow \quad \text{fixed points at } x^* = \pm\sqrt{-r}$

As the control parameter $r$ approaches zero from below, the stable and unstable fixed points approach each other and merge at $r = 0$ (the bifurcation point). The eigenvalue of the linearized system $\lambda = -2\sqrt{-r}$ approaches zero, meaning the system’s recovery rate vanishes.

Observable early warning signals (EWS) include:

Increased autocorrelation: the system “remembers” perturbations longer (lag-1 autocorrelation $\to 1$)
Increased variance: fluctuations grow as the restoring force weakens
Flickering: the system alternates between two states before committing
Increased skewness: distribution becomes asymmetric toward the alternative state

These EWS have been detected in paleoclimate records before abrupt transitions (Dansgaard–Oeschger events, the end of the African Humid Period) and are being monitored in modern observations of the AMOC and Amazon rainfall. For ecological impacts of these transitions, see the Climate & Biodiversity course.

Tipping Point Early Warning Simulation

This simulation demonstrates critical slowing down as a system approaches a fold bifurcation, showing how autocorrelation and variance increase before the tipping point.

Tipping Point: Critical Slowing Down & Early Warning Signals

Python

script.py129 lines

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

np.random.seed(42)

# ======== Simulate system approaching fold bifurcation ========
# dx/dt = -(x - x_eq)^2 + r + noise
# As r approaches 0 from below, the system tips

n_steps = 5000
dt = 0.01
sigma = 0.15  # noise amplitude

# Control parameter r drifts from -2 toward 0 (tipping)
r = np.linspace(-2.0, 0.05, n_steps)

x = np.zeros(n_steps)
x[0] = -np.sqrt(max(-r[0], 0)) if r[0] < 0 else 0  # start at stable eq

for i in range(1, n_steps):
    if r[i] < 0:
        x_eq = -np.sqrt(-r[i])
    else:
        x_eq = 0
    # Potential: V = (x - x_eq)^2 * restoring - r*x
    dxdt = -(x[i-1] - x_eq) * 2 * np.sqrt(max(-r[i], 0.001)) + sigma * np.random.randn() / np.sqrt(dt)
    x[i] = x[i-1] + dxdt * dt

# After tipping (r > 0), system diverges
for i in range(1, n_steps):
    if r[i] > 0 and i > 0:
        x[i] = x[i-1] + (0.5 + x[i-1]**0.3) * dt + sigma * np.random.randn() * np.sqrt(dt)

time = np.arange(n_steps) * dt

# ======== Compute early warning signals ========
window = 200
autocorr = np.zeros(n_steps)
variance = np.zeros(n_steps)

for i in range(window, n_steps):
    segment = x[i-window:i]
    detrended = segment - np.mean(segment)
    if np.std(detrended) > 1e-10:
        autocorr[i] = np.corrcoef(detrended[:-1], detrended[1:])[0, 1]
    variance[i] = np.var(detrended)

autocorr[:window] = autocorr[window]
variance[:window] = variance[window]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# ======== Panel 1: Bifurcation diagram ========
ax1 = axes[0, 0]
r_plot = np.linspace(-2, 0.5, 500)
stable_branch = np.where(r_plot < 0, -np.sqrt(-r_plot), np.nan)
unstable_branch = np.where(r_plot < 0, np.sqrt(-r_plot), np.nan)

ax1.plot(r_plot, stable_branch, color='#4ade80', linewidth=2.5, label='Stable equilibrium')
ax1.plot(r_plot, unstable_branch, color='#f87171', linewidth=2.5, linestyle='--', label='Unstable equilibrium')
ax1.axvline(x=0, color='#fbbf24', linewidth=2, linestyle=':', label='Tipping point (r=0)')
ax1.annotate('BIFURCATION', xy=(0, 0), xytext=(0.15, 0.8),
             color='#fbbf24', fontsize=14, fontweight='bold',
             arrowprops=dict(arrowstyle='->', color='#fbbf24'))
ax1.set_xlabel('Control parameter r', color='white')
ax1.set_ylabel('Equilibrium x*', color='white')
ax1.set_title('Fold Bifurcation Diagram', color='white', fontsize=14)
ax1.legend(fontsize=9)
ax1.set_facecolor('#0a0a1a')
ax1.tick_params(colors='white')
ax1.grid(True, alpha=0.2)

# ======== Panel 2: Time series ========
ax2 = axes[0, 1]
tip_idx = np.argmax(r >= 0)
ax2.plot(time[:tip_idx], x[:tip_idx], color='#60a5fa', linewidth=0.5, alpha=0.8)
ax2.plot(time[tip_idx:], x[tip_idx:], color='#f87171', linewidth=0.5, alpha=0.8)
ax2.axvline(x=time[tip_idx], color='#fbbf24', linewidth=2, linestyle=':', label='Tipping point')
ax2.set_xlabel('Time', color='white')
ax2.set_ylabel('State variable x', color='white')
ax2.set_title('System State Approaching Tipping Point', color='white', fontsize=14)
ax2.legend(fontsize=9)
ax2.set_facecolor('#0a0a1a')
ax2.tick_params(colors='white')
ax2.grid(True, alpha=0.2)

# ======== Panel 3: Autocorrelation (EWS) ========
ax3 = axes[1, 0]
ax3.plot(time, autocorr, color='#c084fc', linewidth=1.5)
ax3.axvline(x=time[tip_idx], color='#fbbf24', linewidth=2, linestyle=':', label='Tipping point')
ax3.set_xlabel('Time', color='white')
ax3.set_ylabel('Lag-1 Autocorrelation', color='white')
ax3.set_title('Early Warning: Autocorrelation', color='white', fontsize=14)
ax3.legend(fontsize=9)
ax3.set_facecolor('#0a0a1a')
ax3.tick_params(colors='white')
ax3.grid(True, alpha=0.2)
ax3.annotate('Rising autocorrelation\n= critical slowing down',
             xy=(time[tip_idx]-5, 0.85), fontsize=10, color='#c084fc',
             bbox=dict(boxstyle='round,pad=0.3', facecolor='#1e1b4b', edgecolor='#c084fc', alpha=0.8))

# ======== Panel 4: Variance (EWS) ========
ax4 = axes[1, 1]
ax4.plot(time, variance, color='#fb923c', linewidth=1.5)
ax4.axvline(x=time[tip_idx], color='#fbbf24', linewidth=2, linestyle=':', label='Tipping point')
ax4.set_xlabel('Time', color='white')
ax4.set_ylabel('Variance (rolling window)', color='white')
ax4.set_title('Early Warning: Variance', color='white', fontsize=14)
ax4.legend(fontsize=9)
ax4.set_facecolor('#0a0a1a')
ax4.tick_params(colors='white')
ax4.grid(True, alpha=0.2)
ax4.annotate('Rising variance\n= weakening restoring force',
             xy=(time[tip_idx]-5, np.max(variance)*0.7), fontsize=10, color='#fb923c',
             bbox=dict(boxstyle='round,pad=0.3', facecolor='#1a0a00', edgecolor='#fb923c', alpha=0.8))

fig.patch.set_facecolor('#0a0a1a')
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')
plt.close()
print("Tipping point simulation complete.")
print(f"Tipping occurs at time = {time[tip_idx]:.1f}")
print(f"Pre-tip autocorrelation: {np.mean(autocorr[window:tip_idx-100]):.3f}")
print(f"Near-tip autocorrelation: {np.mean(autocorr[tip_idx-200:tip_idx]):.3f}")
print(f"Pre-tip variance: {np.mean(variance[window:tip_idx-100]):.4f}")
print(f"Near-tip variance: {np.mean(variance[tip_idx-200:tip_idx]):.4f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

References

Friedlingstein, P. et al. (2023). Global Carbon Budget 2023. Earth System Science Data, 15(12), 5301–5369.
Canadell, J.G. et al. (2007). Contributions to accelerating atmospheric CO$_2$ growth from economic activity, carbon intensity, and efficiency of natural sinks. Proceedings of the National Academy of Sciences, 104(47), 18866–18870.
Revelle, R. & Suess, H.E. (1957). Carbon dioxide exchange between atmosphere and ocean and the question of an increase of atmospheric CO$_2$ during the past decades. Tellus, 9(1), 18–27.
Farquhar, G.D., von Caemmerer, S. & Berry, J.A. (1980). A biochemical model of photosynthetic CO$_2$ assimilation in leaves of C3 species. Planta, 149(1), 78–90.
Norby, R.J. et al. (2005). Forest response to elevated CO$_2$ is conserved across a broad range of productivity. Proceedings of the National Academy of Sciences, 102(50), 18052–18056.
Bony, S. et al. (2006). How well do we understand and evaluate climate change feedback processes? Journal of Climate, 19(15), 3445–3482.
Sherwood, S.C. et al. (2020). An assessment of Earth’s climate sensitivity using multiple lines of evidence. Reviews of Geophysics, 58(4), e2019RG000678.
Lenton, T.M. et al. (2008). Tipping elements in the Earth’s climate system. Proceedings of the National Academy of Sciences, 105(6), 1786–1793.
Scheffer, M. et al. (2009). Early-warning signals for critical transitions. Nature, 461(7260), 53–59.
Armstrong McKay, D.I. et al. (2022). Exceeding 1.5°C global warming could trigger multiple climate tipping points. Science, 377(6611), eabn7950.

M5. Climate Models & GCMs M7. Sea Level & Ice Sheets

Carbon Cycle & Climate Feedbacks

6.1 The Global Carbon Budget

Airborne Fraction

Carbon Reservoirs

6.2 The Ocean Carbon Sink

Ocean Carbon Chemistry

The Revelle Factor

Ocean Acidification

6.3 The Terrestrial Carbon Sink

CO\(_2\) Fertilization

FACE Experiments

Global Carbon Cycle

Carbon Budget Projection

Global Carbon Budget: Historical Fit & Future Projections

6.4 Climate Feedback Mechanisms

Feedback Factor Derivation

Major Feedback Mechanisms

Feedback Analysis Simulation

Climate Feedback Analysis: ECS & Feedback Factor

6.5 Tipping Points & Early Warning Signals

Major Tipping Elements

Critical Slowing Down

Tipping Point Early Warning Simulation

Tipping Point: Critical Slowing Down & Early Warning Signals

References