Validation and Spatiotemporal Analysis of Ceres Surface Net Radiation Product

The Clouds and the Earth's Radiant Energy System (CERES) generates one of the few global satellite radiation products. The CERES ARM Validation Experiment (CAVE) has been providing long-term in situ observations for the validation of the CERES products. However, the number of these sites is low and their distribution is globally sparse, and particularly the surface net radiation product has not been rigorously validated yet. Therefore, additional validation efforts are highly required to determine the accuracy of the CERES radiation products. In this study, global land surface measurements were comprehensively collected for use in the validation of the CERES net radiation (R n) product on a daily (340 sites) and a monthly (260 sites) basis, respectively. The validation results demonstrated that the CERES R n product was, overall, highly accurate. The daily validations had The accuracy was slightly lower for the high latitudes. Following the validation, the monthly CERES R n product, from March 2000 to July 2014, was used for a further analysis. The global spatiotemporal variation of the R n , which occurred during the measurement period, was analyzed. In addition, two hot spot regions, the southern Great Plains and south-central Africa, were then selected for use in determining the driving factors or attribution of the R n variation. We determined that R n over the southern Great Plains decreased by´0.33 W¨m´2 per year, which was mainly driven by changes in surface green vegetation and precipitation. In south-central Africa, R n decreased at a rate of´0.63 W¨m´2 per year, the major driving factor of which was surface green vegetation.


Introduction
The surface radiation budget (SRB) plays an important role in hydrological and biological modeling, and agricultural science and its applications [1,2].The SRB components consist of the downward shortwave irradiance (DSW), upward shortwave irradiance (USW), downward longwave flux (DLW), upward longwave flux (ULW), and net all-wave radiation (R n ).The surface R n is the sum of the four formerly-mentioned radiation components, and can be represented mathematically: where R s Ó denotes the DSW (W¨m ´2), R s Ò denotes the USW (W¨m ´2), R l Ó denotes the DLW (W¨m ´2), and R l Ò denotes the ULW (W¨m ´2).
The surface R n can be obtained by ground observation stations, general circulation model (GCM) simulations, and remote sensing retrieval [3].Comparatively, in situ observations have the highest accuracy; however, they are very expensive to measure at ground sites.Additionally, these ground sites are sparsely distributed, globally, and cannot meet the requirements of various applications.Data simulated by GCMs usually have a coarse spatial resolution and low accuracy [2].Remote sensing retrieval has become one of the most important methods for obtaining the radiation data [4].However, clouds make it impossible to retrieve surface thermal radiation components directly from satellite observations, and the alternative solution is to estimate the all-wave net radiation from the satellite shortwave net radiation product, in conjunction with other information [5,6].
The CERES surface net radiation data are calculated using the radiative transfer model, with the atmospheric and surface properties obtained from a variety of sources used as the model inputs, but constrained by both solar-reflected and Earth-emitted radiation from the top of the atmosphere (TOA) for all-sky and clear-sky conditions [7].Charlock et al. [8] conducted the CERES/ARM/GEWEX Experiment (CAGEX) before the launch of satellites, in order to validate the Fu-Liou radiative transfer model [9].They found that the modeled longwave radiation agrees well with the in situ observations, while the shortwave radiation exhibits a large difference between the modeled and in situ data.Since then, CAGEX has been developed to become the CERES ARM Validation Experiment (CAVE) [10], which provides validation data and model parameters, which are measured on the ground for each version of the CERES radiation product.It includes 84 sites, with limited spatial distribution [11].David Rutan et al. [12] utilized 40 sites in order to validate the Clouds and Radiative Swath (CRS), and found that the RMSEs were related to the land cover type.Gui et al. [13] validated the CERES-Monthly Gridded Radiative Fluxes and Clouds (FSW) in the Tibetan Plateau, and found that the CERES had a lower accuracy for shortwave radiation.Wang et al. [14] also found that FLASHFlux/CERES (Fast Longwave and Shortwave Fluxes Time Interpolated and Spatially Averaged/CERES) satellite products over rugged terrains showed absolute errors of about 31.4W¨m ´2 (12.3%).Pan et al. [15] assessed the CERES R n product by using 50 ground sites in China, and showed that 56% of the net radiation errors come from the net shortwave radiation.Kratz et al. [16] validated the CERES_Ed2b algorithms, and found that the shortwave radiation had systematic errors of approximately 14-21 W¨m ´2, which does not meet the meteorological research.In addition, the CERES Ocean Validation Experiment (COVE) [17] has been collecting data for the validation of CERES on the sea using one site.Charlock et al. [18] validated the CERES data on the sea in 2001, using the Chesapeake Lighthouse and Aircraft Measurements for Satellites (CLAMS) and COVE data, and the results were good.However, all of these previous studies into the validation of the CERES had a low number and distribution of sites, and the accuracies of the CERES R n around the global land surface are still generally unknown.Therefore, the CERES surface net radiation product needs to be validated globally.
Wild et al. [19] obtained a global R n mean value from 2000 to 2004 by calculating the best estimation of shortwave and longwave radiation over the land surface and oceans utilizing 43 Coupled Model Intercomparison Project Phase 5 (CMIP5) model datasets.Sobrino et al. [20] researched the DSW and R n around the globe from 1980 to 2010 by using National Centers for Environmental Prediction (NCEP) and the National Center of Atmospheric Research (NCAR) re-analysis data and detected a decrease in these values in selected areas.Gao et al. [21] detected a decrease in R n over China by using the 1961-2010 data of 699 meteorological stations.Shi et al. [22] fused several radiation data and analyzed the spatiotemporal variation of SRB over the Tibetan Plateau in the past 30 years.However, these previous studies on R n variation were based mainly on data with high uncertainties, and some research focused only on partial areas.Therefore, global spatiotemporal analysis of R n utilizing data with high accuracy is urgently needed.
In this study, the measurements obtained from observation sites around the world were collected to be used for validating the CERES_SYN1deg_Ed3A R n product under all-sky conditions on a daily (340 sites) and a monthly (260 sites) basis, respectively.Considering that most R n variation studies were mainly developed in limited areas and the chosen radiation data had high uncertainties, the CERES monthly R n product, from March 2000 to July 2014, was used for a further spatiotemporal analysis across the global land surface.As the southern Great Plains and south-central Africa are located in special and vital geographical locations, they were chosen as hot spots for a detailed analysis of the relationship between the variations in surface R n and their driving factors.
This paper proceeds as follows: the CERES, in situ observations, driving factor datasets, and experiment methods are introduced in Section 2; the validation and spatiotemporal analysis results and discussions are shown in Section 3; and the conclusions are provided in Section 4.

The CERES Data
The CERES_SYN1deg_Ed3A product is used in this study.The "SYN" (Synoptic Radiative Fluxes and Clouds) means that this version provides radiation data on clear and all sky conditions, the "1deg" means it has a 1 degree spatial resolution, and the "Ed3A" is the version number [23].In addition to the radiation flux, this version also provides MODIS cloud properties and aerosols information.As the three-hourly sampled geostationary (GEO) radiances and cloud properties were used to model the diurnal variability between observations, this version achieves more accurate retrieval results than previous versions [24].The CERES consists of three-hourly, daily, and monthly data, released on Greenwich Mean Time (GMT) [25].The data we used ranged from 1 March 2000 to 30 July 2014.In this study, three-hourly and monthly DSW, USW, DLW, and ULW products, in all-sky conditions, were selected to calculate the R n by us to be used in the validation.
Remote Sens. 2016, 8, 90 3/19 analysis across the global land surface.As the southern Great Plains and south-central Africa are located in special and vital geographical locations, they were chosen as hot spots for a detailed analysis of the relationship between the variations in surface Rn and their driving factors.
This paper proceeds as follows: the CERES, in situ observations, driving factor datasets, and experiment methods are introduced in Section 2; the validation and spatiotemporal analysis results and discussions are shown in Section 3; and the conclusions are provided in section 4.

The CERES Data
The CERES_SYN1deg_Ed3A product is used in this study.The "SYN" (Synoptic Radiative Fluxes and Clouds) means that this version provides radiation data on clear and all sky conditions, the "1deg" means it has a 1 degree spatial resolution, and the "Ed3A" is the version number [23].In addition to the radiation flux, this version also provides MODIS cloud properties and aerosols information.As the three-hourly sampled geostationary (GEO) radiances and cloud properties were used to model the diurnal variability between observations, this version achieves more accurate retrieval results than previous versions [24].The CERES consists of three-hourly, daily, and monthly data, released on Greenwich Mean Time (GMT) [25].The data we used ranged from 1 March 2000 to 30 July 2014.In this study, three-hourly and monthly DSW, USW, DLW, and ULW products, in allsky conditions, were selected to calculate the Rn by us to be used in the validation.

In Situ Observation Data
Observations taken from 15 measurement networks, shown in Figure 1, were collected, and the metadata for these measurement networks are given in Table 1 [26][27][28][29][30][31][32][33][34][35][36].The sites are located around the Earth's surface, representing different land cover types, such as forests (135 sites), scrublands (20 sites), grasslands (76 sites), bare lands (four sites), wetlands (18 sites), tundra (two sites), ice (23 sites), croplands (61 sites), and urban areas (one site).Some ground sites were eliminated from the monthly validation because their monthly in situ observations were less than 12.   Data from all those networks were checked site by site, and some extremely abnormal values were removed.Most networks had long time range observations, whereas CEOP and HiWATER in China had shorter time ranges and discontinuous observation data.All are trustworthy based on network control and continuous project support and have been utilized in many applications.

Land Cover Type Data
Given that the surface R n is closely related to the land cover type, MCD12Q1 [37] land cover type data were used in this study, as shown in Figure 2 for reference.Data from all those networks were checked site by site, and some extremely abnormal values were removed.Most networks had long time range observations, whereas CEOP and HiWATER in China had shorter time ranges and discontinuous observation data.All are trustworthy based on network control and continuous project support and have been utilized in many applications.

Land Cover Type Data
Given that the surface Rn is closely related to the land cover type, MCD12Q1 [37] land cover type data were used in this study, as shown in Figure 2 for reference.

Data for Attribution Analysis
As the sum of the four surface radiation components, Rn can be synthetically affected by many atmospheric and land surface conditions.Therefore, it is necessary to select multiple driving factors for further analysis in order to provide evidence to support the explanations for the variations in Rn.Normalized Difference Vegetation Index (NDVI), which characterizes the vegetation growth and coverage states, affected surface albedo directly [38].Additionally, vegetation also influences surface

Data for Attribution Analysis
As the sum of the four surface radiation components, R n can be synthetically affected by many atmospheric and land surface conditions.Therefore, it is necessary to select multiple driving factors for further analysis in order to provide evidence to support the explanations for the variations in R n .Normalized Difference Vegetation Index (NDVI), which characterizes the vegetation growth and coverage states, affected surface albedo directly [38].Additionally, vegetation also influences surface temperature and diurnal temperature range, and further surface longwave radiation [39].The snow cover affects surface albedo that regulates USW.Albedo largely determines surface shortwave net radiation.Solar radiation significantly affects the Tm that determines the ULW directly as a feedback according to Stefan-Boltzmann law.As for atmospheric factors, both precipitation (Pre) and cloud cover (CC) affect the DSW and DLW [40].In addition, Pre greatly affects vegetation growth.So, all of the chosen factors are helpful to analyze the R n variations.The characteristics of these data are listed in Table 2. CRU datasets are based on interpolated station observations.Harris et al. [47] found that each variable of the CRU datasets agree tightly with the other published datasets at regional and global scales, such as mean temperature data of the University of Delaware (UDEL), precipitation data of the Global Precipitation Climatology Centre (GPCC), DTR from Vose et al. [48].The GLASS products, which are spatiotemporally continuous, have long temporal ranges and higher spatial resolutions than the similar satellite products.They are produced from multiple satellite data but mainly from MODIS data for the period of 2000-2013 [42].

Validation
The three-hourly surface R n values were calculated from the four radiation components using Equation (1).As the CERES data were released in GMT, while ground sites were in local time, the CERES sample data were converted into local time formats and averaged to the daily values for use in the for daily validations.The CERES monthly sample data were used for the monthly validations.
Due to the fact that different measurement networks have different temporal resolutions and units, the observations were handled as follows: for calculating daily in situ observations, measurement data must be recorded in every hour of one day; for obtaining the monthly in situ observations, at least six daily data should be recorded every 10 days in one month (the daily data should be recorded for more than 18 days in February).If the data did not meet these conditions, they would be eliminated before the validation.We unified the units of measurement to be in W¨m ´2 and calculated the average for the daily/monthly in situ observations.This method is trustworthy and controls the data quality strictly.Data recorded at sites that did not have at least 12 samples were also eliminated.In total, there were 340 sites used for the daily validations and 260 sites used for the monthly validations.
In order to express the land cover complexity in a better way, land cover data were re-sampled from 0.5 km to 0.05 ˝.The normalized square deviation of the merging pixels was then taken to be the new pixel value.The new dataset was called the "land cover complexity".This dataset was the base map on which the ground sites were displayed, and used as a reference for comparing the land surface complexity.In additionally, we also made a complexity data set with a 1 ˝resolution, matching the CERES data.The code names for the International Geosphere-Biosphere Program (IGBP) land cover types are shown in Table 3.The Mean Bias Error (MBE), Root-Mean-Square Error (RMSE), and coefficient of determination (R 2 ) were selected to measure the validation results.

Spatiotemporal Analysis
The monthly CERES R n data were chosen for the spatiotemporal analysis, but also aggregated to seasonal and annual scales.The driving factor datasets have different spatial and temporal resolutions, which are higher than those of the CERES.Before beginning the quantitative analysis, we re-sampled all of the datasets at a resolution of 1 ˝, using the following principles.If all the merging pixels have no null values, average them directly.In the case where the number of null merging pixels is more than 60% of the total, this new pixel value is judged a null value.If the number of null pixels is less than 60%, take the average of the remaining non-null values as the new pixel value.All the datasets are re-sampled to obtain monthly mean values (the NDVI re-sampling is based on the maximum pixel value of each month).If the time range of the data recorded in a month is less than 20 days, the data in that month are excluded.
The effects of the driving factors on the R n were analyzed using general linear models (GLMs) [49].First, we obtained the correlation coefficient (R) between the driving factors and the R n .Then, the driving factors entering the model were selected by a forward stepwise procedure by correlation order [50].When adding a new independent variable or driving factor into the GLMs, the explained variance of the new model increases.We subtracted the previous explained variance to calculate the percentage sum of squares explained (%SS).However, the result cannot get rid of factors' collinearity.The %SS of dependent driving factors may have a larger difference as the latter factor provides little information for the GLMs [51].Considering the limitations of the %SS, the partial least squares (PLS) regression [52] was developed and combined with the variable importance in the projection (VIP) scores to be used in a further analysis.The aim of this further analysis was to characterize the relative importance of these driving factors, in terms of their influence on the R n dynamics.All of the independent variables are projected into principal components, and the VIP score of independent variable j can be represented mathematically: where k is the number of independent variable; c h is the one of m principal components; r(Y, c h ) is the correlation index between c h and dependent variable Y; w hj is the weight of independent variable j for c h .Thus, the PLS regression is the combination of multivariable linear regression modeling, correlation analysis, and principal components analysis, showing its strength when the independent variables are interrelated and when there are fewer data samples available [53].correlation analysis, and principal components analysis, showing its strength when the independent variables are interrelated and when there are fewer data samples available [53].The validation results of the CERES Rn, at each site, are shown in Figure 4.The grey scale represents the land cover complexity, to be used as a reference.

Validation of the CERES Rn Product
The MBE and RMSE values, at a number of sites in China, are much larger than most sites in Europe and America.This may be related to the measurement accuracy of different networks, as, for example, the La Thuile and the ARM networks have higher accuracy requirements.Compared with other networks, the CEOP, HiWATER, and CEOPInt, built in China, have poorer data continuity and a lower quantity of data.The CERES accuracy in Greenland appears to be lower.Inamdar et al., also found the CERES net shortwave radiation has large errors in the underlying surface snow [54].The CERES Rn showed bad performance on a site in the central southern region of the United States.By checking the 0.05° land cover complexity results, we determined that this site is located on a more complicated underlying surface.Considering that the CERES Rn showed strong agreement with other sites in this region, the accuracy is still considered trustworthy.The validation results of the CERES R n , at each site, are shown in Figure 4.The grey scale represents the land cover complexity, to be used as a reference.
The MBE and RMSE values, at a number of sites in China, are much larger than most sites in Europe and America.This may be related to the measurement accuracy of different networks, as, for example, the La Thuile and the ARM networks have higher accuracy requirements.Compared with other networks, the CEOP, HiWATER, and CEOPInt, built in China, have poorer data continuity and a lower quantity of data.The CERES accuracy in Greenland appears to be lower.Inamdar et al., also found the CERES net shortwave radiation has large errors in the underlying surface snow [54].The CERES R n showed bad performance on a site in the central southern region of the United States.By checking the 0.05 ˝land cover complexity results, we determined that this site is located on a more complicated underlying surface.Considering that the CERES R n showed strong agreement with other sites in this region, the accuracy is still considered trustworthy.

Monthly Scale
The monthly CERES R n samples were compared with the observations taken from 260 sites.The resulting density of scatter diagram is shown in Figure 5.The validation results, at each ground site, are displayed in Figure 6.
Compared with the daily validation results, the performance of the monthly product is much better.It should be noted that there appears to be a poorer precision at a number of sites near the coast, and in Greenland.This reduction in precision may be caused by the large pixel areas with higher land cover complexity.Both daily validation and monthly validation showed lower accuracy of CERES Rn over Greenland.
The preprocessing of in situ data and CERES data was very strict, and the results are trustworthy.The validation results demonstrate that the CERES Rn product is, overall, highly accurate, although the accuracy is slightly lower for high latitudes.
(a) At a monthly scale, the MBE was 3.40 W¨m ´2, the RMSE was 25.57W¨m ´2, and the R 2 was 0.84.This indicates that the CERES R n data are slightly larger than the in situ observations, regardless of noise.
The validation results, at each ground site, are displayed in Figure 6.
Compared with the daily validation results, the performance of the monthly product is much better.It should be noted that there appears to be a poorer precision at a number of sites near the coast, and in Greenland.This reduction in precision may be caused by the large pixel areas with higher land cover complexity.Both daily validation and monthly validation showed lower accuracy of CERES R n over Greenland.
The preprocessing of in situ data and CERES data was very strict, and the results are trustworthy.The validation results demonstrate that the CERES R n product is, overall, highly accurate, although the accuracy is slightly lower for high latitudes.The validation results, at each ground site, are displayed in Figure 6.
Compared with the daily validation results, the performance of the monthly product is much better.It should be noted that there appears to be a poorer precision at a number of sites near the coast, and in Greenland.This reduction in precision may be caused by the large pixel areas with higher land cover complexity.Both daily validation and monthly validation showed lower accuracy of CERES Rn over Greenland.
The preprocessing of in situ data and CERES data was very strict, and the results are trustworthy.The validation results demonstrate that the CERES Rn product is, overall, highly accurate, although the accuracy is slightly lower for high latitudes. (a)

Global Analysis of Rn
The monthly CERES Rn data, from 2001 to 2013, were used to detect the mean values of the Rn and annual trends in the Rn across the land surface, as showed in Figure 7. Data from years 2000 and 2014 covered less than one full year and were used only for the seasonal analysis.
There are entirely different distributions and trends of the Rn across the land surface over the past 13 years.It is evident that Rn, in some regions, has been increasing, such as in Greenland, the African tropical rainforest regions, the southeast of Brazil, the coast of Antarctica, and the west of Australia.Over the southwest of Greenland, Rn increased the fastest, at a rate of 4.21 W•m −2 per year.Conversely, there are some regions where Rn has decreased over time, such as in the southeast of Argentina, the African savannah, the southern Great Plains, and especially on the eastern side of the Andes Mountains, where Rn has decreased at a rate of 5.10 W•m −2 per year.

Global Analysis of R n
The monthly CERES R n data, from 2001 to 2013, were used to detect the mean values of the R n and annual trends in the R n across the land surface, as showed in Figure 7. Data from years 2000 and 2014 covered less than one full year and were used only for the seasonal analysis.
There are entirely different distributions and trends of the R n across the land surface over the past 13 years.It is evident that R n , in some regions, has been increasing, such as in Greenland, the African tropical rainforest regions, the southeast of Brazil, the coast of Antarctica, and the west of Australia.Over the southwest of Greenland, R n increased the fastest, at a rate of 4.21 W¨m ´2 per year.Conversely, there are some regions where R n has decreased over time, such as in the southeast of Argentina, the African savannah, the southern Great Plains, and especially on the eastern side of the Andes Mountains, where R n has decreased at a rate of 5.10 W¨m ´2 per year.
11/19 In order to understand what caused these significant changes, we selected two geographical hotspot regions (the southern Great Plains and south-central Africa, where significant changes happened, the CERES Rn product and ancillary data for attribution analysis are relatively accurate), marked by the thick lines in Figure 7, for further analysis.

Southern Great Plains
The Great Plains is a well-known geographic area of the United States, and its southern chosen region covers parts of Colorado, Kansas, New Mexico, Oklahoma, and Texas [55].This region is a major cotton and wheat production area of the U.S., and animal husbandry is another important industry [56].The distribution of surface Rn influences its continental arid climate, which determines the agricultural yield and husbandry production in this region.Moreover, this region includes Tornado Alley [57]; the Rn variation affects the wind speed, which causes sand-dust disasters and tornadoes.The southern Great Plains is currently experiencing a remarkable decrease in Rn.
The temporal trend of the Rn over the Southern Great Plains was determined by a linear regression.The annual and seasonal trend plots were generated, and shown in Figure 8.The four seasons are conventionally defined as spring (March, April, and May), summer (June, July, and August), autumn (September, October, and November), and winter (December, January, and February).In order to understand what caused these significant changes, we selected two geographical hotspot regions (the southern Great Plains and south-central Africa, where significant changes happened, the CERES R n product and ancillary data for attribution analysis are relatively accurate), marked by the thick lines in Figure 7, for further analysis.

Southern Great Plains
The Great Plains is a well-known geographic area of the United States, and its southern chosen region covers parts of Colorado, Kansas, New Mexico, Oklahoma, and Texas [55].This region is a major cotton and wheat production area of the U.S., and animal husbandry is another important industry [56].The distribution of surface R n influences its continental arid climate, which determines the agricultural yield and husbandry production in this region.Moreover, this region includes Tornado Alley [57]; the R n variation affects the wind speed, which causes sand-dust disasters and tornadoes.The southern Great Plains is currently experiencing a remarkable decrease in R n .
The temporal trend of the R n over the Southern Great Plains was determined by a linear regression.The annual and seasonal trend plots were generated, and shown in Figure 8.The four seasons are conventionally defined as spring (March, April, and May), summer (June, July, and August), autumn (September, October, and November), and winter (December, January, and February).In the significantly changing areas, the correlation analysis was carried out and %SS values were calculated between the CERES Rn and the driving factors, as shown in Table 4.
Table 4. Statistics for the correlation coefficient (R) and percentage sum of squares explained (%SS) between the Rn and the driving factors for the southern Great Plains.Numbers marked with a * indicates p-values < 0.10, and ** indicates p-values < 0.05, in which case the results are deemed to be significant.R stands for the correlation coefficient, Rn for the net radiation, Pre for precipitation, CC for cloud cover, TΔ for temperature range, Tm for surface temperature, and SC for snow cover.The R n decreased over the entire southern Great Plains by ´0.33 W¨m ´2 per year.The fastest decrease occurred in summer, with a gradient of ´0.58 W¨m ´2 per year, while the slowest decrease in spring was only ´0.09 W¨m ´2 per year.Furthermore, the R n decreased in the central region, and increased in the west edge.In the summer, the central part decreased sharply at a rate of ´2.41 W¨m ´2 per year.However, also in the summer, the R n in the northwest increased with a gradient of 1.62 W¨m ´2 per year.R n in the east flat area did not show any significant trend.

Annual
In the significantly changing areas, the correlation analysis was carried out and %SS values were calculated between the CERES R n and the driving factors, as shown in Table 4.
Table 4. Statistics for the correlation coefficient (R) and percentage sum of squares explained (%SS) between the R n and the driving factors for the southern Great Plains.Numbers marked with a * indicates p-values < 0.10, and ** indicates p-values < 0.05, in which case the results are deemed to be significant.R stands for the correlation coefficient, R n for the net radiation, Pre for precipitation, CC for cloud cover, T∆ for temperature range, Tm for surface temperature, and SC for snow cover.At an annual scale, NDVI explained 51% of the R n variation, and CC explained 16%.In the spring, albedo explained 34% of the R n variation.The lower %SS of the NDVI and SC may be caused by the co-linearity in the GLMs.In summer, CC and Pre played an important role, accounting for 61% of the R n variation together.In autumn, NDVI explained 50% of the R n variation.In winter, SC explained 62% of R n variation.

Annual
However, due to the relevant dependency of certain variables within the GLMs (for example, albedo and SC in the winter), some factors cannot be used individually in an explanatory capacity.Additionally, the majority of correlation coefficients were not significant.Therefore, we conducted a further analysis utilizing the PLS regression to deal with the problems brought on by the interrelated nature of these driving factors.The VIP scores of these driving factors are shown in Figure 9.
Remote Sens. 2016, 8, 90 13/19 At an annual scale, NDVI explained 51% of the Rn variation, and CC explained 16%.In the spring, albedo explained 34% of the Rn variation.The lower %SS of the NDVI and SC may be caused by the co-linearity in the GLMs.In summer, CC and Pre played an important role, accounting for 61% of the Rn variation together.In autumn, NDVI explained 50% of the Rn variation.In winter, SC explained 62% of Rn variation.
However, due to the relevant dependency of certain variables within the GLMs (for example, albedo and SC in the winter), some factors cannot be used individually in an explanatory capacity.Additionally, the majority of correlation coefficients were not significant.Therefore, we conducted a further analysis utilizing the PLS regression to deal with the problems brought on by the interrelated nature of these driving factors.The VIP scores of these driving factors are shown in Figure 9. Considering the %SS and VIP scores, it is clear that NDVI and TΔ were more dominating at an annual scale.Pre and CC effects also performed noticeably well.We think that the decreasing NDVI increased the albedo in the southern Great Plains, which affected the Rn in shortwave net radiation, and declining vegetation caused the larger TΔ and negative correlation with Rn, which affected the longwave net radiation.CC and precipitation affected the DSW and the decreasing precipitation or drought might be one of the reasons that influenced the NDVI.In the spring, precipitation, SC, and albedo were the key factors.SC changed the albedo and affected the USW so they were negatively correlated with Rn, and precipitation in the spring often was in the form of snow, further affecting the USW, especially of the western higher altitude areas.
In the summer, many factors, such as CC, TΔ, Tm, albedo, and NDVI, have large VIP scores.The CC, which reduces the DSW, showed positive correlation with Rn notably.We think that DSW was ample in the summer, as the NDVI decreased and albedo went up, more shortwave radiation was lost, longwave radiation played an important role in Rn variation.The Tm performance provided the proof to support this viewpoint, while CC might just decrease matched with Rn variation due to little rainfall, and fewer CC decreased the DLW especially in the night, and increased the TΔ combined with fewer NDVI due to the heat preservation [58].
In autumn, the key factors were the NDVI and albedo.SC also became an active factor.NDVI and SC both affected the Rn by changing the albedo, so land cover types played a more vital role and atmospheric factors had limited influence.
In the winter, the SC and NDVI had the highest VIP scores.It is clear that snowfall increased the albedo and NDVI was easily affected by the SC and vegetation growth periodicity in the winter, so these two factors both affected the albedo, thus changing the Rn.
Overall, we deduced that the NDVI was the main driving factor in the southern Great Plains, and precipitation played an indispensable role.Considering the %SS and VIP scores, it is clear that NDVI and T∆ were more dominating at an annual scale.Pre and CC effects also performed noticeably well.We think that the decreasing NDVI increased the albedo in the southern Great Plains, which affected the R n in shortwave net radiation, and declining vegetation caused the larger T∆ and negative correlation with R n , which affected the longwave net radiation.CC and precipitation affected the DSW and the decreasing precipitation or drought might be one of the reasons that influenced the NDVI.In the spring, precipitation, SC, and albedo were the key factors.SC changed the albedo and affected the USW so they were negatively correlated with R n , and precipitation in the spring often was in the form of snow, further affecting the USW, especially of the western higher altitude areas.
In the summer, many factors, such as CC, T∆, Tm, albedo, and NDVI, have large VIP scores.The CC, which reduces the DSW, showed positive correlation with R n notably.We think that DSW was ample in the summer, as the NDVI decreased and albedo went up, more shortwave radiation was lost, longwave radiation played an important role in R n variation.The Tm performance provided the proof to support this viewpoint, while CC might just decrease matched with R n variation due to little rainfall, and fewer CC decreased the DLW especially in the night, and increased the T∆ combined with fewer NDVI due to the heat preservation [58].
In autumn, the key factors were the NDVI and albedo.SC also became an active factor.NDVI and SC both affected the R n by changing the albedo, so land cover types played a more vital role and atmospheric factors had limited influence.
In the winter, the SC and NDVI had the highest VIP scores.It is clear that snowfall increased the albedo and NDVI was easily affected by the SC and vegetation growth periodicity in the winter, so these two factors both affected the albedo, thus changing the R n .
Overall, we deduced that the NDVI was the main driving factor in the southern Great Plains, and precipitation played an indispensable role.

South-Central Africa
The second hot spot is situated in south-central Africa, including a majority of the Congo basin and part of the South African Veld.Abundant creatures and mineral resources exist in south-central Africa [59,60].Countries located here are dominated by agriculture and mining industry.Therefore, understanding the R n variation can help guide the population towards developing agriculture in a better way.Abundant flora and fauna exist in south-central Africa; the R n variation could be an indicator of climate change that determines the habitat environment [61].Therefore, research and protection of the biodiversity in this region is crucial.It is noteworthy that a clear decrease in R n was observed.As this area has a tropical savanna climate, we performed the R n variation analysis that is shown in Figure 10, separately, during the dry season (from May to October) and wet season (from the previous November to April).The second hot spot is situated in south-central Africa, including a majority of the Congo basin and part of the South African Veld.Abundant creatures and mineral resources exist in south-central Africa [59,60].Countries located here are dominated by agriculture and mining industry.Therefore, understanding the Rn variation can help guide the population towards developing agriculture in a better way.Abundant flora and fauna exist in south-central Africa; the Rn variation could be an indicator of climate change that determines the habitat environment [61].Therefore, research and protection of the biodiversity in this region is crucial.It is noteworthy that a clear decrease in Rn was observed.As this area has a tropical savanna climate, we performed the Rn variation analysis that is shown in Figure 10, separately, during the dry season (from May to October) and wet season (from the previous November to April).In the annual analysis, the Rn decreased at a rate of −0.63 W•m −2 per year, overall.The dry season decreased at a slower rate of −0.39 W•m −2 per year, whereas the wet season decreased at a faster rate of −0.82 W•m −2 per year.With the exception of the western and eastern coasts, the Rn for the majority of south-central Africa was observed to decrease.The correlation analysis and %SS were calculated, and are shown in Table 5. South-central Africa had very little snow cover, and therefore the SC factor was excluded from the analysis.It is clear that the NDVI was the dominant driving factor.It had the largest explanatory power, accounting for more than 50% of the Rn variation for both the annual and dry season analyses.The albedo showed higher correlation with Rn but limited %SS.The results likely occurred because In the annual analysis, the R n decreased at a rate of ´0.63 W¨m ´2 per year, overall.The dry season decreased at a slower rate of ´0.39 W¨m ´2 per year, whereas the wet season decreased at a faster rate of ´0.82 W¨m ´2 per year.With the exception of the western and eastern coasts, the R n for the majority of south-central Africa was observed to decrease.The correlation analysis and %SS were calculated, and are shown in Table 5. South-central Africa had very little snow cover, and therefore the SC factor was excluded from the analysis.It is clear that the NDVI was the dominant driving factor.It had the largest explanatory power, accounting for more than 50% of the R n variation for both the annual and dry season analyses.The albedo showed higher correlation with R n but limited %SS.The results likely occurred because albedo does not provide more information than NDVI owing to the collinearity with NDVI.Tm showed a positive correlation with R n ; it is widely accepted that Tm determines the ULW but is affected by the shortwave radiation.We conclude that shortwave radiation was the major part of the surface radiation budget here in the wet season.Considering the limitation mentioned previously, the VIP scores of the driving factors in south-central Africa are shown in Figure 11.
Remote Sens. 2016, 8, 90 albedo does not provide more information than NDVI owing to the collinearity with NDVI.Tm showed a positive correlation with Rn; it is widely accepted that Tm determines the ULW but is affected by the shortwave radiation.We conclude that shortwave radiation was the major part of the surface radiation budget here in the wet season.Considering the limitation mentioned previously, the VIP scores of the driving factors in south-central Africa are shown in Figure 11.The VIP scores corresponded very well with the %SS.No matter what period it is, the NDVI is the major driving factor, and we thought that the albedo, here, was also predominantly affected by this factor, as we infered from %SS.Additionally, the Tm had a significant impact during the wet season compared with the dry season.The Pre and CC had limited impact on the variation of the Rn during the annual analysis and dry season.CC in the wet season got higher VIP scores; we inferred that the increase of CC resisted more DSW, which controlled the surface radiation budget as we discussed.CC also influenced the TΔ and they got similar performance in all time, and precipitation had little variation in all time period so it had limited impact on Rn variation.It was hypothesized that the variation of Rn was perhaps caused by the degradation of forests and grassland due to drought in the dry season and human activities.
However, several uncertainties should be discussed.The accuracy of the atmospheric factors obtained from the CRU in Africa is not as high as the data in the southern Great Plains owing to the lower number of sites.Harris et al. [47] compared the CRU mean temperature data with data from UDEL in southern Africa and detected highly similar trends with a correlation coefficient of 0.91.Precipitation data were also compared with GPCC data with a correlation index of 0.89.Therefore, the CRU data are still trustworthy in this region.The reasons behind the high number of non-significant correlation coefficients were analyzed and were outlined here: (1) the R was calculated over areas with distinctly changing Rn, which includes different land cover types.Different underlying surfaces and elevation may lead to the different, or even contradictory, relationships between the Rn and the driving factors; (2) the time range of the CERES is short; and (3) the Rn can be influenced by many driving factors, and the Rn, itself, can influence the land surface or sky conditions though negative feedbacks.

Conclusions
There is a vast potential for developing remote sensing for use in determining the SRB, and the CERES product has become one of the most representative products.However, with a low number of in situ sites and poor distribution of CAVE sites, additional validation efforts are still required in order to determine the global accuracy of the CERES radiation products.This study validated the The VIP scores corresponded very well with the %SS.No matter what period it is, the NDVI is the major driving factor, and we thought that the albedo, here, was also predominantly affected by this factor, as we infered from %SS.Additionally, the Tm had a significant impact during the wet season compared with the dry season.The Pre and CC had limited impact on the variation of the R n during the annual analysis and dry season.CC in the wet season got higher VIP scores; we inferred that the increase of CC resisted more DSW, which controlled the surface radiation budget as we discussed.CC also influenced the T∆ and they got similar performance in all time, and precipitation had little variation in all time period so it had limited impact on R n variation.It was hypothesized that the variation of R n was perhaps caused by the degradation of forests and grassland due to drought in the dry season and human activities.
However, several uncertainties should be discussed.The accuracy of the atmospheric factors obtained from the CRU in Africa is not as high as the data in the southern Great Plains owing to the lower number of sites.Harris et al. [47] compared the CRU mean temperature data with data from UDEL in southern Africa and detected highly similar trends with a correlation coefficient of 0.91.Precipitation data were also compared with GPCC data with a correlation index of 0.89.Therefore, the CRU data are still trustworthy in this region.The reasons behind the high number of non-significant correlation coefficients were analyzed and were outlined here: (1) the R was calculated over areas with distinctly changing R n , which includes different land cover types.Different underlying surfaces and elevation may lead to the different, or even contradictory, relationships between the R n and the driving factors; (2) the time range of the CERES is short; and (3) the R n can be influenced by many driving factors, and the R n , itself, can influence the land surface or sky conditions though negative feedbacks.

Conclusions
There is a vast potential for developing remote sensing for use in determining the SRB, and the CERES product has become one of the most representative products.However, with a low number of in situ sites and poor distribution of CAVE sites, additional validation efforts are still required in order to determine the global accuracy of the CERES radiation products.This study validated the CERES R n product close to the earth's land surface at daily and monthly basis, and analyzed its spatiotemporal variation from 2001 to 2013.Two hot spot regions (the southern Great Plains and south-central Africa) were then selected for use in determining the driving factors.The conclusions of this study are summarized in the following points: 1.
This paper validated the CERES_SYN1deg_Ed3A all-sky surface R n at a daily and monthly basis, respectively.The daily validations had an MBE of 3.43 W¨m ´2, RMSE of 33.56 W¨m ´2, and R 2 of 0.79.The monthly validations had an MBE of 3.40 W¨m ´2, RMSE of 25.57W¨m ´2, and R 2 of 0.84.Considering that the CERES_Ed2b was previously found to have a systematic bias of 14~21 W¨m ´2 [16], the integral accuracy of the CERES_SYN1deg_Ed3A is high enough for use in scientific research.

2.
The accuracy distribution of the daily and monthly CERES R n products were calculated.It was found that areas in the middle and low latitudes had a higher accuracy, and that the R 2 was lower for areas in the high latitudes, such as Greenland.This may be caused by the lower accuracy of CERES inversion over the snow surface.The uncertainty near the coast was caused by the coarse resolution of the CERES products matching with in situ observations.3.
This study analyzed the spatiotemporal variation of the monthly CERES R n , from 2001 to 2013, and found that different regions exhibited vastly differing rates of change.In terms of the hot spot analysis, the R n decreased slightly over the entire southern Great Plains, with a gradient of ´0.33 W¨m ´2 per year.As for the attribution analysis, percentage sum of squares explained and VIP scores were calculated.From the annual and seasonal analysis, we concluded that NDVI, which represents the surface vegetation condition, was the major driving factor of the R n variation over the southern Great Plains and that precipitation played an indispensable role.We determined that little precipitation was one of reasons for the deterioration of the vegetation.Moreover, barren soil increased the T∆ and added more ULW.In the winter, SC was the main factor of R n variation because it easily changed the albedo and controlled shortwave radiation.The R n over south-central Africa decreased at a rate of ´0.63 W¨m ´2 per year, with NDVI as the main driving factor.Regardless of the analysis period, NDVI was always the leading factor, as evidenced by the high VIP scores and %SS.Different from the case of the Great Plains, precipitation showed little correlation with R n variation.This may be attributed to the little variation of precipitation, particularly during the wet season.
To identify the contributing factors that caused the changes in R n in the two hotspots regions, CRU data were used.Although CRU data is based on interpolated ground stations, the data can be considered as, overall, reliable for our analysis in the hotspot regions.
Ridge National Laboratory, University of California-Berkeley, and the University of Virginia.We would also like to thank Yan Li and Ran Li from the Xinjiang Institute of Ecology and Geography, CAS, and Jean Christophe Calvet from CNRM-GAME for sharing their data.

Figure 1 .
Figure 1.Distribution of the 15 measurement networks.Figure 1. Distribution of the 15 measurement networks.

Figure 1 .
Figure 1.Distribution of the 15 measurement networks.Figure 1. Distribution of the 15 measurement networks.

3. 1 .
Validation of the CERES R n Product 3.1.1.Daily Scale The daily CERES R n samples, from March 2000 to July 2014, were compared with the observations taken from 340 sites.The resulting density of scatter diagram is shown in Figure 3. Remote Sens. 2016, 8, 90 7/19 3.1.1.Daily Scale The daily CERES Rn samples, from March 2000 to July 2014, were compared with the observations taken from 340 sites.The resulting density of scatter diagram is shown in Figure 3.

Figure 3 .
Figure 3. Density of scatter plot of the CERES and in situ daily samples.

Figure 3
Figure 3 indicates that the daily CERES Rn data show a strong consistency with the in situ observations.The MBE of the CERES data samples was 3.43 W•m −2 , the RMSE was 33.56 W•m −2 , and the R 2 was 0.79.The MBE was positive, therefore showing that the CERES Rn data are slightly larger than the in situ observations.The validation results of the CERES Rn, at each site, are shown in Figure4.The grey scale represents the land cover complexity, to be used as a reference.The MBE and RMSE values, at a number of sites in China, are much larger than most sites in Europe and America.This may be related to the measurement accuracy of different networks, as, for example, the La Thuile and the ARM networks have higher accuracy requirements.Compared with other networks, the CEOP, HiWATER, and CEOPInt, built in China, have poorer data continuity and a lower quantity of data.The CERES accuracy in Greenland appears to be lower.Inamdar et al., also found the CERES net shortwave radiation has large errors in the underlying surface snow[54].The CERES Rn showed bad performance on a site in the central southern region of the United States.By checking the 0.05° land cover complexity results, we determined that this site is located on a more complicated underlying surface.Considering that the CERES Rn showed strong agreement with other sites in this region, the accuracy is still considered trustworthy.

Figure 3 .
Figure 3. Density of scatter plot of the CERES and in situ daily samples.

Figure 3
Figure3indicates that the daily CERES R n data show a strong consistency with the in situ observations.The MBE of the CERES data samples was 3.43 W¨m ´2, the RMSE was 33.56 W¨m ´2, and the R 2 was 0.79.The MBE was positive, therefore showing that the CERES R n data are slightly larger than the in situ observations.The validation results of the CERES R n , at each site, are shown in Figure4.The grey scale represents the land cover complexity, to be used as a reference.The MBE and RMSE values, at a number of sites in China, are much larger than most sites in Europe and America.This may be related to the measurement accuracy of different networks, as, for example, the La Thuile and the ARM networks have higher accuracy requirements.Compared with other networks, the CEOP, HiWATER, and CEOPInt, built in China, have poorer data continuity and a lower quantity of data.The CERES accuracy in Greenland appears to be lower.Inamdar et al., also found the CERES net shortwave radiation has large errors in the underlying surface snow[54].The CERES R n showed bad performance on a site in the central southern region of the United States.By checking the 0.05 ˝land cover complexity results, we determined that this site is located on a more complicated underlying surface.Considering that the CERES R n showed strong agreement with other sites in this region, the accuracy is still considered trustworthy.

Figure 4 .
Figure 4. Daily validation results of the CERES Rn data at individual sites (a) MBE; (b) RMSE; and (c) R 2 .Figure 4. Daily validation results of the CERES R n data at individual sites (a) MBE; (b) RMSE; and (c) R 2 .

Figure 4 .
Figure 4. Daily validation results of the CERES Rn data at individual sites (a) MBE; (b) RMSE; and (c) R 2 .Figure 4. Daily validation results of the CERES R n data at individual sites (a) MBE; (b) RMSE; and (c) R 2 .
The monthly CERES Rn samples were compared with the observations taken from 260 sites.The resulting density of scatter diagram is shown in Figure5.

Figure 5 .
Figure 5. Density of scatter plot of the CERES and in situ monthly samples.

Figure 5 .
Figure 5. Density of scatter plot of the CERES and in situ monthly samples.
The monthly CERES Rn samples were compared with the observations taken from 260 sites.The resulting density of scatter diagram is shown in Figure5.

Figure 5 .
Figure 5. Density of scatter plot of the CERES and in situ monthly samples.At a monthly scale, the MBE was 3.40 W•m −2 , the RMSE was 25.57W•m −2 , and the R 2 was 0.84.This indicates that the CERES Rn data are slightly larger than the in situ observations, regardless of noise.The validation results, at each ground site, are displayed in Figure6.Compared with the daily validation results, the performance of the monthly product is much better.It should be noted that there appears to be a poorer precision at a number of sites near the coast, and in Greenland.This reduction in precision may be caused by the large pixel areas with higher land cover complexity.Both daily validation and monthly validation showed lower accuracy of CERES Rn over Greenland.The preprocessing of in situ data and CERES data was very strict, and the results are trustworthy.The validation results demonstrate that the CERES Rn product is, overall, highly accurate, although the accuracy is slightly lower for high latitudes.

Figure 6 .
Figure 6.Monthly validation results of the CERES R n data at individual sites (a) MBE; (b) RMSE; and (c) R 2 .

Figure 7 .
Figure 7. Mean values of the Rn from 2001 to 2013 (a); and annual trend of the Rn (b).The pixels in (b) marked with a dot indicate their statistical significance, with p-values < 0.05.

Figure 7 .
Figure 7. Mean values of the R n from 2001 to 2013 (a); and annual trend of the R n (b).The pixels in (b) marked with a dot indicate their statistical significance, with p-values < 0.05.

Figure 8 .
Figure 8. Annual and seasonal trends of Rn over the southern Great Plains.The pixels marked with a dot indicate their statistical significance, with p-values < 0.05.Units: W•m −2 per year.(a) annual trend; (b) trend in spring; (c) trend in summer; (d) trend in autumn; and (e) trend in winter.The Rn decreased over the entire southern Great Plains by −0.33 W•m −2 per year.The fastest decrease occurred in summer, with a gradient of −0.58 W•m −2 per year, while the slowest decrease in spring was only −0.09 W•m −2 per year.Furthermore, the Rn decreased in the central region, and increased in the west edge.In the summer, the central part decreased sharply at a rate of −2.41 W•m −2 per year.However, also in the summer, the Rn in the northwest increased with a gradient of 1.62 W•m −2 per year.Rn in the east flat area did not show any significant trend.In the significantly changing areas, the correlation analysis was carried out and %SS values were calculated between the CERES Rn and the driving factors, as shown in Table4.

Figure 8 .
Figure 8. Annual and seasonal trends of R n over the southern Great Plains.The pixels marked with a dot indicate their statistical significance, with p-values < 0.05.Units: W¨m ´2 per year.(a) annual trend; (b) trend in spring; (c) trend in summer; (d) trend in autumn; and (e) trend in winter.

Figure 9 .
Figure 9. Variable importance in the projection (VIP) scores of the driving factors in the southern Great Plains.

Figure 9 .
Figure 9. Variable importance in the projection (VIP) scores of the driving factors in the southern Great Plains.

Figure 10 .
Figure 10.Annual and seasonal trends of the Rn variation over south-central Africa.The pixels marked with a dot indicate their statistical significance, with p-values < 0.05.Units: W•m −2 per year.(a) annual trend; (b) trend in dry seasons; and (c) trend in wet seasons.

Figure 10 .
Figure 10.Annual and seasonal trends of the R n variation over south-central Africa.The pixels marked with a dot indicate their statistical significance, with p-values < 0.05.Units: W¨m ´2 per year.(a) annual trend; (b) trend in dry seasons; and (c) trend in wet seasons.

Figure 11 .
Figure 11.Variable importance in the projection (VIP) scores of the driving factors in south-central Africa.

Figure 11 .
Figure 11.Variable importance in the projection (VIP) scores of the driving factors in south-central Africa.

Table 1 .
Metadata for the 15 measurement networks.

Table 1 .
Metadata for the 15 measurement networks.

Table 2 .
Metadata for the Driving Factors.

Table 3 .
Code Names of IGBP Land Cover Types.

Table 5 .
Statistics for the R and %SS between the Rn and the driving factors in south-central Africa.Numbers marked with a * indicates p-values < 0.10 and ** indicates p-values < 0.05, in which case the results are deemed to be significant.

Table 5 .
Statistics for the R and %SS between the R n and the driving factors in south-central Africa.Numbers marked with a * indicates p-values < 0.10 and ** indicates p-values < 0.05, in which case the results are deemed to be significant.