Paper deep dive
Missing the Butterfly and Predicting the Past: Features or Bugs of Accurate AI Weather Models?
Pedram Hassanzadeh, Weidong Li, Y. Qiang Sun, Jiangdi Wang, Alexander Wikner, Justin Finkel, Jonathan Q. Weare
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 8/29/2026, 3:30:23 AM
Summary
This paper investigates the surprising accuracy, missing butterfly effect, and skillful backcasting capabilities of AI Weather Prediction (AIWP) models. Through experiments on ERA5 reanalysis, the PlaSim General Circulation Model (GCM), and the Lorenz 96 system, the authors demonstrate that these phenomena stem from the coarse-graining of training data. Coarse-graining removes fast, small-scale dynamics, allowing AI models to predict large scales accurately without inheriting the rapid error growth (butterfly effect) characteristic of physical systems. Reducing coarse-graining makes AI models more physically consistent (restoring the butterfly effect and making backcasting impossible) but reduces forecast accuracy.
Entities (8)
Relation Signals (9)
Coarse-graining → causes → Missing Butterfly Effect
confidence 95% · We trace the surprising forecast accuracy, missing butterfly, and skillful backcasting to a single cause: inevitable coarse-graining of training data
Coarse-graining → causes → High Forecast Accuracy
confidence 95% · Results offer an explanation for AIWP models' forecast skill... they implicitly learn how fast, small scales affect large scales without inheriting their rapid error growth.
AI Weather Prediction (AIWP) → exhibits → Missing Butterfly Effect
confidence 95% · all these forecasting and backcasting models miss the butterfly effect.
Pangu-Weather → isinstanceof → AI Weather Prediction (AIWP)
confidence 95% · From the Lorenz system to official Pangu-Weather models...
Coarse-graining → causes → Backcasting
confidence 90% · We trace the surprising forecast accuracy, missing butterfly, and skillful backcasting to a single cause: inevitable coarse-graining of training data
AI Weather Prediction (AIWP) → exhibits → Backcasting
confidence 90% · AI models can be trained to skillfully predict the past (backcast)
Lorenz-96 → isusedtosimulate → AI Weather Prediction (AIWP)
confidence 90% · Here, across a hierarchy spanning... the multi-scale Lorenz system, we show that AI models can be trained...
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
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.
Tags
Links
- Source: https://arxiv.org/abs/2608.25835v1
- Canonical: https://arxiv.org/abs/2608.25835v1
Trouble viewing inline? Open PDF directly →
Full Text
128,389 characters extracted from source content.
Expand or collapse full text
Missing the Butterfly and Predicting the Past: Features or Bugs of Accurate AI Weather Models? Pedram Hassanzadeh Affiliation: Department of Geophysical Sciences, University of Chicago, Chicago, USA Affiliation: Committee on Computational and Applied Mathematics, University of Chicago, Chicago, USA. Weidong Li Y. Qiang Sun Affiliation: State Key Laboratory of Severe Weather Meteorological Science and Technology, Nanjing University, Nanjing, China. Affiliation: Key Laboratory of Mesoscale Severe Weather/Ministry of Education, Nanjing University, Nanjing, China. Jiangdi Wang Affiliation: Committee on Computational and Applied Mathematics, University of Chicago, Chicago, USA. Alexander Wikner Affiliation: Department of Geophysical Sciences, University of Chicago, Chicago, USA Affiliation: Data Science Institute, University of Chicago, Chicago, USA. Justin Finkel Affiliation: Data Science Institute, University of Chicago, Chicago, USA. Jonathan Q. Weare Affiliation: Courant Institute School of Mathematics, Computing, and Data Science, New York University, New York City, USA. ∗Corresponding authors: pedramh@uchicago.edu, qiangsun@nju.edu.cn†These authors contributed equally to this work. 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. 1 Introduction Artificial intelligence (AI) weather prediction (AIWP) models have transformed short- to medium-range (days to two weeks) forecasting [1, 2, 3]: they outperform the best physics-based numerical weather prediction (NWP) models across a wide range of common skill metrics at orders of magnitude lower computational cost [4, 5, 6, 7]. State-of-the-art AIWP models are trained in a broadly similar fashion: to advance the atmosphere’s 3D global state x from time t to t+Δt+ t using a deep neural network N, i.e., (t+Δt)=θ((t)),x(t+ t)=N_θ (x(t) ), (1) where Δt t is often 6-24 hours, θ represents the learnable parameters, and the training data come from ∼40 40 years of observation-based reanalysis products like ERA5 [8]. Longer forecasts are produced autoregressively. Given their accuracy and computational efficiency, AIWP models are now being operationalized by leading forecasting centers [9, 10] and increasingly used to provide actionable weather information to those most in need [11]. More recently, these models have been pushed to subseasonal-to-seasonal time scales and even to longer-term climate projections, with promising early results [12, 13]. Yet, a growing number of studies have identified several failure modes of AIWP models, e.g., their inability to forecast unprecedented “gray swan” weather extremes [14, 15, 16] and to reproduce the small-scales’ statistics [17, 18]. Most surprisingly, these highly skillful models were found to lack the butterfly effect [19]: the rapid growth of ensemble spread, which accelerates as initial-condition perturbations vanish, in multi-scale chaotic systems [20, 21, 22]. This deficiency appears shared across otherwise very different deterministic and generative AI architectures [23]. While this failure has generated concerns, its root cause and implications remain unclear. Meanwhile, the sources of AIWP models’ surprising accuracy are also poorly understood, hindering trust and further improvements. Indeed, this accuracy comes despite long-standing skepticism about skillful data-driven weather forecasting, e.g., from Lorenz [24] and others, including an estimated requirement of ∼1030 \!10^30 years of training data [25]; see Discussion. Underlying these puzzles, fundamental questions about the chaotic atmosphere’s predictability limits, and the butterfly effect’s practical relevance, remain a subject of extensive research and debate [26, 27, 28, 29, 30]. Here, we investigate the missing butterfly effect and surprising forecast accuracy of AIWP models. First, since nothing in these AIWP models’ design restricts them to the arrow of time, we also train “backcasting” models: neural networks that predict the past (Fig. 1A), i.e., (t−Δt)=ϕ((t)),x(t- t)=N_φ (x(t) ), (2) the Eq. (1) “forecast” model’s time-reversed counterpart, with its own (independent) learnable parameters ϕφ. As shown below, autoregressive backcasting offers a new, useful dimension for exploring AIWP models and fundamentals of atmospheric predictability. Across a hierarchy spanning ERA5 reanalysis, an intermediate-complexity general circulation model (GCM), and multi-scale Lorenz 96, we show that, while skillful, these backcasting AI models consistently behave in ways sharply at odds with fundamental understanding of dynamical systems and their numerical integrations. Simply put, skillful backcasting appears to violate the second law of thermodynamics. We also show that across the hierarchy, the butterfly effect is missing from these AI forecasting and backcasting models alike. We then explicitly connect all three (surprising AI forecast skills, the missing butterfly effect, skillful AI backcasting), showing they are symptoms of the same underlying cause. We examine whether each symptom is a feature or a bug of skillful AIWP models. 2 Results For each system in our hierarchy, the backcasting AI models (Eq. (2)) are architecturally identical to the forecasting models (Eq. (1)) and are trained on the same data, optimizer, and loss function (often the mean absolute or squared error, MAE/MSE); the sole difference is the reversed input–output pairing (Fig. 1A). We have trained a number of deterministic and generative AIWP models based on state-of-the-art architectures (see Methods and Data). 2.1 Skillful AI backcasting and persistent asymmetry: ERA5 and PlaSim (GCM) Figure 1B shows the anomaly correlation coefficient (ACC) of 500-hPa geopotential (Z500) as a function of lead time. On ERA5, the backcasting Transformer (deterministic) remains skillful to 6.5 days into the past while the forecast is skillful for 9.1 days. It might be surprising that the past is harder to predict (asymmetry =9.1/6.5>1=9.1/6.5>1). However, backcasting should be practically impossible for dissipative dynamical systems such as the atmosphere, because, simply put, it leads to an apparent violation of the second law of thermodynamics! As further explained in the Supplementary Materials, time reversal turns dissipation (which is prevalent in the atmosphere across scales) into anti-dissipation (and reverses the sign of Lyapunov exponents), so unavoidable small errors grow explosively, fastest at the smallest scales. Accordingly, backcasting using numerical integrations fails catastrophically; see examples from the heat equation to Lorenz 63 to Lorenz 96 in the Supplementary Materials and Fig. 2A. This is not the case for the AIWP model of ERA5. The asymmetry >1>1 is not an artifact of ERA5, which, as an assimilation product continually updated by observations, is not a closed dynamical system: PlaSim, which is one (a GCM), exhibits even larger asymmetries, including in generative and neural-operator AIWP models (Figs. 1B and S1-S2). The asymmetry is robust across variables, pressure levels, and skill metrics in both ERA5 and PlaSim (Figs. S3-S4). Geographically, the asymmetry is overall larger in the tropics than in the extratropics (Table S1); likely a reflection of the stronger diabatic and dissipative processes (higher entropy production rate) there [31, 32], so that the tropical atmosphere is, in this sense, further from reversibility (see the Supplementary Materials). Two robust findings therefore need explaining: that backcasting works at all, and that the past is nevertheless consistently harder to predict than the future. But first, let us also examine the error growth and butterfly effect in these forecasting and backcasting AIWP models. 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 θN_θ advances the atmospheric state (Eq. (1)); an independently trained backcast network ϕN_φ 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= t= 24 h and the horizontal resolution is ∼100 100 km (ERA5) and ∼300 300 km (PlaSim). Shading marks the ACC >0.6>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 t, and neural-operator and generative architectures. (C) Difference kinetic energy (DKE) for perturbations of different amplitudes, decreasing from a reference value, ϵ0 _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. 2.2 Missing butterflies in AI forecasting and backcasting alike: ERA5 and PlaSim (GCM) In a multi-scale chaotic system like the atmosphere, ensemble spread from vanishingly small initial perturbations should grow rapidly, fed by upscale error transfer from the small scales, and saturate after a finite time. This is the “real” butterfly effect [20, 21, 22, 28], which sets the upper limit of atmospheric predictability. We quantify the spread with the difference kinetic energy (DKE), the area-weighted ensemble variance of the horizontal winds (Methods and Data). Figure 1C shows that numerical models, ICON (NWP) and the PlaSim GCM, reproduce this behavior; however, no AIWP model does, in forecasting or backcasting, for ERA5 or PlaSim. The growth of DKE in the numerical integrations of physics-based forecast models behaves as expected: in both ICON and PlaSim GCM, the smallest-amplitude ensembles grow by orders of magnitude within a day or two, showing diminishing return as the initial-condition errors shrink (Fig. 1C). This is the butterfly effect [21, 22]. For the large-amplitude perturbations, ensemble forecasts from numerical integrations and AIWP models behave similarly, showing exponential growth in DKE over several days (unlike the butterfly effect, the growth rate here is independent of the perturbation amplitude). The forecasting results in ERA5 agree with those of Selz and Craig in their pioneering papers [19, 23]. Here, we show the same behavior in the AIWP models of a GCM, but most strikingly, in backcasting AIWP models. The fundamental understanding of time-reversal physics (Supplementary Materials) and numerical experiments with Lorenz 96 (Fig. 3A) suggest even a more explosive growth of DKE in backcasting. However, we robustly see a near symmetry between the forecasting and backcasting DKE growth across the hierarchy. We are thus left with three puzzles about AIWP models: on one hand, forecast accuracy that long-standing arguments deemed unattainable, and on the other, skillful backcasting and a missing butterfly effect that disagree with the fundamental understanding of multi-scale chaotic dynamics, robustly confirmed by numerical integrations of physics-based models. To further probe these findings, we extend the hierarchy to a canonical testbed of multi-scale, chaotic dynamics, the Lorenz 96 system. 2.3 Training-data coarse-graining dramatically alters forecasts, backcasts, and butterflies: Lorenz 96 The two-scale Lorenz 96 system (Eqs. (S7)-(S8)) couples 8 slow, large-amplitude variables XiX_i (representing large-scale circulation) to 256 fast, small-amplitude variables Yi,jY_i,j whose characteristic time scale is an order of magnitude shorter (representing small-scale, often subgrid-scale processes); see Methods and Data. Its numerical integration displays the time irreversibility that should make any backcasting impractical: Fig. 2A shows that integrated backward, the trajectory blows up rapidly due to anti-dissipation. As demonstrated in Fig. 3, the numerical ensemble forecasts of this system exhibit the butterfly effect, and the overall DKE growth curves show a behavior similar to those from numerical integration of physics-based models like ICON and PlaSim GCM as the perturbation amplitudes are varied (compare to Fig. 1C). The AI forecasting and backcasting models can exhibit very different behaviors depending on the degree of coarse-graining of the training data: How many spatio-temporal scales the training data includes. We first train AI models with data closest to the “real-world” regime: only using the large scales and after removing the fastest variabilities (X-only, Δt=10δt t=10δ t where δtδ t is the time step of the numerical solver/ground truth). As shown in the top-left corners of Figs. 2B-C and 3, the forecasts and backcasts of these AI models have remarkably similar characteristics as those of AIWP models of ERA5 and PlaSim (Figs. 1B-C and S3-S4): Accurate forecasts, skillful backcasts with asymmetry >1>1, amplitude-independent exponential growth in DKE, and missing butterflies. Hovmöller diagrams and DKE curves of backcasts demonstrate behavior that is markedly different from that of numerical integrations. Note that the missing butterfly effect in the real-world-regime AI model of two-scale Lorenz 96 and in the AIWP models of ERA5 and PlaSim (Figs. 1C and 3) does not mean that these models miss chaos. Their predictions are sensitive to initial-condition perturbations across perturbation amplitudes. In fact, their amplitude-independent exponential DKE growth closely resembles that of numerical integrations of Lorenz 63 (Fig. S5), the hallmark example of deterministic chaos, but itself misses the real butterfly effect because it lacks multi-scale dynamics (see Palmer [21, 22] for an illustrative discussion of Lorenz 63 versus Lorenz 96). 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)X_i\,(i=1… 8). (A) The numerical-model reference trajectory alongside forecasts and backcasts from numerical integrations of Eqs. (S7)-(S8) initialized at lead time 00; 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 XiX_i and on both scales (XiX_i and Yi,j(j=1…32CLOSEY_i,j\,(j=1… 32)) with time interval Δt=10δt t=10δ t, 5δt5δ t, and 1δt1δ t (δ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. Next, we develop additional (independent) forecasting and backcasting AI models trained on data with more spatial and/or temporal scales: using X or (X,Y)(X,Y) and Δt=1 t=1, 5, or 10 δtδ t (see Methods and Data for details of extensive model design exploration). At the bottom-right corner of Figs. 2B-C and 3, we show the results of the AI models trained in the “perfect-data” regime: these models have seen all the scales the numerical solution/ground truth contains. Compared with the models trained in the “real-world” regime, the AI models trained in the perfect-data regime have a worse forecast skill (by a factor of ∼7 7) and larger asymmetry (by a factor of ∼4 4) but also blown-up backcasting and rapid DKE growth that closely mimics the butterfly effect. The last two in particular show that the perfect-data-regime AI models are consistent with the fundamental understanding of multi-scale chaotic systems and their numerical integrations. We observe these trends in forecast skill, asymmetry, DKE growth, and Hovmöller diagrams consistently as we move from the “real-world” regime corner to the “perfect-data” regime corner. Therefore, we propose the following three hypotheses supported by these Lorenz 96 experiments: coarse-graining used in producing the training data is the common cause of H1) Unexpected AI forecast skill (for large scales), H2) Missing butterfly effect, and H3) Skillful backcasting. We broadly define coarse-graining, which is inevitable in practice, as filtering of small/fast scales and/or ignoring some state variables. The emergence of H1 from coarse-graining can be further seen in the analysis of leading finite-time Lyapunov exponents of the large and small scales, λXλ^X and λYλ^Y (Table S2). In the ground truth, λXλ^X and λYλ^Y are the same (within the uncertainty range), reflecting the fact that the error growth in the large scales is determined by the small scales. The perfect-data-regime AI model reproduces these Lyapunov exponents fairly well. However, the real-world-regime AI model, which lacks small scales, has λXλ^X around 8 times smaller than the ground truth/perfect-data-regime AI model, leading to its enhanced forecast skills. We note that here, H1 refers to the skill for large scales, which is relevant for practical forecasting and the basis of common verification metrics (see Discussion). As for H2, Fig. 3 shows the dramatic impact of training-data coarse-graining on the DKE growth of the 12 AI forecasting and backcasting models, which share an identical architecture and MSE loss. In particular, when trained on complete or nearly complete data, these MSE-trained models exhibit a butterfly-like effect: rapid, amplitude-dependent DKE growth comparable to, or even faster than, that of the numerical integrations (Fig. 3B). This argues against the suggestion that optimizing toward a conditional mean is the primary cause of the missing butterfly effect (see the Discussion). These results raise the question of whether a butterfly-like effect can likewise emerge in a state-of-the-art AIWP model when the coarse-graining of its training data is reduced, which we address next. 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)X_i\,(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 XiX_i and (B) on both scales (XiX_i and Yi,j(j=1…32CLOSEY_i,j\,(j=1… 32)) with time interval Δt=10δt t=10δ t, 5δt5δ t, and 1δt1δ t (δ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 _0. See Methods and Data for details about the perturbations and metrics. 2.4 As coarse-graining is reduced, a butterfly-like effect emerges and forecast skill degrades: Pangu-Weather We use the official Pangu-Weather models [2]: four independent deterministic neural networks with time steps Δt=24 t=24, 66, 33, and 11 h, all trained with an MAE loss on 0.25∘0.25 ERA5, the native spatial resolution of the dataset. We only perform inference (forecasts) with these models; we have not trained, fine-tuned, or otherwise modified them. Figure 4A shows that as Δt t decreases, the small-amplitude ensembles transition from slow, amplitude-independent exponential DKE growth for the 2424-h Pangu-Weather, consistent with [19, 23] and our 1∘1 Transformer (Fig. 1C), to rapid, amplitude-dependent growth for the 11-h Pangu-Weather that approaches the butterfly-like behavior of the ICON integrations. Meanwhile, the skillful lead time drops from 9.3 to 2.0 days. Consistent with the Lorenz 96 experiments (Figs. 2-3) and hypotheses H1 and H2, reducing (here, temporal) coarse-graining degrades forecast skill and strengthens butterfly-like growth. This behavior has not been reported before because the Δt=1 t=1 h and 33 h Pangu-Weather models had not been run continuously. Selz and Craig [19] used all four networks, but through hierarchical temporal aggregation, in which the Δt=24 t=24 h model takes the longest steps and the shorter-Δt t models only fill the intermediate hours; the drops in DKE at network switches that they reported, largest when the 24-h model takes over, are consistent with Fig. 4A. Their follow-up study [23] ran only the 6- and 24-h versions. The spectra of the ensemble perturbations (Fig. 4B) further characterize the emerging DKE growth. For Δt=24 t=24 and 66 h, the small-amplitude spectra jump within the first step and then grow slowly, with their shape roughly preserved, remaining orders of magnitude below the ERA5 spectrum even at 3 days. For Δt=1 t=1 h, perturbation energy instead grows fastest at the small scales, saturates there first, and progressively fills in the larger scales toward the ERA5 spectrum: the upscale error growth expected in multi-scale chaotic systems [20, 21]; see Fig. S6. However, at the smallest scales, the late-time spread overshoots the ERA5 spectrum, beyond the factor of 2 allowed by physical saturation; the 1-h model generates variance at scales that ERA5 itself barely contains, a warning sign discussed in Methods and Data. We emphasize that we are not claiming that this is the (real) butterfly effect. Even at its native resolution, ERA5 (and especially its subset of variables used for training) is highly coarse-grained: the reanalysis is produced by a ∼ 31-km NWP model with parameterized subgrid processes, it misses state variables, and Δt=1 t=1 h remains far longer than the time step of the NWP model’s numerical solver. The Lorenz 96 hierarchy shows what to expect in this situation: with small spatial scales fully missing (X-only), AI models overall remain in the real-world-regime behavior even when trained at the solver’s own time step (Fig. 3, top-right corner, where Δt=1δt t=1δ t). Even in the perfect-data regime, we call the behavior in Fig. 3 butterfly-like: the DKE grows faster than the numerical model’s near t=0t=0, whereas the real-world-regime AI models match the early numerical growth rates better but never accelerate as time progresses forward or backward. Returning to the Pangu-Weather model with Δt=1 t=1 h, unlike what numerical experiments show [33, 34, 19], the largest early DKE growth is not necessarily collocated with precipitation (Fig. S7). However, we offer a word of caution and argue that assessing the butterfly effect, and more broadly, the physicality of error growth, in AI models requires deeper investigation and a revisiting of what these terms mean, as the current understanding heavily relies on numerical experiments, with their own known shortcomings. Disagreement with a numerical model does not, by itself, make an AIWP model wrong: these data-driven models might be learning, and forecasting, weather differently. One must also be cautious about hardware artifacts [23] and nonphysical forecasts while examining butterflies in AI models; see Methods and Data for how these are controlled for and assessed. What is robust here is the trade-off: within a single family of state-of-the-art AIWP models, reducing temporal coarse-graining results in butterfly-like growth and degrades forecast skill, further supporting H1 and H2. 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 t=24, 6, 3, and 1 h. (A) DKE versus lead time, as in Fig. 1C, for three amplitudes of initial perturbations. DKE is averaged over 60o60^oS-60o60^oN. 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. 1B. (B) Power spectra Sv(n)S_v(n) of the 300-hPa meridional wind perturbation v300′v _300 (each member minus ensemble mean) as a function of total wavenumber n, for the initial-condition perturbations with the smallest (ϵ0/103 _0/10^3, upper row) and largest (ϵ0 _0, lower row) amplitudes. These curves are averaged over the same ensemble members and initial conditions as the ERA5 Transformer in Fig. 1C. Colors denote (nonuniform) lead times from +1+1 to +72+72 h; the gray curve is the imposed initial Perlin noise perturbation at t=0t=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. 3 Discussion Across the hierarchy of ERA5, PlaSim GCM, and Lorenz 96, we find that the coarse-graining of the training data is the common cause of the following (seemingly unrelated) key features of AI prediction models: H1. The surprising forecast accuracy, H2. The missing butterfly effect, H3. The skillful backcasting. In short, the more completely the training data represent the true dynamical system spatially and temporally, the more the AI model behaves like the system, but the worse it forecasts the large scales, which are often the desired target of prediction and most verification metrics. Before discussing H1-H3 and their implications, we note that these results offer a step toward answering a central question about AIWP models: Do they learn atmospheric dynamics, or merely learn the evolution of weather the way one learns a video? A video plays backward as easily as forward. Skillful backcasting (H3), taken alone, might appear to support this video interpretation, but only in the strongly coarse-grained regime characteristic of the AIWP models’ common training datasets. As the training data become more complete, the same architecture, trained with the same loss on pairs one time step apart, becomes markedly more physics-like: backcasting skill disappears, the forecast–backcast asymmetry grows, and a butterfly-like effect emerges. We are not suggesting that AIWP models trained on strongly coarse-grained data (e.g., ERA5 with Δt=24 t=24 h) learn weather much like a video: dynamical tests and other analyses have demonstrated that they learn aspects of atmospheric dynamics [14, 15, 35, 36, 37, 38]. Our findings, however, show that with more complete training data, these models behave more like the physical system, including its arrow of time, although they become less accurate for practical forecasting. A major question then is whether there is a downside to having AIWP models that are more accurate but less physics-like; we return to this later. Regarding H1, coarse-graining removes the fast, small scales through which errors grow upscale to the large scales [20, 28, 33, 34]. Therefore, the effective system that the AI models learn has much slower large-scale error growth than the true system (quantified for Lorenz 96 via Lyapunov exponents; Table S2). That lower spatio-temporal resolution improves forecast accuracy might seem counterintuitive, but this intuition is built on numerical simulations, in which higher resolution improves accuracy. This is because of lower discretization errors and, more importantly, reduced reliance on “explicit” subgrid-scale parameterization schemes, the leading source of structural/epistemic uncertainty in these models. AI models, in contrast, learn subgrid-scale parameterizations “implicitly”: they can account for the averaged effects of the fast, small scales on the large scales’ evolution without explicitly representing the former, which can otherwise drive rapid error growth. Consistent with this implicit parameterization, AI climate emulators, trained on coarse-grained outputs of high-resolution physics-based simulations can closely reproduce key climate statistics of the original simulations [39]. Moreover, even in NWP models, increasing resolution yields diminishing returns in large-scale forecast accuracy [40], in part because the newly resolved fast, small scales drive rapid error growth [29, 34, 41]. Taken together, we argue that this implicit parameterization in AIWP trained with coarse-grained data significantly enhances their forecast skill and is a major reason for their surprising accuracy. We should clarify that the surprising accuracy is not necessarily about comparison with the NWP models. Such a comparison would have to account for the quality of the physics-based models’ subgrid-scale parameterization schemes, a complex question. Rather, H1 concerns the AIWP models’ skill, which is unexpected given earlier studies arguing that data-driven weather forecasting requires an enormous amount of training data [24], with 103010^30 years as one estimate [25]. Our findings suggest that such estimates need to be revisited in light of the effects of AIWP models’ training-data coarse-graining on error growth and the butterfly effect, key contributors to those high data requirements. Other factors likely contribute to the gap between these estimates and practice, including the underlying assumptions, e.g., about the attractor dimension [42, 43], and, perhaps more importantly, the learning algorithm. Analog forecasting in those studies is essentially nearest-neighbor regression, whose scaling suffers from the curse of dimensionality [44, 45]. Current AIWP models instead use deep neural networks (over-parameterized nonlinear regression) to learn a smooth Δt t flow map [46, 47, 48]. Furthermore, AIWP models have been shown to “translocate”, i.e., learn from dynamically similar events in other regions, which significantly enrich their training set [15, 14]. Developing theoretical scaling laws relating training-data size to AIWP model accuracy, beyond emerging empirical estimates [49, 50], remains an important open research direction and can build on our findings. As for the missing butterfly effect (H2), it is a symptom of coarse-grained training data, not of the architecture, the loss function, or deterministic-versus-generative design. The Lorenz 96 experiments isolate the cause cleanly (Fig. 3): all 12 networks share one architecture and one MSE loss, yet the butterfly effect appears or disappears with the content of the training data alone. The Pangu-Weather experiments with varying Δt t (Fig. 4) show the same dependence in a state-of-the-art AIWP model. Spatially, however, ERA5 (and PlaSim) are already highly coarse-grained at their native resolutions, leaving no room to informatively further vary the spatial scales in either direction; future explorations should focus on high-resolution global NWP models and GCMs. These findings argue against optimizing toward a conditional mean or blurring, which is due in part to spectral bias [17], as the primary cause: MSE-trained networks reproduce the butterfly effect in Lorenz 96, while the Diffusion AIWP models of PlaSim studied here and GenCast and FourCastNet2 examined in [23] miss butterflies despite their little spectral bias (Figs. 1C, S2, and S8). Note that after showing the lack of butterfly effect across AIWP models of different architectures and loss functions, [23] suggested a training-data origin, the remaining common factor (although all these AIWP models also share the same learning principle: supervised prediction of the short-Δt t evolution). However, the problem has been attributed to the specific nature of the data-assimilated analysis products used for training and fine-tuning [23, 22]. The hierarchy here, whose PlaSim and Lorenz 96 training sets involve no data assimilation, tests the training-data origin directly and identifies the key property: its spatial and temporal coarse-graining. Our results also suggest that restoring butterflies can come at the cost of forecast accuracy; still, if one views their absence as a bug rather than a feature of accurate AIWP models, the question becomes how best to restore them (see below). Skillful backcasting (H3) is made possible by the same coarse-graining. The learned effective system is smoother, with smaller leading Lyapunov exponents and less of the fast, small scales that dominate the anti-dissipative directions responsible for explosive backward integration (see Supplementary Materials). The apparent violation of the second law of thermodynamics is thereby resolved: the learned effective system is less irreversible than the true one, but still irreversible. This residual irreversibility is manifested in the asymmetry remaining above one, which is largest in the tropics, consistent with the higher entropy production there [31, 32], and in preliminary explorations in which long-term AI backcasting of Lorenz 96, PlaSim, and ERA5 leads to unstable or unphysical states, a sign of irreversibility, though much slower than the true system. This work introduces backcasting as a new, useful lens for investigating the fundamentals of predictability and better understanding AIWP models. Backcasting can also provide practical tools of its own. Speculatively, a single AI model trained both forward and backward might internalize more of the underlying dynamics, perhaps even rudiments of causality. That said, our preliminary bidirectional forecasting-backcasting model (Methods and Data) shows no improvement so far, but the design space is essentially unexplored. Backcast models might also help sample the rarest weather extremes and find their dynamical precursors, tasks traditionally addressed with limited linear tools such as adjoints. AI-based approaches have recently emerged, including AI-boosted rare-event sampling [51], iterative use of AIWP models’ adjoints [52], and sampling emulators like cBottle [53, 54]. Backcasting could augment these as a stable, fully nonlinear, autoregressive counterpart of the adjoints. Growing efforts focus on identifying the physics that AI models miss [14, 16, 17, 18, 19, 23, 55] and on making them more physics-like, e.g., through inductive biases such as constraints in architectures and loss functions [56, 57, 58, 59]. Our results prompt a question: do we actually want AIWP models to fully behave like physics-based models even if restoring the physics costs performance, e.g., accuracy of practical forecasting (Figs. 2 and 4)? This trade-off may prove solvable, and there is precedent for excluding physics for the benefit of practical skill: sound waves, computationally costly to resolve yet largely irrelevant to weather, have been filtered out of the equations solved by many NWP models and GCMs [60]. However, the distinction is that while sound waves carry little of what matters for weather and climate prediction, the fast, small scales filtered by coarse-graining carry much of it. Thus, the deeper question is what these AI models are missing because coarse-graining has removed physics from their training data. Is this absence part of why AIWP models struggle with gray swans [14, 16]? What does it imply for probabilistic forecasts, especially at subseasonal-to-seasonal lead times, which rely on perturbation growth and must be well calibrated [11, 61]? The question might be even more pressing for long-term climate emulators: AI models of the emerging global km-scale, physics-based simulations [62] are typically trained on output first heavily coarse-grained, e.g., to ∼ 100 km [39], to Δt t of hours to a day, and even to daily-averaged state variables [13], which can discard significant physics. Consistent with this concern, idealized data-driven models with partial state vectors (missing variables) fail to reproduce forced responses [63]. Such failures offer a caution for emulators built for sampling internal variability and projecting climate change [56]. Answering these questions requires a deeper understanding of what AI weather and climate models learn, do not learn, and why. To conclude, the development of AI weather (and climate) models has concentrated on architectures and objectives: deeper and larger networks, customized loss functions, and embedded physical constraints. Those directions address the optimization and expressiveness of AI models. Our results point instead to the ingredient that no architecture or loss function can substitute for: the information content of the training data. Two remedies might come to mind immediately: hybrid AI-physics modeling, and foundation modeling, in which heterogeneous datasets are combined through large-scale pretraining. However, state-of-the-art examples of each, NeuralGCM [5] and Aurora [64], miss the butterfly effect [23]; restoring at least this feature appears to require other approaches. Coarse-graining does entail an irreversible loss of information, yet some of what is lost can, in principle, be represented rather than resolved, as memory and stochasticity: the Mori–Zwanzig formalism [65, 66] and its data-driven descendants [67, 68, 69] (see Supplementary Materials) might provide a principled starting point. For predicting the large scales of weather days ahead, coarse-graining is a feature, the source of AIWP skill and of learnable backcasting; for representing the physics of the atmosphere, its error growth, irreversibility, and predictability limits, it is a bug, likely most consequentially for long-term emulation. The content of the training data is therefore not a preprocessing detail but the central design decision of an AI weather or climate model. Acknowledgments and Disclosure of Funding We thank Kerry Emanuel and Fabrizio Falasca for insightful discussions about atmospheric predictability and time reversal, and Boris Bonev and Daniel Boscu for helpful comments on GPU precision. Computational resources were provided by NSF ACCESS (allocation ATM170020), NCAR’s CISL (allocations URIC0009 and UCHI0018), and the University of Chicago’s Research Computing Center and Data Science Institute (DSI). Claude Fable 5 has been used for editing and proofreading the text, and for generating the visualization code for Movie S1. Funding: This work was supported by NSF grant AGS-2531264, Schmidt Sciences LLC (through the InMOS project), and the University of Chicago’s Data Science Institute (DSI) and the Institute for Climate and Sustainable Growth (through the AI for Climate Initiative) to PH. AW and JF are grateful to DSI for an Eric and Wendy Schmidt AI in Science and an AI for Climate postdoctoral fellowship, respectively. JQW was supported by a Pritzker AI+Science Visiting Scholarship from DSI. Author contributions: PH conceived the idea and wrote the paper. PH and YQS supervised research. PH, YQS, and JQW designed the study. WL and YQS (ERA5 and PlaSim) and JW (Lorenz systems) trained AI models, analyzed data, generated the figures, and wrote the Methods and Data section. AW developed AI models and generated the PlaSim data. All authors interpreted the results and edited the paper. Competing interests: There are no competing interests to declare. Data and code availability: We downloaded the ERA5 dataset from Copernicus https://cds.climate.copernicus.eu/ and regridded it with CDO https://code.mpimet.mpg.de/projects/cdo/. The ICON data are provided by [19]. The PlaSim GCM is available at https://github.com/HartmutBorth/PLASIM. The official Pangu-Weather models are obtained from https://github.com/198808xc/Pangu-Weather. Movie S1 is available at https://doi.org/10.5281/zenodo.22062019. We will make all trained AI models and their inference data public upon publication. References [1] J. Pathak, S. Subramanian, P. Harrington, S. Raja, A. Chattopadhyay, M. Mardani, T. Kurth, D. Hall, Z. Li, K. Azizzadenesheli, P. Hassanzadeh, K. Kashinath, and A. Anandkumar (2022) FourCastNet: a global data-driven high-resolution weather model using adaptive Fourier neural operators. arXiv preprint arXiv:2202.11214. External Links: Document Cited by: §B.2, §1. [2] K. Bi, L. Xie, H. Zhang, X. Chen, X. Gu, and Q. Tian (2023) Accurate medium-range global weather forecasting with 3d neural networks. Nature 619 (7970), p. 533–538. External Links: Document Cited by: §A.1, §A.1, §A.1, Figure S10, Figure S8, Figure S9, §B.2, Table S3, §1, Figure 1, Figure 1, Figure 4, Figure 4, §2.4. [3] R. Lam, A. Sanchez-Gonzalez, M. Willson, P. Wirnsberger, M. Fortunato, F. Alet, S. Ravuri, T. Ewalds, Z. Eaton-Rosen, W. Hu, et al. (2023) Learning skillful medium-range global weather forecasting. Science 382 (6677), p. 1416–1421. Cited by: §1. [4] I. Price, A. Sanchez-Gonzalez, F. Alet, T. R. Andersson, A. El-Kadi, D. Masters, T. Ewalds, J. Stott, S. Mohamed, P. Battaglia, et al. (2025) Probabilistic weather forecasting with machine learning. Nature 637 (8044), p. 84–90. Cited by: §1, Figure 1, Figure 1. [5] D. Kochkov, J. Yuval, I. Langmore, P. Norgaard, J. Smith, G. Mooers, M. Klöwer, J. Lottes, S. Rasp, P. Düben, et al. (2024) Neural general circulation models for weather and climate. Nature 632 (8027), p. 1–7. Cited by: §1, §3. [6] S. Rasp, S. Hoyer, A. Merose, I. Langmore, P. Battaglia, T. Russell, A. Sanchez-Gonzalez, V. Yang, R. Carver, S. Agrawal, et al. (2024) WeatherBench 2: a benchmark for the next generation of data-driven global weather models. Journal of Advances in Modeling Earth Systems 16 (6), p. e2023MS004019. Cited by: §1. [7] R. Masiwal, C. Aitken, A. Marchakitus, M. Gupta, K. Kowal, H. A. Pahlavan, T. Yang, Y. Q. Sun, M. Kremer, A. Jina, et al. (2026) Decision-oriented benchmarking to transform ai weather forecast access: application to the indian monsoon. arXiv preprint arXiv:2602.03767. Cited by: §1. [8] H. Hersbach, B. Bell, P. Berrisford, S. Hirahara, A. Horányi, J. Muñoz-Sabater, J. Nicolas, C. Peubey, R. Radu, D. Schepers, A. Simmons, C. Soci, S. Abdalla, X. Abellan, G. Balsamo, P. Bechtold, G. Biavati, J. Bidlot, M. Bonavita, G. D. Chiara, P. Dahlgren, D. Dee, M. Diamantakis, R. Dragani, J. Flemming, R. Forbes, M. Fuentes, A. Geer, L. Haimberger, S. Healy, R. J. Hogan, E. Hólm, M. Janisková, S. Keeley, P. Laloyaux, P. Lopez, C. Lupu, G. Radnoti, P. de Rosnay, I. Rozum, F. Vamborg, S. Villaume, and J. Thépaut (2020) The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological Society 146 (730), p. 1999–2049 (en). External Links: ISSN 1477-870X Cited by: §A.1, §1. [9] E. G. Daub, T. Dunstan, T. Bennett, M. Burnand, J. Chappell, A. Coca-Castro, N. Eftekhari, J. S. Hosking, M. Janmaijaya, J. Lillis, et al. (2025) Technical overview and architecture of the fastnet machine learning weather prediction model, version 1.0. arXiv preprint arXiv:2509.17658. Cited by: §1. [10] G. Moldovan, E. Pinnington, A. Prieto Nemesio, S. Lang, Z. Ben Bouallègue, J. Dramsch, M. Alexe, M. Santa Cruz, S. Hahner, H. Cook, et al. (2026) AIFS single 1.1. 0: an update to ecmwf’s machine-learned weather forecast model aifs. Geoscientific Model Development 19 (10), p. 4703–4724. Cited by: §1. [11] C. Aitken, R. Masiwal, A. Marchakitus, K. Kowal, M. Gupta, T. Yang, A. Jina, P. Hassanzadeh, W. R. Boos, and M. Kremer (2026) Designing probabilistic ai monsoon forecasts to inform agricultural decision-making. arXiv preprint arXiv:2603.07893. Cited by: §1, §3. [12] O. Watt-Meyer, B. Henn, J. McGibbon, S. K. Clark, A. Kwa, W. A. Perkins, E. Wu, L. Harris, and C. S. Bretherton (2025) ACE2: accurately learning subseasonal to decadal atmospheric variability and forced responses. npj Climate and Atmospheric Science 8 (1), p. 205. Cited by: §1. [13] B. Henn, C. S. Bretherton, N. Kodunov, C. Lessig, M. J. Molina, T. Arcomano, O. Watt-Meyer, G. Couairon, R. Singh, R. Brunstein, et al. (2026) AIMIP Phase 1: systematic evaluations of AI weather and climate models. arXiv preprint arXiv:2605.06944. Cited by: §1, §3. [14] Y. Q. Sun, P. Hassanzadeh, M. Zand, A. Chattopadhyay, J. Weare, and D. S. Abbot (2025) Can ai weather models predict out-of-distribution gray swan tropical cyclones?. Proceedings of the National Academy of Sciences 122 (21), p. e2420914122. Cited by: §1, §3, §3, §3. [15] Y. Q. Sun, P. Hassanzadeh, T. Shaw, H. A. Pahlavan, and A. Marchakitus (2026) Predicting regional gray swans via translocation: AI weather models and Dubai’s unprecedented 2024 rainfall. Science Advances (in press). External Links: Document Cited by: §1, §3, §3. [16] Z. Zhang, E. Fischer, J. Zscheischler, and S. Engelke (2026) Physics-based models outperform ai weather forecasts of record-breaking extremes. Science Advances 12 (18), p. eaec1433. Cited by: §1, §3. [17] A. Chattopadhyay, Y. Q. Sun, and P. Hassanzadeh (2023) Challenges of learning multi-scale dynamics with ai weather models: implications for stability and one solution.. arXiv preprint arXiv:2304.07029. Cited by: §1, §3, §3. [18] M. Bonavita (2024) On some limitations of current machine learning weather prediction models. Geophysical Research Letters 51 (12), p. e2023GL107377. Cited by: §1, §3. [19] T. Selz and G. C. Craig (2023) Can artificial intelligence-based weather prediction models simulate the butterfly effect?. Geophysical Research Letters 50 (20), p. e2023GL105747. Cited by: §A.2, Figure S7, §1, Figure 4, Figure 4, §2.2, §2.4, §2.4, §2.4, §3, Data and code availability:. [20] E. N. Lorenz (1969) The predictability of a flow which possesses many scales of motion. Tellus 21 (3), p. 289–307. Cited by: §1, §2.2, §2.4, §3. [21] T. N. Palmer, A. Döring, and G. Seregin (2014) The real butterfly effect. Nonlinearity 27 (9), p. R123–R141. Cited by: Figure S5, §1, §2.2, §2.2, §2.3, §2.4. [22] T. Palmer (2024) The real butterfly effect and maggoty apples. Physics Today 77 (5), p. 30–35. Cited by: Figure S5, §1, §2.2, §2.2, §2.3, §3. [23] T. Selz and G. Craig (2026) Can ai-based weather prediction models simulate the butterfly effect? the role of architecture and implementation. Journal of Geophysical Research: Machine Learning and Computation 3 (3), p. e2025JH001180. Cited by: §A.1, §A.3, §1, §2.2, §2.4, §2.4, §2.4, §3, §3, §3, §3. [24] E. N. Lorenz (1969) Atmospheric predictability as revealed by naturally occurring analogues. Journal of Atmospheric Sciences 26 (4), p. 636–646. Cited by: §1, §3. [25] H. Van den Dool (1994) Searching for analogues, how long must we wait?. Tellus A 46 (3), p. 314–324. Cited by: §1, §3. [26] E. N. Lorenz (1996) Predictability: a problem partly solved. In Proc. Seminar on predictability, Vol. 1, p. 1–18. Cited by: §A.4, §A.4, §1. [27] D. R. Durran and M. Gingrich (2014) Atmospheric predictability: why butterflies are not of practical importance. Journal of the Atmospheric Sciences 71 (7), p. 2476–2488. Cited by: §1. [28] Y. Q. Sun and F. Zhang (2016) Intrinsic versus practical limits of atmospheric predictability and the significance of the butterfly effect. Journal of the Atmospheric Sciences 73 (3), p. 1419–1438. Cited by: §1, §2.2, §3. [29] F. Zhang, Y. Q. Sun, L. Magnusson, R. Buizza, S. Lin, J. Chen, and K. Emanuel (2019) What is the predictability limit of midlatitude weather?. Journal of the Atmospheric Sciences 76 (4), p. 1077–1091. Cited by: §1, §3. [30] P. T. Vonich and G. J. Hakim (2026) Atmospheric predictability beyond 30 days with machine learning. Artificial Intelligence for the Earth Systems 5 (3), p. 260009. Cited by: §1. [31] J. P. Peixoto, A. H. Oort, M. De Almeida, and A. Tomé (1991) Entropy budget of the atmosphere. Journal of Geophysical Research: Atmospheres 96 (D6), p. 10981–10988. Cited by: §B.1, §2.1, §3. [32] O. Pauluis and I. M. Held (2002) Entropy budget of an atmosphere in radiative–convective equilibrium. part i: latent heat transport and moist processes. Journal of the Atmospheric Sciences 59 (2), p. 140–149. Cited by: §B.1, §2.1, §3. [33] F. Zhang, C. Snyder, and R. Rotunno (2003) Effects of moist convection on mesoscale predictability. Journal of the Atmospheric Sciences 60 (9), p. 1173–1185. Cited by: §A.5, Figure S7, §2.4, §3. [34] T. Selz and G. C. Craig (2015) Upscale error growth in a high-resolution simulation of a summertime weather event over europe. Monthly Weather Review 143 (3), p. 813–827. Cited by: Figure S7, §2.4, §3. [35] S. Van Loon, M. Rugenstein, and E. A. Barnes (2025) Reanalysis-based global radiative response to sea surface temperature patterns: evaluating the ai2 climate emulator. Geophysical Research Letters 52 (14), p. e2025GL115432. Cited by: §3. [36] E. Wu, F. Rebassoo, P. Paul, C. Proistosescu, J. Nugent, D. McCoy, P. Caldwell, and C. S. Bretherton (2025) Applying the ace2 emulator to sst green’s functions for the e3smv3 global atmosphere model. Journal of Geophysical Research: Machine Learning and Computation 2 (3), p. e2025JH000774. Cited by: §3. [37] G. J. Hakim and S. Masanam (2024) Dynamical tests of a deep learning weather prediction model. Artificial Intelligence for the Earth Systems 3 (3). Cited by: §3. [38] G. Craig, T. Selz, M. Beylich, and K. I. Tempest (2026) The physics of ai weather models. arXiv preprint arXiv:2605.23778. Cited by: §3. [39] W. A. Perkins, A. Kwa, J. McGibbon, T. Arcomano, S. K. Clark, O. Watt-Meyer, C. S. Bretherton, and L. M. Harris (2025) HiRO-ace: fast and skillful ai emulation and downscaling trained on a 3 km global storm-resolving model. arXiv preprint arXiv:2512.18224. Cited by: §3, §3. [40] R. Buizza (2010) Horizontal resolution impact on short-and long-range forecast error. Quarterly Journal of the Royal Meteorological Society 136 (649), p. 1020–1035. Cited by: §3. [41] F. Judt (2018) Insights into atmospheric predictability through global convection-permitting model simulations. Journal of the Atmospheric Sciences 75 (5), p. 1477–1497. Cited by: §3. [42] C. Nicolis (1998) Atmospheric analogs and recurrence time statistics: toward a dynamical formulation. Journal of the Atmospheric Sciences 55 (3), p. 465–475. Cited by: §3. [43] P. Platzer, P. Yiou, P. Naveau, P. Tandeo, J. Filipot, P. Ailliot, and Y. Zhen (2021) Using local dynamics to explain analog forecasting of chaotic systems. Journal of the Atmospheric Sciences 78 (7), p. 2117–2133. Cited by: §3. [44] C. J. Stone (1982) Optimal global rates of convergence for nonparametric regression. The annals of statistics, p. 1040–1053. Cited by: §3. [45] L. Györfi, M. Kohler, A. Krzyżak, and H. Walk (2002) A distribution-free theory of nonparametric regression. Springer. Cited by: §3. [46] A. R. Barron (1993) Universal approximation bounds for superpositions of a sigmoidal function. IEEE Transactions on Information Theory 39 (3), p. 930–945. Cited by: §3. [47] T. Poggio, H. Mhaskar, L. Rosasco, B. Miranda, and Q. Liao (2017) Why and when can deep-but not shallow-networks avoid the curse of dimensionality: a review. International Journal of Automation and Computing 14 (5), p. 503–519. Cited by: §3. [48] M. Belkin, D. Hsu, S. Ma, and S. Mandal (2019) Reconciling modern machine-learning practice and the classical bias–variance trade-off. Proceedings of the National Academy of Sciences 116 (32), p. 15849–15854. Cited by: §3. [49] Y. Yu, L. Huang, A. Calotoiu, and T. Hoefler (2026) Scaling laws of global weather models. arXiv preprint arXiv:2602.22962. Cited by: §3. [50] S. Subramanian, A. Kiefer, A. Nigmetov, A. Gholami, D. Morozov, and M. W. Mahoney (2026) On neural scaling laws for weather emulation through continual training. arXiv preprint arXiv:2603.25687. Cited by: §3. [51] A. Lancelin, A. Wikner, L. Dubus, C. Le Priol, D. S. Abbot, F. Bouchet, P. Hassanzadeh, and J. Weare (2026) AI-boosted rare event sampling to characterize extreme weather. Physical Review Letters 137 (6), p. 064201. Cited by: §3. [52] P. T. Vonich and G. J. Hakim (2024) Predictability limit of the 2021 pacific northwest heatwave from deep-learning sensitivity analysis. Geophysical Research Letters 51 (19), p. e2024GL110651. Cited by: §3. [53] N. D. Brenowitz, T. Ge, A. Subramaniam, P. Manshausen, A. Gupta, D. M. Hall, M. Mardani, A. Vahdat, K. Kashinath, and M. S. Pritchard (2025) Climate in a bottle: towards a generative foundation model for the kilometer-scale global atmosphere. arXiv preprint arXiv:2505.06474. Cited by: §A.1, §3. [54] J. Lin, M. Chien, M. Sakarvadia, and E. A. Barnes (2026) Extremes on rewind: generating 1,000-member ensembles initialized at a final condition. arXiv preprint arXiv:2608.19008. Cited by: §3. [55] H. Kim, J. Ryu, S. Son, J. Jeong, H. Kim, and J. Yoon (2026) A spectral test of the butterfly effect and physical consistency in the diffusion-based gencast’s ensembles. npj Climate and Atmospheric Science 9 (1), p. 110. Cited by: §3. [56] C. Lai, P. Hassanzadeh, A. Sheshadri, M. Sonnewald, R. Ferrari, and V. Balaji (2025) Machine learning for climate physics and simulations. Annual Review of Condensed Matter Physics 16 (1), p. 343–365. Cited by: §3. [57] Z. Li, H. Zheng, N. Kovachki, D. Jin, H. Chen, B. Liu, K. Azizzadenesheli, and A. Anandkumar (2024) Physics-informed neural operator for learning partial differential equations. ACM/IMS Journal of Data Science 1 (3), p. 1–27. Cited by: §3. [58] Y. Verma, M. Heinonen, and V. Garg (2024) Climode: climate and weather forecasting with physics-informed neural odes. In International Conference on Learning Representations, Vol. 2024, p. 8408–8430. Cited by: §3. [59] Y. Sha, J. S. Schreck, W. Chapman, and D. J. Gagne (2025) Improving ai weather prediction models using global mass and energy conservation schemes. Journal of Advances in Modeling Earth Systems 17 (11), p. e2025MS005138. Cited by: §3. [60] G. K. Vallis (2006) Atmospheric and oceanic fluid dynamics. Cambridge University Press. Cited by: §3. [61] A. Asch, R. Rossellini, P. Hassanzadeh, and R. Willett (2026) Rigorous uncertainty quantification of probabilistic ai weather forecasts with conformal prediction. arXiv preprint arXiv:2606.19642. Cited by: §3. [62] B. Stevens, M. Satoh, L. Auger, J. Biercamp, C. S. Bretherton, X. Chen, P. Düben, F. Judt, M. Khairoutdinov, D. Klocke, et al. (2019) DYAMOND: the dynamics of the atmospheric general circulation modeled on non-hydrostatic domains. Progress in Earth and Planetary Science 6 (1), p. 61. Cited by: §3. [63] F. Falasca (2025) Probing forced responses and causality in data-driven climate emulators: conceptual limitations and the role of reduced-order models. Physical Review Research 7 (4), p. 043314. Cited by: §3. [64] C. Bodnar, W. P. Bruinsma, A. Lucic, M. Stanley, A. Allen, J. Brandstetter, P. Garvan, M. Riechert, J. A. Weyn, H. Dong, et al. (2025) A foundation model for the earth system. Nature 641 (8065), p. 1180–1187. Cited by: §3. [65] H. Mori (1965) Transport, collective motion, and brownian motion. Progress of theoretical physics 33 (3), p. 423–455. Cited by: §B.1, §B.1, §3. [66] R. Zwanzig (1960) Ensemble method in the theory of irreversibility. The Journal of Chemical Physics 33 (5), p. 1338–1341. Cited by: §B.1, §B.1, §3. [67] J. Wouters and V. Lucarini (2013) Multi-level dynamical systems: connecting the ruelle response theory and the mori-zwanzig approach. Journal of Statistical Physics 151 (5), p. 850–860. Cited by: §B.1, §3. [68] D. Kondrashov, M. D. Chekroun, and M. Ghil (2015) Data-driven non-markovian closure models. Physica D: Nonlinear Phenomena 297, p. 33–55. Cited by: §B.1, §3. [69] F. Falasca and L. Zanna (2026) Physics and causally constrained discrete-time neural models of turbulent dynamical systems. arXiv preprint arXiv:2602.13847. Cited by: §B.1, §3. [70] K. Fraedrich, H. Jansen, E. Kirk, U. Luksch, and F. Lunkeit (2005) The planet simulator: towards a user friendly model. Meteorologische Zeitschrift 14 (3), p. 299–304. External Links: Document Cited by: §A.3. [71] B. Bonev, T. Kurth, C. Hundt, J. Pathak, M. Baust, K. Kashinath, and A. Anandkumar (2023) Spherical fourier neural operators: learning stable dynamics on the sphere. In International conference on machine learning, p. 2806–2823. Cited by: §A.3, Figure S1, §B.2. [72] W. Peebles and S. Xie (2023) Scalable diffusion models with transformers. In Proceedings of the IEEE/CVF international conference on computer vision, p. 4195–4205. Cited by: §A.3. [73] F. Ragone, J. Wouters, and F. Bouchet (2018) Computation of extreme heat waves in climate models using a large deviation algorithm. Proceedings of the National Academy of Sciences 115 (1), p. 24–29. Cited by: §A.3. [74] D. S. Wilks (2005) Effects of stochastic parametrizations in the Lorenz ’96 system. Quarterly Journal of the Royal Meteorological Society 131 (606), p. 389–407. External Links: Document Cited by: §A.4. [75] T. Thornes, P. Düben, and T. Palmer (2017) On the use of scale-dependent precision in earth system modelling. Quarterly Journal of the Royal Meteorological Society 143 (703), p. 897–908. External Links: Document Cited by: §A.4, §A.4. [76] T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama (2019) Optuna: a next-generation hyperparameter optimization framework. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, p. 2623–2631. External Links: Document Cited by: §A.4. [77] L. Li, K. Jamieson, G. DeSalvo, A. Rostamizadeh, and A. Talwalkar (2018) Hyperband: a novel bandit-based approach to hyperparameter optimization. Journal of Machine Learning Research 18 (185), p. 1–52. Cited by: §A.4. [78] G. Benettin, L. Galgani, A. Giorgilli, and J. Strelcyn (1980) Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part 1: Theory. Meccanica 15, p. 9–20. External Links: Document Cited by: §A.4. [79] A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano (1985) Determining Lyapunov exponents from a time series. Physica D: Nonlinear Phenomena 16 (3), p. 285–317. External Links: Document Cited by: §A.4. [80] T. Krishnamurti, K. Rajendran, T. Vijaya Kumar, S. Lord, Z. Toth, X. Zou, S. Cocke, J. E. Ahlquist, and I. M. Navon (2003) Improved skill for the anomaly correlation of geopotential heights at 500 hpa. Monthly weather review 131 (6), p. 1082–1102. Cited by: §A.5. [81] J. Hadamard (2014) Lectures on cauchy’s problem in linear partial differential equations. Courier Corporation. Cited by: §B.1. [82] L. E. Payne (1975) Improperly posed problems in partial differential equations. SIAM. Cited by: §B.1. [83] E. N. Lorenz (1963) Deterministic nonperiodic flow. Journal of the Atmospheric Sciences 20 (2), p. 130–141. Cited by: §B.1. [84] D. Viswanath (1998) Lyapunov exponents from random fibonacci sequences to the lorenz equations. Cornell University. Cited by: §B.1. [85] V. I. Oseledets (1968) A multiplicative ergodic theorem. Lyapunov characteristic numbers for dynamical systems. Transactions of the Moscow Mathematical Society 19, p. 197–231. Cited by: §B.1. [86] L. Young (2013) Mathematical theory of lyapunov exponents. Journal of Physics A: Mathematical and Theoretical 46 (25), p. 254001. Cited by: §B.1. [87] F. Ginelli, H. Chaté, R. Livi, and A. Politi (2013) Covariant lyapunov vectors. Journal of Physics A: Mathematical and Theoretical 46 (25), p. 254005. Cited by: §B.1. [88] W. G. Hoover, C. G. Tull, and H. A. Posch (1988) Negative Lyapunov exponents for dissipative systems. Physics Letters A 131 (3), p. 211–215. External Links: Document Cited by: §B.1. [89] M. Ghil and V. Lucarini (2020) The physics of climate variability and climate change. Reviews of Modern Physics 92 (3), p. 035002. Cited by: §B.1. [90] J. Eckmann and D. Ruelle (1985) Ergodic theory of chaos and strange attractors. Reviews of Modern Physics 57 (3), p. 617. Cited by: §B.1. [91] J. L. Kaplan and J. A. Yorke (1979) Chaotic behavior of multidimensional difference equations. In Functional Differential Equations and Approximation of Fixed Points, Lecture Notes in Mathematics, Vol. 730, p. 204–227. Cited by: §B.1. [92] A. Karimi and M. R. Paul (2010) Extensive chaos in the Lorenz-96 model. Chaos 20, p. 043105. External Links: Document Cited by: §B.1. [93] M. Carlu, F. Ginelli, V. Lucarini, and A. Politi (2019) Lyapunov analysis of multiscale dynamics: the slow bundle of the two-scale Lorenz 96 model. Nonlinear Processes in Geophysics 26, p. 73–89. Cited by: §B.1. [94] L. De Cruz, S. Schubert, J. Demaeyer, V. Lucarini, and S. Vannitsem (2018) Exploring the Lyapunov instability properties of high-dimensional atmospheric and climate models. Nonlinear Processes in Geophysics 25, p. 387–412. Cited by: §B.1. [95] S. Vannitsem (2017) Predictability of large-scale atmospheric motions: Lyapunov exponents and error dynamics. Chaos 27, p. 032101. External Links: Document Cited by: §B.1. [96] R. Lattès and J. L. Lions (1969) The method of quasi-reversibility: applications to partial differential equations. American Elsevier. Cited by: §B.1. [97] A. N. Tikhonov and V. Y. Arsenin (1977) Solutions of Ill-Posed Problems. V. H. Winston & Sons. Note: Translated from the Russian; translation editor F. John. Distributed by Halsted Press (Wiley), New York Cited by: §B.1. [98] K. A. Emanuel, J. David Neelin, and C. S. Bretherton (1994) On large-scale circulations in convecting atmospheres. Quarterly Journal of the Royal Meteorological Society 120 (519), p. 1111–1143. Cited by: §B.1. [99] A. Seif, M. Hafezi, and C. Jarzynski (2021) Machine learning the thermodynamic arrow of time. Nature Physics 17 (1), p. 105–113. Cited by: footnote 1. [100] R. Zwanzig (2001) Nonequilibrium statistical mechanics. Oxford university press. Cited by: §B.1, §B.1. [101] J. L. Lebowitz (1993) Boltzmann’s entropy and time’s arrow. Physics Today 46 (9), p. 32–38. Cited by: §B.1. [102] H. D. Zeh (2007) The physical basis of the direction of time. Vol. 4, Springer. Cited by: §B.1. [103] C. Gardiner (2009) Stochastic methods. Vol. 4, Springer Berlin Heidelberg. Cited by: §B.1. [104] B. D. Anderson (1982) Reverse-time diffusion equation models. Stochastic Processes and their Applications 12 (3), p. 313–326. Cited by: §B.1. [105] U. G. Haussmann and E. Pardoux (1986) Time reversal of diffusions. The Annals of Probability, p. 1188–1205. Cited by: §B.1. [106] A. J. Chorin, O. H. Hald, and R. Kupferman (2000) Optimal prediction and the mori–zwanzig representation of irreversible processes. Proceedings of the National Academy of Sciences 97 (7), p. 2968–2973. Cited by: §B.1. [107] A. J. Majda, I. Timofeyev, and E. Vanden Eijnden (2001) A mathematical framework for stochastic climate models. Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences 54 (8), p. 891–974. Cited by: §B.1. [108] C. Franzke, T. J. O’Kane, J. Berner, P. D. Williams, and V. Lucarini (2015) Stochastic climate theory and modelling. Wiley Interdisciplinary Reviews Climate Change 6, p. 63–78. Cited by: §B.1. Appendix A Methods and Data Throughout the Methods and Data, δtδ t denotes the numerical integration time step used to generate training and testing data (in PlaSim and Lorenz 96), Δt t denotes the input–output prediction interval used to train and infer with an AI model, and τ denotes the forecast or backcast lead time. With the exception of Fig. S10C, all results and analyses presented in this paper are from AI models trained and inferred with single precision (FP32) on the NCAR CISL Derecho cluster using A100 GPUs (more details below). We provide details of the following in this section: • The ERA5 dataset and the AI models trained on it, • The ICON NWP model and data, • The PlaSim GCM, the dataset generated from it, and the AI models trained on this data, • The two-scale Lorenz 96 system, the dataset generated from it, and the AI models trained on this data, • The evaluation metrics. The Supplementary Materials includes • Discussion on time reversal in multi-scale, dissipative dynamical systems, • Movie S1 visualizing the 3D phase space of numerical forecasts and backcasts in Lorenz 63, • More details on the Transformer-, SFNO-, and Diffusion-based AI model architectures. A.1 Observation-based data: ERA5 Data ERA5 [8] is the fifth-generation global atmospheric reanalysis of the European Centre for Medium-Range Weather Forecasts (ECMWF) and serves as the observationally constrained atmospheric dataset. We use five upper-air variables on 17 pressure levels and a selection of surface variables as our prognostic state x (see Table S3 for details). Each AIWP model evolves the prognostic variables in time conditioned on prescribed time-varying forcings (t)b(t) and static boundary fields c, which are supplied as inputs but are not evolved. As described below, we have also trained independent, separate forecasting and backcasting AIWP models with a transformer architecture that is based on Pangu-Weather. These deterministic models (Transformers, hereafter), used in Fig. 1, are trained on data remapped conservatively to a uniform 1∘×1∘1 × 1 latitude–longitude grid. We have trained independent, separate models with Δt=6 t=6 and 2424 h. We train on 6-hourly data from 1979–2018 (40 years) and validate and select checkpoints on 2019. The official Pangu-Weather model, which we obtained from [2], is trained on data with 0.25∘×0.25∘0.25 × 0.25 latitude–longitude grid; there are independent, separate versions trained with Δt=1,3,6 t=1,3,6 and 2424 h. We have only performed inference for “forecasting” with the Pangu-Weather model (Fig. 4). These Transformer AIWP models and the Pangu-Weather model are tested using the same initial conditions taken once a month in 2020 and 2021. Forecasting and backcasting AIWP models: Transformer The forecast and backcast AI models are formulated based on Eqs. (1)–(2), though here, the N takes two additional inputs (t)b(t) and c: (t+Δt) (t+ t) = = θ((t),(t),), _θ (x(t),b(t),c ), (S1) (t−Δt) (t- t) = = ϕ((t),(t),), _φ (x(t),b(t),c ), (S2) where (t)b(t) is the time-varying forcing top-of-atmosphere incident solar radiation and c includes the land–sea mask and surface geopotential field (topography) through channel-wise concatenation. Table S3 lists all of the variables used in x, b, and c. We have trained independent, separate models for Δt=6 t=6 and 2424 h, all on the 1∘×1∘1 × 1 horizontal resolution for computational feasibility. The deterministic forecasting and backcasting neural networks (θN_θ and ϕN_φ) have identical architectures, based on the Pangu-Weather architecture [2], hereafter referred to as the “Transformer”. The full architecture description is provided in the Supplementary Materials. The sole difference between θN_θ and ϕN_φ is the reversed input–output pairing, [(t)→(t−Δt)][x(t) (t- t)] instead of [(t)→(t+Δt)][x(t) (t+ t)], drawn from the same training dataset. Note that backcasting here should not be confused with the common meteorological terms “hindcast” or “reforecast”, which refer to a forecast initialized retrospectively for a historical date (see Fig. 1). Note that this autoregressive backcasting is also different from the backward prediction in cBottle [53], a novel non-autoregressive sampling emulator, which can produce a backward trajectory of a pre-defined length at once. For training, each field is standardized by its mean and standard deviation computed from the training set. Both models are trained on single-time-step Δt t predictions; multi-step trajectories (rollouts) are generated autoregressively at inference by feeding each predicted state back in as input. The Transformer minimizes a weighted mean absolute error (MAE) loss. Below, x represents an ERA5 training sample and x x denotes the AI model prediction; the subscripts a, z, and s index upper-air variables, pressure levels, and surface variables, respectively (e.g., xsx_s includes surface variables and x^a,z x_a,z denotes predicted upper-air variables). All losses are latitude weighted, with ⟨f⟩w≡∑φ,λw(φ)f/∑φ,λw(φ) f _w≡ _ ,λw( )\,f/ _ ,λw( ) denoting the latitude-weighted spatial mean and w(φ)=cosφw( )= (here, λ and φ are longitude and latitude, respectively). The MAE loss is: ℒTransformer=1NaNz∑a,z⟨|x^a,z−xa,z|⟩w+αNs∑s⟨|x^s−xs|⟩w,L_Transformer= 1N_aN_z _a,z | x_a,z-x_a,z| _w+ αN_s _s | x_s-x_s| _w, (S3) where α=1/4α=1/4 gives differential weighting between upper-air and surface variables. Nz=17N_z=17, Na=5N_a=5, and Ns=9N_s=9, are the number of vertical levels, upper-air variables, and surface variables, respectively. The ERA5 forecasting and backcasting Transformers are trained with the AdamW optimizer and a OneCycleLR learning-rate schedule. We have explored a range of hyperparameters, e.g., the peak learning rate, learning-rate schedule, and number of epochs, and selected the best-performing configuration based on the one-time-step (1Δt1 t) validation MAE. At the end, we have found the same hyperparameters for the best-performing forecasting and backcasting models, which have, overall, similar 1Δt1 t accuracy, comparable to that of the official Pangu-Weather forecast, at the same Δt t (Fig. S9). The complete architectural and training hyperparameters are listed in Table S4. Ensemble generation Ensembles are generated by adding perturbations to the initial conditions. Perlin noise following [2] is added to all standardized fields of the initial condition (T, u, v, q, Z on all pressure levels, and the surface variables). Each member is drawn from an independent random seed, and the perturbation amplitude is a function of a reference amplitude ϵ0=0.1 _0=0.1. The initial conditions are taken once a month in 2020 and 2021, and a 16-member ensemble is generated for each. Each two-dimensional perturbation field is a 3-octave sum on the 180×360180× 360 grid with periods of (4, 8, 16) cycles across the domain (spatial scales of roughly 90∘90 , 45∘45 , and 22.5∘22.5 per cycle) and octave weights of (0.2, 0.1, 0.05). Hardware (GPU) precision and details can matter in small-amplitude perturbation experiments [23]: apparent butterflies can be hardware artifacts, as the rapid DKE growth of several AIWP models has been shown to be GPU numerical noise that vanishes or shrinks on CPUs. Therefore, we have carefully ensured that all training of and inference with our AI models (for ERA5, PlaSim, and Lorenz 96) use single floating-point precision (FP32). We have disabled TF32 and conducted all experiments on the same cluster (CISL Derecho) and same GPUs (A100). We have also confirmed that setting ϵ0=0 _0=0 leads to DKE =0=0 over time in Pangu-Weather used in Fig. 4. To further see the effect of hardware, Fig. S10 presents experiments with Pangu-Weather and lower precision (TF32), which show noticeably larger DKE growth. This is due to state-dependent random perturbations that are injected into the calculations over time due to lower precision, consistent with the results reported in [23]. We have further verified that the DKE growth in Fig. 4A is robust to the averaging domain (global versus 60∘60 S–60∘60 N) and to the perturbation type (Perlin versus Gaussian); see Fig. S10. Fig. 4A excludes latitudes poleward of 60∘60 , where the freely running Pangu-Weather (with small Δt t) can develop unphysical behaviors. Part of the rapid DKE growth of the Δt=1 t=1 h model at high latitudes (included in total wavenumbers in Figs. 4B and S6) reflects such instabilities rather than learned dynamics. Cleanly disentangling the two requires further study. Forecasting-backcasting AI model: Preliminary exploration We have conducted preliminary exploration of building a single AI model that performs both forecasting and backcasting, as such a model might better learn dynamics. This bidirectional AI model can be represented as (t+aΔt) (t+a t) = = ϑ((t),a,(t),), _ (x(t),a,b(t),c ), (S4) where a=+1a=+1 (forecast) or −1-1 (backcast) determines the direction of the time evolution. This single neural network is trained with all training samples, [(t)→(t+Δt),a=+1][x(t) (t+ t),a=+1] and [(t)→(t−Δt),a=−1][x(t) (t- t),a=-1], combined. The time indicator a can be fed into the model in different ways. For the Transformer, we use a learned directional embedding, a vector of the token dimension added to the patch embedding of every token, as is done for the positional embedding. Apart from this embedding, the architecture, the weighted MAE loss of Eq. (S3), and all training settings are the same as before (Table S4). We have only trained this model for Δt=24 t=24 h. This preliminary investigation has not shown any noticeable improvement in accuracy or any change in DKE growth, though the design space for encoding a remains largely unexplored. As described later, we reach the same conclusion with a preliminary exploration in PlaSim, for which we also tried a cyclic loss. A.2 Numerical weather prediction (NWP) model: ICON As a physics-based NWP model reference, we use ICON, an operational global model from the German Weather Service and Max Planck Institute for Meteorology, which solves the discretized governing differential equations on an icosahedral grid. The ICON data (the DKE in Figs. 1C and 4A) used in this study are provided by [19], where the ICON ensembles are run at both 2.5 km and 20 km horizontal resolution. A.3 Atmosphere-only GCM: PlaSim GCM setup and data PlaSim [70] is an intermediate-complexity global climate model (GCM) that couples a spectral atmospheric dynamical core, which solves the primitive equations for vorticity, divergence, temperature, and humidity, to a simplified land model. Sea surface temperature (SST) and sea ice are prescribed as boundary conditions that repeat identically each year. Numerical integrations are performed at T42 spectral truncation mapped onto a 64×12864× 128 Gaussian grid with 10 native vertical sigma levels. We have performed 110 years of integration. The first 60 years are discarded as spin-up (admittedly, unnecessarily too long). After integration, the 3D atmospheric state is interpolated onto 13 fixed pressure levels (Table S3). Six-hourly data from simulation years 61–100 (40 years) are used for training and after a 4-year gap, year 105 is used for validation and checkpoint selection. Initial conditions taken every 50 days in Years 106–109 are employed for testing. We train separate, independent AIWP models with Δt=6 t=6 and 2424 h for Transformer and Δt=24 t=24 h for SFNO and Diffusion. Forecasting and backcasting AIWP models: Transformer, Neural Operator, and Diffusion For PlaSim, we train forecast and backcast AIWP models based on three architectures: the same Transformer used for ERA5, a spherical Fourier neural operator (SFNO) [71], and a diffusion transformer (Diffusion, hereafter) [72]. Full architecture descriptions are provided in the Supplementary Materials. The same Eqs. (S1)-(S2), single-time-step training, autoregressive inference, and standardization described for ERA5 are used. Beyond the prognostic state x, the Transformer and SFNO receive the time-varying forcings b and static fields c (Table S3) through channel-wise concatenation, whereas the Diffusion uses them as cross-attention context and is additionally conditioned on day of year and hour of day. For the forecasting Diffusion model, Eq. (S1) represents a sample from the learned conditional distribution pθ((t+Δt)∣(t),(t),)p_θ\! (x(t+ t) (t),b(t),c ) rather than a deterministic mapping to a unique future state. The same applies to the backcasting Diffusion model. All three architectures minimize a latitude-weighted loss using the notation defined earlier for ERA5. The Transformer minimizes the weighted MAE loss of Eq. (S3). The SFNO minimizes a latitude-weighted mean squared error (MSE) on the predicted state: ℒSFNO=1NaNz+Ns[∑a,z⟨(x^a,z−xa,z)2⟩w+∑s⟨(x^s−xs)2⟩w],L_SFNO= 1N_aN_z+N_s [ _a,z ( x_a,z-x_a,z)^2 _w+ _s ( x_s-x_s)^2 _w ], (S5) where atmospheric and surface fields share a single normalizer NaNz+NsN_aN_z+N_s. Nz=13N_z=13, Na=5N_a=5, and Ns=4N_s=4, are the number of vertical levels, upper-air variables, and surface variables, respectively. The Diffusion AIWP model is trained by denoising score matching and minimizes the same latitude-weighted MSE, applied to the predicted noise ε : ℒDiffusion=1NaNz+Ns[∑a,z⟨(ε^a,z−εa,z)2⟩w+∑s⟨(ε^s−εs)2⟩w],L_Diffusion= 1N_aN_z+N_s [ _a,z ( _a,z- _a,z)^2 _w+ _s ( _s- _s)^2 _w ], (S6) where ε is the noise injected at diffusion step τd∈1,…,Ndiff _d∈\1,…,N_diff\ (τd _d is the diffusion step index, distinct from the forecast lead time τ), with Ndiff=100N_diff=100 and a cosine noise schedule. All PlaSim models are trained with the AdamW optimizer; the learning-rate schedules and other training hyperparameters are listed in Table S4. As for ERA5, we have explored a range of hyperparameters (e.g., peak learning rate, learning-rate schedule, and number of epochs). For Transformers, the forecasting and backcasting models are tuned independently based on their single-step (1Δt1 t) validation MAE, which yields the same hyperparameters except for the peak learning rate (the 1Δt1 t errors of the forecast and backcast models are overall comparable, especially for Δt=24 t=24 h; see Fig. S9). We also conducted hyperparameter optimization of the Transformer models against the error of the day-3 prediction, but this does not change any key results or conclusions, so all results reported in the paper are from models whose hyperparameters are chosen based on the 1Δt1 t loss. For the SFNO and Diffusion models, hyperparameters are tuned on the forecasting model only, and the backcasting model adopts the same configuration. Ensemble generation For all AIWP models of PlaSim, ensembles are generated by adding Perlin noise perturbations to initial conditions, following what we did for the ERA5 AIWP models. Because the Diffusion model is stochastic and can otherwise generate ensemble spread intrinsically through its sampling process, we run it in a “quasi-deterministic” mode following [23]: the sampling seed is kept constant across ensemble members, so that ensemble spread arises solely from the explicit initial perturbations (with pre-defined amplitude) rather than from random sampling. Initial conditions are taken every 50 days during simulation years 106–109, and a 16-member ensemble is generated for each. Each two-dimensional noise field is constructed as for ERA5. For the numerical integration of the physics-based PlaSim model, we use the PlaSim’s perturbation scheme implemented in [73]. We generate a 16-member ensemble for each initial condition by adding random white noise of small amplitude (10−810^-8) to the spectral coefficients of the logarithmic surface pressure. Each perturbed state is advanced with the δt=20δ t=20 min time step. The first 6 h is discarded as spin-up to allow rapid dynamical adjustment of the unphysical noise; at 6 h, we then subtract the ensemble mean to obtain the perturbations, which are scaled so that the initial ensemble spread (DKE) matches that of the corresponding AIWP models at 300 hPa (as shown in Fig. 1C). Forecasting-backcasting AI model: Preliminary exploration As done for ERA5, we have explored a single bidirectional forecasting-backcasting Transformer for PlaSim, following Eq. (S4). Here, we have additionally experimented with a cycle-consistency constraint in the loss function, which requires that one Δt t forward (a=+1a=+1) followed by one Δt t backward (a=−1a=-1) returns the original state. However, in these preliminary explorations, neither approach resulted in any improvement in accuracy or any change in the DKE growth. A.4 Canonical multi-scale, chaotic system: Two-scale Lorenz 96 Equations and Data We use the two-scale Lorenz 96 model [26], a canonical testbed for multi-scale, chaotic dynamics. Slow- and large-scale variables XiX_i and fast- and small-scale variables Yi,jY_i,j are governed by dXidt dX_idt =Xi−1(Xi+1−Xi−2)−Xi+F−hγβ∑j=1JYi,j, =X_i-1(X_i+1-X_i-2)-X_i+F- hγβ _j=1^JY_i,j, (S7) dYi,jdt dY_i,jdt =−γβYi,j+1(Yi,j+2−Yi,j−1)−γYi,j+hγβXi, =-γβ\,Y_i,j+1(Y_i,j+2-Y_i,j-1)-γ Y_i,j+ hγβX_i, (S8) with periodic boundary conditions (i=1,…,Ki=1,…,K; j=1,…,Jj=1,…,J). We adopt standard parameters [74] K=8K=8, J=32J=32, F=20F=20, h=1h=1, and β=γ=10β=γ=10, giving a 264-dimensional state vector and chaotic dynamics. Time in this system is commonly measured in model time units (MTU). Based on the error doubling time of the large scales, past studies [26, 75] have estimated that with the current parameters, 1MTU≈51~MTU≈ 5 days, which we have adopted in all figures showing results from the Lorenz 96 system. The forward dynamics is chaotic and dissipative. In this system, dissipation is stronger at the smaller scales by a factor of γ=10γ=10. To generate ground truth and training data, Eqs. (S7)-(S8) are integrated with a fourth-order Runge–Kutta scheme at δt=0.005δ t=0.005 MTU in double floating-point precision (FP64). The initial condition is drawn with independent standard normal values for all 264 components; after a 5,000-step (25 MTU) spin-up to discard transients, the subsequent 10710^7 steps are saved at every δtδ t in double precision. The 10710^7-step dataset is sequentially partitioned into 90% training and 10% validation, and all variables are standardized using the training set’s statistics. The validation segment is used for hyperparameter pruning, the learning-rate schedule, early stopping, and model selection (more details below). After model selection, an additional independent test trajectory equal in length to the validation segment is generated. The 100 evaluation initial conditions are drawn from this test trajectory and are spaced 1,000 δtδ t (5 MTU, more than 100 Lyapunov times) apart, far beyond the system’s decorrelation time. For the numerical results in Fig. 2A, we numerically re-integrate the equations from a reference state on the archived trajectory both forward and backward in time, the latter using the same Runge–Kutta scheme with a negated time step −δt-δ t, in both single (FP32) and double (FP64) precision (the former is used to match the precision of AI models’ training and inference). To prevent floating-point overflow, state values in the numerical integration are capped at ±1010± 10^10; the backward trajectory reaches this cap within a fraction of an MTU. In Fig. 2A, the saturated color regions correspond to |Xi||X_i| exceeding the color-scale range as the backcast solution grows without bound; in Fig. 3, crosses mark the last lead time at which the ensemble DKE of the backcasting numerical integration remains below 10310^3 (the plotted range). Forecasting and backcasting AI models: Multilayer perceptron Following Eqs. (1) and (2), we train independent forecast and backcast AI models for the Lorenz 96 system. Each is a multilayer perceptron with ReLU activations, of width W (neurons per hidden layer) and depth D (number of hidden layers), that outputs the full state x at t±Δt± t directly. The AI-model prediction interval is Δt=nδt t=n\,δ t with n∈1,5,10n∈\1,5,10\ (reminder that δt=0.005δ t=0.005 MTU is the numerical solver’s time step). Note that the largest time step of AI models, 10δt=0.0510δ t=0.05 MTU =0.25=0.25 day, is comparable to one Lyapunov time (1/λ≈0.041/λ≈ 0.04 MTU, from the measured leading finite-time Lyapunov exponent, FTLE, in Table S2). The Lorenz 96 results are based on 12 independently trained AI models: • Forecast and backcast models, • Models with =Xix=X_i and =(Xi,Yi,j)x=(X_i,Y_i,j) (from concatenation), • Models with Δt=nδt t=nδ t where n=1,5,or 10n=1,5,or\,10. For each n, training pairs ((t),(t±nδt)) (x(t),\,x(t± nδ t) ) use all saved samples and each model is trained on the full 9×1069× 10^6-sample training set. All models are trained with the Adam optimizer by minimizing the MSE of the predicted state. We have selected the best performing AI models based on extensive hyperparameter optimization, though the key results and conclusions are robust with respect to these choices. We have explored a range of width W∈512,1024,2048W∈\512,1024,2048\ and depth D∈3,5,7D∈\3,5,7\ and four weight initialization seeds. Then for each fixed (W,D)(W,D) configuration and seed, the learning rate and batch size have been tuned independently for each of the twelve models with the Optuna framework [76] and a Hyperband pruner [77] (Table S5). All results are reported using the best-performing AI models with W=1024W=1024 and D=7D=7, selected based on the criterion of using the same (W,D)(W,D) for forecast and backcast that produce similar lowest single-step validation error. Single time-step (one Δt t) forecast and backcast RMSE are comparable for all models (backcast-to-forecast RMSE ratios of 0.90−1.130.90-1.13 except for one case; see Table S6), indicating that the forecast–backcast asymmetry in Figs. 2 and 3 is dominated by autoregressive rollout rather than by differences in the single time-step accuracy of the AI models. Ensemble generation Ensembles are generated by adding random perturbations to initial conditions: independent Gaussian noise of amplitude ϵ0=0.1 _0=0.1 is added to all components of the initial state x to produce a 100-member ensemble. Reduced-amplitude ensembles are obtained by scaling down ϵ0 _0. Computation of leading Lyapunov exponents in Lorenz 96 In multi-scale chaotic systems like Lorenz 96, the leading global Lyapunov exponent is set by the fast, small-scale variables (Y): a generic infinitesimal perturbation asymptotically grows at this fastest rate, so the infinite-time measure is uninformative about the large scales (X). Predictability measures such as the error doubling time [75] are informative about the large scales, but only because they track finite-amplitude errors, whose fast, small-scale components saturate while the large-scale error continues to grow. We instead require the growth rate of infinitesimal perturbations as they impact the large-scale (X) variables alone. Thus, for the numerical and AI models, we estimate a projected FTLE: continuous perturbation rescaling keeps the perturbation infinitesimal, following the classical algorithm of Benettin et al. [78] as implemented by Wolf et al. [79], and the growth rate is measured, over a finite window, from the X components alone. Because the estimator requires only evaluations of the model itself (no tangent-linear equations), it is applied identically to the numerical and AI models. A random initial state 0x_0 and a perturbed state 0′=0+δ00x _0=x_0+ _0v_0 (where 0v_0 is a random unit vector in the full state space and δ0=10−4 _0=10^-4) are integrated through the AI model (time step Δt t) or the numerical model (time step δtδ t). At each step m, the full perturbation vector Δm=m′−m _m=x _m-x_m is recorded, and the perturbed state is re-normalized to magnitude δ0 _0 along the new divergence vector m=Δm/‖Δm‖v_m= _m/\| _m\|. The projected forward FTLE of the large-scale variables, λ(X)λ^(X), is the time average over a finite window of NstepN_step steps: λ(X)=1NstepΔt∑m=1Nstepln(‖Δm(X)‖δ0‖m−1(X)‖),λ^(X)= 1N_step t _m=1^N_step ( \| ^(X)_m\| _0\|v^(X)_m-1\| ), (S9) where, to isolate the predictability of the large scales, the logarithmic growth rate uses exclusively the norm of the X components. The fast-variable exponents λ(Y)λ^(Y) in Table S2 are computed similarly but using the Y components. The averaging window spans 0.10.1–0.30.3 MTU after a 0.10.1-MTU transient, so Nstep=0.2MTU/ΔtN_step=0.2~MTU/ t (40 steps for the numerical model and the 1δt1δ t AI models; 4 steps at Δt=10δt t=10δ t). To capture the variability of predictability across the reference trajectory, reported values represent the mean and standard deviation evaluated over 20 independent initial conditions drawn from the reference trajectory. The λ(X)λ^(X) and λ(Y)λ^(Y) for the ground truth and some of the key AI models are reported in Table S2. A.5 Evaluation metrics Forecast and backcast skills are assessed using two complementary metrics: the anomaly correlation coefficient (ACC) and the root mean square error (RMSE), which measure the accuracy of trajectory predictions. We also compute the difference kinetic energy (DKE), which quantifies the growth of initial perturbations and characterizes the chaotic behavior of the system. ACC and RMSE are evaluated for several predicted variables, whereas DKE is evaluated from the horizontal wind components. For Lorenz 96, RMSE and DKE are both based only on the large-scale variable, X. Let τ denote forecast or backcast lead time, i=1,…,Ni=1,…,N the grid point index, and wi=cosφiw_i= _i the area weight for latitude φi _i. For deterministic verification, x^i(τ) x_i(τ) and xi(τ)x_i(τ) are the predicted and ground truth values, and x¯ic x_i^\,c is the climatological mean at grid point i. For ensemble verification, ui(m)(τ)u_i^(m)(τ) and vi(m)(τ)v_i^(m)(τ) are the zonal and meridional winds of member m=1,…,Mm=1,…,M, with instantaneous ensemble means u~i(τ) u_i(τ) and v~i(τ) v_i(τ). All metrics are averaged over the evaluation cases (inference from many initial conditions in the test set) for each system to yield lead-time-dependent skill curves. Anomaly correlation coefficient (ACC) ACC measures the area-weighted pattern similarity between predicted and target anomalies relative to climatology (defined as the calendar-day mean over 1979–2018 for ERA5 and years 61–100 for PlaSim), ACC(τ)=∑i=1Nwi(x^i(τ)−x¯ic)(xi(τ)−x¯ic)∑i=1Nwi(x^i(τ)−x¯ic)2∑i=1Nwi(xi(τ)−x¯ic)2.ACC(τ)= _i=1^Nw_i ( x_i(τ)- x_i^\,c ) (x_i(τ)- x_i^\,c ) _i=1^Nw_i\, ( x_i(τ)- x_i^\,c )^2\; _i=1^Nw_i\, (x_i(τ)- x_i^\,c )^2. (S10) Note that ACC primarily measures the large-scale forecast accuracy. The conventional skill threshold for weather prediction is ACC>0.6ACC>0.6 [80]. For all ERA5 and PlaSim results, we define the asymmetry as the ratio of forecast to backcast lead times at ACC=0.6ACC=0.6. Root mean square error (RMSE) RMSE quantifies prediction error as RMSE(τ)=∑i=1Nwi(x^i(τ)−xi(τ))2∑i=1Nwi.RMSE(τ)= _i=1^Nw_i\, ( x_i(τ)-x_i(τ) )^2 _i=1^Nw_i. (S11) For Lorenz 96, RMSE is defined as in Eq. (S11) with uniform weights over the K=8K=8 large, slow variables X. The asymmetry in Lorenz 96 is then defined as the backcast/forecast RMSE ratio at a fixed lead (1.5 days; Fig. 2). Difference kinetic energy (DKE) DKE measures ensemble spread rather than error against a target. It is the area-weighted mean of half the ensemble variance of the horizontal wind components, computed about the instantaneous ensemble mean (u~i(τ) u_i(τ), v~i(τ) v_i(τ)), distinct from the calendar-day climatological mean x¯ic x_i^\,c used in ACC. Analogous to the kinetic component of the difference total energy of [33], it quantifies spread directly in dynamical units: DKE(τ)=∑i=1Nwi12[σu,i2(τ)+σv,i2(τ)]∑i=1Nwi,σu,i2(τ)=1M−1∑m=1M(ui(m)(τ)−u~i(τ))2,DKE(τ)= _i=1^Nw_i\, 12 [ _u,i^2(τ)+ _v,i^2(τ) ] _i=1^Nw_i, _u,i^2(τ)= 1M-1 _m=1^M (u_i^(m)(τ)- u_i(τ) )^2, (S12) where σu,i2 _u,i^2 and σv,i2 _v,i^2 are the unbiased ensemble variances of the two horizontal wind components at grid point i, wiw_i is the area weight, M is the ensemble size, and N the number of grid points. Units are m2s−2m^2\,s^-2. Globally averaged DKE at 300 hPa is reported, unless specified otherwise. For the Lorenz 96 system, DKE is analogously defined as half the sum of the ensemble variances of the large-scale variables: DKELorenz 96(τ)=12∑i=1K1M−1∑m=1M(Xi(m)(τ)−X~i(τ))2,DKE_Lorenz\,96(τ)= 12 _i=1^K 1M-1 _m=1^M (X_i^(m)(τ)- X_i(τ) )^2, (S13) where Xi(m)X_i^(m) is the large-scale state of ensemble member m and X~i X_i is the instantaneous ensemble mean. Appendix B Supplementary Materials B.1 Time reversal in multi-scale, dissipative, chaotic systems and the 2nd2^nd law of thermodynamics This section further describes why predicting the past of a multi-scale, dissipative, chaotic system, such as the atmosphere, should be practically impossible. It then explains why backcasting AI models can nevertheless be skillful without violating any of these arguments, or the second law of thermodynamics. We use three examples of increasing relevance: the heat equation, Lorenz 63, and the two-scale Lorenz 96 system (Methods and Data). Multi-scale dynamics: Heat equation and Hadamard ill-posedness The linear, one-dimensional heat equation, ∂T/∂t=κ∂2T/∂x2∂ T/∂ t=κ\,∂^2T/∂ x^2 with diffusivity κ>0κ>0 and spatial coordinate x, evolves each Fourier mode of wavenumber k independently: T^(k,t)=T^(k,0)e−κk2t. T(k,t)= T(k,0)\,e^-κ k^2t. (S14) Forward in time (t>0t>0), perturbations decay, most quickly at the smallest spatial scales; the analytical and numerical solutions are stable to initial-condition errors. Reversing time (t∗=−t>0t^*=-t>0) flips the exponent: a perturbation of amplitude ε at wavenumber k, unavoidable in any real computation or measurement, grows backward as εe+κk2t∗ \,e^+κ k^2t^* and reaches order unity by te∗(k)=ln(1/ε)κk2t^*_e(k)= (1/ )κ k^2. Even at double floating-point precision (ε∼10−16 10^-16, ln(1/ε)≈37 (1/ )≈ 37), te∗t^*_e quickly collapses quadratically with wavenumber: the smallest scales dominate the backcasting solution first, and in a system with a broad range of scales, they do so almost immediately. This is the classical statement that the backward heat problem is ill-posed in the sense of Hadamard [81, 82]. Note that here, the exponents in forecasting and backcasting have the same amplitudes (as a function of length scale) but opposite signs (−κk2-κ k^2 versus +κk2+κ k^2). Chaotic dynamics: Lorenz 63 Chaotic dissipative systems further refine the picture of backcasting. The nonlinear, single-scale Lorenz 63 system [83], dxdt dxdt =σ(y−x), =σ(y-x), (S15) dydt dydt =x(ρ−z)−y, =x(ρ-z)-y, (S16) dzdt dzdt =xy−βz, =xy-β z, (S17) with standard parameters (σ=10σ=10, ρ=28ρ=28, β=8/3β=8/3) has Lyapunov exponents (λ1,λ2,λ3)≈(0.91, 0,−14.57)( _1, _2, _3)≈(0.91,\,0,\,-14.57) [84]. Their sum equals the phase-space divergence: ∑i=13λi=−(1+σ+β)≈−13.67, _i=1^3 _i=-(1+σ+β)≈-13.67, (S18) i.e., volumes contract, and trajectories settle onto an attractor, the Lorenz strange attractor, while nearby states separate exponentially at the modest rate λ1≈+0.91 _1≈+0.91 (Movie S1). Reversing time negates the Lyapunov spectrum, which reorders to ≈(+14.57, 0,−0.91)≈(+14.57,\,0,\,-0.91); see below for more discussion on this reordering. Two major consequences follow, and Movie S1 shows both. First, the reversed system is dramatically more unstable along trajectories: infinitesimal errors e-fold 14.57/0.91≈1614.57/0.91≈ 16 times faster in backcasting compared with forecasting, so an initial error reaches order one 16×16× faster. Second, and more damning, the sum of the reversed exponents is −∑i=13λi=+(1+σ+β)≈+13.67.- _i=1^3 _i=+(1+σ+β)≈+13.67. (S19) Backward in time, phase-space volumes expand, and the attractor of the forward dynamics is a repeller of the backward dynamics. Any perturbation transverse to the attractor, which forward dynamics would harmlessly contract away, is exponentially amplified backward at rates up to e+14.57t∗e^+14.57t^*, and trajectories are expelled from the attractor altogether (see Movie S1). As a result, numerical backcasting of single-scale Lorenz 63 fails for two distinct reasons: errors grow faster and into regions the dynamics never visits. Note that here, unlike for the heat equation, time-reversal not only negates the signs of exponents but also changes the amplitude of the leading exponent, leading to 16×16× faster growth of errors in backcasting. These two failures follow from two distinct properties of Lorenz 63’s Lyapunov spectrum under time reversal: A) The sign reversal of all exponents, and B) |λ3|>λ1| _3|> _1. Next, we discuss the conditions under which (A) and (B) apply to more complex systems such as the atmosphere. As for (A), if the dynamics is invertible, admits an ergodic invariant probability measure, and satisfies a mild integrability condition on the tangent dynamics, the Oseledets multiplicative ergodic theorem guarantees a well-defined Lyapunov spectrum, and the exponents of the time-reversed dynamics are the forward exponents with reversed signs [85, 86, 87, 88]. For dissipative systems such as the atmosphere [89], the invariant measure is supported on the attractor, so the reversal applies to trajectories on the attractor itself [90]; e.g., the initial conditions that we choose from forward numerical integrations of Lorenz 96 and PlaSim, or from ERA5 reanalysis. Off-attractor states diverge under the reversed flow and have no backward Lyapunov exponents [90]; such states are not available from ERA5 reanalysis or any forward numerical integration (Lorenz 96 or PlaSim) and are not of concern here. Therefore, (A) is broadly true. Let us order the L Lyapunov exponents of a system from the most positive λ1 _1 (leading) to the most negative λL _L. Then, property (B) is that |λL|>λ1| _L|> _1, so that upon reversal and reordering, the leading backward exponent exceeds the leading forward one, making backward error growth exponentially faster than forward. In three-dimensional dissipative, chaotic systems such as Lorenz 63 this is automatic: the exponent sum equals the (negative) phase-space divergence and the intermediate exponent vanishes, so |λ3|=λ1+(σ+1+β)>λ1| _3|= _1+(σ+1+β)> _1 (equivalently, the Kaplan–Yorke dimension lies below three [91]). In higher dimensions, however, a negative exponent sum only guarantees that the contracting directions collectively outweigh the expanding ones, so |λL|>λ1| _L|> _1 is not guaranteed in general. But the picture changes once we consider multi-scale systems. In such systems, dissipation acts most strongly on the smallest and fastest scales, leading to |λL|>λ1| _L|> _1, as documented in computed Lyapunov spectra of Lorenz 96 [92, 93] and of atmospheric and coupled climate models [94, 95]. Therefore, both (A) and (B) are expected to apply to the atmosphere. Considering the above discussion on the heat equation and Lorenz 63, backcasting of multi-scale, dissipative, chaotic systems should be practically impossible. Next, we discuss a canonical system that embodies all three characteristics, the two-scale Lorenz 96. Multi-scale, chaotic dynamics: Two-scale Lorenz 96 This system (Eqs. (S7)-(S8)) adds the ingredient the atmosphere has and Lorenz 63 lacks: multi-scale dynamics, specifically, faster dissipation at smaller scales Yi,jY_i,j. Under time reversal, the damping terms −Xi-X_i and −γYi,j-γ Y_i,j become anti-dissipative, with the fast variables, growing γ=10γ=10 times more strongly, blowing up first and dragging the large scales XiX_i with them through the coupling. Consequently, backcasting via numerical integration diverges rapidly at both single and double floating-point precision (Fig. 2A). The second law of thermodynamics, the arrow of time, and the asymmetry It may seem that a skillful AI backcasting conflicts with the second law of thermodynamics: entropy increase selects a direction (arrow) of time, and the atmosphere, a forced-dissipative system full of mixing and irreversible diabatic processes, plainly should follow that arrow. However, coarse-graining removes the spatio-temporal scales at which the strongest dissipation happens, leading to a training set that is closer to reversible. Seen in a different way, coarse-graining mimics the classical remedy to Hadamard ill-posedness: regularization by truncating/damping the high wavenumbers, after which the time-reversed integrations become stable over a finite window [96, 97]. Furthermore, the training data contain only states on (or extremely near) the attractor, so the AI model never learns, and is never asked to represent, the explosive off-attractor directions (and as described above, because of the coarse-graining, the on-attractor instability the training data inherits is milder as well). A learned backcast does not integrate the reversed vector field; it estimates, by regression, the attractor state Δt t earlier that is most compatible with the current coarse state. It operates on the attractor, at the “resolved” scales, or not at all. In the end, the skillful AI backcasting does not violate the second law of thermodynamics; it just operates on a dataset that is closer to reversible due to coarse-graining and on-attractor sampling. The forecast-backcast asymmetry that we see in the hierarchy (Figs. 1B and 2B-C) is then a measure of the residual irreversibility (dynamical dissipation and, for the atmosphere, thermodynamic entropy production) that survives coarse-graining and is retained in the training data. This interpretation makes a testable prediction: the asymmetry should be largest where irreversible processes are strongest in the resolved dynamics of the atmosphere. The larger asymmetry we find in the tropics, where diabatic heating and convective dissipation dominate [31, 32, 98], compared with the extratropics (Table S1) is consistent with this prediction, and warrants further investigation. Finally, it is worth mentioning the double role “coarse-graining” might appear to play, and the distinction between the two operations this term refers to can be instructive. Microscopic dynamics (at the molecular level) are time-reversal symmetric; irreversibility and entropy growth emerge at the level of coarse-grained, macroscopic descriptions11 1 It is worth mentioning that Seif et al. [99] has addressed the question of whether neural networks can learn the “thermodynamic arrow of time” from movies of microscopic processes., when information about the discarded degrees of freedom is given up, as formalized from Boltzmann’s molecular chaos assumption to the projection-operator formalisms of Mori and Zwanzig [66, 65, 100]. Thus, coarse-graining from molecular (microscopic) to continuous (macroscopic) dynamics is the operation through which the latter acquires its arrow of time and time-irreversibility [101, 102]. In the context of the work here, coarse-graining refers to the removal of the fast and small scales of the macroscopic dynamics, which brings the training data closer to time-reversible. The role of stochasticity and the Mori-Zwanzig Formalism All of the arguments above concern deterministic dynamics. Stochastic dynamics behave qualitatively differently under time reversal. Consider the linear ordinary differential equation dx/dt=−axdx/dt=-a\,x with a>0a>0, which, for example, governs the Fourier mode of the heat equation with a=κk2a=κ k^2. The forecast solution is x(t)∼e−atx(t) e^-at, which exponentially decays. The time-reversed equation is dx/dt∗=+axdx/dt^*=+a\,x, where t∗=−t>0t^*=-t>0. Thus, the backcast solution is x(t∗)∼e+at∗x(t^*) e^+at^*, which exponentially blows up. Now let’s consider the stochastic version driven by Gaussian white noise η(t)η(t): dxdt=−ax+η(t), dxdt=-a\,x+η(t), (S20) i.e., an Ornstein–Uhlenbeck process [103]. When the noise is additive and the process is statistically stationary, the equation obeyed by the time-reversed process [104, 105] is: dxdt∗=−ax+η~(t∗), dxdt^*=-a\,x+ η(t^*), (S21) where η~ η is another realization of the same white noise. This is basically the same equation as (S20), showing that the stationary stochastic process is time-reversible. The past is exactly as predictable as the future, which is dramatically different from the behavior of the same equation when the noise vanishes (η=0η=0). The Mori–Zwanzig formalism [65, 66, 100] reconciles the coarse-graining versus stochasticity explanations. By the Mori–Zwanzig formalism, eliminating the fast, small scales of a deterministic multiscale system (coarse-graining) leaves resolved-scale dynamics that contain memory and a stochastic forcing [106]. This observation underlies stochastic climate modeling [107, 108, 69] and stochastic parameterizations [67, 68]. Coarse-grained deterministic data are thus statistically of the same kind as realizations of a stochastic process: the stochastic-process view of ERA5 and the coarse-graining explanation of the main text are one interpretation expressed in two languages. B.2 More details on neural network architectures The Transformer model The Transformer follows the three-dimensional Earth-specific transformer architecture of [2]. A 3D patch embedding groups every 2×2×22×2×2 block in (level, latitude, longitude) into a 240240-dimensional token. These tokens pass through a hierarchical encoder–decoder of Swin-style shifted-window self-attention blocks in four stages of depth 2, 6, 6, and 2, each attending within a local window spanning two pressure levels, with 6, 12, 12, and 6 attention heads per stage. The encoder halves the horizontal resolution at each stage and the decoder restores it; a sub-pixel deconvolution head and a recovery head reconstruct the full-resolution prognostic state. The ERA5 Transformer uses a window size of 2×6×102×6×10. The PlaSim Transformer is identical except for a window size of 2×6×122×6×12, matching the T42 grid; see Table S4 for details. The spherical Fourier Neural Operator (SFNO) model The SFNO [71] uses spherical harmonic transform (SHT), so that spectral filtering respects the Earth’s spherical geometry. In plain terms, the model learns transformations in spherical spectral space rather than directly on the latitude–longitude grid (used in FourCastNet [1]). A point-wise encoder lifts the input to a 384384-dimensional channel space. Twelve spectral blocks each apply a forward SHT, a learned degree-dependent linear filter, and an inverse SHT, followed by a two-layer GELU multilayer perceptron (expansion ratio 2) with a residual connection and per-channel instance normalization. A global residual connection links encoder to decoder. The Diffusion model This is a conditional score-based diffusion model. This generative model learns to reverse an artificial noising process and generates a target state through iterative denoising. A forward process adds Gaussian noise to the target state over Ndiff=100N_diff=100 steps with a cosine noise schedule. The denoiser is a vision transformer that tokenizes the input via a 2×22×2 patch embedding into 10241024-dimensional tokens, to which a learned spherical-harmonics positional embedding is added. Twelve self-attention blocks (16 heads each) operate under adaptive layer normalization conditioned on the diffusion noise level and calendar position (day-of-year and hour-of-day), setting per-block shift, scale, and gate parameters. After each of the first 6 blocks, a cross-attention layer allows state tokens to attend to boundary-forcing context tokens. A noise-level-modulated unpatchify head projects back to the prognostic channels. 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>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 t. For each variable and AI model, the largest value among the three regions is shown in bold. For ERA5, consistent across both Δt 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. 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 Z500 Global 1.40 1.63 1.63 1.76 1.59 1.94 NH 1.35 1.64 1.59 1.73 1.48 1.89 Tropics 1.53 1.97 2.03 2.75 1.82 2.17 SH 1.39 1.65 1.63 1.72 1.70 1.94 U850 Global 1.47 2.03 1.82 2.35 1.77 2.19 NH 1.38 2.07 1.80 2.34 1.75 2.22 Tropics 1.79 2.52 1.79 2.25 1.71 2.27 SH 1.33 1.84 1.80 2.24 1.83 2.06 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 _0=10^-4, computed over the window 0.10.1–0.30.3 MTU (after discarding a 0.10.1-MTU transient). The FTLE values are projected onto the slow (X) or fast (Y) variables. System / AI model Projection λ (MTU-1) Ground truth X 21.9 ± 6.3 Ground truth Y 24.3 ± 4.0 XYXY, 1δt1δ t AI model X 26.8 ± 6.5 XYXY, 1δt1δ t AI model Y 26.8 ± 5.9 X, 1δt1δ t AI model X 3.05 ± 1.84 X, 10δt10δ t AI model X 2.81 ± 1.59 Table S3: Variable inventories for the Pangu-Weather, ERA5, and PlaSim AIWP models. Upper-air prognostic state uax_ua, surface prognostic state sfcx_sfc, 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∘0.25 grid, which we use for inference only. The ERA5 Transformer column describes the data for the 1∘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. 1B-C remain the same. 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 uax_ua 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 ✓ ✓ ✓ 10-m zonal & meridional wind ✓ ✓ — Mean sea-level pressure ✓ ✓ — Surface sfcx_sfc Surface pressure — ✓ ✓ Skin temperature — ✓ ✓ 0–7 cm soil water & soil temp. — ✓ — Soil moisture — — ✓ Sea-surface temperature — ✓ — Time-varying top-of-atmosphere incident solar radiation — ✓ ✓ forcings b Sea-surface temperature — — ✓† Sea-ice concentration — — ✓ Static fields c Land–sea mask — ✓ ✓ Surface geopotential — ✓ ✓ Horizontal grid 721×1440721× 1440 (0.25∘0.25 ) 180×360180× 360 (1∘1 ) 64×12864× 128 (T42) Pressure levels (hPa) 50, 100, 150, 200, 250, 300, 400, 500, 600, 700, 850, 925, 1000 5, 10, 20, 30, 50, 70, 100, 150, 250, 300, 400, 500, 600, 700, 850, 925, 1000 50, 100, 150, 200, 250, 300, 400, 500, 600, 700, 850, 925, 1000 Time interval Δt t 1 h, 3 h, 6 h, 24 h 6 h, 24 h 6 h, 24 h Training period 1979–2017 (39 years) 1979–2018 (40 years) 61–100 (40 years) Validation period 2019 (1 year) 2019 (1 year) 105 (1 year) Test period 2020-2021 (2 years) 2020-2021 (2 years) 106–109 (4 years) 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. ERA5 Transformer PlaSim Transformer PlaSim SFNO PlaSim Diffusion Embedding dim DembD_emb 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 ×102\!×\!6\!×\!10 ×122\!×\!6\!×\!12 — — Patch size ×22\!×\!2\!×\!2 ×22\!×\!2\!×\!2 — ×22\!×\!2 Optimizer AdamW AdamW AdamW AdamW Training Peak learning rate ×10−45\!×\!10^-4 ×10−46\!×\!10^-4 / ×10−42\!×\!10^-4† ×10−46\!×\!10^-4 ×10−55\!×\!10^-5 LR schedule OneCycleLR OneCycleLR ReduceLROnPlateau warmup + cosine decay Weight decay 0.01 0.01 0.01 0.05 Batch size 32 32 32 32 Epochs 120 120 200 800 Loss function weighted MAE weighted MAE spherical MSE weighted MSE Diffusion Steps NdiffN_diff — — — 100 Noise schedule — — — cosine Table S5: Architectural and training hyperparameters for the Lorenz 96 AI models. Set notation indicates discrete values explored during the hyperparameter optimization phase. Phase Hyperparameter Value / Search space Architecture multilayer perceptron, ReLU activations Architecture Width W 512,1024,2048\512,1024,2048\ Depth D 3,5,7\3,5,7\ hidden layers Input dimension 8 (X only) or 264 (XYXY) Optimizer Adam Hyperparameter search Learning rate Log-uniform [10−4,10−2][10^-4,10^-2] Batch size 256,512,1024,2048,4096\256,512,1024,2048,4096\ Pruning Hyperband (150 epochs/trial) Max epochs 500 Final training LR schedule ReduceLROnPlateau (factor 0.5, patience 15) Early stopping patience 50 epochs Training loss mean squared error (MSE) Evaluation metric mean absolute error (MAE) Table S6: One-time-step RMSE of the forecast and backcast AI models. RMSE of a single prediction step of length Δt t, computed over the X variables against the reference trajectory, for the selected AI models with W=1024W=1024 and D=7D=7. Values are the mean ± standard deviation over the same 100 initial conditions as Fig. 2. x Δt t Forecast RMSE Backcast RMSE X 1δt1δ t 0.0056 ± 0.0020 0.0056 ± 0.0021 X 5δt5δ t 0.026 ± 0.0093 0.027 ± 0.011 X 10δt10δ t 0.048 ± 0.016 0.046 ± 0.018 XYXY 1δt1δ t 0.038 ± 0.012 0.043 ± 0.015 XYXY 5δt5δ t 0.16 ± 0.054 0.23 ± 0.067 XYXY 10δt10δ t 0.29 ± 0.073 0.26 ± 0.066 Figure S1: Missing butterfly and forecast-backcast asymmetry in the SFNO-based AIWP models of PlaSim GCM. As in Fig. 1B–C, but for independently trained forecasting and backcasting AIWP models, based on spherical Fourier neural operator (SFNO) architecture [71], with Δt=24 t=24 h (see Methods and Data). ACC of U850U850, U250U250, and Z500Z500 versus lead time are shown, as well as the DKE at 300 hPa for ensembles initialized with perturbations of amplitude ϵ0=0.1 _0=0.1, ϵ0/10 _0/10, and ϵ0/103 _0/10^3. 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 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 t in Transformer-based AIWP models of ERA5 and PlaSim GCM. As in Fig. 1B, but for U850 (top), Z500 (middle), and U250 (bottom) for ERA5 and PlaSim Transformers (Δt=24 t=24 and 6 h). Colors denote Δt 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 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 t=10δ t. Perturbations are Gaussian, with amplitudes decreasing from a reference value, ϵ0 _0. Curves are averaged over 10,00010,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 t. Power spectra E(n)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 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 t of the official Pangu-Weather model. Panels marked “No data” correspond to lead times shorter than the model’s Δt 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 _0/10^3) for one representative initial condition. The Δt=1 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 11-h model instead resembles a mixture of small-scale features amplified in the 66-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)E(n), as a function of total wavenumber n, computed with the spherical harmonic transform, for one-time-step (1Δt1 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-5/3 and −3-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Δt1 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 t=1, 3, 6, and 24 h [2]) and our forecasting (blue) and backcasting (red) Transformers (Δt=6 t=6 and 24 h). Right: the PlaSim forecasting and backcasting Transformers (Δt=6 t=6 and 2424 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 t=24, 6, 3, and 1 h. (A) The same as Fig. 4A, the reference, which uses DKE averaged over 60o60^oS-60o60^oN, 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 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. Caption for Movie S1. Numerical forecasts and backcasts on the Lorenz 63 attractor. Animation of the Lorenz 63 (Eqs. (S15)-(S17)) states numerically integrated forward and backward from 10 initial conditions that are slightly different (same numerical methods as those used for Lorenz 96; see Methods and Data). The Lorenz attractor (thin gray lines) is visualized in the x–y–z phase space. The initial conditions are on the attractor. The forecast trajectories remain on the attractor (due to the phase-space contraction), though they diverge significantly in the long term due to the +0.91+0.91 leading Lyapunov exponent. The backcast trajectories quickly leave the attractor, which is now a repeller (due to the phase-space expansion). See the Supplementary Materials for more discussion. The visualization code is generated by Claude Fable 5. The movie is available at https://doi.org/10.5281/zenodo.22062019.