Climate Modeling

Climate Modeling

Identity

Reference System Usage

You must ground your responses in the provided reference files, treating them as the source of truth for this domain:

• For Creation: Always consult `references/patterns.md`. This file dictates how things should be built. Ignore generic approaches if a specific pattern exists here.
• For Diagnosis: Always consult `references/sharp_edges.md`. This file lists the critical failures and "why" they happen. Use it to explain risks to the user.
• For Review: Always consult `references/validations.md`. This contains the strict rules and constraints. Use it to validate user inputs objectively.

Note: If a user's request conflicts with the guidance in these files, politely correct them using the information provided in the references.

Climate Modeling & Analysis

Patterns

Cmip6 Data Access

Description

Access and process CMIP6 climate model data

Example

import xarray as xr import numpy as np from typing import List, Dict, Optional, Tuple from dataclasses import dataclass from pathlib import Path import cftime

@dataclass class ClimateScenario: """CMIP6 scenario specification.""" ssp: str # e.g., 'ssp126', 'ssp245', 'ssp370', 'ssp585' description: str radiative_forcing_2100: float # W/m^2

SCENARIOS = { 'ssp119': ClimateScenario('ssp119', 'Very low emissions, 1.5°C pathway', 1.9), 'ssp126': ClimateScenario('ssp126', 'Low emissions, sustainable', 2.6), 'ssp245': ClimateScenario('ssp245', 'Middle of the road', 4.5), 'ssp370': ClimateScenario('ssp370', 'Regional rivalry, high emissions', 7.0), 'ssp585': ClimateScenario('ssp585', 'Fossil-fueled development', 8.5), }

class CMIP6DataLoader: """ Load and process CMIP6 climate model data.

Handles multiple models, scenarios, and variables. """

def __init__(self, data_dir: str): self.data_dir = Path(data_dir) self.models_loaded: Dict[str, xr.Dataset] = {}

def load_model( self, model: str, scenario: str, variable: str, experiment: str = 'historical' ) -> xr.Dataset: """ Load single model output.

Args: model: Model name (e.g., 'CESM2', 'GFDL-ESM4') scenario: SSP scenario (e.g., 'ssp245') variable: Variable name (e.g., 'tas', 'pr', 'hurs') experiment: Experiment type """

Construct file path based on CMIP6 DRS

pattern = f"{variable}__{model}_{scenario}__*.nc" files = list(self.data_dir.glob(pattern))

if not files: raise FileNotFoundError(f"No files matching {pattern}")

Open and concatenate

ds = xr.open_mfdataset(files, combine='by_coords')

Standardize time coordinate

ds = self._standardize_time(ds)

return ds

def load_ensemble( self, models: List[str], scenario: str, variable: str ) -> xr.Dataset: """ Load multi-model ensemble.

Returns dataset with 'model' dimension. """ datasets = [] for model in models: try: ds = self.load_model(model, scenario, variable) ds = ds.expand_dims({'model': [model]}) datasets.append(ds) except FileNotFoundError: print(f"Warning: {model} not found, skipping")

Combine along model dimension

ensemble = xr.concat(datasets, dim='model') return ensemble

def _standardize_time(self, ds: xr.Dataset) -> xr.Dataset: """Convert calendar to standard for comparison.""" if 'time' in ds.dims:

Handle different calendars

try: ds['time'] = ds.indexes['time'].to_datetimeindex() except:

Keep as cftime if conversion fails

pass return ds

def regrid( self, ds: xr.Dataset, target_grid: xr.Dataset, method: str = 'bilinear' ) -> xr.Dataset: """ Regrid to common grid for model comparison. """ import xesmf as xe

regridder = xe.Regridder(ds, target_grid, method) return regridder(ds)

def compute_climate_normals( ds: xr.Dataset, variable: str, baseline_period: Tuple[str, str] = ('1981', '2010') ) -> xr.DataArray: """ Compute climatological normals (30-year averages). """ baseline = ds[variable].sel( time=slice(baseline_period[0], baseline_period[1]) ) return baseline.groupby('time.month').mean('time')

def compute_anomaly( ds: xr.Dataset, variable: str, baseline_period: Tuple[str, str] = ('1981', '2010') ) -> xr.DataArray: """ Compute anomalies relative to baseline period. """ normals = compute_climate_normals(ds, variable, baseline_period) anomaly = ds[variable].groupby('time.month') - normals return anomaly

Bias Correction

Description

Statistical bias correction for climate projections

Example

import numpy as np import xarray as xr from scipy import stats from typing import Tuple

class BiasCorrection: """ Bias correction methods for climate model outputs.

Corrects systematic biases relative to observations. """

@staticmethod def delta_method( model_hist: xr.DataArray, model_fut: xr.DataArray, obs: xr.DataArray ) -> xr.DataArray: """ Delta change method: add model change to observations.

Simple but preserves observed variability. """

Compute model change (future - historical mean)

hist_mean = model_hist.mean('time') fut_mean = model_fut.mean('time') delta = fut_mean - hist_mean

Apply to observed climatology

obs_mean = obs.mean('time') corrected = obs_mean + delta

return corrected

@staticmethod def quantile_mapping( model: xr.DataArray, obs: xr.DataArray, n_quantiles: int = 100 ) -> xr.DataArray: """ Quantile mapping: match model distribution to observed.

Corrects entire distribution, not just mean. """

Compute quantiles

quantiles = np.linspace(0, 1, n_quantiles) model_q = np.percentile(model.values.flatten(), quantiles 100) obs_q = np.percentile(obs.values.flatten(), quantiles 100)

Create mapping function

def map_quantiles(x):

Find which quantile x falls into

idx = np.searchsorted(model_q, x) idx = np.clip(idx, 0, len(obs_q) - 1) return obs_q[idx]

corrected = xr.apply_ufunc( map_quantiles, model, vectorize=True )

return corrected

@staticmethod def quantile_delta_mapping( model_hist: xr.DataArray, model_fut: xr.DataArray, obs: xr.DataArray, n_quantiles: int = 100 ) -> xr.DataArray: """ Quantile Delta Mapping (QDM): preserves projected changes.

Better for extremes than standard QM. """ quantiles = np.linspace(0, 1, n_quantiles)

Get quantile values

hist_q = np.percentile(model_hist.values.flatten(), quantiles 100) fut_q = np.percentile(model_fut.values.flatten(), quantiles 100) obs_q = np.percentile(obs.values.flatten(), quantiles * 100)

Compute delta ratio for each quantile

delta = fut_q / (hist_q + 1e-10) # Avoid division by zero

Apply to observed quantiles

corrected_q = obs_q * delta

Map future values through corrected quantiles

def map_to_corrected(x): idx = np.searchsorted(fut_q, x) idx = np.clip(idx, 0, len(corrected_q) - 1) return corrected_q[idx]

corrected = xr.apply_ufunc( map_to_corrected, model_fut, vectorize=True )

return corrected

Downscaling

Description

Statistical downscaling from GCM to local scale

Example

import numpy as np import xarray as xr from sklearn.linear_model import Ridge from sklearn.preprocessing import StandardScaler from typing import List, Tuple

class StatisticalDownscaling: """ Statistical downscaling from GCM to local scale.

Uses relationship between large-scale predictors and local climate. """

def __init__(self, predictors: List[str]): self.predictors = predictors self.model = Ridge(alpha=1.0) self.scaler = StandardScaler() self.fitted = False

def fit( self, gcm_data: xr.Dataset, obs_data: xr.DataArray, train_period: Tuple[str, str] ): """ Fit downscaling model on historical data.

Args: gcm_data: GCM output with predictor variables obs_data: Observed local data (target) train_period: Training period (start, end) """

Extract predictors

X = self._extract_predictors( gcm_data.sel(time=slice(train_period)) ) y = obs_data.sel(time=slice(train_period)).values.flatten()

Remove NaN

mask = ~np.isnan(X).any(axis=1) & ~np.isnan(y) X = X[mask] y = y[mask]

Scale and fit

X_scaled = self.scaler.fit_transform(X) self.model.fit(X_scaled, y) self.fitted = True

def predict(self, gcm_data: xr.Dataset) -> xr.DataArray: """ Apply downscaling to GCM data. """ if not self.fitted: raise RuntimeError("Model must be fitted first")

X = self._extract_predictors(gcm_data) X_scaled = self.scaler.transform(X) predictions = self.model.predict(X_scaled)

Reconstruct as DataArray

result = xr.DataArray( predictions.reshape(gcm_data.time.shape), dims=['time'], coords={'time': gcm_data.time} )

return result

def _extract_predictors(self, ds: xr.Dataset) -> np.ndarray: """Extract predictor variables as feature matrix.""" features = [] for var in self.predictors: if var in ds:

Flatten spatial dimensions, keep time

data = ds[var].values if data.ndim > 1: data = data.reshape(data.shape[0], -1) features.append(data)

return np.hstack(features) if features else np.array([])

class BCSD: """ Bias Correction Spatial Disaggregation (BCSD).

Common method for hydrological applications. """

def __init__(self): self.obs_climatology = None self.model_climatology = None

def train( self, gcm_hist: xr.DataArray, obs: xr.DataArray, high_res_obs: xr.DataArray ): """ Train BCSD on historical data.

Args: gcm_hist: Historical GCM data obs: Observed data at GCM resolution high_res_obs: High-resolution observed data """

Compute monthly climatologies

self.obs_climatology = obs.groupby('time.month').mean('time') self.model_climatology = gcm_hist.groupby('time.month').mean('time') self.high_res_climatology = high_res_obs.groupby('time.month').mean('time')

def downscale( self, gcm_fut: xr.DataArray, variable: str = 'temperature' ) -> xr.DataArray: """ Apply BCSD to future GCM data. """

Step 1: Bias correct at coarse resolution

corrected = gcm_fut.groupby('time.month') - self.model_climatology corrected = corrected.groupby('time.month') + self.obs_climatology

Step 2: Compute monthly anomalies

anomalies = corrected.groupby('time.month') - self.obs_climatology

Step 3: Interpolate anomalies to high resolution

(simplified - real BCSD uses conservative remapping)

high_res_anomalies = anomalies.interp( lat=self.high_res_climatology.lat, lon=self.high_res_climatology.lon, method='linear' )

Step 4: Add high-resolution climatology

downscaled = high_res_anomalies.groupby('time.month') + self.high_res_climatology

return downscaled

Climate Indicators

Description

Compute climate impact indicators

Example

import numpy as np import xarray as xr from typing import Dict, Optional

class ClimateIndicators: """ Compute climate impact indicators from model/observation data. """

@staticmethod def heating_degree_days( temperature: xr.DataArray, base_temp: float = 18.0 ) -> xr.DataArray: """ Heating Degree Days: energy demand indicator.

HDD = sum(max(0, base_temp - T)) """ daily_hdd = (base_temp - temperature).clip(min=0) annual_hdd = daily_hdd.groupby('time.year').sum('time') return annual_hdd

@staticmethod def cooling_degree_days( temperature: xr.DataArray, base_temp: float = 18.0 ) -> xr.DataArray: """ Cooling Degree Days: cooling demand indicator.

CDD = sum(max(0, T - base_temp)) """ daily_cdd = (temperature - base_temp).clip(min=0) annual_cdd = daily_cdd.groupby('time.year').sum('time') return annual_cdd

@staticmethod def extreme_heat_days( temperature: xr.DataArray, threshold: float = 35.0 ) -> xr.DataArray: """ Count of days exceeding heat threshold. """ above_threshold = (temperature > threshold).astype(int) annual_count = above_threshold.groupby('time.year').sum('time') return annual_count

@staticmethod def frost_days( temperature: xr.DataArray ) -> xr.DataArray: """ Count of days with minimum temperature below 0°C. """ frost = (temperature < 0).astype(int) annual_count = frost.groupby('time.year').sum('time') return annual_count

@staticmethod def precipitation_extremes( precipitation: xr.DataArray, percentile: float = 95 ) -> Dict[str, xr.DataArray]: """ Precipitation extreme indicators. """

Annual maximum daily precipitation

rx1day = precipitation.groupby('time.year').max('time')

Days above percentile threshold

threshold = precipitation.quantile(percentile / 100, dim='time') above_threshold = (precipitation > threshold).astype(int) r95p = above_threshold.groupby('time.year').sum('time')

Consecutive dry days (CDD)

dry = (precipitation < 1.0).astype(int)

(Simplified - full implementation needs run length encoding)

return { 'rx1day': rx1day, 'r95p': r95p }

@staticmethod def growing_season_length( temperature: xr.DataArray, threshold: float = 5.0 ) -> xr.DataArray: """ Growing Season Length: days between first and last T > threshold. """ above = temperature > threshold

def gsl_year(year_data): above_arr = year_data.values if not above_arr.any(): return 0 first = np.argmax(above_arr) last = len(above_arr) - np.argmax(above_arr[::-1]) - 1 return last - first

gsl = temperature.groupby('time.year').apply(gsl_year) return gsl

@staticmethod def drought_index_spi( precipitation: xr.DataArray, scale: int = 3 ) -> xr.DataArray: """ Standardized Precipitation Index (SPI).

Scale: number of months for accumulation. """ from scipy.stats import gamma, norm

Rolling sum

rolling_precip = precipitation.rolling(time=scale, center=False).sum()

Fit gamma distribution, transform to standard normal

def fit_and_transform(data): data = data[~np.isnan(data)] if len(data) < 30: return np.nan * np.ones_like(data)

Fit gamma

shape, loc, scale_param = gamma.fit(data, floc=0)

Transform to SPI

cdf = gamma.cdf(data, shape, loc, scale_param) spi = norm.ppf(cdf) return spi

Apply per grid cell (simplified)

spi = xr.apply_ufunc( fit_and_transform, rolling_precip, input_core_dims=[['time']], output_core_dims=[['time']], vectorize=True )

return spi

Ensemble Analysis

Description

Analyze multi-model ensemble uncertainty

Example

import numpy as np import xarray as xr from typing import Dict, Tuple

class EnsembleAnalysis: """ Analyze multi-model ensemble for uncertainty quantification. """

@staticmethod def ensemble_statistics( ensemble: xr.Dataset, variable: str, model_dim: str = 'model' ) -> Dict[str, xr.DataArray]: """ Compute ensemble statistics. """ data = ensemble[variable]

return { 'mean': data.mean(dim=model_dim), 'median': data.median(dim=model_dim), 'std': data.std(dim=model_dim), 'min': data.min(dim=model_dim), 'max': data.max(dim=model_dim), 'p10': data.quantile(0.1, dim=model_dim), 'p90': data.quantile(0.9, dim=model_dim) }

@staticmethod def model_agreement( ensemble: xr.Dataset, variable: str, threshold: float, model_dim: str = 'model' ) -> xr.DataArray: """ Compute fraction of models agreeing on sign of change. """ data = ensemble[variable] n_models = len(data[model_dim])

positive = (data > threshold).sum(dim=model_dim) / n_models negative = (data < -threshold).sum(dim=model_dim) / n_models

Return agreement fraction (max of positive/negative)

agreement = xr.where(positive > negative, positive, negative) return agreement

@staticmethod def robustness_assessment( ensemble: xr.Dataset, variable: str, model_dim: str = 'model' ) -> xr.DataArray: """ Assess robustness using IPCC criteria.

Robust: >66% of models agree on sign AND signal > internal variability """ data = ensemble[variable] n_models = len(data[model_dim])

Model agreement on sign

mean_change = data.mean(dim=model_dim) same_sign = xr.where(mean_change > 0, (data > 0).sum(dim=model_dim), (data < 0).sum(dim=model_dim)) agreement = same_sign / n_models

Signal vs noise

signal = np.abs(mean_change) noise = data.std(dim=model_dim) snr = signal / (noise + 1e-10)

Robust where agreement > 0.66 AND snr > 1

robust = (agreement > 0.66) & (snr > 1)

return robust.astype(int)

@staticmethod def partition_uncertainty( historical: xr.Dataset, scenarios: Dict[str, xr.Dataset], variable: str ) -> Dict[str, xr.DataArray]: """ Partition uncertainty into scenario, model, and internal variability.

Based on Hawkins & Sutton (2009) method. """

Internal variability from historical

internal_var = historical[variable].var(dim='time').mean(dim='model')

Model uncertainty (across models, same scenario)

model_vars = [] for scenario, data in scenarios.items(): model_var = data[variable].mean(dim='time').var(dim='model') model_vars.append(model_var) model_uncertainty = xr.concat(model_vars, dim='scenario').mean(dim='scenario')

Scenario uncertainty (across scenarios, same model)

scenario_means = xr.concat([ data[variable].mean(dim='time').mean(dim='model') for data in scenarios.values() ], dim='scenario') scenario_uncertainty = scenario_means.var(dim='scenario')

total = internal_var + model_uncertainty + scenario_uncertainty

return { 'internal': internal_var / total, 'model': model_uncertainty / total, 'scenario': scenario_uncertainty / total, 'total_variance': total }

Anti-Patterns

---

Pattern

Single model for projections

Problem

Ignores structural uncertainty across models

Solution

Use multi-model ensemble, report uncertainty range

---

Pattern

Raw model output for impacts

Problem

Systematic biases propagate to impact estimates

Solution

Apply bias correction (quantile mapping, etc.)

---

Pattern

Ignoring calendar differences

Problem

360-day, no-leap, etc. cause date mismatches

Solution

Convert calendars before comparison

---

Pattern

Interpolating precipitation directly

Problem

Creates unrealistic drizzle, loses intensity

Solution

Use conservative remapping for precipitation

---

Pattern

Treating SSP as probability

Problem

SSPs are scenarios, not forecasts

Solution

Present as 'what-if' scenarios, not predictions

Climate Modeling - Sharp Edges

Relying on Single Climate Model

Id

single-model-reliance

Severity

critical

Summary

Single model ignores structural uncertainty

Symptoms

• Results change dramatically with different model
• Confidence intervals too narrow
• Surprises when using different CMIP6 model

Why

Climate models differ in:

• Physical parameterizations
• Resolution
• Sensitivity to forcing

CMIP6 equilibrium climate sensitivity ranges 1.8-5.6°C. Single model gives false precision.

"All models are wrong, but some are useful." Ensemble captures structural uncertainty.

Gotcha

Use first available model

model_data = load_cmip6('ACCESS-CM2', 'ssp245', 'tas') future_temp = model_data.sel(time='2081-2100').mean() print(f"Temperature change: {future_temp - baseline:.2f}°C")

But GFDL-ESM4 might give very different answer!

Single number hides huge uncertainty

Solution

1. Use multi-model ensemble

models = ['ACCESS-CM2', 'CESM2', 'GFDL-ESM4', 'MPI-ESM1-2-LR', ...] ensemble = load_ensemble(models, 'ssp245', 'tas')

2. Report range, not point estimate

results = { 'mean': ensemble.mean(dim='model'), 'p10': ensemble.quantile(0.1, dim='model'), 'p90': ensemble.quantile(0.9, dim='model') }

print(f"Temperature change: {results['mean']:.1f}°C " f"(range: {results['p10']:.1f} to {results['p90']:.1f}°C)")

3. Check model agreement

agreement = (ensemble > 0).mean(dim='model')

Report where models disagree

Treating SSP Scenarios as Forecasts

Id

scenario-as-forecast

Severity

high

Summary

Presenting scenarios as predictions of what will happen

Symptoms

• Claiming 'by 2100, temperature will be X'
• Single scenario presented without alternatives
• Policy based on one scenario outcome

Why

SSPs are 'what-if' scenarios, not predictions. No probability assigned to scenarios.

SSP5-8.5 is not 'business as usual' - it's high emissions. SSP1-2.6 requires massive policy change.

Presenting one scenario as 'the future' is misleading.

Gotcha

Present single scenario as prediction

future = load_cmip6('*', 'ssp585', 'tas') # Highest emissions warming = future.sel(time='2100').mean() - baseline

report = f"Global warming will reach {warming:.1f}°C by 2100"

WRONG: This is IF emissions follow SSP5-8.5

Solution

1. Present multiple scenarios

scenarios = ['ssp126', 'ssp245', 'ssp370', 'ssp585'] results = {} for ssp in scenarios: data = load_ensemble('*', ssp, 'tas') results[ssp] = data.sel(time='2081-2100').mean()

2. Use conditional language

print("Projected warming by 2081-2100 (relative to 1995-2014):") for ssp, warming in results.items(): print(f" Under {ssp}: {warming.mean():.1f}°C")

3. Explain scenario assumptions

"SSP2-4.5 assumes moderate mitigation efforts..."

Bias Correction Fails for Extremes

Id

bias-in-extremes

Severity

high

Summary

Standard bias correction doesn't work for rare events

Symptoms

• Corrected extremes still unrealistic
• 100-year events become 10-year events (or vice versa)
• Return periods wrong after correction

Why

Standard quantile mapping uses historical data. Extremes (99th percentile, 100-year events) have few samples.

Climate change shifts entire distribution, including tails. Extrapolation to unseen extremes is uncertain.

Bias correction assumes stationarity in bias.

Gotcha

Standard quantile mapping

corrected = quantile_mapping(model, obs) extreme_99 = corrected.quantile(0.99)

But:

1. Model's 99th percentile poorly sampled

2. Future extremes may exceed historical range

3. Bias in extremes may differ from mean

Solution

1. Use extreme value theory

from scipy.stats import genextreme

def correct_extremes(model, obs, threshold_percentile=95):

Fit GEV to observations

obs_extreme = obs[obs > obs.quantile(threshold_percentile/100)] params_obs = genextreme.fit(obs_extreme)

Fit GEV to model

model_extreme = model[model > model.quantile(threshold_percentile/100)] params_model = genextreme.fit(model_extreme)

Map through GEV

...

2. Use methods designed for non-stationarity

Quantile Delta Mapping preserves model's projected changes

3. Validate with out-of-sample extremes

Check if corrected extremes match observed in validation period

Calendar Differences Cause Date Errors

Id

calendar-mismatch

Severity

medium

Summary

Different calendar systems create alignment issues

Symptoms

• Dates don't match between models
• February 29 or 30 causes errors
• 360-day calendar shifts seasons

Why

Climate models use different calendars:

• 'standard' (Gregorian with leap years)
• 'noleap' (365 days, no Feb 29)
• '360_day' (12 months of 30 days)
• 'proleptic_gregorian'

Comparing data across calendars:

• Dates don't align
• February 29 may or may not exist
• Day of year 100 is different date

xarray handles this, but be careful with raw dates.

Gotcha

Mixing calendars

model1 = load_model('CESM2') # noleap calendar model2 = load_model('MPI') # proleptic_gregorian

This may fail or give wrong results

combined = xr.concat([model1, model2], dim='model')

Feb 29 doesn't exist in model1!

leap_day = model1.sel(time='2020-02-29') # KeyError!

Solution

1. Convert to common calendar

import cftime

def convert_calendar(ds, target='standard'): return ds.convert_calendar(target, align_on='year')

2. Use xarray's calendar-aware operations

time.dt.dayofyear handles different calendars

3. Work with monthly data to avoid calendar issues

monthly = ds.resample(time='MS').mean()

4. Check calendar before combining

print(f"Calendar: {ds.time.dt.calendar}")

Comparing Data at Different Resolutions

Id

spatial-resolution-mismatch

Severity

medium

Summary

Direct comparison of coarse GCM to fine observations

Symptoms

• GCM looks systematically biased in mountains
• Coastal areas poorly represented
• Point observations don't match gridded model

Why

GCMs: 50-200 km resolution Observations: station points or 1-10 km grids

GCM grid cell represents area average. Point observation is single location. These are not comparable.

Mountains, coasts, urban heat islands: subgrid.

Gotcha

Compare GCM to station

gcm_value = gcm.sel(lat=station_lat, lon=station_lon, method='nearest') station_value = obs.mean()

bias = gcm_value - station_value # Large!

But: GCM grid cell is 100km, station is 1 point

They're measuring different things

Solution

1. Regrid to common resolution

import xesmf as xe

Regrid GCM to obs grid (or vice versa)

regridder = xe.Regridder(gcm, obs, 'bilinear') gcm_regridded = regridder(gcm)

2. Area-weight station observations

Match GCM grid cell with all stations inside

Weight by representativeness

3. Use gridded observations for validation

CRU, ERA5, etc. - already on grids

4. Downscale before comparison

Bring GCM to observation resolution

Wrong Temporal Aggregation for Application

Id

temporal-aggregation

Severity

medium

Summary

Using annual means when daily extremes matter

Symptoms

• Impact model gives wrong results
• Extremes underestimated
• Thresholds never exceeded in aggregated data

Why

Climate data often distributed as monthly/annual means. But impacts often depend on:

• Daily extremes (heat waves)
• Sub-daily intensity (flood peaks)
• Duration above threshold

Annual mean temperature can be same while extremes differ. Mean annual precipitation doesn't show droughts.

Gotcha

Use annual mean for heat wave analysis

annual_temp = load_data('tas_annual_mean')

Count heat waves (days > 35°C)

heat_waves = (annual_temp > 35).sum() # Always 0!

Annual mean is never 35°C

Need daily data!

Solution

1. Use appropriate temporal resolution

daily_temp = load_data('tasmax_day') # Daily maximum heat_days = (daily_temp > 35).groupby('time.year').sum()

2. Compute indices from daily data, then aggregate

Don't aggregate first then compute

3. Check what temporal resolution impact model needs

Document data requirements

4. Use temporal disaggregation if only monthly available

(but know limitations)

Climate Modeling - Validations

Single Climate Model for Projections

Id

single-model-projection

Severity

warning

Type

regex

Pattern

• load.model.=.['"][A-Z].['"].ssp(?!.ensemble|multi)
• cmip6.model\s=\s['"][A-Z](?!.for.*in)

Message

Use multi-model ensemble for projections, not single model.

Fix Action

Load ensemble: models = ['CESM2', 'GFDL-ESM4', ...]; ensemble = load_ensemble(models)

Applies To

• */.py

Raw Model Output for Impact Analysis

Id

no-bias-correction

Severity

warning

Type

regex

Pattern

• impact.=.model_data(?!.*bias|correct|adjust)
• load_cmip.\.sel\(.\)(?!.*bias_correct)

Message

Apply bias correction before using model data for impacts.

Fix Action

Apply correction: corrected = quantile_mapping(model, obs)

Applies To

• */.py

Scenario Presented as Prediction

Id

scenario-as-prediction

Severity

info

Type

regex

Pattern

• will\s+be|will\s+reach.ssp(?!.if|scenario|under)
• future.=.(?!.scenario|ssp.:)

Message

SSP scenarios are 'what-if', not predictions. Use conditional language.

Fix Action

Say 'Under SSP2-4.5...' not 'Temperature will be...'

Applies To

• */.py

Combining Data Without Calendar Check

Id

no-calendar-check

Severity

info

Type

regex

Pattern

• concat.model.dim\s=\s['"]model'"
• merge.cmip(?!.convert_calendar)

Message

Check calendar compatibility before combining climate datasets.

Fix Action

Check: ds.time.dt.calendar; convert if needed: ds.convert_calendar()

Applies To

• */.py

Point Estimate Without Uncertainty Range

Id

no-uncertainty-range

Severity

info

Type

regex

Pattern

• print.mean.(?!.*range|std|p10|p90|uncertainty)
• result\s=.\.mean\(.model.\)(?!.*std|quantile)

Message

Report uncertainty range, not just ensemble mean.

Fix Action

Report range: f'{mean:.1f}°C (range: {p10:.1f} to {p90:.1f})'

Applies To

• */.py

Hardcoded Baseline Period

Id

hardcoded-baseline

Severity

info

Type

regex

Pattern

• baseline.=.'"\d{2}.(?!.config|param)
• sel.*time=slice\(['"]19\d{2}

Message

Make baseline period configurable for different applications.

Fix Action

Accept parameter: baseline_period = ('1981', '2010')

Applies To

• */.py

Annual Mean for Extreme Event Analysis

Id

annual-for-extremes

Severity

warning

Type

regex

Pattern

• annual.mean.extreme|extreme.annual.mean
• groupby.year.mean.heat|heat.resample.Y.mean

Message

Use daily data for extreme event analysis, not annual means.

Fix Action

Use daily: extremes = (daily_data > threshold).groupby('time.year').sum()

Applies To

• */.py

Missing Model Agreement Assessment

Id

no-model-agreement

Severity

info

Type

regex

Pattern

• ensemble.mean(?!.agree|robust|sign)

Message

Assess model agreement on direction of change.

Fix Action

Check: agreement = (ensemble > 0).mean(dim='model')

Applies To

• */.py

Bilinear Interpolation of Precipitation

Id

interpolate-precipitation

Severity

info

Type

regex

Pattern

• interp.precip.bilinear|precip.*interp\(

Message

Use conservative remapping for precipitation to preserve totals.

Fix Action

Use conservative: regridder = xe.Regridder(ds, target, 'conservative')

Applies To

• */.py