Forecasting record-breaking weather extremes at least as well as operational NWP

A Science Advances 2026 study (arXiv 2508.15724) evaluates forecasts on record-breaking heat, cold and wind events. It finds that ECMWF HRES consistently outperforms GraphCast, Pangu-Weather and FuXi there, and that the AI models underestimate both the frequency and the magnitude of record exceedance. Average-skill leaderboards (WeatherBench 2) hide this because record events are a tiny fraction of the scored points.

Source: science.org

import Mathlib.Data.Real.Basic
import Mathlib.Algebra.Order.BigOperators.Group.Finset
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.NormNum
import Mathlib.Tactic.Positivity
import Mathlib.Data.Fintype.Prod

/-!
# Forecasting record-breaking weather extremes at least as well as operational NWP

On mean WeatherBench-2 skill, AI forecasts already beat ECMWF IFS. What stays
open is the tail. Evaluated on events where the observed value breaks the prior
record, IFS HRES still outperforms GraphCast, Pangu-Weather and FuXi, and the AI
models underestimate both how often records break and by how much (Science
Advances 2026, arXiv 2508.15724). This file states "solved" on that event set
only.

Concrete here:
* The **record-event set**, `recordEvents`. It holds every (initialisation, cell)
  pair whose verifying truth exceeds that cell's prior record. It is computed
  from truth and records, so neither the submitter nor the host picks the events.
* The **latitude-weighted squared error on that set**, and "beats the baseline
  by margin δ in RMSE". The comparison is stated in squared form, over the same
  events, so the normalisation cancels.
* **Bounded underestimation.** The weighted mean signed error on record events
  must be at least `-biasTol`. This is the failure mode the literature reports.
* **Detection without crying wolf.** The submission must catch strictly more
  record events than the baseline, and raise no more false alarms than the
  baseline. False alarms are forecasts above record where none occurred.
  Without this clause, "always forecast extreme" could score well.

Orientation: records here are maxima, as for heat and wind. For cold records, a
board supplies negated temperature fields; nothing else changes.

Opaque oracle (`Eval`), for one variable at one lead time:
* verification truth (for example ERA5 or station observations over the held-out period);
* the IFS HRES forecasts from the archive;
* latitude weights;
* per-cell prior records, from climatology that ends before the evaluation period starts;
* attestations of no look-ahead, and that model weights and records predate the
  evaluation period.

Who must attest: the benchmark host publishes truth, records, weights and the
baseline archive version. An independent party attests that the submissions are
genuine initial-time forecasts.

Units: fields and `biasTol` are in the variable's own units (K for T2M, m/s for
10 m wind). Weights and `margin` are dimensionless; 0.05 means 5 % lower RMSE.

Checked with Lean 4.33.0 / Mathlib using:
  lake env lean <absolute-path-to-this-file>
-/

namespace OpenSpur.Problems.WeatherRecordExtremes

/-- A field indexed by initialisation time and grid cell, for one variable and lead. -/
abbrev Field (n g : ℕ) := Fin n → Fin g → ℝ

/-- Sealed evaluation data. -/
structure Eval (n g : ℕ) where
  /-- Latitude-cosine grid weights. -/
  weight : Fin g → ℝ
  truth : Field n g
  /-- Operational baseline forecast (IFS HRES) at the same lead. -/
  baseline : Field n g
  /-- Per-cell record maximum from climatology ending before the evaluation period. -/
  record : Fin g → ℝ
  noLookAhead : Prop
  weightsPredateEval : Prop
  recordsPredateEval : Prop

variable {n g : ℕ}

/-- Record-breaking events: the truth exceeds the prior record. -/
noncomputable def recordEvents (E : Eval n g) : Finset (Fin n × Fin g) :=
  Finset.univ.filter (fun x => E.record x.2 < E.truth x.1 x.2)

/-- Pairs where a forecast `f` exceeds the prior record. -/
noncomputable def forecastRecords (E : Eval n g) (f : Field n g) : Finset (Fin n × Fin g) :=
  Finset.univ.filter (fun x => E.record x.2 < f x.1 x.2)

/-- Record events the forecast also puts above record. -/
noncomputable def hits (E : Eval n g) (f : Field n g) : ℕ :=
  (recordEvents E ∩ forecastRecords E f).card

/-- Forecast records where no record occurred. -/
noncomputable def falseAlarms (E : Eval n g) (f : Field n g) : ℕ :=
  (forecastRecords E f \ recordEvents E).card

/-- Weighted squared error of `f` on the record-event set. -/
noncomputable def eventSSE (E : Eval n g) (f : Field n g) : ℝ :=
  ∑ x ∈ recordEvents E, E.weight x.2 * (f x.1 x.2 - E.truth x.1 x.2) ^ 2

/-- Weighted signed error of `f` on the record-event set (negative = underforecast). -/
noncomputable def eventBias (E : Eval n g) (f : Field n g) : ℝ :=
  ∑ x ∈ recordEvents E, E.weight x.2 * (f x.1 x.2 - E.truth x.1 x.2)

/-- Total weight of the record-event set. -/
noncomputable def eventWeight (E : Eval n g) : ℝ :=
  ∑ x ∈ recordEvents E, E.weight x.2

structure Thresholds where
  /-- RMSE improvement over the baseline on record events; must lie in `[0, 1]`. -/
  margin : ℝ
  /-- Largest tolerated weighted mean underforecast on record events. -/
  biasTol : ℝ
  /-- Minimum number of record events, so the test is not decided by a handful. -/
  minEvents : ℕ

/-- Acceptance predicate for a submitted forecast `f`. -/
def Solved (T : Thresholds) (E : Eval n g) (f : Field n g) : Prop :=
  E.noLookAhead ∧ E.weightsPredateEval ∧ E.recordsPredateEval ∧
  (∀ p, 0 < E.weight p) ∧
  0 ≤ T.margin ∧ T.margin ≤ 1 ∧
  T.minEvents ≤ (recordEvents E).card ∧
  eventSSE E f ≤ (1 - T.margin) ^ 2 * eventSSE E E.baseline ∧
  -T.biasTol * eventWeight E ≤ eventBias E f ∧
  hits E E.baseline < hits E f ∧
  falseAlarms E f ≤ falseAlarms E E.baseline

/-- Rung `k`: `5 (k+1)` % lower RMSE on record events, underforecast at most 0.5
units on average, and at least 100 record events. -/
noncomputable def Rung (k : ℕ) : Thresholds where
  margin := 0.05 * ((k : ℝ) + 1)
  biasTol := 0.5
  minEvents := 100

theorem rung_margin_mono {j k : ℕ} (h : j ≤ k) : (Rung j).margin ≤ (Rung k).margin := by
  have hjk : (j : ℝ) ≤ (k : ℝ) := Nat.cast_le.mpr h
  simp only [Rung]
  linarith

theorem eventSSE_nonneg (E : Eval n g) (hw : ∀ p, 0 < E.weight p) (f : Field n g) :
    0 ≤ eventSSE E f :=
  Finset.sum_nonneg fun x _ => mul_nonneg (hw x.2).le (sq_nonneg _)

/-- Solving a higher rung implies solving every lower rung. -/
theorem solved_rung_mono {j k : ℕ} (h : j ≤ k) (E : Eval n g) (f : Field n g) :
    Solved (Rung k) E f → Solved (Rung j) E f := by
  rintro ⟨h1, h2, h3, hw, _, hk1, hev, hsse, hbias, hhit, hfa⟩
  have hm := rung_margin_mono h
  have hj0 : 0 ≤ (Rung j).margin := by
    simp only [Rung]; positivity
  have hsq : (1 - (Rung k).margin) ^ 2 ≤ (1 - (Rung j).margin) ^ 2 := by
    have : 0 ≤ 1 - (Rung k).margin := by linarith
    nlinarith
  refine ⟨h1, h2, h3, hw, hj0, by linarith, hev, ?_, hbias, hhit, hfa⟩
  exact le_trans hsse (mul_le_mul_of_nonneg_right hsq (eventSSE_nonneg E hw _))

end OpenSpur.Problems.WeatherRecordExtremes

• WeatherBench 2 (benchmark + leaderboard vs IFS HRES/ENS)
• ECMWF AI Weather Quest (live sub-seasonal ML competition, 2025-26)
• Record-breaking extremes evaluation (Science Advances / arXiv 2508.15724)
• GenCast code, weights and forecasts (open release)
• NVIDIA Earth2Studio / Earth-2 open models
• Silurian NEXRAD storm catalog (HF dataset)
• S2S AI Challenge (WMO, 2021; benchmark vs recalibrated ECMWF)