Module 8

Climate Projections & Scenarios

SSP scenarios, climate sensitivity, detection & attribution, carbon budgets, and uncertainty

Projecting future climate requires integrating our understanding of climate physics with scenarios of human activity. The Shared Socioeconomic Pathways (SSPs) provide a framework linking socioeconomic development to radiative forcing trajectories. Combined with our knowledge of climate sensitivity and remaining carbon budgets, these projections inform the most consequential policy decisions of the century. Understanding the sources and structure of uncertainty is essential for robust decision-making.

8.1 Shared Socioeconomic Pathways

The SSP-RCP framework combines socioeconomic narratives (SSPs) with radiative forcing levels (RCPs) to create internally consistent scenarios for climate projections. The naming convention SSP$x$-$y$ denotes pathway $x$ with end-of-century forcing $y$ W/m$^2$:

SSP1-2.6 — Sustainability

Low challenges. Rapid clean energy transition, strong international cooperation. CO$_2$ peaks ~2025, net-zero by ~2070. Warming: 1.3–2.4°C by 2100.

SSP2-4.5 — Middle of the Road

Moderate challenges. Current trends continue with slow progress. CO$_2$ peaks ~2050. Warming: 2.1–3.5°C by 2100.

SSP3-7.0 — Regional Rivalry

High challenges. Nationalism, slow technology transfer, high population growth. CO$_2$ doubles by 2100. Warming: 2.8–4.6°C.

SSP5-8.5 — Fossil-fueled Development

High challenges for mitigation. Rapid economic growth powered by fossil fuels. CO$_2$ triples by 2100. Warming: 3.3–5.7°C.

CO$_2$ Concentration from Emission Pathways

Given an emission pathway $E(t)$, the atmospheric CO$_2$ concentration can be computed using an impulse response function (IRF). The IRF describes how a pulse of emissions decays over time:

$\text{IRF}(t) = a_0 + \sum_{i=1}^{3} a_i \exp(-t/\tau_i)$

where $a_0 \approx 0.22$ (permanent airborne fraction),$a_1 \approx 0.28$ with $\tau_1 \approx 394$ yr (deep ocean uptake),$a_2 \approx 0.30$ with $\tau_2 \approx 36$ yr (biosphere + surface ocean),$a_3 \approx 0.20$ with $\tau_3 \approx 4.3$ yr (fast biosphere). The CO$_2$ concentration is then:

$C(t) = C_0 + \frac{1}{2.12} \int_0^t E(t') \cdot \text{IRF}(t - t') \, dt'$

where the factor 2.12 converts Gt CO$_2$ to ppm. The permanent fraction $a_0$ means that ~22% of every CO$_2$ pulse remains in the atmosphere essentially forever (on human timescales), underscoring the importance of emission reductions over carbon removal.

8.2 Equilibrium vs. Transient Climate Sensitivity

Equilibrium Climate Sensitivity (ECS) is the equilibrium warming from a doubling of CO$_2$. First estimated by Charney (1979) as 1.5–4.5°C, the IPCC AR6 assessment narrows this to a likely range of 2.5–4.0°C with best estimate 3.0°C.

Transient Climate Response (TCR) is the warming at the time of CO$_2$ doubling in a 1%/yr CO$_2$ increase experiment (year 70). TCR is lower than ECS because the deep ocean has not reached equilibrium. The IPCC AR6 likely range is 1.4–2.2°C, best estimate 1.8°C.

Relationship: ECS/TCR

The relationship between ECS and TCR can be derived from a two-box energy balance model. The ocean absorbs heat at a rate proportional to the temperature difference between the surface and the deep ocean:

$C_s \frac{dT_s}{dt} = F - \lambda T_s - \kappa(T_s - T_d)$

$C_d \frac{dT_d}{dt} = \kappa(T_s - T_d)$

where $C_s$ and $C_d$ are the heat capacities of the surface and deep ocean layers, $F$ is the radiative forcing,$\lambda$ is the climate feedback parameter, and $\kappa$ is the ocean heat uptake efficiency. At equilibrium ($dT/dt = 0$):

$\text{ECS} = \frac{F_{2\times}}{\lambda} \qquad \text{TCR} = \frac{F_{2\times}}{\lambda + \kappa} \cdot \frac{\kappa + \lambda}{\lambda + \kappa \cdot r}$

To a good approximation, the ratio is:

$\frac{\text{TCR}}{\text{ECS}} \approx \frac{\lambda}{\lambda + \kappa} \approx 0.56$

The ratio depends on the ocean heat uptake efficiency $\kappa$. Higher $\kappa$ (more efficient ocean heat uptake) means a lower TCR/ECS ratio — the transient response is suppressed more by ocean heat storage, but the committed warming is larger.

8.3 Detection & Attribution

How do we know the observed warming is caused by human activities and not natural variability? The answer lies in optimal fingerprinting, a statistical technique that separates the contributions of different forcing agents.

Optimal Fingerprinting Method

The observed climate response $\mathbf{y}$ (a vector of temperatures across space and time) is modeled as a linear combination of forcing patterns (fingerprints) from GCM simulations:

$\mathbf{y} = \sum_i \beta_i \mathbf{x}_i + \boldsymbol{\varepsilon}$

where $\mathbf{x}_i$ are the model-simulated response patterns to individual forcings (greenhouse gases, aerosols, solar, volcanic), $\beta_i$ are scaling factors to be estimated, and $\boldsymbol{\varepsilon}$ is the residual (internal variability).

The scaling factors are estimated using generalized least squares (GLS), which accounts for the correlation structure of internal variability:

$\hat{\boldsymbol{\beta}} = (\mathbf{X}^T \mathbf{C}^{-1} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{C}^{-1} \mathbf{y}$

where $\mathbf{C}$ is the covariance matrix of internal variability (estimated from pre-industrial control simulations). If $\beta_{\text{GHG}}$ is significantly different from zero (and consistent with unity), the greenhouse gas signal is detected. If$\beta_{\text{natural}}$ alone cannot explain the observations, the warming isattributed to human forcing.

The IPCC AR6 states: “It is unequivocal that human influence has warmed the atmosphere, ocean and land.” The human-caused warming of ~1.1°C since pre-industrial is detectable with >99% confidence, while natural forcings contribute approximately zero net warming since 1950.

SSP Scenario Comparison

This diagram compares the four main SSP scenarios across temperature, CO$_2$ concentration, and key socioeconomic characteristics, with an uncertainty cascade from forcing to impacts.

Figure 8.1: SSP scenario comparison and the cascade of uncertainty from emission scenarios to impacts. Near-term projections are dominated by internal variability; long-term by scenario choice.

8.4 Remaining Carbon Budget

The concept of a remaining carbon budget derives from the near-linear relationship between cumulative CO$_2$ emissions and global warming, characterized by the Transient Climate Response to Cumulative Emissions (TCRE):

$\Delta T \approx \text{TCRE} \times \sum E_{\text{CO}_2}$

The TCRE is approximately $1.65 \pm 0.65$ K per 1000 Gt CO$_2$ (IPCC AR6). The remaining carbon budget for a given temperature target is:

$\text{Budget} = \frac{T_{\text{target}} - T_{\text{current}}}{\text{TCRE}} - \text{non-CO}_2 \text{ adjustment}$

For 1.5°C with 50% probability (from start of 2024):

Temperature gap: 1.5 − 1.3 = 0.2°C remaining
Budget from TCRE: ~400 Gt CO$_2$ (after non-CO$_2$ adjustment)
Current emissions: ~40 Gt CO$_2$/yr
Time remaining at current rates: ~10 years

The linearity of TCRE is remarkable and arises from a compensation between carbon cycle and climate feedbacks: as more CO$_2$ is emitted, the ocean and land sinks become less efficient (positive carbon-climate feedback), but the logarithmic forcing-concentration relationship means each additional unit of CO$_2$ produces less forcing. These two effects approximately cancel, producing the near-linear relationship.

Two-Box Climate Model Simulation

This simulation implements a two-box energy balance model that captures both transient (TCR) and equilibrium (ECS) climate sensitivity, computing temperature projections for different SSP scenarios and deriving carbon budgets.

Two-Box Climate Model: TCR, ECS & Carbon Budget

Python

script.py176 lines

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

# ======== Two-box energy balance model ========
# C_s dT_s/dt = F - lambda*T_s - kappa*(T_s - T_d)
# C_d dT_d/dt = kappa*(T_s - T_d)

def two_box_model(forcing, dt=1.0, C_s=8.0, C_d=100.0, lam=1.18, kappa=0.7):
    """Run two-box EBM. Returns surface and deep temperatures."""
    n = len(forcing)
    T_s = np.zeros(n)
    T_d = np.zeros(n)
    for i in range(1, n):
        dTs = (forcing[i] - lam * T_s[i-1] - kappa * (T_s[i-1] - T_d[i-1])) / C_s * dt
        dTd = kappa * (T_s[i-1] - T_d[i-1]) / C_d * dt
        T_s[i] = T_s[i-1] + dTs
        T_d[i] = T_d[i-1] + dTd
    return T_s, T_d

# ======== Forcing scenarios ========
years = np.arange(1850, 2101)
n = len(years)

# Historical forcing (simplified)
F_hist = np.zeros(n)
for i, yr in enumerate(years):
    if yr <= 2020:
        F_hist[i] = 2.5 * (yr - 1850) / 170  # ramp to ~2.5 W/m^2 by 2020
    else:
        F_hist[i] = 2.5  # placeholder

# SSP scenarios (post-2020 forcing)
ssp_params = {
    'SSP1-2.6': {'F_2100': 2.6, 'peak_yr': 2040, 'color': '#4ade80'},
    'SSP2-4.5': {'F_2100': 4.5, 'peak_yr': 2070, 'color': '#fbbf24'},
    'SSP3-7.0': {'F_2100': 7.0, 'peak_yr': 2100, 'color': '#fb923c'},
    'SSP5-8.5': {'F_2100': 8.5, 'peak_yr': 2100, 'color': '#f87171'},
}

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

# ======== Panel 1: Temperature projections ========
ax1 = axes[0, 0]

for name, params in ssp_params.items():
    F = np.copy(F_hist)
    for i, yr in enumerate(years):
        if yr > 2020:
            if yr <= params['peak_yr']:
                F[i] = 2.5 + (params['F_2100'] - 2.5) * (yr - 2020) / (params['peak_yr'] - 2020)
            else:
                F[i] = params['F_2100'] * np.exp(-0.01 * (yr - params['peak_yr']))

T_s, T_d = two_box_model(F)

# Plot from 1950 onward
    mask = years >= 1950
    ax1.plot(years[mask], T_s[mask], color=params['color'], linewidth=2, label=f"{name}: {T_s[-1]:.1f}°C")

ax1.axhline(y=1.5, color='#818cf8', linestyle=':', linewidth=1.5, alpha=0.7, label='1.5°C target')
ax1.axhline(y=2.0, color='#c084fc', linestyle=':', linewidth=1.5, alpha=0.7, label='2.0°C target')
ax1.axvline(x=2024, color='white', linestyle=':', linewidth=1, alpha=0.3)
ax1.set_xlabel('Year', color='white')
ax1.set_ylabel('Global warming (°C above pre-industrial)', color='white')
ax1.set_title('Temperature Projections (Two-Box Model)', color='white', fontsize=14)
ax1.legend(fontsize=8, loc='upper left')
ax1.set_facecolor('#0a0a1a')
ax1.tick_params(colors='white')
ax1.grid(True, alpha=0.2)

# ======== Panel 2: TCR vs ECS ========
ax2 = axes[0, 1]

# 1%/yr CO2 experiment
years_1pct = np.arange(0, 200)
F_1pct = 3.7 * np.log(1.01**years_1pct) / np.log(2)  # forcing from 1%/yr CO2

ecs_values = [2.0, 2.5, 3.0, 3.5, 4.0, 4.5]
colors_ecs = ['#4ade80', '#a3e635', '#fbbf24', '#fb923c', '#f87171', '#ef4444']

for ecs_val, col in zip(ecs_values, colors_ecs):
    lam = 3.7 / ecs_val
    T_s, _ = two_box_model(F_1pct, lam=lam)
    ax2.plot(years_1pct, T_s, color=col, linewidth=1.5, label=f'ECS={ecs_val}°C')
    # Mark TCR (year 70)
    ax2.plot(70, T_s[70], 'o', color=col, markersize=6)

ax2.axvline(x=70, color='white', linestyle='--', linewidth=1, alpha=0.5, label='2xCO₂ (yr 70)')
ax2.set_xlabel('Years of 1%/yr CO$_2$ increase', color='white')
ax2.set_ylabel('Warming (°C)', color='white')
ax2.set_title('TCR vs ECS (1%/yr experiment)', color='white', fontsize=14)
ax2.legend(fontsize=8, ncol=2)
ax2.set_facecolor('#0a0a1a')
ax2.tick_params(colors='white')
ax2.grid(True, alpha=0.2)

# ======== Panel 3: Carbon budget ========
ax3 = axes[1, 0]

# TCRE relationship
cum_emissions = np.linspace(0, 5000, 200)  # Gt CO2
tcre_best = 1.65 / 1000  # K per Gt CO2
tcre_low = 1.0 / 1000
tcre_high = 2.3 / 1000

warming_best = tcre_best * cum_emissions
warming_low = tcre_low * cum_emissions
warming_high = tcre_high * cum_emissions

ax3.fill_between(cum_emissions, warming_low, warming_high, alpha=0.2, color='#818cf8',
                 label='TCRE uncertainty')
ax3.plot(cum_emissions, warming_best, color='#818cf8', linewidth=2.5, label='Best estimate')

# Historical emissions ~2400 Gt CO2
ax3.axvline(x=2400, color='#fbbf24', linestyle=':', linewidth=1.5, label='Historical (2400 Gt CO$_2$)')
ax3.axhline(y=1.5, color='#4ade80', linestyle='--', linewidth=1.5, alpha=0.7)
ax3.axhline(y=2.0, color='#fb923c', linestyle='--', linewidth=1.5, alpha=0.7)

# Budget annotations
budget_15 = (1.5 / tcre_best)
budget_20 = (2.0 / tcre_best)
remaining_15 = budget_15 - 2400
remaining_20 = budget_20 - 2400
ax3.annotate(f'1.5°C budget: {remaining_15:.0f} Gt remaining', xy=(budget_15, 1.5),
             xytext=(3200, 1.2), color='#4ade80', fontsize=10,
             arrowprops=dict(arrowstyle='->', color='#4ade80'))
ax3.annotate(f'2.0°C budget: {remaining_20:.0f} Gt remaining', xy=(budget_20, 2.0),
             xytext=(3200, 2.5), color='#fb923c', fontsize=10,
             arrowprops=dict(arrowstyle='->', color='#fb923c'))

ax3.set_xlabel('Cumulative CO$_2$ emissions (Gt CO$_2$)', color='white')
ax3.set_ylabel('Global warming (°C)', color='white')
ax3.set_title('Carbon Budget: TCRE Relationship', color='white', fontsize=14)
ax3.legend(fontsize=8)
ax3.set_facecolor('#0a0a1a')
ax3.tick_params(colors='white')
ax3.grid(True, alpha=0.2)

# ======== Panel 4: Budget countdown ========
ax4 = axes[1, 1]

current_rate = 40  # Gt CO2/yr
budgets = {
    '1.5°C (50%)': 400,
    '1.5°C (67%)': 275,
    '2.0°C (50%)': 1150,
    '2.0°C (67%)': 900,
}

targets = list(budgets.keys())
budget_vals = list(budgets.values())
years_left = [b / current_rate for b in budget_vals]
colors_budget = ['#4ade80', '#86efac', '#fb923c', '#fdba74']

bars = ax4.barh(targets, years_left, color=colors_budget, edgecolor='white', linewidth=0.5)
for bar, yr, gt in zip(bars, years_left, budget_vals):
    ax4.text(yr + 0.3, bar.get_y() + bar.get_height()/2,
             f'{yr:.0f} yr ({gt} Gt)', va='center', color='white', fontsize=11)

ax4.set_xlabel('Years remaining at current rate (40 Gt CO$_2$/yr)', color='white')
ax4.set_title('Carbon Budget Countdown', color='white', fontsize=14)
ax4.set_facecolor('#0a0a1a')
ax4.tick_params(colors='white')
ax4.grid(True, alpha=0.2, axis='x')

fig.patch.set_facecolor('#0a0a1a')
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')
plt.close()
print("Two-box climate model simulation complete.")
print(f"TCRE best estimate: {tcre_best*1000:.2f} K per 1000 Gt CO2")
print(f"1.5C remaining budget: ~{remaining_15:.0f} Gt CO2 ({remaining_15/current_rate:.0f} years)")
print(f"2.0C remaining budget: ~{remaining_20:.0f} Gt CO2 ({remaining_20/current_rate:.0f} years)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

8.5 Cascade of Uncertainty

Hawkins & Sutton (2009) showed that the relative importance of different uncertainty sources changes with the projection timescale. For global mean temperature:

Near-term (2020–2040): Internal variability dominates. ENSO, PDO, AMO create year-to-year fluctuations of $\pm 0.2$°C that dwarf the inter-scenario differences. All SSPs predict similar warming (~0.5°C above current).
Mid-term (2040–2060): Model uncertainty (climate sensitivity) becomes the dominant factor. The spread in ECS (2.5–4.0°C) produces a range of ~0.5°C even within a single scenario.
Long-term (2060–2100): Scenario uncertainty dominates. The difference between SSP1-2.6 and SSP5-8.5 is ~3°C, far exceeding model uncertainty or internal variability. This means human choices matter most for end-of-century climate.

The total variance can be decomposed as:

$\sigma^2_{\text{total}}(t) = \sigma^2_{\text{scenario}}(t) + \sigma^2_{\text{model}}(t) + \sigma^2_{\text{internal}}(t)$

This framework has profound policy implications: reducing uncertainty in climate sensitivity (through better observations and models) matters most for mid-century projections, while reducing scenario uncertainty (through mitigation commitments) matters most for end-of-century outcomes.

8.6 Overshoot & Carbon Dioxide Removal

Most pathways consistent with 1.5°C involve temperature overshoot — temporarily exceeding the target before net-negative emissions bring temperatures back down. This requires large-scale carbon dioxide removal (CDR).

CDR Technologies

DACCS (Direct Air Carbon Capture & Storage): Chemical sorbents capture CO$_2$ from ambient air. Current energy requirement: ~6–10 GJ per tonne CO$_2$. Cost: $250–600/tonne (declining).
BECCS (Bioenergy with CCS): Biomass absorbs CO$_2$during growth, then combustion with carbon capture creates net-negative emissions. Land requirements are a major constraint.
Enhanced weathering: Spreading crusite minerals (olivine, basalt) accelerates natural CO$_2$ drawdown through silicate weathering:$\text{CaSiO}_3 + \text{CO}_2 \to \text{CaCO}_3 + \text{SiO}_2$.
Ocean alkalinity enhancement: Adding alkaline materials to the ocean increases its capacity to absorb CO$_2$.
Afforestation/reforestation: Most mature and lowest cost, but limited by available land and permanence risks (fire, disease).

Energy Requirements for DACCS

We can derive the thermodynamic minimum energy for separating CO$_2$ from air. The minimum work of separation from a mixture at concentration $c$ is:

$W_{\min} = -RT\left[\ln(c) + \frac{1-c}{c}\ln(1-c)\right]$

For CO$_2$ at 420 ppm ($c = 4.2 \times 10^{-4}$), at $T = 300$ K:

$W_{\min} \approx 20 \text{ kJ/mol CO}_2 \approx 0.45 \text{ GJ/tonne CO}_2$

Real DACCS systems operate at 15–25 times this thermodynamic minimum due to kinetic limitations, sorbent regeneration, and compression for storage. To remove 10 Gt CO$_2$/yr (the scale needed for significant overshoot correction), DACCS would require ~60–100 EJ/yr of energy — roughly 10–17% of current global primary energy supply. This underscores that CDR is a complement to, not a substitute for, emissions reduction.

For ecological and societal implications of climate projections and CDR deployment, see the Climate & Biodiversity course.

Detection, Attribution & Uncertainty Simulation

This simulation demonstrates the optimal fingerprinting method for climate attribution and decomposes the cascade of uncertainty following Hawkins & Sutton (2009).

Detection & Attribution Regression + Uncertainty Cascade

Python

script.py158 lines

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

np.random.seed(42)

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

# ======== Panel 1: Detection & Attribution ========
ax1 = axes[0, 0]

years_da = np.arange(1950, 2024)
n_da = len(years_da)

# True signals (from GCM simulations)
ghg_signal = 0.015 * (years_da - 1950)  # ~1.1 K by 2024
aerosol_signal = -0.004 * (years_da - 1950) * np.exp(-0.01 * (years_da - 1980))
solar_signal = 0.03 * np.sin(2 * np.pi * (years_da - 1950) / 11)
volcanic_signal = np.zeros(n_da)
volcanic_years = [1963, 1982, 1991]  # Agung, El Chichon, Pinatubo
for v_yr in volcanic_years:
    idx = np.argmin(np.abs(years_da - v_yr))
    if idx < n_da:
        volcanic_signal[idx:min(idx+5, n_da)] -= 0.3 * np.exp(-np.arange(min(5, n_da-idx)) / 2)

# Internal variability
internal = 0.12 * np.random.randn(n_da)
# Add low-frequency variability (PDO-like)
internal += 0.08 * np.sin(2 * np.pi * (years_da - 1950) / 30)

# Observations = sum of all
obs = ghg_signal + aerosol_signal + solar_signal + volcanic_signal + internal

# GLS regression (simplified)
X = np.column_stack([ghg_signal, aerosol_signal + solar_signal + volcanic_signal])  # 2 fingerprints
beta = np.linalg.lstsq(X, obs, rcond=None)[0]

ax1.plot(years_da, obs, 'w-', linewidth=1, alpha=0.6, label='Observed')
ax1.plot(years_da, beta[0] * ghg_signal + beta[1] * (aerosol_signal + solar_signal + volcanic_signal),
         color='#f87171', linewidth=2, label=f'All forcings ($\\beta_{{GHG}}$={beta[0]:.2f})')
ax1.plot(years_da, beta[1] * (solar_signal + volcanic_signal),
         color='#60a5fa', linewidth=2, label=f'Natural only ($\\beta_{{nat}}$={beta[1]:.2f})')
ax1.fill_between(years_da, obs - 0.15, obs + 0.15, alpha=0.1, color='white')

ax1.set_xlabel('Year', color='white')
ax1.set_ylabel('Temperature anomaly (°C)', color='white')
ax1.set_title('Detection & Attribution', color='white', fontsize=14)
ax1.legend(fontsize=9, loc='upper left')
ax1.set_facecolor('#0a0a1a')
ax1.tick_params(colors='white')
ax1.grid(True, alpha=0.2)

# ======== Panel 2: Scaling factors with uncertainty ========
ax2 = axes[0, 1]

# Monte Carlo for scaling factor uncertainty
n_mc = 5000
betas_ghg = np.zeros(n_mc)
betas_nat = np.zeros(n_mc)

for i in range(n_mc):
    noise = 0.12 * np.random.randn(n_da)
    obs_mc = ghg_signal + aerosol_signal + solar_signal + volcanic_signal + noise
    b = np.linalg.lstsq(X, obs_mc, rcond=None)[0]
    betas_ghg[i] = b[0]
    betas_nat[i] = b[1]

# Plot as distributions
ax2.hist(betas_ghg, bins=50, alpha=0.6, color='#f87171', density=True, label='$\\beta_{GHG}$')
ax2.hist(betas_nat, bins=50, alpha=0.6, color='#60a5fa', density=True, label='$\\beta_{natural}$')
ax2.axvline(x=1.0, color='white', linestyle='--', linewidth=2, label='Expected = 1.0')
ax2.axvline(x=0.0, color='#fbbf24', linestyle=':', linewidth=1.5, label='Zero')

ax2.set_xlabel('Scaling factor', color='white')
ax2.set_ylabel('Probability density', color='white')
ax2.set_title('Scaling Factor Distributions', color='white', fontsize=14)
ax2.legend(fontsize=9)
ax2.set_facecolor('#0a0a1a')
ax2.tick_params(colors='white')
ax2.grid(True, alpha=0.2)
ax2.annotate(f'GHG detected:\n$\\beta$ = {np.mean(betas_ghg):.2f} ± {np.std(betas_ghg):.2f}',
             xy=(np.mean(betas_ghg), 0), xytext=(1.5, 2.5), fontsize=10, color='#fca5a5',
             arrowprops=dict(arrowstyle='->', color='#f87171'),
             bbox=dict(boxstyle='round,pad=0.3', facecolor='#450a0a', edgecolor='#f87171', alpha=0.8))

# ======== Panel 3: Uncertainty Cascade ========
ax3 = axes[1, 0]

years_proj = np.arange(2020, 2101)
dt = years_proj - 2020

# Variance components (Hawkins & Sutton style)
var_internal = 0.04 * np.ones(len(dt))  # constant ~0.2 K std
var_model = 0.001 * dt**1.2  # grows with time
var_scenario = 0.0001 * dt**2.2  # grows fastest

var_total = var_internal + var_model + var_scenario

# Fractional contributions
f_internal = var_internal / var_total
f_model = var_model / var_total
f_scenario = var_scenario / var_total

ax3.fill_between(years_proj, 0, f_internal, alpha=0.7, color='#4ade80', label='Internal variability')
ax3.fill_between(years_proj, f_internal, f_internal + f_model, alpha=0.7,
                 color='#818cf8', label='Model uncertainty')
ax3.fill_between(years_proj, f_internal + f_model, 1, alpha=0.7,
                 color='#f87171', label='Scenario uncertainty')

ax3.set_xlabel('Year', color='white')
ax3.set_ylabel('Fraction of total variance', color='white')
ax3.set_title('Cascade of Uncertainty (Hawkins & Sutton)', color='white', fontsize=14)
ax3.legend(fontsize=9, loc='center right')
ax3.set_facecolor('#0a0a1a')
ax3.tick_params(colors='white')
ax3.set_ylim(0, 1)
ax3.grid(True, alpha=0.2)

# ======== Panel 4: Temperature fan chart with uncertainty ========
ax4 = axes[1, 1]

# Mean projection (SSP2-4.5 like)
T_mean = 1.3 + 1.5 * (1 - np.exp(-dt / 50))

# Uncertainty bands from three sources
sigma_total = np.sqrt(var_total)

for ci, alpha_val, label in [(1.96, 0.1, '95% CI'), (1.0, 0.25, '68% CI')]:
    ax4.fill_between(years_proj, T_mean - ci * sigma_total, T_mean + ci * sigma_total,
                     alpha=alpha_val, color='#818cf8', label=label)

# Individual scenario lines
for name, dT, color in [('SSP1-2.6', 0.8, '#4ade80'), ('SSP2-4.5', 1.5, '#fbbf24'),
                          ('SSP3-7.0', 2.5, '#fb923c'), ('SSP5-8.5', 3.5, '#f87171')]:
    T_ssp = 1.3 + dT * (1 - np.exp(-dt / 50))
    ax4.plot(years_proj, T_ssp, color=color, linewidth=1.5, alpha=0.8, label=name)

ax4.axhline(y=1.5, color='white', linestyle=':', linewidth=1, alpha=0.5)
ax4.axhline(y=2.0, color='white', linestyle=':', linewidth=1, alpha=0.5)
ax4.set_xlabel('Year', color='white')
ax4.set_ylabel('Global warming (°C)', color='white')
ax4.set_title('Temperature Projections with Uncertainty', color='white', fontsize=14)
ax4.legend(fontsize=7, ncol=2)
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("Detection, attribution & uncertainty analysis complete.")
print(f"GHG scaling factor: {np.mean(betas_ghg):.3f} +/- {np.std(betas_ghg):.3f}")
print(f"Natural scaling factor: {np.mean(betas_nat):.3f} +/- {np.std(betas_nat):.3f}")
print(f"At 2050: internal={f_internal[30]:.0%}, model={f_model[30]:.0%}, scenario={f_scenario[30]:.0%}")
print(f"At 2100: internal={f_internal[-1]:.0%}, model={f_model[-1]:.0%}, scenario={f_scenario[-1]:.0%}")