Module 5

Climate Models & GCMs

From primitive equations to petascale simulations of Earth's future

5.1 The Primitive Equations

General circulation models (GCMs) solve the primitive equations, which are the full Navier-Stokes equations on a rotating sphere simplified by the hydrostatic approximation (valid for large-scale flow where horizontal scales \(\gg\) vertical scales). The system consists of five coupled equations:

Horizontal Momentum (Zonal & Meridional)

\( \frac{Du}{Dt} - fv + \frac{1}{\rho}\frac{\partial p}{\partial x} = \mathcal{F}_x \)

\( \frac{Dv}{Dt} + fu + \frac{1}{\rho}\frac{\partial p}{\partial y} = \mathcal{F}_y \)

where \(D/Dt = \partial/\partial t + u\partial/\partial x + v\partial/\partial y + \omega\partial/\partial p\)is the material derivative, \(f = 2\Omega\sin\phi\) is the Coriolis parameter, and\(\mathcal{F}\) represents subgrid-scale friction and turbulent mixing.

Hydrostatic Equation

\( \frac{\partial p}{\partial z} = -\rho g \quad \Longleftrightarrow \quad \frac{\partial \Phi}{\partial p} = -\frac{RT}{p} \)

The hydrostatic balance replaces the vertical momentum equation, filtering out sound waves and allowing much larger time steps. This is valid for horizontal scales \(\gtrsim\) 10 km.

Continuity Equation

\( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial \omega}{\partial p} = 0 \)

In pressure coordinates, the continuity equation takes this simple divergence-free form (mass conservation). Here \(\omega = Dp/Dt\) is the pressure vertical velocity.

Thermodynamic Energy Equation

\( c_p \frac{DT}{Dt} - \frac{RT\omega}{p} = Q_{\text{rad}} + Q_{\text{latent}} + Q_{\text{diff}} \)

The first law of thermodynamics, where the left-hand side gives the adiabatic temperature change (including work done by expansion), and the right-hand side contains radiative heating, latent heat release from condensation/precipitation, and diffusive mixing.

Moisture Conservation

\( \frac{Dq}{Dt} = E - C + \mathcal{D}_q \)

where \(q\) is specific humidity, \(E\) is evaporation source,\(C\) is condensation sink, and \(\mathcal{D}_q\) is turbulent moisture diffusion. This equation is critical because water vapour is both the dominant greenhouse gas and the source of latent heating.

5.2 Discretisation Methods

The continuous equations must be discretised for numerical solution. Four major approaches are used in modern climate models:

Finite Difference (GFDL FMS, CESM)

Approximate derivatives on a lat-lon grid: \(\partial u/\partial x \approx (u_{i+1} - u_{i-1})/(2\Delta x)\). Simple to implement but suffers from the pole problem (converging meridians create very small grid cells near the poles).

Spectral (ECMWF IFS, older NCAR models)

Expand fields in spherical harmonics: \(\psi(\lambda,\phi) = \sum_{m,n} \psi_n^m Y_n^m(\lambda,\phi)\). Differentiation becomes multiplication in spectral space. Truncation at wavenumber \(N\) gives resolution \(\sim 180°/N\). T255 \(\approx\) 80 km.

Finite Volume (GFDL FV3, MPAS)

Conserves mass, momentum, and energy by integrating fluxes through cell boundaries. MPAS uses a Voronoi mesh on an icosahedral grid, avoiding the pole problem. FV3 uses a cubed-sphere grid.

Finite Element (FESOM, ICON)

Unstructured triangular meshes that can refine locally (e.g., higher resolution near coastlines). Combines flexibility of finite volumes with higher-order accuracy.

CFL Stability Condition

All explicit time-stepping schemes must satisfy the Courant-Friedrichs-Lewy (CFL)condition for numerical stability. For a wave with phase speed \(c\):

\( \text{CFL} = \frac{c \cdot \Delta t}{\Delta x} \leq 1 \quad \Longrightarrow \quad \Delta t \leq \frac{\Delta x}{c} \)

For climate models, the fastest wave is typically the external gravity wave with\(c \approx \sqrt{gH} \approx 300\) m/s. At \(\Delta x = 100\) km, this gives \(\Delta t \leq 333\) s (\(\sim\)5 min). At convection-permitting resolution (\(\Delta x = 3\) km), \(\Delta t \leq 10\) s, making these simulations enormously expensive.

5.3 Subgrid Parameterisation

Processes occurring below the grid scale (\(\lesssim \Delta x\)) cannot be resolved explicitly and must be parameterised: expressed in terms of resolved-scale variables. This is the largest source of uncertainty in climate projections.

Convection

Individual cumulus clouds (\(\sim\)1\(\text{--}\)10 km) are unresolved at typical GCM resolution. Mass-flux schemes (e.g., Zhang-McFarlane, Tiedtke) represent the collective effect of a population of clouds in each grid cell.

Cloud Microphysics

Processes like droplet formation, ice nucleation, coalescence, and precipitation occur at \(\mu\)m\(\text{--}\)mm scales. Bulk schemes track mass mixing ratios of cloud water, ice, rain, snow, and graupel.

Radiation

Radiative transfer through the atmosphere is computed column-by-column using correlated-k methods (e.g., RRTMGP). Cloud-radiation interaction (overlap assumptions) remains a major source of uncertainty.

Boundary Layer

Turbulent mixing in the lowest 1\(\text{--}\)2 km is parameterised using eddy diffusivity or mass-flux approaches. Surface exchange coefficients (\(C_k\), \(C_d\) as in Module 3) are critical for air-sea interaction.

Gravity Wave Drag

Mountain waves and convectively generated gravity waves transport momentum vertically. Unresolved orographic drag decelerates the flow, affecting jet stream position and blocking frequency.

Convective Adjustment Scheme

The simplest convection scheme: if the lapse rate exceeds the moist adiabatic rate, the column is unstable and is adjusted to neutral buoyancy. The buoyancy of a parcel is:

\( B = g\frac{\theta_v' - \bar{\theta}_v}{\bar{\theta}_v} \)

where \(\theta_v\) is the virtual potential temperature. When \(B > 0\)(unstable), the scheme redistributes heat and moisture vertically until \(B = 0\)(neutral), mimicking the effect of convective overturning. The heating rate from convective adjustment is:

\( \frac{\partial T}{\partial t}\bigg|_{\text{conv}} = \frac{T - T_{\text{adj}}}{\tau_{\text{conv}}} \)

where \(\tau_{\text{conv}} \sim 1\text{--}4\) hours is the convective relaxation timescale.

5.4 CMIP6 & Model Intercomparison

The Coupled Model Intercomparison Project Phase 6 (CMIP6)coordinates climate modelling centres worldwide to run standardised experiments with shared forcing scenarios (SSPs). CMIP6 includes \(\sim\)50 models from \(\sim\)30 centres.

Shared Socioeconomic Pathways (SSPs)

Scenario

Forcing (2100)

Warming (2100)

SSP1-1.9

1.9 W/m\(^2\)

1.0\(\text{--}\)1.8°C

SSP1-2.6

2.6 W/m\(^2\)

1.3\(\text{--}\)2.4°C

SSP2-4.5

4.5 W/m\(^2\)

2.1\(\text{--}\)3.5°C

SSP3-7.0

7.0 W/m\(^2\)

2.8\(\text{--}\)4.6°C

SSP5-8.5

8.5 W/m\(^2\)

3.3\(\text{--}\)5.7°C

Multi-Model Ensembles

Model agreement is a key metric for confidence. When >80% of models agree on the sign of a change (e.g., precipitation increase), confidence is high. The multi-model mean often outperforms individual models because random errors partially cancel. However, structural uncertainties (shared biases) do not cancel. As discussed in ourClimate & Biodiversity course, model spread is particularly large for regional precipitation projections that drive ecosystem impacts.

5.5 Why Fortran Dominates Climate Modelling

Nearly all production climate models (CESM, GFDL, ECMWF IFS, UKESM, MPAS) are written in Fortran. This dominance persists because:

Array operations: Fortran's native multi-dimensional array support maps directly to stencil computations on grids
Compiler optimisation: Fortran's strict aliasing rules allow aggressive optimisation (vectorisation, loop unrolling)
Parallelism: MPI + OpenMP + coarray Fortran enable efficient scaling to \(\sim 10^5\) cores
Legacy code: decades of validated physics packages that would be costly and risky to rewrite
Numerical precision: built-in support for various floating-point precisions critical for conservation properties

5.6 Model Evaluation & Taylor Diagrams

Climate models are evaluated against observations using several metrics:

Root Mean Square Error (RMSE)

\( \text{RMSE} = \sqrt{\frac{1}{N}\sum_{i=1}^N (m_i - o_i)^2} \) where \(m_i\) and \(o_i\) are model and observation at point \(i\)

Bias

\( \text{Bias} = \bar{m} - \bar{o} \) — systematic offset between model mean and observed mean

Correlation

\( R = \frac{\sum(m_i - \bar{m})(o_i - \bar{o})}{\sqrt{\sum(m_i - \bar{m})^2 \sum(o_i - \bar{o})^2}} \) — pattern correlation measures spatial structure

Taylor Diagram Geometry

The Taylor diagram (Taylor, 2001) compactly displays three statistics on a single plot using the mathematical identity:

\( \text{RMSE}'^2 = \sigma_m^2 + \sigma_o^2 - 2\sigma_m \sigma_o R \)

where RMSE\('\) is the centred (bias-removed) RMSE, \(\sigma_m\) and\(\sigma_o\) are the standard deviations of model and observations, and \(R\)is the correlation coefficient. This has the form of the law of cosines:

\( c^2 = a^2 + b^2 - 2ab\cos\theta \)

with \(a = \sigma_m\), \(b = \sigma_o\), \(\theta = \arccos(R)\), and \(c = \text{RMSE}'\). This means we can place \(\sigma_m\) as the radial distance, \(R\) as the azimuthal angle, and read RMSE\('\) as the distance from the observation point. A perfect model sits exactly on the observation point.

5.7 Resolution Hierarchy

Climate modelling spans a vast range of resolutions, each suited to different questions:

CMIP6 GCMs: ~100 km

Multi-century projections, parameterised convection. Grid cells cover ~\(10^4\) km\(^2\). Cannot resolve individual storms but captures large-scale circulation.

Regional models / high-res GCMs: ~25 km

Better representation of orography, coastlines, tropical cyclones (see Module 3). Still requires convective parameterisation.

Convection-permitting: ~1\(\text{--}\)4 km

Explicitly resolves deep convective cells. Critical for extreme precipitation (see Module 4) and urban flooding. Typical runs cover small domains for ~10 year periods.

Large Eddy Simulation (LES): ~10\(\text{--}\)100 m

Resolves turbulent eddies in the boundary layer and individual cloud dynamics. Used for process studies and developing parameterisations for coarser models.

Computational cost scales as \(\propto (\Delta x)^{-3}\) in 3D (doubling resolution requires \(2^3 = 8\times\) more grid points) and \(\propto (\Delta x)^{-4}\)when the CFL condition reduces \(\Delta t\) proportionally. Running a 100-year CMIP6 simulation at convection-permitting resolution would require\(\sim (100/3)^4 \approx 10^6\) times more compute than current practice.

5.8 Ocean Component Models

The ocean component solves similar primitive equations but with additional complexities: the equation of state is nonlinear (UNESCO or TEOS-10), requiring both temperature and salinity as prognostic variables. The Boussinesq approximation is typically applied:

\( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u}\cdot\nabla\mathbf{u} + f\hat{k}\times\mathbf{u} = -\frac{1}{\rho_0}\nabla p + \frac{\rho'}{\rho_0}g\hat{k} + \mathcal{F} \)

\( \frac{\partial T}{\partial t} + \mathbf{u}\cdot\nabla T = \nabla\cdot(K\nabla T) + Q_T \)

\( \frac{\partial S}{\partial t} + \mathbf{u}\cdot\nabla S = \nabla\cdot(K\nabla S) + Q_S \)

Major ocean models include MOM6 (GFDL), NEMO (European), POP (CESM), and MPAS-Ocean. A critical parameterisation is the Gent-McWilliams scheme for mesoscale eddies (\(\sim\)10\(\text{--}\)100 km), which are unresolved at typical 1° ocean resolution. These eddies transport heat poleward and mix tracers, playing a key role in Southern Ocean dynamics and the overturning circulation (seeModule 2).

5.9 Data Assimilation & Reanalysis

Data assimilation optimally combines model predictions with observations to produce the best estimate of the atmospheric state. The most widely used approach is 4D-Var (four-dimensional variational assimilation), which minimises a cost function:

\( J(\mathbf{x}_0) = \frac{1}{2}(\mathbf{x}_0 - \mathbf{x}_b)^T \mathbf{B}^{-1}(\mathbf{x}_0 - \mathbf{x}_b) + \frac{1}{2}\sum_i (\mathbf{y}_i - H_i[\mathbf{x}_i])^T \mathbf{R}_i^{-1}(\mathbf{y}_i - H_i[\mathbf{x}_i]) \)

where \(\mathbf{x}_0\) is the initial state, \(\mathbf{x}_b\) is the background (model forecast), \(\mathbf{B}\) is the background error covariance, \(\mathbf{y}_i\) are observations at time \(t_i\),\(H_i\) is the observation operator, and \(\mathbf{R}_i\) is the observation error covariance.

Reanalysis products (ERA5, JRA-55, MERRA-2) apply data assimilation retrospectively to produce gridded, physically consistent records of the atmosphere spanning decades. ERA5 provides hourly fields at 31 km resolution from 1940 to present, ingesting millions of observations from radiosondes, satellites, aircraft, and surface stations. These products are essential for evaluating climate models and studying extreme events (see Module 4).

Climate Model Architecture

A coupled Earth System Model consists of multiple component models exchanging fluxes through a coupler:

Simulation: 1D Radiative-Convective Equilibrium

A minimal climate model: iterate radiative heating and convective adjustment to find the equilibrium temperature profile. This captures the essential physics of Earth's mean temperature structure:

1D Radiative-Convective Equilibrium Model

Python

Iterating to steady state: grey-atmosphere radiation + convective adjustment

script.py185 lines

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

# 1D Radiative-Convective Equilibrium Model
# Grey atmosphere with convective adjustment

# Parameters
nlevels = 30           # number of atmospheric layers
p_surf = 1013.25       # surface pressure (hPa)
p_top = 10.0           # top of atmosphere (hPa)
S0 = 1361.0            # solar constant (W/m^2)
albedo = 0.30          # planetary albedo
sigma = 5.67e-8        # Stefan-Boltzmann
cp = 1005.0            # specific heat (J/kg/K)
g_val = 9.81           # gravity
tau_total = 2.0        # total IR optical depth
lapse_rate = 6.5e-3    # critical lapse rate (K/m)
dt = 86400.0           # time step (1 day)
H_scale = 8500.0       # scale height (m)

# Pressure levels (layer midpoints)
p_levels = np.linspace(p_surf, p_top, nlevels)
dp = np.abs(np.diff(np.linspace(p_surf, p_top, nlevels + 1))) * 100  # Pa

# Optical depth profile (proportional to pressure)
tau = tau_total * (p_levels / p_surf)

# Height (approximate)
z = -H_scale * np.log(p_levels / p_surf)

# Initial temperature profile (isothermal 250 K)
T = np.full(nlevels, 250.0)
T_surf = 288.0  # surface temperature

# Store history
T_history = [T.copy()]
T_surf_history = [T_surf]

def radiative_heating(T, T_surf, tau, dp):
    """Compute radiative heating rates (grey atmosphere)"""
    # Upward longwave flux at each level
    F_up = np.zeros(nlevels + 1)
    F_up[0] = sigma * T_surf**4  # from surface
    for k in range(nlevels):
        dtau = tau[min(k+1, nlevels-1)] - tau[k] if k < nlevels-1 else tau[k] - tau[max(k-1,0)]
        dtau = max(abs(dtau), 0.001)
        transmit = np.exp(-dtau)
        emit = sigma * T[k]**4 * (1 - transmit)
        F_up[k+1] = F_up[k] * transmit + emit
    
    # Downward longwave flux
    F_down = np.zeros(nlevels + 1)
    F_down[nlevels] = 0  # no incoming LW at TOA
    for k in range(nlevels-1, -1, -1):
        dtau = tau[min(k+1, nlevels-1)] - tau[k] if k < nlevels-1 else tau[k] - tau[max(k-1,0)]
        dtau = max(abs(dtau), 0.001)
        transmit = np.exp(-dtau)
        emit = sigma * T[k]**4 * (1 - transmit)
        F_down[k] = F_down[k+1] * transmit + emit
    
    # Solar absorption (simple: all absorbed at surface)
    F_solar = S0/4 * (1 - albedo)
    
    # Heating rate: dT/dt = -g/(cp) * dF_net/dp
    F_net = F_up - F_down
    dT = np.zeros(nlevels)
    for k in range(nlevels):
        dF = F_net[k+1] - F_net[k]
        dT[k] = -g_val / cp * dF / (dp[k]) * dt  # K per time step
    
    # Surface energy balance
    dT_surf = (F_solar + F_down[0] - sigma * T_surf**4) * dt / (cp * dp[0] / g_val * 0.1)
    
    return dT, dT_surf

def convective_adjustment(T, T_surf, z):
    """Adjust unstable layers to critical lapse rate"""
    T_adj = T.copy()
    # Check from surface up
    T_adj[0] = min(T_adj[0], T_surf - lapse_rate * (z[0]))
    
    for k in range(1, nlevels):
        dz = z[k] - z[k-1]
        T_max = T_adj[k-1] - lapse_rate * dz
        if T_adj[k] < T_max:
            pass  # stable, no adjustment
        else:
            # Unstable: adjust to critical lapse rate
            # Mix to conserve energy
            T_mean = 0.5 * (T_adj[k] + T_adj[k-1])
            T_adj[k-1] = T_mean + lapse_rate * dz / 2
            T_adj[k] = T_mean - lapse_rate * dz / 2
    
    return T_adj

# Iterate to equilibrium
n_steps = 800
for step in range(n_steps):
    # Radiative step
    dT, dT_surf = radiative_heating(T, T_surf, tau, dp)
    T = T + dT * 0.05  # damping factor for stability
    T_surf = T_surf + dT_surf * 0.005
    
    # Convective adjustment
    T = convective_adjustment(T, T_surf, z)
    
    # Enforce physical bounds
    T = np.clip(T, 150, 350)
    T_surf = np.clip(T_surf, 200, 340)
    
    if step % 100 == 0:
        T_history.append(T.copy())
        T_surf_history.append(T_surf)

# Plot results
fig, axes = plt.subplots(1, 3, figsize=(16, 6))
fig.patch.set_facecolor('#0a0a1a')

# Panel 1: Final temperature profile
ax = axes[0]
ax.set_facecolor('#0a0a2e')
ax.plot(T, z/1000, color='#ff6644', linewidth=2.5, label='RCE solution')

# Plot observed-like profile for comparison
T_obs = T_surf - lapse_rate * z
T_obs = np.where(z < 11000, T_obs, 216)  # tropopause
ax.plot(T_obs, z/1000, color='#44ccaa', linewidth=2, linestyle='--', label='Observed approx.')

ax.plot(T_surf, 0, 'o', color='#ffaa44', markersize=10, label=f'Surface: {T_surf:.1f} K')
ax.set_xlabel('Temperature (K)', color='white', fontsize=12)
ax.set_ylabel('Height (km)', color='white', fontsize=12)
ax.set_title('Equilibrium Temperature Profile', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=9)

# Panel 2: Evolution
ax = axes[1]
ax.set_facecolor('#0a0a2e')
colors_ev = plt.cm.plasma(np.linspace(0.2, 0.9, len(T_history)))
for i, (Th, c) in enumerate(zip(T_history, colors_ev)):
    label = f'Day {i*100}' if i < 4 or i == len(T_history)-1 else None
    ax.plot(Th, z/1000, color=c, linewidth=1.5, alpha=0.8, label=label)
ax.set_xlabel('Temperature (K)', color='white', fontsize=12)
ax.set_ylabel('Height (km)', color='white', fontsize=12)
ax.set_title('Approach to Equilibrium', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=8)

# Panel 3: Sensitivity to optical depth (CO2 proxy)
ax = axes[2]
ax.set_facecolor('#0a0a2e')
for tau_val, color, label in [(1.0, '#4488ff', 'tau=1.0 (low CO2)'),
                               (2.0, '#ffaa44', 'tau=2.0 (present)'),
                               (3.0, '#ff4444', 'tau=3.0 (2xCO2)'),
                               (4.0, '#ff66aa', 'tau=4.0 (4xCO2)')]:
    T_test = np.full(nlevels, 250.0)
    T_s_test = 288.0
    tau_test = tau_val * (p_levels / p_surf)
    for step in range(600):
        dT_t, dTs_t = radiative_heating(T_test, T_s_test, tau_test, dp)
        T_test = T_test + dT_t * 0.05
        T_s_test = T_s_test + dTs_t * 0.005
        T_test = convective_adjustment(T_test, T_s_test, z)
        T_test = np.clip(T_test, 150, 350)
        T_s_test = np.clip(T_s_test, 200, 340)
    ax.plot(T_test, z/1000, color=color, linewidth=2, label=f'{label} (Ts={T_s_test:.0f}K)')
    ax.plot(T_s_test, 0, 'o', color=color, markersize=6)

ax.set_xlabel('Temperature (K)', color='white', fontsize=12)
ax.set_ylabel('Height (km)', color='white', fontsize=12)
ax.set_title('Sensitivity to Optical Depth', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=8)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')

print(f"Surface temperature at equilibrium: {T_surf:.1f} K ({T_surf-273.15:.1f} °C)")
print(f"Tropopause temperature: {T[-1]:.1f} K")
print(f"Surface lapse rate: {-(T[1]-T_surf)/(z[1])*1000:.1f} K/km")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Simulation: CFL Condition & Numerical Stability

CFL Condition Demonstration

Python

Showing numerical instability when CFL > 1 and stability when CFL < 1

script.py148 lines

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

def advection_ftcs(nx, nt, c, dx, dt, initial):
    """Forward-time central-space advection scheme.
    du/dt + c * du/dx = 0
    """
    u = initial.copy()
    history = [u.copy()]
    cfl = c * dt / dx
    
    for n in range(nt):
        u_new = u.copy()
        for i in range(1, nx-1):
            u_new[i] = u[i] - 0.5 * cfl * (u[i+1] - u[i-1])
        # Periodic BC
        u_new[0] = u[0] - 0.5 * cfl * (u[1] - u[-2])
        u_new[-1] = u_new[0]
        u = u_new
        if n % (nt//10) == 0:
            history.append(u.copy())
    history.append(u.copy())
    return history, cfl

def advection_upwind(nx, nt, c, dx, dt, initial):
    """Upwind scheme (stable for CFL <= 1)."""
    u = initial.copy()
    history = [u.copy()]
    cfl = c * dt / dx
    
    for n in range(nt):
        u_new = u.copy()
        if c > 0:
            for i in range(1, nx):
                u_new[i] = u[i] - cfl * (u[i] - u[i-1])
            u_new[0] = u[0] - cfl * (u[0] - u[-2])
        u_new[-1] = u_new[0]
        u = u_new
        if n % (nt//10) == 0:
            history.append(u.copy())
    history.append(u.copy())
    return history, cfl

# Setup
nx = 200
L = 1.0
dx = L / nx
c = 1.0  # wave speed
x = np.linspace(0, L, nx)

# Initial condition: Gaussian pulse
initial = np.exp(-200 * (x - 0.3)**2)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.patch.set_facecolor('#0a0a1a')

# Panel 1: Stable (CFL = 0.5, upwind)
ax = axes[0, 0]
ax.set_facecolor('#0a0a2e')
dt_stable = 0.5 * dx / c
hist, cfl = advection_upwind(nx, 200, c, dx, dt_stable, initial)
colors = plt.cm.viridis(np.linspace(0.2, 0.9, len(hist)))
for i, (h, col) in enumerate(zip(hist, colors)):
    ax.plot(x, h, color=col, linewidth=1.5, alpha=0.8)
ax.set_title(f'Upwind Scheme (CFL = {cfl:.2f}) - STABLE', color='#44ccaa', fontsize=12, fontweight='bold')
ax.set_xlabel('x', color='white')
ax.set_ylabel('u', color='white')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.set_ylim(-0.5, 1.5)

# Panel 2: Unstable (CFL = 1.2, upwind)
ax = axes[0, 1]
ax.set_facecolor('#0a0a2e')
dt_unstable = 1.2 * dx / c
hist_u, cfl_u = advection_upwind(nx, 100, c, dx, dt_unstable, initial)
for i, (h, col) in enumerate(zip(hist_u, plt.cm.magma(np.linspace(0.2, 0.9, len(hist_u))))):
    ax.plot(x, np.clip(h, -5, 5), color=col, linewidth=1.5, alpha=0.8)
ax.set_title(f'Upwind Scheme (CFL = {cfl_u:.2f}) - UNSTABLE', color='#ff4444', fontsize=12, fontweight='bold')
ax.set_xlabel('x', color='white')
ax.set_ylabel('u', color='white')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.set_ylim(-5, 5)

# Panel 3: CFL number vs resolution (climate model context)
ax = axes[1, 0]
ax.set_facecolor('#0a0a2e')
dx_range = np.array([200, 100, 50, 25, 10, 5, 3, 1])  # km
c_grav = 300  # gravity wave speed m/s
c_sound = 340  # sound wave speed m/s
c_advect = 50  # typical wind speed m/s

for c_val, label, color in [(c_grav, 'Gravity wave (300 m/s)', '#ff6644'),
                              (c_advect, 'Advection (50 m/s)', '#44ccaa')]:
    dt_max = dx_range * 1000 / c_val  # seconds
    ax.semilogy(dx_range, dt_max, 'o-', color=color, linewidth=2, markersize=6, label=label)

# Reference time steps
ax.axhline(y=1800, color='#ffaa44', linewidth=1, linestyle='--', alpha=0.5)
ax.text(150, 2200, 'Typical GCM dt (30 min)', color='#ffaa44', fontsize=9)
ax.axhline(y=10, color='#ff4444', linewidth=1, linestyle='--', alpha=0.5)
ax.text(150, 13, 'LES dt (10 s)', color='#ff4444', fontsize=9)

ax.set_xlabel('Grid Spacing (km)', color='white', fontsize=11)
ax.set_ylabel('Maximum Time Step (s)', color='white', fontsize=11)
ax.set_title('CFL-Limited Time Step vs Resolution', color='white', fontsize=12, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=9)
ax.invert_xaxis()

# Panel 4: Computational cost scaling
ax = axes[1, 1]
ax.set_facecolor('#0a0a2e')
dx_norm = np.array([200, 100, 50, 25, 10, 5, 3, 1])
# Cost ~ (dx)^-3 for 3D, * (dx)^-1 for CFL = (dx)^-4 total
cost_3d = (100 / dx_norm)**3  # relative to 100 km
cost_4d = (100 / dx_norm)**4  # including CFL constraint

ax.semilogy(dx_norm, cost_3d, 's-', color='#6688ff', linewidth=2, markersize=6, label='Grid points only (~dx^-3)')
ax.semilogy(dx_norm, cost_4d, 'o-', color='#ff6644', linewidth=2, markersize=6, label='Including CFL (~dx^-4)')

# Label key resolutions
for dx_val, label, y_off in [(100, 'CMIP6', 1.5), (25, 'HighRes', 1.5), (3, 'CPM', 1.5)]:
    idx = np.argmin(np.abs(dx_norm - dx_val))
    ax.annotate(label, xy=(dx_norm[idx], cost_4d[idx]),
                xytext=(dx_norm[idx]*0.7, cost_4d[idx]*y_off*5),
                color='#ffaa44', fontsize=10, fontweight='bold',
                arrowprops=dict(arrowstyle='->', color='#ffaa44'))

ax.set_xlabel('Grid Spacing (km)', color='white', fontsize=11)
ax.set_ylabel('Relative Computational Cost', color='white', fontsize=11)
ax.set_title('Cost Scaling with Resolution', color='white', fontsize=12, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=9)
ax.invert_xaxis()

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')

print("CFL condition: dt < dx/c")
print(f"At 100 km, gravity wave: dt < {100000/300:.0f} s ({100000/300/60:.0f} min)")
print(f"At 3 km, gravity wave: dt < {3000/300:.0f} s")
print(f"Cost ratio 3km vs 100km: {(100/3)**4:.0f}x")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Simulation: Spectral Methods in Climate Models

Spectral methods expand fields in basis functions (Fourier modes or spherical harmonics). Differentiation in spectral space is exact (just multiplication by wavenumber), eliminating truncation error from finite differences:

\( f(x) = \sum_{k=-N/2}^{N/2} \hat{f}_k e^{ikx} \quad \Longrightarrow \quad \frac{df}{dx} = \sum_k ik\hat{f}_k e^{ikx} \)

Spectral vs Finite Difference Accuracy

Python

Comparing derivative accuracy: spectral (FFT) vs finite difference on a periodic domain

script.py118 lines

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

def spectral_derivative(f, dx):
    """Compute derivative using FFT (spectral method)"""
    N = len(f)
    k = np.fft.fftfreq(N, d=dx/(2*np.pi))
    f_hat = np.fft.fft(f)
    df_hat = 1j * k * f_hat
    return np.real(np.fft.ifft(df_hat))

def fd_derivative(f, dx):
    """Central finite difference derivative"""
    return (np.roll(f, -1) - np.roll(f, 1)) / (2 * dx)

def fd4_derivative(f, dx):
    """4th order finite difference"""
    return (-np.roll(f, -2) + 8*np.roll(f, -1) - 8*np.roll(f, 1) + np.roll(f, 2)) / (12 * dx)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.patch.set_facecolor('#0a0a1a')

# Test function: sum of sines (smooth)
for ax_idx, (N, title) in enumerate([(32, 'N=32 (coarse)'), (128, 'N=128 (fine)')]):
    ax = axes[0, ax_idx]
    ax.set_facecolor('#0a0a2e')
    
    x = np.linspace(0, 2*np.pi, N, endpoint=False)
    dx = x[1] - x[0]
    f = np.sin(x) + 0.5 * np.sin(3*x) + 0.2 * np.sin(7*x)
    df_exact = np.cos(x) + 1.5 * np.cos(3*x) + 1.4 * np.cos(7*x)
    
    df_spectral = spectral_derivative(f, dx)
    df_fd2 = fd_derivative(f, dx)
    df_fd4 = fd4_derivative(f, dx)
    
    ax.plot(x, df_exact, 'k-', linewidth=3, alpha=0.3, label='Exact')
    ax.plot(x, df_spectral, color='#ff6644', linewidth=2, label='Spectral (FFT)')
    ax.plot(x, df_fd2, color='#4488ff', linewidth=2, linestyle='--', label='FD 2nd order')
    ax.plot(x, df_fd4, color='#44ccaa', linewidth=2, linestyle=':', label='FD 4th order')
    
    ax.set_title(title, color='white', fontsize=13, fontweight='bold')
    ax.set_xlabel('x', color='white')
    ax.set_ylabel('df/dx', color='white')
    ax.tick_params(colors='white')
    ax.grid(True, alpha=0.15, color='#4444aa')
    ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=9)

# Panel 3: Error convergence
ax = axes[1, 0]
ax.set_facecolor('#0a0a2e')
Ns = [8, 16, 32, 64, 128, 256, 512]
errors_spectral = []
errors_fd2 = []
errors_fd4 = []

for N in Ns:
    x = np.linspace(0, 2*np.pi, N, endpoint=False)
    dx = x[1] - x[0]
    f = np.sin(x) + 0.5 * np.sin(3*x) + 0.2 * np.sin(7*x)
    df_exact = np.cos(x) + 1.5 * np.cos(3*x) + 1.4 * np.cos(7*x)
    
    err_s = np.max(np.abs(spectral_derivative(f, dx) - df_exact))
    err_2 = np.max(np.abs(fd_derivative(f, dx) - df_exact))
    err_4 = np.max(np.abs(fd4_derivative(f, dx) - df_exact))
    
    errors_spectral.append(max(err_s, 1e-16))
    errors_fd2.append(err_2)
    errors_fd4.append(err_4)

ax.loglog(Ns, errors_fd2, 's-', color='#4488ff', linewidth=2, markersize=6, label='FD 2nd order (~N^-2)')
ax.loglog(Ns, errors_fd4, 'D-', color='#44ccaa', linewidth=2, markersize=6, label='FD 4th order (~N^-4)')
ax.loglog(Ns, errors_spectral, 'o-', color='#ff6644', linewidth=2, markersize=6, label='Spectral (exponential)')

# Reference slopes
ax.loglog([16, 256], [errors_fd2[1], errors_fd2[1]*(16/256)**2], '--', color='#4488ff', alpha=0.3)
ax.loglog([16, 256], [errors_fd4[1], errors_fd4[1]*(16/256)**4], '--', color='#44ccaa', alpha=0.3)

ax.set_xlabel('Number of Points N', color='white', fontsize=11)
ax.set_ylabel('Max Error', color='white', fontsize=11)
ax.set_title('Error Convergence', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=9)

# Panel 4: Spectral truncation (spherical harmonics analogy)
ax = axes[1, 1]
ax.set_facecolor('#0a0a2e')
x_fine = np.linspace(0, 2*np.pi, 1000)
# Original function with many modes
f_orig = np.sin(x_fine) + 0.5*np.sin(3*x_fine) + 0.3*np.sin(7*x_fine) + 0.15*np.sin(15*x_fine) + 0.08*np.sin(31*x_fine)

ax.plot(x_fine, f_orig, color='#888888', linewidth=1, alpha=0.5, label='Full (all modes)')

for N_trunc, color, label in [(4, '#4488ff', 'T4 (~500 km)'), (8, '#44ccaa', 'T8 (~250 km)'), 
                                (16, '#ffaa44', 'T16 (~125 km)'), (32, '#ff6644', 'T32 (~60 km)')]:
    f_trunc = np.zeros_like(x_fine)
    for k in range(1, N_trunc+1):
        # Reconstruct with truncated modes
        coeff = 2/len(x_fine) * np.sum(f_orig * np.sin(k*x_fine)) * (x_fine[1]-x_fine[0]) / (2*np.pi/len(x_fine))
        f_trunc += coeff * np.sin(k * x_fine)
    ax.plot(x_fine, f_trunc, color=color, linewidth=1.5, label=label)

ax.set_xlabel('Longitude (rad)', color='white', fontsize=11)
ax.set_ylabel('Field Value', color='white', fontsize=11)
ax.set_title('Spectral Truncation (resolution analogy)', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=9)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')

print("Spectral methods achieve machine precision for smooth functions.")
print(f"At N=64: Spectral error = {errors_spectral[3]:.2e}, FD2 error = {errors_fd2[3]:.2e}")
print(f"Spectral is {errors_fd2[3]/max(errors_spectral[3],1e-16):.0f}x more accurate at N=64")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Simulation: Taylor Diagram for Model Evaluation

Taylor Diagram: Comparing Climate Models to Observations

Python

Visualising correlation, standard deviation ratio, and centred RMSE for multiple models

script.py124 lines

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

np.random.seed(42)

# Synthetic model evaluation data
# Each model has: (correlation, normalized std dev, name)
# These represent typical CMIP6 model performance for global temperature
models = [
    {'name': 'CESM2', 'R': 0.97, 'sigma_ratio': 1.05, 'color': '#ff4444'},
    {'name': 'GFDL-ESM4', 'R': 0.95, 'sigma_ratio': 0.92, 'color': '#ff8844'},
    {'name': 'UKESM1', 'R': 0.93, 'sigma_ratio': 1.12, 'color': '#ffaa44'},
    {'name': 'MPI-ESM', 'R': 0.96, 'sigma_ratio': 0.98, 'color': '#44cc88'},
    {'name': 'EC-Earth3', 'R': 0.94, 'sigma_ratio': 1.08, 'color': '#44aacc'},
    {'name': 'ACCESS-ESM', 'R': 0.91, 'sigma_ratio': 0.85, 'color': '#4488ff'},
    {'name': 'CanESM5', 'R': 0.89, 'sigma_ratio': 1.20, 'color': '#8866ff'},
    {'name': 'MIROC6', 'R': 0.92, 'sigma_ratio': 0.95, 'color': '#cc44aa'},
    {'name': 'INM-CM5', 'R': 0.88, 'sigma_ratio': 0.78, 'color': '#aa6644'},
    {'name': 'Multi-model', 'R': 0.98, 'sigma_ratio': 0.96, 'color': '#ffffff'},
]

fig = plt.figure(figsize=(10, 10))
fig.patch.set_facecolor('#0a0a1a')
ax = fig.add_subplot(111, polar=True)
ax.set_facecolor('#0a0a2e')

# Taylor diagram: angle = arccos(R), radius = sigma_ratio
# Reference point (observation) at (0, 1)

# Set up angular axis (correlation)
ax.set_thetamin(0)
ax.set_thetamax(90)
ax.set_theta_direction(-1)  # clockwise
ax.set_theta_offset(np.pi/2)  # start at top

# Correlation labels
corr_labels = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99]
for c in corr_labels:
    angle = np.arccos(c)
    ax.plot([angle, angle], [0, 2.0], color='#444466', linewidth=0.5, alpha=0.5)
    ax.text(angle, 2.05, f'{c:.2f}', color='#8888aa', fontsize=8, ha='center')

# Correlation axis label
ax.text(np.arccos(0.5), 2.3, 'Correlation', color='#aaaacc', fontsize=12,
        fontweight='bold', ha='center', rotation=-60)

# Standard deviation circles
for r in [0.5, 1.0, 1.5, 2.0]:
    theta_range = np.linspace(0, np.pi/2, 100)
    ax.plot(theta_range, np.full_like(theta_range, r), color='#444466', linewidth=0.5, alpha=0.5)

# RMSE circles (centred on observation point at (0, 1))
for rmse_val in [0.25, 0.5, 0.75, 1.0, 1.25]:
    theta_rmse = np.linspace(0, np.pi/2, 200)
    # RMSE'^2 = sigma_m^2 + sigma_o^2 - 2*sigma_m*sigma_o*R
    # With sigma_o = 1: RMSE'^2 = sigma_m^2 + 1 - 2*sigma_m*cos(theta)
    # This is a circle centred at (0, 1) with radius rmse_val
    r_rmse = []
    t_rmse = []
    for theta in theta_rmse:
        # Solve: rmse^2 = r^2 + 1 - 2r*cos(theta)
        # r^2 - 2r*cos(theta) + 1 - rmse^2 = 0
        a_coeff = 1
        b_coeff = -2 * np.cos(theta)
        c_coeff = 1 - rmse_val**2
        disc = b_coeff**2 - 4*a_coeff*c_coeff
        if disc >= 0:
            r1 = (-b_coeff + np.sqrt(disc)) / 2
            r2 = (-b_coeff - np.sqrt(disc)) / 2
            for r_sol in [r1, r2]:
                if 0 < r_sol < 2.0:
                    r_rmse.append(r_sol)
                    t_rmse.append(theta)
    if r_rmse:
        # Sort by angle for clean plotting
        order = np.argsort(t_rmse)
        ax.plot(np.array(t_rmse)[order], np.array(r_rmse)[order],
                color='#446644', linewidth=0.8, linestyle='--', alpha=0.5)
        # Label
        if len(t_rmse) > 0:
            mid = len(t_rmse)//2
            ax.text(np.array(t_rmse)[order][mid], np.array(r_rmse)[order][mid] + 0.05,
                    f'{rmse_val}', color='#668866', fontsize=8)

# Plot observation reference point
ax.plot(0, 1.0, 'k*', markersize=15, color='#ffcc44', zorder=10)
ax.text(0.05, 1.08, 'OBS', color='#ffcc44', fontsize=10, fontweight='bold')

# Plot each model
for m in models:
    theta = np.arccos(m['R'])
    r = m['sigma_ratio']
    marker = '*' if m['name'] == 'Multi-model' else 'o'
    size = 200 if m['name'] == 'Multi-model' else 80
    ax.scatter(theta, r, c=m['color'], s=size, marker=marker, zorder=5,
              edgecolors='white', linewidth=0.5)
    # Label offset
    ax.text(theta + 0.02, r + 0.06, m['name'], color=m['color'],
            fontsize=8, fontweight='bold')

# Axis settings
ax.set_rlabel_position(0)
ax.set_rticks([0.5, 1.0, 1.5, 2.0])
ax.set_yticklabels(['0.5', '1.0', '1.5', '2.0'], color='#8888aa', fontsize=9)
ax.set_rlim(0, 2.0)

# Title
ax.set_title('Taylor Diagram: CMIP6 Model Evaluation\n(Global Mean Temperature)', 
             color='white', fontsize=14, fontweight='bold', pad=30)

# Add labels
fig.text(0.5, 0.02, 'Normalised Standard Deviation (radial) | Correlation (angular) | RMSE\' (green dashed)',
         color='#8899aa', fontsize=10, ha='center')

plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')

print("Taylor Diagram Summary:")
print(f"{'Model':<15} {'Corr':>6} {'Std Ratio':>10} {'RMSE\'':>8}")
for m in models:
    rmse_prime = np.sqrt(m['sigma_ratio']**2 + 1 - 2*m['sigma_ratio']*m['R'])
    print(f"{m['name']:<15} {m['R']:>6.3f} {m['sigma_ratio']:>10.3f} {rmse_prime:>8.3f}")
print("\nMulti-model mean typically closest to observations (lowest RMSE').")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5.8 Emergent Constraints & Climate Sensitivity

Emergent constraints use observable relationships in the present climate to narrow projections of future change. The idea is to find a quantity\(X\) that can be observed and a future projection \(Y\) such that across the model ensemble, \(X\) and \(Y\) are strongly correlated:

\( Y = \alpha + \beta X + \epsilon \)

The observed value \(X_{\text{obs}}\) then constrains \(Y\)through the regression. Key examples include:

Equilibrium Climate Sensitivity (ECS): constrained by subtropical low-cloud feedback, seasonal cycle amplitude, and ENSO variability. IPCC AR6 assessed ECS = 2.5\(\text{--}\)4.0°C (likely range)
Arctic amplification: constrained by observed seasonal sea ice retreat
Land carbon sink: constrained by present-day CO\(_2\) seasonal cycle amplitude
Tropical precipitation: constrained by present-day interannual variability

The formal definition of ECS itself comes from the energy balance:

\( \text{ECS} = \frac{\Delta F_{2\times\text{CO}_2}}{\lambda} = \frac{3.7\;\text{W/m}^2}{\lambda} \)

where \(\lambda\) is the total climate feedback parameter (W/m\(^2\)/K) and \(\Delta F_{2\times\text{CO}_2} \approx 3.7\) W/m\(^2\) is the radiative forcing from CO\(_2\) doubling. The feedback parameter decomposes as:

\( \lambda = \lambda_{\text{Planck}} + \lambda_{\text{WV}} + \lambda_{\text{LR}} + \lambda_{\text{albedo}} + \lambda_{\text{cloud}} \)

The Planck response (\(\lambda_{\text{Planck}} \approx -3.2\) W/m\(^2\)/K) is the stabilising blackbody feedback. Water vapour (\(\sim +1.8\)), lapse rate (\(\sim -0.6\)), and albedo (\(\sim +0.4\)) are well constrained. The cloud feedback (\(\sim +0.5 \pm 0.4\)) remains the dominant source of uncertainty in ECS, hence the importance of cloud parameterisation discussed in Section 5.3.

5.9 Machine Learning in Climate Modelling

Machine learning is increasingly used to augment traditional climate modelling:

ML parameterisation: neural networks trained on high-resolution simulations to replace or augment physics-based parameterisations of convection, radiation, and turbulence. Achieved 10\(\text{--}\)100\(\times\) speedup over traditional schemes
Weather foundation models: Pangu-Weather, GraphCast, FourCastNet produce 10-day forecasts in seconds rather than hours, competitive with ECMWF IFS
Downscaling: super-resolution networks map coarse GCM output to high-resolution regional climate, replacing expensive dynamical downscaling
Emulators: surrogate models that approximate the full GCM response surface, enabling rapid exploration of parameter uncertainty
Bias correction: ML learns systematic errors in model output and corrects them using observational constraints

Challenges remain: ML models lack physical conservation guarantees (mass, energy, momentum), may fail under out-of-distribution scenarios (novel climate states), and their black-box nature makes physical interpretation difficult. Hybrid approaches embedding physical constraints into neural network architectures (physics-informed neural networks, equivariant architectures) are an active research frontier.

A particularly promising direction is using ML to accelerate the most expensive components of climate models. Radiation calculations, which can consume 30\(\text{--}\)50% of total compute, have been successfully emulated with neural networks achieving\(R^2 > 0.999\) and \(100\times\) speedup. Similarly, ML-based ocean mixing parameterisations trained on high-resolution simulations show improved fidelity compared to traditional schemes, with implications for long-term climate drift in coupled models.

5.10 Computational Requirements & Exascale

Climate modelling is among the most computationally demanding scientific applications. Representative costs for a 100-year simulation:

Resolution

Grid Points

Core-Hours (100 yr)

100 km

\(\sim 10^5\)

25 km

\(\sim 10^7\)

3 km

\(\sim 10^9\)

\(\sim 10^{11}\)

100 m (LES)

\(\sim 10^{12}\)

\(\sim 10^{16}\)

The exascale computing era (\(10^{18}\) FLOPS) enables the first global convection-permitting simulations covering decades. Projects like DYAMOND and nextGEMS have demonstrated 2.5\(\text{--}\)5 km global simulations for multi-month periods, revealing storm-resolving features impossible at CMIP6 resolution. GPU-accelerated climate codes are achieving order-of-magnitude speedups, potentially making km-scale century-length projections feasible within the next decade.

References

Washington, W. M. & Parkinson, C. L. (2005). An Introduction to Three-Dimensional Climate Modeling. University Science Books, 2nd edition.

Eyring, V., et al. (2016). Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6). Geoscientific Model Development, 9(5), 1937-1958.

Taylor, K. E. (2001). Summarizing multiple aspects of model performance in a single diagram. Journal of Geophysical Research, 106(D7), 7183-7192.

Courant, R., Friedrichs, K., & Lewy, H. (1928). Über die partiellen Differenzengleichungen der mathematischen Physik. Mathematische Annalen, 100(1), 32-74.

Randall, D. A. (2015). An Introduction to the Global Circulation of the Atmosphere. Princeton University Press.

Satoh, M. (2014). Atmospheric Circulation Dynamics and General Circulation Models. Springer, 2nd edition.

O'Neill, B. C., et al. (2016). The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6. Geoscientific Model Development, 9(9), 3461-3482.

Palmer, T. N. & Stevens, B. (2019). The scientific challenge of understanding and estimating climate change. Proceedings of the National Academy of Sciences, 116(49), 24390-24395.

Sherwood, S. C., et al. (2020). An assessment of Earth's climate sensitivity using multiple lines of evidence. Reviews of Geophysics, 58(4), e2019RG000678.

Lam, R., et al. (2023). Learning skillful medium-range global weather forecasting. Science, 382(6677), 1416-1421.

Stevens, B., et al. (2019). DYAMOND: the DYnamics of the Atmospheric general circulation Modeled On Non-hydrostatic Domains. Progress in Earth and Planetary Science, 6, 61.

Hersbach, H., et al. (2020). The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological Society, 146(730), 1999-2049.