Continuous Ensemble Weather Forecasting with Diffusion models

## Abstract

Abstract Weather forecasting has seen a shift in methods from numerical simulations to data-driven systems. While initial research in the area focused on deterministic forecasting, recent works have used diffusion models to produce skillful ensemble forecasts. These models are trained on a single forecasting step and rolled out autoregressively. However, they are computationally expensive and accumulate errors for high temporal resolution due to the many rollout steps. We address these limitations with Continuous Ensemble Forecasting, a novel and flexible method for sampling ensemble forecasts in diffusion models. The method can generate temporally consistent ensemble trajectories completely in parallel, with no autoregressive steps. Continuous Ensemble Forecasting can also be combined with autoregressive rollouts to yield forecasts at an arbitrary fine temporal resolution without sacrificing accuracy. We demonstrate that the method achieves competitive results for global weather forecasting with good probabilistic properties.

## 1 Introduction

Forecasting of physical systems over both space and time is a crucial problem with plenty of real-world applications, including in the earth sciences, transportation, and energy systems.
A prime example of this is weather forecasting, which billions of people depend on daily to plan their activities.
Weather forecasting is also crucial for making informed decisions in areas such as agriculture, renewable energy production, and safeguarding communities against extreme weather events.
Current Numerical weather prediction (NWP) systems predict the weather using complex physical models and large supercomputers (Bauer et al., ).
Recently Machine learning weather prediction (MLWP) models have emerged, rivaling the performance of existing NWP systems.
These models are not physics-based but data-driven and have been made possible thanks to advancements in deep learning.
By analyzing patterns from vast amounts of meteorological data (Hersbach et al., ), MLWP models now predict the weather with the same accuracy as global operational NWP models in a fraction of the time (Kurth et al., ; Lam et al., ; Bi et al., ).

Following the success of deterministic MLWP models, the focus of the field has increasingly shifted towards probabilistic modeling.
The probabilistic models generate samples of possible future weather trajectories.
By drawing many such samples, referred to as ensemble members, it is possible to generate an ensemble forecast for quantifying forecast uncertainty and detecting extreme events (Leutbecher & Palmer, ).
Sampling forecasts from deep generative models also address the blurriness often observed in predictions from deterministic MLWP, yielding forecasts that better preserve the variability of the modeled entities.
A popular class of deep generative models used for probabilistic MLWP are diffusion models (Ho et al., ; Price et al., ; Lang et al., ; Shi et al., ).
While these models generate accurate and realistic looking forecasts, they are computationally expensive due to requiring multiple sequential forward passes through the neural network to generate a sample.
Moreover, they are often applied iteratively to roll out longer forecasts (Price et al., ), which exacerbates the computational issue.
Naively switching out this iterative rollout to directly forecasting each future time step does not result in trajectories that are consistent over time.
The auto-regressive rollout additionally puts some limitations on the temporal resolution of the forecast.
Taking too small timesteps results in large accumulation of error over time (Bi et al., ), which forces existing models to resort to a temporal resolution of 12 h Price et al. (); Lang et al. ().
However, in many situations it is crucial to obtain probabilistic forecasts at a much higher temporal resolution. This is true not least when the forecasts are used as decision support in extreme weather situations and for capturing rapidly changing weather events.

We propose a continuous forecasting diffusion model that takes lead time as input and forecasts the future weather state in a single step, while maintaining a temporally consistent trajectory for each ensemble member.(^1^11Our code will be made openly available at publication. [https://github.com/martinandrae/Continuous-Ensemble-Forecasting](https://github.com/martinandrae/Continuous-Ensemble-Forecasting)) This enables both autoregressive and direct forecasting within a single model, improves the accuracy compared to purely autoregressive models at high temporal resolutions, and enables forecasting at arbitrary (non-equidistant) lead times throughout the forecast trajectory. To generate ensemble forecasts, the model uses a lead-time-dependent deterministic ODE-solver to solve the probability flow ODE starting in different pure noise samples. To ensure temporal consistency, we correlate the driving noises for the different lead times, e.g., by using a single noise sample for all timesteps, as illustrated in fig. 1, which enables generating a continuous trajectory for each member. This enables parallel sampling of individual ensemble members without the use of any multi-step loss functions, significantly enhancing the efficiency of ensemble forecasting.

Figure: Figure 1: Diagram of Continuous Ensemble Forecasting.
Refer to caption: extracted/5907711/images/state-t2m-t0.png

#### Contributions.

We propose a novel method for ensemble weather forecasting built on diffusion that: 1) can generate ensemble member trajectories without iteration, 2) can forecast arbitrary lead-times, 3) can be used together with iteration to improve performance on long rollouts, 4) achieves competitive performance in global weather forecasting.

## 2 Related Work

#### Ensembles from perturbations.

In NWP, ensemble forecasts are created by perturbing initial conditions or model parameterizations. The same idea has been applied in MLWP using initial state perturbations (Bi et al., ; Kurth et al., ; Chen et al., ; Li et al., ), model parameter perturbations (Hu et al., ; Weyn et al., ) and Monte Carlo dropout (Scher & Messori, ; Hu et al., ).
In all these cases the underlying model is still deterministic, trained with a mean squared error loss function to predict the average future weather.
These ensembles can thus be viewed as a mixture of means, inheriting the spatial oversmoothing characteristic of deterministic forecasts.

#### Generative models.

Another way of getting ensembles is by modeling the distribution directly through generative models. Price et al. () train a graph-based diffusion model to produce 15-day global forecasts, at 12-hour steps and 0.25°resolution.
Their GenCast model generates realistic forecasts but is slow to sample from, requiring sequential forward passes through the network both for sampling each time step and to roll out the forecast over time. Latent-variable-based MLWP models (Hu et al., ; Oskarsson et al., ; Zhong et al., ) can perform inference in a single forward pass, but have been shown to suffer from more blurriness compared to diffusion models and require delicate hyperparameter tuning (Oskarsson et al., ).
Apart from forecasting, diffusion models have also been successfully applied to other weather-related tasks, including downscaling (Chen et al., ; Mardani et al., ; Pathak et al., ), enlarging physics-based ensembles (Li et al., ) and generating realistic weather from climate scenarios (Bassetti et al., ). However, it is important to note that since the output from these generative models does not directly influence the forecasting process, there are no guarantees that the resulting dynamics will remain continuous over time.
Hua et al. () explore the possibility to incorporate prior information in MLWP diffusion forecasting, by guidance from existing NWP forecasts or climatology, but they do not consider ensemble forecasting.

#### Temporal resolution.

Unlike in NWP, where stability conditions dictate the time step, MLWP models are free to predict at any temporal resolution. Still, the most common approach is to learn to forecast a single short time-step (6h) and iterate this process until the desired lead time (Lam et al., ; Chen et al., ). Although intuitive, this process can lead to error accumulation and is impossible to parallelize due to its sequential nature (Bi et al., ). Multi-step losses have been shown to reduce the error accumulation for deterministic (Lam et al., ) and latent-variable based models (Oskarsson et al., ), but are not trivially implemented in diffusion models. Taking longer timesteps (24h) has been shown to give better results (Couairon et al., ; Bi et al., ), but comes at the loss of temporal resolution.
Bi et al. () resolves this by training multiple models to forecast different lead times, which are then combined in different ways to reach the lead times of interest.
Nguyen et al. (; ) use a similar setup, but train a single model taking the forecast lead time as an input.
This continuous forecasting parallelizes the prediction of the fine temporal scales, but has so far only been applied to deterministic models.
Other approaches have tried to learn fully time-continuous dynamics by using an ODE to generate forecasts (Verma et al., ; Saleem et al., ; Rühling Cachay et al., ; Kochkov et al., ).

#### Spatio-temporal forecasting with diffusion models.

Outside of MLWP, diffusion models have also been applied to forecasting other spatio-temporal processes (Yang et al., ).
Notable examples include turbulent flow simulation (Kohl et al., ; Rühling Cachay et al., ) and PDE solving (Lippe et al., ).
Yang & Sommer () apply diffusion models conditioned on the prediction lead time to a specific floating-smoke fluid field, but do not consider ensemble forecasting.
Another alternative to autoregressive rollouts is the DYfussion framework (Rühling Cachay et al., ), where stochastic interpolation and deterministic forecasting is combined into a diffusion-like model.
The method allows for forecasting at arbitrary temporal resolution, but still requires sequential computations for sampling the prediction.
DYffusion has been successfully applied to climate modeling (Cachay et al., ), but not weather forecasting.

## 3 Background

#### Problem statement.

This paper targets the global weather forecasting problem, an initial-value problem with intrinsic uncertainty. Consider a weather state space $\mathcal{X}$ containing the grid of target variables. Given information about previous weather states $X(\Omega)\subset\mathcal{X}$ at times $\Omega\subset(-\infty,0]$, the aim is to forecast a trajectory $X(\mathcal{T})$ of future weather states $X:(-\infty,\infty)\to\mathcal{X}$ at times $\mathcal{T}\subset(0,T]$ for some time horizon $T$. In particular, the task is to learn and sample from the conditional distribution $p(X(\mathcal{T})|X(\Omega))$.

#### Autoregressive forecasting.

To simplify the problem, $\mathcal{T}$ is often chosen as a set of discrete equally spaced times $\{k\delta\}_{k=1}^{N}$ for some timestep $\delta$. By choosing $\Omega=\{-k\delta\}_{k=0}^{M}$ and assuming $M$th order Markovian dynamics, the joint distribution of $X[k]:=X(k\delta)$ can be factorized over successive states,

$$ $\displaystyle p(X[1{:}N]|X[-M{:}0])=\prod_{k=0}^{N-1}p(X[k{+}1]|X[k{-}M{:}k]),$ (1) $$

and the forecasts can be sampled autoregressively. This way, the network only has to learn to sample a single step $p(X[k{+}1]|X[k{-}M{:}k])$.

#### Conditional Diffusion Models.

Similarly to Price et al. (), we model $p(X[k{+}1]|X[k{-}M{:}k])$ using a conditional diffusion model (Ho et al., ; Song et al., ; Karras et al., ). Diffusion models generate samples by iteratively transforming noise into data. To forecast a future weather state $X[k{+}1]$ given $X[k{-}M{:}k]$ we sample a latent noise variable $Z\sim p_{\text{noise}}$ and iteratively transform it until it resembles a sample from $p(X[k{+}1]|X[k{-}M{:}k])$. We consider the SDE formulation of diffusion models presented by Karras et al. (), which separates training and sampling. Sampling can then be done by solving the probability flow ODE

$$ $\text{d}\bm{z}=-\sigma S_{\theta}(\bm{z};X[k{-}M{:}k],\sigma)\text{d}\sigma$ (2) $$

starting in $Z=\bm{z}(0)\sim p_{\text{noise}}$ and ending in $X[k{+}1]=\bm{z}(1)$. Here $S_{\theta}$ is the neural network trained to match the score function $S$ through the denoising training objective as presented by Karras et al. (). The noise level $\sigma$ is fed to $S_{\theta}$ in each layer as a Fourier embedding similar to Karras et al. ().
Repeating this process for different noise samples $\bm{z}(0)$ gives an ensemble of forecasts.
For a more detailed description of diffusion models see appendix B.

## 4 Continuous Ensemble Forecasting

Iterative forecasting is simple to train but can suffer from error accumulation at long horizons. Consider again the general problem of sampling from $p(X(\mathcal{T})|X(\Omega))$.
In continuous forecasting, a single-step prediction is given by conditioning on the lead time $t\in\mathcal{T}$ and training a conditional score network $S_{\theta}(\bm{z};X(\Omega),t,\sigma)$ to simulate directly from the marginal distribution $p(X({t})|X(\Omega),t)$.
This does not require setting a fixed $\delta$, making the setup more flexible.
The lead time $t$ is added to the conditioning arguments for clarity and can be passed to the network in the same way as the noise level $\sigma$.
While this allows sampling a distribution of states $X(t)$ at each time $t\in\mathcal{T}$, combining these naively does not result in a trajectory $X(\mathcal{T})$. This is because $X(t)$ are samples from the marginal distributions $p(X(t)|X(\Omega),t)$ and not the joint trajectory distribution $p(X(\mathcal{T})|X(\Omega))$.
We propose Continuous Ensemble Forecasting, a novel method of combining samples from $p(X(t)|X(\Omega),t)$ into forecast trajectories $X(\mathcal{T})$ without resorting to autoregressive predictions.

The core idea in our method is to control the source of randomness. To sample from $p(X(t)|X(\Omega),t)$, we sample some noise $Z\sim p_{\text{noise}}$ and feed it to an ODE-solver that solves the probability flow ODE, giving us a forecast $X(t)$ for lead time $t$. Since the ODE-solver is deterministic, the randomness is limited to the noise initialization $Z$. If we freeze the noise, the ODE-solver becomes a deterministic map $f_{\theta}^{Z}:\mathcal{X}^{|\Omega|}\times\mathcal{T}\to\mathcal{X}$ parameterized by the neural network $S_{\theta}$. Applying this map to previous states $X(\Omega)$ and a time $t\in\mathcal{T}$ gives a forecast $X(t)$. Repeating this for several $t_{1},\dots,t_{N}\subset\mathcal{T}$ gives a sequence of forecasts $X(\{t_{i}\}_{i=1}^{N})$. Extending this to all $t\in\mathcal{T}$, we can construct a trajectory $X(\mathcal{T})=f_{\theta}^{Z}(X(\Omega),\mathcal{T})$. We treat this as a sample from $p(X(\mathcal{T})|X(\Omega))$ and propose to use Algorithm 1 to sample such trajectories.

Figure: Algorithm 1 The Continuous Ensemble Forecasting algorithm

### 4.1 Mathematical motivation

In a deterministic system, the dynamics can be described by a forecasting function $f:\mathcal{X}^{|\Omega|}\times\mathcal{T}\to\mathcal{X}$ that maps previous states $X(\Omega)\in\mathcal{X}^{|\Omega|}$, to future states $X(t)$. While weather is in principle governed by deterministic equations, its chaotic nature, lack of information, and our inability to resolve the dynamics at sufficiently high spatio-temporal resolution gives it an intrinsic uncertainty. This motivates the formulation of weather as a stochastic dynamical system. In such a system, no function can describe the entire dynamics. At each instance, there might be a range of functions $f^{1},f^{2},\dots$ that all describe some possible evolution, some more likely than others. If we consider the space $\mathcal{F}$ of all such functions, we can formalize this by defining a probability density $\mu$ over this space, describing the likelihood of each function. To create an ensemble of possible functions, we sample several $f^{1},\dots,f^{N}$, from $\mu$. Given previous states $X(\Omega)$, evaluating each function gives an ensemble of possible trajectories $X^{i}(\mathcal{T})=f^{i}(X(\Omega),\mathcal{T})$.

Figure: Figure 2: Diagram depicting the connection between the latent noise space and the solution space.

The key insight in motivating our method is the identification of the latent noise space $(\mathcal{X},p_{\text{noise}})$ with the solution space $(\mathcal{F},\mu)$.
In our method, the sampling algorithm can be represented by a parameterized function $f_{\theta}^{Z}$, where we have frozen the noise $Z\sim p_{\text{noise}}$. Under some regularity conditions, this function is uniquely defined by $\theta$ and the ODE-solver for any given $Z$. Thus, as illustrated in fig. 2, our setting mirrors the theoretical setting.
Given sufficient data and model capacity, the neural network $S_{\theta}$ matches the score function $S$. Consequently, the distribution of the functions $f_{\theta}^{Z}$ should mirror that of the solutions $f^{i}$.
Since $f^{i}$ describes a possible evolution of weather, it has to be continuous as a function of time.
To ensure that the generated trajectories are also continuous in $t$, regularity conditions need to be imposed on $S_{\theta}$ to ensure that it depends smoothly on the lead time. While we do not impose any restrictions on the network, we show empirically in sec. 5.1 that this is satisfied in practice, and leave a formal analysis for future work.

### 4.2 Autocorrelated noise

As mentioned above, under regularity conditions, we expect the solution map $f_{\theta}^{Z}(X(\Omega),t)$ to be continuous in $t$. A shortcoming of this approach, however, is that it also constrains the trajectories to be conditionally deterministic, conditionally on the solution $f_{\theta}^{Z}(X(\Omega),t^{\star})$ at any fixed time point $t^{\star}\in\mathcal{T}$. Specifically, consider the conditional distribution $p(X(t)\mid X(t^{\star}),X(\Omega),t)$ for $t,t^{\star}\in\mathcal{T}$ under the proposed model.
Conditionally on $X(t^{\star})=f_{\theta}^{Z}(X(\Omega),t^{\star})$ we can, conceptually, invert the ODE which generated $X(t^{\star})$ to recover the driving noise $Z$. Now, if $X(t)=f_{\theta}^{Z}(X(\Omega),t)$ is generated using the exact same noise variable, we find that $p(X(t)\mid X(t^{\star}),X(\Omega),t)$ is a Dirac point mass concentrated at $f_{\theta}^{Z}(X(\Omega),t)$. This violates the assumed stochasticity of the dynamical process that we are modelling.

To address this shortcoming, a simple extension of the proposed model is to replace the fixed $Z$ with a stochastic process $Z(\mathcal{T})$ with continuous sample trajectories. The stochastic process is chosen to be stationary with marginal distribution $Z(t)\sim p_{\text{noise}},\forall t\in\mathcal{T}$. This ensures that the generative model at any fixed $t$ is probabilistically equivalent with the fixed noise setting. Specifically, no changes to the training algorithm are needed since the training is based solely on time marginals. Allowing temporal stochasticity in the driving noise process results in a non-deterministic relationship between the states at different lead times, while we can keep the temporal consistency by ensuring that the driving process is sufficiently autocorrelated. A simple choice is to select $Z(\mathcal{T})$ as an Ornstein-Uhlenbeck process (Gillespie, ),
$dZ(t)=-\rho Z(t)dt+\sqrt{2\rho}dW(t),$
for some correlation parameter $\rho>0$ and $W(t)$ being Brownian motion. This is what we use in Algorithm 2. Alternatively, we can use a Gaussian process with stronger smoothness properties, e.g., with a squared exponential kernel. This could prove useful for formally proving time continuity of the resulting trajectories, but we leave such an analysis for future work.

Figure: Algorithm 2 The Extended Continuous Ensemble Forecasting algorithm

### 4.3 Autoregressive roll-outs with continuous interpolation

Continuous forecasting is effective for forecasting hours to days but can struggle to forecast longer lead times where the correlation is weaker. Iterative forecasting excels at long lead times when used with longer (24h) timesteps (Bi et al., ), but comes at the loss of temporal resolution. In our framework, it becomes possible to sample both iterative and continuous forecasts with the same model. We propose to leverage this by iterating on a longer timestep and forecasting the intermediate timesteps using Continuous Ensemble Forecasting, as outlined in Alg. 3. We refer to this combined method as Autoregressive Rollouts with Continuous Interpolation (ARCI). Our method limits the error accumulation without sacrificing temporal resolution. This allows for producing forecasts at an arbitrary fine temporal resolution, while retaining the accuracy of the best iterative methods throughout the whole forecast. By limiting the number of iterative steps, more of the forecast also becomes parallelizable, allowing for rapidly generating forecasts on large compute clusters.
Furthermore, we can straightforwardly use different time resolutions during different parts of the forecast trajectories, by for instance forecasting with 1h steps for the first few days and then switching to longer time steps for long lead times.
Note that ARCI can be used with either fixed (with Alg. 1) or stochastic (with Alg. 2) driving noise. However, we emphasize that these two algorithms are probabilistically equivalent for all time marginals, and only differ in the autocorrelation of forecast trajectories.

Figure: Algorithm 3 ARCI (Autoregressive roll-outs with continuous interpolation)

## 5 Experiments

#### Data.

We evaluate our method on global weather forecasting up to 10 days at 1, 6, and 24 hour timesteps.
We use the downsampled ERA5 reanalysis dataset (Hersbach et al., ) at 5.625°resolution and 1-hour increments provided by WeatherBench (Rasp et al., ). The models are trained to forecast 5 variables from the ERA5 dataset: geopotential at 500hPa (z500), temperature at 850hPa (t850), ground temperature (t2m) and the ground wind components (u10, v10). The atmospheric fields z500 and t850 offer a comprehensive view of atmospheric dynamics and thermodynamics, while the surface fields t2m and u10, v10 are important for day-to-day activities. We also evaluate the forecast of ground wind speed ws10, computed from the model outputs as $\text{{ws10}}=\sqrt{\text{{u10}}^{2}+\text{{v10}}^{2}}$. This is useful for evaluating how well the methods model cross-variable dependencies. All variables are standardized by subtracting their mean and dividing by their standard deviation. Together with the previous states we also feed the models with static fields. These include the land-sea mask and orography, both rescaled to $[0,1]$. All models are trained on the period 1979–2015, validated on 2016–2017 and tested on 2018. We consider every hour of each year as forecast initialization times, except for the first 24 h and last 10 days in each subset. This guarantees that all times forecasted or conditioned on lie within the specific years.

#### Metrics.

We evaluate the skill of the forecasting models by computing the Root Mean Squared Error (RMSE) of the ensemble mean.
As a probabilistic metric we also consider Continuous Ranked Probability Score (CRPS) (Gneiting & Raftery, ), which measures how well the the predicted marginal distributions capture the ground truth.
We also evaluate the Spread/Skill-Ratio (SSR), which is a common measure of calibration for ensemble forecasts.
For a model with well calibrated uncertainty estimates the SSR should be close to 1 (Fortin et al., ).
Detailed definitions of all metrics are given in appendix A.

#### Models.

We propose to use the ARCI model described in algorithm 3 referred to as ARCI-24/6h. We train it to forecast $t\in\{6,12,18,24\}$ (hours) and roll it out autoregressively with 24 h steps, hence the name. Training is done on a 40GB NVIDIA A100 GPU and takes roughly 2 days. We emphasize again that the choice of algorithm inside ARCI results in probabilistically equivalent forecasts for all time marginals, and only differ in the autocorrelation of forecast trajectories. Hence, for all metrics below that are computed for specific lead times we only report results for one version of the algorithm. We return to the difference between Alg. 1 and Alg. 2 when studying the temporal difference below.

To evaluate the effectiveness of our approach, we compare it to other MLWP baselines.
Deterministic is a deterministic model trained using MSE-loss on a single 6 hour time step, and unrolled up to 10 days.
AR-6/24h is a diffusion model trained only on forecasting a single fixed $\delta$ ahead, and then autoregressively unrolled up to 10 days. This is the forecasting setup of (Price et al., ).
CI-6h is a diffusion model performing continuous forecasting conditioned on a specific lead time. It is trained on uniformly sampled lead times from $\{k\delta\}_{k=1}^{40}$, with $\delta=\text{6 h}$. This is the method proposed in alg. 1.
Sampling a 10-day forecast with 6h resolution for a single member from AR-6h takes 32 seconds, but by parallelizing the 6h timesteps in ARCI-24/6h this reduces to 8 seconds.

To compare against another family of ensemble forecasting models from the literature we retrain the Graph-EFM model (Oskarsson et al., ) on our exact data setup.
Graph-EFM is a graph-based latent-variable model that produces forecasts by 6 h iterative rollout steps.
For all models, unless otherwise specified, we condition on the two previous timesteps $\Omega=\{0,-\delta\}$ and sample 50 ensemble members at each initialization time.
All models except Graph-EFM use the same architecture based on the U-net in Karras et al. () as presented in appendix B.

### 5.1 Results

#### Quantitative results.

Table 1 and figure 3 show metrics for a selection of lead times and variables. Scores for the remaining variables are listed in appendix C. All probabilistic models show a clear improvement over the deterministic model. CI-6h performs well on short-term forecasting but struggles at longer horizons. This is likely due to the challenge of learning any useful relationships between initial states and later lead times, which are weakly correlated. ARCI-24+6h outperforms all models at 6h resolution, including the external baseline Graph-EFM, and matches the best overall model AR-24h in almost all scores. All diffusion-based models have SSR $<$ 1, indicating some systematic underdispersion.

**Table 1: Selection of results for 5 and 10 day forecasting of geopotential at 500 hPa (z500) and temperature at 850 hPa (t850). For RMSE and CRPS, lower values are better, and SSR should be close to 1. The best values are marked with bold and second best underlined.**
|  |  | Lead time 5 days | Lead time 10 days |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Variable | Model | RMSE | CRPS | SSR | RMSE | CRPS | SSR |
| z500 | Deterministic | 766.7 | 483.9 | - | 1042 | 661.5 | - |
|  | Graph-EFM | 699.1 | 317.5 | 1.13 | 817.1 | 373.6 | 1.1 |
|  | AR-6h | 602.3 | 287.8 | 0.75 | 811.8 | 391.9 | 0.88 |
|  | AR-24h | 544.2 | 242.7 | 0.84 | 750.6 | 335.2 | 0.94 |
|  | CI-6h | 707.8 | 321.2 | 0.59 | 885.7 | 406.6 | 0.6 |
|  | ARCI-24/6h | 560.9 | 256.7 | 0.86 | 765.6 | 355.2 | 0.93 |
| t850 | Deterministic | 3.48 | 2.36 | - | 4.54 | 3.17 | - |
|  | Graph-EFM | 3.12 | 1.56 | 1.11 | 3.51 | 1.77 | 1.12 |
|  | AR-6h | 2.72 | 1.34 | 0.82 | 3.39 | 1.69 | 0.92 |
|  | AR-24h | 2.55 | 1.24 | 0.89 | 3.25 | 1.6 | 0.96 |
|  | CI-6h | 3.06 | 1.5 | 0.74 | 3.68 | 1.85 | 0.71 |
|  | ARCI-24/6h | 2.6 | 1.27 | 0.9 | 3.29 | 1.63 | 0.95 |

Figure: Figure 3: RMSE, CRPS and SSR for temperature at 850 hPa (t850) with a selection of models at 6h resolution.
Refer to caption: x1.png

Figure: Figure 4: Example forecasts for temperature at 850hPa (t850) at lead time 10 days. The forecasts are generated using ARCI-24/6h except for Deterministic which is sampled using the deterministic model. The bottom row shows 4 ensemble members, randomly chosen out of the 50.
Refer to caption: extracted/5907711/images/forecast-t850.png

#### Qualitative results.

Figure 4 shows an example forecast from ARCI-24/6h for temperature at 850 hPa (t850) at 10 days lead time.
The forecasts are rich in detail, resembling the true state more than the ensemble mean. Examples of other variables are given in appendix D.

Figure: Figure 5: Temporal difference for temperature on 850 hPa (t850) for different values of $\rho$ in algorithm 2. Choosing $\rho=0$ fixes the noise, $\rho=\text{ln}10$ allows it to vary and $\rho\to\infty$ gives completely uncorrelated noise. The black line refers to the temporal difference of the data.
Refer to caption: x2.png

#### Temporal Difference.

Autoregressively sampled forecast trajectories are necessarily continuous. Since there is no standard way of measuring the continuity of a forecast, we propose using the mean temporal difference ${\Delta X=|X(t)-X(t-1)|}$ as a measure of forecast continuity. Figure 5 shows $\Delta X$ for a CI-1h continuous model trained up to 24 hours with 1-hour timesteps. Compared to using different noise at each step ($\rho\to\infty)$, the temporal difference of our model $(\rho=0)$ stays close to the temporal difference of the data. This supports our claim that continuous ensemble forecasting produces continuous trajectories.

Figure 5 also shows the temporal difference of the continuous model above. When the noise is fixed ($\rho=0$) the temporal difference decreases with lead time, corresponding to predictions with smaller temporal variations. Letting the noise vary with noise factor $(\rho=\text{ln}10)$ as in alg. 2 stops this from happening. The bias between $\Delta X$ of the data and our model is likely due to the model producing slightly blurrier forecasts, making the differences smaller.

#### Continuous Time Forecasting.

Our ARCI method allows for producing forecasts at arbitrary fine temporal resolution while retaining the accuracy of methods taking longer autoregressive steps.
We here demonstrate this by producing hourly forecasts.
Figure 6 shows the scores of a selection of models for z500 for 10-day forecasts at 1h resolution. The AR-1h model has the same setup as AR-6h but on 1h resolution. ARCI-24/1h and ARCI-24/2h* are both trained with $\Omega=\{0,-24\}$ to have access to the same information at each timestep. The autoregressive AR-1h model performs much worse than on 6 or 24 hours. The continuous model, however, does not lose performance by increasing the temporal resolution, making the 1h timestep forecasts as skillful as the 24 hour ones.
An alternative to directly producing forecasts at a fine temporal resolution would be to linearly interpolate the forecasts sampled using an autoregressive model. In fig. 9 in appendix C we show that linearly interpolated forecasts behave much worse than the ARCI model on both 1 and 6-hour resolution.

The ARCI-24/2h* model is trained only on lead times 2h apart (lead times in $\{2k\}_{k=1}^{12}$), but used to forecast each 1h timestep. This showcases the ability of the method to generalize to lead times not considered during training. It performs similarly to the ARCI-24/1h model trained on all timesteps, indicating that it can generalize beyond its training setup to even finer resolutions. For highly time-dependent fields such as t2m, the model performs worse at the first forecast in each iteration (1h, 25h,…), as seen in figure 11 in appendix C.
For other lead-times $t$ not considered during training, the model has trained on forecasting $t-1$ and $t+1$, thus only having to interpolate to $t$. However, since we do not train on forecasting 0h, the model instead has to extrapolate to $t=\text{1h}$ what was learned for 2h forecasts. This issue could possibly be fixed by letting the network also train on 0h forecasts.

Figure: Figure 6: RMSE, CRPS, and SSR for geopotential at 500 hPa (z500) with a selection of models at 1h resolution.
Refer to caption: x3.png

## 6 Conclusion

We present Continuous Ensemble Forecasting, a novel framework for probabilistic MLWP that increases the efficiency, accuracy, and flexibility of weather forecasts at high temporal resolution. When combined with autoregressive prediction, our ARCI method produces 10-day forecasts with a 1-hour resolution that matches the accuracy of a purely autoregressive model with 24-hour steps.
With this work, we hope to show that the possibilities with generative modeling for spatio-temporal predictions are still largely unexplored and a fruitful area of research.

#### Limitations.

While our proposed framework achieves good results on 5.625°Weatherbench (Rasp et al., ) data, we have yet to show that the method scales to problems with higher spatial resolution. Additionally, as the lead time increases, the correlation between initial and future states becomes weaker, limiting the application of continuous forecasting. While our method parallelizes more of the sampling than previous autoregressive models, solving the probability flow ODE in eq. 2 still requires many sequential forward passes through the network. Sampling is thus still slower than for latent variable models, but the predicted distribution more accurate.
Finally, due to computational limitations, we were not able to retrain our models several times and as such do not report error bars in our results.

#### Future work.

One interesting direction for future work is a further investigation of autocorrelated noise, in particular, how the choice of stochastic process can aid in producing continuous trajectories with a stationary temporal difference.
This includes correlating the noise in the autoregressive steps with the continuous steps, which could help ease the transition between them.
Another idea is to take the direction of DYffusion (Rühling Cachay et al., ) and directly adjust the diffusion objective to better suit temporal data.
While we have demonstrated continuous ensemble forecasting for weather, the idea is generally applicable and it would also be of interest to apply it to other spatio-temporal forecasting problems.

## Acknowledgments

This work was partially supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation.
Our computations were enabled by the Berzelius resource at the National Supercomputer Centre, provided by the Knut and Alice Wallenberg Foundation.

## Appendix A Metrics

We consider the following evaluation metrics used to assess the probabilistic forecasts produced by the diffusion model. Metrics from meteorology and general uncertainty quantification, such as RMSE, Spread/Skill ratio (SSR), and Continuous Ranked Probability Score (CRPS) are employed to measure the effectiveness and reliability of the model outputs. All of our metrics are weighted by latitude-dependent weights. For a particular variable and lead time,

- •
$x^{k}_{i,n}$ represents the value of the $k$-th ensemble member at initialization time indexed by $n=1\dots N$ for grid cell in the latitude and longitude grid indexed by $i\in I$.
- •
$y_{i,n}$ denotes the corresponding ground truth.
- •
$\bar{x}_{i,n}$ denotes the ensemble mean, defined by $\bar{x}_{i,n}=\frac{1}{n_{\text{ens}}}\sum_{k=1}^{n_{\text{ens}}}x_{i,n}^{k}$.
- •
$a_{i}$ denotes the area of the latitude-longitude grid cell, which varies by latitude and is normalized to unit mean over the grid (Rasp et al., ).

RMSE or skill measures the accuracy of the forecast. Following Rasp et al. () we define the RMSE as the mean square root of the ensemble mean:

$$ $\displaystyle\text{RMSE}:=\frac{1}{N}\sum_{n=1}^{N}\sqrt{\frac{1}{|I|}\sum_{i \in I}a_{i}(y_{i,n}-\bar{x}_{i,n})^{2}}.$ (3) $$

In the case of deterministic predictions, the ensemble mean is taken as the deterministic prediction.

Spread represents the variability within the ensemble and is calculated as the root mean square of the ensemble variance:

$$ $\displaystyle\text{Spread}:=\frac{1}{N}\sum_{n=1}^{N}\sqrt{\frac{1}{|I|}\sum_{ i\in I}\frac{1}{n_{\text{ens}}-1}\sum_{k=1}^{n_{\text{ens}}}a_{i}(x_{i,n}^{k}- \bar{x}_{i,n})^{2}}.$ (4) $$

Ideally, the forecast achieves a balance where skill and spread are proportional, leading to an optimal spread/skill ratio (SSR) close to 1, indicating effective uncertainty estimation:

$$ $\displaystyle\text{SSR}:=\sqrt{\frac{n_{\text{ens}}+1}{n_{\text{ens}}}}\frac{ \text{Spread}}{\text{RMSE}}.$ (5) $$

Continuous Ranked Probability Score (CRPS) (Gneiting & Raftery, ) measures the accuracy of probabilistic forecasts by comparing the cumulative distribution functions (CDFs) of the predicted and observed values. It integrates the squared differences between these CDFs, providing a single score that penalizes differences in location, spread, and shape of the distributions. An estimator of the CRPS is given by:

$$ $\text{CRPS}:=\frac{1}{N}\sum_{n=1}^{N}\frac{1}{|I|}\sum_{i\in I}a_{i}\left( \frac{1}{n_{\text{ens}}}\sum_{k=1}^{n_{\text{ens}}}|x_{i,n}^{k}-y_{i,n}|-\frac {1}{2n_{\text{ens}}^{2}}\sum_{k=1}^{n_{\text{ens}}}\sum_{k^{\prime}=1}^{n_{ \text{ens}}}|x_{i,n}^{k}-x_{i,n}^{k^{\prime}}|\right).$ $$

Temporal Difference measures the mean absolute difference between states at consecutive times. It’s used to measure the continuity of a forecast. For forecasts $x_{i,n}^{k},\hat{x}_{i,n}^{k}$ at consecutive lead times, it is given by:

$$ $\displaystyle\Delta X:=\frac{1}{N}\sum_{n=1}^{N}\frac{1}{n_{\text{ens}}}\sum_{ k=1}^{n_{\text{ens}}}\frac{1}{|I|}\sum_{i\in I}a_{i}|x_{i,n}^{k}-\hat{x}_{i,n} ^{k}|.$ (6) $$

## Appendix B Model

#### Preconditioning

The sampling process is based on the denoising neural network $D_{\theta}$ that takes a noisy residual and tries to denoise it. To help in this, it is also given the noise level $\sigma$, the previous state $X(\Omega)$ and the lead time $t$. To simplify learning, $D_{\theta}$ is parameterized by a different network $F_{\theta}$ defined by

$$ $D_{\theta}(X_{t}^{\sigma};\sigma,X(\Omega),t)=c_{\text{skip}}(\sigma)\cdot X_{ t}^{\sigma}+c_{\text{out}}(\sigma)\cdot F_{\theta}\left(c_{\text{in}}(\sigma) \cdot X_{t}^{\sigma};c_{\text{noise}}(\sigma),X(\Omega),t\right),$ $$

where $X_{t}^{\sigma}$ denotes a noisy version of the target time $X({t})$ at noise level $\sigma$, and $c_{\text{skip}}$, $c_{\text{out}}$, $c_{\text{in}}$ and $c_{\text{noise}}$ are scaling functions taken from (Karras et al., ) defined in Tab. 2. These scaling functions $c_{\text{skip}}$, $c_{\text{out}}$, $c_{\text{in}}$ and $c_{\text{noise}}$, are specifically chosen to handle the influence of the noise level within the network, allowing $D_{\theta}$ to adapt dynamically to different noise intensities without the need for adjusting the scale of $\sigma$ externally. Consequently, for consistency with the normalization of the data where $\sigma_{\text{data}}$ is set to 1, the lead time $t$ is also scaled to fit within the range $[0,1]$. This normalization ensures that the network inputs are uniformly scaled, enhancing the efficiency and effectiveness of the denoising process.

**Table 2: Scaling functions.**
| Skip scaling | $c_{\text{skip}}(\sigma)$ | $\sigma^{2}_{\text{data}}/(\sigma^{2}+\sigma^{2}_{\text{data}})$ |
| --- | --- | --- |
| Output scaling | $c_{\text{out}}(\sigma)$ | $\sigma\cdot\sigma_{\text{data}}/\sqrt{\sigma^{2}+\sigma^{2}_{\text{data}}}$ |
| Input scaling | $c_{\text{in}}(\sigma)$ | $1/\sqrt{\sigma^{2}+\sigma^{2}_{\text{data}}}$ |
| Noise scaling | $c_{\text{noise}}(\sigma)$ | $\frac{1}{4}\ln(\sigma)$ |

#### Conditioning

To condition on the initial conditions $X(\Omega)$ and static fields, these are concatenated along the channel dimension with the input to the denoiser, increasing the dimension of the input. To condition on the noise level $\sigma$ and lead time $t$, we use Fourier embedding as specified in (Karras et al., ). Fourier embedding captures periodic patterns in noise and time, enhancing the model’s ability to handle complex time-series dependencies effectively. They work by transforming the time/noise into a vector of sine/cosine features at 32 frequencies with period 16. These vectors are added and then passed through two fully connected layers with SiLU activation to obtain a 128-dimensional encoding.

#### Architecture

The backbone of the diffusion model is a U-Net architecture. Our model is based on the one used in (Karras et al., ), reconfigured for our purposes with 32 filters as the base multiplier. It is built up by blocks configured as in fig. 7. The blocks consist of two convolutional layers and are constructed as in fig. 8. If the block is a down-/up-sample or if the number of input filters is different from the number of output filters, there is an additional skip layer from the input to the output. The time/noise embedding is fed directly into each block and not passed through the network. Unlike the network in Karras et al. (), our convolutions uses zero padding at the poles and periodic padding at the left/right edges. This periodic padding ensures periodicity over longitudes. The model has 3.5M parameters.

Figure: Figure 7: Overview of the U-Net Architecture, detailing layer configurations and the flow of information through convolutional blocks and skip connections.
Refer to caption: x4.png

Figure: Figure 8: Construction of a Diffusion Model Block, showing the sequence of operations and the integration of embeddings with add connections.
Refer to caption: x5.png

#### Sampling

To generate forecasts using our diffusion model, we solve the probability flow ODE as defined in (Karras et al., )

$$ $\text{d}\bm{z}=-\dot{\sigma}(s)\sigma(s)\nabla_{\bm{z}}\log p_{s}\left(\bm{z} \right)\text{d}s.$ (7) $$

We employ the second-order Heun’s method, a deterministic ODE solver, as outlined in Algorithm 4. For the noise parameters, we define the noise level function as $\sigma(s)=s$. Additionally, we set a noise level schedule to lower the noise during sampling from $\sigma_{\max}$ to $\sigma_{\min}$ over $N$ steps:

$$ $s_{i}=\left(\sigma_{\max}^{\frac{1}{\rho}}+\frac{i}{N-1}\left(\sigma_{\min}^{ \frac{1}{\rho}}-\sigma_{\max}^{\frac{1}{\rho}}\right)\right)^{\rho},\quad i\in \{0,\dots,N-1\}.$ $$

The relevant parameters for training and sampling are given in tab. 3.

Figure: Algorithm 4 Deterministic sampling using Heun’s $2^{\text{nd}}$ order method.

**Table 3: Parameters used for sampling and training.**
| Name | Notation | Value, sampling | Value, training |
| --- | --- | --- | --- |
| Maximum noise level | $\sigma_{\text{max}}$ | 80 | 88 |
| Minimum noise level | $\sigma_{\text{min}}$ | 0.03 | 0.02 |
| Shape of noise distribution | $\rho$ | 7 | 7 |
| Number of noise levels | $N$ | 20 |  |

#### Training

The dataset is partitioned into three subsets: training, validation, and testing. The training subset is used for model training, the validation subset for evaluating generalization, and the testing subset to determine final accuracy. The diffusion model is trained using the following training objective

$$ $\mathbb{E}_{t\sim p^{t}}\mathbb{E}_{\sigma\sim p_{\sigma}}\mathbb{E}_{(X( \Omega),X(t))\sim p_{\text{data}}}\mathbb{E}_{\epsilon|\sigma\sim\mathcal{N}(0 ,\sigma^{2}\mathbf{I})}\left[\frac{1}{\sigma^{2}}\sum_{i}\sum_{j}\frac{a_{i}}{ s_{j}(t)}\frac{1}{|I||J|}\left(\hat{X}(t)_{i,j}-X(t)_{i,j}\right)^{2}\right].$ $$

with $\hat{X}(t)=D_{\theta}(X(t)+\epsilon;\sigma,X(\Omega),t)$ and $J$ being the set of variables. Here, $p^{t}$ represents a uniform distribution over the lead times. We have also included a scaling term $s_{j}(t)^{-1}$ which scales the loss by the precomputed standard deviation $s_{j}(t)$ based on lead time $t$ for each variable $j\in J$. This normalization process is designed to weigh short and longer times equally. The noise level distribution $p_{\sigma}$ is chosen to be consistent with the sampling noise level described above. Specifically, its inverse CDF is:

$$ $F^{-1}(u)=\left(\sigma_{\max}^{\frac{1}{\rho}}+u\left(\sigma_{\min}^{\frac{1}{ \rho}}-\sigma_{\max}^{\frac{1}{\rho}}\right)\right)^{\rho},$ $$

and we sample from it by drawing $u\sim U[0,1]$. The training process is executed in Pytorch, with setup and parameters detailed in Tab. 4.

**Table 4: Optimizer Hyperparameters.**
| Optimizer hyperparameters |  |
| --- | --- |
| Optimiser | AdamW (Loshchilov & Hutter, ) |
| Initialization | Xavier Uniform (Glorot & Bengio, ) |
| LR decay schedule | Cosine (Loshchilov & Hutter, ) |
| Peak LR | 5e-4 |
| Weight decay | 0.1 |
| Warmup steps | 1e3 |
| Epochs | 300 |
| Batch size | 256 |
| Dropout probability | 0.1 |

#### Graph-EFM Baseline

For the Graph-EFM baseline we use the same data setup as for the other models.
Since we are working on a coarser resolution than Oskarsson et al. () some adaptation was necessary to the exact graph structure used in the model.
We construct the graph by splitting a global icosahedron 2 times, resulting in 3 hierarchical graph levels.
The training follows the same schedule as in Oskarsson et al. (), with pre-training on single step 6 h prediction and fine-tuning on rollouts.
Details of the training schedule are given in tab. 5.

**Table 5: Training schedule for Graph-EFM, using the notation from Oskarsson et al. ().**
| Epochs | Learning Rate | Unrolling steps | $\lambda_{\text{KL}}$ | $\lambda_{\text{CRPS}}$ |
| --- | --- | --- | --- | --- |
| 20 | $10^{-3}$ | 1 | 0 | 0 |
| 75 | $10^{-3}$ | 1 | 0.1 | 0 |
| 20 | $10^{-4}$ | 4 | 0.1 | 0 |
| 8 | $10^{-4}$ | 8 | 0.1 | $10^{5}$ |

## Appendix C Additional Results

#### Interpolation

An alternative to directly producing forecasts at a fine temporal resolution would be to simply interpolate the forecasts sampled using an autoregressive model. We argue that our ARCI method produces more accurate and realistic predictions than simple linear interpolation. Figure 9 shows the RMSE and Spread/Skill for a linear interpolation of AR-24h, compared to ARCI-24/6h. The interpolated forecasts behave much worse than the ARCI model on the same resolution, which strengthens the role of our model as an advanced interpolator.

Figure: Figure 9: RMSE and SSR for linear interpolations of geopotential at 500 hPa (z500) and ground wind speed (ws10). Linear-1h/Linear-6h are linear interpolations of AR-6h/AR-24h.
Refer to caption: x6.png

Table 6 and figure 10 show the results on 10 day forecasting for all variables. Figure 11 shows the results on forecasting with 1h resolution. Figure 12 shows the temporal difference for different values of $\rho$.

**Table 6: Results on 5 and 10 day forecasting for all variables. For RMSE and CRPS, lower values are better, and SSR should be close to 1.**
|  |  | Lead time 5 days | Lead time 10 days |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Variable | Model | RMSE | CRPS | SSR | RMSE | CRPS | SSR |
| z500 | Deterministic | 766.7 | 483.9 | - | 1042 | 661.5 | - |
|  | Graph-EFM | 699.1 | 317.5 | 1.13 | 817.1 | 373.6 | 1.1 |
|  | AR-6h | 602.3 | 287.8 | 0.75 | 811.8 | 391.9 | 0.88 |
|  | AR-24h | 544.2 | 242.7 | 0.84 | 750.6 | 335.2 | 0.94 |
|  | CI-6h | 707.8 | 321.2 | 0.59 | 885.7 | 406.6 | 0.6 |
|  | ARCI-24/6h | 560.9 | 256.7 | 0.86 | 765.6 | 355.2 | 0.93 |
| t850 | Deterministic | 3.48 | 2.36 | - | 4.54 | 3.17 | - |
|  | Graph-EFM | 3.12 | 1.56 | 1.11 | 3.51 | 1.77 | 1.12 |
|  | AR-6h | 2.72 | 1.34 | 0.82 | 3.39 | 1.69 | 0.92 |
|  | AR-24h | 2.55 | 1.24 | 0.89 | 3.25 | 1.6 | 0.96 |
|  | CI-6h | 3.06 | 1.5 | 0.74 | 3.68 | 1.85 | 0.71 |
|  | ARCI-24/6h | 2.6 | 1.27 | 0.9 | 3.29 | 1.63 | 0.95 |
| t2m | Deterministic | 2.71 | 1.72 | - | 3.56 | 2.28 | - |
|  | Graph-EFM | 2.51 | 1.1 | 1.09 | 2.88 | 1.32 | 1.14 |
|  | AR-6h | 2.13 | 0.96 | 0.83 | 2.62 | 1.18 | 0.89 |
|  | AR-24h | 1.98 | 0.87 | 0.9 | 2.48 | 1.09 | 0.95 |
|  | CI-6h | 2.29 | 1.0 | 0.79 | 2.75 | 1.21 | 0.75 |
|  | ARCI-24/6h | 2.02 | 0.9 | 0.9 | 2.51 | 1.11 | 0.94 |
| u10 | Deterministic | 4.37 | 2.93 | - | 5.14 | 3.53 | - |
|  | Graph-EFM | 3.81 | 1.93 | 0.97 | 4.08 | 2.07 | 0.97 |
|  | AR-6h | 3.47 | 1.71 | 0.86 | 3.95 | 1.97 | 0.94 |
|  | AR-24h | 3.32 | 1.62 | 0.92 | 3.85 | 1.9 | 0.97 |
|  | CI-6h | 3.87 | 1.91 | 0.77 | 4.27 | 2.14 | 0.78 |
|  | ARCI-24/6h | 3.35 | 1.64 | 0.93 | 3.87 | 1.92 | 0.97 |
| v10 | Deterministic | 4.48 | 3.0 | - | 5.23 | 3.58 | - |
|  | Graph-EFM | 3.88 | 1.96 | 0.94 | 4.11 | 2.08 | 0.93 |
|  | AR-6h | 3.55 | 1.76 | 0.86 | 4.02 | 2.0 | 0.96 |
|  | AR-24h | 3.42 | 1.68 | 0.93 | 3.96 | 1.96 | 0.99 |
|  | CI-6h | 4.0 | 1.99 | 0.77 | 4.39 | 2.2 | 0.8 |
|  | ARCI-24/6h | 3.45 | 1.69 | 0.94 | 3.97 | 1.97 | 0.99 |
| ws10 | Deterministic | 3.09 | 2.2 | - | 3.44 | 2.49 | - |
|  | Graph-EFM | 2.56 | 1.35 | 1.01 | 2.65 | 1.4 | 1.01 |
|  | AR-6h | 2.38 | 1.23 | 0.9 | 2.57 | 1.34 | 0.96 |
|  | AR-24h | 2.3 | 1.18 | 0.96 | 2.53 | 1.31 | 0.99 |
|  | CI-6h | 2.54 | 1.32 | 0.88 | 2.68 | 1.4 | 0.89 |
|  | ARCI-24/6h | 2.31 | 1.19 | 0.96 | 2.54 | 1.32 | 0.99 |

Figure: Figure 10: RMSE, SSR, and CRPS for a selection of models at 6h resolution.
Refer to caption: x7.png

Figure: Figure 11: RMSE, SSR and CRPS for a selection of models at 1h resolution.
Refer to caption: x10.png

Figure: Figure 12: Temporal difference for different values of $\epsilon$ in algorithm 2. Choosing $\rho=0$ fixes the noise, $\rho=\text{ln}(10)$ allows it to vary slightly and $\rho\to\infty$ gives completely uncorrelated noise. The black line refers to the temporal difference of the data.
Refer to caption: x13.png

## Appendix D Example forecasts

Figures 13–18 show example forecasts for the remaining variables.

Figure: Figure 13: Example forecasts for geopotential at 500 hPa (z500) at lead time 10 days. The forecasts are generated using ARCI-24/6h except for Deterministic which is sampled using the deterministic model. The bottom row shows 4 ensemble members, randomly chosen out of the 50.
Refer to caption: extracted/5907711/images/forecast-z500.png

Figure: Figure 14: Example forecasts for temperature at 850hPa (t850) at lead time 10 days. The forecasts are generated using ARCI-24/6h except for Deterministic which is sampled using the deterministic model. The bottom row shows 4 ensemble members, randomly chosen out of the 50.
Refer to caption: extracted/5907711/images/forecast-t850.png

Figure: Figure 15: Example forecasts for ground temperature (t2m) at lead time 10 days. The forecasts are generated using ARCI-24/6h except for Deterministic which is sampled using the deterministic model. The bottom row shows 4 ensemble members, randomly chosen out of the 50.
Refer to caption: extracted/5907711/images/forecast-t2m.png

Figure: Figure 16: Example forecasts for u-component of wind at 10m (u10) at lead time 10 days. The forecasts are generated using ARCI-24/6h except for Deterministic which is sampled using the deterministic model. The bottom row shows 4 ensemble members, randomly chosen out of the 50.
Refer to caption: extracted/5907711/images/forecast-u10.png

Figure: Figure 17: Example forecasts for v-component of wind at 10m (v10) at lead time 10 days. The forecasts are generated using ARCI-24/6h except for Deterministic which is sampled using the deterministic model. The bottom row shows 4 ensemble members, randomly chosen out of the 50.
Refer to caption: extracted/5907711/images/forecast-v10.png

Figure: Figure 18: Example forecasts for wind speed at 10m (ws10) at lead time 10 days. The forecasts are generated using ARCI-24/6h except for Deterministic which is sampled using the deterministic model. The bottom row shows 4 ensemble members, randomly chosen out of the 50.
Refer to caption: extracted/5907711/images/forecast-ws10.png

## Peer Reviews

confidence: 2
contribution: 2
initial_rating: 5
presentation: 3
questions: In the weakness section.
soundness: 2
strengths: The proposed method is presented clearly.
summary: This paper uses the diffusion model for continuous ensemble weather forecasting without requiring autoregressive steps to produce future values. As introduced in Section 4.2, the main idea is to replace the fixed Z with a stochastic process Z(T) with continuous sample trajectories.
weaknesses: W1: there have been exsiting works for nonautoregressive sequence forecasting by diffusion models. However, realted investigation and discussions are missing, e.g., 

[1] NeurIPS'21 CSDI Conditional Score-based Diffusion Models for Probabilistic Time Series Imputation

[2] ICML'23 Non-autoregressive Conditional Diffusion Models for Time Series Prediction 

W2: moreover, some questions and details are not provided.
- it would be better to explain what ensemble forecasting is formally
- it is still not very convincing why a stochastic process is better than random Gaussian noises for nonautoregressive generation. Because we can simply represent the forecast outputs as a vector of a length N. Then we can use a Gaussian vector to yield forecasts. 
- is there a significant difference between weather forecasting and general time series forecasting?
- what is the input size (the number of variables) of the model in the experiments? It seems there are many variables in the weather forecasting task. Did you perform univariate forecasting tasks for each variable?   
- how to fuse $n_{ens}$ trajectories to get the final predictions?

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 3
questions: Could you clarify the rationale behind using fixed noise across time steps? Additionally, has any sensitivity analysis or ablation study been conducted to validate this design choice? If so, please provide these findings.

Could you further explain how interpolation using an ODE solver contributes to the model’s performance, and how it compares to alternative interpolation methods?

Only two of the five forecasting variables are included in Table 2. Could you explain the rationale for this selective presentation and consider including results for all five variables to enable a comprehensive evaluation?

Could you provide formal proofs or derivations for key claims, particularly those regarding trajectory continuity and temporal stability?
soundness: 2
strengths: The model developed in this paper presents an approach that eliminates traditional autoregressive steps during sampling from diffusion models, allowing it to achieve high temporal resolution and efficiency in forecasting. This is an exciting approach for improving computational efficiency, a feature that has been challenging to achieve in weather prediction. 

The paper's application is exciting, as diffusion models are typically used in areas like image generation rather than meteorological forecasting.
summary: The "Continuous Ensemble Weather Forecasting with Diffusion Models" paper introduces an innovative framework for probabilistic weather forecasting using diffusion models. It focuses on enhancing computational efficiency by removing the autoregressive nature of sampling from diffusion models.
weaknesses: The paper lacks references to closely related methods, such as Zheng et al. (2023) on fast sampling of diffusion models via operator learning, which shares similarities with the proposed model’s sampling approach. It would be beneficial to clarify how the proposed sampling method diverges from these approaches.

Although not directly used for weather data, recent models like Dyfusion, Alternator, and Mamba (with s-mamba used specifically for weather forecasting) have shown promise in similar domains, such as sea surface temperature forecasting. Including these models as baselines or discussing their relevance could offer a more comprehensive comparison, illustrating where this model stands relative to other contemporary approaches.

The paper mentions “regularity conditions” (e.g., Line 234) without specifying them rigorously. These conditions are fundamental to claims about trajectory continuity and temporal stability of forecasts. A formal mathematical description of these conditions, alongside proof or empirical validation, would be necessary to substantiate these claims.

The claims regarding fixed noise across time steps and ODE-based interpolation are under-supported. No mathematical evidence or ablation study is presented to show the impact of fixing noise on forecasting performance. Including such a study could clarify the advantages of this design choice.

The proposed Algorithm 2 seems tailored to a specific diffusion model setup. It would be valuable to discuss whether this algorithm generalizes to other diffusion models or is fundamentally limited to the framework outlined in the paper. A broader discussion on this topic could enhance the paper’s applicability to related models.

The results lack completeness, as only two of the five variables used for forecasting are presented in Table 2. Including comparisons for all five variables would provide a clearer picture of the model's strengths and weaknesses across different atmospheric variables.

The absence of error bars due to computational limitations weakens the robustness of the reported results. Given the probabilistic nature of the model, showcasing variability over multiple runs would offer greater confidence in the findings. Including error bars or confidence intervals for at least a subset of key results would be helpful, even if computational constraints prevent a complete assessment.

Numerous claims about model performance, such as trajectory continuity, lack rigorous mathematical proof and appear to rely on empirical observations. Addressing this would lend greater theoretical robustness to the proposed model.

There is no mathematical derivation for the training objective and loss function. A step-by-step derivation of these, with explanations of any modifications relative to standard diffusion model objectives and the rationale for these modifications in the weather forecasting context, would clarify the theoretical foundation of the approach.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: **Major:**
- It is a good idea to use the probability flow ODE from the work proposed from https://arxiv.org/pdf/2206.00364 in diffusion models. But what is the assumption behind the customized ODE in Eq. (2)? Why is there a linear noise level $sigma$ in Eq. (2)? What if $\sigma^2$? More explanation would help. 
- From the paper title and section 4 name, "CONTINUOUS" seems to be used to describe $ensemble forecasting$. However, what the authors really want to convey is the temporal continuity, i.e., any arbitrary fine temporal resolution. You may need to think over the rephrasing. Otherwise, it is confusing for readers (at least for me).
- In Figure 1, do you need to include the 1-hour prediction as the input for the 6-hour prediction? I saw you have explanations in sections 3 and 4. But Figure 1 deserves a more precise description within its title since it is the main framework of your method.
- It would be better to provide more comparisons with previous SOTA ML models. Even though I know not all models did not release source codes to the public, but some did, e.g., FourCastNet and ClimaX. 
- The authors mentioned the limitation of diffusion models: computationally expensive. I suggest adding more analysis of computational efficiency.

**Minor:**
- You are using the letter 'k' for different representations (see Eq. 1 and Alg. 1).
- It is customary to denote a time point with the superscript or subscript. It looks like a function if using $X[k+1]$ or $X(\tau)$.
- In the line under Eq. (2), why is $z(1)$ for $X[k+1]$? Isn't it the number of diffusion steps or the noise levels?
- What is $v_i^{k}$ in Algorithm 2?
soundness: 2
strengths: - It is interesting to use probability flow ODE of the noise in diffusion models. Leveraging a lead-time-dependent deterministic ODE-solver to generate a smooth and continuous trajectory is meaningful compared to autoregressive methods.
- Authors provide intensive experiments on their method.
summary: This work proposes an ODE-based method to address the limitations of previous a single forecasting step and rolled out autoregressively methods for weather prediction. It can generate long-range weather predictions without iteration and forecast arbitrary lead times.
weaknesses: - Need more descriptions on how probability flow ODE works.
- Need more comparisons with other SOTAs. 
- No available codes to help readers better understand the paper.

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 2
questions: 1.Were the limitations of previous methods that they were computationally expensive? How does the method proposed in the paper address the limitations of previous approaches?

2.What specific contributions does the paper make that differentiate it from other weather forecasting models?
soundness: 2
strengths: 1. Weather forecasting is an important problem.

2.The method introduced in the paper seems to exhibit a certain level of novelty.

3.Some of the figures in the paper look good.
summary: This paper introduces Continuous Ensemble Forecasting for global weather forecasting.
weaknesses: 1.The paper's presentation seems somewhat disorganized. The introduction does not clearly articulate the limitations of current weather forecasting methods or how the proposed method addresses these limitations. Figure 1, being the first figure, lacks sufficient information; ideally, a single figure should clearly highlight the differences between the proposed method and previous methods, as well as explain how the proposed method overcomes the limitations of prior approaches. In the Methods section, there is a tendency to mix descriptions of existing methods with the new method, which can obscure the exact nature of the contribution made by the proposed method. Figure 2 should also provide a clearer depiction of the methodology.

2.Based on the unclear description provided in the paper, the contributions appear insufficient, and the contributions listed in the introduction seem to be something that any weather forecasting model could achieve.

3.The baseline models compared in the paper are somewhat limited and do not appear to be the latest SOTA models. Moreover, according to Table 1, the ARCI-24/6h model proposed in the paper (if I understand correctly) is not the optimal model under most settings.

Meta-review: The paper introduces a new approach for ensemble weather forecasting using diffusion modeling. The method is not autoregressive and can generate ensemble member trajectories from initial conditions without relying on iterative steps, even though it can be combined with autoregressive methods for long-range forecasts. Uncertainty quantification is provided via repeated sampling. The approach is theoretically sound.

Decision: Accept (Poster)