Feather

Lightweight climate model evaluation framework supporting multiple high-resolution model sets.

Feather compares high-resolution climate models against observations (ERA5, CERES, EN4, ESA-CCI, HadISST, OSI-SAF, PIOMAS/GIOMAS, MSWEP, Berkeley Earth) and CMIP6/HighResMIP multi-model means. It generates diagnostic figures, LLM-analyzed results, a LaTeX report, and a static web dashboard.

Supported model sets

Model set	Models	Resolution	Grid	Data backend
DestinE	IFS-FESOM, IFS-NEMO, ICON	~5 km	HEALPix	intake catalogs
EERIE Ensemble	IFS-FESOM2-SR, IFS-NEMO-ER, ICON-ESM-ER, HadGEM3-GC5	~10-20 km atm, ~5-10 km ocean	0.25° lat/lon	CMOR directory tree
EERIE multi-member	IFS-FESOM2-SR (r1–r3), IFS-NEMO-ER (r1–r3), ICON-ESM-ER, HadGEM3-GC5	~10-20 km atm, ~5-10 km ocean	0.25° lat/lon	CMOR + kerchunk parquet
EERIE 10-member	IFS-FESOM2-SR (r1–r3), IFS-NEMO-ER (r1–r3), ICON-ESM-ER (r1–r3), HadGEM3-GC5	~10-20 km atm, ~5-10 km ocean	0.25° lat/lon	CMOR + kerchunk parquet + ICON kerchunk
IFS-FESOM T319	IFS-FESOM (T319)	~60 km	HEALPix	per-year NetCDF
DestinE GRIB	IFS-FESOM TCO399/TCO319	~25-35 km	lat/lon	GRIB files
TerraDT	IFS-FESOM, IFS-NEMO, ICON	~5 km	HEALPix	intake catalogs
Combined	TCO2999 + TCO399 + TCO319	mixed	mixed	multi-source

The framework is grid-agnostic: diagnostics automatically dispatch between HEALPix and regular lat/lon grids based on per-model configuration. The CompositeModelLoader enables mixing models from different data backends in the same evaluation run.

Installation

From scratch (conda/mamba)

environment.yml uses conda-forge only (nodefaults); no Anaconda defaults channel, so the Anaconda Terms of Service do not apply. A conda-forge-first installer such as Miniforge or micromamba is recommended.

# Create the environment with all dependencies
conda env create -f environment.yml

# Or using mamba/micromamba (faster)
mamba env create -f environment.yml

# Activate and install feather in editable mode
conda activate feather
pip install -e .

Container (Docker / Singularity)

A ready-built image of the latest main is published to ghcr.io/clidyn/feather:latest (every build also as :sha-<commit>):

docker run --rm ghcr.io/clidyn/feather:latest --help                  # Docker
singularity pull docker://ghcr.io/clidyn/feather:latest              # Levante

Build locally with docker build -t feather:latest .. See docker/README.md for bind mounts, API keys and SLURM use.

Into an existing environment

pip install -e .

See environment.yml for the full list of dependencies. Key packages that are best installed via conda-forge: cartopy, healpy, netcdf4, eccodes, cfgrib, regionmask.

regionmask (>= 0.13) supplies the AR6 reference-region polygons used by the per-region Added Value output. It pulls in geopandas and rasterio, which ship their own GDAL, so install it from conda-forge rather than pip to avoid a second GDAL sitting alongside cartopy's.

nereus must be installed directly from GitHub, from main rather than a tag:

pip install git+https://github.com/koldunovn/nereus.git@main

main is required for method="conservative", which feather uses to remap flux and precipitation fields (see Regridding). It was merged upstream after the 0.4.1 tag and is not in the PyPI release. The version string was not bumped, so an installed nereus reports 0.4.1 either way — feather therefore probes the capability at runtime rather than checking the version, and warns and falls back to point interpolation if it is missing. To upgrade an existing environment:

pip install --no-deps --force-reinstall git+https://github.com/koldunovn/nereus.git@main

--no-deps avoids disturbing conda-installed packages; all of nereus's dependencies (including shapely>=2.0, now required) are already present in a working feather environment.

After installation, the feather command is available.

Pipeline

Feather has four pipeline stages that can be run individually or together:

diagnostics  →  analyze  →  report  →  website
(compute+plot)  (Gemini)   (OpenAI)   (HTML)

Run the full pipeline

# DestinE (default config)
feather --config configs/default.yaml -v

# EERIE Ensemble (4 models)
feather --config configs/eerie.yaml -v

# EERIE — 3 IFS-NEMO-ER members + IFS-FESOM2-SR + ICON-ESM-ER + HadGEM3-GC5
feather --config configs/eerie_ifsnemo_members.yaml -v

# EERIE — all 8 members (3×IFS-FESOM2-SR + 3×IFS-NEMO-ER + ICON-ESM-ER + HadGEM3-GC5)
feather --config configs/eerie_all_members.yaml -v

# EERIE — 10 members (as above + ICON-ESM-ER r2/r3), CMIP6 + HighResMIP benchmarks
feather --config configs/eerie_10_mems_cmip6.yaml -v

# EERIE Tropical Nights climate change signal (SSP2-4.5)
feather --config configs/eerie_climchange_tn.yaml \
        --steps diagnostics \
        --diagnostics tropical_nights_change -v

# EERIE Heatwave indices climate change signal (SSP2-4.5)
feather --config configs/eerie_climchange_hw.yaml \
        --steps diagnostics \
        --diagnostics heatwave_change -v

# EERIE Heatwave hotspots — extreme-heat tail-widening vs ERA5 (PNAS Fig 2–4)
feather --config configs/eerie_all_members_extremes.yaml \
        --steps diagnostics \
        --diagnostics heatwave_hotspots --variables tasmax -v

# EERIE daily precipitation extremes (ETCCDI) — 10 members vs daily CMIP6,
# ERA5 and MSWEP; present day + SSP2-4.5 change (derive the daily obs first,
# see "Daily precipitation extremes" below)
feather --config configs/eerie_10_mems_precip_extremes.yaml \
        --steps diagnostics \
        --diagnostics precip_extremes --save-netcdf -v

# EERIE daily temperature extremes (14 ETCCDI indices, land) — 10 members vs
# daily CMIP6 and ERA5; present day + SSP2-4.5 change
feather --config configs/eerie_all_members_indices.yaml \
        --steps diagnostics \
        --diagnostics temp_extremes_indices --save-netcdf -v

# EERIE marine heatwaves (Hobday et al. 2016) in daily SST — 7 members with
# daily SST vs OISST v2.1, ERA5 and daily CMIP6 (derive the daily obs first,
# see "Marine heatwaves" below)
feather --config configs/eerie_10_mems_marine_heatwaves.yaml \
        --steps diagnostics \
        --diagnostics marine_heatwaves --save-netcdf -v

# Subsurface marine heatwaves, shallow (<= 50 m) vs deep, on daily thetao
# model levels (IFS-FESOM2-SR r1, ICON-ESM-ER r1); run after marine_heatwaves
feather --config configs/eerie_10_mems_marine_heatwaves.yaml \
        --steps diagnostics \
        --diagnostics marine_heatwaves_subsurface --save-netcdf -v

# SSH variability (mean + std of zos) vs AVISO — daily (8 members) and
# monthly (10 members + CMIP6); derive the AVISO tree first, see
# "Sea-surface height" below
feather --config configs/eerie_10_mems_ssh.yaml \
        --steps diagnostics \
        --diagnostics ssh_variability --save-netcdf -v

# TerraDT baseline evaluation
feather --config configs/terradt.yaml -v

# DestinE Added Value 1990–2025 (hist + ssp3-7.0 stitched; timeseries to 2049/2044)
feather --config configs/destine_added_value.yaml -v

# Multi-resolution IFS-FESOM comparison (mixed data sources)
feather --config configs/ifs_fesom_combined.yaml -v

Run individual steps

# Generate diagnostic figures (requires compute node for real data)
feather --steps diagnostics -v

# LLM analysis of figures (Gemini)
feather --steps analyze --api-key $GEMINI_API_KEY -v

# Generate LaTeX report (OpenAI)
feather --steps report --openai-api-key $OPENAI_API_KEY -v

# Build static website
feather --steps website -v

# Combine steps
feather --steps analyze report website -v

CLI options

# Filter diagnostics by name
feather --diagnostics global_biases timeseries -v

# Filter by variable (CMOR canonical names)
feather --variables tas pr -v

# Include individual CMIP6 model lines/biases alongside MMM
feather --cmip6-individual -v

# Per-region Added Value. Bare flag = CORDEX-14 bar charts; name a set for
# the others. 'ar6' reads the bias NetCDFs written by global_biases /
# precipitation_mswep, so run those with --save-netcdf first.
feather --diagnostics added_value --added-value-regions -v            # CORDEX-14
feather --diagnostics added_value --added-value-regions ar6 -v        # 58 AR6 regions
feather --diagnostics added_value --added-value-regions cordex14 ar6 -v

# Figures-only website (no LLM analysis)
feather --no-llm -v

# Custom output directory
feather --output ./my_output -v

# Re-run everything (ignore cached results)
feather --no-skip-existing -v

# Also write per-source diagnostic fields to NetCDF (obs, each model,
# CMIP6/HighResMIP MMM). Filenames include the analysis period and
# already-present files are skipped, so it is safe to re-run.
feather --save-netcdf --diagnostics global_biases temperature_berkeley -v

# Custom time period
feather --period 1990 2014 -v

# Compile LaTeX report to PDF
feather --steps report --compile-pdf -v

# Also works as a module
python -m feather --config configs/eerie.yaml -v

NetCDF export (--save-netcdf)

With --save-netcdf, diagnostics also write the per-source fields underlying their figures — the observations, each evaluated model, and each benchmark MMM (CMIP6, HighResMIP) — to NetCDF under {output_dir}/netcdf/{diagnostic}/. Each filename carries the analysis period (e.g. tas_annual_1980-2014.nc), and files already present are skipped, so the flag is incremental and safe to re-run.

Supported by global_biases, temperature_berkeley, precipitation_mswep, global_trends, climate_variability, ocean_sst, sst_hadisst, ocean_en4, radiation_budget, sea_ice, timeseries, seasonal_cycle, and teleconnections. The obs-comparison and extremes/classification diagnostics write NetCDF as part of their normal operation. A new diagnostic can opt in by following the template in NEW_DIAGNOSTIC_SPEC.md.

Ocean Added Value chain. With --save-netcdf (and benchmarks configured), the ocean diagnostics also write model/ensemble/benchmark bias NetCDFs in the shared bias-map schema (ocean_sst→tos vs ESA-CCI, sst_hadisst→tos vs HadISST, ocean_en4→thetao/so vs EN4, sea_ice→siconc vs OSI-SAF), all regridded to a common grid (nereus.resolution). The added_value diagnostic then reuses these to build the Ocean Added Value page — so run the ocean diagnostics first, then added_value last:

CFG=configs/eerie_6_members_cmip6_highresmip.yaml
# 1) Ocean diagnostics write the benchmark-bias NetCDFs
feather --config $CFG --steps diagnostics --save-netcdf -v \
    --diagnostics ocean_sst sst_hadisst ocean_en4 sea_ice
# 2) Added Value reuses them (also runs the atmospheric AV)
feather --config $CFG --steps diagnostics --save-netcdf -v --diagnostics added_value

Statistical significance tests

Bias-map diagnostics attach two area-weighted significance tests to each evaluated source (every model, the CMIP6 MMM, and — with --cmip6-individual — each individual CMIP6 model), comparing its spatial field against the observations:

• Paired t-test — treats each grid point as a paired sample and tests whether the mean bias (model − obs) differs significantly from zero. Cell areas are used as weights so the test matches the area-weighted global-mean bias and does not overrepresent polar cells; the weighted degrees of freedom use Kish's effective sample size. Reported as t_test_statistic / t_test_p_value.
• Variance-ratio F-test — computes F = Var(model) / Var(obs) across grid points (population variance, ddof=0, area-weighted) to test whether the model reproduces the observed spatial variability. F = 1 means equal variability, F > 1 the model is too noisy, F < 1 too smooth; the p-value is two-sided. Reported as variance_ratio / variance_ratio_p_value.

Both statistics are written into each figure's sidecar JSON metadata and surfaced on the dashboard (the website's format_stat filter renders very small p-values in scientific notation). The underlying helpers are spatial_ttest() and spatial_variance_ratio() in feather/util/spatial.py. Currently emitted by global_biases; other bias-map diagnostics can reuse the same helpers.

Python API

from feather import FeatherConfig, run_pipeline

# Load any config
cfg = FeatherConfig.from_yaml("configs/eerie.yaml")

# Full pipeline
result = run_pipeline(cfg, steps="all", api_key="...", openai_api_key="...")

# Individual steps
result = run_pipeline(cfg, steps=["diagnostics"], diagnostics=["global_biases"])
result = run_pipeline(cfg, steps=["report"], openai_api_key="...")

# Result dict
print(result)
# {"figures": 42, "analyses": 42, "syntheses": 12, "report": Path(...), "site_dir": Path(...)}

Available diagnostics

Feather provides 26 registered diagnostics across atmosphere, ocean, cryosphere, extremes, cross-domain evaluation, and model intercomparison:

Atmosphere

Diagnostic	Class	Observation	What it produces
global_biases	GlobalBiases	ERA5	Multi-panel bias maps (obs + all model biases), annual/DJF/JJA
timeseries	TimeseriesDiag	ERA5	Global-mean time series, all models + obs overlaid
seasonal_cycle	SeasonalCycleDiag	ERA5	12-month climatological cycle, all models + obs overlaid
global_trends	GlobalTrends	ERA5	Per-grid-point linear trends over the evaluation period
climate_variability	ClimateVariability	ERA5	STD maps of deseasonalised, detrended monthly fields
radiation_budget	RadiationBudget	CERES EBAF	TOA/surface radiation budget bars, Gregory plot, imbalance time series, bias maps

Ocean

Diagnostic	Class	Observation	What it produces
ocean_sst	OceanSST	ESA-CCI	SST bias maps, time series, seasonal cycle, zonal mean
sst_hadisst	SSTHadISST	HadISST	SST evaluation vs HadISST (1°, full 1980–2014 record vs ESA-CCI's 1990–2014): bias maps, time series, seasonal cycle, zonal mean, warming-trend maps (global + Arctic/Antarctic), and a Taylor diagram, with CMIP6/HighResMIP overlays
ocean_en4	OceanEN4	EN4 v4.2.2	3D ocean T/S: surface bias maps, Hovmoller diagrams, depth-layer time series

Cryosphere

Diagnostic	Class	Observation	What it produces
sea_ice	SeaIceDiag	OSI-SAF, PIOMAS/GIOMAS	Sea ice area/extent/volume time series, seasonal cycles, trends, polar spatial maps

Extremes

Diagnostic	Class	Observation	What it produces
tropical_nights	TropicalNightsDiag	Berkeley Earth Land TMIN	Annual tropical nights count (TN > 20 °C) per grid point; land-only bias map vs Berkeley Earth
tropical_nights_change	TropicalNightsChangeDiag	Berkeley Earth Land TMIN	Climate change signal in TN under SSP2-4.5: [Reference | Future | Change] maps per model, land-mean time series hist+SSP stitched, mean Tmin bias map
heatwave	HeatwaveDiag	Berkeley Earth Land TMAX	Five TX90 heatwave indices (HWN, HWF, HWD, HWM, HWA): climatological maps + annual time series; TMAX bias map
heatwave_change	HeatwaveChangeDiag	Berkeley Earth Land TMAX	Climate change signal in five TX90 heatwave indices under SSP2-4.5: [Reference | Future | Change] maps per model per index, land-mean time series hist+SSP stitched, mean Tmax bias map
heatwave_hotspots	HeatwaveHotspotsDiag	ERA5 (derived daily tasmax)	Extreme-heat tail-widening (trend in yearly 99th − 87.5th percentile of daily tasmax, °C/decade): global trend maps per model with significance stippling and region boxes, regional box-and-whisker, PDF/CDF discrepancy panels (PNAS Fig 2–4)
precip_extremes	PrecipExtremesDiag	ERA5 + MSWEP v2.8 (daily, derived to 0.25°)	Nine ETCCDI daily-precipitation indices (R1mm, R10mm, R20mm, PRCPTOT, SDII, R95p, R99p, Rx1day, Rx5day): reference climatology per model family + daily CMIP6 MMM + obs, SSP2-4.5 change maps, EERIE − CMIP6 difference, bias maps vs each obs dataset, global-mean time series
temp_extremes_change	TempExtremesChangeDiag	— (model change signal)	Climate change signal in mean daily tasmin and tasmax under SSP2-4.5, annual + four seasons (DJF/MAM/JJA/SON): [Reference | Future | Change] maps per model + stitched global-mean time series, per variable
temp_extremes_indices	TempExtremesIndicesDiag	ERA5 (daily tasmax/tasmin, derived to 0.25°)	Fourteen ETCCDI daily-temperature indices over land (TXx, TNx, TXn, TNn, DTR, FD, ID, TX10p, TX90p, TN10p, TN90p, WSDI, CSDI, GSL): reference climatology per model family + daily CMIP6 MMM + ERA5, SSP2-4.5 change maps, EERIE − CMIP6 difference, bias maps vs ERA5, land-mean time series
marine_heatwaves	MarineHeatwavesDiag	NOAA OISST v2.1 + ERA5 (daily SST, derived to 0.25°)	Hobday et al. (2016) marine heatwaves detected per ocean cell in daily SST (de Toma et al. 2026 surface analysis): frequency, MHW days, mean duration, max/mean/cumulative intensity, intensity variability, onset/decline rates, max category — maps per family + CMIP6 MMM, every member, bias vs OISST (grey within 2σ of the member spread), member spread, zonal means, yearly MHW days, violins over the paper's seven areas of interest
marine_heatwaves_subsurface	MarineHeatwavesSubsurfaceDiag	— (no daily subsurface obs; surface OISST/ERA5 as 0 m reference)	Hobday marine heatwaves detected independently on each daily thetao model level down to ~200 m and grouped shallow (≤ 50 m) vs deep (de Toma et al. 2026 §2.4): shallow/deep maps per family, deep − shallow, every member level, depth profiles (global + seven regions, with mixed-layer depth), zonal means, yearly MHW days, shallow/deep violins per region
ssh_variability	SSHVariabilityDiag	AVISO DT2024 allsat L4 adt (daily, derived to 0.25°)	Sea-surface height per grid cell: time mean (minus its area mean) and standard deviation of the raw series, of anomalies (seasonal cycle + trend removed) and of high-pass (< 91 d, mesoscale) anomalies — daily for the 8 members with daily zos, monthly for all 10 + CMIP6 Omon/zos on 1°: maps per family, every member, difference to AVISO, zonal means, regional means over boundary currents

The Tropical Nights Index (TN20) counts nights per year where daily minimum temperature exceeds 20 °C. Requires daily tasmin (CMOR) or kerchunk-parquet mn2t24 store. Per-model NC checkpoints are written to {output_dir}/tropical_nights/.

The Tropical Nights Climate Change Signal (tropical_nights_change) compares mean annual TN between a historical reference period (default 1981–2000) and a future period (default 2031–2050) under SSP2-4.5. It is purpose-built for experiments that span multiple data backends (e.g. CMOR r1 + kerchunk-native r2/r3):

• Group A — Combined map: one row per model × three columns [Reference climatology | Future climatology | Change (Future − Reference)]. Sequential YlOrRd colormap for absolute counts; diverging RdYlBu_r for the change column. Models whose SSP run ends before the future period show a placeholder in columns 2–3.
• Group B — Stitched land-mean time series: historical segment (hist-1950) joined to SSP2-4.5 at 2015 with a dashed vertical line; reference and future windows highlighted; dashed black line shows Berkeley Earth approximate observed TN (months with mean TMIN > 20 °C weighted by days-in-month).
• Group C — Mean daily Tmin bias map vs Berkeley Earth Land TMIN over the reference period (skipped when BE data is unavailable).

Configuration lives under project.climate_change (see configs/eerie_climchange_tn.yaml). Each model declares its hist/future data source (cmor or kerchunk_native) and optionally a future_only_to year. Land masking is applied per-model using the Berkeley Earth land mask (NaN over ocean). Per-model NC checkpoints written to {output_dir}/tropical_nights_change/.

The Heatwave diagnostic computes five interconnected indices using the TX90 method (tasmax > DOY-specific 90th-percentile threshold, runs ≥ 3 consecutive days):

• HWN — heatwave number (events per year)
• HWF — heatwave frequency (heatwave days per year)
• HWD — heatwave duration (length of longest event, days)
• HWM — heatwave magnitude (mean tasmax on heatwave days, °C)
• HWA — heatwave amplitude (peak tasmax on heatwave days, °C)

Summer windows are hemisphere-aware (NH: May–Sep; SH: Nov–Mar). Requires daily tasmax (CMOR day/tasmax table or kerchunk-parquet mx2t24 store). Per-model NC checkpoints stored in {output_dir}/heatwave/.

The Heatwave Climate Change Signal (heatwave_change) compares all five TX90 indices between a historical reference period (default 1981–2000) and a future period (default 2031–2050) under SSP2-4.5. The T90 threshold is computed once from the reference period and applied to both periods. One figure set is produced per index:

• Group A — Combined change map: one row per model × three columns [Reference climatology | Future climatology | Change (Future − Reference)]. Sequential colormap for absolute values; diverging RdYlBu_r for the change column. Models whose SSP run ends before the future period show a placeholder in columns 2–3.
• Group B — Stitched land-mean time series: historical segment joined to SSP2-4.5 at 2015 with a dashed vertical line; reference and future windows highlighted.
• Group C — Mean daily Tmax bias map vs Berkeley Earth Land TMAX over the reference period (skipped when BE data is unavailable).

Configuration lives under project.climate_change (see configs/eerie_climchange_hw.yaml). Each model declares its hist/future data source (cmor or kerchunk_native). Land masking is applied per-model using the Berkeley Earth land mask. Per-model NC checkpoints written to {output_dir}/heatwave_change/.

The Daily Tmin/Tmax Climate Change Signal (temp_extremes_change) compares the mean of daily minimum (tasmin) and maximum (tasmax) temperature between a historical reference period (default project.period) and a future period (default 2031–2050) under SSP2-4.5. Unlike the index-count *_change diagnostics it maps the raw Tmin/Tmax fields themselves, globally (land + ocean), for the annual mean and each of the four seasons. Per variable it produces six figures:

• Change map panels ({var}_{season}_change_panels, ×5: annual + DJF/MAM/JJA/SON) — one row per model × three columns [Reference | Future | Change (Future − Reference)]. cmocean.thermal for the absolute °C climatology columns; diverging RdBu_r for the change column (ΔK ≡ Δ°C). Models with no SSP data in the future window show a placeholder.
• Global-mean time series ({var}_change_timeseries) — historical global-mean annual series joined to SSP2-4.5 at 2015, with reference and future windows shaded and the 2015 boundary marked.

It reuses the multi-backend loader factory of the other *_change diagnostics, so a model can mix CMOR and kerchunk-native segments. To avoid re-reading the daily archive, each model persists its annual- and seasonal-mean series (year, lat, lon) as NetCDF checkpoints {output_dir}/temp_extremes_change/{model}_{var}_{hist,ssp}_{start}_{end}.nc. Configuration lives under project.climate_change (see configs/eerie_all_members_indices.yaml).

The Heatwave Hotspots diagnostic (heatwave_hotspots) reproduces Figures 2–4 of Sambartusek, Kornhuber et al. (2024, PNAS, doi:10.1073/pnas.2411258121) — "A global emergence of regional heatwave hotspots". The core metric is tail-widening: at each land grid point and each year it computes the 99th and 87.5th percentiles of daily tasmax, forms the yearly tail width D = P99 − P87.5, and takes the linear trend of D over the analysis period (°C/decade). A positive trend means the hot tail is widening faster than the bulk of the warm-season distribution — the heatwave-hotspot signature. All fields are regridded to a common 0.25° grid (per-grid interpolator cache), land-masked (>25 %), and masked grey where the P87.5 trend is negative. Five figures are produced:

• Trend map (heatwave_hotspots_trend_map) — global D-trend map per model on Robinson projection, with parametric OLS slope t-test stippling for significance and the ten study-region boxes overlaid.
• Region trend maps (heatwave_hotspots_regions) — zoomed PlateCarree panels for the ten heatwave-hotspot regions (Europe, China, Argentina, Oman, Australia, Japan, Arctic, NW Canada, N Africa, Siberia).
• Box-and-whisker (heatwave_hotspots_boxwhisker) — regional D-trend distributions per region with bootstrap (10 000-sample) confidence intervals on the regional area-weighted mean.
• Discrepancy (heatwave_hotspots_discrepancy) — model-minus-reference cumulative-difference panels highlighting where each model over/under-represents the tail-widening signal.
• PDF (heatwave_hotspots_pdf) — area-weighted probability density of the D-trend across land, model vs ERA5.

The reference is ERA5 (the derived daily tasmax CMOR tree); the model ensemble is the EERIE high-resolution members. CMIP6 is not included in this first pass — daily tasmax is only on the DKRZ pool DRS (/pool/data/CMIP6/.../day/tasmax/), not in the zarr cache — but the config carries a stub for wiring in a yearly-percentile cache later. Per-model yearly-percentile fields are checkpointed to NetCDF in {output_dir}/heatwave_hotspots/ so re-runs skip the heavy daily-percentile pass. Configuration lives in configs/eerie_all_members_extremes.yaml.

With project.ensemble_mode: per_family (set in that config: ICON-ESM-ER runs without a parameterised convection scheme, so its members are not pooled with the IFS runs) the figures that summarise the models as one ensemble are drawn per model family, and their ids gain a _family suffix:

• heatwave_hotspots_regions_family — one thick family-mean line per family over its members as thin lines.
• heatwave_hotspots_discrepancy_family — one column per family; top row: observed trend minus the family's member range (0 inside it), bottom row: observed trend minus the family mean.
• heatwave_hotspots_boxwhisker_family — per region, each member as a dot and the family mean as a bar, one column per family; CMIP6 keeps its box.
• heatwave_hotspots_pdf_family — family-mean densities over the member densities, and one cumulative obs − family curve per family.
• heatwave_hotspots_eerie_vs_cmip6_family — family-mean trend minus the CMIP6 MMM, one panel per family.

Three multi-panel trend-map figures are drawn in both modes, every panel masked like the ERA5 trend map (land, ERA5 P87.5 trend > 0):

• heatwave_hotspots_trend_maps_members — ERA5, each evaluated member and the CMIP6 MMM on one colour scale.
• heatwave_hotspots_trend_maps_family — ERA5, the family means and the CMIP6 MMM on one colour scale.
• heatwave_hotspots_trend_diff_family — the ERA5 trend (own scale), then each family mean minus ERA5 and the CMIP6 MMM minus ERA5 on one difference scale, with mean difference, RMSE, t-test, variance ratio and pattern correlation per panel.

The trend map and the CMIP6-only maps are unchanged. The per-model percentile caches are the same in both modes, so switching mode only replots.

Daily precipitation extremes (precip_extremes)

Computes nine ETCCDI indices from daily pr, following the C3S CMIP6 extreme-indices definitions (wet day: RR ≥ 1 mm):

Index	Meaning	Units
R1mm / R10mm / R20mm	days with RR ≥ 1 / 10 / 20 mm	days/yr
PRCPTOT	total wet-day precipitation	mm/yr
SDII	simple daily intensity (PRCPTOT / R1mm)	mm/day
R95p / R99p	wet-day total above the 95th / 99th wet-day percentile	mm/yr
Rx1day / Rx5day	maximum 1-day / consecutive 5-day precipitation	mm

Periods come from project.climate_change (reference 1981–2000, future SSP2-4.5 2031–2050 in configs/eerie_10_mems_precip_extremes.yaml). precip_extremes.seasons (DJF, MAM, JJA, SON in the EERIE config) adds seasonal versions of every index — the same definitions over the season's days, DJF of year Y being December Y−1 + January–February Y — with the full figure set per season ({idx}_{season}_reference, …). R95p/R99p thresholds are each dataset's own wet-day percentiles over precip_extremes.base_period (default: the reference period — ETCCDI's 1961–1990 is not covered by the 1975-start members) and stay fixed for the future. Indices are computed year by year on each model's native grid; the period means are then interpolated to 0.25° (EERIE family means) or 1° (CMIP6 MMM). The CMIP6 ensemble is auto-discovered from the DKRZ DRS tree (cmip6_daily): one member per model, the first (r1i1p1f1 preferred) that publishes day/pr in both historical and ssp245, under any institution — 35 models as of 2026-10 (of 67 historical models, 16 have no daily pr and 16 no ssp245 daily pr).

Figures per index:

• {idx}_reference — reference-period climatology: ERA5, MSWEP, one panel per EERIE model family (mean of its members), CMIP6 MMM.
• {idx}_change — SSP2-4.5 future minus reference; counts and R95p/R99p as absolute change, PRCPTOT/SDII/Rx1day/Rx5day in percent.
• {idx}_diff_cmip6 — EERIE family minus CMIP6 MMM, reference period.
• {idx}_bias_era5, {idx}_bias_mswep — each EERIE family, the CMIP6 MMM and the other observational dataset minus the obs, with area-weighted mean bias, RMSE, t-test and variance-ratio F-test per panel. The obs-minus-obs panel shows the observational uncertainty, which is large for the intensity indices.
• {idx}_timeseries — global-mean annual index, historical + SSP2-4.5, members thin, family means thick, CMIP6 MMM with min–max range, ERA5/MSWEP lines.

Daily observations. Both references are derived to daily 0.25° (1440×721, kg m⁻² s⁻¹) CMOR trees with area-conservative remapping before the indices are computed, so they are processed exactly like a model:

Dataset	Source on Levante	Script	Output tree
ERA5	/pool/data/ERA5/E5/sf/fc/1D/228 — daily totals, N320 reduced Gaussian	scripts/era5_derive_pr_daily.sh	/work/bm1344/AWI/OBS/era5_derived/CMOR/ECMWF/ERA5/era5/r1i1p1f1/day/pr/gr/v1/
MSWEP v2.8	/pool/data/ICDC/atmosphere/mswep_precipitation/DATA — 0.1°, 3-hourly, one file per step	scripts/mswep_derive_pr_daily.sh	/work/bm1344/AWI/OBS/mswep_derived/CMOR/GloH2O/MSWEP/mswep/r1i1p1f1/day/pr/gr/v1/

The aqua-dvc MSWEP used by precipitation_mswep is monthly only. MSWEP days are the sum of the eight 3-hourly steps 00–21 UTC; a day with a missing step is skipped (treated as missing), never summed short. The datasets are declared under project.precip_extremes.obs (name → data_root, experiment, color); obs_load_period sets the years read (default: the historical load period). Without an obs block the diagnostic runs without observations.

Full workflow:

# 1) Daily 0.25° observations (SLURM, one job per year)
sbatch --array=1981-2014 scripts/era5_derive_pr_daily.sh
sbatch --array=1981-2014%12 scripts/mswep_derive_pr_daily.sh

# 2) Indices and figures (compute node). Daily data are read once; the
#    per-model/per-obs checkpoints make re-runs cheap.
feather --config configs/eerie_10_mems_precip_extremes.yaml --steps diagnostics \
        --diagnostics precip_extremes --save-netcdf -v

# 3) LLM analysis, report, website
feather --config configs/eerie_10_mems_precip_extremes.yaml --steps analyze \
        --diagnostics precip_extremes --api-key $VERTEX_API_KEY -v
feather --config configs/eerie_10_mems_precip_extremes.yaml --steps report \
        --diagnostics precip_extremes --openai-api-key $OPENAI_API_KEY -v
feather --config configs/eerie_10_mems_precip_extremes.yaml --steps website \
        --diagnostics precip_extremes -v

Add --no-skip-existing to analyze and report when figures changed since the last run — otherwise the cached analyses and the cached report structure (publication/) are reused.

Checkpoints (annual index fields, native grid): {output_dir}/precip_extremes/{model}_{hist|ssp}_{start}_{end}.nc, …/cmip6/… and …/obs/{dataset}_hist_{start}_{end}.nc, with the seasonal fields next to each in …_seasons.nc; delete one to force its recomputation. Adding seasons to a finished run keeps the annual files and re-reads the daily data once for the seasons. Open them with decode_timedelta=False — count indices carry units = "days". --save-netcdf writes {output_dir}/netcdf/precip_extremes/{idx}_{ref0}-{fut1}_summary.nc (family/CMIP6/obs fields on 0.25°).

Daily temperature extremes (temp_extremes_indices)

Fourteen ETCCDI indices from daily maximum and minimum 2 m temperature, as defined for the C3S CMIP6 extreme-indices dataset (sis-extreme-indices-cmip6); maths in feather/util/temp_indices.py, diagnostic in feather/diag/temp_extremes_indices.py. It is the temperature counterpart of precip_extremes and shares its periods, family means, CMIP6 MMM, figure set and summary export.

Index	Definition (per year or season)	Units
TXx / TNx	Maximum of daily TX / TN	°C
TXn / TNn	Minimum of daily TX / TN	°C
DTR	Mean diurnal temperature range, mean(TX − TN)	°C
FD / ID	Frost days (TN < 0 °C) / ice days (TX < 0 °C)	days
TX10p / TX90p	Cold / warm days: share of days with TX below its 10th / above its 90th calendar-day percentile	%
TN10p / TN90p	Cold / warm nights, the same for TN	%
WSDI / CSDI	Warm / cold spell duration: days in runs of ≥ 6 days beyond TX90 / TN10 (annual only)	days
GSL	Growing season length: first ≥ 6-day run of TG > 5 °C to the first ≥ 6-day run of TG < 5 °C after 1 July (NH; SH July–June), TG = (TX + TN)/2 (annual only)	days

• Land only. Every index is computed on land cells (land-area fraction > 50 %: the dataset's own sftlf, else the ERA5 mask from obs_datasets.ERA5_SFTLF); maps blank the ocean and the time series are land means. No EERIE member publishes sftlf, so they all use the ERA5 mask.
• Percentiles. Each dataset's own calendar-day 10th/90th percentiles over temp_extremes_indices.base_period (default the reference period, 1981–2000), 5-day window, Hyndman–Fan type 8 as in climdex, no bootstrap, held fixed for the future. Over the base period TX10p…TN90p are therefore ≈ 10 % by construction (slightly less: the in-base inhomogeneity ETCCDI's bootstrap removes), so these four get change maps and time series only — no reference, bias or CMIP6-difference maps.
• Seasons (temp_extremes_indices.seasons): every index except WSDI, CSDI and GSL; DJF of year Y is December Y−1 + January–February Y.
• Members. 10 EERIE members (3 × IFS-FESOM2-SR, 3 × IFS-NEMO-ER, 3 × ICON-ESM-ER, HadGEM3-GC5); HadGEM3-GC5 publishes no daily tasmax/tasmin in either MOHC tree and is skipped. ICON-ESM-ER r2 publishes its daily maximum only on the native R2B8 grid; ICONKerchunkLoader regrids it on the fly with the same CDO conservative weights that produced the published _remap025 stores (they reproduce r2's own remapped daily minimum exactly). A flat model that is also listed under temp_extremes_indices.obs (ERA5 in the EERIE config, where it is a "model" for tropical_nights/heatwave) is treated as observations only.
• CMIP6. cmip6_daily discovery, keeping the first member with day/tasmax and day/tasmin in both historical and ssp245 — 31 models as of 2026-10.
• Observations. ERA5 daily tasmax/tasmin from scripts/era5_derive_tasminmax.sh, processed exactly like a model.

feather --config configs/eerie_all_members_indices.yaml --steps diagnostics \
        --diagnostics temp_extremes_indices --save-netcdf -v

Checkpoints (annual index fields, native grid, ocean NaN): {output_dir}/temp_extremes_indices/{model}_{hist|ssp}_{start}_{end}.nc, …/cmip6/…, …/obs/…, seasonal fields in …_seasons.nc. The percentile thresholds are not stored (≈ 1 GB per 0.25° member); if only an SSP file is missing they are recomputed from the historical run. temp_extremes_indices.percentile_memory_gb (default 12) bounds the base-period block behind them — a 0.25° member needs about 9 GB per variable; above the budget the land cells are processed in groups.

Marine heatwaves (marine_heatwaves)

The surface analysis of de Toma et al. (2026, EGUsphere-2026-5747, code ROHAN), applied to the EERIE members instead of ocean reanalyses. Detection (feather/util/mhw.py) is Hobday et al. (2016) as implemented in Oliver's marineHeatWaves.detect, vectorised over cells and checked against a per-cell port of it: 90th percentile of an 11-day window per calendar day over the climatology period, climatology and threshold smoothed with a 31-day running mean, events ≥ 5 days, gaps ≤ 2 days joined, no detrending. All settings live in project.marine_heatwaves.

• Inputs. Daily tos (table="Oday"): CMOR Oday for the r1 members, the IFS-FESOM2 kerchunk ocean store (avg_tos, now served daily as well as monthly) and ICON-ESM-ER's ocean_native_2d_daily_mean_remap025 store for r2/r3. IFS-NEMO-ER r2/r3 publish monthly ocean output only, and HadGEM3-GC5's daily ocean store holds only zos, so 7 of the 10 members enter. Missing days (FESOM r2/r3 lack about a week) are NaN and never exceed the threshold.
• Observations. scripts/oisst_derive_tos_daily.sh (NOAA OISST v2.1 AVHRR-only from the ICDC copy, 0.25° cell-centred → the models' 721×1440 grid, bilinear) and scripts/era5_derive_tos_daily.sh (ERA5 daily-mean SST, code 34, N320 → 0.25° bilinear, K → °C) write CMOR trees processed exactly like a member. OISST is the reference.
• CMIP6. Daily historical Oday/tos, one member per model (32 models in the 1982–2014 run, 2026-10; r1i1p1f1 where it has Oday/tos), box-averaged from the native ocean grid onto 1° before detection; empty 1° cells take the nearest source point within 1°. Each dataset is thus evaluated at its own resolution.
• Per-cell metrics. Event frequency, % of days in a MHW, mean duration, maximum intensity over all events, mean peak / mean / cumulative intensity, mean within-event intensity standard deviation, mean onset and decline rates, maximum category (intensities relative to the seasonal climatology), plus yearly MHW days and events.
• Masking. Cells where the OISST annual-mean climatology is below ice_mask_sst (−1 °C) are masked in all figures.
• Figures. Per metric: mhw_{m}_maps (obs, family means, CMIP6 MMM), mhw_{m}_members, mhw_{m}_bias (family means, ERA5 and the CMIP6 MMM minus OISST; grey where a multi-member family's bias is within twice its member standard deviation, as in the paper's Fig. 1), mhw_{m}_spread. Plus mhw_zonal_means, mhw_days_timeseries (ocean mean 60°S–60°N) and mhw_regions_{region} (violins over the paper's seven areas of interest).
• Cost. A 0.25° member is ~690k ocean cells × 12,053 days (33 GB of float32). Blocks of latitude rows are read so the block plus workers × cell_chunk temporaries fit memory_gb (default 64 GB, 8 workers, 5000 cells); detection takes about 3 ms per cell for 33 years on one core (≈ 35 min per member, a few minutes on 8 workers), so reading the daily data dominates.

sbatch --array=1982-2014 scripts/oisst_derive_tos_daily.sh
sbatch --array=1982-2014 scripts/era5_derive_tos_daily.sh
feather --config configs/eerie_10_mems_marine_heatwaves.yaml --steps diagnostics \
        --diagnostics marine_heatwaves --save-netcdf -v

Checkpoints (per-cell maps and yearly counts, native grid): {output_dir}/marine_heatwaves/{model}_{start}_{end}.nc, …/obs/…, …/cmip6/… (1°). They store the detection settings and are recomputed when those change.

Subsurface marine heatwaves (marine_heatwaves_subsurface)

The subsurface part of de Toma et al. (2026, §2.4 and §3.3). The same Hobday detection is run independently on every model level of daily thetao (table: HROday) down to max_depth (210 m, so ICON's 201.8 m level is included), and the levels are grouped into shallow (depth ≤ shallow_max_depth, 50 m) and deep. Settings are project.marine_heatwaves overridden by project.marine_heatwaves_subsurface, so detection matches the surface run.

• Inputs. Daily thetao exists for 1982–2014 only for IFS-FESOM2-SR r1 (2.5, 12.5 | 97.5, 195 m) and ICON-ESM-ER r1 (1 | 50.55, 201.8 m). ICON's 50.55 m level falls just below the divide, so the two families sample the column differently; the per-level maps and profiles show this. ICON-ESM-ER r2/r3 load from the new ocean5lev_day kerchunk store (ICONKerchunkLoader, to/so with a depth axis), but it ends in 2000 and r2 lacks 1979, 1981 and 1987, so the 90 % coverage check skips them for 1982–2014. IFS-NEMO, IFS-FESOM2 r2/r3 and HadGEM3 have no daily subsurface temperature.
• No vertical interpolation. Each model is analysed on its own levels.
• Grouping (pool_levels). Event means (duration, intensities, rates) are pooled over all events of the group's levels (per-level means weighted by event frequency); maximum intensity and category take the maximum over levels; frequency, MHW-day share and climatology are the mean over levels.
• Reference. There is no daily subsurface observation product on Levante for 1982–2014 (ICDC ORAS5 is monthly, ARMOR3D covers 2014–15 weekly, GLORYS12 the Arctic only), so this compares the members with each other. The surface run's checkpoints ({output_dir}/marine_heatwaves/[obs/]…, same settings) mark 0 m on the profiles (OISST/ERA5 and each member's tos) and supply the OISST ice mask. Mixed-layer depth comes from Omon/mlotst (mld_table).
• Figures. mhwz_{m}_layers (shallow and deep per family), mhwz_{m}_deep_minus_shallow, mhwz_{m}_levels (every member level), mhwz_profiles_{global|region} (area means vs depth over cells that are ocean on every level of every member, with MLD), mhwz_zonal_means, mhwz_days_timeseries, mhwz_regions_{region} (shallow/deep violins per member, paper Fig. 17). The paper's time–depth sections of named events (Figs. 14–16) are not reproduced: free-running models do not contain the observed events.
• Cost. All requested levels of a latitude-row block are read in one pass (the CMOR files chunk 3–4 levels together). One FESOM member is ~120 GB of float32, hence memory_gb: 200 in the config block — use a ≥ 256 GB node or lower it to read in more blocks.

Checkpoints: {output_dir}/marine_heatwaves_subsurface/{model}_z{depth}m_{start}_{end}.nc per level (only missing levels are recomputed) and {model}_mld_{start}_{end}.nc.

Sea-surface height (ssh_variability)

The notebook /work/bm1344/AWI/EERIE/science/zos_multimodel_std.ipynb (maps of the temporal standard deviation and the mean of daily SSH, models next to AVISO) for the whole EERIE ensemble, over 1993–2014 (AVISO starts in 1993, the historical runs end in 2014). Statistics per grid cell (feather/util/ssh_stats.py):

• mean — time-mean SSH minus its area mean over the cells every dataset covers. Models and AVISO use different reference levels, so only the pattern is compared.
• std_raw — standard deviation of SSH as it is (the notebook's quantity).
• std_anom — after removing the seasonal climatology (calendar-day means smoothed over 31 days; monthly means for monthly data) and a linear trend.
• std_hp — daily only: anomaly minus its centred 91-day running mean (hp_window_days): periods shorter than about 3 months, mostly the mesoscale (but see the ICON caveat below).

Two comparisons, both on the 0.25° grid:

• Daily — 8 members: CMOR HROday for the r1 members, the kerchunk stores for IFS-FESOM2 and ICON-ESM-ER r2/r3, and HadGEM3-GC5 from /work/bm1344/DKRZ/MOHC/HadGEM3-GC5-EERIE-N640-ORCA12 (gr1; the AWI copy has no zos). daily_tables are tried in order (Oday, HROday) and a table returning monthly data is passed over. IFS-NEMO-ER r2/r3 have monthly ocean output only. CMIP6 publishes no daily zos.
• Monthly — all 10 members (Omon) plus the CMIP6 historical Omon/zos ensemble (55 models, one member each), box-averaged from the native ocean grid onto 1°. AVISO is box-averaged onto 1° before its statistics for this comparison: averaging removes eddy variance, so the std of a 1° field is lower than the 1° mean of the 0.25° std.

Observations: scripts/aviso_derive_zos_daily.sh turns the CMEMS DUACS DT2024 allsat L4 files behind the EERIE intake entry dkrz.disk.observations.AVISO.ssh-L4.daily_multi (adt, 0.125°, daily) into Oday and Omon CMOR trees on the models' grid (cdo remapcon, about 17 s per weight build). Cells AVISO covers on fewer than 90 % of the steps (min_valid_fraction; seasonal sea ice) are masked. The gridded AVISO product is smoothed in space and time, which lowers its high-pass variance compared with eddy-rich models. Its absolute dynamic topography also has the dynamic atmospheric correction removed (inverse-barometer response and wind- and pressure-driven barotropic signals shorter than about 20 days); a model whose zos contains these signals shows extra high-frequency variance that is not mesoscale. This is stated in every figure's computation notes. ssh_variability.notes in the config appends free-text caveats to those notes, which the LLM analysis reads.

Figures (f = daily | monthly): ssh_{stat}_{f}_maps (AVISO, family means, and for monthly AVISO at 1° + CMIP6 MMM), ssh_{stat}_{f}_members, ssh_{stat}_{f}_bias (grey where a multi-member family's difference is within twice its member spread; area-mean ratio to AVISO in the statistics), ssh_{f}_zonal_means, ssh_{f}_regions (Gulf Stream, Kuroshio Extension, Agulhas Retroflection, Brazil–Malvinas Confluence, East Australian Current, Drake Passage, equatorial Pacific, ocean 60°S–60°N).

sbatch --array=1993-2014 scripts/aviso_derive_zos_daily.sh
feather --config configs/eerie_10_mems_ssh.yaml --steps diagnostics \
        --diagnostics ssh_variability --save-netcdf -v

Checkpoints: {output_dir}/ssh_variability/[obs/|cmip6/]{name}_{daily|monthly}_{start}_{end}.nc. They hold the statistics on each dataset's own grid, so a rerun after a plotting failure only redraws the figures.

Results of the 2026-10-13 run (area means over the cells all datasets cover, cm):

	AVISO	IFS-FESOM2-SR	IFS-NEMO-ER	ICON-ESM-ER	HadGEM3-GC5	CMIP6 MMM (1°)
daily raw std	7.55	7.22	6.94	9.94	7.90	—
daily anomaly std	6.30	6.08	5.84	8.19	6.67	—
daily high-pass std	3.69	3.90	3.73	6.62	4.33	—
monthly raw std	6.90 (1°: 6.55)	6.28	6.11	7.52	6.86	5.31
monthly anomaly std	5.57 (1°: 5.13)	4.98	4.88	5.35	5.46	3.58

Members of a family agree closely (daily high-pass std: FESOM 3.88/3.91/3.91 cm, ICON 6.60/6.62/6.64 cm), although r1 is read from CMOR and r2/r3 from kerchunk stores.

ICON-ESM-ER caveat. ICON's daily high-pass std is about 1.8 times AVISO's, with the excess spread over whole basins in the Southern Ocean and the high-latitude North Pacific and Atlantic (Drake Passage: 17.3 vs 9.7 cm raw daily std), not concentrated in the boundary currents. Monthly, ICON is close to AVISO. This points to sub-monthly barotropic or inverse-barometer signals in ICON's zos that AVISO removes, not to excess eddy activity. It is not yet verified; the check is the regression of ICON's daily high-pass zos on psl (about −1 cm/hPa for an inverse-barometer response).

HadGEM3-GC5's gr1 latitudes and longitudes are float32, up to 3·10⁻⁵° away from the regular 0.25° values. The shared regridding helper (_to_grid in feather/diag/precip_extremes.py) now gives a field that is on the target grid within float tolerance the exact target coordinates. Before, it kept the field's own coordinates, and xarray arithmetic with another field kept only the points with exactly equal coordinates (655 × 1340 of 721 × 1440 for HadGEM3). That caused a shape mismatch in the ssh_variability bias statistics and could silently shrink fields in the other diagnostics that use the helper (precip_extremes, temp_extremes_indices, both marine heatwave diagnostics).

CMORLoader now reads only the CMOR-named files ({var}_{table}_*.nc) of a directory when there are any: the IFS-NEMO-ER r1 HROday/zos directory also holds hand-merged zos_merged_*.nc copies (one broken) that broke or duplicated the read.

Climate classification

Diagnostic	Class	Observation	What it produces
climate_classification	KTClimateClassification	ERA5, Berkeley Earth HR + MSWEP	Köppen–Trewartha (KT14) land climate type per grid cell; discrete classification maps, % land-area bar chart, and summary table per source

The Köppen–Trewartha Climate Classification (climate_classification) classifies the monthly climatology of tas and pr into the 14 KT types (Ar, Aw, As, BW, BS, Cs, Cw, Cf, Do, Dc, Eo, Ec, Ft, Fi) using the Trewartha dryness threshold R = 2.3·Tann − 0.64·Pw + 41 (cm) and the F→E→B→D→C→A priority order. The classification is restricted to land via the Berkeley Earth land mask and computed on a common 0.25° grid so per-type area fractions are directly comparable. Each model is classified from its own tas+pr; observations come from ERA5 (t2m+tp) and a Berkeley Earth HR + MSWEP combination; the CMIP6 multi-model mean is included when enabled. Four figures are produced:

• Classification maps — one discrete KT panel per source, ordered ERA5, Berkeley Earth HR + MSWEP, the EERIE models, then the CMIP6 multi-model mean.
• Ensemble comparison — the observational datasets alongside the mean and median of each model family. Families are defined by a per-model ensemble label in the config (e.g. EERIE, or CORDEX/CMIP5/CMIP6 for the Africa config), and each ensemble is classified from the per-cell, per-month mean/median of its members' tas and pr climatologies (not by averaging KT codes).
• Bar chart and summary table — % land area per KT type per source, including every family's ensemble mean/median and the CMIP6 MMM (each source sums to 100 %).

Regional mode (Africa). Setting project.region (e.g. africa) restricts the classification to a regional bounding box, builds the common grid on that box (−180..180 longitude convention so western Africa is not dropped), handles rotated-pole (2-D lat/lon) sources, and draws zoomed PlateCarree maps. configs/africa_climate.yaml classifies, on a common 0.25° Africa grid for the 1981–2010 reference: ERA5/obs, CORDEX-CORE AFR-22 (10 GCM-driven members + 4 ERAINT evaluation runs; CORDEXLoader stitches historical 1981–2005 + rcp85 2006–2010 and falls back from mon to day-resampled-to-monthly), CMIP5 (5 driving GCMs, 1st member, historical+rcp85 stitch via CMIP5Loader), and CMIP6 (5 driving-GCM successors, historical, via the NetCDF-tree CMIP6NCLoader). Coarse GCMs use a 250 km regrid influence-radius floor.

Climate shifts at a warming level (climate_shifts, Phase 2). Compares the reference classification with the classification at a target global warming level (default +2 °C above the 1850–1900 pre-industrial baseline, configured under project.gwl). The 30-year window is found per GCM from its own global-mean temperature (feather/util/gwl.py; CMIP5→rcp85, CMIP6→ssp585), and CORDEX regional models inherit the window of their driving GCM (driving_gcm in the config). Per model family it produces: shift maps [reference | future | change], a reference-vs-future area bar, a net gain/loss summary (Δ % land area per type), and a transition matrix (reference type → warming-level type, area-weighted). Reference and future KT-code NetCDFs (per source and per family ensemble mean) and an area-shift CSV are written to {output_dir}/climate_shifts/. Run with --diagnostics climate_shifts on the same africa_climate.yaml.

Outputs written to {output_dir}/climate_classification/ for later regional analysis: a {source}_clim_{period}.nc per source with the monthly tas (°C) and pr (mm/month) climatology — both a full-global field and a land-only version, on the common grid — for ERA5, Berkeley Earth HR, MSWEP, every EERIE model, every CMIP6 model, the EERIE ensemble mean/median, and the CMIP6 MMM; a {source}_kt_{period}.nc with the integer kt_code field per classified source; and a kt_area_percent CSV. See feather/util/koeppen_trewartha.py for the pure, vectorised classifier.

Cross-domain

Diagnostic	Class	Observation	What it produces
precipitation_mswep	PrecipitationMSWEP	MSWEP v2.8	Precipitation bias maps (absolute + relative), time series, seasonal cycle, zonal mean, intensity PDF
temperature_berkeley	TemperatureBerkeley	Berkeley Earth	T2m bias maps, warming trend maps (global + polar), Taylor diagram
teleconnections	TeleconnectionDiag	ERA5	Climate variability modes (ENSO, NAO, SAM, AO, IOD, PDO, QBO): index time series, spatial patterns, power spectra, seasonal variance
obs_comparison	ObsComparisonDiag	ERA5 + Berkeley Earth	ERA5 vs Berkeley Earth T2m trend and bias comparison across two periods (1980–2014, 1980–2024)
precip_obs_comparison	PrecipObsComparisonDiag	ERA5 + MSWEP v2.8	ERA5 vs MSWEP precipitation trend and bias comparison across two periods (1980–2014, 1980–2023)

Model intercomparison

Diagnostic	Class	Observation	What it produces
added_value	AddedValueDiag	ERA5 / Berkeley Earth / MSWEP (atmos); ESA-CCI + HadISST / EN4 / OSI-SAF (ocean)	Dosio et al. (2015) Added Value: EERIE ensemble vs CMIP6/HighResMIP MMM — ensemble summary maps + per-model panels. Includes a dedicated Ocean Added Value page (tos vs ESA-CCI and HadISST, thetao/so vs EN4, siconc vs OSI-SAF) that reuses the ocean diagnostics' benchmark-bias NetCDFs

Added Value (AV) quantifies where the EERIE ensemble outperforms the CMIP6 multi-model mean relative to observations. AV ∈ [-1, 1]: AV > 0 means EERIE reduces squared error vs CMIP6 MMM at that grid point. Two figures per period (annual, DJF, JJA): ensemble mean/median summary and one panel per individual EERIE and CMIP6 model.

Per-region Added Value (--added-value-regions)

Off by default — only the global-domain figures are produced. Two region catalogues are available, and both can run in one invocation.

cordex14 (also what the bare flag selects, preserving the original behaviour) — the 14 CORDEX-CMIP6 *-11 domains from feather/util/regions.py, as summary bar charts per region for annual/DJF/JJA. Written to figures/added_value_regions/ as their own "CORDEX Regions" nav page.

ar6 — the 58 IPCC AR6 reference regions of Iturbide et al. (2020) (46 land + 15 ocean, three belonging to both sets), with polygons from regionmask. Written to figures/added_value_ar6/ as an "AR6 Regions" nav page:

• Tables — region × member heatmaps for every configured season. Three per season: all 58 regions, plus one per keep-set threshold. With two references each cell is split on the diagonal; a black box marks a positive value.
• Maps — annual gridded AV with region outlines, one panel per member per reference; plus ensemble mean/median with the keep-set outlined, one figure per threshold.
• CSV — {var}_ar6_av_per_member.csv, one row per (reference, period, member, region).

AV is computed per grid cell and then area-averaged over each region (cos-lat weighted, NaN-aware), rather than derived from region-mean biases — otherwise a warm half cancelling a cold half would score as skill.

A region is drawn and tabulated when at least N individual members have positive regional AV; the ensemble mean and median are excluded from that count, being summaries of the members rather than further evidence. Both thresholds come from added_value.region_thresholds in the config (default [6, 3]: a majority of a ten-member ensemble, and a looser "helps at least a subset" view). With more than one reference the keep-set is the union across them, so the map outlines always match the table columns.

References per variable are set by added_value.ar6_references, defaulting to tas: [ERA5] and pr: [ERA5, MSWEP]. Berkeley Earth (BE) is available but is not a tas default: it is land-only, so it cannot score the 15 ocean regions.

Unlike the rest of the diagnostic, the ar6 set reads the bias NetCDFs that the bias-map diagnostics already wrote ({output_dir}/netcdf/{global_biases,precipitation_mswep,temperature_berkeley}/) instead of recomputing, which takes the whole regional pass from hours to seconds. Generate them first; a missing file is reported and skipped rather than silently triggering a full recompute:

# 1) bias NetCDFs for the references you want (compute node)
feather --config configs/eerie_10_mems_cmip6.yaml \
        --diagnostics global_biases precipitation_mswep \
        --variables tas pr --benchmarks CMIP6 --save-netcdf \
        --steps diagnostics -v

# 2) AR6 maps + tables (cheap enough for a login node)
feather --config configs/eerie_10_mems_cmip6.yaml \
        --diagnostics added_value --variables tas pr \
        --added-value-regions ar6 --benchmarks CMIP6 \
        --steps diagnostics -v

All diagnostics support:

• variables=["tas", ...] — filter which variables to evaluate (CMOR canonical names)
• experiment="hist-1950" — experiment identifier (from config by default)
• period=("1980", "2014") — time period for climatologies (from config by default)
• cmip6_individual=True — individual CMIP6 model lines/biases alongside MMM

Configuration

Feather uses YAML configuration files. Sixteen configs are provided:

Config	Model set	Data source	Comparison type
configs/default.yaml	DestinE (3 models)	intake catalogs	multi_model
configs/eerie.yaml	EERIE Ensemble (4 models)	CMOR directory tree	multi_model
configs/eerie_ifsnemo_members.yaml	EERIE — 3 IFS-NEMO-ER + 3 other models	CMOR (hist-1950 + hist-1975)	multi_model
configs/eerie_all_members.yaml	EERIE — 8 models (3×FESOM2 + 3×NEMO + 2)	CMOR + kerchunk parquet	multi_model
configs/eerie_climchange_tn.yaml	EERIE — IFS-FESOM2-SR (r1–r3) + ICON-ESM-ER, SSP2-4.5 TN signal	CMOR + kerchunk native	multi_model
configs/eerie_climchange_hw.yaml	EERIE — IFS-FESOM2-SR (r1–r3) + ICON-ESM-ER, SSP2-4.5 heatwave signal	CMOR + kerchunk native	multi_model
configs/eerie_all_members_extremes.yaml	EERIE — ERA5 + the 10 members (3×FESOM2 + 3×NEMO + 3×ICON incl. r2/r3 via icon_kerchunk; HadGEM3 has no daily tasmax and is skipped), heatwave-hotspot tail-widening (PNAS Fig 2–4), per-family summaries	CMOR + kerchunk parquet (ERA5 derived daily + EERIE day table)	multi_model
configs/eerie_10_mems_precip_extremes.yaml	EERIE — 10 members (3×FESOM2 + 3×NEMO + 3×ICON + HadGEM3-GC5E-HH), ETCCDI daily precipitation extremes, 1981–2000 + SSP2-4.5 2031–2050, vs daily CMIP6, ERA5 and MSWEP	CMOR + kerchunk (EERIE day tables), DRS NetCDF (CMIP6), derived daily obs	multi_model
configs/eerie_all_members_indices.yaml	EERIE — ERA5 + 10 members (3×FESOM2 + 3×NEMO + 3×ICON + HadGEM3, which has no daily tasmax/tasmin), daily temperature extremes: tropical nights, heatwaves, temp_extremes_change and the 14 ETCCDI temp_extremes_indices, 1981–2000 + SSP2-4.5 2031–2050, vs ERA5 and daily CMIP6	CMOR + kerchunk (EERIE day tables), DRS NetCDF (CMIP6), derived daily ERA5	multi_model
configs/eerie_10_mems_marine_heatwaves.yaml	EERIE — 10 members, of which 7 publish daily SST (3×FESOM2, NEMO r1, 3×ICON), marine heatwaves (Hobday et al. 2016; de Toma et al. 2026) 1982–2014 vs OISST v2.1, ERA5 and daily CMIP6, per family	CMOR + kerchunk (EERIE Oday), DRS NetCDF (CMIP6 Oday), derived daily obs	multi_model
configs/himansu_319.yaml	IFS-FESOM T319	per-year NetCDF	single_model
configs/tco_grib.yaml	IFS-FESOM TCO399/TCO319	GRIB files	resolution_sensitivity
configs/destine_ifs_fesom.yaml	IFS-FESOM only	intake catalogs	single_model
configs/destine_ifs_nemo.yaml	IFS-NEMO only	intake catalogs	single_model
configs/destine_icon.yaml	ICON only	intake catalogs	single_model
configs/terradt.yaml	TerraDT (3 models)	intake catalogs	baseline_evaluation
configs/destine_added_value.yaml	DestinE (3 models), 1990–2025	intake catalogs (hist + ssp3-7.0 stitched)	multi_model
configs/ifs_fesom_combined.yaml	IFS-FESOM multi-res	mixed (catalog + GRIB)	resolution_sensitivity

Config format

The models key determines the config format:

• List (models: [ifs-fesom, ifs-nemo, icon]) — legacy DestinE format, auto-populates HEALPix grid defaults
• Dict (models: {ModelA: {institution: ..., grids: ...}}) — structured format with per-model configuration

data_source.type selects the loading backend: "cmor", "destine_catalog", "netcdf_healpix", "grib", or "kerchunk_parquet". Individual models can override the global backend with a per-model data_source_type key (used by eerie_all_members.yaml to mix CMOR and kerchunk parquet in the same run).

Structured config example (EERIE-style)

project:
  name: "EERIE"
  description: "EERIE Ensemble evaluation"
  experiment: "hist-1950"
  period: ["1980", "2014"]
  resolution: "high-resolution (~10-20 km atm, ~5-10 km ocean)"
  comparison_type: "multi_model"         # multi_model | resolution_sensitivity | single_model | baseline_evaluation
  comparison_description: ""             # optional free-text for LLM prompt framing

data_source:
  type: "cmor"                           # cmor | destine_catalog | netcdf_healpix | grib | kerchunk_parquet
  root: "/path/to/CMOR/tree"

models:
  IFS-FESOM2-SR:
    institution: AWI
    experiment: hist-1950
    variant: r1i1p1f1
    grids:
      sfc: latlon                        # "healpix" or "latlon" per domain
      o2d: latlon
      o3d: latlon
    color: "#1f77b4"

Multi-source config (CompositeModelLoader)

For evaluations mixing data backends, set per-model data_source_type:

data_source:
  type: "destine_catalog"                # default backend

models:
  TCO2999:
    data_source_type: "destine_catalog"  # per-model override (optional, uses default if omitted)
    catalog_key: "ifs-fesom"             # catalog entry name (if different from model name)
    member: 1                            # ensemble member index
    grids: {sfc: healpix, o2d: healpix, o3d: healpix}
  TCO399:
    data_source_type: "grib"             # different backend
    grids: {sfc: latlon, o2d: latlon, o3d: latlon}

The same pattern is used in eerie_all_members.yaml to mix CMOR (r1 members) with kerchunk parquet (r2/r3 members):

data_source:
  type: "cmor"
  root: "/path/to/CMOR/EERIE/HighResMIP"

models:
  IFS-FESOM2-SR:                         # r1: standard CMOR path
    institution: AWI
    variant: r1i1p1f1
    ...
  IFS-FESOM2-SR-r2:                      # r2: kerchunk parquet reference store
    institution: AWI
    variant: r2i1p1f1
    data_source_type: "kerchunk_parquet"
    data_root: "/path/to/kerchunks/IFS-FESOM2/hist-1950"
    ...

Per-model overrides

ModelConfig supports:

• data_root — custom data path (e.g., HadGEM3 outside shared CMOR root)
• grid_label — CMOR grid label override (default "gr")
• variable_aliases — CMOR var → on-disk name mapping (e.g., thetao: thetao-con)
• scale_factors — per-variable multiplier (e.g., clt: 100 for fraction → %)
• absolute_salinity: true — enable SA→SP salinity conversion (NEMO-based models)

Experiment stitching (extended time periods)

Some evaluation windows are not covered by a single experiment. For DestinE, the historical baseline_hist run ends in 2014 while the projections_ssp3-7.0 run continues to 2049 (IFS-FESOM / IFS-NEMO) or 2044 (ICON). To analyse e.g. 1990–2025, list the experiments to concatenate along the time axis:

project:
  period: ["1990", "2025"]
  # Stitched along time: historical baseline + ssp3-7.0 projection.
  experiments: ["baseline_hist", "projections_ssp3-7.0"]
  # Optional: let the timeseries diagnostic extend to each model's native end
  # (IFS-FESOM/IFS-NEMO → 2049, ICON → 2044) while every other diagnostic
  # uses `period` above. Defaults to `period` when omitted.
  timeseries_period: ["1990", "2050"]

cmip6:
  # Stitch CMIP6 symmetrically so added-value / bias diagnostics cover the
  # same window. Models lacking ssp370 for their variant fall back to
  # historical-only (end at 2014).
  experiments: ["historical", "ssp370"]

Segments are concatenated, sorted, and de-duplicated on overlapping months (the first experiment in the list wins). Models that did not run a listed experiment are skipped gracefully. Period slicing happens after concatenation, so a single stitched series serves both the analysis window and the extended timeseries. Omitting experiments preserves the legacy single-experiment behaviour, so existing configs are unaffected. See configs/destine_added_value.yaml for a complete example.

Note: an ssp370 variant label may differ from its historical counterpart (e.g. CanESM5 ssp370 is r1i1p2f1, not r1i1p1f1); the configured variant must exist for both experiments, or the model reverts to historical-only.

Comparison types

The comparison_type field controls how LLM prompts frame the analysis:

Type	Focus
multi_model	Inter-model differences and common biases across a multi-model ensemble
resolution_sensitivity	How model skill scales with resolution (same model at different grid spacings)
single_model	Detailed evaluation of one model against observations and CMIP6
baseline_evaluation	Whether a model is adequate as a starting point for a downstream project

Common config sections

All configs share:

• Observations: ERA5, CERES, EN4, ESA-CCI, OSI-SAF, MSWEP, Berkeley Earth
• CMIP6: 12 models from zarr archives (configurable ensemble mode)
• LLM analysis: Gemini (GEMINI_API_KEY / VERTEX_API_KEY)
• Report generation: OpenAI (OPENAI_API_KEY)
• nereus: Interpolation settings (method, resolution, influence_radius)

Benchmark ensembles (CMIP6 & HighResMIP)

Feather can overlay two reference benchmark ensembles as multi-model means (MMM): CMIP6 historical and HighResMIP hist-1950. Each is built once from the DKRZ pool into a per-variable zarr cache (scripts/convert_pool_cmip6.py, SLURM wrapper scripts/convert_pool_cmip6.sh) and configured with models: auto, so ensemble membership is discovered at runtime from the cache — see configs/eerie_all_members_cmip6_highresmip.yaml.

Daily cache (day table). The monthly converter deliberately skips the huge day table. A companion pair — scripts/convert_pool_cmip6_daily.py (SLURM wrapper scripts/convert_pool_cmip6_daily.sh) — builds a daily zarr cache of tas, tasmax, tasmin and pr (+ areacella) for the daily extremes / Tmin/Tmax-change diagnostics. It reuses the monthly converter's I/O helpers (atomic, resumable, calendar-aware) and writes stores named {model}_{experiment}_{member}_day_{var}.zarr into a dedicated cache dir (default /work/bk1580/cmip6_zarr); the experiment is encoded in the filename so historical (reference) and ssp245 (SSP2-4.5 future) coexist. Run sbatch scripts/convert_pool_cmip6_daily.sh for both experiments (or … historical / … ssp245), and preview with bash scripts/convert_pool_cmip6_daily.sh both --dry-run.

Selection rule. A model contributes to a benchmark MMM for a variable only when its converted store:

1. is non-empty (the conversion actually wrote the field);
2. is on a regriddable grid — rectilinear, curvilinear (2-D lat/lon), or unstructured with 1-D lat/lon cell-centre coords (e.g. ICON's triangular grid, regridded as scattered points like HEALPix); and
3. covers the full analysis period 1980–2014 (require_full_coverage, default on; opt out with require_full_coverage: false in the benchmark config block).

Membership grows automatically as more stores are converted. The lists below reflect the current cache for atmospheric tas; per-variable membership varies with which stores converted.

CMIP6 historical — 62 of 63 candidate models included:

ACCESS-CM2, ACCESS-ESM1-5, AWI-CM-1-1-MR, AWI-ESM-1-1-LR, AWI-ESM-1-REcoM, BCC-CSM2-MR, BCC-ESM1, CAMS-CSM1-0, CAS-ESM2-0, CESM2, CESM2-FV2, CESM2-WACCM, CESM2-WACCM-FV2, CIESM, CMCC-CM2-HR4, CMCC-CM2-SR5, CMCC-ESM2, CNRM-CM6-1, CNRM-CM6-1-HR, CNRM-ESM2-1, CanESM5, CanESM5-CanOE, E3SM-1-0, E3SM-1-1, E3SM-1-1-ECA, EC-Earth3, EC-Earth3-AerChem, EC-Earth3-CC, EC-Earth3-Veg, EC-Earth3-Veg-LR, FGOALS-f3-L, FGOALS-g3, FIO-ESM-2-0, GFDL-CM4, GFDL-ESM4, GISS-E2-1-G, GISS-E2-1-G-CC, GISS-E2-1-H, HadGEM3-GC31-LL, HadGEM3-GC31-MM, ICON-ESM-LR, IITM-ESM, INM-CM4-8, INM-CM5-0, IPSL-CM5A2-INCA, IPSL-CM6A-LR, KACE-1-0-G, KIOST-ESM, MCM-UA-1-0, MIROC-ES2L, MIROC6, MPI-ESM-1-2-HAM, MPI-ESM1-2-HR, MPI-ESM1-2-LR, MRI-ESM2-0, NESM3, NorCPM1, NorESM2-LM, NorESM2-MM, SAM0-UNICON, TaiESM1, UKESM1-0-LL

Excluded (1) — pending (re)conversion: IPSL-CM6A-LR-INCA — no Amon tas store in the cache (its other variables are present, so this is a per-variable gap rather than a missing model). It rejoins the ensemble automatically once the store is written.

The seven models previously listed here as excluded — AWI-ESM-1-REcoM, CAS-ESM2-0, CNRM-CM6-1-HR, HadGEM3-GC31-MM, IPSL-CM6A-LR, MPI-ESM-1-2-HAM and MPI-ESM1-2-LR — have since been converted and rejoined without any code change, which is the selection rule working as intended.

HighResMIP hist-1950 — 20 of 22 candidate models included:

BCC-CSM2-HR, CESM1-CAM5-SE-HR, CESM1-CAM5-SE-LR, CMCC-CM2-HR4, CMCC-CM2-VHR4, CNRM-CM6-1, CNRM-CM6-1-HR, EC-Earth3P, EC-Earth3P-HR, ECMWF-IFS-HR, ECMWF-IFS-LR, ECMWF-IFS-MR, GFDL-CM4C192, HadGEM3-GC31-HH, HadGEM3-GC31-HM, HadGEM3-GC31-LL, HadGEM3-GC31-MM, INM-CM5-H, MPI-ESM1-2-HR, MPI-ESM1-2-XR

Excluded (2):

• AWI-CM-1-1-HR, AWI-CM-1-1-LR — no usable tas store yet (empty / not yet converted); rejoin once rebuilt.

BCC-CSM2-HR and split version directories. BCC publishes its hist-1950 Amon record as two halves under different DRS versions — v20200822 holds tas_..._195001-200012.nc, v20200921 holds tas_..._200101-201412.nc. Reading only the latest version, the normal CMIP6 rule, therefore yielded 168 of 780 months and the model failed the full-coverage check. The converter now carries a version-union override for it (VERSION_UNION_OVERRIDES in scripts/convert_pool_cmip6.py), which unions non-overlapping segments across versions, and after reconversion tas, pr and ts span the full 1980–2014.

psl, uas and prw remain at 168 months: only the 2001–2014 file was ever published for them, so BCC-CSM2-HR counts towards the tas ensemble above but is absent from those three variables' MMMs. Omon/SImon were never affected — tos/siconc keep both time chunks under a single version.

Extend the behaviour to other publications with --union-versions MODEL[/EXPERIMENT[/TABLE[/VARIABLE]]], or drop the built-in entry with --no-builtin-union. Selection is strictly additive — an older file is taken only when no newer version publishes the same span — so a genuine supersede is never downgraded. A scan of every multi-version Amon/tas directory in the CMIP historical and HighResMIP pools found this to be the only such publication.

Re-converting only the missing stores: the converter's skip_existing treats an empty store as "done", so delete the empties first — the verification sweep writes their paths to …/cmip6_empty_stores.txt and …/highresmip_empty_stores.txt — or pass --no-skip.

Data backends

Feather supports seven data loading backends:

Backend	Class	Config type	Grid
DestinE intake catalogs	MultiCatalogLoader	destine_catalog	HEALPix (zarr)
CMOR directory tree	CMORLoader	cmor	regular lat/lon (NetCDF)
Per-year NetCDF	NetCDFLoader	netcdf_healpix	HEALPix (NetCDF)
GRIB files	GRIBLoader	grib	regular lat/lon (GRIB)
Kerchunk parquet reference stores	KerchunkParquetLoader	kerchunk_parquet	regular lat/lon (via fsspec)
ICON kerchunk stores	ICONKerchunkLoader	icon_kerchunk	regular 0.25° lat/lon (via fsspec)
Multi-source	CompositeModelLoader	mixed	per-model

KerchunkParquetLoader reads IFS-FESOM2 ensemble members stored as kerchunk parquet reference files pointing to raw GRIB/FESOM output. Two store layouts are supported:

• Regridded 0.25° (gr025): {root}/{variant}/atmos/gr025/*.parq — atmosphere 2D/3D monthly; ocean/gr025/2D_daily_avg_*.parq ocean 2D daily. Used by eerie_all_members.yaml.
• Native grid (native): {root}/{variant}/atmos/native/2D_daily_native_atmos_min.parq / …_max.parq — daily min/max on the native HEALPix grid. Selected via data_source_type: kerchunk_native in the model config; used by eerie_climchange_tn.yaml for the tropical_nights_change diagnostic.

Requires kerchunk, fastparquet, and fsspec in the Python environment.

ICONKerchunkLoader reads the ICON-ESM-ER hist-1950 members r2/r3, which the EERIE CMOR tree does not publish (only r1 is CMOR'd). They exist as intake catalogues over kerchunk reference stores; the loader reads those stores directly. Despite the name overlap it shares nothing with KerchunkParquetLoader — the ICON stores are already rectilinear (time, lat, lon) on the standard 0.25° grid with near-CMOR variable names. Used by eerie_10_mems_cmip6.yaml via data_source_type: icon_kerchunk plus member: 2 / member: 3. Note these members have no usable 3-D ocean or atmosphere, so ocean_en4 and the teleconnections QBO mode skip them.

The CompositeModelLoader automatically routes load_var() calls to the correct backend per model based on data_source_type in the model config.

Regridding

Diagnostics that build bias maps put every model, observation and benchmark onto a common lat/lon grid with nereus. Two settings govern how.

Interpolation method

nereus.method (default nearest) sets the scheme for state variables — temperature, pressure, wind, cloud cover, SST, sea ice. linear is recommended and used by the shipped configs.

Flux and precipitation fields are handled differently. Point-interpolation schemes do not preserve an area integral: a target cell takes the value of one source point (nearest) or a distance-weighted blend (linear) instead of the area-weighted mean of the source cells it covers. Coarsening precipitation that way biases the domain total and reports point values as box means, which overstates extremes. So pr, hfss, hfls and the radiation fields are remapped area-conservatively.

Coarsening a spiky precipitation-like field from 0.5° to 2°:

method	integral error	reported maximum
nearest	+0.649 %	50.1
idw	−0.487 %	18.3
linear	−0.011 %	37.9
conservative	+0.004 %	11.4

On a smooth field all four agree to ~0.01 %, so this only matters for intermittent fields — precipitation above all.

nereus:
  method: "linear"                  # state variables
  conservative_fluxes: true         # true | false | require (see below)
  conservative_max_points: 8000000  # cost guard, see below

conservative_fluxes takes three values:

value	behaviour
true (default)	use conservative for fluxes; warn and fall back if unavailable
false	point interpolation everywhere
require	fail rather than fall back

Prefer require for production runs. The capability cannot be read from nereus.__version__ — upstream did not bump it — so an environment built from a stale pin silently produces different numbers instead of an obvious failure, and the warning is easy to lose in a long log.

Classification is registry-driven (is_flux_variable() in feather/data/variables.py), keyed on variable group, so a flux variable added later is picked up automatically.

Cost guard

Conservative weights come from a spherical Voronoi tessellation plus polygon overlap, far costlier than a KD-tree. Measured onto a 0.25° target (1,036,800 cells), single core:

source points	weight build	apply
4,050	220 s	0.07 s
1,039,682 (a 0.25° model)	905 s (15 min)	0.01 s
6,480,000 (MSWEP 0.1°)	3,073 s (51 min), peak RSS 13.6 GB	0.18 s

Even a tiny source costs 220 s, because the fixed cost is set by the target tessellation. Applying cached weights is free, so the cost is one-off per source grid and amortises across every variable, season and period.

Pole-inclusive grids

A regular lat/lon grid spanning −90→90 collapses every longitude to a single point at each pole, which scipy.spatial.SphericalVoronoi rejects as duplicate generators. Since the EERIE, ICON and ERA5 grids are all 721×1440 spanning ±90, feather routes conservative remapping through feather/util/regrid.py, which merges coincident source points (averaging their values) before building weights. This is transparent — nothing in the configs or diagnostics changes — but it is a consumer-side workaround, and the durable fix belongs in nereus.

conservative_max_points budgets source + target together; over budget, feather warns and falls back to nereus.method. Per-dimension limits cannot distinguish 6.5M→1M (affordable) from 6.5M→6.5M (not).

Two practical consequences:

• The budget must clear the observation grid, not just the models. MSWEP is 0.1° = 6,480,000 points. At the 2,000,000 default, the MSWEP→0.25° coarsening — the step that most needs conserving — silently falls back while the near-identity 0.25°→0.25° model regrid gets conservative treatment, i.e. exactly backwards. eerie_10_mems_cmip6.yaml sets 8,000,000 for this reason.
• The 13.6 GB peak is in the main process, not a dask worker, so dask.memory_limit does not bound it. Request at least 32 GB for a job that builds conservative weights at this scale.

Still excluded at any sane budget: anything onto a 0.1° common grid (6.48M target cells on its own) and DestinE nside=1024 sources (12.6M points). DestinE precipitation therefore uses linear.

Precipitation common grid

precipitation_mswep historically built its common grid at the MSWEP native 0.1°, refining every model onto it regardless of the model's own resolution. nereus.precip_resolution now sets it (unset → obs native, the original behaviour):

nereus:
  precip_resolution: 0.25   # EERIE: models are 0.25 deg
  # precip_resolution: 0.1  # DestinE: runs are ~5 km, finer than the obs

At 0.25° the models are compared like with like and MSWEP is coarsened conservatively onto them, rather than the models being refined onto the obs grid. added_value evaluates pr against MSWEP too and honours the same key.

This is deliberately not nereus.resolution: every shipped config declares that as 0.25, so it cannot distinguish the two cases.

Adding a new model set

1. Create a config YAML using the structured format (see configs/eerie.yaml as template)
2. Set data_source.type to one of: "cmor", "destine_catalog", "netcdf_healpix", "grib", "kerchunk_parquet", "icon_kerchunk"
3. Define models as a dict with per-model institution, experiment, variant, grids, color
4. Set project.name, project.period, project.experiment, project.resolution
5. Optionally set project.comparison_type for LLM prompt framing
6. For mixed data sources, set per-model data_source_type and use CompositeModelLoader
7. Run: feather --config configs/my_project.yaml -v

If your data format is not supported, implement a new loader class (see CMORLoader or GRIBLoader as templates) and add a dispatch case in feather/run.py:_create_model_loader().

Output structure

output/
  figures/                          # Diagnostic figures + JSON metadata
    global_biases/
      tas_annual_bias_combined.png
      tas_annual_bias_combined.json
      ...
    timeseries/
    seasonal_cycle/
    radiation_budget/
    sea_ice/
    ocean_sst/
    ocean_en4/
    global_trends/
    climate_variability/
    precipitation_mswep/
    temperature_berkeley/
    teleconnections/
    obs_comparison/
    precip_obs_comparison/
    added_value/
    tropical_nights/
    tropical_nights_change/
    heatwave/
    heatwave_change/
    heatwave_hotspots/
    temp_extremes_change/
  tropical_nights/                    # NC checkpoints (outside figures tree)
    {model}_tropical_nights_tn20_{start}_{end}.nc
  tropical_nights_change/             # NC checkpoints (outside figures tree)
    {model}_tn_hist_{start}_{end}.nc
    {model}_tn_ssp_{start}_{end}.nc
  heatwave/                           # NC checkpoints (outside figures tree)
    {model}_heatwave_tx90_{start}_{end}.nc
  heatwave_change/                    # NC checkpoints (outside figures tree)
    {model}_hw_hist_{start}_{end}.nc
    {model}_hw_ssp_{start}_{end}.nc
  heatwave_hotspots/                  # NC checkpoints (outside figures tree)
    {model}_percs_{start}_{end}.nc
  temp_extremes_change/               # NC checkpoints (outside figures tree)
    {model}_tasmin_hist_{start}_{end}.nc
    {model}_tasmin_ssp_{start}_{end}.nc
    {model}_tasmax_hist_{start}_{end}.nc
    {model}_tasmax_ssp_{start}_{end}.nc
  analysis/                         # LLM analysis (Gemini)
    global_biases/
      tas_annual_bias_combined_analysis.json
      synthesis.json
    ...
  publication/                      # LaTeX report (OpenAI)
    structure.json                  # Stage 1: editorial curation
    sections/                       # Stage 2: per-section prose
    figures/                        # Copied figures for LaTeX
    report.tex                      # Final LaTeX document
    report.pdf                      # (if --compile-pdf)
  site/                             # Static HTML dashboard
    index.html
    global_biases.html
    ...

Project structure

feather/
  cli.py                 # CLI entry point (feather command)
  config.py              # FeatherConfig + ModelConfig dataclasses
  run.py                 # Pipeline orchestration (run_pipeline)
  __main__.py            # python -m feather support
  data/
    loader.py            # MultiCatalogLoader (DestinE intake catalogs)
    cmor_loader.py       # CMORLoader (CMOR directory tree, e.g. EERIE)
    netcdf_loader.py     # NetCDFLoader (per-year NetCDF on HEALPix grid)
    grib_loader.py       # GRIBLoader (GRIB files on regular lat/lon grid)
    kerchunk_loader.py   # KerchunkParquetLoader (kerchunk parquet reference stores)
    composite_loader.py  # CompositeModelLoader (multi-source routing)
    obs.py               # ObsLoader (observations from config)
    cmip6.py             # CMIP6Loader (multi-model mean from zarr)
    variables.py         # VARIABLE_REGISTRY (36 vars, CMOR canonical names)
  util/
    spatial.py           # Zonal/global means, regridding, latlon areas
    temporal.py          # Climatology, anomaly, deseason, detrend
    eof.py               # EOF computation via SVD (teleconnections)
    spectrum.py          # Power spectrum via Welch (teleconnections)
    units.py             # Unit conversions
  plot/
    maps.py              # Bias maps (nereus), combined multi-panel maps
    lines.py             # Time series, seasonal cycle, zonal profiles, Taylor diagram
    styles.py            # Model colors, plot defaults
  diag/
    base.py              # DiagnosticBase ABC + grid-agnostic helpers
    registry.py          # @register decorator, diagnostic discovery
    figure_meta.py       # Figure metadata sidecar system
    global_biases.py     # Climatology bias maps (18 variables)
    timeseries.py        # Global-mean time series
    seasonal_cycle.py    # Monthly climatological cycle
    radiation_budget.py  # Radiation budget analysis (CERES EBAF)
    sea_ice.py           # Sea ice evaluation (OSI-SAF, PIOMAS/GIOMAS)
    ocean_sst.py         # SST evaluation (ESA-CCI)
    ocean_en4.py         # 3D ocean T/S (EN4 v4.2.2)
    global_trends.py     # Per-grid-point linear trends
    climate_variability.py # STD of deseasonalised, detrended fields
    precipitation_mswep.py # Precipitation evaluation (MSWEP v2.8)
    temperature_berkeley.py # T2m evaluation (Berkeley Earth)
    teleconnections.py    # Variability modes (ENSO, NAO, SAM, AO, IOD, PDO, QBO)
    obs_comparison.py     # ERA5 vs Berkeley Earth T2m trends and biases
    precip_obs_comparison.py # ERA5 vs MSWEP precipitation trends and biases
    added_value.py        # Added Value: EERIE ensemble vs CMIP6 MMM (Dosio et al. 2015)
    tropical_nights.py    # Tropical Nights Index (TN > 20 °C, daily tasmin)
    tropical_nights_change.py # Tropical Nights climate change signal (SSP2-4.5)
    heatwave.py           # Heatwave Indices (TX90: HWN, HWF, HWD, HWM, HWA)
    heatwave_change.py    # Heatwave climate change signal (SSP2-4.5)
    heatwave_hotspots.py  # Heatwave hotspots — extreme-heat tail-widening (PNAS Fig 2–4)
    temp_extremes_change.py # Daily Tmin/Tmax climate change signal (SSP2-4.5, annual + seasons)
  llm/
    analyzer.py          # FigureAnalyzer (Gemini, comparison-type aware)
    schemas.py           # FigureAnalysis, DiagnosticSynthesis (Pydantic)
    prompts.py           # Gemini prompts (composable, comparison-type aware)
  export/
    report.py            # ReportGenerator (OpenAI, 3-stage)
    schemas.py           # ReportStructure, WrittenSection
    prompts.py           # OpenAI prompts (composable, comparison-type aware)
    openai_client.py     # Thin OpenAI wrapper with retry
    latex_builder.py     # LaTeX escaping, template rendering, PDF
    templates/           # Jinja2 LaTeX template
  website/
    generator.py         # SiteGenerator (static HTML dashboard)
    templates/           # Jinja2 HTML templates
    static/              # Dark-theme CSS

Adding a new diagnostic

Use CMOR canonical variable names and grid-agnostic base class helpers:

from feather.diag.base import DiagnosticBase
from feather.diag.registry import register

@register
class MyDiagnostic(DiagnosticBase):
    name = "my_diagnostic"
    title = "My Diagnostic"
    domain = "sfc"
    variables = ["tas"]         # CMOR canonical names
    group = "temperature"

    def compute(self):
        for model in self.config.models:
            mdata = self._load_model_var(model, "tas")      # Grid-agnostic
            ts = self._model_global_mean(mdata.data, model)  # HEALPix or latlon
            # ...
        obs = self._load_obs_var("tas")                      # Sign-convention aware
        return {"model_ts": ..., "obs_ts": ...}

    def plot(self, results):
        fig, ax = plt.subplots()
        meta = self._build_metadata(title="...", figure_id="...", models=[...])
        return [(fig, meta)]

Then import the module in feather/diag/__init__.py so @register fires. See NEW_DIAGNOSTIC_SPEC.md for the full implementation guide.

Testing

conda activate feather

# Unit tests (synthetic data, no real data needed)
pytest tests/ -v -m "not integration"

# Integration tests (requires data access on DKRZ Levante)
pytest tests/ -v -m "integration"

# All tests
pytest tests/ -v

2483 tests across 63 test files (2468 unit + 15 integration).

Requirements

Python 3.10+. See environment.yml for the full dependency list, or install via conda env create -f environment.yml.

Key dependencies: xarray, dask, numpy, scipy, matplotlib, cartopy, healpy, nereus, gsw, cfgrib, pydantic, google-genai, openai.