Merging ground-based sunshine duration observations with satellite cloud and aerosol retrievals to produce high-resolution long-term surface solar radiation over China

Although great progress has been made in estimating surface solar radiation (R_s) from meteorological observations, satellite retrieval, and reanalysis, getting best-estimated long-term variations in R_s are sorely needed for climate studies. It has been shown that R_s data derived from sunshine duration (SunDu) can provide reliable long-term variability, but such data are available at sparsely distributed weather stations. Here, we merge SunDu-derived R_s with satellite-derived cloud fraction and aerosol optical depth (AOD) to generate high-spatial-resolution (0.1^∘) R_s over China from 2000 to 2017. The geographically weighted regression (GWR) and ordinary least-squares regression (OLS) merging methods are compared, and GWR is found to perform better. Based on the SunDu-derived R_s from 97 meteorological observation stations, which are co-located with those that direct R_s measurement sites, the GWR incorporated with satellite cloud fraction and AOD data produces monthly R_s with R²=0.97 and standard deviation =11.14 W m⁻², while GWR driven by only cloud fraction produces similar results with R²=0.97 and standard deviation =11.41 W m⁻². This similarity is because SunDu-derived R_s has included the impact of aerosols. This finding can help to build long-term R_s variations based on cloud data, such as Advanced Very High Resolution Radiometer (AVHRR) cloud retrievals, especially before 2000, when satellite AOD retrievals are not unavailable. The merged R_s product at a spatial resolution of 0.1^∘ in this study can be downloaded at https://doi.org/10.1594/PANGAEA.921847 (Feng and Wang, 2020).

1 Introduction

A clear knowledge of variations in surface solar radiation (R_s) is vitally important for an improved understanding of the global climate system and its interaction with human activity (Jia et al., 2013; Myers, 2005; Schwarz et al., 2020; Wang and Dickinson, 2013; Wild, 2009, 2017; Zell et al., 2015). Direct measurements have shown that R_s has significant decadal variability, namely a decrease (global dimming) from the 1950s to the late 1980s and subsequent increase (global brightening) (Wild, 2009). The variation in R_s is closely related to the Earth's water cycle, the whole biosphere, and the amount of available solar energy. This situation emphasizes the urgency to develop reliable R_s products to obtain the variability in R_s.

Great progress has been made in the detection of variability in R_s by meteorological observations, satellite retrieval, and radiation transfer model simulations or reanalysis R_s products in previous studies (Rahman and Zhang, 2019; Wang et al., 2015). However, each estimation has its advantages and disadvantages. Directly observed data provide accurate R_s records at short timescales; however, careful calibration and instrument maintenance are needed to maintain their long-term homogeneity. Previous studies have reported that directly observed R_s over China may have major inhomogeneity problems due to sensitivity drift and instrument replacement (Wang, 2014; Wang et al., 2015; Yang et al., 2018). Before 1990, the imitations of the USSR pyranometers had different degradation rates of the thermopile, resulting in an important sensitivity drift. To overcome the radiometer ageing problem, China replaced its instruments from 1990 to 1993. However, the new solar trackers failed frequently and introduced a high data missing rate for the direct radiation component of R_s (Lu and Bian, 2012; Mo et al., 2008). After 1993, although the instruments were substantially improved, the Chinese-developed pyranometers still had high thermal offset with directional response errors, and the stability of these instruments was also worse than that of the first-class pyranometers recommended by the World Meteorological Organization (WMO) (Lu et al., 2002; Lu and Bian, 2012; Yang et al., 2010). Yang et al. (2018) show that nearly half of observed R_s values (60 out of the 119 R_s observation stations) have inhomogeneity issues. These artificial change points in observed R_s are mainly caused by instrument change (42 shifts), station relocation (34 shifts), and observation schedule change (20 shifts). The remaining 64 change points could not be identified.

Sunshine duration (SunDu) data are relatively widely distributed and have a long-term record (Sanchezlorenzo et al., 2009; Wild, 2009). Existing studies have also confirmed that SunDu-derived R_s data are reliable R_s data, which can capture long-term trends of R_s and reflect the impacts of both aerosols and clouds at timescales ranging from daily to decadal (Feng and Wang, 2019; Manara et al., 2015; Sanchezlorenzo et al., 2013; Sanchezromero et al., 2014; Tang et al., 2011; Y. Wang et al., 2012; Wild, 2016). Even though SunDu data do not provide a direct estimate of R_s and have different sensitivity to atmospheric turbidity changes, compared with R_s observations, they are still a good proxy for variations in R_s (Manara et al., 2017).

Sunshine duration observations collected at weather stations in China have been used to reconstruct long-term R_s (Che et al., 2005; Feng et al., 2019; He et al., 2018; He and Wang, 2020; Jin et al., 2005; Shi et al., 2008; Yang et al., 2006, 2020). Based on the global SunDu-derived R_s records, He et al. (2018) found that SunDu permitted a revisit of global dimming from the 1950s to the 1980s over China, Europe, and the USA, with brightening from 1980 to 2009 in Europe and a declining trend in R_s from 1994 to 2010 in China. Wang et al. (2015) found that the dimming trend from 1961 to 1990 and nearly constant zero trend after 1990 over China, as calculated from the SunDu-derived R_s, was consistent with independent estimates of AOD (Luo et al., 2001); they also observed changes in the diurnal temperature range (K. C. Wang et al., 2012; Wang and Dickinson, 2013) and the observed pan evaporation (Yang et al., 2015). Although direct observations and SunDu-derived R_s can provide accurate long-term variations in R_s, both direct observations and sunshine duration records are often sparsely spatially distributed.

Satellite R_s retrievals and radiation transfer model simulations or reanalysis R_s products can provide R_s estimation with global coverage at high spatial resolution. However, model simulations and reanalysis R_s products have substantial biases due to the deficiency of simulating cloud and aerosol quantities (Feng and Wang, 2019; Zhao et al., 2013). Previous comparative studies have shown that the accuracies of R_s from reanalyses are lower than those of satellite products (Wang et al., 2015; Zhang et al., 2016) due to the good capability of capturing the spatial distribution and dynamic evolution of clouds in satellite remote sensing data.

Table 1 lists the current satellite-based R_s products, which have been widely validated in previous studies. Zhang et al. (2004) found that the monthly International Satellite Cloud Climatology Project-Flux Data (ISCCP-FD) R_s product had a positive bias of 8.8 W m⁻² using Global Energy Balance Archive (GEBA) archived data as a reference. By comparing 1151 global sites, Zhang et al. (2015) evaluated four satellite-based R_s products, including ISCCP-FD, the Global Energy and Water Cycle Experiment-Surface Radiation Budget (GEWEX-SRB), the University of Maryland/Shortwave Radiation Budget (UMD-SRB), and the Earth's Radiant Energy System energy-balanced and filled product (CERES EBAF) and concluded that CERES EBAF shows better agreement with observations than other products. A similar overall good performance of CERES EBAF can also be found (Feng and Wang, 2018; Ma et al., 2015).

Table 1 Current satellite-derived surface solar radiation (R_s) products.

Although CERES EBAF uses more accurate input data to provide R_s data, its spatial resolution is only 1^∘ (Kato et al., 2018). Since 2010, new-generation geostationary satellites have provided opportunities for high-temporal- and high-spatial-resolution R_s data, such as Himawari-8 (Hongrong et al., 2018; Letu et al., 2020). However, the time span of the new-generation satellite-based R_s product is short. The long-term AVHRR records provide the possibility of building long-term radiation datasets. The CLoud, Albedo and RAdiation dataset, the AVHRR-based data-second edition (CLARA-A2), covers a long time period, but the spatial resolution is only 0.25^∘ (Karlsson et al., 2017). Recently, Tang et al. (2019) built a satellite-based R_s (SSR-tang) dataset using ISCCP-HXG cloud data. By using a variety of cloud properties derived from AVHRR, Stengel et al. (2020) presented the Cloud_cci AVHRR-PMv3 dataset generated within the Cloud_cci project.

Validation against the BSRN data indicated that SSR-tang values have the mean bias error (MBE) of −11.5 W m⁻² and root-mean-square error (RMSE) of 113.5 W m⁻² for the instantaneous R_s estimates at a 10 km scale, but Tang et al. (2019) point out that care should be taken when using this dataset for trend analysis due to the absence of realistic aerosol input data. Stengel et al. (2020) also show that R_s derived from Cloud_cci AVHRR-PMv3 reveals a very good agreement against BSRN stations, with low standard deviations of 13.8 W m⁻² and correlation coefficients above 0.98, and a bias for shortwave fluxes that is small (1.9 W m⁻²). However, default aerosol optical depth of 0.05 or data from Aerosol cci Level-2 or NASA MODIS Level-2 are used in the BUGSrad model to calculate clear-sky R_s, indicating that the impact of aerosols is not perfectly parameterized in Cloud_cci AVHRR-PMv3.

On the other hand, the long-term cloud records also contain uncertainties. For example, ISCCP cloud products, which directly combine geostationary and polar orbiter satellite-based cloud data, have large inhomogeneity due to different amounts of data from polar orbit and geostationary satellites and their different capabilities for detecting low-level clouds (Dai et al., 2006; Evan et al., 2007). This inhomogeneity of the cloud products might introduce significant inhomogeneity to the R_s values calculated from the cloud products (Montero-Martín et al., 2020; Pfeifroth et al., 2018b), and R_s long-term variability estimation still needs improvement.

Efforts have been made to further improve R_s products. Merging multisource data has become an effective empirical method for improving the quality of R_s products (Camargo and Dorner, 2016; Feng and Wang, 2018; Hakuba et al., 2014; Journée et al., 2012; Lorenzo et al., 2017; Ruiz-Arias et al., 2015). For instance, to produce spatiotemporally consistent R_s data, multisource satellite data are used in Global LAnd Surface Satellite (GLASS) R_s products (Jin et al., 2013). By merging reanalysis and satellite R_s data by the probability density function-based method, the reanalysis R_s biases can be substantially reduced (Feng and Wang, 2018). This finding suggests that fusion methods are effective ways to improve the estimation of R_s, especially when R_s impact factors are considered (Feng and Wang, 2019). Although linear regression fusion methods can produce R_s data incorporated with R_s impact factors, the stable regression parameters might have negative effects on the final fusion results due to the complex characteristics of R_s spatial–temporal variability.

On the other hand, the spatial resolution of R_s data is crucial for regional meteorology studies, as the minimum requirement of the spatial resolution of R_s data, as suggested by the Observing Systems Capabilities Analysis and Review of WMO OSCAR), is 20 km (Huang et al., 2019). Interpolation methods are often included in R_s fusion methods to further improve the spatial resolutions of R_s data (Loghmari et al., 2018). For example, Zou et al. (2016) estimated global solar radiation using an artificial neural network based on an interpolation technique in southeast China. By integrating R_s data from 13 ground stations with Meteosat Second Generation satellite R_s products, Journée and Bertrand (2010) found that kriging with the external drift interpolation method performed better than mean bias correction, interpolated bias correction, and ordinary kriging with satellite-based correction. However, interpolation results have uncertainties due to the lack of detailed high-spatial-resolution information. Although traditional linear regression fusion methods can incorporate high-spatial-resolution data during the fusion process, the impacts of the stable regression parameters need further investigation.

The performances of different machine learning methods have been evaluated in many previous studies, including simulation R_s at regional scale with support of satellite retrievals (Wei et al., 2019; Yeom et al., 2019) and site scale by using routine meteorological observations (Cornejo-Bueno et al., 2019; Hou et al., 2020). Whatever models or training data are selected, the impacts of the spatial relationship are not taken into account in these machine-learning-based models, and therefore a large number of input data are required to ensure accuracy.

Geographically weighted regression (GWR) is an extension of the traditional regression model by allowing the relationships between dependent and explanatory variables to vary spatially. Researchers have examined and compared the applicability of GWR for the analysis of spatial data relative to that of other regression methods (Ali et al., 2007; Gao et al., 2006; Georganos et al., 2017; LeSage, 2004; Sheehan et al., 2012; Q. Zhou et al., 2019). Due to the large spatial heterogeneity of R_s over China, the GWR method might produce accurate R_s variability estimations with an improved spatial resolution.

This study is established to merge SunDu-derived R_s data with satellite-derived cloud fraction (CF) and AOD data to generate high-spatial-resolution (0.1^∘) R_s over China from 2000 to 2017. The GWR and ordinary least-squares (OLS) regression merging methods are compared. CF and AOD are important R_s impact factors; however, many long-term R_s satellite products use climatology aerosol data as input. Whether much improvement is made in merging SunDu-derived R_s by incorporating AOD is also evaluated in this study, instead of evaluating direct merging current R_s products with SunDu-derived R_s. Since current high-quality R_s such as CERES EBAF have low spatial resolution, the output of this study provides reliable high-resolution grid R_s data to avoid the disadvantage of CERES EBAF for capturing the variability of R_s within a 1^∘ box and provides guidance to merge multisource data to generate long-term R_s data over China.

2 Data and methodology

2.1 Ground-based observations

2.1.1 Direct observations

R_s direct observations from 2000 to 2017 are obtained from the China Meteorological Data Service Centre (CMDC, http://data.cma.cn/data/index.html, last access: 3 March 2021) of the China Meteorological Administration (CMA). TBQ-2 pyranometers and DFY4 pyranometers have been used to measure R_s since 1993. Daily R_s values from 97 R_s stations are collected, and we calculated monthly R_s values by averaging daily R_s values when daily observed data are available for more than 15 d for each month at each radiation station. These monthly R_s values from direct measurements and collocated SunDu-derived R_s are used as independent reference data to investigate the performances of the fusion methods (Fig. 1). The whole area over China is further divided into nine zones by the K-mean cluster method based on geographic locations and multiyear mean R_s using 97 R_s direct observation sites, as shown in Fig. 1. The download instructions of the R_s direct observations can be found in Table 2.

Figure 1 The 2400 sunshine duration (SunDu) merging sites are shown as light brown points, and 97 independent validation sites, including R_s direct measurements and SunDu-derived R_s measurements, are shown as brown red points. The whole region is classified into nine subregions (I to IX) by the K-mean cluster method based on geographic locations and multiyear mean R_s using 97 R_s direct observation sites. The base hillshade map was produced by an elevation map of the study area using the global digital elevation model (DEM) derived from the Shuttle Radar Topography Mission 30 (SRTM30) dataset.

Table 2 Summary of availability information for all source data used in this study. CMDC is the China Meteorological Data Service Centre. SunDu is the sunshine duration data. R_s is surface solar radiation, and AOD is the aerosol optical depth.

2.1.2 SunDu-derived R_s

Sunshine duration observations (SunDu) and other meteorological data (e.g. air temperature, relative humidity, and surface pressure) from 1980 to 2017, which were collected from approximately 2400 meteorological stations (http://data.cma.cn/data/index.html) from the CMA, are used to calculate the SunDu-derived R_s (Fig. 1). R_s values are calculated following the method of the revised Ångström–Prescott equation (Eqs. 1–2) (He et al., 2018; Wang, 2014; Wang et al., 2015; Yang et al., 2006).

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{R_{\mathrm{s}}}{R_{\mathrm{c}}}}=a_{\mathrm{0}}+a_{\mathrm{1}}{\displaystyle \frac{n}{K}}+a_{\mathrm{2}}{\left({\displaystyle \frac{n}{K}}\right)}^{\mathrm{2}}\text{(2)}& {\displaystyle }{R}_{\mathrm{c}}=\int ({\mathit{\tau }}_{\mathrm{c}\mathrm{\_}\mathrm{dir}}+{\mathit{\tau }}_{\mathrm{c}\mathrm{\_}\mathrm{dif}})\times I_{\mathrm{0}}{d}_t\end{array}$$

Here n represents the measured SunDu, and K represents the theoretical value of the SunDu. a₀, a₁, and a₂ are the station-dependent parameters by tuning this equation with measurements of R_s and SunDu, and then the method is applied regionally (Wang, 2014). Instead of using observations from weather stations in Japan (Yang et al., 2006), observations in CMA are used (Wang, 2014). R_c is the daily total solar radiation at the surface under clear-sky conditions (Eq. 2). τ_{c_dir} and τ_{c_dif} represent the direct radiation transmittance and the diffuse radiation transmittance under clear-sky conditions. I₀ is the solar irradiance at the top of the atmosphere (TOA). For the clear-sky R_s, τ_{c_dir} and τ_{c_dif} are calculated using a modified a broadband radiative transfer model by simplifying Leckner's spectral model (Leckner, 1978), where the effect of transmittance functions of permanent gas absorption, Rayleigh scattering, water vapour absorption, ozone absorption, and aerosol extinction are parameterized using the surface air temperature, surface pressure, precipitable water, thickness of the ozone layer, turbidity, and sunshine duration as inputs (Yang et al., 2006). Calculation of R_s also includes impacts of aerosols because SunDu is impacted by changes in both clouds and aerosols (Wang, 2014).

Based on the classified subregions using 97 direct R_s observations in Fig. 1, the intercomparison results in Figs. 2 and 3 show that the agreement between SunDu-derived R_s and CERES EBAF R_s estimates is better than that between the direct observations and SunDu-derived R_s estimates, which is likely due to the inhomogeneity issue of direct R_s observations over China, as mentioned in many previous studies (Wang, 2014; Yang et al., 2018). The satellite R_s retrievals and SunDu-derived R_s are totally independent, but the high agreements of these two datasets indicate that they both are of higher accuracy. Similar results are also reported by Wang et al. (2015) that low agreement between SunDu-derived R_s and direct R_s observation is likely due to the directional response errors of the direct observations of R_s.

The SunDu-derived R_s observations, excluding SunDu observations located at direct observation sites, are used for merging. A total of 10 % merging observations are randomly selected for GWR parameter optimization. The download instructions of the SunDu observations can be found in Table 2.

Figure 2 Statistical summary of annual anomaly R_s from direct observed R_s, SunDu-derived R_s, and CERES EBAF R_s estimates in different subregions. The statistics include the mean absolute bias (MAB), standard deviation (SD), and root-mean-square error (RMSE). We use MAB due to the cancelling out effect of positive bias and negative bias. Nine subregions (I to IX) over China are shown in Fig. 1. S-D represents comparisons between SunDu-derived R_s and directly observed R_s. C-D represents comparisons between CERES EBAF R_s and directly observed R_s. S-C represents comparisons between SunDu-derived R_s and CERES EBAF R_s. The unit of the y axis is W m⁻².

Figure 3 Similar to Fig. 2, but this statistical summary is for R².

2.2 Satellite data

R_s data from the Clouds and Earth's Radiant Energy System energy balanced and filled product (CERES Synoptic (CERES) EBAF) surface product (edition 4.1) (Kato et al., 2018), cloud fraction from MODAL2 M CLD data product (Platnick et al., 2017), and AOD from the CERES SYN1deg edition 4A product (Doelling et al., 2013) are collected in this study. CERES EBAF R_s data are used as reference data. AOD from CERES SYN1deg and cloud fraction from MODAL2 M CLD are used as input data for fusion methods.

CERES is a three-channel radiometer measuring three filtered radiances, including shortwave (0.3–5 µm), total (0.3–200 µm), and window (8–12 µm). R_s values from CERES EBAF are adjusted using radiative kernels, including bias correction and Lagrange multiplier processes. The input data of CERES EBAF are adjusted during the product-generating process constrained by CERES observations at the TOA. The biases in temperature and specific humidity from the Goddard Earth Observing System (GEOS) reanalysis are adjusted by atmospheric infrared sounder (AIRS) data. Cloud properties, such as optical thickness and emissivity, from MODIS and geostationary satellites are constrained by cloud profiling radar; Cloud-Aerosol Lidar, and Infrared Pathfinder Satellite Observations (CALIPSO) detectors, and CloudSat. The uncertainties of CERES EBAF data, reported by Kato et al. (2018), in all sky global annual mean R_s is 4 W m⁻². Previous studies (Feng and Wang, 2019, 2018; Ma et al., 2015; Wang et al., 2015) have shown that the CERES EBAF surface product provides reliable estimations of R_s.

CERES SYN1deg AOD derived from an aerosol transport model, named Atmospheric Transport and Chemistry Modelling (MATCH) (Collins et al., 2001), which assimilates MODIS AOD data, is used to obtain spatiotemporally consistent AOD data. Different aerosol constituents, including small dust (<0.5 µm), large dust (>0.5 µm), stratosphere, sea salt, soot, and soluble, are used to compute the optical thickness for a given constituent optical thickness for a given constituent. We did not use AOD from MODIS, because MODIS AOD conation misses values and cannot meet the requirements of spatio-temporal continuity of AOD input in this study. In addition, MODIS AOD is only available under clear-sky conditions while AOD provided by the assimilation system is averaged under all conditions.

Cloud fraction data from MODAL2 M CLD are collected as input cloud fraction data with a spatial resolution of 0.1^∘ and time span from 2000 to 2017 (Platnick et al., 2017). The MODAL2 M CLD data are synthesized based on the cloud data from MOD06. Cloud fraction data from MOD06 are generated by the cloud mask product of MOD35 with a spatial resolution of 1 km. The MOD35 cloud mask is determined by applying appropriate single field-of-view (FOV) spectral tests to each pixel with a series of visible and infrared threshold and consistency tests. Each land type has different algorithms and thresholds for the tests. For each pixel test, an individual confidence flag is determined and then combined to produce the final cloud mask flag. The three confidence levels included in the cloud mask flag output are (i) high confidence for cloudless pixels (group confidence values >0.95), (ii) low confidence for unobstructed views on the surface (group confidence values Q≤0.66), and (iii) values between 0.66 and 0.95. Spatial and temporal continuity tests are further applied to determine whether the pixel is absolutely cloudless. Then, the cloud fraction is calculated from the 5 km × 5 km cloud mask pixel groupings; i.e. given the 25 pixels in the group, the cloud fraction for the group equals the number of cloudy pixels divided by 25.

2.3 Methods

2.3.1 Fusion models

OLS regression and GWR are used to build fusion methods for estimating R_s data. Cloud fraction and AOD have been important factors that affect variations in R_s. We compare different combinations of input data for the fusion methods, which can be classified into two types. The first type only contains cloud fraction data. The second type contains cloud fraction and AOD (Feng and Wang, 2020).

The OLS regression model is a commonly used model to estimate dependent variables by minimizing the sum-of-square differences between the independent and dependent variables. GWR is a regression model that allows the relationships between the independent and dependent variables to vary by locality (Brunsdon et al., 2010, 1998). GWR deviates from the assumption of homoscedasticity or static variance but calculates a specific variance for data within a zone or search radius of each predictor variable (Brunsdon et al., 1998; Fotheringham et al., 1996; Sheehan et al., 2012). The regression coefficients in GWR are not based on global information; rather, they vary with location, which is generated by a local regression estimation using subsampled data from the nearest neighbouring observations. The principle of GWR is described as follows:

$$\begin{array}{}\text{(3)}& y_i=\mathit{\delta }\left(i\right)+\sum _k{\mathit{\delta }}_k\left(i\right)x_{ik}+{\mathit{\epsilon }}_i,\end{array}$$

where y_i is the value of R_s unit i; $i=\mathrm{1},\mathrm{2},\mathrm{\dots },n,n$ denotes location i; x_ik indicates the value of the x_ik variable, such as cloud fraction and AOD; and ε denotes the residuals. δ_(i) is the regression intercept. δ_k(i) is the vector of regression coefficients determined by spatial weighting function w_(i), which is the weighting function quantifying the proximities of location i to its neighbouring observation sites; X is the variable matrix; and b is the bias vector.

$$\begin{array}{}\text{(4)}& {\mathit{\delta }}_k\left(i\right)={\left({\mathbf{X}}^{T}w\left(i\right)\mathbf{X}\right)}^{-\mathrm{1}}{\mathbf{X}}^{T}w\left(i\right)\mathit{b}\end{array}$$

The weighting functions are generally determined using the threshold method, inverse distance method, Gauss function method, and bi-square method. Due to the irregular distribution of observation sites and computer ability, the adaptive Gaussian function method is selected as a weighting function that varies in extent as a function of R_s observation site density.

$$\begin{array}{}\text{(5)}& w_{ij}=\exp (-(d_{ij}/\mathit{b}{)}^{\mathrm{2}}),\end{array}$$

where w_ij is the weighting function for observation site j that refers to location i; d_ij denotes the Euclidian distance between j and i; and b is the size of the neighbourhood, the maximum distance away from regression location i, called “band width”, which is determined by the number of nearest-neighbour points (NNPs).

2.3.2 GWR parameter comparison

To perform the local regression for every local area, the number of NNPs is required to estimate spatially varying relationships between CF, AOD, and R_s in the GWR-based fused method. To identify the best combination of parameter values, we test the number of NNPs ranging from 29 to 1000. A total of 10 % of merging SunDu-derived R_s data are randomly selected to validate these GWR parameters (Fig. 1). The results show that R² increases and bias decreases when the number of NNPs decreases. However, when the NNP is smaller than 30, the GWR-based fusion method produces spatially incomplete R_s data due to the local collinearity problem with large spatial variability. Therefore, 30 is selected as the NNP parameter (Table 3).

Table 3 Statistical summary of GWR parameter optimization. NNP is the number of nearest-neighbour points. GWR-CF presents the GWR-based fused method using only cloud fraction (CF) input, and GWR-CF-AOD presents that using both CF and aerosol optical depth (AOD) as input. MAB is the mean absolute bias. SD is the standard deviation. RMSE is the root-mean-square error.

3 Results

3.1 Site validation

Based on the independent SunDu validation sites, both the GWR and OLS methods explain 97 %–86 % of R_s variability (Fig. 4). The GWR method generally shows an improved performance compared with the OLS method due to the representativeness of the spatial heterogeneity relationship between R_s and its impact factors in GWR. Both the GWR and OLS methods produce better simulations of R_s if satellite and AOD data are incorporated.

Direct observations from 2000 to 2016 are also used to further evaluate the performance of the fusion methods (Fig. 4). The comparative result shows that both fusion methods show slightly reduced performances when using direct R_s observations rather than the SunDu-derived R_s. Both the GWR and OLS methods explain 91 %–82 % of R_s variability by using direct observations as reference data. Similarly, the GWR method exhibits better performances than the OLS-based fusion method, with an R² of 0.91 and root-mean-square error (RMSE) ranging from 19.89 to 19.97 W m⁻² at the monthly timescale (Table 4).

Figure 4 Comparison of surface solar radiation (R_s) derived from the GWR method and the OLS method. Panels (a), (c), (e), and (g) represent validation results using SunDu-derived R_s data as a reference. Panels (b), (d), (f), and (h) show directly observed R_s data. Panels (a)–(d) denote the GWR validation results, and panels (e)–(h) denote the OLS validation results.

Table 4 Validation of fusion methods driven by cloud fraction (CF) and AOD. GWR-CF and OLS-CF represent the GWR fusion method and OLS fusion method driven only by CF. GWR-CF-AOD and OLS-CF-AOD represent GWR and OLS fusion methods driven by CF and AOD, respectively.

3.2 Seasonal and annual variations in R_s

To analyse the impacts of AOD on the GWR fusion results, the GWR driven with only CF (GWR-CF) and GWR driven with CF and AOD (GWR-CF-AOD) are compared. Two validation sites (Changchun, 43.87^∘ N, 125.33^∘ E and Beihai, 21.72^∘ N, 109.08^∘ E) are randomly selected to evaluate the seasonal and annual variations in R_s derived from the GWR method (Fig. 5). The multiyear mean of AOD from Changchun and Beihai is 0.49 and 0.70, respectively. As shown in panels (a) and (b), both GWR-CF and GWR-CF-AOD produce similar seasonal variation patterns compared with SunDu-derived R_s and CERES EBAF R_s data. Small differences are found in the seasonal variation in R_s derived from GWR regardless of whether AOD was incorporated. Examination of the annual variation in R_s from the GWR-CF and GWR-CF-AOD is shown in panels (c) and (d) of Fig. 5. The two fusion methods also produce similar annual R_s variations. The similar performances of the GWR-CF and GWR-CF-AOD might suggest that the impacts of AOD have already been included in the SunDu-derived R_s site data.

Figure 5 Seasonal and annual variations in R_s at two sites: Changchun (a, c, 43.87^∘ N, 125.33^∘ E) and Beihai (b, d, 23.50^∘ N, 99.72^∘ E). SunDu R_s is the SunDu-derived R_s data, and GWR-CF R_s is R_s produced by the GWR method incorporating only the cloud fraction. GWR-CF-AOD is R_s produced by the GWR method incorporating cloud fraction and AOD. The multiyear mean of AOD from Changchun and Beihai is 0.49 and 0.70, respectively.

We also analyse the performances of fusion methods for different seasons at all validation sites, as shown in Table 4. At seasonal scales, both the GWR-CF and GWR-CF-AOD methods have high R² values ranging from 0.94 to 0.96, compared with direct R_s measurement or SunDu-derived R_s. GWR-CF and GWR-CF-AOD show slight differences, indicating that both fusion methods produce consistent R_s seasonal variation patterns, which might be because the impacts of AOD have already been included in the SunDu-derived R_s site data at seasonal timescales. Comparatively, the GWR methods perform best in autumn, with RMSEs ranging from 9.23 to 9.56 W m⁻² followed by winter, spring, and summer. Both the GWR-CF and GWR-CF-AOD methods produce similar annual variations in R_s from 2000 to 2016, with R² values ranging from 0.57 to 0.58 (Table 4). The statistics indicate that the GWR can produce reasonable seasonal and annual variations in R_s.

3.3 Multiyear mean and long-term variability in R_s

Figure 6 shows the performance of GWR-CF and GWR-CF-AOD in simulating the multiyear mean R_s by using 97 direct R_s observation sites and independent SunDu-derived R_s sites. Based on direct R_s measurements, both GWR-based methods show good performances with high R² (0.89–0.95) and low RMSE (11.03–11.11 W m⁻²), and few differences are found for the GWR merging results, whether or not AOD is taken as input data (Table 4).

Figure 6 Comparison of multiyear mean surface solar radiation (R_s) derived from the GWR method. Panels (a) and (c) represent validation results using SunDu-derived R_s data as a reference, while that of panels (b) and (d) uses directly observed R_s data.

The spatial distributions of the multiyear means of R_s from 2000 to 2017 are shown in Fig. 7. The SunDu sites show that R_s is high in northwest China, ranging from 180 to 300 W m⁻², and low in eastern China, ranging from 120 to 180 W m⁻². Both the GWR-CF and GWR-CF-AOD methods show consistent R_s spatial patterns with SunDu-derived R_s observations and CERES EBAFs, indicating that the relationship between R_s and impact factors is not linearly stable and is closely related to spatial position. The spatial distribution of the R_s trend derived from the GWR method is also consistent with the SunDu-derived R_s trend, especially in western China (Fig. 8). In order to prove that SunDu-derived R_s can add value to the 0.1^∘ product, instead of cloud fraction data alone, we perform a similar GWR analysis but using CERES EBAF interpolated to 0.1^∘ and 0.1^∘ cloud, and we compare the results with those using SunDu-derived R_s and 0.1^∘ cloud (Fig. 9). The results indicate that SunDu-derived R_s can add value to the 0.1^∘ product and the merged R_s by using interpolated CERES EBAF, and 0.1^∘ cloud products are also similar to the original CERES EBAF but not the 0.1^∘ cloud product.

Figure 7 Spatial distribution of multiyear mean monthly surface solar radiation (R_s) from 2000 to 2017. The first line (a, b) shows the observed multiyear mean monthly R_s from SunDu and CERES EBAF; the multiyear mean monthly R_s derived from the GWR method is shown in the second line (c, d), respectively.

Figure 8 Spatial distributions of monthly anomaly trends of surface solar radiation (R_s) from 2000 to 2017. The first line (a, b) shows the SunDu-derived R_s and CERES EBAF R_s; the R_s-derived GWR fusion methods are shown in the second line (c, d). Panel (c) incorporates only CF, and panel (d) incorporates CF and AOD. The black dots on the maps represent significant trends (P<0.05).

Figure 9 Spatial distributions of monthly anomaly trends of surface solar radiation (R_s) from 2000 to 2017. The first line (a, b) shows the SunDu-derived R_s and CERES EBAF R_s; the R_s-derived GWR fusion methods are shown in the second line (c, d). (c) SunDu-derived R_s incorporates only CF, and (d) is the results of GWR analysis using CERES data interpolated to 0.1^∘ and 0.1^∘ cloud. The black dots on the maps represent significant trends (P<0.05).

Based on the classified subregions using 97 direct R_s observations in Fig. 1, the regional means of R_s annual anomaly variation from 2000 to 2016 are shown in Fig. 10. Compared with observations, both the GWR-CF and GWR-CF-AOD methods produce consistent long-term R_s trends with SunDu-derived R_s and CERES EBAF R_s (Figs. 2, 3, and 10), indicating that the GWR-CF and GWR-CF-AOD methods can produce reasonable annual R_s variations over China.

In zones I and II, located in northern arid and semiarid regions, the annual anomaly R_s variation shows small fluctuations ranging from −10 to 10 W m⁻². In contrast, zones IV, V, VIII, and IX covering the Sichuan Basin, Yunnan–Guizhu Plateau, Qinghai–Tibet Plateau, and North China Plain show large R_s variation trends. Li et al. (2018) found a sharply increasing R_s trend over east China, especially in the North China Plain, which is due to controlling air pollution and reducing aerosol loading. However, our results indicate that the increased surface solar radiation in north China is not confirmed by satellite retrieval (CERES) and SunDu-derived R_s.

Figure 10 The regional mean of the annual anomaly of the surface solar radiation (R_s) for different subregions. Nine subregions (I to IX) over China are shown in Fig. 1. Direct R_s observations, SunDu-derived R_s, and CERES EBAF are shown as black lines, green lines, and red lines, respectively. Light and dark blue represent the R_s variation derived from the GWR-CF and the GWR-CF-AOD methods.

4 Discussion

4.1 Impact factors of R_s

In this study, we merged more than 2400 sunshine-duration-derived R_s site data with MODIS CF and AOD data to generate high-spatial-resolution (0.1^∘) R_s over China from 2000 to 2017. The results show that the GWR method incorporated with CF and AOD (GWR-CF-AOD) performs best, indicating the non-neglected role of clouds and aerosols in regulating the variation in R_s over China.

Clouds and aerosols impact the solar radiation reaching the surface by radiative absorption and scattering (Tang et al., 2017). Recent R_s trend studies over Europe suggest that CF may play a key role in the positive trend of R_s since the 1990s (Pfeifroth et al., 2018a). In terms of input data, our results also indicate that the cloud fraction might be a major factor affecting R_s, which is consistent with our previous studies (Feng and Wang, 2019).

Changes in aerosol loading have also been reported to be an important impact factor (Che et al., 2005; Li et al., 2018; Liang and Xia, 2005; Qian et al., 2015; Xia, 2010; Z. Zhou et al., 2019). The atmospheric visibility data show that the slope of the linear variation in surface solar radiation with respect to atmospheric visibility is distinctly different at different stations (Yang et al., 2017), implying that the relationship between R_s and aerosols varies with location.

4.2 Performances of the fusion methods

The good overall performances of the GWR model have been reported in many previous studies, including geography (Chao et al., 2018; Georganos et al., 2017), economics (Ma and Gopal, 2018), meteorology (Li and Meng, 2017; Q. Zhou et al., 2019), and epidemiology (Tsai and Teng, 2016). Chao et al. (2018) used the GWR method to merge satellite precipitation and gauge observations to correct biases in satellite precipitation data and downscale satellite precipitation to a finer spatial resolution at the same time. Q. Zhou et al. (2019) used GWR to analyse haze pollution over China and found that the GWR estimate was better than the OLS estimate, with an improvement in correlation coefficient from 0.20 to 0.75.

Compared with other traditional interpolation methods, such as optimal interpolation (OI), GWR can theoretically integrate geographical location and R_s impact factors for spatial R_s estimations and reflect the non-stationary spatial relationship between R_s and its impact factors. The thin plate spline method can include CF and AOD as covariates to simulate the approximately linear dependence of these impact factors on R_s, but this linear function cannot fully describe the relationship among CF, AOD, and R_s (Hong et al., 2005).

Comparison results from (Wang et al., 2017) also indicate that the GWR method is better than the multiple linear regression method and spline interpolation method for near-surface air temperature. By using the spatial interpolation method, CERES EBAF R_s can also be downscaled to 1 km or 30 m. These interpolated CERES R_s data cannot represent the detailed R_s distributions at spatial resolution of 1 km or 30 m due to the variability of R_s within a 1^∘ box. Without additional high-spatial-resolution data, interpolated data cannot capture more detail variability of R_s.

5 Code availability

The software code used for this article can be accessed at https://doi.org/10.5281/zenodo.4584571 (forgetbear, 2021).

6 Data availability

The merged R_s product by GWR methods with cloud fraction and AOD data as input in this study is available at https://doi.org/10.1594/PANGAEA.921847 (Feng and Wang, 2020). The cloud data can be downloaded from https://neo.sci.gsfc.nasa.gov/view.php?datasetId=MODAL2_M_CLD_FR (last access: 3 March 2021; NASA Earth Observations, 2021). The CERES SYN data can be downloaded from https://ceres.larc.nasa.gov/data/ (last access: 3 March 2021; CERES Science Team, 2021).

7 Conclusions

Accurate estimation of R_s variability is crucially important for regional energy budget, water cycle, and climate change studies. Recent studies have shown that SunDu-derived R_s data can provide reliable long-term R_s series. In this study, we merged SunDu-derived R_s data with satellite-derived cloud fraction (CF) and aerosol optical depth (AOD) data to generate high spatial resolution (0.1^∘) R_s over China from 2000 to 2017 (Feng and Wang, 2020). The GWR and OLS merging methods were also compared.

Our results show that the spatial resolutions of all fusion results are improved to 0.1^∘ by incorporating MODIS cloud fraction data. The GWR shows better performance than OLS, with increases in R² by 9.21 %–12.81 % and RMSEs reduced by 49.56 %–54.68 %, indicating that R_s has complex characteristics of spatial variability over China, which has also indicated the necessity of the high spatial resolution of R_s data. As clouds and aerosols play vital roles in the variability in R_s, apparent improvements in the results of SunDu-derived R_s data-merging are found if both cloud fraction and AOD are incorporated. Based on the merging results incorporating only cloud fraction, cloud fraction is suggested to be the major factor impacting R_s, which explained approximately 86 %–97 % of R_s variability. Generally, SunDu-derived R_s data merging results derived from GWR show more consistent multiyear mean R_s and long-term R_s trends compared with those from OLS. Our results show that the improvement in R_s variability estimation is closely related to R_s impact factors and R_s spatial heterogeneity. The merged R_s products derived from GWR-CF-AOD can be downloaded at https://doi.org/10.1594/PANGAEA.921847.