# The added value of combining solar irradiance data and forecasts: A probabilistic benchmarking exercise

Philippe Lauret, Rodrigo Alonso-Suárez, Rodrigo Amaro E Silva, John

Boland, Mathieu David, Wiebke Herzberg, Josselin Le Gall La Salle, Elke

Lorenz, Lennard Visser, Wilfried van Sark, et al.

► To cite this version:

Philippe Lauret, Rodrigo Alonso-Suárez, Rodrigo Amaro E Silva, John Boland, Mathieu David, et al.. The added value of combining solar irradiance data and forecasts: A probabilistic benchmarking exercise. Renewable Energy, 2024, 237, 10.1016/j.renene.2024.121574. hal-04752368

HAL Id: hal-04752368

https://univ-reunion.hal.science/hal-04752368v1

Submitted on 24 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# The added value of combining solar irradiance data and forecasts: A probabilistic benchmarking exercise

Philippe Lauret<sup>a,∗</sup>, Rodrigo Alonso-Suárez<sup>b</sup>, Rodrigo Amaro e Silva<sup>c</sup>, John Boland<sup>d</sup>, Mathieu David<sup>a</sup>, Wiebke Herzberg<sup>e</sup>, Josselin Le Gall La Salle<sup>a</sup>, Elke Lorenz<sup>e</sup>, Lennard Visser<sup>f</sup>, Wilfried van Sark<sup>f</sup>, Tobias Zech<sup>e</sup>

## ARTICLE INFO

**Keywords:**

Probabilistic solar forecasting

Benchmarking exercise

Blended point forecast

CRPS

IEA PVPS T16

## ABSTRACT

Despite the growing awareness in academia and industry of the importance of solar probabilistic forecasting for further enhancing the integration of variable photovoltaic power generation into electrical power grids, there is still no benchmark study comparing a wide range of solar probabilistic methods across various local climates. Having identified this research gap, experts involved in the activities of IEA PVPS T16¹ agreed to establish a benchmarking exercise to evaluate the quality of intra-hour and intra-day probabilistic irradiance forecasts.

The tested forecasting methodologies are based on different input data including ground measurements, satellite-based forecasts and Numerical Weather Predictions (NWP), and different statistical methods are employed to generate probabilistic forecasts from these. The exercise highlights different forecast quality depending on the method used, and more importantly, on the input data fed into the models.

In particular, the benchmarking procedure reveals that the association of a point forecast that blends ground, satellite and NWP data with a statistical technique generates high-quality probabilistic forecasts. Therefore, in a subsequent step, an additional investigation was conducted to assess the added value of such a blended point forecast on forecast quality. Three new statistical methods were implemented using the blended point forecast as input.

To ensure a fair evaluation of the different methods, we calculate a skill score that measures the performance of the proposed model relative to that of a trivial baseline model. The closer the skill score is to 100%, the more efficient the method is. Overall, skill scores of methods that use the blended point forecast ranges from 42% to 46% for the intra-hour scenario and 27% to 32% for the intra-day scenario. Conversely, methods that do not use the blended point forecast exhibit skill scores ranging from 33% to 43% for intra-hour forecasts and 8% to 16% for intra-day forecasts.

These results suggest that using (a) blended point forecasts that optimally combine different sources of input data and (b) a post-processing with a statistical method to produce the quantile forecasts is an effective and consistent way to generate high-quality intra-hour or intra-day probabilistic forecasts.

<footnote><sup>a</sup> University of La Réunion - PIMENT laboratory, 15, avenue René Cassin, 97715 Saint-Denis, France</footnote>

<footnote><sup>b</sup> Laboratorio de Energía Solar (LES), Departamento de Física, CENUR Litoral Norte, Udelar, Uruguay</footnote>

<footnote><sup>c</sup> O.I.E. Centre Observation, Impacts, Energy, MINES Paris, PSL Research University, 06904 Sophia Antipolis, France</footnote>

<footnote><sup>d</sup> Industrial AI, Centre for Industrial and Applied Mathematics, UniSA STEM, University of South Australia, Mawson Lakes Boulevard, Mawson Lakes, SA 5095, Australia</footnote>

<footnote><sup>e</sup> Fraunhofer Institute for Solar Energy Systems ISE, Heidenhofstr. 2, 79110 Freiburg, Germany</footnote>

<footnote><sup>f</sup> Copernicus Institute of Sustainable Development, Utrecht University, Princetonlaan 8a, 3584 CB, Utrecht, The Netherlands</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: philippe.lauret@univ-reunion.fr (P. Lauret), rodrigoa@fing.edu.uy (R. Alonso-Suárez), rodrigo.amaro_e_silva@minesparis.psl.eu (R. Amaro e Silva), John.Boland@unisa.edu.au (J. Boland), mathieu.david@univ-reunion.fr (M. David), wiebke.herzberg@ise.fraunhofer.de (W. Herzberg), josselin.le-gal-la-salle@univ-reunion.fr (J. Le Gall La Salle), elke.lorenz@ise.fraunhofer.de (E. Lorenz), l.r.visser@uu.nl (L. Visser), w.g.j.h.m.vansark@uu.nl (W. van Sark), tobias.zech@ise.fraunhofer.de (T. Zech).</footnote>

<footnote n="1">¹ International Energy Agency - Photovoltaic Power Systems - Solar Resource for High Penetration and Large Scale Applications.</footnote>
# 1. Introduction

Accurate forecasts of solar energy generation play a crucial role in effectively integrating solar power into existing grids and reducing associated expenses [1]. This is because power output from photovoltaic plants (PV) is greatly influenced by weather conditions, making it highly variable. Consequently, having precise information about future solar power production is essential to minimize the need for costly balancing services and power reserves. Hence, enhancing solar forecasting models to increase the value of solar power generation becomes critically important. This work will focus in Global Horizontal Irradiance (GHI) forecasting since it is deemed one of the main drivers for solar power forecasting [2].

Today, users are faced with a myriad of forecasting methods. They can vary in the nature of the approach, the kind of inputs fed into the models, the outputs they produce or even in the forecasting horizon under consideration. Review publications aim to identify and structure all these elements, often comparing the performance values reported across the literature [3-6]. However, benchmark studies aim to do so using a methodology that ensures a fair comparison, i.e. the same locations and training and test periods. This has been done for many different aspects of solar forecasting: the post-processing of numerical weather prediction models (NWP) [7], baseline approaches [8], autoregressive statistical learning approaches [9,10], cloud motion vector techniques [11], deep learning approaches using sky images as inputs [12], among others. In other words, these benchmark studies ensure a comprehensive understanding and comparable indicators of the benefits associated with a particular approach. For instance, in the frame of IEA SHC Task 46 [13], Lorenz et al. [14] designed a standardized procedure to evaluate the accuracy of day-ahead deterministic irradiance forecasts.

However, all these previous studies focus on deterministic forecasting and dismiss the inherent uncertainty of a forecast. Indeed, when it comes to decision-making for grid operators, utilities, aggregators, balancing responsible parties and others, having not only a point forecast but also an associated uncertainty or prediction interval becomes immensely valuable. In other words, reliable probabilistic predictions can significantly enhance the integration of variable energy sources within the energy network, leading to improved efficiency [15]. Unlike the mature field of wind power forecasting, where probabilistic forecasting is well-established [15-18], probabilistic solar forecasting is still relatively nascent [19-21]. Consequently, there are considerably less benchmark studies focusing on solar probabilistic forecasts.

A literature review restricted to intra-hour/intra-day solar probabilistic forecasts reveals that a few studies have started to address this gap. Among others, one can cite the following works. Grantham et al. [22] used a non parametric bootstrapping method for generating prediction intervals of GHI at a forecast horizon of 1 h. The bootstrap technique requires a point forecast which is, in their work, delivered by a linear auto-regressive (AR) model. With only past ground data, David et al. [23] used a combination of point forecast ARMA model and a parametric GARCH model to generate intra-hour (up to 1 h ahead with a time step of 10 min) and intra-day (up to 6 h ahead with a time step of 1 h) GHI probabilistic forecasts. Lauret et al. [24] evaluated the quality of three probabilistic models for intra-day solar forecasting. A linear quantile regression technique is used to build three models for generating 1 to 6 h ahead probabilistic forecasts. Inputs of the models are either only ground data or ground data with day-ahead forecasts provided by the European Center for Medium-Range Weather Forecast (ECMWF). The results demonstrated that the Numerical Weather Prediction (NWP) exogenous inputs improve the quality of the intra-day probabilistic forecasts. Using only past ground GHI measurements, David et al. [25] set up a combination of 3 points forecasting methods and 7 probabilistic methods to issue intra-day GHI forecasts. None of the model combinations clearly outperformed the others. However, regardless of the point forecasting method used, linear models in quantile

regression, weighted quantile regression and gradient boosting decision trees appear to produce probabilistic forecasts with higher quality than the other proposed methods. In their work, Alonso-Suárez et al. [26] developed three models aimed at generating intra-day probabilistic GHI forecasts, spanning lead times from 10 min to 3 h with a granularity of 10 min. The initial model solely relies on historical ground measurements. The second model enhances the first one by integrating a variability metric derived from these historical ground measurements. The third model incorporates satellite albedo as an additional input. A linear quantile regression is employed to create directly (i.e. without using a point forecast) a range of quantiles summarizing the predictive distributions of the global solar irradiance. The findings demonstrate that the inclusion of satellite data further enhances the quality of the probabilistic forecasts. Mazorra-Aguiar et al. [27] assessed the performance of two approaches for solar probabilistic forecasting to generate intra-day solar forecasts covering time horizons from 1 h to 6 h. The first approach involves a two-step process. Initially, point forecasts are generated for each forecast horizon, followed by the utilization of quantile regression techniques to estimate the prediction intervals. The second methodology directly predicts the quantiles of the predictive distribution using past ground data as input. Yang et al. [28] benchmarked 5 forecasting intrahour/intraday solar probabilistic methods (including notably an Analog Ensemble method and a linear quantile regression technique) on a standardized dataset set up by Pedro et al. [29]. All the proposed methods generate directly the quantile forecasts without resorting to a point deterministic forecast. The findings clearly highlight the significance of exogenous inputs in probabilistic solar forecasting, as all methods demonstrate enhanced results upon the integration of exogenous features computed from sky, satellite images and NWP outputs provided by the NAM the North American Mesoscale (NAM) forecast system.

Finally, it must be noted that specific methods based on Cloud-Motion Vector (CMV) approach or combination of sky and satellite images have been recently proposed in the literature. For instance, Carrière et al. [30] designed a CMV-based probabilistic method which is an extension of the deterministic CMV approach by adding Gaussian noise to the norm and direction of the cloud motion vector estimates. Paletta et al. [31] used a hybrid deep learning method combining sky images, satellite observations and/or past ground irradiance to generate intra-hour solar forecasts.

Following the previous literature review, the following comments can be made. To evaluate the quantile forecasts two methodologies can be distinguished. The first one leverages on a point deterministic forecast to produce with a specific statistical technique the prediction intervals. The second one generates directly (i.e. without resorting to a point forecast) the quantile forecasts. Regarding the first methodology, no work tries to evaluate the impact of a high quality point forecast on the generation of probabilistic forecasts.

Moreover, to the best of our knowledge, no benchmarking exercise has been conducted to compare classical probabilistic techniques like quantile regression or Analog Ensemble with a CMV-based probabilistic approach on multiple sites experiencing different climate conditions.

Therefore, as part of IEA PVPS T16 [32], experts engaged in Activity 3.3 on solar probabilistic forecasts found it essential to complement these previous studies regarding intra-hour and intra-day solar probabilistic forecasts. In other words, it appears important to experts of the IEA PVPS T16 to propose to the solar forecasters community a comprehensive benchmarking exercise related to intra-day and intra-hour solar irradiance forecasts.

To this end, five participants set-up a benchmarking exercise based on a shared ground measurement, satellite and NWP data. Eight European sites with diverse climatic conditions were chosen for this purpose. The proposed benchmarking procedure is implemented to compare 15-min irradiance probabilistic forecasts up to 6 h issued by each participant with their own forecasting methods. In particular, to fill the gaps highlighted by the literature review, it appeared first
<table><caption>**Table 1**

List of participants. The code associated to each participant also identifies the forecasting methodology used by the participant.</caption><thead><tr><th>Participant</th><th>Code/Method</th><th>Forecasting methodology</th><th>Input data</th></tr></thead><tbody><tr><td>Mines Paris (OIE)</td><td>OIE</td><td>CMV-based probabilistic approach</td><td>Satellite data</td></tr><tr><td>University of La Réunion (PIMENT)</td><td>PIMENT</td><td>Parametric method ARMA-GARCH</td><td>Ground data</td></tr><tr><td>Fraunhofer (ISE)</td><td>ISE</td><td>Blended point forecast +Analog Ensemble</td><td>Ground + satellite + NWP data</td></tr><tr><td>Utrecht University (UU)</td><td>UU</td><td>Non-linear Quantile Regression Forest</td><td>Ground data</td></tr><tr><td>Laboratorio de Energía Solar (Udelar)</td><td>LES</td><td>Linear quantile regression</td><td>Ground + satellite data</td></tr></tbody></table>

important to the IEA PVPS experts to jointly evaluate a CMV probabilistic system with traditional quantile forecasting methods. Second, an assessment of the impact of a high quality point forecast on the quality of the generated probabilistic predictions is also conducted in this work.

To evaluate the quality of the probabilistic predictions, different diagnostic tools and scoring rules can be employed [33]. For user convenience, the verification scheme should be kept simple. For that, we propose using the reliability diagram as a visual diagnostic tool and the Continuous Ranked Probability Score (CRPS), as the numerical score. It is commonly adopted by the community in the verification process of solar irradiance probabilistic forecasts.

Further, in this work, unlike most of the previous studies, and in order to better highlight the skill of a forecasting method, we propose the numerical decomposition of the CRPS into the reliability and resolution components as in Lauret et al. [33].

The rest of the paper is organized as follows. Firstly, this paper introduces the benchmarking exercise. Section 3 details the data used in the exercise while Section 4 gives an overview of the diverse forecasting methodologies. The verification framework is presented in Section 5. Section 6 gives the main results of the benchmark and Section 7 discusses the impact of combining a blended point forecast with statistical techniques to generate probabilistic solar forecasts. Finally, Section 8 concludes this paper.

## 2. The benchmarking exercise

In the frame of the IEA PVPS Task 16 [32], five participants agreed to set up a benchmarking exercise related to intra-day and intra-hour solar irradiance probabilistic forecasting. Table 1 lists the participants of this benchmark, the code that will be used to identify them throughout the following sections and plots, as well as the forecasting methods used. For this exercise, each participant submitted their forecasts under the form of quantile forecasts (i.e. the quantiles of the predictive distribution).

Together, the participants designed the framework that would guide this benchmark, namely the type of input data that could be fed into the forecasting models, the forecast horizons to be considered, and the selection of the probability levels of the quantile forecasts. Table 2 gives details regarding these decisions. Thus, each participant was responsible for generating 15-min irradiance probabilistic forecasts up to 2 h ahead (for intra-hour) and up to 6 h (for intra-day) for 8 selected European sites (described in Section 3). The verification of the forecasts was conducted blindly by one of the authors of this paper. The Mean Absolute Error (MAE) of the median of the predictive distribution was also calculated. This allows for a direct comparison between the performance of the pointwise deterministic median and that of the proposed probabilistic methods, as the MAE is equivalent to the CRPS for deterministic counterparts.

## 3. Data for the benchmarking exercise

### 3.1. Ground measurements

It is crucial to use identical data for evaluation when comparing various prediction methods. The selected dataset comprises 15-min measured GHI values from eight locations of Europe. We restricted the

<table><caption>**Table 2**

Parameters of the benchmarking exercise. Details regarding the verification metrics are provided in Section 5.</caption><tr><td>Type of Input Data</td><td>- Ground GHI measurements<br/>- Satellite estimates<br/>- NWP forecasts<br/>- Solar geometry variables (e.g. Solar Zenith Angle (SZA), etc.)</td></tr><tr><td>Forecast horizon</td><td>- Intra-hour: 8 horizons (15 to 120 min, in 15-min steps)<br/>- Intra-day: 16 horizons (135 to 360 min, in 15-min steps)</td></tr><tr><td>Forecasts specification</td><td>15 GHI quantile forecasts with probability levels of [0, 0.025, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.975, 1]</td></tr><tr><td>Verification metrics</td><td>- Reliability diagram<br/>- CRPS<br/>- CRPS skill score with the CSD-CLIM model as baseline — see [34].<br/>- MAE of the median of the predictive distribution</td></tr></table>

evaluation to European sites since the methods of two of the participants used satellite data only covering most of the European domain. The evaluation period spans from January 2017 to December 2018. The year 2017 was chosen for the training set of the different methods described below, while 2018 was used for testing these methods.

The original reference database comprises high temporal resolution GHI data (1 min) that were collected for a benchmarking exercise of modeled solar radiation data [35]. This benchmark exclusively incorporates quality-assured data, meticulously checked using an extensive range of best practices and newly established quality-control procedures [35]. These procedures encompass automated and manual data quality tests along with descriptive quality flagging, conducted by a team of experts of the IEA PVPS T16 subtask 1 [32].

The 15-min dataset results from a downsampling of these original 1-min. More precisely, the 1-min raw GHI were averaged at 15 min resolution. Also, in case of missing raw data, a linear interpolation is done if the gap is below 1 h otherwise the whole day is discarded. Finally, data for which solar elevation ≤ 10° have been filtered out and are consequently not taken into account in the evaluation process.

Table 3 gives all the details related to each site. Note that, except for the TAB site provided by CIEMAT/DLR [36] and the MIL site provided by Ricerca sul Sistema Energetico (RSE) all sites are part of the Baseline Surface Radiation Network (BSRN) [37].

### 3.2. Satellite data

In this study, two data sets of GHI estimates based on satellite data are considered. Both are based on images obtained by the SEVIRI instrument onboard the Meteosat Second Generation (MSG) satellite using the MSG 15-minute visible channel with a spatial resolution of approximately 1 × 2 km at Europe.

The first data set used for the OIE and LES forecasts is derived from the satellite images using the Heliosat-4 model [38]. The GHI estimates in the second data set used for the ISE forecasts are based on a modified version of the Heliosat method [39].

Generally, for forecasting purposes, a sequence of satellite images is used to infer cloud motion vectors (CMV), i.e. vectors that describe cloud advection, which can be extrapolated into the future to make
<caption>**Table 3**</caption>

Locations and key figures of ground measurements used for the benchmark. Column “Köppen C.” lists the Köppen-Geiger climate classification of each site while column “Mean GHI” gives the average GHI of the test dataset.

<table><thead><tr><th>Site</th><th>Country</th><th>Code</th><th>Latitude (°N)</th><th>Longitude (°E)</th><th>Altitude (m)</th><th>Köppen C.</th><th>Source</th><th>Mean GHI (W/m²)</th></tr></thead><tbody><tr><td>Cabauw</td><td>NLD</td><td>CAB</td><td>51.9711</td><td>4.9267</td><td>0</td><td>Cfb</td><td>BSRN</td><td>315.0</td></tr><tr><td>Carpentras</td><td>FRA</td><td>CAR</td><td>44.083</td><td>5.059</td><td>100</td><td>Csa</td><td>BSRN</td><td>411.1</td></tr><tr><td>Cener</td><td>ESP</td><td>CEN</td><td>42.816</td><td>-1.601</td><td>471</td><td>Cfb</td><td>BSRN</td><td>381.6</td></tr><tr><td>Milan</td><td>ITA</td><td>MIL</td><td>45.4761</td><td>9.2545</td><td>150</td><td>Cfa</td><td>RSE</td><td>394.1</td></tr><tr><td>Palaiseau</td><td>FRA</td><td>PAL</td><td>48.713</td><td>2.208</td><td>156</td><td>Cfb</td><td>BSRN</td><td>347.7</td></tr><tr><td>Payerne</td><td>CHE</td><td>PAY</td><td>46.815</td><td>6.944</td><td>491</td><td>Cfb</td><td>BSRN</td><td>368.9</td></tr><tr><td>Plataforma Solar</td><td>ESP</td><td>TAB</td><td>37.0909</td><td>-2.3581</td><td>500</td><td>Bsk</td><td>CIEMAT/DLR</td><td>499.5</td></tr><tr><td>Toravere</td><td>EST</td><td>TOR</td><td>58.254</td><td>26.462</td><td>70</td><td>Dfb</td><td>BSRN</td><td>318.1</td></tr></tbody></table>

a prediction. For a better understanding of the assumptions and limitations of the use of CMV, the reader is directed to Lorenz et al. [2].

### 3.3. NWP forecasts

The blended forecasts results from the combination of three types of inputs namely a persistence forecast, a satellite-based forecast and the deterministic ECMWF IFS forecast (see Section 4.3.1). Thereto, ISE includes the atmospheric model high resolution 10-day forecast (HRES) product of the ECMWF IFS forecast with a spatial resolution of 0.125° and a time resolution of 1 h. The forecast of shortwave solar radiation downwards (ssrd) is used at base-times 0:00 UTC and 12:00 UTC. It is spatially smoothed over 9 × 9 grid points, and upsampled from the original 1 h time resolution to 15 min via clear sky index (see Eq. (1)) interpolation.

## 4. Description of the probabilistic methods

Let us recall that all probabilistic methods depicted in this work generate quantile forecasts with the probability levels given in Table 2. However, each participant utilizes their own method to produce this set of quantiles (listed and briefly described above in Table 1).

Three classes of approaches are employed in this benchmark. The first class extends to the probabilistic domain a framework traditionally used to produce deterministic CMV-based forecasts, mainly by adding Gaussian noise to its inputs (similarly to a Monte Carlo approach). That is the case of Mines Paris OIE. The second class corresponds to a two-step approach where a deterministic forecast is generated and then used as input by a statistical technique to generate the quantile forecasts. This is the case of participants PIMENT and ISE. Conversely, the third class produces directly (in one step) the quantile forecasts from a set of input variables. This is the case of participants UU and LES.

Finally, in the field of solar forecasting, it is a standard procedure to detrend the Global Horizontal Irradiance (GHI) time series due to its non-stationary nature, characterized by daily cycles and annual seasonalities [40]. This detrending process involves utilizing the output of a clear sky model. Specifically, a new deseasonalized series, known as the clear sky index $k_c$ time series, is derived by employing the following data transformation

$$k_c = \frac{I}{I_c}, \qquad (1)$$

where *I* is the measured global horizontal irradiance and $I_c$ is the output of a clear sky model. All the proposed forecasting models described below make use of the clear sky index $k_c$ time series. However, it should be noted that the choice of the clear sky model may vary depending on the participant.

### 4.1. Description of Mines Paris OIE model

This approach, proposed originally in Carrière et al. [30], combines physical and statistical elements and leverages a standard satellite-based solar forecasting framework which is traditionally used for deterministic forecasting. Fig. 1 gives an overview of the principle of the method.

First, a 25 × 25 grid with 0.04° resolution centered in the location of interest is defined. For each grid cell, satellite-derived time series for GHI and its clear-sky expectation ($I_c$) are obtained from the Copernicus Atmospheric Monitoring Services (CAMS) Radiation product [41] using the pvlib Python interface [42]. Note that while the native resolution of this product depends on the distance to nadir, the CAMS Radiation product adjusts to any requested coordinates by means of interpolation. Based on this data, the corresponding clear-sky index $k_c$ grid was calculated according Eq. (1) and its spatial resolution was further increased by a factor of 3 through a 2D linear interpolation. Then, the CMV of each downscaled $k_c$ cell is inferred using an Optical Flow technique following the work of Chow et al. [43] and using an efficient method proposed by Liu [44].

Deterministic CMV-based forecasts are calculated, for example, by extrapolating the $k_c$ grid in space (according to the CMVs and the forecast horizon) and selecting the advected $k_c$ value which is closest to the location of interest. In this approach, an Eulerian spatial extrapolation is considered, where the clouds are assumed to move in a straight-line trajectory.

Here, the probabilistic aspect is enabled by three elements: (i) a Gaussian noise distribution relative to the CMV estimates, namely to the norm and direction of a vector; (ii) a Gaussian noise distribution relative to the $k_c$ estimates; and (iii) a monitoring perimeter. The first two aim to describe the uncertainty associated with the estimation of the satellite-derived variables, whereas the third defines a distance threshold below which an extrapolated grid cell is considered a plausible forecast candidate. This allows the generation of a set of plausible advection scenarios. The third corresponds to a distance threshold below which an extrapolated $k_c$ pixel is considered a plausible forecast candidate. Thus, the combination of the viable candidates from all the generated scenarios constitutes an ensemble of $k_c$ forecasts from which an empirical cumulative distribution function (CDF) is built.

Finally, the $k_c$ CDF is converted back to GHI by multiplying it with the $I_c$ obtained from CAMS Radiation (see Eq. (1)), which considers the McClear clear-sky model, which accounts for water vapor and aerosol effects [45]. To mitigate potential calibration issues, the forecasts are post-processed by first considering the baseline model CSD-CLIM [34] for defining the bounding quantiles (Q0 and Q100) and then implementing a quantile mapping approach for calibration, adjusting the effective probability rate of each quantile to its theoretical rate according to the training data.

More details on this implementation can be found in Carrière et al. [30], including the parameters assumed for the considered sources of Gaussian noise. Concerning the model implementation, a few remarks:

* It is only tested for the locations covered by the CAMS Radiation service (i.e., covered by the Meteosat Second Generation geostationary satellite);
<figure><img src="image_3.png" /><figcaption>Fig. 1. Mechanisms leveraging the probabilistic CMV approach: (a) the generation of advection scenarios by inputting noise to the base CMV; (b) the consideration of all advected grid cells that fall inside a monitoring perimeter.</figcaption></figure>

* It is only tested for the horizon range between 15 min and 2 h;

* In two situations, the baseline approach proposed by Le Gal La Salle et al. [34] is considered instead of the CMV-based one: (i) when for a given day, forecast time, and horizon, there is yet no available satellite image; (ii) when this approach leads to less than 50 $k_c$ candidates, which possibly compromises the representativeness of the produced ensemble.

## 4.2. Description of PIMENT model

The PIMENT model is based on a parametric approach commonly used in the financial domain. It combines a AutoRegressive Moving Average model (ARMA) and Generalized AutoRegressive Conditional Heteroskedasticity (GARCH), which successively generate a point forecast and then its associated uncertainty. This combination in the field of solar energy has been first introduced by David et al. [23]. The model is applied to the clear sky index $k_c$ time series with the McClear model [45] selected as the clear-sky model.

The AutoRegressive Moving Average model (ARMA) stands as a prevalent and widely-applied method in time series prediction. Its extensive utilization in forecasting renewable energy has underscored its competitive edge, owed largely to its parsimonious nature. Notably, its application spectrum encompasses the forecasting of solar irradiance among other domains [25,46]. A general formulation of an ARMA(p,q) model with p autoregressive (AR) terms and q moving average (MA) terms is given by Tsay [47]. Its application to the h-ahead forecast of a variable y is given by the following equation

$$ \hat{y}_{t+h} = \alpha_0 + \sum_{i=1}^{p} \alpha_i \times y_{t-i+1} + \sum_{j=1}^{q} \beta_j \times \epsilon_{t-j+1}, \quad (2) $$

with h = 1, 2, ... the forecast horizon and $\alpha_0, \alpha_1, ..., \alpha_p, \beta_1, ..., \beta_q$ the coefficients to be estimated. The error term $\epsilon$ is the difference between the previous forecasts and observations as defined in the following equation:

$$ \epsilon_t = \hat{y}_t - y_t, \quad (3) $$

The ARCH (AutoRegressive Conditional Heteroskedasticity) models, introduced by Engle [48], is used to model the variance of time series in the financial domain. These models are particularly efficient to predict changes in variance over the time, for instance the error of point forecast generated with an ARMA model [49]. PIMENT applied the Generalized AutoRegressive Conditional Heteroskedasticity (GARCH) model proposed by Bollerslev [49], which gives a more parsimonious formulation than the simple ARCH model. In GARCH models, the conditional variance is a linear function of lagged squared error terms and

also lagged conditional variance terms [50]. The general formulation of a GARCH(r,s) model, with r error terms, s conditional variance terms and an horizon of forecast h, is given by:

$$ \hat{y}_{t+h} = \hat{\sigma}_{t+h} \times \epsilon_t, \quad (4) $$

with $\epsilon$ an uniformly distributed random variable with a null mean and a unitary variance, and $\hat{\sigma}_{t+h}$ the predicted standard deviation given by

$$ \hat{\sigma}_{t+h}^2 = \gamma_0 + \sum_{i=1}^{r} \gamma_i \times \epsilon_{t-i+1} + \sum_{j=1}^{s} \delta_j \times \sigma_{t-j+1}^2, \quad (5) $$

As for the ARMA models, $\gamma_0, \gamma_1, ..., \gamma_r, \delta_1, ..., \delta_s$ are the coefficients to be estimated. There are numerous methods to estimate these coefficients. The two most widely used are the least squares (LS) and the Maximum Likelihood Estimation (MLE) methods. Here, we propose to implement the Recursive Least Squares (RLS) method, which is a variation of the LS method. This method reduces the computational cost with the coefficient of the model being updated in real-time as new data become available. The RLS method is very efficient in an operational context where forecast have to be timely delivered.

To determine the lag parameters p, q of the ARMA model, PIMENT ran the model on the training year for different combinations of the lag parameters with values varying from 1 to 10. The best combination is the one that minimizes the RMSE of the point forecast. For the probabilistic part, PIMENT used a GARCH(1,1), which is appropriate for the error times series of the point forecast.

## 4.3. Description of Fraunhofer ISE model

The approach by Fraunhofer ISE consists of two steps. In a first step, blended point forecasts are derived from different input data (Section 4.3.1). In a second step, quantile forecasts are generated from these blended point forecasts using the Analog Ensemble (AnEn) method (Section 4.3.2).

### 4.3.1. Blended forecasts by Fraunhofer ISE

Deriving blended forecasts from several distinct input forecasts using statistical or machine learning methods is a common approach in deterministic forecasting [see e.g.2].

Here, GHI forecasts are generated by blending three different types of forecasts

* a persistence forecast,

* a satellite-based forecast,

* and the deterministic ECMWF IFS forecast (see Section 3.3).
<figure><img src="image_4.png" /><figcaption>Fig. 2. MAE vs. forecast horizon, computed for all eight sites together for the test year of 2018. The slight deviations of the median of the predictive distribution denoted here in the figure by ensemble median compared to the median of the predictive distribution of the benchmark (see Fig. 4b) result from differences in the filtering process.</figcaption></figure>

The persistence forecast is created by deriving a clear sky index from the latest measurement, which is then extrapolated into the future. Fraunhofer ISE employs the clear-sky model from Dumortier [51] and the turbidity model from Bourges [52] to compute clear sky irradiance and the clear sky index, not only for persistence, but for all modeling steps described in this section. The satellite-based forecasts are based on CMVs derived from MSG satellite images (see section 3.2) following Kühnert et al. [53]. The CMVs are computed using a block matching algorithm. Future images are obtained by repetitive application of the cloud motion vectors. Finally, smoothing filters depending on forecast lead times are applied to the future images.

The three different forecast types are blended using a set of linear regression models, fitted for each forecast horizon and time of the day. This allows to adjust the regression weights to the varying performance of the different input forecasts in dependence of the forecast horizon (see Fig. 2). The resulting blended forecast can be written as

$$I_{\text{blend}}^{h,\tau} = c_{\text{pers}}^{h,\tau} \cdot I_{\text{pers}}^{h,\tau} + c_{\text{cmv}}^{h,\tau} \cdot I_{\text{cmv}}^{h,\tau} + c_{\text{nwp}}^{h,\tau} \cdot I_{\text{nwp}}^{h,\tau} \quad (6)$$

where $I_X$ are the GHI of the blended, persistence, CMV or NWP forecast, $c_X$ are the corresponding regression weights, and $h$ and $\tau$ denote the index of the forecast horizon and time of the day respectively. When determining the regression weights $c_X^{h,\tau}$ for a forecast horizon $h$, data from $h \pm 2$ horizons were included in the fit, to enable the generation of a smooth blended output forecast from the different linear models. The regression weights are trained for the year 2017 and applied to generate the forecasts for 2018. With multiple sites of observation and forecast data, either one set of the blending parameters for all sites or separate sets for each site can be derived. It was decided for the single site training, to obtain close-to-bias-free blended forecasts for each site.

When generating the blended forecasts, missing data is handled in the following way: If one of the component forecasts is missing, a simple mean of the remaining component forecasts is calculated instead of using the regression weights. For horizons or times of the day for which no regression weights could be determined, NWP data is returned. It should be noted here that persistence as well as satellite based forecasts can only be calculated after sunrise and before sunset,

which impacts availability of these forecasts in the early morning hours depending also on the forecast horizon. For example, if the earliest satellite based forecast could be calculated at 7:00 am, forecasts for horizons of 4 h ahead are available only from 11:00 am onwards.

The performance of the different forecasts, quantified in terms of MAE, is shown as a function of forecast horizon in Fig. 2. Here the different availabilities of the forecasts discussed above have to be considered. For reasons of comparability, calculation of MAE includes only data points for which all displayed forecasts are available. The Figure shows that the persistence forecasts are best performing for short horizons up to 45 min, satellite-based forecasts are best performing for intermediate horizons up to about 3 h, and the NWP performs best for even larger horizons. The blended forecasts always exhibit a lower MAE than any of the individual input forecasts. They form the basis to derive probabilistic forecasts in a next step.

### 4.3.2. Analog Ensemble

The Analog Ensemble (AnEn) is a non-parametric ensemble prediction method [54,55]. The method is based on the evaluation of historic observations and deterministic forecasts. Past forecasts are compared to the current forecast and the observations corresponding to the most similar forecast situations, called analogs, then form an ensemble of possible future values. The quantities, which are used to measure similarity, are called predictors. Alessandrini et al. [56] constructed an AnEn for PV power prediction and used the GHI, solar elevation and azimuth, cloud cover, and ambient temperature (T2M) from the deterministic ECMWF IFS forecast as predictors. Here, an AnEn to predict the GHI instead of the PV power is created. Furthermore, just one predictor quantity is used for the AnEn, namely the forecasted clear sky index.

To set up the Analog Ensemble, in a first step, the clear sky index values of the blended forecasts and the measurements are computed. Situations with measurement-based clear sky index values above 1.2 are excluded. To identify the analogs for a current forecast value with a forecast horizon $h_0$, similarity to past forecasts is evaluated using the Euclidean distance of the forecast values of five horizons $\{h_0 - 30 \text{ min}, h_0 - 15 \text{ min}, h_0, h_0 + 15 \text{ min}, h_0 + 30 \text{ min}\}$, centered
around $h_0$. This window helps to reduce fluctuations in the distance measure and to improve meteorological similarity between situations. The measurement-based clear sky indices corresponding to the 40 most similar situations are taken to form the AnEn.

Similar to a k-nearest-neighbor regression, the AnEn has no explicit training phase, instead the analogs are selected from a search space at the time a prediction is made. The analog search is performed separately for each forecast horizon and each time of the day, reflecting different uncertainties in dependence of the forecast horizon (see Fig. 2) and the time of the day. The search space for the analogs is composed using a rolling window of the last 180 days and integrating all eight European sites together, resulting in 1440 historic data points from which the 40 analogs are selected. Integrating the different sites increases the search space and thus the reliability of the ensemble. It is made possible by using the clear sky index as a predictor instead of GHI and close-to-bias-free forecasts for the different sites. The ensemble of clear sky indices is, then, transformed back to GHI values by multiplication with clear sky irradiance. The quantiles are finally obtained by a linear interpolation of the ensemble members, see method 7 of (Hyndman and Fan, 1996).

## 4.4. Description of Utrecht University (UU) model

Quantile regression forest (QRF) is a nonlinear ensemble model that is based on the random forest regression (RF) model [57,58]. Similar to a RF model, QRF is made up of a predefined set of decision trees that exist of a number of layers ($t_n$) and decision nodes ($2t_n$). The trees are constructed independently from each other by considering bootstrap samples from the training dataset in the training stage. The nodes are constructed by selecting a random subset of the predictor variables and optimizing the decision node on a preset loss function, e.g the mean squared error. In contrast to a RF model, the QRF model predicts a conditional distribution function (or weighted distribution of observations). Hence, given a set of predictor variables, each tree in the QRF model predicts the conditional quantiles of the target variable, i.e., GHI. Finally, a post-processing step is added in which each quantile value is set to be equal or higher than zero and equal or lower than the clear sky irradiance estimate ($I_c$).

The UU forecast model follows the approach described in Visser et al. [59] to find the optimal hyperparameter settings. Hence, the training set is split into several training and validation subsets, using k-fold cross-validation with $k = 8$ [60]. Once the optimal hyperparameter settings are found, the QRF model is trained considering the entire training set, i.e., one year (2017) and then applied to predict the GHI for the test year (2018).

The QRF model considered in this study operates endogenously. This implies that the model only relies on historical observations of the target variable, i.e., the GHI ($I$), as well as variables that are available at any time, i.e., the clear sky irradiance $I_c$. Using $I$ and $I_c$, we construct two additional variables: the clear-sky index $k_c$ (see Eq. (1)) and the expected GHI using a clear sky-based smart persistence model, similar as discussed in Section 4.3.

From these four main variables, a large set of predictor variables can be generated by simply considering a multitude of lagged values. In this study, UU optimizes the number of historical values by means of an iterative process. Hence, starting with a base model, at each iteration one lagged value is added, where after an evaluation if the addition leads to a significant performance improvement is made. The final set of variables considers 18 predictor variables, including: the previous eleven GHI measurements ($I(t), (...) , I(t-10)$), the clear-sky irradiance ($I_c(t+k)$), the clear-sky index ($k_c(t)$) and the persistence forecast considering the three most recent irradiance measurements ($I_{pers}(t+k, t), I_{pers}(t+k, t-1), I_{pers}(t+k, t-2)$).

## 4.5. Description of Udelar LES model

The LES forecast is an adaptation of the methodology proposed by Alonso-Suárez et al. [26]. This approach utilizes lagged ground measurements and geostationary satellite data as inputs for a Linear Quantile Regression (LQR) model, as described by Koenker and Bassett [61]. The LQR model is used to predict quantiles of the clear-sky index ($k_c$). These quantiles are then converted to quantiles of the Global Horizontal Irradiance (GHI) using the McClear clear-sky model (see Eq. (1)). While the mathematical formulation is relatively simple, the crucial aspect lies in the predictors' selection. The forecasting model incorporates the present time and the six preceding $k_c$ values, along with four other predictors derived from either the past $k_c$ values or a satellite space cell that surrounds the specific location. The first additional predictor is the local short-term variability ($\sigma_c$), which is calculated as the standard deviation of the last six changes in $k_c$. For more in-depth information on how $\sigma_c$ is calculated, please refer to Alonso-Suárez et al. [26]. The remaining three predictors are derived from the Heliosat-4 satellite estimates in a 25 × 25 px space cell provided by the Copernicus Atmosphere Monitoring Service (CAMS). By employing the McClear model, the clear-sky index is calculated for each pixel in the satellite cell. From this index, the average, standard deviation, and cloud coverage are computed and utilized as input variables in the LQR method. The cloud coverage is estimated as the fraction of pixels in the cell with $k_c < 0.85$. In summary, the inputs for the LQR model are six lagged $k_c$ values, the local short-term variability, and four variables related to the current time. These current time variables consist of the measured $k_c$, the space average and standard deviation of the satellite-derived $k_c$, and the satellite-estimated cloud coverage within the cell. Different LQR parameters are trained for each site, forecast horizon, and quantile.

# 5. Proposed evaluation framework

Visual diagnostic tools and quantitative scores are used to assess the quality of the probabilistic forecasts (i.e. the correspondence between ground truth and the forecasts). Diagnostic tools are used to visually assess the quality of probabilistic forecasts, while numerical scores are used to quantify the skills of a forecasting system and to rank competing prediction methods.

In this study, following the recommendation of Lauret et al. [33], we adopt the CRPS as the scoring rule to assess the overall performance of the forecasting method. Moreover, to gain a deeper understanding of the forecast skill of each forecasting method, we further decompose the CRPS into two components: reliability and resolution. In case of predictive distributions summarized by discrete quantile forecasts, Lauret et al. [33] proposed specific formulae to compute the CRPS and its related decomposition. The interested reader is referred to Lauret et al. [33] for more details regarding the computation of this CRPS decomposition. Another useful assessment is whether a prediction system outperforms a trivial baseline model. To this end, we compute the CRPS skill score with the climatological model CSD-CLIM [34] as the reference model.

Finally, in this work, we use reliability diagrams to visually evaluate reliability of the different forecasts.

### 5.1. Visual assessment with reliability diagrams

The reliability diagram serves as a graphical tool for assessing the reliability of probabilistic forecasts. In this paper, we follow the methodology established by Pinson et al. [62], which is tailored for predictive distributions summarized by quantile forecasts. Specifically, quantile forecasts are considered reliable when their stated probabilities match the observed proportions. In essence, over a sufficiently large evaluation dataset, the disparity between observed and nominal probabilities should be minimized [62].
One of the advantages of this representation is that it allows deviations from perfect reliability, represented by the diagonal line, to be readily visualized [62]. However, it is important to acknowledge that due to the finite sample of observation/forecast pairs and potential serial correlation in the sequence of forecasts and observations, observed proportions may not align precisely along the diagonal, even if the forecasts are perfectly reliable [62]. In other words, reliability diagrams can sometimes be misleading because even for perfectly reliable forecasts, deviations from the ideal diagonal case can be observed.

To address the limitations arising from the finite number of observation/forecast pairs, Bröcker and Smith [63] introduced reliability diagrams with consistency bars. Additionally, Pinson et al. [62] has proposed consistency bars that consider the combined effects of serial correlation and limited data. In this work, consistency bars are calculated according to Pinson et al. [62]. When interpreting reliability diagrams with consistency bars, it becomes clear that one cannot reject the hypothesis of the quantile forecasts being perfectly reliable if the observed proportions fall within the consistency bars. In practice, incorporating consistency bars into reliability diagrams can provide additional support for the user's, possibly subjective, assessment of the reliability of the different models.

## 5.2. Continuous Rank Probability Score (CRPS) and its decomposition

The Continuous Ranked Probability Score (CRPS) is a numerical score that quantifies the difference between the predicted and observed cumulative distribution functions (CDF) [64]. It is formulated as follows:

$$CRPS = \frac{1}{N} \sum_{i=1}^{N} \int_{-\infty}^{+\infty} \left[ \hat{F}_{f_{cst}}^{i}(x) - F_{x_{obs}}^{i}(x) \right]^2 dx, \quad (7)$$

where $\hat{F}_{f_{cst}}(x)$ is the predictive CDF of the variable of interest $x$ (e.g. GHI) and $F_{x_{obs}}(x)$ is a cumulative-probability step function that jumps from 0 to 1 at the point where the forecast variable $x$ equals the observation $x_0$ (i.e. $F_{x_{obs}}(x) = 1_{\{x \ge x_{obs}\}}$). The squared difference between the two CDFs is averaged over the $N$ observation/forecast pairs.

The CRPS score is designed to reward forecasts that concentrate their probabilities around the step function located at the observed value, promoting accuracy and precision in forecast predictions [65]. Put simply, the CRPS serves as a penalty for both insufficient resolution in predictive distributions and biased forecasts. It is worth noting that the CRPS is oriented negatively, meaning that smaller values indicate better performance, and it has the same unit as the forecast variable.

As previously mentioned and in accordance with its nature as a proper scoring rule [66], the CRPS can be decomposed into two fundamental aspects of probabilistic forecasts: reliability and resolution. This decomposition of the CRPS yields the following equation:

$$CRPS = REL+UNC-RES. \quad (8)$$

The reliability REL component of the CRPS provides an assessment of forecast biases, while the resolution RES component quantifies the improvement achieved by issuing case-dependent probability forecasts. The uncertainty component UNC, on the other hand, is inherent to the observations and cannot be influenced by the forecast system; it depends solely on the variability of the observed data [65].

Given that the CRPS is negatively oriented, the objective of a forecast system is to minimize the reliability component as much as possible, while also maximizing the resolution component. By employing this decomposition of the CRPS, a detailed evaluation of the forecast performance of different forecasting methods can be obtained.

Besides, in the case of deterministic forecasts, the CRPS reduces to the Mean Absolute Error (MAE). This characteristic enables a direct comparison between the performance of a probabilistic model and a deterministic one, or equivalently, it allows for assessing the additional value provided by a probabilistic approach [67]. In this study, we

calculate the Mean Absolute Error (MAE) of the forecast distribution's median to evaluate the extent to which the probabilistic approach enhances (or fails to enhance) the overall quality of the forecasts over its deterministic counterpart.

## 5.3. The CSD-CLIM model and the associated CRPS skill score

Probabilistic scores do not allow fair comparisons between different sites or datasets. To do so, it is a good practice to consider the relative performance against reference models [68] through skill scores. Over the past years, several benchmark models for probabilistic forecasting were introduced in the literature [68,69]. For this work, the CSD-CLIM model has been selected. For each site, the measurements of the training dataset are gathered according to a set of Nb bins of clear-sky irradiance. For each bin, CSD-CLIM generates a forecast cumulative distribution function (CDF) using only the historical GHI training data within that bin. As a result, Nb predictive CDFs are produced. In the test period, the clear-sky model is used to select the appropriate bin and the associated forecasting CDF. Thus, the CSD-CLIM approach is climatological in the sense that it only uses historical data and a clear-sky model. More details about theory and implementation are available in Le Gal La Salle et al. [34]. In this work, the McClear clear-sky model [45] has been chosen with 30 clear-sky irradiance bins.

A skill score represents the degree of improvement of a forecasting model compared to the reference baseline model. The CRPS skill score (CRPSS) reads as

$$CRPSS = 1 - \frac{CRPS_{model}}{CRPS_{reference}}. \quad (9)$$

# 6. Results

## 6.1. Reliability diagrams

Reliability diagrams related to each scenario (intra-hour or intra-day) for each participant and for all sites are given in Fig. 3. Note that the computation of reliability diagrams for all sites is performed by aggregating GHI observations and forecasts for each site into two distinct time series. This procedure will be also used to calculate the overall CRPS results in Section 6.2.1. Consistency bars for a 90% confidence level are individually computed for each nominal proportion. Note that comments related to Figs. 3f, 3g and 3h will be made in Section 7.

Irrespective of the scenario (intra-hour or intra-day), the visual analysis shows that only probabilistic forecasts derived from ISE (Fig. 3c) possibly has a high reliability. All the other forecasts are possibly non reliable namely those generated by OIE (Fig. 3a), PIMENT (Fig. 3b), UU (Fig. 3d) and to a lesser extent LES (Fig. 3e). More specifically, forecasts provided by OIE appear to be clearly non reliable and forecasts issued by LES, UU and PIMENT experience high deviations from the ideal line for high nominal proportions.

For the ISE forecasts for which the observed proportions lie within the consistency bars, while this does not confirm the perfect reliability of the quantile forecasts, it also does not allow us to confidently assert their lack of reliability at a 10% significance level.

## 6.2. CRPS and its decomposition

### 6.2.1. Overall results for all sites

Following Yang et al. [28], Tables 4 and 5 report the overall results (i.e. computed from the aggregation of forecasts and data of all sites) obtained by each metric in terms of “mean ± standard deviation”. More precisely, for intra-hour forecasting, the mean and standard deviation are computed from the 8 values of the metric corresponding to the 8 horizons while for intra-day forecasting, the mean and standard deviation are computed from the 16 values related to the 16 horizons (see Table 2). We recall here that intra-hour forecasts correspond to 15–120 min-ahead forecasting at 15 min timesteps while the intra-day
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Reliability diagrams for all sites averaged over all the forecast horizons related to each participant. Consistency bars for a 90% confidence level around the ideal line are individually computed for each nominal proportion. The red curves stand for intra-hour forecasts while the blue one are for intra-day forecasts.</figcaption></figure>

forecasts are for 120–360 min-ahead forecasting at 15 min timesteps. Also, in this section, we will comment on the results obtained by the five participants with their original proposed method (see Table 1). The methods BLEND-GARCH, BLEND-LQR and BLEND-QRF will be presented and discussed later in Section 7.

Regarding intra-hour forecasts (see Table 4), the best performer is ISE regardless of the metric while the worst one is PIMENT. In terms of skill score, the forecast skill of ISE is (in average) 46.6% while PIMENT exhibits a skill score of 32.9%. It appears clearly that the better performance of ISE originates from its better resolution and reliability. In line with the reliability diagrams of Fig. 3a, the CMV-based probabilistic method of OIE leads to poor results notably in terms of reliability. More generally, the quantitative reliability (REL) component of the CRPS confirms the visual diagnosis provided by the reliability diagrams in Fig. 3. Finally, it should be noted that the linear LQR method proposed by LES, fed with ground and satellite data, achieves similar results to the nonlinear QRF method of UU, which uses only ground data.

As shown by Table 5, for intra-day forecasts, the same comments made above for intra-hour forecasts still hold. However, except for ISE, the only participant integrating NWP forecasts, one can observe a strong decrease in forecast skill particularly for PIMENT for which a decrease of 25 points in the mean skill score is noted. Similar to deterministic forecasts (see Fig. 2), with increasing forecast horizon, the positive impact of integrating NWP forecasts in the modeling process is clearly demonstrated.

We complement the above quantitative metrics analysis summarized for intra-day and day-ahead by plotting over all the forecast horizons the numerical scores selected for this benchmarking exercise namely CRPS, CRPS reliability, CRPS resolution, MAE and CRPS skill score (see Fig. 4). It should also be noted that we deliberately use the same Y-scale for the CRPS and MAE plots to emphasize the improvement in quality brought by the probabilistic approach. Indeed, as shown by Figs. 4a and 4b, the CRPS (i.e., the MAE) of the median of the predictive distribution is clearly worse than the CRPS of the entire predictive distribution.

As shown by all the plots, the highest overall skill of the ISE forecasts is clearly demonstrated irrespective of the forecast horizon. In terms of CRPS skill score, over the whole range of forecast horizons, the best performer is ISE with skill scores between 60% and almost 30% while PIMENT forecasts lead to the worst forecasting results with a CRPS skill score ranging between 55% and almost 0%

Again, the specific method developed by OIE does not outperform the other methods. The forecasts issued by this technique are clearly non reliable and confirms the visual inspection of the related reliability diagram.

### 6.2.2. Results for each site

To gain a deeper insight in the performances of the different methods, Fig. 5 focuses on the CRPS values obtained by each contestant for each of the eight sites. As shown, again, the best CRPS values are obtained by ISE whatever the location. Interestingly, the parametric
<caption>**Table 4**</caption>

Intra-hour forecasts overall results. For each method, the metrics are presented as “mean ±standard deviation,” calculated over all forecast horizons. For these overall results, the CRPS of the CSD-CLIM is 82.3 W m⁻² and the uncertainty component UNC of the CRPS is 150.8 W m⁻².

<table><caption>**Table 5**</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>53.8 ± 5.6</td><td>34.6 ± 6.7</td><td>4.5 ± 0.3</td><td>101.8 ± 5.4</td><td>67.9 ± 8.3</td></tr><tr><td>PIMENT</td><td>55.2 ± 10.3</td><td>32.9 ± 12.5</td><td>3.9 ± 0.7</td><td>99.7 ± 9.6</td><td>72.9 ± 15.4</td></tr><tr><td>ISE</td><td>43.9 ± 6.6</td><td>46.6 ± 8.1</td><td>0.8 ± 0.2</td><td>107.9 ± 6.5</td><td>59.4 ± 9.4</td></tr><tr><td>UU</td><td>48.3 ± 8.4</td><td>41.3 ± 10.2</td><td>1.1 ± 0.2</td><td>103.8 ± 8.3</td><td>64.1 ± 12.2</td></tr><tr><td>LES</td><td>47.1 ± 8.6</td><td>42.8 ± 10.4</td><td>1.5 ± 0.3</td><td>105.5 ± 8.3</td><td>62.7 ± 12.1</td></tr><tr><td>BLEND-GARCH</td><td>47.7 ± 6.9</td><td>42.1 ± 8.3</td><td>2.0 ± 0.3</td><td>105.4 ± 6.6</td><td>61.6 ± 9.2</td></tr><tr><td>BLEND-LQR</td><td>47.3 ± 6.4</td><td>42.6 ± 7.8</td><td>2.1 ± 0.2</td><td>105.9 ± 6.2</td><td>61.4 ± 9.2</td></tr><tr><td>BLEND-QRF</td><td>44.2 ± 5.9</td><td>46.3 ± 7.2</td><td>1.1 ± 0.1</td><td>107.9 ± 5.8</td><td>59.5 ± 8.6</td></tr></tbody></table>

Intra-day forecasts overall results. For each method, the metrics are presented as “mean ±standard deviation,” calculated over all forecast horizons. OIE has “NA” values since OIE method is limited to intra-hour forecasting.

<table><caption>**Table 6**</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>PIMENT</td><td>76.0 ± 4.0</td><td>7.7 ± 4.9</td><td>4.7 ± 0.1</td><td>79.8 ± 3.9</td><td>105.3 ± 6.5</td></tr><tr><td>ISE</td><td>56.0 ± 2.0</td><td>31.9 ± 2.5</td><td>1.3 ± 0.2</td><td>96.3 ± 1.9</td><td>76.6 ± 2.7</td></tr><tr><td>UU</td><td>69.3 ± 5.6</td><td>15.8 ± 6.8</td><td>1.6 ± 0.2</td><td>83.3 ± 5.4</td><td>96.3 ± 8.9</td></tr><tr><td>LES</td><td>69.4 ± 5.5</td><td>15.6 ± 6.7</td><td>2.7 ± 0.4</td><td>84.3 ± 5.1</td><td>94.5 ± 8.3</td></tr><tr><td>BLEND-GARCH</td><td>60.4 ± 2.2</td><td>26.6 ± 2.6</td><td>2.9 ± 0.3</td><td>93.6 ± 1.8</td><td>79.4 ± 3.0</td></tr><tr><td>BLEND-LQR</td><td>59.2 ± 2.0</td><td>28.0 ± 2.4</td><td>3.1 ± 0.4</td><td>94.9 ± 1.6</td><td>79.1 ± 2.9</td></tr><tr><td>BLEND-QRF</td><td>55.5 ± 2.0</td><td>32.5 ± 2.5</td><td>1.7 ± 0.2</td><td>97.2 ± 1.8</td><td>75.7 ± 2.8</td></tr></tbody></table>

<table><caption>New proposed forecasting methods.</caption><thead><tr><th>Method</th><th>Forecasting technology</th><th>Input data</th></tr></thead><tbody><tr><td>BLEND-GARCH</td><td>ISE blended point fcst + GARCH</td><td>ground data+SAT+NWP</td></tr><tr><td>BLEND-LQR</td><td>ISE blended point fcst + LQR technique</td><td>ground data+SAT+NWP</td></tr><tr><td>BLEND-QRF</td><td>ISE blended point fcst + QRF technique</td><td>ground data+SAT+NWP</td></tr></tbody></table>

model of PIMENT cannot beat the climatological model CSD-CLIM for some sites namely CAR, MIL and TAB for certain forecast horizons.

For the site TAB which is located in the South of Spain and corresponds to a semi-arid climate (see Koppen-Geiger classification provided in Table 3) and which experiences a high share of clear skies PIMENT is outperformed by CSD-CLIM at 1 h lead time while at 2 h lead time, this is the case of LES and UU methods. Notice again the similar CRPS behavior of LES and UU models.

Finally, the interested reader is directed to Appendix where results related to the 8 sites are tabulated.

## 7. Impact of a high quality point forecast on the skills of the probabilistic methods

In the previous Section 6, we observed that the ISE forecasting methodology clearly outperforms the other methods. We hypothesize that the skill of the ISE method comes from the blended point forecasts which are used by the Analog Ensemble technique and not necessarily by the Analog technique by itself.

In order to confirm our assumption, we use the same blended point forecasts as inputs to three other approaches to generate probabilistic forecasts, including the PIMENT parametric GARCH approach, applied to ARMA point forecasts before. This new forecast is denoted BLEND-GARCH. In addition, we designed two other models based respectively on the LQR and QRF technique that use as input the blended ISE forecasts. These 2 new models are denoted BLEND-LQR and BLEND-QRF. Table 6 lists the new combinations of the ISE blended point forecast with the different techniques employed to generate the quantile forecasts. The results of the newly proposed methods are listed in the last 3 lines of Tables 4 and 5.

The combination BLEND-GARCH clearly improves the original PIMENT method for all the considered metrics. In particular, for the intra-hour scenario, the gain in average skill score is 9 points while for intra-day forecast, the gain in average forecast skill is 19 points.

The decomposition of the CRPS permits to highlight the improvement in resolution brought by the BLEND-GARCH combination.

Regardless of the scenario (intra-hour or intra-day forecasts), we can state that the BLEND-GARCH and BLEND-LQR exhibit similar performances. The same statement is also valid for BLEND-QRF and ISE forecasts.

Specifically, for intra-day forecasts, the BLEND-QRF slightly outperforms the original ISE method in terms of forecast skill.

Again, for a better inspection of the results, we provide the visual display of the metrics. Here also, in Fig. 6 the metrics are computed from the aggregation of forecasts and GHI data of all sites while Fig. 7 plots the CRPS obtained on each site. Notice that, for sake of comparison, the metrics related to the ISE and PIMENT previous methods are also plotted.

As shown by Fig. 6, irrespective of the forecast horizon and metric, a clear improvement is brought by using the blended point forecasts. For instance, Fig. 6e shows that the CRPS skill score of the new BLEND-GARCH now ranges from 58% to 22% instead of 56% to 1% obtained by the previous PIMENT method. The improvement is more pronounced at higher forecasting horizons and at the last forecast horizon the skill score of the previous PIMENT method based on ground measurements only is increased by almost 22 percentage points.

Also, combining the blended forecasts with techniques like LQR or QRF improves the skills of the probabilistic forecasts. From the decomposition of the CRPS into REL and RES, it appears that the BLEND-QRF slightly outperforms the ISE method in terms of resolution. However, the ISE method is still the best performer in terms of reliability.

Fig. 7 displays the CRPS obtained by the new combinations for each site under study. Now, all the new proposed forecasting techniques beat the CSD-CLIM model but for the site TAB (see Fig. 7g) the CSD-CLIM outperforms the BLEND-GARCH from a 2 h forecast horizon. Note that for the site TAB, the BLEND-QRF exhibits the best skill, a considerable improvement compared to the Analog Ensemble is found from two hour on-wards.
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** CRPS and its associated decomposition for all stations. Average GHI for all sites is 382.4 W m⁻² and can be used to calculate the relative counterparts of the different metrics. Notice that the same Y-scale is used for the CRPS and MAE plots to highlight the improvement brought by the probabilistic approach.</figcaption></figure>

From the previous results, we can conclude that the use of the blended point forecasts of participant ISE improves substantially the PIMENT parametric approach forecasting models. Further, the improvements are slightly better when the blended point forecasts are inputted to a nonlinear machine learning technique such as QRF.

In terms of reliability diagrams, the situation is also clearly improved when one compares Fig. 3f against 3b albeit it seems that BLEND-GARCH intra-hour forcecasts still suffer from a lack of reliability at high nominal proportions. Moreover, Figs. 3g and 3h reveal that the new proposed method BLEND-LQR and BLEND-QRF generate reliable forecasts.

## 8. Summary and conclusion

A benchmarking procedure set up by a group of experts of the IEA PVPS Task 16 was implemented to assess the performance of intra-hour and intra-day probabilistic solar irradiance forecasts. This procedure was utilized to evaluate eight distinct forecasting algorithms on eight sites located in Europe.

In the initial stage, the benchmarking exercise involved evaluating probabilistic forecasts submitted directly by five different participants. This initial comparison of forecasts using different input data and methods revealed a significant variation in performance. Particularly, during this first step, a probabilistic forecast that utilized a blended point
<figure><img src="image_7.png" /><figcaption>Fig. 5. CRPS for the different locations. Table 3 lists the average GHI of each site that can be used to compute the relative CRPS.</figcaption></figure>

forecast outperformed the other methods. In a second analysis, to better understand the impact of input data versus methodology on forecast quality, we combined the well-performing blended point forecast from the first comparison with different probabilistic approaches.

As mentioned earlier, the first stage of this work revealed that the initially proposed methods exhibit varying levels of forecast quality. In particular, the satellite-based method recently developed by OIE, which directly generates the set of quantile forecasts from a CMV model,

suffers from a lack of reliability that significantly impacts its overall performance. However, a calibration technique could be employed to enhance this attribute. As expected, a parametric approach, like the one proposed by PIMENT is not suitable to provide high quality probabilistic solar forecast, even with the high performing point forecast as input.

The second stage of this study demonstrated that a high quality point forecast (that blends measurements, satellite-based and NWP
<figure><img src="image_8.png" /><figcaption>(a) CRPS</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Reliability (REL)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(e) CRPS skill score</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) MAE</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(d) Resolution (RES)</figcaption></figure>

<caption>Fig. 6. CRPS and its associated decomposition for all stations and for the new proposed models.</caption>

forecasts) used in combination with a statistical technique is able to generate probabilistic forecasts with high quality. Overall, the skill scores of methods employing the blended point forecast vary between 42% and 46% for the intra-hour scenario and between 27% and 32% for the intra-day scenario. In contrast, methods that do not utilize the blended point forecast but are based on measurements and/or satellite data only exhibit skill scores ranging from 33% to 43% for intra-hour forecasts and from 8% to 16% for intra-day forecasts.

Besides a good forecast skill, the methodology that consists in gen-
erating probabilistic forecasts in a two step approach has the advantage
that it is easy to implement in combination with blended deterministic
forecasts that are well understood and used operationally. It allows
to benefit from high quality deterministic point forecasts with com-
paratively simple probabilistic techniques applied in a post-processing
step.

An alternative to blending the deterministic forecasts before apply-
ing the probabilistic techniques, would be to directly use the three
deterministic forecasts as input to these techniques. LQR or QRF can
be applied to generate quantile forecasts from different inputs in one
step. For the AnEn method the different inputs can be combined with
predictor weighting. Setting up these more complex models will be
subject of future investigations.

Ongoing efforts by members of the IEA PVPS Task 16 involve
the continuing development of solar probabilistic methods. Conse-
quently, the evaluation and comparison of probabilistic forecasts will
persist, and further analysis will be conducted using recent ground
measurement data, satellite or NWP data. Also, it is expected that
solar forecasting community should benefit from the recent advances
provided by the deep learning techniques. This continuing research
aims to enhance the quality of solar irradiance probabilistic forecasts.
<figure><img src="image_13.png" /><figcaption>(a) CAB</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) CAR</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) CEN</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) MIL</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(e) PAL</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(f) PAY</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(g) TAB</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(h) TOR</figcaption></figure>

<caption>Fig. 7. Improved CRPS by using the blending deterministic forecasts of Fraunhofer ISE for the different locations.</caption>
<caption>**Table 7**</caption>

Same as for Table 4 (intra-hour) but for site CAB. For site CAB, the CRPS of the CSD-CLIM is 81.8 W m⁻² and the uncertainty component UNC of the CRPS is 125.9 W m⁻².

<table><caption>**Table 8**</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>48.8 ± 6.1</td><td>40.3 ± 7.5</td><td>4.7 ± 0.4</td><td>81.8 ± 5.7</td><td>61.1 ± 8.5</td></tr><tr><td>PIMENT</td><td>50.9 ± 8.6</td><td>37.8 ± 10.5</td><td>2.1 ± 0.2</td><td>77.1 ± 8.4</td><td>69.2 ± 13.2</td></tr><tr><td>ISE</td><td>40.4 ± 5.3</td><td>50.7 ± 6.5</td><td>0.9 ± 0.2</td><td>86.4 ± 5.1</td><td>55.4 ± 7.6</td></tr><tr><td>UU</td><td>47.4 ± 7.8</td><td>42.1 ± 9.5</td><td>1.8 ± 0.4</td><td>80.3 ± 7.4</td><td>63.7 ± 11.4</td></tr><tr><td>LES</td><td>44.3 ± 7.4</td><td>45.9 ± 9.1</td><td>1.7 ± 0.5</td><td>83.2 ± 7.0</td><td>59.8 ± 10.2</td></tr><tr><td>BLEND-GARCH</td><td>41.9 ± 5.4</td><td>48.8 ± 6.6</td><td>1.9 ± 0.3</td><td>85.8 ± 5.1</td><td>55.0 ± 7.6</td></tr><tr><td>BLEND-LQR</td><td>42.2 ± 5.3</td><td>48.4 ± 6.5</td><td>2.1 ± 0.4</td><td>85.7 ± 4.9</td><td>54.7 ± 7.3</td></tr><tr><td>BLEND-QRF</td><td>39.9 ± 4.9</td><td>51.3 ± 6.0</td><td>1.1 ± 0.2</td><td>87.0 ± 4.7</td><td>54.3 ± 7.0</td></tr></tbody></table>

Same as for Table 5 (intra-day) but for site CAB.

<table><caption>**Table 9**</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>PIMENT</td><td>69.1 ± 3.8</td><td>15.6 ± 4.6</td><td>2.4 ± 0.2</td><td>59.1 ± 3.6</td><td>98.2 ± 6.1</td></tr><tr><td>ISE</td><td>50.2 ± 1.6</td><td>38.7 ± 1.9</td><td>1.0 ± 0.1</td><td>76.6 ± 1.6</td><td>69.4 ± 2.4</td></tr><tr><td>UU</td><td>68.4 ± 5.9</td><td>16.5 ± 7.2</td><td>3.7 ± 0.8</td><td>61.2 ± 5.1</td><td>96.3 ± 9.6</td></tr><tr><td>LES</td><td>65.6 ± 5.3</td><td>19.9 ± 6.4</td><td>3.7 ± 0.7</td><td>64.0 ± 4.6</td><td>89.8 ± 8.2</td></tr><tr><td>BLEND-GARCH</td><td>51.7 ± 1.5</td><td>36.9 ± 1.8</td><td>2.2 ± 0.1</td><td>76.4 ± 1.6</td><td>69.5 ± 2.2</td></tr><tr><td>BLEND-LQR</td><td>49.2 ± 1.6</td><td>39.9 ± 2.0</td><td>1.3 ± 0.1</td><td>78.0 ± 1.7</td><td>67.8 ± 2.5</td></tr><tr><td>BLEND-QRF</td><td>52.0 ± 1.5</td><td>36.5 ± 1.9</td><td>2.7 ± 0.0</td><td>76.6 ± 1.6</td><td>68.5 ± 2.4</td></tr><tr><td>BLEND-QRF</td><td>49.2 ± 1.6</td><td>39.9 ± 2.0</td><td>1.3 ± 0.1</td><td>78.0 ± 1.7</td><td>67.8 ± 2.5</td></tr></tbody></table>

Same as for Table 4 (intra-hour) but for site CAR. For site CAR, the CRPS of the CSD-CLIM is 73.2 W m⁻² and the uncertainty component UNC of the CRPS is 152.3 W m⁻².

<table><thead><tr><th>Method</th><th>CRPS</th><th>CRPSS</th><th>REL</th><th>RES</th><th>MAE</th></tr></thead><tbody><tr><td>OIE</td><td>49.2 ± 5.1</td><td>32.8 ± 7.0</td><td>4.3 ± 0.4</td><td>108.5 ± 5.1</td><td>63.0 ± 7.4</td></tr><tr><td>PIMENT</td><td>53.3 ± 10.9</td><td>27.2 ± 14.9</td><td>5.3 ± 0.7</td><td>105.4 ± 10.2</td><td>68.8 ± 16.3</td></tr><tr><td>ISE</td><td>39.6 ± 7.2</td><td>46.0 ± 9.9</td><td>1.1 ± 0.3</td><td>114.9 ± 7.0</td><td>53.3 ± 10.0</td></tr><tr><td>UU</td><td>43.8 ± 8.9</td><td>40.2 ± 12.2</td><td>1.1 ± 0.2</td><td>110.7 ± 8.8</td><td>58.0 ± 12.9</td></tr><tr><td>LES</td><td>43.1 ± 8.7</td><td>41.1 ± 12.0</td><td>1.7 ± 0.3</td><td>112.0 ± 8.5</td><td>57.3 ± 12.6</td></tr><tr><td>BLEND-GARCH</td><td>44.3 ± 7.4</td><td>39.5 ± 10.2</td><td>2.9 ± 0.4</td><td>112.0 ± 7.0</td><td>55.8 ± 9.5</td></tr><tr><td>BLEND-LQR</td><td>43.1 ± 6.6</td><td>41.2 ± 9.0</td><td>2.9 ± 0.2</td><td>113.2 ± 6.3</td><td>55.0 ± 9.6</td></tr><tr><td>BLEND-QRF</td><td>40.0 ± 6.2</td><td>45.4 ± 8.5</td><td>1.3 ± 0.2</td><td>114.7 ± 6.1</td><td>53.2 ± 8.9</td></tr></tbody></table>

## CRediT authorship contribution statement

**Philippe Lauret:** Writing – review & editing, Writing – original draft, Supervision, Software, Methodology. **Rodrigo Alonso-Suárez:** Writing – original draft, Validation, Software, Investigation. **Rodrigo Amaro e Silva:** Writing – review & editing, Writing – original draft, Validation, Supervision, Software, Investigation. **John Boland:** Writing – review & editing, Writing – original draft. **Mathieu David:** Writing – review & editing, Writing – original draft, Validation, Software, Investigation. **Wiebke Herzberg:** Writing – review & editing, Software, Investigation. **Josselin Le Gall La Salle:** Writing – original draft, Software. **Elke Lorenz:** Writing – review & editing, Writing – original draft, Validation, Supervision. **Lennard Visser:** Writing – review & editing, Writing – original draft, Validation, Software, Investigation. **Wilfried van Sark:** Supervision. **Tobias Zech:** Writing – review & editing, Writing – original draft, Validation, Software, Investigation.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This research paper contributes to the IEA PVPS task 16 “Task 16 Solar resource for high penetration and large-scale applications”, which is a collaborative effort between the International Energy Agency’s. The study also received support from the TwInSolar project funded by the

European Union’s Horizon Europe research and innovation program grant number 101076447.

R. Alonso-Suárez acknowledges financial support from the CSIC Group’s Program, Universidad de la República, Uruguay.

The research contribution of Tobias Zech, Wiebke Herzberg, Elke Lorenz (Fraunhofer ISE) received funding by the German Federal Ministry for Economic Affairs and Climate Action (BMWK) within the SOLREV project (Grant Agreement no. 03EE1010) on the basis of a decision by the German Bundestag.

The authors express their gratitude to Anna Forstinger for providing the ground datasets.

## Appendix A. Intra-hour and intra-day forecasts results site CAB

See Tables 7 and 8.

## Appendix B. Intra-hour and intra-day forecasts results site CAR

See Tables 9 and 10.

## Appendix C. Intra-hour and intra-day forecasts results site CEN

See Tables 11 and 12.

## Appendix D. Intra-hour and intra-day forecasts results site MIL

See Tables 13 and 14.
<table><caption>**Table 10**
Same as for Table 5 (intra-day) but for site CAR.</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>PIMENT</td><td>74.9 ± 4.4</td><td>-2.3 ± 6.0</td><td>6.3 ± 0.2</td><td>84.9 ± 4.2</td><td>102.2 ± 7.0</td></tr><tr><td>ISE</td><td>52.6 ± 2.0</td><td>28.2 ± 2.7</td><td>1.8 ± 0.2</td><td>102.7 ± 1.8</td><td>71.0 ± 2.6</td></tr><tr><td>UU</td><td>65.2 ± 5.6</td><td>10.9 ± 7.7</td><td>1.3 ± 0.1</td><td>89.5 ± 5.7</td><td>90.8 ± 8.9</td></tr><tr><td>LES</td><td>64.7 ± 5.5</td><td>11.6 ± 7.5</td><td>2.5 ± 0.2</td><td>91.1 ± 5.3</td><td>88.1 ± 8.0</td></tr><tr><td>BLEND-GARCH</td><td>57.7 ± 2.0</td><td>21.2 ± 2.8</td><td>3.7 ± 0.2</td><td>99.4 ± 1.8</td><td>73.7 ± 2.4</td></tr><tr><td>BLEND-LQR</td><td>55.3 ± 1.8</td><td>24.4 ± 2.5</td><td>4.1 ± 0.4</td><td>102.2 ± 1.4</td><td>73.6 ± 2.6</td></tr><tr><td>BLEND-QRF</td><td>51.7 ± 1.9</td><td>29.4 ± 2.6</td><td>2.1 ± 0.3</td><td>103.8 ± 1.6</td><td>69.7 ± 2.6</td></tr></tbody></table>

<table><caption>**Table 11**
Same as for Table 4 (intra-hour) but for site CEN. For site CEN, the CRPS of the CSD-CLIM is 87.7 W m⁻² and the uncertainty component UNC of the CRPS is 153.1 W m⁻².</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>60.9 ± 6.4</td><td>30.6 ± 7.3</td><td>5.5 ± 0.3</td><td>96.2 ± 6.5</td><td>78.4 ± 10.1</td></tr><tr><td>PIMENT</td><td>60.9 ± 10.6</td><td>30.5 ± 12.1</td><td>4.0 ± 0.7</td><td>94.6 ± 9.9</td><td>82.5 ± 16.2</td></tr><tr><td>ISE</td><td>50.6 ± 7.4</td><td>42.3 ± 8.5</td><td>1.0 ± 0.2</td><td>102.0 ± 7.2</td><td>69.4 ± 10.9</td></tr><tr><td>UU</td><td>53.7 ± 8.7</td><td>38.7 ± 9.9</td><td>0.9 ± 0.1</td><td>98.7 ± 8.6</td><td>73.1 ± 13.2</td></tr><tr><td>LES</td><td>53.1 ± 9.0</td><td>39.5 ± 10.2</td><td>1.2 ± 0.3</td><td>99.7 ± 8.7</td><td>72.2 ± 13.0</td></tr><tr><td>BLEND-GARCH</td><td>53.8 ± 7.8</td><td>38.7 ± 8.9</td><td>1.9 ± 0.3</td><td>99.7 ± 7.5</td><td>72.0 ± 11.0</td></tr><tr><td>BLEND-LQR</td><td>53.8 ± 7.5</td><td>38.6 ± 8.6</td><td>2.4 ± 0.5</td><td>100.1 ± 7.1</td><td>71.4 ± 11.3</td></tr><tr><td>BLEND-QRF</td><td>51.0 ± 6.9</td><td>41.8 ± 7.9</td><td>0.8 ± 0.1</td><td>101.4 ± 6.8</td><td>70.3 ± 10.4</td></tr></tbody></table>

<table><caption>**Table 12**
Same as for Table 5 (intra-day) but for site CEN.</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>PIMENT</td><td>82.8 ± 4.2</td><td>5.6 ± 4.8</td><td>5.3 ± 0.1</td><td>74.1 ± 4.1</td><td>116.6 ± 6.6</td></tr><tr><td>ISE</td><td>65.2 ± 2.6</td><td>25.6 ± 2.9</td><td>1.7 ± 0.2</td><td>88.1 ± 2.4</td><td>90.6 ± 3.5</td></tr><tr><td>UU</td><td>75.1 ± 5.6</td><td>14.4 ± 6.4</td><td>0.8 ± 0.1</td><td>77.3 ± 5.7</td><td>107.9 ± 9.4</td></tr><tr><td>LES</td><td>75.4 ± 5.3</td><td>14.1 ± 6.0</td><td>2.1 ± 0.2</td><td>78.2 ± 5.0</td><td>105.8 ± 8.6</td></tr><tr><td>BLEND-GARCH</td><td>69.2 ± 2.6</td><td>21.1 ± 3.0</td><td>2.7 ± 0.2</td><td>85.0 ± 2.4</td><td>94.0 ± 3.5</td></tr><tr><td>BLEND-LQR</td><td>68.3 ± 2.3</td><td>22.2 ± 2.7</td><td>3.2 ± 0.2</td><td>86.5 ± 2.1</td><td>94.0 ± 3.7</td></tr><tr><td>BLEND-QRF</td><td>65.2 ± 2.6</td><td>25.7 ± 3.0</td><td>1.1 ± 0.1</td><td>87.4 ± 2.6</td><td>91.9 ± 3.9</td></tr></tbody></table>

<table><caption>**Table 13**
Same as for Table 4 (intra-hour) but for site MIL. For site MIL, the CRPS of the CSD-CLIM is 78.4 W m⁻² and the uncertainty component UNC of the CRPS is 154.6 W m⁻².</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>50.2 ± 4.9</td><td>36.0 ± 6.2</td><td>4.2 ± 0.4</td><td>109.6 ± 4.6</td><td>63.8 ± 7.4</td></tr><tr><td>PIMENT</td><td>53.7 ± 10.5</td><td>31.5 ± 13.3</td><td>6.3 ± 1.1</td><td>108.2 ± 9.4</td><td>66.5 ± 15.2</td></tr><tr><td>ISE</td><td>41.1 ± 7.0</td><td>47.5 ± 8.9</td><td>1.4 ± 0.4</td><td>115.8 ± 6.6</td><td>54.6 ± 9.7</td></tr><tr><td>UU</td><td>43.9 ± 8.4</td><td>44.0 ± 10.7</td><td>0.9 ± 0.1</td><td>112.6 ± 8.3</td><td>58.3 ± 12.2</td></tr><tr><td>LES</td><td>44.8 ± 8.2</td><td>42.8 ± 10.4</td><td>3.8 ± 0.5</td><td>114.6 ± 7.7</td><td>56.1 ± 11.3</td></tr><tr><td>BLEND-GARCH</td><td>46.8 ± 7.0</td><td>40.2 ± 9.0</td><td>4.3 ± 0.6</td><td>113.1 ± 6.4</td><td>57.2 ± 9.2</td></tr><tr><td>BLEND-LQR</td><td>46.0 ± 6.7</td><td>41.3 ± 8.6</td><td>4.4 ± 0.8</td><td>114.0 ± 5.9</td><td>58.8 ± 9.3</td></tr><tr><td>BLEND-QRF</td><td>42.9 ± 6.0</td><td>45.2 ± 7.6</td><td>2.6 ± 0.5</td><td>115.3 ± 5.5</td><td>56.3 ± 8.2</td></tr></tbody></table>

<table><caption>**Table 14**
Same as for Table 5 (intra-day) but for site MIL.</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>PIMENT</td><td>75.7 ± 4.6</td><td>3.4 ± 5.8</td><td>7.8 ± 0.1</td><td>87.8 ± 4.5</td><td>100.2 ± 7.2</td></tr><tr><td>ISE</td><td>54.9 ± 2.6</td><td>30.0 ± 3.3</td><td>3.2 ± 0.6</td><td>104.0 ± 2.0</td><td>73.4 ± 3.2</td></tr><tr><td>UU</td><td>65.4 ± 6.4</td><td>16.5 ± 8.2</td><td>0.9 ± 0.1</td><td>91.1 ± 6.3</td><td>91.0 ± 10.0</td></tr><tr><td>LES</td><td>67.0 ± 6.2</td><td>14.5 ± 7.9</td><td>5.5 ± 0.6</td><td>94.1 ± 5.6</td><td>86.5 ± 9.0</td></tr><tr><td>BLEND-GARCH</td><td>61.0 ± 2.8</td><td>22.2 ± 3.5</td><td>6.5 ± 0.7</td><td>101.1 ± 2.1</td><td>77.4 ± 4.3</td></tr><tr><td>BLEND-LQR</td><td>60.2 ± 2.8</td><td>23.2 ± 3.6</td><td>8.0 ± 1.3</td><td>103.4 ± 1.6</td><td>79.9 ± 4.6</td></tr><tr><td>BLEND-QRF</td><td>55.1 ± 2.4</td><td>29.8 ± 3.0</td><td>4.3 ± 0.4</td><td>104.8 ± 2.0</td><td>72.7 ± 3.1</td></tr></tbody></table>

## Appendix E. Intra-hour and intra-day forecasts results site PAL

See Tables 15 and 16.

## Appendix F. Intra-hour and intra-day forecasts results site PAY

See Tables 17 and 18.

## Appendix G. Intra-hour and intra-day forecasts results site TAB

See Tables 19 and 20.

## Appendix H. Intra-hour and intra-day forecasts results site TOR

See Tables 21 and 22.
<caption>**Table 15**</caption>

Same as for Table 4 (intra-hour) but for site PAL. For site PAL, the CRPS of the CSD-CLIM is 88.0 W m⁻² and the uncertainty component UNC of the CRPS is 140.5 W m⁻².

<table><caption>**Table 16**</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>55.7 ± 6.5</td><td>36.7 ± 7.3</td><td>5.3 ± 0.5</td><td>92.4 ± 6.0</td><td>69.3 ± 9.1</td></tr><tr><td>PIMENT</td><td>56.5 ± 9.2</td><td>35.8 ± 10.5</td><td>2.9 ± 0.4</td><td>89.3 ± 8.9</td><td>76.9 ± 14.3</td></tr><tr><td>ISE</td><td>45.0 ± 5.7</td><td>48.9 ± 6.5</td><td>1.2 ± 0.2</td><td>99.0 ± 5.5</td><td>61.3 ± 8.2</td></tr><tr><td>UU</td><td>52.9 ± 7.8</td><td>39.9 ± 8.9</td><td>2.0 ± 0.3</td><td>92.0 ± 7.6</td><td>70.3 ± 11.6</td></tr><tr><td>LES</td><td>50.0 ± 8.5</td><td>43.2 ± 9.6</td><td>1.7 ± 0.4</td><td>94.6 ± 8.1</td><td>68.0 ± 11.9</td></tr><tr><td>BLEND-GARCH</td><td>47.3 ± 5.8</td><td>46.3 ± 6.6</td><td>1.8 ± 0.2</td><td>97.4 ± 5.7</td><td>62.3 ± 8.1</td></tr><tr><td>BLEND-LQR</td><td>47.9 ± 5.7</td><td>45.5 ± 6.5</td><td>2.2 ± 0.3</td><td>97.1 ± 5.4</td><td>62.2 ± 8.1</td></tr><tr><td>BLEND-QRF</td><td>45.2 ± 5.4</td><td>48.6 ± 6.1</td><td>1.3 ± 0.2</td><td>98.9 ± 5.2</td><td>61.4 ± 7.8</td></tr></tbody></table>

Same as for Table 5 (intra-day) but for site PAL.

<table><caption>**Table 17**</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>PIMENT</td><td>76.9 ± 4.5</td><td>12.6 ± 5.1</td><td>3.6 ± 0.2</td><td>69.6 ± 4.3</td><td>109.4 ± 7.3</td></tr><tr><td>ISE</td><td>55.1 ± 1.4</td><td>37.4 ± 1.6</td><td>1.4 ± 0.1</td><td>89.2 ± 1.5</td><td>75.9 ± 2.0</td></tr><tr><td>UU</td><td>73.1 ± 5.8</td><td>16.9 ± 6.6</td><td>2.7 ± 0.3</td><td>72.5 ± 5.5</td><td>102.6 ± 9.8</td></tr><tr><td>LES</td><td>73.7 ± 6.0</td><td>16.2 ± 6.8</td><td>3.4 ± 0.6</td><td>72.5 ± 5.4</td><td>102.3 ± 9.4</td></tr><tr><td>BLEND-GARCH</td><td>57.3 ± 1.6</td><td>34.9 ± 1.8</td><td>2.2 ± 0.2</td><td>87.8 ± 1.4</td><td>76.9 ± 2.3</td></tr><tr><td>BLEND-LQR</td><td>58.0 ± 1.5</td><td>34.1 ± 1.8</td><td>3.2 ± 0.3</td><td>88.0 ± 1.2</td><td>76.9 ± 2.3</td></tr><tr><td>BLEND-QRF</td><td>54.9 ± 1.7</td><td>37.6 ± 1.9</td><td>1.8 ± 0.2</td><td>89.7 ± 1.5</td><td>75.4 ± 2.3</td></tr></tbody></table>

Same as for Table 4 (intra-hour) but for site PAY. For site PAY, the CRPS of the CSD-CLIM is 89.0 W m⁻² and the uncertainty component UNC of the CRPS is 151.6 W m⁻².

<table><caption>**Table 18**</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>53.3 ± 5.7</td><td>40.1 ± 6.4</td><td>4.7 ± 0.3</td><td>105.3 ± 5.8</td><td>68.4 ± 9.0</td></tr><tr><td>PIMENT</td><td>57.4 ± 11.7</td><td>35.5 ± 13.1</td><td>4.6 ± 0.6</td><td>101.2 ± 11.1</td><td>75.2 ± 17.7</td></tr><tr><td>ISE</td><td>44.7 ± 7.6</td><td>49.8 ± 8.5</td><td>0.9 ± 0.1</td><td>110.1 ± 7.4</td><td>61.1 ± 11.0</td></tr><tr><td>UU</td><td>48.0 ± 9.3</td><td>46.1 ± 10.5</td><td>1.3 ± 0.3</td><td>107.2 ± 9.1</td><td>64.0 ± 13.3</td></tr><tr><td>LES</td><td>48.2 ± 9.7</td><td>45.8 ± 10.9</td><td>2.0 ± 0.5</td><td>107.7 ± 9.2</td><td>64.3 ± 13.8</td></tr><tr><td>BLEND-GARCH</td><td>50.1 ± 8.0</td><td>43.7 ± 9.0</td><td>3.0 ± 0.3</td><td>106.8 ± 7.8</td><td>64.9 ± 11.2</td></tr><tr><td>BLEND-LQR</td><td>49.5 ± 7.7</td><td>44.4 ± 8.6</td><td>3.4 ± 0.6</td><td>107.8 ± 7.1</td><td>64.6 ± 11.2</td></tr><tr><td>BLEND-QRF</td><td>45.6 ± 6.8</td><td>48.7 ± 7.6</td><td>1.8 ± 0.2</td><td>110.1 ± 6.6</td><td>61.7 ± 10.1</td></tr></tbody></table>

Same as for Table 5 (intra-day) but for site PAY.

<table><caption>**Table 19**</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>PIMENT</td><td>81.1 ± 4.3</td><td>8.8 ± 4.8</td><td>5.5 ± 0.1</td><td>78.2 ± 4.2</td><td>112.0 ± 6.7</td></tr><tr><td>ISE</td><td>60.7 ± 3.2</td><td>31.8 ± 3.6</td><td>1.4 ± 0.2</td><td>94.6 ± 3.0</td><td>84.7 ± 4.5</td></tr><tr><td>UU</td><td>71.9 ± 6.0</td><td>19.2 ± 6.7</td><td>1.9 ± 0.2</td><td>83.9 ± 5.8</td><td>100.0 ± 9.4</td></tr><tr><td>LES</td><td>73.8 ± 6.3</td><td>17.1 ± 7.1</td><td>4.6 ± 1.1</td><td>84.7 ± 5.2</td><td>101.2 ± 9.4</td></tr><tr><td>BLEND-GARCH</td><td>66.8 ± 3.2</td><td>24.9 ± 3.6</td><td>4.1 ± 0.4</td><td>91.2 ± 2.8</td><td>88.4 ± 4.3</td></tr><tr><td>BLEND-LQR</td><td>65.8 ± 3.1</td><td>26.1 ± 3.5</td><td>5.4 ± 0.6</td><td>93.5 ± 2.5</td><td>88.4 ± 4.4</td></tr><tr><td>BLEND-QRF</td><td>61.2 ± 3.4</td><td>31.2 ± 3.9</td><td>2.7 ± 0.3</td><td>95.4 ± 3.2</td><td>84.3 ± 4.8</td></tr></tbody></table>

Same as for Table 4 (intra-hour) but for site TAB. For site TAB, the CRPS of the CSD-CLIM is 58.7 W m⁻² and the uncertainty component UNC of the CRPS is 160.0 W m⁻².

<table><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>56.3 ± 4.9</td><td>4.1 ± 8.4</td><td>6.2 ± 0.8</td><td>111.3 ± 4.2</td><td>70.2 ± 7.1</td></tr><tr><td>PIMENT</td><td>55.9 ± 10.6</td><td>4.9 ± 18.1</td><td>6.0 ± 1.2</td><td>111.5 ± 9.5</td><td>72.7 ± 16.0</td></tr><tr><td>ISE</td><td>45.9 ± 7.7</td><td>21.9 ± 13.1</td><td>1.6 ± 0.5</td><td>117.2 ± 7.2</td><td>59.9 ± 10.5</td></tr><tr><td>UU</td><td>47.7 ± 8.3</td><td>18.7 ± 14.2</td><td>1.4 ± 0.3</td><td>115.1 ± 8.0</td><td>61.2 ± 11.5</td></tr><tr><td>LES</td><td>46.2 ± 8.9</td><td>21.4 ± 15.1</td><td>1.2 ± 0.2</td><td>116.4 ± 8.7</td><td>61.2 ± 12.6</td></tr><tr><td>BLEND-GARCH</td><td>51.5 ± 7.8</td><td>12.3 ± 13.3</td><td>3.7 ± 0.9</td><td>113.7 ± 6.9</td><td>63.9 ± 9.6</td></tr><tr><td>BLEND-LQR</td><td>48.7 ± 6.3</td><td>17.1 ± 10.7</td><td>3.1 ± 0.3</td><td>115.8 ± 6.0</td><td>63.6 ± 9.1</td></tr><tr><td>BLEND-QRF</td><td>46.1 ± 6.3</td><td>21.5 ± 10.8</td><td>1.6 ± 0.3</td><td>116.9 ± 6.1</td><td>60.5 ± 8.8</td></tr></tbody></table>
<table><caption>**Table 20**
Same as for Table 5 (intra-day) but for site TAB.</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>PIMENT</td><td>73.8 ± 2.6</td><td>-25.7 ± 4.3</td><td>7.8 ± 0.4</td><td>95.5 ± 2.2</td><td>101.2 ± 4.2</td></tr><tr><td>ISE</td><td>56.2 ± 1.2</td><td>4.2 ± 2.0</td><td>3.4 ± 0.6</td><td>108.6 ± 0.5</td><td>73.9 ± 1.2</td></tr><tr><td>UU</td><td>64.9 ± 3.3</td><td>-10.5 ± 5.7</td><td>1.7 ± 0.1</td><td>98.2 ± 3.4</td><td>84.1 ± 4.2</td></tr><tr><td>LES</td><td>65.5 ± 3.7</td><td>-11.6 ± 6.4</td><td>2.1 ± 0.3</td><td>98.0 ± 3.4</td><td>87.4 ± 4.6</td></tr><tr><td>BLEND-GARCH</td><td>63.0 ± 1.6</td><td>-7.2 ± 2.7</td><td>6.4 ± 0.9</td><td>104.9 ± 0.7</td><td>80.1 ± 2.5</td></tr><tr><td>BLEND-LQR</td><td>57.0 ± 0.8</td><td>2.9 ± 1.4</td><td>4.1 ± 0.5</td><td>108.5 ± 0.4</td><td>76.6 ± 1.3</td></tr><tr><td>BLEND-QRF</td><td>53.8 ± 0.6</td><td>8.3 ± 1.1</td><td>2.3 ± 0.2</td><td>109.9 ± 0.5</td><td>71.7 ± 0.8</td></tr></tbody></table>

<table><caption>**Table 21**
Same as for Table 4 (intra-hour) but for site TOR. For site TOR, the CRPS of the CSD-CLIM is 84.5 W m⁻² and the uncertainty component UNC of the CRPS is 125.9 W m⁻².</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>55.4 ± 5.3</td><td>34.5 ± 6.3</td><td>7.3 ± 0.4</td><td>79.9 ± 4.9</td><td>67.3 ± 7.6</td></tr><tr><td>PIMENT</td><td>52.4 ± 9.0</td><td>38.0 ± 10.7</td><td>2.6 ± 0.2</td><td>78.1 ± 8.8</td><td>70.0 ± 13.4</td></tr><tr><td>ISE</td><td>43.0 ± 4.8</td><td>49.1 ± 5.6</td><td>0.9 ± 0.0</td><td>85.8 ± 4.8</td><td>59.3 ± 7.0</td></tr><tr><td>UU</td><td>48.8 ± 7.8</td><td>42.2 ± 9.3</td><td>2.3 ± 0.6</td><td>81.4 ± 7.3</td><td>64.4 ± 11.4</td></tr><tr><td>LES</td><td>45.9 ± 7.9</td><td>45.7 ± 9.3</td><td>1.7 ± 0.5</td><td>83.7 ± 7.4</td><td>62.1 ± 10.7</td></tr><tr><td>BLEND-GARCH</td><td>44.8 ± 5.2</td><td>47.0 ± 6.2</td><td>1.7 ± 0.1</td><td>84.8 ± 5.1</td><td>60.0 ± 7.2</td></tr><tr><td>BLEND-LQR</td><td>46.1 ± 5.2</td><td>45.5 ± 6.2</td><td>2.8 ± 0.4</td><td>84.7 ± 4.9</td><td>59.4 ± 7.3</td></tr><tr><td>BLEND-QRF</td><td>41.9 ± 4.6</td><td>50.4 ± 5.4</td><td>1.4 ± 0.1</td><td>87.4 ± 4.5</td><td>56.9 ± 6.9</td></tr></tbody></table>

<table><caption>**Table 22**
Same as for Table 5 (intra-day) but for site TOR.</caption><thead><tr><th>Method</th><th>CRPS (W/m²)</th><th>CRPSS (%)</th><th>REL (W/m²)</th><th>RES (W/m²)</th><th>MAE (W/m²)</th></tr></thead><tbody><tr><td>OIE</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>PIMENT</td><td>72.0 ± 4.1</td><td>14.8 ± 4.8</td><td>2.8 ± 0.1</td><td>58.8 ± 4.0</td><td>100.9 ± 6.8</td></tr><tr><td>ISE</td><td>51.8 ± 1.6</td><td>38.7 ± 1.9</td><td>1.1 ± 0.1</td><td>77.2 ± 1.6</td><td>72.0 ± 2.3</td></tr><tr><td>UU</td><td>70.7 ± 6.3</td><td>16.3 ± 7.4</td><td>5.1 ± 1.0</td><td>62.3 ± 5.3</td><td>98.3 ± 10.9</td></tr><tr><td>LES</td><td>69.1 ± 6.1</td><td>18.3 ± 7.2</td><td>4.6 ± 1.0</td><td>63.5 ± 5.1</td><td>93.7 ± 9.4</td></tr><tr><td>BLEND-GARCH</td><td>54.7 ± 1.9</td><td>35.3 ± 2.2</td><td>2.5 ± 0.3</td><td>75.7 ± 1.6</td><td>72.6 ± 2.3</td></tr><tr><td>BLEND-LQR</td><td>56.2 ± 1.9</td><td>33.4 ± 2.2</td><td>4.2 ± 0.4</td><td>75.9 ± 1.5</td><td>72.4 ± 2.2</td></tr><tr><td>BLEND-QRF</td><td>51.6 ± 2.0</td><td>39.0 ± 2.4</td><td>2.2 ± 0.2</td><td>78.5 ± 1.8</td><td>70.2 ± 2.6</td></tr></tbody></table>

## References

[1] G. Notton, M.-L. Nivet, C. Voyant, C. Paoli, C. Darras, F. Motte, A. Fouilloy, Intermittent and stochastic character of renewable energy sources: Consequences, cost of intermittence and benefit of forecasting, Renew. Sustain. Energy Rev. 87 (2018) 96–105.

[2] E. Lorenz, J.A. Ruiz-Arias, L. Martin, S. Wilbert, C. Köhler, R. Fritz, A. Betti, P. Lauret, M. David, J. Huang, R. Perez, A. Kazantzidis, P. Wang, Y.-M. Saint-Drenan, Forecasting solar radiation and photovoltaic power, in: Best Practices Handbook for the Collection and Use of Solar Resource Data for Solar Energy Applications: Third Edition, (NREL/TP-5D00-77635) National Renewable Energy Laboratory, Golden, CO, 2021.

[3] J. Antonanzas, N. Osorio, R. Escobar, R. Urraca, F. Martinez-de Pison, F. Antonanzas-Torres, Review of photovoltaic power forecasting, Sol. Energy 136 (2016) 78–111.

[4] S. Sobri, S. Koohi-Kamali, N.A. Rahim, Solar photovoltaic generation forecasting methods: A review, Energy Convers. Manage. 156 (2018) 459–497.

[5] R. Blaga, A. Sabadus, N. Stefu, C. Dughir, M. Paulescu, V. Badescu, A current perspective on the accuracy of incoming solar energy forecasting, Prog. Energy Combust. Sci. 70 (2019) 119–144.

[6] D. Yang, W. Wang, C.A. Gueymard, T. Hong, J. Kleissl, J. Huang, M.J. Perez, R. Perez, J.M. Bright, Á.X. Xiang, D. van der Meer, I.M. Peters, A review of solar forecasting, its dependence on atmospheric sciences and implications for grid integration: Towards carbon neutrality, Renew. Sustain. Energy Rev. 161 (2022) 112348.

[7] H. Verbois, Y.-M. Saint-Drenan, A. Thiery, P. Blanc, Statistical learning for NWP post-processing: A benchmark for solar irradiance forecasting, Sol. Energy 238 (2022) 132–149, URL https://www.sciencedirect.com/science/article/pii/S0038092X22001839.

[8] R. Alonso-Suárez, D. Aicardi, F. Marchesoni-Acland, Analysis of persistence-based solar irradiance forecasting benchmarks, 2022, arXiv:2203.13819.

[9] H.T. Pedro, C.F. Coimbra, Assessment of forecasting techniques for solar power production with no exogenous inputs, Sol. Energy 86 (7) (2012) 2017–2028.

[10] P. Lauret, C. Voyant, T. Soubdhan, M. David, P. Poggi, A benchmarking of machine learning techniques for solar radiation forecasting in an insular context, Sol. Energy 112 (2015) 446–457.

[11] D. Aicardi, P. Musé, R. Alonso-Suárez, A comparison of satellite cloud motion vectors techniques to forecast intra-day hourly solar global horizontal irradiation, Sol. Energy 233 (2022) 46–60.

[12] Q. Paletta, G. Arbod, J. Lasenby, Benchmarking of deep learning irradiance forecasting models from sky images, Aì an an in-depth analysis, Sol. Energy 224 (2021) 855–867.

[13] IEA-SHC-T46, Solar resource assessment and forecasting, 2024, https://task46.iea-shc.org. (Accessed: 22 January 2024).

[14] E. Lorenz, J. Remund, S.C. Müller, W. Traunmüller, G. Steinmaurer, D. Pozo, J.A. Ruiz-Arias, V.L. Fanego, L. Ramirez, M.G. Romeo, C. Kurz, L.M. Pomares, C. Geijo Guerrero, Benchmarking of different approaches to forecast solar irradiance, in: 24th European Photovoltaic Solar Energy Conference, Hamburg Germany, 2009, pp. 21–25.

[15] J.M. Morales, A.J. Conejo, H. Madsen, P. Pinson, M. Zugno, Integrating renewables in electricity markets, in: International Series in Operations Research & Management Science, vol. 205, Springer US, 2014.

[16] E.B. Iversen, J.M. Morales, J.K. Møller, H. Madsen, Short-term probabilistic forecasting of wind speed using stochastic differential equations, Int. J. Forecast. (2015).

[17] J. Jung, R.P. Broadwater, Current status and future advances for wind speed and power forecasting, Renew. Sustain. Energy Rev. 31 (2014) 762–777.

[18] P. Pinson, H.A. Nielsen, J.K. Møller, H. Madsen, G.N. Kariniotakis, Non-parametric probabilistic forecasts of wind power: required properties and evaluation, Wind Energy 10 (6) (2007) 497–516.

[19] T. Hong, P. Pinson, S. Fan, H. Zareipour, A. Troccoli, R.J. Hyndman, Probabilistic energy forecasting: Global energy forecasting competition 2014 and beyond, Int. J. Forecast. 32 (3) (2016) 896–913.

[20] D. van der Meer, J. Widén, J. Munkhammar, Review on probabilistic forecasting of photovoltaic power production and electricity consumption, Renew. Sustain. Energy Rev. 81 (2018) 1484–1512, URL https://www.sciencedirect.com/science/article/pii/S1364032117308523.

[21] T. Hong, P. Pinson, Y. Wang, R. Weron, D. Yang, H. Zareipour, Energy forecasting: A review and outlook, IEEE Open Access J. Power Energy 7 (2020) 376–388.

[22] A. Grantham, Y.R. Gel, J. Boland, Nonparametric short-term probabilistic forecasting for solar radiation, Sol. Energy 133 (2016) 465–475.
[23] M. David, F. Ramahatana, P. Trombe, P. Lauret, Probabilistic forecasting of the solar irradiance with recursive ARMA and GARCH models, Sol. Energy 133 (2016) 55-72.

[24] P. Lauret, M. David, H. Pedro, Probabilistic solar forecasting using quantile regression Models, Energies 10 (10) (2017) 1591.

[25] M. David, L. Mazorra Aguiar, P. Lauret, Comparison of intraday probabilistic forecasting of solar irradiance using only endogenous data, Int. J. Forecast. 34 (3) (2018) 529-547.

[26] R. Alonso-Suárez, M. David, V. Branco, P. Lauret, Intra-day solar probabilistic forecasts including local short-term variability and satellite information, Renew. Energy 158 (2020) 554-573.

[27] L. Mazorra-Aguiar, P. Lauret, M. David, A. Oliver, G. Montero, Comparison of two solar probabilistic forecasting methodologies for microgrids energy efficiency, Energies 14 (6) (2021) 1679.

[28] D. Yang, D. van der Meer, J. Munkhammar, Probabilistic solar forecasting benchmarks on a standardized dataset at folsom, California, Sol. Energy 206 (2020) 628-639.

[29] H.T.C. Pedro, D.P. Larson, C.F.M. Coimbra, A comprehensive dataset for the accelerated development and benchmarking of solar forecasting methods, J. Renew. Sustain. Energy 11 (3) (2019) 036102.

[30] T. Carrière, R. Amaro E Silva, F. Zhuang, Y.-M. Saint-Drenan, P. Blanc, A new approach for satellite-based probabilistic solar forecasting with cloud motion vectors, Energies 14 (16) (2021) 4951.

[31] Q. Paletta, G. Arbod, J. Lasenby, Omnivision forecasting: Combining satellite and sky images for improved deterministic and probabilistic intra-hour solar energy predictions, Appl. Energy 336 (2023) 120818.

[32] IEA-PVPS-T16, Solar resource for high penetration and large scale applications, 2024, https://iea-pvps.org/research-tasks/solar-resource-for-high-penetration-and-large-scale-applications/. (Accessed: 22 January 2024).

[33] P. Lauret, M. David, P. Pinson, Verification of solar irradiance probabilistic forecasts, Sol. Energy 194 (2019) 254-271.

[34] J. Le Gal La Salle, M. David, P. Lauret, A new climatology reference model to benchmark probabilistic solar forecasts, Sol. Energy 223 (2021) 398-414.

[35] A. Forstinger, S. Wilbert, A.R. Jensen, B. Kraas, C.F. Peruchena, C.A. Gueymard, D. Ronzio, D. Yang, E. Collino, J.P. Martinez, J.A. Ruiz-Arias, N. Hanrieder, P. Blanc, Y.M. Saint-Drenan, Expert quality control of solar radiation ground data sets, in: Proceedings of SWC 2021: ISES Solar World Congress, International Solar Energy Society, 2021, pp. 1037-1048, SWC 2021: ISES Solar World Congress ; Conference date: 25-10-2021 Through 29-10-2021.

[36] CIEMAT, CIEMAT plataforma solar de almería - europe's biggest test center for concentrating solar power (CSP), 2024, https://www.dlr.de/sf/en/desktopdefault.aspx/. (Accessed: 22 January 2024).

[37] BSRN, World radiation monitoring center (WRMC), the central archive of the baseline surface radiation network (BSRN), 2024, https://bsrn.awi.de. (Accessed: 22 January 2024).

[38] Z. Qu, A. Oumbe, P. Blanc, B. Espinar, G. Gesell, B. Gschwind, L. Klüser, M. Lefèvre, L. Saboret, M. Schroedter-Homscheidt, L. Wald, Fast radiative transfer parameterisation for assessing the surface solar irradiance: The heliosat24 method, Meteorol. Z. 26 (1) (2017) 33-57.

[39] A. Hammer, D. Heinemann, C. Hoyer-Klick, R. Kuhlemann, E. Lorenz, R. Müller, H.G. Beyer, Solar energy assessment using remote sensing technologies, Remote Sens. Environ. 86 (2003) 423-432.

[40] P. Lauret, R. Alonso-Suárez, J. Le Gal La Salle, M. David, Solar forecasts based on the clear sky index or the clearness index: Which is better? Solar 2 (4) (2022) 432-444.

[41] Z. Qu, A. Oumbe, P. Blanc, B. Espinar, G. Gesell, B. Gschwind, L. Klüser, M. Lefèvre, L. Saboret, M. Schroedter-Homscheidt, L. Wald, Fast radiative transfer parameterisation for assessing the surface solar irradiance: The heliosat24 method, Meteorol. Z. 26 (1) (2017) 33-57.

[42] A.R. Jensen, K.S. Anderson, W.F. Holmgren, M.A. Mikofski, C.W. Hansen, L.J. Boeman, R. Loonen, Pvlib iotools, open-source Python functions for seamless access to solar irradiance data, Sol. Energy 266 (2023) 112092.

[43] C.W. Chow, S. Belongie, J. Kleissl, Cloud motion and stability estimation for intra-hour solar forecasting, Sol. Energy 115 (2015) 645-655.

[44] C. Liu, Beyond Pixels: Exploring New Representations and Applications for Motion Analysis (Ph.D. thesis), Massachusetts Institute of Technology, Cambridge, MA, USA, 2009.

[45] M. Lefèvre, A. Oumbe, P. Blanc, B. Espinar, B. Gschwind, Z. Qu, L. Wald, M. Schroedter-Homscheidt, C. Hoyer-Klick, A. Arola, A. Benedetti, J.W. Kaiser, J.-J. Morcrette, McClear: a new model estimating downwelling solar radiation at ground level in clear-sky conditions, Atmos. Meas. Tech. 6 (9) (2013) 2403-2418.

[46] P. Bacher, H. Madsen, H.A. Nielsen, Online short-term solar power forecasting, Sol. Energy 83 (10) (2009) 1772-1783.

[47] R.S. Tsay, Analysis of financial time series, first ed., in: Wiley Series in Probability and Statistics, Wiley, 2010, URL https://onlinelibrary.wiley.com/doi/book/10.1002/9780470644560.

[48] R.F. Engle, Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation, Econometrica 50 (4) (1982) 987-1007, Publisher: [Wiley, Econometric Society].

[49] T. Bollerslev, Generalized autoregressive conditional heteroskedasticity, J. Econometrics 31 (3) (1986-04) 307-327.

[50] J.W. Taylor, Volatility forecasting with smooth transition exponential smoothing, Int. J. Forecast. 20 (2) (2004) 273-286.

[51] D. Dumortier, Modelling global and diffuse horizontal irradiances under cloudless skies with different turbidities, Final Report Vol. 2, Daylight II, JOU2-CT92-0144, 1995.

[52] B. Bourges (Ed.), Climatic data handbook for europe, Kluwer Acad. Publ, Europäische Kommission, Dordrecht, 1992, Available online at http://bookshop.europa.eu/en/-pbEUNA13537/.

[53] J. Kühnert, E. Lorenz, D. Heinemann, Satellite-based irradiance and power forecasting for the german energy market, Solar Energy Forecast. Resour. Assessment (2013) 267-295.

[54] L. Delle Monache, F. Eckel, D. Rife, B. Nagarajan, K. Searight, Probabilistic weather prediction with an analog ensemble, Mon. Weather Rev. 141 (2013) 3498-3516.

[55] C. Junk, L. Delle Monache, S. Alessandrini, G. Cervone, L. Bremen, Predictor-weighting strategies for probabilistic wind power forecasting with an analog ensemble, Meteorol. Z. 24 (2015) 361-379.

[56] S. Alessandrini, L. Delle Monache, S. Sperati, G. Cervone, An analog ensemble for short-term probabilistic solar power forecast, Appl. Energy 157 (2015) 95-110.

[57] R. Koenker, Quantile regression, in: Econometric Society Monographs, vol. 38, Cambridge University Press, Cambridge, 2005.

[58] N. Meinshausen, G. Ridgeway, Quantile regression forests, J. Mach. Learn. Res. 7 (6) (2006) 983-999.

[59] L. Visser, T. AlSkaif, J. Hu, A. Louwen, W. van Sark, On the value of expert knowledge in estimation and forecasting of solar photovoltaic power generation, Sol. Energy 251 (2023) 86-105.

[60] S. Raschka, V. Mirjalili, Python Machine Learning: Machine Learning and Deep Learning with Python, scikit-learn, and TensorFlow 2, third ed., Packt Publishing Ltd, Birmingham, 2019.

[61] R. Koenker, G. Bassett, Regression Quantiles, Econometrica 46 (1) (1978) 33-50.

[62] P. Pinson, P. McSharry, H. Madsen, Reliability diagrams for non-parametric density forecasts of continuous variables: Accounting for serial correlation, Q. J. R. Meteorol. Soc. 136 (646) (2010) 77-90.

[63] J. Bröcker, L.A. Smith, Increasing the reliability of reliability diagrams, Weather Forecast. 22 (3) (2007) 651-661.

[64] H. Hersbach, Decomposition of the continuous ranked probability score for ensemble prediction systems, Weather Forecast. 15 (5) (2000) 559-570.

[65] D.S. Wilks, Statistical methods in the atmospheric sciences, second ed., [Nachdr.], in: International geophysics series, (91) Elsevier [u.a.], Amsterdam, 2009.

[66] T. Gneiting, A.E. Raftery, Strictly proper scoring rules, prediction, and estimation, J. Amer. Statist. Assoc. 102 (477) (2007) 359-378.

[67] Z. Ben Boualègue, Assessment and added value estimation of an ensemble approach with a focus on global radiation forecasts, MAUSAN (2015) 541-550.

[68] T. Gneiting, S. Lerch, B. Schulz, Probabilistic solar forecasting: Benchmarks, post-processing, verification, Sol. Energy 252 (2023) 72-80.

[69] K. Doubleday, V. Van Scyoc Hernandez, B. Hodge, Benchmark probabilistic solar forecasts: Characteristics and recommendations, Sol. Energy 206 (2020) 52-67.