Determining Factors of Monthly to Decadal Variability in Surface Solar Radiation in China: Evidences From Current Reanalyses

Clouds and aerosols play essential roles in regulating surface incident solar radiation (Rs). It has been suggested that the increased aerosol loading over China is a key factor for the decadal variability in Rs and can explain the bias in its trend from reanalyses because the reanalyses do not include the interannual variability of aerosols. In this study, we compare the Rs derived from sunshine duration at 2,400 weather stations in China and that from five reanalyses from 1980 to 2014. The determining factors for the biases in the mean values and trends of Rs from the reanalyses are examined, with the help of Rs and the cloud fraction (CF), from satellite and 2,400 weather stations. Our results show that all reanalyses overestimate the multiyear Rs by 24.10–40.00 W/m2 due to their underestimations of CF, which is more obvious in southern China. The biases in the simulated CF in the reanalyses can explain the biases in Rs by 55–41%, and the bias in clear‐sky surface solar radiation (Rc), which is primarily due to biases in aerosol loading, can explain 32–9% of the bias in Rs. The errors in the trend of the simulated CF can explain the errors in the Rs trends in the reanalyses by 73–12%, and the trend errors in the Rc can explain 43–30% of the trend error in Rs. Our study suggests that more work is needed to improve the simulation of aerosols, clouds, and aerosol‐cloud‐radiation interactions in the reanalyses.


Introduction
Surface incident solar radiation (R s ) is a key component of the surface energy budget. It drives the global climate system and the hydrological and carbon cycles (C. Dorno, 1920;Roderick & Farquhar, 2002;Sedlar, 2018;Sellers et al., 1990;Wang et al., 2017). Widespread solar radiation measurements have shown that R s has significant decadal variability, which is known as global dimming (from the 1950s to 1980s) and brightening (since the mid-1980s, Wild, 2009). The variation in R s is closely related to climate change and global warming (Ruosteenoja & Raisanen, 2013). This has been regarded as an important driver of the observed decadal variability of hemisphere and global mean surface temperature during the last century (Ramanathan et al., 2005;Ruckstuhl et al., 2008;Wang & Dickinson, 2013;Wild, 2009).
Clouds and aerosols have been considered the main contributors to the variability in R s (Bodas-Salcedo et al., 2014;Folini & Wild, 2011;Ghan et al., 2012;Pyrina et al., 2015;Tang et al., 2012;. It is argued that the reduction in aerosol loading is the primary factor of R s brightening from 1980 to 2014 over Europe (Nabat et al., 2015). However, studies based on satellite retrievals have found that the reason for the increasing trend in R s from 1992 to 2015 over Europe is the decrease in the cloud fraction (CF, Pfeifroth et al., 2018). Augustine and Dutton (2013) found that the increasing trend in R s is mainly due to a decrease in CF over the United States from 1996 through 2011.
Changes in aerosol loading have been reported to be the primary cause of variations in R s over China (Che et al., 2005;J. Li et al., 2018;Liang & Xia, 2005;Qian et al., 2015;Xia, 2010). For example, Li et al. (2018) find that aerosol loading was the most likely major cause of variations in R s from 2005 to 2015 over China, especially for the northeast and southern parts of China. In contrast, Tang et al. (2017) conclude that clouds and the interactions between clouds and aerosols are the main reasons for the variation in R s .
The above-mentioned studies on the attributions of the trend in R s to clouds and aerosols are primarily qualitative due to limited numbers of surface measurements of R s , CF, and aerosol loading. Observations of R s are sparsely distributed and have limited spatial and temporal coverage (Wild, 2016). Moreover, the solar radiation measurements over China suffer from problems such as sensitivity drift and instrument replacement (Wang, 2014;Wang et al., 2015;Yang et al., 2018;Zhang & Lu, 1988. Sunshine duration (SunDu) records have been suggested that can be used to reconstruct long-term R s with comparatively large spatial coverage (He et al., 2018;Matuszko, 2014;Sanchez-Lorenzo et al., 2008;Wang et al., 2015b;. Satellite-derived R s products also provide an estimate of R s (Feng & Wang, 2018b;Ma et al., 2015). More importantly, satellite data provide CF in addition to the ground-based manual observations of CF.
Compared with observations and satellite data, the R s reanalysis data have complete spatial and temporal coverage by combining observations and forecast models. The radiative fluxes from the reanalyses are widely used as forcing data in many climate models (Essou et al., 2017). However, R s reanalyses contain substantial biases (Slater, 2016;Wang et al., 2015;Wu et al., 2015) due to the uncertainties in clouds and aerosols in the reanalyses (Fujiwara et al., 2017). Complete knowledge of the interactions among clouds, aerosols and radiation, and the related parameterizations in the climate models or reanalyses can help to reduce the uncertainty in predicting potential future climate changes, especially at regional scales (Loew et al., 2016;Zadra et al., 2018).
In this study, we compare the R s , CF from five current reanalyses with field observations, satellite retrievals, and the satellite and reanalyses-derived clear-sky surface solar radiation (R c ), which is closely related to aerosols, by analyzing the climatologies, spatial patterns, seasonal variations, and trends in R s , CF, and R c from these reanalyses. Second, we evaluate the relationship between the bias in R s and those in the CF and R c . These analyses help to figure out the determining factors of monthly to decadal variability in surface solar radiation in the reanalysis system and observations in China from 1980 to 2014.
where n represents the measured sunshine duration and K represents the theoretical value of the sunshine duration. a0, a1, and a2 are determined following the method of Wang (2014). R c is the daily total solar radiation under clear sky. I is solar irradiance at the top of the atmosphere. T b and T d denote atmospheric transmittance for direct solar radiation and diffuse solar radiation. h is the solar elevation.
Existing studies (Manara et al., 2015;Sanchezlorenzo et al., 2013;Tang et al., 2011;Yang et al., 2018) have shown that SunDu-derived R s is a reliable R s proxy data set at time scales ranging from monthly to decadal that can reflect the impacts of aerosols and clouds on R s over China. SunDu data are relatively widely distributed and have a long time record, which extends from the late nineteenth century to the present (Sanchez-Lorenzo et al., 2009;Wild, 2009).
Moreover, R s measurements such as pyranometers require careful calibrations (Wood et al., 2015;Yang et al., 2018) due to the thermal offsets (Philipona, 2002;Zo et al., 2017) and the directional response errors (Myers et al., 2002). Moreover, R s measurements are sparsely distributed in China. In spite of uncertainty in short time scales, SunDu-derived R s data have their advantage in quantifying annual to decadal variability of R s (Feng & Wang, 2018a;Wang, 2014).
We also collected synoptic observations of the total CF from 2,400 weather stations, which are collocated with the SunDu measurement sites. Daily observational total cloud amounts (in tenths), which are observed by human eyes 4 times a day (2:00, 8:00, 14:00, and 20:00) based on the specifications for surface meteorological observation, are collected and then averaged to get monthly mean observed cloud amount and are linearly converted into the percentage scale to enable comparison with the reanalysis data.

Satellite Data
R s data from the latest release of the Clouds and the Earth's Radiant Energy System energy balanced and filled product (CERES EBAF) surface product (edition 4) and cloud amounts from the CERES Synoptic (CERES SYN1deg) edition 4 product are collected in this study (Kato et al., 2018). CERES is a radiometer on board the Tropical Rainfall Measuring Mission, Terra, Aqua, and Suomi National Polar-orbiting Partnership satellites and National Oceanic and Atmospheric Administration-20. CERES measures three filtered radiances, including shortwave (0.3-5 μm), total (0.3-200 μm), and window (8-12 μm).
The CERES-EBAF product provides global monthly mean R s data. The input cloud properties for CERES EBAF, such as optical thickness and emissivity from Moderate Resolution Imaging Spectradiometer (MODIS) and geostationary satellites, are constrained by a cloud profiling radar, Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations detectors, and CloudSat. MODIS collection 5-derived aerosol data are used as input data for the CERES EBAF product. R s of CERES EBAF is produced with a radiative transfer model constrained by CERES observations at the top of the atmosphere (Kato et al., 2018).
Compared with the previous version of CERES EBAF (Ed2.8), the newly released product (version 4.0) has extensive improvements in the inputs and algorithms. For example, the newly released version 4.0 product eliminated the discontinuities in the temperature and humidity time series at the beginning of 2008. Furthermore, the biases in temperature and specific humidity from the Goddard Earth Observing System (GEOS) version 5.4.1 reanalysis are corrected using atmospheric infrared sounder data. All the improvements help to improve its capability in quantifying long-term trend of R s . The uncertainties of CERES EBAF data, reported by Kato et al. (2018), in all-sky global annual mean R s is 4 W/m 2 . The new CERES products are supplemented with geostationary satellite data between the twice-per-day passes to document diurnal variation of solar and IR irradiance overcoming a major flaw in polar orbiter (i.e., CERES) data.
Existing studies (Feng & Wang, 2018b;Ma et al., 2015;Wang et al., 2015) also show that the CERES EBAF surface product provides reliable estimations of monthly R s . Feng and Wang (2018a) find that CERES EBAF has a bias of 7.94 W/m 2 compared with observation and 7.53 W/m 2 compared with SunDu-derived R s from 2000 to 2014. The comparison results of CERES EBAF, GEWES-SRB, and ISCCP-FD-SRF from Li, Xin, and Peng (2017) show that CERES EBAF have lowest root-mean-square error (RMSE, 14.73 W/m 2 ) compare with solar radiation measurements in China.

Reanalysis Data
Five of the latest global reanalysis products, including the European Centre for Medium Range Weather Forecasting Reanalysis-Interim (ERAI), Japanese 55-year reanalysis (JRA55), Climate Forecast System Reanalysis (CFSR), Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA2) and MERRA, were collected in this study. The brief descriptions of these five reanalyses are summarized in Table 1. It is necessary to note that the time span of the CFSR is from 1979 to 2010. We use its successor, version 2 (CFSR2), to provide the subsequent CFSR estimates from 2010 to 2014 according to previous studies, in which CFSR2 has been used as a continuation of CFSR due to few changes in the physical models (Suranjana Saha et al., 2014). ERAI does not provide a R c product, and we calculate the R c from its albedo and surface net solar radiation products for a clear sky.
These five reanalyses do not assimilate any R s data from conventional or satellite observations. R s derived from these reanalyses is calculated using radiation transfer models. Specifically, the ERAI used a six-band shortwave parameterization model developed by Fouquart and Bonnel (1980). JRA55 used a shortwave parameterization model developed by Briegleb (1992). The shortwave parameterization of the CFSR is based on the Rapid Radiative Transfer Model for General Circulation Models developed by Clough et al. (2005). Both MERRA2 and MERRA applied the Goddard climate and radiation package (CLIRAD) for shortwave parameterization (Fujiwara et al., 2017).
The simulation of the CF in reanalyses falls into two categories: (1) the diagnostic scheme and (2) the prognostic scheme. CF in the diagnostic scheme is parameterized by an empirical function of relative humidity or a probability density function (PDF) of moisture variability. In the prognostic scheme, the CF is calculated using a prognostic equation of the sources and sinks of cloud areas. For the reanalyses, the CF from the ERAI is calculated using a prognostic scheme (Berrisford et al., 2009). CFSR uses a diagnosed scheme for CF computation (Saha et al., 2006). JRA55 uses a PDF-based diagnosed cloud scheme (Smith, 2010). MERRA applies a prognostic cloud scheme (Bacmeister et al., 2006), and MERRA2 updates its atmospheric model, which results in increasing cloud condensation that has a substantial impact on CF simulation (Molod et al., 2015).
ERAI, CFSR, MERRA2, and MERRA use a maximum-random overlap scheme to depict different vertical layers with cloud overlap, whereas the JRA55 uses a random overlap scheme. Climatological aerosol values are used in the ERAI, CFSR, JRA55, and MERRA. MERRA2 uses analyzed aerosol values (Randles et al., 2017), which are produced by the Global Ozone Chemistry Aerosol Radiation and Transport model, and assimilates satellite-derived Aerosol Optical Depth (AOD) values and Aerosol Robotic Network AOD values.

Methods
We assess the R s values and CF, from the reanalyses by comparing them with surface observations and satellite retrievals. Both satellite retrievals and field observations are used as all-sky R s reference. Clear-sky surface solar radiation (R c ), which is closely related to atmospheric aerosol loading, can be calculated with input from reanalyses and satellite. R c from reanalyses were also compared with satellite retrievals to illustrate the impact of aerosols.   To eliminate the impact of surface observations distribution inhomogeneities, we interpolate the field observations into 1°× 1°grid cells using the area-weighted averaging method (Du et al., 2018). To be specific, we divide the study region into 1°× 1°grids covering China and assigned all sites to the grids. For a certain grid contained more than one site, the area-weighted average of these sites are calculated as this grid value. For consistency, all data are transformed into 1°× 1°grids based on the bilinear interpolation method.
The spatial distributions of the multiyear mean R s , CF, and R c over China, along with the biases in these parameters from the reanalyses, are analyzed in this study. The linear trends in R s , CF, satellite-derived R c , and R c from reanalyses are calculated by the least squares method. The trends are evaluated using observations and satellite retrievals. To ensure consistent comparisons, all the trends are also calculated by using 1°× 1°grids data including sites observation. To quantify the impacts of CF on all-sky R s error between the reanalysis data and observations, correlation coefficients (R) of linear regression plots between CF error and all-sky R s error are calculated. It is essential to assess the sensitivity of error in simulated R s to CF and R c . Therefore, we calculated a polynomial regression model to analyze the sensitivity of R s to the biases in the CF and R c where β i represents the sensitivity, which is the slope of the linear regression line between the bias in R s plotted against the bias in x i . Δx i represents the bias in the annual anomaly for the related variable (e.g., CF and R c ). ΔRs represents the bias in the annual anomaly R s , and m and ε represent the constant and residual, respectively. The coefficient of determination (R 2 ) of those linear plots is used to determine the degree to which Δx i explain the variability of ΔRs. The R s trend errors are calculated as follows: where δ represents the sensitivity in the trend error and ΔT xi represents the trend error in the annual anomaly for the related variable (e.g., CF and R c ). ΔT Rs represents the trend error in the annual anomaly R s , and β i is calculated using equation (2). The national mean sensitivity in the trend error is calculated by the mean of absolute values.

Results
3.1. Multiyear Mean R s , CF, and R c Over China   larger than that over other regions of China (180-260 W/m 2 , Figure 2), and the variation in R c in the Tibetan Plateau is small due to its high altitude, low aerosol loading and water vapor concentration, not to mention the low value center for total column ozone over the Tibetan Plateau (Zhou et al., 2006).   Figure 2.
Reanalyses R s display larger relative biases in winter and early spring than those in other seasons ( Figure 3). The zonal mean relative biases in R s , CF, and R c are shown in Figure 4. The spatial patterns of the relative bias in R s further indicate that large biases in winter and spring mainly occur in the area between 25 and 30°N. The relative bias in the CF shows an opposite spatial pattern compared with that in R s , especially via MERRA and MERRA2. The R c estimated by the reanalyses exhibits an overall low relative bias for different seasons.

Trends of R s , CF, and R c Over China
The trends in R s , CF, and R c over China from 2000 to 2014 and the corresponding biases via the satellite and reanalyses are shown in Figure 5. Generally, observed R s shows a national mean trend of −0.89 W/m 2 per decade ( Table 3). The observed CF shows an overall increasing trend (0.02 per decade). The national mean decreasing trend in R s and the national mean increasing trend in CF can also be seen during the period of 1980-2014. In spatial terms, most of China shows a slight decreasing trend in R s ( Figure 5). The increasing trend in R s is comparatively high in southwest China. Conversely, the increasing trend in the CF can be seen in many parts of China, except in southwest China. Based on the CERES EBAF data, the R c shows a decreasing trend in eastern China and an increasing trend in north China. Moreover, the distribution diagrams of Figure 5 indicate that reanalyses have more positive R s trend bias and negative CF trend bias. Based on the CERES EBAF data, the national mean trend in R c from 2000 to 2014 is −0.88 W/m 2 per decade.
Compared with the reanalyses, the CERES EBAF shows a more consistent national mean R s trend and CF trend with those of the ground observations (Table 3). All reanalyses show poor performance in simulating Table 3 Statistical Summary of the Linear Trends (10 years

10.1029/2018JD030214
Journal of Geophysical Research: Atmospheres CF trends and R c trends compared with the CERES EBAF data (Table 3 and Figure 5). ERAI and MERRA produce large, positive R s trend biases ranging from 10 to 20 W/m 2 per decade and correspondingly large negative CF trend biases (0 to −0.3 per decade) in east China. For R c simulation results, CERES EBAF produced a decreasing trend in R c in the North China Plain (32-40°N, 114-121°E) and east China. This might be attributed to the increasing trend in aerosol loading in these areas that is consistent with previous studies (Tang et al., 2017;Yang et al., 2009). ERAI, JRA55, CFSR, and MERRA show overestimated R c trends in the North China Plain because these reanalyses do not include the increase in aerosols that corresponds a decreasing trend in R c . ERAI exhibits negative biases in the R c trend in northern China, and therefore, the uncertainties in the water vapor trend from ERAI might not be ruled out (Ning et al., 2013). Figure 6 illustrates the correlation between the annual anomalies of CF, R c , and annual anomalies of R s . Based on observations, reanalyses fail to produce the weak negative correlation between the annual Figure 7. Maps of the correlation coefficient (R) between the bias in the annual anomaly cloud fraction (CF) and the bias in the annual anomaly surface solar radiation (R s ). In the left column (a, d, g, j, and m), the observed CF and the observed R s represent the CF reference data derived from the observations, and the R s reference data are the SunDu-derived R s data. In the middle column (b, e, h, k, and n), CERES CF and CERES R s represent the CF reference data and R s reference data are calculated by the CERES EBAF data. In the right column (c, f, i, l, and o), CERES CF and observed R s represent the CF reference data and the CERES EBAF data, and the R s reference data are the SunDu-derived R s data. The time span is from 2000 to 2014. The probability density plot embedded in each subplot is the same with Figure 2 but stands for the distribution of correlation coefficient. CERES EBAF = Clouds and the Earth's Radiant Energy System energy balanced and filled product.

10.1029/2018JD030214
Journal of Geophysical Research: Atmospheres anomalies of CF and R s in northern China because most reanalyses ignore the interannual variability in dust, which leaves the interannual variability in R s from reanalyses highly dependent on the CF (Fujiwara et al., 2017;Randles et al., 2017;Zhang et al., 2019). The distribution diagrams of Figure 6 further show that reanalyses have more negative correlation between the CF and R s than that of observations.
Positive correlations between annual anomalies in R c and annual anomalies in R s are found in northern China based on the observations and CERES EBAF data ( Figure 6). However, only MERRA2 shows consistent correlation spatial patterns, which is likely because of its improvement in aerosol simulation. Figure 7 shows the correlation between the bias in annual anomaly CF and the bias in the annual anomaly R s . Using field observations as reference, reanalyses produce negative correlation between the bias in the annual anomaly CF and that in the annual anomaly R s , with the national mean determination coefficient (R 2 ) ranging from 0.36 to 0.43 (Table 4). Figure 7 also show that reanalyses have stronger negative correlation when both CF and R s reference data use CERES EBAF data. Because R s and CF are independently observed for the field observations, the CERES EBAF R s is calculated from its CF. In spatial terms, the negative correlation coefficient between the annual anomaly CF bias and annual anomaly R s bias is high in southern China (−0.50 to −1.00) and low in northern China (−0.40 to 0.00), compared with field observation. Similar spatial patterns can also be found by using the CERES EBAF data as a reference, excluding MERRA2, which shows a high negative correlation in northeast China.
The sensitivity of the annual anomaly bias in R s to that in the CF from reanalyses (Figure 8) correspondingly exhibits an almost similar spatial pattern compared with the correlation results from Figure 7. MERRA2 has less sensitivity of the annual bias in R s to that in the CF than MERRA, ERAI, JRA55, and CFSR because these reanalyes except MERRA2 do not include the impact of interannual variability in aerosols on R s . Figure 9 shows the correlation between the bias in annual anomaly R c and that in the annual anomaly R s . All reanalyses show an almost positive correlation between the bias in annual anomaly R c and that in the annual anomaly R s with national means of R 2 ranging from 0.08 to 0.32 (Table 5). In spatial terms, a high positive correlation coefficient between the bias in the annual anomaly R c and the bias in the annual anomaly R s over China can be seen in northern China based on JRA55, CFSR, MERRA2, and MERRA, which is probably because the low CF in northern China has a small impact on the variation in R s , while aerosols and water vapor have large contributions to the modification of the variation in the annual anomaly R s in northern China. Comparatively, the ERAI shows a low correlation coefficient between the bias in the annual anomaly R c and the bias in the annual anomaly R s , which might be due to the calculation method for R c in the ERAI. We calculate the R c using the ERAI clear-sky net R s and its forecast albedo product. However, the calculated R c may be a little different from the R c data used in the original ERAI radiation transfer calculation (Hogan, 2014).
For the impacts of R s mean bias, the correlation results (Table 4) show that the biases in the CF in MERRA can largely explain the biases in R s by 55%, followed by ERAI (45%), CFSR (43%), JRA55 (42%), and MERRA2 (41%) using satellite data as a reference. The biases in R c can explain the biases in R s by 32% in JRA55, 29% in MERRA, 22% in CFSR, 14% in MERRA2, and 9% in ERAI using satellite data as a reference (Table 5). The correlation results suggest that the bias in the CF has more impacts on the bias in R s than that in R c . Figure 10 shows the sensitivity of trend error in R s to that in CF, and Figure 11 illustrate the sensitivity of trend error in R s to that in R c . In spatial terms, the trend error in the CF in all five reanalyses (except Figure 8. Maps of sensitivity between the bias in the annual anomaly cloud fraction (CF) and the bias in the annual anomaly of surface solar radiation (R s ) for the five reanalyses using various sources as the reference. In the left column (a, d, g, j, and m), the observed CF and the SunDu R s represent the reference data. In the middle column (b, e, h, k, and n), the CERES CF and CERES R s represent the reference data. In the right column (c, f, i, l, and o), the CERES CF and the SunDu R s represent the reference data. The time span is from 2000 to 2014. The probability density plot embedded in each subplot is the same with Figure 2 but stands for the distribution of sensitivity between the bias in the annual anomaly CF and the bias in the annual anomaly of R s . CERES EBAF = Clouds and the Earth's Radiant Energy System energy balanced and filled product. MERRA2) can largely explain the trend error in R s, especially in east China ( Figure 10). However, the trend error in R c in all five reanalyses, except MERRA2 and JRA55, can only slightly explain the trend error in R s (Figure 11). The distribution diagrams of Figures 10 and 11 further show that the trend error in R s is more sensitive to that of CF than the trend error in R c .
In summary, the trend error in the CF from 2000 to 2014 in MERRA can explain the trend error in R s by 73%, followed by ERAI (64%), CFSR (44%), MERRA2 (35%), and JRA55 (12%) using the CERES EBAF data as a reference ( Table 6). The trend error in R c from 2000 to 2014 can explain the trend error in R s by 36% in ERAI, 32% in JRA55, 34% in CFSR, 43% in MERRA2, and 30% MERRA using the CERES EBAF data as a reference.

Performances of the Reanalyses
Using the best records of historic surface radiation data and trends by Sun duration observed data, this study presents a comprehensive analysis of the performance of five reanalyses commonly used in atmospheric research with regard to geographically distributed surface solar radiation trends over China. Our results show that all five reanalyses are wrong in the same direction with regard to R s , regardless of the different methods and parameterizations used for CF and aerosols, which are consistent with previous findings (Jia et al., 2013;Wang et al., 2015;Wu et al., 2015;Xia et al., 2006). This is largely because the current climate models have difficulty simulating low-level clouds, such as stratus clouds in southern China (Yu et al., 2001;Yu et al., 2004). High AOD in winter and aerosol-cloud interaction make it worse (Li et al., 2016;. Results also demonstrate the dangers of using only radiative transfer and climatological aerosol information without data assimilation in the reanalysis packages Our sensitivity analyses illustrate that the biases in CF in the reanalyses can explain the biases in R s by 55-41%, and the biases in R c can explain the biases in R s by 32-9%. These two factors together can explain the biases in R s by 87-50%. The biases in R s are more sensitive to the CF than those in the R c in the reanalyses. The improvement in aerosol simulation of MERRA2 recently draws much attention (Randles et al., 2017).
Our previous study Wang, 2014) also shows MERRA2 has a better performance in the simulation of the R s trend and R c than MERRA, and the correlation between CF and R s in northern China from MERRA2 is smaller than MERRA. However, the R s trend error in MERRA2 is still very large because the simulation error of clouds still exists.
CERES EBAF R s retrievals have been used to compare with the reanalyses. Our previous study also suggests that CERES EBAF R s retrievals have 9.9 W/m 2 mean absolute biases comparing with BSRN sites data (Feng  & Wang, 2018b). Existing study suggest that changing the MODIS AOD from Collections 4 to 5 resulted in the discontinuity in the CERES EBAF global land R s in 2006 (Jia et al., 2018). Moreover, satellite AOD retrievals are only available under cloud-free conditions, which cannot fully represent AOD trend in the all-weather situation. Although upgrades have been made in CERES EBAF data recently, the uncertainties of MODIS AOD trend used in CERES EBAF data might not be ignored. Likewise, the limited number of observations in the west part of China might have potential impacts on the application of the area-weighted method in data-void regions and its uncertainties cannot be ruled out. Moreover, SunDu data are derived from human read of a burned signal. Alternation of human observer should produce inhomogeneous time series of SunDu data. Although the SunDu-derived R s used in this study is processed carefully with quality control (He et al., 2018;Wang et al., 2015), the uncertainties of SunDu data also cannot be ignored.

Impact Factors of Variability in R s
Previous studies (Boers et al., 2017;Pfeifroth et al., 2018;Wang, Dickinson, et al., 2012) have confirmed that clouds and aerosols are the dominant factors in the variation in R s . Comparatively, radiatively active gases Figure 10. The sensitivity of trend error in annual anomaly surface solar radiation (R s ) to that in the annual anomaly cloud fraction (CF). The left column is calculated based on field observations, the middle column is calculated based on CERES EBAF data, and the right column shows the cross terms. The time span is from 2000 to 2014. The probability density plot embedded in each subplot is the same with Figure 2 but stands for the distribution of the sensitivity of trend error in annual anomaly R s to that in the annual anomaly CF. CERES EBAF = Clouds and the Earth's Radiant Energy System energy balanced and filled product. and water vapor have minor effects on the variation in R s (Kvalevåg & Myhre, 2007;Mateos et al., 2013;Wild, 2016) and indicated that aerosols are the main factor affecting the variation in R c . Large aerosol emissions due to rapid development in China are expected to be the main factor for the variation in R s over China (Liang & Xia, 2005;Qian et al., 2015;Wang, Yang, et al., 2012). Li et al. (2018) further stated that the variation in single-scattering albedo can also affect the trend in R s , especially for diffuse radiation. However, the study conducted by Tang et al. (2017) suggests that clouds are likely to be the main contributor to the variation Figure 11. The sensitivity of trend error in annual anomaly surface solar radiation (R s ) to that in the annual anomaly surface solar radiation under clear-sky condition (R c ). In the left column (a-e), the R s trend bias is calculated based on field observations, and the right column (f-j) is calculated based on CERES EBAF data. The time span is from 2000 to 2014. CERES EBAF = Clouds and the Earth's Radiant Energy System energy balanced and filled product.
in R s . The controlling factors of variability in R s is difficult to qualify, systematic studies including all impact factors may conduct only at site scale or using climate models. All these factors jointly lead to the different conclusions from earlier studies. Therefore, improving the R s simulation model such as reanalyses is helpful for further studies on this issue.
Atmospheric aerosol loading is fixed in most reanalyses except MERRA2; the major driver of R c produced by most reanalyses are water vapor. Compared with CERES EBAF, the main error sources of R c data from reanalyses are mainly aerosol, although water vapor cannot be rule out (Dolinar et al., 2016). This is also true for MERRA2.
Cloud overlap schemes have close relationships with the vertical resolutions in the reanalyses (Oreopoulos & Khairoutdinov, 2003;Stephens et al., 2004), which can also introduce uncertainties in the simulation of CF. Wang et al. (2016) found that there were nonnegligible biases in simulating clouds and their shortwave radiation effects by using a maximum and random overlap scheme. Recently, the newly developed Monte Carlo independent column approximation scheme is expected to have an overwhelming advantage in simulating clouds compared with other cloud overlap schemes.
This study discusses the impacts of CF on R s in the reanalyses and the contribution of aerosols to R s variability. Other cloud and aerosol properties might also be important factors, for example, changes in cloud radiative properties, such as cloud optical thickness (Deneke et al., 2008;McCoy et al., 2014). Cloud height, which is relative to the distributions of low clouds, also impacts the variation in R s Matuszko, 2012). Li et al. (2018) show that the variation in single-scattering albedo has a substantial impact on the trend in R s .
Moreover, elevation also impacts R s . For example, comparison results of observed R s , GEWEX-SRB, and ISCCP-FD in Tibet from  show the differences in elevations and these products might result in nonnegligible uncertainties. Du et al. (2018) also find that the elevation bias in reanalyses has important impacts on temperature biases of reanalyses.

Conclusions
Our results show that all reanalyses overestimate the multiyear mean R s over China (24.10-40.00 W/m 2 ) and correspondingly underestimate the multiyear mean CF (−0.18 to −0.03) for the highest and lowest of reanalyses, especially in southern China. ERAI and CFSR have a larger bias in simulating R c compared with those of JRA55, MERRA2, and MERRA. ERAI and CFSR show high positive biases in R c in south China and negative biases in northwest China. All reanalyses have larger R s relative biases in winter and spring than other seasons.
We partite the error source between clouds and aerosols and find the biases in the CF can explain the biases in R s by 41-55%, and the bias in R c , which is primarily due to errors in atmospheric aerosol loading, can explain 9-32% of the bias in R s in the reanalyses. The trend error in the CF in the reanalyses can largely explain approximately 12-73% of the trend error in R s , while the trend error in R c can explain 30-43% of the trend error in R s using CERES EBAF as a reference. In east China, the trend error in R s can be largely explained by those in the CF in the five reanalyses (excluding MERRA2). The trend error in R c in all five reanalyses have a weak relationship with those in R s (excluding MERRA2 and JRA55). The biases in R s are more sensitive to the biases in the CF than those in the R c . Reanalysis data are commonly used as truth and credible support for atmospheric research; our study suggests that improvement of the cloud and aerosol representation in reanalyses is needed, especially for aerosol-cloud interactions parameterization.