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

Abstract

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/m², −5.9 W/m², and 14.9 W/m² for clear-sky, cloudy-sky, and all-sky conditions, respectively. The overall root mean square error (RMSE) of instantaneous SSI was 80.0 W/m² (16.8%), 127.6 W/m² (55.1%), and 99.5 W/m² (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/m². 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.

1. Introduction

Surface Solar Irradiance (SSI) is commonly referred to as the amount of downward solar energy incident to a horizontal surface, and is a major component of the surface energy balance that governs the exchange processes of energy between Earth’s land surface and atmosphere [1,2]. SSI is required by land-surface models, hydrological models, and ecological models to simulate land–atmosphere interactions [3,4]. Accurate observation and estimation of global energy spatial-temporal distribution is essential for climate change monitoring and forecasting [5].

A non-uniform spatial and temporal distribution of SSI has large effects on regional and global climates. However, sparse networks of ground SSI measurements are insufficient for modeling land-surface processes and Earth radiation budget research. Furthermore, fewer surface stations are located in mountainous areas, yet SSI is highly dependent on topography and features larger temporal and spatial variations than horizontal surfaces [6]. Numerous attempts have been made at estimating SSI from satellite data on local or regional scales with multi-scale temporal resolutions in order to overcome the limitations of in situ records [7,8,9,10,11,12,13,14,15]. Perez et al. (1997) demonstrated that satellite-derived irradiation is more accurate compared to interpolation techniques obtained from station measurements if the distance from the station exceeds 34 km for hourly irradiation and 50 km for daily irradiances [16].

Global SSI datasets have been available since the 1990s at different spatiotemporal resolutions based on multi-source remotely sensed data. These include the International Satellite Cloud Climatology Project (ISCCP) [17], the Earth Radiation Budget Experiment (ERBE) [18], the National Centers for Environmental Prediction and National Center for Atmospheric Research Reanalysis Project [19], the Global Energy and Water Cycle Experiment Surface Radiation Budget (GEWEX-SRB), the Clouds and the Earth’s Radiant Energy System (CERES) [20], Satellite Application Facility on Climate Monitoring Solar Surface Radiation Heliosat (CM SAF SARAH) [21], and Global Land Surface Satellite (GLASS) products [22]. The above satellites and their parameter-based meteorological products provide long-term, multiple time-scale global SSI data, but are generally associated with coarser spatial resolutions (e.g., >1°), excluding GLASS and CM SAF SARAH, which have a 5-km resolution and bias estimation. The majority of these products cannot meet the requirements for studying land-surface processes and fail to describe the spatial changes with sufficient accuracy due to their coarse spatial resolutions [14]. Zhang et al. (2015) evaluated four products using 1151 ground sites and found that SSI was generally overestimated by approximately 10 W/m², while the averaged global annual mean SSI from the ground-measured-calibrated value was 180.6 W/m² [23]. Differences range from 10 to 30 W/m², with maximum discrepancies in areas of high cloud cover in the tropics between the ISCCP and ERBE datasets [24]. These differences may be partly due to spatial resolution; in fact, Pinker and Laszlo found an average difference of about 8–9% in daily surface irradiance when adjusting the resolution from 8 to 50 km [25]. Another possible explanation for the bias found in these products can be attributed to the cloud fractional cover and aerosol optical depth [6].

Satellite-based SSI products are useful for historical global SSI analysis, while the general circulation model (GCM) and the numerical weather prediction (NWP) model can be used to estimate SSI at timespans ranging from days to decades using projection scenarios of emissions and land use. The errors obtained from NWP models are generally less than 50 W/m² or exceed 200 W/m² for clear-sky and cloudy-sky conditions, respectively [26]. Lara-Fanego et al. (2012) found the forecast errors produced by Weather Research and Forecasting (WRF) to be 2% under clear-sky conditions and 18% for cloudy skies [27]. Due to the coarse resolution of most WRF models and the GCM, detailed cloud properties and Earth’s energy budget have not been clearly demonstrated. Accurately estimated “kilometer-level” SSI datasets are necessary to overcome the limitations of cloud representations in the climate model.

Besides “single point” ground observations and “kilometer-level” SSI datasets, SSI can be directly retrieved using the relationship between SSI and the top-of-atmosphere (TOA) radiance measured by satellite sensors [4,11,28] or indirectly retrieved through rigorous radiative transfer models (RTMs). However, RTMs are disadvantageous since they are generally time-consuming and require a substantial amount of unavailable detailed atmospheric profile data, and are thus not convenient for applications to large areas with a fine-resolution grid. Semi-empirical models have been developed that contain meteorological variable inputs and feature a parameterized hybrid model with simplified atmospheric transmittance. Some of these models include the pre-computed lookup tables (LUT) method based on RTMs, which has a reduced computational time at the expense of accuracy. These studies considered the extinction and absorption of solar radiation caused by aerosols, water vapor (PW), ozone, and clouds [4,7,14,29].

SSI estimation under cloudy skies is much more complex compared to clear-sky models. The performance of physical SSI models under cloudy skies is largely dominated by cloud macrophysical and microphysical properties, such as cloud fractional cover (CFC), cloud optical thickness (COT), and cloud effective particle radius (ER) with high variability in space and time [30,31]. The Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the Terra and Aqua satellites provides detailed and consistent atmospheric, terrestrial, and oceanic products, and studies have been developed for SSI mapping from pairs of MODIS products [2,32]. Barzin et al. (2017) proposed a combination of principal components analysis (PCA) and regression models for estimating daily average downward solar radiation using MODIS data for ten synoptic stations in Fars Province, Iran, with a root mean square error (RMSE) of 0.9–2.04 MJ/(m²·d) [33]. The largest uncertainties resulting from SSI retrieval arise from inadequate information on cloud properties. Many studies have taken advantage of the fine spatial resolution and higher temporal resolution (5–30 min) of geostationary satellites to derive inhomogeneous and rapidly changing atmospheric parameters. The HELIOSAT algorithm uses a simplified parameterization for cloud transmission, also denoted the cloud index, derived from geostationary satellite measurements and a clear-sky model to calculate sky SSI [34]. Newly improved HELIOSAT-based models have been proposed [35]. An Artificial Neural Network (ANN) can be used to predict SSI simulation from meteorological parameters and satellite images using training data [12,36]. Few studies have focused on SSI estimation for different cloud phases despite the thermodynamic effects of cloud processes being significantly different [14].

The bispectral solar reflectance method has been widely used for retrieving COT and ER from passive satellite multispectral imagers [37]. It has been employed in cloud property retrieval for MODIS [38], the Advanced Very High Resolution Radiometer (AVHRR) [39], and the Spinning Enhanced Visible and Infrared Imager (SEVIRI) [40]. However, the passive optical remote-sensing-based pixel-level COT retrieval uncertainty will be larger than 10% for clouds having COT >70 (MODIS cloud optical properties product Algorithm Theoretical Basis Document for collection 6, MOD06/MYD06-ATBD. See details from https://modis-atmos.gsfc.nasa.gov/sites/default/files/ModAtmo/C6MOD06OPUser-Guide.pdf). The millimeter-wavelength cloud-profiling radar (CPR) on CloudSat and the cloud-aerosol lidar with orthogonal polarization (CALIOP) on Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) provides an opportunity for detailing the structures of clouds, but the temporal resolution of them is too low to monitor the rapid changes of clouds.

This paper presents a lookup-table-based method for all-sky SSI retrieval from combined polar-orbiting and geostationary satellites with rapid retrieval of changing cloud micro-physical properties. The cloud properties retrieval was based on an assumption of a homogeneous one-layer cloud model, and the method is valid for pixels with a solar zenith angle less than 81.4 to better match the “daylight” region as referenced by MOD06/MYD06-ATBD. The SSI estimation approach we proposed should be applicable to clear sky and cloudy sky having COT <70 with less COT retrieval uncertainty as indicated by MOD06/MYD06-ATBD. This paper is organized as follows: Section 2 describes our method and the datasets used in our study to estimate SSI. Validation and a comparison of the results are provided in Section 3 and Section 4. Conclusions are provided in Section 5.

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² (7%) and 75 W/m² (60%) for upper cloud having COT of 70 when lower cloud COT changes from 1 to 100 (Figure 1 and Figure 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.

2. Materials and Methods

2.1. Materials

2.1.1. 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 (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.

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

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

2.2. 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 detection is briefly described in Section 2.2.1. The details for retrieving AOD and cloud microphysical parameters are described in Section 2.2.2 and Section 2.2.3.

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

2.2.2. 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. LUT dimensions for AOD retrieval.

Input Variable	Value Range	Increment
Solar zenith angle	0–89°	5°
Viewing zenith angle	0–89°	5°
Relative azimuth angle	0–180°	30°
Aerosol horizontal visibility	5, 10, 20, 30, 40, 50, 70, 100 km	-
Aerosol type	Rural	-
Water vapor	0.01–5.0 g/cm2	0.5
Surface altitude	0–6 km	1 km
Surface reflectance	0–1.0	0.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):

$$AO{D}_{(0.55}{}_{\mathsf{\mu }\text{m})}=3.912\times \frac{1.05\times W+1.51\times (1-W)}{VIS}$$

(1)

where W is a weighting factor, which is a piecewise function depending on the value of VIS and is given by the following equation:

$$\{\begin{array}{c}\mathrm{W}=\left(\frac{1/\mathrm{VIS}-1/23}{1/5-1/23}\right),\hspace{1em}5<\mathrm{VIS}<23\\ \mathrm{W}=1,\hspace{1em}\mathrm{VIS}<5\\ \mathrm{W}=0,\hspace{1em}\mathrm{VIS}>23\end{array}$$

(2)

2.2.3. 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_{3.9}^{obs}$) 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:

$$L_{0.65}^{obs}=L_{0.65}^{cloud}+L_{0.65}^{sr}$$

(3)

$$L_{3.9}^{obs}=L_{3.9}^{cloud}+L_{3.9}^{sr}+L_{3.9}^{th(cloud)}+L_{3.9}^{th(sr)}$$

(4)

where $L_{0.65}^{cloud}$ and $L_{3.9}^{cloud}$ are the cloud-reflected radiance at the VIS and IR3.9 channels, respectively, $L_{0.65}^{sr}$ and $L_{3.9}^{sr}$ are the ground-reflected radiance at the VIS and IR3.9 channels, respectively. $L_{3.9}^{th(cloud)}$ and $L_{3.9}^{th(sr)}$ are the cloud and ground thermal radiance, respectively. $L_{3.9}^{sr}$ and $L_{3.9}^{th(sr)}$ were assumed to be 0 for thick clouds (COT >16). $L_{3.9}^{sr}$ and $L_{3.9}^{th(sr)}$ 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. Characteristics of LUT for the retrieval of cloud microphysical parameters.

Input Variable	Value Range	Increment
Solar zenith angle	0–89°	5°
Viewing zenith angle	0–89°	5°
Relative azimuth angle	0–180°	30°
COT	0.5, 1, 2, 5, 8, 11, 15, 20, 30, 50, 70, 100	-
ER (μm)	Water cloud: 2, 4, 8, 16, 32 Ice cloud: 2, 4, 8, 16, 32, 64	-
Surface albedo	0–1.0	0.1
Surface temperature (K)	280–320	2
Cloud-top temperature	195–300	5
Cloud phase	Water, ice	-

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

2.2.4. 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 Section 2.2.2 and Section 2.2.3). CO₂ 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² 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 Table 1 and Table 2. The instantaneous SSI was estimated by linear interpolation from the lookup table once the above input data were known.

3. 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²), RMSE (in W/m²), and correlation coefficient (R²).

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. Overview of statistics comparing estimated SSI and SURFRAD measurements for the year of 2012.

Site	Latitude	Longitude	Clear Sky			Cloudy Sky			All Sky
			R2	BIAS W/m2 (%)	RMSE W/m2 (%)	R2	BIAS W/m2 (%)	RMSE W/m2 (%)	R2	BIAS W/m2 (%)	RMSE W/m2 (%)
BON	40.06°N	88.37°W	0.96	6.4 (1.3)	61.0 (12.6)	0.72	4.97 (2.4)	111.4 (53.6)	0.92	4.5 (1.2)	83.3 (21.7)
DRA	36.63°N	116.02°W	0.97	19.0 (3.6)	61.4 (11.8)	0.63	−54.5 (−18.6)	178.2 (60.9)	0.92	5.3 (1.1)	94.8 (19.8)
FPK	48.31°N	105.10°W	0.94	43.8 (10.9)	79.9 (19.9)	0.60	−0.1 (−0.04)	141.5 (65.2)	0.87	29.3 (8.6)	104.4 (30.7)
GWN	34.25°N	89.87°W	0.96	2.1 (0.4)	62.2 (12.1)	0.76	−2.2 (−0.9)	122.3 (47.9)	0.91	0.4 (0.1)	90.7 (22.1)
PSU	40.72°N	77.93°W	0.93	25.7 (5.7)	81.5 (18.2)	0.80	6.0 (2.8)	98.5 (45.2)	0.90	16.1 (4.8)	90.2 (26.9)
SXF	43.73°N	96.62°W	0.92	20.6 (4.3)	81.5 (17.1)	0.69	5.7 (2.7)	111.5 (53.8)	0.90	15.0 (4.0)	93.9 (25.1)
TBL	48.31°N	105.24°W	0.90	65.2 (13.6)	118.6 (24.7)	0.60	−28.7 (−11.1)	155.6 (60.0)	0.84	33.8 (8.3)	132.1 (32.5)
All			0.94	26.4 (5.5)	80.0 (16.8)	0.69	−5.9 (−2.6)	127.6 (55.1)	0.89	14.9 (3.8)	99.5 (25.5)

. The overall root mean square error (RMSE) values were 99.5 W/m² (25.5%), 80.0 W/m² (16.8%), and 127.6 W/m² (55.1%) for all-sky, clear-sky, and cloudy-sky conditions, respectively. The validation revealed a positive bias of 26.4 W/m² (5.5%) and a negative bias of −5.9 W/m² (−2.6%) for clear and cloudy skies. The RMSE values for all-sky ranged from 83.3 W/m² (21.7%) to 132.1 W/m² (32.5%), the RMSE values for clear skies ranged from 61.4 W/m² (11.8%) to 118.6 W/m² (24.7%), and the RMSE values for cloudy skies ranged from 98.5 W/m² (45.2%) to 141.5 W/m² (65.2%).

These statistics indicate that the quality of the retrieval was better for clear skies, which had a correlation coefficient (R²) 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.

4. Discussion

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

Table 4. Overview of error statistics for all-sky SSI for the year of 2008 (Zhang et al., 2014). RMSE: root mean square error.

Site	R2	BIAS (W/m2)	RMSE (W/m2)
BON	0.86	20	100
DRA	0.88	−55	119
FPK	0.82	5.5	111
GWN	0.92	1.7	86
PSU	0.87	12	100
SXF	0.86	14	102
TBL	0.77	−8.7	140

and

Table 5. Overview of error statistics for all-sky SSI derived from MODIS products for the years 2006–2008 (Qin et al., 2015).

Site	Clear Sky (W/m2)				All Sky (W/m2)
	Terra		Aqua		Terra		Aqua
	BIAS	RMSE	BIAS	RMSE	BIAS	RMSE	BIAS	RMSE
BON	11.5	41.1	15.3	54.4	4.0	86.3	7.6	95.0
DRA	−11.3	41.9	8.4	34.4	−11.0	55.0	8.1	69.7
FPK	20.2	43.8	29.9	49.0	−4.3	105.6	7.8	95.1
GWN	21.7	47.2	24.7	56.7	17.1	72.4	22.0	92.8
PSU	27.0	57.8	25.1	59.7	21.8	101.7	15.1	99.0
SXF	17.7	43.3	19.1	47.0	−5.3	101.0	−2.2	98.8
TBL	2.1	37.6	7.1	42.9	−17.7	113.6	−1.9	123.0

.

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² and 86 W/m², 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² at TBL, compared with −55 W/m² produced by Zhang’s method at DRA.

Qin et al. (2015) developed a physical parameterization to estimate SSI from MODIS atmospheric and land products and the retrievals were validated against in situ measurements at SURFRAD for three years (2006–2008) (Table 5) [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² 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² 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² (25.5%) is comparable to the MODIS-products-driven Breathing Earth System Simulator (BESS) shortwave products with an RMSE of 111.1 W/m² (22.6%) and 137.1 W/m² (31.7%) for temperate and continental climate zones, respectively [49].

As indicated by Yeom, a reduced RMSE of about 10 W/m² 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.

4.2. 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).

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², 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² and 90 W/m² 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². The SSI will change by about 548 and 600 W/m² 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.

An overall underestimation of SSI with a maximum bias of 62.5 W/m² and a minimum of 1 W/m² for water clouds is indicated in Figure 7 and

Table 6. Validation results at seven SURFRAD stations for water and ice clouds.

Site	Water Clouds				Ice Clouds				Mixed Clouds	Undetected Clouds
	R2	Bias W/m2 (%)	RMSE W/m2 (%)	NO.	R2	Bias W/m2 (%)	RMSE W/m2 (%)	NO.	NO.	NO.
BON	0.74	−8 (−3.5)	106 (49.5)	280	0.66	27 (19.8)	93 (68.7)	121	222	31
DRA	0.70	−63 (−21.0)	178 (59.6)	149	0.38	−44 (−15.7)	211 (75.3)	78	74	20
FPK	0.62	−25 (−10.0)	157 (63.7)	187	0.54	15 (6.6%)	150 (67.2)	133	220	34
GWN	0.78	−31 (−10.6)	136 (46.3)	326	0.75	45 (32.0)	101 (71.2)	130	181	44
PSU	0.81	−1 (−0.42)	105 (43.7)	576	0.81	44 (35.8)	83 (67.7)	83	217	20
SXF	0.80	−5 (−2.6)	80 (43.0)	233	0.49	29 (16.5)	135 (78.1)	109	243	47
TBL	0.64	−38 (−12.9)	169 (57.5)	197	0.61	−0.2 (−0.1)	118 (61.1)	186	173	20

. 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² at DRA, thereby exceeding the values for clear skies (61.4 W/m²) by a factor of three. The minimum RMSE value is 80.2 W/m² 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², 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.

5. 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² (16.8%), 127.6 W/m² (55.1%), and 99.5 W/m² (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.