Postprocessing of Ensemble Weather Forecasts Using Permutation-invariant
  Neural Networks
Statistical postprocessing is used to translate ensembles of raw numerical
weather forecasts into reliable probabilistic forecast distributions. In this
study, we examine the use of permutation-invariant neural networks for this
task. In contrast to previous approaches, which often operate on ensemble
summary statistics and dismiss details of the ensemble distribution, we propose
networks that treat forecast ensembles as a set of unordered member forecasts
and learn link functions that are by design invariant to permutations of the
member ordering. We evaluate the quality of the obtained forecast distributions
in terms of calibration and sharpness and compare the models against classical
and neural network-based benchmark methods. In case studies addressing the
postprocessing of surface temperature and wind gust forecasts, we demonstrate
state-of-the-art prediction quality. To deepen the understanding of the learned
inference process, we further propose a permutation-based importance analysis
for ensemble-valued predictors, which highlights specific aspects of the
ensemble forecast that are considered important by the trained postprocessing
models. Our results suggest that most of the relevant information is contained
in a few ensemble-internal degrees of freedom, which may impact the design of
future ensemble forecasting and postprocessing systems.
Kevin Höhlein, Benedikt Schulz, Rüdiger Westermann, Sebastian Lerch
# Postprocessing of Ensemble Weather Forecasts

###### Abstract

Statistical postprocessing is used to translate ensembles of raw numerical weather forecasts into reliable probabilistic forecast distributions. In this study, we examine the use of permutation-invariant neural networks for this task. In contrast to previous approaches, which often operate on ensemble summary statistics and dismiss details of the ensemble distribution, we propose networks which treat forecast ensembles as a set of unordered member forecasts and learn link functions that are by design invariant to permutations of the member ordering. We evaluate the quality of the obtained forecast distributions in terms of calibration and sharpness, and compare the models against classical and neural network-based benchmark methods. In case studies addressing the postprocessing of surface temperature and wind gust forecasts, we demonstrate state-of-the-art prediction quality. To deepen the understanding of the learned inference process, we further propose a permutation-based importance analysis for ensemble-valued predictors, which highlights specific aspects of the ensemble forecast that are considered important by the trained postprocessing models. Our results suggest that most of the relevant information is contained in few ensemble-internal degrees of freedom, which may impact the design of future ensemble forecasting and postprocessing systems.

Key words and phrases:

## 1 Introduction

Operational weather forecasting relies on numerical weather prediction (NWP) models. Since such models are subject to multiple sources of uncertainty, such as uncertainty in the initial conditions or model parameterizations, a quantification of the forecast uncertainty is indispensable. To achieve this, NWP models generate a set of deterministic forecasts, so-called ensemble forecasts, based on different initial conditions and variations of the underlying physical models. Since these forecasts are subject to systematic errors such as biases and dispersion errors, statistical postprocessing is used to enhance their reliability (see, e.g., Vannitsem et al.2018).

Recently, machine learning (ML) approaches for statistical postprocessing have shown superior performance over classical methods. For instance, Rasp and Lerch (2018) propose a distribution regression network (DRN) which predicts the parameters of a temperature forecast distribution from a suitable family of parametric distributions. In subsequent work, Schulz and Lerch (2022b) found that shallow multi-layer perceptrons (MLPs) with forecast distributions of different flexibility achieve state-of-the-art results in postprocessing wind gust ensemble forecasts.

An ensemble forecast consists of multiple separate member forecasts, which are generated by repeatedly running NWP simulations with different model parameterizations and initial conditions. Typically, the configurations of different runs are sampled randomly from an underlying distribution of plausible simulation conditions, obtained, e.g., from uncertainty-aware data assimilation. The member forecasts can then be seen as identically distributed and interchangeable random samples from a distribution of possible future weather states. In this setting, statistical postprocessing of ensemble forecasts can be phrased as a prediction task on unordered predictor vectors and requires solutions that are tuned to match the predictor format. Specifically, member interchangeability demands that the predictions of a well-designed postprocessing system should not be affected by permutations, i.e. shuffling, of the ensemble members. Systems that satisfy this requirement are called permutation invariant. Established postprocessing methods rely on basic summary statistics of the raw ensemble forecast to inform the estimation of the postprocessed distribution and are thus permutation invariant by design. However, especially in large ensembles, the details of the distribution may carry valuable information for postprocessing, and a more elaborate treatment of the inner structure of the raw forecast ensembles may be advisable.

Alleviating such restrictions, Bremnes (2020) employs MLPs for postprocessing of wind speed forecasts, which receive information about the full state of a univariate ensemble forecast. Yet, the size of the ensemble, i.e. the number of member forecasts, acts as a multiplier to the dimension of the predictors in the proposed models, resulting in an increased model complexity and a higher tendency to overfitting. Only recently, studies have started to explore how more dedicated model architectures can help to improve postprocessing (Mlakar et al.2023; Ben-Bouallegue et al.2023), and ML provides a variety of further approaches
to enforcing permutation invariance in data-driven learning (e.g., Ravanbakhsh et al., 2016; Vaswani et al., 2017; Zaheer et al., 2017; Lee et al., 2019; Sannai et al., 2019; Zhang et al., 2019). The increasing adoption of permutation-invariant statistical models in postprocessing thus raises the question of how capable different model architectures are in extracting information from the ensemble forecasts and how much value is added by considering ensemble-valued predictors instead of summary statistics.

### Contribution

In this study, we investigate the capabilities of different permutation-invariant NN architectures for univariate postprocessing of station predictions. We evaluate the proposed models on two exemplary station-wise postprocessing tasks with different characteristics. The ensemble-based network models are compared to classical methods and basic NNs which operate only on ensemble summary statistics but are trained under identical predictor conditions otherwise. We further assess how much of the predictive information is carried within the details of the ensemble distribution, and how much of the model skill arises from other factors. To shed light on the sources of model skill, we propose an ensemble-oriented importance analysis and study the effect of ensemble-internal degrees of freedom using conditional feature permutation.

## 2 Related work

### Statistical postprocessing of ensemble forecasts

Two of the first methods for statistical postprocessing of ensemble forecasts are ensemble model output statistics (EMOS; Gneiting et al., 2005) and Bayesian model averaging (BMA; Raftery et al., 2005). While EMOS performs a distributional regression based on a suitable family of parametric distributions and summary statistics of the ensemble, BMA generates a mixture distribution based on the individual ensemble members. Due to its simplicity, EMOS has been applied to a wide range of weather variables including temperature (Gneiting et al., 2005), wind gusts (Pantillon et al., 2018), precipitation (Scheuerer, 2014) and solar radiation (Schulz et al., 2021). Following the simple statistical approaches, ML approaches such as quantile regression forests (Taillardat et al., 2016) or a gradient boosting extension of EMOS (Messner et al., 2017) have been introduced. First NN-based approaches are the DRN approach (Rasp and Lerch, 2018) as an extension of the EMOS framework, and the Bernstein quantile network (BQN; Bremnes, 2020) that provides a more flexible forecast distribution. In Schulz and Lerch (2022b), NN-based approaches were adapted towards the prediction of wind gusts and outperformed classical methods. Recently, research has shifted towards the use of more sophisticated network architectures. Examples include convolutional NNs that incorporate spatial NWP output fields (Scheuerer et al., 2020; Gronquist et al., 2021; Veldkamp et al., 2021; Horat and Lerch, 2023), and generative models to produce multivariate forecast distributions (Dai and Hemri, 2021; Chen et al., 2022). Only recently, Mlakar et al. (2023) propose NN models that explicitly address the ensemble structure of the inputs by employing a dynamic attention mechanism. This model performs best in the benchmark study of Demaeyer et al. (2023). In orthogonal work, Ben-Bouallegue et al. (2023) postprocess each member individually with hierarchical ensemble transformers, and Orlova et al. (2022) found that exploiting the ensemble structure enhances the predictive performance in the context of sub-seasonal forecasting.

For a general review of statistical postprocessing of weather forecasts, we refer to Vannitsem et al. (2018), a review of recent developments and challenges can be found in Vannitsem et al. (2021); Haupt et al. (2021).

### Neural network architectures for regression on set-structured data

Vinyals et al. (2015) compare a sequential method for processing unordered predictors against a permutation-invariant alternative model and demonstrate that the lack of built-in permutation invariance may substantially affect prediction quality. Ravanbakhsh et al. (2016) introduce a permutation-equivariant layer for NNs that operate on set-structured data and combine these to design permutation-invariant networks. Similar layers were later used by Zaheer et al. (2017), who propose the framework _DeepSets_, which is discussed in more detail in section 4a. Murphy et al. (2018) propose Janossy pooling, which obtains a permutation-invariant mapping as the average of a permutation-sensitive function applied to all possible reorderings of the set. The authors propose methods to lessen the computational burden of the method, but the resulting mappings are subject to constraints or achieve permutation invariance only approximately so we do not consider the approach in our comparison. Limitations of representing functions on sets have been discussed by Wagstaff et al. (2019).

Lyle et al. (2020) demonstrate further that algorithmically enforced permutation invariance is favorable in a variety of tasks compared to alternative approaches, such as data augmentation. Pooling-type network architectures were introduced by Edwards and Storkey (2016) and investigated in more detail by Zaheer et al. (2017) and Sannai et al. (2019), who prove that pooling architectures with additive pooling are universal approximators of functions on sets. Yet, Soelch et al. (2019) highlight that the use of more expressive pooling functions may enhance model performance. A different approach is considered
by Lee et al. (2019), who use (multi-head) attention functions (Vaswani et al., 2017) for permutation-invariant inference on set-valued data. Attention-based models, also known as transformers, have proven powerful in a variety of computer vision tasks (e.g., Khan et al., 2022) and have more recently also found meteorological applications, such as data-driven weather forecasting (e.g. Pathak et al., 2022) and postprocessing (e.g. Finn, 2021; Ben-Bouallegue et al., 2023).

### Machine learning explainability and feature importance

ML explainability has attracted substantial interest throughout the last decade (Guidotti et al., 2018; Linardatos et al., 2020; Sahakyan et al., 2021; Burkart and Huber, 2021; Zhang et al., 2021). Model explanation aims at understanding black-box algorithms by assessing the general logic of the algorithm, often using sample-based explanations. Many variants exist (e.g., Bach et al., 2015; Ribeiro et al., 2016; Shrikumar et al., 2017; Lundberg and Lee, 2017) and are increasingly adopted in the earth-system sciences (e.g., Labe and Barnes, 2021; Farokhmanesh et al., 2023). Such techniques explain model predictions by assigning sample-specific relevance or attribution scores to the model inputs (i.e., the predictors), and deriving the effective strength of certain predictors. While attribution-based approaches are well suited for an in-depth investigation of the model inference they usually come at high computational cost and provide information that is too fine-grained for a comparative evaluation of algorithms.

For these higher-level tasks, averaged importance scores of certain predictors, as obtained, e.g. through feature permutation importance (FPI; Breiman, 2001), are more informative. FPI is commonly implemented as a post-training step, in which relevance scores are assigned to the predictors based on the accuracy loss after permuting the predictor values within the test dataset. In this work, we propose a conditional permutation importance measure for ensemble-valued predictors, which allows attributing importance values to different aspects of the ensemble-internal variability. Conditional perturbation measures have been considered in earlier works (e.g., Strobl et al., 2008; Molnar et al., 2023), yet there the importance of specific predictors is evaluated in the context of the remaining predictors, whereas our approach addresses specifically the ensemble structure of the raw forecasts encountered in postprocessing.

## 3 Benchmark methods and forecast distributions

### Assessing predictive performance

We evaluate probabilistic forecasts based on the paradigm of Gneiting et al. (2007), i.e., a forecast should maximize sharpness subject to calibration. Both sharpness and calibration can be assessed quantitatively using proper scoring rules (Gneiting and Raftery, 2007). A popular choice is the continuous ranked probability score (CRPS; Matheson and Winkler, 1976)

\[\text{CRPS}(F,y)=\int_{-\infty}^{\infty}\left(F(z)-\mathbb{1}\left\{y\leq z \right\}\right)^{2}dz,\]

wherein \(y\in\mathbb{R}\) is the observed value, \(F\) the cumulative distribution function (CDF) of the forecast distribution, and \(\mathbb{1}\) the indicator function. The CRPS can be computed analytically for a wide range of distributions including the truncated logistic distribution and probabilistic forecasts in ensemble form (Jordan et al., 2019).

In addition to the CRPS, we assess calibration based on the empirical coverage of prediction intervals (PIs) derived from the forecast distribution, and sharpness on the corresponding length. Under the assumption of calibration, the observed coverage of a PI should match the nominal level, and a forecast is sharper the smaller the length of the PI. In line with Schulz and Lerch (2022), we choose the PI level based on the size of the underlying ensemble. For an ensemble of size \(M\), this gives rise to a PI with nominal level \((M-1)/(M+1)\).

Further, we qualitatively assess calibration based on (unified) probability integral transform (PIT) histograms (Gneiting and Katzfuss, 2014; Vogel et al., 2018). While a flat histogram indicates that the forecasts are calibrated, systematic deviations indicate miscalibration. For more details on the evaluation of probabilistic forecasts, we refer to Gneiting and Katzfuss (2014).

### Distributional regression with a parametric forecast distribution

In this study, we consider postprocessing of the ensemble forecast distribution of a real-valued random variable \(Y\) as a distributional regression task on ensemble-structured predictors. We consider the case of station-wise forecasts, which are given as prediction vectors \(\mathbf{x}\in\mathcal{P}\subseteq\mathbb{R}^{P}\), each comprising the predictions of \(p\) scalar-valued meteorological variables, such as surface temperature or 10-m wind speed at a station site. For \(M\in\mathbb{N}\), an \(M\)-member ensemble forecast \(X\in\left\{\mathcal{P}\right\}_{M}\), with \(\left[\mathcal{P}\right]_{M}:=\left\{\left\{\mathbf{x}_{1},...,\mathbf{x}_{M}\right\} :\mathbf{x}_{m}\in\mathcal{P},m=1,...,M\right\}\), is a set of separate prediction vectors, which all concern the same forecasting task. Within the (parametric) distributional regression framework, the parameter vector \(\mathbf{\theta}\in\Theta\subseteq\mathbb{R}^{D}\), \(D\in\mathbb{N}\), of a parametric distribution \(\mathcal{F}_{\theta}\) is linked to the predictors via a function that is estimated by minimizing a proper scoring rule. The underlying model can be written as

\[Y\mid X\sim\mathcal{F}_{\mathbf{\theta}},\quad\mathbf{\theta}=g(X)\in\Theta, \tag{1}\]

where \(g:\left[\mathcal{P}\right]_{M}\to\Theta\) is called the link function.

For EMOS, the link function is typically a generalized affine-linear function of ensemble summary statistics \(s:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}\), where \(D\) is the dimension of the model, \(\mathbf{\theta}\) is the mean of the prediction vector, and \(\mathbf{\theta}\) is the mean of the prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{2}\]

where \(\mathbf{\theta}\) is the mean of the prediction vector, and \(\mathbf{\theta}\) is the mean of the prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{3}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{4}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{5}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{6}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{7}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{8}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{9}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{10}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{11}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{12}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{13}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{14}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{15}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{16}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{17}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{18}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{19}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{20}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{21}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{22}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{23}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as

\[\mathbf{\theta}=\mathbf{\theta}+\mathbf{\theta}, \tag{24}\]

where \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector, and \(\mathbf{\theta}\) is the mean of prediction vector. The model is then defined as
\(\left[\mathcal{P}\right]_{M}\rightarrow\mathbb{R}\), such as ensemble mean or standard deviation. I.e., given \(F\in\mathbb{N}\) summary features \(s_{f}\), for \(f=1,...,F\), the link function reads

\[g(X)=\rho(\Gamma\mathbf{s}_{X}+\mathbf{\gamma})\,, \tag{2}\]

wherein \(\mathbf{s}_{X}:=(s_{1}(X),...,s_{F}(X))\in\mathbb{R}^{F}\), and \(\Gamma\in\mathbb{R}^{D\times F}\) and \(\mathbf{\gamma}\in\mathbb{R}^{D}\) denote the parameters of the optimized affine-linear transform. The function \(\rho:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}\) indicates a combination of element-wise activation functions, such as \(\exp(\cdot)\) or the identity.

DRN (Rasp and Lerch, 2018; Schulz and Lerch, 2022) overcomes the need to pre-define the detailed structure of \(g\) by admitting the data-driven estimation of arbitrary link functions using NNs, i.e.,

\[g(X)=\phi_{\mathbf{\beta}}\left(\mathbf{s}_{X}\right) \tag{3}\]

with \(\phi_{\mathbf{\beta}}\) denoting an MLP with parameters \(\mathbf{\beta}\). The forecast distribution as well as the underlying proper scoring rule used for optimization are two implementation choices. We note that Schulz and Lerch (2022) employ a learned station embedding to reuse the same set of model parameters for multiple weather stations. Since this design choice does not affect the treatment of the forecast ensemble, we subsume this embedding within the model parameters \(\mathbf{\beta}\) for brevity.

### Flexible distribution estimator

Distributional regression methods based on a parametric forecast distribution are robust but lack flexibility as they are bound to the parametric distribution family of choice. Typical choices of forecast distributions include the normal (Gneiting et al., 2005; Rasp and Lerch, 2018), logistic (Schulz and Lerch, 2022) or generalized extreme value distribution (Lerch and Thorarinsdottir, 2013; Scheuerer, 2014). They all lack the ability to express multi-modalities that are required when different weather patterns occur. Hence, methods that do not rely on parametric assumptions have been proposed in the postprocessing literature. Examples are the direct adjustment of the ensemble members (van Schaeybroeck and Vannitsem, 2015) or quantile regression forests (Taillardat et al., 2016).

To incorporate the full ensemble structure, we consider flexible forecast distributions that are able to represent the distributional structure of the ensemble forecast in more detail. BQN (Bremnes, 2020) models the forecast distribution as a quantile function

\[Q(p|\mathbf{\alpha})=\sum_{\nu=0}^{d}\alpha_{\nu}\,B_{\nu d}(p)\,,\quad p\in[0,1]\,, \tag{4}\]

which is a linear combination of Bernstein (basis-)polynomials \(B_{\nu d}(p)\) of degree \(d\in\mathbb{N}\), \(\nu=0,\ldots,d\), with mixing coefficients \(\mathbf{\alpha}=(\alpha_{0},...,\alpha_{d})\) such that \(\alpha_{0}\leq\alpha_{1}\leq...\leq\alpha_{d}\). The inference network is designed to output parameters \(\mathbf{\theta}\) that parameterize the mixing coefficients, i.e., \(\mathbf{\alpha}=\mathbf{\alpha}(\mathbf{\theta})\). In contrast to DRN, this formulation offers increased flexibility for modeling multi-modality, while requiring hard upper and lower bounds for the values of the forecast variable. For BQN models, the optimization is guided by an average of quantile scores (Koenker and Bassett, 1978; Gneiting and Raftery, 2007), which can be seen as a discrete approximation of the CRPS (Gneiting and Ranjan, 2011).

Note that our BQN implementation differs from that of Bremnes (2020) and Schulz and Lerch (2022) in how the ensemble predictors are used. Bremnes (2020) specifies the link function as

\[g(X)=\phi_{\mathbf{\beta}}\left(\text{sort}(X)\right), \tag{5}\]

wherein \(\phi_{\mathbf{\beta}}\) is a NN and sort\((\cdot)\) indicates a sorting operation of the ensemble members wrt. a pre-determined reference quantity. While Bremnes (2020) restrict their models to univariate ensemble predictors, Schulz and Lerch (2022) use ensemble-valued predictors only for the predictor of the target variable and ensemble means for additional auxiliary predictor variables. The sorting imposes a fixed ordering of the members and thus ensures permutation invariance of the model predictions. Yet, this increases the number of weights in the initial layer of the network and restricts the trained model to ensembles of fixed size. Hence, we incorporate the ensemble as described in Eq. (3), analogous to DRN using ensemble summary statistics. A comparison of both variants is conducted in the supplementary materials, demonstrating the equivalence of both approaches.

### Usage of auxiliary predictors

In addition to the predictions of the postprocessing target variable, most algorithms use auxiliary information to improve the prediction performance (see Table 1. We distinguish between ensemble-valued and scalar-valued predictors, where ensemble-valued predictors vary between different members and scalar-valued predictors do not. In the ensemble-valued case, we further differentiate the prediction of the postprocessed quantity, termed the primary prediction, from auxiliary predictions of other meteorological variables. For either of these, postprocessing models can have access to the full set of ensemble values or only to summary statistics. Scalar-valued predictors refer to contextual information, such as station-specific coordinates and orography details (cf. Table 1: station pred.), as well as to temporal information, such as the day of the year. We consider only models that are trained on predictions for specific initialization and lead times, such that information about the diurnal cycle is not required. While most approaches include the scalar predictors explicitly as features in the regression process, EMOS takes advantage of categorical location and time information implicitly by fitting
independent models for each station and month (Schulz and Lerch, 2022b). Notably, the permutation-invariant models (cf. Table 1: perm.-inv.) considered in this study have access to the richest predictor pool.

## 4 Permutation-invariant neural network architectures

From the variety of permutation-invariant model architectures, we select two representative approaches, _set pooling architectures_ and _set transformers_, which we adapt for distributional regression. Compared with the benchmark methods of section 3, the proposed networks replace the summary-based ensemble processing while the parameterization of the forecast distributions remains unchanged. A schematic comparison of both permutation-invariant architectures is shown in Fig. 1.

### Set pooling architectures

Set pooling architectures (Zaheer et al., 2017), also known as _DeepSets_, achieve permutation invariance via extraction and permutation-invariant summarization of learned latent features. The features are obtained by applying an encoder network to all ensemble members separately. The resulting link function can be expressed as

\[g(X) =\phi_{\mathbf{\beta}}^{\text{Dec}}(\text{pool}(\tilde{X})), \tag{6}\] \[\text{where}\ \tilde{X} =\{\phi_{\mathbf{\beta}}^{\text{Enc}}(\mathbf{x}):\mathbf{x}\in X\}\}. \tag{7}\]

Therein, pool is a permutation-invariant pooling function, and \(\phi_{\mathbf{\beta}}^{\text{Enc}}\) and \(\phi_{\mathbf{\beta}}^{\text{Dec}}\) are trainable MLPs, acting as encoder and decoder, respectively. We will thus use the names _set pooling_ and _encoder-decoder_ (ED) architecture synonymously.

We consider different variants of ensemble summarization based on average and extremum pooling, as well as adaptive pooling functions based on an attention mechanism (Lee et al., 2019; Soelch et al., 2019), discussed below. Overall, we find that the pooling mechanism is of minor importance. Detailed comparisons are thus deferred to the supplementary materials. In all subsequent experiments, we use attention-based pooling.

### Set Transformer

Set transformers (Lee et al., 2019) are NNs, which model interactions between set members via self-attention. _Attention_ is a form of nonlinear activation function, in which the relevance of the inputs is determined via a matching of input-specific key and query vectors. _Multi-head attention_ allows the model to attend to multiple key patterns in parallel (Vaswani et al., 2017). Lee et al. (2019) combine multi-head attention with member-wise NN to build a permutation-invariant set-attention block, from which a set transformer is built by stacking multiple instances. Set transformers apply straight-forwardly to ensemble data and can exploit all aspects of the available ensemble dataset by allowing for information exchange between ensemble members early in the inference process (cf. (1)). We construct a set transformer by using three set-attention blocks with 8 attention heads (Vaswani et al., 2017; Lee et al., 2019). Each block comprises a separate MLP with two hidden layers. Additionally, the first set-attention block is preceded by a linear layer to align the channel number of the ensemble input with the hidden dimension of the set-attention blocks. To construct vector-valued predictions from set-valued inputs, Lee et al. (2019) propose attention-based pooling, in which the output query vectors are implemented as learnable parameters. After pooling, the final prediction is obtained by applying another two-layer MLP.

## 5 Data

We evaluate the performance of the proposed models in two postprocessing tasks using the datasets described in Table 2.

### Wind gust prediction in Germany

In the first case study, we employ our methods for station-wise postprocessing of wind gust forecasts using a dataset that has previously been used in Pantillon et al. (2018) and Schulz and Lerch (2022b). The ensemble forecasts are based on the COSMO-DE (Baldauf et al., 2011) ensemble prediction system (EPS) and consist of 20 members forecasts, simulated with a horizontal resolution of 2.8 km. The forecasts are initialized at 00 UTC, and we consider the lead times 6, 12 and 18h. Other than wind gusts, the

Figure 1: Set pooling architecture (top), consisting of encoder and decoder MLPs, and set transformer (bottom), featuring attention blocks and intermediate MLPs with residual connections. While the encoder-decoder architecture admits interactions between members only inside the pooling step, the set transformer admits information transfer between the members in each attention step.
dataset comprises ensemble forecasts of several meteorological variables, such as temperature, pressure, precipitation and radiation. The predictions are verified against observations measured at 175 stations of the German weather service (Deutscher Wetterdienst; DWD). Forecasts for the individual weather stations are obtained from the closest grid point. The time period of the forecast and observation data starts on 9 December 2010 and ends at 31 December 2016. The models use the data from 2010-2015 for model estimation, using 2010-2014 as training and 2015 as validation period. The forecasts are then verified in 2016. As in Schulz and Lerch (2022b), each lead time is processed separately.

As detailed in Schulz and Lerch (2022b), a minor caveat is caused by a non-trivial substructure of the forecast ensembles. The 20-member ensembles constitute a conglomerate of four sub-ensembles, which are generated with slightly different model configurations. While this formally violates the assumption of statistical interchangeability of the members, the sub-ensembles are sufficiently similar to justify the application of permutation-invariant models.

For the benchmark methods EMOS and DRN, we use the exact same forecasts as in Schulz and Lerch (2022b), both estimating the parameters of a truncated logistic distribution by minimizing the CRPS, see their section 3 for details. BQN is adapted as described in section 3 and Table 1.

### Temperature forecasts from the EUPPBench dataset

In a second example, we postprocess ensemble forecasts of surface temperature using a subset of the EUPPBench postprocessing benchmark dataset (Demaeyer et al., 2023). EUPPBench provides paired forecast and observation data from two sets of samples. The first part consists of 20 years of reforecast data (1997 - 2016) from the Integrated Forecasting System (IFS) of the ECMWF with 11 ensemble members. Mimicking typical operational approaches, the reforecast dataset is used as training data, complemented by additional two years (2017 and 2018) of 51-member forecasts as test data. EUPPBench comprises sample data from multiple European countries - Austria, Belgium, France, Germany and the Netherlands - which are publicly accessible via the CliMetLab API (ECMWF, 2013). Additional data for Switzerland can be requested from the Swiss weather service, but is not used in this study. EUPPBench constitutes a comprehensive dataset of samples over a long time period. In contrast to the wind gust forecasts, the EUPPBench ensemble members are exchangeable, so that permutation-invariant model architectures are optimally suited.

\begin{table}
\begin{tabular}{l c c c} \hline \hline Predictors & Ensemble-valued & Scalar-valued \\ \hline Method & Primary prediction & Auxiliary predictions & Spatial & Temporal \\ \hline EMOS (Schulz and Lerch, 2022b, ours) & mean + std. dev. & – & different models per station and month \\ \hline BQN (Bremnes, 2020) & ensemble (sorted) & – & station embed. & – \\ BQN (Schulz and Lerch, 2022b) & ensemble (sorted) & mean & station pred. + embed. & day of year \\ BQN (ours) & mean + std. dev. & mean & station pred. + embed. & day of year \\ DRN (Schulz and Lerch, 2022b, ours) & mean + std. dev. & mean & station pred. + embed. & day of year \\ \hline Perm.-inv. DRN + BQN (ours) & ensemble (perm.-inv.) & ensemble (perm.-inv.) & station pred. + embed. & day of year \\ \hline \hline \end{tabular}
\end{table}
Table 1: Predictor utilization by postprocessing methods. Methods used in this study are indicated by _ours_.

\begin{table}
\begin{tabular}{l r r} \hline \hline Dataset & Wind gust forecasts & EUPPbench (re)forecasts \\ \hline Underlying NWP model & COSMO-DE-EPS & ECMWF-IFS \\ Initialization time & 00 UTC & 00 UTC \\ Ensemble size \(M\) & 20 & Reforecasts: 51 \\  & & Forecasts: 11 \\ Predicted ensemble forecast quantities \(p\) & 61 & 28 \\ Region & Germany & Central Europe \\ Stations & 175 & 117 \\ Lead times considered in h & 6, 12, 18 & 24, 72, 120 \\ Training samples & 315,000 & 374,000 \\ Test samples & 63,000 & Reforecasts: 97,000 \\  & & Forecasts: 85,000 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Overview of the data used in the postprocessing applications described in section 5.
Deviating from the EUPPBench convention, models are tested on the 51-member forecasts and the last 4 years of the reforecast dataset are considered as an independent test set of 11-member forecast samples. This allows us to assess the generalization capabilities of the ensemble-based postprocessing models on data equivalent to the training data, as well as on data with larger ensemble size. Furthermore, we use the full set of available surface- and pressure-level predictor variables, whereas the original EUPPBench task is restricted to using only surface temperature data. While this design choice hinders the direct comparison of the evaluation metrics in this paper with original EUPPBench models, it enables a more comprehensive assessment of the relative benefits of using summary-based vs. ensemble predictors. From the pool of available forecast lead times, we select 24h, 72h and 120h for a closer analysis.

Unlike previous postprocessing applications for temperature (e.g., Gneiting et al. 2005; Rasp and Lerch 2018), we employ a zero-truncated logistic distribution as parametric forecast distribution in Eq. (1), instead of a zero-truncated normal, as preliminary tests showed a slightly superior predictive performance (see supplementary material for details). Note that the zero-censoring arises from temperatures being measured in Kelvin here. Further, this allows using the same configuration as for the wind gust predictions. In particular, both the EMOS and DRN benchmark approaches are identical for both data sets.

## 6 Performance evaluation

For each of the postprocessing methods, we generated a pool of 20 networks in each forecast scenario. To ensure a fair comparison to the benchmark methods, we follow the approach from Schulz and Lerch (2022a,b), who build an ensemble of 10 networks and aggregate the forecasts via quantile aggregation. Hence, we draw 10 members from the pool and repeat this procedure 50 times to quantify the uncertainty of sampling from the general pool. For all model variants and resamples, we select those configurations as the final forecast that yield the lowest CRPS on the validation set. Details on hyperparameter tuning are discussed in the supplementary materials. For both datasets, we compute the average CRPS, PI length and PI coverage for the different forecast lead times based on the respective test datasets. The average is calculated over the resamples of the aggregated network ensembles. In what follows, the prefixes ED and ST refer to pooling-based encoder-decoder models and set transformers, respectively, and suffixes DRN and BQN indicate the parameterization of the forecast distribution. The model categories DRN and BQN without additional prefix refer to the benchmark models based on summary statistics.

### Wind gust forecasts

Table 3 shows the quantitative evaluation for lead times 6h, 12h and 18h. All permutation-invariant model architectures perform similar to the DRN and BQN benchmarks and outperform both the EPS and conventional postprocessing via EMOS, thus achieving state-of-the-art performance for all lead times. Further, the PI lengths and coverages are similar to that of the benchmark methods with the same forecast distribution, indicating that the ensemble-based models achieve approximately the same level of sharpness as the benchmark networks while being well-calibrated. Note that the underlying distribution type should be taken into account when comparing the sharpness of different postprocessing models based on the PI length, as the DRN and BQN forecast distributions exhibit different tail behavior, which affects the PI lengths for different nominal levels (see supplementary materials for details). A noticeable difference between the different network classes is that the ED models result in sharper PIs than the ST models. This coincides with the empirical PI coverages of the methods in that wider PIs typically result in a higher coverage. Fig. 2 shows the PIT histograms of the postprocessed forecasts. While differences are seen between DRN-type and BQN-type models, all DRN-type and all BQN-type models show very similar patterns. While all models are well calibrated, DRN-type models reveal limitations in the resolution of gusts in the lower segment of the distribution. BQN-type models all yield very uniform calibration histograms.

### EUPPBench surface temperature reforecasts

As shown in Table 4, also for the EUPPBench dataset both ED and ST models show significant advantages compared to the EPS and EMOS in terms of CRPS and PI length. Differences between the network variants arise mainly due to the use of different forecast distribution types. Note that the lead times of the wind gust dataset are in the short range with a maximum of 18h, whereas the lead times considered in the EUPPBench dataset range from one to five days. Hence, the differences between the lead times in the effects of postprocessing are more pronounced. E.g., for a lead time of 120h, the improvement of the network-based postprocessing methods over the conventional EMOS approach is much smaller than for shorter lead times. In particular, ST models perform the best for lead time 24h and all newly proposed models result in the smallest CRPS for lead time 120h. In terms of the PI length and coverage, we find that the ED and ST models tend to generate slightly sharper predictions. A more detailed discussion of the differences in the PI lengths due to the choice of the underlying distribution is provided in the supplementary material. The PIT histograms in Fig. 2 show that the BQN models struggle to set accurate upper and lower bounds for the predicted distribution, whereas
DRN distributions do not show such issues. Instead, they face the problem that the tail is too heavy. Overall, all postprocessing methods result in calibrated forecasts, while the DRN forecasts appear slightly better calibrated than the DRN forecasts.

\begin{table}
\begin{tabular}{c c c c c c c c c c} \hline \hline Lead Time & \multicolumn{3}{c}{6h} & \multicolumn{3}{c}{12h} & \multicolumn{3}{c}{18h} \\ \hline Method & CRPS & PI length & PI coverage & CRPS & PI length & PI coverage & CRPS & PI length & PI coverage \\ \hline EPS & 1.31 & 2.37 & 43.18 & 1.26 & 3.31 & 56.32 & 1.32 & 3.80 & 59.78 \\ EMOS & 0.88 & 5.58 & 92.83 & 0.97 & 6.01 & 91.92 & 1.04 & 6.43 & 92.46 \\ \hline BQN & 0.79 & 4.60 & 90.23 & 0.85 & 4.90 & 89.65 & 0.95 & 5.56 & 90.70 \\ DRN & 0.79 & 4.75 & 91.43 & 0.85 & 5.11 & 91.08 & 0.95 & 5.68 & 91.78 \\ \hline ED-BQN & 0.80 & 4.56 & 89.83 & 0.86 & 4.92 & 89.56 & 0.95 & 5.55 & 90.55 \\ ED-DRN & 0.79 & 4.70 & 91.17 & 0.86 & 5.15 & 91.13 & 0.95 & 5.76 & 92.07 \\ ST-BQN & 0.80 & 4.67 & 90.20 & 0.87 & 5.01 & 89.94 & 0.96 & 5.61 & 90.70 \\ ST-DRN & 0.80 & 4.77 & 91.34 & 0.86 & 5.17 & 91.13 & 0.96 & 5.83 & 92.24 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Mean CRPS in m/s, PI length in m/s and PI coverage in % of the postprocessing methods for the different lead times of the wind gust data. Recall that the nominal level of the PIs is approximately 90.48%.

Figure 2: Calibration histograms of the postprocessing models on 20-member wind gust forecasts (bottom), 11-member EUPPBench reforecast ensembles (top, left) and 51-member forecast ensembles (top, right).
BQN forecasts, yielding PIT histograms with a wave-like structure.

### Generalization to 51-member forecast ensembles

As before, postprocessing outperforms the EPS forecasts and results in calibrated and accurate forecasts (cf. Table 5 and Fig. 2). Notably, all models have been trained purely on 11-member reforecasts and are not fine-tuned to the 51-member forecast ensembles. The CRPS scores are similar with almost identical values for all models, except EMOS, for all lead times. The ST models again perform the best for the shortest lead time. For the DRN forecasts, we find that the ensemble-based networks tend to reduce the PI length, as it is smaller for all cases except lead time 120h. The corresponding PI coverages are closely connected to the length of the PIs and indicate that the PIs are too large for most postprocessing models, as the observed coverages are above the nominal level.

The calibration of the methods is not as good as in the other case studies, as indicated by the PIT histograms in Fig. 2, which may be a consequence of the large learning rate used in training the models (cf. supplementary materials). All BQN forecasts have problems in the tails, where the lower and upper bound are too extreme, such that not sufficiently many observations fall into the outer bins. DRN yields similar results as for the reforecast data with too heavy-tailed forecast distribution, as indicated by the least frequent last bin. The differences between the methods itself are again minor. Still, all postprocessing methods generate reasonably well-calibrated forecasts. Overall, the ensemble-based models result in state-of-the-art performance for generalization on 51-member forecasts or offer advantages over the summary-based benchmark methods.

## 7 Analysis of predictor importance

We analyse how the different model types distill relevant information out of the ensemble predictors. For this, we propose an ensemble-oriented permutation feature importance (PFI) analysis to assess which distribution properties of the ensemble-valued predictors have the most effect on the final prediction. In its original form, PFI (e.g., Breiman 2001; Rasp and Lerch 2018; Schulz and Lerch 2022b) is used to assign relevance scores to scalar-valued predictors by randomly shuffling the values of a single predictor across the dataset. While the idea of shuffling predictor samples translates identically from scalar-valued to ensemble-valued predictors, ensemble predictors possess internal degrees of freedom (DOFs), such as ensemble mean and ensemble range, which may affect the prediction differently. In addition to ensemble-internal DOFs, the perturbed predictor ensemble is embedded in the context of the remaining multivariate ensemble predictors, such that covariances, copulas or the rank order of the ensemble members may carry information. To account for such effects, we introduce a conditional permutation strategy that singles out the effects of different ensemble properties.

### Importance of the ensemble information

Following the notation of section 3, let \(g\) denote a postprocessing system that translates a raw ensemble forecast \(X=\left\{\mathbf{x}_{m}\in\mathcal{P}:m=1,...,M\right\}\in\left[\mathcal{P}\right] _{M}\) into a postprocessed distribution descriptor \(\theta\), and let for each member forecast \(\mathbf{x}_{m}\) and predictor channel \(i\), \(1\leq i\leq p\), \(\mathbf{x}_{m}^{(i)}\) denote the value of predictor \(i\) in the respective member. Let further \(\mathcal{D}=\left\{\left(X(t),y(t)\right):1\leq t\leq T\right\}\) be a test dataset for evaluation, consisting of \(T\in\mathbb{N}\) known raw forecast-observation pairs. Given a (negatively oriented) accuracy measure \(\bar{S}\), we write \(\bar{S}_{\mathcal{D}}\left[g\right]\) to denote the accuracy score of \(g\), subject to data \(\mathcal{D}\). From here, we choose \(\bar{S}\) to be the expected CRPS, and assume that all scores are computed based on the same test dataset, allowing us to drop the dataset index, i.e., \(\bar{S}_{\mathcal{D}}\left[g\right]\equiv\bar{S}\left[g\right]\). In this notation, the relative FPI, as used in Schulz and Lerch (2022b), can be written as

\[\Delta_{0}(P):=\frac{\bar{S}[g\circ P]-\bar{S}[g]}{\bar{S}[g]}, \tag{8}\]

wherein \(P\) indicates a perturbation operator that alters parts of the predictor data, and \(\circ\) denotes function composition. For the classical FPI, we denote the permutation operator as \(\Pi_{\pi}^{(i)}\), which shuffles the \(i\)-th predictor channel of the raw ensembles according to a permutation \(\pi\) of the dataset \(\mathcal{D}\) (omitted in the notation for brevity).

For ensemble-valued predictors, we consider two generalizations of this operator. We refer to these as the fully-random permutation, \(\Pi_{\pi}^{(i)}\), and the rank-aware random permutation, \(\tilde{\Pi}_{\pi}^{(i)}\). The former acts as a direct analog of the scalar-valued permutation case, i.e., for all \(1\leq t\leq T\) and \(1\leq m\leq M\), it replaces the values \(\mathbf{x}_{m}^{(i)}\) of the ensemble \(X(t)\) with arbitrary values \(\mathbf{x}_{m^{\prime}}^{(i)}\), \(1\leq m^{\prime}\leq M\), from the ensemble \(X(\pi(t))\), without replacement. Thus, it destroys all information of the original ensemble. The latter ranks the member values \(\mathbf{x}_{m}^{(i)}\) in \(X(t)\) and replaces them with values \(\mathbf{x}_{m^{\prime}}^{(i)}\) from \(X(\pi(t))\), where \(m^{\prime}\) are chosen such that all members are used exactly once and the perturbed ensemble possesses the same ranking order as the original one. It thus preserves the ordering of the perturbed predictors in the context of the remaining predictors. In practice, we note that the differences in feature importance for both variants are very minor, such that we select only the rank-aware variant for further analysis.

To probe the importance of ensemble-internal DOFs, we consider additional perturbation operators, which rely on conditional shuffling of the ensemble predictors. For this, let \(s:\left[\mathbb{R}\right]_{M}\rightarrow\mathbb{R}\) be a summary function, which translates an ensemble of scalar predictor values into a real
valued summary statistic, such as ensemble mean or standard deviation. Then an \(s\)-conditional shuffling operator \(\Pi_{(\pi_{b})\mid s}^{(i)}\) is defined as follows. For all raw predictions \(X(t)\) in the dataset, the predictor ensemble for the \(i\)-th predictor, \(X^{(i)}(t)=\{\mathbf{x}_{m}^{(i)}:\mathbf{x}_{m}\in X(t)\}\), is extracted and summary statistics \(s(X^{(i)}(t))\) are computed. The observed summary statistics are ranked from 1 to \(T\) and the corresponding ensembles \(X(t)\) are distributed into \(B\in\mathbb{N}\) evenly spaced bins, according to these ranks. For each bin \(b\), \(0\leq b<B\), a permutation \(\pi_{b}\) is sampled randomly and the values of the \(i\)-th predictor are shuffled bin-wise according to these permutations. For suitably sized bins, the shuffling preserves information about \(s\) and erases information about other DOFs, thereby ensuring that each of the bins contains an approximately equal number of samples, independent of the details of the predictor distribution. In our experiments, \(B=100\) bins yielded a good balance between information preservation and randomization. Results for larger and smaller bin sizes were qualitatively similar. Note that for predictors in which certain values appear with large multiplicity, such as zero in censored variables like precipitation, the ranking is computed on the unique values of the summary statistics. This enforces a small amount of variation even in bins with degenerate values. In analogy to the rank-aware (unconditional) shuffling, the rank-aware \(s\)-conditional shuffling is denoted as \(\tilde{\Pi}_{(\pi_{b})\mid s}^{(i)}\). For the conditional FPI analysis, we suggest the computation of importance ratios,

\[\chi(P|R):=\frac{\tilde{S}[g\circ P]-\tilde{S}[g]}{\tilde{S}[g\circ R]-\tilde{ S}[g]}, \tag{9}\]

which measure the fraction of skill restored (or destroyed) by applying a shuffling operation \(P\) instead of a reference operation \(R\). The ratios of interest are \(\chi(\tilde{\Pi}_{(\pi_{b})\mid s}^{(i)},\tilde{\Pi}_{\pi}^{(i)})\), which measure how much of the prediction skill deficit due to randomized shuffling of predictor \(i\) is restored by preserving information about the summary statistic \(s\). In absence of sampling errors due to finite data, \(\chi\left(\tilde{\Pi}_{(\pi_{b})\mid s}^{(i)},\tilde{\Pi}_{\pi}^{(i)}\right)\) yields values between 0 and 1, with 0 indicating uninformative summary statistics, and 1 suggesting that knowledge of \(s\) is sufficient to restore the original model skill entirely. Empirically, we find that the theoretical bounds are preserved well for predictors with sufficiently large FPI.

\begin{table}
\begin{tabular}{c c c c c c c c c c} \hline \hline Lead Time & \multicolumn{2}{c}{24h} & \multicolumn{2}{c}{72h} & \multicolumn{2}{c}{120h} \\ \hline Method & CRPS & PI length & PI coverage & CRPS & PI length & PI coverage & CRPS & PI length & PI coverage \\ \hline EPS & 1.21 & 1.81 & 39.85 & 1.28 & 3.29 & 56.62 & 1.54 & 4.89 & 63.44 \\ EMOS & 0.82 & 3.85 & 82.56 & 0.96 & 4.72 & 83.96 & 1.25 & 6.05 & 83.12 \\ \hline BQN & 0.67 & 3.32 & 84.73 & 0.87 & 4.44 & 86.08 & 1.20 & 6.24 & 86.09 \\ DRN & 0.67 & 3.28 & 84.16 & 0.86 & 4.27 & 84.58 & 1.19 & 5.70 & 83.09 \\ \hline ED-BQN & 0.67 & 3.29 & 84.57 & 0.87 & 4.45 & 86.05 & 1.19 & 6.03 & 85.45 \\ ED-DRN & 0.67 & 3.19 & 83.39 & 0.87 & 4.25 & 84.11 & 1.19 & 5.64 & 82.60 \\ ST-BQN & 0.66 & 3.16 & 84.01 & 0.87 & 4.31 & 84.85 & 1.19 & 6.07 & 85.15 \\ ST-DRN & 0.66 & 3.06 & 82.67 & 0.87 & 4.18 & 83.44 & 1.19 & 5.77 & 83.17 \\ \hline \hline \end{tabular}
\end{table}
Table 4: Mean CRPS in K, PI length in K and PI coverage in % of the postprocessing methods for the different lead times for the EUPPBench reforecast data. Recall that the nominal level of the PIs is approximately 83.33 %.

\begin{table}
\begin{tabular}{c c c c c c c c c c} \hline \hline Lead Time & \multicolumn{2}{c}{24h} & \multicolumn{2}{c}{72h} & \multicolumn{2}{c}{120h} \\ \hline Method & CRPS & PI length & PI coverage & CRPS & PI length & PI coverage & CRPS & PI length & PI coverage \\ \hline EPS & 1.21 & 2.65 & 57.54 & 1.18 & 4.71 & 74.78 & 1.38 & 7.14 & 83.26 \\ EMOS & 0.79 & 6.31 & 96.26 & 0.90 & 7.74 & 97.49 & 1.16 & 9.92 & 97.47 \\ \hline BQN & 0.64 & 4.32 & 94.13 & 0.80 & 6.52 & 97.23 & 1.13 & 9.18 & 97.58 \\ DRN & 0.64 & 5.48 & 97.92 & 0.80 & 7.21 & 98.37 & 1.13 & 9.58 & 98.28 \\ \hline ED-BQN & 0.64 & 4.74 & 96.30 & 0.81 & 6.49 & 97.42 & 1.12 & 8.81 & 97.15 \\ ED-DRN & 0.64 & 5.31 & 97.62 & 0.81 & 7.09 & 98.19 & 1.12 & 9.61 & 97.90 \\ ST-BQN & 0.62 & 4.61 & 95.96 & 0.80 & 6.18 & 96.31 & 1.13 & 8.68 & 96.05 \\ ST-DRN & 0.62 & 5.10 & 97.40 & 0.81 & 6.88 & 97.55 & 1.13 & 9.43 & 97.11 \\ \hline \hline \end{tabular}
\end{table}
Table 5: Mean CRPS in K, PI length in K and PI coverage in % of the postprocessing methods for the different lead times for EUPPBench forecast data. Recall that the nominal level of the PIs is approximately 96.15%.
### _Results_

We compute FPI scores \(\Delta_{0}\left(\Pi_{\pi}^{(i)}\right)\) for all ensemble predictors and model variants. Fig. 3 depicts a selection of the FPI scores of the most important ensemble-valued predictors in both tasks. A figure with all ensemble-valued predictors is shown in the supplementary materials. Scalar-valued predictors (cf. section 3d for the terminology) are omitted for easier comparison with the conditional importance measures. The boxplots indicate statistics obtained from 20 separate model runs, which have been evaluated independently, the bars show the mean importance. The accuracy of the wind gust models is dominated by VMAX-10M, and supplemented by additional predictors with lower importance. Temperature-like predictors obtain similar or higher scores than, e.g. winds at 850 hPa and 950 hPa pressure levels. Note that for each lead time, the importance highlights different temperature predictors, which may be attributed to the diurnal cycle. Similar arguments can explain the increasing importance of ASOB-S (short-wavelength radiation balance at the surface) with increasing lead time. In a direct comparison of the model variants, we find that the differences between BQN-type and DRN-type models are very minor. However, ED-type models attribute higher importance to the most relevant predictors (VMAX-10M, T1000, T-2M), whereas ST-type models distribute the importance more evenly and use more diverse predictor information.

In the EUPPBench case, the models focus mainly on temperature-like predictors as well as surface radiation balances. Notably, for the summary-based models, mn2t6 and mx2t6 tend to be more important than the primary predictor t2m up to lead time 72h. Since the diurnal cycle does not cause variations between the lead times here, differences in the predictor utilization must be due to the increasing uncertainty at longer lead times. The ensemble-based models rely relatively more strongly on the t2m predictor for the shorter lead time, whereas for longer lead times, the information utilization is more diverse. Qualitative differences between ED- and ST-type models are observed with respect to the humidity-related predictors tcw and tcw. Only ST models recognize the value in these predictors, which may explain in parts the different generalization properties of ED and ST models on the EUPPBench reforecast and forecast datasets.

Figs. 4 and 5 investigate the importance of ensemble-internal DOFs of selected ensemble predictors for the permutation-invariant model architectures. For both datasets, we choose a set of representative high-importance predictors, and display the DOF importance for the ensemble-based models. For all predictors and lead times, we compute importance ratios \(\chi\left(\tilde{\Pi}_{(\pi_{b})|s}^{(i)},\tilde{\Pi}_{\pi}^{(i)}\right)\) for a selection of commonly used ensemble summary statistics. Specifically, we consider the ensemble mean as a proxy for the location of the distribution, ensemble maximum and minimum to assess the impact of extreme values, standard deviation, inter-quartile range and full range (difference between maximum and minimum) to quantify the scale of the distribution, as well as skewness and kurtosis as higher-order summary statistics. Due to the pairwise similarity of some of the measures, it is to be expected that conditional shuffling with respect to one of the measures preserves information also about others. To assess the information overlap between shuffling patterns with different reference statistics, Spearman rank correlations are computed between the shuffled statistics and the original statistics. The resulting correlation matrices illustrate how accurately the rank order for one statistic is preserved if the data is conditionally shuffled with respect to another. Rank correlations are chosen to minimize the effect of the marginal distribution of the respective statistics values, since these may vary considerably between different predictors and summary statistics. The results are depicted as heatmaps in Figs. 4 and 5.

For wind gust postprocessing (Fig. 4), the importance ratios suggest in many cases that virtually all of the predictor information can be restored by conditioning the shuffling procedure on the ensemble mean. Notably, this is the case for VMAX-10M, consistently for all lead times, as well as for T-G and FI850. The interaction plots suggest that the mean-conditioning preserves information about extrema to a high degree, whereas ensemble range and higher-order statistic information are mixed up. These findings are supported by observations in Schulz and Lerch (2022b), who note that omitting the standard deviation of the auxiliary ensemble predictors helps to improve the quality of the network predictions.

ASOB-S and WIND850 are interesting corner cases, in which the mean-conditioning restores substantial amounts of the model skill, but fails to reproduce it completely. This indicates that, while the ensemble mean is an important predictor, the remaining DOFs deliver complementary information that modulates the interpretation of the mean value. Note that the information overlap between location-like and scale-like metrics for ASOB-S predictors at 6h lead time is again an artifact due to the diurnal cycle. At 6h lead time, a substantial fraction of the ASOB-S predictor ensembles fall zero mean and no variance due to the lack of solar irradiation, which impacts the correlation values. The relevance scores for FI850 statistics suggest possible higher-order interactions. For all lead times, it is possible to restore half of the lost model skill by conditioning on skewness and kurtosis, and the corresponding correlation plots suggest that this cannot be attributed to information overlap with other predictors, because correlations are consistently close to zero. An explanation of this observation could be that a fully randomized shuffling of the predictor ensembles destroys the information coherence with other predictors.
In surface temperature postprocessing, t2m is an interesting case, in which for all lead times neither of the summary statistics alone is sufficient to restore the unperturbed model performance. This indicates that both ED- and ST-type models learn to attend to the details of the ensemble distribution and marks a difference to the wind gust case study, where most of the information is conveyed in the ensemble means. With increasing lead time, the mean-conditional shuffling becomes more effective in restoring the model skill. This may be due to the decreasing reliability of the EPS prediction system with increasing lead time. Similar patterns are observed also in the remaining predictors. While the model skill cannot be restored with mean-only conditioning for 24h lead time, the mean appears to become more informative for longer lead times. The radiation parameter ssrd6 sticks out visually with high correlations between location-related predictors, which occurs due to the same reasons as for the ASOB-S parameter discussed before.

## 8 Discussion and Conclusion

We have introduced permutation-invariant NN architectures for postprocessing ensemble forecasts by selecting two exemplary model families and adapting them to postprocessing. In two case studies, using datasets for wind gusts and surface temperature postprocessing, we have validated the model performance and compared the permutation-invariant models against benchmark models from prior work. Our results show that permutation-invariant postprocessing networks achieve state-of-the-art performance in both applications. All permutation-invariant architectures outperform both the raw ensemble forecast and conventional postprocessing via EMOS by a large margin, but no systematic differences can be observed between the (more complex) permutation-invariant models and existing NN-based solutions.

Based on a subsequent assessment of the permutation importance of ensemble-internal DOFs, we have seen that for many auxiliary ensemble predictors, preserving information about the ensemble mean is sufficient to maintain almost the complete information about the postprocessing target, while more detailed information is required about the primary predictors. These findings are consistent with prior work and are more comprehensive due to the larger variety of summary statistics considered in the analysis.

A striking advantage of the permutation-invariant models lies in the generality of the approach, i.e., the models possess the flexibility to attend to the important features in

Figure 3: Permutation feature importance for summary-based networks (top) and permutation-invariant models (bottom) for EUPPBench and wind gust postprocessing. Predictors named _ens_ in the top figure correspond to the primary predictors t2m and VMAX-10M, respectively. The suffix _sd_ indicates the ensemble standard deviation of the predictor.
the predictor ensembles and they are capable of identifying those during training (as shown in our feature analysis). As the added flexibility comes with a surplus of computational complexity, the benefits of the respective methods should be weighed carefully. In operational settings, it may be reasonable to consider permutation-invariant models, as proposed here, as a tool for identifying relevant aspects of the input data. The gained knowledge can then be used for data reduction and to train reduced models with a more favorable accuracy-complexity trade-off.

Despite these advantages, the apparent similarity between the performance of the ensemble-based and summary-based models remains baffling and requires further clarification. Supposing capable ensemble predictions, it seems reasonable, from a meteorological perspective, to expect that postprocessing models that operate on the entire ensemble can learn more complex patterns and relationships than models that operate on simple summary statistics. The lack of substantial improvements, as seen in this study, admits different explanations. One possibility would be that the available datasets are insufficient to establish statistically relevant connections between higher-order ensemble-internal patterns and the predicted variables. Problems could arise, e.g. due to insufficient sample counts of the overall datasets or due to ensemble sizes being too low to provide reliable representations of the forecast distribution. Yet, another reason could lie in the fact that the generation mechanisms underlying the NWP
ensemble forecasts fail to achieve meaningful representations of such higher-order distribution information, which would raise follow-up questions regarding the design of future ensemble prediction systems. Given the impact and potential implications of the latter alternative, future work should examine the information content of raw ensemble predictions in more detail. The proposed permutation-invariant model architectures may help to achieve this, e.g., by conducting postprocessing experiments with dynamical toy systems that are cheap to simulate and simple to understand, such that large datasets can be generated and evidence for both hypotheses can be distinguished.

Acknowledgments.This research was funded by the subprojects B5 and C5 of the Transregional Collaborative Research Center SFB/TRR 165 "Waves to Weather" (www.wavestoweather.de) funded by the German Research Foundation (DFG). Sebastian Lerch gratefully acknowledges support by the Vector Stiftung through the Young Investigator Group "Artificial Intelligence for Probabilistic Weather Forecasting".

Data availability statement.The case study on surface temperature postprocessing is based on the EUPPBench dataset, which is publicly available. See Demaeyer et al. (2023) for details. The wind gust dataset is proprietary but can be obtained from the DWD for research purposes. Code with implementations of all methods is publically available (Hohlein, 2023).

Figure 5: Importance of ensemble-internal DOFs for temperature postprocessing. Same as Fig. 4.
## Appendix A Description of predictors

The descriptions of the ensemble-valued predictor variables used in both case studies are shown in Tables A1 and A2 for wind-gust and surface-temperature postprocessing, respectively. The predictors listed in Table A3 are not ensemble-valued and are used equally in both cases.
\begin{table}
\begin{tabular}{l l l l} \hline \hline Short name & Units & Full name & Levels \\ \hline VMAX & m/s & Maximum wind, i.e. wind gusts & 10 m \\ U & m/s & U-component of wind & 10 m, 1000 hPa, 950 hPa, 850 hPa, 700 hPa, 500 hPa \\ V & m/s & V-component of wind & 10 m, 1000 hPa, 950 hPa, 850 hPa, 700 hPa, 500 hPa \\ WIND & m/s & Wind speed, derived from U and V via \(\sqrt{\text{U}^{2}+\text{V}^{2}}\) & 10 m, 1000 hPa, 950 hPa, 850 hPa, 700 hPa, 500 hPa \\ OMEGA & Pa/s & Vertical velocity (Pressure) & 1000 hPa, 950 hPa, 850 hPa, 700 hPa, 500 hPa \\ \hline T & K & Temperature & Ground-level, 2 m, \\  & & 1000 hPa, 950 hPa, 850 hPa, 700 hPa, 500 hPa \\ T-D & K & Dew point temperature & 2 m \\ \hline RELIUM & \% & Relative humidity & 1000 hPa, 950 hPa, 850 hPa, 700 hPa, 500 hPa \\ TOT-PREC & kg/m\({}^{2}\) & Total precipitation (Acc.) & – \\ RAIN-GSP & kg/m\({}^{2}\) & Large scale rain (Acc.) & – \\ SNOW-GSP & kg/m\({}^{2}\) & Large scale snowfall - water equivalent (Acc.) & – \\ W-SNOW & kg/m\({}^{2}\) & Snow depth water equivalent & – \\ W-SO & kg/m\({}^{2}\) & Column integrated soil moisture & multilayers: 1, 2, 6, 18, 54 \\ CLC & \% & Cloud cover & T: total; \\  & & & L: soil to 800 hPa; M: 800 to 400 hPa; H: 400 to 0 hPa \\ HBAS-SC & m & Cloud base above mean sea level, shallow connection & – \\ HTOP-SC & m & Cloud top above mean sea level, shallow connection & – \\ \hline ASOB-S & W/m\({}^{2}\) & Net short wave radiation flux & surface \\ ATHB-S & W/m\({}^{2}\) & Net long wave radiation flux (m) & surface \\ ALB-RAD & \% & Albedo (in short-wave) & – \\ \hline PMSL & Pa & Pressure reduced to mean sea level & – \\ FI & m\({}^{2}\)/s\({}^{2}\) & Geopotential & 1000 hPa, 950 hPa, 850 hPa, 700 hPa, 500 hPa \\ \hline \hline \end{tabular}
\end{table}
Table 10: Description of meteorological parameters for wind-gust postprocessing (cf. Schulz and Lerch 2022b). Target variable: wind speed of gust (observations). Primary predictor: VMAX-10m (ensemble forecast).
\begin{table}
\begin{tabular}{l l l} \hline \hline Predictor & Type & Description \\ \hline yday & Temporal & Cosine transformed day of the year \\ lat & Spatial & Latitude of the station \\ lon & Spatial & Longitude of the station \\ alt & Spatial & Altitude of the station \\ orog & Spatial & Difference of station altitude and model surface height of nearest grid point. \\ loc-bias & Spatial & Mean bias of ensemble forecasts, computed from the training data. \\ loc-cover & Spatial & Mean coverage of ensemble forecasts, computed from the training data. \\ \hline \hline \end{tabular}
\end{table}
Table 10: Experiment of meteorological parameters for surface temperature postprocessing (EUPPBench, cf. Demaeyer et al. 2023). Target variable: t2m (observations). Primary predictor: t2m (ensemble forecast).