Using Copernicus Atmosphere Monitoring Service (CAMS) Products to Assess Illuminances at Ground Level under Cloudless Conditions

: Natural daylight is recognized as an important variable in the energy performance of buildings. A method that estimates the global illuminance received on a horizontal surface at ground level and its direct component at normal incidence under cloudless conditions is presented. The method uses the k -distribution method and the correlated- k approximation to compute a set of clearness indices integrated over 13 spectral bands covering the range 380–780 nm. A spectral resampling technique, including a spectral disaggregation and a spectral linear interpolation, is applied to these indices for providing a detailed set of solar irradiances at 1 nm in spectral resolution over the whole range. Then, these are weighted by the standardized CIE action spectrum for human eye for assessing the illuminance. Inputs to the method include the total column contents of ozone and water vapor as well as aerosol optical properties produced by the Copernicus Atmosphere Monitoring Service. Estimates of illuminance were compared to high-quality 1 min measurements of illuminance that were collected from two experimental sites located in two different climatic zones. A slight overestimation is observed for the global illuminance: the bias is between +1 klx and +3 klx, i.e., between +1% and +4% in relative value. The root mean square error varies between 5 klx (8%) and 6 klx (9%). The squared correlation coefﬁcient ranges between 0.95 and 0.97. At the site providing the direct illuminance at normal incidence, the performance of the method is lower compared to global illuminance with a lower squared correlation coefﬁcient of 0.53. The bias, relative bias, RMSE, and rRMSE are +7 klx, +9%, 12 klx, and 15%, respectively. The uncertainty of the method is of the order of the uncertainty of the measurements. The method offers accurate estimates of illuminance in cloudless conditions at high spatial and temporal resolutions useful for construction industries and operators as well as thermal simulation tools for optimal building design strategies.


Introduction
Natural daylight radiation reaching the surface of the Earth has become an essential climatic variable in the human daily life. For instance, daylight radiation is well known as the primary source of our visual and thermal comfort. The maximal usage of daylight radiation in building design strategies contributes to the reduction of energy consumption Kato et al. [12] scheme is accurate and useful for representing clearness indices in each of the 13 KBs covering the daylight radiation under cloudless and cloudy sky conditions. Hence, the combination of the Kato et al. [12] scheme and inputs from CAMS is a good candidate for an accurate estimation of the illuminance.
Nevertheless, the KBs do not precisely overlay the daylight spectral range. A technique has been developed to overcome this issue and to take into account the CIE standardized action spectrum. This paper describes the method to derive the illuminance from the clearness indices in the 13 KBs, which are themselves calculated from runs of libRadtran with inputs from CAMS. The method is not entirely new, since it has already been experimented in the case of UV fluxes using the most improved version of absorption parametrization included in libRadtran and named katoandwandji ( [16,17]) and photosynthetically active radiation (PAR) [18,19], but it has never been tested in the case of illuminance nor in combination with a spectral response. This study is part of an extensive project whose overarching goal is to freely and publicly offer datasets of illuminance at any time and any location in an operational way similar to the McClear model [13,14] with inputs on atmosphere variables originating from CAMS. To achieve this project, a step-by-step approach was adopted in which the performance of the clear-sky procedure is separately assessed to be better able to later understand the performance of the modelling of the cloud effects. The accuracy of the presented method is evaluated at two sites providing high-quality 1-min measurements of illuminance under cloudless conditions. Overall, this study contributes to answering the growing need of construction industries and investors on the spatial and temporal availability of illuminance at ground level for optimal building design strategies.

Description of Measurements Used for the Validation
Measurements of illuminance used for validation were collected from two stations: Vaulx-en-Velin and a station of National Renewable Energy Laboratory (NREL) in Golden, Colorado, located in France and the United States, respectively. Table 1 gives the geographical coordinates, the altitude of the CAMS cell containing the station from which atmospheric inputs are extracted, and the temporal period of measurements. At both stations, measurements are illuminances provided every minute. The station of Vaulx-en-Velin is located in the eastern part of Lyon urban community and is part of the International Daylight Measurement Program network of the CIE. The environment surrounding the station is made up of 70% urban housing and 30% of cultivated fields and parks. The climate is temperate with a maritime influence. Diffuse and global illuminances on horizontal surface were measured by two independent LMT BAP 30 FCT photometric sensors. One of these sensors is equipped with a shadow ring for measuring the diffuse component. The relative uncertainty of the global or diffuse illuminance is estimated at 5% ( [20]). These measurements are freely/publicly accessible and can be downloaded from the website http://idmp.entpe.fr/mesfr.htm, last access: 1 December 2020. The direct illuminance at normal incidence L BN is derived from the global and diffuse illuminances L G and L D measured on a horizontal surface as follows: where θ s is the solar zenith angle, which is calculated with the SG2 algorithm [21]. The station at Golden is a high-altitude site experiencing a quite dry climate located in the Colorado state of the United States. The station field is flat and is covered by natural grasses without trees. Only measurements of the global illuminance on a horizontal surface were available and measured by a Global Photometric sensor model LI-210. These measurements have been downloaded from the website midcdmz.nrel.gov/apps/sitehome. pl?site=BMS, accessed on 1 December 2020. The relative uncertainty of the measurements is approximately 8% ([22]).
In addition, measurements of broadband diffuse and global irradiances received on a horizontal surface and of direct irradiance at normal incidence were collected at both stations with 1 min of temporal resolution. With the assumption that cloudless instants identified by scrutinizing broadband irradiances are also cloudless instants in the illuminance measurements, the three datasets of broadband irradiances have served to select the cloudless instants instead of a single dataset of global illuminance at Golden or two datasets at Vaulx-en-Velin. When applied to the time series of the broadband direct, diffuse, and global irradiances at each site, the algorithm of Lefèvre et al. [13] yields a series of detected cloudless instants. The first of the two filters in this algorithm only retains those values for which the ratio of the diffuse to the global irradiance is under 0.3. The second filter uses the product of the clearness index KT by a typical air mass and inspects the temporal variability of this quantity, which should be steady for three hours under cloudless conditions. It may happen that in certain circumstances, the illuminance may be affected by scattered cloudiness, which may go undetected in the broadband measurements. Consequently, the retained series of cloudless periods may include cloudy periods for the illuminance. Given the high selectivity of the Lefèvre et al. [13] algorithm, the authors believe that such cases are rare and that the conclusions are unaffected as a whole. An additional test rejects cloudless instants for which the measurements of the global or diffuse illuminance are less than or equal to 0.

Description of the Method
In a nutshell, the method applies a spectral resampling technique, including a spectral disaggregation and a spectral linear interpolation, to the 13 KB-integrated clearness indices for providing a detailed set of clearness indices at 1 nm of spectral resolution over the whole range 380-780 nm. These detailed clearness indices are converted into spectral irradiances and then weighted by the standardized CIE action spectrum for human eye for assessing the global illuminance on horizontal surface and the direct illuminance at normal incidence.

Data Exploited by Libradtran
Under cloudless conditions, illuminances at the surface level depend mostly on the solar zenithal angle θ s , aerosol optical properties that may be characterized by the Angström exponent, aerosol type and aerosol optical depth (AOD), total amount of water vapor (TWV) and ozone (TOC), ground albedo, vertical profiles of temperature, pressure, density, volume mixing ratio for gases as a function of altitude, and altitude of the ground above mean sea level. As previously mentioned, this study participates in the development of an operational method, and the sources of data to be input to libRadtran should be chosen to ease the estimation of illuminance at any location and any time. A convenient way to fulfill this operational constraint is the exploitation of the atmospheric products delivered by CAMS. The SoDa service (http://www.soda-pro.com/, accessed on 1 December 2020) provides an access to the CAMS aerosol optical properties together with TOC and TWV. θ s is calculated with the SG2 algorithm [21]. libRadtran offers several solar spectra; that of Gueymard [23] was selected for computing the solar spectral irradiance Io Nλ received at the top of atmosphere at normal incidence at any time, λ being the wavelength. Air Force Geophysics Laboratory (AFGL) vertical profiles are used and selected at any location using the map of Gschwind et al. [14]. The Shuttle Radar Topography Mission dataset is used to extract the ground altitude.
The albedo of the ground in the daylight spectral band denotes the portion of the incident sunlight reflected by the ground in this band and is unknown in most cases. To solve this problem, it was assumed that the "daylight" albedo is equal to the PAR albedo because their spectral ranges are fairly similar. The acronym PAR stands for Photosynthetically Active Radiation, which is the part of solar radiation that lies in the wavelength range of 400-700 nm. The PAR albedo itself is assumed to be 0.47 times the broadband albedo [24]. The broadband albedo is defined as the integral of the bidirectional reflectance distribution function (BRDF), depending on the surface type and its roughness. Here, the broadband albedo is given by the series of maps of Blanc et al. [25] that provide the MODIS-derived BRDF parameters for each calendar month with no missing values at a spatial resolution of 0.05 • . One-min values of the inputs listed above are conveniently retrieved by machine-to-machine requests made in the verbose mode to the McClear web service on the Soda website (Gschwind et al. [26], http://www.soda-pro.com/, accessed on 1 December 2020).

Spectral Resampling Technique
Let L G and G λ be the global illuminance and the global spectral irradiance received on a horizontal surface at ground level, where λ is the wavelength (in nm). L G is given by: where K max is a constant defined as the maximum luminous efficacy and equal to 683 lm W −1 , and S λ is the standardized CIE action spectrum for human eye. In the same way, the direct illuminance at normal incidence L BN is given by: where B Nλ is the direct spectral irradiance at normal incidence. Let KT λ and KT Bλ be the spectral clearness index and the direct clearness index. They are given by: The libRadtran runs provide the clearness indices KT KBi and KT B_KBi in each of the 32 spectral bands. Altogether, the thirteen KBs, from KB 6 to KB 18 , do not precisely overlay the daylight spectral range 380-780 nm, since the KB 6 363-408 nm and KB 18 743-791 nm only partially cover the daylight spectral range.
The method that has been developed overcomes this issue and takes into account the CIE standardized action spectrum. The driving idea is to find a reasonable small number of spectral bands of 1 nm in width, noted FB j for fine band, whose spectral clearness indices KT FBj and KT B_FBj are calculated from the integrated ones in the KB i by affine functions. This is the spectral disaggregation step. Then, a complete and detailed set of clearness indices at 1 nm of spectral resolution over the whole range 380-780 nm is obtained by a linear interpolation of these FB j clearness indices. The 1 nm clearness indices are converted into 1 nm irradiances G λ and B Nλ that are then weighted by the standardized CIE action spectrum for human eye for assessing the global illuminance L G on horizontal surface and the direct illuminance at normal incidence L BN .
By assuming that the spectral clearness indices and the irradiances are constant over 1 nm, the integrals in Equations (2) and (3) can be approximated by finite sums, called Riemann sums. For instance, L G and L BN can be computed as follows: or using the clearness index KT n and the direct clearness index KT Bn : The spectral resampling technique has been described in Wandji Nyamsi et al. [16,18,19] for the UV and PAR ranges respectively and is outlined here for a better understanding. The proposed method is a pure modeling concept with radiative transfer simulations made with libRadtran. No measurements have been used for its development. First, a set of 60,000 atmospheric condition parameters under cloudless conditions has been built with Monte Carlo draws, following the statistical distribution of each input, as reported in Table 2 in [16]. For each condition, libRadtran is run twice for both the direct and global irradiances: one with the detailed spectral calculations every 1 nm and the other with the Kato et al. [12] scheme. Then, the irradiances are converted into clearness indices in order to eliminate the influence of the daily and annual variations in θ s and the solar spectrum. For both the global and direct irradiances, this ensemble of runs provides two sets of clearness indices: the detailed indices at 1 nm resolution and the indices spectrally integrated over each KBi.
For a given KBi, 2D histograms were drawn for the 1 nm clearness indices KT n , respectively KT Bn , and KT KBi , respectively KT B_KBi , for the range 380-780 nm. A visual inspection of each 2D histogram clearly shows a straight line with a squared correlation coefficient greater than 0.99 in all cases. Therefore, affine functions were established between the clearness indices by a least-square fitting technique. There is a considerable number of affine functions, and for operational purposes, a limited set of 29 intervals of 1 nm in width, the fine bands FBj, was selected and then used in a linear interpolation process to obtain the clearness indices at each 1 nm without losing accuracy to compute the illuminance. The current approach is empirical with no guarantee that the selected set of FBj is the optimum. It could have been possible to use some mathematical optimization tools. For each FBj in a given KBi, two affine functions have been established: where the parameters a FBj , b FBj , c FBj , and d FBj are given in Table 2. These two sets of affine functions are obtained once for all. The operational method is as follows. For a given set of inputs, libRadtran with the Kato et al. [12] scheme is run to provide the set of thirteen KT KBi and KT B _ KBi from which the set of 29 KT FBj and KT B_FBj is calculated using the affine functions. Then, KT n and KT Bn are estimated at each 1 nm between 380 and 780 nm using spectral linear interpolation of KT FBj and KT B_FBj . Then, this complete set of 1 nm clearness indices is converted into a set of spectral irradiances G n and B Nn , which in turn are multiplied by the standardized CIE action spectrum for human eye, yielding the illuminances L G and L BN .   16 , KB 17 , and KB 18 , where KT n exhibits a nonlinear behavior that cannot be accounted for with a single KT FBj . That is the reason for selecting more than one FB j (magenta crosses). The linear interpolation of KT FBj (in blue line) provides a fairly accurate estimate of KT n . Overall, the graph illustrates the high capability of the spectral resampling technique to reproduce the spectral variation of the clearness index throughout the daylight spectral band.

Results and Discussion
The validation of the presented method was carried out by comparing the estimates to the 1 min measurements of illuminance. The errors were as follows: estimates minus measurements, were synthetized by usual statistical quantities such as the bias (mean error), the root mean square error (RMSE), and their relative values rBias and rRMSE with respect to the average of the measurements. The squared correlation coefficient also known as the coefficient of determination (R 2 ) was also computed. Figure 2 exhibits the 2D histogram between 1 min measured and estimated global illuminances on horizontal surface at Golden. The 2D histogram displays the number of data pairs, also known as a count, between measured and estimated illuminances within a given bin of area 1.5 klx × 1.5 klx. Counts increase from blue to yellow. In other words, the yellower the color, the greater the count; the bluer the color, the lower the count. The cloud of points is well elongated along the 1:1 line with a limited scattering, although the spread of points increases with increasing magnitude of global illuminance. The coefficient of determination is 0.95, denoting that more than 95% of the variability in measurements, expressed by the variance, is reproduced by the presented method. The method reveals a small positive bias and rbias of 1 klx and 1% respectively meaning a small overestimation. The RMSE is 6 klx; i.e., 9% of the average value of measurements. Statistics of errors for each station are reported in Table 3, together with the number of cloudless instants Ndata. With respect of the accuracy of the instrument, the method shows a good performance in assessing global illuminance on a horizontal surface. In addition, one observes that the vast majority of the points drops within the area delimited by the two thin magenta lines representing the uncertainty of the Global Photometric sensor model LI-210 estimated here to be approximately ±8% from the measurements. In order to investigate the sources of errors, the ratio of estimated and measured global illuminances was plotted as a function of the CAMS inputs and the solar zenithal angle (θ s ). The investigation has been done on deviations as well. Figure 3 shows the ratio (top) and deviation (bottom) as function of θ s , daylight albedo, TOC, TWV, and AOD for Golden. The dependence of the errors is graphically shown with the boxplot aiming to understand the distribution of errors and to measure the degree of dispersion and skewness of the errors for a specific interval based on a given CAMS input. The boxplot exhibits both upper and lower quartiles and the median. The results are displayed on the panels referring to one CAMS input. A boxplot is shown for a given range of values.

Validation on Global Illuminance on Horizontal Surface
Overall, there is no obvious influence of the investigated inputs on the errors. An exception is made on AOD where the ratios and deviations tend to become respectively less than 1 and more negative when the AOD increases beyond 0.6. For the ratios and deviations, the magnitudes of boxes in the panel are quite small, denoting a very limited scattering of errors. These magnitudes are fairly identical in the panel of each investigated input denoting a weak influence of inputs on errors except at high AOD. Nevertheless, it is interesting to explore these strong underestimations for which AOD is greater than 0.7 in somewhat more detail.
Most of these cases occurred on 15 April 2017 between 20:08 and 22:00 UTC. Figure 4 illustrates the diurnal variability of measured and CAMS AOD at Golden over a longer period i.e., from 18:00 to 22:00 UTC. The graph provides a closer look on the temporal matching between measured and CAMS AOD. As used in this study, the original 3-h CAMS AOD has been interpolated for providing 1 min CAMS AOD estimates. For the selected cloudless instants (SCI) in this longer period, CAMS AOD (blue dot) is always much greater than the measured AOD (green dot), clearly explaining the observed strong underestimations. In addition, the graph reveals the issue of interpolated CAMS AOD, which might not adequately describe the high temporal variability of measured AOD.  Similarly to Figure 2, Figure 5 exhibits the 2D histogram at Vaulx-en-Velin. The points follow quite well the 1:1 line with a limited spread. The plot shows an overestimation close to a systematic bias. The spread of points is more limited compared to the previous case of Golden. R 2 is very large, up to 0.98, meaning that 98% of the variability contained in the measurements is well captured by the method. The bias is still small and is 3 kx, i.e., +4% of the average of the measurements. The RMSE is small with a value of 5 klx (8%).
At this station, the investigation of the influence of CAMS inputs on the errors yields observations that are similar to those made at Golden. In addition, the dependency of statistical quantities with the month were explored. It was found that the relative bias as well as relative RMSE are greater around the wintertime i.e., between November and February compared to other months. The sources of errors are discussed later in more details.

Validation on Direct Illuminance at Normal Incidence
The validation of the direct illuminance at normal incidence L BN is only carried out at Vaulx-en-Velin. Similarly to Figure 2, Figure 6 shows the 2D histogram between measurements of L BN and estimates for Vaulx-en-Velin while the statistics are reported in Table 4. R 2 is low, with a value of 0.53 denoting that the variability in L BN is only fairly well reflected by the estimates. This small coefficient of determination may be partly explained by the restrained range of variations of the direct illuminance at normal incidence in case of cloudless conditions, and the small values do not strictly imply weaknesses in the method. The bias is 7 klx, i.e., 9% of the average of the measurements. The RMSE and rRMSE are 12 klx and 15%, respectively. Figure 6. Two-dimensional (2D) histogram between measurements and estimates of direct illuminance at normal incidence at Vaulx-en-Velin. The colorbar depicts the number of points in the area of 1.5 klx × 1.5 klx. The two thin magenta lines delimit the area of relative errors of ±5% from the measurements. The investigations on the influence of the inputs on the errors (not shown) do not reveal an obvious dependency except for AOD at 550 nm, where the influence of aerosols is more visible in direct illuminance at normal incidence than in global illuminance. It is clearly seen that both ratio and deviation decrease from greater to less than 1 and from positive to negative respectively with increasing AOD.
The statistical quantities show that the performance of the method is worse during the period from November to February compared to other months. These observations are similar from one year to another.

Error Analysis on Estimated Illuminances
Despite data being available at two stations only, it is worth investigating the sources of errors resulting in an overall overestimation by the presented method: slightly positive bias on the global illuminance at horizontal surface and higher positive bias on direct illuminance at normal incidence at Vaux-en-Velin. We recall that there was no measured illuminance at normal incidence. The values were derived from the measured global and diffuse illuminances at horizontal surface following Equation (1). As a positive bias is observed on validations of global illuminance, it is somehow understandable to see an overestimation on direct illuminance at normal incidence for this site. In general, four sources of errors are possible: either the main CAMS inputs, or ground albedo, or the method itself, or errors in connection with the ground-based measurements.
The accuracy of CAMS AOD and TWV were explored in more details by validating the CAMS AOD and TWV estimates to AOD and TWV measured at the closest AErosol RObotic NETwork stations. The comparison does not reveal any noticeable bias, even at the closest site Aubiere_LAMP located at 140 km. This finding is exactly in line with the worldwide performance of CAMS AOD, including Europe ( [27]). The ground albedo used here reasonably represents the type of surface over year and does not exhibit any seasonality. This is confirmed by a visual analysis of images of the surface throughout the year obtained from the website: https://worldview.earthdata.nasa.gov/, accessed on 24 April 2021). Furthermore, an additional sensitivity analysis performed with libRadtran by modifying the daylight albedo does not remarkably reduce the bias. For aerosol type, "urban" and "continental average" types from OPAC (Optical Properties of Aerosols and Clouds: The Software Package OPAC, [28]) were used at Vaux-en-Velin and Golden respectively. This could affect the estimates and particularly the RMSE, as the local atmosphere may experiment other aerosol types.
An additional numerical validation was carried out by comparing the outputs of the presented method against those from detailed spectral calculations by libRadtran serving as reference to assess the accuracy of the method itself on both global illuminance at horizontal surface and direct illuminance at normal incidence. An independent sample of 10,000 atmospheric states has been randomly built following the marginal distribution variables described earlier in Section 3.2. Overall, squared correlation coefficient of 1.00 was found. In addition, both bias and RMSE (rbias and rRMSE) are less than 0.3 klx (less than 0.3%). These statistics denote the very good level of accuracy of the presented method and therefore indicate that the errors do not originate from the method itself.
It is possible that the illuminance measurements face some issues, especially at Vauxen-Velin. However, measurements seem to be of good quality, since similar errors were found when validating both global irradiance and illuminance on the horizontal surface at this site, and since a positive bias on global irradiance on horizontal surface was also found at the close site of Payerne with the McClear model ( [13]). Therefore, it might be possible that the CAMS products may not capture all atmospheric features at the regional scale particularly around Vaux-en-Velin.

Conclusions
A novel method for deriving the global illuminance on the horizontal surface and its direct component at normal incidence under cloudless conditions has been described here, and its accuracy has been assessed against high quality 1 min measurements. It is a preliminary and necessary step to complete an operational method for deriving illuminances in any sky conditions at any location and any time.
For global illuminance on a horizontal surface, the method exhibits a coefficient of determination greater than 0.95 at both stations, demonstrating that more than 95% of the temporal variability contained in the measurements is very well represented by the method. In general, the method exhibits a small overestimation characterized by a relative positive bias less than 5% of the average of the measurements. The relative RMSE is close to the relative bias, denoting a small standard deviation of errors. This proves a good level of performance of the presented method in assessing the global illuminance with CAMS products as inputs. The level of accuracy from the method is quite similar to the uncertainty of the measurements.
For direct illuminance at normal incidence and only limited to one station, the performance of the method is still good but less when compared to the performance for global illuminance on horizontal surface. Nevertheless, the relative bias is small. The relative RMSE is close to the relative bias denoting a limited scattering of cloud of points. As for global illuminance on a horizontal surface, this also proves the good level of accuracy of the method to estimate direct illuminance at normal incidence.
In their extensive numerical study with two RTMs, Oumbe et al. [29] observed that the influence of the cloud effects on the solar radiation received at ground level may be accurately approximated by the product of the irradiance under cloudless conditions by a cloud modification factor, which is also known as the clear-sky index. This approximation has been used at several occasions for operational computations of the broadband irradiances [30,31] or UV fluxes [32][33][34] for example. Since the proposed method delivers accurate illuminance estimates under cloudless conditions, a clear advantage is that any model assessing the influence due to clouds can be associated to this method to estimate illuminance in all-sky conditions. Data Availability Statement: All used measurements at NREL, Golden, CO and Vaulx-en-Velin were freely available and downloaded through the following websites idmp.entpe.fr/mesfr.htm and midcdmz.nrel.gov/apps/sitehome.pl?site=BMS respectively, accessed on: 1 Decmeber 2020. The McClear atmospheric products are accessible after registration and were downloaded from: http: //www.soda-pro.com/, accessed on 1 December 2020. Ground surface images can be downloaded from the following website: https://worldview.earthdata.nasa.gov/, accessed on 24 April 2021).