Significant discrepancies of land surface daily net radiation among ten remotely sensed and reanalysis products

ABSTRACT

Land surface all-wave net radiation (R_n) is crucial in determining Earth's climate by contributing to the surface radiation budget. This study evaluated seven satellite and three reanalysis long-term land surface R_n products under different spatial scales, spatial and temporal variations, and different conditions. The results showed that during 2000-2018, Global Land Surface Satellite Product (GLASS)-Moderate Resolution Imaging Spectroradiometer (MODIS) performed the best (RMSE=25.54 Wm^-2, bias=−1.26 Wm^-2), followed by ERA5 (the fifth-generation of European Centre for Medium-Range Weather Forecast Reanalysis) (RMSE=32.17 Wm^-2, bias=−4.88 Wm^-2) and GLASS-AVHRR (Advanced Very-High-Resolution Radiometer) (RMSE=33.10 Wm^-2, bias=4.03 Wm^-2). During 1983-2018, GLASS-AVHRR and ERA5 ranked top and performed similarly, with RMSE values of 31.70 and 33.08 Wm^-2 and biases of −4.56 and 3.48 Wm^-2, respectively. The averaged multi-annual mean R_n over the global land surface of satellite products was higher than that of reanalysis products by about 10∼30 Wm^-2. These products differed remarkably in long-term trends variations, particularly pre-2000, but no significant trends were observed. Discrepancies were more frequent in satellite data, while reanalysis products showed smoother variations. Large discrepancies were found in regions with high latitudes, reflectance, and elevation which could be attributed to input radiative components, meteorological variables (e.g., cloud properties, aerosol optical thickness), and applicability of the algorithms used. While further research is needed for detailed insights.

KEYWORDS:

• All-wave net radiation
• remote sensing
• reanalysis
• evaluation
• spatio-temporal variation
• product

1. Introduction

All-wave net surface radiation (R_n) at the land surface is the algebraic sum of the radiative energy budget of the underlying surface from the shortwave (0.3 to 3.0 μm) to the longwave (3.0 to 100.0 μm) portions of the electromagnetic spectrum (Bisht and Bras 2010; Jiang et al. 2016). Mathematically, R_n is the sum of the net shortwave and net longwave radiation and is expressed as: (1) $\begin{array}{rl}& R_n=R_{ns}+R_{nl}\\ & R_{ns}=R_{si}-R_{so}=\left(1-\alpha \right)R_{si}\\ & R_{nl}=R_{li}-R_{lo}\end{array}$(1) where $R_{ns}$ is the net shortwave radiation; $R_{nl}$ is the net longwave radiation; $R_{si}$ and $R_{so}$ are the downward and upward shortwave radiations; α is the broadband shortwave albedo; $R_{li}$ and $R_{lo}$ are the downward and upward longwave radiations. In this study, the downward direction was defined as positive, and the unit for all radiative components was Wm⁻².

R_n characterizes the land surface energy balance, and it also plays a critical role in ecological, agricultural, biogeochemical, and hydrological processes, such as evapotranspiration (ET) and some energy-consuming processes (i.e. photosynthesis) (Dickinson et al. 2006; Guo et al. 2020). Hence, R_n is an essential input for various ecological models, and weather forecast models (Al-Riahi, Al-Jumaily, and Kamies 2003; Liu and Yu 2006). High-quality long-term land surface R_n is of great significance for a variety of scientific studies and applications (Gao et al. 2013).

Due to unreliable data continuity, poor spatial representativeness, high maintenance costs, and uneven distribution of ground measurements, R_n is usually obtained from model reanalysis and satellite data in practical applications. Model reanalysis products are usually generated by assimilating various types of observations into an atmospheric model to obtain the best estimates of the state of the atmosphere and surface (Decker et al. 2012). Many reanalysis products provide radiative components to derive R_n, such as National Center for Atmospheric Research (NCAR) Reanalysis Ⅰ and its updated version namely NCAR II (Kanamitsu et al. 2002) from the National Center for Environmental Prediction (NCEP), Climate Forecast System Reanalysis (CFSR) (Decker et al. 2012; Saha et al. 2010); the forty-years Reanalysis (ERA40) (Uppala et al. 2005), the Interim Reanalysis (ERA-Interim) (Simmons et al. 2006), and their latest version ERA5 (Hersbach et al. 2019) from the European Centre for Medium-Range Weather Forecast (ECMWF), note that ERA5 has been extended to 1950s (Bell et al. 2021); Modern-Era Retrospective Analysis for Research and Applications version (MERRA) (Bosilovich et al. 2011) and its updated version MERRA2 (Gelaro et al. 2017a) from the National Aeronautics and Space Administration (NASA); the Japanese 25-year Reanalysis (JRA25) (Onogi et al. 2007) and its latest version the Japanese 55-year Reanalysis (JRA55) (Kobayashi et al. 2015b) from the Japan Meteorological Agency (JMA), and so on. Remotely sensed products are mainly generated using satellite observations, and can be roughly divided into two categories: generated from Radiative Transfer Model (RTM) and from parameterization or statistical models. Specifically, RTM-based products are produced from RTMs taking satellite observations and atmospheric and surface parameters derived from satellite data (e.g. cloud optical thickness) as inputs, such as the Clouds and Earths Radiant Energy System (CERES) Synoptic Radiation Fluxes and Clouds 1 degree series products (namely CERES_SYN1deg_Ed3A <CERES3A for short> and CERES_SYN1deg_Ed4A <CERES4A for short>) (Kato et al. 2015a), International Satellite Cloud Climatology Project (ISCCP) flux data (namely ISCCP-FD) (Zhang et al. 2004) and its high-resolution version ISCCP-FH (Young et al. 2018), and the Global Energy and Water Cycle Experiment Surface Radiation Budget (GEWEX-SRB) (Fu et al. 1997; Pinker and Laszlo 1992). In contrast, some remotely sensed products, most of which were released recently, can provide radiative components or R_n generated from parameterization models (i.e. TERRA (Verma et al. 2016)) or statistical models (i.e. FLUXCOM (Jung et al. 2019) and Global Land Surface Satellite <GLASS> (Jiang et al. 2016)). The statistical models were usually built by linking the radiative parameters with other ancillary information using various machine learning methods, and some of them can provide R_n as an integral based on the relationship between R_n and shortwave radiation (i.e. GLASS), or between R_n and top-of-atmosphere (TOA) satellite observations (Xu, Liang, and Jiang 2022). Table 1 lists the characteristics of these global R_n product.

Table 1. The existing global R_n products (not all).

Products	Spatial resolution	Temporal resolution	Period	References or link
Reanalysis data products
NCAR/ NCAR II	T62 (200 km)	6 hourly	1979-present	Kanamitsu et al. (2002)
CFSR	T382 (38 km)	6 hourly	1979-2010	Saha et al. (2010) and Decker et al. (2012)
MERRA	0.5° × 0.625°	1 hourly	1979-present	Bosilovich et al. (2011)
MERRA2	0.5° × 0.625°	1 hourly	1979-present	Gelaro et al. (2017a)
ERA40	2.5°/1.125°	Monthly/6 hourly	1957-2002	Allan (2004)
ERA-Interim	T255 (80 km)	3 hourly	1979-2019	Simmons et al. (2006)
ERA5	0.25° × 0.25°	1 hourly	1950-present	Hersbach et al. (2019)
JRA25	2.5°/1.125°	Monthly/6 hourly	1979-2004	Xu and Randall (1996)
JRA55	T319 (∼55 km)	3 hourly	1958-present	Kobayashi et al. (2015a)
Remotely sensed products
CERES SYN1deg_Ed3A	1° × 1°	3 hourly	2000.03-present	Wielicki et al. (1998)
CERES SYN1deg_Ed4A	1° × 1°	1 hourly	2000.03-present	Kato et al. (2015a)
GEWEX-SRB	1° × 1°	3 hourly	1983-2007	Pinker and Laszlo (1992), Fu et al. (1997)
ISCCP-FD	280 km	3 hourly	1983-2011	Zhang et al. (2004)
ISCCP-FH	1° × 1°	3 hourly	1983-2017	Young et al. (2018)
FLUXCOM_RS/FLUXCOM_RS + METEO	0.0833° × 0.0833° / 0.5° × 0.05°	8 daily/daily	2001-2010/2001-2015	Jung et al. (2019)
GLASS-MODIS daytime	0.05° × 0.05°	daytime	2000-2018	Jiang et al. (2016)
GLASS-MODIS	0.05° × 0.05°	daily	2000-2020	https://www.glass.umd.edu/
GLASS-AVHRR	0.05° × 0.05°	daily	1982-2018	https://glass-product.bnu.edu.cn/

Although the performance of different R_n products was evaluated and claimed satisfactory by their producers, their validation procedures, including validation samples, study regions, or assessment periods, differed considerably. Hence, several studies conducted inter-comparisons of different R_n products for objective evaluations. For example, Decker et al. (2012) evaluated R_n from five reanalysis products (NCAR, CFSR, ERA40, ERA-Interim, and MERRA) against the in situ measurements from 33 global flux tower sites. The results indicated that nearly all these products had high correlation with each other and their R_n variabilities were very similar to that of solar radiation. The MERRA and ERA-Interim products performed the best, with lower root-mean-square-error (RMSE) values ranging from 50 to 70 Wm⁻², but all five products tended to overestimate R_n, mostly owing to their overestimation of shortwave radiation. With the release of the updated version of reanalysis products, such as NCAR II, MERRA2, and ERA5, Urraca et al. (2018) concluded that the new version of reanalysis products provided better radiative estimates, with higher accuracy and finer spatio-temporal resolutions than their previous versions (NCAR, MERRA, and ERA-Interim). In addition, some studies have also conducted for remotely sensed R_n products evaluation. Jia et al. (2016) validated surface R_n from CERES3A against comprehensive measurements collected from hundreds of stations at daily and monthly scales. And then, Jia et al. (2018) conducted the comparison between CERES3A to three reanalysis R_n products, including JRA55, ERA-Interim, and MERRA2, also using these validation samples. The overall validation accuracy of CERES3A R_n was higher than that of other three reanalysis products, with the RMSEs of 33.56 and 25.57 Wm⁻² at daily and monthly resolutions, respectively, which agreed with other studies showing that remotely sensed products usually perform better than reanalysis products in radiative estimates of the land (Liang et al. 2010) and sea surface (Liang et al. 2022). However, Jia et al. (2018) also pointed out an inconsistency occurring in the long term variations in CERES3A R_n around 2007, possibly due to a different version of the model inputs. Recently, Seo et al. (2020) evaluated R_n from ERA5 and NCAR and two remotely sensed radiative products, namely CERES_EBAF (Energy Balanced and Filled) and GEWEX-SRB, in the Arctic region based on ground measurements from six sites, and indicated that ERA5 R_n had the highest accuracy over this region, while the satellite-based products had the tendency to be overestimated at high latitudes. However, there are still only few R_n products evaluation studies, while most studies focused on evaluating the radiative components from different products, such as solar radiation (Damiani et al. 2018; Du et al. 2022; Jiang et al. 2020; Yu et al. 2019; Zhang et al. 2016), longwave radiation (Qin et al. 2020; Sianturi and Sartika 2020; Tang et al. 2021), and albedo (Corbea-Pérez et al. 2021; Liu et al. 2013; Qin et al. 2011). The major conclusions made by these studies all concluded that there are large discrepancies among various products in each radiative component, confirming that the reliability of the R_n obtained by combining the radiative components of different products needs to be carefully assessed. Moreover, even in the case of these R_n evaluation studies, there were limitations in terms of the comprehensiveness of the evaluated products and validation samples, the study region, the evaluation periods, and so on. Particularly, several R_n products listed in Table 1, especially the recently released ones, including the new version reanalysis (i.e. MERRA2, and JRA55) and remotely sensed products (e.g. CERES4A, ISCCP-FH, FLUXCOM, GLASS), have not been thoroughly assessed. Therefore, the performance of global R_n products needs to be validated and evaluated more rigorously and comprehensively, so that their advantages can help to improve the R_n estimation accuracy in the future and provide a valuable reference for product selection in practical applications.

The main aim of this study was to objectively evaluate the performance of ten R_n products, including seven remotely sensed products (GLASS-MODIS, GLASS-AVHRR, CERES4A, ISCCP-FH, GEWEX-SRB, FLUXCOM_RS, and FLUXCOM_RS + METEO <referred to as FLUXCOM_ALL>), and three reanalysis products (ERA5, MERRA2, and JRA55), in accuracies and spatio-temporal variations at daily time scales over the land surface using the most comprehensive R_n ground measurements collected from nearly 600 globally distributed sites during 1991–2018. The organization of this paper is as follows: Section 2 gives a detailed information about the R_n ground measurements for validation and the ten evaluated R_n products. Section 3 describes the methods. Section 4 provides the evaluation results and their analysis. Finally, major conclusions and discussions are provided in Section 5.

2. Data and pre-processing

2.1. In situ R_n measurements

In this study, the ground radiative observations from 584 sites in 26 measurement networks around the world (Figure 1) from 1991 to 2018 were collected for products evaluation. The spatial distribution of the 584 sites was shown in Figure 1(a). All sites were globally distributed, from the Antarctic to the Arctic, but were mostly located in the northern hemisphere and ranged from −9 m to 5038.6 m above the sea level in elevation. These sites were characterized by various land-cover types and ecosystem conditions within various climatic zones. The land covers of these sites were classified into sixteen types defined by the International Geosphere-Biosphere Programme (IGBP). The proportions of these sites for all sixteen land cover types and five elevation ranges were presented in Figure 1(b,c), respectively. In addition, the homogeneousness of these sites has been examined by referring to Li et al. (2021) by using the MODIS Land Cover Climate Modeling Grid (MCD12C1) product (https://lpdaac.usgs.gov/products/mcd12c1v006). Based on the results, within the 25 × 25 km² regions surrounding each site, the land cover types for more than 79% of the sites were the same or similar. Hence, this demonstrates the spatial representativeness and the comprehensiveness of the measurements used for Rn products evaluation.

Figure 1. (a) Spatial distribution of 584 stations from 26 measuring networks. The proportions of all sites for (b) sixteen land-cover types defined by IGBP: and (c) five elevation ranges.

The detailed information about the 26 measuring networks/projects is listed in Table 2. Note that only the site at which both R_n and R_si observations provided was applied in this study as R_si was used for sky condition determination (see Equation (2)). As shown in Table 2, these radiative measurements were sampled at different frequencies at different sites and were recorded in various formats. Hence, at first, all quality-controlled R_n measurements or the ones obtained by adding radiative components measurements (see Equation (1)) were transferred into the local time at first, and then, the measurements whose sampling frequency was less than hourly were aggregated into hourly means with no missing data; thirdly, the hourly values were aggregated into daily means as long as at least 22 hourly values were avaliable in one day. Finally, all calculated daily values were checked manually again to remove any unreasonable values. The daily R_si samples were obtained at the same way. Therefore, there were 513,486 daily R_n and R_si samples during 1991–2018 were collected in this study.

Table 2. Details of the 26 observing networks.

Network Abbr.	No. of sites	Temporal resolution	Period	Instruments
LaThuile/Fluxnet	310	30 min	1991–2008	Kipp & Zonen CNR-1
AsiaFlux	31	30 min	1999–2008	Kipp & Zonen, CM-6F
ARM^a	22	1 min	2002–2013	Kipp & Zonen Pyrgeometer
BSRN^b	21	1 or 3 min	1992–2012	Eppley, PIR/Kipp & Zonen CG4
SURFRAD^c	7	3 min	1995–2012	Eppley, PIR
GAME.AAN	2	30 min	1997–2003	EKO MS0202F
BOREAS^d	5	30 min	1993–1996	Kipp & Zonen CM-5
GC-Net^e	11	1 min	1995–2012	Li Cor Photodiode & REBS Q* 7
CEOP^f	13	30 min	2007–2009	Eppley PIR, CG4
CEOPInt	8	30 min	2002–2009	QMN101
SMOSREX^g	1	30 min	2005–2010	Kipp & Zonen CNR-1
CERN^h	1	30 min	2007	/
FIFE^i	8	30 min	1987–1989	REBS net radiometer Model Q*3
LBA-ECO^j	6	1 h	1999–2006	REBS Q*7.1 CG2
Hi-Water	32	10 min	2012	CNR-4, CNR-1
HAPEX-Sahel	19	30 min	1992	Net radiometer
SAFARI.2000^k	2	30 min	2000	Kipp & Zonen Pyrgeometer
AmeriFlux	2	30 min	2000–2017	Kipp & Zonen, CNR-4
ASIAQ	2	30 min	2008–2011	Kipp & Zonen CNR-1
ChinaFlux	4	30 min	2003–2015	Kipp & Zonen CNR-1
EOL^l	10	30 min	2005–2009	/
EFDC^m	19	30 min	2000–2017	Kipp & Zonen CNR1-
IMAU^n	6	1 h	2000–2017	Kipp & Zonen, CNR-1
PROMICE^o	21	10 min	2000–2017	Kipp & Zonen CNR-1/CNR-2
FGI^p	9	10 min	2009–2019	/

^a ARM: Atmospheric Radiation Measurement (Phillips et al. 2017), ^bBSRN: Baseline Surface Radiation Network (Zo et al. 2017),^cSURFRAD: Surface Radiation Network(Augustine, DeLuisi, and Long 2000), ^dBOREAS: Boreal Ecosystem-Atmosphere Study, ^eGC-Net: Greenland climate network (Ni et al. 2020), ^fCEOP: Coordinated Enhanced Observing Period (Yao et al. 2014), ^gSMOSREX: Surface Monitoring Of Soil Reservoir Experiment (de Rosnay et al. 2006), ^hCERN: Chinese Ecosystem Research Network (Fu et al. 2010), ⁱFIFE: International Satellite Land Surface Climatology Project, ^jLBA-ECO: Large Scale Biosphere-Atmosphere Experiment (Tota et al. 2008), ^kSAFARI.2000: Southern African Regional Science Initiative Project 2000 (Swap et al. 2003),^lEOL: Earth Observing Laboratory, ^mEFDC: European Fluxes Database Cluster, ⁿIMAU: Institute for Marine and Atmospheric research Utrecht, ^oPROMICE: Programme for Monitoring of the Greenland Ice Sheet, ^pFGI: Finnish Meteorological Institute.

By considering the generation of some products, such as FLUXCOM and two GLASS R_n which were generated from the empirical models trained with numerous in-situ measurements, all validation samples were further carefully checked to guarantee their independence. However, for a more objective evaluation, the validation samples from the sites whose measurements have not been used for any R_n products production were screened and used for further validation. Figure 2 displays the spatial distribution of the sites providing only validation samples (pink circles). The figure shows that these sites are distributed unevenly, mainly in North America and Eastern China, and they account for less than 20% of all sites, providing 84,952 samples during 1991–2018. Therefore, two validation results, including the one using all validation samples and the other using only the spatially and temporally independent validation samples, were obtained in this study.

Figure 2. Spatial distribution of validation-only sites and sites for both training and validation.

In order to evaluate the R_n product under different sky conditions, the clearness index (CI) which is commonly used to represent atmospheric transmittance was used to classify the sky condition. CI is defined as the ratio of the daily $R_{si}$ to extraterrestrial radiation ($R_{se},\mathrm{W}{\mathrm{m}}^{-2}$) as Equations (2) and (3) (Li et al. 2021; Liu and Jordan 1960) shows, and the daily $R_{si}$ samples from sites were used for CI calculation. (2) $\begin{array}{c}CI={\displaystyle \frac{R_{si}}{R_{se}}}\end{array}$(2) where $R_{se}$ is calculated using Equation (3a) according to Irmak et al. (2003): (3a) $R_{se}={\displaystyle \frac{1440G_{sc}{d}_r}{\pi }\left(w_s\sin \left(\phi \right)\sin \left(\delta \right)+\cos \left(\delta \right)\sin \left(w_s\right)\right)}$(3a) (3b) $d_r=1+0.033\cos \left({\displaystyle \frac{2\pi doy}{365}}\right)$(3b) (3c) $\delta =0.409\sin \left({\displaystyle \frac{2\pi doy}{365}-1.39}\right)$(3c) (3d) $w_s=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\left(-\tan \left(\phi \right)\tan \delta \right)$(3d) (3e) $\mathrm{r}\mathrm{a}\mathrm{d}={\displaystyle \frac{\pi }{180}\ast \left(decimaldeg\right)}$(3e) where $G_{sc}$ is the solar constant (0.0820 MJm⁻² min⁻¹), $d_r$ is the inverse relative distance from the Earth to the Sun, $w_s$ is the sunset hour angle (rad), $\phi$ is the latitude (rad), δ is the solar declination (rad), and doy is the day of the year. According to a previous study (Iziomon, Mayer, and Matzarakis 2000), the sky was considered clear when CI was larger than 0.7, cloudy when CI was between 0.3 and 0.7, and overcast when CI was smaller than 0.3.

2.2. R_n datasets

Ten daily R_n products, including seven remotely sensed products, and three reanalysis products were evaluated in this study, and the characteristics of these products are listed in Table 1. For an objective evaluation, all R_n products should be pre-processed into the same format. Firstly, these products were converted into the local time as the in situ measurements. Then, they were aggregated into daily means. Products (e.g. GEWEX-SRB, ISCCP-FH, and JRA55) with temporal resolution larger than 1 h should be interpolated into an hourly scale at first. Finally, the daily R_n values of the ten products at each of the 584 validation sites were extracted according to the locations of these sites, one by one. Note that all products were extracted at their original spatial resolution for evaluation because the influence of their different spatial scales was small as long as the spatial scale was larger than 5 km according to Yan et al. (2016); additionally, the spatial representativeness of the validation samples was satisfactory.

Figure 3 shows that the ten evaluated products were of different duration. After the year 2000, remotely sensed R_n products increased as more satellite observations became available. Hence, two periods during which the evaluated products were the most and typical were considered for evaluation: Period I for years 1983–2018 (No. of samples = 513,486) and Period Ⅱ for years 2000–2018 (No. of samples = 470,203). The validation samples for some products were fewer because of their shorter durations (i.e. ISCCP-FH, GEWEX-SRB, FLUXCOM_RS, and FLUXCOM_ALL) or coarser temporal resolution (>1 day) (i.e. FLUXCOM_RS and FLUXCOM_ALL). A brief description of each R_n product follows below.

Figure 3. Gantt plot of the available time ranges of each of the ten products. Period Ⅰ (1983–2018) and Period Ⅱ (2000–2018) were the two evaluation periods considered in this study.

2.2.1. Remotely sensed products

1. GLASS-AVHRR and GLASS-MODIS

GLASS daily R_n products at 0.05° include GLASS-MODIS from 2000 to 2020 and GLASS-AVHRR from 1982 to 2018, named after the satellites Moderate-resolution Imaging Spectro-radiometer (MODIS) and Advanced Very-High-Resolution Radiometer (AVHRR) used for production, which explains the difference in the two datasets’ durations. The production flowchart is shown in Figure 4.

Figure 4. Workflow of GLASS-MODIS/GLASS-AVHRR daily R_n algorithm and production. Detailed information about the DSR-based Algorithm and LRD Algorithm are shown in the right panel.

Based on the idea for generating the GLASS daytime R_n product (Jiang et al. 2016; 2018), GLASS-MODIS and GLASS-AVHRR daily R_n were both estimated using the Downward Shortwave Radiation (DSR)-based algorithm, which was built based on the relationship between daily R_n and $R_{si}$ as well as its variations under various conditions by combining other ancillary information (e.g. meteorological factors, Normalized Difference Vegetation Index <NDVI>, broadband albedo <α>, etc.). However, in comparison to the GLASS daytime R_n estimation, it was found that the variations in the relationship between daily R_n and $R_{si}$ were also influenced by the length of daytime (dt) in hours besides α and NDVI, which were combined to characterize the land surface. Hence, dt was added to the daily R_n estimation in the DSR-based algorithm. As shown in Figure 4(b), the DSR-based algorithm consists of four conditional models (case1-case4) classified by the combinations of α, NDVI, and dt. Specifically, the case2-case4 models are applied when the land surface is determined as non-vegetated (NDVI < 0.2 & α > 0.7), while the case1 model is used for vegetated surfaces. Particularly, the case4 model was built for snow-covered or high-latitude areas with high albedo. All the four models were constructed using the Multivariate Adaptive Regression Splines (MARS) method (Jiang et al. 2016). However, the DSR-based algorithm is effective as long as $R_{si}$ is available or it dominates R_n (mainly at low-mid latitudes), while in the other cases (i.e. high-latitudes or short dt) the daily R_n in GLASS-MODIS and GLASS-AVHRR is generated in different ways. For GLASS-MODIS, the daily R_n was calculated by the LRD (length ratio of daytime) algorithm (Figure 4(c)) developed by Chen et al. (2020), which estimates daily R_n directly from MODIS TOA observations with a machine learning method. However, for GLASS-AVHRR, the daily R_n was filled with the linearly calibrated ERA5 daily R_n according to the one from GLASS-MODIS since the AVHRR TOA observations had too many missing data to apply the LRD algorithm. Moreover, for the overlapping regions where the daily R_n estimates derived from different algorithms or from ERA5, the final R_n values were weight-averaged with the weights determined by the in situ measurements within this region. The data used for GLASS-MODIS and GLASS-AVHRR daily R_n production are shown in Figure 4(a), including in situ R_n/$R_{si}$ measurements, GLASS high-level products including broadband albedo (Liu et al. 2013), NDVI (Xiao et al. 2017), and DSR (Zhang et al. 2019), meteorological factors (i.e. air temperature <Ta>, wind speed <W>, surface air pressure <PS>, etc.) from MERRA2, ERA5 daily R_n, and MODIS TOA bands (MOD/MYD02) and the corresponding geometry information (MOD/MYD03). Although the GLASS daily R_n products have not been evaluated thoroughly, the study conducted by Guo et al. (2020) indicated that the estimated ET was more accurate when using daily R_n from GLASS-MODIS than when using MERRA2.

According to Figure 3, GLASS-MODIS was only evaluated in period I, and GLASS-AVHRR was evaluated for both periods I and II.

• (2) CERES4A

The CERES SYN1deg series radiative products (Wielicki et al. 1996) were generated based on a radiative transfer model, with the atmospheric and surface properties obtained from a variety of sources used as the model inputs, but constrained by both solar-reflected and Earth-emitted radiation from the TOA for all-sky and clear-sky conditions (Gupta et al. 2010; Kato et al. 2015b; Rutan et al. 2015). Recently, CERES4A, which is the updated version of CERES3A by altering input data such as aerosol optical thickness, atmospheric profiles, and cloud properties, was released. Compared to CERES3A, the temporal resolution of CERES4A increased to hourly. While the inconsistency issue was reported to be well addressed (Doelling et al. 2013; 2016), the performance of CERES4A daily R_n has not been assessed comprehensively. The daily R_n of CERES4A obtained by adding four radiative components during period II (2000–2018) was validated in this study. Note that CERES4A provides data only since March, 2000.

• (3) ISCCP-FH

The ISCCP-FH climate data record (CDR) is the third-generation and high-resolution version of the ISCCP radiative flux profile product. Compared with its predecessor ISCCP-FD (Qin et al. 2011; Zhang 2004), ISCCP-FH contains a suite of level 2 and 3 products for monitoring the distribution and variation of cloud and surface properties to better understand the effects of clouds on climate and the radiation budget. The radiative estimates in ISCCP-FH were obtained every three hours on a 110-km equal-area (EQ) global grid covering the period July 1983 to June 2017 (Wang, Zhou, and Duan 2011; Young et al. 2018), and its R_n was also the sum of four radiative components. In this study, the daily R_n from ISCCP-FH were validated for periods 1983 to 2017 and from 2000 to 2017, which were shorter that period I and II.

(4)

GEWEX-SRB

The GEWEX-SRB project is a major component of GEWEX radiation research, which aims to offer global surface and atmospheric shortwave and longwave fluxes at 1° every 3 h to analyze climate variations and decadal trends from July 1983 to 2007. The GEWEX-SRB shortwave products were generated from the algorithm developed by Pinker and Laszlo (1992) based on radiative transfer theory, and the longwave radiation algorithm used the thermal radiative transfer code developed by Fu et al. (1997). Therefore, the GEWEX-SRB R_n is also the sum of four radiative components and was validated between 1983 and 2007.

(5)

FLUXCOM

The FLUXCOM initiative aims to improve understanding of the multiple sources and facets of uncertainties in empirical upscaling and, ultimately, to provide an integration of machine learning-based global flux products to the scientific community. It comprises two sets of products. FLUXCOM_RS was estimated exclusively from MODIS satellite data and does not require global climate forcing datasets as inputs. It also excludes potentially important information on meteorological conditions for biosphere-atmosphere fluxes and limits the temporal coverage to the MODIS era. For this reason, such datasets are themselves subject to uncertainty and are limited in spatial resolution. The second product is FLUCOM_RS + METEO (referred to FLUCOM_ALL), which was estimated from daily meteorological data and mean seasonal cycles of satellite data. It makes use of daily meteorological conditions and mean seasonal cycles of satellite-derived input drivers allowing for estimating fluxes beyond the satellite era. The FLUXCOM_RS gridded products at 0.0833° every 8 days were generated with nine machine learning methods, while FLUXCOM_ALL used three machine learning techniques with four global climate forcing data sets that yielded products with daily temporal and 0.5° spatial resolution (Jung et al. 2019). The R_n of FLUXCOM_RS and FLUXCOM_ALL is the sum of R_ns and R_nl. In this study, FLUXCOM_RS and FLUXCOM_ALL daily R_n were validated during 2001–2010 and 2001–2015, respectively.

2.2.2. Reanalysis product

1. ERA5

ERA5 is the 5th generation climate reanalysis dataset from ECMWF. The most substantial upgrades compared to the ERA-Interim are the finer spatial grid from 79 to 31 km, the higher time resolution from 3-hourly to hourly, the higher number of vertical levels from 60 to 137, a new numerical weather product model (IFS <integrated forecasting system> Cycle 41r2) and the increase of the amount of data assimilated. The data assimilation model is also a 12-hourly 4DVar (Hersbach et al. 2020; Martens et al. 2020; Seo et al. 2020). In addition, this dataset contains records from 1979 to the present. The daily R_n in ERA5 is also the sum of the radiative components and it was evaluated in both period I and II.

(2)

MERRA2

MERRA2 is a new modern satellite era atmospheric reanalysis from the NASA GMAO (Global Modeling and Assimilation Office) (Reichle et al. 2017). It replaced its predecessor MERRA with an upgraded version of the Goddard Earth Observing System model accompanied by advanced assimilation technology, microwave observations, and GPS-Radio Occultation datasets (Zhang et al. 2020). It also has reduced some water cycle deviations and imbalances and improved in aerosol data assimilation, cryosphere, and stratosphere (Gelaro et al. 2017b). The MERRA2 longwave radiative processes are described in Chou (1994) and the shortwave radiative processes are based on Suárez and Chou (2013). MERRA2 directly provides hourly net radiation at 0.5° × 0.625° since 1979. The MERRA2 daily R_n was evaluated in both period I and II.

(3)

JRA55

JRA55 is the second global atmospheric reanalysis product and the first comprehensive reanalysis developed by JMA covering the period from 1958 to the present. JRA55 calculated longwave radiation fluxes based on the band emissivity method with diffusivity approximation and shortwave radiation fluxes by a two-stream method using the Delta-Eddington approximation (Joseph, Wiscombe, and Weinman 1976). It is considerably better than its predecessor JRA25, as it is based on a new data assimilation and prediction system (DA) that improves many deficiencies found in JRA25 and added several additional observational data sources. JRA55 provides four radiative components and R_n is their sum. This product was evaluated in both the period I and II.

3. Methods

In this study, the performance of daily R_n from ten products was first evaluated at site scale, followed by an inter-comparison of their spatio-temporal variations. For clarification, the evaluation products for the two periods are listed in Table 3. Hence, there are six products evaluated in period I and ten products in period II. Note that the lengths of some products were shorter than period I or II (see Section 2.2).

Table 3. The R_n products evaluated in the two periods.

	Products for evaluation
Period I (1983–2018)	GLASS-AVHRR, ISCCP-FH*, GEWEX-SRB*, ERA5, JRA55, MERRA2
Period II (2000–2018)	GLASS-MODIS, GLASS-AVHRR, ISCCP-FH*, GEWEX-SRB*, CERES4A*, FLUXCOM_ALL*, FLUXCOM_RS*, ERA5, JRA55, MERRA2

*Indicates that the length of this product is shorter than the corresponding evaluation period. Detailed information is given in Section 2.2.

3.1. Validation accuracy calculation

Validation accuracy was assessed with three common statistical measures: Bias, RMSE, and R², which depict the sample’s overall bias, errors, and overall goodness of fit with a 1:1 line, respectively, calculated as: (4) $Bias={\displaystyle \frac{{\sum }_{i=1}^n\left(E_i-O_i\right)}{n}}$(4) (5) $RMSE=\sqrt{{\displaystyle \frac{{{\sum }_{i=1}^n\left(E_i-O_i\right)}^2}{n}}}$(5) (6) $R^2=1-{\displaystyle \frac{{{\sum }_{i=1}^n\left(E_i-\overline{O}\right)}^2}{{{\sum }_{i=1}^n\left(O_i-\overline{O}\right)}^2}}$(6) where $E_i$ is the R_n estimate from product, $O_i$ is the corresponding R_n measurement, and $n$ is the number of samples.

3.2. Spatio-temporal variations calculation

In this study, the spatio-temporal variations were represented by the global mean and the long-term trend.

1. Annual mean R_n over global land surface

The annual mean R_n over the global land surface for each year in the evaluation period of each product was calculated. By considering the change in grid area at different latitudes, the latitude weighted average method was applied for each product according to Equation (7): (7) $\begin{array}{c}\overline{R_{ny}}={\displaystyle \frac{\sum \left(R_{ny,k}\cos \left({\theta }_k\right)\right)}{\sum \cos \left({\theta }_k\right)}}\end{array}$(7) where $\overline{R_{ny}}$ is the averaged annual mean R_n over global land surface, $R_{ny,k}$ is the annual mean R_n of the kth pixel by averaging all available daily R_n values in this year for this pixel, ${\theta }_k$ is the corresponding latitude of the kth pixel, and cos (${\theta }_k$) is the weight of this pixel.

(2)	Long-term trend

To avoid the inclusion of parametric deviations (Fensholt and Proud 2012; Liu et al. 2015), the widely used nonparametric Mann-Kendall (M–K) test (Kendall 1948) was employed for calculating the trend in this study. M–K is a rank-based method that identifies trend of a time-series dataset, which does not require the samples to conform to a certain type of normal distribution as the linear regression method, and is not easily affected by missing values and outliers (Gao et al. 2013). Its calculation is as follows: (8a) $\begin{array}{c}Z=\left\{\begin{array}{cc}{\displaystyle \frac{S-1}{\sqrt{Var\left(S\right)}},}& \mathrm{i}\mathrm{f}S>0\\ 0,& \mathrm{i}\mathrm{f}S=0\\ {\displaystyle \frac{S+1}{\sqrt{Var\left(S\right)}},}& \mathrm{i}\mathrm{f}S<0\end{array}\right\}\end{array}$(8a) (8b) $S=\sum _{i=1}^{n-1}\sum _{j=i+1}^{n}sgn\left(R{n}_i-R{n}_j\right)$(8b)

(8c) $\begin{array}{c}sgn\left(R{n}_i-R{n}_j\right)=\left\{\begin{array}{cc}1,& \mathrm{i}\mathrm{f}R{n}_j-R{n}_i>0\\ 0,& \mathrm{i}\mathrm{f}R{n}_j-R{n}_i=0\\ -1,& \mathrm{i}\mathrm{f}R{n}_j-R{n}_i<0\end{array}\right\}\end{array}$(8c) (8d) $\begin{array}{c}Var\left(S\right)={\displaystyle \frac{n\left(n-1\right)\left(2n+5\right)}{18}}\end{array}$(8d) where $Z$ is the standardized normal statistic with a positive, negative and zero value indicating an increasing, decreasing, and no trend, respectively; $S$ is the statistic of the test; $sgn\left(\right)$ is a symbolic function; $R{n}_j$ and $R{n}_i$ are the $t_j$ and $t_i$ ($t_j>t_i$) values of the data, respectively; $n$ is the number of samples. In this study, $\left|Z\right|$ ≥ 1.96 indicates a significant increasing (decreasing) trend at the 95% significance level. No significant trend is present if $\left|Z\right|$ < 1.96.

4. Results and analysis

4.1. Evaluation of validation accuracy

The performance of the ten R_n products at site scales was evaluated in terms of the validation accuracy and its robustness under various conditions in terms of land cover, elevation, and sky condition.

4.1.1. Site scale validation accuracy

The ten daily R_n products were validated against the corresponding in situ measurements during Period Ⅰ and Period Ⅱ, respectively. Table 4 presents the overall validation accuracy calculated by using all validation samples and the other one by using spatially and temporally independent validation samples (in italic font) for the two periods, respectively.

Table 4. The validation accuracy of each product in periods I and II.

Product	Time range	R^2	RMSE (Wm^−2)	Bias (Wm^−2)	No. of samples
Period Ⅰ (1983–2018)
GLASS-AVHRR	1983–2018	0.79 0.79	33.08 33.15	3.48 6.33	513,486 84,952
ISCCP-FH	1983–2017	0.71 0.69	46.69 43.77	−16.13 −10.12	440,520 67,435
GEWEX-SRB	1983–2007	0.72 0.64	37.42 42.50	7.84 16.00	132,099 37,209
MERRA2	1983–2018	0.74 0.75	37.40 36.41	1.69 4.01	513,486 84,952
ERA5	1983–2018	0.80 0.82	31.70 29.26	−4.56 −1.54	513,486 84,952
JRA55	1983–2018	0.72 0.68	42.35 43.49	−7.61 −0.27	513,486 84,952
Period Ⅱ (2000–2018)
GLASS-MODIS	2000–2018	0.88 0.85	25.54 27.18	−1.26 −1.94	470,203 67,832
GLASS-AVHRR	2000–2018	0.80 0.78	33.10 34.02	4.03 3.40	470,203 67,832
ISCCP-FH	2000–2017	0.68 0.68	48.84 45.75	−17.21 −12.57	451,354 51,911
GEWEX-SRB	2000–2007	0.75 0.61	35.34 41.88	4.70 8.06	88,193 6,371
CERES4A	2000.03–2018	0.76 0.75	35.49 36.59	3.50 5.05	470,203 67,832
FLUXCOM_ALL	2001–2015	0.78 0.75	31.41 33.73	2.62 6.75	193,388 15,791
FLUXCOM_RS	2001–2010	0.71 0.69	35.10 36.73	0.11 3.53	37,224 5,106
ERA5	2000–2018	0.80 0.82	32.17 30.73	−4.88 −4.71	470,203 67,832
MERRA2	2000–2018	0.76 0.75	37.61 36.92	1.04 0.01	470,203 67,832
JRA55	2000–2018	0.74 0.66	42.29 46.40	−8.24 −1.33	470,203 67,832

Notes: For each product, the validation results in the first row were calculated using all validation samples (overall validation accuracy), whereas the ones in the second row in italic font were calculated using the samples from the sites providing validation sample only.

From the overall validation accuracy calculated by using all validation samples, daily R_n of the remotely sensed products was generally superior to the reanalysis data products for both periods I and II, except in the case of ISCCP-FH, with the validated RMSE ranging from 25.54 to 37.42 Wm⁻², the bias ranging from −1.26 to 7.84 Wm⁻², and R² ranging from 0.71 to 0.88. In terms of the worse accuracy of ISCCP-FH, it is assumed that might be caused by the significant errors in its longwave radiation (Qin et al. 2011). Most remotely sensed products slightly overestimate daily R_n which is consistent with the study of Bojanowski, Vrieling, and Skidmore (2014), while the reanalysis products have the tendency to underestimate daily R_n with negative biases, except for MERRA2. Besides, there was little difference in the validation accuracies between Period Ⅰ and Period Ⅱ (i.e. GLASS-AVHRR, ISCCP-FH, ERA5, MERRA2, and JRA55), which indicated the comprehensiveness of the validation samples and the reliability of the validation results. Specifically, for Period Ⅱ, the validation accuracy in daily R_n of GLASS-MODIS was the best among all products evaluated with the validated RMSE of 25.54 Wm⁻², bias of −1.26 Wm⁻² and R² of 0.88, followed by ERA5 (RMSE = 32.17 Wm⁻², bias = −4.88 Wm⁻², R² = 0.80) and GLASS-AVHRR (RMSE = 33.10 Wm⁻², bias = 4.03 Wm⁻², R² = 0.80) with their RMSE values close to each other and all less than 34 Wm⁻². CERES4A, whose previous version CERES3A was acknowledged as one of the best R_n products (Jia et al. 2018), was at the fifth place (R² = 0.76, RMSE = 35.49 Wm⁻², bias = 3.50 Wm⁻²). Surprisingly, the daily ISCCP-FH R_n has poor accuracy with the largest RMSE value (48.44 Wm⁻²) and the largest bias (−17.21 Wm⁻²), and it even worse than that of MERRA2 (RMSE = 37.61 Wm⁻², bias = 1.04 Wm⁻²) and JRA55 (RMSE = 42.29 Wm⁻², bias = −8.24 Wm⁻²). In terms of the other two new remotely sensed products FLUXCOM_ALL and FLUXCOM_RS, their accuracy ranked ahead, with their validated RMSE values of 31.41 and 35.34 Wm⁻² and biases of 2.62 and 0.11 Wm⁻², which somewhat differed from the conclusions made by Jung et al. (2019), which used a very different validation samples size. For inter-comparison, GLASS-MODIS and GLASS-AVHRR R_n were also validated using the same validation samples as FLUXCOM_ALL, and their accuracies were still better with the RMSEs of 23.48 and 28.78 Wm⁻² and biases of 1.97 and 7.99 Wm⁻², respectively. For Period Ⅰ, ERA5 daily R_n ranked the best and GLASS-AVHRR performed similarly, followed by GEWEX-SRB, MERRA2, and JRA55, while ISCCP-FH daily R_n was still ranked last with the largest RMSE (46.69 Wm⁻²) and bias (−16.13 Wm⁻²). Regarding the validation from sites providing validation samples only, the rank of the evaluation results remained nearly unchanged. Although the accuracy of most products decreased slightly, only a few reanalysis data products improved somewhat. In particular, the accuracy for ERA5 daily R_n improved for both period I and II by reducing the RMSEs by 1.56 and 1.44 Wm⁻², respectively. Compared with the results in Figure 2, we speculate that the reason is the lack of these validation samples at high latitudes, where ERA5 daily R_n performed relatively poorly. Overall, the two validation results indicated that the validation results were reliable. In summary, GLASS-MODIS and GLASS-AVHRR daily R_n were the best performing satellite products and GLASS-MODIS was even better; ERA5 performed the best among all three reanalysis products in accordance with Wang and Clow (2021), and comparable to or even better than GLASS-AHVRR.

The spatial distribution of the overall validation accuracy at sites of the ten products for Period Ⅱ in terms of RMSE is shown in Figure 5. The variations in the validated bias values with latitude (30° S to 90° N) for each product are also presented, and were calculated by averaging all available bias values at the sites per 10° latitude.

Figure 5. Spatial distribution of site overall validation accuracy (represented by RMSE) for period Ⅱ of the ten evaluated products: (a) GLASS-MODIS, (b) GLASS-AVHRR, (c) ERA5, (d) CERES4A, (e) MERRA2, (f) JRA55, (g) ISCCP-FH, (h) GEWEX-SRB, (i) FLUXCOM_RS and (j) FLUXCOM_ALL. The left graph in each plot shows the variations in the validated bias with latitude, which was calculated by averaging all bias values per 10° latitude from 30°S to 90°N. The number of sites shown in (i) FLUXCOM_RS and (j) FLUXCOM_ALL was less due to missing data.

Figure 5 depicts that the ten daily R_n products showed worse accuracy in tropical areas (30°S to 30°N) than other regions. They almost all tended to be overestimated (positive bias) at low-middle latitudes and underestimated (negative bias) at high latitudes, except for the two FLUXCOM products (Figure 5(i–j)), which estimated R_n simply at high latitudes (Jung et al. 2019). Comparatively, the sites with the lower RMSE values (<25 Wm⁻², in blue) were the most for GLASS-MODIS (Figure 5(a)) and mainly distributed in eastern Asia, eastern American, European, and even Greenland. The bias magnitudes were below 5 Wm⁻² and decreased with the increases of latitude; whereas GLASS-AVHRR (Figure 5(b)) showed a similar pattern as GLASS-MODIS but with more sites with larger RMSE values (>40 Wm⁻², in warm red). The spatial pattern of ERA5 R_n (Figure 5(c)) was much better than other reanalysis products and was similar to GLASS-AVHRR, but its accuracy was worse than that of GLASS-AVHRR at high latitudes, and this demonstrated the necessity to calibrate ERA5 R_n referring to GLASS-MODIS at high-latitudes for GLASS-AVHRR. For other products, the sites with RMSE values larger than 40 Wm⁻² and magnitude of bias larger than 25 Wm⁻² were the most in ISCCP-FH (Figure 5(e)), the maximum magnitude of bias shown in GEWEX-SRB (Figure 5(h)) at around 32° N due to the overestimation of albedo (Qin et al. 2011), two FLUXCOM products performed comparatively with the GLASS products, although FLUXCOM_ALL was a little better, whereas the results of any other products were consistent with previous results. However, the evaluation of these products in the Southern hemisphere was insufficient due to too few ground measurements available.

4.1.2. Product performance under various conditions

The accuracy robustness of each product under different conditions varies due to differences in the input data for different products, as well as the applicability of the algorithms used (Liang et al. 2010). Hence, the validation accuracy of the ten R_n products (in terms of RMSE) for each land-cover type, elevation range, clear-/cloudy-sky conditions, and month was further examined and compared for Period Ⅱ (Figure 6(a–d)). The RMSE was used here because it is usually given more weight in accuracy evaluations, although a small RMSE does not always mean a small uncertainty when the average value of the parameter is small, such as in the case of R_n for land surfaces with high albedo (i.e. ICE) or at high-latitudes. Note that all validation samples were applied in this analysis.

Figure 6. Comparison of the validation accuracy in daily R_n (represented by RMSE) of the ten evaluated products under various conditions during Period Ⅱ: (a) land-cover types, (b) elevation ranges (200 m, 200–500 m, 500–1000 m, 1000–1500 m, and >1500 m), (c) sky conditions defined by CI (clear: CI > 0.7, mixed cloudy: 0.3≤CI≤0.7, overcast: CI < 0.3), and (d) twelve months. The accuracy for ICE and TUN of FLUXCOM_ALL and FLUXCOM_RS products was not available in (a) due to missing data (shown in oblique line).

Figure 6(a) shows the RMSE values of the ten products in daily R_n for sixteen land-cover types (see Section 2.1). The accuracy of all ten R_n products varied with different land-cover types, and almost all products estimated daily R_n poorly for forest (including ENF, EBF, DBF, and MF), WET, URB, and ICE, which was consistent with previous studies (Jiang et al. 2019). GLASS-MODIS and ERA5 yielded a similar and much better performance across all land-cover types than other products, with their RMSE mostly around 25 Wm⁻², followed by GLASS-AVHRR. Regarding the two FLUXCOM products, their performance was relatively robust for most land-cover types especially FLUXCOM_ALL, but their results were not comprehensive enough as much fewer samples were available for ICE and TUN. However, the fewer available ground measurements for some land-cover types, such as DNF, URB, and WET, influence the validation results. The accuracy comparisons for five elevation ranges is shown in Figure 6(b). Overall, all products performed worse at high elevation (>1500 m) regions. GLASS-MODIS R_n showed a robust and consistent performance for elevations less than 1500 m, with RMSE less than 25 Wm⁻², followed by FLUXCOM_ALL. GLASS-AVHRR and ERA5 R_n were the third-best performing products, with the RMSE of GLASS-AVHRR R_n increasing with elevation and vice versa in the case of ERA5. Though the performance of CERES4A R_n was also robust for all evaluation ranges, its RMSE values were mostly around 35 Wm⁻². Figure 6(c) shows the performance of these products under clear-sky (CI > 0.7), mixed cloudy (0.3 ≤ CI ≤ 0.7), and overcast (CI < 0.3). GLASS-MODIS R_n was also the best either for clear or cloudy sky, and the uncertainty in R_n from GLASS-AVHRR increased as cloudiness decreased; in contrast, that of ERA5 was different with the largest uncertainty in R_n appearing in the presence of mixed cloudy sky. GEWEX-SRB and ISCCP-FH show poor accuracy in cloudy sky condition as they mistake clouds for snow (Qin et al. 2011). Figure 6(d) shows obvious seasonal variations in the R_n estimates accuracy of all products, with the smaller RMSE values appearing in winter months and the larger ones appearing in summer months, which was possibly due to the uncertainty in radiative estimates, such as R_si (Zhang et al. 2016), caused by the difference in cloud coverage and water vapor content of the two seasons Zhang et al. (2020). Relatively, the GLASS and FLUXCOM series products, as well as ERA5 show a stable performance, with their RMSEs ranging from 15 to 30 Wm⁻² across twelve months. And, the R_n from ISCCP-FH performed similarly but was still least accurate under the four conditions discussed, and this was consistent with the results described above.

Moreover, we also found that GLASS-MODIS R_n performed particularly better when the land surface was covered by ice or snow, likely because of the application of the case4 model of the DSR-based algorithm and LRD algorithm, which were designed specifically for this case. For better illustration, the temporal variations in daily R_n from seven products and ground measurements at three sites with ICE land-cover (PM-SCO_L <72.223°N, 26.8182°W>, PM-EGP <75.62°N, 35.97°W>, and PM-QAS_M <61.09°N, 46.83°W>) during the randomly selected period were shown in Figure 7. The R_n of GEWEX-SRB and two FLUXCOM series products were not included in Figure 7 because of unavailable data. Figure 7 shows that all products could capture the general variations in R_n at these ICE sites but most with the tendency to underestimate R_n. Relatively, the R_n values from GLASS-MODIS (red line) and GLASS-AVHRR (orange line) were the closest to the in situ measurements at the three sites, especially during summer days, but they varied differently due to different algorithms under ice cover land at high latitudes. CERES4A R_n was the only one to be overestimated often, such as R_n at PM-EGP (Figure 7(b)), particularly during the doy of 170–220. Surprisingly, R_n estimates from ERA5 were heavily underestimated though its site scale validation accuracy was superiority, and the underestimation could be more than 100 Wm⁻² (Figure 7(a)), and the results were in accordance with the results in Table 4.

Figure 7. Temporal variations in daily R_n from GLASS-MODIS, GLASS-AVHRR, ISCCP-FH, CERES4A, ERA5, MERRA2, JRA55 and, ground measurements (black dots) at three validation sites of ICE land-cover: (a) PM-SCO_L, (b) PM-EGP, and (c)PM-QAS_M.

In general, GLASS-MODIS was the best performing among the ten evaluated products, and ERA5 R_n was the best one among the three new-generation reanalysis products and performed similarly to the GLASS series products. GLASS-MODIS and GLASS-AVHRR R_n products were satisfactory in validation accuracy and robustness. Besides, the uncertainty in R_n estimates was generally large for nearly all products in some regions, such as high-latitudes, tropical areas, high elevation (>1500 m), with a high albedo, even coastal regions, and so on. Hence, more efforts should be made in these regions to improve R_n estimation accuracy.

4.2. Inter-comparison of spatio-temporal variations in R_n

According to the results described above, R_n from six products, including three remotely sensed products (GLASS-MODIS, GLASS-AVHRR, and CERES4A) and three reanalysis products (ERA5, MERRA2, and JRA55), were further inter-compared in their spatio-temporal variations during Period II (2000–2018) based on their relative better validation accuracies, global continuous coverage, and long durations.

4.2.1. Spatial distribution in multi-annual mean R_n

At first, the global spatial distributions of the multi-annual mean R_n during Period II for each of the six products were calculated separately and presented in the diagonal of Figure 8. Generally, the spatial patterns of the multi-annual mean R_n estimates from the six products were similar and reasonable, with their values decreasing from tropic to the polar area. The highest values in the multi-annual mean R_n (>100 Wm⁻²) appeared in central and southern Africa, northern south America, and southeast Asia regions, whereas the lowest values showed in Antarctica, Arctic (e.g. Greenland), and high elevation areas (e.g. Mount Qomolangma).

Figure 8. Comparison of the global spatial distributions of multi-annual mean R_n between six products (GLASS-MODIS, GLASS-AVHRR, CERES4A, ERA5, MERRA2, and JRA55) for Period Ⅱ.

In Figure 8, the maps below the diagonal present the difference in the multi-annual mean R_n between any two products; they were computed by using the multi-annual mean R_n of the product in every row on the diagonal minus each of the other ones. For example, the plot at the right bottom corner is the difference in multi-annual mean R_n by GLASS-MODIS minus JRA55. Overall, there is a large discrepancy between these six products. Specifically, the GLASS-MODIS multi-annual mean R_n is very close to that of GLASS-AVHRR with only a little bit larger over northern Africa (20°–30° N), where they are both extremely larger than that of all other products with the difference even larger than 40 Wm⁻² (the first two columns from right to left), and this is possibly caused by the overestimated R_si as pointed by Zhang et al. (2019). For other products in multi-annual mean R_n, CERES4A was larger than others especially over the Qinghai-Tibet Plateau and Greenland which can be contributed to overestimation of downward longwave radiation according to Wang, Wang, and Xue (2021), an evaluation in three poles regions, however, it suggested that CERES4A has the highest accuracy in Qinghai-Tibet Plateau; ERA5 was similar to the three remotely sensed products within a smaller magnitude; MERRA2 showed a large overestimation in the south and southwest China, possibly due to the uncertainty in cloud fraction estimates (Du et al. 2022). JRA55 was lower than others, especially in high latitude such as Antarctica region and Greenland, whereas it showed an overestimation over the tropics which may due to unrealistic trends in analyzed tropical cyclone strength (Kobayashi et al. 2015b).

4.2.2. Spatio-temporal trend in annual mean R_n

The averaged multi-annual mean R_n over the global land surface ($\overline{R_{n_{multi}}}$), which is used in numerous studies, was calculated for each of the evaluated eight products (the two FLUXCOM products with missing data were excluded) and the results were list in Table 5.

Table 5. Multi-annual mean R_n over the global land surface ($\overline{R_{n_{multi}}}$) of eight products for Periods I and II.

	Product	$\overline{R_{n_{multi}}}$ (Unit: Wm^−2)
		Period Ⅱ	Period Ⅰ
Remotely Sensed Products	GLASS-MODIS	92.729	/
	GLASS-AVHRR	95.527	95.861
	CERES4A	98.162	/
	GEWEX-SRB	87.232 (2000–2007)	90.155 (1984–2007)
	ISCCP-FH	65.619 (2000–2017)	66.261 (1984–2017)
Reanalysis Products	ERA5	76.267	75.816
	MERRA2	86.277	86.398
	JRA55	68.098	66.918

Theoretically, the $\overline{R_{n_{multi}}}$ value should be robust regardless of the length of the calculated years as long as the time period is long enough. The calculated $\overline{R_{n_{multi}}}$ values for the products in both period I and II were very similar (i.e. GLASS-AVHRR, GEWEX-SRB, ISCCP-FH, ERA5, MERRA2, and JRA55), hence, demonstrating the reliability of the results. At a global perspective, the $\overline{R_{n_{multi}}}$ of nearly all remotely sensed products (except ISCCP-FH) was larger than that of three reanalysis products, with values ranging from ∼87 to ∼98 Wm⁻² and from ∼68 to ∼86 Wm⁻² for remotely sensed products and reanalysis products, respectively. The overall underestimation of ISCCP-FH daily R_n may due to the overestimation of surface albedo in caused by mistakenly treating the cloud as the snow when the cloud detection algorithms are conducted (Qin et al. 2011). Taking period II as an example, $\overline{R_{n_{multi}}}$ of CERES4A was the largest, with 98.162 Wm⁻², followed by that of GLASS-AVHRR and GLASS-MODIS with 95.527 and 92.729 Wm⁻², respectively; their values were similar and larger than those of ERA5 by about 20 Wm⁻², though their site validation accuracies were very close to each other. But the value of ERA5 (76.467 Wm⁻²) was close to the result from Wild et al. (2014), which calculated $\overline{R_{n_{multi}}}$ from CMIP5 simulations equal to 72.5 Wm⁻². Hence, combining the validation results leads to speculate that the R_n generated from models was possibly an underestimation.

The climatology of the anomaly in $\overline{R_{n_y}}$ for each of the eight products during their own periods from 1983 to 2018 was inter-compared in Figure 9.

Figure 9. Inter-annual variations in the anomalies of the annual mean R_n averaged over the global land surface for eight products (as Table 5) during their own duration periods from 1983 to 2018. The gray area indicates Period II (2000–2018).

Overall, the discrepancy in the variations in the anomalies of $\overline{R_{n_y}}$ of all eight products was large, especially before the year 2000. For all remotely sensed products, the fluctuations in the $\overline{R_{n_y}}$ from GEWEX-SRB and ISCCP-FH were dramatic and different from others, which can be attributed to errors in the cloud detection algorithm, which may reflect inadequacies in simulating the effects of volcanic aerosols on longwave radiative transfer or unrealistic cloud structure and stratospheric thermal layers, as they both used the cloud properties derived from the ISCCP-DX data (Qin et al. 2011). Hence, GEWEX-SRB and ISCCP-FH present homologous systematic errors that could explain the similarity in the fluctuating trends shown in Figure 9 (Zhang et al. 2015). Besides, Qin et al. (2011) implied that neglecting the elevation effect influenced the shortwave of GEWEX-SRB as well. While the $\overline{R_{n_y}}$ of other three products (GLASS-MODIS, GLASS-AVHRR, and CERES4A) were relatively stable after 2000 but still in different trends, even in the case of the two GLASS products; whereas a sudden drop around 2000 was observed in $\overline{R_{n_y}}$ for all four remotely sensed products (GLASS-AVHRR, GLASS-MODIS, GEWEX-SRB, and ISCCP-FH) that contained data of this year, which was possibly caused by the satellite data inputs. A sudden jump in 2015 and a sharp decrease since 2017 were observed in $\overline{R_{n_y}}$ of CERES4A, which differed from other products and seemed unreasonable. The overestimated aerosol optical depths (AODs) in 2015 of CERES4A used for clear-sky irradiance computations (Fillmore et al. 2022), the significant errors in the upward longwave radiation estimates pointed out by Liu et al. (2020), and a shift in the retrieved aerosol optical thickness over land caused by the switch from MODIS Collection 5 to Collection 6 for the inputs in 2017 might be the possible reasons. The three reanalysis products varied smoothly during the whole period but were still inconsistent in trend. A sudden drop in $\overline{R_{n_y}}$ during 1991–1993, which was influenced by the increase in aerosols caused by the eruption of Mount Pinatubo in 1991 (Frame and Stone 2013), was only captured by GLASS-AVHRR, ERA5, and MERRA2. What is striking is that JRA55 could not capture the significant plunge due to eruption, although Kobayashi et al. (2015b) has stressed that numerous meteorological are not necessarily of sufficient quality for evaluation purposes, the reason behind this merits further investigation.

Moreover, the consistency in the spatio-temporal variations in annual mean R_n of each product calculated with the M–K method was also examined. By considering the reasonability and the number of products, four products during Period Ⅰ (GLASS-AVHRR, ERA5, MERRA2, and JRA55), and six products during Period II, including three remotely sensed products (GLASS-AVHRR, GLASS-MODIS, and CERES4A) and three reanalysis products (ERA5, MERRA2, and JRA55), were analyzed respectively. The spatial sizes of the six products were all processed to 5 km before the calculation of the trend. Then, the number of products with the same trends at each pixel was counted and represented by different colors and shown in Figure 10.

Figure 10. Statistical results on the spatio-temporal trends in annual mean R_n from (a) four products (GLASS-AVHRR, ERA5, MERRA2, and JRA55) during Period Ⅰ, and (b) six products (GLASS-AVHRR, GLASS-MODIS, CERES4A, ERA5, MERRA2, and JRA55) during Period II. The triangle legend represents the consistence among these products in different trends: blue for significant increase (‘I’), green for significant decrease (‘D’), and pink for non-significant trend (‘N’); the tone of the three colors indicates the number of products having the same trend; the darker the color, the greater the number of products.

From the results in Figure 10, the hotspots with a significant trend in annual mean R_n that agreed across most products were mainly located in tropical and high-latitude areas either in Period I or II. Specifically, for both period I and II, a significant increase in annual mean R_n in northwestern Africa, southeastern Australia, and some sparsely distributed regions at high latitudes, and a significant decrease over northern South America, particularly over the Amazon rainforest region, southern Africa, central oceanic, most Southeast Asia, and most areas of Antarctic close to the ocean were observed, which indicates that these hotspots which are very sensitive to radiation balance should be paid more attention and the consistence trend in R_n over these hotspots is in a high probability. As a matter of fact, some of these hotspots have also been discussed in previous studies. For instance, it was speculated that the reduced available R_n in the Amazon rainforest region was caused partly by the altered albedo and radiation partitioning due to deforestation (Oliveira et al. 2019; Restrepo-Coupe et al. 2013). However, more hotspots need to be further studied. In addition, the discrepancy in the global R_n trend between the two periods was more significant, especially at mid-high latitudes in the Northern hemisphere (e.g. east southern China, central Europe, south Asia, and so on), where the significant trend in R_n during period I seems mainly occurring before 2000. However, by considering the dramatically drop around 1991 in R_n, the large uncertainty in most products before 2000, an inconsistent variations in R_n around 2000, as well as only one remotely sensed product used in calculating the results in Figure 10(a). Consequently, those regions with a different trend in R_n between the two periods should be analyzed more cautiously.

Therefore, the results point out that the discrepancy in the variations in R_n among different products is large, and all the questions raised above require further exploration and investigation in future research.

5. Conclusions and discussions

As one of the most critical parameters in land-surface interactions, the surface radiation budget R_n is usually obtained from various products in practical use due to many limitations in ground measurements. Hence, the evaluation of R_n products is of great significance. In this study, ten R_n products, including seven remotely sensed products (GLASS-MODIS, GLASS-AVHRR, CERES4A, ISCCP-FH, GEWEX-SRB, FLUXCOM_RS, and FLUXCOM_ALL), and three reanalysis products (ERA5, MERRA2, and JRA55), have been compared and evaluated against the most comprehensive in situ measurements collected from 584 globally distributed sites during 1991 to 2018. The analysis found that large discrepancies exist in various R_n products in terms of accuracy, robustness, consistency, and spatio-temporal variations. Several major conclusions have been drawn as follows:

1. The performances of these ten R_n products varied differently from each other, but nearly all products performed poor in some special regions, particularly at high-latitudes, within high elevation zones, and with high reflectance. Hence, more efforts should be made in these regions to improve the accuracy of R_n estimation. Moreover, better site scale validation precision does not mean reasonable spatio-temporal trends (e.g. ERA5, GLASS-MODIS, and GLASS-AVHRR).

2. Remotely sensed products match the in situ land surface R_n much better than those of most reanalysis products, and they are also better at characterizing the spatio-temporal variations in R_n. However, inconsistencies appear in most remotely sensed products, mainly due to their satellite data inputs, and some of them show the dramatic and unreasonable trends (i.e. ISCCP-FH and GEWEX-SRB). Regarding the reanalysis products, their R_n estimates tend to be underestimated, and their temporally variations are too smooth though their time series are usually longer and more consistent.

3. Among these products, GLASS-MODIS and GLASS-AVHRR R_n performed satisfactorily, showing strong potential to be widely used because of their finer spatial resolution, long time series, and global continuous coverage without missing. Particularly GLASS-MODIS, which had the best validation accuracy yielding R² of 0.88, RMSE of 25.54 Wm⁻², and bias of −1.26 Wm⁻² during the Period II, and performed robust and almost the best under various conditions, especially for ICE land cover. However, the overestimation of the two GLASS products in North Africa should be paid sufficient attention in use.

4. The values of the averaged multi-annual mean R_n over global land surface from nearly all remotely sensed products (except ISCCP-FH) were larger than that of the reanalysis products by about 10–30 Wm⁻². The top three values were from CERES4A (2000.03–2018), GLASS-AVHRR (1983–2018), and GLASS-MODIS (2000–2018) with 98.162, 95.527, and 92.729 Wm⁻², respectively.

5. Globally, the climatology of R_n from different products differed from each other, especially before the year of 2000, and no significant trend was observed during the whole past decades. However, some hotspots with significant trends in R_n agreed with most products were found mainly in the Amazon rainforest region, western Africa, Greenland, and most of the Antarctic, which are sensitive to climate change and human activity.

Overall, large discrepancies exist in all evaluated R_n products in performance, temporal mean, and temporal variations, and each of the ten R_n products has its advantages and disadvantages. Hence, it is necessary to take the characteristics and performance into account to select the right R_n product in practical use. However, due to the few available ground measurements in the Southern Hemisphere and tropical areas, the evaluation results might be impacted. In summary, more efforts should be made to improve the evaluation work and to conduct the uncertainty analysis to help to generate a better R_n product, and further to explore the factors to drive the variations in R_n particularly in the hotspots.

Authorship contribution statement

Xiuwan Yin: Investigation, Methodology, Software, Data curation, Visualization, Writing-original draft. Bo Jiang: Conceptualization, Methodology, Resources, Investigation, Data curation, Writing-review & editing, Funding acquisition, Supervision. Shunlin Liang: Conceptualization, Investigation, Supervision. Shaopeng Li, Qian Wang: Data curation. Xiang Zhao, Jianglei Xu, Jiakun Han, Hui Liang, Xiaotong Zhang, Qiang Liu, Yunjun Yao, Kun Jia, and Xianhong Xie: Discussion, review.

Acknowledgements

We appreciate the radiation measurements provided by the FLUXNET community, in particular by the following networks (listed in Table 2): LaThuile, AsiaFlux, ARM, BSRN, SURFRAD, GAME.AAN, BOREAS, GC-Net, CEOP, CEOPInt, SMOSREX, CERN, FIFE, LBA-ECO, Hi-Water, HAPEX-Sahel, SAFARI.2000, AmeriFlux (U.S. Department of Energy, Biological and Environmental Research, Terrestrial Carbon Program (DE-FG02-04ER63917)), ASIAQ, ChinaFlux, EOL, EFDC, IMAU, PROMICE, and FGI. The authors acknowledge the financial support to the eddy covariance data harmonization provided by CarboEuropeIP, FAO-GTOS-TCO, Ileaps, Max Planck Institute for Biogeochemistry, National Science Foundation, University of Tuscia, Université Laval, Environment Canada and US Department of Energy and the datasets development and technical support from Berkeley Water Center, Lawrence Berkeley National Laboratory, Microsoft Researche Science, Oak Ridge National Laboratory, University of California-Berkeley and the University of Virginia. We thank the data support from ‘National Earth System Science Data Center, National Science and Technology Infrastructure of China’ (https://www.geodata.cn, accessed on 31 January 2023), and the valuable suggestions made by the reviewers to help improve the quality of this paper.

Disclosure statement

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

Additional information

Funding

This work was funded by the National Natural Science Foundation of China [grant numbers 42090012 and 41971291].