Probabilistic Weather Forecasting with Hierarchical Graph Neural Networks
Authors: Joel Oskarsson
Tomas Landelius
Marc Peter Deisenroth
Fredrik Lindsten
Venue: NeurIPS 2024 spotlight (2024)
Keywords: weather forecasting
graph neural network
probabilistic
ensemble forecasting
latent variable model
earth system modeling
Primary area: machine_learning_for_physical_sciences
TL;DR: We introduce a probabilistic graph neural network model for weather forecasting, capturing uncertainty by generating ensemble forecasts.

Abstract: In recent years, machine learning has established itself as a powerful tool for high-resolution weather forecasting. While most current machine learning models focus on deterministic forecasts, accurately capturing the uncertainty in the chaotic weather system calls for probabilistic modeling. We propose a probabilistic weather forecasting model called Graph-EFM, combining a flexible latent-variable formulation with the successful graph-based forecasting framework. The use of a hierarchical graph construction allows for efficient sampling of spatially coherent forecasts. Requiring only a single forward pass per time step, Graph-EFM allows for fast generation of arbitrarily large ensembles. We experiment with the model on both global and limited area forecasting. Ensemble forecasts from Graph-EFM achieve equivalent or lower errors than comparable deterministic models, with the added benefit of accurately capturing forecast uncertainty.

Probabilistic Weather Forecasting with Hierarchical Graph Neural Networks
Joel Oskarsson
Tomas Landelius Swedish Meteorological and Hydrological Institute tomas. ins@smhi
Linköping University joel oswarsson@lu...ss
Mare Peter Deisenroth University College London
Fredrik Lindsten Linköping University fredrik. innstee@lli..s
deissenooth@ucll..... ac. uk
Abstract
In recent years, machine learning has established itself as a powerful tool for high-resolution weather forecasting. While most current machine learning mod- els focus on deterministic forecasts, accurately capturing the uncertainty in the chaotic weather system calls probabilistic modeling. We propose probabilistic weather forecasting model called Graph-EFM, combining flexible antent-variable formulation with the successful graph-basen forecasting framework. The use of hierarchical graph construction allows for efficient sampling of spatially coherent forecasts. Requiring only a single forward pass per time step, Graph-EFT EFM allows for fast generation of arbitrarily large ensembles. We experiment with the model on both global and limited area forecasting. Ensemble forecasts from Graph-EFM achieve equivalent lower errors than comparable deterministic models, with the added benefit accurately capturing forecast uncertainty.
1 Introduction
Forecasting the dynamics of Earth's atmosphere is scientific problem of utmost importance. Society is dependent on fast and informative weather forecasts for planning in areas such as transportation and agriculture and for balancing the energy system [3. Especially important the use of forecasts to Learning-based Weather Prediction (MLWP) have enabled models that produce accurate forecasts in a fraction the time of traditional physics-based systems (37 23 So far these developments have largely been focused on deterministic modeling. However, forecasting only one likely weather scenario ignores the many uncertainties in predicting future weather.
issue warnings for extreme weather events [I]. Recent advances in
Weather chaotic system, resulting in high forecast uncertainty [2] This uncertainty comes from both imperfect representations of initial states and inaccurate descriptions of the function mapping from one time step to the next (26]. Accurately modeling this uncertainty significantly increases value of weather forecasts. Such uncertainty can be communicated to end-users to improve decision making or be used in downstream products, for example to compute a distribution over solar power generation. Capturing the full forecast uncertainty requires trajectory, but a collection of possible future weather states. Due to the complexity and dimensionality of the weather system the feasible way achieve this is by generating samples from modeled distribution. Such ensemble forecasting is today performed using physics- methods, where a number of ensemble members simulated samples from this distribution. The computational cost of this however massive, often limiting the spatial resolution or size the ensemble [3]
predict not just single likely state
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

MLWP is a promising approach for addressing this limitation and enabling large ensemble forecasts. However, for the ensemble add value the machine learning model needs to accurately represent the distribution. Initial attempts at MLWP ensemble forecasting either rely on ad-hoc initial state perturbations
10 37 or have not been scaled to spatial resolutions of interest [9]. Also diffusion models [8] have been applied to the problem, but sampling forecasts from these is computationally
expensive and can be prohibitively slow [39]. We propose a Graph-based Ensemble Forecasting
Model (Graph-EFM), enabling efficient sampling of ensemble members with only one forward-pass per time step. The method builds on graph-based MLWP 22 231, which is a flexible framework that can be adapted to different geometries and state grid representations [24] By combining latent-variable formulation with a hierarchical Graph Neural Network (GNN) the distribution is modeled in swwerddmeensioal space and sampled forecasts are spatially coherent.
MLWP models are typically trained for and evaluated on global weather forecasting Another common forecasting setup in practice is the use of Limited Area Models (LAMs) to produce high-resolution regional forecasts [ Such L.MSs are for example used by local weather services in order to provide forecasts tailored to the geographical properties and societal needs of the region [38 7 32). These high-resolution models are also invaluable to various industrial sectors, including energy forecasters, who rely precise weather predictions to manage supply and demand. This motivates research into also constructing MLWP LAMs, which brings new challenges related to the high resolution and boundary conditions of the limited area. In this work we experiment not just with global forecasting, but consider also how probabilistic LAMs can be trained to produce forecasts for the Nordic region.
40
23
Our main contributions are: 1) We develop hierarchical GNN framework for both deterministic and probabilistic MLWP. The hierarchical construction encourages spatially coherent fields in forecasts. 2) We use this framework to define the probabilistic weather forecasting model Graph-EFM capable of efficient sampling of arbitrarily large ensemble forecasts. 3) We develop training method targeting both forecast quality and ensemble calibration. 4) We experiment with both global forecasting on 1.5" resolution and novel limited-area modeling task at 10 km resolution.
2 Related Work
Deterministic MLWP Multiple machine learning methods have been successfully applied to
large-scale weather forecasting. These include graph-based models 120 23 24], transformers
25 34 and neural operators (37 []. While large neural network models learn weather dynamics purely from data, there are also parallel developments in building hybrid physics-MLWP
4 10
33
models 22 50
Ensembles from perturbations Most existing methods for MLWP ensemble forecasting follow closely the physics- based methods, where initial states and model parameters perturbed to create ensemble diversity. A number of MLWP works create ensembles by ad-hoc perturbing initial states
with random noise [10 37
4 16, 6 More informed perturbations have been re- re-used from physics based ensembles 39 and created based on model-informed singular vectors [43] Others try to perturb the forecast model itself, rolling ensemble members using different neural network
parameters [51] 43 Such multi-model approaches require training, or at least fine-tuning, a pre-
defined number of MLWP models. Graubner et al. [16] use the SWAG method [30] to allow for constructing multi-model ensembles of arbitrary size.
Generative modeling Probabilistic machine learning approaches aim to directly learn generative models producing ensemble members. Similar to our approach, the Swin model [19 uses a latent variable formulation, but combined with Swin Transformer architecture (29). SwinVRNN is developed for global forecasting 5" resolution and scales poorly higher spatial resolutions. Also building on the graph- based framework, Price et al. 339 train a diffusion model [18] 46) to sample each time step. Their Gencast model produces ensemble forecasts of 0.25 global data with 12 h time steps. Diffusion models produce realistic-lookin samples, but typically require solving an ordinary differential equation involving multiple passes through the neural network to sample each time step. For GenCast, this results in a sampling time of 8 minutes for single 15 day forecast on TPUv5 device (39]. Other works use diffusion models to increase the size of physics-based ensembles 27 or stochastically downscale deterministic forecasts 9 31]

---

Hierarchical GNNs
Motivated by capturing multiple spatial scales, hierarchical GNNs have been
used for modeling general partial differential equations [8]. The overall hierarchical framework shares much of its structure with the popular U-Net architecture 41 for computer vision tasks, but extended to a general graph setting.
Background
3.1 Problem Definition
The weather forecasting problem can summarized as mapping from a set of initial states X-1:0 (X-1,xx) to the sequence of future states X = (X .x."....."""" X²). A table of notation is provided in appendix A Each weather state X E RNdx here contains dz variables modeled at N different locations. Geospatial data often represented as regular grids, in which case these locations correspond to the grid cells. The d, variables can include both atmospheric variables, modeled at multiple vertical levels, and surface variables. As is common in MLWP we assume the initial states to consist of two time steps, which allows for capturing first-order state dynamics. To produce forecast, a set of forcing inputs FL. are also available. These contain known quantities, such as the time of day. There are also static features associated with the grid cells, such as the orography, which we here consider part of the forcing.
Many variables impact the chaotic weather system, all of which are not fully captured in initial states represented on finite grids. This induces forecast uncertainty, which we view as a distribution p(x1:1 1:7 X x FI.7). In deterministic forecasting we seek model that minimizes the Mean modeling only the mean of the distribution. In probabilistic forecasting we instead aim to model the full distribution. Note that we here specifically model the conditional distribution p(X1:I' X -1:0, 1::0, p1:²), rather than p(X I:T FI.T). Hence we do not marginalize over uncertainty in initial states.
Squared Error (MSE) to the future weather states [23 B3] 34 This is equivalent
3.2 Graph-based Weather Forecasting
Graph-based MLWP models use an autoregressive mapping X = f(xt-2tt1111 consisting of sequence of GNNs (20 23 24). Starting from the initial states, this mapping can be iteratively applied to roll out a full forecast X 1:2", Central to the graph-based framework is the idea of mapping from the original N grid locations to the grid locations as a set Ve of grid nodes. By choosing [Var < [Vel = N it becomes efficient to perform the majority of computations forecasting setting, for example to respect the spherical geometry in global forecasting [20]. The mapping / realizes single-step prediction by passing Xt-22t-1 and F through a series of GNN layers. In sequence, these layers: 1) map grid inputs to representations on the mesh graph; 2) perform number of processing steps on the mesh; 3) map back to the grid to produce the prediction for X Steps and use bipartite graphs GGZM = (Ve U Vvr εG2M) and GM2G = (Va U V with edges connecting the grid and mesh nodes. The GNN layers in each step compute updates for node representations H R [V]xd and edge representations E R (|xd, in the graphs. For simplicity all representation vectors have dimensionality d2.
mesh graph GM (VM, EM). In the graph-context
refer
the mesh. Such mesh graph can also be tailored to the
Interaction Networks The specific GNN layers used in previous works are Interaction Networks 23) The layers in these networks pass messages from a set of sender nodes along directed graph edges to set of receiver nodes. Based on these messages the edge and receiver node representations are then updated. For a graph g = (V, () let Ca-g E R.d. be the row of E corresponding to the edge (a, 3) E. Let H be the matrix with rows containing sender node representations and / the corresponding matrix for receiver nodes. Interaction Networks then implement the representation update H E GNN(G, H E, HR)
H H3) Ea- MLP (Eo- MLP(E
(Ia)
ea+g t- e ea+g + 2.49
H + - + MLP HH. EnENe(S)
(1b)
where Ne(33) = fa: (a, B) E E} are the incoming neighbors of node 3. Parameters in Multi-Layer Perceptrons (MLPs) are shared across nodes and edges in the graph, but not between GNN layers.
3

---

Sample Z
x F*)
Residual
connection
9
GNN Layer
X1-21-1 F
X
91
P [x' zz, (Xx'("="--_))
/raph--MM
Figure 1: Overview of our Graph-EFM model, with example data and graphs for a Limited Area Model. The corresponding overview for the global setting is given in fig. 6 in appendix C
Global mesh graphs Keisler 20 proposed to construct a mesh graph for global MLWP as an icosahedral grid covering the globe. This approach was extended in the GraphCast model (33 by introducing a multi-scale mesh graph with edges of varying length. Such multi-scale edges are capable of propagating information and capturing statistical dependencies both locally and over long distances the graph. The multi-scale mesh graph is created by sequentially splitting the faces of an icosahedron into sequence of graphs .., 91 with node sets satisfying Vz C... V1 by construction. The original icosahedron 9L has the longest edges &L, stretching far across the globe, whereas the final graph 9i has short edges ε, only connecting nodes locally. The final multi-scale mesh graph is constructed as GMS &r U... E taking the nodes from the final graph but connecting these using edges of all different lengths 23
4 Weather Forecasting with Hierarchical Graph Neural Networks
Two great challenges in weather forecasting is to accurately capture processes unfolding over different spatial scales and modeling the uncertainty in the chaotic system [52]. To tackle these challenges, we propose to construct hierarchical mesh graph, working with different length scales at each level in the hierarchy. We use a sequence .1...., 9L of graphs as the different levels in the hicrarchy, additionally adding connections between the nodes of adjacent levels. This construction is also highly suitable as a basis for building probabilistic forecasting models, as discussed below. Figure I shows an overview of the hierarchical mesh used in our model. See figs. 12 and 14 in the appendix for illustrations of how this differs from the multi-scale graph.
There are multiple benefits to such hierarchical mesh construction for MLWP. By keeping the graphs at different levels separate, we can define GNN layers with independent parametrizations at each level. This adds flexibility by allowing the model to learn different representation updates for edges of different spatial scales. A hierarchical mesh graph also offers a natural, spatially-aware dimensionality reduction, as the state in the grid is encoded into few nodes at the top level. Such a representation can capture the general structure of each weather state, with finer details added as this is propagated down through the hierarchy. We leverage this property construct probabilistic model by imposing a distribution over these lower-dimensional representations at the top level. This allows for efficiently drawing spatially coherent samples from the distribution of future weather states.
4.1 Hierarchical Graph
Our hicrarchical mesh graph consists of graph levels 9I,..., Gr with G. = EV, Only level of the hierarchy is connected to the grid, so re-define GGIM (Va U EG2M) and GMDG = (Va U V,, EM2G). The number of nodes [i] decreases with the level l. The smallest set of nodes are found at the top level L.
To pass information between the levels of the hierarchy we introduce additional graphs connecting the different levels. Let Gi, (V) U V1+1. E1,+1) be graph containing directed edges from mesh level to level /+1 1. We make use of graph sequence G1,2, Gr.-1.F. I, to propagate information up through the hierarchy and similarly sequence Gr. ... G2, in the downward direction. The exact layout of nodes and edges at and in-between levels are design choices that should be tailored to the
specific forecasting setting, Examples for global and limited-area forecasting given in section 5

---

4.2 Graph-FM: Deterministic Forecasting
The hierarchical graph allows for defining GNN layers both on and in-between the different levels. By sequentially updating node and edge representations at different levels in the hierarchy, information can be propagated up from the grid to the different levels. As these levels have edges of different lengths, the processing at each level happens on different spatial scales, Note that this differs from the multi-scale graph approach, where information processing over all different spatial scales happen in the same GNN layer (23] As a step towards probabilistic model, we define alternative deterministic Graph- based Forecasting Model Graph-FM operating on the hierarchical graph.
In Graph-FM one processing step on the mesh graph is defined as complete sweep through the hierarchy. GNNs are applied sequentially to the inter-level and intra-level graphs in the order updating edge and node representations at the different levels. Process- ing steps going up the hierarchy are alternated with similar steps going down from level L to 1. The single step mapping f consists of multiple such sweeps up and down (see appendix C.2
G1 21,2, , 1.1.1,1,,11,
4.3 Graph-EFM: Probabilistic Forecasting
76-1
Z
To capture the uncertainty in the chaotic weather system we next aim to construct a probabilistic model from the ground up to capture the full
distribution p(X - -x00 FI:T). We start by assuming the weather
system to satisfy a second-order Markov assumption, decomposing
X
x pI:T) [X'x F²).
Figure 2:
Graphical
Factoring the distribution over time steps allows us to work with forecasts of varying length. Specifying the model for single-step prediction avoids having to learn separate parameters for different lead times. Next, we seek a flexible, but computa- tionally efficient parametrization for the distribution p( xx'''" 2:t-1 F². This can be achieved by introducing latent random variable Z', and letting
model for eq. 3
p(x*|x(-2111 F²) = / ,( 21-2111 ( 201-1 F')dZ'.
Here the stochasticity in Z' should capture the uncertainty over X at each time step. The corre- sponding graphical model is shown in fig. We impose a spatial structure over the latent variable by letting Z' be |vr.| d: matrix-valued, with each row a dg_ -i-dennnsoo vector associated with one node in the top level Gr. of the mesh graph.
The single-step model consists of two components, a latent map p(Z' X FF) and predictor Z'
( 2 F²). The latent map is parametrized using GNNs, mapping the conditioning variables to parameters of Gaussian distribution. We consider the predictor to be concentrated around its mean, and realize p(X' pX'1 xt-2:t11, F²) as a deterministic mapping of a similar form as Graph-FM. By sampling Z' and passing this through the predictor we can draw a sample of X from eq. 3 This sample can then be conditioned on at the next time step, continuing this sampling process to roll out a forecast following eq. (2). This forecast constitutes one ensemble member, and the process can be repeated to sample ensemble of arbitrary size. We call our Graph-based Ensemble Forecasting Model Graph Full details about the model are given in appendix C
Latent map
We let the latent map be an isotropic Gaussian
p x Fr) = IL NZ (Zzz)/z(( ,F)),,))
(4)
with the mean as a function of the conditioning variables. The variance is fixed, imposing a fixed scale for the learned latent space. The mean function μzz consists of a sequence of GNNs. These take the inputs at the grid, propagate representations up through the hierarchical mesh graph, and finally predicts the mean of Z at each node a at level L. In appendix L.2 we verify empirically the importance of using the latent map over static distribution for Z'
'The deterministic Graph-FN FM model was first proposed in a preliminary version of this work [35 but there only for the LAM setting under the name Hi- HAM.
5

---

Predictor The predictor is deterministic mapping
X = g(Z'. x 2:1-1 _F*) x
- x2t22t-1 .F").
(5)
With the small time steps used in MLWP, X does not change dramatically in a single step. We thus follow the common practice of including a skip connection to the previous state 23 5 The predictor takes both inputs X²-2:--1, F² at the grid and Z' at the top of the mesh graph. To incorporate both we design g similar to Graph-FM, performing sweeps up and down through the mesh hierarchy. At the top of the hierarchy Z' is added to node representations H through the residual connections in the GNN layers. A sampled value of Z' then affects the prediction X through the downward sweep. While multiple such sweeps are possible, we found one be sufficient in practice.
19
Spatial dependencies We want each sample of X' to contain spatially coherent atmospheric fields. One approach would be to impose spatial dependencies in the joint distribution over Z". However, learning and sampling from such a distribution typically comes with computational challenges [19]. Instead, we impose spatial dependencies by integrating the latent variable formulation with the hierarchical graph. We argue that as the independent components of Z' are propagated down through the mesh graph, gradually increasing the spatial resolution, spatial dependencies are introduced by the model in the GNN layers. The hierarchical graph is key to this property, as the stochasticity in Z necessarily spread out over the forecast region, rather than only affecting the output locally.
4.4 Training Objective
Deterministic forecasting models can be straightforwardly trained by minimizing a weighted MSE [23 or Negative Log-Likell Log-Likelihood (NLL) loss [8] for rolled out forecasts. To train Graph- EFM we instead leverage the fact that the single-step model has structure similar to (conditional) Variational AutoEncoder (VAE) 45), allowing us to use variational objective. We introduce at each time step, approximating the true poste- rior pl (Z) X X²,F*) over Z This variational distribution is parametrized in a similar way as the latent map, with GNN layers mapping to a Gaussian over Z', Note however that & also depends on Xt, since approximates the posterior. Using q. we can then define
variational approximation q( (Z X (x1-2:11,, x F²))
Lvaa Lx x F') DKL (1 X X F²) ( 2:0-1 FF))
(6)
- Enevo Σ log. N (xas (z x F²) an 225]
which is equal to the (negative) Evidence Lower Bound (ELBO) when the weighting is Aki. 1. While the predictor is a deterministic mapping, we introduce Gaussian likelihood in eq. 6 to
get a well-defined learning problem. This setup corresponds to the common practice in VA of assuming Gaussian observation noise, but not adding this to samples from the model 42 The standard deviation Jo.j can either be a second output from the predictor or manually chosen (see appendix D for details). As with deterministic models 23 20 g 34, we found it crucial to fine-tune on rolled out forecasts of multiple time steps. This improves stability and performance for longer lead times. In the final fine-tuning
include also a Continuous Ranked Probability Score (CRPS) loss term LCRPS [15] 22). The full objective function is then L = L + ACRPS CCRPS, with ACRPS a weighting hyperparameter. Including this CRPS loss improves the calibration of ensemble forecasts.
4.5 Improved GNN Layers: Propagation Networks
In Graph-EFM there is a large amount of information that needs be propagated between the grid and z However, the Interaction Network GNNs are biased towards keeping old representations of receiver nodes, rather than updating this with new information from incoming edges. Note in eq. (T) that if the MLPs are initialized to give outputs close to 0, there will be no change to ca and H
In practice the model has a hard time learning to propagate useful information up from the grid to Zt. Even when trained purely as an auto-encoder (Aąı. = 0), Z' easily ends up being ignored. To remedy this we propose an alternative GNN formulation that we call Propagation Network, defined by
Ca-B f- H + MLP (Ca-.g H H
c 3
Ca-18 E
(7a)
1
i +- Ne(BB EaeNe(n)ea--p
R t i n - MLP H ) H
(7b)

---

For MLPs initialized with outputs close to 0, Propagation Networks reduce to averaging the values of neighboring nodes. This encourages the propagation of information from H to H by construction. Propagation Networks were found to perform better also in the deterministic model (see comparison in appendix I..T. so we employ these in both Graph-FM and Graph-EFM.
5 Experiments
To evaluate our models we conduct experiments on both global and limited area forecasting. The mod- els are implemented in Py PyTorc and trained on 8 A100 80 GB GPUs in data- tata--arallll configuration. Training takes 700-1400 total GPU-hours for the global models, and around half of that for the limited area models. The computational demands prevent from re-training multiple models for statistical analysis. Once trained, sampling from Graph-EFM is highly efficient. Using batched sampling on single GPU, 80 ensemble members are produced in s00 (2.5s s per member) for global forecasting.
Metries We measure the skill of deterministic models by Root Mean Squared Error (RMSE). For probabilistic models we compute the RMSE for the ensemble mean. Good skill terms of RMSE is however not enough for ensemble forecasts, where we want to capture the full distribution. For these we also assess the ensemble calibration by computing the Frread-Skill-aee Ratio (SpSkR). Calibrated uncertainty corresponds to SpSkR = [13]. We additionally use CRPS to measure how well the marginal distributions of the model matches the data. For deterministic models the CRPS reduces to Mean Absolute Error (MAE). Complete definitions of all metrics are given in appendix E
Models Achieving a fair comparisons of the actual machine learning methodology in MLWP is challenging due to models using different spatial resolution, variables and initial states. We here train an illustrative set of models on the same data and with comparable training setups. Our full Graph-EFM model is compared to: 1) Graph-FFM (ms), a version of Graph-EFM using multi- scale mesh graph instead of the hierarchical one. 2) Graph-FM, our deterministic model using the hierarchical graph. 3) GraphCast reimplementation of GraphCast [23], adapted and trained on our datasets. 4) GraphCast*
multi-model ensemble created by applying Stochastic Weight Averaging Gaussian (SWAG) (30 to GraphCast Inspired by Graubner et al. [16], this represents a simple way to augment deterministic model to perform ensemble forecasting. Further details about the baseline models are given in appendix C.5 For ensemble models we sample 80 members for the glohal experiments and 100 members for limited area forecasting. In appendix [..3 L..] we investigate the impact of ensemble size on the evaluation. We find that improvements in metric values quickly saturate when increasing the ensemble size. This shows that sampling even more members would have negligible impact on the results of our experiments.
5.1 Global Forecasting with ERAS
Data and graphs We experiment on global weather forecasting up to 10 days with h time steps. The dataset used for training and evaluation is 1.5" version of the global ERA5 reanalysis 3 [17], provided through the WeatherBench 2 benchmark [40]. The models forecast de 83 different variables in total, including both surface-level variables and atmospheric variables at 13 different pressure levels. We use the years 1959-2017 7 for training, 2018-201 2002 for validation and 2020 as a test set. Forecasts are always started from initial conditions taken directly from ERA5, both during training and evaluation. For global forecasting we use the graph generation process from GraphCast [23]. The multi-scale graph GMs is created by refining the icosahedron 4 times. The hierarchical graph contains 4 levels of such icosahedral grids. More details on the global experiments are given in appendix H
Results As the models forecast many different variables we present only a selection of results the main paper. Metric values for geopotential (z500) and 2 m temperature (2t) are listed in table and results for mean sea level pressure (msl) plotted in fig. 3 Line plots for all metrics and a large number of variables are given in appendix J.1. In the appendix we also show comparisons to additional models from the literature, trained on different data, as well as the physics-based IFS-ENS model 12]. The
"Our code available at (glohal forecasting) and andhttt::: / //t///gtttuub htam/mllam/neurraa http://eeeerraalaaaerreerrr.abba wttrre//roomooeeeee_eeee 1ammem_mmmmmm____ (L.MM.
http://ttw..cccmmmmmm//aaaaaaaaa http://wwwwaaaeeeaeaaaaaaaaaaaaa wwa.mmeeeeoaoaallllllllllleaaaaa global
Provided by the Copernicus Climate Change Service under the ECMWF Copernicus License.
7

---

Table Selection results for global forecasting, including geopotential at 500 hPa (z500) and m temperature (2t). For RMSE and CRPS lower values are better, and SpSkR should be close to 1 for calibrated ensemble. The best metric values are marked with bold and second best underlined.
Lead time days
Lead time 10 days
Variable
Model
RMSE
CRPS
SpSkR
RMSE
CRPS
SpSkR
z500
GraphCast*
387
236
808
498
Graph-FM
363
223
825
510
GSWP
437
269
0.07
960
590
0.12
Graph-EFM (ms)
472
211
333
0.83
Graph-EFM
399
169
1.18
695
299
1.15
GraphCast*
1.65
1.00
2.82
1.69
Graph-FM
1.57
0.94
2.82
1.66
GraphCast*+SWAG
2.03
1.20
0.06
3.58
2.04
0.13
Graph-EFM (ms)
1.76
0.77
0.75
2.55
1.09
0.82
Graph-EFM
1.64
0.71
2.32
1.00
0.99
Gmphhass
Graph
Graph (ms)
Graph-EPM
GraphCast*
Calibrated
750
1.0
100
500
200
0.5
250
0.0
8
10
5
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) RMSE
(b) CRPS
(c) SpSkR
Figure 3: Results for global forecasting of mean sea level pressure (msl) at all lead times.
ensemble mean from Graph-EFM often shows improvements in RMSE over the deterministic models, especially for longer lead times. Across the ensemble models, Graph-EFM achieves lower CRPS values, better capturing the distribution of the weather data. Without any perturbations to initial states Graph-EFM reaches a SpSkR close to We note that GraphCast* ensemble forecasts, these poorly calibrated and in general do not lead improved forecast errors. Figure 4 shows an example forecast from Graph-EFi for specific humidity (q700) at 10 days lead time. Examples for other variables are given in appendix [.2
does not produce useful
Ground Truth
Ens. Mean
Ens. Std.-Dev.
2003
2.002
Ens. Members
Figure 4: Example Graph-EFe ensemble forecast for specific humidity at 700 (q700), for lead time 10 days. The bottom row shows 3 ensemble members, randomly chosen out of the 80.

---

Table 2: Selection of results for LAM forecasting, including geopotential 500 hPa (2500) and integrated column of water vapor (wvint).
Lead time h
Lead time 57 h
Variable
Model
RMSE
CRPS
SpSkR
RMSE
CRPS
SpSkR
z500
GraphCast*
153
108
201
138
Graph-FM
230
162
354
238
GSSWG
219
136
0.08
376
206
0.10
Graph-EFM (ms)
400
261
0.22
711
470
0,23
Graph-EFM
172
91
0.84
219
115
0.75
wvint
GraphCast*
1.51
1.01
2,06
1.32
Graph-FM
1.64
1.08
2.48
1.58
SSWAG
1.78
1.17
2.34
1.50
0.05
Graph-EFM (ms)
2.39
1.43
0.16
3.51
2.12
0.13
Graph-EFM
1.61
0.79
0.57
2.08
1.00
0.53
Extreme weather case study An important use case for ensemble forecasting is modeling extreme weather events. While higher resolutions than 1.5% are generally desirable for accurately capturing such extremes, we conduct one case study on using Graph-EFM for forecasting hurricane Laura. The full case study with visualized forecasts available appendix F For this example we show that there exists ensemble members accurately predicting the landfall location of the hurricane at 7 days lead time, while the deterministic models still show no sign of the hurricane in the region. Closer to the land- fall event the ensemble forecast from Graph-EFM indicates uncertainties associated with the landfall location and wind intensity. This demonstrates the added value of probabilistic forecasting model.
5.2 Limited Area Modeling with MEPS Data
In LAMs weather forecasts are produced for a bounded region of the globe. LAM forecasting allows for higher resolution modeling and regionally tailored model configurations [1I), properties that can be inherited by MLWP models by training on LAM data. To model weather over a limited domain, boundary conditions need be taken into account. In physics-based LAMs these are typically given by a global forecast (38 7 32). We adapt a similar approach for MLWP LAMs, by taking boundary conditions as additional forcing along the boundary of the forecast area. The problem of LAM forecasting is thus about simulating physics not just based on the initial state, but also consistent with these boundary inputs. In the models we introduce N, additional grid nodes along the area boundary, for the boundary forcing B' E IR Nod, Boundary forcing B is always fed together with X² the model. Grid nodes on the boundary and within the area are treated identically by the GNN layers. We perform this adaptation to all models in our experiment.
Data and graphs We experiment with a dataset containing 6069 forecasts from the MetCoOp Ensemble Prediction System (MEPS) LAM. Training on forecasts, the goal is here to learn a fast surrogate model for MEPS. We use forecasts started during April 2021 Jun 2022 for training and validation, and forecasts from July 2022 March 2023 as a test set. The data is laid out in a 238 x 268 grid with spatial resolution 10 km, covering the Nordic region. This dataset contains in total d, 17 weather variables, some repeated on multiple vertical levels. Forecasts are rolled out with 3 time steps up to lead time 57 h. In this experiment we take also the boundary forcing directly from the MEPS dataset. We define the boundary as the outermost 10 grid positions. Using the same dataset for the area and boundary allows us to investigate the modeling choices in controlled experimental setup. In an operational scenario the boundary forcing would instead come from a re-gridded global forecast. the LAM setting we define our graphs as regular quadrilateral meshes covering the MEPS forecasting area, but with far fewer nodes than the original grid. The graph hierarchy G1, ......, Gr is created by constructing such meshes at different resolutions. By placing each node in Gi at the center 3 x 3 nodes in Gi-1, we merge 4 such graph levels to create GMS. In the hierarchical graph we instead introduce edges from each node in Gi to the x 3 nodes in the level below. More details about the MEPS data and experiment can be found in appendix !!
Results A selection of metrics are shown in table 2 and full results given in appendix K At these shorter lead times there is no clear benefit of probabilistic modeling in terms of RMSE. Still,

---

Graph-EFM
Graph-EFM (ms) Ens. Member
Graph-EFM
Ground Truth
Ens. Member
10
Ens. Std.-Dev.
33000
-15
130000
Figure Example forecasts for net solar longwave radiation (nlwrs) at lead time 57
exemplified by the coandadd-evvatiio plotted in fig. 5 probabilistic modeling provides useful information about the forecast uncertainty. Comparing the ensemble members in fig. 5 highhights the improved spatial coherency of the hierarchical graph in Graph-FFM. In contrast, the Graph-EFM (ms) forecast looks patchy and lacks physically intuitive features. There are also clear visual artifacts, that can be traced the multi-scale graph structure. We discuss this more in-depth appendix G In the LAM setting all models are under-dispersed, with SpSkR < 1. One explanation for this is that the boundary forcing constrains the space of plausible forecasts, hindering the ensemble spread.
Discussion
In this paper we have explored MLWP ensemble weather forecasting using graph-based latent variable models. Our Graph-EFM model is capable of efficiently producing accurate ensemble forecasts. This paves the way for large-scale MLWP ensemble forecasting both in operational use and research settings. In appendix B we further discuss the societal impact of this research. With this work hope to emphasize that MLWP models are not just deterministic mappings, but parametrize distributions of weather states. It follows that ensemble forecasting should not be achieved by perturbing models, but by directly modeling the distribution of interest.
Limitations The training process comes with some complications in terms of choosing a training schedule and byperparameters ÅKL and ACRPS. While the CRPS fine-tuning is an important training step. we have found that choosing the Graph-EFM (ms) model (see appendix G). While Graph-EFN EFM produces diverse and physically plausible ensemble members, the forecasts still suffer from some of the blurriness common to deterministic models [23] 40). We here trade off some of the visual fidelity achieved for example by diffusion models
high ACRPS can introduce visual artifacts, especially for
for more efficient sampling of ensemble members.
Future work Interesting avenues for future work include learning probabilistic weather models based on other types of autoencoders [48 49), or by directly optimizing scoring rules [36 22 Another approach for achieving efficient ensemble forecasting is to explore techniques for speeding up diffusion model sampling [47].
Acknowledgments and Disclosure of Funding
This research is financially supported by the Swedish Research Council via the project Handling Uncertainty Machine Learning Systems (contract number: 2020-04122), the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation, the Excellence Center Linköping Lund in Information Technology (ELLITT), and the project OWGRE, funded by partners of the ERA-Net Smart Energy Systems and Mission Innovation through the Joint Call 2020. As such, this project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 883973. Our computa- tions were enabled by the Berzelius resource the National Supercomputer Centre, provided by the Knut and Alice Wallenberg Foundation.
10

---

References
[1] M. Astitha and Nikolopoulos. Definition of extreme weather events (subchapter 1.1). In Extreme Weather Forecasting, pages 1-7. Elsevier, 2023.
[2] P. Battaglia, R. Pascanu, M. Lai, D. Jimenez Rezende, and K. Kavukcuoglu. Interaction networks for learning about objects, relations and physics. In Advances in Neural Information Processing Systems, volume 29, 2016.
[3] Bauer, A. Thorpe, and G. Brunet. The quiet revolution of numerical weather prediction. Nature, 2015.
[4] K. Bi, L. Xie, H. Zhang, X. Chen, X. Gu, and Q. Tian. Accurate medium-range range global weather forecasting with 3d neural networks. Nature, 2023.
[5] B. Bonev, T. Kurth, C. Hundt, J. Pathak, M. Baust, K. Kashinath, and A. Anandkumar. Spherical Fourier neural operators: Learning stable dynamies on the sphere. In Proceedings of the 40th International Conference on Machine Learning, 2023.
[6] C. Bülte, N. Horat, J. Quinting, and S. Lerch. Uncertainty quantification for data-driven weather models. arXiv preprint arXiv:2403. 1ate24403333458, 2024.
[7] M. Bush, I. Boutle, J. Edwards, A. Finnenkoetter, C. Franklin, K. Hanley, A. Jayakumar, H. Lewis, A. Lock, M. Mittermaier, S. Mohandas, R. North, A. Porson, B. Roux, S. Webster, and M. Weeks. The second met office unified model-JULES regional atmosphere and land configuration, RAL2. Geoscientific Model Development, 2023.
[8] K. Chen, T. Han, Gong, L. Bai, Ling, J.-J. Luo, X. Chen. Ma, T. Zhang, R. Su, Y. Ci. B. Li, X. Yang, and W. Ouyang. FengWu: Pushing the skillful global medium-range weather forecast beyond 10 days lead. arXiv preprint arXiv:2304. 0300:23004.0248,, 2023.
[9] I.. Chen, F. Du, Y. Hu, Z. Wang, and F. Wang. Swinrdm: integrate swinrnn with diffusion model towards high-resolution and high-quality weather forecasting. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 37, 2023.
[10] L. Chen, X. Zhong, F. Zhang, Y. Cheng, Y. Xu, Y. Qi, and H. Li. Fuxi: a cascade machine learning forecasting system for 15-day global weather forecast. npj Climate and Atmospheric Science, 2023.
[11] J. Coiffier. Fundamentals of numerical weather prediction. Cambridge University Press, 2011.
[12] ECMWF. Ifs documentation cy46rl part Ensemble prediction system, 2019, URL https:/ h/pp:://wwww http::/wwwwwwww.w...mm//m/mm inttiinne//iiieiieiiioooooo 19309
[13] V. Fortin, M. Abaza, F. Anetil, and R. Turcotte. Why should ensemble spread match the RMSE of the ensemble mean? Journal of Hydrometeorology, 2014.
[14] M. Fortunato, T. Pfaff, Wimsserger, A. Pritzel, and Battaglia, Multiscale meshgraphnets. In ICML 2022 2nd for Science Workshop, 2022.
[15) T. Gneiting and A. E. Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 2007.
[16] A. Graubner, K. Kamyar Azizzadenesheli, J. Pathak, M. Mardani, M. Pritchard, K. Kashinath, and A. Anandkumar. Calibration of large neural weather models. In NeurlPS 2022 Workshop on Tackling Climate Change with Machine Learning, 2022.
[17] H. Hersbach, B. Bell, P. Berrisford, S. Hirahara, A. Horányi, J. Muñoz-Sabater, J. Nicolas, C. Peubey, Radu, D. Schepers, Simmons, C. Soci, Abdalla, X. Abellan, G. Balsamo, P. Bechtold, G. Biavati, J. Bidlot, M. Bonavita, G. De Chiara, P. Dahlgren. D. Dee, M. Dia- mantakis, R. Dragani, Flemming, R. Forbes, M. Fuentes, A. Geer, L. Haimberger, S. Healy, R. J. Hogan, E. H6lm, M. Janisková, S. Keeley, P. Laloyaux, P. Lopez, C. Lupu, G. Radnoti, de Rosnay, Rozum, Vamborg, Villaume, and J.-N. Thépaut. The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological Society, 2020.
11

---

[18] J. Ho. A. Jain, and P. Abbeel. Denoising diffusion probabilistic models. In Advances in Neural Information Processing Systems, volume 33, 2020.
[19] Y. Hu, L. Chen, Z. Wang, and H. Li. SwinVRNN: A data-driven ensemble forecasting model via learned distribution perturbation. Journal of Advances in Modeling Earth Systems, 2023.
[20] R. Keisler. Forecasting global weather with graph neural networks. arXiv preprint arXiv:2202.07575, 2022.
[21] D. P. Kingma and M. Welling, Auto-encoding variational Bayes. In International Conference on Learning Representations, 2014.
[22] D. Kochkov, J. Yuval, I. Langmore, P. Norgaard, J. Smith, G. Mooers, M. Klöwer, J. Lottes, Rasp, P. Düben, Hatfield, P. Battaglia, A. Baachez-Gonzalez, M. Willson, M. P. Brenner, and S. Hoyer. Neural general circulation models for weather and climate. Nature, 2024.
[23] R. Lam, A. Sanzhez-Gonzalez, M. Willson, Wirnsberger, M. Fortunato, F, Alet, S. Ravuri, T. Ewalds, Z. Eaton-Rosen, W. Hu, A. Merose, S. Hoyer, G. Holland, O. Vinyals, J. Stott, A. Pritzel, S. Mohamed, and P. Battaglia. Learning skillful medium-range global weather forecasting, Science, 2023.
[24] S. Lang, M. Alexe, M. Chantry, J. Dramsch, F. Pinault, B. Raoult, M. C. A. Clare, C. Lessig, M. Maier-Gerber, L. Magnusson, Z. B. Bouallègue, A. P. Nemesio, P. D. Dueben, A. Brown, F. Pappenberger, and F. Rabier. AIFS-ECMWF's data-driven forecasting system. arXiv preprint 26X1::2406.01665, 2024.
[25] C. Lessigg I. Luise, B. Gong, M. Langguth, S. Stadler, and M. Schultz. AtmoRep: A stochastic model of atmosphere dynamics using large scale representation learning. arXiv preprint arXiv:2308.. 13280, 2023.
[26] M. Leutbecher and T. Palmer. Ensemble forecasting. Journal of Computational Physics, 2008
[27] L. Li, R. Carver, I. Lopez-Gomez, F. Sha, and J. Anderson. Generative emulation of weather forecast ensembles with diffusion models. Science Advances, 2024.
[28] M. Lino, C. Cantwell, A. A. Bharath, and S. Fotiadis. Simulating continuum mechanics with multi-scale graph neural networks. arXiv preprint arXiv:2 2021.
[29] Z. Liu, Y. Lin, Y. Cao, H. Hu, Y. Wei, Z. Zhang, S. Lin, and B. Guo. Swin transformer: Hierarchieal vision transformer using shifted windows. In Proceedings of the IEEE/CVF International Conference on Computer Vision, 2021.
[30] W. J. Maddox, P. Izmailov, T. Garipov, D. P. Vetrov, and A. G. Wilson. A simple baseline for bayesian uncertainty in deep learning. In Advances in Neural Information Processing Systems, volume 32, 2019.
[31] M. Mardani, N. Brenowitz, Y. Cohen, J. Pathak, C.-Y. Chen, C.-C. Liu, A. Vahdat, K. Kashinath, J. Kautz, and M. Pritchard. Residual diffusion modeling for km-scale atmospheric downscaling. arXiv preprint arXiv:2309. 15214, 2023.
[32] M. Müller, M. Homleid, K.-I. Ivarsson, M. A. Ø. Koltzow, M. Lindskog, K. H. Midtbo, U. An- drae, T. Aspelien, L. Berggren, D. Bjørge, P. Dahlgren, J. Kristiansen, R. Randriamanpianina M. Ridal, and O. Vignes. AROME- MetMoo:etCoOp: A nordic convective-scale operational weather prediction model. Weather and Forecusting, 2017.
[33] Nguyen, J. Brandstetter, A. Kapoor, J. K. Gupta, and A. Grover. ClimaX: A foundation model for weather and climate. In Proceedings of the 40th International Conference on Machine Learning, 2023.
[34] T. Nguyen, R. Shah, H. Bansal, T. Arcomano, S. Madireddy, R. Maulik, V. Kotamarthi, I. Foster, and A. Grover. Scaling transformer neural networks for skillful and reliable medium-range weather forecasting, arXiv preprint arXiv:2312.03887 2023.
12

---

[35] J. Oskarsson, T. Landelius, and F. Lindsten. Graph-based neural weather prediction for limited area modeling, In NeurlPS 2023 Workshop on Tackling Climate Change with Machine Learning. 2023.
[36] L. Pacchiardi, R. A. Adewoyin, P. Dueben, and R. Dutta, Probabilistic forecasting with generative networks via scoring rule minimization. Journal of Machine Learning Research, 2024.
[37] J. Pathak, S. Subramanian, P. Harrington, S. Raja, A. Chattopadhyay, M. Mardani, T. Kurth, D. Hall, Z. Li, K. Azizzadenesheli, Hassanzadeh, K. Kashinath, and A. Anandkumar. Four- castnet: A global data-driven high-resolution weather model using adaptive fourier neural operators. arXiv preprint arXiv:2202. 1214, 2022.
[38] T. V. Pham. Steger, B. Rockel, K. Keuler, I. Kirehner, M. Mertens, D. Rieger, G. Zängl, and B. Früh. ICON in climate limited-area mode (ICON release version 2.6.1): a new regional climate model. Geoscientific Model Development, 2021.
[39] I. Price, A. Sanche7-Gonzalez, F. Alet, T. Ewalds, A. EI-Kadi, J. Stott, S. Mohamed, P. Battaglia, R. Lam, and M. Willson. Gencast: Diffusion- baasd ensemble forecasting for medium-range weather. arXiv preprint arXiv:2312.157%6 2023.
[40] S. Rasp, S. Hoyer, A. Merose, I. Langmore, P. Battaglia, T. Russel, A. Sanzaez-Gonzalez, V. Yang. R. Carver, S. Agrawal, M. Chantry, Z. B. Bouallegue, P. Dueben, C. Bromberg, J. Sisk, L.. Barrington, A. Bell, and F. Sha. WeatherBench 2: A benchmark for the next generation of data-driven global weather models. arXiv preprint arXiv:2308. 15560, 2023.
[41] O. Ronneberger, P. Fischer, and T. Brox. U-net: Convolutional networks for biomedical image segmentation. In Medical Image Computing and Computer-Assisted Intervention. Springer International Publishing, 2015.
[42] O. Rybkin, K. Danillidis, and S. Levine. Simple and effective vae training with calibrated decoders. In Proceedings of the 38th International Conference Machine Learning, 2021.
[43] S. Scher and G. Messori. Ensemble methods for neural network-based weather forecasts. Journal of Advances in Modeling Earth Systems, 2021.
[44] Y. Seity, P. Brousseau, S. Malardel, G. Hello, P. Bénard, F. Bouttier, C. Lac, and V. Masson. The AROME-france convective-scale operational model. Monthly Weather Review, 2011 1.
[45] K. Sohn, H. Lee, and X. Yan. Learning structured output representation using deep conditional generative models. In Advances in Neural Information Processing Systems, volume 28, 2015.
[46] Y. Song and S. Ermon. Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems, 32, 2019.
[47] Y. Song, P. Dhariwal, M. Chen, and I. Sutskever. Consistency models. In Proceedings of the 40th International Conference on Machine Learning, 2023.
[48] I. Tolstikhin, O. Bousquet, S. Gelly, and B. Schoelkopf. Wasserstein auto-encoders. In International Conference on Learning Representations, 2018.
[49] A. van den Oord, O. Vinyals, and k. kavukcuoglu. Neural discrete representation learning. In Advances in Neural Information Processing Systems, volume 30, 2017.
[50] Y. Verma, M. Heinonen, and V. Garg. ClimODE: Climate forecasting with physics-informed
neural ODEs. In International Conference on Learning Representations, 2024.
[51] A. Weyn, D. R. Durran, R. Caruana, and N. Cresswell-Clay: Sub-seasonal forecasting with a large ensemble of deep-learning weather prediction models. Journal of Advances in Modeling Earth Systems, 2021.
[52] J.-I. Yano, M. Z. Ziemiański, M. Cullen, P. Termonia, Onvlee, L. Bengtsson, A. Carrassi, R. Davy, A. Deluca, S. L. Gray, V. Homar, M. Köhler, S. Krichak, S. Michaelides, V. T. J. Phillips, P. M. M. Soares, and A. A. Wyszogrodrki. Scientific challenges of convective-scale numerical weather prediction. Bulletin of the American Meteorological Society, 2018.
13

---

Appendix Table of Contents
Table of Notation
14
B Societal Impact
14
Model Details
16
D Training Objectives
21
Metrie Definitions
23
Extreme Weather Case Study: Hurricane Laura
24
Spatial Coherency of Forecasts
24
H Experiment Details: Global Forecasting with ERAS
28
Experiment Details: Limited Area Modeling with MEPS Data
32
J Additional Results: Global Forecasting with ERAS
36
Additional Results: Limited Area Modeling with MEPS Data
50
L Additional experiments
65
A Table of Notation
Notation used throughout the paper is listed in table 3
B Societal Impact
Extreme weather Due climate change the prevalence and severity of extreme weather events is expected to increase substantially, endangering both property and human life [56]. These events are also getting harder to predict and the cost of damages per event has increased nearly 77% over the past five decades 622. In order to detect extreme events, there is a need for ensemble forecasting to better capture the full distribution of possible weather states. While ensemble forecasting historically has been limited by computational costs 33, efficient MLWP ensemble models have the potential to vastly improve our ability to model extreme weather. In appendix we present case study showing how our Graph-EFM could be used for forecasting Hurricane Laura. More extensive evaluation of the abilities of MLWP models capture extreme events is complex, but important issue [44 55 We view these capabilities as one of the main motivations for further developments of MLWP ensemble forecasting models.
Forecast failures There will always be cases where weather forecasting systems fail, and produce inaccurate predictions for future weather. Depending on the weather event, such forecast errors can have disastrous consequences. In these cases the lack of interpretability of black-box MLWP systems can be a problem, making it hard to understand why the model was wrong. Traditional physics- based systems can also be tough to interpret due to their complexity. but at their core are physical equations understood by researchers. With the rapid progress of MLWP, it is likely that we will soon see landscape where the physical models take the back seat to plethora of skillful MLWP forecasting systems. In such scenario the question of how to investigate forecast failures becomes pressing. seems desirable to be able to fall back on physical models for understanding impactful events poorly
14

---

Table 3: Table of notation
Weather state at time step
X
Predicted weather state at time step t
((m)
Ground truth weather state at time step in sample from dataset
Predicted weather state time step f in sample from dataset
(.m)).00)
Prediction of ensemble member It time step in sample from dataset Forcing inputs time step
pt
Zt
Latent random variable for time step Boundary forcing input at time step 1, for limited area modeling
B²
H
Matrix with node representation vectors as rows Matrix with edge representation vectors as rows Representation vector for the edge going from node a to node
a
Autoregressive mapping predicting the next state in deterministic models Predictor part of Graph-EFM, predicting the next state conditioned sample of Z Predictor output before skip connection
g
Pz
Mean function of latent map, mapping to the mean of Z
Length of one forecast, in discrete time steps
N
Number of grid nodes (grid cells)
No
Number of grid nodes in the boundary region, for limited area modeling Number of levels in hierarchical graph
L
S
Number of samples (forecasts) in dataset
K
Ensemble size, number of members in ensemble forecast
dg
Weather state dimensionality (total number of variables) in each grid node
d,
Latent dimensionality, dimensionality of each node or edge representation vector
Gm
Mesh graph
Va
Nodes of mesh graph
Ex
Edges between the nodes of the mesh graph
Va
Grid nodes, corresponding to grid cells in the data Graph connecting grid nodes to the mesh
GGzM
EOzm
Edges of GuzM- going from grid nodes to mesh nodes Graph connecting the mesh to grid nodes
GMzG
EmG
Edges of GMzG. going from mesh nodes to grid nodes Multi-scale graph
GMs
Gi
Intra-level graph at level / in hierarchical mesh graph Nodes at level of hierarchical mesh graph
VI
&
Edges of G. intra-level edges at level I in hierarchical mesh graph
G
Inter-level graph with edges from level to level + in hierarchical mesh graph Edges of Gr,41,, from level / level ( in hicrarchical mesh graph Incoming neighbors of node 3, {a (a. B) E E}
E
Ne(3)
Area weighting for node a, proportional to the area of corresponding grid cell Weighting for KL-term in variational training objective
AKL.
AcRPS
Weighting for CRPS term in fine-tuning objective of Graph-EFM Standard deviation used in NLL or variational training objective
15

---

Algorithm Single-step prediction / for graph-based MLWP
1: HG MT.P(X1-22-1, Ft)
Embedd grid inputs dy-dimensional vectors >> Map grid representation to mesh
2: HM .. GNNNGGGM. HC ECEM, H
3: for all mesh processing steps do
4: 11 MMM, + GNN(GS H EMM HA
Update mesh representations
5: end for
6: H t- H + MLL(H
7: HO. .. GNNNGMZZ. HM, , HG)
b Map mesh representation back to grid
8: if outputs then
9: [Y.]] - MLP(HG)
10: return X Softplus(U))
> Return prediction of X and a for loss
11: else
12:
return X MML((H))
D Return prediction of X
13:
end if
forecast by MLWP. However, to allow this we can not do away with all infrastructure and expertise related to physical modeling. These are important considerations as weather forecasting moves into an era operational MLWP.
Renewable energy The production of renewable energy, such as solar and wind power, can be highly volatile [60]. This creates challenges for including these sources in the larger energy system. Accurate weather forecasts, translated into forecasts of energy production, fill an important role as enablers of these energy sources by making their output predictable. Detailed probability estimates from ensemble forecasting can additionally allow for improved cost-loss decision making these systems.
The energy footprint of weather forecasting Traditional weather forecasting systems utilize massive computing clusters [II], resulting in substantial energy footprint of the forecasting process. In comparison, MLWP models are highly energy efficient, even when taking into account their initial training process. The total energy required to train the large global MLWP model FourCastNet is comparable to running 10 day forecast with 50 ensemble members using traditional forecasting system 377 Producing forecast using the same model uses four orders of magnitude less energy than physics- thssics-based model. A similar reduction in energy footprint is to be expected for our models, if applied at the same scale. For ensemble forecasting the total energy saving even greater, as is multiplied by the number of ensemble members. However, one has to beware of rebound effects, where efficiency improvements result in more extensive use of resources. If ensemble forecasting becomes 1000 times more energy efficient and we run 1000 times as many ensemble members there no energy saved in the end.
C Model Details
We here give more details about the different models discussed the main paper.
C.1 Deterministic graph-based models [23 20
Deterministic graph-based MLWP models represent the single-step prediction function X = f(x 2:t-1 Ft) as sequence of MLPs and GNN layers. Algorithm describes the full prediction process. Note that for some of the GNN layers the edge representations are not updated, as these edges are not used in later steps of the process. We denote this by 11, - GNN in algorithm MLPs always act on the last dimension of its concatenated inputs. Mesh node and edge representations are initialized by embedding related static features using additional MLPs. These static features contain information about node and edge positions (see appendices H. and I.1). D). The number of processing steps to run on mesh is a hyperparameter. Each such step contains one GNN layer, and these GNNs do not share parameters. The mapping f can optionally output also standard deviations a, which we use when training with NL. loss (see appendix D.Tp. We restrict to be positive by applying the Softplus function.
16

---

Algorithm 2 Single-step prediction f Graph-FM
1: HG MT.P(X Ft)
Embedd grid inputs dy-dimensional vectors
Map from the grid p...and and up through the whole mesh hierarchy 4: 11',"-11 t GNN(G H² EE H1)
2: H¹ . .GNNNGGaaa H EG2M )
3: for t= :....... do
5: end for
6: for number of mesh processing steps 2 do
> Each iteration contains two steps, down and up
7:
H² E t GNN(GL, H E H
8:
for - (L 1)....11 do
> Downward pass of sweep
9:
1',E"+ii + GNNNG H , E HI) 11' E' + GNN H1, E 11) end for
10:
12:
H', E +- GNN H E H
for I ...... do 1 E'
Upward pass of sweep
14:
H' J",E"S ",,E-11 H1, E. < GNN(G, H E","'')
15:
end for
17: end for for
= (L ......,,dl. do H',. GNN(G++1. 1141 El+11, H²)
19:
Map down through the mesh hierarchy
20: end for
H H + ML)(HH))
22: H GNN(GM2G. H EMG H.)
D and finally to the grid
if outputs a then
24: [Y,U] < MLP
25: return X Y, Softpllus(UU)
D. Return prediction of X and a for loss
else
27: return X +ML(H))
Return prediction of X
28:
end if
C.2 Graph-FM
In the hierarchical graph there are node representations associated with every level and edge represen- tations associated with each subset of edges. We let H' be the representations associated with nodes Vi level / the mesh hierarchy. Similarly, let E contain representations of intra-level edges & El,1+1 representations of upwards edges Er,+1, and E representations of downward edges E As in the non-hierarchical case, all representations associated with the mesh graph are initialized by MLPs applied to static features. Independent MLPs are used for the different levels of the hierarchy.
Forecasting using Graph-FM follows similar structure as previous works, but replaces the mappings between the grid and mesh to encompass all levels and changes the processing steps into sweeps through the hierarchy. Recall that for the hierarchical graph we re-define GG2m (VGUV, EGTM) and GM2G (Ve U V1, EM2G) The full prediction mapping f for Graph-FM is described in algorithm 2 We define one processing step Graph-FM as one update all node and edge representations in the graph, keeping similar interpretation as processing steps in non-hieree non-hierarchical graphs. This means that one sweep down and up (lines 7 16 in algorithm ( is counted two processing steps. follows that only even numbers of processing steps are reasonable in Graph-FM. Note also here that GNNs only share parameters across nodes and edges the specific graph they operate on. There no parameter sharing across processing steps nor levels in the hierarchy. We include Propagation Networks in Graph-FM, but not all parts the model can expected to benefit from this change. For some steps of algorithm we want to keep the inductive bias of Interaction Networks to retain information. Specifically, we use Propagation Networks for the mappings directly between the grid and lowest mesh level (lines 2 and 22 and for the upward processing steps (line
17

---

Sample Z
g3
pl 2: pr)
Residual
connection
G2
GNN Layer
X 2:21111 Ft
1
p X PP)
Graph-FM
Figure 6: Overview of our Graph-EFM model, with example data and graphs for global forecasting. The corresponding overview for the LAM setting is given fig. [ Note that the predictor part of the model has the same structure as one sweep through the hierarchy in Graph-FM.
Algorithm 3 Sampling of ensemble member from p( 1:77 |X FI:T)
1: X --10 f- X- -1:0
2: for = 1...., , do
3:
~ --p Z x1-22-11 F()
D Sample Z' from latent map
4: X + x'+8-1 + t 81--1122-1111
>> Compute next state using predictor
5: end for
6: return X
C.3
Graph-EFM
An overview of the Graph-EFM model with global graphs is shown in fig. 6 At each tim step Graph-EFM takes input X Fs and produces a sample of X By feeding the sample from the model back at the next time step, an ensemble member can be rolled out. This process described in algorithm 3 where we let X an initial state for < 1 and then assign sampled predictions for the actual forecast time steps 1....,
When creating an ensemble of size K from deterministic initial conditions some care has to be taken as to how many members to sample. If we at each time would sample K values of X for each realization of Xt-2:t-1 we would end up with K" members at time T. Instead, we sample K realizations of X given the initial conditions and after that only one realization of X for each of the members. This equivalent to doing K independent runs of algorithm 3 Note that since we can never restricted in the ensemble size.
always draw new samples of Zt
Latent map
The latent map is defined as
p(z' x F*)
II aEVi
NZZz)uz(X' :..)..1)))))
(8)
where the mean function μzz is realized hierarchy. The exact process is described in algorithm 4 Since only one upward pass is performed we do update any of the edge representations. The final MLP at minpss from node representations at mesh level L to the mean of p Z' X x
sequence of GNN layers mapping up through the
2:t-1 F In the latent map we use exclusively Propagation Networks, encouraging the model to encode useful information from the grid into Z'.
Predictor
The predictor p( I 2 F²) is chosen as dirac measure, with all probability mass concentrated in one point. It thus takes the form of deterministic mapping
X g( Z X F*) x + 8 (Z' - (
This choice emphasizes that all randomness p ( F²) should come from the latent variable Z¹. An alternative would be parametrize the predictor as a diagonal Gaussian. We found this inferior in practice, as that would require also sampling from this Gaussian when rolling out
18

---

Algorithm 4 Latent map mean function
1: HG MT.P(X Ft)
Embedd grid inputs dj-dimensional vectors
2: H¹ . GNNNGGGM. H EG2M )²²)
H², -GNN(G,, H E" 111))
4: for | = ,.....LA do
> Propagate up the hierarchy
5: H' GNN(9-111 H E HH)
6: H'.- - GN H E')
7: end for
8: return ML
Return mean of p x 2 ,F)
Algorithm 5 Predictor function
1: Hº MLP Ft)
Embedd grid inputs d dy-dimensional vectors
2: H' <- CNN((G2M) H EG2M H
3: for -....., --1....-11
Upward pass of sweep
4:
H' E' - GNN(G), H', H')
5:
if I. then
6:
H < GNN(GL- H4-1 FL-1LL z
> Incorporate Z' through GNN
7: else
8:
H + GNN(G H Eut+1 H1+1)
9: end
10: end for
11: H4, +- GNN(9L, EL, H²)
12: for = (L-11 1)......, do
Downward pass of sweep
13:
H' + GNN(G H44 E++LI H))
14:
H'.-t GNN(G, H E H")
end for
16: HG H + MLP(H()
17: HC. HG,.+ +. GNN(GM2C, H., M2G HG)
18: if outputs then
19: [YU]= + MLPP (H)
20: return X Y, Softplus(U)
b Return prediction of X and a for loss
21: else
22:
return 1--1 + MLP(HO))
D Return prediction of X
23:
end if
forecasts. This sampling would simply entail adding independent Gaussian noise the forecasts, which only degrades their quality. We still re-introduce the Gaussian form in the variational objective (eq. (6). as some choice of likelihood is necessary for well-defined learning problem. The interpretation of this is that the forecast is rolled out as a latent process with all randomness coming from Z Observation noise is then assumed added this forecast independently each time step. We use the latent, noise-free free forecast as our prediction. Connecting back to the VAE analogue of our model, this setup corresponds to the common practice in the VAE literature using a Gaussian likelihood but treating the mean the generated sample [2).
The exact form of the predictor function is described in algorithm It has similar form Graph- FM, doing sweeps up and down the mesh hierarchy. We found it sufficient to perform one sweep up and down in the predictor, but it can easily be generalized multiple sweeps as in Graph-FM The predictor function optionally returns also the to be used in our variational objective. In the predictor we use Interaction Networks in the upward direction and Propagation Networks the downward direction. This is to guarantee that randomness and useful information encoded in Z' at the top level reaches all the way down to the grid. Note that on line 6 Z' combined with H through the Interaction Network layer mapping from level L happens, see eq. [. where Z' takes the role of the receiver node representations H R Note that edge and node representations updated on line in the upward sweep then re-appear in the downward pass. This constitutes residual connections between the upward and downward passes, as illustrated in fig.
to L. For the details of how this combination
19

---

Algorithm 6 Parameter mappings Mas Og of variational approximation
1: HG <- -T.P(X Xt, Ft)
b Embedd grid inputs (including Xt) to dg- vectors
2: Ht + -GNNGGGMM, H EGGMM _
3: H², +t -GNN(G,, H E .1!))
4: for w = .2....,La.
Propagate up the hierarchy
5: H', + GNN(G--11 H 10-,,EE--_ E H1')
6: H', - GNN E", )
7: end for
8: [
9: return Y, Softplus(U)
>> Return mean and standard deviation of g( ( x
x F
Variational approximation Defining the ELBO for the model requires a variational approximation ( X X F²) at each time step, approximating the true posterior p( Z X²-2:1-1, _X'', F²) over Zt. Due to the autoregressive structure of the model, it is sufficient to condition q states at time steps 2, t - and t. This conditioning removes all dependence between Z' and other variables (see fig. 2) simplifying the parametrization and training of q. We choose the variational
approximation to have a similar Gaussian form as the latent map
q(Z 2:t-1 x X F*)F)=
de II II Zas (Xt x²,F') (. a.j og (x (x1-23-11 x xx,, (2) aEVL
(10)
For added flexibility we also model the variance of q. Here both B4y and og are functions mapping from the grid inputs to |Vrl x dz matrices. These functions are implemented jointly as described in algorithm This mapping similar to the latent map, with the differences being that X is taken as an input and both mean and standard deviation returned. Similarly to the latent map, all GNNs used in the variational approximation are Propagation Networks.
C.4 MLP parametrization
Following Lam et al. (23), all MI.Ps have one hidden layer with Swish activation functions (59). With the exception of MLPs outputting predictions distribution parameters, apply LayerNorm [53 to all other MLP outputs.
While MLPs in different GNN layers never share parameters, we use some parameter sharing across the embedding MLPs in the different components of Graph-EFM. Specifically, the latent map and predictor share parameters for the MLP embedding the grid input (line I in algorithm and line in algorithm (5). The MLPs that embedd static graph features, for initializing graph representation vectors, are also shared across all components of the Graph-EFM model.
C.5 Baseline models
GraphCast* is reimplementation of the GraphCast model 23 in our codebase. This was done to allow for a fair comparison by using the same data and multi-scale graphs as other models in our experiments.
GraphCast
Graph-EFM (ms) The Graph-EFN, (ms) model is a version of Graph-EFM, but using a multi- scale mesh graph GMs [23], instead of the hierarchical one. We replace the sweeps up and down the hierarchy with multiple processing steps performed on the same mesh graph, similar to the deterministic described in algorithm In the predictor of this model, the latent Zt is integrated already in the mapping from the grid to the mesh as
MM, GNN(GGMM. HC (EE.,M, 22))
(11)
Note that as Z' here associated with the nodes GMS. it has shape |i] d, as opposed to [Vrl] x d₂ for Graph-EFM. As Vi has many more nodes, the total dimensionality of z' is higher for Graph-FFM (ms). This is a consequence of the multi-scale graph offering no further dimensionality reduction beyond mapping from the grid to the mesh.
20

---

GraphCas One simple way to create an MLWP ensemble is by re-training multiple deterministic models from random initializations. Training many large weather models from scratch is however infeasible in practice. It also limits the number of ensemble members to the number of models trained. Inspired by Graubner al. [66, we create multi-mode ensemble baseline using SWAG [S)] SWAG is technique for approximate Bayesian inference based on running Stochastic Gradient Descent (SGD) training with a constant high learning rate. At regular intervals throughout this optimization process the model parameters are saved. The full SWAG process includes estimating a high-dimensiooal Gaussian over these parameters and drawing samples to create the ensemble.
This gets around the limitation of needing to decide the maximum ensemble size when training the model. As we never require more ensemble members than saved during the SGD training we skip estimating this Gaussian and resampling. Instead, we directly use subset of the saved model parameters to create our ensemble members. Graubner et al. [16 combine SWAG with initial state perturbations, but as model the conditional distribution p(X X F1:2) we do not include any such perturbations to the initial states.
To create GraphCast*
ensemble we start from the trained GraphCast* model and run SGD training for T = 8 unrolled steps and with fixed learning rate (10 for the global experiment and 5 10-5 for MEPS). Even higher learning rates led to numerical issues or the model training diverging. Model parameters are then saved every 100 gradient descent steps. During evaluation these parameters are loaded one after another and used to produce the forecast for each ensemble member. The error and spread of the G ensemble is highly dependent on the correlation of consecutively saved model parameters. Given the poor spread of GraphCast SSWG tested also saving parameters with 1000 steps in-between, but this led to similar results.
D Training Objectives
We here give detailed definitions of the objectives used to train our models.
D.1 Deterministic models
Weighted MSE Deterministic forecasting models can be trained by minimizing weighted MSE
23 of rolled-out forecasts
Wnisjss 8 X(ad) )²
LWMSE
(12)
TN
tm aaVa
where we recall that N Vall and:
W is the normalized area weighting for grid cell a. These are computed as described by
Rasp al. (40
Xj is the inverse variance of time differences for variable j. These can be computed directly from iterating over the data.
wy is weight associated with the vertical level of variable j. For global experiments we use the same weights wy as Lam al. (23).
Negative Log-Like Log-Likelihood loss For the MEPS data however, choosing weights wy is challenging due to the many different variables and their irregular vertical levels. In such cases, an alternative approach is
NLL loss, where the model itself determines the trade-offs between different variables [ As the MEPS area has grid cells of approximately equal size we can also simply set W, 1. For the MEPS experiments we thus use
T
de Wa log (x (
ENLL = TN
Σ 1=1
&a. (2.))
(13)
NEVE
The standard deviation Øaj is here second output from the mapping f. When unrolling forecasts we do not sample from this Gaussian, but feed the mean X to the next time steps. Note that eq. (13 decomposes into weighted squared error terms and log-determinant terms preventing too large
21

---

outputs. The weighted MSE in eq. (12) can therefore also be viewed as a special case of the NLL loss, up an additive constant. This relationship corresponds eq. (33) with constant variances given by
La.t 2wgd
(14)
D.2 Probabilistic model
For the probabilistic model we can not directly apply the training objectives above, as they only seek to match the first moments of the model distribution to the data. To train the probabilistic model we instead leverage the fact that the single-step model has structure similar (conditional) VAE [2145] 45. We can therefore optimize variational objective, the ELBO, to match the distribution of the data.
Variational objective
For a single time step, our variational objective is
Xx', F²) [lp(2 X (F)) Ever (xt=2t-11 2:t-1 X F*) DKL X ",',
(15)
Bu F() Evere Ej=ra wa log N (xxx (0 2-----1 F*) F')a.) oa.)]
We include weight ÅKL. in front of the KL-regularizer to allow for some tuning between the two parts of the objective function [6]. practice we found this useful to make sure the model does not collapse pure deterministic predictions. Note that when ΛKL this is exactly equal to the (negative) ELBO. The NLL from eq. (13) here shows up again in the second term. Similarly to the deterministic models, øa.j can either be output by the model (the predictor in this case) chosen constant. For global experiments we choose Ja.j as in eq. (14), reducing the term to a weighted MSE loss. For the MEPS data we again let each Oa.j be output by the model. Equation (15) can in practice be computed efficiently by using single-sample Monte Carlo estimate for the expectation and the reparametrization trick (11] The KL-divergence is available in closed form as both distributions are Gaussian. As with deterministic models (23 20 g 34], it is crucial to fine-tune on rolled out forecasts of multiple time steps. This induces stability and improves performance for longer lead times. We thus define our fine-tuning objective as
LNar Et= Lvar ( (81-2:11 X F²)
(16)
where X is an initial state from the dataset for t 1 and otherwise sampled using the variatonal approximation for Z". In eq. (16) X is always from the training data.
CRPS fine-tuning In practice we found that fine-tuning the model with only C tends to result in underdispersed ensemble forecasts, underestimating the variance of the distribution. To remedy this,
we include an additional fine-tuning objective based the CRPS [15] 63 22]. The CRPS measures how well a univariate distribution matches the observed data and is defined as
CRPS(p(x).v) = E| vl
E(Ix - z'1]
(17)
2
where v, x' are independent copies of random variables with distribution p(z) and the observed data. The CRPS is minimized only when the predicted distribution matches that of the data [15]. making it suitable for calibrating our ensemble forecasting model. We define our CRPS fine-tuning loss an unbiased two-sample estimator, summed over all dimensions the predictive distribution
di
CCRPS = = 1=1
Σ
Wa
(xxa)
- a.j +| - x - - -x(,
X(,
(18)
LEVE j=1
with X and X coming from two independent ensemble members sampled from the model.
Note that the CRPS objective targets only the marginal distributions, while spatio-temporally coherent samples requires matching the joint distribution of the data. To avoid degenerate solutions we weight LCRPS by ACRES and combine it with the variational objective for the final fine-tuning loss = Lvar + >CRPS Computing L requires rolling out three forecasts, one for L (using samples from q) and two for LERPS (using samples from the latent map). This can have substantial
22

---

memory cost. It is however sufficient to include the CRPS term only in the final steps of the training process, making this less of an issue in practice.
Empirically we found that including the CRPS objective improves both ensemble calibration and forecast accuracy. The accuracy improvement can be attributed to the MAE terms in eq. (T8). Note that in Var forecasts are rolled out by sampling Z' from the variational approximation. while in CCRPS these samples are from the latent map. This matches how forecasting is carried As consequence the CRPS loss has an additional effect of bridging the gap between the variational training and the final forecasting.
time.
E
Metric Definitions
We here give full definitions for the metrics used to evaluate our models.
E.1
Deterministic Metrics
We define evaluation metrics based on set S forecasts, each a sequence ( Xl.m),..., (.mmm of predicted weather states. Metries computed per lead-time and variable. For deterministic forecasting we define the RMSE variable j at lead time as
S Σ
Σ Wa (XL(m)
(mm) X a.
RMSE,
(19)
SN
VoVo
In all metrics we include area weighting Wa to handle the fact that the grid cells for the global data has different area. These weights are computed described in Rasp al. [40]. For the MEPS data
these are all set to Wa = 1. Note that the square root is applied after sample averaging, following standard convention and the WeatherBench 2 benchmark [40]
E.2 Probabilistic Metrics
In the probabilistic setting we evaluate an ensemble forecast of K members. Let (t)( be the prediction of ensemble member k. For ensemble forecasts we compute the RMSE of the ensemble mean
gt..(m) 0.d
x-f((m) 2 or.
RMSE, =
Wa
(20)
SN
m=
VEVo
K
gl(mm)
(R))
(21)
K
k=1
To measure the calibration of uncertainty expressed by the ensemble the (bias-corrected) Spread-Skill-Ratio
K +
Spread
SpSkR,,, ..J
(22a)
K
RMSE,
S
Spread
SKN
Σ
Wa
(xt((m) a.j
(t)(m))(()) 2
a.
(22b)
the ensemble members represent realistic forecasts, and observe exchangeability with the ground truth, SpSkR,,, 13)
As probabilistic metric we also use the CRPS 15 here computed as finite sample estimate 63 over all the ensemble members
CRPS,,j
1 SN m=loeVe
Σ Σ
Wa
Σ
(m))
as
x-f(mm)
K
k=
(23)
K
K
2K(K
-1
0.j
km
23

---

Note that this is spatial average over marginal CRPS values and does not take any covariance structure of the data or model distributions into account. To avoid the quadratic sum in eq. (23) we in practice use the idea of Zamo and Naveau 63 to compute this with linear memory by sorting the ensemble members. For deterministic models the CRPS reduces to MAE, See this for example from eq. (23) by letting all ensemble members be the same predicted value.
F Extreme Weather Case Study: Hurricane Laura
One key motivations for probabilistic weather forecasting is to allow for better predictions of extreme events and improved estimates of uncertainties associated with these. To exemplify this we include here a case study looking at forecasts for Hurricane Laura, using the global models.
In August of 2020 Laura developed in the Atlantic Ocean and reached hurricane levels in the Gulf of Mexico, eventually making landfall in Louisiana and causing major damages [58]. We study here forecasts for 2020-08-2780 28-2712 UTC, which is 6 hours after the hurricane hit land. All forecasts are for this exact time point, but initialized a varying number of days before. We here run 50 member ensemble forecasts using Graph-EFM and deterministic forecasts using the GraphCast* and Graph-FM models. Figure shows 10 m wind speeds in ERA5 and the forecasts. For Graph-EFN -FM we plot both forecasts from randomly sampled ensemble members and a cherry- that was deemed to most closely match ERA5, Note that the 1.5 resolution that we work with makes determining similarity exact positions somewhat challenging. The resolution also puts some limitations on how well the most extreme wind speeds can be captured.
best member
We see in fig. that at 7 days lead time there is great uncertainty, and the deterministic models do not show the hurricane at all. In the ensemble forecast from Graph-EFM there exists however already members indicating the possibility of hurricane making landfall week ahead (see for example the cherry-picked "Best member"). While the ensemble includes many possible scenarios, a total of 7 members show the development of a hurricane in the area. Having information of such possible scenarios a long time ahead allows for planning and readying disaster response efforts that might be needed. Note that discovering these scenarios is only possible through ensemble forecast, as the deterministic models do not indicate such an event. At 5 days lead time all models are indicating the development of hurricane. While the deterministic models do a good job here, they are indicating landfall location slightly too far eastward. The ensemble members from Graph-EFM however show range of different positions for the hurricane, indicating the uncertainty in the landfall location. At 3 days ahead 42 out of 50 ensemble members show the hurricane making landfall, but still with some uncertainty about the exact location. At day lead time all models give an accurate forecast of the position of the hurricane. Apart from position, it is also interesting to consider how the models capture the intensity in terms of wind speeds. At 1 day ahead the deterministic models are somewhat underestimating the wind speed. The Graph-EFM ensemble shows a range of possible values, indicating the uncertainty in the exact wind intensity. Overall this study exemplifies how ensemble forecasts from a machine learning model can be used to discover possible extreme weather scenarios long lead times and uncover the uncertainties associated with them.
G Spatial Coherency of Forecasts
Apart from giving a low forecast error, it is also important for forecasts to look realistic. Specifically, we want forecasts to be physical, containing features that are possible under laws of atmospheric physics. This is necessary for forecasts to be interpretable by meteorologists and builds trust in the MLWP model. One key property for realistic forecasts is that of spatial coherency. Forecasts should take realistic values not just locally, in each grid cell, but show larger-scale features that are consistent over the forecasting area. This is important since the weather system contains spatial dependencies on multiple different scales. Note that in the probabilistic setting this connects closely to capturing the correct joint distribution spatially, rather than just the marginal distributions in each grid cell. As the metrics commonly used to evaluate MLWP models only consider the marginal distributions, it is important to also inspect sample forecasts to gauge the spatial coherency.
This modeling aspect is especially relevant for probabilistic models, since being able to visualize and interpret individual ensemble members is an important capability of ensemble forecasting. Other uncertainty quantification techniques can be applied to estimate forecast uncertainty (6), but these
24

---

Lead Time
2020-08-27T12
7 days
days
days
day
Best Member
18
ERAS
16
Wind Speed (m/s)
12
Member
Graph-EFM
Member 2
Member 3
Member 4
Graph-FM GraphCast*
Figure 7: Example forecasts of 10 wind speeds during Hurricane Laura, all for 2020-08-27T12 UTC. The first column shows the wind speeds in ERA5, with the red arrow indicating the location of the hurricane. Remaining columns show model forecasts for the specific time, initialized different lead times ahead. For Graph-EFM we plot both 4 randomly sampled ensemble members and one cherry-picked member that deemed to best match BRA5.
25

---

(a) nlars Graph-FFM
(b) nlwrs Graph-FFM (ms)
(c) Graph-EFM
Graph-FFM (ms)
(c) u65 Graph-EFM
(f) u65
(g) v65 Graph-EFM
(h) v65 Graph-RFM (ms)
Graph-FFM (ms)
Figure 8: Example forecasts from Graph-EFM and Graph-EFM (ms) trained on the MEPS data. All forecasts from single ensemble members, plotted for lead time 57 h.
generally do not generate samples of possible future weather states. Having samples to inspect can build trust and understanding of the model. If one member is predicting extreme event, it is possible to inspect this forecast and see what sequence of events in the atmosphere could lead to this outcome.
For the final Graph-EFa -FM model forecasts are generally spatially coherent and represent physically plausible scenarios. Examples of this is shown in appendices].2]andPp this with some the shortcomings of baseline models. also discuss some of the challenges that we have encountered with achieving spatially coherent forecasts and ways to tackle these.
In this appendix we contrast
G.1 Limited Area Modeling with MEPS Data
In the MEPS experiment the benefits of the hierarchical graph structure become the most clear when visually inspecting forecasts. Figure 8 shows comparison between a few forecasts from Graph-EFM and Graph-EFM (ms). In the Graph-EFM (ms) model samples tend to be more patchy and poorly reproduce spatial features. We connect this to the fact that the randomness in Z' is more local, as is associated with the nodes in GMS which directly connect to the grid. While the randomness introduced from sampling Z' can spread to the full forecast area using the multi-scale graph, this is something the model has to explicitly learn. For matching marginal distributions grid cells it can be enough for Graph-EFM (ms) to just have the randomness in each mesh node impact the connected grid nodes locally. We can contrast this to Graph-EFM, where top of the hierarchical graph. For this randomness to affect the prediction it must necessarily be propagated down the hierarchy, which disperses it spatially by construction.
associated with nodes at the
We have observed that the CRPS fine-tuning (see appendix D.2) can be central to the performance of Graph-EFM. This aligns the model distribution with the distribution the data. However, since this objective only encourages matching of the marginal distributions it is not sufficient for making model capture full joint distribution. This has to be achieved by constraints to the model or through additional parts of the training objective. We have observed that the weighting AcRPS used for the CRPS term can have a large impact on the spatial coherency forecasts. When AcRPS chosen too high the models trade off the ability to generate meaningful large-scale patterns for capturing
26

---

(a) nlars
2r
(c) z1000
Figure 9: Examples of artifacts appearing in Graph-EFM (ms) forecasts when trained with too high ACRPS the MEPS data. The checkerboard-lik pattern can be related to the layout of the multi-scale graph used.
(a) presOs Graph-FM
(b) presos GraphCast*
(c) u65
(d) 2r GraphCast*
Graph-FM
Figure 10: Examples of artifacts in smaller versions the deterministic models. The circular artifacts in GraphCast become much more prevalent when the model capacity is limited. These models use processing steps, d2 64 and are trained using NLL loss.
local variations. This is especially severe for the Graph-EFM (ms) model, forcing us to use lower ACRPS values in order to still get physical-looking forecasts. Figure shows example of the type of local noise-patterns that can appear when training Graph-EFM (ms) with a too high ACRPS-
Another problem observed with the multi-scale graph GMs in the LAM setting is the appearance of circular artifacts over the forecasting area. These appear in samples from Graph-EFM (ms) (see figs. 8b and (d). but also be found in deterministic forecasts from GraphCast*. These artifacts can be traced to the heterogeneous structure of the multi-scale graph. Due to how GMS constructed, mesh nodes can have different number of neighbors. The artifact positions match the positions of mesh nodes with many neighbors. Note that the hierarchical graph does not have this issue, as G1 (that connects to the grid) has more uniform structure. We noticed that these artifacts even more prevalent for smaller models, as shown in fig. I0 This points to the fact that with increased capacity, the GraphCast* model
some extent learn compensate for this problem.
G.2 Global Forecasting with ERA5
We have found it less challenging to achieve spatially coherent forecasts in the global setting, Both Graph-EFM and Graph-EFM (ms) generally produce samples with realistic physical features. Early in the training process observe some hexagonal patterns for longer lead times (also noted by Keisler (20), but these disappear as we train the models on longer rollouts.
One challenge that we encountered with the global probabilistic models was the appearance of small areas of instabilities when rolling out the models to longer lead times. In fig. [T we show some examples of this in forecasts from preliminary models. These unphysical artifacts typically covered
27

---

(a) msl
(b) msl
Graph-EFM
Graph-EFM (ms)
(c) u700
(d) u700
Graph-EFM
Graph-EFM (ms)
Examples of artifacts in preliminary version the global models unrolled to 10 days lead time. Small patches of deviating values appear and impact multiple variables. While these were more prevalent in Graph-EFM (ms), we did also observe this issue in Graph-EFM.
Figure 11:
only a few grid-cells, but are of course undesirable in a forecasting model. To remedy this for Graph-EFM found it sufficient train longer rollouts. The stabilizing effect of unrolling up to T 12 time steps during training did suppress these issues. This was however not enough to solve the problem for Graph-EFM (ms), for which additionally needed further lower ACRPS.
H Experiment Details: Global Forecasting with ERA5
In this appendix we give further detail about the graphs, data and experimental setup used in our global forecasting experiment.
H.1 Graph Construction
Multi-scale graph The global multi-scale graph is created in the same way in GraphCast 23
by recursively splitting the faces of an icosahedron. As we work with data on coarser resolution, we perform only 4 steps of such splitting. resulting in the graphs Is..... G1. These are all then merged to create GMs, used by GraphCast* and Graph-EFM (ms). By splitting 4 times we end up with ratio of N Rs 11 between grid nodes and mesh nodes. This can be compared with the 6 splitting [ steps of GraphCast, resulting in ratio with splitting only 3 times (resulting in N Rv. 45), but models using these graphs showed inferior
[PM] 25 for their 0.25" We initially experimented also
performance.
Hierarchical graph As the original icosahedron 95 contains very few nodes, we do not use this for the hierarchical mesh graph, but rather construct the hierarchy using 94...., G1. Edges E between levels are constructed by connecting each mesh node at level with nodes at level + 1 within a distance of 1.1 times the edge length in Gr. This guarantees that each node at level / has or 2 connection to the level above (see fig. 12g) Downward edges E are created by simply flipping the edge directions E1,41. We visualize global mesh graphs in fig. 12 and list the exact number of nodes and edges in table
28

---

(a) 91
(b) Gz
(d)
(e) Multi-scale mesh graph GMs
(I) Hierarchical mesh graph
Inter-level graph Ga.4
Figure 12: Mesh graphs used in the global experiment. Note that the vertical positioning (away from earth's surface) is purely for visualization purposes.
29

---

Table 4: Global graph statistics.
Graph
Nodes
Edges
G1
2562
15 360
G2
642
3840
G3
162
960
G4
42
240
Hierarchical
9
4482
G2.3/9022
1122
93,4/94,3
282
Total
3408
32 172
GMS
2562
20.60
GG2M
46 158
GM2G
120
Grid
29
Connecting the grid and mesh Edges connecting the grid and mesh graph also constructed following Lam et al. [23]. The edges EG2m are created by connecting each grid node to mesh nodes within distance of 0.6 times the edge length in Gi. Edges EM2G are constructed by for each grid node finding the triangle in G1 containing it, and adding three edges from the corner mesh nodes of that triangle. Note that the edge sets EGZM and EM2 identical for the multi-scale and hierarchical mesh graphs, since GMS and 91 contain the same nodes. Since the grid nodes are laid out in latitude-longee longitude grid and the mesh is icosahedral all mesh nodes will not be connected to the same number of grid nodes. Specifically, mesh nodes close to the poles will have many more grid connections than mesh nodes around the equator.
Static features All mesh nodes are associated with static features related to their position, specifi- cally cos(Latitud)), sin (Longiitd)))) and cos(gggudddd)), Edges (in the mesh, &G2M and MMG. have static features containing the edge length and vector difference between its endpoints (see Lam et al. for details). Edge features are normalized by the length of the longest edge.
23
H.2
Dataset Details
We version of ERA5 17 re-gridded to 1.5 latitude-longitude gridding, provided through the
WeatherBench 2 benchmark 40). Data from the period 1959-01-01T12. to 2017-12-31TI2 is used for training. 2017-12-31T18 to 2019-12-31T12 for validation and 2019-12-31T18 2021-01-10T18 for testing. These exact time stamps guarantee that there is no overlap between the subsets. For evaluation consider forecasts initialized times 00 and 12 UTC each day of 2020, following WeatherBench 2. We perform 10 day forecasts using h time steps. During training we also start forecasts at times 06 and 18 UTC.
The exact set of forecast variables, forcing and static fields is listed in table 3 We use the same inputs as Lam et al. [23]. The forcing is windowed over three consecutive time steps. meaning that each F contains forcing from times 1, and t + as well as the the static fields. For training we rescale the values of each variable to zero mean and unit variance.
H.3 Model and Training Contigurations
We train all models using the AdamW optimizer 57 and utilize BFloat16 mixed precision to save GPU memory. The training costs for the models makes extensive hyperparameter tuning unfeasible. We choose hyperparameters based on initial experimentation with smaller models. For Graph-EFI the important weightings Åkr. and ACRPS L can chosen based on monitoring the model behavior during training. This was done by plotting example forecasts throughout the training process and monitoring metrics on the validation set. The weight ÅKI. was tuned to prevent the model from collapsing to deterministic predictions, and to still achieve useful predictions when Z' was sampled
30

---

Table Variables, forcing and static fields from ERA5 used in our global forecasting experiments. kg of water vapour per kg of air.
Abbreviation
Unit
Vertical Levels
Variables
Geopotential
m²ls
Specific humidity
kg/kg'
50 hPa, 100 hPa, 150 hPa, 200 hPa, 250 hPa, 300 hPa, 400 hPa,
Temperature
K
t-component of wind t-component of wind Vertical velocity
m/s
500 hPa, 600 hPa, 700 hPa,
850 hPa, 925 hPa, 1000 hPa
m/s
Pa/s
Temperature
K
m above surface
u-component wind e-component of wind
m/s
10 m above surface
m/s
m above surface
Mcan sea level pressure
msl
Pa
Sea level
Total precipitation (Acc. over h)
tp
m
Surface
Forcing
Top of atm. solar radiation flux Sine-encoded time of day Cosine-encoded time of day Sine-encoded time of year Cosine-encoded time of year
toa
W/m Top of atmosphere
sin_tod tod
cod
sin_toy
cos_toy
Static Fields
Land-sea mask
water topography cos_lat
[0, mřl²
Surface Surface
Surface geopotential (topography) cos cos(Latiuud))
sin (Longitude)
sin_lon
cos (Longiitte))
cos lon
Table 6: Details of model architectures and training times for global forecasting.
Model
d₂
Processing steps
Parameters
Training time (GPU-hours)
GraphCast*
256
8
5.2 106
716
Graph-FM
256
30.9 106
1040
Graph-EFM (ms)
256
Pz: 2, g: 4, q: 4
7.7 106
1372
Graph-EFM
256
16.3 105
1264
from the latent map. The CRPS weight AcRpS tuned to achieve good ensemble spread while avoiding artifact issues, as discussed in appendix G
The exact model configurations and training times for our global experiments are listed in table 6 Parameter counts for probabilistic models include parameters in variational approximation, although this component does not play role during forecasting. We report training times as hours of active computations on a single GPU, In practice we use 80GB NVIDIA A.00 GPUs in parallel, meaning that the wall-clock time of our training is given by dividing the numbers in table by 8. For the Graph-FM model one processing step on the mesh graph constitutes a full pass through the hierarchy (either from the grid up or from level L down). In Graph-EFM (ms) multiple GNN processing steps operate on GMs- We list the number of such steps in table 6 separately for the latent map (μz). predictor (g) and variational approximation (q).
train all models by sequence training step, starting from single-step prediction and then unrolling predictions 7 steps. For the probabilistic models we first train in pure auto-encoder setup with AgI. = 0, encouraging y to encode useful information in the distribution over Zt. We find this useful for preventing the model from ignoring Z' and collapsing purcly deterministic forecasting.
31

---

Table 7: Training schedule for deterministic models (GraphCast* and Graph-FM) on global data.
Epochs
Learning Rate
Unrolling T
70
10-3
20
10
4
20
10
8
Table 8: Training schedule for Graph-EFM on global data. For Graph-EFM (ms) we use the same
schedule but with different constants (Agı. 1, ACRPS 10%).
Epochs
Learning Rate
Unrolling T
AKL
ACRPS
50
10-3
1
0
0
75
10-1
1
0.1
0
10-4
4
0.1
10
10
8
0.1
2
10-4
12
0.1
104
The probabilistic models additionally include the fine-tuning using CRPS as a final training step. The full training schedule for the deterministic models is given in table Note that we here define one epoch as initializing the model at each possible time in the training set (00, 06, 12 and UTC in each day) such that all unrolled lead times are still within the training set period. The reason for probabilistic models being unrolled to longer lead time T is mainly to combat the artifacts discussed in appendix G.2
for the probabilistic models
I Experiment Details: Limited Area Modeling with MEPS Data
In this appendix we give further detail about the graphs, data and experimental setup used in our experiments with the MEPS data.
L1 Graph Construction
Multi-scale graph In this limited-area setting we construct graphs as regular quadrilateral meshes covering the rectangular MEPS forecasting area. To create these we lay out mesh nodes in regular rows and columns over the area. Each node then connected with bidirectional edges its neighbors horizontally, vertically and diagonally (see fig. 1 This results in all nodes (except those at the edge of the area) having neighbors. The procedure repeated at 4 different resolutions, tripling the distance between nodes at each resolution. This means that a node at resolution level is positioned the center of group of 3 x nodes at resolution level 1, sharing its exact position with the center node of the group (illustrated in fig. 13b). To create the multi multi-scale mesh graph GMS we then merge the graphs at different resolutions, combining any nodes that sit the same coordinates into one node. Note that this is possible due to how nodes align across the resolution levels.
Hierarchical graph For the hierarchical model the graphs of different resolution are not merged, but used the different levels of the hierarchy. We include only the 3 finest meshes G1, G2 and Gs. as G4 contains only 9 nodes. Additional inter-level edge sets Σ1,41 and E1+,,1 are then created for all adjacent levels. Each set ε1,4.1 upwards edges created by connecting each node on level with the closest node on level 1 + This means that each node at levels >>>>>>>>>>>>> 1 will have 9 incoming edges from the level below. The downward edges 1, are copy of ε with the direction of each edge flipped. The mesh graphs used for the MEPS experiment are visualized in fig. 14 and the number of nodes and edges each graph listed table
Connecting the grid and mesh To form εGZM, each grid node connected to mesh nodes closer than 0.67 times the distance between nodes in 9. All distances are here 2-dimensional euclidean, computed in the MEPS Lambert projection coordinates. The set EMpg constructed by iterating over the grid nodes, and at each creating edges from the closest mesh nodes to the node in the grid.
32

---

(a) Each mesh node a is connected to its neighbors horizontally, vertically and diagonally.
(b) Alignment of mesh nodes G/ with mesh nodes in Gi-1.
Figure 13: Illustration of mesh graph construction in the LAM setting.
(a) Multi-scale mesh graph Gus
(b) Hierarchical mesh graph
(c) Mesh levels G1, Ga, Go
Figure 14: Mesh graphs used in the MEPS experiment. Note that the vertical positioning and size of nodes is purely for visualization purposes.
33

---

Table 9: MEPS graph statistics.
Graph
Nodes
Edges
91
6561
51
G2
729
5512
G3
81
544
Hierarchical
g1,2/92,1
6561
92,299,2
729
Total
7371
156
GMS
6561
67616
9G2M
150
GM2G
255 136
Grid
63 784
2022
Apr May June
Sep
Oct
Nov
Dec
Jan
Mar
Apr May June July Aug Sep
Oct
Nov
Dec
Jan
Feb
Mar
Training
Validation
Test
Figure 15: Overview of the period covered the dataset and the training/validationness split.
Static features The static features associated with mesh nodes their 2-dimensional coordinates in the MEPS Lambert projection, normalized by the maximum coordinate value. All edges have static features including the length of the edge and the vector difference between the source and target nodes (using the projection coordinates). The edge features are normalized by the length the longest edge in the whole mesh graph.
1.2
Dataset Details
The MEPS dataset consists of archived forecasts from the operational MEPS system during the period April 2021 March 2023. This period was chosen due to the system configuration being reasonably stable, preventing distributional shifts within the dataset. From the chosen period we extract the forecasts started at 00 and 12 UTC each day. At each initialization time there are ensemble forecasts, started from slightly different initial conditions. This results in 10 forecasts per day (assuming no ensemble members missing due to operational issues). When retrieving the data we additionally downsample the spatial resolution from the original 2.5 km to 10 km. This results in dataset of 6069 forecasts of length 66 h with h time steps. We split the forecasts into training, validation and test sets according to fig. 15) The specific validation months were chosen to reasonably cover the seasonal variations. Note that the dataset contains around 46 years of individual time steps. However, since there are obvious correlations between ensemble members and successively started forecasts the actual information content far lower than in 46 year reanalysis dataset.
Dataset
Variables and forcing At each grid cell we model 17 weather variables, including broad range of different quantities and different vertical levels in the atmosphere. All variables, forcing and static fields are described in table TO The particular choice of variables was motivated by combination of metcorological relevance. data availability and striving for diverse set of variables to evaluate the model on. We use the same type of windowing for the forcing and standardization of variables as in the global experiment. For solar radiation (nlwrs and nswrs) level, aggregated over the past 3 hours (since the last time step). Apart from the solar radiation all other variables are instantaneous. For the MEPS data we use the fraction of open water in the grid cell as a forcing input. We assume this to be constant over the forecast period and take the value from time of the initial state. the global experiment the land-sea mask is static, forcing could be useful for taking into account seasonal fluctuations of the ice cover in the Nordic region. The boundary forcing Bd consists of the same variables as listed in table 10 We include static binary indicator variable describing if a node is in the boundary or forecast area.
consider the flux at ground
treating this as
34

---

Table 10: Variables, forcing and static fields in the MEPS dataset. 'In the MEPS system 65 vertical model levels defined from the ground to the top of the atmosphere [32] The lowest MEPS level (Ly165) sits at approximately 12.5 m above ground.
Abbreviation
Unit
Vertical Levels
Variables
Atmospheric pressure
pres
Pa Wm3
Ground level (0g), Sea level (0s)
Net longwave solar radiation flux Net shortwave solar radiation flux Relative humidity
nlwrs
Surface
nswrs
W/m²
Surface
r
[0, I]
2 m, Lv165' (65)
Temperature
t
K
m. Lv165 (65), 500 hPa, 850 hPa
u-component of wind
u
m/s
Lv165 (65), 850 hPa
v-component of wind
V
m/s
Lv165 (65), 850 hPa
Water vapor, integrated column Geopotential
wvint
kg/m²
Full column
2
m
500 hPa, 1000 hPa
Forcing
Top of atm. solar radiation flux Fraction of open water Sine-encoded time of day Cosine-encoded time of day Sinc-encoded time of year Cosine-encoded time of year
toa
W/m
Top of atmosphere
water
[0,
Surface
ain_tod
cos_tod
sin_toy
cos_toy
Statie Fields
Surface geopotential (topography) u-coordinate in MEPS projection y-coordinate in MEPS projection On-boundary binary indicator
topography x_oor yooord boundary
m²ls
Surface
0/1
Forecast steps and length The original data uses h time steps, but our MLWP models predict in 3h steps. Because of this we can extract 3 training samples from each forecast in the dataset (i.c. original time steps {1, 4, ...... {2, 5, ... and (3, 6,,,..... ...,..... .). As we train only T < 8 rollout steps, we also only need a subset each such time series for each training iteration. To make maximum of the data during training randomly sample which time step for unrolling the model. The combination of 1) the 3 time steps, 2) using the last two weather states model inputs, 3) including forcing from multiple times input, means that reduce the effective length of our ground truth forecasts in pre- pre-processing. This explains why we predict for 57 rather than the original 66 h. Note however that nothing prevents forecasts for 66 h or beyond. Although this is possible, we do not have ground truth data to compare against past 57 and therefore view such experiments to be of limited interest. There is still reason to expect the model performance to become drastically worse specifically past 57 h, as all models are anyhow fine-tuned on shorter rollouts than this.
as initial state
from unrolling the models create
1.3 Boundary Forcing
In all models we use boundary forcing Be in order to include information about the surrounding area. At each time step boundary forcing Bt-2:t11 passed to the model together with Xt-2:t- , as described in fig. [6] There is a separate set of grid nodes used for the boundary input. These nodes treated identically as grid nodes within the forecasting area by the GNN layers. As the boundary forcing is only a change to the model inputs this does not require any substantial re-design and we can adapt all models in our experiments to this LAM setting. Note that we use a boundary area that lays inside the original MEPS area, allowing us to use parts of the ground truth forecasts boundary forcing. We specifically define the boundary as the 10 grid cells closest the edge of the limited area. In the operational system the boundary area is defined along the edge outvide the MEPS area. There is no major conceptual difference between these and one could easily re-define the different areas to match the operational MEPS system.
35

---

Xt
Model
Xt-2:t-1 B
AB B²
F
X B
B
Figure 16: Schematic showing how we feed boundary forcing as input to the model in each autore- gressive step.
Table 11: Details of model architectures and training times for LAM forecasting.
Model
d2
Processing steps
Parameters
Training time (GPU-hours)
GraphCast*
128
1.1 x 106
336
Graph-FM
128
6
4.5 106
432
Graph-EFM (ms)
128
1zz 2, g: 4, 9:2
1.7 x 106
456
Graph-EFM
128
3.1 x 106
556
1.4 Model and Training Configurations
For the MEPS experiment we again use the AdamW optimizer [57], but not mixed precision compu- tations. Hyperparameter tuning follows a similar strategy as for the global experiment. The exact model configurations and training times for our MEPS experiments are listed in table The training schedule for the deterministic models is given in table 12 and for the probabilistic models in table [3 These tables follow the same formats as tables
and we refer to global experiment details in H.3 for further explanations. In the MEPS case epoch means training one sub-sample (see description of sub-sampling in appendix I.2) of each forecast in the training set once.
appendix H.3
Additional Results: Global Forecasting with ERA5
In this appendix present additional results for global forecasting.
J.1 Metrics
Comparisons between existing MLWP models are challenging, due the many factors that impact forecast quality. While most global models proposed in the literature are trained ERA5, the
Table 12: Training schedule for deterministic models (GraphCast* and Graph-FM) on MEPS data.
Epochs Learning Rate
Unrolling T
600
10-3
300
10-4
200
10-4
8
36

---

Table 13: Training schedule for Graph-EFM on MEPS data. For Graph-EFM (ms) we use the same schedule but with different constants (AKı. = 100. AcRPS = 10", fine-tuning learning rate 10 +).
Epochs Learning Rate
Unrolling
XKL
ACRPS
300
10-3
1
0
0
100
10-3
1
1
200
x 10
4
50
5 x 10
4
1
106
100
5 x 10-4
6
106
exact choice of train/test split, spatial resolution and forecasted variables can vary (compare for example Keisler (20), Pathak et al. [37 and Hu et al. [1]). Models that operate on a higher spatial resolution and with more variables get expected to produce better forecasts. When it comes to MLWP ensemble forecasting, an even greater challenge is that different initial conditions are used [39 6 [6]. Given this situation, it is hard to disentangle the proposed machine learning methods from their surrounding design choices and make fair comparisons of model architectures.
information in initial states, and should therefore be
We approach these complications by re-training a smaller set of models on the same data and experimental setup. Our focus is graph-based, non-hybrid MLWP models, this is the regime of Graph-FM and Graph-EFM. The baseline models are described in detail in appendix [C.5) In this appendix also include additional metrics taken directly from the WeatherBench 2 Benchmark [40]
These are for the original GraphCast model
[23]. KeislerNet (20] and the IFS-ENS operational
ensemble from the European Centre for Mediumm-aange Weather Forecasts (ECMWF) [12. The evaluation setup for these results match ours (evaluation against 2020 from ERA5, identical metrics computed on 1.5 data). However, as discussed above there are differences in the exact training data, variables and operating resolutions compared to our models.
Metric values from our global experiments presented in figs. 19 We here showcase values for surface variables and atmospheric variables at pressure levels 500, 700 and 850 hPa. The results taken directly from WeatherBench 2 shown in gray, emphasize that any comparison with these comes with caveats. Note that results for all variables are not available for all models.
J.2 Example Forecasts
In fig. 20 we showcase example forecasts from different models at lead time 10 days. Note that forecasting this far into the future challenging, and sampled ensemble members show only one possible scenario predicted by the models.
Figure 21 shows example ensemble forecasts from Graph-EFM for surface variables, atmospheric variables at 700 hPa, z500 and t850. All forecasts are for 10 days lead time. For practical reasons we use only 20 members estimate the ensemble mean and standard deviation these plots.
37

---

GraphCast*
Graph EFM (ms))
GraphCast* +SWA
GraphCast
KeislerNet
ENS
Graph FM
Graph EFM
3
2
3
9
10
4
7
10
5
8
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) 2t
(b) 10u
(c) 10v
750
2F-3
500
26-3
300
250
1B-3
5B-4
10
7
10
2
9
Lead time (Days)
Lead time (Days)
Lead time (Days)
(d) msl
(e) tp
(f) z500
88:4
7.5
2
5B-4
5.0
3F.4
0
3
6
9
10
2
4
7
8
10
3
5
7
9
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
(g) t500
(h) 9500
(i) u500
10
750
0.175
500
250
0.125
0
10
5
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
(i v500
(k) w500
(1) 2700
2B-3
16:3
2
5B-4
0
3
4
6
7
9
10
2
4
5
7
8
10
2
4
5
7
8
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
(m) t700
(n) 9700
(o) u700
38

---

GraphCast
Graph-EFM
GraphCast*
GraphCast
KeislerNet
IFS-ENS
Graph-FM
Graph-EFM
600
0.175
400
0.150
200
0.125
1
2
4
5
10
4
8
2
6
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
(p) v700
(q) w700
(r) 2850
26-3
6
1E.3
2
9
9
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
(s) t850
(0) q850
(u) u850
0.16
0.12
0.10
3
9
Lead time (Days)
Lead time (Days)
(v) v850
(w) w850
Figure 17: RMSE results for global experiment.
39

---

GraphCast*
Graph EFM (ms))
GraphCast*+S
GraphCast
KeislerNet
ENS
Graph FM
Graph EFM
2:0
1.5
1.0
0.5
2
6
8
10
4
7
10
5
8
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) 2t
(b) 10u
(c) 10v
600
8F-4
6B-4
400
200
4E-4
200
2B-4
0
10
7
10
7
Lead time (Days)
Lead time (Days)
Lead time (Days)
(d) msl
(e) tp
(f) z500
6E-4
2
4b-4
0
3
6
9
10
2
4
7
10
2
5
7
8
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
(g) t500
(h) 9500
(i) u500
0.12
6
400
0.10
1
0.08
200
2
0.06
0
3
10
10
2
7
Lead time (Days)
Lead time (Days)
Lead time (Days)
(i) v500
(k) w500
(1) 2700
2
1B-3
5E-4
0
3
4
6
7
9
10
2
4
5
7
8
10
2
4
5
7
8
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
(m) t700
(n) 9700
(o) u700
40

---

GraphCast
Graph-EFM Graph-EFM
GraphCast*
GraphCast
KeislerNet
IFS-ENS
Graph-FM
0.12
400
0.10
200
0.08
0.06
1
2
4
5
10
3
6
9
2
6
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
(p) v700
(q) w700
(r) 2850
4
1E.3
5E-4
9
9
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
(s) t850
(0) q850
(u) u850
0.06
3
9
Lead time (Days)
Lead time (Days)
(V) v850
(w) w850
Figure 18: CRPS results for global experiment.
41

---

Graph-EFM (ms)
Graph-EFM
GSWA
Calibrated
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
2
4
6
8
9
10
2
3
6
8
10
1
3
5
6
7
8
9
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) 2t
(b) 10u
(c) 10v
1.0
1.0
0.5
0.5
0.5
0.0
1
0.0
0.0
2
6
8
10
10
3
5
Lead time (Days)
Lead time (Days)
Lead time (Days)
(d) msl
(c) tp
(f) 2500
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
2
4
6
8
10
2
3
6
8
10
3
5
7
8
9
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
(g) t500
(h) 9500
(i) u500
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
2
4
6
8
10
8
10
3
5
7
Lead time (Days)
Lead time (Days)
Lead time (Days)
v500
(k) w500
(1) 2700
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
2
4
6
8
10
2
6
10
1
3
5
7
10
Lead time (Days)
Lead time (Days)
I Lead time (Days)
(m) t70
(n) q700
(o) u700
42

---

Gruph- BP
Graph- EFM
GSSA
Calibrated
1.0
0.5
0.5
0.5
0.0
0.0
Lead time (Days)
Lead time (Days)
Lead time (Days)
(p) v700
(q) w700
z850
1.0
1.0
1.0
0.5
0.5
0.0
0.0
0.0
3
5
10
3
5
7
9
10
2
5
9
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
t850
(1) q850
(u) u850
0.5
0.5
0.0
3
5
Lead time (Days)
Lead time (Days)
(v) v850
(W) w850
Figure 19: SpSkR results for global experiment.
43

---

(a) 10u
(b) 10u GraphCast*
(c) 10u Graph-FM
Ground Truth
(d)
Graph-EFM (ms), Ens. Member
Graph-EFM Ens. Member
(1) 9700 Ground Truth
(g) 9700 GraphCast*
(h) q700 Graph-FM
(i) q700 Graph-EFM (ms), Ens. Member
(i) q700 Graph-EFF Ens. Member
(k) z500 Ground Truth
(1) z500 GraphCast*
(m) Graph-FM
(n) z500 Graph-EFM (ms), Ens. Member
(o) z500
Graph-EFM, Ens. Member
Figure 20: Comparison of global model forecasts at 10 days lead time for u-component of 10 m wind (10u), specific humidity 700 hPa (q700) and geopotential at 500 hPa (z500). For probabilistic models we show sampled ensemble members.
44

---

Ground Truth
Ens. Mean
Ens. Std.-Dev.
10.0
Ens. Members
(a) 2t
Ground Truth
Ens. Mean
Ens. Std. Dev.
96 100
Ens. Members
(b) 10u
Ground Truth
Ens. Mean
Ens. Std.-Dev.
100
Ens. Members
-30
(c) 10v
45

---

Ground Truth
Ens. Mean
Ens. Std.-Dar
X00
-180
1:40
-5.0
Ens. Members
102000
100000
130000
(d) msl
Ground Truth
Ens. Mean
Ens. Std.-Dev.
0.005
Ens. Members
(e) tp
Ground Truth
Ens. Mean
Ens. Std. Mev.
200)
-000
2003
Ens. Members
23000
42000
(l) z500
46

---

Ground Truth
Ens. Dev.
Ens. Members
Sa sapuer
t850
Ground Truth
Ens. Mean
Ens. Std.-Det
3.015
Ens. Members
3016
2.000
(h) q700
Ground Truth
Ens. Mean
Ens. Std. Mev.
(200)
BAR
2000
300
Ens. Members
25000
2200
(i 2700
47

---

Ground Truth
Ens. Mean
Ens, Std.-Dev.
Ens. Members
Raz
280
200
ii) t700
Ground Truth
Ens. Mean
Ens. Std. Ael.DDev
Std Dec
Ens. Members
(k) u700
Ground Truth
Ens. Mean
Ens. Mev.
Sod. DWN.
Ens. Members
31
20
(l) v700
48

---

Ground Truth
Ens. Mean
Ens. Std. Dev.Dev.
0.
Ens. Members
95
-25
(m) w700
Figure 21: Example Graph-EFM global ensemble forecasts at lead time 10 days.
49

---

K
Additional Results: Limited Area Modeling with MEPS Data
In this appendix we show additional results from our experiment with the MEPS data.
K.1
Metrics
Figures
22 to 24 show metric values for all variables and lead times in the MEPS dataset.
The poor performance of the Graph-EFM (ms) model is noteworthy, both in terms of forecast accuracy and ensemble calibration. This can to a large extent be attributed to the exact training objective (see table 13), particularly lower weighting ACRPS- Using a lower value for ACRPS in the multi-scale model was necessary to avoid the artifacts discussed appendix G The does however mean that Graph-EFM (ms) does not gain much of the benefits that come with the CRPS fine-tuning.
K.2 Example Forecasts
comparison between example forecasts from different models is given fig. 25 Note that ensemble members sampled from Graph-EFM show more detailed features than deterministic forecasts.
In fig. 26 we plot example ensemble forecasts from the Graph-EFM model for all variables in the MEPS data. Note that these plots include the boundary area, which is not forecast, as a faded border in each plot. All forecasts for lead time 57 Rolling out LAM ensemble forecast from Graph-EFM to 57 h is even faster than generating the global ensemble. Using batched sampling on single GPU, we can produce 100 ensemble members 140 (1.4 per member).
GraphCast
Graph-FM
Graph-EPM (ms)
Graph-EFM
GraphCast*
600
4F+5
400
200
2F+5
2
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) pres0g
(b) presOs
(c) nlwrs
1B+6
0.10
0.100
0.08
0.075
5B+5
0.06
0.050
0.04
Lead time (Days)
Lead time (Days)
Lead time (Days)
(d) nswrs
(e) 2r
(f) r65
50

---

CruphCast
Gruph-PM
Graph-EFM (ms)
EHM
GraphClst**SSW.
2.0
1.5
2
1.0
Lead time (Days)
Lead time (Days)
Lead time (Days)
(g) 2t
(h) t65
(i) t500
Lead time (Days)
Lead time (Days)
Lead time (Days)
(j) t850
(k) u65
(I)
Lead time (Days)
Lead time (Days)
Lead time (Days)
(m) u850
(n) v850
(0) wvint
600
100
100
200
200
2
Lead time (Days)
Lead time (Days)
(p)
(q) z1000
Figure 22: RMSE for MEPS experiment.
51

---

GraphCast
Graph-FM
Graph-EPM (ms)
Graph-EFM
GraphCast-SSW.
400
100
4E+5
200
200
2E+5
1B+5
2
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) pres0g
(b) presOs
(c) nlwrs
6B+5
0.06
4B+5
0.01
0.04
2B+5
0.02
0.02
2
Lead time (Days)
Lead time (Days)
Lead time (Days)
(d) nswrs
(c) 2r
(f) r65
1.0
1.0
1.0
0.5
0.5
Lead time (Days)
Lead time (Days)
Lead time (Days)
(g) 2t
(h) t65
(i) t500
1.5
2.0
2.0
1.0
1.5
1.5
0.5
0.5
0.5
Lead time (Days)
Lead time (Days)
Lead time (Days)
(j) t850
(k) u65
(l) v65
2.0
3-
1.5
2
1.0
0.5
2
Lead time (Days)
Lead time (Days)
Lead time (Days)
(m) u850
(n) v850
(o) wvint
52

---

GraphCast
Gruph FM
Graph EFM (ms)
Graph EFM
FFM
GraphCast*-SAA
300
200
100
Lead time (Days)
Lead time (Days)
(p) z500
(q) z1000
Figure 23: CRPS for MEPS experiment.
Graph-EFM (ms)
Graph-EFM
GraphCast SSWA
Calibrated
1.0
1.0
1.0
0.5
0.5
0.0
0.0
0.0
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) pres0g
(b) presOs
(c) nlvrs
1.0
0.5
0.5
0.5
0.0
0.0
0.0
Lead time (Days)
Lead time (Days)
Lead time (Days)
(d) nswrs
(c) 2r
(f) r65
53

---

Graph-EFM (ms)
Gruph-EFM
GraphiCast* 4SWA
Calibrated
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
Lead time (Days)
Lead time (Days)
Lead time (Days)
(g) 2t
(h) t65
t500
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
Lead time (Days)
Lead time (Days)
Lead time (Days)
(j) t850
(k) u65
(I) v65
1.0
1.0
1.0
0.5
0.5
0.0
0.0
0.0
2
Lead time (Days)
Lead time (Days)
Lead time (Days)
(m) u850
(n) v850
(0) wvint
1.0
1.0
0.5
0.5
Lead time (Days)
Lead time (Days)
(p) z500
(q) z1000
Figure 24: SpSkR for MEPS experiment.
54

---

(a) u850
(b) u850 GraphCast*
(c) u850 Graph-FM
Ground Truth
(d) u850 Graph-EFM (ms) Ens. Member
(e) u850 Graph-EFM Ens. Member
(f) 2r Ground Truth
(g) 2r GraphCast*
(h) 2r Graph-FM
(i) 2r
(j) Graph-EFM Ens. Member
Graph-EFM (ms) Ens. Member
Figure 25: Comparison of LAM forecasts for u-component of wind 8P0 (u850) and m relative
humidity (2r) at 57 h lead time. For probabilistic models we show sampled ensemble members.
55

---

Ground Truth
Ens. Mean
Ens. Std.-Dev.
200
Ponk AR
150
Std. Dev.
-50
Ens. Members
1000000
95000
Phi
30000I
85000
S0000
(a) presOg
Ground Truth
Ens. Mean
Ens. Std.-Dev.
200
150
Sld- Dec.
100
500
Ens. Members
1022000
101500
1010000
1005500
1000000
94500
(b) presOs
56

---

Ground Truth
Ens. Mean
Ens. Std.-Dev.
1000000
3000000
Sid. Dew
2000000
1000000
Ens. Members
100
0.5
Whm
-15
(c) nlwrs
Ground Truth
Ens. Mean
Ens. Std.-Dev.
1.75
1.50
1.25
Std. Dev.
0.75
0.50
Rict
0.00
Ens. Members
10
2
nswrs
57

---

Ground Truth
Ens. Mean
Ens. Std.-Dev.
0.177
0..50
0.125
Std. Dev.
0.100
0.075
0.050
0.025
0.0000
Ens. Members
0.8
0.6
0.4
(c) 2r
Ground Truth
Ens. Mean
Ens. Std.-Dev.
0.175
0.150
0.125
Std. DeV.
0.100
0.075
0.050
0.025
0.000
Ens. Members
1.0
0.9
0.8
0.7
0.6
D.4
0.3
(f) x65
58

---

Ground Truth
Ens. Mean
Ens. Std.-Dev.
2.5
Std. Dev.
1.5
1.0
0.5
0.0
Ens. Members
900
295
290
285
280
275
270
205
(g) 2t
Ground Truth
Ens. Mean
Ens. Std.-Dev.
2.0
Dev.
1..0
0.5
0.0
Ens. Members
295
290
285
280
275
270
(h)
59

---

Ground Truth
Ens. Mean
Ens. Std.-Dev.
1.75
1.50
1.25
0.50
0.25
0.00
Ens. Members
200
255
250
245
ii +500
Ground Truth
Ens. Mean
Ens. Std.-Dev.
2.0
Dex
1.0
0.5
0.0
Ens. Members
200
285
280
275
270
265
i) t850
60

---

Ground Truth
Ens. Mean
Ens. Std.-Dev.
"Sdd-Dev Dev.Dev
Std
Ens. Members
15
10
5
-0
15
(k) u65
Ground Truth
Ens. Mean
Ens. Std.-Dev.
5
Std. Dev.
Ens. Members
TOR Ro PE
10
m/s
(1) u850
61

---

Ground Truth
Ens. Mean
Ens. Std.-Dev.
6
Std J.-Dev.
Ens. Members
15
5
-10
(m) v65
Ground Truth
Ens. Mean
Ens. Std.-Dev.
Std. Sev.
-2
Ens. Members
20
-0
-10
20
(n) v850
62

---

Ground Truth
Ens. Mean
Ens. Std.-Dev.
7
6
2
0
Ens. Members
30
(0) wvint
Ground Truth
Ens. Mean
Ens. Std.-Dev.
200
150
3dd Dev
100
Ens. Members
56000
55000
mn//
540000
53000
52000
(p) 2500
63

---

Ground Truth
Ens. Mean
Ens. Std.-Dev.
175
150
125
100
50
25
Ens. Members
1500
Ro M
1000
500
500
(g) z1000
Figure 26: Example Graph-EFM LAM ensemble forecasts at lead time 57 h.
64

---

L
Additional experiments
L.1
Comparing Interaction and Propagation Networks in Graph-FM
Given the usefulness of Propagation Networks in Graph-EFM it is reasonable to ask if these could be beneficial to use also in the deterministic Graph-FM model. Also Graph-FM uses a hierarchical mesh graph, so the propagation of information between levels is important also for this model. We test this empirically, by training versions of Graph-FM using Interaction and Propagation networks. As the retention of information is also important for some parts of the architecture, we do not replace all GNNs with Propagation Networks (see appendix (C.2). b. We compare models on both global and LAM forecasting, using the same experimental setups as described in appendices H and
Global forecasting with ERAS Figure 27 shows RMSEs for highlighted variables from ERAS and Graph-FM using Interaction and Propagation Networks. The models have very similar errors for almost all variables. However, for variables in the upper atmosphere the Propagation Network model shows lower errors.
Limited area modeling with MEPS data In fig. 28 we compare RMSEs for Graph-FM models using the different GNN layers on MEPS data. Here we see a greater advantage of the Propagation Networks. We hypothesize that this relates to the boundary forcing. Using Propagation Networks the information from boundary nodes should easier reach the top graph level, where it can faster spread throughout the forecasting area. Motivated by these results we Propagation Networks in our final Graph-FM architecture, both for global and LAM forecasting.
L.2 Importance of Latent Map
In Graph-EFM we define the latent map p pZZ [Z F²) using a neural network with explicit dependence on - Fr. Given that also the predictor takes xt-2:t-1, F as inputs, it is not immediately clear that conditioning on this also in the latent map is necessary. Indeed, due to the ability of deep neural networks to well approximate arbitrary functions, the predictor network should internally be able to transform simple Z' ~ (00 I) to introduce the dependence on Xt-2:t11, F This does however assume infinite flexibility in the predictor, which might be far from the situation in practice.
To investigate the importance of learnable latent map we compare our Graph-EFM model on the MEPS dataset with one where Z' is sampled from a static distribution N 1). This experiment was carried out using an earlier version of Graph-EFM with the same architecture but a slightly different training schedule. To save on computations we here only sample 16 ensemble members from each model. RMSE, CRPS and SpSkR are shown in figs. 29 to 3T In terms of RMSE and CRPS there is clear benefit to letting the distribution over Zt be a learnable mapping. The SpSkR fig. shows no clear trends between the models. For practical network architectures the latent map does add flexibility, changing the mean of Z' that enters the predictor. The learnable latent map should also simplify the inference problem solved when optimizing our variational objective. With the dependence on previous states we can expect smaller discrepancy between p F²) and p( Z' 21-2:t1, X F²). simplifying the optimization of the variational
approximation q X | 2:-1 empirical benefit of using the latent map.
X F*). We believe that this is an important aspect of the observed
L.3 Impact of Ensemble Size
We here study how the performance of the model, in terms of metric values, varies when sampling different numbers of ensemble members. To investigate this we ran the evaluation of Graph-EFM with 5 80 members for the global model and 5-100 for the LAM model. Results for a selection of variables are shown in fig. ( (global, ERA5) and fig. 33 (LAM, MEPS). As expected the RMSE of the ensemble mean decreases when sampling more members. However, already when sampling 20 or 25 members the results are fairly close to the full ensemble. For SpSkR the differences are even smaller. As the CRPS is property of the distribution of the model forecast, its true value does not depend on the number of samples drawn. In practice we compute CRPS using an unbiased estimator, and the variance of this estimator decreases with ensemble size. When averaged spatially and over
65

---

Propagation Networks
Interaction Networks
750
500
1.0
0.5
3
6
8
10
2
6
8
10
2
4
5
7
8
10
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) z500
(b) 9700
(c) t:50
2
0.018
Lead time (Days)
Lead time (Days)
Lead time (Days)
(d) 2t
(e) t100
(I) w100
Figure 27: RMSE of Graph-FM models with Propagation and Interaction Networks, evaluated on the
ERA5
Propagation Networks
Interaction Networks
400
200
2645
2
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) z500
(b) wvint
(c) nlwrs
1.5
1.0
0.5
Lead time (Days)
Lead time (Days)
Lead time (Days)
(d) 2t
(c) u65
( v65
Figure 28: RMSE of Graph-FM models with Propagation and Interaction Networks, evaluated on the MEPS test set.
66

---

Latent map
Static
1.5
1.0
0.5
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) 2t
(b) wvint
(c) u850
Figure 29: RMSE for Graph-EFM models with static distribution for Z or learnable latent map, evaluted on the MEPS validation dataset.
Latent map
Static
0.8
1.0
0.6
0.4
1.0
0.5
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) 2t
(b) wvint
(c) u850
Figure 30: CRPS for Graph-EFM. models with a static distribution for Zt or a learnable latent map. evaluted the MEPS validation dataset.
Latent map
Static
Calibrated
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
2
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) 2t
(b) wvint
(c) u850
Figure 31: SpSkR for Graph-EFM models with static distribution for Z learnable latent map, evaluted on the MEPS validation dataset.
67

---

-.....
5
10
20
40
80
750
300
1.0
300
0.5
0.0
10
3
10
3
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) 2500, RMSE
(h) 2500, CRPS
(c) 2500, SpSkR
1.0
1.0
2
0.6
0.5
3
10
10
5
Lead time (Days)
Lead time (Days)
Lead time (Days)
(d) 2t. RMSE
(e) 2t, CRPS
(l) 2t, SpSkR
Figure 32: Metric values for z500 and 2t from the global experiment when sampling different
numbers of ensemble members from Graph-EFM.
the whole we do however not see any difference in CRPS for different ensemble sizes. All these trends hold consistently for all variables in both the ERAS and MEPS datasets.
In main experiments in section 5 improvements to metrics saturate, we would not expect results to meaningfully change from sampling members than this. It should however be noted that the motivation for sampling very large ensembles is mainly not to improve on metrics such as these. More important motivations for large ensembles include estimating probabilities of rare events studying different possible scenarios of extreme weather.
the full 80/100 member ensembles. Given that any
68

---

5
10
25
50
100
1845
2B+5
1.0
81++4
2B+5
0.5
6F44
1E+5
0.0
Lead time (Days)
Lead time (Days)
Lead time (Days)
(a) nlwra, RMSE
(h) nlwrs, CRPS
(c) nlure, SpSkR
1.5
0.6
1.0
1.0
0.5
0.0
Lead time (Days)
Lead time (Days)
Lead time (Days)
(d) 2t. RMSE
(e) 2t, CRPS
(I) 2t, SpSkR
Figure 33: Metric values for nlwrs and 2t from the LAM experiment when sampling different numbers of ensemble members from Graph-EFM.
Appendix References
[53] J. L. Ba, J. R. Kiros, and G. E. Hinton. Layer normalization. arXiv preprint arXiv:1607.06450 2016.
[54] Z. B. Bouallègue, M. C. A. Clare, L. Magnusson, E. Gascón, M. Maier Gerber, M. Janoušek, M. Rodwell, Pinault, J. Dramsch, S. Lang, B. Raoult, F. Rabier, M. Chevallier, I. Sandu, P. Dueben, M. Chantry, and F. Pappenherger. The rise of data-driven weather forecasting: first statistical assessment of machine learning- ceanning-based based weather forecasts in an operational-lin like context. Bulletin of the American Meteorological Society, 2024.
[55] A. J. Charlton-Pere! H. F. Dacre, S. Driscoll, S. L.. Gray, B. Harvey, N. J. Harvey, K. M. R. Hunt, R. W. Lee, R. Swaminathan, R. Vandaele, and A. Volonté. Do Al models produce better weather forecasts than physics-based models? quantitative evaluation study of storm ciarán. npj Climate and Atmospheric Science, 2024.
[56) Core writing team, H. Lee, and J. R. (eds.). IPCC, 2023: Climate change 2023: Synthesis report. contribution of working groups I, II and to the sixth assessment report of the intergovern- mental panel on climate change. Technical report, Intergovemmentall Panel Climate Change (IPCC), 2023.
[57] I. Loshchilov and F. Hutter. Decoupled weight decay regularization. In International Conference on Learning Representations, 2019.
[58] National Weather Service. Hurricane laura, 2020. URL 2020Laura
[59] Ramachandran, B. Zoph, and Q. V. Le. Searching for activation functions. arXiv preprint urXiv:1710.0553 2017.
[60] R. Sims, P. Mereado, W. Krewitt, G. Bhuyan, D. Flynn, H. Holttinen, G. Jannuzzi, S. Khennas, Liu, L. Nilsson, J. Ogden, K. Ogimoto, M. O' Mallley H. Outhred, Ø. Ulleberg, and v. Hulle. Integration of renewable energy into present and future energy systems. In IPCC Special Report on Renewable Energy Sources and Climate Change Mitigation. Cambridge University Press, 2012.
69

---

[61] C. K. Sønderby, T. Raiko, L. Maalće. S. K. Sønderby, and O. Winther. Ladder variational autoencoders. Advances in neural information processing systems, 29, 2016.
[62] Whitt and S. Gordon. This the economic cost of extreme weather. In World Eco- nomic Forum Annual Meeting, 2023, URL httpp://wwii weforum argg/aenda//2223, 01/ extreme- weather- www.oorreccooic---o cost- weett-wef3///
[63] M. Zamo and P. Naveau. Estimation the continuous ranked probability score with limited information and applications to ensemble weather forecasts. Mathematical Geosciences, 2018.
70

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We outline the scope and contributions of our work in the introduction. See specifically the last paragraph where we list our main contributions.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: Limitations are discussed in section i
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [NA]
Justification: The paper does not present theoretical results.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
We provide enough details about the model architecture and experimental setups to reproduce our results. Apart from the descriptions in the main paper, model details are given in appendix Ceandixxannnd and experimental details in appendices H and Our code is also openly available.
Justification:
5. Open access to data and code
Question: Does the paper provide open access the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: Our code is openly available, together with instructions on how to run it. Along with details given in the appendices this allows for reproducing our results. The ERA5 dataset is openly available. The MEPS data is available to other researchers upon request.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification: See appendices H and
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [No]
The main sources of randomness in our experiments relate to the random seed used for parameter initialization and sampling. A proper statistical analysis would require training and evaluating multiple models. However, the computational cost associated makes it unfeasible to train enough models draw any well founded conclusions.
Justification:
8. Experiments Compute Resources
71

---

For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: See section tables 6 and
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics http://tww..pprr....pppo/opptttt
Answer: [Yes]
Justification:
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: See section and appendix B
11. Safeguards
Question:
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification: We do not judge the shared models have a high risk of misuse.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We explicitly refer to the ERA5 license. The MEPS dataset is used with permission from the the Swedish Meteorological and Hydrological Institute.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We are releasing our code and models with the paper. The code comes with documentation to the degree that other researchers can reproduce and extend our work.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
72


Reviews:
[{"review_id": "PuFl1Ux901", "reviewer": "Reviewer_i1dL", "summary": "The paper proposes Graph-EFM, a method that combines a hierarchical multi-scale graph neural network with a variational objective for probabilistic weather forecasting. The method performs on par with Graphcast on deterministic metrics with the extra benefit of uncertainty estimation.", "questions": "- How important do you think the architecture is to the performance versus the objective function? The proposed architecture has an intuition similar to UNet, i.e., multi-scale features and the lowest layer can be used to parameterize the hidden variable.\n- Diffusion models are considered the best family of models for generative modeling, surpassing GANs and latent variable models for other fields such as computer vision. What is the reason to believe latent variable models are the way to go for probabilistic weather forecasting?\n- Why does the paper compare with Graphcast+SWAG but not the perturbed version of Graphcast and Gencast?\n- How does the performance vary w.r.t. the number of ensemble samples? Given that sampling from a latent variable is fast, have the authors tried using more ensemble members?\n- Is there an explanation why Graphcast is better than Graph-FM, but Graph-EFM is better than Graph-EFM (ms)?\n- Why in LAM, the Graphcast architecture is doing better than the proposed architecture?", "limitations": "See above.", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- The paper is well-written and easy to follow.\n- The paper tackles probabilistic weather forecasting, which is an important problem in the field.\n- The proposed method is intuitive and makes sense. Overall, generative modeling is a potential direction for probabilistic weather forecasting. People have used GANs and diffusion, so a latent variable model is a natural addition to the literature.\n- The performance looks promising, and it is more efficient than existing methods using diffusion.", "weaknesses": "- The authors should replace Table 1 with a line graph figure instead, as it allows comparison across different variables and lead times.\n- Please see my questions below.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer i1dL for useful comments. See our response below: \n\n1. The authors should replace Table 1 with a line graph figure instead, as it allows comparison across different variables and lead times. \n\nGiven the limited space in the main paper we did not find a way to fit line plots for all metrics, multiple variables and both datasets. We found that shrinking the line plots to make this feasible made them hard to read, so instead opted for Table 1 and keeping all line plots in the appendix. If we have the space in the camera ready version we would be happy to bring back some line plots also to the main paper. \n\n2. How important do you think the architecture is to the performance versus the objective function? \n\nBoth of these contributions are indeed important. The objective function, including the CRPS fine-tuning, is important for learning a calibrated probabilistic model. The hierarchical architecture does also fill an important role in Graph-EFM, spreading out the stochasticity from the latent random variable over the forecast region. Since Graph-EFM (ms) uses the same objective function, but not the hierarchical architecture, it can be viewed as an ablation of this part of the model. The poor performance of Graph-EFM (ms) confirms the importance of both of these components in order to learn a useful probabilistic model. \n\n3. What is the reason to believe latent variable models are the way to go for probabilistic weather forecasting? (as opposed to diffusion models) \n\nWhile we believe there is a place also for diffusion models in weather forecasting, one strong reason to explore latent variable models is the difference in sampling speed. Being able to quickly sample large ensembles is a desirable property of probabilistic weather forecasting models. As diffusion models require multiple forward passes through a neural network to produce one sample, this can be a slow process. This issue is aggravated by the common practice of producing forecasts auto-regressively, meaning that each forecast time step must be sampled sequentially. In Graph-EFM this only requires one pass through the network per time step. However, speeding up diffusion model sampling is a very active area of research. We view this as another useful research direction also for weather forecasting, complementary to our approach. \n\n4. Why does the paper compare with Graphcast+SWAG but not the perturbed version of Graphcast and Gencast? \n\nRegarding GenCast comparison see the global author rebuttal. \n\nWe assume here that the perturbed version of GraphCast refers to the one considered in the GenCast paper, where GraphCast is initialized from perturbed initial conditions. We generally consider ensemble forecasting methods that start from multiple (possibly perturbed) initial conditions to be outside of our scope, and thus do not include perturbed GraphCast as a baseline. As outlines in section 3.1 we define our problem setting as modeling the distribution of future weather states conditioned on one specific set of initial conditions. The reason for this is that creating ensembles only based on multiple initial conditions typically limits the ensemble size or requires the use of ad-hoc perturbations without any guarantees of modeling the correct distribution. \n\n5. How does the performance vary w.r.t. the number of ensemble samples? \n\nSee the “global” author rebuttal for an investigation regarding the impact of the ensemble size. \n\n6. Why in LAM, the Graphcast architecture is doing better than the proposed architecture? \n\nThere is some more nuance when comparing GraphCast* and Graph-FM than saying one of these is strictly doing better in any of the experiments. Which out of these two that performs the best depends on what variable and lead time is considered, as can be seen from the full results in Appendix H and I. However, it is true that for longer lead times Graph-FM performs overall worse then GraphCast* in LAM. This is likely related to the structure of the problem in the LAM setting, that features the boundary forcing and other length scales than the global case. We would like to emphasize that we consider the probabilistic Graph-EFM model the main contribution and focus of the paper, and Graph-FM mainly a step in getting there. While further investigation into the possible failure modes of Graph-FM would be interesting, we prefer this to not be the focus of this work. \n\n7. Is there an explanation why Graphcast is better than Graph-FM, but Graph-EFM is better than Graph-EFM (ms)? \n\nAs noted above, the difference between Graphcast* and Graph-FM has some nuance. Still, what should be remembered here is that the graph partly fills a different purpose in the deterministic and probabilistic models. While in the deterministic model it is used to spread information from the input and produce the best possible forecast for the next step, in the probabilistic models the graph has the key role of augmenting the sampled latent variables into a coherent realization of the weather state. In particular, as we associate the latent variable with the mesh nodes, the stochasticity enters in different ways when using the multiscale and hierarchical graphs. Since the multiscale graph does not offer any further dimensionality reduction than mapping from the grid to the mesh, it also means that the latent variable $Z^t$ will have a higher total dimensionality ($|\\mathcal{V}_1| \\times d_z$) in Graph-EFM (ms) than when using the hierarchical graph in Graph-EFM ($|\\mathcal{V}_L| \\times d_z$). For these reasons one should not assume that any performance differences between using the multiscale or hierarchical graphs in deterministic models straightforwardly translates to the probabilistic case.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "G83KG9HWjI", "reviewer": "Reviewer_LjpL", "summary": "The authors propose a graph-based ensemble forecasting model (Graph-EFM) to provide weather prediction with a hirearchical GNN framework. They used a hierarchical mesh graph to handle the challenges of capturing processes unfolding over different spatial scales and modeling the uncertainty in the chaotic system. The Graph-EFM provides a probabilistic weather forecasting with the benefit of capturing forecast uncertainity. The experiment results show the effectiveness and advantages of Graph-EFM compared to other deterministic models.", "questions": "1. The authors could better explain the underlying meaning of meaures, RMSE, CRPS, and SpSkR, for the results of weather forecast. Basically, I would like to ask authors to show how close their Graph-EFM's forecast is to the ground truth weather.", "limitations": "The authors have adequately addressed the limitations.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The hierarchical mesh graph provides a reasonable idea to handle different spatial scales for weather forecasting, which could inspire other researchers to handle problems in different domains.\nThe spatial dependencies are considered and handled within GNN layers.\nUsing ensumble-based model could capture the uncetainty of weather system.", "weaknesses": "In Figure 3, it seems like the selected ensemble members vary a lot, and how close is your forecast to the ground truth. Possibly, explaining a little bit of the underlying meaning of the measures in table 1 & 2 in the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer LjpL for useful comments. See our response below: \n\n1. In Figure 3, it seems like the selected ensemble members vary a lot, and how close is your forecast to the ground truth? \n\nNote that Figure 3 shows the forecasts for 10 days in the future. At such lead times there is indeed a lot of uncertainty in any forecast that can be produced. In Figure 3 we want to show that the ensemble forecast does indeed capture this variability by producing varying ensemble members that each still represent a realistic weather state. The closest forecast (in terms of RMSE) to the ground truth is instead achieved by the ensemble mean. But as expected, due to the large variability at 10 days lead time, this ensemble mean is blurry and not a realistic weather state. We can however still see that some of the large-scale features present in the ground truth are present also in this prediction. \n\n2. The authors could better explain the underlying meaning of measures, RMSE, CRPS, and SpSkR, for the results of weather forecast. \n\nWe agree that the presentation of these metrics in the main part of the paper is quite short. Due to space restrictions we chose to move the full definitions and descriptions of these metrics to appendix D. We encourage anyone unfamiliar with these metrics in the context of weather forecasting to read also this appendix. To give a short, more high-level explanation of these metrics in the weather forecasting context: \n\n**RMSE** measures the average deviation of the forecast from the ground truth data. For probabilistic models, that produce forecasts as samples from a distribution, we compute the RMSE using a forecast representing the mean of this distribution. This is because RMSE is minimized by predicting the mean. \n\n**CRPS** is a probabilistic metric, and measures how well the ground truth value is captured by the distribution specified by the model. In this case this corresponds to how well the ground truth weather is captured by the probability distribution over possible future weather states specified by the model. One useful interpretation of the CRPS is that it measures the difference between the Cumulative Distribution Function (CDF) of the model distribution and the CDF of a Dirac distribution centered at the ground truth value. One can compute a version of CRPS also for deterministic models, in which case the model distribution also becomes a Dirac distribution, but centered at the predicted value. The CRPS then reduces to the Mean Absolute Error (MAE), and this is indeed what we report for the deterministic models in our paper. CRPS is thus useful as a metric that can compare both deterministic and probabilistic forecasts. See [1] for more about CRPS.  \n\n**SpSkR** measures how calibrated the uncertainty of the model is. At SpSkR = 1 the ensemble forecasts are exchangeable with the ground truth, meaning that they exhibit similar levels of variability. This means that we can accurately use the ensemble spread as an indicator for forecast uncertainty. Each forecast will necessarily have some error, and a SpSkR close to 1 indicates that the model expresses the uncertainty about this error correctly. \n\nAll the metrics are computed for each variable and lead time separately, but spatially averaged for all grid points. \n\n[1] T. Gneiting and A. E. Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 2007.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "szSQI3M6qy", "reviewer": "Reviewer_568u", "summary": "This paper introduces a new method for predicting weather using advanced deep learning models. The approach, called Graph-EFM, improves accuracy and better handles uncertainties in weather forecasts. It uses a 1.5 degree version of ERA5 and making weather predictions more reliable and useful for real-world applications.", "questions": "I need to some figures during extreme events, e.g cyclones like Yaku.\n\nCan the model deal with higher resolution data?\n\nOr ERA6 when available?\n\nOr more localised higher resolution data?\n\nDo the results change a lot if it's only trained from 1980 onwards?", "limitations": "Not many, this is a fantastic piece of work", "rating": 10, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The paper's strengths include the innovative use of Graph-EFM for accurate probabilistic weather forecasting, detailed experiments on large datasets, and clear presentation of methods. Graph-EFM significantly enhances uncertainty estimation and forecast reliability adding value to both research and practical weather prediction applications.", "weaknesses": "What happened if 0.25 degree ERA5 is used?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer 568u for useful comments. See our response below: \n\n1. I need to some figures during extreme events, e.g cyclones like Yaku. \n\nWe have now included such a case study for Hurricane Laura, that can be found in the global author rebuttal. \n\n2. About higher resolution data and ERA6 \n\nThere is indeed nothing that technically prevents us from applying these methods to higher resolution data. Some minor adaptations have to be done to the exact graph structures used, but the overall framework is directly applicable. The method should scale well to even higher resolutions, although naturally with an increasing memory requirement. Our choice to focus on 1.5° data is mainly due to the computational needs. \n\nAssuming that ERA6 will follow a similar format as ERA5, although higher resolution, there is no reason to believe our methods would not be applicable also to that dataset. \n\nRegarding higher resolution localized data, note that the MEPS data is exactly an example of this. The 10 km spatial resolution of this data is far higher than the 1.5° ERA5 data for this area (spatial resolution $\\approx$ 167 x 76 km) and even higher than 0.25° ERA5 (spatial resolution $\\approx$ 28 x 13 km). Also for limited area modeling we expect the method to be readily applicable also to data at even higher resolutions. \n\n3. Do the results change a lot if it's only trained from 1980 onwards? \n\nWhile we have not had the time and resources to train such a model now, we would not expect substantial differences in the results from using only data from 1980 onwards. We opted for using as much data as possible, to make sure the model was trained on a set of weather scenarios as diverse as possible. However, in the GraphCast paper it was reported that including data from before 1980 had only minor impact on model performance [1]. Due to variations in climate one might even expect an improvement in performance when considering only more recent data, as the climate in 2020 would be more similar to that of 1980-2017 than that of 1959-1979. Using only more recent years would thus create a smaller shift in the data distribution between training and testing. \n\n[1] R. Lam et al. Learning skillful medium-range global weather forecasting. Science, 2023.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "14MxjFkCoU", "reviewer": "Reviewer_PLfK", "summary": "This work introduces a VAE variant of GraphCast for global medium-range weather forecasting and a VAE variant of a UNet (that is formulated as a GNN) for limited area modeling over Scandinavia. For this, they adapt GraphCast to have a similar hierarchical structure to UNets, and then treat the coarsest hierarchical layer (the bottleneck) as a latent variable representing the mean of isotropic Gaussians. The ensemble predictions from the model are similarly fast as a single deterministic prediction, achieved through batching. Their calibration can be good for some variables (e.g. global t2m for 10 day lead time has a spread/skill ratio of 0.99), while poorer for others (e.g. local wvint has a spread/skill ratio of 0.57 for 24h lead time).", "questions": "Why does the original GraphCast paper not report the visual artifacts that you found for Graph-EFM (ms)? Could it be that your models have simply not been trained sufficiently or that there is a bug in your implementation?\n\nWhy do you use a fixed variance of the latent variable for the global predictions? It would be interesting to see this ablation.", "limitations": "The section on limitations is somewhat short. I believe the key limitation of this work is in the evaluation, i.e. it is unclear to me after reading the work how well the models perform. E.g., are these models robust for longer rollouts or if applied at prediction time further away from training time?\nMoreover, the work does transparently communicate a limitation of their multi-scale Graph-EFM: visual artifacts. But, this is strange to me, as it could mean two things: a) the VAE formulation and multiscale edges simply don't work well together or b) since GraphCast did not observe such issues, the presented Graph-EFM (ms) is simply not trained sufficiently or suffers from a buggy implementation. While a) would be an interesting finding, i believe b) can not be ruled out given the presented results in this paper.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "1. The proposed VAE extension to GraphCast is significantly faster than diffusion-based approaches (like the GenCast model).\n2. The work does not limit itself to just global weather forecasting, but also presents results for limited area modeling, which is the class of models used by many national weather services.\n3. The paper is reasonably well written, keeping a good amount of detail in the main paper, and presenting many additional details in the appendix.", "weaknesses": "Major points:\n1. Questionable baselines: I am unsure if the chosen baselines are very strong, let me name a few reasons for this:\n    - Tab 1 presents performance for GraphCast, e.g. RMSE=387 for z500, 5 day leadtime. However, if I check the headline scores in the WeatherBench 2 (https://sites.research.google/weatherbench/) for GraphCast, i see RMSE=274 for z500, 5 day leadtime, which is significantly higher and beats all models presented in this study.\n    - Both Tab 1 and Tab 2 do not include scores for the conventional weather models. I would expect Tab 1 to include IFS & IFS-ENS scores and Tab 2 to include MEPS scores.\n    - Since this work introduces a probabilistic weather model, i would expect comparison with other recent works on probabilistic weather models, like the ones cited in this paper (e.g. GenCast).\n    - Graph-FM has almost 6x the parameters compared to GraphCast* (Tab 5) - which quite possibly could be the major reason for its improved performance, and not the introduced architectural feature of hierarchical layers.\n2. Overlooked connection to UNets: The Graph-FM that was introduced for the LAM setting looks to me as equivalent to a UNet:\n    - The input data comes on a regular grid with 236 x 268 pixels. Which is subsequently downsampled using 3x3 windows. Processing at each depth level is done with a locally connected layer (in other words: a local convolutional filter). A semantically simpler description of such a model would be a UNet with 3x3 pooling (learned in this case) and 3x3 conv filters at each stage. Possibly, implementing it as a UNet could also be computationally advantageous, making use of the highly optimized kernels for 2d convolutions and pooling operations, instead of GNN layers that rely on scatter sums.\n    - UNets have been previously used for Weather forecasting: e.g. https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018GL080704 & https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2020MS002203\n3. Proposed VAE implementation physically not meaningful? Your VAE is implemented with a latent variable at the coarsest level. This is supposed to capture epistemic uncertainty related to the forward model (and not due to initial state). However, one may argue for atmospheric models most model uncertainty comes from the subgrid-scale parametrizations and not from the coarse-scale representation of atmospheric dynamics. Hence, to me it seems far more intuitive to introduce the stochasticity at the finest level, representing the small scales. I assume you chose the hierarchical description mostly for computational reasons, but given a lack of physical basis, i would at least expect a more thorough investigation of potential errors introduced by this, e.g. is the ensemble variability too smooth?\n4. Missing reference to previously published work? A workshop paper at last years NeurIPS has introduced both the hierarchical GNN and the MEPS dataset https://arxiv.org/abs/2309.17370 , if I am not misstaken. I am not really sure about NeurIPS policy here, but even if this work is a direct extension of the previous work and the previous work is to be considered as non-archival, I still believe you should at least cite the workshop paper.\n\nMinor points:\n1. GraphCast + SWAG: This is a baseline with poor performance, that is somewhat arbitrarily picked from many possible approaches to obtain ensemble predictions from neural networks. I see two options here: Either you keep it, but also introduce many other such baselines, to make clear that you did not cherrypick a particularly weak one. Other approaches  that should not be prohibitively expensive to run could e.g. be MC-Dropout or Laplace Approximation. Or, you simply drop it, as is, it does not  add much to the paper.\n2. Introduction lacks motivation for LAM: This is an ML conference that you are submitting to. It would probably be good to briefly motivate why doing LAM is even necessary (i.e., why can't we just rely on global models instead)?\n3. Extreme Weather evaluation / case study: One key reason for ensemble prediction is capturing the tails, i.e. the extremes. You state in Appendix A that this is out-of-scope for the work. I would argue you are making your life too easy here. Since the presented models are likely not useful unless they display decent performance also for extreme weather, it would be important to evaluate just that. It may be enough for this paper to e.g. study a single extreme event as a case study.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer PLfK for useful comments. See our response below: \n\n1. Scores for GraphCast*\n\nAs we state clearly when introducing this baseline, this is a version of GraphCast trained on the same 1.5° dataset as the other models. It thus has different scores than the original GraphCast model evaluated in Weatherbench. As we emphasize in the paper, we focus here on a fair comparison of models operating at the same resolution and trained on the same data. In the paper we refer to appendix H for comparisons with other state-of-the-art models, trained on higher resolution data.\n\n2. Comparison to conventional models\n\nResults for IFS-ENS are included in appendix H. We realize that this might not be mentioned clearly in the main paper and will clarify. As the paper is mainly about probabilistic modeling, we find it more relevant to include IFS-ENS rather than IFS. \n\nSince the limited area modeling experiment is a pure surrogate modeling task, with data coming from MEPS forecasts, it does not make sense to compare to the MEPS model itself.\n\n3. Comparison to GenCast\n\nSee the global author rebuttal.\n\n4. Graph-FM parameter count\n\nWhen instantiating models of similar size we use the dimensionality $d_z$ of representation vectors as the main scaling factor for model size. The hierarchical architecture in Graph-FM is what enables more flexibility at the same value of $d_z$, resulting also in more parameters. Note that this flexibility does not introduce a significant increase in computational cost, which would be the case if one for example added 6 times as many layers to GraphCast.\n\n5. U-Net connection\n\nWe disagree that this connection is overlooked and refer the reviewer to the related work section where we clearly note this similarity. It should however be emphasized that our GNN model can not be seen as *equivalent* to a U-net. While some hierarchical structure is shared, the GNN layers can not be reduced to simple convolutions, even in the LAM setting. For example, these feature edge representation updates that have no correspondence in convolutional layers.\n\nWe would also like to emphasize that we present a general framework that is applicable to different forecasting regions and graph constructions. The exact similarity to a U-net model will depend on these choices. In the specific LAM case from section 5.2 this connection is stronger (although the models are not equivalent), whereas in the global forecasting setting this connection is more conceptual.\n\n6. Physical meaningfulness of VAE implementation\n\nWe would like to clarify that while the latent variables are associated with the top level of the graph, this does not mean that they only influence the coarsest spatial scales. As the GNN layers of the predictor are applied through the graph hierarchy the randomness in the latent variables can spread to all parts of the graph and introduce variations on different scales. However, explicitly constraining different latent variables to control variation at different spatial scales is an interesting idea for future work.\n\nWe have found that introducing the latent variables closer to the grid points in the model can lead to less spatially coherent forecasts. Graph-EFM (ms) can act as an ablation of this, as the latent variables are there directly associated with the nodes of the single level multiscale graph. We found Graph-EFM (ms) to be harder to train and produce worse forecasts.\n\n7. Reference to workshop paper\n\nBy the NeurIPS policy, workshop papers not published in proceedings are not to be counted as publications. Therefore we did not explicitly cite the mentioned workshop paper. However, if it is deemed advisable to do so we would not mind including a reference to it. Perhaps the AC can help clarify the situation.\n\n8. GraphCast* + SWAG baseline\n\nThis baseline was not picked arbitrarily, but as described it was inspired by the usage of SWAG in Graubner et al. The point of this baseline is to include a comparison to multi-model ensembles, which we see can produce too little variability. While there are indeed other ways to achieve this, we wanted to focus on a method that has been proposed in the ML weather prediction literature, rather than choosing a method arbitrarily.\n\n9. Motivation for LAM\n\nWe currently give some motivation for this in section 5.2, but we do agree that this can be expanded on already in the introduction of the paper. We will change this.\n\n10. Extreme weather case study\n\nWe have now included such a case study for Hurricane Laura, that can be found in the global author rebuttal. \n\n11. Fixed variance of latent variable for global predictions \n\nThere might be some misunderstanding here, as the variance of the latent variable (Eq. 4) is fixed for both the global and LAM models. However, the variance $\\sigma^2_{\\alpha, j}$ of the likelihood term (Eq. 6) is fixed in the global setting while output by the model in LAM. See Appendix C for more discussion about this. \n\n12. Robustness for longer rollouts \n\nNote that all models are unrolled for much fewer time steps during training than at inference time. In the global case for example, Graph-EFM is unrolled 40 steps for evaluation at 10 days lead time, but only trained by unrolling up to 12 steps. As we focus here on weather forecasting rather than climate modeling, the stability of unrolling the models for months or years is of limited interest in this work.  \n\n13. Visual artefacts in Graph-EFM (ms) \n\nWe would like to emphasize that when considering the same setting as in the GraphCast paper (global, deterministic modeling) we do not see any artefacts. The artefacts in Graph-EFM (ms) are mainly observed in the LAM setting (e.g. Figure 6.d) and for the global probabilistic model, none of which were considered in the GraphCast paper. Note however that the final global Graph-EFM (ms) model does not have these artefacts, as the artefacts mentioned as limitations only appear when training with a poor choice of $\\lambda_{\\text{CRPS}}$.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank all reviewers for valuable comments and questions that we are sure will improve the overall quality of our paper. We have responded to the points raised by each reviewer separately, but also include this general rebuttal with a few points that we think could be relevant to all. 

### Extreme Weather Case Study

Multiple reviewers requested a case study on ensemble forecasting for extreme weather events. We agree that this would be a valuable addition to the paper, as this is one important motivation for the work. 

We include here a case study for Hurricane Laura, that will be added to the paper. In August of 2020 Laura developed in the Atlantic Ocean and reached Hurricane levels in the Gulf of Mexico, eventually making landfall in Louisiana and causing major damages [1]. We study here forecasts for 2020-08-27T12 UTC, which is 6 hours after the Hurricane hit land. All forecasts are for this exact time point, but initialized a varying number of days before. We here run 50 member ensemble forecasts using Graph-EFM and deterministic forecasts using the GraphCast* and Graph-FM models. Figure 1 in the attached pdf shows 10 m wind speeds in ERA5 and the forecasts. For Graph-EFM we plot both forecasts from randomly sampled ensemble members and a cherry-picked best member that was deemed to most closely match ERA5. Note that the 1.5° resolution that we work with makes determining similarity or exact positions somewhat challenging. 

We see in Figure 1 that at 7 days lead time there is great uncertainty, and the deterministic models do not show the hurricane at all. In the ensemble forecast from Graph-EFM there exists however already members indicating the possibility of a Hurricane making landfall a week ahead (see for example the cherry-picked “Best member”). While the ensemble includes many possible scenarios, a total of 7 members show the development of a hurricane in the area. Having information of such possible scenarios a long time ahead allows for planning and readying disaster response efforts that might be needed. Note that discovering these scenarios is only possible through an ensemble forecast, as the deterministic models do not indicate such an event. At 5 days lead time all models are indicating the development of a hurricane. While the deterministic models do a good job here, they are indicating a landfall location slightly too far eastward. The ensemble members from Graph-EFM however show a range of different positions for the hurricane, indicating the uncertainty in the landfall location. At 3 days ahead 42 out of 50 ensemble members show the hurricane making landfall, but still with some uncertainty about the exact location. At 1 day lead time all models give an accurate forecast of the position of the hurricane. Apart from position, it is also interesting to consider how the models capture the intensity in terms of wind speeds. At 1 day ahead the deterministic models are somewhat underestimating the wind speed. The Graph-EFM ensemble shows a range of possible values, indicating the uncertainty in the exact wind intensity. Overall this study exemplifies how a large ensemble forecast from a machine learning model can be used to discover possible extreme weather scenarios at long lead times and uncover the uncertainties associated with them. 

[1] Hurricane Laura 2020. National Weather Service. https://www.weather.gov/lch/2020Laura  

### Impact of Ensemble Size on Results

As pointed out by reviewer i1dL it would be interesting to know how the performance of the model, in terms of metric values, varies when sampling different number of ensemble members. To investigate this we ran the evaluation of Graph-EFM with 5-80 members for the global data and 5-100 for the limited area model.  

Results for a selection of variables and metrics are shown in Figure 2 in the attached pdf. As expected the RMSE of the ensemble mean decreases when sampling more members. However, already when sampling 20 or 25 members the results are fairly close to the full ensemble. For SpSkR the differences are even smaller. As the CRPS is a property of the distribution of the model forecast its true value does not depend on the number of samples drawn. In practice we compute CRPS using an unbiased estimator, and the variance of this estimator decreases with ensemble size. When averaged spatially and over the whole test set we do however not see any difference in CRPS for different ensemble sizes. All these trends hold consistently for all variables in both the global and limited area datasets. We intend to add the results of this investigation to the paper. 

Given that any improvements to metrics saturate, we would not expect the results to meaningfully change from sampling more than 80/100 members. It should however be noted that the motivation for sampling very large ensembles is mainly not to improve on metrics such as these. More important motivations for large ensembles include estimating probabilities of rare events or studying different possible scenarios of extreme weather. 

### Comparison to GenCast

The reason for not comparing against GenCast is the lack of any openly available code, pre-trained models or produced forecasts (at the time of writing). Such a comparison would require reimplementing and retraining a GenCast model from scratch, which we deem outside the scope of this work. Note that doing this for GraphCast was much more feasible, as there is code openly available and the model is more similar to ours.

Meta-review: The authors propose an autoencoder version of GraphCast and UNet for probabilistic global medium-range weather forecasting.

- **Reviews**: While the paper received (very) high scores from some reviewers, some reviewers were initially more critical.
- **Rebuttal**: The authors successfully addressed the concerns of the more critical reviewers by discussing computational challenges and justifying their approach through experiments.
- **Decision**: I recommend acceptance.
Decision: Accept (spotlight)