Missing the Butterfly and Predicting the Past: Features or Bugs of Accurate AI Weather Models?
Authors: Pedram Hassanzadeh, Weidong Li, Y. Qiang Sun, Jiangdi Wang, Alexander Wikner, Justin Finkel, Jonathan Q. Weare
Organizations: Department of Geophysical Sciences, University of Chicago, Chicago, USA · Committee on Computational and Applied Mathematics, University of Chicago, Chicago, USA. · Data Science Institute, University of Chicago, Chicago, USA. · State Key Laboratory of Severe Weather Meteorological Science and Technology, Nanjing University, Nanjing, China. · Key Laboratory of Mesoscale Severe Weather/Ministry of Education, Nanjing University, Nanjing, China. · Courant Institute School of Mathematics, Computing, and Data Science, New York University, New York City, USA.
AI weather prediction (AIWP) models rival physics-based models, yet the sources of their unexpected forecast accuracy and the degree of their physical fidelity remain unclear. Here, across a hierarchy spanning observation-based reanalysis, a general circulation model, and the multi-scale Lorenz system, we show that AI models can be trained to skillfully predict the past (backcast), though backcasts are systematically less accurate than forecasts. However, skillful backcasting appears to violate the second law of thermodynamics, and all these forecasting and backcasting models miss the butterfly effect. We trace the surprising forecast accuracy, missing butterfly, and skillful backcasting to a single cause: inevitable coarse-graining of training data, which removes fast, small scales and/or some variables. From the Lorenz system to official Pangu-Weather models, reducing coarse-graining makes AI predictions more physics-like (arrow of time and butterfly-like effects emerge), but forecast accuracy declines. Results offer an explanation for AIWP models' forecast skill: unlike physics-based models, they implicitly learn how fast, small scales affect large scales without inheriting their rapid error growth. Broader implications are that AI models' proliferation calls for revisiting predictability theories and long-term climate emulation strategies, and backcasting offers a useful, new lens for such analyses.
Figures & tables
Figure 1: Missing butterfly and forecast-backcast asymmetry in ERA5 and GCM AIWP models. (A) Schematic of the autoregressive forecasting and backcasting models. The forecast network NNθ advances the atmospheric state (Eq. ( 1 )); an independently trained backcast network NNϕ steps it backward (Eq. ( 2 )). (B) Anomaly correlation coefficient (ACC) of Z500 versus lead time for ERA5 (left; Transformer (deterministic), based on Pangu-Weather [ 2 ] ) and the PlaSim GCM (right; Transformer and conditional Diffusion (stochastic), similar to GenCast [ 4 ] ), averaged over all initial conditions from the test sets. Δt= 24 h and the horizontal resolution is ∼100 km (ERA5) and ∼300 km (PlaSim). Shading marks the ACC >0.6 skill window; annotations give the crossing times and their forecast/backcast ratio, the asymmetry. See Figs. S1 - S4 for other variables, metrics, Δt , and neural-operator and generative architectures. (C) Difference kinetic energy (DKE) for perturbations of different amplitudes, decreasing from a reference value, ϵ0 . We also show DKE from numerical integrations of physics-based models (ICON 20 km and PlaSim GCM). See Methods and Data for details about the AIWP models, physics-based models, perturbations, and metrics.
Figure 2: Effects of training data coarse-graining on forecasting and backcasting AI models of Lorenz 96. All results (RMSE, asymmetry, and Hovmöller diagrams) are based on the large-scale variable, Xi(i=1…8) . ( A ) The numerical-model reference trajectory alongside forecasts and backcasts from numerical integrations of Eqs. ( S7 )-( S8 ) initialized at lead time 0 ; saturated dark red/blue regions mark numerical blow-up in backcasts. Single floating-point precision calculations are shown to match the precision used in AI models’ inference; it does not make any noticeable difference. ( B ) and ( C ) show RMSE curves (averaged over 100 initial conditions from the test set) and Hovmöller diagrams (from the same initial condition as in (A)) for 12 independent AI models trained on only the large-scale Xi and on both scales ( Xi and Yi,j(j=1…32 )) with time interval Δt=10δt , 5δt , and 1δt ( δt is the time step of the numerical solver). Shading indicates the RMSE interquartile range. Annotation boxes report forecast RMSE at 4.5 days and the asymmetry (backcast/forecast RMSE ratio) at 1.5 days. Note the difference between the RMSE ranges in rows (B) and (C), as well as between the lead-time ranges, chosen for clearer illustration. See Methods and Data for more details of the system and AI models.
Figure 3: Effects of training-data coarse-graining on ensemble growth and the butterfly effect in AI models of Lorenz 96. Curves show DKE of the large-scale variable, Xi(i=1…8) for 12 independent, separate AI models that have an identical architecture (including the same MSE loss) trained on ( A ) only the large-scale Xi and ( B ) on both scales ( Xi and Yi,j(j=1…32 )) with time interval Δt=10δt , 5δt , and 1δt ( δt is the time step of the numerical solver). Gray curves show forecasts and backcasts from numerical integrations of the two-scale Lorenz 96 equations ( S7 )–( S8 ); crosses mark numerical blow-up in backcasts. Curves are averaged over the same 100 initial conditions from the test set. Perturbations are of different amplitudes, decreasing from a reference value, ϵ0 . See Methods and Data for details about the perturbations and metrics.
Figure 4: Effects of temporal coarse-graining of training data on ensemble growth and the butterfly effect in the Pangu-Weather forecasts. Columns show, from left to right, forecasts from the independent, separate official Pangu-Weather models [ 2 ] with time step Δt=24 , 6, 3, and 1 h. ( A ) DKE versus lead time, as in Fig. 1 C, for three amplitudes of initial perturbations. DKE is averaged over 60o S- 60o N. See Fig. S10 for results with global averaging and lower GPU precision. Gray curves show the DKE from ICON forecasts (from [ 19 ] ). The number in each panel is the skillful lead time, i.e., when ACC drops to 0.6, as in Fig. 1 B. ( B ) Power spectra Sv(n) of the 300-hPa meridional wind perturbation v300′ (each member minus ensemble mean) as a function of total wavenumber n , for the initial-condition perturbations with the smallest ( ϵ0/103 , upper row) and largest ( ϵ0 , lower row) amplitudes. These curves are averaged over the same ensemble members and initial conditions as the ERA5 Transformer in Fig. 1 C. Colors denote (nonuniform) lead times from +1 to +72 h; the gray curve is the imposed initial Perlin noise perturbation at t=0 , and the black curve is the spectrum of the full ERA5 field for reference. The top axis denotes the corresponding wavelength. Figure S6 shows the spectra of the forecast errors as a function of lead time. See Methods and Data for more information.
Appendix figures & tables16 assets
Supplementary material from the paper’s appendix.
Appendix
ERA5
PlaSim
Architecture
Transformer
Transformer
SFNO
Diffusion
Variable
Region
24 h
6 h
24 h
6 h
24 h
24 h
U250
Global
1.50
1.81
1.78
1.98
1.77
2.03
NH
1.36
1.66
1.68
1.84
1.68
1.94
Tropics
1.91
2.23
2.04
2.15
2.00
2.32
SH
1.42
1.74
1.72
1.87
1.79
1.99
Appendix
Table S1 : Regional variation of the forecast-backcast asymmetry . As defined in the main text, the asymmetry is the ratio of the lead times at which the forecast and backcast ACC drops to 0.6. Asymmetry >1 means that backcasting is less skillful than forecasting. Asymmetry values are reported globally and separately for the extratropical Northern Hemisphere (NH: 30 ∘ N–70 ∘ N), Tropics: 30 ∘ S–30 ∘ N, and extratropical Southern Hemisphere (SH: 70 ∘ S–30 ∘ S) for both ERA5 and PlaSim and for a variety of architectures and Δt . For each variable and AI model, the largest value among the three regions is shown in bold. For ERA5, consistent across both Δt values and variables at 3 vertical levels, the asymmetry is noticeably largest in the tropics. For PlaSim, the same trend is seen for Z500 and U250, but the trend is not robust for U850.
System / AI model
Projection
λ (MTU -1 )
Ground truth
X
21.9 ± 6.3
Ground truth
Y
24.3 ± 4.0
XY , 1δt AI model
X
26.8 ± 6.5
XY , 1δt AI model
Y
26.8 ± 5.9
X , 1δt AI model
X
3.05 ± 1.84
X , 10δt AI model
X
2.81 ± 1.59
Appendix
Table S2 : Leading Lyapunov exponents of Lorenz 96 numerical and AI models. The leading finite-time Lyapunov exponents (FTLE) λ are estimated using the rescaling method (see Methods and Data). The reported values are the mean and standard deviation over 20 independent initial conditions on the reference trajectory, each with an independent random initial perturbation of amplitude δ0=10−4 , computed over the window 0.1 – 0.3 MTU (after discarding a 0.1 -MTU transient). The FTLE values are projected onto the slow ( X ) or fast ( Y ) variables.
Category
Variable
Pangu-Weather [ 2 ]
ERA5 Transformer
PlaSim AIWP
Temperature ( T )
13 levels
17 levels
13 levels
Zonal wind ( u )
13 levels
17 levels
13 levels
Upper-air xua
Meridional wind ( v )
13 levels
17 levels
13 levels
Specific humidity ( q )
13 levels
17 levels
13 levels
Geopotential ( Z )
13 levels
17 levels
13 levels
2-m temperature
✓
✓
✓
Appendix
Table S3 : Variable inventories for the Pangu-Weather, ERA5, and PlaSim AIWP models. Upper-air prognostic state xua , surface prognostic state xsfc , time-varying boundary forcings b , and static fields c are listed. Checkmarks indicate inclusion; dashes indicate not applicable. The Pangu-Weather column describes the pretrained model on its native 0.25∘ grid, which we use for inference only. The ERA5 Transformer column describes the data for the 1∘ ERA5 Transformer. The PlaSim AIWP column describes data for the Transformer, SFNO, and Diffusion AIWP models. † We have trained versions of the PlaSim Transformer with sea-surface temperature as a prognostic variable (in x ) and the key observations and conclusions from Fig. 1 B-C remain the same.
ERA5 Transformer
PlaSim Transformer
PlaSim SFNO
PlaSim Diffusion
Embedding dim Demb
240
240
384
1024
Depths / blocks
2+6+6+2
2+6+6+2
12
12 self-attn. + 6 cross-attn.
Architecture
Attention heads
6/12/12/6
6/12/12/6
—
16
Window size
2×6×10
2×6×12
—
—
Patch size
2×2×2
2×2×2
—
2×2
Optimizer
AdamW
AdamW
AdamW
AdamW
Appendix
Table S4 : Architectural and training hyperparameters for ERA5 and PlaSim AIWP models. Entries marked with a dash are not applicable. † Forecast / backcast peak learning rates.
Phase
Hyperparameter
Value / Search space
Architecture
multilayer perceptron, ReLU activations
Architecture
Width W
{512,1024,2048}
Depth D
{3,5,7} hidden layers
Input dimension
8 ( X only) or 264 ( XY )
Optimizer
Adam
Hyperparameter search
Learning rate
Log-uniform [10−4,10−2]
Appendix
Table S5 : Architectural and training hyperparameters for the Lorenz 96 AI models. Set notation indicates discrete values explored during the hyperparameter optimization phase.
x
Δt
Forecast RMSE
Backcast RMSE
X
1δt
0.0056 ± 0.0020
0.0056 ± 0.0021
X
5δt
0.026 ± 0.0093
0.027 ± 0.011
X
10δt
0.048 ± 0.016
0.046 ± 0.018
XY
1δt
0.038 ± 0.012
0.043 ± 0.015
XY
5δt
0.16 ± 0.054
0.23 ± 0.067
XY
10δt
0.29 ± 0.073
0.26 ± 0.066
Appendix
Table S6 : One-time-step RMSE of the forecast and backcast AI models. RMSE of a single prediction step of length Δt , computed over the X variables against the reference trajectory, for the selected AI models with W=1024 and D=7 . Values are the mean ± standard deviation over the same 100 initial conditions as Fig. 2 .
Figure S1 : Missing butterfly and forecast-backcast asymmetry in the SFNO-based AIWP models of PlaSim GCM. As in Fig. 1 B–C, but for independently trained forecasting and backcasting AIWP models, based on spherical Fourier neural operator (SFNO) architecture [ 71 ] , with Δt=24 h (see Methods and Data). ACC of U850 , U250 , and Z500 versus lead time are shown, as well as the DKE at 300 hPa for ensembles initialized with perturbations of amplitude ϵ0=0.1 , ϵ0/10 , and ϵ0/103 . All curves are averaged over the same initial conditions from the test set as in Fig. 1 .
Figure S2 : Missing butterfly and forecast-backcast asymmetry in the Diffusion-based AIWP models of PlaSim GCM. As in Fig. S1 , but for independently trained forecasting and backcasting AIWP models, based on conditional diffusion models, with Δt=24 h (see Methods and Data). ACC curves are shown for the 16-member ensemble mean (solid) and for a single member (dashed). In the DKE panel, solid curves are obtained in the quasi-deterministic mode used throughout the paper (see Methods and Data), whereas the dashed curve is the spread generated by stochastic sampling alone, with unperturbed initial conditions.
Figure S3 : Forecast-backcast asymmetry across variables and Δt in Transformer-based AIWP models of ERA5 and PlaSim GCM. As in Fig. 1 B, but for U850 (top), Z500 (middle), and U250 (bottom) for ERA5 and PlaSim Transformers ( Δt=24 and 6 h). Colors denote Δt (legend at the bottom); the color-matched annotations show the asymmetry. Curves are averaged over the test set initial conditions (same as in Fig. 1 ).
Figure S4 : RMSE-based forecast and backcast accuracy and asymmetry across variables and Δt in Transformer-based AIWP models of ERA5 and PlaSim GCM. As in Fig. S3 , but showing RMSE rather than ACC for the same models, variables, and initial conditions.
Figure S5 : DKE growth in the numerical and AI models of Lorenz 63. Curves show DKE, computed over all three variables, of forecast ensembles for the numerical model (based on integration of Eqs. ( S15 )–( S17 )) and the AI model. The AI model is trained similarly to the one used for Lorenz 96 (Methods and Data): a multilayer perceptron trained on pairs of all three variables with Δt=10δt . Perturbations are Gaussian, with amplitudes decreasing from a reference value, ϵ0 . Curves are averaged over 10,000 initial conditions drawn from a held-out test trajectory, each with a 100-member ensemble; the same initial conditions and perturbations are used for both models. The AI model closely tracks the numerical model’s DKE curves at all amplitudes. They both lack the “real” butterfly effect [ 21 , 22 ] .
Figure S6 : Scale dependence of the forecast errors in the Pangu-Weather models with different Δt . Power spectra E(n) of the 300-hPa error field (prediction minus ERA5) as a function of total wavenumber n and lead time (colors, nonuniformly from 1 h to 3 days; see legend), compared with the spectrum of the full field (black). The top axis gives the corresponding wavelength. Error saturates first at small scales and progressively fills in the larger scales; at a fixed lead time, the error is larger at all scales for the models trained with smaller Δt . Curves are averaged over the test set initial conditions.
Figure S7 : Spatial distribution of Pangu-Weather forecasts’ early DKE growth and the ERA5 precipitation field. For the rotated figure, rows show lead time from 1 h (first row) to 24 h (last row). The first four columns from the left show the Δt of the official Pangu-Weather model. Panels marked “No data” correspond to lead times shorter than the model’s Δt . As an example, maps of the 300-hPa DKE over parts of the tropical western Pacific are shown for ensembles with the smallest perturbation amplitude ( ϵ0/103 ) for one representative initial condition. The Δt=1 h ensembles show coherent DKE patterns rather than unphysical noise. However, unlike what numerical experiments suggest [ 19 , 33 , 34 ] , the regions of largest early DKE growth in the 1-h Pangu-Weather forecasts do not necessarily collocate with precipitation. Taking the 6-h lead time as an example, the DKE pattern in the 1 -h model instead resembles a mixture of small-scale features amplified in the 6 -h model under the influence of precipitation.
Figure S8 : Spectra of one-time-step forecasts in the AIWP models of ERA5 and PlaSim GCM. Power spectra of the 300-hPa kinetic energy, E(n) , as a function of total wavenumber n , computed with the spherical harmonic transform, for one-time-step ( 1Δt ) forecasts. Left: The four official Pangu-Weather models trained on ERA5 (0.25 ∘ ) from [ 2 ] . Right: Transformer-, SFNO-, and Diffusion-based AIWP models of PlaSim GCM, which has a numerical resolution of T42. Curves are averaged over the test set initial conditions. Bottom row zooms into the high wavenumber end for better visualization. Black curves are the reference spectra of the full ERA5 and PlaSim fields. Vertical dashed lines mark each dataset’s spectral truncation (T639 for ERA5, T42 for PlaSim); gray dotted and dashed lines are −5/3 and −3 slope references.
Figure S9 : One-time-step error of the forecasting and backcasting AIWP models of ERA5 and PlaSim GCM. RMSE of one-time-step prediction ( 1Δt ) for U850 (top), Z500, and U250, evaluated over the test set initial conditions (see Methods and Data). Left: ERA5, which also shows forecast accuracy of the official Pangu-Weather models (gray; Δt=1 , 3, 6, and 24 h [ 2 ] ) and our forecasting (blue) and backcasting (red) Transformers ( Δt=6 and 24 h). Right: the PlaSim forecasting and backcasting Transformers ( Δt=6 and 24 h). Boxes span the interquartile range, whiskers the 5th–95th percentiles, and the horizontal line the median.
Figure S10 : Sensitivity of the Pangu-Weather forecasts’ DKE growth to the averaging domain, GPU precision, and initial-perturbation structure. Each panel shows the ensemble DKE at 300 hPa. Columns correspond to the four official Pangu-Weather models [ 2 ] with Δt=24 , 6, 3, and 1 h. (A) The same as Fig. 4 A, the reference, which uses DKE averaged over 60o S- 60o N, single precision (FP32) inference, and Perlin-noise initial perturbations. (B) as in (A) but averaged globally. (C) as in (A) but without strict FP32 inference, so that the GPU carries out the internal matrix operations in lower precision (TF32). (D) as in (A) but with Gaussian noise instead of Perlin noise. Colors denote the initial-perturbation amplitudes. Gray curves are the ICON reference simulations, solid for 2.5 km and dashed for 20 km resolutions. Panels (B) and (C) show that the faster DKE growth for small-amplitude perturbations might be observed for the small- Δt AIWP models, but they can be due to high-latitude unphysical instabilities or noise from lower-precision GPU calculations. (D) shows little sensitivity to the structure of the initial-condition perturbation.
Weather forecasting predicts future atmospheric conditions, from local storms to global patterns, using physical models, observations, and machine learning.
Could it be that AI weather models are solving physical equations, although they may not be the equations used by conventional NWP models? We compute correlations of forecast skill and Centered Kernel Alignment, providing evidence that different AI weather models represent the atmosphere in similar ways, despite differences in architecture and capacity. We argue that the architecture and training of the AI models constrains the form of the physical laws that they might simulate. In particular, we propose that the models implement a particle description of the atmosphere, where the latent variables at each mesh point correspond to the position of a particle in the high dimensional latent space. We hypothesize that the movement of the particles follows a gradient flow in the latent space towards a minimum of a learned free energy functional. Analysis of the GraphCast and Aurora models show that they make changes on large spatial scales in the early processor layers and move to smaller scale with increasing layer depth, consistent with the gradient flow hypothesis.
George Craig, Tobias Selz, Matthias Beylich +1
Meteorological Institute, LMU Munich · Karlsruhe Institute of Technology
The subseasonal-to-seasonal (S2S) timescale, roughly from two weeks to two months ahead, is a critical forecast window for sectors such as agriculture, energy, and water management. Yet, it is widely known as the `predictability desert'. Recent AI weather models excel up to two weeks ahead but deteriorate beyond, largely because they are trained to predict fine-scale details that are neither predictable nor essential at S2S timescales. We argue that a more physically grounded objective is to forecast only the slowly varying components that remain predictable. Computer vision reached the same conclusion with the Joint-Embedding Predictive Architecture (JEPA), which predicts in latent space, discarding unpredictable details. In this work, we introduce S2S-JEPA, which brings the JEPA paradigm to S2S forecasting. It is tailored to this task through design elements from state-of-the-art AI weather models. S2S-JEPA achieves comparable skill to the gold-standard ECMWF physics-based ensemble and surpasses it on multiple metrics at weeks 5 to 6.
Chenyu Dong, Gianmarco Mengaldo
Department of Mechanical Engineering, College of Design and Engineering, National University of Singapore, 9 Engineering Drive 1, 117575, Singapore.
While AI weather models excel at short-to-medium range forecasts (up to 15 days), they frequently suffer from ill-defined "instabilities" when rolled out over longer horizons. This work addresses the lack of a formal taxonomy by categorizing these failures into three distinct regimes: blow-up, drift, and loss of seasonality, through year-long rollouts of nine state-of-the-art AI weather models. Our analysis reveals that stability hinges on the treatment of small spatio-temporal scales: unstable models amplify high-frequency energy, while stable models act as denoisers when noise is added to their inputs. Far from reducing these models to mere stochastic parrots, our findings highlight that stable models generate unique weather trajectories, conditioned on the initial state. We verify our findings through ablation studies on architectural design choices, conducted using state-of-the-art Vision Transformer (ViT) AI weather model architectures.
Fanny Lehmann, Firat Ozdemir, Yun Cheng +4
ETH AI Center ETH Zurich · Swiss Data Science Center ETH Zurich & EPFL · Scalable Parallel Computing Lab ETH Zurich +3