6.4 Climate Models

Climate models are mathematical representations of Earth's climate system that solve the governing equations of fluid dynamics, thermodynamics, and radiative transfer on computational grids. They range from simple zero-dimensional energy balance models to comprehensive Earth System Models coupling atmosphere, ocean, land surface, ice sheets, and biogeochemical cycles. Climate models are indispensable for understanding past climate, attributing observed changes, and projecting future trajectories under different emission scenarios. The Coupled Model Intercomparison Project (CMIP) coordinates international experiments that directly inform IPCC assessments.

Model Hierarchy: EBMs to ESMs

Climate models span a hierarchy of complexity. Simpler models provide conceptual understanding and rapid exploration; comprehensive models capture emergent feedbacks and regional detail.

Energy Balance Models (EBMs)

0D or 1D (latitude). Balance incoming solar against outgoing longwave radiation:

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

where $C$ is heat capacity, $S_0 = 1361$ W/m$^2$, $\alpha \approx 0.30$, and $\epsilon \approx 0.62$. Equilibrium: $T_{eq} \approx 288$ K.

Radiative-Convective Models (RCMs)

1D column models resolving vertical temperature structure. Radiative equilibrium is unstable, so convective adjustment is applied:

$$\frac{\partial T}{\partial t} = -\frac{1}{\rho c_p}\frac{\partial F_{net}}{\partial z} + Q_{conv}$$

$Q_{conv}$ adjusts the lapse rate to the moist adiabat ($\Gamma_m \approx 6.5$ K/km) when static instability occurs. First developed by Manabe and Strickler (1964).

General Circulation Models (GCMs)

3D models solving the primitive equations on a rotating sphere. Coupled atmosphere-ocean GCMs simulate dynamical response to forcing, including ENSO, AMO, and jet stream variability. Typical resolution: 50-100 km horizontal, 30-90 vertical levels. Time steps of 10-30 minutes. Multi-century runs require millions of CPU-hours.

Earth System Models (ESMs)

GCMs extended with interactive carbon cycle, atmospheric chemistry, aerosol microphysics, ice sheet dynamics, and dynamic vegetation. Capture critical carbon-climate feedbacks ($\gamma_L \approx -30$ GtC/K for land, $\gamma_O \approx -8$ GtC/K for ocean). Examples: CESM2, GFDL-ESM4, UKESM1, MPI-ESM1-2.

Governing Equations on the Sphere

GCMs solve the hydrostatic primitive equations in spherical coordinates ($\lambda$ longitude, $\phi$ latitude) with a generalized vertical coordinate $\eta$:

$$\frac{\partial u}{\partial t} = -\mathbf{v}\cdot\nabla u + \left(f + \frac{u\tan\phi}{a}\right)v - \frac{1}{a\cos\phi}\frac{\partial \Phi}{\partial \lambda} + F_u$$

$$\frac{\partial v}{\partial t} = -\mathbf{v}\cdot\nabla v - \left(f + \frac{u\tan\phi}{a}\right)u - \frac{1}{a}\frac{\partial \Phi}{\partial \phi} + F_v$$

$$\frac{\partial T}{\partial t} = -\mathbf{v}\cdot\nabla T + \kappa T\frac{\omega}{p} + \frac{\dot{Q}}{c_p}$$

where $a$ is Earth's radius, $f = 2\Omega\sin\phi$, $\Phi$ is geopotential, $\kappa = R_d/c_p$

Sigma and Hybrid Coordinates

$$p(\eta, \lambda, \phi, t) = A(\eta)\,p_0 + B(\eta)\,p_s(\lambda, \phi, t)$$

Near surface $B \to 1$ (terrain-following); at model top $B \to 0$ (pure pressure). ESMs use 40-90 hybrid levels to ~0.01 hPa.

Spectral vs. Finite-Volume Methods

Spectral Transform

Expand fields in spherical harmonics $Y_n^m(\lambda,\phi)$. Truncation at N gives resolution $\sim 40000/(2N)$ km. Exact differentiation, no pole problem. Used by ECMWF IFS, CESM CAM-SE. Disadvantage: global communication, not locally conservative.

Finite-Volume / Finite-Element

Discretize on cubed-sphere (GFDL FV3), icosahedral (ICON, MPAS), or lat-lon grids. Locally conservative. Scale well on parallel architectures. MPAS uses Voronoi meshes for variable resolution. Growing preference for next-generation models.

Parameterization of Sub-Grid Processes

Convection Schemes

Arakawa-Schubert (1974)

Spectrum of cloud types parameterized by entrainment rate $\lambda$. Quasi-equilibrium closure: $dA/dt \approx 0$ where $A$ is cloud work function. Mass flux: $M_b(\lambda) = -\frac{dA_L/dt}{\partial A/\partial M_b}$.

Zhang-McFarlane (1995)

Bulk mass-flux with CAPE closure: $M_b = \text{CAPE}/\tau_c$ where $\tau_c \sim 2$ hours. Used in CESM/CAM. Organized entrainment: $\epsilon = \epsilon_0(1 + \delta \cdot \text{RH}_{env})$.

Radiation: Two-Stream and Correlated-k

$$\frac{dF^{\uparrow}}{d\tau} = F^{\uparrow} - (1 - \tilde{\omega})B(\tau) - \frac{\tilde{\omega}}{2}(F^{\uparrow} + F^{\downarrow})$$

Correlated-k: Groups spectral lines with similar absorption into k-bins, reducing calculations to ~16 LW and ~14 SW bands (RRTMG). Accuracy within 1 W/m$^2$ of line-by-line.

Cloud Microphysics

Autoconversion: $(\partial q_r/\partial t)_{auto} = \alpha \cdot q_c^{a} \cdot N_c^{b}$. Double-moment schemes (Morrison, Thompson) predict mass and number concentration for each hydrometeor species. Cloud fraction predicted prognostically in modern ESMs.

Grid Resolution and Numerical Stability

$$\Delta t \leq \frac{\Delta x}{c_{max}} \quad \text{(Courant-Friedrichs-Lewy condition)}$$

With $c_{max} \approx 340$ m/s and $\Delta x = 100$ km: $\Delta t \leq 294$ s. Semi-implicit treatment relaxes this to the advective CFL ($\Delta t \sim 900$ s).

Coarse (~100 km)

CMIP6 standard. Resolves synoptic eddies, planetary waves. Parameterizes convection. Multi-century ensembles feasible. Effective resolution ~5-8x grid spacing.

High-Res (~25 km)

HighResMIP. Resolves tropical cyclones, western boundary currents. ~50x compute of coarse. CESM HR, HadGEM3-GC3.1-HH.

Storm-Resolving (~3 km)

DYAMOND project. No cumulus parameterization. ICON 2.5 km, IFS 4 km. Limited to weeks-months. Transformative for cloud feedbacks.

CMIP6 and SSP Scenarios

The Shared Socioeconomic Pathways define future emission trajectories:

SSP1-2.6 (Sustainability)

CO$_2$ peaks ~2030 then declines. Forcing 2.6 W/m$^2$ by 2100. Warming: 1.3-2.4 K. Net-zero by ~2075. Consistent with Paris 2°C target.

SSP2-4.5 (Middle of the Road)

CO$_2$ peaks ~2050 at ~600 ppm. Forcing 4.5 W/m$^2$. Warming: 2.1-3.5 K. Close to current policy trajectory.

SSP3-7.0 (Regional Rivalry)

CO$_2$ ~850 ppm by 2100. Forcing 7.0 W/m$^2$. Warming: 2.8-4.6 K. Fragmented world, high aerosol emissions.

SSP5-8.5 (Fossil-Fueled)

CO$_2$ exceeds 1100 ppm. Forcing 8.5 W/m$^2$. Warming: 3.3-5.7 K. Used for risk assessment; increasingly unlikely baseline.

Equilibrium Climate Sensitivity: $\Delta T_{eq} = -\Delta F / \lambda$ where $\lambda \approx -1.0 \pm 0.5$ W m$^{-2}$ K$^{-1}$. CMIP6 range: 1.8-5.6 K; likely 2.5-4.0 K (AR6).

Budyko-Sellers 1D Energy Balance Model

The 1D EBM (Budyko 1969, Sellers 1969) resolves latitude dependence with meridional heat transport and ice-albedo feedback:

$$C\frac{\partial T}{\partial t} = S(\phi)(1-\alpha(T)) - [A + BT] + D\nabla^2 T$$

$S(\phi)$ is latitude-dependent insolation, $A + BT$ is linearized OLR, $D\nabla^2 T$ is diffusive heat transport

$$\alpha(T) = \begin{cases} \alpha_i = 0.62 & \text{if } T < T_{crit} = -10°\text{C (ice-covered)} \\ \alpha_0 = 0.30 & \text{if } T \geq T_{crit} \text{ (ice-free)} \end{cases}$$

Parameters

A = 202 W/m$^2$ (OLR constant), B = 1.9 W/m$^2$/K (OLR sensitivity), D = 0.44 W/m$^2$/K (diffusion). Empirically calibrated from satellite observations of OLR vs. surface temperature.

Key Results

Multiple equilibria exist (present climate, Snowball Earth, ice-free). Small ice cap instability: if ice line reaches ~30° latitude, runaway glaciation occurs. Hysteresis between entering and exiting Snowball state requires different solar thresholds.

Linearized Climate Response

Linearizing the 0D EBM around equilibrium $T_0$ with $T = T_0 + T'$ gives the transient response to a forcing perturbation:

$$C\frac{dT'}{dt} = -\lambda T' + \Delta F(t), \quad \text{where } \lambda = 4\epsilon\sigma T_0^3 \approx 3.3 \text{ W m}^{-2}\text{ K}^{-1}$$

For step forcing: $T'(t) = \frac{\Delta F}{\lambda}(1 - e^{-t/\tau})$ with e-folding time $\tau = C/\lambda$

Mixed Layer (50 m)

$C \approx 4 \times 10^8$ J/m$^2$/K

$\tau \approx 4$ years

Deep Ocean

$C \approx 4 \times 10^{10}$ J/m$^2$/K

$\tau \approx 400$ years

2xCO$_2$ Forcing

$\Delta F \approx 3.7$ W/m$^2$

$\Delta T_{eq} \approx 1.1$ K (no feedbacks)

Coupled System Components

Modern ESMs couple multiple component models through a coupler framework (OASIS, ESMF/NUOPC). At each coupling timestep (1-6 hours), fluxes of heat, momentum, freshwater, and chemical species are exchanged:

Atmosphere

CAM, IFS, ECHAM, UM. Dynamics + physics parameterizations. Provides wind stress, heat flux, precipitation to ocean/land. Receives SST, sea ice fraction, land surface fluxes.

Ocean

POP, MOM, NEMO, MPAS-Ocean. Resolves thermohaline circulation, mesoscale eddies (at high res). 50-75 vertical levels. Provides SST and heat uptake. Slow component with centennial memory.

Land & Ice

CLM, JSBACH, JULES (land). CICE, SI3 (sea ice). CISM (ice sheets). Land model handles soil moisture, vegetation, carbon. Sea ice model tracks concentration, thickness, dynamics.

Downscaling Methods

Statistical Downscaling

Empirical relationships between large-scale predictors and local variables: $y_{local} = f(\mathbf{X}_{GCM}) + \epsilon$. Methods: regression, weather typing, stochastic generators. Fast but assumes stationarity.

Dynamical Downscaling (RCMs)

High-resolution regional models (WRF, RegCM, COSMO-CLM) nested within GCM output. 10-25 km typical; convection-permitting at 1-4 km. CORDEX coordinates international RCM experiments.

Model Evaluation and Emergent Constraints

Taylor Diagrams

Display correlation, standard deviation ratio, and centered RMSE: $E'^2 = \sigma_f^2 + \sigma_r^2 - 2\sigma_f\sigma_r R$. Compact visualization of multi-model performance against observations.

Emergent Constraints

Observable across-ensemble relationships constraining projections. E.g., tropical low-cloud variability vs. ECS, seasonal CO$_2$ cycle vs. carbon sensitivity. Leverage present-day observables to narrow future uncertainty.

Skill Scores: BSS = $1 - BS/BS_{ref}$, ranked probability skill score, Climate Prediction Index. Models evaluated against ERA5, CERES, GPCP for mean state, variability, trends, and extremes.

Fortran: 2D Shallow Water Model on a Beta-Plane

Solves the linearized shallow water equations on a beta-plane ($f = f_0 + \beta y$) to simulate Rossby wave propagation, demonstrating the key mechanism for planetary wave dynamics in climate models.

program shallow_water_beta_plane
  implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: nx = 128, ny = 64, nsteps = 5000
  real(dp), parameter :: Lx = 4.0e6_dp, Ly = 2.0e6_dp
  real(dp), parameter :: g = 9.81_dp, H0 = 500.0_dp
  real(dp), parameter :: f0 = 1.0e-4_dp, beta = 1.6e-11_dp
  real(dp) :: dx, dy, dt, c_gw
  real(dp) :: h(nx,ny), u(nx,ny), v(nx,ny)
  real(dp) :: h_new(nx,ny), u_new(nx,ny), v_new(nx,ny)
  real(dp) :: f_loc, dhdx, dhdy, dudx, dvdy, x, y, ke, pe
  integer  :: i, j, n, ip, im

  dx = Lx/real(nx,dp); dy = Ly/real(ny,dp)
  c_gw = sqrt(g*H0); dt = 0.25_dp*min(dx,dy)/c_gw
  write(*,'(A,F8.2,A)') 'Rossby radius = ', c_gw/f0/1000.0_dp, ' km'

  u = 0.0_dp; v = 0.0_dp; h = 0.0_dp
  do j = 1, ny
    do i = 1, nx
      x = (i-1)*dx; y = (j-1)*dy - Ly/2.0_dp
      h(i,j) = H0 + 5.0_dp*cos(2.0_dp*3.14159_dp*x/Lx) &
               * exp(-y**2/(2.0e5_dp)**2)
    end do
  end do

  do n = 1, nsteps
    do j = 2, ny-1
      do i = 1, nx
        ip = mod(i,nx)+1; im = mod(i-2+nx,nx)+1
        y = (j-1)*dy - Ly/2.0_dp; f_loc = f0 + beta*y
        dhdx = (h(ip,j)-h(im,j))/(2.0_dp*dx)
        dhdy = (h(i,j+1)-h(i,j-1))/(2.0_dp*dy)
        u_new(i,j) = u(i,j) + dt*(f_loc*v(i,j) - g*dhdx)
        v_new(i,j) = v(i,j) + dt*(-f_loc*u(i,j) - g*dhdy)
      end do
    end do
    u_new(:,1) = u_new(:,2); u_new(:,ny) = u_new(:,ny-1)
    v_new(:,1) = 0.0_dp; v_new(:,ny) = 0.0_dp
    do j = 2, ny-1
      do i = 1, nx
        ip = mod(i,nx)+1; im = mod(i-2+nx,nx)+1
        dudx = (u_new(ip,j)-u_new(im,j))/(2.0_dp*dx)
        dvdy = (v_new(i,j+1)-v_new(i,j-1))/(2.0_dp*dy)
        h_new(i,j) = h(i,j) - dt*H0*(dudx + dvdy)
      end do
    end do
    h_new(:,1) = h_new(:,2); h_new(:,ny) = h_new(:,ny-1)
    u = u_new; v = v_new; h = h_new
    if (mod(n,1000)==0) then
      ke = 0.0_dp; pe = 0.0_dp
      do j = 2, ny-1
        do i = 1, nx
          ke = ke + 0.5_dp*H0*(u(i,j)**2+v(i,j)**2)*dx*dy
          pe = pe + 0.5_dp*g*(h(i,j)-H0)**2*dx*dy
        end do
      end do
      write(*,'(A,I6,A,ES12.5,A,ES12.5)') 'Step ',n,' KE=',ke,' PE=',pe
    end if
  end do
  write(*,'(A)') '--- Height at y=0 cross-section ---'
  do i = 1, nx, 4
    x = (i-1)*dx/1000.0_dp
    write(*,'(I5,F10.1,F12.4)') i, x, h(i,ny/2)
  end do
end program shallow_water_beta_plane

Interactive Simulation: Energy Balance Model (EBM)

Python

Simple 0-D energy balance model: C dT/dt = (1-alpha)S/4 - eps*sigma*T^4 + dF. Solve for temperature response to step CO2 doubling with different ocean heat capacities (fast/slow response). Shows transient warming and ocean heat uptake.

energy_balance_model.py109 lines

import numpy as np
import matplotlib.pyplot as plt

# 0-D Energy Balance Model: C dT/dt = (1-alpha)S/4 - eps*sigma*T^4 + dF
sigma_sb = 5.670e-8  # Stefan-Boltzmann constant
S0 = 1361.0          # Solar constant (W/m^2)
alpha = 0.30         # Planetary albedo
eps = 0.62           # Effective emissivity (accounts for greenhouse effect)
dF_2xCO2 = 3.7       # CO2 doubling forcing (W/m^2)
yr_to_s = 365.25 * 24 * 3600  # Seconds per year

# Total feedback parameter (Planck + all feedbacks)
lambda_total = 1.3   # W/m^2/K (net restoring coefficient)

# Ocean heat capacities (J/m^2/K)
# Mixed layer depth: 50m -> C = rho_w * c_w * h = 1025 * 3850 * 50
C_configs = [
    {'name': 'Fast (50m mixed layer)', 'C': 1025 * 3850 * 50, 'color': '#f87171', 'ls': '-'},
    {'name': 'Medium (200m)', 'C': 1025 * 3850 * 200, 'color': '#fbbf24', 'ls': '--'},
    {'name': 'Slow (500m deep ocean)', 'C': 1025 * 3850 * 500, 'color': '#38bdf8', 'ls': ':'},
    {'name': 'Very slow (1000m)', 'C': 1025 * 3850 * 1000, 'color': '#a78bfa', 'ls': '-.'},
]

# ECS (equilibrium)
T_eq_pre = ((S0 * (1 - alpha) / 4) / (eps * sigma_sb))**0.25
ECS = dF_2xCO2 / lambda_total

# Time integration: 500 years
n_years = 500
dt = 0.1  # years
n_steps = int(n_years / dt)
t = np.linspace(0, n_years, n_steps)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0f172a')
for ax in [ax1, ax2]:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')
    ax.grid(True, color='#334155', alpha=0.4)

print("0-D Energy Balance Model: Step CO2 Doubling Response")
print("=" * 55)
print(f"CO2 doubling forcing: {dF_2xCO2} W/m^2")
print(f"Net feedback parameter: {lambda_total} W/m^2/K")
print(f"Equilibrium warming (ECS): {ECS:.2f} K")
print()
print(f"{'Configuration':<30} {'tau (yr)':>10} {'T(50yr)':>10} {'T(200yr)':>10}")
print("-" * 60)

for cfg in C_configs:
    C = cfg['C']
    T_anom = np.zeros(n_steps)
    # Forward Euler: C dT/dt = dF - lambda * T
    for i in range(1, n_steps):
        dTdt = (dF_2xCO2 - lambda_total * T_anom[i-1]) / C
        T_anom[i] = T_anom[i-1] + dTdt * dt * yr_to_s

# e-folding time
    tau = C / (lambda_total * yr_to_s)  # years
    T_50 = T_anom[int(50 / dt)]
    T_200 = T_anom[int(200 / dt)]
    print(f"{cfg['name']:<30} {tau:>10.1f} {T_50:>10.2f} {T_200:>10.2f}")

ax1.plot(t, T_anom, color=cfg['color'], lw=2.5, ls=cfg['ls'], label=cfg['name'])

# Analytical solution overlay
tau_fast = C_configs[0]['C'] / (lambda_total * yr_to_s)
T_analytical = ECS * (1 - np.exp(-t / tau_fast))
ax1.plot(t, T_analytical, color='#f87171', lw=1, alpha=0.4, ls='-')

ax1.axhline(y=ECS, color='#2dd4bf', ls='--', lw=1.5, label=f'ECS = {ECS:.1f} K')
ax1.set_xlabel('Time (years)', color='#cbd5e1', fontsize=12)
ax1.set_ylabel('Temperature Anomaly \u0394T (K)', color='#cbd5e1', fontsize=12)
ax1.set_title('Transient Response to Step CO\u2082 Doubling', color='#cbd5e1', fontsize=13)
ax1.legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
ax1.set_xlim(0, 500)
ax1.set_ylim(0, ECS * 1.1)

# Panel 2: Energy imbalance and ocean heat uptake
C_med = C_configs[1]['C']
T_med = np.zeros(n_steps)
for i in range(1, n_steps):
    dTdt = (dF_2xCO2 - lambda_total * T_med[i-1]) / C_med
    T_med[i] = T_med[i-1] + dTdt * dt * yr_to_s

N_imbalance = dF_2xCO2 - lambda_total * T_med  # Net TOA imbalance
ocean_heat = np.cumsum(N_imbalance * dt * yr_to_s) / 1e22  # in 10^22 J

ax2.plot(t, N_imbalance, color='#f97316', lw=2.5, label='TOA imbalance N(t)')
ax2b = ax2.twinx()
ax2b.plot(t, ocean_heat, color='#22d3ee', lw=2.5, ls='--', label='Cumulative heat uptake')
ax2b.set_ylabel('Ocean Heat Uptake (10\u00b2\u00b2 J/m\u00b2)', color='#22d3ee', fontsize=11)
ax2b.tick_params(colors='#22d3ee')

ax2.set_xlabel('Time (years)', color='#cbd5e1', fontsize=12)
ax2.set_ylabel('TOA Imbalance N (W/m\u00b2)', color='#f97316', fontsize=11)
ax2.set_title('Energy Imbalance & Ocean Heat Uptake (200m)', color='#cbd5e1', fontsize=13)
ax2.legend(fontsize=9, loc='upper right', facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
ax2b.legend(fontsize=9, loc='center right', facecolor='#1e293b', edgecolor='#334155', labelcolor='#22d3ee')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor=fig.get_facecolor())
print()
print(f"Analytical e-folding time (50m): {tau_fast:.1f} years")
print(f"T(t) = ECS * (1 - exp(-t/tau)) for single-layer model")
print(f"TCR/ECS ratio depends on ocean heat uptake efficiency.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← 6.3 Climate Feedbacks Part 7: Observation & Prediction →