physics.ao-phAug 26, 2026
SaveMissing the Butterfly and Predicting the Past: Features or Bugs of Accurate AI Weather Models?
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. ·
Abstract
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 advances the atmospheric state (Eq. ( 1 )); an independently trained backcast network 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. 24 h and the horizontal resolution is km (ERA5) and km (PlaSim). Shading marks the ACC skill window; annotations give the crossing times and their forecast/backcast ratio, the asymmetry. See Figs. S1 - S4 for other variables, metrics, , and neural-operator and generative architectures. (C) Difference kinetic energy (DKE) for perturbations of different amplitudes, decreasing from a reference value, . 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, . ( A ) The numerical-model reference trajectory alongside forecasts and backcasts from numerical integrations of Eqs. ( S7 )-( S8 ) initialized at lead time ; 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 and on both scales ( and )) with time interval , , and ( 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, for 12 independent, separate AI models that have an identical architecture (including the same MSE loss) trained on ( A ) only the large-scale and ( B ) on both scales ( and )) with time interval , , and ( 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, . 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 , 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 S- 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 of the 300-hPa meridional wind perturbation (each member minus ensemble mean) as a function of total wavenumber , for the initial-condition perturbations with the smallest ( , upper row) and largest ( , 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 to h; the gray curve is the imposed initial Perlin noise perturbation at , 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 | |
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 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 . For each variable and AI model, the largest value among the three regions is shown in bold. For ERA5, consistent across both 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 | 21.9 6.3 | |
| Ground truth | 24.3 4.0 | |
| , AI model | 26.8 6.5 | |
| , AI model | 26.8 5.9 | |
| , AI model | 3.05 1.84 | |
| , AI model | 2.81 1.59 |
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 , computed over the window – MTU (after discarding a -MTU transient). The FTLE values are projected onto the slow ( ) or fast ( ) variables.
| Category | Variable | Pangu-Weather [ 2 ] | ERA5 Transformer | PlaSim AIWP |
| Temperature ( ) | 13 levels | 17 levels | 13 levels | |
| Zonal wind ( ) | 13 levels | 17 levels | 13 levels | |
| Upper-air | Meridional wind ( ) | 13 levels | 17 levels | 13 levels |
| Specific humidity ( ) | 13 levels | 17 levels | 13 levels | |
| Geopotential ( ) | 13 levels | 17 levels | 13 levels | |
| 2-m temperature |
Table S3 : Variable inventories for the Pangu-Weather, ERA5, and PlaSim AIWP models. Upper-air prognostic state , surface prognostic state , time-varying boundary forcings , and static fields are listed. Checkmarks indicate inclusion; dashes indicate not applicable. The Pangu-Weather column describes the pretrained model on its native grid, which we use for inference only. The ERA5 Transformer column describes the data for the 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 ) and the key observations and conclusions from Fig. 1 B-C remain the same.
| ERA5 Transformer | PlaSim Transformer | PlaSim SFNO | PlaSim Diffusion | ||
| Embedding dim | 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 | — | — | |||
| Patch size | — | ||||
| Optimizer | AdamW | AdamW | AdamW | AdamW |
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 | |
| Depth | hidden layers | |
| Input dimension | 8 ( only) or 264 ( ) | |
| Optimizer | Adam | |
| Hyperparameter search | Learning rate | Log-uniform |
Table S5 : Architectural and training hyperparameters for the Lorenz 96 AI models. Set notation indicates discrete values explored during the hyperparameter optimization phase.
| Forecast RMSE | Backcast RMSE | ||
|---|---|---|---|
| 0.0056 0.0020 | 0.0056 0.0021 | ||
| 0.026 0.0093 | 0.027 0.011 | ||
| 0.048 0.016 | 0.046 0.018 | ||
| 0.038 0.012 | 0.043 0.015 | ||
| 0.16 0.054 | 0.23 0.067 | ||
| 0.29 0.073 | 0.26 0.066 |
Table S6 : One-time-step RMSE of the forecast and backcast AI models. RMSE of a single prediction step of length , computed over the variables against the reference trajectory, for the selected AI models with and . 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 h (see Methods and Data). ACC of , , and versus lead time are shown, as well as the DKE at 300 hPa for ensembles initialized with perturbations of amplitude , , and . 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 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 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 ( and 6 h). Colors denote (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 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 . Perturbations are Gaussian, with amplitudes decreasing from a reference value, . Curves are averaged over 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 . Power spectra of the 300-hPa error field (prediction minus ERA5) as a function of total wavenumber 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 . 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 of the official Pangu-Weather model. Panels marked “No data” correspond to lead times shorter than the model’s . 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 ( ) for one representative initial condition. The 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 -h model instead resembles a mixture of small-scale features amplified in the -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, , as a function of total wavenumber , computed with the spherical harmonic transform, for one-time-step ( ) 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 and 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 ( ) 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; , 3, 6, and 24 h [ 2 ] ) and our forecasting (blue) and backcasting (red) Transformers ( and 24 h). Right: the PlaSim forecasting and backcasting Transformers ( and 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 , 6, 3, and 1 h. (A) The same as Fig. 4 A, the reference, which uses DKE averaged over S- 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- 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.
More in this topic
Weather Forecasting
Weather forecasting predicts future atmospheric conditions, from local storms to global patterns, using physical models, observations, and machine learning.
See all 176 papers