Advancing subseasonal forecasting with machine learning – Nature

Probabilistic subseasonal forecasting

We begin with a formal description of the probabilistic subseasonal forecasting problem, following the conventions of the ECMWF 2025 AI Weather Quest competition¹⁸. For each week-long period beginning on a date t, let \({{\bf{y}}}_{t}\in {{\mathbb{R}}}^{G}\) represent a vector of ground-truth weather measurements y_t,g for each grid cell g in a global latitude–longitude grid. We used 1979–2026 ERA5 global reanalysis³⁰ regridded to 1.5° × 1.5° latitude–longitude resolution as our measure of ground truth and focus on predicting average 2-m temperature, total precipitation and average mean sea-level pressure across each week-long period. Both temperature and precipitation are evaluated only over land-dominated regions (that is, grid cells with at least 50% land coverage).

For each grid cell g, our aim is to predict a distribution over weather at time t, expressed as a CDF F_t,g over K = 5 historical climatological quantile bins. To this end, we set q_t,g(K) = ∞ and, for each k < K, define the kth climatological quintile q_t,g(k) for (t,g) as the \(\frac{k}{K}\)th quantile of the observations y_s,g falling 0 days, 2 days or 4 days away from the target month–day combination in each of the previous 20 years. For each k, the probabilistic forecast F_t,g(k) represents the predicted probability that the measurement y_t,g will fall below the kth quintile threshold q_t,g(k).

Given a test set \({{\mathcal{T}}}^{\star }\) of target dates, we measure the global quality of these probabilistic forecasts by means of the RPSS²⁸,

where \({O}_{t,g}(k)={\mathbb{I}}[{y}_{t,g} < {q}_{t,g}(k)]\) is the observed indicator of whether y_t,g falls below the kth quintile threshold, \({w}_{g}=\cos ({\rm{latitude}}(g))\) adjusts for the relative size of each grid cell, RPS_t,g is termed the ranked probability score (RPS) of the forecast and \({\overline{{\rm{RPS}}}}_{t,g}\) is the RPS of a baseline climatological forecast that assigns equal probability to each quintile bin. For precipitation, arid grid cells g with no climatological precipitation (that is, with q_t,g(K − 1) = 0) are excluded from this evaluation.

Similarly, for a given grid cell g, we measure forecast quality over a test set \({{\mathcal{T}}}^{\star }\) by means of the spatial RPSS

Throughout, we evaluate RPSS at two subseasonal forecasting horizons, days 19–25 (week 3) and days 26–32 (week 4). Unless otherwise noted, we use all non-leap-day Monday and Friday target dates in the years 2016–2024 as our default test set.

Baselines

Extended Data Table 1 provides a summary of the dynamical, AI, hybrid (dynamical + AI) and post-processing baselines compared with PBC in this work. Here we detail the processing steps taken to evaluate each baseline.

Raw ECMWF

At each subseasonal lead time ℓ, the operational dynamical model of ECMWF¹⁷ outputs an ensemble of deterministic forecasts, \({({\widehat{y}}_{m,{\ell },t,g})}_{m=1}^{M}\), each representing a possible future weather state. For each target weather variable and subseasonal horizon, we downloaded all weekly-aggregated ECMWF ensemble forecasts, spanning the years 2015–2026, from the subseasonal to seasonal (S2S) prediction project database⁵⁴ through the International Research Institute for Climate and Society Data Library (IRIDL)⁵¹. Forecasts with M = 51 members (1 control and 50 perturbed) were issued every Monday and Thursday from 14 May 2015 to 26 June 2023 and forecasts with M = 101 members were issued daily starting on 28 June 2023. We converted each ensemble into a CDF forecast \({F}_{{\ell },t,g}(k)=\frac{1}{M}{\sum }_{m=1}^{M}{\mathbb{I}}[{\widehat{y}}_{m,{\ell },t,g} < {q}_{t,g}(k)]\) using the climatological quintile thresholds q_t,g.

Debiased ECMWF

To form the debiased ECMWF baseline, we followed the operational post-processing protocol of ECMWF³³ for each weather variable and subseasonal horizon. First, we downloaded all weekly-aggregated ECMWF ensemble hindcasts from the S2S database⁵² through the IRIDL⁵¹. Each ensemble hindcast consists of M = 11 members, \({({\widehat{y}}_{m,{\ell },\delta ,t,g})}_{m=1}^{M}\), indexed by a forecast target date t and the number of years δ ∈ {1, …, 20} between the forecast and corresponding hindcast target dates. Next, for each forecast target date t, we identified the set \({{\mathcal{S}}}_{t}\) of model climatology dates: the closest date no later than t with available hindcasts along with the next four and previous four available dates. We then formed model-based quintile thresholds for each forecast target date t, lead time ℓ and grid cell g by setting \({q}_{{\ell },t,g}^{{\rm{ECMWF}}}(K)=\infty \) and, for each k < K, defining the kth model quintile \({q}_{{\ell },t,g}^{{\rm{ecmwf}}}(k)\) as the \(\frac{k}{K}\)th quantile of all hindcasts \({\widehat{y}}_{m,{\ell },\delta ,s,g}\) with \(s\in {{\mathcal{S}}}_{t}\). Finally, we constructed the debiased ECMWF CDF forecast \({F}_{{\ell },t,g}(k)=\frac{1}{M}{\sum }_{m=1}^{M}{\mathbb{I}}[{\hat{y}}_{m,{\ell },t,g} < {q}_{{\ell },t,g}^{{\rm{ECMWF}}}(k)]\).

FuXi-S2S

For each target weather variable and subseasonal horizon, we downloaded 51-member ensemble forecasts, spanning the test years 2017–2021, and 51-member ensemble hindcasts, spanning the training years 2002–2016, from the official FuXi-S2S forecast dataset⁵³. Forecasts were subsequently aggregated to a weekly resolution, with precipitation accumulated and temperature and mean sea-level pressure averaged. Following ref. ¹³, we form model-based quintile thresholds for each forecast target date t, lead time ℓ and grid cell g by setting \({q}_{{\ell },t,g}^{{\rm{FuXi}}\,}(K)=\infty \) and, for each k < K, defining the kth model quintile \({q}_{{\ell },t,g}^{{\rm{FuXi}}\,}(k)\) as the \(\frac{k}{K}\)th quantile of all hindcasts matching the target month–day combination. Finally, we constructed the CDF forecast \({F}_{{\ell },t,g}(k)=\frac{1}{M}{\sum }_{m=1}^{M}{\mathbb{I}}[{\hat{y}}_{m,{\ell },t,g} < {q}_{{\ell },t,g}^{{\rm{FuXi}}\,}(k)]\) from each FuXi-S2S ensemble forecast \({({\widehat{y}}_{m,{\ell },t,g})}_{m=1}^{M}\) with M = 51.

Because FuXi-S2S and ECMWF historically issued their forecasts on different schedules, we evaluated all available forecasts from FuXi-S2S and PBC-ECMWF with target dates in the period 2017–2021. This corresponds to all Monday and Friday target dates for PBC-ECMWF and, for FuXi-S2S, all Tuesday and Friday target dates in 2017, all Wednesday and Saturday target dates in 2018, all Thursday and Sunday target dates in 2019, all Monday and Friday target dates from 3 January 2020 to 24 February 2020, all Tuesday and Saturday target dates from 29 February 2020 to 29 December 2020 and all Wednesday and Sunday target dates from 3 January 2021 to 26 December 2021.

AIFS-SUBS

For each target weather variable, subseasonal horizon and odd issuance date in 2025 (except for the 31st of any month), we generated weekly-aggregated 100-member ensemble forecasts using the AIFS-SUBS model¹⁶. AIFS-SUBS follows the architecture introduced in ref. ⁵⁴, combining a transformer-based graph neural network encoder and decoder with a sliding-window transformer processor in which attention is computed along spiral longitudinal bands. Ensemble members are generated by conditioning the processor on independently sampled noise, injected by means of conditional layer normalization in each processor layer and the model is trained with the fair continuous RPS loss computed over four ensemble members against ERA5 targets⁵⁵. AIFS-SUBS takes two consecutive states (t − 24 h, t) as input and predicts the atmospheric state at t + 24 h. Inputs and outputs are represented on the reduced Gaussian o96 grid (approximately 1° resolution). The variant of AIFS-SUBS used in this work includes a representation of the stratosphere and was trained exclusively on ERA5 over the periods 1979–2006 and 2012–2024, leaving a 5-year validation period.

Following ref. ¹⁶, for each forecast target date t, lead time ℓ and grid cell g, we generated 10-member ensemble hindcasts for the corresponding month–day combination in each of the previous 20 years, set \({q}_{{\ell },t,g}^{{\rm{AIFS}}\,}(K)=\infty \) and, for each k < K, defined the kth model quintile \({q}_{{\ell },t,g}^{{\rm{AIFS}}\,}(k)\) as the \(\frac{k}{K}\)th quantile of the aforementioned hindcasts. Finally, we constructed the CDF forecast \({F}_{{\ell },t,g}(k)=\frac{1}{M}{\sum }_{m=1}^{M}{\mathbb{I}}[{\hat{y}}_{m,{\ell },t,g} < {q}_{{\ell },t,g}^{{\rm{AIFS}}\,}(k)]\) from each AIFS ensemble forecast \({({\widehat{y}}_{m,{\ell },t,g})}_{m=1}^{M}\) with M = 100.

POET

The POET method²⁹ is a member-by-member post-processing technique that has been shown to significantly improve the medium-range weather forecasting skill of ECMWF dynamical forecasts and to outperform other statistical and AI-based post-processing techniques, including the member-by-member model output statistics post-processing technique³⁴ used operationally at the Royal Meteorological Institute of Belgium³⁵, ensemble model output statistics³⁶, LeNet neural network post-processing³⁷ and pure transformer-based post-processing³⁸. POET includes a self-attention mechanism across the ensemble member dimension in its encoder-decoder convolutional U-Net architecture, allowing the model to use context from the full input forecast distribution to minimize probabilistic forecasting error. We extended POET to subseasonal horizons by making a few minor modifications to its data-processing pipeline and architecture.

First, we applied POET to weekly aggregates of weather variables rather than the instantaneous values appropriate for medium-range forecasting. As input features, we used all available variables from the S2S database⁵² and a collection of 36 constant fields from ERA5³⁰, including orography, a land–sea mask and soil and vegetation types. Following common practice for feature embeddings, metadata, including the constant fields, lead time and weekly-averaged top-of-atmosphere incoming solar radiation, were embedded with a learned convolutional layer and used to update the highest-level features additively. To accommodate the larger set of input and output features, the model size was increased to approximately 27 million trainable parameters by increasing the number of latent features and the number of ensemble transformer layers.

POET is configured by default to predict several targets, including geopotential heights, wind and humidity on a few pressure levels and surface winds and humidity, in addition to the target variables temperature, mean sea-level pressure and precipitation. Although predicting several targets improves consistency between variables, we found that the single-target model produced higher-quality outputs for precipitation on our validation set. Hence, the final results reported for POET are from the multi-target model for temperature and pressure and from the single-target model for precipitation. The POET models were developed using the 20 years of 10-member ECMWF hindcasts issued in 2024, with the hindcasts for the years 2004–2022 used for training and the 2023 hindcasts used for validation.

For each target weather variable, subseasonal horizon and ECMWF forecast issuance date in 2024 (every Monday and Thursday to 21 November and every odd date thereafter), we generated weekly-aggregated ensemble forecasts, \({({\widehat{y}}_{m,{\ell },t,g})}_{m=1}^{M}\), by applying POET to the M = 100 perturbed ECMWF ensemble members. We then converted each ensemble into a CDF forecast \({F}_{{\ell },t,g}(k)=\frac{1}{M}{\sum }_{m=1}^{M}{\mathbb{I}}[{\widehat{y}}_{m,{\ell },t,g} < {q}_{t,g}(k)]\) using the climatological quintile thresholds q_t,g.

PBC

PBC (Fig. 1) corrects the systematic errors in probabilistic forecasts by training and applying two complementary machine learning corrections (Debias++ and Persistence++), projecting the output of each onto the space of valid distributions and finally averaging the resulting distributions. For a given weather variable, target date t^⋆ and lead time ℓ^⋆, each machine learning algorithm is given access to a model’s uncorrected CDF forecasts \({{\bf{F}}}_{{{\ell }}^{\star },{t}^{\star }}={({F}_{{{\ell }}^{\star },{t}^{\star },g})}_{g=1}^{G}\) for the target, a dataset of historical CDF forecasts \({{\bf{F}}}_{{\ell },t}={({F}_{{\ell },t,g})}_{g=1}^{G}\) from the model with training set target dates \(t\in {\mathcal{T}}\) and lead times \({\ell }\in {\mathcal{L}}\) and observed cumulative quantile bin indicators \({{\bf{O}}}_{t}={({O}_{t,g})}_{g=1}^{G}\) for each \(t\in {\mathcal{T}}\). To mimic operational deployment, each algorithm is only given access to training observations that would have been fully observable ℓ^⋆ − 1 days previous to the target date t^⋆.

Debias++

For each quantile bin k, Debias++ (Algorithm 1) performs probabilistic forecast correction in four steps. First, to account for heterogeneity in forecast errors over time, Debias++ adaptively selects a span s of observations around the target day of year for probabilistic debiasing and a range d^⋆ of forecast issuance dates for ensembling. Second, Debias++ forms ensemble mean probabilistic forecasts \({\overline{{\bf{F}}}}_{t}(k)\) by averaging its input CDF forecasts over the selected range of issuance dates. Third, to directly counteract systematic errors, Debias++ corrects the target CDF forecast \({\overline{{\bf{F}}}}_{{t}^{\star }}(k)\) by subtracting the probabilistic prediction bias—the average forecast value \({\overline{{\bf{F}}}}_{t}(k)\) minus the average observed value O_t(k)—over the selected debiasing window across the previous Y years. Finally, the results are projected into the valid probability range [0, 1].

In every application of Debias++, we set Y = 20, selected the configuration (s, d^⋆) ∈ {(14, 1), (28, 1), (35, 1)} that minimized mean RPS over the 3 years preceding the target date t^⋆ and converted each available ensemble forecast \({({\widehat{y}}_{m,{\ell },t,g})}_{m=1}^{M}\) and hindcast \({({\widehat{y}}_{m,{\ell },\delta ,t,g})}_{m=1}^{M}\) into an input CDF forecast using the observed climatological quintile thresholds q_t,g:

Algorithm 1

Debias++

input target date t^⋆; quantile bin k; lead time ℓ^⋆; span s; # issuance dates d^⋆; # training years Y; training set of observed cumulative quantile bin indicators \({({{\bf{O}}}_{t})}_{t\in {\mathcal{T}}}\) and cumulative probability forecasts \({({{\bf{F}}}_{{\ell },t})}_{{\ell }\in {\mathcal{L}},t\in {\mathcal{T}}\cup \{{t}^{\star }\}}\)

initialize days per year D = 365.242199

\({\mathcal{S}}=\left\{t\in {\mathcal{T}}:\mathrm{year}\_\mathrm{diff}:= \left\lfloor \frac{{t}^{\star }-t}{D}\right\rfloor \le Y\,{\rm{and}}\,\mathrm{day}\_\mathrm{diff}:= \frac{365}{2}-\left|\lfloor ({t}^{\star }-t)\,{\rm{mod}}\,D\rfloor -\frac{365}{2}\right|\le s\right\}\)

// Form ensemble mean probabilistic forecasts across issuance dates

for target and training dates \(t\in {\mathcal{S}}\cup \{{t}^{\star }\}\) do

\({\bar{{\bf{F}}}}_{t}(k)\) mean(\({({{\bf{F}}}_{{{\ell }}^{\star }-d+1,t}(k))}_{1\le d\le {d}^{\star }}\))

output clip (\({\bar{{\bf{F}}}}_{{t}^{\star }}(k)\) + mean(\({({{\bf{O}}}_{t}(k)-{\bar{{\bf{F}}}}_{t}(k))}_{t\in {\mathcal{S}}}\)), [0,1])

Persistence++

For each quantile bin k and grid cell g, Persistence++ (Algorithm 2) accounts for both recent weather trends and climatological shifts by regressing the observed indicator O_t,g(k) onto the historical model forecast \({F}_{{{\ell }}^{\star },t,g}(k)\), two recent and fully observable lagged indicators \({O}_{t-{{\ell }}^{\star }+6,g}(k)\) and \({O}_{t-2{{\ell }}^{\star }+5,g}(k)\) and a probabilistic climatology formed from the preceding 20 years of observed indicators with matching month and day. In the climatology computation, any missing indicators are imputed with the prior value \(\frac{k}{K}\). The learned regression weights are applied to the target forecast \({F}_{{{\ell }}^{\star },{t}^{\star },g}(k)\), lagged observations and climatology to produce a corrected output that is then projected into the valid probability range [0, 1].

When applying Persistence++ to POET forecasts or ECMWF precipitation forecasts, we again converted each available ensemble forecast and hindcast into an input CDF forecast using the observed climatological quintile thresholds. In all other applications of Persistence++, we used model-based quintile thresholds \({q}_{{\ell },t,g}^{{\rm{model}}}\) for a given forecast target date t to convert both the ensemble forecast \({({\widehat{y}}_{m,{\ell },t,g})}_{m=1}^{M}\) and the associated ensemble hindcasts \({({\widehat{y}}_{m,{\ell },\delta ,t,g})}_{m=1}^{M}\) into input CDF forecasts:

Algorithm 2

Persistence++

input target date t^⋆; quantile bin k; lead time ℓ^⋆; training set of observed cumulative quantile bin indicators \({({{\bf{O}}}_{t})}_{t\in {\mathcal{T}}}\) and cumulative probability forecasts \({({{\bf{F}}}_{{{\ell }}^{\star },t})}_{t\in {\mathcal{T}}\cup \{{t}^{\star }\}}\)

initialize forecast period length L = 7, # climatology years Y = 20

// Form rolling probabilistic climatology from observed indicators

for target and training dates \(t\in {\mathcal{T}}\cup \{{t}^{\star }\}\) do

C_t(k) = mean({O_s(k) : monthday(s) = monthday(t), 1 ≤ year(t) − year(s) ≤ Y})

// Combine probabilistic forecast, climatology and lagged observations

for grid points g = 1 to G do

\(\begin{array}{l}{\widehat{{\boldsymbol{\beta }}}}_{g}\in {{\rm{argmin}}}_{{\boldsymbol{\beta }}}{\sum }_{t\in {\mathcal{T}}}({O}_{t,g}(k)\\ \,-\,{{\boldsymbol{\beta }}}^{\top }[1,{C}_{t,g}(k),{O}_{t-{{\ell }}^{\star }-L+1,g}(k),{O}_{t-2{{\ell }}^{\star }-L+2,g}(k),{F}_{{{\ell }}^{\star },t,g}(k)]{)}^{2}\end{array}\)

output clip (\(({\widehat{{\boldsymbol{\beta }}}}_{g}^{\top }[1,{C}_{{t}^{\star },g}(k),{O}_{{t}^{\star }-{{\ell }}^{\star }-L+1,g}(k),{O}_{{t}^{\star }-2{{\ell }}^{\star }-L+2,g}(k),\)

    \({F}_{{{\ell }}^{\star },{t}^{\star },g}(k)]{)}_{g=1}^{G},\) [0, 1])

Projection and ensembling

It is possible for the quantile bin predictions of Debias++ and Persistence++ to violate the constraint that cumulative probabilities must be non-decreasing across quantile thresholds. In such cases, we use isotonic regression⁵⁶ to project their output forecasts onto the space of valid CDFs. This operation is guaranteed to improve forecast skill as, by minimizing squared error to the set of valid CDFs, it increases RPSS for all possible observations simultaneously. The final PBC forecast averages the two projected corrections to capture the complementary benefits of each.

MicroDuet

To form the MicroDuet model used in the AI Weather Quest, we considered three combinations of the PBC-ECMWF and PBC-POET forecasts—an equal-weighted average, full weight on PBC-ECMWF and full weight on PBC-POET—and selected the combination that maximized 2024 RPSS for each variable and lead time. The resulting selection gave equal weight to PBC-ECMWF and PBC-POET for temperature, full weight to PBC-ECMWF for precipitation and full weight to PBC-POET for mean sea-level pressure.

AI Weather Quest evaluation

For each target weather variable and model, Fig. 4 reports the AI Weather Quest leaderboard scores for the September–October–November contest season⁴¹. The Quest scored each model against ERA5-T ground truth³⁰ using the RPSS across all Monday target dates \({{\mathcal{T}}}^{\star }\) from 1 September 2025 to 24 November 2025 (inclusive) for week 3 and from 8 September 2025 to 1 December 2025 (inclusive) for week 4.

Forecasting extreme events

When forecasting extreme highs or extreme lows, each model is used to predict the probability of the top semidecile (measurements verifying above the 95th percentile) or the bottom semidecile (measurements verifying below the 5th percentile) respectively. Given a test set \({{\mathcal{T}}}^{\star }\) of target dates, we measured extreme forecast quality by means of the BSS⁵⁷,

where K = 20, k = 1 for extreme lows, k = K − 1 for extreme highs, BS_t,g is termed the Brier score (BS) of the forecast and \({\overline{{\rm{BS}}}}_{t,g}\) is the BS of a baseline climatological forecast that assigns equal probability to each semidecile bin.

Flood forecasting

To assess forecast skill during observed flood events, we constructed a global catalogue of flood occurrences from the GDACS (https://www.gdacs.org). Flood events (event type FL) at all alert levels (green, orange, red) were retrieved through the GDACS REST API for the date range 1 January 2016 to 29 May 2026. Each GDACS event record provides a centroid (latitude, longitude) coordinate and an event start date; from these we derived for every event (1) a geographic bounding box of ±10° around the centroid (clamped to [−90°, 90°] in latitude) and (2) a 7-day target period beginning on the event start date. Only events with start dates matching the target date of an ECMWF dynamical forecast (that is, Monday or Friday start dates until 15 July 2023 for week 3 or until 22 July 2023 for week 4 and all start dates thereafter) were evaluated. This yielded 2,435 flood events for week 3 and 2,428 for week 4. The severe flood analysis filtered these events to floods with orange or red alert levels yielding 88 events at each lead time.

For a given test set \({{\mathcal{T}}}^{\star }\) of flood event target dates, we measured flood forecast quality by means of the BSS (1) for the top precipitation semidecile bin (k = K − 1 for K = 20) with g now enumerating the bounding box grid cells.

Statistical methods

Paired Wilcoxon signed-rank tests were used to test for significant RPSS differences across the T = 26 weeks of the September–October–November and December–January–February AI Weather Quest seasons⁴². Because the exact issuance dates of the FuXi-S2S and ECMWF dynamical model seldom coincided in the 2017–2021 test period of Extended Data Fig. 6, we reported (two-sided) stationary block bootstrap⁵⁸ confidence intervals based on the bias-corrected and accelerated bootstrap⁵⁹ with automatic block length selection⁶⁰, 5,000 bootstrap replicates, T = 514 target dates for FuXi-S2S and T = 522 target dates for PBC-ECMWF.

For significance testing in Fig. 5c, we used paired one-sided bootstrap tests of BSS improvement using the bias-corrected and accelerated bootstrap with 5,000 bootstrap replicates, T = 88 severe flood events, T = 2,435 total week-3 flood events and T = 2,428 total week-4 flood events. For all other significance tests, we used paired one-sided stationary block bootstrap tests of skill score improvement using the bias-corrected and accelerated bootstrap with automatic block length selection, 5,000 bootstrap replicates and T target dates for T = 939 (Figs. 2 and 5a and Extended Data Figs. 2, 7a and 9a), T = 179 (Fig. 3a), T = 104 (Fig. 3b and Extended Data Fig. 7b, week 3), T = 107 (Fig. 3b and Extended Data Fig. 7b, week 4), T = 232 (Extended Data Fig. 1, December–January–February), T = 237 (Extended Data Fig. 1, March–April–May), T = 236 (Extended Data Fig. 1, June–July–August) and T = 234 (Extended Data Fig. 1, September–October–November).