Ensemble weather forecast post-processing with a flexible probabilistic
  neural network approach
Ensemble forecast post-processing is a necessary step in producing accurate
probabilistic forecasts. Conventional post-processing methods operate by
estimating the parameters of a parametric distribution, frequently on a
per-location or per-lead-time basis. We propose a novel, neural network-based
method, which produces forecasts for all locations and lead times, jointly. To
relax the distributional assumption of many post-processing methods, our
approach incorporates normalizing flows as flexible parametric distribution
estimators. This enables us to model varying forecast distributions in a
mathematically exact way. We demonstrate the effectiveness of our method in the
context of the EUPPBench benchmark, where we conduct temperature forecast
post-processing for stations in a sub-region of western Europe. We show that
our novel method exhibits state-of-the-art performance on the benchmark,
outclassing our previous, well-performing entry. Additionally, by providing a
detailed comparison of three variants of our novel post-processing method, we
elucidate the reasons why our method outperforms per-lead-time-based approaches
and approaches with distributional assumptions.
Peter Mlakar, Janko Merše, Jana Faganeli Pucer
# Ensemble weather forecast post-processing with a flexible probabilistic neural network approach

###### Abstract

Ensemble forecast post-processing is a necessary step in producing accurate probabilistic forecasts. Conventional post-processing methods operate by estimating the parameters of a parametric distribution, frequently on a per-location or per-lead-time basis, which limits their expressive power. We propose a novel, neural network-based method, which produces forecasts for all locations and lead times, jointly. To relax the distributional assumption made by many post-processing methods, our approach incorporates normalizing flows as flexible parametric distribution estimators. This enables us to model varying forecast distributions in a mathematically exact way. We demonstrate the effectiveness of our method in the context of the EUPPBench benchmark, where we conduct temperature forecast post-processing for stations in a sub-region of western Europe. We show that our novel method exhibit state-of-the-art performance on the benchmark, outclassing our previous, well-performing entry. Additionally, by providing a detailed comparison of three variants of our novel post-processing method, we elucidate the reasons why our method outperforms per-lead-time-based approaches and approaches with distributional assumptions.

## 1 Introduction

Forecast post-processing is a crucial task when constructing skillful weather forecasts. This is due to the inherent biases of the numerical weather predictions (NWP) which are the results of initial condition errors, computational simplifications, and sub-grid parametrizations [3, 40, 18]. All these factors combine and result in errors that compound as the forecast horizon increases. Therefore, due to these issues, the forecast is uncertain by nature. To quantify this uncertainty ensemble forecasts are issued by centers for weather forecasts, such as the European Center for Medium Range Weather Forecasts (ECMWF) [11]. ECMWF generates ensemble forecasts by perturbing the initial atmospheric conditions and varying the NWP parametrizations. This results in 50 different ensemble members, each expressing a potentially unique future weather situation. However, biases are still present in ensemble forecasts resulting in a disconnect between the specific observations and the described ensemble distribution.

To mitigate the effects of biases in forecasts with the aim of constructing calibrated forecast probability distributions [13] forecast providers frequently employ forecast post-processing techniques [40]. These statistical approaches transform raw ensemble forecasts to better concur with actual weather variable observations. Post-processing methods vary from simple to complex, from rigid [15, 33, 36, 30] to flexible [24, 21, 29, 39], and encompass a large spectra of statistical approaches [40]. Another way of distinguishing between different post-processing methods is to characterize them based on their distributional assumptions or lack thereof. Since the nature of ensemble forecasts is inherently probabilistic (ensemble forecasts quantify the weather forecast uncertainty) the post-processing output should reflect that. Methods that assume the target weather variable distribution are usually simpler to implement however, they suffer from the pitfall of assuming the wrong distribution. This can result in additional biases. Likewise, different weather variables can exhibit different distributions therefore, tuning and exploration of the best target distribution is required before constructing the model. Methods that do not make such assumptions are usually more complex to implement and require more data for parameter estimation.

Recently neural networks started showing promising results in this field [27, 26, 5, 37, 41, 4, 17, 35, 34, 9, 27, 14]. Applications such as those by [34] exhibited state of the art results compared to conventional post-processing
techniques. By using the ensemble mean and variance with the combination of station embeddings [34] construct a neural network that estimates the parameters of a normal distribution for each lead time. This approach is further extended in [37] by applying different distribution estimators in place of the normal distribution. To be more specific, they apply the same neural network architecture as a base parameter predictor to three different models: a distributional regression model, a Bernstein quantile regression model [4], and a histogram estimation model. Indeed, their aim to reduce the distributional assumption is apparent and also required as many weather variables can not be effectively described using simple probability distributions. We can make a similar observation by investigating the work by [41] where the authors construct a convolutional neural network for post-processing wind speeds. Similarly to [37], they tried different probability models as the output of their convolutional neural network: a quantized softmax approach, kernel mixture networks [1], and truncated normal distribution fitting approach. [41] concluded that the most flexible distribution assessment method exhibited the best performance. This further bolsters the need for sophisticated distribution estimation techniques as these increase the performance of post-processing algorithms.

To further improve upon the existing post-processing techniques we propose two novel neural network approaches, Atmosphere NETwork 1 (ANET1) and Atmosphere NETwork 2 (ANET2). The novelty of the ANET methods is threefold and can be summarized thusly:

* Two novel neural network architectures for post-processing ensemble weather forecasts
* Joint probabilistic forecasting for all lead times and locations
* Flexible parametric distribution estimation based on normalizing spline flows, implemented in ANET2

Both ANET variants are constructed such that there is only one model for the entire post-processing region and all lead times. This is often not the case with frequently used post-processing methods [6] which typically construct individual models per-lead-time and sometimes even per-station. Our joint location-lead-time method ANET1 exhibits state-of-the-art performance on the EUPPBench benchmark [6] compared to other submitted methods, which are implemented on a per-lead time or per station basis.

ANET2 is the second iteration of our post-processing approach and features an optimized neural network architecture (similar to [34]) and training procedure. Furthermore, ANET2 improves upon ANET1 by relaxing the distribution assumption of ANET1, which assumes that the target weather variable is distributed according to a parametric distribution which, in the case of temperature, is a normal distribution. ANET2 relaxes this assumption by modeling the target distribution in terms of a normalizing spline flow [10, 8, 25, 22]. This enables us to model mathematically exact distributions without specifying a concrete target distribution family. We demonstrate that ANET2 further improves upon ANET1 in a suite of performance metrics tailored for probabilistic forecast evaluation. Additionally, we provide intuition as to why such a type of joint forecasting and model construction leads to better performance by analyzing the feature importance encoded in the ANET2 model.

To demonstrate the performance of the ANET variants we compare four different novel methods for probabilistic post-processing of ensemble temperature forecasts: ANET1, ANET2, ANET2\({}_{\text{NORM}}\), and ANET2\({}_{\text{BERN}}\). This helps us quantify the impact different neural network architectures and distributional models have on the final post-processing performance of our approaches. The description of the aforementioned methods and training procedures is available in Section 2. We follow this with the evaluation results in Section 3 which we further elaborate on in Section 4.

## 2 Methods

### Anet1

ANET1 is a neural network based approach for probabilistic ensemble forecast post-processing. The ANET1 architecture, displayed in Figure 1, is tailored for processing ensemble forecasts with a varying number of ensemble members. ANET1 achieves this by first processing ensemble forecasts for the entire lead time individually, by passing them through a shared forecast encoder structure. We concatenate per-station predictors to each ensemble forecast to allow ANET1 to adapt to individual station conditions. These predictors include station and model altitude, longitude, latitude, and land usage, with the addition of a seasonal time encoding, defined as \(\cos(\frac{2\pi d}{365})\), where \(d\) denotes the day of year the individual input forecast was issued. The shared forecast encoder transforms each forecast-predictor pair into a high-dimensional latent encoding, which are then weighed by a dynamic attention block and averaged into a single, mean ensemble encoding. This dynamic attention mechanism is implemented with the goal of determining the importance of individual ensemble members in a given weather situation. The mean
ensemble encoding is then passed through a regression block, the outputs of which are the per lead time additive corrections to the ensemble mean and standard deviation. We use the corrected mean and standard deviation as parameters of the predictive normal distribution.

### Anet2

ANET1 suffers from the distributional assumption drawback which mitigates its expressive power. To improve upon ANET1 we developed ANET2 which utilizes normalizing spline flows [10] in place of the normal distribution. Normalizing flows are designed for density estimation and can approximate complex distributions in a tractable and mathematically exact way. ANET2 embraces this methodology, albeit in a modified manner which better suits our context of application, by conditioning the final distributional model on the provided raw ensemble forecast.

#### 2.2.1 Flexible parametric distribution estimation using modified normalizing flows

We begin our overview of the normalizing spline flow density estimation procedure by defining a rational-quadratic spline transformation \(\tilde{T}(x;\theta)\), parameterized by \(\theta\), with \(x\) as the target temperature realization. The transformation \(\tilde{T}_{\theta}\) is a strictly increasing and differentiable function and is the key behind the expressive power of normalizing spline flows. The parameter set \(\theta\) includes the spline knots and the spline's values at those knots. Since our goal is to model the target temperature distribution \(F_{temp}\) we can use the change of variable approach [25] to express this unknown distribution in terms of the transformation \(\tilde{T}_{\theta}\) and a base distribution \(F_{norm}\) which, in our case, is a univariate normal distribution. The target variable density can then be defined as

\[p_{temp}(x;\theta)=p_{norm}(\tilde{T}_{\theta}(x))\frac{\partial\tilde{T}_{ \theta}(x)}{\partial x}.\]

Since \(p_{norm}\) is the density of a normal distribution with zero mean and standard deviation of one, the final loss function minimized (the negative log-likelihood) for a given sample \(x_{i}\) is

\[L(p_{temp}(x_{i}))=\frac{\tilde{T}_{\theta}(x_{i})^{2}}{2}-\ln(\frac{ \partial\tilde{T}_{\theta}(x_{i})}{\partial x}).\]

In practice we replace \(\tilde{T}_{\theta}\) with a composition of rational-quadratic spline transformations \(T_{\Theta}\), where

\[T_{\Theta}=\tilde{T}_{\theta_{i}}(\tilde{T}_{\theta_{i-1}}(...)),\]

Figure 1: The neural network architecture of ANET1. The parameter f denotes the number of features in a dense layer. The parameters t and p denote the lead time and the number of per-lead-time parameters required by the normal distribution model. The circular blocks denote the following operations: \(\mathbf{M}\) and \(\mathbf{S}\) the ensemble forecast mean and standard deviation, \(\mathbf{X}\) the element-wise product, \(\mathbf{C}\) concatenation operation, and \(\mathbf{F}\) the softplus activation [32]. The square blocks denote the following operations: \(\mathbf{E}\) corresponds to the shared forecast encoder block, \(\mathbf{A}\) corresponds to the dynamic attention block, and \(\mathbf{R}\) corresponds to the regression block. The variable \(\mathrm{e_{i}}\) denotes the i-th ensemble member in the forecast.
and \(\Theta\) denotes the entire set of individual spline parameters \(\theta_{l}\). ANET2 uses a composition of four rational-quadratic spline transformations, each described by four knot-value pairs. Since ANET2 generates a distribution for each lead time this entails that the parametric distribution described by this model contains 840 parameters (\(21\times 10\times 4\): 21 for each lead time, 10 for each of the 4 spline transformations). To fully describe a rational-quadratic spline we require the spline's knot-value pairs and the value of the spline's derivative at those knots. The knot-value pairs are estimated by the ANET2 neural network (described in more detail in Section 2.2.4). For a specific ensemble forecast with a lead time of \(t\) steps, ANET2 computes \(t\) parameter sets \(\Theta_{j}\), where \(j\in[1,t]\) (in our case, \(t=21\)).

Each \(\theta_{l}\) is comprised of two vectors, where the first contains the knots and the second the values of the spline. More formally,

\[\theta_{l}:=\{\mathbf{k}_{l},\mathbf{v}_{l}\},\]

where the vector \(\mathbf{k}_{l}\) denotes the spline knots and \(\mathbf{v}_{l}\) the values. We use a neural network to estimate these values based on the input ensemble forecast. The monotonicity of the spline knot-value pairs is ensured using the Softplus function [32]. Therefore, the final vectors of knots \(\mathbf{k}_{l}\) and values \(\mathbf{v}_{l}\) are defined as

\[\mathbf{k}_{l} :=\text{CumSum}([\mathbf{k}^{\prime}_{l,1},1\text{e}^{-3}+\text{ Softplus}(\mathbf{k}^{\prime}_{l,[2,5]})]),\] \[\mathbf{v}_{l} :=\text{CumSum}([\mathbf{v}^{\prime}_{l,1},1\text{e}^{-3}+\text{ Softplus}(\mathbf{v}^{\prime}_{l,[2,5]})]),\]

where \(\mathbf{k}^{\prime}_{l}\) and \(\mathbf{v}^{\prime}_{l}\) refer to the raw neural network outputs for the knot-value pairs, \(\mathbf{k}^{\prime}_{l,1}\) refers to the first element of the vector \(\mathbf{k}^{\prime}_{l}\), while \(\mathbf{k}^{\prime}_{l,[2,5]}\) contains the remaining elements (the same holds true for \(\mathbf{v}^{\prime}_{l}\). The index into \(\mathbf{k}^{\prime}_{l}\) and \(\mathbf{v}^{\prime}_{l}\) goes up to five as we only have five knot-value pairs per spline. Additionally, we limit the minimal distance between two consecutive knot-value pairs to be \(1\text{e}^{-3}\). We implement this restriction to ensure numerical stability and in our testing this does not impact the regression performance. The function CumSum denotes the cumulative sum operation, where the first element is kept identical to that of the input vector. This in combination with the Softplus function ensures that the knots-value pairs increase monotonically, in turn guaranteeing the monotonicity of the rational quadratic spline transformation [16].

However, to provide a full estimate of the transformation we also require the positive derivative values at the knots, denoted as \(d_{l,i}\), where \(i\in[1,...,5]\). [10] propose that we determine these derivatives from the neural network in much the same way as the knot-value pairs. This definition of the derivative can yield discontinuities in the splines and therefore in the final estimated density. This is not practical for our application context as we frequently inspect the probability density to determine the most likely weather outcome. A density full of discontinuities conveys unnatural properties and impedes the required analysis. To rectify this we look at the work by [16], where they propose two derivative estimation schemes that compute the required derivatives from the knot-value pairs. We adapt the second approach in our work, which is expressed as

\[d_{l,i}=\frac{\Delta_{l,i}\cdot\Delta_{l,i-1}}{\delta_{l,j}}\]

for \(i\in[2,4]\). The derivative values at the edge knots are defined as

\[d_{l,1}=\frac{\Delta_{l,1}^{2}}{\delta_{l,2}},\;\;d_{l,5}=\frac{\Delta_{l,4}^{ 2}}{\delta_{l,4}},\]

where \(\Delta_{l,j}\) and \(\delta_{l,j}\) are computed by

\[\Delta_{l,j}=\frac{\mathbf{v}_{l,j+1}-\mathbf{v}_{l,j}}{\mathbf{k}_{l,j+1}-\mathbf{k}_{l,j}}, \;\;\delta_{l,j}=\frac{\mathbf{v}_{l,j+1}-\mathbf{v}_{l,j-1}}{\mathbf{k}_{l,j+1}-\mathbf{k}_{l,j-1}}.\]

Due to our constraint on the minimal differences between consecutive knots and values we do not have to concern ourselves with specific edge cases outlined in [16] where the differences in above equalities would be zero.

Additionally, we do not restrict the spline knots to a predetermined interval. This is contrary to [10], however, in our testing, this constraint relaxation works well and provides no tangible difference in performance. It does, however, eliminate the need for parameter tuning as the knot ranges are determined automatically during optimization.

#### 2.2.2 Anet2norm

ANET2NORM combines the ANET2 neural network model described in Section 2.2.4 with a normal distribution acting as its probability distribution model. Therefore, we model the mean and standard deviation vectors of the
distribution, each containing 21 values corresponding to the lead times. Let us denote the raw ensemble mean and standard deviation vectors as \(\mathbf{\mu}^{\rm e},\mathbf{\sigma}^{\rm e}\). We compute the final model mean and standard deviation vectors \(\mathbf{\mu},\mathbf{\sigma}\) as

\[\mathbf{\mu} =\mathbf{\mu}^{\rm e}+\mathbf{\mu}^{\prime}, \tag{1}\] \[\mathbf{\sigma} =\mathbf{\sigma}^{\rm e}+\text{Softplus}(\mathbf{\sigma}^{\prime}), \tag{2}\]

where \(\mathbf{\mu}^{\prime},\mathbf{\sigma}^{\prime}\) denote the raw neural network estimates for those parameters. We again use the Softplus function to enforce the positivity constraint on the standard deviation. We found that if we use the neural network to predict the mean and standard deviation correction residual (expressed as an additive correction term to the raw ensemble statistics) the model performs better than it would if it were to directly predict the target mean and standard deviation without the residual.

#### 2.2.3 Anet2bern

Similarly to ANET2\({}_{\text{NORM}}\), we form ANET2\({}_{\text{BERN}}\) by combining the neural network parameter estimation architecture of ANET2 with the quantile regression framework described by [4] as its probability distribution model. In this case, the output of the neural network model is a set of 21 Bernstein polynomial coefficient vectors, each containing 13 values. Therefore, the degree of the Bernstein polynomial we fit is 12 (the degree is one less than the number of parameters). As per the suggestion of [4] we train the model by minimizing the quantile loss on 100 equidistant quantiles. To ensure that no quantile crossing can occur, [4] suggest one might restrict the polynomial coefficient to the positive reals. However, in our testing, this limits the expressive power of the model. Therefore, we let the coefficients be unconstrained as the quantile crossing event is rare in practice [4] and we did not encounter it in our testing (however, one must not neglect this correctness concession as it can be a source of potential issues if left unchecked).

#### 2.2.4 Parameter estimation neural network architecture

To determine the values of these spline parameters (in the case of ANET2\({}_{NORM}\) the parameters of a normal distribution and in the case of ANET2\({}_{BERN}\) the coefficients of the Bernstein polynomial) we use a dense neural network, whose architecture is displayed in Figure 2. ANET2 conducts post-processing jointly for the whole lead time and all spatial locations. Therefore, the input to the network is an ensemble forecast with \(m\) ensemble members, each containing 21 forecasts (we will discuss the training dataset and training procedure in the following section). ANET2 first computes the ensemble forecast mean and variance. These two vectors, each containing 21 elements, are then concatenated with the per-station predictors and seasonal time encodings to form the input to the neural network. The per-station predictors include station and model altitude, longitude, latitude, and land usage. The seasonal time encoding is defined as \(\cos(\frac{2\pi d}{365})\), where \(d\) denotes the day of year the forecast was issued. There are a total of six dense layers in the ANET2 neural network. All dense layers, except for the last one, are followed by a SiLU [20] activation, and each of the layers, except for the first and last, are preceded by a dropout layer [38] with a dropout probability of 0.2. Each of these dense layers is a residual layer implying that the transformation done by the dense layer and the consequent activation is added to the input of that layer, which then forms the output of that residual block (inspect Figure 2 for more details). The final layer produces 21 sets of parameters \(\theta\), representing the parametric distribution for each forecast time.

Figure 2: The neural network architecture of ANET2. The parameter f denotes the number of features in a dense layer. The parameters t and p denote the lead time and the number of per-lead-time parameters required by the normalizing flow. The circular blocks denote the following operations: \(\mathbf{M}\) denotes the computation of the ensemble forecast mean, \(\mathbf{S}\) block denotes the computation of the standard deviation of the ensemble forecast, and \(\mathbf{C}\) denotes a concatenation operation.
### Training procedure

The training procedure is different between ANET1 and the remaining ANET2 variants. The reason for this discrepancy was the fact that ANET1 preceded ANET2 in development. We describe both approaches and give our arguments as to why we modified the training procedure with ANET2.

In both cases, we used the datasets provided by the EUPPBench benchmark [6] to construct the model. The pilot study focused on the post-processing of temperature forecasts for stations in a limited region in Europe. This study included two datasets. The first dataset, denoted as \(D_{11}\), consisted of 20 year _re-forecasts_[6] for the years \(2017,2018\), with 209 time samples for both years. The ensemble for this dataset counted \(m=11\) members. The \(D_{11}\) dataset was the designated training dataset. The second dataset, denoted as \(D_{51}\), included _forecasts_ for \(2017,2018\), with 730 time samples for both years. The ensemble for this dataset counted \(m=51\) members. The \(D_{51}\) dataset was the designated test dataset. We performed post-processing for 229 stations which were equal for both datasets. Each training batch consisted of 256 randomly selected samples across all stations, years, and time samples.

The datasets described above and their designated uses were equal across ANET1 and all ANET2 variants. The main difference between the training procedures of ANET1 and the ANET2 variants lies in how we split the \(D_{11}\) dataset into the training and validation subsets.

#### 2.3.1 ANET1 training procedure

To create the training-validation split we randomly partitioned the \(D_{11}\) dataset across all stations, years, and time samples. The first partition was the final training dataset subset, consisting of 80 percent of \(D_{11}\). The remaining 20 percent formed the validation dataset. We select the model that exhibits the lowest validation loss (the loss function being the negative log-likelihood) as the final candidate model. Additionally, the model did not use the validation or test datasets for parameter estimation during training. This way we select the model that minimizes the unseen data loss in a particular training setup. The random partitioning required us to implement early stopping conditions as otherwise the training quickly lead to overfitting due to the training-validation dataset distribution similarity. Gradient descent was our optimization procedure of choice in conjunction with the Adam optimizer [23] in the Pytorch [31] framework, with its default parameters, a batch size of 128 samples, learning rate of \(10^{-3}\), and weight decay of \(10^{-9}\). If the validation loss did not improve for more than 20 epochs the training was terminated. We reduced the learning rate by a factor of 0.9 after the validation loss plateaued for 10 epochs to increase numerical stability.

Even though ANET1 managed to perform well in the final evaluation against competing methods [6], this training set-up poorly reflected the nature of the train-test relationship. It also lead to stability issues. We rectified this in the next iteration of the training procedure for ANET2.

#### 2.3.2 ANET2 variants training procedure

To rectify the issues of the previous training approach we forwent the random split approach. Therefore, we formed the validation subset of \(D_{11}\) with the re-forecasts corresponding to the year 2016 while we used all the remaining data to train the model. This better reflects the train-test dynamic of \(D_{11}\) and \(D_{51}\), measuring the "out-of-sample" predictive power of the model, as the validation subset is not sampled from the same years as the training one. Similarly to ANET1, we select the model that exhibits the lowest validation loss (ANET2 and ANET2\({}_{\text{NORM}}\) minimize the negative log-likelihood, while ANET2\({}_{\text{BERN}}\) is trained minimizing the quantile loss) as the final candidate model, with not using the validation and test datasets for parameter estimation during training. This way we select the model that minimizes the unseen data loss in a particular training setup however, the training procedure is now more stable and less prone to overfitting. The loss functions were minimized using gradient descent and the Adam optimizer in the Pytorch framework with its default parameters, a batch size of 256 samples, learning rate of \(10^{-3}\), and weight decay of \(10^{-6}\). A bigger weight decay helped with the stability of the learning procedure as ANET2's distribution modeling approach is more flexible compared to ANET1's and therefore required additional regularization. Finally, we also reduced the learning rate by a factor of 0.9 after the validation loss plateaued for 10 epochs to increase numerical stability.

## 3 Results

In this Section we present the evaluation of ANET1, ANET2, ANET2\({}_{\text{NORM}}\), ANET2\({}_{\text{BERN}}\) on the \(D_{51}\) test dataset. We evaluate the performance of the above methods using the continuous ranked probability score (CRPS) [13], bias, quantile loss (QL), and quantile skill score (QSS) [4].
We forwent the direct comparison of ANET variants against conventional post-processing approaches. ANET1 was already evaluated in detail against other methods in [6]. The conclusion from that comparison was that, while most post-processing methods performed similarly, ANET1 achieved the lowest CRPS and further attenuated the error variability in the day-night cycle. ANET1's advantage compared to other approaches in terms of CRPS is most prominent at high-altitude stations (stations with altitudes bigger than 1000 meters). For more details, about this comparison and included method, please refer to [6].

Our evaluation of the ANET variants is displayed in Figure 3 and Table 1. In the case of QSS, we use ANET1 as the reference model, quantifying the performance gain of ANET2 variants relative to ANET1. We compute the CRPS and bias by averaging them across all stations and time samples, resulting in a per-lead-time performance report. The performance values in Table 1 are a result of averaging the corresponding metrics across the entire test dataset.

### Impact of improved network architecture and training procedure

To quantify the contributions of the ANET2 neural network architecture and training procedure, we look at the performance differences between ANET1 and the ANET2 variants. All ANET2 variants perform better on average than ANET1 across all lead times in terms of CRPS and QSS. The difference between individual ANET2 variants is smaller, with ANET2 leading the pack. We notice that ANET2\({}_{\text{NORM}}\) constantly outperforms ANET1 both in terms of average CRPS and QSS, which is further corroborated by the results in Table 1. However, both predict the same target distribution, that being the normal distribution. These performance enhancements of ANET2\({}_{\text{NORM}}\) can be attributed to the architecture and training protocol changes from ANET1. ANET2\({}_{\text{NORM}}\) also exhibits the lowest bias amongst all methods. This could be due to the symmetric nature of the target normal distribution ANET2\({}_{\text{NORM}}\). Conversely, ANET1 outputs the same distribution as ANET2\({}_{\text{NORM}}\) but has a higher bias. We believe that this difference in bias is due to the less stable and less aligned training procedure of ANET1, which results in a sub-optimal convergence.

### Forecast calibration

To further quantify the probabilistic calibration of individual methods we use rank histograms [19], displayed in Figure 3. The results are aggregated across all stations, time samples, and lead times, for each quantile. This gives us an estimate of how well the individual methods describe the distribution of the entire \(D_{51}\) dataset. ANET1 and ANET2\({}_{\text{NORM}}\) exhibit similar rank histograms due to both using a normal distribution as their target, implying similar deficiencies. The standout models are ANET2, which uses the modified normalizing flow approach, and ANET2\({}_{\text{BERN}}\), which is based on Bernstein polynomial quantile regression. Both methods are much more uniform compared to the remaining alternatives. For example, ANET2\({}_{\text{NORM}}\) and ANET1 seem to suffer from over-dispersion as is implied by the hump between the 20-th and 40-th quantiles. We can observe similar phenomena with other post-processing methods that fit a normal distribution to temperature forecasts [6]. ANET2 and ANET2\({}_{\text{BERN}}\) do not exhibit this specific over-dispersion, albeit with a small over-dispersion between the 10-th and 30-th quantile. Additionally, both methods under-predict upper quantiles, albeit to a lesser degree than ANET1 and ANET2\({}_{\text{NORM}}\). We also found that the ANET2\({}_{\text{BERN}}\) method exhibits artefacts that are most pronounced at the low and high quantiles. We believe that this is due to the nature of the Bernstein quantile regression method where the fixed degree of the polynomial affects its expressive value, especially on the quantile edges. Overall, ANET2 displays the most uniform rank histogram out of all evaluated methods.

### Per-altitude and per-station performance

Station altitude correlates with the decrease in model predictive performance, as we show in Figure 5. Therefore, we investigate the altitude-related performance of each model in terms of QSS over three station altitude intervals: \((-5,800]\), \((800,2000]\), and \((2000,3600)\), displayed in Figure 4. We compute the per-altitude QSS by aggregating stations whose altitudes fall into a specific interval. Then we average the QL across all lead times, time samples, and corresponding stations, producing the QL on a per-quantile level. We use ANET1 as the reference model for the QSS. Looking at Figure 4, we can see that all ANET2 variants improve upon ANET1. The only exceptions are ANET2\({}_{BERN}\) and ANET2\({}_{NORM}\) that exhibit roughly equal QL to ANET1 for certain quantile levels at the altitude intervals of \((-5,800]\) and \((2000,3600)\) respectively. ANET2 exhibiting the highest QSS of all methods. Indeed, ANET2 outperforms all the remaining methods over all quantile levels and all altitudes, with the exception of lower quantile levels at the stations with the lowest altitude where ANET2\({}_{\text{NORM}}\) performs the best. The closest method, exhibiting similar trends in improvement to ANET2, is ANET2\({}_{\text{BERN}}\). In our opinion this and the results we show
in Figure 3 further bolster the value of flexible distribution estimators. Similarly, by observing the average CRPS ranking per station, shown in Figure 5, we can see that ANET2 exhibits the best CRPS across 204 stations from a total of 229 stations. ANET2\({}_{\text{BERN}}\) is the most performant method at 18 stations, with ANET2\({}_{\text{NORM}}\) leading in CRPS only at 7 locations. ANET1 is not displayed as its CRPS does not outperform any of the ANET2 variants at any of the stations.

## 4 Discussion

The results we presented provide a clear argument for the use of neural networks in conjunction with flexible distribution estimation methods. Even in the case of the initial run of the EUPPBench dataset with a limited number of predictors and one target variable, temperature, neural networks can tangibly increase the state-of-the-art performance in post-processing. The comparative study we conducted between ANET1, the improved ANET2, and its variants, revealed useful guidelines for future method research and development.

\begin{table}
\begin{tabular}{l|l|l|l}  & CRPS & Bias & QL \\ \hline ANET2\({}_{\text{NORM}}\) & 0.940 & **0.038** & 0.373 \\ ANET2\({}_{\text{BERN}}\) & 0.935 & 0.076 & 0.367 \\ ANET2 & **0.923** & 0.069 & **0.363** \\ ANET1 & 0.988 & 0.092 & 0.386 \\ \end{tabular}
\end{table}
Table 1: CRPS, bias, and QL (quantile loss) for all evaluated methods averaged over all stations, time samples, and lead times for the \(D_{51}\) test dataset. The values in bold represent the best results for each metric.

Figure 3: (Top row): The CRPS, bias, and QSS for all ANET variants. A lower CRPS is better, a bias value closer to zero is better, and a higher QSS is preferable. (Bottom row): The rank histogram of all ANET variants. A more uniform histogram (column heights closer to the black dashed line) is better. Histograms with central humps imply over-dispersion while histograms with outliers on the edges represent under-dispersion.
### ANET1 versus ANET2

First, our results corroborate the findings and practices outlined by [34]. ANET1 uses a per-ensemble-member dynamic attention mechanism with the idea of determining the individual ensemble member importance, conditioned on the predictors and weather situation. While the idea is sound in our opinion, further testing with ANET2 revealed that performing regression on the ensemble mean and variance resulted in equal or even slightly better performance, with a similar number of model parameters. This result is similar to many existing post-processing approaches [34]. We also evaluated the effects of including additional forecast statistics into the input, such as the minimal, maximal, and median temperature. We found that these do not impact the model's predictive performance. We conclude that for the EUPPBench dataset v1.0 the ensemble statistics such as the mean and variance contain enough information for the formation of skillful probabilistic forecasts of temperature. However, it likely that we have not yet found an efficient way of extracting information from individual ensemble members.

### ANET2 versus ANET2\({}_{\text{NORM}}\) and ANET2\({}_{\text{Bern}}\)

ANET2\({}_{\text{NORM}}\) operates by assuming the target weather distribution. However, this method is outperformed by the remaining two, more flexible methods when observing the mean QSS and CRPS. Similar observations were found by other researchers [36, 41] further underlining the need for flexible distribution estimators. From the two flexible approaches we implemented, ANET2 outperforms ANET2\({}_{\text{BERN}}\) in all metrics. The ability to model the probability density through a cascade of spline transformations seems to offer greater flexibility compared to the Bernstein polynomial approach. Additionally, ANET2 does not suffer from the quantile crossing issue that can occur with Bernstein quantile regression [4]. ANET2 produces a probabilistic model from which we can estimate the exact density, distribution, and samples from the distribution.

### Joint lead time forecasting

Where we believe that ANET2 and ANET1 innovate compared to previously applied neural network approaches is in the post-processed forecast lead time. All ANET variants perform predictions for the entire lead time. They are single models for all lead times and all stations, taking as input the entire ensemble forecast (whole lead time). This enables ANET to model dependencies between different forecast times. These ideas are further corroborated by investigating the feature importance generated using input feature permutation [12]. We investigate how different lead times in the input help form the post-processed output for other lead times.

We display the input variable importance relative to a target lead time in Figure 6. Each row signifies a specific lead time in the post-processing output while a column denotes a specific lead time in the ensemble forecast that forms the model input. In many cases, off diagonal elements of each row exhibit high importance implying that ensemble forecasts for different lead times contribute to the correction for a specific target lead time. For example, when ANET2 predicts the distribution for the first target lead time (row denoted with zero) it heavily relies on the first three input forecast lead times (the first three columns). We can also observe that past forecasts exhibit higher importance compared to future forecasts, relative to a specific target lead time. This is implied by the darker

Figure 4: QSS for all compared methods relative to stations with different altitudes. ANET1 is used as the reference model, therefore, the \(y\)-axis denotes the percentage improvement compared to ANET1.
lower triangle of the importance matrix (area under the blue dashed line). Likewise, we can identify an interesting, periodic trend when observing the columns corresponding to lead times at noon We can see that these forecasts non-trivially impact post-processing any target lead time. We hypothesise that this is due to the daily temperature, which frequently reaches its maximum value at noon (this is due to the six-hour resolution), playing an important role in the amount of heat retained throughout the night (void of all other predictors). Since the first noon forecast is the most accurate compared to the rest it acts as an estimator, with consecutive noon forecasts decreasing in importance due to forecast errors.

Additionally, looking at the CRPS comparison of ANET1 with other state-of-the-art post-processing techniques evaluate on the EUPPBench dataset (page 17, Figure 3 ion [6]) we can observe that ANET1 exhibits better correction for temperatures at midnight. The remaining methods suffer from stronger periodic spikes in error when post-processing those lead times. However, the majority of the implemented methods in [6] (except for reliance calibration) operate on a single lead time. If we turn our attention again to Figure 6 and look at each labeled row (corresponding to midnight forecasts), we can see that the importance is more spread out between multiple input lead times. Therefore, it seems that information important for post-processing forecast at midnight is not solely concentrated at that time but is spread out. This might explain why our method can better suppress forecast errors at these lead times.

### Modified derivative estimator

In this section we present our modification to the derivative estimator of ANET2's distribution model which is based on normalizing spline flows [10]. The implementation described by [10] requires the estimation of three sets of parameters for the spline transformations: spline knots, spline values at the knots, and spline derivatives at the knots. [10] suggest that these parameters be estimated by a neural network. However, in our testing, this resulted in sharp discontinuities in the probability density of the final estimated distribution, as can be seen in Figure 7. This is because the values of the derivatives are selected independently of the knots and values, resulting in potential discontinuities. [10] state that this induces multi-modality into the density. However, this might hamper future efforts for determining the most likely predicted weather outcome as sharp jumps in the density could introduce noise around the modes. To produce smoother densities, we modified the derivative estimation procedure such that it now entirely depends on the spline knots and values. This modified derivative estimator was described

Figure 5: (Left): Average per-station CRPS relative to the absolute difference between the model and station altitude. We used median filtering with a kernel size of 15 to suppress outliers. An increase in CRPS correlates with an increase in the absolute difference in altitude. (Right): Average per-station CRPS. The number of cases in which a specific method outperformed the remaining is denoted in the legend. ANET1 is not displayed since it never performed better versus the ANET2 variants.
by [16] and offers certain spline continuity guarantees. The resulting probability density is much smoother (right panel in Figure 7) and in our tests, the modified distribution estimator performs equally well relative to the default derivative estimator in terms of predictive power in our test scenario. Another positive attribute of the modified derivative estimation approach is the fact that it reduced the required number of parameters our neural network has to estimate by a third. This is because in our case the derivatives are estimated from the existing spline knots and values, while the default implementation required the explicit estimation of all three sets of parameters by the neural network. For example, ANET2 issues a probabilistic forecast for each lead time. If we use a normalizing spline flow with \(s\) consecutive spline transformations, each containing \(k\) knots, this results in \(21\times(s\times k\times 3)\) parameters for the default derivative estimation implementation (21 for each lead time, times \(s\) for each spline, times \(k\) for the number of knots, times 3 for each set of parameters - knots, values, derivatives). Contrary to that, the modified derivative approach requires \(21\times(s\times k\times 2)\) parameters, as the derivatives are estimated from the remaining parameters.

### ANET2 drawbacks

The enhanced distributional flexibility of ANET2 comes with two drawbacks: an increased number of parameters in the distribution model, and increased execution time. The execution time and distribution model parameter count comparison is displayed in Table 2. We can see that ANET2 is the slowest of all three variants, exhibiting a roughly five times slower inference time for the \(D_{51}\) dataset compared to ANET2BERN. ANET2NORM has the lowest execution speed, followed by ANET2BERN. This difference in time is mainly due to the sequential nature of the normalizing spline flow: input data has to first pass through a cascade of spline transformations before the density can be estimated. Additionally, when a data sample passes through a spline transformation it has to locate the appropriate spline bin or interval which incurs additional computational burdens of searching for that bin ([10] suggest the use of bisection to speed up this search, however, we do not implement this procedure). The execution times are still small in an absolute sense, as all methods require less than a minute to post-process the entire \(D_{51}\) test dataset, which contains 167170 forecasts. Nevertheless, we must keep this in mind as these computational burdens might compound in a more complex setting with multiple forecast variables in a spatial context.

We can also see that ANET2 requires the highest number of parameters for determining the distribution model. Although ANET2 outperforms ANET2BERN in all forecast evaluation metrics, we have reached the area of diminishing returns. ANET2 requires three times as many parameters compared to ANET2BERN to produce these improvements.

Still, even in the face of these drawbacks, we believe that ANET2 is the most compelling of the three compared approaches in the EUPPBench setting due to its superior distribution estimation performance, and exact probability distribution estimation with no potential quantile crossing. Perhaps these would require attention in a more complex forecast post-processing scenario, however, this requires additional research and testing.

Figure 6: Input forecast lead times (\(x\) axis) importance relative to ANET2’s post-processed lead time outputs (\(y\) axis). Values on the diagonal (blue dashed line) represent the importance of each input lead time to its corresponding post-processed lead time output. Example: the column labeled 24 contains the importance measures of the forecast issued at that lead time on all lead times in the post-processed output.
## 5 Conclusion

In this work we introduced the ANET2 forecast post-processing model to combat two main issues plaguing current state-of-the-art techniques: the lack of expressive power in the parameter estimation and distribution models. We addressed the first by developing a novel neural network architecture and training procedure, processing the entire forecast lead time at once. By implementing a modified normalizing flow approach we tackled the second issue, which resulted in the creation of a flexible, mathematically exact distribution model. We demonstrated the effectiveness of ANET2 against the current state-of-the-art method ANET1, showing that ANET2 outperforms ANET1 in all tested metrics.

This is not to say that ANET2 is without flaws. From all the compared methods, ANET2 is the slowest in terms of inference time (about 10 times slower compared to ANET\({}_{\text{2NORM}}\), requiring 36 seconds for the processing of the D\({}_{51}\) dataset). Likewise, the distribution model described by the normalizing flow requires the most parameters out of all evaluated models. While this is not an issue in the context of the initial EUPPBench experiment, these drawbacks have to be kept in mind if one wishes to scale ANET2 to a spatial domain with more than one weather

Figure 7: (Top): ANET2 with the default derivative estimator for the normalizing spline flow, as described by [10]. (Bottom): ANET2 with the derivative estimation modification based on [16]. QF denotes the plot containing the quantile function (inverse cumulative distribution function). The ANET2 variant with the modified derivative estimator on the bottom suppresses sharp non-linearities, resulting in more ”natural” probability densities and a reduced amount of parameters that the model needs to predict. The observations are generated by sampling two normal distributions with equal variance at different mean temperatures.
variable to be post-processed at the same time.

Still, we consider ANET2 to be a significant contribution to the field of forecast post-processing. To the best of our knowledge, ANET2 is the first method that utilizes normalizing flows to construct weather probability distributions based on weather forecasts. When combined with the novel neural network parameter estimation model, ANET2 achieves state-of-the-art post-processing results. Our intention in introducing this new tool for weather forecast post-processing is to empower national weather forecast providers to generate higher quality probabilistic models with increased confidence.

### Data Availability

The data we used in this study is part of the EUPPBench dataset. However, a subset of the dataset pertaining to the Swiss station data is not freely available. For more information about the dataset please refer to [6]. The EUPPBench dataset is available at [7].

### Code Availability

All ANET variants are open source and available on GitHub, with ANET1 accessible at [2], and ANET2 at [28].

## Acknowledgments

We would like to thank EUMETNET for providing the necessary resources, enabling the creation of the EUPPBench benchmark and to our colleagues behind the benchmark, whose dedicated work lead to its realization.

This work was supported by the Slovenian Research Agency (ARRS) research core funding P2-0209 (Jana Faganeli Pucer).
