Evaluation and intercomparison of multiple satellite-derived and reanalysis downward shortwave radiation products in China

ABSTRACT

Downward shortwave radiation (DSR) is a critical variable in energy balance driving Earth’s surface processes. Satellite-derived and reanalysis DSR products have been developed and continuously improved during the last decades. However, as those products have different temporal resolutions, their performances in different time scales have not been well-documented, particularly in China. This study intended to evaluate several DSR products across multiple time scales (i.e. instantaneous, 1-hourly, daily, and monthly average) and ecosystems in China. Six DSR products, including GLASS, BESS, CLARA-A2, MCD18A1, ERA5 and MERRA-2, were evaluated against ground measurements at Chinese Ecosystem Research Network (CERN) and integrated land-atmosphere interaction observation (TPDC) sites from 2009 to 2012. The instantaneous DSR of MCD18 showed a root mean square error (RMSE) of 146.02 W/m². The hourly RMSE of ERA5 (155.52 W/m²) was largely smaller than MERRA-2 (188.53 W/m²). On the daily and monthly scale, BESS had the most optimized accuracy among the six products (RMSE of 36.82 W/m²). For the satellite-derived DSR products, the monthly accuracy at CERN can meet the threshold accuracy requirement set by World Meteorological Organization (WMO) for Global Numerical Weather Prediction (20 W/m²).

KEYWORDS:

• Downward shortwave radiation
• product validation
• remote sensing products
• reanalysis products

1. Introduction

Surface downward shortwave radiation (DSR) is the total solar irradiance from 0.3–3.0 μm reaching the Earth's surface. As a critical component of the surface radiation budget (SRB), DSR is essential for biological, geophysical, and ecological processes (Liang et al. 2019). Moreover, it is an important input for land surface numerical models such as irrigation scheduling and crop growth simulation models (Shook and Pomeroy 2011; Jiang et al. 2019). Therefore, accurate DSR estimation and characterization of its regional and global variabilities are crucial for climate monitoring and ecosystem research.

Ground measurements provide the most accurate DSR observations (Hartmann et al. 2013); however, they are limited by sparse spatial distribution and coverage owing to the high costs of building and maintaining an observation platform (Jahani et al. 2017). Alternatively, satellite-derived and reanalysis products are available for global and regional DSR monitoring (Yang and Bright 2020). Satellite-derived DSR products play a critical role in accurate DSR estimations because satellite observations provide a unique perspective on radiation at the top of the atmosphere (TOA; Hilburn et al. 2015). Reanalysis products can provide long-term historical DSR estimates from data assimilation that combines geophysical process dynamics models with ground and satellite measurements (Besharat, Dehghan, and Faghih 2013; Zhang et al. 2016). Unfortunately, reanalysis and satellite-derived DSR products are subject to uncertainties such as inaccurate retrieval of clouds and aerosols (Huang et al. 2019), topographic effects (Ma et al. 2023), geolocation, radiometric errors in sensor measurements (Liang et al. 2010), and simplified retrieval algorithms (Letu et al. 2020). Therefore, comprehensive evaluations should be performed to elucidate the uncertainties and deviations in DSR estimates.

In 2011, the World Meteorological Organization (WMO) published accuracy requirements/recommendations for DSR products for four specific applications (WMO 2022; Table 1). The requirements are divided into three levels: goals, breakthroughs, and thresholds. In 2022, the Global Climate Observing System (GCOS) program published the 2022 GCOS Implementation Plan and the 2022 ECV Requirements for DSR (GCOS 2022); 1, 5, and 10 W/m² were set as the goal, breakthrough, and threshold, respectively. Huang et al. (2019) summarized the typical accuracy of DSR products for different time scales. For the instantaneous scale, the root means square error (RMSE) for a clear sky condition is below 50 W/m², while that of a cloudy sky condition is 60–140 W/m². The RMSE for daily DSR of high-resolution kilometer-scale products is ∼35 W/m² (∼10% of the mean DSR), and the monthly RMSE has been reported to be below 10 W/m². DSR records at different timescales are suitable for different applications and have distinct requirements (Table 1). Instantaneous DSR represents the inversion capability of the satellite. Hourly DSR records detail solar radiation changes within a day, which are required for nowcasting and numerical weather prediction (WMO 2022), diurnal variabilities of solar-induced fluorescence (Damm et al. 2010), and solar energy utilization in photovoltaic power plants (Voyant et al. 2012). The daily scale is the most common temporal scale for current DSR products (Slater 2016). Daily DSR estimates are required for agricultural meteorology and climate monitoring (Abraha and Savage 2008). Although monthly DSR records are not practical for precise climate monitoring, validation results on a monthly scale can disclose systematic biases in the algorithm (Hilburn et al. 2015) by reducing random errors. Monthly DSR records are also capable of detecting long-term climate trends and characterizing solar dimming and brightening (SDB; Yang et al. 2022). Annual DSR records are related to photovoltaic power, and their trends are frequently used to capture SDB and climate change (Tang et al. 2017; Fang et al. 2021).

Table 1. Quantitative requirements for DSR collected by WMO in weather, water, and climate-related areas (WMO 2022).

Applications	Spatial resolution			Temporal resolution			Uncertainty (RMSE)
	Goal	Break	Thres	Goal	Break	Thres	Goal	Break	Thres
Global Numerical Weather Prediction	10km	30km	100km	60m	3h	12h	1 W/m^2	10 W/m^2	20 W/m^2
Agricultural Meteorology	1km	5km	20km	24h	2d	7d	—	—	—
Nowcasting and Very Short-Range Forecasting	5km	15km	50km	60s	10m	60 m	1 W/m^2	10 W/m^2	20 W/m^2
Climate monitoring	25km	50km	100km	24h	2d	5d	5 W/m^2	6.5 W/m^2	10 W/m^2

China necessitates an examination of DSR product performance because of its heavy aerosol load accompanied by rapid urbanization and economic development (Li et al. 2018; Guo et al. 2019; Lou et al. 2019). The Qinghai–Tibet Plateau (QTP), located in western China, is characterized by varied topography and complex atmospheric conditions that substantially impact atmospheric circulation and the energy budget (Yang et al. 2008; Wang et al. 2022). Therefore, DSR product validation on the QTP should be given additional consideration. Numerous validations and assessments of DSR products have been conducted using ground measurements from China Meteorological Administration (CMA) network stations. For the DSR reanalysis products, Zhang et al. (2020) evaluated monthly records from ERA-Interim and MERRA-2 using surface measurements from the CMA network station from 2000 to 2009. They pointed out that both ERA-Interim and MERRA-2 products overestimate DSR in China, with MERRA-2 overestimation being more pronounced (18.28 W/m² bias for ERA-Interim and 43.84 W/m² bias for MERRA-2). Jiang et al. (2020) evaluated global, direct, and diffuse DSR from the ERA5 across China using CMA stations. They concluded that ERA5 underestimated hourly diffuse radiation by 43.08 W/m² but overestimated total solar and direct radiation by 30.87 and 73.95 W/m², respectively. Less work has been done for satellite-derived products. Zhang et al. (2022) evaluated daily and monthly DSR estimates of CERES using CMA stations and documented RMSE and MBE of 21.2 (31.6) W/m² and 4.8 (6.2) W/m² at the monthly (daily) scale, respectively. Wang et al. (2018) reported that CLARA-A2 and SARAH-E overestimated monthly global solar radiation by 10.0 and 7.5 W/m² compared with ground measurements from CMA stations for 1993–2015 and 1999–2015, respectively. Nevertheless, most of these validation studies have been conducted at certain temporal scales and with little focus on ecosystems.

However, determining DSR accuracy across different temporal scales is crucial for meeting the requirements of various applications and research studies (Huang et al. 2019). In addition, potential changes in the land-atmosphere feedback of ecosystems may create major uncertainties in the energy budget (Heffernan et al. 2014) and exhibit different DSR features in various ecosystems. Therefore, it is necessary to conduct product evaluations in different ecosystems.

With a 1 h temporal resolution for DSR records and ecosystem monitoring, the Chinese Ecosystem Research Network (CERN) is an ideal alternative to the CMA network to conduct multi-temporal-scale DSR product evaluation. As an ecosystem research network, CERN provides valuable long-term ecological research by answering the question, ‘How can ecosystem structures be both a cause and consequence of ecological fluxes of energy?’(Li et al. 2015; Mirtl et al. 2018) Integrated land-atmosphere interaction observations (hereafter referred to as TPDC) can provide ground measurements in the QTP and elucidate the performance of DSR products.

Using CERN and TPDC data from 2009 to 2012, we aimed to validate and compare the multi-temporal-scale (hourly, daily, and monthly) and ecosystem-based performance of widely used surface parameter products and their newest DSR versions (GLASS, BESS, MCD18A1, CLARA-A2, ERA5, and MERRA-2) in China. Additionally, we studied their performance across different seasons and compared their interannual variation trends. Moreover, we investigated the spatial patterns of the products and accessed their quantitative discrepancies across China. The remainder of this paper is organized as follows. The DSR products and ground measurements are introduced in Section 2. The processing procedures, accuracy statistical metrics, and interannual analysis method used in this study are presented in Section 3. The results and analysis are presented in Section 4. The discussion and conclusions are presented in sections 5 and 6, respectively.

2. Materials

The latest versions of the most widely used DSR products were selected for this study. Table 2 lists the temporal and spatial resolutions, temporal ranges, and estimation methods for the DSR products. Details of each DSR product are presented in Section 2.1.

Table 2. Summary of the data sources used in the evaluation.

Category	Dataset	Temporal res.	Spatial res.	Temporal ext.	Method	Reference
Satellite-derived	GLASS	Daily	0.05°	2000 −2020	Direct estimation	(Zhang et al. 2019)
	BESS	Daily/monthly	0.05°/0.5°	2001–2020	Neural network (ANN)	(Ryu et al. 2018)
	CLARA-A2	Daily/monthly	0.25°	1982–2015	LUT	(Karlsson et al. 2017)
	MCD18A1	Instantaneous/3h	1 km	2000-present	LUT	(Wang 2021)
Reanalysis	ERA5	Hourly/monthly	0.25°/0.1°	1940-present	Assimilation (Ensemble of 4DVAR)	(Muñoz Sabater 2019)
	MERRA-2	Hourly/monthly	0.5°×0.625°	1980-present	Assimilation (3DVAR SI)	((GMAO) 2015)

2.1. Satellite-derived and reanalysis of DSR products

2.1.1. Global land surface satellite-DSR

Production of the Global Land Surface Satellite (GLASS) product suite began in 2009 (Liang et al. 2013a; Liang et al. 2013), and subsequent development enhanced the accuracy and quality of the products (Liang et al. 2021). The old version of GLASS DSR was enhanced by improving the look-up table (LUT) method using visible-band TOA radiance from multi-source polar-orbiting and geostationary satellite data (Zhang et al. 2014). However, this development was limited by (i) poor performance among bright objects, such as snow – and ice-covered areas, and (ii) high time consumption for processing and constructing long-term DSR products. Therefore, the GLASS DSR product (V60) was created from Moderate Resolution Imaging Spectroradiometer (MODIS) TOA spectral reflectance using a direct-estimation method (Zhang et al. 2019). GLASS DSR outputs daily average data with a spatial resolution of 0.05°, which can be obtained from http://www.glass.umd.edu/DSR/ (last accessed May 3, 2023).

2.1.2. Breathing earth system simulator radiation

The Breathing Earth System Simulator (BESS) is a simplified process-based model that couples atmospheric and canopy radiative transfers, canopy photosynthesis, transpiration, and energy balance. The BESS DSR products were first published in 2017 at 0.05°daily and 0.5° monthly resolutions (Ryu et al. 2018). DSR products were generated using a combination of an atmospheric radiative transfer model and an artificial neural network (ANN). The ANN was fed with various land and atmospheric data from MODIS. BESS DSR data from 2000 to 2019 can now be obtained from https://www.environment.snu.ac.kr/bess-rad (last accessed May 3, 2023).

2.1.3. CM SAF cloud, albedo and surface radiation-A2

The CLARA-A2 record provides cloud properties, surface albedo, and surface radiation parameters derived from the Advanced Very High Resolution Radiometer (AVHRR) sensor onboard the polar-orbiting NOAA and METOP satellites (Karlsson et al. 2017). Surface radiation products included daily and monthly average DSR estimates. The retrieval of DSR is an LUT-based method (Mesoscale Atmospheric Global Irradiance Code, MAGIC) proposed by Mueller et al. (2009) that uses satellite-derived radiances in the visible and near-infrared AVHRR satellite channels. The nowcasting SAF software (PPSv2014) was used to distinguish between cloudy and clear-sky pixels. Monthly and daily averages are available on a 0.25° × 0.25° global grid and cover 34 years (1982–2015). The data can be accessed free from https://wui.cmsaf.eu/safira/action/viewDoiDetails?acronym = CLARA_AVHRR_V002 (last accessed May 3, 2023).

2.1.4. MCD18A1

MCD18A1 is a high-resolution global DSR product founded by NASA and distributed by the NASA EOSDIS Land Processes DAAC. The first version of MODIS DSR was released in 2017 as Collection 6 and provided a spatial resolution of 5 km and a temporal resolution of 3 h from 2000 to 2017 (Wang 2017). In the 2021 version of MCD18A1 (Version 6.1), the spatial resolution was improved to 1 km, and certain programming errors in the earlier version were resolved (Wang 2021). An LUT-based retrieval strategy was selected as the primary operational algorithm to generate products, mainly from TOA satellite observations based on MODIS Terra and Aqua. The direct use of TOA data reduces data gaps and error accumulation compared with atmospheric-parameter-based (e.g. aerosol and cloud) DSR estimation strategies (Wang et al. 2020). Users can access the data at https://lpdaac.usgs.gov/products/mcd18a1v061/ (last accessed May 3, 2023).

2.1.5. ERA5

ERA5 is a reanalysis dataset obtained from the European Center for Medium-Range Weather Forecasts (ECMWF; Muñoz Sabater 2019). Based on the Integrated Forecasting System (IFS) Cy41r2, ERA5 benefits from 10 years of modeling and assimilation techniques (such as the Rapid Radiative Transfer Model shortwave scheme; Hersbach et al. 2020). The most notable improvement in the ERA5 dataset is its significantly higher spatial (0.25°) and temporal (hourly) resolutions {Dee et al. 2011 #314} compared with its precursors (ERA-15, ERA-40, and ERA-Interim). We used the ‘Shortwave Solar Radiation Downward’ variable from ‘ERA5 hourly data on single levels from 1940 to the present’ (hereafter referred to as ‘ERA5 DSR’). The ERA5 hourly DSR data were accessed from the Climate Data Store (CDS) at https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab = overview (last accessed May 3, 2023).

2.1.6. Modern-Era retrospective analysis for research and applications version 2

The Modern-Era Retrospective Analysis for Research and Applications Version 2 (MERRA-2) was produced by the Global Modeling and Assimilation Office (GMAO; Gelaro et al. 2017). MERRA-2 is a timely replacement for MERRA and has made significant progress toward GMAO’s goal of Earth System reanalysis (Kleist et al. 2009; Molod et al. 2015). MERRA-2 implements version 5.12.4, of the GEOS atmospheric data assimilation system (including the GEOS atmospheric model and GSI analysis scheme). The spatial resolution of the MERRA-2 DSR data is 0.5° × 0.625°, and the temporal resolutions are 1 h and monthly. The variable used in our experiment, Surface Downward Shortwave Flux (henceforth referred to as ‘MERRA-2 DSR’), can be downloaded from https://disc.gsfc.nasa.gov/datasets/M2T1NXRAD_5.12.4/summary?keywords = MERRA-2%20tavg1_2d_rad_Nx (last accessed May 3, 2023).

2.2. Ground measurements

We gathered 4 years (2009-2012) of ground measurements from 48 sites, including 31 from CERN (Fu et al. 2010), 6 sites from the TPDC (Ma et al. 2020), and 11 sites from the Global Energy Balance Archive (GEBA) network. The distribution of these sites is shown in Figure 1.

Figure 1. Geographical distribution of the observation stations used in this study from CERN (stations represented by circles, different colors indicate different ecosystems), TPDC (stations represented by orange triangles), and GEBA (stations represented by cyan rectangles).

Supported by the Chinese Academy of Sciences (CAS), CERN was established in 1988 (Fu et al. 2010). Each ground observation in CERN collects hourly total incoming radiation (305–2800 nm) using a pyranometer CM11 (Kipp & Zonen, Delft, The Netherlands), with an accuracy of 5%. After registration and application, we collected data for 31 sites from CERN, including 14 agricultural sites, 8 forest stations, 2 grassland sites, 4 desert sites, 1 marsh sites, and 2 lake ecosystem sites. For data privacy for the 1-h scale records of CERN, nine sites (ASA, CWA, DHF, LCA, LZD, NMD, YGA, AKA, and HJA) provided available hourly scale records; the other 22 sites were recorded on a daily scale.

In addition to CERN, TPDC was implemented to characterize the accuracy of the DSR product on the Tibetan Plateau. Most of the TPDC stations used four-component radiometers (Hukseflux NR01 and Kipp & Zonen CNR 1) with an accuracy of $\pm$10% and spectral range of 285−3000 nm (Ma et al. 2020). BJ station was an exception, as it uses CM21 (Kipp & Zonen, Delft, The Netherlands) with an accuracy of $\pm$2% and spectral range of 300−2800 nm. In the trend analysis section, GEBA (http://www.geba.ethz.ch; Wild et al. 2017) was included to illustrate DSR variance. Maintained by ETH Zürich, GEBA collects monthly means of DSR measurements globally. Estévez et al. (2022) estimated that the accuracy of monthly mean values was 5% in terms of relative random measurement error.

3. Methodology

3.1. Quality control of ground measurements

As ground measurements were used as true values to gauge the accuracy of the gridded products, their quality was crucial for the correctness of the validation results. Previous studies have shown that ground measurements obtained in China are limited with respect to reliability owing to technical failures and operation-related problems (Xu et al. 2008; Tang et al. 2010). Therefore, we applied a quality control (QC) procedure to hourly measurements (nine sites described in Section 2.2) based on the study performed by Espinar et al. (2012). The identifiable faulty data and inaccurate measurements were eliminated using the following three filters:

• $\mathrm{SEA}\le -3^{\circ },\mathrm{GH}{\mathrm{I}}_{\mathrm{GRD}}=0$,

• $-3^{\circ }\le \mathrm{SEA}\le 0^{\circ }\mathrm{\&}\mathrm{GH}{\mathrm{I}}_{\mathrm{GRD}}<0,\mathrm{GH}{\mathrm{I}}_{\mathrm{GRD}}=0$,

• $\mathrm{ERL}:0.03\mathrm{GH}{\mathrm{I}}_{\mathrm{TOA}}<\mathrm{GH}{\mathrm{I}}_{\mathrm{GRD}}<1.2\mathrm{BN}{\mathrm{I}}_{\mathrm{TOA}}\cos (\mathrm{SZA}{)}^{1.2}+50$,

where SEA is the solar elevation angle, $\mathrm{GH}{\mathrm{I}}_{\mathrm{GRD}}$ is the global horizontal irradiance on the ground, $\mathrm{GH}{\mathrm{I}}_{\mathrm{TOA}}$ is the global horizontal irradiance at TOA, and $\mathrm{BN}{\mathrm{I}}_{\mathrm{TOA}}$ is the normal beam irradiance at TOA.

In the first two filters, GHI_GRD night values or negative values at sunrise and sunset were identified and set to zero. Subsequently, the QC procedure for the BSRN proposed by the WMO was implemented. The extremely rare limit (ERL; De Miguel et al. 2001; Muneer and Fairooz 2002; Long and Dutton 2010) excludes records that cannot guarantee logical minimum and maximum values. Therefore, only data that passed the three limits were retained. Figure 2 shows the data before and after quality control, taking 1,095 days of data (from 2009 to 2011) from the CERN CWA site as an example. Certain unstable values at the day–night threshold (i.e. at sunrise and sunset) and anomalous records above the set threshold were removed.

Figure 2. 1095 days data from CWA site before and after the quality control (The colors indicate the value of hourly DSR, the vertical coordinates indicate the variation in radiation values over the 24 h of the day, and the horizontal coordinates indicate the variation in data for each day over the three years).

3.2. Spatio-temporal matching of products and ground measurements

The investigated DSR products (see Section 2.1) have spatial resolutions ranging from 1 km to 0.5° and temporal resolutions of instantaneous, 1-hourly, daily, and monthly. Therefore, the spatiotemporal matching of different products is a prerequisite for ensuring that all validated data are comparable.

For the spatial scale, the products were validated and compared at a uniform spatial resolution of 0.5° × 0.5°, as recommended by Huang et al. (2016). If the spatial resolution of the DSR product was finer than 0.5° × 0.5°, the station location (data point) was regarded as the window center, and the data within 0.5° × 0.5° windows were resampled into the target resolution. The results presented in this study refer to those aggregated to a grid of 0.5° × 0.5° unless otherwise mentioned.

For the temporal scale, one of the greatest limitations of sub-daily scale validation is the time mismatch between ground measurements and DSR products. In this study, the hourly ground measurements nearest to the instantaneous/hourly stamps of the DSR products were used as the ground reference. Daily and monthly observation records were generated using the following steps. On a daily scale, nine quality-controlled hourly measurements (see Section 2.2) were averaged into a daily scale if there were more than 20 valid observations per day. Moreover, the 22 stations that were initially provided on a daily scale were combined for validation. In other words, 9 plus 22 stations were included in the daily scale validation. On a monthly scale, at least six daily data points should be recorded every 10 days for 1 month (in February, daily data should be recorded for more than 18 days).

3.3. Statistical metrics

3.3.1. Evaluation with ground measurements

Following the aforementioned processing procedures, the DSR products were evaluated by comparing the surface measurements with their collocated grid values. R², MBE, relative MBE (rMBE), RMSE, and relative RMSE (rRMSE) were used to characterize the accuracy of the DSR estimates. The correlation coefficient is an important index; the closer it is to the ideal value of 1, the more ‘in phase’ the ground measurements and product estimates are. MBE indicates the average over- or underestimation of the product estimates compared with reference datasets, and RMSE indicates the uncertainty of the product estimates. The formulas for these indicators are as follows: (1) $\begin{array}{c}{R}^2={\displaystyle \frac{{\left(n{\sum }_{i=1}^n{e}_i\times m_i-{\sum }_{i=1}^n{e}_i\times {\sum }_{i=1}^n{m}_i\right)}^2}{\left[n{\sum }_{i=1}^n{e}_i^2-{({\sum }_{i=1}^n{e}_i)}^2\right]\left[n{\sum }_{i=1}^n{m}_i^2-{({\sum }_{i=1}^n{m}_i)}^2\right]}}\end{array}$(1) (2) $\begin{array}{c}MBE={\displaystyle \frac{1}{n}\sum _{i=1}^n(e_i-m_i),rMBE=\left|{\displaystyle \frac{bias}{\overline{m}}}\right|\times 100\text{\%}}\end{array}$(2) (3) $\begin{array}{c}RMSE=\sqrt{{\displaystyle \frac{1}{n}\sum _{i=1}^n{(e_i-m_i)}^2}},rRMSE=\left|{\displaystyle \frac{RMSE}{\overline{m}}}\right|\times 100\text{\%}\end{array}$(3) where $e_i$ is the product estimate of the DSR, $m_i$ is the ground radiation measurement, $n$ is the number of samples, and $\overline{m}$is the average of the ground measurements.

3.3.2. Product intercomparison

To explore the differences between products directly and continuously, we calculated the 2009 annual mean values of the DSR products instead of using discrete ground measurements as an intermediary. The difference was represented by calculating the deviation between the collocated annual mean values using R (linear correlation), rRMSD (relative root mean square deviation), and rMBD (relative mean bias deviation). The formulas for these indicators are as follows: (4) $\begin{array}{c}R={\displaystyle \frac{{\sum }_{i=1}^n(x_i-\overline{x})(\mathrm{y}-\overline{y})}{\sqrt{{\sum }_{i=1}^n(x_i-\overline{x})}\sqrt{{\sum }_{i=1}^n(y_i-\overline{y})}}}\end{array}$(4) (5) $\begin{array}{c}MBD={\displaystyle \frac{1}{n}\sum _{i=1}^n(x_i-y_i),rMBD=\left|{\displaystyle \frac{MBD}{\overline{m}}}\right|\times 100\text{\%}}\end{array}$(5) (6) $\begin{array}{c}RMSD=\sqrt{{\displaystyle \frac{1}{n}\sum _{i=1}^n{(x_i-y_i)}^2}},rRMSD=\left|{\displaystyle \frac{RMSD}{\overline{m}}}\right|\times 100\text{\%}\end{array}$(6)

Note that (r)RMSD and (r)MBD were calculated using formulas similar to those for (r)RMSE and (r)MBE listed in (equations 2(2) $\begin{array}{c}MBE={\displaystyle \frac{1}{n}\sum _{i=1}^n(e_i-m_i),rMBE=\left|{\displaystyle \frac{bias}{\overline{m}}}\right|\times 100\text{\%}}\end{array}$(2) and 3(3) $\begin{array}{c}RMSE=\sqrt{{\displaystyle \frac{1}{n}\sum _{i=1}^n{(e_i-m_i)}^2}},rRMSE=\left|{\displaystyle \frac{RMSE}{\overline{m}}}\right|\times 100\text{\%}\end{array}$(3) ). Herein, ‘Deviation (D)’ was used instead of ‘error (E)’ because in this case, one product ($x_i$)$$was compared with a product ($y_i$) other than the ‘true value'.

3.4. Interannual variation analysis

We used the Mann–Kendall (M–k) test and Sen's slope estimator to assess the trend of DSR in the time series data. The M–k test is a non-parametric method suitable for environmental, meteorological, and hydrological time series data. Sen's slope estimator is a robust method for estimating the magnitude and direction of trends. Statistical significance at the 95% confidence level was used to ensure high confidence in trend estimates. We employed deseasonalized monthly anomalies [DA(m)], defined as follows: (7) $\begin{array}{c}DA(m)=DSR(m)-\overline{DSR(m)}\end{array}$(7) where $DSR(m)$ is the monthly mean DSR in a year and $\overline{DSR(m)}$ is the corresponding averaged value. DSR anomalies were preferred to absolute DSR in the computation of the trends of the time series because the use of deseasonalized anomalies minimizes the seasonal variability of the component.

4. Results and analysis

4.1. Multiple-temporal scale evaluation

4.1.1. Sub-daily scale evaluation

At sub-daily scales, we separated the products into clear, cloudy, and all-sky conditions to evaluate their performances under different sky conditions. The division was based on the clearness index (difference between observed DSR and TOA). If the clearness index was > 0.5, the threshold used by Cao et al. (2016), the hourly DSR value was considered to indicate clear-sky conditions. Otherwise, it was treated as a cloudy-sky observation.

First, we compared the quality-controlled hourly DSR ground measurements of CERN with instantaneous estimates of MCD18A1. The evaluation results of MCD18A1 under different sky conditions are shown in Figure 3. Under all-sky conditions, R² was 0.75, RMSE was143.42 W/m², and MBE was 18.11 W/m². The results under clear-sky conditions show that (r)RMSEs were (21.39%) 139.11 W/m² with a slight underestimation (MBE of −3.03 W/m²). Under cloudy-sky conditions, the performance deteriorated substantially: (r)RMSEs were (75.82%) 152.17 W/m² with (r)MBE of (31.36%) 62.94 W/m².

Figure 3. Evaluation results of instantaneous DSR estimates of MCD18A1 against nine sites from CERN under different conditions from 2009–2012: all-sky, clear-sky, and cloudy-sky conditions. The warm colors represent the high density of samples, while the cool colors represent low density, and the following density scatter plots have the same color scheme.

Figures 4 and 5 show the validation results under different sky conditions for hourly DSR estimations from ERA5 and MERRA-2, respectively. Under all-sky conditions, the overall RMSE (MBE) for ERA5 was 155.52 (32.14) W/m². In terms of MERRA-2, the overall RMSE (MBE) was 188.53 (72.7) W/m². These results are consistent with those of previous studies (Liang et al. 2010; Huang et al. 2019; Jiang et al. 2020), thereby indicating that both the ERA5 and MERRA-2 products overestimate DSR in China, with the overestimation of MERRA-2 being more apparent.

Figure 4. As Figure 3 but for hourly estimates of ERA5.

Figure 5. As Figure 3 but for hourly estimates of MERRA-2.

For clear-sky conditions, the reanalysis products showed good agreement between the DSR estimates and ground-based measurements. The (r)RMSEs were (27.66%) 134.43, and (26.88%) 130.65 W/m² for ERA5 and MERRA-2, respectively. In terms of MBE, ERA5 tended to underestimate DSR with an MBE of −18.75 W/m²: MERRA-2 slightly underestimated DSR with an MBE of 3.65 W/m².

Under cloudy-sky conditions, performance deteriorated substantially. The (r)RMSEs were (126.33%) 184.85, and (175.50%) 256.07 W/m² for ERA5 and MERRA-2, respectively. A different degree of overestimation was observed for the reanalysis products; (r)MBEs were (78.73%) 115.20 and (126.73%) 184.91 W/m², for ERA5 and MERRA-2, respectively.

4.1.2. Daily scale evaluation

MERRA-2, ERA5, and MCD18A1 do not directly provide daily DSR values. Therefore, we aggregated the hourly DSR estimates of MERRA-2 and ERA5 and the 3-hourly DSR estimates of MCD18A1 to the daily DSR using a simple averaging approach. The daily DSR estimates calculated from the products were validated against CERN sites (Figure 6), TPDC sites (Figure 7), and all sites (Figure 8). The MBEs of ERA5 were 16.87 W/m² for all sites, and 16.58 and 18.60 W/m² for CERN and TPDC sites, respectively. The MBEs of MERRA-2 were 33.04 W/m² for all sites, and 33.91 and 27.96 W/m² for CERN and TPDC sites, respectively. The four satellite-derived products showed an apparent underestimation of the TPDC records. In CERN, GLASS and CLARA tended to overestimate DSR with an MBE of approximately 7 W/m², whereas BESS and MCD18A1 tended to underestimate DSR with an MBE of approximately −2 W/m². The discrepancy over CERN finally led to differences in the overall MBEs, which were 1.70, −1.42, 2.65, and −6.58 W/m² for GLASS, BESS, CLARA, and MCD18A1, respectively. The overall RMSE values ranging from small to large were 30.24 W/m² (BESS), 31.70 W/m² (CLARA), 35.70 W/m² (GLASS), 37.11 W/m² (MCD18A1), 44.22 W/m² (ERA5), and 60.29 W/m² (MERRA-2). Based on the results, the satellite-derived products performed better than the reanalysis products overall (Babar, Graversen, and Boström 2019). However, neither met the thresholds set by the WMO.

Figure 6. Evaluation results of daily DSR values against CERN network: (a) GLASS; (b) BESS; (c) CLARA; (d) MCD18A1; (e) ERA5; (f) MERRA-2.

Figure 7. Evaluation results of daily DSR values against TPDC network: (a) GLASS; (b) BESS; (c) CLARA; (d) MCD18A1; (e) ERA5; (f) MERRA-2.

Figure 8. Evaluation results of daily DSR values against ground measurements combing CERN and TPDC together: (a) GLASS; (b) BESS; (c) CLARA; (d) MCD18A1; (e) ERA5; (f) MERRA-2.

Figure 9 shows the geographical distribution of the statistical metrics. The size of the dots varies from small to large, corresponding to the value of rRMSEs, and the color of the dots varies from blue (the greater the underestimation, the bluer) to red (the greater the overestimation, the redder), corresponding to the value of MBE. The reanalysis products showed an overestimation for almost all stations, and ERA5 dominated MERRA-2 except for AKA, NMG, and MAWORS. Consistent with the hourly timescale validation, the performances of sites in southern China were poorer than those of sites in northern China. Satellite-derived products perform differently in various locations and ecosystems. The optimal retrieval regions were desert sites in Northwest China (MBEs < 10 W/m²) The poorest retrieval was observed over the QTP, with passive MBE of < −20 W/m² except for BESS, which outperformed the other products with smaller MBEs by approximately $\pm$10 W/m². GLASS and CLARA yielded different accuracies on the two sides of the Hu Huanyong Line (Hu Line, population density boundary; dotted red lines in Figure 9a and d). The sites west of the Hu Line (23 of 26 for GLASS and 21 of 26 for CLARA) mainly yielded positive MBEs, whereas MBEs east of the Hu Line tended to be negative. This coincidence may result from increasing aerosol loading accompanied by dramatic urbanization in southeastern China (Wang et al. 2022). As suggested by previous studies(Xia et al. 2006; Wang et al. 2018), increasing aerosol loading from pollutant emissions is responsible for reduced DSR and increased diffuse radiation in cloud-free skies. It is highly probable that the GLASS and CLARA algorithms do not adequately account for changes in aerosols, resulting in a systematic overestimation(Zhang et al. 2019).

Figure 9. Spatial distribution of accuracy of product DSR estimates at daily time scale at ground measurements (a) GLASS; (b) BESS; (c) ERA5; (d) CLARA; (e) MCD18A1; (f) MERRA-2.

4.1.3. Monthly scale evaluation

Among the six products used in this study, GLASS and MCD18A1 did not archive monthly DSRs; therefore, their monthly DSR data were calculated using the daily integrated DSR. Consistent with the daily scale evaluation, the monthly mean DSR of the product estimates was validated using all the selected ground measurements from CERN and TPDC. Figures 10–12 show the evaluation results for CERN, TPDC, and all the sites, respectively. The evaluation results were consistent with those obtained on a daily timescale, indicating that the DSR estimates from the reanalysis products were significantly overestimated. However, for satellite-derived product data, BESS and MCD18A1 tended to underestimate DSR, whereas GLASS and CLARA tended to overestimate it in China. In particular, the (RMSEs) MBEs of monthly DSR estimates were (20.54) 1.02 W/m², (15.94) −2.32 W/m², (20.16) 1.85 W/m², (23.44) −7.31 W/m², (24.32) 15.53 W/m², and (41.39) 32.57 W/m² for GLASS, BESS, CLARA, MCD18A1, ERA5, and MERRA-2, respectively. BESS maintains the optimal estimates in the monthly scale, particularly in TPDC, showing an MBE of −3.89 W/m² and RMSE of 20.06 W/m². At CERN, BESS did not significantly outperform the other three satellite-derived products (RMSE of BESS was approximately 2–3 W/m² lower than that of the other three satellite-derived products). Except for the MERRA-2 data, the RMSEs of the products were < 25 W/m². Notably, monthly DSR estimates of the satellite-derived products meet the WMO's ‘threshold’ standards for numerical weather prediction.

Figure 10. Evaluation results of monthly DSR values against CERN network: (a) GLASS; (b) BESS; (c) CLARA; (d) MCD18A1; (e) ERA5; (f) MERRA-2.

Figure 11. Evaluation results of monthly DSR values against TPDC network: (a) GLASS; (b) BESS; (c) CLARA; (d) MCD18A1; (e) ERA5; (f) MERRA-2.

Figure 12. Evaluation results of monthly DSR values against ground measurements combining CERN and TPDC together: (a) GLASS; (b) BESS; (c) CLARA; (d) MCD18A1; (e) ERA5; (f) MERRA-2.

The geographic distribution of the monthly mean values for each site is shown in Figure 13. The overall geographical distribution pattern is consistent with the daily average values shown in Figure 9. BESS maintains the best performance in QTP. Owing to the shortage of cloud detection algorithms, MERRA-2 was still overestimated compared with sites in southern China. The monthly estimates of most products are acceptable because random errors have been significantly eliminated.

Figure 13. Spatial distribution of accuracy of DSR product estimates at ground sites. (The same as Figure 9, but at a monthly time scale).

4.2. Ecosystem-based variation

DSR drives canopy photosynthesis, which is fundamental to regional energy and ecological processes (Ryu, Berry, and Baldocchi 2019). A previous study reported that radiation sources could result in a 20% uncertainty in yield estimates (Zhang et al. 2022). In this study, we conducted a preliminary evaluation of the performance of DSR products, focusing on ecosystem characteristics, by averaging over stations from diverse ecosystems (agriculture, forests, lakes, grasslands, and marshes). As shown in Figure 14, a significant overestimation of the desert ecosystem was observed; DSR in agriculture, forest, and lake ecosystems tended to be overestimated, while grasslands and marshes were underestimated. In addition, the performance of CLARA was optimal in certain ecosystems (deserts and marshes), whereas BESS yielded the most optimized estimations for agriculture, forests, lakes, and grasslands.

Figure 14. Comparison of the MBD and rRMSE values obtained from different sets of ecosystems.

Although CERN provides valuable insights into ecosystems, only a limited number of sites from grasslands, marshes, and lake ecosystems are available, which may cause a lack of representation of these ecosystems. Therefore, we made another attempt to characterize the influence of ecosystem heterogeneity on DSR estimation using GLASS, BESS, and MCD18A1 (chosen because of their fine resolution) as examples. We divided the selected DSR products’ 2009 annual average DSR estimates into nine primary land cover types (evergreen forest, deciduous forest, mixed forest, shrub, savannas, grasslands, cropland, bare land, and water) according to MCD12C1. Figure 15 shows the dispersion and outliers of the 2009 annual average DSR of different land-cover classes as boxplots, with whiskers representing the 25th (lower) and 75th (upper) percentiles. The selected products have diverse dispersions over different land covers but show similarities within a certain type of land cover. For example, the DSR estimates for the deciduous forest, cropland, and bare land were spread out on a relatively wide scale compared with the boxplots corresponding to evergreen forest, shrubland, and water, without much separation between the different products.

Figure 15. Box plots of the 2009 annual mean of different products over different land cover types; the red line within the box represents the median value; the blue markers represent the outliers.

4.3. Seasonal and interannual variation

4.3.1. Seasonal variation

We investigated the seasonal dependency on the accuracy of the products by dividing monthly DSR estimates into spring (MAM), summer (JJA), autumn (SON), and winter (DJF) seasons. As shown in Figure 16, BESS yields optimal estimates during spring, summer, and winter, whereas CLARA slightly outperforms BESS during autumn. Figure 16 also shows their consistent seasonal accuracy: the lowest R² and highest rMBE and rRMSE in winter and a relatively higher R² and lower rMBE and rRMSE in the other seasons. Zhang et al. (2016) arrived at a similar conclusion. Winter is distinguished from other seasons by (i) inaccuracy in DSR retrievals over snow-covered surfaces, (ii) degraded retrieval quality under the low sun and extensive viewing angle conditions in the Northern Hemisphere during winter, and (iii) low absolute DSR levels. Furthermore, disregarding that sub-monthly aerosol variability in satellite retrievals may contribute to the substantial relative bias during winter when China is assumed to have significant aerosol concentrations owing to domestic heating.

Figure 16. variations of R²; (b) seasonal variations of rMBE; (c) seasonal variations of rRMSE. DJF represents December, January, and February. MAM means March, April, and May. JJA represents June, July, and August. SON represents September, October, and November.

4.3.2. Interannual variation

First, we utilized surface measurements from 43 sites from CERN, TPDC, and GEBA between January 2009 and December 2012 to analyze site-averaged time trends of all products. Figure 17 shows the site-averaged time series of site-measured and product-estimated DSR monthly anomalies. Note that, with only 4 years of data, our purpose was not to accurately represent interannual variation in DSR, but instead to examine whether the products could replicate the anomalies of the given ground measurements. The correlation coefficients of the calculated and observed anomalies and the standard deviations of the calculated and observed anomaly differences are shown at the bottom left of each plot.

Figure 17. Time series of monthly mean DSR anomalies from January 2009 to December 2012. The solid blue line shows the site-averaged DSR anomalies; The solid red line shows the averaged DSR anomalies of the grids overlying the sites for each product. Sen's slopes of the site-measured and product-estimated time series as dark blue and red dashed lines (** marked as the significance level p < 0.05). The product ensemble mean is denoted by the light blue dashed line with 2 standard deviation uncertainties represented by the shadings. The bottom left of each plot shows the correlation coefficient of calculated and observed anomalies and the standard deviation of their anomaly differences.

Five out of six products have correlation coefficients larger than 0.90. CLARA and MCD18A1 DSR had the highest correlations and lowest standard deviations relative to the observations. Standard deviations between estimated and observed MERRA-2 DSR estimates were 3.0 W/m², achieving the largest fluctuations among these six products. Regarding linear trends, the observed and estimated DSR anomalies exhibited similar slope values with the same negative sign. In conclusion, the time series of the estimated DSR anomalies were generally consistent with the observed DSR anomalies, displaying the capability of DSR products to reproduce monthly and annual variability.

Furthermore, grid-area-averaged DSR data from the six products were compared for long-term trends over China, both annually and across seasons. To maintain consistency among the products, we divide the period (2001-2020) into two segments: Period I, spanning from 2001 to 2015, and Period II, spanning from 2016 to 2020. By doing so, we examine the temporal changes of all six products over the longer time frame of Period I and compare the five products (excluding CLARA, which ends in 2015) in the more recent and shared time frame of Period II.

The annual results are presented in Figure 18. During Period I, CLARA exhibited a slightly increasing trend of 0.07 W/m² per year, while the other products showed decreasing trends ranging from −0.45 W/m² per year to −0.02 W/m² per year, resulting in a six-product mean of 0.3 W/m⁻² per year. Notably, the products still exhibited significant consistency from 2009 to 2012, which supports our previous conclusion of their consistency. However, over the longer term, the deviation between products could be significant. For example, large standard deviations between products were observed in 2001 (3.59 W/m²), 2002 (2.64 W/m²), and 2012 (1.48 W/m²). In Period II, all five products showed statistically insignificant trends in DSR, with GLASS showing a flat trend, BESS and MCD18A1 showing slightly increasing trends, as well as MCD18A1 and MERRA-2 displaying slightly decreasing trends.

Figure 18. Time series of grid-area-averaged DSR anomalies in China for different products; the product ensemble means with 2 standard deviation uncertainties represented by the shadings. In all panels, the period was divided into 2001–2015 (Period I) and 2016–2020 (Period II). Dash lines represent the Sen's slopes (** marked as the significance level p < 0.05) of the product estimated time series.

Decadal DSR trends for individual seasons are depicted in Figure 19, revealing trends varying from season to season in both Periods I and II. During Period I, a consistent negative trend (from −0.41 W/m² per year to −0.26 W/m² per year) was observed in autumn for all six products. However, in the other seasons, the trends of each product were less consistent and mostly nonsignificant. In Period II, consistently increasing trends were observed in spring and autumn while decreasing trends were observed in summer and winter. These opposing trends resulted in a corresponding overall nonsignificant annual trend in Period II.

Figure 19. As Figure 18, but for the seasonal trend.

Multiple studies have disclosed that neither satellite nor reanalysis retrievals can accurately reproduce the decadal variability of DSR (Yang et al. 2023). As reported by He, Wang, and Feng (2021), ERA-Interim and ERA5 failed to reproduce trends in DSR variations from 1979 to 2014. Additionally, Wang et al. (2015) found that both CERES EABF and MERRA recorded invalid long-term trends in China. CERES EABF's recorded decreasing trend from 1983 to 2007 was caused by inhomogeneity due to imperfect usage of cloud observations, while MERRA's predicted increasing trend after 1990 was found to be invalid due to the products ignoring the canceling effect of increased aerosol loading. Li, Wang, and Liang (2021) have found that CLARA and GLASS products exhibited inaccurate long-term global trends due to sensor degradation over the years. Our study revealed inconsistency in the long-term trend in China across these products, which casts doubt on their ability to capture the long-term variation in the region. However, additional research using ground truth is necessary to quantify the long-term trend's accuracy and identify the source of error.

4.4. Spatial distribution and intercomparison

4.4.1. Annul mean DSR distribution

Analyzes based on ground measurements showed that DSR products have varying accuracies and systematic biases across different networks; therefore, a detailed spatial comparison is necessary. In this study, the average annual mean values of all DSR products were calculated. As shown in Figure 20a–f, the DSR estimates for the six products displayed a similar spatial distribution in China. In particular, DSR was higher in the non-monsoon region than in the monsoon region and higher in areas with higher elevations than in areas with lower elevations. Nationwide, maximum DSR occurs on the QTP, which is related to increased atmospheric transparency. At the same time, the minimum is in southwestern China (Szechwan Basin and Guizhou), which is associated with more clouds and higher aerosol loadings.

Figure 20. Spatial distribution of annual mean DSR (W/m²) of (a) GLASS; (b) BESS, (c) CLARA; (d) MCD18A1; (e) ERA5; (f) MERRA-2 in China.

The datasets were compared directly by computing their differences. Owing to its optimal performance, BESS was regarded as a benchmark for subtracting other DSR products, as shown in the validation section. Figures 21a–f compare the BESS estimates with the other five product estimates for the 2009 annual mean DSR. The red (blue) colors of various tints indicate that BESS is higher (lower) than the other products by varying degrees. A striking discrepancy can be observed in Figure 21. Most areas are red for ERA5 (d) and MERRA-2 (e), indicating that their DSR values are more systematically significant than those of BESS. In contrast, the differences between the BESS and the three other satellite-derived products are highly region-related, being negative (blue) in high elevation areas (the QTP and Yunnan–Guizhou Plateau in the west) and positive (red) in low elevation areas (plains and hills in the east).

Figure 21. Spatial distribution of annual mean DSR difference(W/m²) of (a) CLARA - BESS; (b) GLASS-BESS; (c) MCD18A1-BESS; (d) ERA5-BESS; (e) MERRA-2-BESS in China.

4.4.2. Product intercomparison

The geographical distribution of products provides qualitative verification results, such as whether they are spatially homogenous or heterogeneous. The quantitative discrepancies between products were further displayed by comparing the collocated grid values at 0.5° spatial resolution and illustrating them in terms of R, (r)MBD, and (r)RMSD.

As shown in Figure 22, R ranged from 0.68–0.93, rMBD from 1.17% to 21.55%, and rRMSD from 5.222% to 23.90%. GLASS and CLARA yielded the most optimized consistencies (R, rMBD, and rRMSD were 0.93, 1.17%, and 5.22%, respectively). Generally, the DSR estimation errors of the four satellite-derived products were strongly correlated with each other, with R > 0.86, rMBD < 5%, and rRMSE < 11%. Differences between the satellite-derived and reanalysis products were significantly relatively large (ERA5 > 10%; MERRA-2 > 20%). The correlation coefficient between BESS and MERRA-2 was high (0.88), but rMBD was low (21.54%). This high correlation is probably attributable to the use of MERRA-2 reanalysis data to compute filled data gaps under cloudy-sky radiation (Ryu et al. 2018). A large rMBD reflects the intrinsic difference between satellite-derived and reanalysis products. These results enhance the findings of previous studies that reported a large discrepancy between satellite-derived and reanalsis products (Liang et al. 2010; Huang et al. 2019; Jiang et al. 2020). This study goes a step further by providing quantitative data.

Figure 22. Comparison for the six DSR products (a) R; (b) rMBD; (c) rRMSE.

5. Discussion

5.1. Spatial representativeness error analysis

Representativeness errors result in discrepancies between in situ measurements and collocated DSR product estimates (Hakuba et al. 2013; Huang et al. 2019). The current study illustrates that representativeness errors are strongly contingent on the spatial scale of satellite-derived products (Huang et al. 2016). However, these inferences regarding particular products are controversial. For example, Wang, Wang, and Xue (2021) validated the 3 h DSR at the original resolution (1 × 1 km) and 1° × 1° scale and concluded that the accuracy was almost unchanged. In contrast, Wang et al. (2020) showed that when instantaneous MCD18A1 was aggregated from 1 to 50 km, the RMSE decreased from 119.2–101.0 W/m². To determine this contradiction, we evaluated the MCD18A1 product at multiple temporal (instantaneous, daily, and monthly) and spatial (50, 25, 5, and 1 km) scales. Consistent with Wang, Wang, and Xue (2021), for the instantaneous scale of MCD18, RMSE was strongly correlated with the spatial scale, ranging from 160.43–146.02 W/m² (Table 3). When the spatial scale changed from 5 to 25 km, RMSE decreased from 155.43–147.02 W/m², indicating that the accuracy of MCD18 is possibly subject to the representativeness of the surface observation. On a daily timescale, the results showed a 7.8 W/m² discrepancy between the 1 and 50 km scales. Conversely, on a temporal scale larger than or equal to 1 month, the spatial scale had no impact on data accuracy.

Table 3. Comparison between mean DSR of MCD18A1 with the ground measurements collected from CERN stations from 2009 to 2012 at multi-spatiotemporal scales.

Temporal scales	Accuracy	Spatial scales
		50km	25km	5 km	1km
Instantaneous	R^2	0.74	0.74	0.72	0.71
	RMSE (W/m^2)	146.02	147.96	155.45	160.43
	MBE (W/m^2)	16.94	17.44	19.35	20.12
Daily mean	R^2	0.87	0.86	0.83	0.81
	RMSE (W/m^2)	33.34	35.2	39.18	41.14
	MBE (W/m^2)	−1.95	−2.19	−2.32	−2.22
Monthly mean	R^2	0.9	0.89	0.88	0.88
	RMSE (W/m^2)	18.98	19.45	20.32	20.54
	MBE (W/m^2)	−2.63	−2.84	−2.96	−2.86

As the spatial resolution is an important factor that determines the value of statistical indicators, it is difficult to conclude that BESS outperforms GLASS based merely on accuracy at a scale of 0.5° × 0.5°. Therefore, we performed further validation on the original scale (0.05°). GLASS yielded an R² of 0.85, RMSE of 35.03 W/m², and MBE of 7.83 W/m², while BESS yielded an R² of 0.86, RMSE of 33.24 W/m², and MBE of −1.05 W/m². From this perspective, GLASS is comparable to BESS in terms of RMSE and R². A similar conclusion was drawn by Li, Wang, and Liang (2021).

5.2. Cloud effect on accuracy at sub-daily scales

Previous research has found that clouds significantly affect both satellite-derived and reanalysis DSR retrievals, making them unsuitable in cloudy regions or seasons (Li, Wang, and Liang 2021; Zhang et al. 2022). The sub-daily evaluation results of this study under different sky conditions (Section 4.1) reinforce this conclusion. Under cloudy-sky conditions, MCD18A1, ERA5, and MERRA-2 showed acceptable accuracies with MBEs of −3.03, −18.75, and 3.65 W/m², respectively. Under cloudy-sky conditions, a varying degree of overestimation was observed for the MCD18A1, ERA5, and MERRA-2 products: the MBEs were 62.94, 115.20, and 184.91 W/m², respectively.

However, the main reasons for the overestimation were different for the satellite-derived product (MCD18A1) and the reanalysis products (ERA5 and MERRA-2). For MCD18A1, the 3D effects of clouds may have significantly affected the validation results. Huang et al. (2016) demonstrated that uncertainty originating from the 3D effects of clouds is particularly large for instantaneous high-resolution satellite-derived DSR products. For ERA5 and MERRA-2, the main reason for their positive MBE may be inaccurate cloud characteristics (especially underestimation of cloudiness; Du et al. 2022; Stamatis et al. 2022). This is consistent with evaluation at the daily scale, where sites in southern China (more affected by clouds and precipitation) generally performed worse than those in northern China (see Figure 9e and f).

5.3. Differences between BESS and other satellite products

BESS did not correlate with the three other satellite-derived products (rMBDs between BESS and GLASS, CLARA-A2, and MCD18A1 were 10.06%, 10.06%, and 10.59%, respectively). The differences were highly region-related, being lower in high elevation areas (i.e. the QTP and Yunnan–Guizhou Plateau in the west) and higher in low elevation regions (i.e. plains and hills in the east). Previous studies have indicated that the choice of elevation model may introduce large uncertainties in DSR modeling (Ma et al. 2022). In fact, BESS offers the most reliable estimation over QTP among the products investigated, with an MBE of −1.75 W/m² and RMSE of 36.82 W/m². However, GLASS, CLARA, and MCD18A1 showed a significant underestimation over the QTP, with MBEs of −30.47, −23.93, and −36.86 W/m², respectively. A previous study reported that DSR estimates of satellite-derived products over the QTP are usually less accurate than those for other regions (Xia et al. 2006; Wu and Fu 2011). Wang, Wang, and Xue (2021) noted that GLASS and CERES-SYN yield fluctuations across the QTP. Wang et al. (2018) attributed the underestimation by CLARA to neglecting reduced atmospheric scattering at higher elevations. However, we observed that BESS performed remarkably across the QTP. This may be because the neural network-based algorithm is distinct from other LUT-based approaches. This phenomenon could reveal differences in the inherent physical mechanisms between the neural network approach, LUT, direct estimation (Ryu et al. 2018), and different elevation treatments.

6. Conclusions and outlook

We evaluated the DSR estimates from six gridded products (GLASS, BESS, MCD18A1, CLARA-A2, ERA5, and MERRA-2) on a 0.5° × 0.5° spatial scale and multiple temporal scales based on quality-controlled ground-measured records of CERN and TPDC from 2009 to 2012. Under all-sky conditions, the instantaneous DSR estimates of MCD18 exhibited an RMSE of 146.02 W/m² and MBE of 16.94 W/m². The hourly DSR of ERA5 had more optimized accuracy (RMSE of 155.52 W/m²) than that of MERRA-2 (RMSE of 188.53 W/m².) The products showed a clear overestimation under cloudy-sky conditions (MBEs of 62.94, 115.20, and 184.91 W/m² for MCD18A1, ERA5, and MERRA-2, respectively). The validation results for the daily and monthly data were consistent. Among the six products, BESS showed the highest overall accuracy on daily and monthly mean scales, particularly at sites from the TPDC (MBE of −1.75 W/m²). However, under certain circumstances, CLARA and MCD18 may outperform BESS. For instance, CLARA showed the best performance in particular ecosystems (desert and marsh) and during a specific season (autumn), whereas MCDA1 showed the best consistency for interannual variation with ground measurements. For the current satellite-derived DSR products, the monthly accuracy at CERN meets the threshold accuracy requirement suggested by the WMO for the application of global numerical weather prediction (20 W/m²). In terms of spatial accuracy distribution, reanalysis products greatly overestimated the DSR at most sites; GLASS and CLARA overestimated DSR at sites east of the ‘Hu Line,’ indicating they were subjected to an urbanization effect accompanied by increased aerosol loading. A preliminary evaluation of ecological dependence revealed an overestimated tendency in desert agriculture, forest, and lake ecosystems, and an underestimation in grassland and marsh ecosystems.

The performance of the DSR products exhibited seasonal dependence (more optimized accuracy during winter than during other seasons). The selected products displayed consistency and the capability of replicating monthly DSR variability in China in the short term (2009–2012). However, inconsistencies were observed in the long-term trends (2001-2020) of DSR among the gridded products. Seasonal variations were also found in their long-term trends, in which the products only showed a consistent dimming trend in autumn, while differed in the other seasons in terms of direction and magnitude. These discrepancies in long-term trends can significantly impact Earth system simulations, potentially leading to differences in energy and carbon budget estimation (Chakraborty, Lee, and Lawrence 2022). In addition, they can affect solar power potential estimation and determine the effectiveness and confidence in designing the layout of solar photovoltaic systems (He et al. 2023). Therefore, we highlight the need to pay attention to dataset selection when conducting the long-term DSR trend analysis over China.

Furthermore, the 2009 annual mean values of the six products were compared by characterizing their spatial distributions and quantitative statistical metrics. The differences between BESS and the three other satellite-derived products were highly region-related, being relatively low in high elevation areas (QTP and the Yunnan–Guizhou Plateau in the west) and relatively high in low elevation regions (i.e. plains and hills in the east). The results further established significant discrepancies between satellite-derived products and reanalysis products (rMBD > 8.5%, rRMSD > 12.3%). The GLASS, CLARA, and MCD18A1 products were reasonably correlated, and GLASS and CLARA yielded the most optimized consistencies (R, rMBD, and rRMSD were 0.93, 1.17%, and 5.22%, respectively). At the same time, BESS was less correlated with the results of the three other satellite-derived products.

This validation can elucidate the performance of DSR products and provide guidelines for the selection of appropriate products for specific applications in China. Future studies using additional ground measurement data are required to ensure representativeness and reinforce our conclusions. Moreover, surface measurement data should be collected over a more extensive time range to account for intra-annual variability and long-term trends and provide users with more comprehensive and representative information. In addition to global DSR components, direct and diffuse components are required for various research applications such as concentrated solar power plants (Abal et al. 2017) and monitoring ecosystem carbon uptake (Zhou et al. 2021). Radiation products containing these components are increasingly available; however, few related studies have been conducted (Jiang et al. 2020). Therefore, future evaluations should focus more on direct and diffuse components.

Acknowledgment

We would like to thank the satellite and reanalysis data developers and providers: ERA5 from ECMWF, MERRA-2 from NASA GMAO, MCD18A1 from NASA’s LP DAAC, CLARA-A2 from CM SAF, BESS from Ecological Sensing AI Lab at Seoul National University and the GLASS team. We are grateful to the Chinese National Ecosystem Research Network for their support. We also extend our special thanks to the 31 CERN stations represented by circles in Figure 1 for providing the necessary data. TPDC data are provided by National Tibetan Plateau Data Center. GEBA is co-funded by the Federal Office of Meteorology and Climatology MeteoSwiss within the framework of GCOS Switzerland. We appreciated anonymous reviewers’ and editors’ constructive comments in improving our paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by National Natural Science Foundation of China Grant (42090012), the Hubei Provincial Natural Science Foundation (2021CFA082), National Key Research and Development Program of China (2020YFA0608704), and the Fundamental Research Funds for the Central Universities through Wuhan University under Grant 2042022dx0001.