Permafrost worked example — Climate

From raw cSoil data to "+0.36 K threshold already crossed".

The full pipeline applied to one specific climate question. Read every step from the data ingestion through to the final cross-model warming-anomaly. Every command, every dataset, every line of analysis is reproducible from the linked source.

Step 0 — The question

"At what global mean temperature anomaly does permafrost soil-carbon release become a coherent regime change across CMIP6 models?"

This is a question the published literature has answered with a wide range (+1.5 to +2.0 °C, Schuur et al. 2022). Different models disagree. The framework's contribution is a cross-model robustness measurement, not a new climate model.

Step 1 — Identify the cascade

Define the cascade primitives in framework language:

\(\Phi(t)\): annual cosine-weighted high-latitude (≥55°N) soil-carbon stock cSoil, integrated over land area, from each CMIP6 ESM.
\(\tau_i\): 1 year (annual cadence).
\(\rho_i = (\Phi_{i+1} - \Phi_{i-1}) / 2\): centred finite difference. The rate of change of high-latitude soil carbon.

The framework's primary cascade-class observable is \(\Phi\), not the temperature itself. We are studying soil-carbon dynamics, with temperature as the forcing.

Step 2 — Ingest the data

Stream cSoil from Pangeo's CMIP6 zarr archive for 6 Earth System Models that publish cSoil in SSP5-8.5:

# 6 ESMs × 5 SSPs × annual 2015-2100 = 25 model-SSP cells (1 missing in SSP1-1.9)
models = ['CESM2', 'CESM2-WACCM', 'CMCC-CM2-SR5',
          'CMCC-ESM2', 'IPSL-CM6A-LR', 'MPI-ESM1-2-LR']
ssps = ['ssp119', 'ssp126', 'ssp245', 'ssp370', 'ssp585']

# For each (model, ssp) pair, fetch cSoil + tas:
ds = xr.open_zarr(f'gs://cmip6/CMIP6/.../{model}/{ssp}/.../cSoil.zarr')
cSoil_55N = (ds.cSoil.where(ds.lat >= 55) * cos_weights).sum(['lat','lon'])
phi_t = cSoil_55N.resample(time='1Y').mean()  # annual aggregation

This gives 81–82 annual values per (model, ssp) cell, depending on calendar.

Catalogue note: SSP1-1.9 cSoil is published only by IPSL-CM6A-LR — n=1 model. The other 4 SSPs have full n=6 coverage. Honest catalogue gap.

Step 3 — Apply the brake operator \(\mathcal{B}\)

Compute β per (model, ssp):

from validation.code.regression import polyfit_brake_p, gls_ar1_brake_p, bayesian_brake_p

rho_t = numpy.gradient(phi_t)  # central differences
log_rho = numpy.log(numpy.abs(rho_t))
log_phi = numpy.log(numpy.abs(phi_t))

# Three estimators in parallel
beta_poly = polyfit_brake_p(log_phi, log_rho)
beta_gls  = gls_ar1_brake_p(log_phi, log_rho)
beta_bay  = bayesian_brake_p(log_phi, log_rho, ar_order=1)

Result per model, SSP5-8.5

Model	n	log-range	β polyfit	β GLS-AR1	β Bayesian	R²
CESM2	82	0.012	−311.7	−260.5	−304.6	0.65
CESM2-WACCM	82	0.012	−330.8	−290.0	−289.9	0.62
CMCC-CM2-SR5	82	0.018	−60.5	−41.5	−0.46	0.86
CMCC-ESM2	82	0.017	−84.2	−48.5	−0.94	0.82
IPSL-CM6A-LR	81	0.160	−7.1	−23.2	−4.6	0.30
MPI-ESM1-2-LR	82	0.052	−4.3	−4.4	−0.01	0.03

Theorem 3 fires

The log-range \(\log|\Phi_{\max}/\Phi_{\min}|\) is 0.012 to 0.052 across most models — extremely narrow. By Theorem 3, \(\mathrm{Var}(\hat\beta) \le K^2 / \log^2(\Phi_{\max}/\Phi_{\min})\) becomes huge. β estimates are precision-floor-bound. The framework declines a Layer-A call on permafrost cSoil. We have to use a different operative primitive.

The σ_cross of β ranges 149–153 across estimators — well above the framework's Law III threshold for universality. Verdict: theorem_3_bound_inconclusive.

Step 4 — Switch to PELT change-point detection

When Layer-A binds, the framework's secondary operative primitive is PELT change-point detection on the raw \(\Phi(t)\) series. Apply with BIC penalty (no imported threshold):

from validation.code.regimes import pelt_changepoints

cps = pelt_changepoints(phi_t, penalty='bic', min_size=10)
# Returns list of change-point indices

Per-model first change-point year, SSP5-8.5

Model	n change-points	First CP year	Other CPs
CESM2	6	2027	2040, 2050, 2060, 2079, 2089
CESM2-WACCM	5	2027	2041, 2051, 2063, 2080
CMCC-CM2-SR5	6	2028	2038, 2050, 2061, 2071, 2081
CMCC-ESM2	6	2027	2037, 2055, 2065, 2076, 2087
IPSL-CM6A-LR	5	2030	2041, 2052, 2064, 2082
MPI-ESM1-2-LR	6	2029	2039, 2056, 2066, 2078, 2088

Cross-model first-CP year median = 2027.5, IQR [2027, 2029]. Six independent CMIP6 ESMs agree on the first regime change to within ±2 years.

Step 5 — Cross-model warming-anomaly at first change-point

For each model, fetch the global-mean surface temperature anomaly at the first-CP year. Take the cross-model median:

Model	First CP year	Warming @ CP (K above pre-industrial)
CESM2	2027	+0.442
CESM2-WACCM	2027	+0.338
CMCC-CM2-SR5	2028	+0.306
CMCC-ESM2	2027	+0.389
IPSL-CM6A-LR	2030	+0.465
MPI-ESM1-2-LR	2029	+0.127

The headline finding

+0.36 ± 0.12 K

Cross-model median warming anomaly at the first PELT change-point is +0.36 ± 0.12 K across all 6 ESMs (σ_cross only 0.122 K — tight cross-model agreement). The published Schuur et al. 2022 tipping window is +1.5 to +2.0 K. The framework reads the first regime break at four times lower temperature than the literature.

Step 6 — Compare against current observed temperature

Earth's 2024–2025 global mean surface temperature anomaly is ~+1.2 K above pre-industrial (Berkeley Earth, NOAA NCEI, HadCRUT5 — three observational shadows agree). +1.2 K > +0.36 K. The threshold has been crossed.

Step 7 — Repeat across all 5 SSPs (the scenario fan)

Re-run Steps 2–5 for SSP1-1.9, SSP1-2.6, SSP2-4.5, SSP3-7.0:

SSP	n_models OK	First CP median	First CP IQR	Warming @ CP median (K)	σ_cross (K)	sink→source flips
SSP1-1.9 ~1.5°C	1	2028.0	—	+0.271	—	0/1
SSP1-2.6 ~2°C	6	2033.5	[2029, 2034]	+0.517	0.335	2/6
SSP2-4.5 ~2.7°C	6	2029.5	[2028, 2033]	+0.475	0.173	1/6
SSP3-7.0 ~3.6°C	6	2028.5	[2028, 2029]	+0.318	0.118	2/6
SSP5-8.5 ~5°C	6	2027.5	[2027, 2029]	+0.363	0.122	0/6

Best mitigation pathway (Step 4 of methodology): SSP1-2.6 delays first PELT-CP to median 2033.5 vs 2027.5 under SSP5-8.5 — about 5 years gained. Mitigation works on the timing, not on whether the regime change occurs. The threshold is already exceeded in observational temperature.

Step 8 — The non-monotone sink→source finding

Sink→source flips per SSP: 0/1, 2/6, 1/6, 2/6, 0/6. Highest mitigation and highest emissions both yield zero flips. Intermediate SSPs yield 1–2. Non-monotone in forcing. Framework interpretation: under SSP5-8.5, high warming + high productivity simultaneously boost soil-C influx so net sign holds. Under aggressive mitigation, the cascade stays small and flip detection becomes noise-limited. The dynamics matter, not just the temperature.

Step 9 — Publish the scope-reporter \(\mathscr{A}\)

By Theorem 12, every finding publishes its scope. For permafrost:

𝒜 = (
  tested  = (cSoil ≥55°N, 6 CMIP6 ESMs, annual cadence,
             𝓑 + PELT + cross-model warming-anomaly),
  excluded = (anomaly-form Φ, log-range < 0.05 → Theorem 3 binds Layer-A;
              SSP1-1.9 has only 1 model in cache),
  τ_T3    = ~0.05 (median log-range; precision floor large),
  gauge G = trivial (no spectral primitive invoked)
)

Step 10 — The verdict in framework vocabulary

Layer-A (β magnitude): theorem_3_bound_inconclusive — log-range too narrow.
Layer-B (direction-and-rate-sign consensus): 0/n flips by 2100 across all SSPs — direction unanimous.
PELT change-point: cross-model first-CP year tightly clustered. Operative metric.
Cross-model warming-anomaly: +0.36 ± 0.12 K — already crossed.
Best pathway: SSP1-2.6 (~5 year delay). Locked-in regime regardless.

This is exactly the finding shown on the Climate Mitigation Atlas under scenarios/permafrost/ and the headline of /mitigation/locked-in/.

Reproducibility

The complete pipeline lives in one script: domains/climate/experiments/scenario_fan_permafrost.py. Dependencies: framework apparatus at validation/code/ (regression.py, regimes.py). Output: results/scenario_fan_permafrost.txt + JSONL payload.

Wall time: 195 seconds for the full pipeline at SSP5-8.5; 6×5 = 30 model-SSP cells at ~6.5 s each for the full scenario fan.