brycewang-stanford/climate-modeling-guide

Climate Modeling Guide

A skill for working with climate models and climate data in research contexts. Covers accessing CMIP archives, processing NetCDF data, running idealized climate simulations, statistical downscaling, and analyzing climate projections with Python tools.

Climate Data Standards

NetCDF and CF Conventions

Climate data is stored in NetCDF (Network Common Data Form) files following CF (Climate and Forecast) conventions:

import xarray as xr
import numpy as np

# Open a CMIP6 temperature dataset
ds = xr.open_dataset("tas_Amon_CESM2_ssp585_r1i1p1f1_gn_201501-210012.nc")

print(ds)
# Dimensions:  (time: 1032, lat: 192, lon: 288)
# Variables:   tas (surface air temperature, K)
# Attributes:  CF-1.6 compliant, CMIP6 metadata

# Basic inspection
print(f"Variable: {ds.tas.long_name}")
print(f"Units: {ds.tas.units}")
print(f"Time range: {ds.time.values[0]} to {ds.time.values[-1]}")
print(f"Spatial resolution: {np.diff(ds.lat.values[:2])[0]:.2f} deg")

CMIP6 Data Access

The Coupled Model Intercomparison Project Phase 6 provides standardized multi-model climate projections:

# Using intake-esm to search the CMIP6 catalog
import intake

# Open the Pangeo CMIP6 catalog (cloud-hosted on Google Cloud)
url = "https://storage.googleapis.com/cmip6/pangeo-cmip6.json"
col = intake.open_esm_datastore(url)

# Search for monthly surface temperature under SSP5-8.5
query = col.search(
    experiment_id="ssp585",
    variable_id="tas",
    table_id="Amon",
    source_id=["CESM2", "GFDL-ESM4", "UKESM1-0-LL", "MPI-ESM1-2-HR"],
    member_id="r1i1p1f1",
)
print(f"Found {len(query)} datasets from {query.nunique()['source_id']} models")

# Load as xarray datasets (lazy, Zarr-backed)
dsets = query.to_dataset_dict(zarr_kwargs={"consolidated": True})

Climate Analysis Techniques

Global Mean Temperature Anomaly

def compute_global_mean_anomaly(ds, baseline_start="1850-01-01",
                                  baseline_end="1900-12-31"):
    """
    Compute area-weighted global mean temperature anomaly
    relative to a baseline period.
    """
    # Area weighting by cosine of latitude
    weights = np.cos(np.deg2rad(ds.lat))
    weights.name = "weights"

    # Weighted global mean time series
    global_mean = ds.tas.weighted(weights).mean(dim=["lat", "lon"])

    # Compute baseline climatology
    baseline = global_mean.sel(time=slice(baseline_start, baseline_end))
    climatology = baseline.groupby("time.month").mean("time")

    # Compute anomalies
    anomaly = global_mean.groupby("time.month") - climatology

    # Annual mean anomaly
    annual_anomaly = anomaly.resample(time="YE").mean()
    return annual_anomaly


def multi_model_ensemble(datasets: dict, baseline_period: tuple):
    """
    Compute multi-model ensemble mean and spread for temperature projections.
    datasets: dict of {model_name: xarray.Dataset}
    Returns ensemble mean and 5th/95th percentile bounds.
    """
    anomalies = []
    for name, ds in datasets.items():
        anom = compute_global_mean_anomaly(ds, *baseline_period)
        anom = anom.assign_coords(model=name)
        anomalies.append(anom)

    ensemble = xr.concat(anomalies, dim="model")
    return {
        "mean": ensemble.mean(dim="model"),
        "p05": ensemble.quantile(0.05, dim="model"),
        "p95": ensemble.quantile(0.95, dim="model"),
    }

Climate Indices

Standard indices used in climate research:

| Index | Full Name | Definition | |-------|-----------|-----------| | ENSO (Nino3.4) | El Nino Southern Oscillation | SST anomaly in 5S-5N, 170W-120W | | NAO | North Atlantic Oscillation | SLP difference Iceland - Azores | | PDO | Pacific Decadal Oscillation | Leading PC of North Pacific SST | | AMO | Atlantic Multidecadal Oscillation | Detrended North Atlantic SST | | IOD | Indian Ocean Dipole | SST difference western - eastern Indian Ocean |

def compute_nino34(sst_dataset, baseline="1991-01-01/2020-12-31"):
    """Compute Nino 3.4 index from SST data."""
    # Select Nino 3.4 region
    nino34_region = sst_dataset.tos.sel(
        lat=slice(-5, 5), lon=slice(190, 240)
    )
    # Area-weighted mean
    weights = np.cos(np.deg2rad(nino34_region.lat))
    nino34_ts = nino34_region.weighted(weights).mean(dim=["lat", "lon"])

    # Remove monthly climatology
    clim = nino34_ts.sel(time=slice(*baseline.split("/"))).groupby("time.month").mean()
    nino34_index = nino34_ts.groupby("time.month") - clim

    # 5-month running mean for standard definition
    nino34_smoothed = nino34_index.rolling(time=5, center=True).mean()
    return nino34_smoothed

Statistical Downscaling

Bias Correction and Spatial Disaggregation

Global climate models (GCMs) typically have 50-200 km resolution, too coarse for impact studies. Statistical downscaling bridges this gap:

def quantile_mapping(obs: np.ndarray, model_hist: np.ndarray,
                     model_future: np.ndarray, n_quantiles: int = 100):
    """
    Quantile mapping bias correction.
    Maps model quantiles to observed quantiles for bias correction.
    """
    quantiles = np.linspace(0, 1, n_quantiles + 1)
    obs_q = np.quantile(obs, quantiles)
    hist_q = np.quantile(model_hist, quantiles)

    # For each future value, find its quantile in historical distribution
    # then map to corresponding observed quantile
    corrected = np.interp(model_future, hist_q, obs_q)
    return corrected

Downscaling Methods Comparison

| Method | Type | Advantages | Limitations | |--------|------|-----------|-------------| | Quantile mapping | Statistical | Simple, preserves distribution | Assumes stationarity | | BCSD | Statistical | Preserves spatial patterns | Limited for extremes | | Delta method | Statistical | Very simple | Only shifts mean | | WRF (dynamical) | Physical | Physically consistent | Computationally expensive | | DeepSD (deep learning) | Hybrid | Learns complex patterns | Requires large training data |

Running Climate Models

Simple Energy Balance Model

def energy_balance_model(S0=1361, albedo=0.30, emissivity=0.612):
    """
    Zero-dimensional energy balance model.
    S0: solar constant (W/m2)
    albedo: planetary albedo
    emissivity: effective atmospheric emissivity
    Returns equilibrium surface temperature (K).
    """
    sigma = 5.67e-8  # Stefan-Boltzmann constant
    # Absorbed solar radiation
    absorbed = S0 * (1 - albedo) / 4
    # Surface temperature with greenhouse effect
    T_surface = (absorbed / (emissivity * sigma)) ** 0.25
    return T_surface

T_eq = energy_balance_model()
print(f"Equilibrium surface temperature: {T_eq:.1f} K ({T_eq - 273.15:.1f} C)")

Tools and Resources

• xarray + dask: Scalable multi-dimensional climate data analysis
• CDO (Climate Data Operators): Command-line NetCDF processing
• NCO (NetCDF Operators): File manipulation and arithmetic
• CESM (Community Earth System Model): Full-complexity coupled GCM
• Pangeo: Cloud-native geoscience data analysis ecosystem
• ESMValTool: Community diagnostic and performance metrics for ESMs
• ClimateData.ca / Copernicus CDS: Processed climate projection portals