# ChaosBench: A Multi-Channel, Physics-Based Benchmark for Subseasonal-to-
Seasonal Climate Prediction

Juan Nathaniel
Columbia University
<EMAIL_ADDRESS>Yongquan Qu
Columbia University
<EMAIL_ADDRESS>Tung Nguyen
UCLA
<EMAIL_ADDRESS>Sungduk Yu
UCI
<EMAIL_ADDRESS>Julius Busecke
LDEO, Columbia University
<EMAIL_ADDRESS>Aditya Grover
UCLA
<EMAIL_ADDRESS>Pierre Gentine
Columbia University
<EMAIL_ADDRESS>
###### Abstract

Accurate prediction of climate in the subseasonal-to-seasonal scale is crucial
for disaster readiness, reduced economic risk, and improved policy-making
amidst climate change. Yet, S2S prediction remains challenging due to the
chaotic nature of such system. At present, existing benchmarks for weather and
climate applications, tend to (1) have shorter forecasting range of up-to 14
days, (2) do not include a wide range of operational baseline forecasts, and
(3) lack physics-based constraints for explainability. Thus, we propose
ChaosBench, a large-scale, multi-channel, physics-based benchmark for S2S
prediction. ChaosBench has over 460K frames of real-world observations and
simulations, each with 60 variable-channels and spanning for up-to 45 years.
We also propose several physics-based, in addition to vision-based metrics,
that enables for a more physically-consistent model. Furthermore, we include a
diverse set of physics-based forecasts from 4 national weather agencies as
baselines to our data-driven counterpart. We establish two tasks that vary in
complexity: full and sparse dynamics prediction. Our benchmark is one of the
first to perform large-scale evaluation on existing models including
PanguWeather, FourCastNetV2, GraphCast, and ClimaX, and finds methods
originally developed for weather-scale applications fails on S2S task: they
perform much worse than just simply taking the long-term climatological
averages. We release our benchmark code and datasets at https://leap-
stc.github.io/ChaosBench.

_https://leap-stc.github.io/ChaosBench_

Figure 1: We propose ChaosBench, a large-scale, multi-channel, physics-based
benchmark for subseasonal-to-seasonal (S2S) climate prediction. It is framed
as a high-dimensional video regression task that consists of 45-year,
60-channel observations for validating physics-based and data-driven models,
and training the latter. Physics-based forecasts are generated from 4 national
weather agencies with 44-day lead-time and serve as baselines to data-driven
forecasts. Our benchmark is one of the first to incorporate physics-based
metrics to ensure physically-consistent and explainable models. We establish
two tasks: full and sparse dynamics prediction. The blurred image at $\Delta
t=44$ represents a common challenge of long-term forecasting: the decay of
high-frequency spectral signal.

##  1 Introduction

Datasets |  # input channels |  # target channels |  # benchmark agencies |  forecast lead (days) |  simulations |  observations |  physics-based metrics |  spatial extent |  temporal extent (input intersect)
---|---|---|---|---|---|---|---|---|---
WeatherBench [43] | 110 | 110 | 2 | 14 | ✓ | ✓ | ✗ | global | 1979-2018
SubseasonalRodeo [15] | 11 | 2 | 1 | 44 | ✓ | ✓ | ✗ | western US | 1979-2018
ClimateLab [56] | 20 | 2 | 2 | 44 | ✓ | ✓ | ✗ | global | 2000-2019
SubseasonalClimateUSA [32] | 25 | 2 | 2 | 44 | ✓ | ✓ | ✗ | contiguous US | 1981-2023
ChaosBench (ours) | 60 | 48 | 4 | 44 | ✓ | ✓ | ✓ | global | 1979-2023
Table 1: Comparison with other benchmark datasets: WeatherBench is focused on
short and medium-range forecasting task of up-to 14 days, SubseasonalRodeo and
SubseasonalClimateUSA are limited to the US, and ClimateLab is only evaluated
against 3 variable-channels. So far, ChaosBench (ours) is evaluated on the
largest set of variables (60 variable-channels), comprised of data with the
longest temporal extent (45 years), benchmarked against the largest number of
operational weather models (4 national agencies across the US, Europe, UK, and
Asia), and incorporates physics-based metrics for a more physically-faithful
and explainable data-driven forecast.

Although critical for economic planning, disaster preparedness, and policy-
making, subseasonal-to-seasonal (S2S) prediction has only recently begun to
gain attention, lagging behind the more established fields of short/medium-
range weather, or long-range climate predictions. For instance, many natural
hazards tend to manifest in the S2S scale, including the slow-onset of
droughts that lead to wildfire [38, 5], heavy precipitations that lead to
flooding [48], and persistent weather anomalies that lead to extremes [65].

So far, current approaches to weather and climate prediction are heavily
reliant on physics-based models in the form of Numerical Weather Prediction
(NWP). Many NWPs are based on the discretization of governing equations that
describe thermodynamics, fluid flows, _etc_. However, these models are
expensive to run especially in high-resolution setting. For example, there are
massive computational overheads to perform numerical integration at fine
spatiotemporal resolutions that are operationally-useful [47]. Furthermore,
their relative inaccessibility to non-experts is a major roadblock to the
broader community. As a result, there is a growing interest to apply data-
driven models to emulate NWPs, as they tend to have faster inference speed,
are less resource-hungry, and more accessible [4, 21, 33, 37, 40, 64].

Nevertheless, many data-driven benchmarks have so far been focused on the
short (1-5 days), medium (5-14 days), and long (months to years) forecasting
ranges. In this work, we include the S2S scale as a more challenging task:
being in between two extreme scales, the S2S system is doubly sensitive to (1)
_initial conditions_ as in the case for short/medium-range weather system, and
(2) _boundary conditions_ (i.e., coupled interaction between different
components of the Earth system, as illustrated in Figure 2) as in the case for
long-range climate system [39, 63, 24].

We propose ChaosBench (Figure 1) to bridge these gaps. We process over 45
years of observations for validating physics-based and data-driven models, and
training the latter. Each observation frame of size $121\times 240$ is
comprised of 60 variable-channels (e.g., temperature at 1000-hPa pressure
level). We also provide 44-day ahead physics-based forecasts as baselines,
generated from 4 national weather agencies over the last 8 years with 48
variable-channels. In addition, we introduce physics-based, in addition to
vision-based metrics, for a more physically-consistent and explainable model
[7, 35]. Also, we establish two tasks based on different research needs: full
and sparse dynamics prediction. As far as we know, ChaosBench is one of the
first to systematically evaluate several state-of-the-art data-driven models
including PanguWeather [3], FourCastNetV2 [37], and GraphCast [21] using both
vision-based and physics-based metrics on long-term S2S forecasting scale.

|  |  | Vertical Pressure Levels (hPa)
---|---|---|---
Variables | Notation | Unit | 1000 | 925 | 850 | 700 | 500 | 300 | 200 | 100 | 50 | 10
Geopotential height | $z$ | $gpm$ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
Specific humidity | $q$ | $kg\ kg^{-1}$ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | - | - | -
Temperature | $t$ | $K$ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
U component of wind | $u$ | $m\ s^{-1}$ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
V component of wind | $v$ | $m\ s^{-1}$ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
Vertical velocity | $w$ | $Pa\ s^{-1}$ | - | - | - | - | ✓ | - | - | - | - | -
Table 2: List of variables and their corresponding pressure levels (hPa)
available for physics-based model (48 variable-channels). The input
observations consist of all variable/level combinations (60 variable-
channels).

In this work, we demonstrate that existing physics-based and data-driven
models perform much worse than climatology as the forecasting range approaches
the S2S scale. The high spectral divergence observed in many state-of-the-art
models suggests the lost of predictive accuracy of small-scale structures.
This leads to significant blurring and a tendency towards smoother
predictions. For one, such averaging is of little use when one attempts to
identify extreme events requiring high-fidelity forecasts on the S2S scale
(e.g., regional droughts, hurricanes, _etc_). Also, performing worse than
climatology renders them _operationally unusable_. This highlights the urgent
need for a robust and unified data-driven S2S intercomparison project.

Figure 2: Generating forecasts from physics-based models that couple different
parts of the Earth system, including the atmosphere, cryosphere, ocean, and
terrestrial ecosystems. Then several operational choices are made, including
model initialization, perturbation of ensemble member, data assimilation to
correct for errors, numerical integration schemes, _etc_. The models then
iteratively generate simulations/forecast for $T$ timesteps.

##  2 Related Work

Building accurate, skillful, and explainable data-driven models for S2S
prediction is useful to reduce computational constraints and inaccessibility
of physics-based models to non-experts. Over recent years, multiple benchmarks
have been introduced to push the field of data-driven weather and climate
prediction [15, 18, 41, 34, 43, 56, 58]. We analyze the limitations of
existing works, and propose how ChaosBench bridges these gaps (Table 1,
Appendix A).

Gap in forecast lead-time. Many existing benchmarks are built for
short/medium-range weather (up to 14 days) [43, 18, 41], and long-term climate
(monthly to yearly scale) [58]. As discussed earlier, these problems tend to
be easier due to the lack of combined sensitivities to initial and boundary
conditions [39, 63].

Limited spatiotemporal extent. Many S2S benchmarks are focused on regional
forecasts, such as the US [15, 32], which ignore significant space-time
correlation in S2S-scale forecasting (e.g., teleconnection; butterfly effect
[24]). In addition, the temporal extent of observation with common interval is
more varied, with some less than 20 years [18, 56]. ChaosBench has the most
extensive overlapping temporal coverage yet, extending to 45 full years of
inputs.

Limited diversity of baseline models. Having a large set of physics-based
forecasts as baselines is key to reduce bias and diversify the target goal-
posts. Some previous benchmarks are mostly concerned with increasing the
number of data-driven models for baselines [43, 15]. In contrast, ChaosBench
also focuses on expanding the diversity of physics-based models, including
those developed by national weather agencies in the US, Europe, UK, and Asia.

Lack of physics-based constraints. So far, limited number of benchmarks have
explicitly incorporated physical principles to improve or constrain forecasts.
The second version of WeatherBench [44] attempts to incorporate elements of
physics, albeit primarily for analytical purposes. ChaosBench introduces
physics-based metrics that can be easily used for comparison (_scalar_) and
integrated into the training/inference pipeline (_differentiable_).

##  3 ChaosBench Dataset

(a) Normalized temperature@850-hpa label, prediction, and residual at the
first ($t=1$) and final ($t=44$) prediction step

(b) Power spectrum $S(k)$ _vs._ wavenumber $k$ plot as a function of
prediction step of normalized temperature @850-hpa

Figure 3: Motivating problem for long-range prediction: as we iteratively
unroll our prediction steps, (a) the residual error (truth - prediction)
grows, and the prediction noticeably becomes more blurry. This behavior is
captured in the Fourier frequency domain, (b) where the power spectra $S(k)$
at low wavenumber $k$ (i.e., low frequency signal) remains constant with
larger prediction step. However, $S(k)$ at higher $k$ (i.e., high frequency
signal) becomes convoluted and less structured. This phenomenon partly
explains why long-term forecasts tend to capture _large-scale pattern_ but
fail to incorporate _fine-grained detail_.

ChaosBench provides dataset and tasks that compare physics-based and data-
driven models for S2S prediction. The primary objective is to _fully_ or
_partially_ emulate physics-based models with data-driven ones, eliminating
computation bottleneck and inaccessibility of the former.

###  3.1 Input Observation for Data-Driven Models

We processed hourly ERA5 reanalysis data from 1979 to 2023 [13]. ERA5
reanalysis provides a comprehensive record of the global atmosphere which
combines physics and historical observations for correction. In order to
generate daily inputs, we select measurements at the 00UTC step. We also
perform re-gridding to 1.5-degree resolution to match that of physics-based
forecasts. The list of variables available for the inputs is provided in Table
2, including entries with the (✓) and (-) signs, totalling 60 channels.

###  3.2 Baseline Forecast from Physics-Based Models

We briefly describe the forecast generation process from physics-based models
(refer to Figure 2 for an illustration), including details on forecast
frequency and the number of ensemble members. More details are provided in
Appendix B. The list of available variables for physics-based forecast is
provided in Table 2, as indicated by the (✓) sign, for a total of 48 channels.
In all, we process forecasts between 2016 and 2023 from 4 national weather
agencies.

  * •

UKMO. The UK Meteorological Office uses the Global Seasonal forecast system
version 6 (GloSea6) model [60] to generate daily forecast for 60-day lead
time.

  * •

NCEP. The National Centers for Environmental Prediction uses the Climate
Forecast System version 2 (CFSv2) model [46] to generate daily forecast for
45-day lead time.

  * •

CMA. The China Meteorological Administration uses the Beijing Climate Center
(BCC) fully-coupled BCC-CSM2-HR model [62] to generate forecast at 3-day
interval for 60-day lead time.

  * •

ECMWF. The European Centre for Medium-Range Weather Forecasts uses the
operational Integrated Forecasting System (IFS) that includes advanced data
assimilation strategies and global numerical model of the Earth system [8]. In
particular, we use the CY41R1 version of the IFS which involves 51 ensemble
members and generates forecast twice every week for 46-day lead time.

###  3.3 Auxiliary Baselines

In addition to baseline forecasts from physics-based and data-driven models,
we provide additional auxiliary baselines. This includes climatology, the
long-term average and variability of weather-states, and persistence, which
uses current observation as next-step prediction.

##  4 Evaluation

(a) RMSE ($\downarrow$ is better)

(b) Bias ($\rightarrow 0$ is better)

(c) MS-SSIM ($\uparrow$ is better)

(d) SDiv ($\downarrow$ is better)

Figure 4: Evaluation results between baseline climatology (black lines) and
physics-based models. Overall, we observe that most models have similar metric
trajectory. Furthermore, at longer forecasting horizon, most physics-based
models collapse away from climatology while maintaining fine-grained
structures as evidenced from their low SpecDiv (barring NCEP). This suggests
knowledge gap in the modeling process, such as poorly represented sub-grid
phenomena.

This section describes the evaluation metrics to standardize the benchmarking
process for different models and tasks. In general, we divide our criteria
into: (1) vision-based, and (2) physics-based metrics. For each metric, unless
otherwise specified, we apply a weighting scheme at each latitude $\theta_{i}$
as defined by Equation 1.

$w(\theta_{i})=\frac{cos(\theta_{i})}{\frac{1}{|\boldsymbol{\theta}|}\sum_{a=1}^{|\boldsymbol{\theta}|}cos(\theta_{a})}$
(1)

where $\boldsymbol{\theta}$ is the set of all latitudes in our data, and
$|\boldsymbol{\theta}|$ is its cardinality. The latitude weighting is
necessary since a grid cell in the equator covers larger areas than that at
the poles because of projection. We denote the input at time $t$ as
$\mathbf{X}_{t}\in\mathbb{R}^{h\times w\times p\times l}$, where $h,w,p,l$
represent the height (i.e., latitude), width (i.e., longitude), parameter
(e.g., temperature), and its associated vertical level (e.g., 1000-hpa). In
addition, we denote
$\\{\mathbf{Y}_{t},\mathbf{\hat{Y}}_{t}\\}\in\mathbb{R}^{h\times w\times
p\times l}$ as the ground-truth label and prediction respectively. Finally, we
denote each element of latitude and longitude as
$\theta_{i}\in\boldsymbol{\theta}$ and $\gamma_{j}\in\boldsymbol{\gamma}$.

###  4.1 Vision-Based Metrics

We first propose several classic vision-based metrics, including RMSE, Bias,
ACC, and MS-SSIM.

Root Mean Squared Error (RMSE) is useful to penalize outliers, which are
especially critical for weather and climate applications such as extreme event
prediction. We apply a latitude-adjusted RMSE as defined in Equation S1.

Bias assists us to identify misspecification and any clues of systematic
errors present in the model. We apply a latitude-adjusted Bias as defined in
Equation S2.

Anomaly Correlation Coefficient (ACC) measures the correlation between
predicted and observed anomalies. This metric is especially useful in weather
and climate applications, where deviations from the norm (e.g., temperature
anomalies) often reveal interesting insights. We apply a latitude-adjusted ACC
as defined in Equation S3.

Multi-Scale Structural Similarity (MS-SSIM) [57] compares structural
similarity between forecast and ground-truth label across scales (refer to
Appendix D.4 for more details). This is especially useful in weather systems
because they occur at multiple scales, from large systems like cyclones, to
smaller features like localized rain thunderstorms.

(a) RMSE ($\downarrow$ is better)

(b) Bias ($\rightarrow 0$ is better)

(c) MS-SSIM ($\uparrow$ is better)

(d) SDiv ($\downarrow$ is better)

Figure 5: Evaluation results between baseline climatology (black line) and
data-driven models (PanguWeather, FourCastNetV2, GraphCast). Overall, we
observe that state-of-the-art models perform significantly worse than
climatology on S2S timescale. They also perform poorly on physics-based
metrics indicating the lack of predictive power on small-scale structures.

###  4.2 Physics-Based Metrics

Our motivation of incorporating physics-based metrics is illustrated in Figure
3. We notice that at longer lead-time, $\Delta t$, not only does the residual
error grow, but our forecast becomes more blurry (Figure 3(a)). This
phenomenon is pervasive and has been widely reported by many working on data-
driven weather and climate forecasting [36, 37, 4]. Several root causes
include the tendency of RMSE toward the mean and power decay in spectral
domain.

When we transform our predictions into the spatial Fourier frequency domain
(details in Appendix E), Figure 3(b) highlights the relatively uniform power
spectra $S(k)$ at low wavenumber $k$ (i.e., low-frequency signal) even with
increasing $\Delta t$. On the other hand, $S(k)$ at higher $k$ (i.e., high-
frequency signal) becomes less structured across $\Delta t$. This phenomenon
partly explains why long-term forecasts tend to capture large-scale pattern
but fail to incorporate fine-grained details, which is reflected by the
qualitative blurriness [44]. In order to improve fidelity, we need to preserve
high-frequency signals to be as close as observations.

This motivates us to propose two physics-based metrics that measure the
deviation and difference between the power spectra of prediction
${\hat{S}}(k)$ and observation $S(k)$, where $k\in\mathbf{K}$, and
$\mathbf{K}$ is the set of all scalar wavenumbers from 2D Fourier transform.
Focusing on high-frequency components, we introduce
$\mathbf{K}_{q}=\\{k\in\mathbf{K}\mid k\geq Q(q)\\}$, where $Q$ is the
quantile function of $\mathbf{K}$ and $q\in[0,1]$. We set $q=0$ or $q=0.9$ for
training and evaluation respectively. However, our implementation allows for
variable $q$. We denote $S_{q}=\\{S(k)\mid k\in\mathbf{K}_{q}\\}$ as the
corresponding power spectra on $\mathbf{K}_{q}$, and we normalize the
distribution to $S^{\prime}(k)$ such that it sums up to 1 and can be likened
to a probability density function. Similarly we use $\hat{S}^{\prime}(k)$ to
denote the normalized power for predictions.

Spectral Divergence (SpecDiv) follows principles from Kullback–Leibler (KL)
divergence [20] where we compute the expectation of the log ratio between
target $S^{\prime}(k)$ and prediction $\hat{S}^{\prime}(k)$ spectra, and is
defined in Equation 2 (see Listing LABEL:lst:specdiv for PyTorch psuedocode).

$\mathcal{M}_{SpecDiv}=\sum_{k}S^{\prime}(k)\cdot\log(S^{\prime}(k)/\hat{S}^{\prime}(k))$
(2)

Spectral Residual (SpecRes) follows principles from RMSE and adapted from [52]
where we compute the root of the expected squared residual, and is defined in
Equation 3 (see Listing LABEL:lst:specres for PyTorch psuedocode).

$\mathcal{M}_{SpecRes}=\sqrt{\mathbb{E}_{k}[(\hat{{S}}^{\prime}(k)-{S}^{\prime}(k))^{2}]}$
(3)

The expectations in Equations 2 and 3 are calculated over $\mathbf{K}_{q}$.
For both physics-based metrics, the value will be zero if the distribution of
high-frequency power of the forecast is identical to that of observations, but
will increase as the former loses more of it. Essentially, both metrics
measure how well the forecasts _preserve_ high-frequency signals relative to
observations.

##  5 Task & Experimental Setup

We propose two tasks: full and sparse dynamics prediction. Both tasks simulate
the full states: $\mathbf{\hat{Y}}_{t+1}=f_{\phi}(\mathbf{X}_{t})$, but the
optimization is performed either over _all_ or a _subset_ of variables, and
$f_{\phi}$ is data-driven model parameterized by $\phi$.

###  5.1 Task 1: Full Dynamics Prediction

In many cases, we seek to train data-driven models so that they are able to
perform well on _all_ atmospheric states by optimizing for the full state (60
variable-channels) for the next forecasting timestep $t+1$, that is,

$\phi^{\ast}=\operatorname*{argmin}_{\phi}\mathcal{L}(\mathbf{\hat{Y}}_{t+1},\mathbf{Y}_{t+1})$

where $\mathcal{L}$ is any loss function. This task is especially useful for
building climate emulators that act as surrogates for the more expensive
physics-based NWP models [36].

###  5.2 Task 2: Sparse Dynamics Prediction

Although the first task is useful for learning the full complex interaction
between variables, it is relatively difficult due to the intrinsic high-
dimensionality of the data. As a result, we introduce a second task that
allows for the optimization on a subset of variables of interest
($\mathbf{Y^{\prime}}\in\mathbf{Y}$):

$\phi^{\ast}=\operatorname*{argmin}_{\phi}\mathcal{L}(\mathbf{\hat{Y}^{\prime}}_{t+1},\mathbf{Y^{\prime}}_{t+1})$

###  5.3 Data-Driven Baseline Models

We use 8 baseline models for comparison, including Lagged Autoencoder [27],
ResNet [42], UNet [43], and state-of-the-art Fourier Neural Operator (FNO)
[37], ViT-ClimaX [36], PanguWeather [3], FourCastNetV2 [37], and GraphCast
[21]. The full details on the architectures and final hyperparameter choices
are discussed in Appendix C.

(a) RMSE ($\downarrow$ is better)

(b) Bias ($\rightarrow 0$ is better)

(c) MS-SSIM ($\uparrow$ is better)

(d) SDiv ($\downarrow$ is better)

Figure 6: Evaluation results between baseline climatology (black lines) and
data-driven models trained using either the direct (ClimaX) or the
autoregressive approach (PanguWeather, FourCastNetV2, GraphCast). The y-axis
represent the metric scores and the x-axis the number of forecasting lead-
time, $\Delta t$. Overall, our findings indicate that there is potential for
improvement in S2S forecasting as evidenced by improvements in models trained
using the direct approach, which serves as an upper-bound to models trained
using the more difficult, yet realistic autoregressive approach.

###  5.4 Experimental Setup

For both tasks described in 5.1 and 5.2, we report results for
$\hat{\mathbf{X}}\in\\{T850,Z500,Q700\\}$, following Weatherbench 2 [44]. The
full benchmark scores are available at https://leap-stc.github.io/ChaosBench.

We select the years 1979-2015 for training, 2016-2021 for validation, and
2022-2023 for testing. All benchmark scores for physics-based and data-driven
models are computed for the year 2022-2023. The reason for not using the
entire 8 years of available simulations for testing (i.e., 2016-2023) is
because we allow for the possibilities to combine physics-based forecasts with
actual observations to build an increasingly popular form of hybrid physics-
informed models [17, 49].

##  6 Results and Ablation Studies

We first evaluate baseline models for S2S forecasting. We then sketch several
challenges and plausible ways to improve S2S prediction, including the
different strategies to incorporate temporal context and the interplay between
physics and machine learning.

###  6.1 Collapse in Predictive Skills

Referring to Figure 4, we find that the ensemble forecast across all
operational centers exhibit similar performance, barring NCEP for physics-
based SpecDiv metric. One notable observation is how these models tend to
produce forecasts worse than climatology at S2S-scale beyond 14 days ahead.
Regardless, physics-based forecasts have comparatively low spectral divergence
when compared to data-driven models (discussed next). This highlights the dual
problems of error accumulation in autoregressive forecasting steps and the
inherent knowledge gaps in the modeling process [39, 63, 24]. Hence, the
phenomenon of skill collapse further emphasizes the notoriously difficult
problem of S2S forecasting.

Similar phenomenon of skill collapse is also observed across state-of-the-art
models as illustrated in Figure 5. However, different from their physics-based
counterparts, these forecasts have significantly higher spectral divergence.
This indicates the lack of predictive skills of small-scale structures at
longer forecasting horizon, resulted from the blurring or smoothing artifacts
showcased throughout this and previous works [44].

Overall, we highlight the different problems apparent in both physics-based
and data-driven models in S2S forecasting tasks. Nonetheless, since data-
driven models are often trained directly on model outputs + observations –
which to some extent captures true physical knowledge – improving fidelity by
preserving high-frequency signals over long-term iterative steps could be a
way to narrow the gap to robust S2S forecasting.

###  6.2 Effect of Different Training Strategies

The next question of interest would be the effect of different training
strategies commonly adopted in the weather and climate forecasting literature.
We consider two increasingly popular temporal-aware training schemes:
autoregressive and direct approaches [36, 43]

(a) ClimaX task 1 using direct approach

(b) ClimaX task 2 using direct approach

Figure 7: Normalized humidity@700-hpa qualitative results for direct approach
of ClimaX. We observe that despite performing well on vision-based metrics at
S2S scale ($\Delta t=44$), the forecast is blurred, implying less than optimal
performance on physics-based metrics, such as SpecDiv.

  * •

Autoregressive approach uses the output of current prediction as the input to
forecast future timestep, and computes the average iterative loss as in
Equation 4. Here $S=\\{1,\cdots,s\\}$ and $s$ is the autoregressive steps. For
this work, we set $s=5$.

$\begin{cases}\text{Task
1:}\quad\frac{1}{|S|}\sum_{i=1}^{s}\mathcal{L}(\mathbf{\hat{Y}_{t+s_{i}}}\mathbf{Y_{t+s_{i}}}),\forall
s_{i}\in S\\\ \text{Task
2:}\quad\frac{1}{|S|}\sum_{i=1}^{s}\mathcal{L}(\mathbf{\hat{Y}^{\prime}_{t+s_{i}}}\mathbf{Y^{\prime}_{t+s_{i}}}),\forall
s_{i}\in S\end{cases}$ (4)

  * •

Direct approach allows us to train models at arbitrary lead-time, $\Delta t$.
For this work, we set $\Delta t\in\\{1,5,10,15,20,25,30,35,40,44\\}$.

Autoregressive approach is the more popular of the two since it is similar to
how physics-based dynamical models generate forecasts. It is widely adopted in
existing models including PanguWeather, FourCastNetV2, and GraphCast.
Nonetheless, a direct approach provides an upper bound to the skillfulness of
an autoregressive model since there is trivial extent of error propagation and
accumulation.

Upper-Bound to Long-Term Forecasting Skill. As observed in Figure 6, direct
approach is better than models trained autoregressively. This solidifies our
hypothesis that a direct approach provides an upper-bound to long-term
forecasting skills as the latter will accumulate and propagate error as
$\Delta t\rightarrow T$. However, a direct approach is more costly since it
requires tuning at each $\Delta t$. On the other hand, autoregressive approach
is closer to how dynamical systems evolve. Even though it is more challenging,
a robust autoregressive model signifies a promising step towards emulating the
right physics.

Disentangling Error Types. Direct approach also distinguishes between _error
types_. For instance, the performance gaps between direct and autoregressive
methods can be attributed to (1) _iterative error propagation_ , while the
rest points to (2) _model deficiencies_. By systematically targeting different
types of errors, we can accelerate progress toward robust S2S forecasts.

###  6.3 Toward Physically-Consistent S2S Prediction

Despite the promise of data-driven models for weather and climate prediction,
especially in the short/medium range, the task remains challenging for the S2S
case.

(a) $\hat{S}(k)$ of ClimaX-direct for Task 1

(b) $\hat{S}(k)$ of ClimaX-direct for Task 2

Figure 8: $\hat{S}(k)$ of ClimaX-direct across tasks where blue and orange
lines indicate forecasts at $\Delta t=1$ and $\Delta t=44$ respectively. There
is an overall decay in $\hat{S}(k)$ (y-axis), especially at high-frequency $k$
(x-axis), as $\Delta t$ increases. Note the logarithmic scale on the y-axis
i.e., gaps in orders of magnitude.

Challenge: Consider again Figure 5 which benchmarks results from existing
state-of-the-art models. Although they are able to perform comparably with
baseline climatology and physics-based models in terms of vision-based metrics
(e.g., RMSE, Bias, MS-SSIM), they still massively underperform in physics-
based metrics (e.g., SpecDiv), even at $\Delta t\geq 1$. This is again
highlighted visually in Figure 7 where we notice substantial blurring, and in
Figure 8 where we observe power decay $\hat{S}(k)$ at $\Delta t=44$ of ClimaX
as one of the benchmark models.

##  7 Conclusion

We introduce ChaosBench, a large-scale, multi-channel, physics-based benchmark
dataset for S2S prediction. ChaosBench consists of 45 years of input
observations with 16K snapshopts of 60 variable-channels, and 8 years of
physics-based forecasts from 4 national weather agencies with 440K frames of
48 variable-channels. We establish two challenges that vary in complexity and
use cases: full and sparse dynamics prediction. Through various experiments,
we demonstrate that existing physics-based and data-driven models perform
significantly worse than climatology as the forecasting range approaches the
S2S scale. The high spectral divergence observed in many state-of-the-art
models suggests the lost of predictive accuracy of small-scale structures.
This leads to significant blurring and a tendency towards smoother forecasts.
Consequently, there is a clear need to integrate physics for improved and
interpretable S2S forecasts. Furthermore, our work has shown the potential to
improve S2S forecasts through a direct training method, which serves as an
upper-bound to the more challenging, yet inherently more realistic,
autoregressive approach. Therefore, we invite the computer vision, and the
broader machine learning community, to take part in solving one of the more
challenging and pressing scientific problems of our time: _predicting chaos_.

## Acknowledgements

The authors would like to acknowledge funding from the NSF LEAP Science and
Technology Center award #2019625, Department of Energy grant. We also
acknowledge the computing and storage resources provided by the NSF Science
and Technology Center (STC) Learning the Earth with Artificial intelligence
and Physics (LEAP) (Award #2019625). AG would like to acknowledge support from
Google and Schmidt Sciences.

## References

  * Arakawa [1997] Akio Arakawa.  Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow. part i.  _Journal of computational physics_ , 135(2):103–114, 1997.
  * Best et al. [2011] Martin J Best, M Pryor, DB Clark, Gabriel G Rooney, R Essery, CB Ménard, JM Edwards, MA Hendry, A Porson, N Gedney, et al.  The joint uk land environment simulator (jules), model description–part 1: energy and water fluxes.  _Geoscientific Model Development_ , 4(3):677–699, 2011.
  * Bi et al. [2022] Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian.  Pangu-weather: A 3d high-resolution model for fast and accurate global weather forecast.  _arXiv preprint arXiv:2211.02556_ , 2022.
  * 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.
  * Buch et al. [2023] Jatan Buch, A Park Williams, Caroline S Juang, Winslow D Hansen, and Pierre Gentine.  Smlfire1. 0: a stochastic machine learning (sml) model for wildfire activity in the western united states.  _Geoscientific Model Development_ , 16(12):3407–3433, 2023.
  * Chen et al. [2022] Boyuan Chen, Kuang Huang, Sunand Raghupathi, Ishaan Chandratreya, Qiang Du, and Hod Lipson.  Automated discovery of fundamental variables hidden in experimental data.  _Nature Computational Science_ , 2(7):433–442, 2022.
  * da Silva et al. [2022] Ademir Ferreira da Silva, Juan Nathaniel, Ken CL Wong, Campbell Watson, Hongzhi Wang, Jitendra Singh, Alexandre Alkmim Chamon, and Levente Klein.  Netzeroco 2, an ai framework for accelerated nature-based carbon sequestration.  In _2022 IEEE International Conference on Big Data (Big Data)_ , pages 4881–4887. IEEE, 2022.
  * [8] IFS DOCUMENTATION.  Part v: The ensemble prediction.
  * Dosovitskiy et al. [2020] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  _arXiv preprint arXiv:2010.11929_ , 2020.
  * Ek et al. [2003] MB Ek, KE Mitchell, Ying Lin, Eric Rogers, Pablo Grunmann, Victor Koren, George Gayno, and JD Tarpley.  Implementation of noah land surface model advances in the national centers for environmental prediction operational mesoscale eta model.  _Journal of Geophysical Research: Atmospheres_ , 108(D22), 2003.
  * Griffies et al. [2004] Stephen M Griffies, Matthew J Harrison, Ronald C Pacanowski, Anthony Rosati, et al.  A technical guide to mom4.  _GFDL Ocean Group Tech. Rep_ , 5(5):371, 2004.
  * Hendrycks and Gimpel [2016] Dan Hendrycks and Kevin Gimpel.  Gaussian error linear units (gelus).  _arXiv preprint arXiv:1606.08415_ , 2016.
  * 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.
  * Hortal and Simmons [1991] Mariano Hortal and AJ Simmons.  Use of reduced gaussian grids in spectral models.  _Monthly Weather Review_ , 119(4):1057–1074, 1991.
  * Hwang et al. [2019] Jessica Hwang, Paulo Orenstein, Judah Cohen, Karl Pfeiffer, and Lester Mackey.  Improving subseasonal forecasting in the western us with machine learning.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pages 2325–2335, 2019.
  * Janssen et al. [2005] Peter Janssen, Jean-Raymond Bidlot, Saleh Abdalla, and Hans Hersbach.  _Progress in ocean wave forecasting at ECMWF_.  ECMWF Reading, UK, 2005.
  * Karniadakis et al. [2021] George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang.  Physics-informed machine learning.  _Nature Reviews Physics_ , 3(6):422–440, 2021.
  * Kashinath et al. [2021] Karthik Kashinath, Mayur Mudigonda, Sol Kim, Lukas Kapp-Schwoerer, Andre Graubner, Ege Karaismailoglu, Leo Von Kleist, Thorsten Kurth, Annette Greiner, Ankur Mahesh, et al.  Climatenet: An expert-labeled open dataset and deep learning architecture for enabling high-precision analyses of extreme weather.  _Geoscientific Model Development_ , 14(1):107–124, 2021.
  * Kobayashi et al. [2015] Shinya Kobayashi, Yukinari Ota, Yayoi Harada, Ayataka Ebita, Masami Moriya, Hirokatsu Onoda, Kazutoshi Onogi, Hirotaka Kamahori, Chiaki Kobayashi, Hirokazu Endo, et al.  The jra-55 reanalysis: General specifications and basic characteristics.  _Journal of the Meteorological Society of Japan. Ser. II_ , 93(1):5–48, 2015.
  * Kullback and Leibler [1951] Solomon Kullback and Richard A Leibler.  On information and sufficiency.  _The annals of mathematical statistics_ , 22(1):79–86, 1951.
  * Lam et al. [2022] Remi Lam, Alvaro Sanchez-Gonzalez, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Alexander Pritzel, Suman Ravuri, Timo Ewalds, Ferran Alet, Zach Eaton-Rosen, et al.  Graphcast: Learning skillful medium-range global weather forecasting.  _arXiv preprint arXiv:2212.12794_ , 2022.
  * Lawrence and Chase [2007] Peter J Lawrence and Thomas N Chase.  Representing a new modis consistent land surface in the community land model (clm 3.0).  _Journal of Geophysical Research: Biogeosciences_ , 112(G1), 2007.
  * Li et al. [2020] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Fourier neural operator for parametric partial differential equations.  _arXiv preprint arXiv:2010.08895_ , 2020.
  * Lorenz [1963] Edward N Lorenz.  Deterministic nonperiodic flow.  _Journal of atmospheric sciences_ , 20(2):130–141, 1963.
  * Loshchilov and Hutter [2017] Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization.  _arXiv preprint arXiv:1711.05101_ , 2017.
  * Loveland et al. [2000] Thomas R Loveland, Bradley C Reed, Jesslyn F Brown, Donald O Ohlen, Zhiliang Zhu, LWMJ Yang, and James W Merchant.  Development of a global land cover characteristics database and igbp discover from 1 km avhrr data.  _International journal of remote sensing_ , 21(6-7):1303–1330, 2000.
  * Lusch et al. [2018] Bethany Lusch, J Nathan Kutz, and Steven L Brunton.  Deep learning for universal linear embeddings of nonlinear dynamics.  _Nature communications_ , 9(1):4950, 2018.
  * Malardel et al. [2016] Sylvie Malardel, Nils Wedi, Willem Deconinck, Michail Diamantakis, Christian Kühnlein, George Mozdzynski, Mats Hamrud, and Piotr Smolarkiewicz.  A new grid for the ifs.  _ECMWF newsletter_ , 146(23-28):321, 2016.
  * Mathiot et al. [2017] Pierre Mathiot, Adrian Jenkins, Christopher Harris, and Gurvan Madec.  Explicit representation and parametrised impacts of under ice shelf seas in the z coordinate ocean model nemo 3.6.  _Geoscientific Model Development_ , 10(7):2849–2874, 2017.
  * Meng et al. [2012] Jesse Meng, Rongqian Yang, Helin Wei, Michael Ek, George Gayno, Pingping Xie, and Kenneth Mitchell.  The land surface analysis in the ncep climate forecast system reanalysis.  _Journal of Hydrometeorology_ , 13(5):1621–1630, 2012.
  * Mogensen et al. [2012] Kristian Mogensen, Magdalena Alonso Balmaseda, Anthony Weaver, et al.  The nemovar ocean data assimilation system as implemented in the ecmwf ocean analysis for system 4.  2012\.
  * Mouatadid et al. [2023] Soukayna Mouatadid, Paulo Orenstein, Genevieve Elaine Flaspohler, Miruna Oprescu, Judah Cohen, Franklyn Wang, Sean Edward Knight, Maria Geogdzhayeva, Samuel James Levang, Ernest Fraenkel, et al.  Subseasonalclimateusa: A dataset for subseasonal forecasting and benchmarking.  In _Thirty-seventh Conference on Neural Information Processing Systems Datasets and Benchmarks Track_ , 2023.
  * Mukkavilli et al. [2023] S Karthik Mukkavilli, Daniel Salles Civitarese, Johannes Schmude, Johannes Jakubik, Anne Jones, Nam Nguyen, Christopher Phillips, Sujit Roy, Shraddha Singh, Campbell Watson, et al.  Ai foundation models for weather and climate: Applications, design, and implementation.  _arXiv preprint arXiv:2309.10808_ , 2023.
  * Nathaniel et al. [2023a] Juan Nathaniel, Jiangong Liu, and Pierre Gentine.  Metaflux: Meta-learning global carbon fluxes from sparse spatiotemporal observations.  _Scientific Data_ , 10(1):440, 2023a.
  * Nathaniel et al. [2023b] Juan Nathaniel, Gabrielle Nyirjesy, Campbell D Watson, Conrad M Albrecht, and Levente J Klein.  Above ground carbon biomass estimate with physics-informed deep network.  In _IGARSS 2023-2023 IEEE International Geoscience and Remote Sensing Symposium_ , pages 1297–1300. IEEE, 2023b.
  * 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.
  * 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.
  * Pendergrass et al. [2020] Angeline G Pendergrass, Gerald A Meehl, Roger Pulwarty, Mike Hobbins, Andrew Hoell, Amir AghaKouchak, Céline JW Bonfils, Ailie JE Gallant, Martin Hoerling, David Hoffmann, et al.  Flash droughts present a new challenge for subseasonal-to-seasonal prediction.  _Nature Climate Change_ , 10(3):191–199, 2020.
  * Privé and Errico [2013] Nikki C Privé and Ronald M Errico.  The role of model and initial condition error in numerical weather forecasting investigated with an observing system simulation experiment.  _Tellus A: Dynamic Meteorology and Oceanography_ , 65(1):21740, 2013.
  * Qu and Shi [2023] Yongquan Qu and Xiaoming Shi.  Can a machine learning–enabled numerical model help extend effective forecast range through consistently trained subgrid-scale models?  _Artificial Intelligence for the Earth Systems_ , 2(1):e220050, 2023.
  * Racah et al. [2017] Evan Racah, Christopher Beckham, Tegan Maharaj, Samira Ebrahimi Kahou, Mr Prabhat, and Chris Pal.  Extremeweather: A large-scale climate dataset for semi-supervised detection, localization, and understanding of extreme weather events.  _Advances in neural information processing systems_ , 30, 2017.
  * 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. [2020] Stephan Rasp, Peter D Dueben, Sebastian Scher, Jonathan A Weyn, Soukayna Mouatadid, and Nils Thuerey.  Weatherbench: a benchmark data set for data-driven weather forecasting.  _Journal of Advances in Modeling Earth Systems_ , 12(11):e2020MS002203, 2020.
  * Rasp et al. [2023] 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_ , 2023.
  * Saha et al. [2010] Suranjana Saha, Shrinivas Moorthi, Hua-Lu Pan, Xingren Wu, Jie Wang, Sudhir Nadiga, Patrick Tripp, Robert Kistler, John Woollen, David Behringer, et al.  Ncep climate forecast system reanalysis (cfsr) monthly products, january 1979 to december 2010.  2010\.
  * Saha et al. [2014] Suranjana Saha, Shrinivas Moorthi, Xingren Wu, Jiande Wang, Sudhir Nadiga, Patrick Tripp, David Behringer, Yu-Tai Hou, Hui-ya Chuang, Mark Iredell, et al.  The ncep climate forecast system version 2.  _Journal of climate_ , 27(6):2185–2208, 2014.
  * Schneider et al. [2023] Tapio Schneider, Swadhin Behera, Giulio Boccaletti, Clara Deser, Kerry Emanuel, Raffaele Ferrari, L Ruby Leung, Ning Lin, Thomas Müller, Antonio Navarra, et al.  Harnessing ai and computing to advance climate modelling and prediction.  _Nature Climate Change_ , 13(9):887–889, 2023.
  * Shamekh et al. [2023] Sara Shamekh, Kara D Lamb, Yu Huang, and Pierre Gentine.  Implicit learning of convective organization explains precipitation stochasticity.  _Proceedings of the National Academy of Sciences_ , 120(20):e2216158120, 2023.
  * Shen et al. [2023] Chaopeng Shen, Alison P Appling, Pierre Gentine, Toshiyuki Bandai, Hoshin Gupta, Alexandre Tartakovsky, Marco Baity-Jesi, Fabrizio Fenicia, Daniel Kifer, Li Li, et al.  Differentiable modelling to unify machine learning and physical models for geosciences.  _Nature Reviews Earth & Environment_, 4(8):552–567, 2023.
  * Sønderby et al. [2020] Casper Kaae Sønderby, Lasse Espeholt, Jonathan Heek, Mostafa Dehghani, Avital Oliver, Tim Salimans, Shreya Agrawal, Jason Hickey, and Nal Kalchbrenner.  Metnet: A neural weather model for precipitation forecasting.  _arXiv preprint arXiv:2003.12140_ , 2020.
  * Storkey et al. [2018] David Storkey, Adam T Blaker, Pierre Mathiot, Alex Megann, Yevgeny Aksenov, Edward W Blockley, Daley Calvert, Tim Graham, Helene T Hewitt, Patrick Hyder, et al.  Uk global ocean go6 and go7: A traceable hierarchy of model resolutions.  _Geoscientific Model Development_ , 11(8):3187–3213, 2018.
  * Takamoto et al. [2022] Makoto Takamoto, Timothy Praditia, Raphael Leiteritz, Daniel MacKinlay, Francesco Alesiani, Dirk Pflüger, and Mathias Niepert.  Pdebench: An extensive benchmark for scientific machine learning.  _Advances in Neural Information Processing Systems_ , 35:1596–1611, 2022.
  * Tian et al. [2002] Yuhong Tian, Curtis E Woodcock, Yujie Wang, Jeff L Privette, Nikolay V Shabanov, Liming Zhou, Yu Zhang, Wolfgang Buermann, Jiarui Dong, Brita Veikkanen, et al.  Multiscale analysis and validation of the modis lai product: I. uncertainty assessment.  _Remote Sensing of Environment_ , 83(3):414–430, 2002.
  * Tsujino et al. [2020] Hiroyuki Tsujino, L Shogo Urakawa, Stephen M Griffies, Gokhan Danabasoglu, Alistair J Adcroft, Arthur E Amaral, Thomas Arsouze, Mats Bentsen, Raffaele Bernardello, Claus W Böning, et al.  Evaluation of global ocean–sea-ice model simulations based on the experimental protocols of the ocean model intercomparison project phase 2 (omip-2).  _Geoscientific Model Development_ , 13(8):3643–3708, 2020.
  * Van Genuchten [1980] M Th Van Genuchten.  A closed-form equation for predicting the hydraulic conductivity of unsaturated soils.  _Soil science society of America journal_ , 44(5):892–898, 1980.
  * Vitart et al. [2022] Frederic Vitart, Andrew W Robertson, Aaron Spring, Florian Pinault, Rok Roškar, W Cao, S Bech, A Bienkowski, N Caltabiano, E De Coning, et al.  Outcomes of the wmo prize challenge to improve subseasonal to seasonal predictions using artificial intelligence.  _Bulletin of the American Meteorological Society_ , 103(12):E2878–E2886, 2022.
  * Wang et al. [2003] Zhou Wang, Eero P Simoncelli, and Alan C Bovik.  Multiscale structural similarity for image quality assessment.  In _The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers, 2003_, pages 1398–1402. Ieee, 2003.
  * Watson-Parris et al. [2022] Duncan Watson-Parris, Yuhan Rao, Dirk Olivié, Øyvind Seland, Peer Nowack, Gustau Camps-Valls, Philip Stier, Shahine Bouabid, Maura Dewey, Emilie Fons, et al.  Climatebench v1. 0: A benchmark for data-driven climate projections.  _Journal of Advances in Modeling Earth Systems_ , 14(10):e2021MS002954, 2022.
  * Wieder et al. [2014] WR Wieder, J Boehnert, GB Bonan, and M Langseth.  Regridded harmonized world soil database v1. 2.  _ORNL DAAC_ , 2014.
  * Williams et al. [2015] KD Williams, CM Harris, A Bodas-Salcedo, J Camp, RE Comer, D Copsey, D Fereday, T Graham, R Hill, T Hinton, et al.  The met office global coupled model 2.0 (gc2) configuration.  _Geoscientific Model Development_ , 88(55):1509–1524, 2015.
  * Wu et al. [2014] Tongwen Wu, Lianchun Song, Weiping Li, Zaizhi Wang, Hua Zhang, Xiaoge Xin, Yanwu Zhang, Li Zhang, Jianglong Li, Fanghua Wu, et al.  An overview of bcc climate system model development and application for climate change studies.  _Journal of Meteorological Research_ , 28:34–56, 2014.
  * Wu et al. [2019] Tongwen Wu, Yixiong Lu, Yongjie Fang, Xiaoge Xin, Laurent Li, Weiping Li, Weihua Jie, Jie Zhang, Yiming Liu, Li Zhang, et al.  The beijing climate center climate system model (bcc-csm): The main progress from cmip5 to cmip6.  _Geoscientific Model Development_ , 12(4):1573–1600, 2019.
  * Wu et al. [2005] Wanli Wu, Amanda H Lynch, and Aaron Rivers.  Estimating the uncertainty in a regional climate model related to initial and lateral boundary conditions.  _Journal of climate_ , 18(7):917–933, 2005.
  * Yu et al. [2023] Sungduk Yu, Walter M Hannah, Liran Peng, Mohamed Aziz Bhouri, Ritwik Gupta, Jerry Lin, Björn Lütjens, Justus C Will, Tom Beucler, Bryce E Harrop, et al.  Climsim: An open large-scale dataset for training high-resolution physics emulators in hybrid multi-scale climate simulators.  _arXiv preprint arXiv:2306.08754_ , 2023.
  * Zeppetello et al. [2022] Lucas R Vargas Zeppetello, David S Battisti, and Marcia B Baker.  The physics of heat waves: What causes extremely high summertime temperatures?  _Journal of Climate_ , 35(7):2231–2251, 2022.
  * Zobler [1986] Leonard Zobler.  A world soil file grobal climate modeling.  _NASA Tech. memo_ , 32, 1986.

ChaosBench: A Multi-Channel, Physics-Based Benchmark for Subseasonal-to-
Seasonal Climate Prediction
Supplementary Material

##  A Related Work

Here we discuss the criteria used to compare different S2S benchmark. This
list is by no means exhaustive and there exists many ways to interpret the
different contribution, strength, and scope of each. We refer interested
reader to the respective benchmark paper and website.

On Input Channels. The number of input channels indicates the number of
_unique_ variables used for training data-driven models. For instance, in the
case of SubseasonalClimateUSA, these include tmin, tmax, tmean, precip_agg,
precip_mean, SST, SIC, z-10, z100, z500, z850, u-250, u-925, v-250, v-925,
surface_P, RH, SSP, precipitable water, PE, DEM, KG, MJO-phase, MJO-amp,
ENSO-I, despite them having similar variable across data sources.

On Target Channels and Agencies. Similarly, the number of target channels
represent the variables these benchmarks are aiming for. This is closely
related to the number of benchmark agencies, which refers to the number of
physics-based simulations used as target, rather than inputs. In the case for
SubseasonalClimateUSA, for instance, the number of target channels correspond
to two: precipitation and surface temperature, while the number of benchmark
agencies is also two: CFSv2 and ECMWF, despite them using multiple other
simulations generated from agencies but as inputs. This is important because
if our aim is to emulate physics-based models, then having a diverse set of
independent simulations as target is key to a robust S2S intercomparison
project.

On Physics Metrics. The flag for physics-based metrics indicates whether these
benchmarks incorporate not just physical explanation, but also formulate them
as scalar and differentiable metrics for future optimization problem.

On Spatial Extent. Furthermore, the spatial extent indicates the extent of the
target benchmark, rather than of the input dataset. This is because some of
the more challenging S2S forecasting task is to get the correct global space-
time correlation, and having a full global coverage provides a more complete
evaluation.

On Temporal Extent. Lastly the temporal extent measures the greatest
overlapping window of all input dataset. For instance, in the case of
SubseasonalClimateUSA where they have different temporal resolution (25 input
variables from 1981-2023, 23 from 1979-2023 (all except SST and SIC), and 6
from 1948-2023), we measure the most common complete years (i.e., 1981 - 2023)
as these represent dense observations that can reliably constrain data-driven
emulators.

##  B Physics-Based Baseline Models

In this section, we describe in detail the physics-based models used as
baselines in ChaosBench. Wherever possible, we discuss specific strategies
regarding coupling to the ocean, sea ice, wave, land, initialization and
perturbation strategies, specifications of initial/boundary conditions, as
well as other numerical considerations to generate forecast.

###  B.1 The UK Meteorological Office (UKMO) [60]

  * •

Initialization and Ensemble. The UKMO model employs the lagged initialization
strategy to generate an ensemble of forecasts (4 in this case) at different
initialization time to improve prediction stability.

  * •

Coupling with ocean is performed with the Global Ocean 6.0 model [51], based
on NEMO3.6 [29] with 0.25 degree horizontal resolution and 75 vertical
pressure levels. The ocean model is initialized and calibrated using Nonlinear
Evolutionary Model VARiation (NEMOVAR) [31], a specific data assimilation
strategy that uses temperature, salinity profiles, altimeter-derived sea level
anomalies to calibrate forecasts. Frequency of coupling is 1-hourly.

  * •

Coupling with sea ice is performed with the Global Sea Ice 8.1 (CICE5.1.2)
model [54], and again initialized from NEMOVAR.

  * •

Coupling with wave model is not yet operational.

  * •

Coupling with land surface is performed with the Joint UK Land Environment
Simulator (JULES) [2]. Soil moisture, soil temperature, and snow are
initialized using JULES and forced using the the Japanese 55-year Reanalysis
(JRA-55) data [19]. The land surface model is paramaterized by land cover type
from a combination of satellite (e.g., MODIS LAI [53]) and radiometer data
(e.g., AVHRR [26]). In addition, another parameterization in the form of soil
characteristics is derived from the Harmonized World Soil Database [59].

  * •

Model grid uses the Arakawa C-grid [1] to solve partial differential equations
on a spherical surface. In particular, the velocity components (such as zonal
and meridional wind) are defined at the center of each face of the grid cells
(in the case of a rectilinear grid) or along cell edges (in the case of a
curvilinear grid). The scalar quantities such as pressure or temperature are
computed at the corners of the grid cells.

  * •

Large-scale dynamics uses the Semi-Lagrangian approach. It does not strictly
follow fluid parcels (i.e., Lagrangian), but it does calculate the value of a
field, such as temperature (i.e., Eulerian) by tracing back along the
trajectory that a fluid parcel would have taken to reach a specific point at
the current time step. This backward trajectory is used to find the origin of
the fluid parcel and determine its properties, which are then used to update
the model fields. This hybrid approach is therefore termed Semi-Lagrangian.

###  B.2 National Centers for Environmental Prediction (NCEP) [45]

  * •

Initialization and Ensemble. The NCEP model adds small perturbation to the
atmospheric, oceanic and land analysis at each cycle across 4 ensemble to
reduce sensitivity to initial conditions.

  * •

Coupling with ocean is performed with the GFDL Modular Ocean Model version 4
(MOM4) model that has a spatial resolution of 0.5-degree and 0.25-degree in
the longitude-latitude directions [11]. There are 40 vertical pressure levels.

  * •

Coupling with sea ice is also performed with the GFDL Sea Ice Simulator (SIS),
which models the thermodynamics and overall dynamics of sea ice [11].

  * •

Coupling with wave model is not yet operational.

  * •

Coupling with land surface is performed with 4-layer Noah Land surface model
2.7.1 [10]. Soil moisture, soil temperature, and snow are initialized using
Noah and forced using the Climate Forecast System [45] and the Global Land
Data Assimilation System [30] reanalysis data. The land surface model is
parameterized by land cover type AVHRR data. In addition, another
paramaterization in the form of soil characteristics is derived from the world
soil climate database [66].

  * •

Model grid uses the Gaussian grid [14], where the longitude (x-axis) are
evenly spaced while the latitudes (y-axis) are not. Instead, they are
determined by the roots of the associated Legendre polynomials, which
correspond to the Gaussian quadrature points for the sphere. This ensures that
the actual area represented by each grid cell is more uniform.

  * •

Large-scale dynamics uses the Spectral approach. It solves partial
differential equations by transforming them from the physical space into the
spectral domain. In the latter case, the equations are transformed into a
series of coefficients that represent the amplitude of waves across scales.
The transformations are usually done using Fourier series for periodic domains
or spherical harmonics when dealing with the whole Earth’s surface [14]. This
method is especially beneficial for smooth functions and for representing
large-scale wave phenomena, such as the Rossby waves, which are important for
understanding weather and climate.

###  B.3 China Meteorological Administration (CMA) [62]

  * •

Initialization and Ensemble. The CMA model uses the lagged average forecasting
(LAF) method across 4 ensemble members to ensure that the mean forecast is
less sensitive to initial conditions.

  * •

Coupling with ocean is performed with the GFDL MOM4 model, which has 40
vertical pressure levels [11]. Frequency of coupling is 2-hourly.

  * •

Coupling with sea ice is performed with the GFDL Sea Ice Simulator (SIS),
similar to that used by NCEP [11].

  * •

Coupling with wave model is not yet operational.

  * •

Coupling with land surface is performed with the Atmosphere-Vegetation
Interaction Model version 2 (AVIM2) model [61] and the NCAR NCAR Community
Land Model version 3.0 (CLMv3) [22]. Soil moisture, soil temperature, and snow
are not initialized directly using reanalysis data, as used by other land
surface models. Rather, air-sea-land-ice coupled model is forced by near-
surface atmospheric and ocean reanalysis in a long-term integration, and the
land initial conditions are produced as a by-product. As a result, the
parameterization of land cover type is done by this process, while soil
characteristics is derived from the Harmonized World Soil Database [59].

  * •

Model grid uses the Gaussian grid [14], similar to that used by the NCEP.

  * •

Large-scale dynamics uses a mixture of Spectral approach for the vorticity,
temperature, and surface pressure, as well as Semi-Lagrangian for specific
humidity and cloud waters other tracers.

###  B.4 European Center for Medium-Range Weather Forecasts (ECMWF) [8]

  * •

Initialization and Ensemble. The operational IFS forecast is generated through
Singular Vectors (SV) method: it creates a variety of initial conditions by
adjusting certain parameters slightly, thus generating different starting
points.

  * •

Coupling with ocean is performed with NEMO3.4.1 with 1-degree resolution and
42 vertical pressure levels. Frequency of coupling is 3-hourly.

  * •

Coupling with sea ice is not operational for this model’s version (but it is
in the newer generation, though the forecast start-date is much later than
2016). As a result, sea ice initial conditions are persisted up to day 15 and
then relaxed to climatology up to day 45.

  * •

Coupling with wave model is performed with ECMWF wave model with 0.5-degree
resolution [16].

  * •

Coupling with land surface is relatively more complex than the rest, and we
refer readers to their documentation. Regardless, it is based on Land Data
Assimilation System (LDAS) that combines heterogenous high-quality dataset
from satellite to ground sensors, and integrated with the operational IFS
model. The parameterization for land cover type is primarily based on MODIS
collection 5 [53] and soil characteristics from the FAO dominant soil texture
class [55].

  * •

Model grid uses the Cubic Octohedral grid [28], where the Earth’s surface is
projected onto a cube. Then, the cube is further subdivided to form an
octahedron, where the faces represent finer grid cells. This multi-scale
gridding scheme allows for parallelization where processes at different scales
could be solved simultaneously.

  * •

Large-scale dynamics uses a mixture of Spectral and Semi-Lagrangian approach,
similar to that used by CMA.

##  C Data-Driven Baseline Models

In this section, we describe in detail implementation and hyperparemeter
selections of our data-driven models used as baselines to ChaosBench. Most of
the choices are based on the original works that are adapted to weather and
climate applications using similar input dataset.

###  C.1 Lagged Autoencoder (AE)

We implement lagged AE from [6] with 5 $encoder$ blocks and 5 $decoder$ block,
with detailed specification in Table S1. Each encoder block is comprised of
$\textsc{MaxPool2d}\circ(\textsc{Conv2d}\rightarrow\textsc{BatchNorm2d}\rightarrow\textsc{ReLU}\rightarrow\textsc{Conv2d}\rightarrow\textsc{BatchNorm2d}\rightarrow\textsc{ReLU})$.
Similarly, the decoder block is comprised of
$\textsc{ConvTranspose2d}\rightarrow\textsc{BactNorm2d}\rightarrow\textsc{ReLU})\bigoplus(\textsc{ConvTranspose2d}\rightarrow\textsc{BactNorm2d}\rightarrow\textsc{Sigmoid})\circ(\textsc{Conv2d})$.

Hyperparameters | Values
---|---
Channels | [64, 128, 256, 512, 1024]
Encoder Kernel | $3\times 3$
Decoder Kernel | $2\times 2$
Max Pooling Window | $2\times 2$
Batch Normalization | True
Optimizer |  AdamW [25]
Learning Rate | $\textsc{CosineAnnealing}(10^{-2}\rightarrow 10^{-3})$
Batch Size | 32
Epochs | 500
Tmax | 500
Table S1: Hyperparameters for Lagged AE

###  C.2 ResNet

We adapt ResNet implementation from [42] using ResNet-50 as feature extractor
and 5 $decoder$ blocks, following specification in Table S2. Each decoder
block is composed of
$\textsc{ConvTranspose2d}\rightarrow\textsc{BactNorm2d}\rightarrow\textsc{LeakyReLU}$.

Hyperparameters | Values
---|---
Backbone | ResNet-50
Decoder Channels | [1024, 512, 256, 128, 64]
Decoder Activation | LeakyReLU(0.15)
Optimizer | AdamW
Learning Rate | $\textsc{CosineAnnealing}(10^{-2}\rightarrow 10^{-3})$
Batch Size | 32
Epochs | 500
Tmax | 500
Table S2: Hyperparameters for ResNet

###  C.3 Fourier Neural Operator (FNO)

We adapt FNO implementation from [23], following specification in Table S3 and
illustrated in S1. We implement the encoder-decoder structure, where we (1)
first transform our input $\mathbf{X}_{t}$ by convolutional layers both in the
Fourier (applying fast fourier transform; FFT) and physical domains, before we
concatenate both (applying inverse FFT for the former convolved features), and
apply non-linear GeLU activation function [12]. We select only the first 4
main Fourier modes to make the number of trainable parameters comparable with
the other data-driven baseline models. The (2) decoder block then applies
deconvolutional operation to the latent features to generate output $Y_{t}$.

Hyperparameters | Values
---|---
Non-Spectral Channels | [64, 128, 256, 512, 1024]
Spectral Channel | [64, 128, 256, 512, 1024]
Activation | GeLU
Fourier Modes | (4,4)
Optimizer | AdamW
Learning Rate | $\textsc{CosineAnnealing}(10^{-2}\rightarrow 10^{-3})$
Batch Size | 32
Epochs | 500
Tmax | 500
Table S3: Hyperparameters for FNO Figure S1: FNO architecture: (1) in the
encoder block, we transform our input $X_{t}$ by convolutional layers both in
the Fourier and physical domains, before we concatenate and apply non-linear
GeLU activation function. The (2) decoder block is then applying
deconvolutional operation to the latent features to generate forecast $Y_{t}$.

###  C.4 UNet

We adapt UNet implementation from [43] using 5 $encoder$ and 5 $decoder$
blocks, with skip connections, following specification in Table S4. The
composition of the encoder and decoder components are similar to those
described for Lagged Autoencoder, with the addition of skip connection between
each corresponding contracting-expansive path.

Hyperparameters | Values
---|---
Channels | [64, 128, 256, 512, 1024]
Activation | LeakyReLU(0.15)
Encoder Kernel | $3\times 3$
Decoder Kernel | $2\times 2$
Max Pooling Window | $2\times 2$
Optimizer | AdamW
Learning Rate | $\textsc{CosineAnnealing}(10^{-2}\rightarrow 10^{-3})$
Batch Size | 32
Epochs | 500
Tmax | 500
Table S4: Hyperparameters for UNet

###  C.5 ClimaX

ClimaX is based on the ViT model [9] with variational positional embedding in
variable-time space. We use ClimaX model as is described and implemented in
the original paper and is pre-trained using CMIP6 [36]. We fine-tune the
original pre-trained model given our training setup, for both Task 1 and 2.

###  C.6 PanguWeather, FourCastNetV2, GraphCast

We perform inference using their latest checkpoints using the API provided
here: https://github.com/ecmwf-lab/ai-models.

##  D Vision-Based Metrics

We describe in detail the four primary vision-based metrics used for this
benchmark, including RMSE, Bias, ACC, and MS-SSIM.

###  D.1 Root Mean-Squared Error (RMSE)

As described in the main text, we apply latitude-adjustment to RMSE
computation.

$\mathcal{M}_{RMSE}=\sqrt{\frac{1}{|\boldsymbol{\theta}||\boldsymbol{\gamma}|}\sum_{i=1}^{|\boldsymbol{\theta}|}\sum_{j=1}^{|\boldsymbol{\gamma}|}w(\theta_{i})(\hat{\mathbf{Y}}_{i,j}-\mathbf{Y}_{i,j})^{2}}$
(S1)

###  D.2 Bias

Similarly, we apply latitude-adjustment to Bias computation.

$\mathcal{M}_{Bias}=\frac{1}{|\boldsymbol{\theta}||\boldsymbol{\gamma}|}\sum_{i=1}^{|\boldsymbol{\theta}|}\sum_{j=1}^{|\boldsymbol{\gamma}|}w(\theta_{i})(\hat{\mathbf{Y}}_{i,j}-\mathbf{Y}_{i,j})$
(S2)

###  D.3 Anomaly Correlation Coefficient (ACC)

We remove the indexing for a more compact representation where the summation
is performed over each grid cell $(i,j)$. The predicted and observed anomalies
at each grid-cell are denoted by
$A_{\hat{\mathbf{Y}}_{i,j}}=\hat{\mathbf{Y}}_{i,j}-C$ and
$A_{\mathbf{Y}_{i,j}}=\mathbf{Y}_{i,j}-C$, where $C$ is the observational
climatology. We apply latitude-adjustment to ACC computation.

$\mathcal{M}_{ACC}=\frac{\sum w(\theta)[A_{\hat{\mathbf{Y}}}\cdot
A_{\mathbf{Y}}]}{\sqrt{\sum w(\theta)A^{2}_{\hat{\mathbf{Y}}}\sum
w(\theta)A^{2}_{\mathbf{Y}}}}$ (S3)

###  D.4 Multiscale Structural Similarity Index Measure (MS-SSIM)

Let $\mathbf{Y}$ and $\hat{\mathbf{Y}}$ be two images to be compared, and let
$\mu_{\mathbf{Y}}$, $\sigma^{2}_{\mathbf{Y}}$ and
$\sigma_{\mathbf{Y}\hat{\mathbf{Y}}}$ be the mean of $\mathbf{Y}$, the
variance of $\mathbf{Y}$, and the covariance of $\mathbf{Y}$ and
$\hat{\mathbf{Y}}$, respectively. The luminance, contrast and structure
comparison measures are defined as follows:

$l(\mathbf{Y},\hat{\mathbf{Y}})=\frac{2\mu_{\mathbf{Y}}\mu_{\hat{\mathbf{Y}}}+C_{1}}{\mu_{\mathbf{Y}}^{2}+\mu_{\hat{\mathbf{Y}}}^{2}+C_{1}},$
(S4)
$c(\mathbf{Y},\hat{\mathbf{Y}})=\frac{2\sigma_{\mathbf{Y}}\sigma_{\hat{\mathbf{Y}}}+C_{2}}{\sigma_{\mathbf{Y}}^{2}+\sigma_{\hat{\mathbf{Y}}}^{2}+C_{2}},$
(S5)
$s(\mathbf{Y},\hat{\mathbf{Y}})=\frac{\sigma_{\mathbf{Y}\hat{\mathbf{Y}}}+C_{3}}{\sigma_{\mathbf{Y}}\sigma_{\hat{\mathbf{Y}}}+C_{3}},$
(S6)

where $C_{1}$, $C_{2}$ and $C_{3}$ are constants given by

$C_{1}=(K_{1}L)^{2},C_{2}=(K_{2}L)^{2},\text{ and }C_{3}=C_{2}/2.$ (S7)

$L=255$ is the dynamic range of the gray scale images, and $K_{1}\ll 1$ and
$K_{2}\ll 1$ are two small constants. To compute the MS-SSIM metric across
multiple scales, the images are successively low-pass filtered and down-
sampled by a factor of 2. We index the original image as scale 1, and the
desired highest scale as scale $M$. At each scale, the contrast comparison and
structure comparison are computed and denoted as
$c_{j}(\mathbf{Y},\hat{\mathbf{Y}})$ and $s_{j}(\mathbf{Y},\hat{\mathbf{Y}})$
respectively. The luminance comparison is only calculated at the last scale
$M$, denoted by $l_{M}(\mathbf{Y},\hat{\mathbf{Y}})$. Then, the MS-SSIM metric
is defined by

$\mathcal{M}_{MS-
SSIM}=[l_{M}(\mathbf{Y},\hat{\mathbf{Y}})]^{\alpha_{M}}\cdot\prod_{j=1}^{M}[c_{j}(\mathbf{Y},\hat{\mathbf{Y}})]^{\beta_{j}}[s_{j}(\mathbf{Y},\hat{\mathbf{Y}})]^{\gamma_{j}}$
(S8)

where $\alpha_{M}$, $\beta_{j}$ and $\gamma_{j}$ are parameters. We use the
same set of parameters as in [57]: $K_{1}=0.01$, $K_{2}=0.03$, $M=5$,
$\alpha_{5}=\beta_{5}=\gamma_{5}=0.1333$, $\beta_{4}=\gamma_{4}=0.2363$,
$\beta_{3}=\gamma_{3}=0.3001$, $\beta_{2}=\gamma_{2}=0.2856$,
$\beta_{1}=\gamma=0.0448$. The predicted and ground truth images of physical
variables are re-scaled to 0-255 prior to the calculation of their MS-SSIM
values.

##  E Physics-Based Metrics

In this section, we describe in detail the definition and implementation of
our physics-based metrics, including PyTorch psuedocode implementation.

Let $\mathbf{Y}$ be a 2D image of size $h\times w$ for a physical variables at
a specific time, variable, and level. Let $f(x,y)$ be the intensity of the
pixel at position $(x,y).$ First, we compute the 2D Fourier transform of the
image by

$F(k_{x},k_{y})=\sum_{x=0}^{w-1}\sum_{y=0}^{h-1}f(x,y)\cdot e^{-2\pi
i\left(k_{x}x/w+k_{y}y/h\right)},$ (S9)

where $k_{x}$ and $k_{y}$ correspond to the wavenumber components in the
horizontal and vertical directions, respectively, and $i$ is the imaginary
unit. The power at each wavenumber component $(k_{x},k_{y})$ is given by the
square of the magnitude spectrum of $F(k_{x},k_{y})$, that is,

$S(k_{x},k_{y})=|F(k_{x},k_{y})|^{2}=\text{Re}[F(k_{x},k_{y})]^{2}+\text{Im}[F(k_{x},k_{y})]^{2}.$
(S10)

The scalar wavenumber is defined as:

$k=\sqrt{k_{x}^{2}+k_{y}^{2}},$ (S11)

which represents the magnitude of the spatial frequency vector, indicating how
rapidly features change spatially regardless of direction. Then, the energy
distribution at a spatial frequency corresponding to k is defined as

$S(k)=\sum_{(k_{x},k_{y}):\sqrt{k_{x}^{2}+k_{y}^{2}}=k}S(k_{x},k_{y}).$ (S12)

Given the spatial energy frequency distribution for observations $E(k)$ and
predictions $\hat{S}(k)$ , we perform normalization for each over
$\mathbf{K}_{q}$, the set of wavenumbers corresponding to high-frequency
components of energy distribution, as defined in Equation S13. This is to
ensure that the sum of the component sums up to 1 which exhibits pdf-like
property.

$S^{\prime}(k)=\frac{S(k)}{\sum_{k\in\mathbf{K}_{q}}S(k)},\quad\hat{S}^{\prime}(k)=\frac{\hat{S}(k)}{\sum_{k\in\mathbf{K}_{q}}\hat{S}(k)},\quad
k\in\mathbf{K}_{q}$ (S13)

⬇

1import torch

2import torch.nn as nn

3

4class SpectralDiv(nn.Module):

5 """

6 Compute Spectral divergence given the top-k percentile wavenumber (higher k
means higher frequency)

7 """

8 def __init__(

9 self,

10 percentile=0.9,

11 input_shape=(121,240)

12 ):

13 super(SpectralDiv, self).__init__()

14

15 self.percentile = percentile

16

17 # Compute the discrete Fourier Transform sample frequencies for a signal of
size

18 nx, ny = input_shape

19 kx = torch.fft.fftfreq(nx) * nx

20 ky = torch.fft.fftfreq(ny) * ny

21 kx, ky = torch.meshgrid(kx, ky)

22

23 # Construct discretized k-bins

24 self.k = specify_k_bins(...)

25

26 # Get k-percentile index

27 self.k_percentile_idx = int(len(self.k) * self.percentile)

28

29 def forward(self, predictions, targets):

30

31 # Preprocess data, including handling of missing values, etc

32 predictions = preprocess_data(...)

33 targets = preprocess_data(...)

34

35 # Compute along mini-batch

36 predictions, targets = torch.nanmean(predictions, dim=0),
torch.nanmean(targets, dim=0)

37

38 # Transform prediction and targets onto the Fourier space and compute the
power

39 predictions_power = torch.fft.fft2(predictions)

40 predictions_power = torch.abs(predictions_power)**2

41

42 targets_power = torch.fft.fft2(targets)

43 targets_power = torch.abs(targets_power)**2

44

45 # Normalize as pdf

46 predictions_Sk = predictions_power / torch.nansum(predictions_power)

47 targets_Sk = targets_power / torch.nansum(targets_power)

48

49 # Compute spectral Sk divergence

50 div = torch.nansum(targets_Sk * torch.log(torch.clamp(targets_Sk /
predictions_Sk, min=1e-9)))

51

52 return div

Listing S1: Psuedocode for computing SpecDiv using PyTorch

⬇

1import torch

2import torch.nn as nn

3

4class SpectralRes(nn.Module):

5 """

6 Compute Spectral residual given the top-k percentile wavenumber (higher k
means higher frequency)

7 """

8 def __init__(

9 self,

10 percentile=0.9,

11 input_shape=(121,240)

12 ):

13 super(SpectralRes, self).__init__()

14

15 self.percentile = percentile

16

17 # Compute the discrete Fourier Transform sample frequencies for a signal of
size

18 nx, ny = input_shape

19 kx = torch.fft.fftfreq(nx) * nx

20 ky = torch.fft.fftfreq(ny) * ny

21 kx, ky = torch.meshgrid(kx, ky)

22

23 # Construct discretized k-bins

24 self.k = specify_k_bins(...)

25

26 # Get k-percentile index

27 self.k_percentile_idx = int(len(self.k) * self.percentile)

28

29 def forward(self, predictions, targets):

30

31 # Preprocess data, including handling of missing values, etc

32 predictions = preprocess_data(...)

33 targets = preprocess_data(...)

34

35 # Compute along mini-batch

36 predictions, targets = torch.nanmean(predictions, dim=0),
torch.nanmean(targets, dim=0)

37

38 # Transform prediction and targets onto the Fourier space and compute the
power

39 predictions_power = torch.fft.fft2(predictions)

40 predictions_power = torch.abs(predictions_power)**2

41

42 targets_power = torch.fft.fft2(targets)

43 targets_power = torch.abs(targets_power)**2

44

45 # Normalize as pdf

46 predictions_Sk = predictions_power / torch.nansum(predictions_power)

47 targets_Sk = targets_power / torch.nansum(targets_power)

48

49 # Compute spectral Sk residual

50 res = torch.sqrt(torch.nanmean(torch.square(predictions_Sk - targets_Sk)))

51

52 return res

Listing S2: Psuedocode for computing SpecRes using PyTorch

##  F Extended Result and Discussion

###  F.1 Results on Short-Range Forecasting

Despite its high predictability, nowcasting is a popular task due to its
perceived urgency [50]. Nowcasting results are summarized in Table S5.

| RMSE $\downarrow$ | Bias $\rightarrow 0$ | MS-SSIM $\uparrow$ | SpecDiv $\downarrow$
---|---|---|---|---
Models |  T850 ($K$) |  Z500 ($gpm$) |  Q700 ($\times 10^{-3}$) | T850 | Z500 | Q700 | T850 | Z500 | Q700 | T850 | Z500 | Q700
Climatology* | 3.38 | 81.03 | 1.63 | -0.19 | -5.08 | +0.18 | 0.85 | 0.82 | 0.62 | 0.01 | 0.01 | 0.03
Persistence* | 3.09 | 61.88 | 1.74 | -0.01 | -0.02 | -0.01 | 0.88 | 0.91 | 0.65 | 0.01 | 0.01 | 0.01
UKMO | 1.17 | 6.58 | 0.98 | +0.05 | +1.75 | +0.20 | 0.87 | 0.99 | 0.93 | 0.07 | 0.09 | 0.05
NCEP | 1.52 | 10.88 | 1.24 | -0.26 | -1.06 | +0.32 | 0.96 | 0.99 | 0.90 | 0.50 | 0.55 | 0.14
CMA | 1.66 | 7.94 | 1.41 | +0.11 | +1.51 | +0.04 | 0.96 | 0.99 | 0.90 | 0.04 | 0.04 | 0.04
ECMWF | 0.92 | 7.85 | 0.80 | -0.08 | +6.17 | +3.23 | 0.99 | 0.99 | 0.96 | 0.05 | 0.06 | 0.05
| Task 1: Full Dynamics Prediction
Lagged AE | 4.07 | 90.44 | 1.76 | -0.52 | -11.63 | -0.11 | 0.81 | 0.79 | 0.56 | 0.83 | 0.83 | 0.31
FNO | 3.63 | 68.88 | 1.71 | -0.25 | -11.36 | -0.03 | 0.81 | 0.82 | 0.56 | 0.19 | 0.09 | 0.16
ResNet | 2.21 | 30.93 | 1.27 | -0.39 | -9.78 | -0.04 | 0.94 | 0.98 | 0.83 | 0.34 | 0.63 | 0.19
UNet | 2.00 | 28.28 | 1.02 | -0.36 | -10.50 | -5.56 | 0.96 | 0.98 | 0.89 | 0.80 | 1.44 | 0.49
| Task 2: Sparse Dynamics Prediction
Lagged AE | 1.83 | 30.67 | 1.15 | -0.63 | -14.35 | -0.06 | 0.96 | 0.98 | 0.87 | 0.51 | 0.98 | 0.10
FNO | 2.82 | 52.00 | 1.52 | -0.37 | -5.03 | +0.01 | 0.88 | 0.89 | 0.65 | 0.26 | 0.26 | 0.23
ResNet | 1.85 | 30.39 | 1.15 | -0.49 | -10.09 | -0.07 | 0.96 | 0.98 | 0.88 | 0.05 | 0.69 | 0.07
UNet | 1.48 | 24.86 | 0.95 | -0.51 | -12.22 | -0.04 | 0.98 | 0.99 | 0.93 | 0.02 | 0.71 | 0.03
Table S5: Short-range nowcasting ($\Delta t=1$) results on select metrics and
target variables between physics-based and data-driven models. _(*) Baseline
model that uses privileged information (target observations) to make
prediction_.

Observation 1: First, we find that the best data-driven model in terms of
vision-based metrics, UNet, does not consistently outperform on physics-based
metrics.

Observation 2: On the contrary, some data-driven models, especially FNO,
outperform others in physics-based metrics while underperform on vision-based
ones. This is by construction since FNO inherently learns the underlying
spectral profiles, unlike other data-driven models.

Observation 3: Baseline models with privileged information, such as
climatology and persistence, tend to underperform on vision-based metrics, but
outperform on physics-based metrics. This is natural since they have access to
the physics embedded in actual observations.

Observation 4: Barring climatology and persistence, most physics-based models
tend to consistently outperform data-driven models on both vision and physics-
based metrics (e.g., T850 and Z500).

Observation 5: For sanity check, we include Lagged AE because of its capacity
to generate compact representation at the expense of losing high-frequency
signal [27]. As expected, Lagged AE has the worse SpecDiv scores across.

Overall, we find models that perform well on vision-based metrics tend to be
physically-inconsistent, and vice versa.

###  F.2 Results on S2S-Range Forecasting

| RMSE $\downarrow$ | Bias $\rightarrow 0$ | MS-SSIM $\uparrow$ | SpecDiv $\downarrow$
---|---|---|---|---
Models |  T850 ($K$) |  Z500 ($gpm$) |  Q700 ($\times 10^{-3}$) | T850 | Z500 | Q700 | T850 | Z500 | Q700 | T850 | Z500 | Q700
Climatology* | 3.38 | 81.03 | 1.63 | -0.19 | -5.08 | +0.18 | 0.85 | 0.82 | 0.62 | 0.01 | 0.01 | 0.03
Persistence | 5.88 | 127.83 | 2.47 | -0.02 | -0.17 | -0.01 | 0.71 | 0.69 | 0.41 | 0.02 | 0.03 | 0.05
UKMO | 5.00 | 116.2 | 2.32 | +0.18 | +1.82 | +1.89 | 0.64 | 0.71 | 0.43 | 0.06 | 0.09 | 0.07
NCEP | 4.90 | 116.7 | 2.30 | -0.94 | -16.3 | +0.25 | 0.75 | 0.71 | 0.43 | 0.53 | 0.55 | 0.10
CMA | 5.08 | 118.7 | 2.49 | +0.58 | +14.2 | +0.31 | 0.75 | 0.72 | 0.45 | 0.05 | 0.04 | 0.06
ECMWF | 4.72 | 115.1 | 2.30 | -0.11 | +6.40 | -2.55 | 0.75 | 0.72 | 0.44 | 0.06 | 0.07 | 0.06
UNet | 11.8 | 248.1 | 3.17 | -1.75 | -4.61 | -0.24 | 0.40 | 0.46 | 0.01 | 0.29 | 0.66 | 0.95
UNet${\dagger}^{a}$ | 9.56 | 189.5 | 2.89 | +0.18 | +26.10 | +0.03 | 0.61 | 0.66 | 0.05 | 1.05 | 0.70 | 0.93
UNet${\dagger}^{b}$ | 4.20 | 99.5 | 1.89 | -0.27 | -7.21 | -0.70 | 0.80 | 0.76 | 0.52 | 0.49 | 0.67 | 0.57
PanguWeather${\dagger}^{a}$ | 5.85 | 120.9 | 2.34 | -1.56 | -30.50 | -0.45 | 0.70 | 0.68 | 0.43 | 0.25 | 0.35 | 0.23
FourCastNetV2${\dagger}^{a}$ | 5.11 | 112.4 | - | -0.24 | -0.74 | - | 0.74 | 0.72 | - | 0.28 | 0.11 | -
GraphCast${\dagger}^{a}$ | 5.87 | 136.1 | 2.28 | -0.17 | +7.46 | +0.06 | 0.70 | 0.66 | 0.45 | 0.05 | 0.03 | 0.27
ClimaX${\dagger}^{b}$ | 3.56 | 83.1 | 1.66 | -0.26 | -7.58 | -0.01 | 0.83 | 0.81 | 0.59 | 0.20 | 0.13 | 0.28
Table S6: Long-range forecasting ($\Delta t=44$) results on select metrics and
target variables between physics-based and data-driven models. Results are
specifically for Task 1. _(*) Baseline model that uses privileged information
(label) to make prediction; ( ${\dagger}^{a}$) Data-driven model with
autoregressive training approach; (${\dagger}^{b}$) Data-driven model with
direct training approach_.

Overall, we find models perform significantly worse on all metrics, visual or
physics-based, when compared to climatology.

###  F.3 Effects of Autoregressive Training Scheme

We showcased more results for autoregressive training strategy. In this case,
we performed autoregressive training using 5 iterative steps (S=5). As
illustrated in Figure S2, we observe that incorporating temporal information
improve the vision-based metrics even at longer forecasting timesteps, with
lower RMSE, higher ACC/MS-SSIM. However, the converse trend is true
incorporating temporal context makes S2S forecast worse off in physics-based
scores, especially in Task 1.

(a) RMSE for Task 1 ($\downarrow$ is better)

(b) RMSE for Task 2 ($\downarrow$ is better)

(c) ACC for Task 1 ($\uparrow$ is better)

(d) ACC for Task 2 ($\uparrow$ is better)

(e) MS-SSIM for Task 1 ($\uparrow$ is better)

(f) MS-SSIM for Task 2 ($\uparrow$ is better)

(g) SDiv for Task 1 ($\downarrow$ is better)

(h) SDiv for Task 2 ($\downarrow$ is better)

Figure S2: Ablation results for incorporating temporal information in an
autoregressive scheme for long-range forecast using UNet models. Panels
(a,c,e,g) represent RMSE, ACC, MS-SSIM, and SDiv for Task 1; Panels (b,d,f,h)
represent RMSE, ACC, MS-SSIM, and SDiv for Task 2. The x-axis represents the
number of forecasting days for T850, Z500, Q700 representative tasks. Blue and
orange lines illustrate autoregressive scheme with $S=1$ and $S=5$
respectively. Overall we observe that incorporating temporal information
improve the vision-based metrics even at longer forecasting timesteps, with
lower RMSE, higher ACC/MS-SSIM. However, the converse trend is true where
incorporating temporal context makes S2S forecast worse off in physics-based
scores, especially in Task 1.

###  F.4 Vision-Based S2S Results

We present extended results that demonstrate the forecasts at $t=1$ and $t=44$
for the best-performing data-driven baseline model across tasks, iterative
schemes, and representative variables. We highlight results for both
autoregressive (S=1, S=5) and direct temporal schemes for UNet and ClimaX.

####  F.4.1 UNet Autoregressive Approach

(a) UNet task 1 using autoregressive approach (S=1)

(b) UNet task 2 using autoregressive approach (S=1)

Figure S3: Normalized temperature@850-hpa qualitative results for
autoregressive approach (S=1) of UNet.

(a) UNet task 1 using autoregressive approach (S=1)

(b) UNet task 2 using autoregressive approach (S=1)

Figure S4: Normalized temperature@850-hpa qualitative results for
autoregressive approach (S=5) of UNet.

(a) UNet task 1 using autoregressive approach (S=1)

(b) UNet task 2 using autoregressive approach (S=1)

Figure S5: Normalized geopotential@500-hpa qualitative results for
autoregressive approach (S=1) of UNet.

(a) UNet task 1 using autoregressive approach (S=5)

(b) UNet task 2 using autoregressive approach (S=5)

Figure S6: Normalized geopotential@500-hpa qualitative results for
autoregressive approach (S=5) of UNet.

(a) UNet task 1 using autoregressive approach (S=1)

(b) UNet task 2 using autoregressive approach (S=1)

Figure S7: Normalized humidity@700-hpa qualitative results for autoregressive
approach (S=1) of UNet.

(a) UNet task 1 using autoregressive approach (S=5)

(b) UNet task 2 using autoregressive approach (S=5)

Figure S8: Normalized humidity@700-hpa qualitative results for autoregressive
approach (S=5) of UNet.

####  F.4.2 UNet Direct Approach

(a) UNet task 1 using direct approach

(b) UNet task 2 using direct approach

Figure S9: Normalized temperature@850-hpa qualitative results for direct
approach for UNet.

(a) UNet task 1 using direct approach

(b) UNet task 2 using direct approach

Figure S10: Normalized geopotential@500-hpa qualitative results for direct
approach for UNet.

(a) UNet task 1 using direct approach

(b) UNet task 2 using direct approach

Figure S11: Normalized humidity@700-hpa qualitative results for direct
approach for UNet.

####  F.4.3 ClimaX Direct Approach

(a) ClimaX task 1 using direct approach

(b) ClimaX task 2 using direct approach

Figure S12: Normalized temperature@850-hpa qualitative results for direct
approach for ClimaX.

(a) ClimaX task 1 using direct approach

(b) ClimaX task 2 using direct approach

Figure S13: Normalized geopotential@500-hpa qualitative results for direct
approach for ClimaX.

###  F.5 Physics-Based S2S Results

We provide extended results on the power spectrum $S(k)$ at $t=1$ and $t=44$
across baseline models (Figures S14-S17). In general we find that high-
frequency signal is lost (less structured $S_{k}$) as we unroll longer into
the future, and that incorporating temporal information does not guarantee
improvement. In addition, from Figures S22-S25 for autoregressive approach and
Figures S26-S27 for direct approach, we observe relative stability of low-
frequency signals across models even at long rollouts. However, high-frequency
signals become unstructured at longer lead-time (except for FNOs which exhibit
relative stability). These results further emphasize the need to incorporate
physics-based metrics into the training/inference pipeline or bake in physical
principles during the modeling stage.

(a) Spectral profile of Lagged AE for S=1 and S=5 in Task 1

(b) Spectral profile of Lagged AE for S=1 and S=5 in Task 2

Figure S14: Spectral profile for Lagged AE at t=1 and t=44 using
autoregressive approach

(a) Spectral profile of FNO for S=1 and S=5 in Task 1

(b) Spectral profile of FNO for S=1 and S=5 in Task 2

Figure S15: Spectral profile for FNO at t=1 and t=44 using autoregressive
approach

(a) Spectral profile of ResNet for S=1 and S=5 in Task 1

(b) Spectral profile of ResNet for S=1 and S=5 in Task 2

Figure S16: Spectral profile for ResNet at t=1 and t=44 using autoregressive
approach

(a) Spectral profile of UNet for S=1 and S=5 in Task 1

(b) Spectral profile of UNet for S=1 and S=5 in Task 2

Figure S17: Spectral profile for UNet at t=1 and t=44 using autoregressive
approach

(a) Spectral profile of UNet-direct in Task 1

(b) Spectral profile of UNet-direct in Task 2

Figure S18: Spectral profile for UNet at t=1 and t=44 using direct approach

Figure S19: Spectral profile for PanguWeather at t=1 and t=44 Figure S20:
Spectral profile for FourCastNetV2 at t=1 and t=44 Figure S21: Spectral
profile for GraphCast at t=1 and t=44

(a) 3D Spectral profile for Lagged AE for Task 1

(b) 3D Spectral profile for Lagged AE for Task 2

Figure S22: 3D Spectral profile for Lagged AE across different tasks and using
autoregressive approach

(a) 3D Spectral profile for FNO for Task 1

(b) 3D Spectral profile for FNO for Task 2

Figure S23: 3D Spectral profile for FNO across different tasks and using
autoregressive approach

(a) 3D Spectral profile for ResNet for Task 1

(b) 3D Spectral profile for ResNet for Task 2

Figure S24: 3D Spectral profile for ResNet across different tasks and using
autoregressive approach

(a) 3D Spectral profile for UNet for Task 1

(b) 3D Spectral profile for UNet for Task 2

Figure S25: 3D Spectral profile for UNet across different tasks and using
autoregressive approach

(a) 3D Spectral profile for UNet-direct for Task 1

(b) 3D Spectral profile for UNet-direct for Task 2

Figure S26: 3D Spectral profile for UNet across different tasks and using
direct approach

(a) 3D Spectral profile for ClimaX-direct for Task 1

(b) 3D Spectral profile for ClimaX-direct for Task 2

Figure S27: 3D Spectral profile for ClimaX across different tasks and using
direct approach Figure S28: 3D Spectral profile for PanguWeather Figure S29:
3D Spectral profile for FourCastNetV2 Figure S30: 3D Spectral profile for
GraphCast