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

Abstract

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 insufficient 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² and 66.80 W/m², respectively. The corresponding daily RMSE values were 27.78 W/m² and 23.51 W/m². The RMSE at BSRN sites was 105.87 W/m² for instantaneous ISR and 32.76 W/m² for daily ISR. The accuracy is similar to the estimation from MODIS data at SURFRAD sites but the computational efficiency 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.

1. 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 Japanese 55-year Reanalysis (JRA-55) [8], the ERA-5 [9], the Modern-Era Retrospective analysis for Research and Applications (MERRA) [10], the National Centers for Environmental Prediction (NCEP) [11,12], and the Climate Forecast System Reanalysis (CFSR) [13].

Many satellite-derived radiation products have been published and are widely used, such as the Clouds and the Earth’s Radiant Energy System (CERES) [14], the International Satellite Cloud Climatology Project (ISCCP) [15], the Global Energy and Water Exchanges (GEWEX) [16], Moderate Resolution Imaging Spectroradiometer (MODIS) MCD18 [17], the Global LAnd Surface Satellite (GLASS) [18], and the Satellite Application Facility on Climate Monitoring (CM SAF) [19]. However, most of the existing products have a relatively coarse spatial resolution, usually one or one-half of a degree, which limits the capacity to quantify regional and local changes. The World Meteorological Organization (WMO) requires a spatial resolution of one kilometer in specific applications such as agricultural meteorology. Moreover, the accuracy of the existing products is still far from the uncertainty requirements. Many studies have assessed the widely used satellite products, finding RMSEs are usually 80–150 W/m² at hourly/3-hourly temporal scales [20,21,22]. WMO requires a daily uncertainty (Root Mean Square Error, RMSE) of less than 20 W/m² for the Numerical Weather Prediction (NWP) application (https://www.wmo-sat.info/oscar/variables/view/50), while all four widely used global radiation products have greater uncertainty than 20 W/m² from extensive validation [21,23].

Satellite-derived ISR products usually are produced using four groups of algorithms: parameterization methods [24,25,26,27,28,29], look-up-table (LUT)-based methods [30,31,32,33], direct radiative transfer calculation methods [14,34], and machine-learning methods [35,36,37,38,39,40,41,42,43]. Most of the NSR estimation algorithms link top of atmosphere (TOA) reflectance directly with NSR using radiative transfer simulations [44,45,46,47], while some algorithms calculate NSR using downward and upward components [48] or ISR and land surface albedo. Many of the ISR and NSR algorithms inherit uncertainties from the input data processing and model which still cannot meet the accuracy requirements.

Most algorithms provide instantaneous ISR results, which are calculated at the satellite overpass time. In energy balance studies, daily integrated ISR results are more frequently required. Kim and Liang [46] calculated daily integrated NSR using a sinusoidal model from instantaneous values. Wang et al. [47] extended the model to Landsat data with the assumption of constant atmospheric conditions during one day. Wang and Liang [49] then refined the method for NASA’s Earth Observing System’s MODIS data without the requirement for additional atmospheric water vapor data input. Mateos et al. [50] validated daily UV surface irradiance from the Ozone Monitoring Instrument (OMI) at 14 ground-based stations.

We recently developed an optimization algorithm for estimating instantaneous ISR from MODIS data [51], which is based on an earlier version for estimating land surface albedo in clear-sky cases [52]. The previous MODIS algorithm requires multiple observations as input for the optimization, which makes the algorithm very time-consuming and hardly practical for global operational production. The previous version of the cloud look-up table only considered water cloud, which introduces errors where ice cloud is present. To improve the efficiency and accuracy, here we modified the algorithm and extended it to estimate instantaneous and daily ISR and NSR from the Visible Infrared Imaging Radiometer Suite (VIIRS) data and conducted a comprehensive global validation. VIIRS is onboard Suomi National Polar-orbiting Partnership (Suomi NPP) and the Joint Polar Satellite System (JPSS) polar satellites, and there is no radiation product publicly available currently.

Section 2 describes the data and methodology used in this research. Section 3 and Section 4 analyzes the validation results from, and Section 5 provides a summary.

2. Data and Methodology

2.1. Datasets Used in This Study

The VIIRS instrument provides observation continuity with MODIS. The VIIRS team provides a suite of operational products, termed sensor data records (SDR). In this study, we used VIIRS SDR archive sets from National Oceanic and Atmospheric Administration (NOAA) Comprehensive Large Array-data Stewardship System (CLASS) to estimate surface ISR in 2013. The VIIRS data from bands 1−11 were used in this study except for the two absorption bands (band 6 and 9) at the spatial resolution of 750-m.

2.2. Validation Sites

Ground measurements from seven SURFRAD and BSRN sites during 2013 were used in this study to validate ISR estimation. The recording cycle of SURFRAD and BSRN are both one minute. 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². 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.

Table 1. SURFRAD sites used in the validation.

Site Name	Latitude	Longitude	Elevation (m)
Fort Peck	48.30798	−105.10177	634
Sioux Falls	43.73431	−96.62334	473
Penn State	40.72033	−77.931	376
Bondville	40.05155	−88.37325	230
Boulder	40.12557	−105.23775	1689
Desert Rock	36.6232	−116.01962	1007
Goodwin Creek	34.2547	−89.8729	98

and

Table 2. BSRN sites used in the validation.

Area	Site Name	Latitude	Longitude	Elevation (m)
North America	BIL	36.605	−97.516	317
	BON	40.0667	−88.3667	213
	BOS	40.125	−105.237	1689
	BOU	40.05	−105.007	1577
	CLH	36.905	−75.713	37
	DRA	36.626	−116.018	1007
	E13	36.605	−97.485	318
	FPE	48.3167	−105.1	634
	GCR	34.2547	−89.8729	98
	PSU	40.72	−77.9333	376
	SXF	43.73	−96.62	473
Europe	CAB	51.9711	4.9267	0
	CAM	50.2167	−5.3167	88
	CAR	44.083	5.059	100
	CNR	42.816	−1.601	471
	LER	60.1389	−1.1847	80
	LIN	52.21	14.122	2862
	PAL	48.713	2.208	156
South America	BRB	−15.601	−47.713	1023
	FLO	−27.6047	−48.5227	11
	PTR	−9.068	−40.319	387
	SMS	−29.4428	−53.8231	489
Oceania	ASP	−23.798	133.888	547
	DAR	−12.425	130.891	30
	DWN	−12.424	130.8925	32
	LAU	−45.045	169.689	350
	MAN	−2.058	147.425	6
Asia	FUA	33.5822	130.3764	3
	ISH	24.3367	124.1644	5.7
	XIA	39.754	116.962	32
Africa	GOB	−23.5614	15.042	407
	TAM	22.7903	5.5292	1385
Greenland	ALE	82.49	−62.42	127

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.

2.3. 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)):

$$J\left(X\right)=\left(R^{est}\left(X\right)-R^{obs}\left(X\right)\right)O^{-1}\left(R^{est}\left(X\right)-R^{obs}\left(X\right)\right)+J_c+\left[\left(A(X\right)-A^{clm}\right)B^{-1}\left(A(X\right)-A^{clm})$$

(1)

$$X_{clear}= {\left[BR{F}_1,BR{F}_2,\dots ,BR{F}_{NB},AO{D}_1,AO{D}_2,\dots ,AO{D}_{NO}\right]}^T$$

(2)

$$X_{cloudy}={\left[CO{D}_1,CO{D}_2,\dots ,CO{D}_{NO},CE{R}_1,CE{R}_2,\dots ,CE{R}_{NO}\right]}^T$$

(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\left({\mu }_0\right)$ and NSR $F_n\left({\mu }_0\right)$ are then calculated using a radiative transfer model according to Equations (4)–(6):

$$F\left({\mu }_0\right)= F_0\left({\mu }_0\right)+\frac{r_s\overline{\rho }}{1-r_s\overline{\rho }}{\mu }_0{E}_0\gamma \left({\mu }_0\right)$$

(4)

$$F_0\left({\mu }_0\right)=F_{dir}\left({\mu }_0\right)+F_{dif}\left({\mu }_0\right)$$

(5)

$$F_n\left({\mu }_0\right)=F\left({\mu }_0\right)\left(1-\alpha \right)$$

(6)

Here, $F_0\left({\mu }_0\right)$ is the radiation without any contribution from the surface, while $F_{dir}\left({\mu }_0\right)$ and $F_{dif}\left({\mu }_0\right)$ denote the direct and diffuse parts, respectively. $r_s$ is the surface reflectance, $\overline{\rho }$ is the spherical albedo, ${\mu }_0$ is the cosine of the solar zenith angle, $E_0$ is the extraterrestrial solar radiation, and $\gamma \left({\mu }_0\right)$ is the total transmittance. For each combination of geometry and optical depth, the $F_0\left({\mu }_0\right)$, $\overline{\rho }$, and ${\mu }_0{E}_0\gamma \left({\mu }_0\right)$ were pre-calculated by radiative transfer simulation in the ISR spectrum range and stored in a look-up table. $\alpha$ denotes the land surface broadband albedo, $F\left({\mu }_0\right)$ and $F_n\left({\mu }_0\right)$ 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.

2.4. 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. The following values were used as entries in the radiative transfer simulations: solar zenith angle (0°–80°, at 10° intervals), viewing zenith angle (0°–80°, at 20° intervals), relative azimuth angle (0°–180°, at 30° intervals), AOD at 550 nm (0.01, 0.025, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0), surface elevation (0, 1, 2, 3 4, 5 km), and water vapor (0, 15, 30, 45, 60, 75, 90, 105 mm). COD (1, 3, 5, 10, 20, 40, 60, 80) and CER (3, 6, 9, 12 um) are set for water cloud while CER (20, 30, 40, 50, 60, 70, 80 um) are set for ice cloud.

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 cloudy-sky 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.

2.5. 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 $\left[\left(A(X\right)-A^{clm}\right)B^{-1}\left(A(X\right)-A^{clm})$] in Equation (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². 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.

3. Results over SURFRAD Sites

3.1. Validation Results of Instantaneous ISR

The validation results of instantaneous ISR are shown in Figure 4 and

Table 3. Validation results of instantaneous ISR at SURFRAD sites.

Site Name	R2	Bias	RMSE	RMSE Origin
	Instantaneous Rs Unit: W/m2
Fort Peck	0.89	4.24	80.62	82.17
Sioux Falls	0.90	2.80	82.00	84.28
Penn State	0.88	9.21	94.66	98.03
Bondville	0.93	3.61	77.39	79.31
Boulder	0.92	−24.17	80.76	92.57
Desert Rock	0.92	−29.36	75.32	73.85
Goodwin Creek	0.89	11.88	93.01	96.18
All	0.91	−3.49	83.76	87.30

. The ISR RMSE ranges from 75.32 W/m² to 94.66 W/m² at the different sites while the bias ranges from −29.36 W/m² to 9.21 W/m². The overall RMSE and bias are 83.76 W/m² and −3.49 W/m², 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.

3.2. Validation Results of Instantaneous NSR

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

Table 4. Validation results of instantaneous NSR at SURFRAD sites.

Site Name	R2	RMSE	Bias
Instantaneous Net Rs Unit: W/m2
Fort Peck	0.94	58.63	2.47
Sioux Falls	0.93	65.5	4.83
Penn State	0.92	73.83	10.36
Bondville	0.96	56.8	4.71
Boulder	0.93	68.57	−17.18
Desert Rock	0.92	62.49	−20.48
Goodwin Creek	0.91	80.87	14.6
All	0.93	66.8	0.62

. The RMSE ranges from 56.8 W/m² to 80.87 W/m² at the different sites while the bias ranges from −20.48 W/m² to 14.6 W/m². The overall RMSE and bias are 66.80 W/m² and 0.62 W/m², respectively.

3.3. Validation Results of Daily ISR

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

Table 5. Validation results of daily ISR at SURFRAD sites.

Site Name	R2	RMSE	Bias
Daily Rs Unit: W/m2
Fort Peck	0.97	24.4	−0.06
Sioux Falls	0.97	24.7	−0.41
Penn State	0.95	31.13	−5.41
Bondville	0.97	26.5	−4.55
Boulder	0.94	31.37	7.49
Desert Rock	0.98	20.64	10.59
Goodwin Creek	0.95	34.11	−8.6
All	0.96	27.78	−0.16

. The ISR RMSE ranges from 24.40 W/m² to 34.11 W/m² at the different sites while the bias ranges from −8.60 to 10.59 W/m². The overall daily RMSE and bias are 27.78 W/m² and −0.16 W/m², respectively.

3.4. Validation Results of Daily NSR

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

Table 6. Validation results of daily NSR at SURFRAD sites.

Site Name	R2	RMSE	Bias
Daily Net Rs Unit: W/m2
Fort Peck	0.94	23.03	4.2
Sioux Falls	0.93	22.82	3.65
Penn State	0.93	22.14	2.77
Bondville	0.93	21.08	2.83
Boulder	0.86	30.31	1.49
Desert Rock	0.94	18.07	−0.99
Goodwin Creek	0.9	25.5	4.51
All	0.93	23.51	2.66

. The RMSE ranges from 18.07 W/m² to 30.31 W/m² at the different sites while the bias ranges from –0.99 W/m² to 4.51 W/m². The overall RMSE and bias are 23.51 W/m² and 2.66 W/m², respectively.

4. Results over BSRN sites

4.1. Validation Results of Instantaneous ISR

The validation results of instantaneous ISR are shown in Figure 8,

Table 7. Validation results of instantaneous ISR at all BSRN sites.

Area	Site Name	RMSE (W/m2)	Bias (W/m2)
North America	BIL	105.66	−6.08
	BON	114.55	−13.23
	BOS	127.65	−26.89
	BOU	112.43	−14.37
	CLH	113.15	−10.31
	DRA	95.41	−24.46
	E13	111.5	−7.8
	FPE	109.86	−19.06
	GCR	120.5	−1.55
	PSU	110.28	−11.96
	SXF	112.64	−16.17
Europe	CAB	93.37	12.49
	CAM	90.77	5.46
	CAR	87.86	−5.5
	CNR	102.18	13.52
	LER	86.72	15.41
	LIN	104.91	15.67
	PAL	106.09	6.54
South America	BRB	148.95	−4.62
	FLO	156.99	13.25
	PTR	111.01	26.04
	SMS	115.36	−10.06
Oceania	ASP	107.41	−34.06
	DAR	123.81	17.92
	DWN	121.35	30.29
	LAU	89.48	−35.15
	MAN	147.77	47.89
Asia	FUA	115.18	30.12
	ISH	132.63	29.55
	XIA	109.98	43.23
Africa	GOB	60.2	−9.8
	TAM	109.39	−53.23
Greenland	ALE	74.28	−40.71
All		105.87	−3.79

and

Table 8. Validation results of instantaneous ISR in different areas.

Area	RMSE (W/m2)	Bias (W/m2)
North America	112.39	−14.05
Europe	95.9	9.48
South America	129	4.29
Oceania	115.16	−1.49
Asia	119.48	34.09
Africa	83.36	−27.11
Greenland	74.28	−40.71
All	105.87	−3.79

. The ISR RMSE in North America ranges from 95.41 W/m² to 127.65 W/m² at the different sites while the bias ranges from −26.89 W/m² to −1.55 W/m². The overall RMSE and bias are 112.39 W/m² and −14.05 W/m², respectively. The RMSE at BSRN sites is generally higher than that at SURFRAD sites even they are close to each other.

The ISR RMSE in Europe ranges from 86.72 W/m² to 106.09 W/m² at the different sites while the bias ranges from −26.89 W/m² to −1.55 W/m². The overall RMSE and bias are 95.90 W/m² and 9.48 W/m², respectively.

The ISR RMSE in South America ranges from 111.01 W/m² to 156.99 W/m² at the different sites while the bias ranges from −10.06 W/m² to 26.04 W/m². The overall RMSE and bias are 129.00 W/m² and 4.29 W/m², respectively.

The ISR RMSE in Oceania ranges from 89.48 W/m² to 147.77 W/m² at the different sites while the bias ranges from −35.15 W/m² to 47.89 W/m². The overall RMSE and bias are 115.16 W/m² and −1.49 W/m², respectively.

The ISR RMSE in Asia ranges from 109.98 W/m² to 132.63 W/m² at the different sites while the bias ranges from 29.55 W/m² to 43.23 W/m². The overall RMSE and bias are 119.48 W/m² and 34.09 W/m², respectively.

The ISR RMSE in Africa ranges from 60.2 W/m² to 109.39 W/m² at the different sites while the bias ranges from −53.23 W/m² to −9.8 W/m². The overall RMSE and bias are 83.36 W/m² and −27.11 W/m², respectively.

There is only one site located in this Greenland area. The RMSE and bias are 74.28 W/m² and −40.71 W/m², respectively.

The validation results of instantaneous ISR are shown in Figure 8, Table 7 and Table 8. The ISR RMSE in each area ranges from 83.36 W/m² to 129 W/m² at the different sites while the bias ranges from −40.71 W/m² to 9.48 W/m². 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² and −3.79 W/m², respectively.

4.2. Validation Results of Daily ISR

The validation results of daily ISR are shown in Figure 9,

Table 9. Validation results of daily ISR at all BSRN sites.

Area	Site Name	RMSE (W/m2)	Bias (W/m2)
North America	BIL	32.2	1.74
	BON	35.68	−2.35
	BOS	39.96	−7.59
	BOU	33.88	−3.87
	CLH	33.88	−4.48
	DRA	26.41	−8.89
	E13	33.24	−0.45
	FPE	34.95	−5.05
	GCR	38.47	2.2
	PSU	33.13	−2.61
	SXF	36.03	−3.76
Europe	CAB	27.93	6.12
	CAM	26.95	4.46
	CAR	27.81	−0.35
	CNR	30.96	5.02
	LER	25.4	6.97
	LIN	32.69	6.86
	PAL	33.07	4.32
South America	BRB	51.55	0.4
	FLO	55.09	4.68
	PTR	37.4	9.4
	SMS	34.13	−3.05
Oceania	ASP	25.66	−7.44
	DAR	30.33	3.61
	DWN	29.41	5.38
	LAU	17.41	−6.49
	MAN	34.63	11.02
Asia	FUA	27.8	9.13
	ISH	34.04	7.78
	XIA	30.24	14.53
Africa	GOB	18.3	−3.52
	TAM	40.56	−17.53
Greenland	ALE	39.75	−27.06
All		32.76	0.29

and

Table 10. Validation results of daily ISR in different areas.

Area	RMSE (W/m2)	Bias (W/m2)
North America	34.57	−3.08
Europe	29.36	4.74
South America	42.98	2.52
Oceania	27.41	0.52
Asia	30.74	10.24
Africa	29.08	−8.98
Greenland	39.75	−27.06
All	32.76	0.29

. The ISR in North America RMSE ranges from 26.41 W/m² to 39.96 W/m² at the different sites while the bias ranges from −8.89 W/m² to 2.2 W/m². The overall RMSE and bias are 34.57 W/m² and −3.08 W/m², respectively.

The ISR in Europe RMSE ranges from 25.4 W/m² to 33.07 W/m² at the different sites while the bias ranges from −0.35 W/m² to 6.97 W/m². The overall RMSE and bias are 29.36 W/m² and 4.74 W/m², respectively.

The ISR RMSE in South America ranges from 34.13 W/m² to 55.09 W/m² at the different sites while the bias ranges from −3.05 W/m² to 9.4 W/m². The overall RMSE and bias are 42.98 W/m² and 2.52 W/m², respectively.

The ISR RMSE in Oceania ranges from 17.41 W/m² to 34.63 W/m² at the different sites while the bias ranges from −7.44 W/m² to 11.02 W/m². The overall RMSE and bias are 27.41 W/m² and 0.52 W/m², respectively.

The ISR RMSE in Asia ranges from 27.8 W/m² to 34.04 W/m² at the different sites while the bias ranges from 7.78 W/m² to 14.53 W/m². The overall RMSE and bias are 30.74 W/m² and 10.24 W/m², respectively.

The ISR RMSE in Africa ranges from 18.3 W/m² to 40.56 W/m² at the different sites while the bias ranges from −17.53 W/m² to−3.52 W/m². The overall RMSE and bias are 29.08 W/m² and −8.98 W/m², respectively.

The RMSE and bias in the Greenland area site are 39.75 W/m² and −27.06 W/m², respectively.

The validation results of daily ISR are shown in Figure 9, Table 9 and Table 10. The ISR RMSE in each area ranges from 28.49 W/m² to 39.75 W/m² at the different sites while the bias ranges from −27.06 W/m² to 10.24 W/m². The overall RMSE and bias are 32.75 W/m² and 0.29 W/m², respectively.

5. Discussion and Conclusions

5.1. Spatial and Temporal Analysis

Figure 10 shows the seasonal variation of relative biases and RMSEs over the different regions. Most of the northern hemisphere sites, including North America, Europe, and Asia sites overestimate ISR in summer and underestimate in winter. Tropical and sub-tropical sites, including South America and Oceania sites, show no significant seasonal changes on biases. Greenland site shows significant underestimation as the algorithm underestimates the cloud optical depth due to the snow land surface, which should be improved in future research. Most of the sites have larger relative RMSE in winter, as the absolute ISR values are usually smaller in winter.

Figure 11 and Figure 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.

5.2. 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.

5.3. Conclusions

This research modified the optimization based ISR estimation algorithm and applied it to VIIRS data. We validated instantaneous and daily ISR globally over SURFRAD and BSRN sites and performed the spatial and temporal analysis.

The ISR products are still not available from VIIRS observation. Here we present a framework for estimating instantaneous and daily ISR products from VIIRS observation. The algorithm was used for MODIS instantaneous ISR estimation. We adapted and extended the algorithm for VIIRS data. We validated the results at seven SURFRAD and 33 BSRN sites. The RMSE over SURFRAD sites were 83.76 W/m² for instantaneous ISR and 66.80 W/m² and instantaneous NSR. Daily results showed 27.78 W/m² and 23.51 W/m² for ISR and NSR, respectively. Validation results at BSRN sites present RMSEs of 105.87 W/m² for instantaneous ISR and 32.76 W/m² for daily ISR.

The algorithm can estimate instantaneous, and daily ISR from VIIRS data at a similar or better accuracy than existing products and can provide much higher spatial resolution. However, the optimization-based method generated greater uncertainty compared to that of the previous MODIS algorithm, mainly because the latter combines observations from the local morning and local afternoon. Further research can be undertaken to improve the cost function and the optimization framework in the case of a lack of high-level products as constraints and apply the algorithm in the global operational production of 375/750 m resolution.