Analysis of Spatial and Temporal Variability of the PAR / GHI Ratio and PAR Modeling Based on Two Satellite Estimates

: The main objectives of this work are to address the analysis of the spatial and temporal variability of the ratio between photosynthetically active radiation ( PAR ) and global horizontal irradiance ( GHI ), as well as to develop PAR models. The analysis was carried out using data from three stations located in mainland Spain covering three climates: oceanic, standard Mediterranean, and continental Mediterranean. The results of this analysis showed a clear dependence between the PAR / GHI ratio and the location; the oceanic climate showed higher values of PAR / GHI compared with Mediterranean climates. Further, the temporal variability of PAR / GHI was conditioned by the variability of clearness index, so it was also higher in oceanic than in Mediterranean climates. On the other hand, Climate Monitoring Satellite Facility (CM-SAF) and Moderate-Resolution Imaging Spectroradiometer (MODIS) data were used to estimate PAR as a function of GHI over the whole territory. The validation with ground measurements showed better performance of the MODIS-estimates-derived model for the oceanic climate (root-mean-square error (RMSE) around 5%), while the model obtained from CM-SAF ﬁtted better for Mediterranean climates (RMSEs around 2%).


Introduction
Photosynthetically active radiation (PAR) is the part of the solar spectrum included in the 400-700 nm range. PAR is a source of energy responsible for photosynthesis. Depending on the number of photons that plants can absorb, photosynthesis rate is affected, that is, specifically the quanta of energy necessary to create the products of photochemical reactions used in the process of CO 2 assimilation [1].
The radiation intensity changes along the daytime and plants try to keep a balance between converting radiation energy, and protecting the photosynthetic system against photoinhibition or fixing possible plant damage.
At lower radiation intensities, the photosynthetic rate increases proportionally to the light intensity, until reaching a maximum value. At higher radiation conditions, the mechanism responsible for light collecting is supplied with an excessive number of photons, producing photoinhibition over light-dependent photosynthesis.
Measurement and characterization of PAR are interesting for different applications, such as gross primary production (GPP) [2], terrestrial net primary production (NPP) [3], and PAR use efficiency [4]. Due to its importance, and because of the lack of PAR surface direct measurements, its estimation from satellite measurements is usually required [5,6]; this estimation is usually based on empirical correlation models between satellite measurements and global horizontal radiation (GHI) surface direct measurements [7], as GHI is usually available at quite common radiometric stations. Then, modeling of PAR/GHI ratio is also required, as the portion of global irradiance that corresponds to the PAR spectrum.
However, different features affect PAR/GHI ratio. GHI is the total solar irradiance that reaches the Earth's surface measured on a horizontal plane, as the sum of the diffuse and direct horizontal irradiance. Diffuse horizontal irradiance is the solar irradiance scattered and reflected by the atmosphere, while direct horizontal irradiance is the solar irradiance coming directly from the sun; therefore, atmospheric processes affecting both irradiances are different, and those processes effects depend on each irradiance spectrum and the atmospheric conditions. The PAR/GHI ratio represents the portion of global irradiance that corresponds to the PAR spectrum, but diffuse and direct global irradiance contribution to PAR depends on the atmospheric conditions in a region.
Typically, the PAR spectrum can be estimated from the Kato bands included in the corresponding part of the solar spectrum [8]. Indeed, the spectrum can be divided into 32 bands (Kato bands) [9] of different widths, in each of which the absorption coefficient of different gases is almost constant. The solar spectrum region corresponding to PAR includes the Kato bands interval from 7 to 16, as shown in Table 1. Table 1. Kato bands in the photosynthetically active radiation (PAR) solar spectrum region.

Kato Band
Wavelength Region (µm) As PAR and GHI spectrums are different, the PAR/GHI ratio is affected by various atmospheric parameters [10][11][12], with regional dependence: relative humidity, water vapor pressure, and, especially, the clearness index (k t ) [13]. This index is defined as the quotient between the GHI that reaches the Earth's surface and the extraterrestrial global irradiance at the top of the atmosphere; its range goes from 0 to 1, so that clear atmospheric conditions corresponded to values close to 1. In all PAR/GHI ratio models, k t is a recurrent input variable, as it allows removing the deterministic component of the solar radiation due to astronomical influences and, also because of its significant cloudiness dependence, like the PAR/GHI ratio. Indeed, spectral attenuation due to clouds follows two different processes: scattering and absorption. Scattering attenuation is independent of the spectrum band (nonselective attenuation) and consequently the PAR/GHI ratio is not significantly affected by scattering in clouds. However, clouds' water absorption is strong in the infrared range, which increases GHI attenuation but not that of PAR. Accordingly, typical PAR/GHI ratios should be higher in cloudy areas (with a lower k t ) than in low cloudiness areas.
Therefore, both the PAR/GHI ratio and the index k t are clearly influenced by the location and along the time as, they are strongly associated with the cloudiness, which is linked to the meteorological conditions. In PAR/GHI ratio estimation, the climate of a region is relevant.
PAR modeling is a common approach. Many works [34][35][36][37][38][39] have previously addressed PAR modeling in different regions (from Costa Rica, China, Thailand, Arabian Sea, Spain, and also Spain, respectively) from satellite estimates. [38,39] used a clustering analysis from the datasets provided by the Climate Monitoring Satellite Application Facility (CM-SAF). While in [38] only GHI was used as input, in [39] the PAR results were improved adding more variables such as global extraterrestrial irradiance, relative humidity and temperature. In this spatial modeling, linear regression provided very similar results than more complex techniques, such as neural networks. Moreover, a site-adaptation technique (using both surface and satellite data) was used to achieve better results at the oceanic stations. Additionally, the models developed in [39] with Spain datasets were compared at several stations in Spain with other models from the bibliography, developed for other locations: the results from Spain datasets model were clearly better, as expected.
The use of different surface datasets, including different PAR-related parameters, was deeply considered, but PAR spatial modeling is also dependent on the satellite data accuracy. Therefore in this work, PAR/GHI models obtained from two different satellite-derived datasets (Climate Monitoring Satellite Application Facility (CM-SAF) and Moderate-Resolution Imaging Spectroradiometer (MODIS)) were applied to develop and apply PAR/GHI models, in order to estimate the satellite data accuracy impact in PAR estimation performance. It is assumed that PAR/GHI satellite data estimations can be applied to surface measurements, in order to spatially extend surface PAR estimations. With this assumption, GHI surface measurements can be used to obtain PAR, providing more accurate results. Following this approach, in [38,39] only CM-SAF data were applied, while in this work an alternative satellite dataset (MODIS) was applied, in order to compare both results.
The possible use of MODIS datasets is particularly interesting due to the different schemes applied for obtaining MODIS irradiance products [40]. Indeed, MODIS GHI and PAR data are generated from a look-up table (LUT)-based shortwave irradiance retrieval algorithm. For this LUT, atmospheric radiative properties are estimated by directly taking the satellite top-of-atmosphere (TOA) signature, so the LUT approach is performed at the native resolution of the satellite data, and fewer data gaps are produced. However, this scheme does not take advantage of the multispectral information of the MODIS sensor, so when the difference between the surface reflectance and TOA reflectance is small, MODIS GHI and PAR data accuracy can be compromised.
In this work, PAR models' results are presented and compared over a geographical territory (mainland Spain) with different climates [41][42][43][44]. According to Figure 1, four main climates (oceanic, continental Mediterranean, standard Mediterranean and mountain climates) can be distinguished in mainland Spain, though the first three ones cover most of the territory. Thus, PAR radiometric station measurements located at those three climatic zones are analyzed. Geographical locations and main characteristics of these climatic zones follow:

•
The oceanic climate is mainly present in the northern and northwestern parts of mainland Spain. This climate is characterized by mild temperatures and the westerlies (prevailing winds from the west in middle latitudes) that bring frequent clouds and precipitation, with precipitation events quite uniformly distributed throughout the year. • Mediterranean climates are common in the rest of mainland Spain, including eastern, southern and central areas of the Iberian Peninsula. However, there are no strict boundaries between oceanic and Mediterranean climates, while the transition between them is smooth: the climate gradually shifts as moving toward southern and southeastern Spain. Mediterranean climates have hot summers, and precipitation events are rare, usually as storms. In this climate, two subtypes can be distinguished: 1. Standard Mediterranean climate. This subtype is present in southern and eastern mainland Spain, and it is characterized by dry and hot summers and mild winters due to sea influence.  • The oceanic climate is mainly present in the northern and northwestern parts of mainland Spain. This climate is characterized by mild temperatures and the westerlies (prevailing winds from the west in middle latitudes) that bring frequent clouds and precipitation, with precipitation events quite uniformly distributed throughout the year. • Mediterranean climates are common in the rest of mainland Spain, including eastern, southern and central areas of the Iberian Peninsula. However, there are no strict boundaries between oceanic and Mediterranean climates, while the transition between them is smooth: the climate gradually shifts as moving toward southern and southeastern Spain. Mediterranean climates have hot summers, and precipitation events are rare, usually as storms. In this climate, two subtypes can be distinguished: 1. Standard Mediterranean climate. This subtype is present in southern and eastern mainland Spain, and it is characterized by dry and hot summers and mild winters due to sea influence. 2. Continental Mediterranean climate. This subtype is typical in central areas of the Iberian Peninsula. In this subtype, the influence of the sea is weaker and the winters are colder. Thus, seasonal variations are more pronounced than in the Standard Mediterranean climate. As a summary, the main objectives of this work follow: • Analyzing the spatial and temporal variability of the PAR/GHI ratio, considering the main different climatic regions in mainland Spain. • Comparing PAR/GHI estimates provided by two satellite-derived datasets to determine their best variability. As a summary, the main objectives of this work follow: • Analyzing the spatial and temporal variability of the PAR/GHI ratio, considering the main different climatic regions in mainland Spain.

•
Comparing PAR/GHI estimates provided by two satellite-derived datasets to determine their best variability.

•
Developing PAR estimation models, using GHI as input, and based on each satellite-derived datasets.
Therefore, the influence of different climate types and different satellite-derived datasets as input data in PAR estimation over a wide territory is investigated.

Data
The datasets used in this work were daily measurements of PAR and GHI measured at three locations spread across the Iberian Peninsula (Tables 2 and 3):   June 22, 2017. Located in the southeast of Spain, its climate is the typical Mediterranean, with warm summers. Additionally, the specific location of this station presents the particular characteristic of being over a desert zone, where the presence of aerosols from the sand is very significant. Here, the mean annual temperature is 17.5 • C and the annual precipitation is 227 mm [45].
At the first two stations, the PAR was measured with an ML-020P (EKO Instruments) and the GHI with a pyranometer CM21 (Kipp and Zonen). At the third station, the PAR was measured with a GUV-2511 multi-band radiometer (Biospherical Instruments Inc.) and the GHI data were taken from the measurements made by MeteoGalicia using a pyrradiometer Type 8111 (PH.SCHENK). The three sites were chosen because they are located in zones with different climatic characteristics [42].
Satellite data were provided by two different sources. First, mean daily data were collected from CM-SAF) [46], which is part of the EUMETSAT Satellite Application Facility network. CM-SAF downloaded data corresponded to the years from 1999 to 2011. As CM-SAF does not provide PAR data directly, the Spectral Resolved Irradiance (SRI) product was used in this work, which contains the Kato bands [9]. From these bands, the GHI was also obtained. MODIS products provided GHI and PAR data as well, in particular data from January 1, 2018, to May 31, 2019, from MCD18A1-MODIS/Terra+Aqua Surface Radiation Daily/3-Hour L3 Global 5 km SIN Grid [47] and MCD18A2-MODIS/Terra+Aqua Photosynthetically Active Radiation Daily/3-Hour L3 Global 5 km SIN Grid [48].
The data downloaded from both CM-SAF and MODIS covered a grid from 35.3 • N to 44.0 • N in latitude and from 9.5 • W to 3.5 • E in longitude. The grid from CM-SAF had a resolution of 0.1 • × 0.1 • (11.1 km × 11.1 km), while the grid of MODIS had a resolution of 5 km (0.045 • ). These grids corresponded to the location of the Iberian Peninsula where mainland Spain is located.

Methodology
The study was carried out in the following steps: 1. Obtaining the clearness index: The clearness index (k t ) was obtained at the three study stations from ground measurements. For that, GHI data measured at these stations were converted into a Remote Sens. 2020, 12, 1262 7 of 22 normalized index by dividing GHI by the extraterrestrial radiation, that is, the radiation at the top of the atmosphere. Thus, where G 0 is the extraterrestrial radiation, which is obtained from where G 0n is the normal extra-terrestrial radiation and θ z is the zenith angle. G 0n was obtained, at the same time, from the solar constant and the Julian day, by the following expression: where G sc is the solar constant and n is the Julian day. 2. Analyzing the spatial variability: To analyze the spatial variability, the PAR/GHI ratio was plotted vs. k t at the three study stations from ground measurements and both satellite-derived datasets (CM-SAF and MODIS).
3. Analyzing the temporal variability: To study the temporal variability, the daily evolution of PAR/GHI and k t from ground measurements at the three stations was plotted. The evolution of the monthly means of these variables was also plotted, using ground measurements and satellite estimates. In this section is assumed that for both variables each month of the year constitutes a congruence class of data (modular arithmetic mod12).
4. Collecting satellite estimates: On the other hand, satellite data from two different sources (CM-SAF and MODIS) were collected for the study area (mainland Spain), which included the three analyzed stations. In the case of CM-SAF, it does not provide PAR data directly, but it does provide the SRI product, which contains the irradiance corresponding to each Kato band. Thus, to obtain PAR estimates, we added the irradiance of the bands corresponding to the spectral region of 400-700 nm (from band 7 to 16, according to Table 1). Likewise, though CM-SAF directly provides GHI estimates, the same methodology of adding bands was followed to obtain both products in the same way. This process was not necessary in the case of MODIS data, as it directly provides PAR and GHI estimates.
5. Obtaining the PAR models from CM-SAF and MODIS estimates: PAR and GHI estimates from CM-SAF and MODIS were used to develop linear models that provide PAR estimates from GHI for each point of the study area and each month as well. Indeed, in [39], linear regressions provided very similar results to the use of more complex techniques such as neural networks, consuming less time in computation than other possible PAR models. So, the models used have the following expression: Thus, the coefficients α and β are the regression coefficients obtained from PAR and GHI data collected in each grid point for each month, so these coefficients provide the particular expression of the model which must be used to obtain PAR according to the site and month.
In addition, the determination coefficients for each regression were also calculated to show the correlation degree between PAR and GHI according to the linear regression, following these expressions: where σ PAR Modeled PAR Measured is the covariance of (PAR Modeled , PAR Measured ), σ PAR Modeled is the variance of PAR Modeled , σ PAR Measured is the variance of PAR Measured .
6. Validation: PAR estimated from CM-SAF-and MODIS-derived models were plotted vs. PAR measured at the three study stations. The accuracy obtained for the models in the different places was assessed by the following statistics: the slope of the regression line, the intercept, the determination coefficient (R 2 ), the mean bias error (MBE), the root-mean-square error (RMSE), and the standard error of the regression (SE). Errors are given as Finally, the statistics obtained at the three stations were compared with the following models proposed in [43].

1.
Considering only GHI as input.

2.
Adding more variables, such as GHI, global extraterrestrial irradiance, relative humidity, and temperature.

3.
Applying site-adaptation technique (only used for Santiago-EOAS station). Figure 2 shows the relationship between PAR/GHI and k t [49], from both ground measurements and satellite estimates. Focusing on the ground measurements, at the three stations considered, there was a reduction in the relationship of PAR/GHI when k t increased. This decrease was very mild, especially for a higher clearness index, fitting the PAR/GHI ratio from ground measurements to a logarithmic function of k t as proposed by [10,24,26,29]. The fitting expressions obtained for the different stations are included in Figure 2. According to these expressions, the logarithmic term was higher at Santiago-EOAS and, thus, the rate of decrease of the PAR/GHI ratio was higher for this station (the slope of Santiago-EOAS: −0.042; the slope of PSA-CIEMAT: −0.035; the slope of CEDER-CIEMAT: −0.017). Furthermore, the Mediterranean stations (PSA-CIEMAT and CEDER-CIEMAT) had less dispersion than the oceanic one; that is, the PAR/GHI ratio always presented more similar values for a given value of the clearness index. The means and standard deviations as error bars in bins of 0.1 k t are also plotted in the figure to illustrate better the spread of the data. On the other hand, the mean values of PAR/GHI were generally higher at Santiago-EOAS than in the other two stations, especially for lower k t values. Further, high clearness indexes were more abundant in these Mediterranean stations. In addition, the Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 23 (c)

Temporal Variability
Due to the known relation between the PAR/GHI ratio and kt, the temporal variability of this ratio is also influenced by the temporal variability (by the minute, hour, day, etc.) of kt. In Figure 3 the temporal evolution of measured daily data of kt and PAR/GHI is plotted for the three study stations. The graphs show the high temporal variability of kt at all stations, while the PAR/GHI ratio remained more stable at the Mediterranean stations. These results can be explained according to the climatic characteristics of the zones in which the different stations are located. PSA-CIEMAT and CEDER-CIEMAT presented similar behavior as both are located in zones with a Mediterranean climate (although of different varieties). In the Mediterranean climate, which is less cloudy than the oceanic, the cloudiness indexes were higher, justifying as well the lower values of the PAR/GHI ratio, according to the greater capacity of absorption of water vapor in the infrared spectral area.
On the other hand, the smaller dispersion in this relation, observed in the Mediterranean stations, can be explained by the greater stability in terms of cloudiness, which is characteristic of the Mediterranean climate. This would explain that high levels of clearness were associated with cloudiness situations and lower PAR/GHI ratios. On the contrary, the greater instability of the oceanic climate is advantageous for greater variability of the PAR/GHI ratio. Indeed, the continuous entry of fronts loaded with cloudiness and precipitation generates greater variability of atmospheric conditions and PAR/GHI values, which can give rise to similar situations of clearness. All this would explain the greater dispersion observed at the Santiago-EOAS station.

Temporal Variability
Due to the known relation between the PAR/GHI ratio and k t , the temporal variability of this ratio is also influenced by the temporal variability (by the minute, hour, day, etc.) of k t . In Figure 3 the temporal evolution of measured daily data of k t and PAR/GHI is plotted for the three study stations. The graphs show the high temporal variability of k t at all stations, while the PAR/GHI ratio remained more stable at the Mediterranean stations. In Table 4 the mean values and the coefficients of variation (CV) for PAR/GHI and kt series were shown. CV is defined as the quotient between the standard deviation and the mean value; thus, it is a normalized index of the degree of variability. According to Table 4, the differences between the Santiago-EOAS PAR/GHI series and the other two series are very clear. Indeed, regarding PAR/GHI mean values, the Santiago-EOAS value was higher than those corresponding to the other two series. Similarly, the highest CV of the PAR/GHI series was also obtained at Santiago-EOAS station.
In the case of kt, the CV is very similar in the three series, though the mean value is clearly higher in Santiago-EOAS and lower in PSA-CIEMAT. Focusing on the two Mediterranean stations, PAR/GHI mean values were slightly higher (43.44 vs. 42.17) at the station of the standard Mediterranean climate (PSA-CIEMAT) than at the station of the continental Mediterranean climate (CEDER-CIEMAT). This difference can be explained by the greater quantity of aerosols, from the desert and the sea, which is characteristic of the zone where the PSA-CIEMAT station is situated, as it can be seen in [50]. According to this work, the aerosol optical In Table 4 the mean values and the coefficients of variation (CV) for PAR/GHI and k t series were shown. CV is defined as the quotient between the standard deviation and the mean value; thus, it is a normalized index of the degree of variability. According to Table 4, the differences between the Santiago-EOAS PAR/GHI series and the other two series are very clear. Indeed, regarding PAR/GHI mean values, the Santiago-EOAS value was higher than those corresponding to the other two series. Similarly, the highest CV of the PAR/GHI series was also obtained at Santiago-EOAS station. In the case of k t , the CV is very similar in the three series, though the mean value is clearly higher in Santiago-EOAS and lower in PSA-CIEMAT.
Focusing on the two Mediterranean stations, PAR/GHI mean values were slightly higher (43.44 vs. 42.17) at the station of the standard Mediterranean climate (PSA-CIEMAT) than at the station of the continental Mediterranean climate (CEDER-CIEMAT). This difference can be explained by the greater quantity of aerosols, from the desert and the sea, which is characteristic of the zone where the PSA-CIEMAT station is situated, as it can be seen in [50]. According to this work, the aerosol optical depth (at 550 nm) derived from the MODIS satellite for the 2003-2012 period was clearly higher in the zone of PSA-CIEMAT than in the other two study zones. Thus, this higher level of aerosols at the PSA-CIEMAT station would increase, as we have explained, the PAR/GHI ratio at this location compared with the other Mediterranean station (CEDER-CIEMAT). This higher value of PAR/GHI at the PSA-CIEMAT station located in the Mediterranean coastal zone was observed even though the lower cloudiness in this zone produced clearness indexes slightly higher than those of the continental Mediterranean climate (0.64 vs. 0.58). Indeed, since clear atmospheric conditions corresponding to k t values close to 1, the lower presence of clouds at PSA-CIEMAT entailed higher k t at this station than at CEDER-CIEMAT. Moreover, the higher variability of k t at the CEDER-CIEMAT station (0.33) compared with the PSA-CIEMAT station (0.24) could be counteracted by the influence of aerosols at the PSA-CIEMAT station, bringing less variability of PAR/GHI at CEDER-CIEMAT (0.03 vs. 0.04).
The relationship between the PAR/GHI ratio and k t , again with both ground measurements and satellite estimates, is also shown by the annual evolution of the monthly means (Figure 4), where the features formerly explained can be seen. As expected, the values of k t were higher in the summer months (with more clearness) than in the winter months. In addition, the monthly evolutions of both variables were, in general, complementary at all stations. Indeed, when the clearness index increased/decreased, the PAR/GHI ratio, in general, decreased/increased. The relationship between the PAR/GHI ratio and kt, again with both ground measurements and satellite estimates, is also shown by the annual evolution of the monthly means (Figure 4), where the features formerly explained can be seen. As expected, the values of kt were higher in the summer months (with more clearness) than in the winter months. In addition, the monthly evolutions of both variables were, in general, complementary at all stations. Indeed, when the clearness index increased/decreased, the PAR/GHI ratio, in general, decreased/increased.

Performance of the Photosynthetically Active Radiation (PAR) Models
The previous sections showed that PAR models must take into account spatial and temporal variability in their performances. Other previous works [21,24,25,31] conducted PAR modeling for specific locations and, thus, those models did not take into account spatial variability. In this sense, satellite-derived data provide a very suitable tool to develop these models, as they are available with high resolution for wide regions and long periods. As mentioned, two satellite-derived datasets (CM-SAF and MODIS) were used to develop models to estimate PAR as a function of GHI for a region with different climates (oceanic and Mediterranean).
We assumed that the relationships between PAR and GHI obtained from satellite estimations can be applied for ground measurements, so these relationships can be used as models to estimate PAR from GHI ground measurements, as in [38,39]. Thus, PAR and GHI estimates from CM-SAF and MODIS were used to establish linear models that provide PAR from GHI for each spatial point and

Performance of the Photosynthetically Active Radiation (PAR) Models
The previous sections showed that PAR models must take into account spatial and temporal variability in their performances. Other previous works [21,24,25,31] conducted PAR modeling for specific locations and, thus, those models did not take into account spatial variability. In this sense, satellite-derived data provide a very suitable tool to develop these models, as they are available with high resolution for wide regions and long periods. As mentioned, two satellite-derived datasets (CM-SAF and MODIS) were used to develop models to estimate PAR as a function of GHI for a region with different climates (oceanic and Mediterranean).
We assumed that the relationships between PAR and GHI obtained from satellite estimations can be applied for ground measurements, so these relationships can be used as models to estimate PAR from GHI ground measurements, as in [38,39]. Thus, PAR and GHI estimates from CM-SAF and MODIS were used to establish linear models that provide PAR from GHI for each spatial point and month, as noted in Equation (5). Thus, this methodology provides two models for each grid point: one from CM-SAF and the other from MODIS. Figures 5-10 show the coefficients α and β, that is, the model corresponding to each point of the study area and each month, from both datasets (CM-SAF and MODIS). Regarding the physical reading of α and β coefficients, α is the percentage of PAR in GHI and thus, it should be higher as k t decreases. While the β coefficient is an adjustment correction term. Figures 7 and 10 show the determination coefficients (R 2 ) for each point, indicating the degree of correlation existing between the satellite-derived PAR estimations and the model for this point. According to these figures, most of the determination coefficients were higher than 0.95, except in July for the MODIS-estimate-derived model, which had lower coefficients in a region of the center and south of the study area. Indeed, in these cases, the correlation between PAR and GHI (measured by the determination coefficient) is quite lower than in the rest of the cases. Thus, the coefficients α and β also showed very different values. MODIS performance characteristics could explain these anomalous results. As it was mentioned, MODIS accuracy to estimate radiation may be compromised in clear sky conditions, so the relation between PAR and GHI could not be captured properly during July in the south and center of mainland Spain.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 23 anomalous results. As it was mentioned, MODIS accuracy to estimate radiation may be compromised in clear sky conditions, so the relation between PAR and GHI could not be captured properly during July in the south and center of mainland Spain. Thus, by using satellite estimations, these linear models can be suitable tools to estimate PAR from GHI measurements.  Thus, by using satellite estimations, these linear models can be suitable tools to estimate PAR from GHI measurements.
To check the degree of accuracy of the proposed models, they were applied at the three stations used in this work, for which we used PAR and GHI ground measurements. Figure 2 (PAR/GHI vs. k t ) shows that we added the values provided by CM-SAF and MODIS for comparison with the values obtained from pyranometer measurements, as well as the plots of the error bars in bins of 0.1 k t . In the case of Santiago-EOAS, MODIS estimates showed behavior more similar to the ground measurements than CM-SAF estimates, as these last estimates more poorly captured the high dispersion observed in the ground data. On the other hand, the mean values were very similar for both satellite-derived datasets, although MODIS provided a better estimate for lower k t . At CEDER-CIEMAT, according to the mean values, MODIS usually overestimated more than CM-SAF, while CM-SAF was in good agreement with the ground measurements for values of k t higher than 0.5. According to the error bars, the spread obtained from CM-SAF estimates was less than that from MODIS estimates and more similar to that obtained from measurements. At PSA-CIEMAT, MODIS overestimated the dispersion of the measured data more than CM-SAF, the mean values of which were, in general, more similar to those of the ground measurements. The LUT approach used to obtain MODIS products can explain its worse performance for the Meditteranean case. Indeed, as explained, for this LUT, the atmospheric radiative properties were estimated by directly taking the satellite TOA signature. Therefore, the performance at the native resolution of the satellite data was obtained and, thus, fewer data gaps were produced, which improved the estimates. However, in zones where the contrast between TOA reflectance and surface reflectance is more limited, as in the Mediterranean case, the lack of discrimination can result in worse estimates.        Additionally, the monthly evolutions of k t and PAR/GHI at the study stations determined by using the models obtained from both datasets are included in Figure 4, along with the corresponding evolutions of measured data. According to this figure, at the Mediterranean stations (CEDER-CIEMAT and PSA-CIEMAT) the monthly evolution of both variables was better captured by CM-SAF estimates than those from MODIS. However, in the case of the Santiago-EOAS station (oceanic climate), regarding the PAR/GHI ratio, MODIS behaved more similarly to ground measurements than CM-SAF did, due to the LUT approach used, which provided greater performance accuracy.
In addition, to establish more clearly the accuracy of the two models, PAR values estimated from CM-SAF-and MODIS-data-derived models are plotted vs. PAR measured in Figures 11-13. The statistics used to assess the accuracy were the slope, the intercept, the R 2 , the MBE, the RMSE, and the SE. According to the results, for the oceanic station, the MODIS-data-derived model provided better results than the model developed from CM-SAF estimates. On the contrary, the model developed from CM-SAF provided better results at the two Mediterranean stations. Indeed, the statistics obtained for the best models can be considered to be of high quality, especially in the Mediterranean stations, since the slopes of the scatterplots were 0.99 and 1 and the RMSE was around 2% for both stations. In the case of the oceanic station, the statistics were slightly worse, which was expected due to the high variability of this climate, and yet, its MBE was similar to that of PSA-CIEMAT, its R 2 was higher than 0.99, and its slope was 0.96. Regarding the standard errors of the regression, its value for the model from CM-SAF is lower clearly at Mediterranean stations (8.22 and 12.41 µmol m −2 s −1 respectively) than in the oceanic one (31.80 µmol m −2 s −1 ). The opposite occurred with the model from MODIS, where the oceanic station obtained lower SE (16.81 µmol m −2 s −1 ) than the Mediterranean ones (29.05 and 25.48 µmol m −2 s −1 ). This fact is in agreement with the previous results according to which MODIS works better at the oceanic station and CM-SAF at the Mediterranean ones. In this sense, these errors can help to know more exactly the accuracy of the different models.  To check the degree of accuracy of the proposed models, they were applied at the three stations used in this work, for which we used PAR and GHI ground measurements. Figure 2 (PAR/GHI vs. kt) shows that we added the values provided by CM-SAF and MODIS for comparison with the values obtained from pyranometer measurements, as well as the plots of the error bars in bins of 0.1 kt. In the case of Santiago-EOAS, MODIS estimates showed behavior more similar to the ground measurements than CM-SAF estimates, as these last estimates more poorly captured the high dispersion observed in the ground data. On the other hand, the mean values were very similar for both satellite-derived datasets, although MODIS provided a better estimate for lower kt. At CEDER-CIEMAT, according to the mean values, MODIS usually overestimated more than CM-SAF, while CM-SAF was in good agreement with the ground measurements for values of kt higher than 0.5. According to the error bars, the spread obtained from CM-SAF estimates was less than that from MODIS estimates and more similar to that obtained from measurements. At PSA-CIEMAT, MODIS overestimated the dispersion of the measured data more than CM-SAF, the mean values of which were, in general, more similar to those of the ground measurements. The LUT approach used to obtain MODIS products can explain its worse performance for the Meditteranean case. Indeed, as explained, for this LUT, the atmospheric radiative properties were estimated by directly taking the satellite TOA signature. Therefore, the performance at the native resolution of the satellite data was obtained and, thus, fewer data gaps were produced, which improved the estimates. However, in zones where the contrast between TOA reflectance and surface reflectance is more limited, as in the Mediterranean case, the lack of discrimination can result in worse estimates. Additionally, the monthly evolutions of kt and PAR/GHI at the study stations determined by using the models obtained from both datasets are included in Figure 4, along with the corresponding evolutions of measured data. According to this figure, at the Mediterranean stations (CEDER-CIEMAT and PSA-CIEMAT) the monthly evolution of both variables was better captured by CM-SAF estimates than those from MODIS. However, in the case of the Santiago-EOAS station (oceanic climate), regarding the PAR/GHI ratio, MODIS behaved more similarly to ground measurements than CM-SAF did, due to the LUT approach used, which provided greater performance accuracy.
In addition, to establish more clearly the accuracy of the two models, PAR values estimated from CM-SAF-and MODIS-data-derived models are plotted vs. PAR measured in Figures 11-13. The statistics used to assess the accuracy were the slope, the intercept, the R 2 , the MBE, the RMSE, and the SE. According to the results, for the oceanic station, the MODIS-data-derived model provided better results than the model developed from CM-SAF estimates. On the contrary, the model developed from CM-SAF provided better results at the two Mediterranean stations. Indeed, the statistics obtained for the best models can be considered to be of high quality, especially in the Mediterranean stations, since the slopes of the scatterplots were 0.99 and 1 and the RMSE was around 2% for both stations. In the case of the oceanic station, the statistics were slightly worse, which was expected due to the high variability of this climate, and yet, its MBE was similar to that of PSA-CIEMAT, its R 2 was higher than 0.99, and its slope was 0.96. Regarding the standard errors of the regression, its value for the model from CM-SAF is lower clearly at Mediterranean stations (8.  Finally, in Table 5, these statistics at the three stations can be seen and compared with three of the proposed models in [39]. 1. Considering only GHI as input (model F-C1) 2. Adding more variables, such as GHI, global extraterrestrial irradiance, relative humidity, and temperature (model F-C2) 3. Applying site-adaptation technique (only used for Santiago-EOAS station) (model F-C3) Finally, in Table 5, these statistics at the three stations can be seen and compared with three of the proposed models in [39].
1. Considering only GHI as input (model F-C1) 2. Adding more variables, such as GHI, global extraterrestrial irradiance, relative humidity, and temperature (model F-C2) 3. Applying site-adaptation technique (only used for Santiago-EOAS station) (model F-C3) Finally, in Table 5, these statistics at the three stations can be seen and compared with three of the proposed models in [39].
Applying site-adaptation technique (only used for Santiago-EOAS station) (model F-C3) For Mediterranean stations, F-C1 and CM-SAF obtained similar results and both were clearly better than those of MODIS. It was necessary to add more variables as inputs (F-C2) to slightly improve the results. In the case of the oceanic station, MODIS provides much better results than CM-SAF, F-C1, and F-C2 models. It was necessary to apply a site-adaptation technique, using ground data in addition to satellite data, to improve the performances. Table 5. Statistics of the models proposed in this work (CM-SAF and MODIS) and of the proposed models in [43] (F-C1, F-C2, and F-C3) for each validation station.

Conclusions
The spatial and temporal variability of the PAR/GHI ratio is conditioned, in part, by the k t index. Indeed, the PAR/GHI ratio decreased, in general, as k t increased. The dispersion of PAR/GHI for any given value of k t and the decrease of this ratio as k t increased were both higher in the oceanic case than in the Mediterranean ones. Regarding the differences within the zones with the Mediterranean climate, the continental variety showed, in general, lower values of the PAR/GHI ratio.
The temporal evolutions of PAR/GHI showed a higher mean value and variability of this ratio for the oceanic case. Within Mediterranean climates, the standard variety can have higher mean values of PAR/GHI than the continental, in spite of the standard variety usually also having higher clearness indexes due to the presence of aerosols at stations situated in deserts or close to the sea, as in the case of the PSA-CIEMAT station.
Spatial and temporal variability must be taken into account in PAR model performances; so, satellite-derived data are a suitable tool to develop these models.
The validation at the three study stations of the models proposed, demonstrated the better performance of the MODIS-data-derived model for the case of Santiago-EOAS, while the model obtained from CM-SAF was better for the Mediterranean stations.