Combination of Models to Generate the First PAR Maps for Spain

: This work addresses the development of a PAR model in the entire territory of mainland Spain. Thus, a speciﬁc model is developed for each location of the study ﬁeld. The new PAR model consists of a combination of the estimates of two previous models that had unequal performances in different climates. In fact, one of them showed better results with Mediterranean climate, whereas the other obtained better results under oceanic climate. Interestingly, the new PAR model showed similar performance when validated at seven stations in mainland Spain with Mediterranean or oceanic climate. Furthermore, all validation slopes ranged from 0.99 to 1.00; the intercepts were less than 3.70 µ mol m − 2 s − 1 ; the R 2 were greater than 0.988, while MBE was closer to zero percent than − 0.39%; and RMSE were less than 6.21%. The estimates of the PAR model introduced in this work were then used to develop PAR maps over mainland Spain that represent daily PAR averages of each month and a full year at all locations in the study ﬁeld. A.A.N., L.F.Z. O.W.; writing— review editing, F.F.-C., J.M.V., R.X.V., A.A.N., L.F.Z. and O.W.; visualization, R.X.V. and A.A.N.; supervision, J.M.V.; project administration, J.M.V. and R.X.V.; funding acquisition, J.M.V. and R.X.V.


Introduction
Photosynthetically active radiation (PAR) is the portion of the solar spectrum that ranges from 400 to 700 nm [1,2]. This range corresponds to the wavelengths that plants use to perform photosynthesis. In recent years, there has been a growing interest in measuring this variable due to its numerous applications in different fields of study, such as calculations involving gross primary production (GPP) or terrestrial net primary production (NPP) [3,4]. PAR is also an important variable for estimating the growth of biomass and microalgae [5][6][7][8][9][10][11]. Another example of PAR use is that when designing or choosing a greenhouse cover, the optical properties of the cover in the PAR range must be taken into account [12,13]. Therefore, PAR is an interesting variable when assessing the potential growth of a crop or when determining, in terms of agricultural productivity, which is the most suitable place.
Radiometric stations that provide PAR measurements are scarce and not always available at the required location. Thus, PAR estimates provided by satellite products. such as the Climate Monitoring Satellite Facility (CM-SAF) and the Moderate Resolution Imaging Spectroradiometer (MODIS). or by empirical models, are often used. For example, some authors have used Kato bands [14] to develop PAR models [15,16]. Others have developed multilinear PAR models that include global horizontal irradiance (GHI), global extraterrestrial irradiance, air temperature, relative humidity, clearness index, skylight brightness, solar elevation angle, solar zenith angle, dew point temperature, aerosol optical depth, air mass, total ozone column, total precipitable water, water vapor pressure, or saturation water vapor pressure as input variables [17][18][19][20][21][22][23][24][25][26]. Non-linear PAR models [27][28][29] have also been carried out in previous works.
The variability of solar irradiance and the spatial and temporal variability of the PAR/GHI ratio have previously been addressed [17,[30][31][32]. Solar irradiance depends on cloudiness and the presence of aerosols [33]. According to Ångström's law [34], the extinction of irradiance due to aerosols is higher for shorter wavelengths. Therefore, the ultraviolet and visible bands are the bands most affected by the presence of aerosols [35]. Consequently, PAR models are highly dependent on the climatic and atmospheric conditions in which they were developed. Therefore, PAR models only obtain accurate estimates in locations where the climate and atmospheric conditions are similar to those in which they were developed.
The present work focuses in PAR modeling in Spain. Previous studies have addressed this topic [17][18][19]36,37]. The main climates present in mainland Spain according to the Köppen-Geiger classification [38][39][40][41] are oceanic varieties in northern and northwestern Spain, Mediterranean varieties in central, eastern and southern areas, arid climates in the southeast, and mountainous climates in the main mountain ranges of the Iberian Peninsula.
In a previous study [17], two PAR models for mainland Spain were developed using two different datasets; the first model was developed using the estimates provided by CM-SAF, while the second was developed using the data set supplied by MODIS. According to [17], the results of the model that used CM-SAF estimates for its development were better in the Mediterranean and arid areas. By contrast, the best results in the oceanic areas were obtained with the model developed from the MODIS estimates.
In the present work, a new daily PAR model is presented that covers mainland Spain. This new model is developed using a linear combination of the estimates provided by the two PAR models elaborated in [17]. The model is validated at seven radiometric stations located on mainland Spain. Next, the estimates supplied by the new model are used to create the first monthly and annual PAR maps in mainland Spain.
This study is structured as follows. First, the data set employed and the mathematical tools used to develop the PAR model are described. Next, the results of the validation of the model with data provided by seven stations located throughout the study territory are presented. The monthly and annual PAR maps of mainland Spain developed using the estimates provided by the new PAR model are then shown. The results are discussed subsequently and, finally, the most significant conclusions are discussed.

Materials and Methods
The two models elaborated in [17] used estimates from two satellite datasets-CM-SAF and MODIS-that cover the location of the Iberian Peninsula, where mainland Spain is located. The grid of both data sets ranges from 35.3 • N to 44.0 • N in latitude and from 9.5 • W to 3.5 • E in longitude. The CM-SAF grid had a resolution of 0.1 • × 0.1 • , while the MODIS grid had a resolution of 5 km. In particular, the CM-SAF data set corresponded to the years 1999 to 2011 of the Spectral Resolved Irradiance (SRI) product, which belongs to the EUMETSAT Satellite Application Facility network [42,43]. On the other hand, the MODIS data set corresponded to the dates from 1 January 2018, to 31 May 2019 of the MCD18A1-MODIS/Terra+Aqua Surface Radiation Daily/3 Hour L3 Global 5 km SIN Grid [44] and MCD18A2-MODIS/Terra+Aqua Photosynthetically Active Radiation Daily/3-Hour L3 Global 5 km SIN Grid [45,46] products.
Both models developed in [17] have the following mathematical structure: where α and β are specific coefficients for each point of the grid corresponding to mainland Spain. Thus, both models are valid for the whole study territory, giving specific estimates for each point. For the present work, data from seven radiometric stations belonging to the GEOPAR Project (Project CGL2016-79284-P AEI/FEDER/UE) were used. Figure 1 shows the loca-tions of the stations that provided data on GHI, PAR, air temperature, and relative humidity. Further details of each station can be found in Table 1. where α and β are specific coefficients for each point of the grid corresponding to mainland Spain. Thus, both models are valid for the whole study territory, giving specific estimates for each point.
For the present work, data from seven radiometric stations belonging to the GEOPAR Project (Project CGL2016-79284-P AEI/FEDER/UE) were used. Figure 1 shows the locations of the stations that provided data on GHI, PAR, air temperature, and relative humidity. Further details of each station can be found in Table 1.  Depending on their climate, the stations can be classified into two groups: those that feature a humid climate (oceanic varieties) and those with a dry climate (Mediterranean  Depending on their climate, the stations can be classified into two groups: those that feature a humid climate (oceanic varieties) and those with a dry climate (Mediterranean varieties). In the first group are Álava-NEIKER, Asturias-SERIDA, and Lugo-USC. Albacete-ITAP, Córdoba-IFAPA, Salamanca-CIALE, and Zaragoza-Aula Dei belong to the second group.
As none of the GEOPAR project stations were located in the Spanish archipelagos, the study field was set to mainland Spain because the climate on the islands is quite particular and there was no way to validate the estimates of the PAR model in those places.
Daily averages of PAR and GHI data from radiometric stations were collected over a period of more than two years, specifically from 7 May 2019 to 30 June 2021. This data Remote Sens. 2021, 13, 4950 4 of 14 set was randomly divided into two subsets with the same number of recordings. The first subset was used to train the PAR model, whereas the second subset was used to validate the model.
To develop the PAR model, an analogous technique to that described in [47,48] was applied. This methodology consisted of a linear fit between the estimates of the models developed in [17] from MODIS and CM-SAF and the measurements at the radiometric stations, as Equation (2) indicates.
where the PAR MODIS Model is the estimates from the MODIS model, the PAR CM-SAF Model is the estimates from the CM-SAF model and a, b, and c are the fitting coefficients.
To obtain estimates from the MODIS and CM-SAF models, the coefficients of these models were linearly interpolated to the coordinates of the location of each station and then both models were fed with the GHI data measured at the stations to calculate the PAR estimates. In this way, it is possible to obtain the coefficients a, b, and c that correspond to each station. Therefore, the fitting coefficients obtained for the location of the stations were linearly extrapolated on a least-squares approximation of the gradient at the boundary to each location of the original data grid, in order to obtain a PAR model that covers the entire territory of mainland Spain.
To validate the new model, its estimates were compared with PAR measurements in the seven stations. Statistics such as the slope and intercept of the scatterplot between the PAR estimates and the PAR measurements, the coefficient of determination (R 2 ), the mean bias error (MBE) and the root mean square error (RMSE) were used to address the goodness of the models.
where, PAR Model are the estimates of the model, PAR Measured are the PAR measurements, n is the number of recordings, σ 2 ModelPAR MeasuredPAR is the covariance of the measured and modeled PAR, σ 2 ModelPAR is the variance of the modeled PAR, and σ 2 MeasurePAR is the variance of the measured PAR.
The daily PAR average of each month was then calculated using the estimates of the new model for every point of the grid corresponding to mainland Spain. Similarly, annual averages of PAR estimates were calculated. In order to make these calculations, monthly and annual averages of the CM-SAF and MODIS PAR estimates were used to feed the PAR model.
These monthly and annual PAR averages, estimated using the PAR model, were subsequently used to develop PAR maps over mainland Spain. To calculate these estimates from the PAR model estimates, it was necessary to also use estimates from the CM-SAF and MODIS datasets.
The daily data collected from the seven stations were filtered, eliminating any data that did not meet any of the following conditions. The PAR/GHI ratio (both variables in W/m 2 so that the ratio was dimensionless) was between 0.3 and 0.6, relative humidity between Remote Sens. 2021, 13, 4950 5 of 14 0 and 100%, air temperatures between −40 and 60 • C, and clearness indexes (k t ) between 0 and 1. These criteria were used to discard any data recorded under environmental conditions outside the tolerance range of the instruments and to discard any data whose values do not make physical sense (for example, PAR values in µmol m −2 s −1 should be higher than GHI values in W m −2 ).

Results
After applying the rejection criteria to all the data, the number of remaining records for each station is shown in Table 2.  Table 2 reveals that in Salamanca-CIALE, all the original recordings passed filtering. By contrast, the highest number of rejected recordings was in Lugo-USC. The number of original recordings was 785 in all stations, except in Asturias-SERIDA, where one day of data was lost.
The filtered data set for each station was randomly divided into two subsets, each containing the same number of recordings. The first subset was used to develop the PAR model, whereas the second subset was used to validate the model. The first subset was also used to calculate the estimates of the CM-SAF and MODIS models of [17] for the location of each station.
The next step was to develop the PAR model itself. Therefore, the coefficients a, b, and c were calculated for each station by multilinear fitting according to Equation (2), where the PAR estimates of the CM-SAF and MODIS models are fitted to real PAR measurements from the seven stations. Table 3 shows the fitting coefficients for each station. The fitting coefficients vary from station to station. Therefore, to develop a model that covers the entire field of study, these coefficients were linearly extrapolated to the surface of mainland Spain on a grid with a resolution of 5 Km from 35.3 • N to 44.0 • N in latitude and from 9.5 • W to 3.5 • E in longitude. The coefficients for each point on the grid are illustrated in Figures 2-4, respectively.
Zaragoza-Aula Dei 0.29 0.79 −9.77 The fitting coefficients vary from station to station. Therefore, to develop a m that covers the entire field of study, these coefficients were linearly extrapolated to surface of mainland Spain on a grid with a resolution of 5 Km from 35.3°N to 44.0° latitude and from 9.5°W to 3.5°E in longitude. The coefficients for each point on the are illustrated in Figures 2-4, respectively.      These coefficients needed to be validated with another data set before developing the PAR model for the entire mainland of mainland Spain. Therefore, a validation test was carried out using the validation data subset of the stations.
For this reason, the PAR estimates obtained using the validation data subset were compared to the real measurements from the stations. That is, the GHI data from the validation subset were utilized to calculate the assessments of the CM-SAF and MODIS models; then these assessments along with the corresponding coefficients a, b and c were used to calculate the estimates of the new PAR model for each station. Finally, the estimates of the new PAR model were compared with the PAR data measured at each station, as illustrated by the following scatterplots. Table 4 summarizes all the validation metrics at each station.  Figure 5 reveals the linear fit between the estimates of the PAR model and the measured PAR data, where the red line represents the linear fit. Indeed, the slope is close to the unit in all the stations, and the intercepts are close to zero, being 3.70 the highest intercept (Salamanca-CIALE). The R 2 obtained is greater than 0.99 in all cases except Lugo-USC and Salamanca-CIALE, where the R 2 is 0.988. The highest RMSE values were also reached at these two stations, with 6.21% and 5.75%, respectively. Regarding the MBE results, in all stations, the MBE was close to zero, −0.39% being the highest value (Albacete-ITAP). ured PAR data, where the red line represents the linear fit. Indeed, the slope is close to the unit in all the stations, and the intercepts are close to zero, being 3.70 the highest intercept (Salamanca-CIALE). The R 2 obtained is greater than 0.99 in all cases except Lugo-USC and Salamanca-CIALE, where the R 2 is 0.988. The highest RMSE values were also reached at these two stations, with 6.21% and 5.75%, respectively. Regarding the MBE results, in all stations, the MBE was close to zero, −0.39% being the highest value (Albacete-ITAP).

Developing PAR Maps
The PAR model developed in this study was then used to develop, for the first time, monthly and annual PAR models in mainland Spain. First, the monthly and annual aver-

Developing PAR Maps
The PAR model developed in this study was then used to develop, for the first time, monthly and annual PAR models in mainland Spain. First, the monthly and annual averages of the daily PAR estimates were calculated using the PAR model. To do that, the monthly and annual averages of the CM-SAF and MODIS datasets were calculated. As the resolution of both datasets is different (0.1 • × 0.1 • in the CM-SAF case and 5 km × 5 km in the MODIS case), the CM-SAF grid was resized to have the same resolution in both grids. This resolution makes this model suitable to conduct generalized or large studies, but not adequate for more detailed studies. In this way, it was possible to obtain the average daily PAR estimate for each point on the grid, for every calendar month and for a year. Next, these estimates were represented on the surface of mainland Spain, as illustrated in Figures 6 and 7. year. Next, these estimates were represented on the surface of mainland Spain, as illustrated in Figures 6 and 7.

Discussion
The PAR model introduced in this work has been validated in seven stations located throughout mainland Spain with different types of climate (see Table 1), according to the Köppen-Geiger classification. This model was elaborated as a linear combination of the two models developed in [17]. These two previous models exhibited different performances depending on the climate of the location. Thus, while the model developed from CM-SAF demonstrated good results in places with Mediterranean or dry climates, the model developed from MODIS obtained its better results in locations with oceanic or humid climates. The new proposed PAR model has been also compared with three models proposed by other authors [20,49].
The model proposed in [20] is described in the following equation: Different PAR models were proposed by [49] depending of the land use of the site. One of them was developed for pasture lands and is shown in Equation (9), while the expression for forest lands is shown in Equation (10). PAR = 2.023·GHI + 28.557 (9) PAR = 1.922·GHI + 3.630 (10) Table 5 illustrates the comparison between the results of the previous models and the new model, carried out with the validation data set of each station.

Discussion
The PAR model introduced in this work has been validated in seven stations located throughout mainland Spain with different types of climate (see Table 1), according to the Köppen-Geiger classification. This model was elaborated as a linear combination of the two models developed in [17]. These two previous models exhibited different performances depending on the climate of the location. Thus, while the model developed from CM-SAF demonstrated good results in places with Mediterranean or dry climates, the model developed from MODIS obtained its better results in locations with oceanic or humid climates. The new proposed PAR model has been also compared with three models proposed by other authors [20,49].
The model proposed in [20] is described in the following equation: Different PAR models were proposed by [49] depending of the land use of the site. One of them was developed for pasture lands and is shown in Equation (9), while the expression for forest lands is shown in Equation (10)   Furthermore, the results of the new model were similar at all stations: all the slopes were close to the unit; all the intercepts were close to zero, with 3.70 µmol m −2 s −1 being the highest value (at Salamanca-CIALE); similarly, all the MBEs were close to zero, with −0.39% being the highest value (at Albacete-ITAP); and the highest RMSE was 6.21% (at Lugo-USC). According to these results, the new model showed no clear tendency as evaluated with these datasets and its performance was similar on every station, regardless of the climate. This result is a consequence of the combination of two previous models. One of them performed well in Mediterranean climates, while the second model obtained good results for oceanic climates.
With regard to the PAR maps, they were developed using the new PAR model introduced in this study. As expected, PAR irradiance levels increase during the summer months and decrease during the winter months. However, some trends and features remain the same every month. For example, PAR maximums are always reached in the southeast and central areas of the Iberian Peninsula along with the Guadalquivir and Ebro valleys. However, minimums of PAR irradiance were reached in northern areas of the Iberian Peninsula and mountain ranges, such as the Pyrenees, the Central Range, the Iberian Range, and the Betic Range. Furthermore, the same features are noted on the annual PAR map, where the minimums and maximums of PAR irradiance are reached in the same areas. Surprisingly, there are locations on mountain ranges, particularly in the Pyrenees in winter months, that the model significantly underestimates, producing a result near zero.