# A Practical Probabilistic Benchmark for AI Weather Models

Noah D. Brenowitz  NVIDIA, Santa Clara, CA.  Corresponding author:
<EMAIL_ADDRESS>Yair Cohen  NVIDIA, Santa Clara, CA.  Jaideep Pathak
NVIDIA, Santa Clara, CA.  Ankur Mahesh  University of California, Berkeley, CA
Boris Bonev  NVIDIA, Santa Clara, CA.  Thorsten Kurth  NVIDIA, Santa Clara,
CA.  Dale R. Durran  NVIDIA, Santa Clara, CA.  Peter Harrington  Lawrence
Berkeley National Laboratory, Berkeley, CA  Michael S. Pritchard  NVIDIA,
Santa Clara, CA.

###### Abstract

Since the weather is chaotic, forecasts aim to predict the distribution of
future states rather than make a single prediction. Recently, multiple data
driven weather models have emerged claiming breakthroughs in skill. However,
these have mostly been benchmarked using deterministic skill scores, and
little is known about their probabilistic skill. Unfortunately, it is hard to
fairly compare AI weather models in a probabilistic sense, since variations in
choice of ensemble initialization, definition of state, and noise injection
methodology become confounding. Moreover, even obtaining ensemble forecast
baselines is a substantial engineering challenge given the data volumes
involved. We sidestep both problems by applying a decades-old idea—lagged
ensembles—whereby an ensemble can be constructed from a moderately-sized
library of deterministic forecasts. This allows the first parameter-free
intercomparison of leading AI weather models’ probabilistic skill against an
operational baseline. The results reveal that two leading AI weather models,
i.e. GraphCast and Pangu, are tied on the probabilistic CRPS metric even
though the former outperforms the latter in deterministic scoring. We also
reveal how multiple time-step loss functions, which many data-driven weather
models have employed, are counter-productive: they improve deterministic
metrics at the cost of increased dissipation, deteriorating probabilistic
skill. This is confirmed through ablations applied to a spherical Fourier
Neural Operator (SFNO) approach to AI weather forecasting. Separate SFNO
ablations modulating effective resolution reveal it has a useful effect on
ensemble dispersion relevant to achieving good ensemble calibration. We hope
these and forthcoming insights from lagged ensembles can help guide the
development of AI weather forecasts and have thus shared the diagnostic
code111http://github.com/NVIDIA/earth2mip.

##  1 Introduction

In recent years, several data-driven weather models have claimed deterministic
skill that are comparable (and even outperform) state of the art numerical
weather models. AI weather models are very computationally efficient,
suggesting that they could be ideal for probabilistic forecasting with large
ensembles.

Unfortunately, little is reported about the probabilistic skill of recent AI
weather models. In contrast, deterministic skill scores have been widely
reported – indeed, most AI forecasting systems have been optimized using
either mean squared error (Bonev et al., 2023; Lam et al., 2023) or mean
absolute error (Bi et al., 2023; Chen et al., 2023b) as a key part of the loss
function, and correspondingly their performance has often been compared to
state-of-the-art deterministic forecasts (Bi et al., 2023; Bonev et al., 2023;
Chen et al., 2023b; Lam et al., 2023). Yet this is unsatisfying as the Root
Mean Squared Error (RMSE) for a forecast from a perfect deterministic model
saturates at a value $\sqrt{2}$ higher than that for a perfect ensemble
forecast as derived theoretically by Leith (1974) and demonstrated for an
ensemble of AI weather models by Weyn et al. (2021). Any model claiming better
deterministic skill at long lead times does so by artificially reducing the
variance of the forecast, which will lead it to underestimate the chances of
unlikely events.

The purpose of this paper is to present a consistent and easy-to-use method
for evaluating deterministic weather forecasts from the lens of probabilistic
skill; to show how this clarifies the pitfalls of using only deterministic
scoring for model comparison, and to add complementary perspective to ongoing
discussions about how best to optimize AI weather predictions. We speculate
the lack of such a tool underlies the dearth of probabilistic score
intercomparisons in the literature, and stems from three key challenges.

The first is that fair comparison against an operational ensemble is an
overwhelming technical challenge, at least if one considers the typical
requirements of this baseline. The common operational practice is to run an
ensemble of deterministic forecasts to estimate the probability distribution
of future weather states. In such ensemble forecasts, the initial state (and
sometimes the model) are perturbed to account for unavoidable uncertainties.
Model-aware initialization techniques tend to require some knowledge of the
forecast model to come up with initial conditions at a time $t$, potentially
from the past history. For instance, some ensemble techniques like vector
breeding Toth and Kalnay (1997) require running the model while others like
singular value decomposition require knowledge of the model’s adjoint (Molteni
et al., 1996; Leutbecher, 2005). Model-dependent noise-injection techniques
are further used to enhance dispersion but to varying degrees across
operational centers Palmer et al. (2009). Perturbation techniques for
traditional weather models cannot be straightforwardly applied to AI weather
models. In short, fair comparison of deterministic weather predictions against
such an ensemble forecasting system is difficult as the perturbation
techniques and their tunable parameters are often model specific. The ensemble
methods used in operational centers have been tuned continuously over decades,
and there is no uniform ensemble perturbation method across numerical weather
prediction centers. That is, a parameter free intercomparison method would be
attractive.

The second challenge is that scoring ensemble forecasts over many lead times
requires large data volumes that are not easily transferred between various
groups that develop models. The amount of data of a modern ensemble of 50-100
members is prohibitive. Assuming we verify at $0.25^{\circ}$ resolution, with
12-hour sampling up to 14 days, and 600 initial times, then the baseline data
needed for scoring for a single ensemble member of one scalar field (e.g.
geopotential height at 500 hPa) is 70GB. For 50 members, this becomes TBs,
whereas the operational archives often have a peak download speed of around 10
MB/s. Thus, especially with the increasing demand from the growing community
experimenting with AI weather prediction, simply attaining the operational
baseline is a challenge. This data volume barrier will only grow higher as the
resolution of global forecasts increases.

The final challenge is ease and speed of ensemble scoring itself. Once
attained, a full operational ensemble baseline becomes a challenging amount of
data to use for the purpose of AI model validation. The score that guides AI
optimizations should be fast to compute, and ideally in an automated fashion
to provide rapid feedback to the ML model developer. Running large ensembles
is still slow for high-throughput scoring, even though data-driven models are
dramatically faster.

In summary, a solution having the following requirements would enable more
systematic probabilistic scoring of forecasts (either data-driven or
traditional).

  1. 1.

Use the same ensembling method for all models. Remove confounding factors.

  2. 2.

Minimize data requirements.

  3. 3.

Minimize tuning knobs. Keep it simple.

  4. 4.

Allow rapid feedback to the AI model developer.

We suggest such a solution already exists: A decades-old idea - lagged
ensemble forecasts (LEF) - can be used as a strategic ensemble of convenience.

Hoffman and Kalnay (1983) introduced LEF as an early ensemble forecast
technique, which treats all forecasts initialized within some time window of
the current time $t$ as separate ensemble members. Theoretically, it exploits
the fact that two deterministic forecasts initialized at different times
diverge from one another, similar to an initial condition perturbation, taking
advantage of the notion that a consistent perturbation of the initial state is
one which perturbs the atmospheric state while respecting the correlations in
space and between variables that are present in the data; as such, the
previous (or next) unperturbed state of the atmosphere can also be viewed as a
valid candidate of the perturbed state at the current time. Despite the fact
that it has been supplanted by more modern techniques, LEF is still useful in
various forecasting problems (Trenary et al., 2018; Vitart and Takaya, 2021;
Chen et al., 2013). Although the number of ensemble members is limited,
because they are drawn directly from the data distribution, no new tuning
knobs are added.

Practically, LEF have several properties that make them ideal for validating
and intercomparing the deterministic forecast models that form the core of an
operational ensemble forecast system. LEF provides a parameter free ensemble
forecasting system (Branković et al., 1990). As such it can be used to fairly
compare between different numerical and machine learning weather models.
Lagged ensembles allow one to compute a probabilistic score from deterministic
forecast archives, which are not prohibitive in size and can reasonably be
moved around the internet or computed on the fly from data-driven models. The
method is extremely simple, and most importantly the same for all models we
consider, though there may be some differences in the analyses used to start
the deterministic forecasts. The technique is especially valuable for
comparing ML weather models consistently against physical model baselines.

In the following sections, we apply LEF to intercompare the IFS and leading AI
weather models for the first time, to draw several insights. First, we
validate the relevance of LEF by showing that, applied to IFS, it can explain
most of the temporal variance of the actual operational ensemble (Sec. 3.1).
Then we apply LEF to re-examine the relative skill of two leading AI weather
models, to illustrate that a particular AI intervention – multi-step fine-
tuning - that leads to apparent gains in deterministic skill yields no
associated gains in probabilistic skill, consistent with the speculation of
(Lam et al., 2023) (Sec. 3.2). We then confirm this hypothesis by testing the
effects of multi-step fine-tuning across a detailed set of ablations with the
spherical Fourier neural operator (SFNO) AI weather model Bonev et al. (2023)
(Sec. 3.3). We then turn to assessing the role of effective resolution in the
context of SFNO architectures where we find that, similarly to fine tuning,
systematically decreasing the effective resolution also reduces the spread of
the lagged ensemble. Together, these results support our overall thesis – that
LEF provides efficient, practical guidance on probabilistic skill, with
potential to accelerate the development of state of the art probabilistic
weather forecasting systems.

##  2 Lagged Ensemble Benchmark

Figure 1: Overview of lagged ensemble forecasting. Each ensemble member
(color) is initialized at a different initial time (dots). The true time
series (-) and the lagged ensemble average (- -) are also shown.

Figure 1 shows a schematic of the lagged ensemble approach on a simple chaotic
dynamical system (see Sec A.1 for more information). The right-hand side shows
how forecasts initialized from different initial times soon diverge from one
another, when evaluated at common validation time. The ensemble mean has a
lower variance, but is clearly too smooth and unlike the truth.

To be more precise, let $x(t,t^{\prime})$ be the state of a simulation valid
at time $t$ and initialized at $t^{\prime}$. The lead time is $t-t^{\prime}$.
Let $h$ be the increment of the forecast interval (often 6 hours in the recent
deep-learning weather prediction literature). Then, we develop our scores
based on the centered $M$-lagged ensemble given by

$x_{m}(t,t^{\prime})=x(t,t^{\prime}-mh),\quad-M\leq m\leq M.$

This ensemble has $2M+1$ members and the initial times are centered around
$t^{\prime}$. While this ensemble cannot be used for a real weather forecast
since it requires initial conditions from the future, it is suitable for
comparing hindcasts. Assuming $t$ and $t^{\prime}$ are uniformly sampled with
$h$ simplifies the implementation. In practice a lagged ensemble can be
generated from any archive of deterministic hindcasts (Rasp et al., 2023a;
Swinbank et al., 2016), allowing probabilistic scoring on data volumes far
smaller than raw ensemble outputs, which Price et al. (2023) pointed out are
prohibitive to obtain.

We have implemented a distributed routine for scoring such ensembles using
https://github.com/NVIDIA/earth2mip, an open source software framework. This
tool supports scoring archives on disk, as well as online scoring of data-
driven weather models.

The lagged ensemble (and deterministic) forecasts are compared with the
corresponding atmospheric states given by ERA5 reanalysis Hersbach et al.
(2020) at $0.25^{\circ}$ resolution, which is also the training data for all
the ML models evaluated in this study. Hindcasts are made with initial times
taken every 12 hours from 2018-01-02 0Z through 2018-10-29 12Z, which is
outside the training set for all the models we consider. Like IFS HRES, we
consider 10-day forecasts.

There are some important caveats about our choice of ERA5 as verification.
First, to fairly evaluate models with different spatial resolution, we should
down-sample to the coarsest available resolution (Rasp et al., 2023b). Also,
ERA5 dataset provides fictitious numerical values on pressure levels that lie
below the surface, based on extrapolating a field downwards from the surface
in a field-dependent way. It does not seem fair to evaluate models on how well
they predict these unphysical values. We sidestep these concerns in the
present work by only considering models trained with these same data so we
leave it to future iterations to addresses these concerns.

We compare M-lagged ensembles with $h=$12\text{\,}\mathrm{h}$$ and $M=4$. The
choice of $M=4$ is somewhat arbitrary and $h$ is out of convenience since the
HRES data we had on premises was sub-sampled to this time-resolution. This
corresponds to a centered 2-day window around any given initial time with 9
ensemble members. We only analyze complete ensembles which start at day 2 and
end at day 8.

We use standard metrics to evaluate the performance of a given lagged
ensemble. First, a summary metric of the probability distribution penalizing
both over-confidence and inaccuracy is measured using the continuous-ranked-
probability score (CRPS)(Gneiting and Raftery, 2007). Second, we evaluate the
accuracy of the ensemble mean using the root-mean-square-error (eRMSE, i.e.
ensemble-RMSE). Third, we define the deterministic RMSE (dRMSE) as the RMSE of
the $m=0$ ensemble member. Finally, we diagnose how well calibrated a lagged
ensemble is by comparing the growth of the ensemble standard deviation
(spread) to the growth of the eRMSE. One ensemble is more dispersive than
another if its spread is higher for a given level of error (eRMSE). All
metrics are globally averaged and area-weighted. Section A.4 contains more
details.

For simplicity we assume all ensemble members are equally weighted. For short
lead time forecast, assigning larger weights to more recent lags is preferable
as in (Hoffman and Kalnay, 1983), but for the purpose of having a parameter-
free benchmark, we propose using uniform weights across the ensemble. This
somewhat complicates the spread-error and other calibration analyses since the
famous 1:1 spread-error relationship assumes ensemble members are
interchangeable with each other and the observations (Fortin et al., 2014).
Nonetheless, we do use the spread-error relationship to diagnose if an
ensemble is more or less dispersive than another, but will take care to avoid
value judgements like “under” and “over” dispersive. We will sometimes use the
IFS results as a target since the ECWMF operational forecast ensemble built on
top of IFS is well-calibrated. Section A.5 elaborates how we could extend this
analysis to weighted ensembles.

##  3 Results

###  3.1 Comparing a lagged and operational ensemble

It is natural to wonder to what extent the LEF method can be used as a valid
proxy of the operational ensemble. To validate this idea, we compare the
lagged ensemble of IFS against the ECMWF operational ensemble ECMWF (2019a),
see also A.2. By comparing the best deterministic and ensemble forecasting
systems from within a single modeling center, we remove many confounding
factors in the analysis, although it is worth noting that the IFS and HRES
ensembles were generated at different resolutions and with different
configurations.

Figure 2 confirms the validity of LEF as a proxy with meaningful operational
relevance. At a lead time of 5 days, i.e. sufficiently long that no member of
the LEF contains the observations, the CRPS of the lagged and operational
ensembles are well-correlated across initial dates: much of the day-to-day
variation in global CRPS of the operational ensemble can be predicted from
that of the lagged ensemble. As expected, the lagged ensemble is more skillful
on average, since it has the artificial advantage that some of its members
have shorter lead times than the operational ensemble. We have confirmed that
the relationship between LEF CRPS and operational CRPS is similarly monotonic
for other variables and shorter lead times (not shown).

In short, even though there are many differences between these ensemble
methods, Figure 2 supports our working hypothesis that LEF is a meaningful
probabilistic benchmark of operational relevance. We thus proceed to employ it
to perform the first clean probabilistic intercomparison of AI weather models
against the IFS LEF baseline.

Figure 2: Comparison between the global mean CRPS of
$\mathrm{5}\mathrm{0}\mathrm{0}$hpa height forecasts from a lagged ensemble
and the IFS operational ensemble at a lead time of 5 days. Each dot shows a
different initial time. The units of each axes are m^2/s^2

.

###  3.2 Comparing data-driven and NWP models with the same ensemble method

Figure 3 compares the ensemble metrics for a set of forecast models including
GraphCast, Pangu Weather, and IFS HRES (see Sec. A.3 for details). From the
vantage of deterministic skill scores (left column) both of these AI weather
models appear to outperform the IFS, and GraphCast especially and increasingly
so as a function of lead time. However, the ensemble skill scores – CRPS and
RMSE of the ensemble mean – tell a different story. All models are virtually
tied, subtle differences between the models’ ensemble RMSE shrink rather than
grow with lead time, and (bottom right) Pangu Weather actually has the lowest
Z500 CRPS for many lead times. Overall, these results demonstrate a non-
monotonic relationship between deterministic skill and the performance of an
ensemble. Also we have shown that some data-driven models, despite being
optimized on deterministic loss functions, can outperform physics-based models
such as IFS when the same ensemble technique is used.

Note that results for lead times shorter than two days are omitted from our
analysis because IFS is initialized with analysis (i.e. knows only about past
weather) but we evaluate it against reanalysis (which is inform also by future
weather). This leads to an initial Z500 error of $30\text{\,}\mathrm{m}$, but
for lead times of 2 days and longer this effect is an order of magnitude
smaller (Rasp et al., 2023b, Fig S6).

Figure 3: Z500 and T850 skill comparison for several deterministic forecast
systems. From left to right, deterministic RMSE (dRMSE), ensemble RMSE
(eRMSE), ensemble spread and CRPS. Ensemble scores are only valid between 2
and 8 days.

###  3.3 The effect of long-lead time training

Many of the existing ML weather models employ some method of aggregating loss
over time longer than 6h, though the exact methods differ, see Sec. A.3 and
Table 1. Focusing on the two leading models (GraphCast and Pangu), LEF scoring
shows that GraphCast is less dispersive than the other models (Figure 3c,g).
That is, for a given level of ensemble mean skill, the ensemble has less
variance.

We hypothesize that this is related to the multi-step training protocol and
loss function used to train it, wherein successive amounts of auto-regressive
time-stepping are unfolded within the optimizer (Brenowitz and Bretherton,
2018; Weyn et al., 2020; Lam et al., 2023), out to lead times of 3 days.
Indeed, the minimizer of the MSE loss is known to be the conditional mean, so
an ML model trained for long lead times will remove all variance that is not
predictable on that time scale.

LEF confirms this hypothesis. We begin by training SFNO models with varying
amounts of long-lead time fine-tuning of the autoregressive training. The base
model is trained for 80 epochs with a context length of a single 6-hour
timestep. We also include results with 2, 3, 5, 9, and 11 steps of
autoregressive timestepping in the MSE loss (12h-3 days) to illustrate the
effect. The pre-trained base model is fine-tuned with a smaller learning rate
for 20 epochs with a context length of two 6-hour timesteps. The output from a
single forward pass serves as the input for the model to generate a 12-hour
forecast. The total accumulated loss from the 6 and 12-hour forecasts is then
used as the training objective. This fine-tuning procedure is sequentially
extended, with the context length progressively increasing from 18 hours (3
timesteps) to 3 days (11 timesteps). For each increment in context length, the
model is trained for 6 epochs. NVIDIA has previously released the 11-step
model with the name “FourCastNet v2-small” 222This model is available on
NVIDIA’s NGC Catalog and also on ECMWF’s ai-models platform.

Figure 4 depicts the results of our LEF analysis applied to this sensitivity
test for the 850 hPa temperature field. The deterministic RMSE decreases
significantly for models trained with more autoregressive timesteps, even
appearing to approach the threshold of IFS deterministic skill limits (panel
a). However, the CRPS and ensemble mean RMSE improve much less than the
deterministic RMSE. The reason is due to plummeting ensemble spread (see
Figure 5a-d) where for a given eRMSE the SFNO checkpoints optimized with long-
lead time fine-tuning result in monotonic drops in ensemble variance,
deviating from the approximately 1:1 error-spread calibration of the LEF IFS
ensemble.

Overall, these results show that autoregressive training behaves similarly to
ensemble averaging – it improves skill by reducing variance. This view is
consistent with Lam et al. (2023) who also noted that long-lead time training
came at the expense of a spatially blurred forecast, speculating it optimizes
for something in between a deterministic and ensemble mean forecast, a point
that has since been elucidated in detail by Price et al. (2023) through
comparison against more directly probabilistic diffusion-based forecasts. We
emphasize that unlike the latter, LEF provides a convenient parameter-free way
to trap such issues without the complications of varying approaches to initial
condition perturbation.

We also note a relationship between ensemble and time averaging. Since chaotic
dynamical systems like the weather forget their initial condition
$\lim_{T\rightarrow\infty}E[x(t)|x(t-T)]=E[x(t)]$. This implies that the
climate mean $E[X]$ is the best possible forecast for long lead times. If we
further assume ergodicity, one can use time averages of the observed
record—known as climatology—to approximate $E[X]$. To sum up, time averaging,
long-lead-time training, and ensemble averaging are roughly equivalent in
their ability to reduce the variance of a forecast and thus improve its
deterministic skill.

Figure 4: Using LEF to explore the effect of fine-tuning. Deterministic RMSE
(a), ensemble mean RMSE (b) and CRPS (c) as a function of lead time. Ensemble
standard deviation as a function of ensemble mean RMSE (d). The color lines
show a range of fine tuning steps and the dashed line shows the ifs. The field
being scored is temperature at 850 hPa (T850).

###  3.4 Sensitivity to effective resolution

We now show that it is possible to make a data-driven forecast more dispersive
by tuning the architecture. For computational reasons, many ML models (see
Sec. 1) project the 0.25∘ inputs into a coarse latent space early in the
network architecture. We suspect the resolution of this latent space limits
the resolutions the architecture can faithfully generate in autoregressive
forecast.

We explore this hypothesis by varying the user-specified scale-factor
parameter of SFNO. In the SFNO architecture, the first layer downsamples the
input by the scale-factor; larger scale factors correspond to lower spatial
resolutions. We train models with scale factors of 6, 4, 3, 2, and 1. Each
model was trained on 64 NVIDIA A100 GPUs. SFNOs with scale factor 6, 4, 3, 2,
and 1 took 6.9, 7.1, 7.1, 8, and 89 hours to train, respectively. Figure 5e-h
shows that SFNOs trained with lower scale factors are more dispersive for
several fields (higher spread-error ratio). While these changes did not
significantly improve the skill of the model, we suspect that further tuning
of the hyper-parameters, especially the embedding dimension, could improve the
results.

Figure 5: Sensitivity of the ensemble spread to changes in hyperparameters.
(a-d) spread/skill results for various fields and differing amounts of auto-
regressive fine tuning. (e-h) shows the same, but with a varying scale factor.

##  4 Discussion and Conclusions

For 60 years weather modelers have known that we can reduce the RMSE of target
variables in a deterministic weather forecast by reducing their variance
(Leith, 1974). Yet, the developers of physics-based models typically seek to
retain sharpness in single-model forecasts and to increase the variance in
ensembles of their model forecasts because those ensembles are typically
undispersive. As part of this effort, benchmarks are devised for the dynamical
core alone (Jablonowski and Williamson, 2006; Thatcher and Jablonowski, 2016)
to select the numerical methods with the least numerical diffusion, and the
resolution is increased in-order to represent physical processes with
increasing detail. Model developers also examine case studies of extreme
storms to ensure they produce realistic winds and rainfall (Shahrban et al.,
2016).

On the other hand, much of the data-driven forecast literature has focused
primarily on minimizing the error in deterministic forecasts. The easiest way
to do this is to reduce the variance of a forecast. It is easy to recognize
ensemble averaging or space-time filtering as variance-reduction techniques.
However, there are a large number of other ML design choices whose effects on
ensemble dispersion are less obvious. For instance, we have identified that
effective resolution can be a helpful lever on dispersion separately from
determinstic skill.

In this work, we devise a parsimonious framework that allows one to compare
physics-based and data-driven models with the same ensemble technique. Within
this framework, we highlight that the choice of loss function can produce
artificial dissipation akin to ensemble averaging. Many recent papers have
used some form of training involving multi-day lead times including, for
example, auto-regressive MSE losses(Lam et al., 2023; Pathak et al., 2022),
replay buffers (Chen et al., 2023a, b), or simply predicting very long lead
times directly Nguyen et al. (2023). Our results show that, when taken too
far, such loss functions lead to a poorly calibrated ensemble.

While lagged ensemble scoring is a practical framework for comparing
deterministic forecasts, several improvements are possible. First of all,
lagged ensemble members are not interchangeable so we know that a weighted
ensemble will outperform our unweighted ensembles here. The analysis of
ensemble calibration is made more ambiguous by this subtlety; we haved use the
LEF of IFS as a baseline to circumvent this issue. Future work can also focus
on applying this technique to more models and working to integrate it into
existing community benchmarks which provide a uniform protocol for e.g. output
resolution (Rasp et al., 2023b).

Finally, this work tests a deterministic forecast in isolation from other
components. This is most useful if we assume that the probabilistic forecast
systems of the future will be built as they are today: from ensembles of
deterministic time steppers, but recent work on intrinsically stochastic
models suggests an alternative path (Price et al., 2023; Kochkov et al.,
2023). While we could score such models with lagged ensembles, this benchmark
may give an unfair advantage to deterministic models that have not yet solved
the more complex end-to-end probabilistic prediction task.

## Acknowledgments

The computation for this paper was supported in part by the DOE Advanced
Scientific Computing Research (ASCR) Leadership Computing Challenge (ALCC)
2023-2024 award ‘Huge Ensembles of Weather Extremes using the Fourier
Forecasting Neural Network’ to William Collins (LBNL). This research was
supported by the Director, Office of Science, Office of Biological and
Environmental Research of the U.S. Department of Energy under Contract No. DE-
AC02-05CH11231 and by the Regional and Global Model Analysis Program area
within the Earth and Environmental Systems Modeling Program. The research used
resources of the National Energy Research Scientific Computing Center (NERSC),
also supported by the Office of Science of the U.S. Department of Energy,
under Contract No. DE-AC02-05CH11231. We further acknowledge Anima Anandkumar
and Kamyar Azizzadenesheli for their helpful discussions.

## References

  * Bi et al. (2023) Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian.  Accurate medium-range global weather forecasting with 3d neural networks.  _Nature_ , 619(7970):533–538, 2023.
  * Bonev et al. (2023) Boris Bonev, Thorsten Kurth, Christian Hundt, Jaideep Pathak, Maximilian Baust, Karthik Kashinath, and Anima Anandkumar.  Spherical fourier neural operators: Learning stable dynamics on the sphere.  _arXiv preprint arXiv:2306.03838_ , 2023.
  * Branković et al. (1990) Č Branković, TN Palmer, F Molteni, S Tibaldi, and U Cubasch.  Extended-range predictions with ecmwf models: Time-lagged ensemble forecasting.  _Quarterly Journal of the Royal Meteorological Society_ , 116(494):867–912, 1990.
  * Brenowitz and Bretherton (2018) N D Brenowitz and C S Bretherton.  Prognostic validation of a neural network unified physics parameterization.  _Geophys. Res. Lett._ , 17:2493, June 2018.  ISSN 0094-8276.  doi: 10.1029/2018GL078510.
  * Chen et al. (2023a) Kang Chen, Tao Han, Junchao Gong, Lei Bai, Fenghua Ling, Jing-Jia Luo, Xi Chen, Leiming Ma, Tianning Zhang, Rui Su, et al.  Fengwu: Pushing the skillful global medium-range weather forecast beyond 10 days lead.  _arXiv preprint arXiv:2304.02948_ , 2023a.
  * Chen et al. (2023b) Lei Chen, Xiaohui Zhong, Feng Zhang, Yuan Cheng, Yinghui Xu, Yuan Qi, and Hao Li.  Fuxi: A cascade machine learning forecasting system for 15-day global weather forecast.  _npj Climate and Atmospheric Science_ , 6:190, 2023b.  doi: 10.1038/s41612-023-00325-5.
  * Chen et al. (2013) Mingyue Chen, Wanqiu Wang, and Arun Kumar.  Lagged ensembles, forecast configuration, and seasonal predictions.  _Monthly weather review_ , 141(10):3477–3497, 2013.
  * ECMWF (2019a) ECMWF.  _IFS Documentation CY46R1 - Part V: Ensemble Prediction System_.  Number 5. 2019 2019a.  doi: 10.21957/38yug0cev.  URL https://www.ecmwf.int/node/19309.
  * ECMWF (2019b) ECMWF.  _IFS Documentation CY46R1 - Part V: Ensemble Prediction System_.  Number 5. 2019 2019b.  doi: 10.21957/38yug0cev.  URL https://www.ecmwf.int/node/19309.
  * Fortin et al. (2014) V Fortin, M Abaza, F Anctil, and R Turcotte.  Why should ensemble spread match the RMSE of the ensemble mean?  _J. Hydrometeorol._ , 15(4):1708–1713, August 2014.  ISSN 1525-755X, 1525-7541.  doi: 10.1175/JHM-D-14-0008.1.
  * Gneiting and Raftery (2007) Tilmann Gneiting and Adrian E Raftery.  Strictly proper scoring rules, prediction, and estimation.  _J. Am. Stat. Assoc._ , 102(477):359–378, March 2007.  ISSN 0162-1459.  doi: 10.1198/016214506000001437.
  * Hersbach et al. (2020) Hans Hersbach, Bill Bell, Paul Berrisford, Shoji Hirahara, András Horányi, Joaquín Muñoz-Sabater, Julien Nicolas, Carole Peubey, Raluca Radu, Dinand Schepers, et al.  The era5 global reanalysis.  _Quarterly Journal of the Royal Meteorological Society_ , 146(730):1999–2049, 2020.
  * Hoffman and Kalnay (1983) Ross N Hoffman and Eugenia Kalnay.  Lagged average forecasting, an alternative to monte carlo forecasting.  _Tellus A_ , 35A(2):100–118, March 1983.  ISSN 0280-6495, 1600-0870.  doi: 10.1111/j.1600-0870.1983.tb00189.x.
  * Jablonowski and Williamson (2006) C Jablonowski and DL Williamson.  A baroclinic wave test case for dynamical cores of general circulation models: Model intercomparisons. ncar tech.  Technical report, Note NCAR/TN-469+ STR, National Center for Atmospheric Research, Boulder …, 2006.
  * Keisler (2022) Ryan Keisler.  Forecasting global weather with graph neural networks.  _arXiv preprint arXiv:2202.07575_ , 2022.
  * Kochkov et al. (2023) Dmitrii Kochkov, Janni Yuval, Ian Langmore, Peter Norgaard, Jamie Smith, Griffin Mooers, James Lottes, Stephan Rasp, Peter Düben, Milan Klöwer, et al.  Neural general circulation models.  _arXiv preprint arXiv:2311.07222_ , 2023.
  * Lam et al. (2023) Remi Lam, Alvaro Sanchez-Gonzalez, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Ferran Alet, Suman Ravuri, Timo Ewalds, Zach Eaton-Rosen, Weihua Hu, et al.  Learning skillful medium-range global weather forecasting.  _Science_ , page eadi2336, 2023.
  * Leith (1974) C. E. Leith.  Theoretical Skill of Monte Carlo Forecasts.  _Monthly Weather Review_ , 102:409–418, 1974.  doi: 10.1175/1520-0493(1974)102<0409:TSOMCF>2.0.CO;2.
  * Leutbecher (2005) Martin Leutbecher.  On ensemble prediction using singular vectors started from forecasts.  _Monthly weather review_ , 133(10):3038–3046, 2005.
  * Lorenz (1963) Edward N Lorenz.  Deterministic nonperiodic flow.  _J. Atmos. Sci._ , 20(2):130–141, March 1963\.  ISSN 0022-4928.  doi: 10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2.
  * Molteni et al. (1996) F Molteni, R Buizza, T N Palmer, and T Petroliagis.  The ECMWF ensemble prediction system: Methodology and validation.  _Quart. J. Roy. Meteor. Soc._ , 122(529):73–119, January 1996.  ISSN 0035-9009, 1477-870X.  doi: 10.1002/qj.49712252905.
  * Nguyen et al. (2023) Tung Nguyen, Johannes Brandstetter, Ashish Kapoor, Jayesh K Gupta, and Aditya Grover.  Climax: A foundation model for weather and climate.  _arXiv preprint arXiv:2301.10343_ , 2023.
  * Palmer et al. (2009) Tim N Palmer, Roberto Buizza, F Doblas-Reyes, Thomas Jung, Martin Leutbecher, Glenn J Shutts, Martin Steinheimer, and Antje Weisheimer.  Stochastic parametrization and model uncertainty.  2009\.
  * Pathak et al. (2022) Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, et al.  Fourcastnet: A global data-driven high-resolution weather model using adaptive fourier neural operators.  _arXiv preprint arXiv:2202.11214_ , 2022.
  * Price et al. (2023) Ilan Price, Alvaro Sanchez-Gonzalez, Ferran Alet, Timo Ewalds, Andrew El-Kadi, Jacklynn Stott, Shakir Mohamed, Peter Battaglia, Remi Lam, and Matthew Willson.  Gencast: Diffusion-based ensemble forecasting for medium-range weather.  _arXiv preprint arXiv:2312.15796_ , 2023.
  * Raftery et al. (2005) Adrian E Raftery, Tilmann Gneiting, Fadoua Balabdaoui, and Michael Polakowski.  Using bayesian model averaging to calibrate forecast ensembles.  _Mon. Weather Rev._ , 133(5):1155–1174, May 2005\.  ISSN 0027-0644.  doi: 10.1175/MWR2906.1.
  * Rasp and Thuerey (2021) Stephan Rasp and Nils Thuerey.  Data-driven medium-range weather prediction with a resnet pretrained on climate simulations: A new model for weatherbench.  _Journal of Advances in Modeling Earth Systems_ , 13(2):e2020MS002405, 2021.
  * Rasp et al. (2023a) Stephan Rasp, Stephan Hoyer, Alexander Merose, Ian Langmore, Peter Battaglia, Tyler Russel, Alvaro Sanchez-Gonzalez, Vivian Yang, Rob Carver, Shreya Agrawal, Matthew Chantry, Zied Ben Bouallegue, Peter Dueben, Carla Bromberg, Jared Sisk, Luke Barrington, Aaron Bell, and Fei Sha.  Weatherbench 2: A benchmark for the next generation of data-driven global weather models, 2023a.
  * Rasp et al. (2023b) Stephan Rasp, Stephan Hoyer, Alexander Merose, Ian Langmore, Peter Battaglia, Tyler Russel, Alvaro Sanchez-Gonzalez, Vivian Yang, Rob Carver, Shreya Agrawal, et al.  Weatherbench 2: A benchmark for the next generation of data-driven global weather models.  _arXiv preprint arXiv:2308.15560_ , 2023b.
  * Shahrban et al. (2016) Mahshid Shahrban, Jeffrey P Walker, QJ Wang, Alan Seed, and Peter Steinle.  An evaluation of numerical weather prediction based rainfall forecasts.  _Hydrological Sciences Journal_ , 61(15):2704–2717, 2016.
  * Swinbank et al. (2016) Richard Swinbank, Masayuki Kyouda, Piers Buchanan, Lizzie Froude, Thomas M Hamill, Tim D Hewson, Julia H Keller, Mio Matsueda, John Methven, Florian Pappenberger, et al.  The tigge project and its achievements.  _Bulletin of the American Meteorological Society_ , 97(1):49–67, 2016.
  * Thatcher and Jablonowski (2016) Diana R Thatcher and Christiane Jablonowski.  A moist aquaplanet variant of the held–suarez test for atmospheric model dynamical cores.  _Geoscientific Model Development_ , 9(4):1263–1292, 2016.
  * Toth and Kalnay (1997) Zoltan Toth and Eugenia Kalnay.  Ensemble forecasting at ncep and the breeding method.  _Monthly Weather Review_ , 125(12):3297–3319, 1997.
  * Trenary et al. (2018) L Trenary, T DelSole, MK Tippett, and K Pegion.  Monthly enso forecast skill and lagged ensemble size.  _Journal of Advances in Modeling Earth Systems_ , 10(4):1074–1086, 2018.
  * Virtanen et al. (2020) Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17:261–272, 2020.  doi: 10.1038/s41592-019-0686-2.
  * Vitart and Takaya (2021) Frédéric Vitart and Yuhei Takaya.  Lagged ensembles in sub-seasonal predictions.  _Quarterly Journal of the Royal Meteorological Society_ , 147(739):3227–3242, 2021.
  * Weyn et al. (2019) Jonathan A Weyn, Dale R Durran, and Rich Caruana.  Can machines learn to predict weather? using deep learning to predict gridded 500-hpa geopotential height from historical weather data.  _Journal of Advances in Modeling Earth Systems_ , 11(8):2680–2693, 2019.
  * Weyn et al. (2020) Jonathan A Weyn, Dale R Durran, and Rich Caruana.  Improving data-driven global weather prediction using deep convolutional neural networks on a cubed sphere.  _J. Adv. Model. Earth Syst._ , 12(9), September 2020.  ISSN 1942-2466.  doi: 10.1029/2020ms002109.
  * Weyn et al. (2021) Jonathan A Weyn, Dale R Durran, Rich Caruana, and Nathaniel Cresswell-Clay.  Sub-seasonal forecasting with a large ensemble of deep-learning weather prediction models.  _Journal of Advances in Modeling Earth Systems_ , 13(7):e2021MS002502, 2021.

##  Appendix A Methods

###  A.1 Lorenz 1963 Dynamical System

Lorenz [1963] introduced a simple differential equation model for chaos in
atmospheric convection. The dynamics are given by

$\displaystyle\frac{dx}{dt}$ $\displaystyle=\sigma(y-x),$
$\displaystyle\frac{dy}{dt}$ $\displaystyle=x(\rho-z)-y,$
$\displaystyle\frac{dz}{dt}$ $\displaystyle=xy-\beta z.$

In Fig. 1 the truth time series is generated using the
`scipy.integreate.odeint` routine from SciPy [Virtanen et al., 2020] with the
parameters set to $\rho=28$, $\sigma=10.0$, $\beta=8/3$. The lagged ensemble
is generated by forecasts initialized from this true time series without noise
perturbation but by setting $\rho=28.14$. We interpret this as a model error
similar to that made by an imperfect weather forecast model. The trajectories
of the lagged ensemble would not diverge without this model error since the
dynamics are deterministic and the initial conditions are unperturbed.

###  A.2 Numerical models

The European Centre Medium Range Forecast’s (ECMWF) hosts the world’s leading
global weather forecasting systems. Its deterministic model, the High
RESolution forecast (HRES) is a configuration of the Integrated Forecast
System (IFS) which produces global 10-day forecasts at about $0.1^{\circ}$
resolution ECMWF [2019b]. It is used operationally as the high resolution
deterministic member of an ensemble forecasting system (ENS) where additional
50 members at $0.25^{\circ}$ are used to generate the operational forecast.

###  A.3 Data-driven models

In the past three years a large number of data driven weather models have been
published. Table 1 provides some key information about these models. These
models share the training data of ERA5 [Hersbach et al., 2020] and for the
most part employ either MAE or MSE loss functions (except FangWu [Chen et al.,
2023a]). Earlier models use a lower resolution [Weyn et al., 2019, Rasp and
Thuerey, 2021, Keisler, 2022], while more recent models explore the full ERA5
resolution of $0.25^{\circ}$. To save memory and compute, many architectures
reduce the resolution of this data internally by embedding it into a coarser
latent space. A common practice is to coarsen the ERA5 input data by a factor
(typically 4) while predicting at full resolution (see the "scale factor")
column in Table 1. The details (and terminology) of such coarsening depends on
the architecture, e.g. a spherical harmonics based model like SFNO uses a
scale factor, effectively limiting the maximum wave number, GraphCast uses the
embedded spatial resolution on its icosahedral grid while a ViT based model
like PanguWeather uses patch embedding.

Many models are trained with autoregressive loss accumulation of different
kinds (see Table 1). Brenowitz and Bretherton [2018] introduced this idea to
stabilize machine learning parameterizations, and then Weyn et al. [2021]
applied it to global weather prediction. Since then, severals works have also
used autoregressive loss, including SFNO [Bonev et al., 2023] and GraphCast
[Lam et al., 2023]. The latter reported a comparison of the effect of number
of steps on the RMSE (see figure S34 there) and also note that auto-regressive
training encourages GraphCast to blur the forecast. PanguWeather [Bi et al.,
2023] took a different approach by training several models at different time-
steps (namely 1,3,6,24 hours) and interleaving them in the forecast. FuXi Chen
et al. [2023b] combines both methods in the so-called cascade approach where
separate models are trained and fine-tuned for different parts of the medium
range (0-14 days lead time). FengWu [Chen et al., 2023a] used a “replay
buffer” in which the model is trained to predict the weather at time $t$ from
both target at time $t-1$ and from stored model forecasts at time $t-1$;
effectively learning to veer back to the target data when provided with
initial condition with a lead time error accumulated over several steps. This
in practice could encourage the model to learn an approximation of the
ensemble mean as well.

Furthermore, most models are evaluated with deterministic metrics (RMSE and
ACC), exceptions are reported in Bi et al. [2023], Chen et al. [2023a].

Table 1: A list of the most relevant global ML weather models. Here NLL used
in Chen et al. [2023a] is a negative log-likelihood loss.

Feature |  DLWP |  FCN-SFNO |  GraphCast |  Pangu |  FuXi |  Keisler |  FengWu
---|---|---|---|---|---|---|---
Architecture |  CNN |  SFNO |  GNN |  swin |  Transformer |  Enc-Dec |  Transformer
Parameters |  $0.7\times 10^{6}$ |  $829\times 10^{6}$ |  $36\times 10^{6}$ |  $64\times 10^{6}$ |  – |  $6.7\times 10^{6}$ |  –
Resolution |  $2^{\circ}$ |  $0.25^{\circ}$ |  $0.25^{\circ}$ |  $0.25^{\circ}$ |  $0.25^{\circ}$ |  $1^{\circ}$ |  $0.25^{\circ}$
Loss |  MSE |  MSE |  MSE |  MAE |  MAE |  MSE |  NLL
Autoregressive training |  4 $\times$ 6h steps |  12 $\times$ 6h steps |  12 $\times$ 6h steps |  interleaving 1,3,6,24h |  12 $\times$ 6h steps + interleaving |  10 $\times$ 6h steps |  Replay Buffer
Scale factor |  4 |  6 |  4 |  4 |  – |  – |  –
Metrics |  RMSE, ACC |  RMSE, ACC |  RMSE, ACC |  RMSE, ACC, CRPS, STD/RMSE |  RMSE, ACC, CRPS, STD, SSR |  RMSE |  ACC, RMSE

#### Spherical Fourier Neural Operator (SFNO)

Bonev et al. [2023] introduced a spherical Neural Operator, which respects the
underlying geometry of the sphere and its inherent symmetries. The resulting
architecture has remarkable stability properties with reported autoregressive
rollouts of up to a full year Bonev et al. [2023].

Fourier Neural Operators are closely related to convolutional networks, where
convolutions are computed in the continuous domain via the spectral theorem.
This makes Neural Operators fundamentally agnostic of the underlying
discretisation of input and output signals and enables them to be trained and
evaluated on any discretisation that permits the computation of the
operations.

In the following, we utilize a range of SFNO models, which differ in model
resolution via a scale factor. This scale factor effectively controls the
downsampling of the signal and the resolution of the internal representation.
A scale factor of one corresponds to no downsampling, whereas a scale factor
of two corresponds to a 2x downsampling in the internal representation. Higher
downsampling ratios therefore lead to lower effective resolutions and a higher
smoothing of the predicted fields.

#### Pangu Weather

Pangu Weather Bi et al. [2023] introduced a hierarchy of vision transformers
trained to take different time steps. The publicly available models include

  * •

3 hour time step,

  * •

6 hour time step, and

  * •

24-6 interleaved time step.

The interleaved time stepper is effectively a 24-hour model, with 6-hour steps
used to interleave in between. For example, to make a 72 hour forecast,
requires 3 24 hour steps and no 6 hour steps. A 36 hour forecast requires a
single 24 time step and two 6 hour steps.

###  A.4 Ensemble Verification Metrics

Probabilistic predictions aim to maximize sharpness subject to calibration
Raftery et al. [2005]. Qualitatively, calibration means that the likelihood of
observing the true value is the same as observing a member drawn from the
ensemble. A necessary condition for calibration is that the spread-error
relationship be 1-to-1 when averaged over a sufficient sample [Fortin et al.,
2014]. With this condition in hand, reducing RMSE of the ensemble mean is the
same as reducing the ensemble spread.

The calibrate-sharpness perspective suggests that ensemble-mean RMSE is
equivalent to the spread of an ensemble calibrated by variance inflation. In
the limit of many ensemble members artificially increasing the spread won’t
change the ensemble mean much, since the variance of the ensemble mean is
given by $\sigma^{2}/n=o(1/n)$. Therefore the spread (i.e. sharpness) of the
calibrated ensemble is the same as the ensemble-mean RMSE. Thus, the ensemble
with the lowest ensemble-mean RMSE regardless of spread is the better ensemble
when viewed from the calibrate-sharpness perspective.

On the other hand, not all ensembles are calibrated. For these, we can compare
the performance using proper scoring rules like continuous-ranked-probability
score (CRPS)[Gneiting and Raftery, 2007]. CRPS is given by

$\int(F(x)-H(x-y))^{2}dx$

where $H(x)$ is the heaviside function and $F(x)=n^{-1}\sum_{k}H(x-x_{k})$ is
the empirical cumulative distribution function, and $y$ is the observation.

We present area-weighted global averages of this metric, where the area
weighting is given by $\cos(lat)$. Unlike Rasp et al. [2023b], for RMSE and
spread, initial time averages are taken across the lagged ensemble at fixed
lead time before the square root, so that the square of RMSE is a time and
space-averaged MSE.

###  A.5 Spread-error correspondence for a weighted ensemble

[Fortin et al., 2014] derived the famous 1-1 spread-error relationship
expected assuming the ensemble members and observations are interchangeable.
This is true for initial conditions ensembles, but the members of a lagged
ensemble are not interchangeable since the more recent members are more
skillful.

Bayesian model averaging (BMA) provides a more general formalism for such
cases [Raftery et al., 2005]. Specifically, let $f_{m}(t_{j})$ be the forecast
valid at time $t_{j}$ but initialized with a lag $m$. The initial time of this
forecast is $t_{j}-h(\ell+m)$ where $\ell h$ is the lead time of the forecast.
Let $S_{j}=\bigcup_{m}\\{f_{m}(t_{j})\\}$ be the set of such forecasts. BMA
assumes the distribution of a variable $Y$ drawn from this ensemble is given
by

$p(y|S_{j})=\sum_{m}w_{m}p(y|f_{m}),$ (1)

where $w_{m}$ is the weight of the ensemble member $m$, which can be
interpreted as a Bayesian posterior probability that $f_{m}$ is the correct
ensemble member.

As an example, the “tempering” approach of Hoffman and Kalnay [1983] can be
cast into the BMA framework by letting $p(y|f_{m})=\delta(y-f_{m})$ where
$\delta$ is Dirac’s delta function, and the $w_{m}$ are fit using least
squares. Under this model, the mean and variance of this distribution are

$E(Y|S_{j})=\sum_{m}w_{m}f_{m}$

and

$Var(Y|S_{j})=E((Y-E[Y|S_{j}])^{2})=\sum_{m}w_{m}\left(f_{m}-\sum_{m}w_{m}f_{m}\right)^{2}.$
(2)

We can check the consistency of this relationship empirically by further
averaging this equation in space and time, and comparing the left-hand-side
(error) to the right (spread). While we do not use this derivation in our
results it could be useful for future work.