A Lookup-Table-Based Approach to Estimating Surface Solar Irradiance from Geostationary and Polar-Orbiting Satellite Data

Incoming surface solar irradiance (SSI) is essential for calculating Earth’s surface radiation budget and is a key parameter for terrestrial ecological modeling and climate change research. Remote sensing images from geostationary and polar-orbiting satellites provide an opportunity for SSI estimation through directly retrieving atmospheric and land-surface parameters. This paper presents a new scheme for estimating SSI from the visible and infrared channels of geostationary meteorological and polar-orbiting satellite data. Aerosol optical thickness and cloud microphysical parameters were retrieved from Geostationary Operational Environmental Satellite (GOES) system images by interpolating lookup tables of clear and cloudy skies, respectively. SSI was estimated using pre-calculated offline lookup tables with different atmospheric input data of clear and cloudy skies. The lookup tables were created via the comprehensive radiative transfer model, Santa Barbara Discrete Ordinate Radiative Transfer (SBDART), to balance computational efficiency and accuracy. The atmospheric attenuation effects considered in our approach were water vapor absorption and aerosol extinction for clear skies, while cloud parameters were the only atmospheric input for cloudy-sky SSI estimation. The approach was validated using one-year pyranometer measurements from seven stations in the SURFRAD (SURFace RADiation budget network). The results of the comparison for 2012 showed that the estimated SSI agreed with ground measurements with correlation coefficients of 0.94, 0.69, and 0.89 with a bias of 26.4 W/m2, −5.9 W/m2, and 14.9 W/m2 for clear-sky, cloudy-sky, and all-sky conditions, respectively. The overall root mean square error (RMSE) of instantaneous SSI was 80.0 W/m2 (16.8%), 127.6 W/m2 (55.1%), and 99.5 W/m2 (25.5%) for clear-sky, cloudy-sky (overcast sky and partly cloudy sky), and all-sky (clear-sky and cloudy-sky) conditions, respectively. A comparison with other state-of-the-art studies suggests that our proposed method can successfully estimate SSI with a maximum improvement of an RMSE of 24 W/m2. The clear-sky SSI retrieval was sensitive to aerosol optical thickness, which was largely dependent on the diurnal surface reflectance accuracy. Uncertainty in the pre-defined horizontal visibility for ‘clearest sky’ will eventually lead to considerable SSI retrieval error. Compared to cloud effective radius, the retrieval error of cloud optical thickness was a primary factor that determined the SSI estimation accuracy for cloudy skies. Our proposed method can be used to estimate SSI for clear and one-layer cloud sky, but is not suitable for multi-layer clouds overlap conditions as a lower-level cloud cannot be detected by the optical sensor when a higher-level cloud has a higher optical thickness.

The variations in cloud vertical structures and morphology affect the atmospheric circulation, radiation budget, and satellite-retrieved cloud properties. Most of the sensor-received radiance came from the top of cloud for conditions in which upper optical thick cloud overlaps the lower optical thin cloud. The maximum difference of sensor-received radiance and surface-received radiance is about 4 W/m 2 (7%) and 75 W/m 2 (60%) for upper cloud having COT of 70 when lower cloud COT changes from 1 to 100 (Figures 1 and 2. Sensor and surface-received radiance were estimated using the following parameters: solar zenith angle is 30°, surface albedo is 0.2, cloud phase is water cloud, upper cloud top height is 8 km, upper cloud base height is 6 km, lower cloud top height is 3 km, lower cloud base height is 1 km, cloud particle effective radius is 6 μm). For upper cloud having higher COT, lower-level cloud has a minor impact on satellite-retrieved cloud properties but a larger impact on SSI. Our proposed method for SSI estimation is suitable for one-layer cloud sky, and larger errors may be introduced for multi-layer clouds overlap conditions.

Geostationary Images
The data used in this study was acquired from the third-generation GOES-13 (Geostationary Operational Environmental Satellite System) satellite operated by the national environmental satellite, data, and information service of the National Oceanic and Atmospheric Administration (NOAA). GOES-13 was used for weather forecasting, severe storm tracking, and meteorology research. GOES-13 was launched on 24 May 2006, and is positioned at 75°W, 35,786 km over the Equator. The imager on-board GOES-13 scans Earth's surface every 30 min and provides five spectral channels. The nadir spatial resolution is 1 km for the visible channel (0.65 μm), 4 km for The variations in cloud vertical structures and morphology affect the atmospheric circulation, radiation budget, and satellite-retrieved cloud properties. Most of the sensor-received radiance came from the top of cloud for conditions in which upper optical thick cloud overlaps the lower optical thin cloud. The maximum difference of sensor-received radiance and surface-received radiance is about 4 W/m 2 (7%) and 75 W/m 2 (60%) for upper cloud having COT of 70 when lower cloud COT changes from 1 to 100 (Figures 1 and 2. Sensor and surface-received radiance were estimated using the following parameters: solar zenith angle is 30°, surface albedo is 0.2, cloud phase is water cloud, upper cloud top height is 8 km, upper cloud base height is 6 km, lower cloud top height is 3 km, lower cloud base height is 1 km, cloud particle effective radius is 6 μm). For upper cloud having higher COT, lower-level cloud has a minor impact on satellite-retrieved cloud properties but a larger impact on SSI. Our proposed method for SSI estimation is suitable for one-layer cloud sky, and larger errors may be introduced for multi-layer clouds overlap conditions.

Geostationary Images
The data used in this study was acquired from the third-generation GOES-13 (Geostationary Operational Environmental Satellite System) satellite operated by the national environmental satellite, data, and information service of the National Oceanic and Atmospheric Administration (NOAA). GOES-13 was used for weather forecasting, severe storm tracking, and meteorology research. GOES-13 was launched on 24 May 2006, and is positioned at 75°W, 35,786 km over the Equator. The imager on-board GOES-13 scans Earth's surface every 30 min and provides five spectral channels. The nadir spatial resolution is 1 km for the visible channel (0.65 μm), 4 km for

Geostationary Images
The data used in this study was acquired from the third-generation GOES-13 (Geostationary Operational Environmental Satellite System) satellite operated by the national environmental satellite, data, and information service of the National Oceanic and Atmospheric Administration (NOAA). GOES-13 was used for weather forecasting, severe storm tracking, and meteorology research. GOES-13 was launched on 24 May 2006, and is positioned at 75 • W, 35,786 km over the Equator. The imager on-board GOES-13 scans Earth's surface every 30 min and provides five spectral channels. The nadir spatial resolution is 1 km for the visible channel (0.65 µm), 4 km for three thermal infrared channels Remote Sens. 2018, 10, 411 5 of 19 (3.9 µm, 6.48 µm, and 10.7 µm), and 8 km for channel 6 (13.3 µm). Details can be seen from http://www.ssec.wisc.edu/datacenter/standard_GOES8-15.html, and data can be downloaded freely from http://www.class.ncdc.noaa.gov/saa/products/welcome.

Ancillary Input Data
The MCD43D (V006) surface albedo product derived from the combined Terra and Aqua satellites was used in this study. The bidirectional reflectance distribution function (BRDF) was estimated from all cloud-free observations during a 16-day period. The MCD43D product incorporates the Climate Modeling Grid (CMG) structure and the pixel resolution is 1000 m. The broadband (0.2-4.0 µm) surface albedo for clear skies was calculated as the interpolation between the white-sky and black-sky albedo values dependent on the aerosol optical depth and solar zenith. Only the white-sky albedo product was used for cloudy skies due to minor differences discovered when introducing black-sky albedo for direct beam reflection [41].
The NCEP Climate Forecast System Reanalysis (CFSR) data created using the National Centers for Environmental Prediction (NCEP) Climate Forecast System version 2 (CFSv2) (https://rda.ucar. edu/datasets/ds094.1/#!description) was used in our study. The files in this dataset were grouped by month. The grid spacing was 0.205~0.204 • from 0 • E to 359.795 • E, and 89.843 • N to 89.843 • S (1760 × 880 Longitude/Gaussian Latitude) [42]. Ground surface temperature was selected to be an ancillary input for the cloud effective radius retrieval method for the GEOS-13 infrared channel. The precipitable water of the entire atmosphere was selected to drive the SSI retrieval algorithm, since the 12.0-µm channel was replaced by a 13.3-µm channel. Furthermore, the GOES-13 satellite and the retrieval of precipitable water from the "split window" method (using channels 11.0 µm and 12.0 µm) was inapplicable.
The global 1-km Shuttle Radar Topography Mission with 30 arc-second resolution data (SRTM30) was used to represent the surface elevation (http://vterrain.org/Elevation/SRTM/) required for the retrieval of SSI.

Pyranometer Data for Validation
The Surface Radiation Budget Network (SURFRAD) was established in 1993 with the support of NOAA's Office of Global Programs. Its primary mission was to support climate research using accurate, continuous, and long-term measurements of the surface radiation budget over the United States. Seven SURFRAD stations are currently operating in climatologically diverse regions in the United States, including Fort Peck, Montana (FPK), Table Mountain, Colorado (TBL), Bondville, Illinois (BON), Goodwin Creek, Mississippi (GWN), Penn State, Pennsylvania (PSU), Desert Rock, Nevada (DRA), and Sioux Falls, South Dakota (SXF). The downwelling global solar irradiance (0.28-3 µm) is measured by a pyranometer (model SpectroSun SR-75) with reported uncertainties of ±2% to ±5% [43]. SURFRAD data are provided daily with a sample rate of 1 min (https://www.esrl.noaa.gov/gmd/grad/surfrad/). The maintenance and quality control of these measurements follow World Meteorological Organization (WMO) standards.

Methods
SSI is retrievable by assuming a homogeneous and plane-parallel atmospheric layer without considering the three-dimensional effects. The discrimination of clear and cloudy conditions was implemented by a cloud detection procedure, and the cloud thermodynamic phase was retrieved using IR channels. The cloud parameters (cloud optical thickness and effective particle radius) were inversed from the visible channel and IR channels based on the previous work of Nakajima [37,39]. SSI was estimated using a LUT-based method with the atmospheric and land-surface parameters derived above. The proposed SSI retrieval scheme is given in Figure 3. Clear and cloudy skies were first labeled using the cloud detection procedure. Aerosol optical depth (AOD) and precipitable water were retrieved for clear skies and cloud microphysical parameters were derived for cloudy skies using the pre-calculated LUT. SSI was calculated using the LUT for both clear and cloudy skies. Cloud

Pre-Processing the Images
Surface reflectance can be estimated from the visible band's at-sensor spectral radiance under clear-sky conditions through atmospheric radiative transfer models, such as the Santa Barbara DISORT Atmospheric Radiative Transfer (SBDART) [44]. The top-of-atmosphere reflectance is converted into surface reflectance given the solar-sensor geostationary viewing geometry, Rayleigh scattering, well-mixed gaseous absorption, ozone and water vapor absorption, and aerosol extinction through atmospheric correction. An aerosol visibility of 100 km and a rural model is used to represent clear atmospheric conditions. The water vapor and other trace gases are initialized with the default values of the SBDART model. Surface reflectance is determined by the minimum reflectance retrieved from the visible band taken at the same local time per daylight hour over a temporal period of one month for cloud-free detection due to the difficulty of discriminating the "clearest" atmospheric conditions. Details of the proposal have been discussed by Liang et al. (2006) [28] and Zhang et al. (2014) [4]. A 30° threshold on the glint cone angle was applied to avoid sun-glint affecting water surfaces, and a lower reflectance threshold of 0.005 was applied for the land surface to exclude cloud shadow pixels [45].
Cloud detection was performed pixel-by-pixel using the coupled Cloud Depiction and Forecast System model using the reflectance of visible bands and the brightness temperature of infrared bands [46]. This procedure incorporated temporal differencing, dynamic thresholding, and spectral discrimination to detect clouds with the appropriate optical thickness.

Pre-Processing the Images
Surface reflectance can be estimated from the visible band's at-sensor spectral radiance under clear-sky conditions through atmospheric radiative transfer models, such as the Santa Barbara DISORT Atmospheric Radiative Transfer (SBDART) [44]. The top-of-atmosphere reflectance is converted into surface reflectance given the solar-sensor geostationary viewing geometry, Rayleigh scattering, well-mixed gaseous absorption, ozone and water vapor absorption, and aerosol extinction through atmospheric correction. An aerosol visibility of 100 km and a rural model is used to represent clear atmospheric conditions. The water vapor and other trace gases are initialized with the default values of the SBDART model. Surface reflectance is determined by the minimum reflectance retrieved from the visible band taken at the same local time per daylight hour over a temporal period of one month for cloud-free detection due to the difficulty of discriminating the "clearest" atmospheric conditions. Details of the proposal have been discussed by Liang et al. (2006) [28] and Zhang et al. (2014) [4]. A 30 • threshold on the glint cone angle was applied to avoid sun-glint affecting water surfaces, and a lower reflectance threshold of 0.005 was applied for the land surface to exclude cloud shadow pixels [45].
Cloud detection was performed pixel-by-pixel using the coupled Cloud Depiction and Forecast System model using the reflectance of visible bands and the brightness temperature of infrared bands [46]. This procedure incorporated temporal differencing, dynamic thresholding, and spectral discrimination to detect clouds with the appropriate optical thickness.

Aerosol Optical Depth Estimation
Aerosol plays a key role in Earth's radiation budget by scattering and absorbing solar and terrestrial radiation. The single broadband visible channel of most geostationary satellites is not sufficient to retrieve the aerosol size and single scattering albedo, although they are important for radiation extinction. The AOD was retrieved using the visible band of GOES with a pre-calculated LUT. The dimensions of the LUT are summarized in Table 1. The rural type was defined as incorporated in the SBDART radiative transfer code with a single scattering albedo at 0.55 µm of 0.9558 and an asymmetry factor of 0.6891. The standard atmospheric profile of the midlatitude summer model was used as the default. The LUT was pre-generated using the SBDART model for a range of discrete atmospheric and land-surface values to improve the calculation efficiency without reducing accuracy. SBDART was numerically integrated with Discrete Ordinate Radiative Transfer (DISORT), which assumes a plane-parallel radiative transfer in a vertically inhomogeneous atmosphere. The number of streams for radiance computations was 20 for the zenith angle and azimuth angle. The surface reflectance and cloud mask was determined for each pixel as described in Section 2.2.1. Aerosol horizontal visibility (VIS) was computed for every cloud-free pixel using the rural aerosol model and the given solar position, satellite position, amount of water vapor, and surface altitude. A linear interpolation of the lookup table entries to the actual aerosol visibility was used in this study. Once the VIS was known, the AOD (at 550 nm) was estimated using the following equation [44] (http://www.ncgia.ucsb.edu/projects/metadata/standard/uses/sbdart.htm): where W is a weighting factor, which is a piecewise function depending on the value of VIS and is given by the following equation:

Retrieving Cloud Microphysical Properties
The cloud thermodynamic phase, cloud optical thickness (COT), and effective particle radius were used to describe the radiative properties of clouds in the solar spectral region. The thermodynamic phases of the cloud were classified as: water clouds, ice clouds, mixed clouds, and undetected clouds following the cloud phase determination proposed by Choi et al. (2007) [47]. The retrieval method was based on the theory that cloud reflectance at non-absorbing wavelengths of the visible band is strongly related to COT, while the reflection at the absorbing wavelengths of the near infrared bands is primarily a function of ER [37]. In this study, the visible channel was used to derive COT and the IR3.9 channel was chosen to obtain ER. The radiance received by the sensor at 3.9 µm (L obs 3.9 ) was composed of solar reflection, cloud thermal radiance, and ground thermal radiance for thin clouds. The radiance for 0.65 µm and 3.9 µm is given simply as follows: where L cloud 0.65 and L cloud 3.9 are the cloud-reflected radiance at the VIS and IR3.9 channels, respectively, L sr 0.65 and L sr 3.9 are the ground-reflected radiance at the VIS and IR3.9 channels, respectively. L th(cloud) 3.9 and L th(sr) 3.9 are the cloud and ground thermal radiance, respectively. L sr 3.9 and L th(sr) 3.9 were assumed to be 0 for thick clouds (COT >16). L sr 3.9 and L th(sr) 3.9 were simulated based on the Planck function of ground temperature (T g ) and cloud-top temperature (T c ) to remove the thermal effects of ER retrieval for thin clouds (COT <16). T g data were derived from the NCEP-CFSv2 dataset, and T c was approximated by the cloud-top brightness temperature given by the IR channel at 10.7 µm.
COT and ER were retrieved using LUTs generated from the SBDART one-dimensional radiative transfer code. The LUTs were calculated for different values of COT, ER, solar zenith angle, satellite viewing angle, relative azimuth angle, surface albedo, surface temperature, and cloud-top temperature using the different spectral response functions of the visible and infrared bands ( Table 2). These calculations were carried out under the following assumptions: (1) there is only one single layer of clouds for every pixel, (2) the clouds are homogeneous, plane-parallel, and cover the whole pixel, and (3) ice clouds are composed of spherical particles. COT and ER were assigned to be the average value of water and ice clouds for mixed-phase clouds.

All-Sky SSI Estimation
SSI was estimated separately for clear and cloudy skies with different input data using AOD data and cloud physical parameters efficiently derived from geostationary images (as discussed in Sections 2.2.2 and 2.2.3). CO 2 and ozone were set to default values in SSI estimation since they had a negligible impact. PW and aerosol had a considerable influence on SSI in cloud-free conditions. Clouds played a dominant role in SSI during cloudy-sky conditions, and PW was set at 2.9 g/cm 2 as defined in the standard atmospheric profile of the midlatitude summer model. Aerosol horizontal visibility was set to 100 km for the SSI estimation of cloudy skies since AOD was insignificant compared to clouds and difficult to derive under cloudy conditions.
The all-sky SSI estimation was derived using LUTs generated for clear and cloudy skies. The common variables used for the LUTs were the solar zenith angle, surface altitude, and surface albedo. The LUT atmospheric variables for clear skies were PW and aerosol visibility, while the LUT for cloudy skies contained cloud phase, COT, and ER. The SSI for "mixed-phase clouds" was assigned to be the averaged SSI estimation for water and ice clouds. The SSI for "undetected cloud phase" pixels was calculated using the LUT of water clouds. The range of values and the increments of the above variables were the same as in Tables 1 and 2. The instantaneous SSI was estimated by linear interpolation from the lookup table once the above input data were known.

Results
In this section, the algorithm discussed above is evaluated using the data from seven SURFRAD stations during the entire year of 2012 and a comparison is performed with other SSI estimates. The performance of the SSI estimate is evaluated using three metrics: the mean bias error (MBE, in W/m 2 ), RMSE (in W/m 2 ), and correlation coefficient (R 2 ). Huang et al. (2016) [48] demonstrated that the observed SSI averaged over 30 min was optimal for a comparison with kilometer-level satellite-based SSI estimation. Therefore, we adopted half-hour averaged SSI observations centered at the acquired time of the satellite images to evaluate the satellite-derived instantaneous SSI estimation. The validation results gathered from seven SURFRAD stations in 2012 under clear-and cloudy-sky conditions are displayed in Figure 4 and the statistics are compared in phase" pixels was calculated using the LUT of water clouds. The range of values and the increments of the above variables were the same as in Tables 1 and 2. The instantaneous SSI was estimated by linear interpolation from the lookup table once the above input data were known.

Results
In this section, the algorithm discussed above is evaluated using the data from seven SURFRAD stations during the entire year of 2012 and a comparison is performed with other SSI estimates. The performance of the SSI estimate is evaluated using three metrics: the mean bias error (MBE, in W/m 2 ), RMSE (in W/m 2 ), and correlation coefficient (R 2 ). Huang et al. (2016) [48] demonstrated that the observed SSI averaged over 30 min was optimal for a comparison with kilometer-level satellite-based SSI estimation. Therefore, we adopted half-hour averaged SSI observations centered at the acquired time of the satellite images to evaluate the satellite-derived instantaneous SSI estimation. The validation results gathered from seven SURFRAD stations in 2012 under clear-and cloudy-sky conditions are displayed in Figure 4 and the statistics are compared in Table 3 These statistics indicate that the quality of the retrieval was better for clear skies, which had a correlation coefficient (R 2 ) ranging from 0.9 to 0.96, in comparison to cloudy skies for all stations, which featured a correlation coefficient ranging from 0.60 to 0.80; this was true for both systematic bias and scatter ( Figure 4). The largest RMSE values for clear-and all-sky conditions both occurred at Table Mountain,     These statistics indicate that the quality of the retrieval was better for clear skies, which had a correlation coefficient (R 2 ) ranging from 0.9 to 0.96, in comparison to cloudy skies for all stations, which featured a correlation coefficient ranging from 0.60 to 0.80; this was true for both systematic bias and scatter (Figure 4). The largest RMSE values for clear-and all-sky conditions both occurred at Table  Mountain, while the smallest RMSE values for clear-and all-sky conditions occurred at Desert Rock. All stations exhibited a positive bias for clear-and all-sky conditions.
Further investigation was carried out in our study due to a larger positive bias being discovered in Table Mountain (TBL) compared to other stations with clear skies. The surface of the TBL station in Colorado was mixed by rocks, sparse grasses, desert shrubs, and small cactus, and the surface altitude was 1689 m. The positive bias was partially due to the errors of cloud detection for a mixed surface with a higher altitude. Some pixels covered by thin clouds or haze were classified as "clear sky", and thus resulted in an overestimation of SSI. Nevertheless, as is well-known, the aerosol "dark target" approach is only valid for a dense dark vegetation (DDV) surface, and it is inappropriate for the TBL station with a lower vegetation fractional cover, and thus will generally lead to substantial errors in the retrieved AOD.
On the "station observation" scale, clouds generally deviate much more under the horizontal/ vertical homogeneity assumption of the SSI estimation approach than other atmospheric variables, such as aerosol and total water vapor. The inhomogeneous properties of clouds may cause substantial errors in retrieving cloud optical thickness from visible channels with a 1-km resolution and an effective particle radius from an infrared channel with approximately 4 km from satellite data. The larger discrepancies for cloudy-sky SSI estimation may be attributed to the horizontal/vertical inhomogeneity of clouds and the spatial observing scale mismatches in sensor footprints between ground-observed and satellite-retrieved data. The negative effects of the mismatches will be enlarged for a lower solar zenith and viewing zenith, resulting in poorer SSI estimation and evaluation accuracy, especially for partially covered clouds or broken clouds.

Comparison with Other SSI Estimates
SSI estimation with in situ observations at SURFRAD sites were collected in order to compare the accuracy of our proposed algorithm with previous studies that estimate SSI from geostationary and polar-orbiting satellite data. The results are listed in Tables 4 and 5. Zhang et al. (2014) used a LUT-based method from geostationary satellite images to estimate incident shortwave radiation at 5-km resolution, which was validated with observation data at seven SURFRAD sites of 2008 (Table 4) [4]. The results revealed that the RMSE values produced by our proposed method were less than the values provided by Zhang's estimation, apart from the validation at GWN, which had RMSE values of 90.7 W/m 2 and 86 W/m 2 , respectively. Our proposed model exhibited an overall positive bias at all seven sites, while Zhang's model provided a negative bias at DRA and TBL. The largest bias in our model was 33.8 W/m 2 at TBL, compared with −55 W/m 2 produced by Zhang's method at DRA.  [29]. Different validation results can be examined between our model and Qin's method. Qin's method yielded a better performance for clear sky at all sites. The unfavorable comparison results may attribute to the inaccurate input data of our model with total precipitable water at approximately 20-km resolution and the uncertainty of retrieved AOD which will be discussed later. Our model provided an improved performance with a lesser RMSE of 1-12 W/m 2 compared to Qin's method for all-sky conditions at BON, FPK, PSU, SXF, and GWN, which used input data from Aqua, while Qin's method yielded values in agreement at DRA, TBL, and GWN with input data from Terra. All three methods yielded a poorer validation at TBL; this might have been caused by pyranometer calibration accuracy, climatic conditions, and mixed ground cover. Furthermore, our proposed model indicated a lesser RMSE of about 2-4 W/m 2 compared to the method provided by Tang et al. (2016) [14], which combined an artificial neural network and parameterization model for SSI estimation from multifunctional transport satellite (MTSAT) geostationary satellite images and MODIS atmospheric and land products. The overall accuracy of our model with an RMSE of 99.5 W/m 2 (25.5%) is comparable to the MODIS-products-driven Breathing Earth System Simulator (BESS) shortwave products with an RMSE of 111.1 W/m 2 (22.6%) and 137.1 W/m 2 (31.7%) for temperate and continental climate zones, respectively [49].
As indicated by Yeom, a reduced RMSE of about 10 W/m 2 can be found with the spatial resolution changed from 1 km to 5 km [13]. Considering the estimated SSI of our model with 1-km resolution and the referenced studies with 5-km resolution, we can draw a conclusion that the performance of our proposed scheme was comparable with or even more accurate than state-of-the-art satellite-based SSI retrieval models.

Error Analysis in SSI Retrieval
Aerosol and clouds are the primary atmospheric parameters (besides the solar zenith and surface altitude) that affect SSI for clear and cloudy skies. The retrieval uncertainty of these two parameters will be discussed in this section.
The diurnal change of the underlying surface reflectance is a key parameter for AOD retrieval, and it is gathered by searching for the minimum value of surface reflectance within a 30-day period. The surface reflectance was inversed using a lookup-table-based method and a horizontal visibility set to 100 km (which was approximated to be 0.06 of the AOD value using the relationship between the VIS and AOD as indicated by Equations (1) and (2)). However, the assumption will inevitably introduce some uncertainty since a great spatial and temporal variation of aerosol has been discovered. The AOD data from SURFRAD sites generated from visible Multi-Filter Rotating Shadow band Radiometers (MFRSR) were collected in our study to further investigate the changes in AOD. The statistical results of observed AOD gathered from six SURFRAD sites are displayed in Figure 5 (except for PSU since there was no observed AOD data in 2012). Cloud microphysical parameters, such as cloud optical thickness, thermal phase, and effective particle radius, are important variables in estimating SSI. Large negative effects on the performance of the SSI retrievals were possibly caused by cloud parameter retrieval errors due to the inhomogeneity and spatiotemporal variation of clouds. The sensitivity of SSI to cloud optical thickness and the effective radius of water clouds is presented in Figure 6. SSI was estimated using the rigorous radiative transfer model (SBDART) with the input variables set as: a midlatitude summer atmospheric profile with total precipitable water of 2.92 g/cm 2 , a total ozone column of 320 DU, a solar zenith angle of 30°, a surface elevation of 0 km, a surface albedo of 0.2, and a horizontal visibility of 100 km. A positive relationship between the effective particle radius and SSI can be seen, while a negative relationship between the SSI and cloud optical thickness can be observed. SSI will change about 88 W/m 2 and 90 W/m 2 with the effective particle radius range of 2 μm to 30 μm for water and ice clouds, respectively. SSI is not dependent on the effective radius when the cloud particle radius is greater than 20 μm since the variation of SSI does not exceed 5 W/m 2 . The SSI will change by about 548 and 600 W/m 2 for a cloud optical thickness range of 2 to 30 for water and ice clouds, respectively. It can be concluded that the error in SSI estimation for cloudy skies is primarily affected by the uncertainty of the cloud optical thickness retrieval. There may be more than one solution for optical thickness and effective radius retrieval for optically thin clouds. The retrieved cloud optical thickness only represents 20-40% of the total optical thickness of the total cloud layer for clouds having COT ≥8, as indicated by Nakajima [37]. This situation can partially explain the positive bias of SSI estimation for cloudy skies. A high yearly temporal variation can be found at all six sites. A maximum AOD value was found in summer and a minimum value in winter, and a maximum variability of AOD occurred during the summer months. The pre-defined AOD of 0.06 for clearest-sky conditions had: an overall underestimation at BON, SXF, and GWN for the entire year of 2012; an overall underestimation in summer and overestimation in winter at TBL and DRA; and an overall underestimation in summer at FPK. A lower assumed AOD may have been responsible for the SSI overestimation in our proposed method. The surface reflectance was overestimated due to the lower predefined AOD value for clearest-sky conditions, which in turn resulted in underestimations of the retrieved AOD and the overall overestimation of SSI for clear skies. The surface reflectance was invariable for every month, and larger uncertainty had arisen for snow seasons with higher reflectance and vegetation-growing seasons with a lower reflectance. Furthermore, an overestimation of surface reflectance caused an underestimation of cloud optical thickness and an underestimation of cloud transmittance for transparency or semi-transparency. Besides the uncertainty of surface reflectance, the aerosol attenuation effects were influenced by a large solar zenith angle and bright surface.
Cloud microphysical parameters, such as cloud optical thickness, thermal phase, and effective particle radius, are important variables in estimating SSI. Large negative effects on the performance of the SSI retrievals were possibly caused by cloud parameter retrieval errors due to the inhomogeneity and spatiotemporal variation of clouds. The sensitivity of SSI to cloud optical thickness and the effective radius of water clouds is presented in Figure 6. SSI was estimated using the rigorous radiative transfer model (SBDART) with the input variables set as: a midlatitude summer atmospheric profile with total precipitable water of 2.92 g/cm 2 , a total ozone column of 320 DU, a solar zenith angle of 30 • , a surface elevation of 0 km, a surface albedo of 0.2, and a horizontal visibility of 100 km. A positive relationship between the effective particle radius and SSI can be seen, while a negative relationship between the SSI and cloud optical thickness can be observed. SSI will change about 88 W/m 2 and 90 W/m 2 with the effective particle radius range of 2 µm to 30 µm for water and ice clouds, respectively. SSI is not dependent on the effective radius when the cloud particle radius is greater than 20 µm since the variation of SSI does not exceed 5 W/m 2 . The SSI will change by about 548 and 600 W/m 2 for a cloud optical thickness range of 2 to 30 for water and ice clouds, respectively. It can be concluded that the error in SSI estimation for cloudy skies is primarily affected by the uncertainty of the cloud optical thickness retrieval. There may be more than one solution for optical thickness and effective radius retrieval for optically thin clouds. The retrieved cloud optical thickness only represents 20-40% of the total optical thickness of the total cloud layer for clouds having COT ≥8, as indicated by Nakajima [37]. This situation can partially explain the positive bias of SSI estimation for cloudy skies.
while a negative relationship between the SSI and cloud optical thickness can be observed. SSI will change about 88 W/m 2 and 90 W/m 2 with the effective particle radius range of 2 μm to 30 μm for water and ice clouds, respectively. SSI is not dependent on the effective radius when the cloud particle radius is greater than 20 μm since the variation of SSI does not exceed 5 W/m 2 . The SSI will change by about 548 and 600 W/m 2 for a cloud optical thickness range of 2 to 30 for water and ice clouds, respectively. It can be concluded that the error in SSI estimation for cloudy skies is primarily affected by the uncertainty of the cloud optical thickness retrieval. There may be more than one solution for optical thickness and effective radius retrieval for optically thin clouds. The retrieved cloud optical thickness only represents 20-40% of the total optical thickness of the total cloud layer for clouds having COT ≥8, as indicated by Nakajima [37]. This situation can partially explain the positive bias of SSI estimation for cloudy skies. Figure 6. The sensitivity of SSI to cloud optical thickness (given an effective particle radius of 20 µm) and effective particle radius (given a cloud optical thickness of 20) for water clouds (left) and ice clouds (right).
An overall underestimation of SSI with a maximum bias of 62.5 W/m 2 and a minimum of 1 W/m 2 for water clouds is indicated in Figure 7 and Table 6. All scatterplots for water and ice clouds are uniformly distributed between the line of 1:1. The RMSE has a maximum value of 177.6 W/m 2 at DRA, thereby exceeding the values for clear skies (61.4 W/m 2 ) by a factor of three. The minimum RMSE value is 80.2 W/m 2 at SXF, which is nearly the same as clear skies. Compared with water clouds, ice clouds tend to have a larger RMSE and a lower correlation coefficient ( Figure 8, Table 6). The validation results for ice cloud cases reveal a negative bias at DRA and TBL, which have an elevation greater than 1000 m. The largest RMSE for ice clouds was 211.3 W/m 2 , which is about 75.3% of the systematic error. The relative accuracy of the modeled SSI for ice cloud cases is lower than 40%. This might be caused by the ice crystal density, particle size, shape, or direction, which are difficult to derive and describe for accurate scattering computation.
Besides the uncertainty of the input variables derived from satellite data, large uncertainty may arise from the assumption of plane-parallel, homogeneous atmospheric conditions. Furthermore, a one-dimensional atmospheric transfer model cannot deal with the geometrical effects of scattering from a higher solar zenith. The model is less reliable when the sub-pixel is partially cloudy, or when a rapid change in the atmospheric profile occurs during satellite observations.  Figure 7. Validation of SSI for water cloud cases at SURFRAD sites (mixed and undetected clouds were not included in the comparison).

Conclusions
The paper describes a novel approach to estimating surface solar irradiance (SSI) from a combined geostationary satellite image, MODIS land-surface albedo, and NCEP CFSR data. Aerosol optical thickness and cloud parameters (cloud phase, effective particle radius, and cloud optical thickness) were directly retrieved from the visible channel of geostationary satellite images for clear and cloudy skies. Total precipitable water was derived from the NCEP data, and other atmospheric

Conclusions
The paper describes a novel approach to estimating surface solar irradiance (SSI) from a combined geostationary satellite image, MODIS land-surface albedo, and NCEP CFSR data. Aerosol optical thickness and cloud parameters (cloud phase, effective particle radius, and cloud optical thickness) were directly retrieved from the visible channel of geostationary satellite images for clear and cloudy skies. Total precipitable water was derived from the NCEP data, and other atmospheric variables, such as ozone, carbon dioxide, and trace gases, were not considered in our SSI estimation. The SSI was obtained by searching for and linearly interpolating a pre-calculated lookup table, which was created using the one-dimensional radiative transfer model for computational efficiency at the cost of calculation accuracy.
The validation was performed via station observations at SURFRAD and other developed algorithms were used with input from satellite data on an instantaneous basis to evaluate the performance of the estimates. The results demonstrated that our method could effectively retrieve instantaneous SSI with correlation coefficients of 0.94, 0.69, and 0.89, and an overall RMSE of 80.0 W/m 2 (16.8%), 127.6 W/m 2 (55.1%), and 99.5 W/m 2 (25.5%) for clear-sky, cloudy-sky, and all-sky conditions, respectively. Our algorithm generally overestimated the SSI for clear-and all-sky conditions. Uncertainty analysis revealed that the accuracy of AOD retrieval was largely dependent upon diurnal surface reflectance. An overestimation of surface reflectance resulted in an underestimation of AOD and led to an overestimation of SSI. Large uncertainty may arise for optically thin clouds due to the ambiguous solutions for cloud optical thickness and effective radius. The RMSE for ice clouds is generally larger than water clouds since the radiative transfer process for ice clouds is mainly affected by ice crystal shape and particle size, which are difficult to directly retrieve with acceptable accuracy. In summary, our proposed method holds great promise for accurately estimating regional or global SSI and conducting research on Earth's energy budget using products from geostationary satellites, such as FY2, Himawari-8, MTG, and MODIS.