Estimation of Land Surface Incident and Net Shortwave Radiation from Visible Infrared Imaging Radiometer Suite (VIIRS) Using an Optimization Method

: Incident surface shortwave radiation (ISR) is a key parameter in Earth’s surface radiation budget. Many reanalysis and satellite-based ISR products have been developed, but they often have insu ﬃ cient accuracy and resolution for many applications. In this study, we extended our optimization method developed earlier for the MODIS data with several major improvements for estimating instantaneous and daily ISR and net shortwave radiation (NSR) from Visible Infrared Imaging Radiometer Suite observations (VIIRS), including (1) an integrated framework that combines look-up table and parameter optimization; (2) enabling the calculation of net shortwave radiation (NSR) as well as daily values; and (3) extensive global validation. We validated the estimated ISR values using measurements at seven Surface Radiation Budget Network (SURFRAD) sites and 33 Baseline Surface Radiation Network (BSRN) sites during 2013. The root mean square errors (RMSE) over SURFRAD sites for instantaneous ISR and NSR were 83.76 W / m 2 and 66.80 W / m 2 , respectively. The corresponding daily RMSE values were 27.78 W / m 2 and 23.51 W / m 2 . The RMSE at BSRN sites was 105.87 W / m 2 for instantaneous ISR and 32.76 W / m 2 for daily ISR. The accuracy is similar to the estimation from MODIS data at SURFRAD sites but the computational e ﬃ ciency has improved by approximately 50%. We also produced global maps that demonstrate the potential of this algorithms to generate global ISR and NSR products from the VIIRS data.


Introduction
Incident surface shortwave radiation (ISR) is a critical parameter in Earth's surface radiation budget. It determines the incoming energy source for the Earth's surface and drives energy, ecological, and hydrological dynamics [1][2][3][4][5][6]. Surface net shortwave radiation (NSR) largely determines the total net irradiance at the Earth's surface, which regulates most biological and physical processes at the surface [2,3]. Because of its importance, many regional and global observation networks have been established, such as the Surface Radiation Budget Network (SURFRAD) [7]. However, because of the limited spatial coverage and representativeness, site-based radiation observations have drawbacks when used in many regional and global applications. Reanalysis products are another source of radiation information, and nearly all reanalysis data include radiation data at the Earth's surface, such as The SURFRAD ISR was measured with Spectrolab and Eppley Pyranometer, with a documented error of 5%. BSRN sites used CG4 (Kipp & Zonen) or PIR (Eppley) pyrgeometers for ISR measurements and presented an error of about 5 W/m 2 . We averaged the observations of thirty minutes around the time of the satellite overpass for instantaneous ISR validation to enhance the spatial representativeness of each site [53] and averaged observations of each day for the daily ISR validation. Tables 1 and 2 show the SURFRAD and 33 BSRN sites used in the validation. When validating against BSRN measurements, we divided all sites into seven different areas, namely North America, Europe, South America, Oceania, Asia, Africa, and Greenland. Figure 1 shows the spatial distribution of the BSRN sites.

Algorithm Framework
The algorithm firstly optimizes surface and atmospheric parameters from TOA reflectance and subsequently estimates ISR. Then we used the Wang and Liang method [49] to calculate daily integrated ISR from VIIRS observations. We validated the results against field measurements from seven SURFRAD sites and 33 BSRN sites globally. There are several major improvements over the

Algorithm Framework
The algorithm firstly optimizes surface and atmospheric parameters from TOA reflectance and subsequently estimates ISR. Then we used the Wang and Liang method [49] to calculate daily integrated ISR from VIIRS observations. We validated the results against field measurements from seven SURFRAD sites and 33 BSRN sites globally. There are several major improvements over the earlier version of the algorithm, including (1) integrated look-up table and optimization framework, (2) calculation of net shortwave radiation as well, (3) adding the estimation of daily values, and (4) extensive global validation. We produced global daily integrated ISR maps on 1 January, April, July, and October 2018 as examples. Figure 2 illustrates the framework of the algorithm. An optimization method has been used [52,54] to estimate surface reflectance and broadband albedo. In our previous research, we developed a similar approach for incident shortwave radiation estimation from MODIS data by revising the cost function considering both satellite observations and optional constraints, including aerosol optical depth (AOD), cloud optical depth (COD), surface reflectance products, and albedo climatology [51]. In this study, we adapted the algorithm for the estimation of ISR from VIIRS data by revising the band configuration and the spectrum of radiation transfer simulation. An assumption is made here that the surface reflectance is invariant during a short period (8 days in this research). Under cloudy sky conditions, we used the surface reflectance calculated by the nearest previous clear observations as surface input. The COD can then be optimized using radiative transfer models.  In this study, atmospheric optical parameters such as spherical albedo, atmospheric downward/upward transmittance, and path reflectance are required to calculate the cost function. We used libRadtran [55] to simulate the parameters. We used surface and atmospheric parameters to implement a forward simulation of TOA reflectance using the radiative transfer model. The cost function is based on the difference between simulated and observed TOA reflectance (Equations (1)-(3)): Here, and refer to satellite-observed TOA reflectance and forward-simulated TOA reflectance from the radiative transfer model. and are the calculated albedo and the albedo climatology, respectively. and are the parameters to be optimized under clear and cloudy sky cases, respectively. O and B are the error matrices for the TOA reflectance and the climatology, respectively. The albedo climatology covariance matrix B is determined from the multiyear satellite albedo climatology uncertainty, while the TOA reflectance covariance matrix O is calculated from the narrowband albedo's contribution to the broadband albedo and the spectral reflectance magnitude. Jc is a punishment component and set to a 100 if the TOA reflectance is an invalid value (negative or greater than one). The shuffled complex evolution (SCE) method [56] was used to search for the optimum.
When the surface and atmospheric parameters reach an optimum, ISR ( ) and NSR ( ) are then calculated using a radiative transfer model according to Equations (4)-(6): In this study, atmospheric optical parameters such as spherical albedo, atmospheric downward/ upward transmittance, and path reflectance are required to calculate the cost function. We used libRadtran [55] to simulate the parameters.
We used surface and atmospheric parameters to implement a forward simulation of TOA reflectance using the radiative transfer model. The cost function is based on the difference between simulated and observed TOA reflectance (Equations (1)-(3)): Here, R obs and R est refer to satellite-observed TOA reflectance and forward-simulated TOA reflectance from the radiative transfer model. A and A clm are the calculated albedo and the albedo climatology, respectively. X clear and X cloudy are the parameters to be optimized under clear and cloudy sky cases, respectively. O and B are the error matrices for the TOA reflectance and the climatology, respectively. The albedo climatology covariance matrix B is determined from the multiyear satellite albedo climatology uncertainty, while the TOA reflectance covariance matrix O is calculated from the narrowband albedo's contribution to the broadband albedo and the spectral reflectance magnitude. J c is a punishment component and set to a 100 if the TOA reflectance is an invalid value (negative or greater than one). The shuffled complex evolution (SCE) method [56] was used to search for the optimum.
When the surface and atmospheric parameters reach an optimum, ISR F(µ 0 ) and NSR F n (µ 0 ) are then calculated using a radiative transfer model according to Equations (4)-(6): Remote Sens. 2020, 12, 4153 6 of 23 Here, F 0 (µ 0 ) is the radiation without any contribution from the surface, while F dir (µ 0 ) and F di f (µ 0 ) denote the direct and diffuse parts, respectively. r s is the surface reflectance, ρ is the spherical albedo, µ 0 is the cosine of the solar zenith angle, E 0 is the extraterrestrial solar radiation, and γ(µ 0 ) is the total transmittance. For each combination of geometry and optical depth, the F 0 (µ 0 ), ρ and µ 0 E 0 γ(µ 0 ) were pre-calculated by radiative transfer simulation in the ISR spectrum range and stored in a look-up table. α denotes the land surface broadband albedo, F(µ 0 ) and F n (µ 0 ) denote surface ISR and NSR, respectively.
Equations (4)-(6) are based on instantaneous observations. The daily ISR were then estimated using the LUT-based algorithm by Wang et al. [57]. We calculated ISR every 30 min according to the interpolation algorithm. Then, we calculated the average ISR for each day.

Atmospheric Look-Up Tables Improvements
Atmospheric optical parameters such as spherical albedo, atmospheric downward/upward transmittance, and path reflectance are required to implement a forward simulation. To make the algorithm more efficient, all the parameters were pre-calculated in representative geometries and atmospheric conditions (AOD, cloud optical depth [COD] cloud effective radius [CER]). LibRadtran [55] software was used for the generation of the LUT. We used the continental-clean model to estimate ISR. For each specific solar/viewing geometry and atmospheric parameter (AOD at 550 nm for clear-sky conditions, COD and CER for cloudysky conditions), radiative transfer simulations generated path reflectance, upward/downward transmittances, and spherical albedo for each of the seven VIIRS bands. We used actual site elevation to estimate ISR at the SURFRAD and BSRN sites, and the Global 30 Arc-Second Elevation (GTOPO30) [58,59] for the global ISR map. With the atmospheric LUT, we calculated the surface broadband albedo and atmospheric index (AOD for clear-sky conditions, COD and CER for cloudy-sky conditions) from the optimization process. For cloudy-sky cases, each observation was optimized using both the water and ice cloud LUTs, of which the one with a smaller cost function result was chosen as the result. ISR could then be calculated under certain geometries using the surface radiation LUT. In this paper, we calculated the ISR for the spectral range of 280-2800 nm to match the field measurements.

Improvement of the Optimization Framework
In the previous MODIS algorithm, most of the unknown variables in the optimization process were surface parameters while only one was for the atmospheric condition (AOD/COD). We prioritized the unknown variables in this research to increase more atmospheric parameters.
The influence of ISR is much more associated with the atmospheric conditions, especially the cloud conditions. Surface parameters also count for multiple scattering which also contributes to the ISR. In this research, we tried to minimize the unknown surface variables to increase both the efficiency and accuracy of the model.
We used principal component analysis (PCA) to analyze the surface spectral reflectance of the nine VIIRS bands using the data from all the site measurements extracted from reflectance products. Results show that the first two components explain more than 98 percent of the variations. In the updated framework, we only use two free variables for the surface condition, which reduce the total unknown variables to three in clear-sky conditions and four in cloudy-sky conditions.
The PCA could increase the efficiency with low errors introduced. The ISR is mainly determined by the path flux (the first component in Equation (4)) and the multiple scattering part (the second component in Equation (4)). The albedo climatology ( (1) constrains the surface reflectance and albedo error. Over dark surfaces, such as forests, the multiple scattering part is minimal in determining the overall ISR due to the low surface albedo. On the other hand, over bright surfaces such as deserts, the principal components could explain much of the overall variations due to relatively simple surface conditions. Figure 3a illustrates the radiative transfer simulated maximum error over dark and bright surfaces with the first two components at nadir (solar zenith angle = 0). The maximum possible differences at nadir are calculated under different visibility conditions. The error is more massive when the cloud is dense (optical depth greater than 80), and when dense aerosol (optical depth between 1~2, usually results from a mixture of aerosol and cloud). In all cases, the induced error is less than 1.6 W/m 2 . Figure 3b,c show the original spectral reflectance and the PCs from multiple sites over bright surfaces. The PC1 explains the variation for spectral bands at a shorter wavelength (<1000 nm), while the PC2 explains the other bands.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 23 multiple scattering part is minimal in determining the overall ISR due to the low surface albedo. On the other hand, over bright surfaces such as deserts, the principal components could explain much of the overall variations due to relatively simple surface conditions. Figure 3a illustrates the radiative transfer simulated maximum error over dark and bright surfaces with the first two components at nadir (solar zenith angle = 0). The maximum possible differences at nadir are calculated under different visibility conditions. The error is more massive when the cloud is dense (optical depth greater than 80), and when dense aerosol (optical depth between 1 ~ 2, usually results from a mixture of aerosol and cloud). In all cases, the induced error is less than 1.6 W/m 2 . Figure 3b,c show the original spectral reflectance and the PCs from multiple sites over bright surfaces. The PC1 explains the variation for spectral bands at a shorter wavelength (<1000 nm), while the PC2 explains the other bands.

Validation Results of Instantaneous ISR
The validation results of instantaneous ISR are shown in Figure 4 and Table 3. The ISR RMSE ranges from 75.32 W/m 2 to 94.66 W/m 2 at the different sites while the bias ranges from −29.36 W/m 2 to 9.21 W/m 2 . The overall RMSE and bias are 83.76 W/m 2 and −3.49 W/m 2 , respectively. RMSE Origin denotes the results from the original optimization framework without the improvements. The updated algorithm provides accuracy improvement at most of the SURFRAD sites and the overall result.

Validation results of instantaneous ISR
The validation results of instantaneous ISR are shown in Figure 4 and Table 3. The ISR RMSE ranges from 75.32 W/m 2 to 94.66 W/m 2 at the different sites while the bias ranges from −29.36 W/m 2 to 9.21 W/m 2 . The overall RMSE and bias are 83.76 W/m 2 and −3.49 W/m 2 , respectively. RMSEOrigin denotes the results from the original optimization framework without the improvements. The updated algorithm provides accuracy improvement at most of the SURFRAD sites and the overall result. .

Validation Results of Instantaneous NSR
The validation results of instantaneous NSR are shown in Figure 5 and

Validation Results of Instantaneous NSR
The validation results of instantaneous NSR are shown in Figure 5 and

Validation Results of Daily ISR
The validation results of daily ISR are shown in Figure 6 and

Validation Results of Daily ISR
The validation results of daily ISR are shown in Figure 6 and

Validation Results of Daily ISR
The validation results of daily ISR are shown in Figure 6 and Table 5. The ISR RMSE ranges from 24.40 W/m 2 to 34.11 W/m 2 at the different sites while the bias ranges from −8.60 to 10.59 W/m 2 . The overall daily RMSE and bias are 27.78 W/m 2 and −0.16 W/m 2 , respectively. . Figure 6. Validation results of daily ISR at SURFRAD sites Figure 6. Validation results of daily ISR at SURFRAD sites.

Validation Results of Daily NSR
The validation results of daily NSR are shown in Figure 7 and

Validation Results of Daily NSR
The validation results of daily NSR are shown in Figure 7 and

Validation Results of Instantaneous ISR
The validation results of instantaneous ISR are shown in Figure 8   respectively. The RMSE at BSRN sites is generally higher than that at SURFRAD sites even they are close to each other.     There is only one site located in this Greenland area. The RMSE and bias are 74.28 W/m 2 and −40.71 W/m 2 , respectively.
The validation results of instantaneous ISR are shown in Figure 8, Tables 7 and 8. The ISR RMSE in each area ranges from 83.36 W/m 2 to 129 W/m 2 at the different sites while the bias ranges from −40.71 W/m 2 to 9.48 W/m 2 . South America sites have the largest overall RMSE as there is more cloud in this area. Greenland site has the smallest RMSE due to the high latitude and the low absolute value of ISR. The overall RMSE and bias are 105.87 W/m 2 and −3.79 W/m 2 , respectively.

Validation Results of Daily ISR
The validation results of daily ISR are shown in Figure 9, Tables Figure 9, Tables 9 and 10. The ISR RMSE in each area ranges from 28.49 W/m 2 to 39.75 W/m 2 at the different sites while the bias ranges from −27.06 W/m 2 to 10.24 W/m 2 . The overall RMSE and bias are 32.75 W/m 2 and 0.29 W/m 2 , respectively.  Figure 11 and 12 show the spatial patterns of the validation results for instantaneous and daily ISR, respectively. Red circles denote sites with positive biases while blue circles denote sites with negative biases. The larger the diameters of the circles, the larger the RMSEs are and vice versa. Overall, sites located in lower latitudes and located in islands have larger RMSEs compared to other sites. The algorithm tends to underestimate ISR over North America and overestimate ISR over Europe and Asia. Figure 13. shows the VIIRS daily ISR map on 1 Jan. 2018.   Figures 11 and 12 show the spatial patterns of the validation results for instantaneous and daily ISR, respectively. Red circles denote sites with positive biases while blue circles denote sites with negative biases. The larger the diameters of the circles, the larger the RMSEs are and vice versa. Overall, sites located in lower latitudes and located in islands have larger RMSEs compared to other sites. The algorithm tends to underestimate ISR over North America and overestimate ISR over Europe and Asia. Figure 13. shows the VIIRS daily ISR map on 1 Jan. 2018.  Figure 11 and 12 show the spatial patterns of the validation results for instantaneous and daily ISR, respectively. Red circles denote sites with positive biases while blue circles denote sites with negative biases. The larger the diameters of the circles, the larger the RMSEs are and vice versa. Overall, sites located in lower latitudes and located in islands have larger RMSEs compared to other sites. The algorithm tends to underestimate ISR over North America and overestimate ISR over Europe and Asia. Figure 13. shows the VIIRS daily ISR map on 1 Jan. 2018.

Discussion
The PCA performed to the surface reflectance could reduce the optimization framework's unknown variables, providing a significant improvement in efficiency and slight inaccuracy. Many previous studies used surface reflectance from single or multiple spectral bands to estimate surface ISR. Traditional single-band LUT algorithms tend to assume a single condition that the band used (often the blue band) itself can represent the surface condition and determine the atmosphere's clearness. On the other hand, the optimization algorithm using all shortwave spectral bands is timeconsuming. The PCA can extract critical information for the surface condition. Apart from improving the efficiency, the approach can also improve the accuracy slightly, as fewer free variables bring less uncertainty to the optimization process.
The most significant improvement is at the Table Mountain site at Boulder (CO, USA). The cloud condition in the mountainous area here is the most complicated among the SURFRAD sites. Overall, sites with dark surfaces have slight improvement while not at the sites with bright surfaces, such as Desert Rock. Further improvements are required to apply this algorithm over the bright surfaces.

Discussion
The PCA performed to the surface reflectance could reduce the optimization framework's unknown variables, providing a significant improvement in efficiency and slight inaccuracy. Many previous studies used surface reflectance from single or multiple spectral bands to estimate surface ISR. Traditional single-band LUT algorithms tend to assume a single condition that the band used (often the blue band) itself can represent the surface condition and determine the atmosphere's clearness. On the other hand, the optimization algorithm using all shortwave spectral bands is timeconsuming. The PCA can extract critical information for the surface condition. Apart from improving the efficiency, the approach can also improve the accuracy slightly, as fewer free variables bring less uncertainty to the optimization process.
The most significant improvement is at the Table Mountain site at Boulder (CO, USA). The cloud condition in the mountainous area here is the most complicated among the SURFRAD sites. Overall, sites with dark surfaces have slight improvement while not at the sites with bright surfaces, such as Desert Rock. Further improvements are required to apply this algorithm over the bright surfaces.

Discussion
The PCA performed to the surface reflectance could reduce the optimization framework's unknown variables, providing a significant improvement in efficiency and slight inaccuracy. Many previous studies used surface reflectance from single or multiple spectral bands to estimate surface ISR. Traditional single-band LUT algorithms tend to assume a single condition that the band used (often the blue band) itself can represent the surface condition and determine the atmosphere's clearness. On the other hand, the optimization algorithm using all shortwave spectral bands is time-consuming. The PCA can extract critical information for the surface condition. Apart from improving the efficiency, the approach can also improve the accuracy slightly, as fewer free variables bring less uncertainty to the optimization process.
The most significant improvement is at the Table Mountain site at Boulder (CO, USA). The cloud condition in the mountainous area here is the most complicated among the SURFRAD sites. Overall, sites with dark surfaces have slight improvement while not at the sites with bright surfaces, such as Desert Rock. Further improvements are required to apply this algorithm over the bright surfaces.
Boundaries between different orbits over southern Africa and Australia in the global daily map indicate that estimating daily averaged ISR from instantaneous observations is difficult. Further work is required in generating global operational products.