# A barycenter-based approach for the multi-model ensembling of subseasonal
forecasts

Camille Le Coz  LMD/IPSL, École Polytechnique, Institut Polytechnique de
Paris, ENS, PSL Research University, Sorbonne Université, CNRS, Palaiseau
France  Alexis Tantet  LMD/IPSL, École Polytechnique, Institut Polytechnique
de Paris, ENS, PSL Research University, Sorbonne Université, CNRS, Palaiseau
France  Rémi Flamary  CMAP, Ecole Polytechnique, Institut Polytechnique de
Paris, Palaiseau, France  Riwal Plougonven  LMD/IPSL, École Polytechnique,
Institut Polytechnique de Paris, ENS, PSL Research University, Sorbonne
Université, CNRS, Palaiseau France

###### Abstract

Ensemble forecasts and their combination are explored from the perspective of
a probability space. Manipulating ensemble forecasts as discrete probability
distributions, multi-model ensembles (MMEs) are reformulated as barycenters of
these distributions. Barycenters are defined with respect to a given distance.
The barycenter with respect to the L2-distance is shown to be equivalent to
the pooling method. Then, the barycenter-based approach is extended to a
different distance with interesting properties in the distribution space: the
Wasserstein distance. Another interesting feature of the barycenter approach
is the possibility to give different weights to the ensembles and so to
naturally build weighted MME.

As a proof of concept, the L2- and the Wasserstein-barycenters are applied to
combine two models from the S2S database, namely the European Centre Medium-
Range Weather Forecasts (ECMWF) and the National Centers for Environmental
Prediction (NCEP) models. The performance of the two (weighted-) MMEs are
evaluated for the prediction of weekly 2m-temperature over Europe for seven
winters. The weights given to the models in the barycenters are optimized with
respect to two metrics, the CRPS and the proportion of skilful forecasts.
These weights have an important impact on the skill of the two barycenter-
based MMEs. Although the ECMWF model has an overall better performance than
NCEP, the barycenter-ensembles are generally able to outperform both. However,
the best MME method, but also the weights, are dependent on the metric. These
results constitute a promising first implementation of this methodology before
moving to combination of more models.

##  1 Introduction

##### Multi-model ensemble methods (MME)

Multi-model ensemble (MME) methods have been shown to improve forecast skill
for different time scales, from short- (Heizenreder et al.,, 2006; Casanova
and Ahrens,, 2009) and medium-range (Hamill,, 2012; Hagedorn et al.,, 2012) to
seasonal forecasting (Palmer et al.,, 2004; Alessandri et al.,, 2011; Kirtman
et al.,, 2014). The added-values of MME over single-model ensemble (SME)
forecasts have been attributed to several factors. First, there is in general
no “best” single-model (Hagedorn et al.,, 2005), the relative performances of
the single-models vary depending on the considered target (i.e. region and
variable of interest, metrics, …). The MME can take advantage of the
complementary skill of the single-models and is able to perform better in
average. Second, Hagedorn et al., (2005) identified error cancellation and
non-linearity of the metric as the main reason for the MME performance being
better than the average performance of the single-models. Third, MMEs allow us
to explore a new dimension of uncertainty which remains unexplored by SMEs.
Indeed, combining different models allows us to explore the uncertainty due to
model formulation: by construction, a SME takes into account uncertainties in
the model initialization and can introduce variations in parameterizations to
sample some of the uncertainty due to parameterizations, but they can not take
into account uncertainty due to model formulation. Besides, Weigel et al.,
(2008) showed that MME can improve the predictive skill if and only if the
single-models ‘fail to capture the full amount of forecast uncertainty’.
Finally, in addition, MME benefit from their larger ensemble size, which
implies a better sampling of the real probability distribution.

##### MME for sub-seasonal to seasonal (S2S) forecasts

Sub-seasonal to seasonal forecasting bridges the gap between weather (medium
range) and seasonal forecast. It corresponds to the time range between two
weeks and up to two months. Predictions at this time scale are the focus of
the Sub-seasonal to Seasonal (S2S) prediction project, whose objectives are to
improve their skill and to promote their use (Vitart et al.,, 2017). As part
of this project, a database containing S2S forecasts from twelve (originally
eleven) operational centers has been made available to the research community.
One of the main research questions the S2S project aims to answer is “what is
the benefit of a multimodel forecast for subseasonal to seasonal prediction
and how can it be constructed and implemented?” (Vitart et al.,, 2017).

Several studies have been investigating the potential benefits of MME for S2S
forecasts (Vigaud et al.,, 2017, 2020; Specq et al.,, 2020; Wang et al.,,
2020; Materia et al.,, 2020; Zheng et al.,, 2019; Pegion et al.,, 2019). These
studies use different MME methods, variables and evaluation criteria, but they
all concur that MMEs are generally performing as well or better than SMEs.
Moreover, Vigaud et al., (2017) and Specq et al., (2020) suggest that MME
improve not only the skill but also the reliability of the probabilistic
forecasts. However, the studies using pooling method point out that the better
performance of the MME is also related to its larger number of members (Zheng
et al.,, 2019; Specq et al.,, 2020; Wang et al.,, 2020). Among other studies,
Karpechko et al., (2018) investigate a specific Sudden Stratospheric Warming
event with a MME, while Ferrone et al., (2017) evaluates several MME methods,
but neither of them compare the skill of the MME to the ones of the SMEs.
Thus, while the potential benefits of MME for the sub-seasonal scales have
clearly been established, the question of how to best combine ensemble
forecasts from different models has received little attention. The present
study focuses on this question.

The most direct and often used method for multi-model combination is the
“pooling method”. This consists in simply concatenating the ensemble members
from the different models. The members of the new multi-model ensemble can
have the same weights or be given different weights based on the model’s
skills (e.g. Weigel et al., (2008) for seasonal scale and Wanders and Wood,
(2016) for sub-seasonal scale). In the above mentioned studies, all used some
variation of the pooling method except Vigaud et al., (2017, 2020) and Ferrone
et al., (2017). They focused on the prediction of terciles, and so couple
multi-model combination with other post-processing methods (e.g. model output
statistics method such as extended logarithmic regression) to predict the
terciles directly. Also focusing on quantiles, Gonzalez et al., (2021) used a
sequential learning algorithm to linearly combine predictors (from the SMEs
but also from the climatology and persistence), their weights being updated at
each step depending on the previous performances. Some other methods for
weighted-MME have been explored, but have not been applied yet to the sub-
seasonal scale. At the seasonal scale, Rajagopalan et al., (2002) developed a
Bayesian methodology to combine ensemble forecasts for categorical predictands
(also used by Robertson et al., (2004) and Weigel et al., (2008) at the
seasonal scale to obtain the terciles of the MME). Other approaches, such as
the Ensemble model output statistics (EMOS) and Bayesian model averaging
(BMA), adopt a probability distribution perspective and aim at building the
PDF of the MME. In the EMOS method, an assumption is made on the shape of the
PDF of the MME and its parameters are then optimized with respect to a chosen
score (e.g. the CRPS) on a training period (Gneiting et al.,, 2005). In
contrast, in the BMA method, an assumption is made on the shape of the input
distributions. The MME’s PDF is then their weighted average, the weights being
equal to the posterior probabilities of the input models (Raftery et al.,,
2005). The following framework of barycenter contains the BMA as will be
discussed later.

##### MME as barycenters of discrete distributions : metrics and weights

In this study, we propose to revisit the combination of multiple model
ensembles from a different point of view. We consider each of the ensemble
forecasts as a discrete probability distribution and reformulate the multi-
model ensemble as a barycenter of those distributions. We show that, in this
framework, the pooled MME is actually the (weighted)-barycenter with respect
to the $L_{2}$-distance. As we work with distributions, instead of a
collection of ensemble members, the notion of barycenter can then be extended
to other metrics in the space of distributions. In particular, a natural
distance in the distribution space is the Wasserstein distance that stems from
the optimal transport theory (Villani,, 2003). The Wasserstein distance is
defined as the cost of the optimal transport between these two distributions.

Optimal transport and the Wasserstein distance have been used in diverse
applications for climate and weather. It has been used to measure the response
of climate attractors to different forcings (Robin et al.,, 2017) and to
evaluate the performance of different climate models (Vissio et al.,, 2020) or
different parametrization (Vissio and Lucarini,, 2018). Moreover, Papayiannis
et al., (2018) use Wasserstein barycenter to point-downscale wind speed from
an atmospheric model. Robin et al., (2019) develop multivariate bias
correction method based on optimal transport, while Ning et al., (2014) use it
in the framework of data assimilation to deal with structural error in
forecasts.

Here, we use the Wasserstein barycenters as a tool to build multi-model
ensembles and compare it to the more traditional $L_{2}$-barycenter (i.e.
pooling method). We investigate the impact of this change of metric on the
MME’s performances. The two barycenters (Wasserstein and $L_{2}$) are applied
to the combination of two models from the S2S database. Focusing on two models
allow us to also explore the importance of the weights given to the models in
the barycenters. We indeed use weighted barycenter to take into account the
single-models performance. The weights are learnt from the data and the two
barycenter-based MMEs are compared for their respective optimal weights.

The paper is organized as follow. The use of the two barycenters as MME
methods is presented in Section 2. We first link the pooling method to the
$L_{2}$-barycenter (2.2.1) before introducing the Wasserstein distance and its
barycenter (2.2.2). The case study is described in Section 3, including the
datasets and the evaluation metrics. The skill of the two MMEs and two SMEs
are evaluated and compared in Section 4. The results are discussed in Section
5, and the main conclusions are highlighted in Section 6.

##  2 Multi-model ensemble methods

###  2.1 Ensemble forecasts as discrete probability distributions

At the S2S time scale, it is necessary to move from a deterministic to a
probabilistic approach using ensemble forecasting (e.g. Kalnay, (2003)). The
ensemble is typically generated by perturbing the initial conditions and
running the model for each of these perturbed states. The set of perturbed
initial conditions represent the initial uncertainty associated with possible
errors in the initial state of the atmosphere. This initial uncertainty is
then transferred in time by the model. Thus, an ensemble forecast is a set of
$N$ perturbed forecasts representing the evolution in time of the probability
density of atmospheric variables according to the model formulation.

An ensemble forecast aims at sampling the probability distribution of the
forecasted variable. Now, this can also be considered as a discrete
probability distribution $\mu$ such that

$\mu=\sum\limits_{i=1}^{N}a_{i}\delta_{\mathbf{x}_{i}}$

where $N$ is the number of members in the ensemble,
$\mathbf{x}_{i}\in\mathbb{R}^{n_{t}}$ is the position of the Dirac
corresponding to the time-series of the $i$-th member, $n_{t}$ is the number
of time steps, and $a_{i}$ is the weight of the $i$-th member (such that
$\forall i\leq n,a_{i}\geq 0$ and $\sum_{i=1}^{N}a_{i}=1$). Thus, here,
$\mathbf{x}_{i}$ does not represent an instantaneous state, but a time series.
In a standard ensemble forecast, all the members are equi-probable, so they
have equal weights $a_{i}=\frac{1}{N}$. In the remainder of this section, we
will look at ensemble forecast time-series from the angle of discrete
distributions.

###  2.2 Barycenters for multi-model combination

The goal of multi-model ensemble methods can be rephrased as combining the
imperfect information from these distributions to obtain a new discrete
distribution representing better the true probability distribution function of
the forecasted variable. A way to summarize a collection of distributions
$(\mu_{1},\mu_{2},...,\mu_{d})$ is to compute their barycenter. The barycenter
is found by solving the following minimization problem

$\displaystyle\operatorname*{arg\,min}_{\mu}\quad\sum_{k=1}^{d}\lambda_{k}.d(\mu_{k},\mu)^{2}$
(1)

where $\lambda_{k}$ represents the weights given to distributions, and
$d(.,.)$ is a distance between distributions. The barycenter, also known as
the Fréchet mean, is effectively the distribution that best represents the
input distributions with respect to a criterion given by the chosen distance
$d$.

As a first step, to both demonstrate the case of feasibility and to
investigate the properties and skills of the different barycenters, we
restrict ourselves to the barycenter of two distributions (or ensemble
forecasts). We use the following notations for the two discrete probability
distributions $\mu_{1}$ and $\mu_{2}$ on the space
$\Omega=\mathbb{R}^{n_{t}}$:

$\displaystyle\mu_{1}=\sum\limits_{i=1}^{N_{1}}a_{i}\delta_{\mathbf{x}_{i}}\quad\text{and}\quad\mu_{2}=\sum\limits_{j=1}^{N_{2}}b_{j}\delta_{\mathbf{y}_{j}}$
(2)

with $X=(\mathbf{x}_{1},\ldots,\mathbf{x}_{N_{1}})\in\Omega^{N_{1}}$,
$Y=(\mathbf{y}_{1},\ldots,\mathbf{y}_{N_{2}})\in\Omega^{N_{2}}$, and where
$\mathbf{a}=(a_{1},\ldots,a_{N_{1}})\in\Sigma_{N_{1}}$ and
$\mathbf{b}=(b_{1},\ldots,b_{N_{2}})\in\Sigma_{N_{2}}$ are probability
vectors. The barycenter of $\mu_{1}$ and $\mu_{2}$ with respect to a distance
$d$ is then given by

$\displaystyle\operatorname*{arg\,min}_{\mu}\quad\alpha.d(\mu_{1},\mu)^{2}+(1-\alpha).d(\mu_{2},\mu)^{2}$
(3)

where $\alpha\in\left[0,1\right]$ represents the weight given to the first
distribution, the second distribution having a weight of $1-\alpha$. These
weights have an important impact on the barycenter and will allow us to take
into account the fact that one distribution has generally better skill. They
can be set from a priori knowledge, e.g. when one model is expected to
represent better the variable than the other one and so it is given more
weight. They can also be learned a posteriori from the data (e.g. from past
forecasts). We choose this second approach. More details on the estimation of
the weights are given in Section 3.3.

####  2.2.1 $L_{2}$ barycenter

The $L_{2}$ distance between two distributions $\mu_{1}$ and $\mu_{2}$ is
given by
$\left\|\mu_{1}-\mu_{2}\right\|_{2}=\left(\int\left(\mu_{1}(z)-\mu_{2}(z)\right)^{2}dz\right)^{1/2}$.
When using this distance in the barycenter equation (3), one can find the
analytical formula for the $L_{2}$-barycenter:

$\displaystyle\mu_{L_{2}}^{\alpha}$
$\displaystyle=\alpha.\mu_{1}+(1-\alpha).\mu_{2}=\alpha\sum_{i=1}^{N_{1}}a_{i}\delta_{\mathbf{x}_{i}}+(1-\alpha)\sum_{j=1}^{N_{2}}b_{j}\delta_{\mathbf{y}_{j}}$
(4)

The first part of Equation (4) shows that the $L_{2}$-barycenter is a weighted
average of distributions. This is similar to the BMA method, which can then be
seen as a $L_{2}$-barycenter with specific weights derived from the Bayesian
framework. Another difference is the assumptions on the shape of the input PDF
made by BMA while we use the ensemble directly as discrete probability
distributions (second part of (4)).

From the second part of Equation (4), one can see that the L2 barycenter
corresponds to the concatenation of the members of the two ensembles with the
re-scaling of the member’s weights using $\alpha$. This can be seen as
“pooling” together the ensemble’s members from the two models. The pooling
method is a simple and well-established MME method. It has been used for the
seasonal scale (Hagedorn et al.,, 2005; Weigel et al.,, 2008; Robertson et
al.,, 2004; Becker et al.,, 2014), for the decadal scale (Smith et al.,, 2013)
and for the medium-range scale (Hamill,, 2012; Hagedorn et al.,, 2012). More
recently, it has also been applied on ensemble sub-seasonal forecasts by
(Karpechko et al.,, 2018; Pegion et al.,, 2019; Zheng et al.,, 2019; Specq et
al.,, 2020; Materia et al.,, 2020). It is interesting to note that the
majority of MME methods for S2S forecasts use pooling and implicitly the
$L_{2}$ distance. However, the $L_{2}$ distance is not the only possible
distance in the distribution space. In the following, we introduce another
distance, the Wasserstein distance, and its associated barycenter.

####  2.2.2 Wasserstein barycenter

The Wasserstein distance stems from the optimal transport theory and can be
seen as the cost of transportation between two distributions $\mu_{1}$ and
$\mu_{2}$. It can be defined on discrete distributions as:

$\displaystyle W_{2}^{2}(\mu_{1},\mu_{2})$ $\displaystyle=\ \min_{T\in
U(\mathbf{a},\mathbf{b})}\langle\mathbf{T},\mathbf{C}\rangle=\min_{T\in
U(\mathbf{a},\mathbf{b})}\sum\limits_{i,j}t_{i,j}\|\mathbf{x}_{i}-\mathbf{y}_{j}\|^{2}$
(5)

where
$U(\mathbf{a},\mathbf{b})=\left\\{\mathbf{T}\in\mathbb{R}_{+}^{N_{1},N_{2}}:\mathbf{T}\mathbf{1}_{N_{2}}=\mathbf{a}~{}~{}\text{and}~{}~{}\mathbf{T}^{T}\mathbf{1}_{N_{1}}=\mathbf{b}\right\\}$
is the set of all the feasible transport matrix $\mathbf{T}$ between the
probability vectors $\mathbf{a}$ and $\mathbf{b}$ (with $\mathbf{1}_{N}$
standing for the all-ones vector of size $N$), and
$\mathbf{C}\in\mathbb{R}^{N_{1}\times N_{2}}$ is a distance matrix whose
elements are the pairwise squared euclidean distances between the elements of
$X$ and $Y$. Here, $\left(\langle\mathbf{T},\mathbf{C}\rangle\right)^{1/2}$ is
the cost associated with the transport $\mathbf{T}$ (with $\langle.,.\rangle$
being the element-wise matrix multiplication operator). The elements $t_{i,j}$
of a transport matrix $\mathbf{T}$ describe the amount of mass going from
$x_{i}$ to $y_{j}$, while $d_{i,j}$ is the cost of moving one mass unit from
$\mathbf{x}_{i}$ to $\mathbf{y}_{j}$. The minimization problem consists in
searching for the optimal transport among all the feasible ones in
$U(\mathbf{a},\mathbf{b})$, i.e. the transport associated with the lowest cost
denoted as the 2-Wasserstein distance. This is a short description of the
Wasserstein distance for discrete distributions, for more information see
Peyré and Cuturi, (2020), Santambrogio, (2015) or Villani, (2003).

The barycenter of two distributions is a special case for which there is a
closed-form expression depending on the optimal plan $\mathbf{T}^{*}$:

$\displaystyle\mu_{W_{2}}^{\alpha}(\mu_{1},\mu_{2})$
$\displaystyle=\sum_{i=1}^{N_{1}}\sum_{j=1}^{N_{2}}t^{*}_{i,j}\delta_{\alpha\mathbf{x}_{i}+(1-\alpha)\mathbf{y}_{j}}$
(6)

where the $\mathbf{T}^{*}=(t^{*}_{i,j})$ is the optimal transport matrix
between $\mu_{1}$ and $\mu_{2}$, i.e. the solution of the minimization problem
(5) (Santambrogio,, 2015, theorem 5.27). For a linear problem on a convex
polytope such as $U(\mathbf{a},\mathbf{b})$, the optimum is always achieved on
a vertex of the feasible region. Thus, the minimum is found on a vertex of the
feasible set $U(\mathbf{a},\mathbf{b})$, which means that the optimal matrix
$\mathbf{T}^{*}$ has at most $N_{1}+N_{2}-1$ non-zeros elements. Thus, the
Wasserstein barycenter $\mu_{W_{2}}$ in (6) has a maximum of $N_{1}+N_{2}-1$
points, or members in terms of ensemble forecast (and not $N_{1}\times N_{2}$
as the formula suggests).

Figure 1: Illustration of the $W_{2}$-barycenter between two 2D discrete
distributions $\mu_{1}$ and $\mu_{2}$. The weights of the points are indicated
by the size of the markers. The lines represent the optimal transport between
$\mu_{1}$ and $\mu_{2}$, with the intensity of the line indicating the mass
transferred to one point to the other. The $W_{2}$-barycenter’s points for
$\alpha=0.5$ are located in the middle of these lines.

An illustration of the Wasserstein barycenter between two distributions is
shown in Figure 1. One can see an example with two 2D discrete distributions,
their optimal transport plan $\mathbf{T}^{*}$ and their Wasserstein
barycenter. In this example, all the points of a distribution have the same
weights. However, the first distribution has more points than the second one
and so its points have a smaller weights individually than the ones of the
second distribution. This is representative of ensemble forecasts: all their
members are usually equi-probable (i.e. they have the same weights) but the
number of members per forecast vary from one model to the other. The optimal
transport between the two input distributions is represented here by lines
between their points. As shown by Equation (6), for $\alpha=0.5$, the points
of the Wasserstein barycenter are located in the center of these lines. The
weights of the points in the barycenter distribution are equal to the mass
($t^{*}_{i,j}$) transported along the lines.

####  2.2.3 Illustration of the application of barycenters to ensemble
forecasts

The $L_{2}$ and $W_{2}$ distances have different properties in the
distribution space, and so lead to different multi-model ensembles. An
interesting property of the $W_{2}$-distance is that it captures the proximity
of two discrete distributions as a whole, whereas the $L_{2}$-distance treats
all the Dirac distributions independently. For example, contrary to the
behavior of the $W_{2}$-distance illustrated in Fig. 1, the $L_{2}$-distance
between two distributions with disjoint supports is equal to the sum of their
$L_{2}$-norms, no matter how distant their supports. These properties will
reflect on the corresponding barycenters.

Figure 2: Illustration of multi-model ensembles using barycenters for a
synthetic variable. The ensemble can be seen as a discrete distribution whose
points are the time series of each member. The weight of the points in the
distribution is indicated here by the thickness of the line.

Figure 2 gives an illustration of the two barycenters applied to synthetic
ensemble forecasts. The $L_{2}$-barycenter in Figure 2b is the concatenation
of the two input ensembles in Figure 2a. The weight of each member of the
$L_{2}$-barycenter is given by the weight of the corresponding member in the
input ensemble ($1/N_{1}$ or $1/N_{2}$) multiplied by the corresponding
model’s weight ($\alpha$ or $1-\alpha$). The $W_{2}$-barycenter shown in Fig.
2c has a different structure compared to the $L_{2}$-barycenter. One major
difference is that none of its members belongs to the input ensembles. In
other words, while the support of the $L_{2}$-barycenter is the union of the
input’s supports, the support of the $W_{2}$-barycenter is on the path of the
optimal transport and so can be different. This also allows the
$W_{2}$-barycenter to retain some characteristics of the input distributions.
For example, in this illustrative case, both input distributions are unimodal
and so is the $W_{2}$-barycenter. This is however not the case for the
$L_{2}$-barycenter which has two modes, each associated with the mode of one
of the input distributions.

These properties of the two barycenters may be attractive for different types
of forecast error: the artificial example of Figure 2 is constructed to
illustrate fundamental differences between the two barycenters. It may be
tempting to interpret the $W_{2}$-barycenter as a mean to correct biases, or
to ponder the advantage of the $L_{2}$-barycenter in preserving bi-modality,
which may be of interest (Bertossa et al.,, 2023). Actual forecasts (Figure 3)
are more complex and an assessment of different barycenters will require
extensive testing with different criteria. The purpose of Figure 2 is to
emphasize the (hitherto largely untapped) variety of ways to build multi-model
ensemble forecasts.

Figure 3 shows a similar example as Figure 2 but for two real S2S forecasts.
In Fig. 3a, one can see the forecasted 2m-temperature in Paris according to
the ECMWF (in blue) and to the NCEP (in orange) models from the S2S database.
Their $L_{2}$\- and $W_{2}$-barycenters are shown in Fig. 3b and 3c
respectively. It is interesting to note that, for a given parameter $\alpha$,
the two barycenters have different ensemble members but their ensemble means
are the same. The difference between them is the way they represent the
forecast uncertainty. This raises the question: how different are the
$L_{2}$\- and $W_{2}$-distributions and which one captures better the forecast
uncertainty?

Figure 3: Same as Fig. 2 but for two real sub-seasonal ensemble forecasts of
2m-temperature over Paris (initialized at 2018-02-01 00:00:00). The reference
2m-temperature is added in black.

##  3 Data and methodology

In this study, we focus on the sub-seasonal forecasting of 2m-temperature
during boreal winter in Europe. The energy demand is higher in winter, and is
particularly dependent on the temperature (due to the use of electrical
heating). Moreover, temperature is a large-scale field for which forecast
models generally have good skillmaking it a good starting point to evaluate
the different MME methods.

###  3.1 Data

To make predictions for this case and to validate these predictions, we need
both a dataset of ensemble forecasts from multiple dynamical models to which
the barycenters will be applied as well as a reference dataset against which
to evaluate the skills of the forecasts.

####  3.1.1 S2S data

For this first implementation, we selected two models from the S2S database
(Vitart et al.,, 2017). The first one is the European Centre for Medium-Range
Weather Forecasts (ECMWF) model, which has been shown to have some skill for
European winters at the sub-seasonal scale (Monhart et al.,, 2018; Büeler et
al.,, 2020; Goutham et al.,, 2022). For the second one, we chose the National
Centers for Environmental Prediction (NCEP) model, since its development was
essentially independent from the one of ECMWF. As explained earlier, one of
the advantage of doing multi-model ensemble is to sample and compensate better
the models error. If the models have similar construction, this error is not
sampled properly. Moreover, both models have a long time range and a large
ensemble size. Their main characteristics are summarized in Table 1.

Table 1: Description of the two models from the S2S database. See Vitart et al., (2017) for more details. | ECMWF | NCEP
---|---|---
Forecasts |  |
Time range |  d 0–46 |  d 0–44
Resolution |  Tco639/Tco319 L137 |  T126 L64
Ensemble size |  50+1 |  15+1
Frequency |  Twice a week
(Monday, Thursday) |  Daily
Reforecasts |  |
Method |  On the fly |  Fixed
Period |  Past 20 years |  1999–2010
Frequency |  Twice a week |  Daily
Ensemble size |  10+1 |  3+1

The 2m-temperature forecasts and reforecasts from both models were retrieved
from the S2S database through the ECMWF’s Meteorological Archival and
Retrieval System (MARS). We retrieve the forecasts of the temperature daily
mean directly on a $$$\times$$$ grid for our study domain shown in Figure 4,
that is Europe ($34^{\circ}N-74^{\circ}N$, $13^{\circ}W-40^{\circ}E$). We
select the forecasts started during the months of December-January-February
(DJF) for the 2015–2022 period, for which ECMWF and NCEP have matching
starting dates. At the end, we obtain a total of 180 simulations . Here, we
are only using the perturbed forecast members of the forecasts to build the
discrete distributions. For consistency, we also discard the control member of
the reforecasts for the calibration.

Figure 4: Study domain with the $1.5^{\circ}\times 1.5^{\circ}$ lat/lon grid
from the forecasts. Only the land points (in grey) are used in this case
study.

Due to model errors, models tend to drift away from the observation climate
toward the model climatology as lead time increases (Takaya,, 2019). It is
thus important to calibrate extended-range forecasts. The model climatology
can be estimated from the reforecasts, which are “retrospective” forecasts run
with the same model than the forecast but in the past, while the observed
climatology is derived from a reference data set (Sect. 3.1.2 below). Then,
calibration methods can be applied to statistically correct the forecasts.
Here, we use the mean and variance adjustment (MVA, Leung et al., (1999);
Manzanas et al., (2019)) method to calibrate the forecasts as in Goutham et
al., (2022).

The reforecast production method differs from center to center. At the ECMWF,
they are produced “on the fly”, that is each forecast is provided with its set
of corresponding reforecasts (initialized for the same day of the year over
the last previous 20 years). On the contrary, NCEP produced the reforecasts at
once for all the days of a fixed period. From this “fixed” set, we select all
the reforecasts initialized the same day of the year than the forecast’s
initialization date. This calibration approach is more favorable to the ECMWF
forecasts since the climatological statistics are computed over a longer
period (20 year and 10 members versus 12 years and 3 members) and their
evolution is taken into account thanks to the rolling average, contrary to the
NCEP forecasts for which a fixed time period is used. However, this approach
was chosen to have a similar construction for both models and in accordance
with the reforecasts availability.

####  3.1.2 Reference data

The Modern-Era Retrospective analysis for Research and Applications, version 2
(MERRA-2, Gelaro et al., (2017)) reanalysis is used as reference for the
calibration and the validation of the forecasts. We use a reanalysis as
reference in order to have a spatially and temporally complete dataset. We
choose the MERRA-2 reanalysis because it is a recent reanalysis covering both
the calibration and validation period, and it is based on a different global
circulation model than both the ECMWF and the NCEP S2S forecasts. This would
not the case, for example, with the ERA5 reanalysis which is also produced by
the ECMWF with a similar model as the S2S forecasts (Hersbach et al.,, 2020).
The 2m-temperature daily means were retrieve from NASA Goddard Earth Sciences
(GES) Data and Information Services Center (DISC) center on MERRA-2’s native
grid, i.e. $None$ lat $\times~{}$$$ lon grid. The data was re-gridded on the
same $$$\times$$$ lat/lon grid as the forecasts using bi-linear interpolation
with the Climate Data Operators (CDO, Schulzweida, (2022)).

We also use the MERRA-2 reanalysis to build a 30-year rolling climatology. The
climatology is a common benchmark for forecast validation, and is used to
compute skill scores (see Section 3.2). For a given day of the year, the
climatology corresponds to the MERRA-2 values for the same day over the
previous 30 years. Using a rolling climatology allows us to take into account
the climatic trend of the temperature. The climatology is thus an ensemble
with 30 members, each one corresponding to a different year.

###  3.2 Skill metrics

The performance of the forecasts are evaluated for weekly averages for weeks 3
to 6, which corresponds to the sub-seasonal time scale. We only evaluate the
forecasts’ performance over the land points of our study domain (that is a
total of 465 grid-points indicated in grey in Fig. 4), and, when spatial
averaging is needed, the scores are weighted with the cosine latitude in order
to account for the spherical geometry of the Earth. We consider several
metrics to evaluate and compare the performances of the different ensemble
forecasts.

####  3.2.1 Continuous Ranked Probability Score (CRPS)

The CRPS is a widely used score for probabilistic forecasts of continuous
variable (Matheson and Winkler,, 1976; Hersbach,, 2000; Wilks,, 2019). The
CRPS is actually the $L_{2}$ distance between the Cumulative Density Function
(CDF) of the ensemble forecast and the CDF of the observation:

$\displaystyle\text{CRPS}_{n}$
$\displaystyle=\int_{-\infty}^{+\infty}\left[F_{\text{fc},n}(y)-F_{\text{ref},n}(y)\right]^{2}dy$
(7)

where $F_{\text{fc}}$ is the CDF of the forecast , $F_{\text{ref}}$ is the CDF
of the observation and $n\in[\\![1,N]\\!]$ is the simulation’s number. The CDF
are computed empirically from the ensembles. In general the reference is
deterministic, and its CDF is then a step function

$F_{\text{ref}}(y)=\begin{cases}0&\text{if }y<x_{\text{ref}},\\\ 1&\text{if
}y\geq x_{\text{ref}}.\end{cases}$

where $x_{\text{ref}}$ is the deterministic value of the reference.

In the reminder, we evaluate the models with respect to the their CRPS
averaged over the simulations. We denote this mean-CRPS by CRPSm (to avoid any
confusion with $\text{CRPS}_{n}$ above). The CRPSm is negatively oriented with
0 being the perfect score.

####  3.2.2 Proportions of skillful forecasts (CRPSp)

The second performance score we consider for the evaluation of the models is
the proportion of skillful forecasts (Goutham et al.,, 2022). It corresponds
to the percentage of simulations for which the model has a better CRPS than
the climatology:

CRPSp
$\displaystyle=\frac{100}{N}\sum_{n=1}^{N}\left[\text{CRPS}_{\text{fc},n}<\text{CRPS}_{\text{clim},n}\right]$
(8)

where $N$ is the number of simulations. The proportion of skillful forecasts
(CRPSp) is positively oriented, and we consider that the model has some skill
if $50\%$ of the forecasts are skillful. It is also more robust to outliers
than the CRPSm because it will not be affected much by a few simulations very
far from the ground truth.

####  3.2.3 Brier Score

In order to investigate different attributes of the ensemble forecasts, we
also compute their Brier score and associated decomposition for the
“2m-temperature being below normal”. The Brier score is a measure of accuracy
for probability forecasts of dichotomous events (similarly to the CRPS for a
continuous predictand, Wilks, (2019)). It is equal to the mean square error of
such probabilistic forecast:

Brier score $\displaystyle=\frac{1}{N}\sum_{n=1}^{N}(y_{n}-o_{n})^{2}$ (9)

where $N$ is the number of simulations (here $N=180$ starting-dates), $y_{n}$
is the forecasted probability of the event and $o_{n}$ the observed event
(i.e. 1 if it occurs, 0 otherwise). If the probability forecasts $y_{n}$ are
only allowed to take a finite number of values (here
$\left[0,0.2,0.4,0.6,0.8,1.0\right]$), it can be decomposed in three terms
including the reliability and the resolution, which quantify different
attributes of the forecasts (Murphy,, 1973; Wilks,, 2019). The reliability
characterizes the correspondence between the predicted probability of the
event with respect to the ensemble forecast and the relative event’s frequency
conditioned on the forecasted probability (and is positively oriented). The
resolution describes whether different predictions lead to different outcomes
(and is negatively oriented).

The ensemble forecast for the temperature is transformed into probabilistic
forecasts $y$ for this dichotomous event by counting the proportion of
ensemble members that are in the lower tercile of the climatology and rounding
it to the nearest allowed values.

####  3.2.4 Relationships between skill metrics and barycenters

  * •

The Wasserstein distance could also be used as a score for forecast
validation. However, in the case of a deterministic observation, the
Wasserstein distance is equivalent to the RMSE (see appendix A). This means
that it does not take into account the uncertainty information of the
ensemble. Thus, we do not use it for the evaluation of the forecasts in this
study.

  * •

It is interesting to note that for univariate distribution, the
$L_{2}$-barycenter is also the barycenter with respect to the energy distance
(see appendix B). Thus, since the CRPS is identical to the energy distance in
1D (see Wilks,, 2019), we can expect to have good CRPS performance for the
$L_{2}$-barycenter.

  * •

One can show that the CRPS of the $L_{2}$-barycenter
$\mu_{L_{2}}^{\alpha}(\mu_{1},\mu_{2})$ can be expressed as a function of the
CRPS of $\mu_{1}$ and $\mu_{2}$:

$\displaystyle\text{CRPS}(\mu_{L2},\nu_{obs})$
$\displaystyle=\alpha.\text{CRPS}(\mu_{1},\nu_{obs})+(1-\alpha).\text{CRPS}(\mu_{2},\nu_{obs})-\alpha.(1-\alpha).\text{CRPS}(\mu_{1},\mu_{2})$
(10)

where $\nu_{obs}$ is the observation (see appendix C for proof).

###  3.3 Model combinations with barycenters

####  3.3.1 Estimation of the weight $\alpha$

The ECMWF and NCEP ensembles are combined following the barycenter-based MME
method described in Section 2. We consider two multi-model ensembles, one
based on the $L_{2}$-barycenter and one based on the $W_{2}$-barycenter. The
weights given to the models into the new multi-model ensemble are represented
by the parameter $\alpha\in[0,1]$ (see Eq. (3)). That is, the weight $\alpha$
is given to the ECMWF ensemble and the weight $1-\alpha$ is given to the NCEP
ensemble. We estimate this parameter $\alpha$ from the data.

We perform cross-validation to validate the multi-model forecasts. To account
for the non-stationarity associated with the seasonal cycle, the 180
simulations are divided in seven folds corresponding to seven winters (with 25
or 26 simulations each). Simulations from six winters are used to determine
the optimal value of $\alpha$, then this value is used to compute the
barycenter ensembles for the last winter. To find the optimal $\alpha$, we use
a grid search in the interval $[0,1]$ with a step of 0.02.

##### Remark:

We mentioned earlier that the BMA method is a particular case of the
$L_{2}$-barycenter. In the case of the BMA, the weights are the posterior
distributions of the input models. They are representing the probability that
the associated model is the best model.

####  3.3.2 Optimal $\alpha$ on the data

In fact, for each barycenter we estimate two optimal $\alpha$ values, one for
the mean CRPSm and one for the mean CRPSp averaged over weeks 3 to 6 and the
whole domain. These optimal $\alpha$ values are shown in Figure 5. The narrow
spreads of the distributions of the optimal $\alpha$ for a given metric and
barycenter indicate that the weights are consistent accross the different
folds of the cross-validation. Our training period of six winters seems to be
sufficient to derive stable weights.

One can see that for both barycenters and both scores, the optimal $\alpha$
values are above $0.5$. That means that more weight is given to the ECMWF
ensemble, which is in agreement with the known skill of the models. In the
next section, we will see that the ECMWF ensemble indeed tends to show better
performance than the NCEP ensemble for all the scores. For the
$L_{2}$-barycenter, both scores are associated with similar optimal values for
$\alpha$, around 0.8. For the $W_{2}$ barycenter, the scores tend to be more
sensitive to the value of $\alpha$. The optimal value of $\alpha$ is around
0.85 for the CRPS, and around 0.7 for the CRPSp.

Figure 5: Distributions of the optimal parameter $\alpha$ obtained by cross-
validation for (a) the $L_{2}$-barycenter and (b) the $W_{2}$-barycenter. The
optimal $\alpha$ across the seven folds are represented as boxplots and their
mean is represented by a square marker.

##  4 Results

###  4.1 Models evaluation

In this section, we validate and compare the weekly 2m-temperature from the
two multi-model (barycenter) and the two single-model forecasts. All the
scores shown here for the former are obtained with their respective optimal
weight $\alpha$.

Figure 6: Mean of (a) the spatial average of the CRPSm and (b) the spatial
proportion of skillful forecasts CRPSp of the weekly 2m-temperature for the
different models and as a function of the lead-time.

##### Spatial performances

Figure 6 shows an inter-comparison of the CRPSm and CRPSp of the different
ensembles. In the first panel 6a, one can see the average of the CRPSm over
the domain. The distribution over the starting dates is represented by its
mean in the plot, and is used to test if the model’s scores are significantly
different from each other at a $5\%$ significance level. The significance is
inferred using the Wilcoxon signed-rank test, a non-parametric paired
statistical test (Wilcoxon,, 1945). A first surprising observation is that the
CRPS of the climatology decreases with lead-times. This can be explained by
the effect of the seasonality on the performances of the climatology. Indeed,
the CRPS of the climatology has a strong annual cycle that goes from about
$0.9^{\circ}C$ in summer to $1.5^{\circ}C$ in winter (values obtained for the
period 2016-2022 and averaged over Europe). Due to our choice of forecast’s
selection (i.e. forecasts initialized in DJF), the six weeks lead-time
contains more weeks of spring while the four weeks lead-time contains mostly
weeks of winter. We also observe that the NCEP ensemble has a significantly
lower performance than the other ensemble forecasts, and even than the
climatology in average. The three other ensembles have very similar mean
CRPSm, with the $L_{2}$-barycenter being slightly better followed by the
$W_{2}$-barycenter. However, the skills of the two barycenter ensembles are
not significantly different from that of the ECMWF ensemble at the $5\%$
significance level, except at week 5 (see Table 2).

Similarly, panel 6b shows the mean of the (spatially averaged) CRPSp
distribution over the starting-dates. It is interesting to note that, despite
being worse than the climatology on average, NCEP is better than the
climatology in terms of proportion for weeks 3 to 5. Indeed, its CRPSp is
above $50\%$, meaning that NCEP performs better than the climatology in terms
of CRPS at more than half the grid-points). However, the three other ensembles
perform significantly better than the NCEP one. This time, the
$W_{2}$-barycenter ensemble has a significantly better proportion of skillful
forecasts than both ECMWF and $L_{2}$-barycenter. The $L_{2}$-barycenter also
performs better than ECMWF, but the difference is significant at the $5\%$
level only for week 4 (see Table 3).

##### Spatial distribution of the best models

Figure 7: Best-model maps with respect to the CRPSm of the 2m-temperature for
different lead-times. The color of the pixels indicates which of the ensemble
forecasts has the best score at this location. The crosses indicate grid-
points for which the best model performs significantly better than the ECMWF
ensemble at the $5\%$ significance level (according to the Wilcoxon signed-
rank test).

The performance of the ensembles also varies across the domain. The spatial
patterns are similar for all four ensembles, in the sense that variations
across the domain for each model tend to dominate over differences between
models (not shown here). An exception is the NCEP ensemble that shows clearly
worse performance than for the three other ensembles. Thus, for easier
comparison, we build best model maps that show for each grid-point which of
the four ensembles perform best with respect to the chosen score. The best-
model maps with respect to the CRPSm for the different weeks is shown in
Figure 7. The best model varies across the domain. There are no grid-points
for which the NCEP ensemble performs significantly better than the ECMWF
ensemble. However, in agreement with the average score from Figure 6a, the
$L_{2}$-barycenter has the best CRPSm at a majority of pixels, but is
significantly different from the ECMWF ensemble only in specific areas, mostly
in the North-East of Europe. This shows that, around this area, the
information provided by the NCEP ensemble to the $L_{2}$-barycenter is helpful
to forecasting, even though the performance of the NCEP ensemble alone tends
to be worse than that of the ECMWF ensemble.

Figure 8: Best model maps with respect to the proportion of skillful forecast
CRPSp of the 2m-temperature for different lead-times. The color of the pixels
indicate which of the ensemble forecasts has the best score at this location.
The crosses indicate grid-points for which the best model performs
significantly better than the ECMWF ensemble at the $5\%$ significance level
(according to the McNemar test).

Similarly, best-model maps with respect to the proportion of skillful
forecasts (CRPSp) are shown in Figure 8. Since this time we have a
distribution of nominal values (true when the forecast is better than
climatology, false otherwise), we use the McNemar test to test whether the
best model is significantly better than the ECMWF ensemble (McNemar,, 1947).
The best-model maps for the CRPSp tend to be less smooth than for the CRPSm
(in terms of spatial correlation). We can observe that the $W_{2}$-barycenter
performs best at a majority of grid-points, but also that the significance of
this result is limited to some areas, in northern and eastern Europe,
depending on the lead-time. There are also a few grid-points where the
$L_{2}$-barycenter performs better than the ECMWF ensemble. Thus, like for the
CRPSm (Fig. 7), a barycenter performs significantly better than the ECMWF
ensemble in northern Europe, but the barycenter performing best for the CRPSp
is the $W_{2}$ rather the $L_{2}$ one.

###  4.2 Forecast attributes

So far, the choice of scores has put the focus on the accuracy of the
forecasts, without distinguishing the different attributes responsible for
that accuracy (Wilks,, 2019). Moreover, the focus was on forecasting the the
whole distributions. However, in applications, one may be interested in
predicting statistics about specific parts of the distributions, e.g. cold
spells for the energy sector. To this end, we now focus on predicting the
probability that the temperature will be in the lower tercile of the
distribution compared to the climatology. In particular, we want to
investigate the impact of the weight $\alpha$ on different attributes of such
forecasts. In order to do that, we use the Brier score and its decomposition
into resolution and reliability.

Figure 9: Reliability-resolution diagrams showing the evolution of the
forecasts’ attributes for the Brier score of the lower tercile of the
temperature with respect to the barycentric weight $\alpha$ for (a) week 3 and
(b) week 4. The optimal values of $\alpha$ with respect to the different
attributes are indicated in the inserts. The dashed lines represent isolines
of the Brier score, with largest values of the Brier score to the top left of
the plots.

Figure 9 shows the reliability and resolution for weeks 3 and 4 of the
different ensembles: the ECMWF and NCEP ensembles, and the two barycenter
ensembles for values of $\alpha$ varying from 0 (equal to NCEP) to 1 (equal to
ECMWF). This allows to trace the variation of the Brier score and of the
forecasts’ attributes as more weight is given to one or the other initial
ensemble. For both lead-times, the ECMWF ensemble has a better accuracy (i.e.
better Brier score) but also better reliability and resolution than the NCEP
ensemble. While the starting and end points of the two barycenter curves for
varying $\alpha$ are the same, the paths followed differ. For both
barycenters, the progression of the two attributes is relatively smooth
between $\alpha=0$ and $\alpha=0.6$, but becomes more nonlinear as more weight
is given to the ECMWF ensemble. For $\alpha<0.6$, the $L_{2}$-barycenter has a
better reliability than the $W_{2}$-barycenter for a given resolution.
Conversely, the $W_{2}$-barycenter has a better resolution for a given
reliability. However, the overall Brier score is larger for both barycenters
than for the ECMWF ensemble. To the contrary, for $\alpha>0.6$, the Brier
scores of the barycenters become larger than the ECMWF scores. This is mostly
due to the better resolution of both barycenters for week 3, while both
reliability and resolution are improved for week 4. Overall, for a given value
of $\alpha$, the $L_{2}$-barycenter has a better Brier score and reliability
and the $W_{2}$-barycenter has a better resolution (not shown here). It is
interesting to note that the best reliability and the best resolution are not
reached for the same value of $\alpha$ and that the optimal values of $\alpha$
depend on the choice of barycenter. It also does not necessarily correspond to
the best Brier score.

##  5 Discussion

###  5.1 Performance of the MME

We evaluated the two (barycenter-based) multi-model and the two single-model
ensembles of 2m-temperature with respect to different metrics: the mean CRPS
(CRPSm), the proportion of skillful forecasts (CRPSp), as well as their Brier
score decomposition into resolution and reliability. In general, the multi-
model ensembles improve forecast skill compared to the single-model ones. On
average, the $W_{2}$-barycenter has a significantly better CRPSp. It is also
the best barycenter for this score at most of the locations where a barycenter
performs significantly better than ECMWF. On the other hand, the
$L_{2}$-barycenter has a better average CRPSm and tends to be the best
barycenter at locations where the results are significant. In other words, the
$L_{2}$-barycenter is better in average but the $W_{2}$-barycenter is better
more often with respect to the CRPS. However, the best model also depends on
the location. In particular, the ECMWF ensemble outperforms the other
ensembles over (most of) the Iberian peninsula for both metrics. Regarding the
prediction of the lowest tercile of the temperature distribution, it is also
the $L_{2}$-barycenter that performs best in terms of accuracy (i.e. Brier
score). The decomposition of the Brier score into reliability and resolution
allows us to investigate how these attributes contribute to the accuracy of
the ensemble forecasts. The $L_{2}$-barycenter has an overall better
reliability than the $W_{2}$-barycenter, but the $W_{2}$-barycenter tends to
have a better resolution.

There can be several reasons why merging the ECMWF ensemble with the less
skillful NCEP ensemble, is leading to barycenters with improved performances.
First, NCEP has generally less skill than ECMWF, but can be punctually better
for some given locations, lead times or initialization dates. The barycenters
can exploit this information to improve skill thanks to error cancellation and
to the non-linearity of the skill metrics Hagedorn et al., (2005). Second, the
ECMWF ensemble may have good performances but be overconfident. In that case,
adding a model with lower skills increases the spread of the ensemble and can
move the ensemble mean towards the truth, as shown by Weigel et al., (2008)
for seasonal ensemble forecasts. This can be seen from Equation (10) in the
case of the $L_{2}$-barycenter’s CRPS. This equation shows that the CRPS of
the $L_{2}$-barycenter is composed of two parts: the weighted average of the
CRPS of the single-model ensemble and the (weighted) CRPS between them. Thus,
if one model has a worse CRPS than the other, it can still improve the CRPS of
the barycenter if the CRPS between the models is large enough to compensate
its CRPS.

An example of an interesting case for which NCEP performed best is the cold-
wave event that occurred in February 2018 in France. Figure 3 shows the ECMWF
and NCEP forecasts initialized the 2018-02-01 (ensembles 1 and 2 respectively
in the legend) as well as the daily 2m-temperature according to MERRA-2
(reference in black). One can observe that ECMWF seems to miss the cold-wave
in week 4. On the other hand, some members of NCEP do predict an important
temperature decrease. The peak of the temperature drop is shifted by one or
two days in NCEP but remains within the same week, which is the the typical
resolution for sub-seasonal timescale. Moreover, the large spread of NCEP’s
members during that week translate well the large uncertainty of the forecast.

###  5.2 Importance of the model’s weights

Equal weighting of the models in the pooling method is the simpler and most
used approach. However, a few studies investigate the use weighted multi-model
ensembles with divergent results (at weather or seasonal scale in Weigel et
al., (2008); Casanova and Ahrens, (2009); Kharin and Zwiers, (2002) and
climate in Haughton et al., (2015)). At the sub-seasonal scale, Wanders and
Wood, (2016) used muti-variable linear regression on the ensemble means to
derive the model weights. They show that the weighted muti-model ensemble have
a better deterministic performance but also probabilistic one (in terms of
Brier score) than the non-weighted multi-model. In this paper, we also derive
the weights on previous model’s performances but using different criteria: the
CRPSm and CRPSp (both related to the CRPS). Thus, we optimize the weight
taking into account the whole distribution instead of focusing on the ensemble
mean. In agreement with Wanders and Wood, (2016), we observe that the weights
have a large impact on the performance of the barycenter-ensemble. In fact,
the superior skill of the barycenters were obtained for an optimal weighting
learned from past data. This result does not hold for all values of the
weight. The two barycenters do outperform the NCEP ensemble for all weights,
but the ECMWF ensemble outperforms them for low values of the parameter
$\alpha$ (equal to the weight on ECMWF). It is coherent that, due to the lower
skill of the NCEP ensemble compared to the ECMWF ensemble, more weight should
be put on the latter. This difference in performance is also reflected in the
optimal model’s weights, with the optimal values of $\alpha$ for both
barycenters and metrics being above $0.5$.

The optimal weight depends on the barycenter method and the metric. The
optimal values of the weight are similar for both domain-average metrics for
the $L_{2}$-barycenter (around $0.8$), while they are more sensitive to the
metric for the $W_{2}$-barycenter (around $0.9$ for the CRPSm and $0.7$ for
the CRPSp). Moreover, in the case of the Brier score for the lowest tercile of
the temperature and its decomposition, the reliability and the resolution do
not have the same optimal value for the weight. Thus, the choice of the weight
has to be done in function of the targeted application.

In this work, we assume that there is no best model and that the models are
complementary as they sample different part of the forecast uncertainty. Our
aim is thus to find the best combination of models, and the weights represent
the contribution of each model in this combination. This is different from the
BMA method which assumes that there is a best model but is uncertain about
which one (Raftery et al.,, 2005). The BMA is a weighted average of
distributions, similar to the $L_{2}$-barycenter, but its weights are
representing the probability of the associated model to be the best one.
However, it would be possible to combine our barycenter approach with a
Bayesian framework but with uncertainty on the weights (instead of on the
models as in BMA). We could in principle assign priors to the weights of our
barycenters and follow a Bayesian approach to derive the corresponding
posterior distributions, instead of using (estimated) optimal deterministic
weights. To our knowledge, however, such an approach has not been developed
yet and it remains to be shown whether or not it would be feasible
analytically, or at least tractable computationally.

##  6 Conclusion

We have explored methods to combine ensemble forecasts from multiple models
based on barycenters of ensembles. Building on the recognition of the
relevance of probabilistic forecasts for S2S prediction, we work directly in
the probability distribution space. That is, the ensemble forecasts are
manipulated as discrete probability distributions. This allow us to use
existing tools from this space and in particular the notion of barycenter. The
barycenter of distributions is the probability distribution that best
represents the collection of input distributions (with respect to a given
metric). The barycenter can thus be seen as the combination of these
distributions and so can be used to build a MME. The barycenter is defined
with respect to a metric (i.e. a distance between distributions). Here, we
explore two barycenters based on different metrics: the $L_{2}$-distance and
the Wasserstein ($W_{2}$)-distance. We show that the $L_{2}$-barycenter is in
fact equivalent to the well-known pooling method and compare it to the new
$W_{2}$-barycenter based method.

This first application of this framework to S2S prediction is illustrated for
the combination of two single-model ensembles to predict the winter surface
temperature over Europe. We show that despite the superior skill of one of the
single-model ensemble over the other, it is still advantageous to combine them
into a barycenter. This reconfirms the interest of multi-model methods shown
by previous studies. However, the comparison of the two barycenter-based MMEs
does not single out clearly a better method. The best method depend on the
chosen metric, with the $L_{2}$-barycenter generally performing better with
respect to the mean Continuous Ranked Probability Score (CRPSm), and the
$W_{2}$-barycenter with respect to the Proportions of skillful forecasts
(CRPSp).

Moreover, we also highlight the importance of weighing the models within the
MME. The model’s weights have a significant impact on both MME’s performance.
They are particularly important in our case study where the single-model
ensembles have contrasting skill. In order to optimize the performance of the
MMEs, we learn the weights from past forecasts (using cross-validation here).
The weights are selected such as maximizing the MME’s skill with respect to
the chosen metrics: the CRPSm and CRPSp. This approach can easily be extended
to other metrics.

This study is a proof of concept to develop the framework and investigate the
properties of the barycenter-based MMEs. These results constitute a promising
first step towards improving S2S predictions using barycenters to merge
ensemble forecasts. A next step would be to implement the barycenter-based
MMEs for the combination of more than two models with weights estimated from
the data.

## Acknowledgment

The authors thank Naveen Goutham for code and advice on forecast calibration.
The authors thank the institut de Mathématiques pour la Planète Terre for
(partial) funding (iMPT 2021). This research was produced within the framework
of Energy4Climate Interdisciplinary Center (E4C) of IP Paris and Ecole des
Ponts ParisTech. This research was supported by 3rd Programme
d’Investissements d’Avenir [ANR-18-EUR-0006-02].

## References

  * Alessandri et al., (2011) Alessandri, A., Borrelli, A., Navarra, A., Arribas, A., Déqué, M., Rogel, P., and Weisheimer, A. (2011).  Evaluation of Probabilistic Quality and Value of the ENSEMBLES Multimodel Seasonal Forecasts: Comparison with DEMETER.  Monthly Weather Review, 139(2):581 – 607.
  * Becker et al., (2014) Becker, E., van den Dool, H., and Zhang, Q. (2014).  Predictability and Forecast Skill in NMME.  Journal of Climate, 27(15):5891 – 5906.
  * Bertossa et al., (2023) Bertossa, C., Hitchcock, P., DeGaetano, A., and Plougonven, R. (2023).  Coherent Bimodal Events in Ensemble Forecasts of 2-m Temperature.  Weather and Forecasting.
  * Büeler et al., (2020) Büeler, D., Beerli, R., Wernli, H., and Grams, C. M. (2020).  Stratospheric influence on ECMWF sub-seasonal forecast skill for energy-industry-relevant surface weather in European countries.  Quarterly Journal of the Royal Meteorological Society, 146(733):3675–3694.
  * Casanova and Ahrens, (2009) Casanova, S. and Ahrens, B. (2009).  On the Weighting of Multimodel Ensembles in Seasonal and Short-Range Weather Forecasting.  Monthly Weather Review, 137(11):3811 – 3822.
  * Ferrone et al., (2017) Ferrone, A., Mastrangelo, D., and Malguzzi, P. (2017).  Multimodel probabilistic prediction of 2 m-temperature anomalies on the monthly timescale.  Advances in Science and Research, 14:123–129.
  * Gelaro et al., (2017) Gelaro, R., McCarty, W., Suárez, M. J., Todling, R., Molod, A., Takacs, L., Randles, C. A., Darmenov, A., Bosilovich, M. G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C., Akella, S., Buchard, V., Conaty, A., da Silva, A. M., Gu, W., Kim, G.-K., Koster, R., Lucchesi, R., Merkova, D., Nielsen, J. E., Partyka, G., Pawson, S., Putman, W., Rienecker, M., Schubert, S. D., Sienkiewicz, M., and Zhao, B. (2017).  The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2).  Journal of Climate, 30(14):5419 – 5454.
  * Gneiting et al., (2005) Gneiting, T., Raftery, A. E., Westveld, A. H., and Goldman, T. (2005).  Calibrated probabilistic forecasting using ensemble model output statistics and minimum crps estimation.  Monthly Weather Review, 133(5):1098 – 1118.
  * Gonzalez et al., (2021) Gonzalez, P. L. M., Brayshaw, D. J., and Ziel, F. (2021).  A new approach to extended-range multimodel forecasting: Sequential learning algorithms.  Quarterly Journal of the Royal Meteorological Society, 147(741):4269–4282.
  * Goutham et al., (2022) Goutham, N., Plougonven, R., Omrani, H., Parey, S., Tankov, P., Tantet, A., Hitchcock, P., and Drobinski, P. (2022).  How Skillful Are the European Subseasonal Predictions of Wind Speed and Surface Temperature?  Monthly Weather Review, 150(7):1621 – 1637.
  * Hagedorn et al., (2012) Hagedorn, R., Buizza, R., Hamill, T. M., Leutbecher, M., and Palmer, T. N. (2012).  Comparing TIGGE multimodel forecasts with reforecast-calibrated ECMWF ensemble forecasts.  Quarterly Journal of the Royal Meteorological Society, 138(668):1814–1827.
  * Hagedorn et al., (2005) Hagedorn, R., Doblas-Reyes, F. J., and Palmer, T. (2005).  The rationale behind the success of multi-model ensembles in seasonal forecasting — I. Basic concept.  Tellus A: Dynamic Meteorology and Oceanography, 57(3):219–233.
  * Hamill, (2012) Hamill, T. M. (2012).  Verification of TIGGE Multimodel and ECMWF Reforecast-Calibrated Probabilistic Precipitation Forecasts over the Contiguous United States.  Monthly Weather Review, 140(7):2232 – 2252.
  * Haughton et al., (2015) Haughton, N., Abramowitz, G., Pitman, A., and Phipps, S. J. (2015).  Weighting climate model ensembles for mean and variance estimates.  Climate Dynamics, 45(11):3169–3181.
  * Heizenreder et al., (2006) Heizenreder, D., Trepte, S., and Denhard, M. (2006).  SRNWP-PEPS: A regional multi-model ensemble in Europe.  The European Forecaster: Newsletter of the WGCEF, 11.
  * Hersbach, (2000) Hersbach, H. (2000).  Decomposition of the Continuous Ranked Probability Score for Ensemble Prediction Systems.  Weather and Forecasting, 15(5):559 – 570.
  * Hersbach et al., (2020) Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.-N. (2020).  The ERA5 global reanalysis.  Quarterly Journal of the Royal Meteorological Society, 146(730):1999–2049.
  * Kalnay, (2003) Kalnay, E. (2003).  Atmospheric modeling, data assimilation and predictability.  Cambridge University Press.
  * Karpechko et al., (2018) Karpechko, A. Y., Charlton-Perez, A., Balmaseda, M., Tyrrell, N., and Vitart, F. (2018).  Predicting Sudden Stratospheric Warming 2018 and Its Climate Impacts With a Multimodel Ensemble.  Geophysical Research Letters, 45(24):13,538–13,546.
  * Kharin and Zwiers, (2002) Kharin, V. V. and Zwiers, F. W. (2002).  Climate Predictions with Multimodel Ensembles.  Journal of Climate, 15(7):793 – 799.
  * Kirtman et al., (2014) Kirtman, B. P., Min, D., Infanti, J. M., Kinter, J. L., Paolino, D. A., Zhang, Q., van den Dool, H., Saha, S., Mendez, M. P., Becker, E., Peng, P., Tripp, P., Huang, J., DeWitt, D. G., Tippett, M. K., Barnston, A. G., Li, S., Rosati, A., Schubert, S. D., Rienecker, M., Suarez, M., Li, Z. E., Marshak, J., Lim, Y.-K., Tribbia, J., Pegion, K., Merryfield, W. J., Denis, B., and Wood, E. F. (2014).  The North American Multimodel Ensemble: Phase-1 Seasonal-to-Interannual Prediction; Phase-2 toward Developing Intraseasonal Prediction.  Bulletin of the American Meteorological Society, 95(4):585 – 601\.
  * Leung et al., (1999) Leung, L. R., Hamlet, A. F., Lettenmaier, D. P., and Kumar, A. (1999).  Simulations of the ENSO Hydroclimate Signals in the Pacific Northwest Columbia River Basin.  Bulletin of the American Meteorological Society, 80(11):2313 – 2330\.
  * Manzanas et al., (2019) Manzanas, R., Gutiérrez, J. M., Bhend, J., Hemri, S., Doblas-Reyes, F. J., Torralba, V., Penabad, E., and Brookshaw, A. (2019).  Bias adjustment and ensemble recalibration methods for seasonal forecasting: a comprehensive intercomparison using the C3S dataset.  Climate Dynamics, 53(3):1287–1305.
  * Materia et al., (2020) Materia, S., Ángel G. Muñoz, Álvarez Castro, M. C., Mason, S. J., Vitart, F., and Gualdi, S. (2020).  Multimodel Subseasonal Forecasts of Spring Cold Spells: Potential Value for the Hazelnut Agribusiness.  Weather and Forecasting, 35(1):237 – 254.
  * Matheson and Winkler, (1976) Matheson, J. E. and Winkler, R. L. (1976).  Scoring Rules for Continuous Probability Distributions.  Management Science, 22(10):1087–1096.
  * McNemar, (1947) McNemar, Q. (1947).  Note on the sampling error of the difference between correlated proportions or percentages.  Psychometrika, 12(2):153–157.
  * Monhart et al., (2018) Monhart, S., Spirig, C., Bhend, J., Bogner, K., Schär, C., and Liniger, M. A. (2018).  Skill of Subseasonal Forecasts in Europe: Effect of Bias Correction and Downscaling Using Surface Observations.  Journal of Geophysical Research: Atmospheres, 123(15):7999–8016.
  * Murphy, (1973) Murphy, A. H. (1973).  A New Vector Partition of the Probability Score.  Journal of Applied Meteorology and Climatology, 12(4):595 – 600\.
  * Ning et al., (2014) Ning, L., Carli, F. P., Ebtehaj, A. M., Foufoula-Georgiou, E., and Georgiou, T. T. (2014).  Coping with model error in variational data assimilation using optimal mass transport.  Water Resources Research, 50(7):5817–5830.
  * Palmer et al., (2004) Palmer, T. N., Alessandri, A., Andersen, U., Cantelaube, P., Davey, M., Délécluse, P., Déqué, M., Díez, E., Doblas-Reyes, F. J., Feddersen, H., Graham, R., Gualdi, S., Guérémy, J.-F., Hagedorn, R., Hoshen, M., Keenlyside, N., Latif, M., Lazar, A., Maisonnave, E., Marletto, V., Morse, A. P., Orfila, B., Rogel, P., Terres, J.-M., and Thomson, M. C. (2004).  DEVELOPMENT OF A EUROPEAN MULTIMODEL ENSEMBLE SYSTEM FOR SEASONAL-TO-INTERANNUAL PREDICTION (DEMETER).  Bulletin of the American Meteorological Society, 85(6):853 – 872\.
  * Papayiannis et al., (2018) Papayiannis, G. I., Galanis, G. N., and Yannacopoulos, A. N. (2018).  Model aggregation using optimal transport and applications in wind speed forecasting.  Environmetrics, 29(8):e2531.
  * Pegion et al., (2019) Pegion, K., Kirtman, B. P., Becker, E., Collins, D. C., LaJoie, E., Burgman, R., Bell, R., DelSole, T., Min, D., Zhu, Y., Li, W., Sinsky, E., Guan, H., Gottschalck, J., Metzger, E. J., Barton, N. P., Achuthavarier, D., Marshak, J., Koster, R. D., Lin, H., Gagnon, N., Bell, M., Tippett, M. K., Robertson, A. W., Sun, S., Benjamin, S. G., Green, B. W., Bleck, R., and Kim, H. (2019).  The Subseasonal Experiment (SubX): A Multimodel Subseasonal Prediction Experiment.  Bulletin of the American Meteorological Society, 100(10):2043 – 2060.
  * Peyré and Cuturi, (2020) Peyré, G. and Cuturi, M. (2020).  Computational Optimal Transport.
  * Raftery et al., (2005) Raftery, A. E., Gneiting, T., Balabdaoui, F., and Polakowski, M. (2005).  Using bayesian model averaging to calibrate forecast ensembles.  Monthly Weather Review, 133(5):1155 – 1174.
  * Rajagopalan et al., (2002) Rajagopalan, B., Lall, U., and Zebiak, S. E. (2002).  Categorical climate forecasts through regularization and optimal combination of multiple gcm ensembles.  Monthly Weather Review, 130(7):1792 – 1811.
  * Robertson et al., (2004) Robertson, A. W., Lall, U., Zebiak, S. E., and Goddard, L. (2004).  Improved Combination of Multiple Atmospheric GCM Ensembles for Seasonal Prediction.  Monthly Weather Review, 132(12):2732 – 2744.
  * Robin et al., (2019) Robin, Y., Vrac, M., Naveau, P., and Yiou, P. (2019).  Multivariate stochastic bias corrections with optimal transport.  Hydrology and Earth System Sciences, 23(2):773–786.
  * Robin et al., (2017) Robin, Y., Yiou, P., and Naveau, P. (2017).  Detecting changes in forced climate attractors with Wasserstein distance.  Nonlinear Processes in Geophysics, 24(3):393–405.
  * Santambrogio, (2015) Santambrogio, F. (2015).  Progress in nonlinear differential equations and their applications.  In Optimal Transport for Applied Mathematicians: Calculus of Variations, PDEs, and Modeling, volume 87. Birkhäuser Cham.
  * Schulzweida, (2022) Schulzweida, U. (2022).  CDO User Guide.
  * Smith et al., (2013) Smith, D. M., Scaife, A. A., Boer, G. J., Caian, M., Doblas-Reyes, F. J., Guemas, V., Hawkins, E., Hazeleger, W., Hermanson, L., Ho, C. K., Ishii, M., Kharin, V., Kimoto, M., Kirtman, B., Lean, J., Matei, D., Merryfield, W. J., Müller, W. A., Pohlmann, H., Rosati, A., Wouters, B., and Wyser, K. (2013).  Real-time multi-model decadal climate predictions.  Climate Dynamics, 41(11):2875–2888.
  * Specq et al., (2020) Specq, D., Batté, L., Déqué, M., and Ardilouze, C. (2020).  Multimodel Forecasting of Precipitation at Subseasonal Timescales Over the Southwest Tropical Pacific.  Earth and Space Science, 7(9):e2019EA001003.
  * Takaya, (2019) Takaya, Y. (2019).  Chapter 12 - Forecast System Design, Configuration, and Complexity.  In Robertson, A. W. and Vitart, F., editors, Sub-Seasonal to Seasonal Prediction, pages 245–259. Elsevier.
  * Vigaud et al., (2017) Vigaud, N., Robertson, A. W., and Tippett, M. K. (2017).  Multimodel Ensembling of Subseasonal Precipitation Forecasts over North America.  Monthly Weather Review, 145(10):3913 – 3928.
  * Vigaud et al., (2020) Vigaud, N., Tippett, M. K., Yuan, J., Robertson, A. W., and Acharya, N. (2020).  Spatial Correction of Multimodel Ensemble Subseasonal Precipitation Forecasts over North America Using Local Laplacian Eigenfunctions.  Monthly Weather Review, 148(2):523 – 539.
  * Villani, (2003) Villani, C. (2003).  Topics in Optimal Transportation.  American Mathematical Society.
  * Vissio et al., (2020) Vissio, G., Lembo, V., Lucarini, V., and Ghil, M. (2020).  Evaluating the Performance of Climate Models Based on Wasserstein Distance.  Geophysical Research Letters, 47(21):e2020GL089385.
  * Vissio and Lucarini, (2018) Vissio, G. and Lucarini, V. (2018).  Evaluating a stochastic parametrization for a fast–slow system using the Wasserstein distance.  Nonlinear Processes in Geophysics, 25(2):413–427.
  * Vitart et al., (2017) Vitart, F., Ardilouze, C., Bonet, A., Brookshaw, A., Chen, M., Codorean, C., Déqué, M., Ferranti, L., Fucile, E., Fuentes, M., Hendon, H., Hodgson, J., Kang, H.-S., Kumar, A., Lin, H., Liu, G., Liu, X., Malguzzi, P., Mallas, I., Manoussakis, M., Mastrangelo, D., MacLachlan, C., McLean, P., Minami, A., Mladek, R., Nakazawa, T., Najm, S., Nie, Y., Rixen, M., Robertson, A. W., Ruti, P., Sun, C., Takaya, Y., Tolstykh, M., Venuti, F., Waliser, D., Woolnough, S., Wu, T., Won, D.-J., Xiao, H., Zaripov, R., and Zhang, L. (2017).  The Subseasonal to Seasonal (S2S) Prediction Project Database.  Bulletin of the American Meteorological Society, 98(1):163 – 173\.
  * Wanders and Wood, (2016) Wanders, N. and Wood, E. F. (2016).  Improved sub-seasonal meteorological forecast skill using weighted multi-model ensemble simulations.  Environmental Research Letters, 11(9):094007.
  * Wang et al., (2020) Wang, Y., Ren, H.-L., Zhou, F., Fu, J.-X., Chen, Q.-L., Wu, J., Jie, W.-H., and Zhang, P.-Q. (2020).  Multi-Model Ensemble Sub-Seasonal Forecasting of Precipitation over the Maritime Continent in Boreal Summer.  Atmosphere, 11(5).
  * Weigel et al., (2008) Weigel, A. P., Liniger, M. A., and Appenzeller, C. (2008).  Can multi-model combination really enhance the prediction skill of probabilistic ensemble forecasts?  Quarterly Journal of the Royal Meteorological Society, 134(630):241–260.
  * Wilcoxon, (1945) Wilcoxon, F. (1945).  Individual Comparisons by Ranking Methods.  Biometrics Bulletin, 1(6):80–83.
  * Wilks, (2019) Wilks, D. S. (2019).  Chapter 9 - Forecast Verification.  In Wilks, D. S., editor, Statistical Methods in the Atmospheric Sciences (Fourth Edition), pages 369–483. Elsevier, fourth edition edition.
  * Zheng et al., (2019) Zheng, C., Chang, E. K.-M., Kim, H., Zhang, M., and Wang, W. (2019).  Subseasonal to Seasonal Prediction of Wintertime Northern Hemisphere Extratropical Cyclone Activity by S2S and NMME Models.  Journal of Geophysical Research: Atmospheres, 124(22):12057–12077.

##  Appendix A Wasserstein distance as a measure of performance

Let consider an ensemble forecast with $N$ members and and the corresponding
deterministic observation, both written as discrete probability distributions
on the space $\Omega=\mathbb{R}^{n_{t}}$, respectively $\mu$ and $\nu$.

$\displaystyle\mu=\sum\limits_{i=1}^{N}a_{i}\delta_{\mathbf{x}_{i}}\quad$
$\displaystyle\nu=\delta_{\mathbf{y}}$ (11)

with $X=(\mathbf{x}_{1},\ldots,\mathbf{x}_{N})\in\Omega^{N}$ being the
ensemble members, $\mathbf{a}=(a_{1},\ldots,a_{N})\in\Sigma_{N}$ their
weights, and $\mathbf{y}$ the observed time-series.

In that case, there is only one feasible transport matrix in
$U(\mathbf{a},\mathbf{b})$: $\mathbf{T}=\mathbf{a}$. That is, the mass of the
different members all go to the observation. The equation for the squared
2-Wasserstein distance becomes

$\displaystyle W_{2}^{2}(\mu,\nu)$ $\displaystyle=\min_{T\in
U(\mathbf{a},\mathbf{b})}\sum\limits_{i=1}^{N}t_{i,j}.\|\mathbf{x}_{i}-\mathbf{y}\|^{2}=\sum\limits_{i=1}^{N}a_{i}.\sum\limits_{k=1}^{nt}\left(x_{i,k}-y_{k}\right)^{2}$
$\displaystyle W_{2}^{2}(\mu,\nu)$
$\displaystyle=\sqrt{\sum\limits_{i=1}^{N}\sum\limits_{k=1}^{nt}a_{i}.\left(x_{i,k}-y_{k}\right)^{2}}$

Thus, the $W_{2}$-distance is equal to the RMSE over the time-steps and the
ensemble members (up to a multiplicative factor $1/\sqrt{n_{t}}$). It does not
take into account the information on the forecast uncertainty carried by the
ensemble spread.

##  Appendix B Energy distance and its associated barycenter

Let $\mu_{1}$ and $\mu_{2}$ be two distributions, and $F_{1}$ and $F_{2}$ be
their CDF. That is, $\forall x\in\mathbb{R}$

$\displaystyle F_{1}(x)=\int_{-\infty}^{x}\mu_{1}(t)dt\quad$
$\displaystyle\text{and}\quad F_{2}(x)=\int_{-\infty}^{x}\mu_{2}(t)dt$

The squared energy distance between $\mu_{1}$ and $\mu_{2}$ is

$\displaystyle\mathcal{E}^{2}(\mu_{1},\mu_{2})$
$\displaystyle=\int_{-\infty}^{+\infty}\left[F_{1}(x)-F_{2}(x)\right]^{2}dx$
$\displaystyle=\int_{-\infty}^{+\infty}\left[\int_{-\infty}^{x}\mu_{1}(t)dt-\int_{-\infty}^{x}\mu_{2}(t)dt\right]^{2}dx$

The energy barycenter of $\mu_{1}$ and $\mu_{2}$ is the solution of the
following minimisation problem

$\displaystyle\mu_{\mathcal{E}}^{\alpha}=\operatorname*{arg\,min}_{\mu}\underbrace{\alpha.\mathcal{E}^{2}(\mu,\mu_{1})+(1-\alpha).\mathcal{E}^{2}(\mu,\mu_{2})}_{\mathcal{B}(\mu)}$

We have:

$\displaystyle\frac{d}{d\mu}\mathcal{B}(\mu)$
$\displaystyle=\frac{d}{d\mu}\left[\alpha\int_{-\infty}^{+\infty}\left[F(x)-F_{1}(x)\right]^{2}dx+(1-\alpha)\int_{-\infty}^{+\infty}\left[F(x)-F_{2}(x)\right]^{2}dx\right]$
$\displaystyle=\frac{d}{d\mu}\int_{-\infty}^{+\infty}\left(\alpha\left[\int_{-\infty}^{x}\mu(t)dt-\int_{-\infty}^{x}\mu_{1}(t)dt\right]^{2}+(1-\alpha)\left[\int_{-\infty}^{x}\mu(t)dt-\int_{-\infty}^{x}\mu_{2}(t)dt\right]^{2}\right)dx$
$\displaystyle=\int_{-\infty}^{+\infty}\frac{d}{d\mu}\left(\alpha\left[\int_{-\infty}^{x}\left(\mu(t)-\mu_{1}(t)\right)dt\right]^{2}+(1-\alpha)\left[\int_{-\infty}^{x}\left(\mu(t)-\mu_{2}(t)\right)dt\right]^{2}\right)dx$
$\displaystyle=\int_{-\infty}^{+\infty}\left(2\alpha\int_{-\infty}^{x}\left[\mu(t)-\mu_{1}(t)\right]dt\mu(x)\right)dx+\int_{-\infty}^{+\infty}\left(2(1-\alpha)\int_{-\infty}^{x}\left[\mu(t)-\mu_{2}(t)\right]dt\mu(x)\right)dx$
$\displaystyle=2\int_{-\infty}^{+\infty}\left(\int_{-\infty}^{x}\left(\mu(t)-\left[\alpha\mu_{1}(t)+(1-\alpha)\mu_{2}(t)\right]\right)dt\right)\mu(x)dx$

Thus,

$\displaystyle\frac{d}{d\mu}\mathcal{B}(\mu_{L_{2}}^{\alpha})$
$\displaystyle=0$

where $\mu_{L_{2}}^{\alpha}=\alpha\mu_{1}+(1-\alpha)\mu_{2}$ is the $L_{2}$
barycenter. The $L_{2}$ barycenter is also a barycenter for the energy
distance.

##  Appendix C CRPS and $L_{2}$-barycenter

Let $\mu_{1}$ and $\mu_{2}$ be two distributions corresponding to two ensemble
forecasts, and $\mu_{L_{2}}^{\alpha}$ be their $L_{2}$-barycenter. The CRPS of
the barycenter $\mu_{L_{2}}^{\alpha}$ with respect to the observation
$\nu_{obs}$ can be computed as follow

$\displaystyle\text{CRPS}(\mu_{L2},\nu_{obs})$
$\displaystyle=\alpha\text{CRPS}(\mu_{1},\nu_{obs})+(1-\alpha)\text{CRPS}(\mu_{2},\nu_{obs})-\alpha(1-\alpha)\text{CRPS}(\mu_{1},\mu_{2})$

##### Proof:

The CRPS can be applied to any pair of distributions $\mu_{1}$ and $\mu_{2}$.
Let consider two discrete probability distributions on the space
$\Omega=\mathbb{R}^{n_{t}}$ such that

$\displaystyle\mu_{1}=\sum\limits_{i=1}^{N_{1}}a_{i}\delta_{\mathbf{x}_{i}}\quad\text{and}\quad\mu_{2}=\sum\limits_{j=1}^{N_{2}}b_{j}\delta_{\mathbf{y}_{j}}$

with $X=(\mathbf{x}_{1},\ldots,\mathbf{x}_{N_{1}})\in\Omega^{N_{1}}$,
$Y=(\mathbf{y}_{1},\ldots,\mathbf{y}_{N_{2}})\in\Omega^{N_{2}}$, and where
$\mathbf{a}=(a_{1},\ldots,a_{N_{1}})\in\Sigma_{N_{1}}$ and
$\mathbf{b}=(b_{1},\ldots,b_{N_{2}})\in\Sigma_{N_{2}}$ are probability
vectors.

If we substitute the expressions for $\mu_{1}$ and $\mu_{2}$ in the formula of
the CRPS from Eq. 7, we obtain

$\displaystyle\text{CRPS}(\mu_{1},\mu_{2})$
$\displaystyle=\int_{-\infty}^{+\infty}\left[F_{\mu_{1}}(u)-F_{\mu_{2}}(u)\right]^{2}du$
$\displaystyle=\int_{-\infty}^{+\infty}\left[\sum_{i=1}^{N_{1}}a_{i}I(x_{i}\leq
u)-\sum_{j=1}^{N_{2}}b_{j}I(y_{j}\leq u)\right]^{2}du$
$\displaystyle=\sum_{i=1}^{N_{1}}\sum_{j=1}^{N_{2}}a_{i}b_{j}\lvert
x_{i}-y_{j}\rvert-\sum_{i=1}^{N-1}\sum_{k=i+1}^{N_{1}}a_{i}a_{k}\lvert
x_{i}-x_{k}\rvert-\sum_{j=1}^{N_{2}-1}\sum_{l=j+1}^{N_{2}}b_{j}b_{l}\lvert
y_{j}-y_{l}\rvert$ (12)

Then, we substitute $\nu_{obs}=\sum_{k=1}^{N_{o}}d_{k}\delta_{\mathbf{y_{k}}}$
and
$\mu_{L2}=\sum_{k=1}^{N_{1}+N_{2}}c_{k}\delta_{\mathbf{x_{L2}}_{k}}=\alpha\sum_{i=1}^{N_{1}}a_{i}\delta_{\mathbf{x_{1}}_{i}}+(1-\alpha)\sum_{j=1}^{N_{2}}b_{j}\delta_{\mathbf{x_{2}}_{j}}$
in the formula of the CRPS from Eq. 12.

$\displaystyle\text{CRPS}(\mu_{L2},\nu_{obs})$
$\displaystyle=\sum_{i=1}^{N_{1}+N_{2}}\sum_{j=1}^{N_{o}}c_{i}d_{j}\lvert{x_{L2}}_{i}-y_{j}\rvert-\sum_{i=1}^{N_{1}+N_{2}-1}\sum_{j=i+1}^{N_{1}+N_{2}}c_{i}c_{j}\lvert{x_{L2}}_{i}-{x_{L2}}_{j}\rvert-\sum_{i=1}^{N_{o}-1}\sum_{j=i+1}^{N_{o}}d_{i}d_{j}\lvert
y_{i}-y_{j}\rvert$ $\displaystyle=\sum_{i=1}^{N_{1}}\sum_{j=1}^{N_{o}}\alpha
a_{i}d_{j}\lvert{x_{1}}_{i}-y_{j}\rvert+\sum_{i=1}^{N_{2}}\sum_{j=1}^{N_{o}}(1-\alpha)b_{i}d_{j}\lvert{x_{2}}_{i}-y_{j}\rvert$
$\displaystyle~{}~{}~{}~{}-\sum_{i=1}^{N_{1}-1}\sum_{j=i+1}^{N_{1}}\alpha^{2}a_{i}a_{j}\lvert{x_{1}}_{i}-{x_{1}}_{j}\rvert-\sum_{i=1}^{N_{1}}\sum_{j=1}^{N_{2}}\alpha(1-\alpha)a_{i}b_{j}\lvert{x_{1}}_{i}-{x_{2}}_{j}\rvert$
$\displaystyle~{}~{}~{}~{}-\sum_{i=1}^{N_{2}-1}\sum_{j=i+1}^{N_{2}}(1-\alpha)^{2}b_{i}b_{j}\lvert{x_{2}}_{i}-{x_{2}}_{j}\rvert-\sum_{i=1}^{N_{o}-1}\sum_{j=i+1}^{N_{o}}d_{i}d_{j}\lvert
y_{i}-y_{j}\rvert$
$\displaystyle=\alpha\left[\sum_{i=1}^{N_{1}}\sum_{j=1}^{N_{o}}a_{i}d_{j}\lvert{x_{1}}_{i}-y_{j}\rvert-\sum_{i=1}^{N_{1}-1}\sum_{j=i+1}^{N_{1}}a_{i}a_{j}\lvert{x_{1}}_{i}-{x_{1}}_{j}\rvert-\sum_{i=1}^{N_{o}-1}\sum_{j=i+1}^{N_{o}}d_{i}d_{j}\lvert
y_{i}-y_{j}\rvert\right]$
$\displaystyle~{}~{}~{}~{}+(1-\alpha)\left[\sum_{i=1}^{N_{2}}\sum_{j=1}^{N_{o}}b_{i}d_{j}\lvert{x_{2}}_{i}-y_{j}\rvert-\sum_{i=1}^{N_{2}-1}\sum_{j=i+1}^{N_{2}}b_{i}b_{j}\lvert{x_{2}}_{i}-{x_{2}}_{j}\rvert-\sum_{i=1}^{N_{o}-1}\sum_{j=i+1}^{N_{o}}d_{i}d_{j}\lvert
y_{i}-y_{j}\rvert\right]$
$\displaystyle~{}~{}~{}~{}+\alpha(1-\alpha)\sum_{i=1}^{N_{1}-1}\sum_{j=i+1}^{N_{1}}a_{i}a_{j}\lvert{x_{1}}_{i}-{x_{1}}_{j}\rvert+(1-\alpha)\alpha\sum_{i=1}^{N_{2}-1}\sum_{j=i+1}^{N_{2}}b_{i}b_{j}\lvert{x_{2}}_{i}-{x_{2}}_{j}\rvert$
$\displaystyle~{}~{}~{}~{}-\alpha(1-\alpha)\sum_{i=1}^{N_{1}}\sum_{j=1}^{N_{2}}a_{i}b_{j}\lvert{x_{1}}_{i}-{x_{2}}_{j}\rvert$
$\displaystyle=\alpha\text{CRPS}(\mu_{1},\nu_{obs})+(1-\alpha)\text{CRPS}(\mu_{2},\nu_{obs})$
$\displaystyle~{}~{}~{}~{}-\alpha(1-\alpha)\left[\sum_{i=1}^{N_{1}}\sum_{j=1}^{N_{2}}a_{i}b_{j}\lvert{x_{1}}_{i}-{x_{2}}_{j}\rvert-\sum_{i=1}^{N_{1}-1}\sum_{j=i+1}^{N_{1}}a_{i}a_{j}\lvert{x_{1}}_{i}-{x_{1}}_{j}\rvert-\sum_{i=1}^{N_{2}-1}\sum_{j=i+1}^{N_{2}}b_{i}b_{j}\lvert{x_{2}}_{i}-{x_{2}}_{j}\rvert\right]$
$\displaystyle=\alpha.\text{CRPS}(\mu_{1},\nu_{obs})+(1-\alpha).\text{CRPS}(\mu_{2},\nu_{obs})-\alpha.(1-\alpha).\text{CRPS}(\mu_{1},\mu_{2})$

##  Appendix D Spatial performance’s significance

Table 2: Two-sided Wilcoxon test’s p-value for the spatial CRPSm. The green color indicate the weeks at which the two models are significantly different at the $5\%$ significance level (i.e. p-values$<0.05$). | Week 3 | Week 4 | Week 5 | Week 6
---|---|---|---|---
ECMWF - NCEP |  green!50 3.2e-08 |  green!50 3.6e-06 |  green!50 6.3e-04 |  green!50 8.1e-09
W2 - NCEP |  green!50 2.0e-10 |  green!50 2.4e-08 |  green!50 1.3e-05 |  green!50 1.4e-10
L2 - NCEP |  green!50 6.4e-12 |  green!50 9.9e-10 |  green!50 8.4e-07 |  green!50 3.9e-12
W2 - ECMWF |  red!50 4.5e-01 |  red!50 8.5e-02 |  green!50 1.0e-02 |  red!50 8.6e-01
L2 - ECMWF |  red!50 4.9e-01 |  red!50 9.6e-02 |  green!50 4.6e-03 |  red!50 7.5e-01
W2 - L2 |  red!50 1.8e-01 |  red!50 7.7e-02 |  green!50 5.0e-05 |  green!50 3.1e-02
Table 3: Two-sided Wilcoxon test’s p-value for the spatial CRPSp. The green color indicate the weeks at which the two models are significantly different at the $5\%$ significance level (i.e. p-values$<0.05$). | Week 3 | Week 4 | Week 5 | Week 6
---|---|---|---|---
ECMWF - NCEP |  green!50 1.4e-08 |  green!50 2.0e-04 |  green!50 3.7e-04 |  green!50 8.5e-08
W2 - NCEP |  green!50 1.1e-16 |  green!50 2.1e-12 |  green!50 1.3e-12 |  green!50 1.1e-14
L2 - NCEP |  green!50 2.0e-13 |  green!50 9.3e-08 |  green!50 1.6e-08 |  green!50 2.2e-11
W2 - ECMWF |  green!50 1.5e-02 |  green!50 2.8e-04 |  green!50 2.7e-03 |  green!50 4.9e-02
L2 - ECMWF |  red!50 2.1e-01 |  green!50 1.5e-02 |  red!50 1.0e-01 |  red!50 6.5e-01
W2 - L2 |  green!50 3.6e-04 |  green!50 1.1e-06 |  green!50 2.5e-06 |  green!50 5.7e-06