Merging Satellite Retrievals and Reanalyses to Produce Global Long-Term and Consistent Surface Incident Solar Radiation Datasets

Surface incident solar radiation (Rs) is a key parameter in many climatic and ecological processes. The data from satellites and reanalysis have been widely used. However, for reanalysis, Rs data has been shown to have substantial spatial bias, and the time span of reliable satellite Rs is too short for climatic and ecological studies. Combining reanalysis and satellite data would be an effective method for generating long-term and consistent Rs datasets. Here, we apply a cumulative probability density function-based (CPDF) method to merge eight reanalyses with the latest available satellite Rs data from Clouds and Earth’s Radiant Energy System Energy Balanced and Filled (CERES EBAF) surface retrievals. The CPDF method not only reduces the spatial bias of the reanalysis Rs data, but also makes the Rs datasets in a global, long-term and consistent way. The observed Rs data collected at 54 Baseline Surface Radiation Network (BSRN) stations from 1992 to 2016 are used to evaluate the method. Results show that the CPDF method could reduce the mean absolute biases (MAB) of the reanalysis Rs effectively by 21.24–64.36%. The European Centre for Medium-Range Weather Forecasts Re-Analysis interim (ERA-interim) reanalysis Rs data, which are available for 1979 onward, perform the best before MAB = 13.20 W·m−2 and after MAB = 10.40 W·m−2 merging. This small post-merging MAB of the ERA-interim reanalysis is caused by the MAB of 9.90 W·m−2 in the satellite Rs retrievals. The Japanese 55-year reanalysis provides Rs values back to 1958, and CPDF can reduce its MAB by 32.87%, to 11.17 W·m−2. The National Oceanic and Atmospheric Administration (NOAA)-CIRES twentieth-century reanalysis (CIRES) and the ECMWF twentieth-century reanalysis (ERA20CM) provide century-long Rs estimates. CIRES performs better after merging. The MAB of CIRES can be reduced by 32.10%, to 12.99 W·m−2, while ERA20CM’s can be reduced by 12.51%, to 16.40 W·m−2.


Introduction
Surface Incident solar radiation (R s ), which is also often referred to as the downward solar irradiance, is a key parameter in many climate and ecological processes, such as evapotranspiration [1][2][3], canopy photosynthesis [4], net primary production [5,6], crop growth management [7], and so on.Long-term R s datasets with global coverage and reasonable accuracy have a great value these days.Globally-distributed ground observations of R s began in 1958, and provide solid evidence for global dimming and brightening [8][9][10].However, there are recent studies showing that such observations have poor performance in quantifying the decadal variability in R s [10][11][12], indicating that accuracy estimation of long-term trends should depend on the measurement methods and instrument maintenance [13].Therefore, several long-term observation networks were developed under the high-quality control.The Baseline Surface Radiation Network (BSRN) [14] which provides high-quality R s observations, has been widely used in the assessment of R s variation at local sites [15] and the evaluation of R s derived by satellite data and reanalyses [16,17].The Surface Radiation Budget Network (SURFRAD), the U.S. component of BSRN, was developed by the National Oceanic and Atmospheric Administration (NOAA) to provide accurate, continuous, and long-term R s data for the United States [18].There are also studies [19][20][21][22][23][24][25] suggesting that sunshine duration records can be used to reconstruct long-term R s with reliable accuracy.However, both ground-based observations and sunshine duration records are spatially sparse.
Satellite R s retrievals provide global coverage.Many satellite-based R s products have been developed, such as the Global Energy and Water Cycle Experiment Surface Radiation Budget (GEWEX-SRB) [26].R s values calculated based on satellite cloud observations, i.e., the International Satellite Cloud Climatology Project (ISCCP) [27,28], have been shown to be inhomogeneous in depicting long-term variability [29,30].Inhomogeneities in cloud products, caused by the different capabilities of geostationary and polar-orbiting satellites for the detection of low-level clouds from ISCCP, have introduced significant inhomogeneity in the R s data from the GEWEX-SRB.Another satellite R s product, Clouds and Earth's Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) product, is also widely used.Different from the GEWEX-SRB product, inhomogeneity issues have been corrected in the CERES EBAF data, by merging the cloud amount data from polar-orbiting satellites and the mean diurnal cycle of cloud cover from the geostationary satellite data [31].Studies have shown that the CERES EBAF surface product has a better performance than other R s products [10,32,33] due to the top of the atmosphere (TOA) observation constraint.Nevertheless, the temporal coverage of the recently developed high-accuracy satellite R s product is relatively short, and limits the application of these products in climate change studies.
Global R s data is provided in the reanalysis systems as well.Current reanalysis systems calculate cloud properties based on the assimilation of water vapor and temperature profiles from radiosonde and satellite observations.However, spatial biases in these data have been found, by comparing R s data derived from reanalysis datasets with field observations and satellite data in numerous validation studies [10,[34][35][36].Wang et al. [10] find that the R s data derived from the Modern Era Retrospective-analysis for Research and Applications (MERRA) reanalysis data and ERA-interim (ERAI) products have a large spatial bias of 29.8 W•m −2 and 19.1 W•m −2 in China, respectively.Boilley and Wald [35] conclude that overestimations of R s by MERRA are probably due to the underestimation of cloudy conditions.Jia et al. [37] find that R s from the National Centers for Environmental Prediction and National Center for Atmospheric Research Reanalysis (NCEP1) has a poor correlation with field observations, showing a root mean square error of 56.5 W•m −2 over East Asia.
The latest satellite R s retrievals from the CERES EBAF surface product can be used to improve the reanalysis R s estimates.However, the CERES EABF-surface product is only available after the year 2000.Direct bias correction can be performed for the reanalysis R s data for this period.To produce long-term and consistent values of R s , a new method is needed to correct the reanalysis R s data before 2000.
Here, we apply a cumulative probability density function-based method (CPDF) to merge R s satellite retrievals and reanalysis data.One can use the probability density distribution curves of R s from the CERES EBAF surface to adjust the R s values from reanalyses after 2000.Accordingly, R s data from reanalyses before 2000 can also be adjusted, using the derived distribution curves.We used the CPDF method to merge the monthly distribution of satellite R s data and the reanalysis R s data.After that, merged R s data before and after the year 2000 were evaluated by R s observations from 1992 to 2016, collected at globally-distributed BSRN stations [14].

Ground Observations of Surface Solar Radiation
The BSRN was established by the World Climate Research Programme (WCRP), to provide high-quality surface radiation measurements with high temporal resolution [14].R s data from the BSRN are available for 1992 onward, and cover the major global climate zones.Each BSRN site has an independent principle investigator, and the data from each BSRN site undergo strict quality control.The R s measurement instruments at BSRN sites are calibrated every six months, regularly maintaining an R s uncertainty of less than 5% [12,38,39].In this study, we used R s data from 54 BSRN stations (Figure 1).We strictly followed the criterion described in [40], that all 24-hourly values must be present in order to compute a monthly mean.

Ground Observations of Surface Solar Radiation
The BSRN was established by the World Climate Research Programme (WCRP), to provide highquality surface radiation measurements with high temporal resolution [14].Rs data from the BSRN are available for 1992 onward, and cover the major global climate zones.Each BSRN site has an independent principle investigator, and the data from each BSRN site undergo strict quality control.The Rs measurement instruments at BSRN sites are calibrated every six months, regularly maintaining an Rs uncertainty of less than 5% [12,38,39].In this study, we used Rs data from 54 BSRN stations (Figure 1).We strictly followed the criterion described in [40], that all 24-hourly values must be present in order to compute a monthly mean.

Satellite Surface Solar Radiation Retrievals
Here we used monthly mean satellite Rs retrievals from the CERES EBAF surface data.CERES is a three-channel (0.3-5 µm, 8-12 µm, and 0.3-200 µm) radiometer that is equipped on the Tropical Rainfall Measuring Mission satellite (TRMM), Terra, Aqua, and Suomi National Polar-orbiting Partnership (S-NPP) satellites.The goal of the CERES instrument is to measure solar-reflected and Earth-emitted radiation from the top atmosphere to the surface Earth.TOA fluxes from CERES are energy-balanced, and serve as a constraint to the surface fluxes.To achieve this, the input datasets were adjusted within their range of uncertainty.To be specific, the cloud property data were then adjusted, using cloud profiling radar and Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) detectors [41], and the temperature and humidity profiles are adjusted using the Atmospheric Infrared Sounder.Based on a modified Lagrange multiplier method, these adjustments of input parameters were then used to calculate the Rs adjustment and add these values to the original estimations of Rs to produce a final monthly mean Rs dataset [42].The CERES EBAF surface product provided accurate Rs data, using adjusted input data.The cloud property retrievals primarily came from observations of the Moderate Resolution Imaging Spectroradiometer (MODIS)

Satellite Surface Solar Radiation Retrievals
Here we used monthly mean satellite R s retrievals from the CERES EBAF surface data.CERES is a three-channel (0.3-5 µm, 8-12 µm, and 0.3-200 µm) radiometer that is equipped on the Tropical Rainfall Measuring Mission satellite (TRMM), Terra, Aqua, and Suomi National Polar-orbiting Partnership (S-NPP) satellites.The goal of the CERES instrument is to measure solar-reflected and Earth-emitted radiation from the top atmosphere to the surface Earth.TOA fluxes from CERES are energy-balanced, and serve as a constraint to the surface fluxes.To achieve this, the input datasets were adjusted within their range of uncertainty.To be specific, the cloud property data were then adjusted, using cloud profiling radar and Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) detectors [41], and the temperature and humidity profiles are adjusted using the Atmospheric Infrared Sounder.Based on a modified Lagrange multiplier method, these adjustments of input parameters were then used to calculate the R s adjustment and add these values to the original estimations of R s to produce a final monthly mean R s dataset [42].The CERES EBAF Remote Sens. 2018, 10, 115 4 of 20 surface product provided accurate R s data, using adjusted input data.The cloud property retrievals primarily came from observations of the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument.Atmospheric parameters, such as temperature and humidity profiles, were derived from the MERRA reanalysis.
Previous studies [10,36] have shown that the CERES EBAF surface product provides relatively accurate estimations of R s .Ma et al. [12] find that the CERES EBAF surface R s had a homogeneous distribution of bias of R s , compared with 446 field observations.These findings are confirmed by Figure 1 in this study.
Compared with the observations at 54 BSRN stations, the CERES EBAF surface R s data have an MAB of 9.9 W•m −2 (5.4%).Previous studies have shown that the CERES EBAF surface R s data product is among the best of the currently available global satellite R s retrievals [10,12,36].Hence, CERES EBAF surface R s data were used as references to correct the reanalyses in this study.Linear nearest-neighbor interpolations were applied, to make the spatial resolution of the reanalyses consistent with that of the CERES EBAF surface satellite retrievals with 1 • spatial resolution.

Reanalysis
The reanalysis products used in this study include the ERAI [43], MERRA2 (MERRA Version 2) [44], Climate Forecast System Reanalysis (CFSR) [45], Japanese 55-year Reanalysis (JRA55) [46], National Centers for Environmental Prediction/Department of Energy reanalysis (NCEP2) [47], National Center for Environmental Prediction/National Center for Atmospheric Research reanalysis (NCEP1) [48], European Center for Medium-Range Weather Forecasting 20th-Century Reanalysis dataset (ERA20CM) [49], and the National Oceanic and Atmospheric Administration/Cooperative Institute for Research in Environmental Sciences (CIRES) Twentieth Century Global reanalysis [50] (see Table 1).These reanalyses include four third-generation reanalysis datasets (ERAI, MERRA2, CFSR, and JRA55), one second-generation reanalysis dataset (NCEP2), one first-generation reanalysis dataset (NCEP1), and two reanalysis datasets for the 20th century (ERA20CM and CIRES).The shortwave radiation parameterization schemes of the eight reanalyses are different from one another, except for CFSR and CIRES, which are both based on a modified version of the Rapid Radiative Transfer Model for GCMs (RRTMG) (Table 1).However, most of the uncertainties in the calculated R s stem from model inputs, including cloud and aerosol parameters, rather than the radiative transfer models themselves [51].The maximum-random overlap method is used to describe the cloud overlap between the different vertical layers for the eight reanalyses, except for NCEP1 and NCEP2, which both used the random overlap method.Climatological aerosol values are used in most reanalyses except for NCEP1 and NCEP2.Among these reanalyses, only the MERRA2 assimilates transient aerosol information [52], including the inter-annual variations of atmospheric aerosols.
Conventional data, such as pressure, temperature, humidity profiles from weather stations, radiosondes, balloons, aircraft, and buoys were assimilated, except for ERA20CM and CIRES.ERA20CM only incorporates sea-surface temperature and sea-ice cover data in its boundary conditions.CIRES assimilates sea-ice concentrations, surface pressure, and sea-surface temperatures.The satellite data, including information of atmospheric temperature and humidity profiles, were assimilated in the eight reanalyses, except for ERA20CM and CIRES.
Satellite cloud motion winds, derived from geostationary meteorological satellites (Meteosat) operated by EUMETSAT, were only assimilated in ERAI, JRA55, and MERRA2.The TRMM Microwave Imager (TMI) radiance and its Total Column Water Vapor (TCWV) and rain rate satellite data were only used in JRA55, although JRA55 does not assimilate the data from the Atmospheric Infrared Sounder (AIRS), which are used by most reanalyses.

Method
The time span of the CERES satellite EBAF-surface estimated R s data is short, but can provide realistic R s values.The long-term R s records in the reanalyses exhibit substantial spatial bias.Merging satellite and reanalysis R s would be an effective way to generate long-term R s data with high accuracy.Here, we applied a CPDF method to the reanalysis R s data, a method that has been proposed to correct the bias in satellite precipitation retrievals using ground-based gauge data [53].The basic idea of the CPDF is to match the distribution of the reference data.Its fusion method combines the value distribution of R s data from satellite and R s records from the reanalyses [53].For example, when the 65 percentile values at the cumulative probability density function curve for the reanalysis estimates is 25 W•m −2 , the two nearest-neighbor satellite estimates for the 65 percentile values are 13 W•m −2 and 14 W•m −2 .Reanalysis estimates of 25 W•m −2 will be reduced to 13.5 W•m −2 , so that the cumulative probability density function of the bias-corrected reanalysis estimates will be the same as for the satellite data.
Figure 2 illustrates how the CPDF methods result in a specific grid point of the ERAI reanalysis.The data from the CERES EBAF surface product and the ERAI products over the same period were used to draw curves of the cumulative probability density function (PDF) for each grid point.Afterwards, a merging procedure was performed, by adjusting the PDF of the ERAI for the whole time period based on the CERES EBAF surface data after 2000.This procedure was conducted for each grid point.The PDF of monthly mean R s from satellite data can be used to correct the bias of reanalysis R s for the whole period, including the period with reference data available and that without reference data.This is the key advantage of the cumulative probability density function based method, compared to the simple correction method.
The performance of the CPDF method was evaluated using the R s observations collected at 54 BSRN stations from 1992 to 2016.The coefficient of determination (R 2 ), MAB, and standard deviation (Std) were also calculated, to assess the merge results.The MAB values were calculated as the mean of the absolute value of the difference between the reference values and the obtained data.MAB was used instead of bias because its mean bias, which is commonly used, has the disadvantage of the positive and negative values cancelling each other out when summarizing all the sites results.The Std here means the standard deviation of the difference between the reference values and the obtained data.
Afterwards, a merging procedure was performed, by adjusting the PDF of the ERAI for the whole time period based on the CERES EBAF surface data after 2000.This procedure was conducted for each grid point.The PDF of monthly mean Rs from satellite data can be used to correct the bias of reanalysis Rs for the whole period, including the period with reference data available and that without reference data.This is the key advantage of the cumulative probability density function based method, compared to the simple correction method.

Comparisons of the Reanalyses and CERES EBAF Surface Product
Figure 3 illustrates the mean bias of the R s from the eight reanalyses from 2000 to 2010, compared with the CERES EBAF surface satellite R s .The period of 2000 to 2010 was selected to cover a common period of all the estimates of R s (see Table 1).All the reanalyses contain significant negative and positive biases on a global scale, except for NCEP1, which exhibits a systematic bias of up to 50 W•m −2 for most areas.
Regionally, all reanalyses showed a positive bias in North America (2.82~39.11W•m −2 ) and East Asia (15.21~53.92W•m −2 ).On the western coast of South America, all reanalyses showed a positive bias (0.90~37.10 W•m −2 ), except for CFSR (−36.51W•m −2 ).MERRA2 (−27.00W•m −2 ), JRA55 (−2.53 W•m −2 ) and NCEP2 (−13.11W•m −2 ) had a negative bias over Indonesia.CFSR and CIRES showed a positive bias (29.63 W•m −2 and 22.74 W•m −2 ).Substantial biases in the reanalysis R s data were found at most marine strato-cumulus regions, particularly off the western coasts of the continents, which is likely due to bias in the cloud properties from the reanalysis systems.Moreover, the reanalyses applied a convection parameterization scheme to estimate marine strato-cumulus clouds, because of their coarse model spatial resolution [54,55].In contrast, the CERES EBAF surface product integrated the cloud observation data from CALIPSO-CloudSat merged cloud profiles and MODIS data, to improve its cloud detection accuracy.

Comparisons of the Reanalyses and CERES EBAF Surface Product
Figure 3 illustrates the mean bias of the Rs from the eight reanalyses from 2000 to 2010, compared with the CERES EBAF surface satellite Rs.The period of 2000 to 2010 was selected to cover a common period of all the estimates of Rs (see Table 1).All the reanalyses contain significant negative and positive biases on a global scale, except for NCEP1, which exhibits a systematic bias of up to 50 W•m −2 for most areas.

Evaluation of the Merged Surface Solar Radiation Data
Figure 4 shows a single site example of the comparison of the merged and raw NCEP1 data.This BSRN Billings site (latitude: 36.6N,longitude: 97.51W) is located in the Southern Great Plains of the United States, and the observation period was from 1993 to 2012.CERES EBAF surface data, raw NCEP1, and merged NCEP1 data are all highly correlated with the BSRN field observations, with an R 2 value of 0.95, due to the monthly records including seasonality well-simulated in the data.Compared with the CERES EBAF surface data, the raw NCEP1 data showed poor performance, with an MAB of 54.39 W•m −2 and a Std of 16.87 W•m −2 .However, the merged NCEP1 data performed better, with a lower MAB (12.30-12.56W•m −2 ) and Std (15.00-15.91W•m −2 ), for both the period after the year 2000 with the reference data and the period before 2000 without the reference data.
Merged reanalyses generally are better quality than raw reanalyses when comparing CERES EBAF surface data and BSRN data (Table 2).The MAB in merged renalyses were reduced by 21.62-72.74%using CERES EBAF surface data as reference, and 12.51-64.36%using BSRN as reference, compared with the raw reanalyses.The mean bias was significantly reduced in merged renalyses, and the Std was also improved by 13.14-56.22%,using CERES EBAF surface data as reference, and 6.12-42.92%using BSRN data as reference.Further analyzing the performance of CPDF for the time period without CERES EBAF surface data, Figure 5 presents histograms of the statistics (i.e., MAB, R 2 and Std) of eight reanalyses, using all 54 BSRN stations as reference.These statistical results are divided into two time periods (the period from 1992 to 1999 and the period from 2000 to 2016).For periods 2000-2016 and 1992-1999, the merged reanalyses all show lower MAB than the raw reanalyses.MAB less than 10 W•m −2 are increased by 6.00-50.30%after 2000 and 37.01-51.90%before 2000 in merged reanalyses, compared with raw reanalyses.What is more, the merged reanalyses have almost same R 2 and Std compared with the raw reanalyses.Although satellite data for the period of 1992-1999 are not available, the CPDF method improves reanalysis data based on the satellite data available after 2000.       1 The prefix 'C-' means CPDF merged reanalysis data; MAB the mean of the absolute value of the differences.

Seasonal Variations in R s
In this section, we evaluate whether the CPDF method could improve the seasonal variations in the R s data.An example site with seasonal variations in R s in both the raw and merged reanalyses is shown in Figure 6.The BSRN Tateno site, used as an example to illustrate the effect of the fusion procedure on the seasonal variation in R s data, has data available from 1996 to 2016, and is located in Japan.Before the merging procedure, all reanalyses (except for CFSR) show positive R s biases, especially the NCEP2 and NCEP1 reanalyses in May, while positive R s biases from these reanalyses were eliminated after the merging procedure.The high R s values in summer from the NCEP1 and NCEP2 indicated that the R s variations in these reanalyses had inaccurate estimations of clouds.The biases from these reanalyses were reduced and the seasonal variations of R s were more consistent with the observations after the merging procedure.Moreover, we calculated the mean absolute bias (MAB) of the ratio of each seasonal R s mean to annual mean R s (Table 3), using BSRN as a reference; the results showed that the merged reanalyses have lower MAB of the ratio of each seasonal mean R s to annual mean R s , and indicated that CPDF can mostly reduce relative bias of seasonal mean R s to annual mean R s , except for the ERA20CM, whose MAB remained the same.
Here we used standard deviation (Std) to evaluate R s reanalyses data's ability to simulate the seasonal variations in the reference of the CERES EBAF surface monthly data.Using R s data from all 54 BSRN sites, we further compared the Std of the multiyear monthly mean R s for the raw and merged reanalyses.From Figure 7, one can infer that the coefficient of determination of the seasonal variations of R s between reanalysis and CERES EBAF surface data improved by 1.0-8.9%.This means the merged reanalysis are consistent with the CERES EBAF surface data in Std, which indicates that the CPDF method improves the seasonal variations in R s .The prefix 'C-' means CPDF merged reanalysis data; MAB the mean of the absolute value of the differences.The prefix 'C-' means CPDF merged reanalysis data; MAB the mean of the absolute value of the differences.

Annual Variations in R s
Further examination of the annual variation in R s from the merged and raw reanalyses in Figure 8 shows the time series at Tateno for each of the datasets.The raw CFSR data agree well with the observation data for annual variation of R s .We found that the merging procedure almost entirely maintains the raw CFSR variation pattern, and that the overestimations of R s (bias) in the reanalyses are eliminated by the CPDF method, especially for NCEP1 and NCEP2.

Annual Variations in Rs
Further examination of the annual variation in Rs from the merged and raw reanalyses in Figure 8 shows the time series at Tateno for each of the datasets.The raw CFSR data agree well with the observation data for annual variation of Rs.We found that the merging procedure almost entirely maintains the raw CFSR variation pattern, and that the overestimations of Rs (bias) in the reanalyses are eliminated by the CPDF method, especially for NCEP1 and NCEP2.To validate the performance of the CPDF method for long-term annual Rs data, a pixel-by-pixel validation was conducted for each reanalysis product at the annual scale (Figure 9).The MAB and standard deviation were calculated for each pixel, using the CERES EBAF surface data as a reference.After the merging procedure, the MAB generally decreased by 59.10-85.10%,and the standard deviations remain nearly the same.This finding shows that the merged reanalysis reduced the bias at the annual scale, with the annual variation for each merged reanalysis remaining nearly unchanged.To validate the performance of the CPDF method for long-term annual R s data, a pixel-by-pixel validation was conducted for each reanalysis product at the annual scale (Figure 9).The MAB and standard deviation were calculated for each pixel, using the CERES EBAF surface data as a reference.After the merging procedure, the MAB generally decreased by 59.10-85.10%,and the standard deviations remain nearly the same.This finding shows that the merged reanalysis reduced the bias at the annual scale, with the annual variation for each merged reanalysis remaining nearly unchanged.

Different Performance of Reanalyses
Different spatial biases in the reanalyses are mainly attributed to the quality and amount of cloud data, aerosol data, and atmospheric data in each reanalysis data.The reasons for the high biases in the first and second generations of reanalyses, such as NCEP1 and NCEP2 [56], are complex.Li [57] notice that the Lacis and Hansen parameterization used in the NCEP1 dataset produced an obvious error in water vapor absorption that resulted in an error in Rs of approximately 5-40 W•m −2 .Moreover, the directly-assimilated satellite data used in NCEP1 and NCEP2 introduced a large uncertainty in Rs for the error propagation and the underutilization of the observational data in the early assimilated system.Additionally, the impacts of aerosol are not taken into account in either the NCEP1 or NCEP2 products.These factors above lead to a high MAB in the NCEP1 and NCEP2.
Given these pitfalls, the third generation of reanalyses improves the shortwave parameterization methods, and provides accurate input data to generate realistic estimations of clouds properties,

Different Performance of Reanalyses
Different spatial biases in the reanalyses are mainly attributed to the quality and amount of cloud data, aerosol data, and atmospheric data in each reanalysis data.The reasons for the high biases in the first and second generations of reanalyses, such as NCEP1 and NCEP2 [56], are complex.Li [57] notice that the Lacis and Hansen parameterization used in the NCEP1 dataset produced an obvious error in water vapor absorption that resulted in an error in R s of approximately 5-40 W•m −2 .Moreover, the directly-assimilated satellite data used in NCEP1 and NCEP2 introduced a large uncertainty in R s for the error propagation and the underutilization of the observational data in the early assimilated system.Additionally, the impacts of aerosol are not taken into account in either the NCEP1 or NCEP2 products.These factors above lead to a high MAB in the NCEP1 and NCEP2.Given these pitfalls, the third generation of reanalyses improves the shortwave parameterization methods, and provides accurate input data to generate realistic estimations of clouds properties, aerosols, and atmospheric conditions.For instance, the ERAI added an improved absorption coefficient, increased the number of spectral bands to 6, and used an hourly calculation of radiance.CFSR adopted a rapid radiative transfer model (RRTM).JRA55 updated O 2 , O 3 , and CO 2 absorption coefficients.CFSR, JRA55, ERAI, and MERRA2 all apply variational satellite bias correction (VarBC) to reduce biases in the assimilated satellite observations.Climatological aerosol data are also used in these reanalyses (Table 1).In addition, the improved assimilated system is used in these third-generation reanalyses.The joint efforts of these improvements produce low MAB values in these latest reanalyses.The differences among the third generation of reanalyses' MAB might be due to the different assimilated data and the shortwave schemes.The reason for the outperformance of ERAI might be related to the improvements in the cloud scheme, with the relatively accurate water vapor data by 1D + 4D-Var rain assimilation methods [43].These water vapor data provide relatively accurate cloud values in ERAI, compared to those of other reanalyses.A few observations, which only use sea surface temperature and sea ice cover data as boundary conditions and forcing, assimilated in ERA20CM and CIRES might be attributed to the poor constraints in the R s parameterization to obtain accurate cloud data.
The differences of post-merging MAB values are mainly attributed to the ability of each reanalysis to simulate the spatial-temporal distribution of clouds, and the impact of the variation of aerosols [58] on R s .For instance, the ERA20CM data assimilated fewer observations and produced a high Std, which indicates a limited ability to simulate the spatial-temporal patterns of clouds.Therefore, the post-merging MAB values are large for ERA20CM.NCEP2, which refines the cloud algorithms that are produced, has a lower post-merging MAB than NCEP1 [47].The current reanalyses, except for MERRA2, primarily used climatology aerosols.Thus, another reason for the post-merging MAB results might be the impact of aerosols on long-term trends in R s values [59].Focusing on the impact of aerosols in this issue will be our subject of future work.

Comparison with Other Fusion Methods
Over the past few decades, researchers have focused on generating high-quality R s datasets, especially using the fusion method.Journée and Bertrand [60] merge the surface solar radiation product from the Meteosat Second Generation (MSG) and the solar measurements network of the Royal Meteorological Institute of Belgium.These authors find that kriging with the external drift method performs best compared with the mean bias correction, the interpolated bias correction, and ordinary kriging.However, long-term trends in R s cannot be derived by these methods, because of the short time span of the observation and satellite data.Moreover, these methods depend heavily on the selection of training observations, and have inhomogeneity problems due to different instruments in the record.Using these methods may be a good choice on a regional scale.
Another commonly used approach is regression [61][62][63].However, the regression method suffers from the regression dilution effect, which means that the regression slope tends to be biased toward zero, which causes an underestimated scaling factor [64].The underestimated scaling factor will in turn produce a "flat" time series.Erickson et al. [65] find that the CPDF method performs better than the linear regression method, because the latter overcorrects the warmest temperatures, due to the lower weight applied to the scaling factor and the higher weight applied to the intercept parameter.The regression dilution effect causes the model data to shift closer to the average value (i.e., reduces the variance), which is also found in [66].Moreover, Alvarez-Garreton et al. [67] find that the CPDF method performs better than the regression method, because the latter ignores extreme values and underestimates the peak values, while the CPDF method can maintain the information provided by extreme values in the observations by assuming a non-linear relationship between the datasets.
Previous studies [33,62] have developed bias correction methods for R s , which are primarily dependent on the cloudiness to clear-sky ratio.Zhao et al. [68] provides a bias correction method using the error ratio.However, the bias in clear-sky radiation would be high, and the empirical coefficients in these algorithms only work well in the small regions where the R s observations are acquired.Additionally, the stability of the coefficients is unclear without the support of field observations.Thus, the CPDF method can provide reliable R s values, without empirical coefficients and other additional data that reduce error propagation.

Advantages and Uncertainties
The proposed CPDF can generate long-term R s data, by merging long-term reanalysis R s and high-quality satellite data.In this study, we assumed that the distribution of R s from satellite is true and that the distribution of R s did not change both before 2000 and after 2000.Therefore the biases in reanalyses can be reflected by the different R s distributions, compared with that of the satellite.The results show that the merged reanalyses show lower MAB than the raw renalyses.The merged R s reanalyses before 2000 were generated by using the distribution of R s from CERES EBAF surface data.CPDF also showed a good performance for the time period before 2000.This is because the CPDF method reduces the system bias in the raw reanalyses, and the trend in the R s data is maintained during the fusion process.The CPDF only changed the distribution of R s values, but kept the sequence of R s values from the original dataset.The results also showed that seasonal variability in the raw reanalysis is improved by the CPDF method; however, the improvement in annual variability is not significant compared with the results of the seasonal variability.Moreover, the MAB of the CERES EBAF surface product, compared with the BSRN stations, was 9.9 W•m −2 , and this bias still exists in the merged results.More work needs to be done to further reduce the bias of the satellite and merged R s data.
The validation is based on a direct comparison with field observations, which might lead to uncertainties.The site-specific observations used to validate the coarse R s grid would lead to potential errors.Hakuba et al. [69] notice that either a site-specific correction factor or increasing the number of sites in a grid box would reduce the sampling bias.Huang et al. [70] find that the sampling error for a 5 × 5 km 2 area ranged from 1.4 to 8.1% on monthly to "instantaneous" timescales, respectively.This finding indicates that the sampling error would be small on a monthly timescale.Li et al. [71] also noted that the sampling error decreases when averaged over periods greater than 5 days.Although the observation representativeness error is small on a monthly time scale, it cannot be entirely neglected.Applying the CPDF method at on different temporal and spatial scales needs to be further investigated.

Conclusions
In this study, we compared the R s values from eight reanalyses, including ERAI, MERRA2, CFSR, JRA55, NCEP2, NCEP1, ERA20CM, and CIRES, to those of the CERES EBAF surface product.A CPDF method was developed to merge the eight reanalyses with the CERES EBAF surface product.Comparisons of the raw and merged reanalyses were performed for the periods before 2000 and after 2000, using BSRN surface observations.Our results show that the "new generation" of reanalysis data, including ERAI, MERRA2, CFSR, and JRA55 perform better, especially ERAI.Previous studies also confirm that these latest-generation reanalyses [17,36,[72][73][74] have good performances.The reason can be mainly attributed to the assimilation of new satellite data and atmospheric sounding data.Moreover, enhanced radiation transfer models have improved cloud parameterization schemes, which can reduce the cloud estimation random errors in these latest generation reanalyses.Substantial spatial R s biases exist in the eight reanalyses when compared to the CERES EBAF surface data.After using the CPDF method, the MAB values in theses eight reanalyses were reduced by 21.24-64.36%,for both after 2000 and before 2000.For the period of 1979 to the present, the merged ERAI data showed the best results (MAB = 10.40 W•m −2 , standard deviation = 14.46 W•m −2 and R 2 = 0.97).For the period of 1958 to present, the merged JRA55 data showed a reduction in the MAB of 32.87% to 11.17 W•m −2 .For the century-long period, the merged CIRES would be a better choice.The merged R s products can be obtained by contacting the corresponding author (Kaicun Wang, Email: kcwang@bnu.edu.cn).

Figure 1 .
Figure 1.Mean surface incident solar radiation (Rs) bias value among data for 54 Baseline Surface Radiation Network (BSRN) and Clouds and Earth's Radiant Energy System Energy Balanced and Filled (CERES EBAF) surface sites.Unit: W•m −2 .The biases exhibit a nearly homogeneous spatial distribution ranging from −10 to 10 W•m −2 .

Figure 1 .
Figure 1.Mean surface incident solar radiation (R s ) bias value among data for 54 Baseline Surface Radiation Network (BSRN) and Clouds and Earth's Radiant Energy System Energy Balanced and Filled (CERES EBAF) surface sites.Unit: W•m −2 .The biases exhibit a nearly homogeneous spatial distribution ranging from −10 to 10 W•m −2 .

Figure 2 .
Figure 2. European Reanalysis Interim (ERAI) data for one grid box, merged by the cumulative probability density function-based (CPDF) method.The grid box was randomly selected.The probability density function curves of original ERAI reanalysis, merged ERAI reanalysis, and CERES EBAF surface (CERES) are in blue, red, and green, respectively.A comparison of raw ERAI (blue), merged ERAI (red) and CERES EBAF surface data is on the right.MAB values were calculated using the mean of the absolute value of difference between the reference and observed data.The Std represents the standard deviation of the difference between the reference and data.The reference data is CERES EBAF surface data.Each point in the right subplot is monthly surface solar radiation values.

Figure 2 .
Figure 2. European Reanalysis Interim (ERAI) data for one grid box, merged by the cumulative probability density function-based (CPDF) method.The grid box was randomly selected.The probability density function curves of original ERAI reanalysis, merged ERAI reanalysis, and CERES EBAF surface (CERES) are in blue, red, and green, respectively.A comparison of raw ERAI (blue), merged ERAI (red) and CERES EBAF surface data is on the right.MAB values were calculated using the mean of the absolute value of difference between the reference and observed data.The Std represents the standard deviation of the difference between the reference and data.The reference data is CERES EBAF surface data.Each point in the right subplot is monthly surface solar radiation values.

Figure 3 .
Figure 3. Differences between the reanalysis data and the CERES EBAF surface data.Rs data were computed by averaging the monthly mean Rs for 2000-2010, at a spatial resolution of 1° for each dataset.

Figure 3 .
Figure 3. Differences between the reanalysis data and the CERES EBAF surface data.R s data were computed by averaging the monthly mean R s for 2000-2010, at a spatial resolution of 1 • for each dataset.

Figure 4 .
Figure 4. Validation of CERES EBAF surface data (CERES), raw reanalysis National Center for Atmospheric Research Reanalysis (NCEP1) data, and merged reanalysis (merged NCEP1) data using BSRN data at the Billings site (latitude: 36.6N,longitude: 97.51W).(a,b) raw CERES EBAF surface and NCEP1 data; (c,d) merged NCEP1 data before 2000 (without reference) and after 2000 (with reference).MAB values are calculated using the mean of the absolute value of the difference between the reference and observed data.The Std represents the standard deviation of the differences between the reference and observed data.

Figure 4 .
Figure 4. Validation of CERES EBAF surface data (CERES), raw reanalysis National Center for Atmospheric Research Reanalysis (NCEP1) data, and merged reanalysis (merged NCEP1) data using BSRN data at the Billings site (latitude: 36.6N,longitude: 97.51W).(a,b) raw CERES EBAF surface and NCEP1 data; (c,d) merged NCEP1 data before 2000 (without reference) and after 2000 (with reference).MAB values are calculated using the mean of the absolute value of the difference between the reference and observed data.The Std represents the standard deviation of the differences between the reference and observed data.

Figure 5 .
Figure 5. Histograms of the statistical parameters of the comparisons between merged and original reanalyses.In each sub-histogram, the horizontal axis represents statistical values, and the number of station is recorded on the vertical scale.The hollow red column represents the results merged using CPDF, and the light blue column is original results for each reanalyses.The MAB, R 2 (determination coefficient), and Std (stand deviation) at each station are calculated from monthly Rs data.Each column represents different reanalyses.The first and second row represent histograms of MAB with the reference (after 2000) and without the reference (before 2000) data.The third and fourth rows are the same, but for the histograms of R 2 at each site.The fifth and sixth rows are the histograms of Std at each site.

Figure 5 .
Figure 5. Histograms of the statistical parameters of the comparisons between merged and original reanalyses.In each sub-histogram, the horizontal axis represents statistical values, and the number of station is recorded on the vertical scale.The hollow red column represents the results merged using CPDF, and the light blue column is original results for each reanalyses.The MAB, R 2 (determination coefficient), and Std (stand deviation) at each station are calculated from monthly R s data.Each column represents different reanalyses.The first and second row represent histograms of MAB with the reference (after 2000) and without the reference (before 2000) data.The third and fourth rows are the same, but for the histograms of R 2 at each site.The fifth and sixth rows are the histograms of Std at each site.

Figure 6 .
Figure 6.Seasonal variation of the Rs at Tateno (36.06N, 140.13E).Obs is the BSRN estimation, and CERES represents the CERES EBAF surface estimation.Seasonal variations of Rs were computed by the multiyear average of the Rs for each month.

Figure 6 .
Figure 6.Seasonal variation of the R s at Tateno (36.06N, 140.13E).Obs is the BSRN estimation, and CERES represents the CERES EBAF surface estimation.Seasonal variations of R s were computed by the multiyear average of the R s for each month.

Figure 7 .
Figure 7.Comparison of the seasonal variability in the merged and raw reanalyses.Seasonal variability is described in using standard deviation (Std) of the multiyear mean monthly Rs from merged and raw reanalyses.In each scatterplot, the horizontal axis represents the Std from the CERES EBAF surface product, and the vertical axis represents the Std from the merged and original reanalyses.The red dots denote merged results using CPDF, and the blue dots denote original reanalyses.Each point in the figures corresponds to each BSRN site.

Figure 7 .
Figure 7.Comparison of the seasonal variability in the merged and raw reanalyses.Seasonal variability is described in using standard deviation (Std) of the multiyear mean monthly R s from merged and raw reanalyses.In each scatterplot, the horizontal axis represents the Std from the CERES EBAF surface product, and the vertical axis represents the Std from the merged and original reanalyses.The red dots denote merged results using CPDF, and the blue dots denote original reanalyses.Each point in the figures corresponds to each BSRN site.

Figure 8 .
Figure 8. Annual variation in Rs at Tateno (36.06N, 140.13E).Obs represents the BSRN estimation, and CERES represents the CERES EBAF-surface estimation.Annual variations in Rs were computed using the average of the Rs for each year.The left subplot represents the raw reanalyses before merging, and the right subplot represents the merged reanalyses.

Figure 8 .
Figure 8. Annual variation in R s at Tateno (36.06N, 140.13E).Obs represents the BSRN estimation, and CERES represents the CERES EBAF-surface estimation.Annual variations in R s were computed using the average of the R s for each year.The left subplot represents the raw reanalyses before merging (a), and the right subplot represents the merged reanalyses (b).

Figure 9 .
Figure 9. Mean absolute bias (MAB) and standard deviation (Std) between each reanalysis and the CERES EBAF surface product for annual mean Rs.In each subplot, the horizontal axis represents statistical values, and the vertical axis represents the normal probability density distribution function (Normpdf).The red line denotes the merged results using CPDF, and the blue line denotes the raw reanalysis.

Figure 9 .
Figure 9. Mean absolute bias (MAB) and standard deviation (Std) between each reanalysis and the CERES EBAF surface product for annual mean R s .In each subplot, the horizontal axis represents statistical values, and the vertical axis represents the normal probability density distribution function (Normpdf).The red line denotes the merged results using CPDF, and the blue line denotes the raw reanalysis.

Table 1 .
Reanalysis data used in this study.

Table 2 .
Summary of statistics of all 54 BSRN sites using CERES EBAF surface (CERES) and BSRN as reference.

Table 3 .
Summary of the mean absolute bias (MAB) of the ratio of each seasonal Rs mean to annual mean Rs of all 54 BSRN sites, using BSRN as reference.

Table 3 .
Summary of the mean absolute bias (MAB) of the ratio of each seasonal R s mean to annual mean R s of all 54 BSRN sites, using BSRN as reference.