Paper deep dive

Do AI Forecast Ensembles Sample the Correct Conditional Distribution?

Lucas J. Howard, Elizabeth A. Barnes

Year: 2026Venue: arXiv preprintArea: physics.ao-phType: PreprintEmbeddings: 42

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%

Last extracted: 8/12/2026, 2:59:37 AM

Summary

This study investigates whether AI forecast ensembles, specifically diffusion models, correctly sample the conditional distribution of outcomes. Using coastal sea level forecasts and Lorenz-96 experiments, the authors find that while AI models exhibit positive marginal skill, they fail to capture joint spatial structures, performing worse than climatological draws. This failure is detected by the variogram score but masked by the energy score, indicating a structural inadequacy in learned emulators that persists regardless of training data volume.

Entities (8)

Diffusion Model · model · 98%Energy Score · metric · 97%Variogram Score · metric · 97%Lorenz-96 · system · 96%Coastal Sea Level · physical-variable · 95%Deterministic Emulator · model · 94%Dynamical Ensemble · model · 94%US East Coast · location · 92%

Relation Signals (7)

Diffusion Model → evaluatedby → Energy Score

confidence 97% · failure is invisible to the energy score but detected by the variogram score

Diffusion Model → evaluatedby → Variogram Score

confidence 97% · failure is invisible to the energy score but detected by the variogram score

Diffusion Model → exhibitsfailurein → Joint Spatial Structure

confidence 96% · joint spatial structure is worse than climatological draws

Diffusion Model → produces → Coastal Sea Level

confidence 95% · We train a diffusion model for probabilistic subseasonal coastal sea level forecasts

Dynamical Ensemble → doesnotexhibit → Joint Spatial Structure Failure

confidence 94% · A dynamical ensemble does not replicate the failure

Deterministic Emulator → exhibits → Joint Spatial Structure Failure

confidence 94% · a deterministic emulator does, suggesting it is specific to learned emulators

Lorenz-96 → usedfor → Sensitivity Analysis

confidence 93% · Lorenz-96 experiments across 0.7-170 equivalent years show the gap persists

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Ensemble forecasting aims to sample the conditional distribution of outcomes; whether AI forecast ensembles do this correctly in a joint sense remains largely untested. We train a diffusion model for probabilistic subseasonal coastal sea level forecasts at eight US East Coast tide gauge stations, with sea level derived from reanalysis, and find that marginal and joint forecast quality decouple: positive skill at every station and lead time marginally, while joint spatial structure is worse than climatological draws. A shuffle-based permutation decomposition reveals this failure is invisible to the energy score but detected by the variogram score. Lorenz-96 experiments across 0.7-170 equivalent years show the gap persists regardless of training volume and is reproduced by a linear baseline, indicating structural inadequacy of the learned distribution. A dynamical ensemble does not replicate the failure while a deterministic emulator does, suggesting it is specific to learned emulators rather than ensemble forecasting generally.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

40,363 characters extracted from source content.

Expand or collapse full text

manuscript submitted to Geophysical Research Letters Do AI Forecast Ensembles Sample the Correct Conditional Distribution? Lucas J. Howard 1 , Elizabeth A. Barnes 1 1 Boston University, Faculty of Computing and Data Sciences Key points • A trained probabilistic AI S2S coastal SSH forecast has positive marginal skill but negative joint skill. • The skill gap persists and is insensitive to training dataset volume in idealized ex- periments. • A deterministic emulator reproduces the failure mode while a dynamical model does not. –1– arXiv:2608.08954v1 [physics.ao-ph] 9 Aug 2026 manuscript submitted to Geophysical Research Letters Abstract Ensemble forecasting aims to sample the conditional distribution of outcomes; whether AI forecast ensembles do this correctly in a joint sense remains largely untested. We train a diffusion model for probabilistic subseasonal coastal sea level forecasts at eight US East Coast tide gauge stations, with sea level derived from reanalysis, and find that marginal and joint forecast quality decouple: positive skill at every station and lead time marginally, while joint spatial structure is worse than climatological draws. A shuffle-based permu- tation decomposition reveals this failure is invisible to the energy score but detected by the variogram score. Lorenz-96 experiments across 0.7-170 equivalent years show the gap persists regardless of training volume and is reproduced by a linear baseline, indicating structural inadequacy of the learned distribution. A dynamical ensemble does not repli- cate the failure while a deterministic emulator does, suggesting it is specific to learned emulators rather than ensemble forecasting generally. Plain Language Summary Forecasts of Earth’s oceans and atmosphere are inherently uncertain. These systems are chaotic, our estimates of their current state are imprecise, and our models’ representations of the underlying physics are imperfect. Reliably quantifying the uncertainty of a prediction is therefore a key component of modern forecasting. Ensemble forecasts, in which multiple plausible predictions are independently generated, are a common tool for this purpose. The rapid advancement of AI forecast models, which require far fewer computational resources than traditional models, allows for much larger ensembles and potentially better uncertainty estimates. But a larger ensemble is only useful if it samples the right distribution, i.e. if its spread honestly reflects true uncertainty. We train an AI model to forecast weekly sea levels at eight tide gauge stations along the US East Coast and find that while the model produces skillful forecasts at each individual station, its predictions of how stations are correlated with one another are worse than random historical draws. Experiments in a simplified representation of the atmosphere show that adding more training data does not fix this gap, and that it appears in AI-based models but not traditional physics-based models suggesting a structural limitation of data-driven forecasting. –2– manuscript submitted to Geophysical Research Letters 1 Introduction The fundamental problem that ensemble forecasting attempts to (approximately) solve is to generate samples from a conditional distribution that represents the true un- certainty in outcomes (Leutbecher & Palmer, 2008). Examples of conditions could be the initial state, observations, or climate mode – but in all cases the final ensemble ide- ally contains information about the ways and extent to which forecasts are constrained by available information (Gneiting & Katzfuss, 2014). With ensemble size generally lim- ited by available computational resources, the proliferation of efficient AI forecast tools offers the potential of much larger ensembles and higher fidelity estimates of uncertainty (Price et al., 2023; Weyn et al., 2021). In this work, we explore a foundational question for generative AI forecast ensembles, namely: what distribution are they truly represent- ing, and how does it relate to actual uncertainty? Earth system forecast models provide specific and actionable information to stake- holders on time scales ranging from weather to climate (Bauer et al., 2015). They un- dergird early warning systems that save lives and protect property, they inform guidance for key sectors of the global economy including transportation and agriculture, and they project future climate scenarios to support policy development and implementation. These models have traditionally been numerical solutions of the relevant governing equations solved on a discrete grid, which on the spatial and temporal scales involved are highly computationally expensive (Palmer, 2019). Rigorous uncertainty estimation via ensem- bles has therefore been challenging, with ensemble sizes limited by the cost of a single deterministic simulation. Rapid advances in the use of Artificial Intelligence (AI) meth- ods in earth system science have changed this. Full AI weather emulators and hybrid sys- tems in which only subgrid processes are AI-based have both demonstrated impressive forecast accuracy, competitive with or even surpassing traditional models on commonly used benchmarks (Bi et al., 2022; Lam et al., 2022; Price et al., 2023). Operational fore- cast agencies now run AI forecast models in parallel with their traditional dynamical mod- els (Lang et al., 2024). Even as AI emulators have demonstrated impressive skill (Bi et al., 2022; Lam et al., 2022; Price et al., 2023; Lang et al., 2024), the underlying question of whether ensemble members are sampled from the correct joint distribution remains open. –3– manuscript submitted to Geophysical Research Letters In this work, we show that a generative model trained using reanalysis to predict weekly coastal sea surface height anomalies has positive marginal skill while simultane- ously failing to capture the joint structure of the true uncertainty. Due to the limited size of the reanalysis dataset used for training, we use an idealized system to isolate the impact of training data volume on the observed gap between CRPS skill and variogram score skill. The results of this experiment show that this gap does not close with increas- ing training data volume out to an equivalent of 2X the ERA5 record, suggesting the fail- ure is due to structural inadequacy of the distribution represented by the models out- put rather than insufficient training data. A dynamical ensemble does not replicate the failure mode, while a deterministic emulator using identical initial conditions does, sug- gesting that it is a vulnerability of data-driven approaches not shared by dynamical fore- casts. 2 Methods 2.1 Marginal and Joint Forecast Evaluation To assess whether forecast ensembles are representing the correct distribution in both a marginal and joint sense, we use several evaluation metrics including the contin- uous ranked probability score (CRPS) (Gneiting & Raftery, 2007), energy score (ES) (Gneiting & Katzfuss, 2014), and variogram score (VS) (Scheuerer & Hamill, 2015). CRPS is a proper scoring rule that assesses the marginal distribution of a forecast ensemble at a single lo- cation and lead time; a forecast minimizing CRPS is not required to correctly represent the joint distribution across locations. ES is a strictly proper multivariate extension of CRPS that in principle penalizes joint as well as marginal errors, but has been shown to be relatively insensitive to joint contributions when marginal skill dominates (Scheuerer & Hamill, 2015). VS is an alternative multivariate proper scoring rule (not derived from CRPS) designed to be more sensitive to the dependence structure of a forecast. We use VS with power parameter p=0.5 and equal weights, as recommended by Scheuerer and Hamill (2015) as a convenient default. All scores are expressed as skill scores of the form 1−S f orecast/S r ef , where S f orecast denotes the mean score over the evaluation pe- riod and S r ef is the score of a climatological reference forecast (described in Section 2.2). To isolate the contribution of the joint spatial structure to each score, we apply a shuffle-based permutation decomposition. For a given forecast ensemble, the joint struc- –4– manuscript submitted to Geophysical Research Letters ture is destroyed by independently permuting the ensemble values at each station, which preserves each stations marginal distribution while eliminating cross-station correlations. We apply this shuffle separately to the forecast ensemble and to the climatological ref- erence ensemble. The joint contribution is then quantified as ∆ = S(original)− S(shuffled),(1) where a larger ∆ indicates greater sensitivity of the score to the joint forecast struc- ture. Comparing ∆ forecast to ∆ ref isolates the additional joint contribution of the DDPM above what climatology provides under the same score. 2.2 Diffusion Model SSH Forecast We train a denoising diffusion probabilistic model (DDPM) (Ho et al., 2020) to gen- erate probabilistic SSH anomaly forecasts at eight NOAA National Water Level Obser- vation Network tide gauge stations along the US East Coast: Eastport ME, Boston MA, Atlantic City NJ, The Battery NY, Sewells Point VA, Charleston SC, Mayport FL, and Key West FL. The model outputs the SSH anomaly simultaneously at all eight locations, so the spatial joint distribution is determined entirely by the learned generative model. Data and Preprocessing Ocean state at initialization is represented using GLORYS SSH and sea surface tem- perature (SST) over the western North Atlantic (5-50 ◦ N, 98-50 ◦ W), resampled to weekly means (Lellouche et al., 2021). Atmospheric variability is represented by the North At- lantic Oscillation (NAO) indexdefined as the leading principal component of ERA5 mean sea level pressure over 20 -80 N, 90 W-40 Eand by the 10-m horizontal wind components from ERA5 over the same domain (Kenigson et al., 2018; Thompson & Wallace, 2000). A Ni ̃no3.4 index (mean SST anomaly in 5 ◦ S-5 ◦ N, 170-120 ◦ W) derived from GLORYS is included to utilize potential ENSO teleconnections (Arcodia et al., 2024). All fields are converted to anomalies by a) subtracting the weekly mean computed over the train- ing period (1993-2018) and b) removing the seasonal cycle. Spatial variability in SSH, SST, and wind is compressed using empirical orthogonal functions (EOFs) computed on the training period; we retain the 15 leading SSH modes, 15 SST modes, and 10 wind –5– manuscript submitted to Geophysical Research Letters modes as predictors. SST and SSH components each explain 75% of the variance in the dataset while wind components explain 71%. Model The DDPM conditioning vector is 131-dimensional and includes: the 15 leading SSH principal components at four initialization lags (weeks 0,-1,-2,-4); the 15 leading SST prin- cipal components at three lags (weeks 0, -1,-2); the 10 leading wind principal components at initialization; the NAO index at three lags (weeks 0,-1,-2); the Ni ̃no 3.4 index at four lags (weeks 0, -4,-8,-13); the week 0 GLORYS SSH anomaly at the nearest valid cell to each of the 8 gauge locations; and the forecast lead time which, for training, is one of 2, 4, 8, or 12 weeks normalized to [0,1]. The denoising network is a three-hidden-layer fully connected neural network with layer widths (512, 512, 256). The diffusion process uses a linear noise schedule over T=1000 steps and ensemble members are generated by sampling the learned conditional distribution using the reverse diffusion process. Full train- ing results and hyperparameter tuning can be found in Figure S1 and full architecture details in Table S1. Training and Evaluation The model is trained on weekly mean initialization times from 1993-2018, validated on 2019-2021, and evaluated on a locked test set spanning 2022-2025 (184 initialization times). Training targets are GLORYS12v1 weekly mean SSH anomalies interpolated to each gauge location. An ensemble of 500 members is drawn independently from the con- ditional distribution for each test initialization at lead times of 2, 4, 8, and 12 weeks. The reference forecast for all skill scores is a climatological draw: 500 samples drawn uniformly from the training-period SSH anomaly pool, stratified by calendar week. Skill scores and confidence intervals are computed as described in Section 2.1. Confidence intervals on all skill scores and ∆ values are obtained by a block boot- strap with block length 3 weeks (1000 resamples), which accounts for temporal autocor- relation between successive weekly initialization times. –6– manuscript submitted to Geophysical Research Letters 2.3 Lorenz-96 Experiments The Lorenz-96 (L96) system is a single-variable N-dimensional idealized system de- fined by N coupled ordinary differential equations with cyclic boundary conditions (Lorenz, 2005). It is widely used as a simplified test system in earth system research (Arnold et al., 2013; Brajard et al., 2020; Gagne et al., 2020; L. J. Howard et al., 2024). We use it here to isolate the impact of training data volume and forecast model type on the CRPS/VS skill gap observed in the SSH experiments. We use N=40 and forcing F=8 which are stan- dard choices for which the system is known to be chaotic (Lorenz, 2005). Based on the magnitude of the leading Lyapunov exponent, 1 model time unit (TU) is equivalent to approximately 5 Earth days. The equations are integrated forward using a fourth-order Runge-Kutta scheme with time step dt=0.05 (20 steps per TU). After discarding 100 TU (2,000 steps) as spinup, we save a trajectory of 650,000 steps (32,500 TU). The first 250,000 steps form the train- ing pool, followed by a 100-step gap, 100,000 validation steps, another 100-step gap, and 299,800 test steps. The gaps are included to ensure that each set is fully independent with no data leakage between training, testing, and validation. 500 test cases are drawn randomly without replacement for evaluation from the test set, ensuring that sampled states are well-separated in time and approximately independent. A DDPM is trained to predict the full 40-dimensional system state at lead times of 20, 40, and 60 steps (1, 2, and 3 TU), conditioned on the initial state. Because the L96 domain is cyclic, the denoising network is a convolutional neural network with cir- cular padding (4 residual blocks, 32 hidden channels). Training is repeated for eight train- ing set sizes drawn from the training pool: 1,000; 5,000; 10,000; 20,000; 50,000; 100,000; 200,000; and 249,000 steps (approximately 0.7 to 170 equivalent years; the ERA5 record since 1940 corresponds to approximately 86 equivalent years). For each trained model, an ensemble of 20 members is generated for each of the 500 test initializations. Full ar- chitecture details can be found in Table S1. As a linear baseline, we fit an ordinary least-squares regression with full residual covariance (multivariate normal, MVN) to each training set using the same condition- ing (initial state) and evaluate it with CRPS and VS. Confidence intervals are obtained by paired bootstrap with 1,000 resamples of the 500 test initializations, applied identi- cally to the DDPM and MVN. –7– manuscript submitted to Geophysical Research Letters We additionally compare the DDPM against two reference ensemble forecast types at lead=20 steps (1 TU). The first is a dynamical ensemble generated by integrating the L96 governing equations forward from perturbed initial conditions, where perturbations are drawn from an isotropic Gaussian with standard deviation σ ic reported in units nor- malized to the climatological standard deviation. The second is a deterministic emula- tor trained on the same L96 trajectories and with the same convolutional architecture as the DDPM. Ensemble members are generated by applying this fixed forward map to independently perturbed initial conditions. Both ensembles use 500 members and the same 500 test initializations as the DDPM. We evaluate each at two perturbation lev- els: σ ic = 0.01 to minimize initial-condition uncertainty and isolate the effect of model structure, and a spread-matched σ ic calibrated by grid search to match the DDPM en- semble average spread at lead=20, making the absolute skill comparison fair. 3 Results 3.1 Diffusion Model SSH Forecast The DDPM achieves statistically significant positive CRPS skill at all four lead times (Figure 1a), with station-mean skill dropping monotonically with increasing lead time. Energy Score skill is positive for short leads and slightly negative for longer leads (Fig- ure 1b). The negative skill at 8 and 12 weeks is not statistically significant. In contrast, Variogram Score skill is negative at all lead times (Figure 1c), indicating that the joint spatial structure of the DDPM ensemble is worse than a climatological draw despite the positive marginal skill. The shuffle decomposition in Figure 1b)-(d) isolates the contribution of joint en- semble structure to each score. Figure 1b shows the ES decomposition. The joint con- tribution of the DDPM is statistically significant (asterisk) but small and only modestly larger than the joint contribution of the climatological reference, i.e. the DDPMs joint structure provides little additional information beyond what a random climatological draw provides. Figure 1c shows the VS decomposition. Using VS, the joint contribution to DDPM skill is substantially and significantly larger than the skill of a climatological fore- cast at every lead despite the fact that overall skill is negative. Figure 1d shows the lead- pooled joint contributions. Under ES, the DDPM’s joint contribution barely exceeds that of a climatological draw – the red bar is near zero. Under VS, the DDPM’s joint con- –8– manuscript submitted to Geophysical Research Letters tribution substantially exceeds climatology, yet overall VS skill remains negative (Fig- ure 1c), indicating that the DDPM produces joint structure that is detectably non-climatological but still worse than the observed joint distribution. ES and VS detect fundamentally dif- ferent aspects of the joint forecast distribution: ES is dominated by marginal contribu- tions and is largely insensitive to the multivariate structure of the forecast, while VS clearly detects a structural failure of the distribution represented by the ensemble. To understand this result within the context of the underlying correlation struc- ture of the system, we next examine the pairwise correlations from the training set (Fig- ure 2a) as well as the magnitude of the negative VS skill for that pair averaged over the test set (Figure 2b). There are strong observed correlations concentrated among geograph- ically adjacent northern stations (Eastport through Sewells Point) and among the two southernmost stations (Mayport and Key West). Most station pairs have negative VS skill, indicating the DDPM ensemble is worse than climatology for nearly all pairs. The skill deficit is strongly heterogeneous: the Atlantic City-Battery pair is an extreme out- lier, with the highest observed correlation (r=0.98) and the largest VS deficit, consis- tent with the DDPM over-correlating this already highly correlated pair. Southern sta- tion pairs (Mayport, Key West) show smaller deficits and, in several cases, near-zero fail- ure, hinting that large scale drivers of joint (i.e. spatially correlated) variability may be better captured at these locations. This heterogeneity suggests that the failure mode is not caused by the DDPM simply uniformly amplifying or suppressing correlation struc- ture. Straightforward correction using postprocessing is therefore non-trivial. 3.2 Lorenz 96 Results 3.2.1 Sensitivity to Training Volume We next explore the Lorenz-96 experiments to determine dependence on training data volume. At a lead of 5 days (1 TU, approximately 1.7 Lyapunov times), both CRPS and VS skill increase with training volume, but CRPS skill is consistently higher than VS skill and a statistically significant gap between CRPS and VS skill is present at all training volumes tested (Figure 3a,b). At a lead of 10 days (2 TU, approximately 3.3 Lyapunov times), CRPS skill continues to grow with training volume while VS skill re- mains near zero across all training set sizes, demonstrating a near-complete decoupling of marginal and joint forecast improvement (Figure 3c,d). At 15 days (3 TU, approx- –9– manuscript submitted to Geophysical Research Letters 2wk4wk8wk12wk Lead time 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 CRPS skill (a) Marginal CRPS skill 2wk4wk8wk12wk Lead time 0.10 0.05 0.00 0.05 0.10 ES skill (b) Energy score skill DDPM DDPM shuffled GLORYS draw shuffled 2wk4wk8wk12wk Lead time 1.75 1.50 1.25 1.00 0.75 0.50 0.25 0.00 VS skill ( p =0.5) (c) Variogram score skill DDPMDDPM shuffledGLORYS draw shuffled ESVS 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 skill (pooled) * * * (d) Joint structure contribution 0.00 0.02 0.04 * * * ES skill DDPM GLORYS DDPMGLORYS p < 0.05 (block bootstrap) Figure 1.DDPM SSH forecast skill evaluated using three probabilistic metrics with 95% confidence intervals indicated with shading and error bars. (a) CRPS skill score; positive values indicate improvement over a climatological forecast. (b) Energy Score (ES) skill score of the DDPM ensemble (filled circles), the spatially shuffled DDPM ensemble (open squares), and a shuffled GLORYS climatological draw (triangles) (c) As in (b) but for Variogram Score (VS). (d) Lead-pooled joint structure contribution (∆ skill) for ES and VS; inset shows ES on an expanded scale. The DDPM achieves positive CRPS and ES skill but negative VS skill at all leads, demon- strating that standard metrics can mask a substantial joint forecast failure. –10– manuscript submitted to Geophysical Research Letters Eastport, ME Boston, MA The Battery, NY Atlantic City, NJ Sewells Point, VA Charleston, SC Mayport, FL Key West, FL NS (a) 12345678 0.30 0.25 0.20 0.15 0.10 0.05 0.00 0.05 VS skill deficit (b) 0.20.40.60.81.0 Observed Pearson r 0.775 0.750 0.725 0.700 0.00.51.0 Pearson r (GLORYS 1993 2018) 1 2 3 4 5 6 7 8 Figure 2.Pairwise SSH correlation structure and per-pair VS skill deficit across the eight tide gauge stations. (a) Arc diagram of pairwise Pearson r between stations ordered north (left) to south (right); arc color and thickness indicate correlation strength (GLORYS 1993-2018 training period). Inset shows station locations along the US East Coast. (b) The multivariate skill deficit vs. the observed correlation between stations. Skill deficit is highly heterogeneous across station correlation strength. –11– manuscript submitted to Geophysical Research Letters imately 5.0 Lyapunov times), both scores are near zero at all training volumes, consis- tent with this lead exceeding the effective predictability horizon of the system (Figure 3e,f). The MVN baselinean OLS regression with full residual covariance fit to the same training dataexhibits the same qualitative pattern at every lead and training volume: CRPS skill tracks closely with the DDPM while VS skill remains near zero or below. The ab- sence of learnable joint signal is not specific to the nonlinear DDPM architecture; the much simpler linear regression recovers the same marginal skill compared to climatol- ogy without recovering joint calibration with increasing data volume. The CRPS-VS de- coupling observed in the SSH forecast therefore should not be dismissed as a consequence of insufficient training data. 3.2.2 Dynamical Model and Deterministic AI Emulator To identify the structural source of the CRPS-VS gap, we compare the DDPM against dynamical and deterministic emulator ensembles at lead=20 steps (1 TU, Figure 4). Ad- ditional results for both forecasts are shown in Figure S2 and S3. In the quasi-linear dy- namical regime (σ ic = 0.01), the VS/CRPS ratio is approximately 1.0 (ratio = 1.02), con- firming that chaotic divergence alone does not produce a joint calibration deficit. When σ ic is tuned to match the DDPM ensemble spread (σ ic = 0.158), the dynamical ratio falls modestly to 0.91, indicating that nonlinear trajectory divergence introduces a small joint deficit but cannot account for the full gap observed in the DDPM (ratio = 0.75). The deterministic emulator performs worst at both noise levels (VS/CRPS 0.42-0.58), con- sistent with a deterministic forward map collapsing ensemble diversity in a way that dis- proportionately harms joint relative to marginal calibration. Both AI-based methods show a substantially lower VS/CRPS ratio than the dynamical ensemble at matched spread, suggesting the joint calibration gap is a property of data-driven emulation rather than a consequence of ensemble spread or nonlinear divergence alone. 4 Discussion and Conclusions We have shown that a diffusion model trained for subseasonal coastal sea level pre- diction achieves positive CRPS skill while producing joint spatial forecast structure that is systematically worse than climatology, and that energy score alone does not reveal this –12– manuscript submitted to Geophysical Research Letters 0.0 0.1 0.2 0.3 0.4 Skill score (a) 0.02 0.00 0.02 0.04 0.06 0.08 (c) 0.08 0.06 0.04 0.02 0.00 0.02 (e) 1 yr3 yr7 yr 10 yr30 yr70 yr 100 yr 0.02 0.00 0.02 0.04 0.06 0.08 0.10 0.12 CRPS VS skill (b) 1 yr3 yr7 yr 10 yr30 yr70 yr 100 yr Equivalent training record (yr) 0.02 0.00 0.02 0.04 0.06 0.08 0.10 0.12 (d) 1 yr3 yr7 yr 10 yr30 yr70 yr 100 yr 0.02 0.00 0.02 0.04 0.06 0.08 0.10 0.12 (f) 1k5k10k20k50k100k200k1k5k10k20k50k100k200k Training samples N 1k5k10k20k50k100k200k Lead = 1.0 TU / 5 daysLead = 2.0 TU / 10 daysLead = 3.0 TU / 15 days CRPS skill (DDPM) CRPS skill (MVN) VS skill (DDPM) VS skill (MVN) ERA5 (1940present) Sig. 0 (95%) Not sig. Figure 3.CRPS and Variogram Score (VS) skill as a function of training data volume for a DDPM trained on the Lorenz-96 system at three forecast lead times. (a, c, e) CRPS skill (blue) and VS skill (orange) scores with 95% confidence intervals as a function of training period length (x-axis, in equivalent years) for lead times of 5, 10, and 15 days, respectively. The ERA5 equiv- alent record length is indicated by a solid green vertical line; the skill of a multivariate normal (MVN) baseline is shown for comparison (dotted line). (b, d, f) Violin plots show distributions of the CRPS-VS skill difference at each training set size for the same lead times; dark fill indicates the gap is significantly different from zero (95% bootstrap CI). A statistically significant positive gap (CRPS> VS) is present at leads of 5 and 10 days across all training volumes and shows no systematic decline as training data increases to 2× the ERA5 equivalent record, indicating that the marginal-joint skill discrepancy is not resolved by additional training data under the condi- tions tested. –13– manuscript submitted to Geophysical Research Letters DDPMDyn. ic =0.01 Dyn. ic =0.158 Det. ic =0.01 Det. ic =0.158 0.0 0.2 0.4 0.6 0.8 1.0 Skill (a) CRPS skill VS skill DDPMDyn. ic =0.01 Dyn. ic =0.158 Det. ic =0.01 Det. ic =0.158 0.0 0.2 0.4 0.6 0.8 1.0 1.2 VS/CRPS ratio 0.75 1.02 0.91 0.42 0.58 VS/CRPS = 1 (b) Figure 4.CRPS (blue) and VS (orange) skill scores for each of the ensemble modeling ap- proaches (a). The DDPM is natively probabilistic, while the dynamical and deterministic emu- lators use two sets of perturbed initial conditions to generate ensemble forecasts. The first uses independently and identically distributed Gaussian noise with a standard deviation equal to 0.01 times the climatological standard deviation. The second uses a standard deviation tuned so that the dynamical and DDPM models have approximately equal average ensemble spread at the chosen forecast lead time of 1 TU (5 days). 95% confidence intervals are included for all bars using block bootstrapping. (b) The ratio of VS to CRPS skill for all five experiments. A hori- zontal dashed line at y=1 corresponds to no gap in marginal and joint skill. The VS/CRPS gap is present in both AI-based methods at both noise levels but absent or substantially reduced in the dynamical forecast, indicating the failure is a property of data-driven emulation rather than ensemble forecasting in chaotic systems generally. –14– manuscript submitted to Geophysical Research Letters failure. A shuffle-based permutation decomposition exposes the discrepancy: variogram score skill is negative at all leads while energy score skill is positive, and the per-pair var- iogram score attribution shows highly spatially heterogeneous decoupling of marginal and joint skill. The joint contribution to ES skill is small and only modestly above what a climatological draw provides, while VS detects a substantially larger joint contribution, but one that is still miscalibrated. Lorenz-96 experiments across training volumes equiv- alent to 0.7-170 years of ERA5 demonstrate that the CRPSvariogram score gap does not close with additional data and is shared by linear and nonlinear forecast models alike; the MVN (which by design captures only the linear signal) exhibits the same insensitiv- ity to training data quantity, establishing that the decoupling is unlikely to be a data- volume artifact. Last, a direct comparison against ensembles produced with a dynamical model and a deterministic emulator using perturbed initial conditions shows the skill gap exists in the emulator but not in the dynamical model. Taken together, these experiments demon- strate that 1) probabilistic generative AI forecast tools may sample from a distribution that does not properly represent the joint conditional structure of the true uncertainty, 2) additional training data cannot be assumed to remedy this failure mode, 3) CRPS and ES used as evaluation or training metrics (Perkins et al., 2026) will not reliably identify or mitigate this failure mode, and 4) the failure mode is reproduced for a deterministic AI emulator but not a dynamical model. Supporting the relevance of our findings, recent work on integrating AI forecast mod- els into data assimilation systems has identified covariance errors as a fundamental lim- itation. Tian et al. (2026) showed that the tangent linear trajectories of deterministic AI models have exaggerated sensitivities and predicted that this would manifest as dis- torted background error covariances. Slivinski et al. (2025) confirmed this prediction em- pirically, finding that cycling EnKF with AI models produced covariance amplification sufficient to cause forecast divergence and instabilities. Our results are consistent with this overall picture of AI emulators producing ensembles that are miscalibrated in a mul- tivariate sense even as they compare favorably on standard skill metrics. AI forecast ensembles may appear skillful while simultaneously representing the mul- tivariate distribution incorrectly with no guarantee this failure will improve with increased training data quantity. This strongly supports the continued use of dynamical and physics- –15– manuscript submitted to Geophysical Research Letters based models to complement AI forecasts, and suggests that changes to both the prac- tice and techniques of model evaluation and training are needed to fully realize the po- tential of AI forecast tools. AI Disclosure Claude Code (Anthropic) was used to assist in generating the code for the exper- iments described in this project, as well as scripts for creating publication figures. AI tools, including Claude (Anthropic) and Chat GPT (OpenAI) were also used to assist with manuscript editing. The authors reviewed and verified the content of both the as- sociated code repository and manuscript text and take responsibility for their accuracy and correctness. Acknowledgments This work was supported, in part, by NOAA grants NA24OARX431C0022 and NA22OAR4310621. The computational work reported in this paper was performed using the Shared Com- puting Cluster which is administered by Boston Universitys Research Computing Ser- vices. Open Research Statement Raw GLORYS SSH and SST and ERA5 wind and SLP data are available from their originating agencies. All code for analysis, Lorenz 96 experiments, and plotting along with processed GLORYS and ERA5 data used in this paper are publicly available (L. Howard, 2026). References Arcodia, M. C., Becker, E., & Kirtman, B. P.(2024, February).Subseasonal Variability of U.S. Coastal Sea Level from MJO and ENSO Teleconnection Interference.Weather and Forecasting , 39 (2), 441–458.Retrieved 2026- 08-03, from https://journals.ametsoc.org/view/journals/wefo/39/2/ WAF-D-23-0002.1.xml doi: 10.1175/WAF-D-23-0002.1 Arnold, H. M., Moroz, I. M., & Palmer, T. N.(2013, May).Stochastic parametrizations and model uncertainty in the Lorenz 96 system.Philo- –16– manuscript submitted to Geophysical Research Letters sophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371 (1991), 20110479.Retrieved 2021-10-12, from https://royalsocietypublishing.org/doi/10.1098/rsta.2011.0479 doi: 10.1098/rsta.2011.0479 Bauer, P., Thorpe, A., & Brunet, G. (2015, September). The quiet revolution of nu- merical weather prediction. Nature, 525 (7567), 47–55. Retrieved 2022-12-15, from https://w.nature.com/articles/nature14956 (Number: 7567) doi: 10.1038/nature14956 Bi, K., Xie, L., Zhang, H., Chen, X., Gu, X., & Tian, Q.(2022).PanguWeather: A 3D highresolution model for fast and accurate global weather forecast. arXiv preprint arXiv:2211.02556 . Brajard, J., Carrassi, A., Bocquet, M., & Bertino, L. (2020, July). Combining data assimilation and machine learning to emulate a dynamical model from sparse and noisy observations: A case study with the Lorenz 96 model.Journal of Computational Science, 44 , 101171.Retrieved 2022-10-11, from https:// w.sciencedirect.com/science/article/pii/S1877750320304725doi: 10.1016/j.jocs.2020.101171 Gagne, D. J., Christensen, H. M., Subramanian, A. C., & Monahan, A. H.(2020). Machine Learning for Stochastic Parameterization: Generative Adversar- ial Networks in the Lorenz ’96 Model.Journal of Advances in Modeling Earth Systems, 12 (3), e2019MS001896. Retrieved 2021-09-03, from https://agupubs .onlinelibrary.wiley.com/doi/abs/10.1029/2019MS001896( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019MS001896) doi: 10.1029/2019MS001896 Gneiting, T., & Katzfuss, M.(2014).Probabilistic Forecasting.Annual Re- view of Statistics and Its Application, 1 (1), 125–151.Retrieved 2023-02-23, from https://doi.org/10.1146/annurev-statistics-062713-085831 (eprint: https://doi.org/10.1146/annurev-statistics-062713-085831)doi: 10.1146/annurev-statistics-062713-085831 Gneiting, T., & Raftery, A. E.(2007, March).Strictly Proper Scoring Rules, Pre- diction, and Estimation.Journal of the American Statistical Association, 102 (477), 359–378.Retrieved 2022-12-26, from https://doi.org/10.1198/ 016214506000001437( eprint: https://doi.org/10.1198/016214506000001437) –17– manuscript submitted to Geophysical Research Letters doi: 10.1198/016214506000001437 Ho, J., Jain, A., & Abbeel, P.(2020, December).Denoising Diffusion Probabilis- tic Models. arXiv. Retrieved 2026-06-10, from http://arxiv.org/abs/2006 .11239 (arXiv:2006.11239 [cs.LG]) doi: 10.48550/arXiv.2006.11239 Howard, L. (2026, August). Do AI forecast ensembles sample the correct conditional distribution? code repository. Zenodo. Retrieved from https://zenodo.org/ records/21864536 doi: 10.5281/zenodo.21864536 Howard, L. J., Subramanian, A., & Hoteit, I.(2024).A Machine Learn- ing Augmented Data Assimilation Method for High-Resolution Ob- servations.Journal of Advances in Modeling Earth Systems, 16 (1), e2023MS003774.Retrieved 2024-01-31, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2023MS003774( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2023MS003774) doi: 10.1029/2023MS003774 Kenigson, J. S., Han, W., Rajagopalan, B., Yanto, & Jasinski, M.(2018, July). Decadal Shift of NAO-Linked Interannual Sea Level Variability along the U.S. Northeast Coast.Journal of Climate, 31 (13), 4981–4989.Retrieved 2026- 08-03, from https://journals.ametsoc.org/view/journals/clim/31/13/ jcli-d-17-0403.1.xml doi: 10.1175/JCLI-D-17-0403.1 Lam, R., SanchezGonzalez, ., Willson, M., Wirnsberger, P., Fortunato, M., Alet, F., . . . others (2022). GraphCast: Learning skillful medium-range global weather forecasting. arXiv preprint arXiv:2212.12794 . Lang, S., Alexe, M., Chantry, M., Dramsch, J., Pinault, F., Raoult, B., . . . oth- ers (2024). AIFS ECMWFs datadriven forecasting system.arXiv preprint arXiv:2406.01465 . Lellouche, J.-M., Greiner, E., Bourdall-Badie, R., Garric, G., Melet, A., Drvil- lon, M., . . . Le Traon, P.-Y.(2021, July).The Copernicus Global 1/12 Oceanic and Sea Ice GLORYS12 Reanalysis.Frontiers in Earth Science, 9 .Retrieved 2024-08-28, from https://w.frontiersin.org/journals/ earth-science/articles/10.3389/feart.2021.698876/fulldoi: 10.3389/feart.2021.698876 Leutbecher, M., & Palmer, T. N.(2008, March).Ensemble forecasting.Jour- nal of Computational Physics, 227 (7), 3515–3539.Retrieved 2023-02- –18– manuscript submitted to Geophysical Research Letters 23, from https://w.sciencedirect.com/science/article/pii/ S0021999107000812 doi: 10.1016/j.jcp.2007.02.014 Lorenz, E. N.(2005, May).Designing Chaotic Models.Journal of the Atmo- spheric Sciences, 62 (5), 1574–1587.Retrieved 2022-06-07, from https:// journals.ametsoc.org/view/journals/atsc/62/5/jas3430.1.xmldoi: 10.1175/JAS3430.1 Palmer, T.(2019).The ECMWF ensemble prediction system: Looking back (more than) 25 years and projecting forward 25 years.Quarterly Journal of the Royal Meteorological Society , 145 (S1), 12–24.Retrieved 2026-08- 03, from https://onlinelibrary.wiley.com/doi/abs/10.1002/qj.3383 ( eprint: https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1002/qj.3383)doi: 10.1002/qj.3383 Perkins, W. A., Kwa, A., McGibbon, J., Arcomano, T., Clark, S. K., Watt-Meyer, O., . . . Harris, L. M.(2026, February).HiRO-ACE: Fast and skillful AI emulation and downscaling trained on a 3 km global storm-resolving model. arXiv.Retrieved 2026-08-03, from http://arxiv.org/abs/2512.18224 (arXiv:2512.18224 [physics.ao-ph]) doi: 10.48550/arXiv.2512.18224 Price, I., SanchezGonzalez, ., Alet, F., Andersson, T. R., ElKadi, A., Masters, D., . . . others (2023). GenCast: Diffusionbased ensemble forecasting for mediumrange weather. arXiv preprint arXiv:2312.15796 . Scheuerer, M., & Hamill, T. M.(2015, April).Variogram-Based Proper Scor- ing Rules for Probabilistic Forecasts of Multivariate Quantities.Monthly Weather Review , 143 (4), 1321–1334.Retrieved 2026-05-04, from https:// journals.ametsoc.org/view/journals/mwre/143/4/mwr-d-14-00269.1.xml doi: 10.1175/MWR-D-14-00269.1 Slivinski, L. C., Whitaker, J. S., Frolov, S., Smith, T. A., & Agar- wal, N.(2025).Assimilating Observed Surface Pressure Into ML Weather Prediction Models.Geophysical Research Letters, 52 (6), e2024GL114396.Retrieved 2025-07-13, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2024GL114396( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2024GL114396) doi: 10.1029/2024GL114396 Thompson, D. W. J., & Wallace, J. M. (2000, March). Annular Modes in the Extra- –19– manuscript submitted to Geophysical Research Letters tropical Circulation. Part I: Month-to-Month Variability.Journal of Climate, 13 (5), 1000–1016.Retrieved 2026-08-03, from https://journals.ametsoc .org/view/journals/clim/13/5/1520-044220000131000amitec2.0.co2 .xml doi: 10.1175/1520-0442(2000)013⟨1000:AMITEC⟩2.0.CO;2 Tian, X., Holdaway, D., & Kleist, D.(2026).Evaluating Machine Learn- ing Weather Models for Data Assimilation: Fundamental Limitations in Tangent Linear and Adjoint Properties.Geophysical Research Let- ters, 53 (2), e2025GL119402.Retrieved 2026-06-10, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2025GL119402( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2025GL119402) doi: 10.1029/2025GL119402 Weyn, J. A., Durran, D. R., Caruana, R., & Cresswell-Clay, N.(2021).Sub- Seasonal Forecasting With a Large Ensemble of Deep-Learning Weather Prediction Models.Journal of Advances in Modeling Earth Systems, 13 (7), e2021MS002502.Retrieved 2026-08-04, from https://onlinelibrary.wiley .com/doi/abs/10.1029/2021MS002502 doi: 10.1029/2021MS002502 –20–

Instant research discovery