Assessment of surface solar irradiance derived from real-time modelling techniques and verification with ground-based measurements

This study focuses on the assessment of surface solar radiation (SSR) based on operational neural network (NN) and multi-regression function (MRF) modelling techniques that produce instantaneous (in less than 1 min) outputs. Using real-time cloud and aerosol optical properties inputs from the Spinning Enhanced Visible and Infrared Imager (SEVIRI) on board the Meteosat Second Generation (MSG) satellite and the Copernicus Atmosphere Monitoring Service (CAMS), respectively, these models are capable of calculating SSR in high resolution (1 nm, 0.05^∘, 15 min) that can be used for spectrally integrated irradiance maps, databases and various applications related to energy exploitation. The real-time models are validated against ground-based measurements of the Baseline Surface Radiation Network (BSRN) in a temporal range varying from 15 min to monthly means, while a sensitivity analysis of the cloud and aerosol effects on SSR is performed to ensure reliability under different sky and climatological conditions. The simulated outputs, compared to their common training dataset created by the radiative transfer model (RTM) libRadtran, showed median error values in the range −15 to 15 % for the NN that produces spectral irradiances (NNS), 5–6 % underestimation for the integrated NN and close to zero errors for the MRF technique. The verification against BSRN revealed that the real-time calculation uncertainty ranges from −100 to 40 and −20 to 20 W m⁻², for the 15 min and monthly mean global horizontal irradiance (GHI) averages, respectively, while the accuracy of the input parameters, in terms of aerosol and cloud optical thickness (AOD and COT), and their impact on GHI, was of the order of 10 % as compared to the ground-based measurements. The proposed system aims to be utilized through studies and real-time applications which are related to solar energy production planning and use.

1 Introduction

Solar energy exploitation is a cornerstone for sustainable development, through efficient energy planning, towards the goal of gradual independence from fossil fuels. In this direction, the European Union (EU), the Middle East and North Africa (MENA) and numerous neighbouring regions and countries have laid out specific technology roadmaps aiming at the integration of low carbon energy technologies linked with the deployment of photovoltaic (PV) installations in the energy market (IPCC, 2012; NREL, 2016; IRENA, 2016; Jager-Waldau, 2016; REN21, 2017; UN, 2017). In addition, the United Nations (2017) has set as its main sustainable development goal by 2030 to ensure universal access to affordable, reliable, and modern energy services. The International Energy Agency (2007) has estimated that the global primary energy demand will increase by 40–50 % from 2003 to 2030. Since energy production, transportation and consumption put considerable pressure on the environment, there is serious concern regarding the sustainability of energy consumption.

Earth observation (EO)-based systems and relevant services already play an important role in the solar energy industry, as well as in human-health-related emerging technologies, but there is still significant potential in increasing their efficiency and exploitation (Schroedter-Homscheidt et al., 2006; Wald et al., 2011; Lefevre et al., 2014). EO from space is already triggering services and applications that can deliver benefits throughout all the phases of energy production and supply. Their contribution ranges from identifying reservoirs and locations with solar energy potential to controlling and monitoring of the distribution networks across Europe, Africa and the Middle East, while providing support to energy policy formulation and enforcement (EU, 2011; IEA, 2010).

The need for improved EO-based surface solar irradiance assessment is increasing as more solar farms are included in national electricity grids worldwide (EC, 2013). Solar-energy-related installations have been increasing their share on the total energy demand as defined by the distribution and transmission system operators (DSOs and TSOs, respectively). As a result, accurate, real-time and short-term forecasting estimations of surface solar radiation (SSR) and, more specifically, global horizontal irradiance (GHI) related to the operation principles of PV installations are vital. The real-time GHI estimations are required at local and regional scales, as well as high temporal frequency (every 5–15 min), in order to be used for near-real-time decisions, linked with the PV-related contribution to the electricity grid.

Since the launch of EO satellites, such as Meteosat Second Generation (MSG), and Sentinel satellite series, real-time image processing techniques have been developed (Suárez and Nesmachnow, 2012). The main advantage of these techniques is the possibility to monitor numerous meteorological variables in almost real time (Derrien et al., 2005; MeteoFrance, 2013). A comprehensive intercomparison of radiation products, codes, algorithms, models and independent databases has been performed by many researchers (Oreopoulos and Mlawer, 2010; Oreopoulos et al., 2012; Ellingson et al., 1991; Ineichen, 2006; Beyer et al., 2009; Cahalan, et al., 2005). Solid steps in estimating the surface GHI were taken by Deneke and Feijt (2008), Schulz et al. (2009), Mueller et al. (2009), Huang et al. (2011) and Qu et al. (2017), who developed GHI retrieval methodologies based on the use of discrete pre-calculated look-up tables (LUTs), while Dorvlo et al. (2002), Zarzalejo et al. (2005), Lopez et al. (2001) and Takenaka et al. (2011) developed solutions based on neural network (NN) models. The validation of most of the above mentioned methodologies was performed against radiative transfer model (RTM) simulations and ground-based measurements, from various networks around the globe. However, from the validation results it was shown that accuracy was inversely proportional to calculation speed under all-sky and terrain conditions. The magnitude of the GHI uncertainty due to the effect of aerosols and clouds is significant and has motivated numerous related studies (Federico et al., 2017; Kosmopoulos et al., 2015; Lara-Fanego et al., 2012; Tegen et al., 1996; Lindfors et al., 2013). Under high aerosol loads the SSR can be reduced by 20–50 % (Eck et al., 1998; Gleeson et al., 2016; Kosmopoulos et al., 2017), while under cloudy conditions the impact was up to 60–90 % for overcast conditions and cloud coverage of 8 octas (Aebi, et al., 2017; Kosmopoulos et al., 2015; Zygmuntowska et al., 2012), highlighting the significant effect of these atmospheric parameters (clouds and aerosols) on the GHI calculations and in the performance of PV installations and energy production.

In the present study, we report on (i) the assessment of the surface solar irradiance calculated in real time, which is defined as the product with a time delay of 1 min or less from an actual atmospheric situation, by developing and using two NN-based techniques and a multi-regression-function-based technique and (ii) the validation of these techniques against ground-based measurements from the Baseline Surface Radiation Network (BSRN). Section 2 presents data, methods and techniques used. Section 3 describes the validation results including a sensitivity analysis of related atmospheric parameters and in Sect. 4 we present our conclusions on the proposed techniques.

Table 1 Coordinates (degrees) and height (metres above sea level) of the BSRN stations used for the validation.

Figure 1 Flowchart of the modelling technique scheme. The initial pixel classification followed by the clear- or cloudy-sky inputs to the real-time solver result the spectral (NNS) and integrated (MRF and NN) SSR-related outputs.

2 Data and methodology

2.1 Data

2.1.1 Ground-based measurements

The verification of the applied SSR real-time modelling techniques was performed against ground-based measurements from nine stations (Table 1) of the Baseline Surface Radiation Network (BSRN; Hegner et al., 1998) equipped with Kipp and Zonen pyranometers (GHI measurements) and a precision filter radiometer (PFR) at Izaña, Spain. BSRN consists of high-quality ground-based measurements of SSR and for the purposes of the comparison we used the dataset from July 2014 to June 2015. Table 1 presents the location and description of the nine BSRN stations used for the validation of the SSR estimations calculated with the modelling techniques. The temporal resolution of the ground-based measurements is 1 min, so in order to match the 15 min resolution of the MSG cloud data (and hence the SSR outputs) we used 15 min averages of all the BSRN and PFR measurements used. The selected BSRN stations represent a variety of different climates, altitudes and aerosol sources in the field of view of MSG and thus provide an opportunity to study the models' performance under various atmospheric conditions.

2.1.2 Real-time cloud observations

The most important inputs to our real-time modelling techniques were the satellite cloud data products from the Spinning Enhanced Visible and Infrared Imager (SEVIRI) on board the Meteosat Second Generation (MSG) satellite. We obtained the cloud type (CT), the cloud phase (CP) and the cloud optical thickness (COT) products so as to efficiently quantify the effect of clouds on SSR. COT depends on the moisture density as well as the vertical thickness of the cloud. The cloud reflectance at channel at 0.6 µm in the visible part of the electromagnetic spectrum is directly related to COT (Roebeling et al., 2006). The MSG geostationary satellite, because of its orbit height (36 000 km above the Equator), allows the continuous monitoring over Europe, Africa and parts of South America at high temporal and spatial resolution (15 min and 0.05^∘, respectively). The operational MSG-SEVIRI data were acquired by the EUMETCast station operated by the Institute for Astronomy, Astrophysics, Space Applications and Remote Sensing of the National Observatory of Athens. The cloud properties are extracted operationally and in real time using the Satellite Application Facility for Nowcasting Weather Conditions software (SAFNWC) installed in-house. CT and CP are standard output products of the SAFNWC computational procedure, while COT is a tailor-made product and as a result its extraction required an additional intervention in the process chain. The cloud product identification is described in Derrien and Gléau (2005) and the MétéoFrance (2013) technical report. In the current implementation, cloud products are provided operationally for the entire Earth disc view area of MSG. We extracted products at specific pixels corresponding to locations of the BSRN stations, which were used as inputs to the SSR modelling techniques.

2.1.3 Aerosol forecasts

For the real-time assessment of the SSR we additionally incorporated the aerosol 1-day forecast data from the Copernicus Atmospheric Monitoring Service (CAMS) as the basic input parameter. These forecasts are based on the Monitoring Atmospheric Composition and Climate (MACC) reanalysis tools, and include validated modelling of aerosol and satellite data assimilation (Eskes et al., 2015). They are able to provide operationally accurate data of aerosol optical depth (AOD) at 550 nm, at 1 h time steps and 0.4^∘ spatial resolution. The estimation of the aerosol sources is extracted from the Emission Database for Global Atmospheric Research (EDGAR) and the Speciated Particulate Emission Wizard (SPEW), while the reliability of the product is supported by continuous assimilation into the model of the MODIS AOD data, applying a bias correction from multiple data sources (Dee and Uppala, 2009). For the purposes of our SSR estimations, the CAMS AOD forecasts with the MSG COT data described above constitute the most important input parameters, together with solar elevation, for the SSR retrieval modelling tools.

2.2 Methodology

In this section we present the SSR real-time modelling techniques, the methodology used for developing operational products and the validation statistics against ground-based measurements. The techniques are the multi-regression function (MRF), the neural network that produces spectral irradiances (NNS) and which is presented in detail in Taylor et al. (2015), and a variant version of the NN that produces integrated irradiances. All three techniques have been optimized based on LUTs that are described in the Sect. 2.2.1 and produce instantaneous (with less than 1 min delay from the time that the MSG image is produced) SSR. The number of outputs depends on the region under study and can be of the order of 10⁶ simulations simultaneously. In this study we used the CAMS AOD and the MSG COT as operational inputs, in conjunction with the solar elevation angle, as they are the major attenuators of the GHI. Since the comparison of real-time modelling techniques with ground-based measurements are performed from southern Africa to northern Europe, the verification will be focused on GHI. Utilization of the direct normal irradiance (DNI) by concentrated solar power (CSP) plants is limited at places with high amounts of DNI (Green et al., 2015), and hence CSPs are de facto outside of energy planning for the majority of the countries represented by the nine BSRN stations and the MSG view. Large-scale, high temporal and spatial resolution EO-based assessment of the SSR seems to be an emerging market prospect (ITA, 2016). The potential application fields of the methodology proposed in this study include the production planning support on large-scale solar farm projects and the efficient control of the electricity balancing and distribution (in support of the TSOs and DSOs) by incorporating the produced energy of the solar farms into the electricity grid. At the same time, SSR in different spectral regions highlight spectrally weighted outputs like the UV index (linked with skin cancer, eye cataracts, DNA damage etc.), vitamin D efficiency (related to pregnancy) and a number of agricultural and oceanographical related processes (plant photosynthesis, crop production, phytoplankton growth etc.). As a result, the developed real-time modelling techniques are able to assist public authorities in energy-planning policies; support the work of various scientific communities dealing with health protection, energy production and consumption and solar energy exploitation; and finally enable the solar industry to better plan clean energies, its transmission and distribution, which in turn will boost the relative contribution to national portfolios. Figure 1 illustrates the procedural flows of the three developed real-time modelling techniques for operational use. Starting with the MSG cloud flags (0 = clear sky and 1, 2, 3 = cloudy sky in terms of water, ice and mixed clouds, respectively), we identify the clear-sky and cloudy-sky pixels. For the cloudy pixels we incorporate the optical properties (COT) and types of clouds (CT), while for clear-sky pixels we take into account the aerosol effect (AOD) and the total ozone column (TOC), which was derived using Ozone Monitoring Instrument (OMI) retrievals (Wandji-Nyamsi et al., 2015). Then, for all-sky conditions we generate the input files to the real-time techniques and, depending on their special characteristics, we produce spectral or spectrally weighted products (see following subsections) at high spectral, spatial and temporal resolution (1 nm, 0.05^∘, 15 min). The actual outputs can be SSR time series, local and regional maps or Earth disc view maps (Fig. 2).

Figure 2 An example of the output maps based on the real-time SSR techniques. Here is the GHI for 15 April 2015 at 12:00 UTC together with the BSRN station locations.

The performance of the real-time techniques was evaluated by comparing the GHI outputs with (i) the initial RTM simulation LUTs and (ii) the BSRN ground-based measurements and with respect to the aerosol and cloud effects. The evaluation was based on the bias and mean bias error (MBE), the root mean square error (RMSE) and their relative components (rMBE and rRMSE, respectively):

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathrm{MBE}=\overline{\mathit{\epsilon }}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\mathit{\epsilon }}_i,\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N{\mathit{\epsilon }}_i^{\mathrm{2}}}.\end{array}$$

The residuals (estimation errors), ${\mathit{\epsilon }}_i=x_e-x_{\mathrm{0}}$, are calculated as the difference between the estimated values by the real-time techniques (x_e) and the measured values (x₀) by BSRN, where N is the total number of data. MBE measures the overall bias and detects the model's overestimation (MBE > 0) or underestimation (MBE < 0). RMSE quantifies the spread in the distribution of errors. Concerning the rMBE and rRMSE error measures, the normalization is done with respect to the mean ground measurement irradiance in the considered station and period. In addition, for the various tests performed in this study we calculated the slope, the correlation coefficient (r), the coefficient of determination (r²), the percentage difference (%), the mean absolute difference and the standard deviation.

Table 2 RTM-simulated GHI at 15 min time intervals as compared to the BSRN ground-based measurements in terms of correlation coefficient (r) and slope.

Figure 3 Percentage difference (%) of the real-time modelling techniques as compared to the RTM simulations for all ground stations. The box charts highlight the more precise estimation approach of the MRF technique as compared to the NN-based techniques.

Figure 4 The altitude correction of GHI for various SZAs as a function of the SSR ratio (SSR at height h compared to SSR at sea level).

2.2.1 Radiative transfer model

All modelling techniques presented in this paper for the real-time assessment of the SSR are based on LUTs, calculated with the radiative transfer model (RTM) libRadtran (Mayer and Kylling, 2005; Emde et al., 2016). These LUTs are described in detail in Taylor et al. (2015) and consist of more than 2.5 million RTM simulations with atmospheric inputs and 1 nm spectral resolution GHI outputs. The interoperable exchange of similar GHI databases is studied by Ménard et al. (2015), highlighting the usefulness and necessity of such LUT-based approaches (Lefevre et al., 2014). Under clear-sky conditions the simulated by libRadtran input parameters were the solar zenith angle (SZA), the AOD, the Ångström exponent (AE), the single-scattering albedo (SSA), TOC and the columnar water vapour (WV), while under cloudy conditions except from SZA and TOC, we also used the optical thicknesses of water and ice clouds (WCOT and ICOT, respectively) as inputs. The AOD is not used for cloudy conditions when COT > 1, as the effects of aerosols are much weaker compared to thick clouds. For the model versus BSRN station comparison, in order to take into account the station altitude, an altitude correction on the solar energy output of the different model simulations has been applied based on RTM (Libradtran) calculations. The outputs are high-resolution spectral irradiances (1 nm) covering the wavelength region between 285 and 2700 nm. In brief, we used the SDISORT radiative transfer solver (Dahlback and Stamnes, 1991) with pseudospherical approximation to produce valid outputs from 0 to 90^∘ SZA; the simulations were calculated using a band parameterization method based on the correlated-K approximation (Kato et al., 1999), while the aerosol and cloud determination was performed based on the default aerosol model described by Shettle (1989) and typical cases for the height of water and ice clouds, the effective radius (R_eff) and the liquid water path (Hess et al., 1998). All the technical and structural information about the RTM simulations, the input parameters and the construction of the LUTs is presented in Taylor et al. (2015). Table 2 presents the slope and the correlation coefficient between the RTM simulations of GHI and the BSRN ground-based measurements for the whole datasets and period. The overall accuracy in terms of slope ranges from 0.866 (CAB) to almost 1 (0.999 at TOR), while the r values range between 0.93 and 0.97.

Figure 5 The mean GHI in W m⁻² of the real-time modelling techniques as compared to the RTM simulations for all ground stations (a), and the mean GHI of all models as compared to the BSRN measurements (b).

Figure 6 Mean GHI differences in W m⁻² derived by MRF as compared to the BSRN stations for each time horizon. The boxes represent the 25th and 75th percentiles, while the in-box lines represent the median of the difference of each station. The upper and lower whiskers represent the minimum and maximum error values.

2.2.2 Multi-regression function

The multi-regression function (MRF) technique was developed as an analytical methodology using the RTM outputs, with the aim to provide results as close as possible to the initial (training set) RTM outputs. The advantage in the use of these functions is that they can be executed very rapidly and can be used for real-time SSR determination. In order to achieve that, analytical functions for the SSR should be constructed. In general, SSR is a function of SZA, COT, AOD, AE, SSA, WV and TOC (Appendix A presents the complete list of nomenclature and abbreviations). For the AE and SSA we used monthly climatological values in order to bridge the gap between the operational input availability and the SSR accuracy. However, a preliminary investigation has been performed for the sensitivity of GHI to WV column and TOC. We compared integrated spectral GHI over the entire spectrum for different TOC values and we found a mean difference of only 0.5 % for TOC ranging between 300 and 400 DU. For WV columns ranging between 0.5 and 2 cm we found a mean difference of 3.2 %, although for SZA < 15^∘ this difference was higher, up to 5 %. Therefore, we chose to neglect these variables in the first place and use TOC = 350 DU and WV = 0.5 cm for further calculations, considering the differences mentioned above as a scale of error introduced by this approach.

Table 3 Values of parameters used for the polynomial function Eq. (3) of the MRF technique for GHI calculations under clear-sky and cloudy-sky conditions.

Then, we constructed different polynomial functions according to Gasca and Sauer (2000) for cloudy and clear-sky conditions, to be applied into the scheme presented in Fig. 1. For cloudy cases the irradiance is expressed as f_cloud(SZA, COT) and for clear-sky cases as f_clear(SZA, AOD). We tested different orders of two-variable polynomials to conclude on the best regression (multi-regression analysis) and we found that the estimates closest to the RTM results were achieved using fifth- and fourth-degree polynomials (Sauer and Xu, 1995), as follows:

$$\begin{array}{ll}{\displaystyle }f\left(x,y\right)& {\displaystyle }=p_{\mathrm{00}}+p_{\mathrm{10}}x+p_{\mathrm{01}}y+p_{\mathrm{20}}{x}^{\mathrm{2}}+p_{\mathrm{11}}xy+p_{\mathrm{02}}{y}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }+p_{\mathrm{30}}{x}^{\mathrm{3}}+p_{\mathrm{21}}{x}^{\mathrm{2}}y+p_{\mathrm{12}}x{y}^{\mathrm{2}}+p_{\mathrm{03}}{y}^{\mathrm{3}}+p_{\mathrm{40}}{x}^{\mathrm{4}}\\ {\displaystyle }& {\displaystyle }+p_{\mathrm{31}}{x}^{\mathrm{3}}y+p_{\mathrm{22}}{x}^{\mathrm{2}}{y}^{\mathrm{2}}+p_{\mathrm{13}}x{y}^{\mathrm{3}}+p_{\mathrm{04}}{y}^{\mathrm{4}}+p_{\mathrm{41}}{x}^{\mathrm{4}}y\\ \text{(3)}& {\displaystyle }& {\displaystyle }+p_{\mathrm{32}}{x}^{\mathrm{3}}{y}^{\mathrm{2}}+p_{\mathrm{23}}{x}^{\mathrm{2}}{y}^{\mathrm{3}}+p_{\mathrm{14}}x{y}^{\mathrm{4}}+p_{\mathrm{05}}{y}^{\mathrm{5}}\end{array}$$

where x is SZA and y is AOD and COT accordingly (clear- or cloudy-sky pixels). Table 3 presents the analytical values of p_xx for the purposes of this study (GHI) under clear- and cloudy-sky conditions. By this approach RTM simulations of SSR are derived in computational times that can be applied in any real-time application.

2.2.3 Neural network

As presented in Taylor et al. (2015), the LUT approach, despite its large size, still provides estimates at discrete input values. The interpolation techniques to correct the input-output parameter intervals are computationally more costly than a continuous function-approximating model, or a NN model, which is more preferable for producing real-time outputs (Hornik et al., 1989). Indicatively, using a test set of 1000 RTM simulations from the developed LUT, we applied an interpolating function to adjacent/nearest value and found that each interpolation calculation required a time in excess (in total ≈ 21 h) of each single run of RTM used to generate the LUT in the first place (≈ 12 h for 1000 RTM simulation outputs with spectral resolution of 1 nm in the range 285–2700 nm), while for the same test set the NN needed almost 0.144 s to generate the 1000 output spectra. Takenaka et al. (2011) have pointed out that the inclusion of many parameters (we incorporated six for the clear-sky and four for the cloudy-sky simulations) and small step sizes (we produced more than 2.5 million RTM simulations in total) can dramatically increase the LUT volume, while Sauer and Xu (1995) and Gasca and Sauer (2000) noted that the multidimensional nature of the dataset requires interpolation/extrapolation procedures that impact strongly on calculation speed. Hence, based on the developed LUT, we trained two sets of NNs, each one consisting of a clear-sky- and a cloudy-sky-specific NN. For multivariate input–output data, feed-forward NNs with a minimum of one layer of “hidden” neurons whose activation functions are nonlinear hyperbolic tangent functions or other general nonlinear sigmoidal functions have been shown in the literature to be a universal function approximator (Cybenko, 1989; Hornik et al., 1989). The input–output vectors used in this study were connected via two network layers – the first containing hidden neurons with tanh activation functions and the second containing output neurons with linear activation functions. The exact mathematical equation relating the NN outputs to the NN inputs for this type of NN is given in the following matrix equation described analytically in Taylor et al. (2014):

$$\begin{array}{}\text{(4)}& {\displaystyle }Y=f^{\mathrm{2}}\left(L{W}^{\mathrm{2},\mathrm{1}}{f}^{\mathrm{1}}\left({\mathbf{IW}}^{\mathrm{1},\mathrm{1}}\mathit{X}+b^{\mathrm{1}}\right)+b^{\mathrm{2}}\right).\end{array}$$

The multiplication of the matrix IW^1,1 and the vector X is a dot product equivalent to the summation of all input connections to each neuron in the hidden layer. This equation is the continuous and nonlinear functional approximation that relates the output vector to the input vector. This NN approach, its training procedure and all the technical details are described analytically in Taylor et al. (2015).

In the first NN set, we produced instantaneous SSR spectra of the order of 1 million in less than 1 min, using as operational inputs the CAMS AOD 1-day forecasts, the MSG COT and real-time calculations of SZA. The output resolution is high in terms of spectral (1 nm), spatial (0.05^∘) and temporal (15 min) components (Taylor et al., 2015), and, operationally speaking, this spectrally based NN (NNS) can incorporate additional inputs as described in Sect. 2.2.1. Similar studies on the temporal variability of SSR by means of spectral representations and the wavelengths absorption parameterization applied to satellite channels and spectral bands were performed by Gasteiger et al. (2014) and Belgulescu et al. (2016). For the purposes of this study we used monthly climatological values for the rest of the input parameters, more specifically TOC from OMI (2007–2016), WV from the Medium Resolution Imaging Spectrometer (MERIS) on board ESA's Environmental Satellite (ENVISAT), and AE and SSA from the AeroCom database (Kinne et al., 2006). The second NN set was trained using integrated SSR over the whole wavelength range using the LUT's spectral data. The SSR results of this technique (called hereafter NN) are more accurate in terms of GHI, DNI and diffuse horizontal irradiance (DHI), as will be discussed in the following section. On the other hand, spectrally weighted products like the UV index, photosynthetically active radiation (PAR) or the vitamin D effective dose (VDED) cannot be produced with this approach, as only the NNS is able to produce the spectral irradiance needed for such applications.

Since the proposed modelling techniques (MRF, NN and NNS) operate in real time, the potential applicability for short-term forecasting purposes for the next few hours is feasible. In this direction, the CAMS AOD is already an operational forecast input (Benedetti et al., 2009), with accurate predictions every 1 h even under high aerosol load conditions (Kosmopoulos et al., 2017). On the other hand, the MSG COT short-term forecasting requires the employment of a cloud motion vector analysis (e.g. Hammer et al., 1999) in high spatial and temporal resolution (5 km × 5 km and 15 min, which is the MSG/SEVIRI resolution) in order to predict the impact of clouds on SSR for the next 2–3 h, while under cloudless conditions the SZA and AOD are the main solar irradiance attenuators and hence are available as input information to the models.

3 Results

3.1 Performance of real-time techniques

3.1.1 Comparison with RTM

This section initially summarizes the performance of all the real-time modelling techniques against the RTM simulations for all BSRN stations. Figure 3 presents the percentage difference between the RTM simulations and the MRF, NN and NNS techniques. All data presented here are GHI model outputs with a 15 min temporal resolution. The box plots represent the interquartile range between the 25th and 75th percentiles with the in-box line to show the median and the upper and lower whiskers to represent the maximum and minimum error values that are within 1.5× the interquartile range of the box edges. The largest differences for all techniques occur for LER and TOR stations followed by CAB, indicating higher introduced uncertainties over highest latitudes, as observed on the MSG Earth view edges. However, differences for MRF are much smaller for these three stations. For all stations MRF shows differences around zero, showing quite efficient representation of the LUT-based RTM simulations. For the altitude correction (described in Sect. 2.2.1) we included corrections of 4.2 % per kilometre at 20^∘ SZA up to 12 % per kilometre at 80^∘ (Fig. 4), based on libRadtran model sensitivity analysis.

Figure 7 Seasonal relative RMSE values of the GHI estimations produced by the real-time techniques as compared to the BSRN measurements.

The NN and NNS approaches showed a systematic underestimation for the NN of ∼ 8 %, while NNS had comparatively the worst performance, with differences in the range of −15 to 80 % (LER station – for the interquartile ranges). The median differences range from −15 to 15 % for NNS and ∼ 5 to 6 % for NN, while they are less than 1 % for MRF. It is clear that spectral output methods (NNS) provide more detailed information (e.g. for specialized studies on spectral impacts on the yield of different PV technologies; Dirnberger et al., 2015; Ishii et al., 2013), but they are more uncertain than the NN and MRF that produce integrated SSR.

3.1.2 Verification with BSRN

The model accuracy was verified against nine BSRN stations. We calculated the regression of the mean GHI between the ground measurements and the model outputs, shown in Fig. 5. We also show the intra-model regression compared to the initial RTM simulations (Fig. 5a) in order to assess the NN and NNS included interpolations of the LUT outputs and the MRF performance. We found that the MRF technique presents identical values to the RTM, for all ground stations and under all climatological conditions. The NN and NNS show quite good agreement too in terms of absolute values, as under all conditions mean GHIs are less than 5 % different from the BSRN measurements. In Fig. 5b we confirmed the similarity of MRF with RTM and in some cases with the NN models, indicating the overall efficiency of all interpolation and multi-function techniques used. A slightly better performance was observed for higher mean GHIs, proving the usefulness under high solar energy potential conditions.

Table 4 GHI evaluation results as a function of season and real-time techniques for all stations. The model MBE and RMSE statistical scores are shown in absolute units (W m⁻²) and as relative magnitude (percentages in brackets).

Figure 8 Scatterplots of real-time (MRF and NN) and RTM-simulated GHI in W m⁻² as compared to the BSRN measurements for all stations under clear-sky (a) and all-sky (clear-sky and cloudy) conditions (b).

Figure 9 Mean percentage difference (blue) and standard deviation (red) of the 15 min GHI produced from the MRF technique as compared to ground-based measurements from all stations as a function of the COT.

Figure 6 shows the accuracy of MRF, being the most reliable technique as presented in Fig. 3, with respect to the ground-based measurements for various temporal integrations, starting from the (actual derived) 15 min to hourly, daily and monthly averages. The uncertainty range of the MRF simulations given as mean interquartile GHI differences is highest (from −100 to 40 W m⁻², depending on the station) for the 15 min resolution. It is reduced for hourly and daily averages (−70 to 40 and −40 to 30 W m⁻², respectively) and is minimized for the monthly averages (−20 to 20 W m⁻²). In particular, IZA and TAM showed the highest differences for all temporal retrievals, while LER and TOR presented minimum differences down to ±20 W m⁻² for the interquartile range of the 15 min averages. The median values are within 10 W m⁻² for the 15 min and hourly resolutions, while the corresponding minimum and maximum error values (represented in Fig. 6 as the upper and lower whiskers) extend from −200 to 100 W m⁻² for the aforementioned resolutions and are reduced to ±60 and ±40 W m⁻² for the daily and monthly averages, respectively. These results are comparable with similar model verification approaches and studies (Riihela et al., 2015; Muller et al., 2015; Thomas et al., 2016; Eissa et al., 2015a, b). Indicatively, Muller et al. (2015) and Riihela et al. (2015) discussed the CM-SAF SARAH (Solar surfAce RAdiation Heliosat) data record, which consists of post-processed data. They calculated a mean monthly error for GHI of 5.5 W m⁻² and a mean daily error of 12.1 W m⁻², with additional uncertainties in terms of spatial representativeness and measurement quality of about ±12 W m⁻², while they did not provide relevant information about the hourly or even higher time resolution. The overall accuracy of all models was also evaluated with respect to seasonality. In Fig. 7 we present the seasonal rRMSE values of the GHI estimations produced by the MRF, NN and NNS models as compared to the BSRN 15 min interval measurements. The rRMSE for MRF ranges from 5 to 48 % for GOB and TOR stations, respectively; for NN the range increases to 6–60 %; and for NNS the corresponding range is 7–87 %. We need to note that the aforementioned large differences correspond to significantly low absolute GHI values, indicating the impact under cloudy conditions mainly in the winter season, at stations with high mean cloudiness (LER, TOR, CAM and CAB). In the summer, results are better for all stations and for all models (5–29 %), while GOB, IZA and TAM stations showed the lowest rRMSE values (5–12 % in summer and 12–17 % in winter), linked with their lower cloudiness. Eissa et al. (2015a, b) validated the HelioClim-3 database and the McClear model in Egypt and in the United Arab Emirates, and they found RMSE of 68.4–151.7 and 22–47 W m⁻², respectively (we found a range of 58.2–70.8 W m⁻²). Thomas et al. (2016) validated the latest version of HelioClim-3 (v5) against BSRN and found rRMSE of 14.1–37.2 % for the 15 min averages, which is directly comparable to our 15 min results (12–35.7 %). In particular, for the LER, TOR, CAB, CAM, CAR and TAM stations they found rRMSE of 37.2, 33, 29.4, 25.9, 16.3 and 15.8 %, while from our MRF performance evaluation results we observe 35.7, 35.6, 29.9, 30.3, 20.2 and 12.2 % for the same stations. This indicates that the use of the suggested real-time modelling techniques enables the production of instantaneous, high-resolution and quite accurate (as compared to the post-processed databases) GHI outputs that can be used for solar-energy-related applications and studies. A detailed presentation of results for all metrics and stations can be found in Table 4.

3.2 Sensitivity analysis

3.2.1 Cloud effect

The cloud effect via the radiative transfer of solar radiation in the atmosphere represents the greatest source of uncertainty in the simulation of SSR, while several models do not have the capability to deal with clouds coexisting with a radiatively active atmosphere (Cahalan et al., 2005). Small changes in cloudiness and its optical properties can impact on GHI. The magnitude of the cloud effects on the model to BSRN comparison can be seen in Fig. 8. Under clear-sky conditions (Fig. 8a), the regression of the 15 min modelled GHI values, in terms of coefficient of determination (r²), shows very good agreement when compared with the BSRN measurements for both MRF (0.952) and NN-based (0.924) techniques. We plotted the RTM simulations as well in order to depict the corresponding regression (0.958). The distinct scatter shown under all-sky conditions (Fig. 8b) with the cloud cases linked with an underestimation of the modelled GHI in comparison to the BSRN values, while the corresponding r² decreased to 0.887 and 0.867 for the MRF and NN techniques, respectively (the RTM was almost identical to the MRF, i.e. r²= 0.889). This effect has to do with the MSG COT uncertainties and hence introduces errors into the outputs of the SSR techniques (Derrien and Le Gléau, 2005; Pfeifroth et al., 2016). In addition, comparison principles of (point) station GHI measurements with a 0.05^∘ MSG cloud “picture” are responsible for part of the observed deviations. As an example, for instants in which the MSG 0.05^∘ grid is partly cloudy, the BSRN GHI measurements could fluctuate more than 100 %, depending on whether the sun is visible or whether clouds attenuate the direct component of the solar irradiance. As a result, in the case of a partly covered 0.05^∘ pixel and in the absence of clouds between the BSRN instrument and the sun, BSRN-measured GHI would be much higher than the modelled one. Of course, the opposite situation is feasible as well, consequently causing an overestimation of the modelled GHI (Koren et al., 2007).

Figure 9 illustrates the mean percentage difference and standard deviation of the 15 min GHI produced by the MRF and the measured values by the BSRN stations (only instances with cloudy conditions were used for all stations) as a function of COT. For COT < 2, the MRF technique results in higher GHI values than those actually measured, 1–12 %, while as the COT values increase the MRF underestimates the measurements by up to −60 % for COT around 35. We note that under such high COT values the mean radiation values are much lower than 50 W m⁻². The standard deviation reaches its highest value of 43 W m⁻² for COT 14–16, while its lowest value of 32 W m⁻² is found for COT 2.6.

3.2.2 Aerosol effect

In addition to the clouds, aerosols play an important role in the solar radiation transfer in the atmosphere. Especially in places with high solar energy potential, where cloud-free conditions prevail during the greater part of the year, significant aerosol sources could exist (Gkikas et al., 2012). The aerosols effect is closely related to the aerosol optical properties and mostly AOD, and as a consequence the uncertainty in the model AOD input could result in significant errors in the assessment of SSR (Oumbe et al., 2015; Kosmopoulos et al., 2017). For the purposes of this study we used the Global Atmosphere Watch (GAW-PFR) station of IZA, which is an internationally recognized test bed for aerosol remote sensing instruments (Cuevas et al., 2016), to quantify the AOD difference between the operational input from CAMS and a PFR instrument, under high-altitude conditions (Garcia et al., 2013). In Fig. 10 we present the yearly frequency distribution of the differences between CAMS and PFR values for cloudless sky conditions. The majority of the AOD differences are lower than 0.2, with the maximum frequency encountered at zero AOD differences, indicating the overall good accuracy of CAMS-derived 1-day forecasts of AOD. The mean absolute difference was found equal to 0.1075 ± 0.1038 (1σ). This shows an overestimation of 0.1 for CAMS that could lead to a small MRF GHI underestimation of 2 % compared with BSRN measured GHI. Finally, in Fig. 11 a scatterplot of the CAMS-PFR differences in AOD is shown as a function of absolute differences in GHI derived between the MRF technique and the IZA measurements. The GHI differences are spread around zero independently of the AOD difference, showing the negligible dependence of such small AOD differences on the GHI model calculations.

Figure 10 Frequency histogram of differences between the CAMS and the PFR AOD at the Izaña station together with the mean absolute difference and standard deviation metrics.

Figure 11 Absolute differences in GHI (in W m⁻²) derived by the MRF technique from the ground-based measurements at Izaña (BSRN pyranometer), as a function of differences in AOD from CAMS and PFR.

4 Summary and conclusions

This study proposed state-of-the-art modelling techniques (NNS, NN, MRF) for the real-time estimation of SSR, which have been validated against ground-based BSRN measurements. The determination and understanding of the input parameter effects on radiative transfer revealed that the accuracy of simulations depends on the quality and resolution of the atmospheric inputs to the models (mostly COT and AOD), while increasing the calculation speed and including spectral GHI information decreases the model accuracy.

We firstly described the developed modelling techniques, which are based on large LUTs for clear-sky and cloudy conditions. Verification of these models was performed for the GHI against ground measurements at nine stations, with variable geographical, atmospheric and altitudinal conditions. The comparison showed a dependence on seasonal variability, with summer rRMSE values below 30 % for all models and under all conditions, and revealed largest errors for the NNS technique because of the spectral special characteristics, as well as for LER and TOR stations. The NN presented a slight underestimation of 8 % against its training RTM simulations, while against BSRN stations it achieved MBE and RMSE values less than 30 and 80 W m⁻², respectively, for the annual period, indicating relatively good agreement under various conditions. The technique with the most accurate results, which was almost identical to the RTM simulations, was MRF. Under different temporal scales the mean GHI differences in terms of 25th to 75th interquartiles, compared to the nine stations, were found to range from −100 to 40 W m⁻² for the 15 min intervals, −70 to 40 W m⁻² for the hourly means, and −40 to 30 and −20 to 20 W m⁻² for the daily and monthly averages, and they were almost 10 W m⁻² for the median of difference of each station.

The results presented here show the potential use of such techniques for solar-energy-related applications and electricity grid support services (IRENA, 2015). Comparison of the proposed real-time models with existing databases (e.g. SARAH), which in most cases are post-processed data using past data series, showed similar results. Finally, we tested the impact of cloud and aerosol inputs to the models in order to reveal the AOD forecast accuracy of CAMS, which turned out to be ∼ 0.1 in absolute terms as compared to ground-based sun-photometric measurements in Izaña. The CAMS AOD performance has also been tested under high aerosol loads (Kosmopoulos et al., 2017) in different regions (e.g. the eastern Mediterranean), showing similar results compared with MODIS. However, its accuracy should be checked in the case of application of the methodology to different regions (e.g. Middle East). The MSG COT is related to MRF underestimation of the order of 60 % under highly cloudy conditions (COT > 30) and negligible GHI levels (< 50 W m⁻²). As a result, the presented real-time models based on the synergy of satellite products, RTM and NN or MRF techniques, are a promising tool to be used within the solar-energy-related community. Improvements on satellite-based model inputs from latest and future satellite missions (e.g. Sentinel missions) could be implemented in the future in the existing system in order to improve spatial and temporal resolution and GHI accuracy.

Appendix A: Nomenclature and abbreviations

AE	Ångström exponent
AOD	Aerosol optical depth
BSRN	Baseline Surface Radiation Network
CAMS	Copernicus Atmosphere Monitoring Service
CM SAF	Satellite Application Facility on Climate Monitoring
COT	Cloud optical thickness
CP	Cloud phase
CSP	Concentrated solar power
CT	Cloud type
DHI	Diffuse horizontal irradiance
DNI	Direct normal irradiance
DU	Dobson unit
EDGAR	Emission Database for Global Atmospheric Research
ENVISAT	Environmental Satellite
EO	Earth observation
ESA	European Space Agency
EU	European Union
GAW	Global Atmosphere Watch
GHI	Global horizontal irradiance
ICOT	Ice cloud optical thickness
libRadtran	Library for Radiative transfer
LUT	Look-up table
MACC	Monitoring Atmospheric Composition and Climate
MBE	Mean bias error
MENA	Middle East and North African countries
MERIS	Medium Resolution Imaging Spectrometer
MODIS	Moderate Resolution Imaging Spectroradiometer
MRF	Multi-regression function
MSG	Meteosat Second Generation
NN	Neural network
NNS	Neural network spectral
OMI	Ozone Monitoring Instrument
PAR	Photosynthetically active radiation
PFR	Precision filter radiometer
PV	Photovoltaic
rMBE	Relative mean bias error
RMSE	Root mean square error
rRMSE	Relative root mean square error
RTM	Radiative transfer model
SAFNWC	Satellite Application Facilities for NoWCasting
SARAH	Solar surfAce RAdiation Heliosat
SEVIRI	Spinning Enhanced Visible and InfRared Imager
SPEW	Speciated Particulate Emission Wizard
SSA	Single-scattering albedo
SSR	Surface solar radiation
SZA	Solar zenith angle
TOC	Total ozone column
UV	Ultraviolet
VDED	Vitamin D effective dose
WCOT	Water cloud optical thickness
WV	Water vapour