Energy-conservation datasets of global land surface radiation and heat fluxes from 2000–2020 generated by CoSEB

Accurately estimating global land surface radiation [including downward shortwave radiation (SW_IN), downward longwave radiation (LW_IN), upward shortwave radiation (SW_OUT), upward longwave radiation (LW_OUT) and net radiation (Rn)] and heat fluxes [including latent heat flux (LE), soil heat flux (G) and sensible heat flux (H)] is essential for quantifying the exchange of radiation, heat and water between the land and atmosphere under global climate change. This study presents the first data-driven energy-conservation datasets of global land surface radiation and heat fluxes from 2000 to 2020, generated by our model of Coordinated estimates of land Surface Energy Balance components (CoSEB). The model integrates GLASS and MODIS remote sensing data, ERA5-Land reanalysis datasets, topographic data, CO₂ concentration data as independent variables and in situ radiation and heat flux observations at 258 eddy covariance sites worldwide as dependent variables within a multivariate random forest technique to effectively learn the physics of energy conservation. The developed CoSEB-based datasets are strikingly advantageous in that [1] they are the first data-driven global datasets that satisfy both surface radiation balance and heat balance among the eight fluxes, as demonstrated by both the radiation imbalance ratio [RIR, defined as $\mathrm{100}\mathit{\%}\times ({\text{SW}}_{\text{IN}}-{\text{SW}}_{\text{OUT}}+{\text{LW}}_{\text{IN}}-{\text{LW}}_{\text{OUT}}-\text{Rn})/\text{Rn}$] and energy imbalance ratio [EIR, defined as $\mathrm{100}\mathit{\%}\times (\text{Rn}-G-LE-H)/\text{Rn}$] of 0, [2] the radiation and heat fluxes are characterized by high accuracies, where (1) the RMSEs (R²) for daily estimates of SW_IN, SW_OUT, LW_IN, LW_OUT, Rn, LE, H and G from the CoSEB-based datasets at 44 independent test sites were 37.52 W m⁻² (0.81), 14.20 W m⁻² (0.42), 22.47 W m⁻² (0.90), 13.78 W m⁻² (0.95), 29.66 W m⁻² (0.77), 30.87 W m⁻² (0.60), 29.75 W m⁻² (0.44) and 5.69 W m⁻² (0.44), respectively, (2) the CoSEB-based datasets, in comparison to the mainstream products/datasets (i.e. GLASS, BESS-Rad, BESSV2.0, FLUXCOM, MOD16A2, PML_V2 and ETMonitor) that generally separately estimated subsets of the eight flux components, better agreed with the in situ observations. Our developed datasets hold significant potential for application across diverse fields such as agriculture, forestry, hydrology, meteorology, ecology, and environmental science, which can facilitate comprehensive studies on the variability, impacts, responses, adaptation strategies, and mitigation measures of global and regional land surface radiation and heat fluxes under the influences of climate change and human activities. The CoSEB-based datasets are open access and available through the National Tibetan Plateau Data Center (TPDC) at https://doi.org/10.11888/Terre.tpdc.302559 (Tang et al., 2025a) and through the Science Data Bank (ScienceDB) at https://doi.org/10.57760/sciencedb.27228 (Tang et al., 2025b).

1 Introduction

Land surface radiation balance and heat balance play important roles in Earth's climate system, representing the physical processes by which the surface-atmosphere absorbs and redistributes radiation and heat fluxes (Berbery et al., 1999; Betts et al., 1996; Mueller et al., 2009; Sellers et al., 1997; Xu et al., 2022a), and facilitating the exchange of water, energy, carbon, and other agents essential to climatic and ecological systems and human society (Jia et al., 2013; Wang et al., 2012; Wild, 2009; Wild et al., 2012; Xia et al., 2006). Accurately monitoring the spatial and temporal variations of global land surface radiation [including downward shortwave radiation (SW_IN), downward longwave radiation (LW_IN), upward shortwave radiation (SW_OUT), upward longwave radiation (LW_OUT) and net radiation (Rn)] and heat fluxes [including latent heat flux (LE), soil heat flux (G) and sensible heat flux (H)] is indispensable for quantifying the exchange of radiation, heat and water between the land and atmosphere under global climate change (Ersi et al., 2024; Liang et al., 2019; Rios and Ramamurthy, 2022; Tang et al., 2024a; Wang et al., 2021), and for studying solar energy utilization (Tang et al., 2024b; Zhang et al., 2017), hydrological cycle (Huang et al., 2015; Wild and Liepert, 2010), ecosystem productivity (Nemani et al., 2003), agricultural management (De Wit et al., 2005) and ecological protection (Tang et al., 2023). Remote sensing (RS) technology, with its high spatial-temporal resolution and applicability over large areas, is considered to be the most effective and economical means for obtaining global land surface radiation and heat fluxes (Liu et al., 2016; Van Der Tol, 2012; Zhang et al., 2010).

In past decades, numerous RS-based products/datasets of global surface radiation and heat fluxes have significantly advanced, which were generally generated by physical (Li et al., 2023; Mu et al., 2011; Yu et al., 2022) or statistical methods (Jiao et al., 2023; Jung et al., 2019; Peng et al., 2020). However, two key limitations still exist in these products. Firstly, most available products provide only a single component of land surface radiation or heat fluxes, e.g. ETMonitor (Zheng et al., 2022) and MOD16A2 (Mu et al., 2011) only estimating LE, leading to the failure to satisfy surface radiation balance and heat balance when the single radiation or heat flux is utilized in conjunction with products containing other radiation and heat components (Wang et al., 2025), and further posing significant uncertainties to understand the interactions and redistributions of surface radiation and energy in the Earth-atmosphere system. Secondly, a few products, e.g., FLUXCOM (Jung et al., 2019) and GLASS (Jiang et al., 2015; Zhang et al., 2014), generated datasets for multiple components of surface radiation and heat fluxes by using separate estimates from the uncoordinated models, which make them difficult to abide by surface radiation and heat conservation. These energy-imbalanced and radiation-imbalanced estimates among multiple components from previous products/datasets severely limit their in-depth applications in analyzing the spatial and temporal trends, simulating the physical processes of radiation, heat and water cycles as well as revealing the attributions and mechanisms in Earth-surface system under global climate change. It was imperative to develop global datasets of land surface radiation and heat fluxes characterized by high accuracy, radiation balance as well as heat balance, to better meet the requirements in practical applications of various fields.

Our proposed data-driven model/framework of Coordinated estimates of land Surface Energy Balance components (CoSEB) (Wang et al., 2025), which effectively learns the underlying physical interrelations (i.e., surface energy conservation law) among multiple targeted variables, provides an unprecedented opportunity to develop global datasets of land surface radiation and heat fluxes that can not only simultaneously provide high-accuracy estimates of these components but also adhere to surface radiation- and heat-conservation laws.

The objectives of this study are twofold: (1) to develop high-accuracy datasets of global land surface radiation and heat fluxes, which comply with the principles of radiation balance and heat balance, using our CoSEB model renewed based on in situ observations, remote sensing data and reanalysis datasets; (2) to validate the datasets/model estimates against data from in situ observations, mainstream products as well as estimates from uncoordinated random forest (RF) techniques. Section 2 introduces the data resources used in this study. Section 3 briefly describes the method we used to estimate global surface radiation and heat fluxes. Section 4 presents the evaluation of the datasets/model estimates generated by our renewed CoSEB model. Section 5 discusses the superiority, potential applications and uncertainties of the developed datasets. Data availability is given in Sect. 6, and a summary and conclusion is provided in Sect. 7.

2 Data

2.1 Ground-based observations

In this study, the in situ observations of land surface radiation and heat fluxes at 302 eddy covariance (EC) sites from the networks of AmeriFlux (174 sites, 2000–2020, https://AmeriFlux.lbl.gov/Data/, last access: 6 August 2024), EuroFlux (72 sites, 2000–2020, http://www.europe-fluxdata.eu/, last access: 6 August 2024), OzFlux (5 sites, 2007–2012, https://data.ozflux.org.au/, last access: 6 August 2024), FLUXNET (108 sites, 2000–2014, https://FLUXNET.org/Data/download-Data/, last access: 6 August 2024), JapanFlux (15 sites, 2001–2020, https://ads.nipr.ac.jp/japan-flux2024/, last access: 10 October 2025), ChinaFLUX (5 sites, 2005–2020, http://www.chinaflux.org/, last access: 6 August 2024) and National Tibetan Plateau/Third Pole Environment Data Center (TPDC, 13 sites, 2012–2020, https://Data.tpdc.ac.cn/en/Data, last access: 6 August 2024) were used (Fig. 1), where 37, 48 and 5 sites in FLUXNET were also shared in AmeriFlux, EuroFlux and OzFlux, respectively. These 302 sites were filtered out from all collected 1098 sites by following the quality-assurance and quality-control steps, including: (1) any site with a missing component of any of the SW_IN, SW_OUT, LW_IN, LW_OUT, LE, H and G was excluded, reducing the 1098 sites to 472 sites for further analysis; (2) any half-hour period with missing data for any of these components was excluded; (3) the half-hourly ground-based observations with quality-control flag of 2 or 3 (bad quality) were removed but quality-control flag of 0 and 1 (good quality) were maintained; (4) a daily average of the half-hour observations was calculated for each day with greater than 80 % good-quality data, further reducing the 472 sites to 355 sites; (5) the aggregated daily LE and H were corrected for energy imbalance using the Bowen ratio method when the daily energy balance closure [defined as $(LE+H)/(\text{Rn}-G)$] varied between 0.2 and 1.8 following Wang et al. (2025) to exclude physically implausible measurements; (6) extreme outliers in the daily evaporative fraction were further removed by excluding values outside the 1st–99th percentile range, a common practice in flux and remote sensing studies (Bartkowiak et al., 2024; Wang et al., 2023), further reducing the 355 sites to 337 sites. Besides, the RS products/datasets involved in this study collocated at the sites should not be missing, finally reducing the 337 sites to 302 sites for analysis. Note that the Rn at these sites used in this study was calculated from the sum of net longwave radiation (LW_IN minus LW_OUT) and net shortwave radiation (SW_IN minus SW_OUT), rather than using the observed Rn directly, to ensure surface radiation balance in training datasets.

Figure 1 Spatial distribution of the 302 eddy covariance sites from AmeriFlux, FLUXNET, EuroFlux, OzFlux, JapanFlux, ChinaFLUX and TPDC, and nine radiation sites from SURFRAD involved for analysis in this study.

These 302 sites used in this study cover a wide range of global climate regimes across 14 land cover types, including (1) evergreen needleleaf forests (ENF, 55 sites); (2) evergreen broadleaf forests (EBF, 12 sites); (3) deciduous needleleaf forests (DNF, 7 sites); (4) deciduous broadleaf forests (DBF, 40 sites); (5) mixed forests (MF, 8 sites); (6) closed shrublands (CSH, 5 sites); (7) open shrublands (OSH, 11 sites); (8) woody savannas (WSA, 6 sites); (9) savannas (SAV, 10 sites); (10) grasslands (GRA, 62 sites); (11) permanent wetlands (WET, 22 sites); (12) croplands (CRO, 59 sites); (13) water bodies (WAT, 1 sites); (14) cropland/natural vegetation mosaics (CVM, 4 sites). Among them, 44 sites (∼15 % of the total, see Table S1 in the Supplement) were isolated to serve as spatially independent sites to test the generated datasets and they did not participate in the development of the model/datasets.

Furthermore, ground-based radiation observations from nine sites that are located in large flat agricultural areas covered by crops and grasses from SURFRAD (https://gml.noaa.gov/, last access: 30 December 2025) were also introduced to validate land surface radiation estimates. Similar to the preprocessing performed on the observations of the 302 EC sites, the SW_IN, SW_OUT, LW_IN, LW_OUT and Rn from the SURFRAD were also quality-controlled and aggregated to daily data. Spatial distribution of the 302 EC sites and nine radiation sites from SURFRAD are shown in Fig. 1, with site details (latitude, longitude, land cover types, digital elevation model and temporal coverage) provided in Tables S1 and S2.

2.2 Climate/meteorology and remote sensing data

To generate global datasets of land surface radiation and heat fluxes from 2000 to 2020, five types of climate/meteorology and remote sensing data were used in this study, including:

1. ERA5-Land reanalysis datasets (https://cds.climate.copernicus.eu/, last access: 6 August 2024) with the spatial resolution of ∼9 km from 1950 (Muñoz-Sabater et al., 2021). Following our previous work (Wang et al., 2025), this study used variables from the ERA5-Land datasets to drive the model, including near-surface 2 m air temperature (T_a), soil temperature in layer 1 (0–7 cm, T_S1), soil volumetric moisture content in layer 1 (0–7 cm, SM1), solar radiation reaching the surface of the earth (${\text{SW}}_{\text{IN}}^{\text{ERA5}}$), net thermal radiation at the surface (LW_net), pressure of the atmosphere (PA), 10 m wind speed (WS), precipitation (P_r) and the 2 m dewpoint temperature, daily minimum and maximum air temperature [for calculating relative air humidity (RH)].

2. GLASS datasets (https://glass.bnu.edu.cn/, last access: 6 August 2024), which provide the 500 m 8 d leaf area index (LAI) and fractional vegetation cover (FVC) from February 2000 to December 2021.

3. MOD44B product (https://lpdaac.usgs.gov/, last access: 6 August 2024), which offers yearly 250 m percent tree cover (PTC) since 2000, representing the percentage (0 %–100 %) of a pixel covered by tree canopy.

4. NOAA/GML atmospheric carbon dioxide (CO₂) concentration data, providing monthly global marine surface mean data since 1958 (ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_gl.txt, last access: 6 August 2024).

5. GMTED2010 topographic data (https://topotools.cr.usgs.gov/gmted_viewer/gmted2010_global_grids.php, last access: 6 August 2024), providing 500 m digital elevation model (DEM), slope, and aspect.

The ∼9 km ERA5-Land datasets were spatially interpolated to 500 m using the cubic convolution method, and the 250 m PTC was resampled to 500 m using the arithmetic averaging method.

2.3 Mainstream datasets/products for inter-comparison

Mainstream RS-based datasets/products of moderate-resolution global land surface radiation and heat fluxes were collected for inter-comparison (Table 1), including (1) the daily 0.05° GLASS SW_IN, LW_IN, LW_OUT and Rn products from 2000 to 2018 (https://glass.bnu.edu.cn/, last access: 6 August 2024), (2) the daily 0.05° Breathing Earth System Simulator Radiation (BESS-Rad) SW_IN products from 2000 to 2020 (https://www.environment.snu.ac.kr/bess-rad, last access: 30 December 2025), (3) the daily 0.05° BESS Version2.0 (BESSV2.0) Rn and LE products from 2000 to 2020 (https://www.environment.snu.ac.kr/bessv2, last access: 30 December 2025), (4) the 8 d 0.0833° FLUXCOM Rn, LE and H products from 2001 to 2020 (https://fluxcom.org/, last access: 6 August 2024), (5) the daily 1 km ETMonitor LE product from 2000 to 2020 (https://data.casearth.cn/, last access: 6 August 2024), (6) the 8 d 500 m Penman-Monteith-Leuning Version2 (PML_V2, https://www.tpdc.ac.cn/, last access: 6 August 2024) LE product from 2000 to 2020; and (7) the 8 d 500 m MOD16A2 (https://lpdaac.usgs.gov/, last access: 6 August 2024) LE product from 2000 to 2020.

Table 1 Summary of mainstream datasets/products for inter-comparison used in this study.

The GLASS SW_IN products are derived from a combination of the GLASS broadband albedo product and the surface shortwave net radiation estimates, where the surface shortwave net radiation is estimated using linear regression with MODIS top-of-atmosphere (TOA) spectral reflectance (Wang et al., 2015). The GLASS LW_IN and LW_OUT products are generated using densely connected convolutional neural networks, incorporating Advanced Very High-Resolution Radiometer (AVHRR) TOA reflectance and ERA5 near-surface meteorological data (Xu et al., 2022b). The GLASS Rn products are estimated from the meteorological variables from MERRA2 and surface variables from GLASS using the multivariate adaptive regression splines model (Jiang et al., 2015). The BESS-Rad and BESSV2.0 estimate SW_IN and Rn using a radiative transfer model (i.e., Forest Light Environmental Simulator, FLiES) with an artificial neural network based on MODIS and MERRA2 reanalysis datasets, and using FLiES based on MODIS products and NCEP/NCAR reanalysis data, respectively (Li et al., 2023; Ryu et al., 2018). Moreover, the BESSV2.0 (Li et al., 2023), MOD16A2 (Mu et al., 2011), PML_V2 (Zhang et al., 2019) and ETMonitor (Zheng et al., 2022) generated global LE by physical models, such as Penman-Monteith equation, Priestley-Taylor equation and/or Shuttleworth-Wallace two-source scheme. The FLUXCOM Rn, LE and H datasets are obtained through multiple machine learning methods based on in situ observations from FLUXNET and remote sensing and meteorological data (Jung et al., 2019). For better consistency, RF-based 8 d 0.0833° Rn and Bowen ratio-corrected LE and H for the periods of 2000 to 2020 from the FLUXCOM were used in this study.

3 Methods

The method used to generate global datasets of land surface radiation and heat fluxes is based on the CoSEB model/framework, which was developed by our previously published work (Wang et al., 2025), to coordinately estimate global land surface energy balance components (including Rn, LE, H and G) using the multivariate random forest technique, with a combination of MODIS and GLASS products, ERA5-Land reanalysis datasets, and in situ observations at 336 EC sites. The CoSEB model was demonstrated to be able to produce high-accuracy estimates of land surface energy components, with the RMSE of <17 W m⁻² and R² of >0.83 for estimating 4 d Rn, LE and H, and the RMSE of <5 W m⁻² and R² of 0.55 for estimating 4 d G. The most praiseworthy superiority of the CoSEB model lies in its ability to balance the land surface energy components, with an energy imbalance ratio [EIR, defined as $\mathrm{100}\mathit{\%}\times (\text{Rn}-G-LE-H)/\text{Rn}$] of 0.

To coordinately estimate land surface radiation and heat fluxes that comply with both radiation balance and heat balance, one of the key procedures in the construction of the CoSEB model was to prepare training datasets that satisfy surface radiation and heat balance. For this purpose, the energy-imbalance corrections on daily in situ observed LE and H were conducted by the most widely applied Bowen ratio method [$H^{\text{corr}}=\frac{H}{H+LE}\times (\text{Rn}-G)$, $L{E}^{\text{corr}}=\frac{LE}{H+LE}\times (\text{Rn}-G)$, where H^corr and LE^corr represent the sensible heat flux and latent heat flux after energy-imbalance correction, respectively] with the aid of Rn and G observations, and the in situ Rn was calculated from the sum of in situ observed net longwave radiation (LW_IN minus LW_OUT) and net shortwave radiation (SW_IN minus SW_OUT). The input variables to renew the CoSEB model include: (1) climate/meteorology: T_a, ${\text{SW}}_{\text{IN}}^{\text{ERA5}}$, LW_net, WS, PA, P_r, RH, CO₂ concentration; (2) vegetation and soil: LAI, FVC, PTC, T_S1, SM1; (3) topography data: DEM, Slope and Aspect, in addition to longitude (Lon), latitude (Lat), and inverse relative distance from the Earth to the Sun (dr), in which the dr was calculated as $\text{dr}=\mathrm{1}+\mathrm{0.033}\times \cos \left(\frac{\mathrm{2}\mathit{\pi }\times \text{DOY}}{\mathrm{365}}\right)$, where DOY represents the day of year. Considering that the footprint of the site-based measurements of turbulent heat fluxes is generally at a scale of hundreds of meters, to reduce the effect of differences of spatial scales between ground-based measurements (dependent variables) and remotely sensed/reanalysis datasets (independent variables), we renewed the CoSEB model at a spatial scale of 500 m for coordinately estimating global daily land surface radiation and heat fluxes, which can be expressed as follows:

$$\begin{array}{}\text{(1)}& \left(\begin{array}{l}{\text{SW}}_{\text{IN}},{\text{SW}}_{\text{OUT}},\\ {\text{LW}}_{\text{IN}},{\text{LW}}_{\text{OUT}},\\ \text{Rn},LE,H,G\end{array}\right)=f\left(\begin{array}{l}\text{Lon},\text{Lat},T_{\mathrm{a}},T_{\text{S1}},\\ \text{SM1},{\text{SW}}_{\text{IN}}^{\text{ERA5}},\\ {\text{LW}}_{\text{net}},\text{PA},\text{WS},P_r,\\ \text{dr},\text{RH},\text{LAI},\text{FVC},\\ \text{PTC},\text{DEM},\text{Slope},\\ \text{Aspect},{\mathrm{CO}}_{\mathrm{2}}\end{array}\right)\end{array}$$

To enhance model generalization, the renewed CoSEB model was reoptimized using random and grid search methods, resulting in different hyperparameters of 281 decision trees, a maximum depth of 21, and minimum samples split and leaf of 8 from those of Wang et al. (2025). Site-based 10-fold cross-validation was employed to evaluate the transferability and generalization of the CoSEB model by randomly dividing all sites into ten folds, where the samples from each fold of sites in turn served as validation datasets while the remaining folds were used as training datasets, ensuring that the validation was conducted on sites spatially independent from the training data. Furthermore, to benchmark the coordinated estimates from the renewed CoSEB model, eight RF-based uncoordinated models were constructed, each separately estimating one of SW_IN, SW_OUT, LW_IN, LW_OUT, Rn, LE, H or G using the same inputs as those in the renewed CoSEB model. Figure 2 illustrates the flowchart for generating global datasets of land surface radiation and heat fluxes by the CoSEB model.

Figure 2 Flowchart for generating energy-conservation datasets of global land surface radiation [including downward shortwave radiation (SW_IN), downward longwave radiation (LW_IN), upward shortwave radiation (SW_OUT), upward longwave radiation (LW_OUT) and net radiation (Rn)] and heat fluxes [including latent heat flux (LE), soil heat flux (G) and sensible heat flux (H)] by the CoSEB model renewed from in situ observations at 258 sites worldwide and collocated remote sensing and reanalysis datasets.

4 Results

4.1 Validation of the CoSEB model

4.1.1 Site-based 10-fold cross-validations at 258 EC sites

Figures 3 and 4 present the scatter density plots of the site-based 10-fold cross-validation of daily SW_IN, LW_IN, SW_OUT, LW_OUT, Rn, LE, H and G estimated from the renewed CoSEB model and the RF-based uncoordinated models, respectively, by using the validation datasets collected at 258 EC sites worldwide. Results indicated that the estimates from both the CoSEB model and the RF-based uncoordinated models agreed well with the in situ observations, with the coefficient of determination (R²) varying between 0.80 and 0.95 for SW_IN, LW_IN, LW_OUT and Rn, and between 0.59 and 0.67 for SW_OUT, LE and H. The CoSEB model, with the root mean square error (RMSE) of 26.82 to 34.25 W m⁻² and mean absolute error (MAE) of 18.83 to 24.49 W m⁻² for SW_IN, Rn, LE and H, the RMSE of 12.24 to 17.75 W m⁻² and the MAE of 8.39 to 13.70 W m⁻² for SW_OUT, LW_IN and LW_OUT, demonstrated comparable accuracies to the RF-based models, with the RMSE of 27.07 to 33.34 W m⁻² and MAE of 19.29 to 23.64 W m⁻² for SW_IN, Rn, LE and H, the RMSE of 12.12 to 16.93 W m⁻² and the MAE of 8.68 to 12.99 W m⁻² for SW_OUT, LW_IN and LW_OUT. In the validation of daily G, both the CoSEB and RF-based models yielded RMSEs below 7 W m⁻². Comparisons with the corresponding training results (Table S3) indicated that although the CoSEB model performed better on the training datasets, its overall performance remained stable, suggesting that the CoSEB model was not affected by overfitting.

Figure 3 Scatter density plots of the site-based 10-fold cross-validation of daily downward shortwave and longwave radiation (SW_IN and LW_IN), upward shortwave and longwave radiation (SW_OUT and LW_OUT), net radiation (Rn), soil heat flux (G), latent heat flux (LE) and sensible heat flux (H) derived by the CoSEB model against in situ observed SW_IN, LW_IN, SW_OUT, LW_OUT, Rn, G, and energy imbalance-corrected LE ($L{E}_{\text{daily}}^{\text{corr}}$) and H ($H_{\text{daily}}^{\text{corr}}$). The EIR and RIR in the subfigure (i) represent the energy imbalance ratio and radiation imbalance ratio, which are defined as $\mathrm{100}\mathit{\%}\times (\text{Rn}-G-LE-H)/\text{Rn}$ and $\mathrm{100}\mathit{\%}\times ({\text{SW}}_{\text{IN}}-{\text{SW}}_{\text{OUT}}+{\text{LW}}_{\text{IN}}-{\text{LW}}_{\text{OUT}}-\text{Rn})/\text{Rn}$, respectively. The colorbar represents the normalized density of data points.

Figure 4 Same as Fig. 3, but for estimates from RF-based uncoordinated models.

Strikingly, the CoSEB model exhibited large superiority in balancing the surface radiation and heat fluxes, with the radiation imbalance ratio [RIR, defined as $\mathrm{100}\mathit{\%}\times ({\text{SW}}_{\text{IN}}-{\text{SW}}_{\text{OUT}}+{\text{LW}}_{\text{IN}}-{\text{LW}}_{\text{OUT}}-\text{Rn})/\text{Rn}$] and energy imbalance ratio [EIR, defined as $\mathrm{100}\mathit{\%}\times (\text{Rn}-G-LE-H)/\text{Rn}$] of 0, while the RF-based uncoordinated models showed substantial imbalances of the surface radiation and heat fluxes, with RIR and EIR that were approximately normally distributed, having absolute mean values of 38.84 % and 31.22 %, respectively, and reaching as high as 50 % in some cases. Furthermore, the RIR as well as EIR tended to be higher under lower solar radiation, air temperature, or FVC, with more frequent low values of these three variables leading to a broader and less peaked distribution of RIR and EIR (see Fig. S1).

4.1.2 Validation at nine radiation sites from SURFRAD

To further illustrate the generality and transferability of the renewed CoSEB model, the validation of estimates of the five radiation components (including SW_IN, SW_OUT, LW_IN, LW_OUT, Rn) derived from both the CoSEB model and RF-based uncoordinated models against observations at nine radiation sites from SURFRAD was performed, as shown in Fig. 5. The results showed that both the CoSEB model and the RF-based models achieved high accuracy in estimating daily SW_IN, SW_OUT, LW_IN, LW_OUT and Rn, with the RMSE of ∼30 W m⁻² for SW_IN, ∼14 W m⁻² for SW_OUT and LW_IN, ∼12 W m⁻² for LW_OUT and ∼24 W m⁻² for Rn, with the R²>0.9 for SW_IN, LW_IN and LW_OUT, ∼0.65 for SW_OUT and ∼0.85 for Rn. Compared to the results of the site-based 10-fold cross-validation at 258 EC sites, the performances at nine radiation sites showed slight improvements, with the RMSE decreasing by 0.74 to 4.54 W m⁻² for SW_IN, LW_IN, LW_OUT and Rn in the CoSEB model, but a slight degradation with the RMSE increasing by ∼1.05 W m⁻² for SW_OUT, suggesting the robust performance of the CoSEB model. Furthermore, the CoSEB model demonstrated a large superiority in maintaining surface radiation balance among the five radiation components, with the RIR of 0, in contrast to the RF-based models, which failed to meet this balance, exhibiting significant RIR exceeding 50 %.

Figure 5 Scatter density plots of the validation of daily downward shortwave and longwave radiation (SW_IN and LW_IN), upward shortwave and longwave radiation (SW_OUT and LW_OUT) and net radiation (Rn) from the renewed CoSEB model (upper two rows) and RF-based uncoordinated models (lower two rows) against in situ observations at nine radiation sites from SURFRAD. The RIR represents the radiation imbalance ratio, defined as $\mathrm{100}\mathit{\%}\times ({\text{SW}}_{\text{IN}}-{\text{SW}}_{\text{OUT}}+{\text{LW}}_{\text{IN}}-{\text{LW}}_{\text{OUT}}-\text{Rn})/\text{Rn}$. The colorbar represents the normalized density of data points.

4.2 Validation and inter-comparisons of the CoSEB-based datasets

As demonstrated in Sect. 4.1, the renewed CoSEB model with a spatial scale of 500 m achieved comparable accuracies to the RF-based uncoordinated models but outperformed them in balancing surface radiation and heat fluxes. Evidenced by the validation for its superiority, the renewed CoSEB model was then applied to the spatially aggregated input datasets to generate our developed global daily datasets with a spatial resolution of 0.05°. To further assess the performance of the developed CoSEB-based datasets, in situ observations from another 44 spatially independent test sites (see Sect. 2.1), which were not involved in model construction and datasets generation, were used for validation. Mainstream products (i.e. GLASS, BESS-Rad, BESSV2.0, FLUXCOM, PML_V2, MOD16A2 and ETMonitor) were also involved for inter-comparison at the 44 test sites.

Note that due to the lack of moderate-resolution global RS-based products/datasets of daily and/or 8 d SW_OUT, H and G, the intercomparison between different products/datasets was impossible. Instead, we conducted a validation of these components from the CoSEB-based datasets against in situ observations at 44 test sites, as shown in Figs. S2 and S3. Results indicated that the CoSEB-based datasets could provide good estimates of SW_OUT, H and G, with the RMSEs (R²) of 14.20 W m⁻² (0.42), 29.75 W m⁻² (0.44) and 5.69 W m⁻² (0.44) at daily scale, respectively, and the RMSE (R²) of 12.19 W m⁻² (0.39) and 4.60 W m⁻² (0.47) for 8 d SW_OUT and G, respectively.

Figures 6 and 7 present the comparison of daily SW_IN, LW_IN and LW_OUT, as well as Rn and LE from the CoSEB-based datasets and mainstream products/datasets (including GLASS, BESS-Rad, BESSV2.0 and ETMonitor), with in situ observations at 44 test sites, respectively. Overall, the estimates from the CoSEB-based datasets exhibited a closer agreement with in situ observations than those from mainstream products/datasets, where the CoSEB-based datasets reduced the RMSE by 0.01 to 4.58 W m⁻² and increased the R² by 0.01 to 0.09 compared to mainstream products. Specifically, the RMSE for the SW_IN, LW_IN, and LW_OUT increased from 37.52, 22.47 and 13.78 W m⁻² in the CoSEB-based datasets to 37.53, 23.37 and 16.46 W m⁻² in the GLASS, respectively, and for SW_IN from 37.52 W m⁻² in the CoSEB-based datasets to 40.87 W m⁻² in the BESS-Rad. Likewise, the RMSEs for daily Rn and LE were 29.66 and 30.87 W m⁻² in the CoSEB-based datasets, which were lower than those of 34.24 and 34.36 W m⁻² in BESSV2.0, respectively, as well as those of 30.60 W m⁻² for Rn in GLASS and 33.62 W m⁻² for LE in ETMonitor.

Figure 6 Comparison of the daily downward shortwave radiation (SW_IN, the first column), downward longwave radiation (LW_IN, the second column) and upward longwave radiation (LW_OUT, the third column) from the CoSEB-based datasets, GLASS and BESS-Rad with the in situ observed SW_IN, LW_IN and LW_OUT at 44 test sites. The colorbar represents the normalized density of data points.

Figure 7 Comparison of the daily net radiation (Rn, the upper row) and latent heat flux (LE, the lower row) from the CoSEB-based datasets, BESSV2.0, GLASS and ETMonitor with the in situ observed Rn, and energy imbalance-corrected LE ($L{E}_{\text{daily}}^{\text{corr}}$) at 44 test sites. The colorbar represents the normalized density of data points.

Figures 8–10 compare the 8 d SW_IN, LW_IN and LW_OUT, Rn and LE, as well as H from the CoSEB-based datasets and mainstream products, with in situ observations at 44 test sites, respectively. Overall, the CoSEB-based datasets outperformed the mainstream products/datasets for all surface radiation and heat fluxes, where the CoSEB-based datasets reduced the RMSE by 0.24 to 10.48 W m⁻² and increased the R² by 0.01 to 0.38 compared to mainstream products. Specifically, for SW_IN, LW_IN and LW_OUT, the RMSE increased from 18.54, 18.50 and 9.41 W m⁻² in the CoSEB-based datasets to 21.35, 20.39 and 14.48 W m⁻² in the GLASS, respectively, and for SW_IN from 18.54 W m⁻² in the CoSEB-based datasets to 18.78 W m⁻² in the BESS-Rad. For Rn, the RMSE increased from 19.12 W m⁻² in the CoSEB-based datasets to ∼23 W m⁻² in the FLUXCOM and GLASS and to >27 W m⁻² in the BESSV2.0, while the R² decreased from 0.82 in the CoSEB-based datasets to 0.75 in the FLUXCOM and GLASS and to 0.62 in the BESSV2.0. Likewise, for LE, the RMSE increased from 22.31 W m⁻² in the CoSEB-based datasets to ∼25 W m⁻² in the FLUXCOM, PML_V2, BESSV2.0 and ETMonitor, and to >32 W m⁻² in MOD16A2, while the R² decreased from 0.67 in the CoSEB-based datasets to ∼0.60 in the FLUXCOM, PML_V2, BESSV2.0 and ETMonitor, and to <0.3 in the MOD16A1. For H, the RMSE increased from 21.63 W m⁻² in the CoSEB-based datasets to 22.64 W m⁻² in the FLUXCOM.

Figure 8 Same as Fig. 6, but for the comparison at 8 d scale.

Figure 9 Comparison of the 8 d net radiation (Rn, the upper two rows) and latent heat flux (LE, the lower three rows) from the CoSEB-based datasets, FLUXCOM, BESSV2.0, GLASS, MOD16A2, PML_V2 and ETMonitor with in situ observed Rn, and energy imbalance-corrected LE ($L{E}_{\text{8\hspace{0.17em}d}}^{\text{corr}}$) at 44 test sites. The colorbar represents the normalized density of data points.

Figure 10 Comparison of the 8 d sensible heat flux (H) from the CoSEB-based datasets and the FLUXCOM with the in situ energy imbalance-corrected H ($H_{\text{8\hspace{0.17em}d}}^{\text{corr}}$) at 44 test sites. The colorbar represents the normalized density of data points.

The differences between the estimates from the CoSEB-based datasets and mainstream datasets are likely multifactorial, arising from the simplification and parameterization uncertainties in physics-based models, as well as the lack of physical constraints, limited training samples, and incomplete consideration of influencing factors in other machine-learning-based models.

4.3 Spatial-temporal patterns of global land surface radiation and heat fluxes

In addition to the validation and inter-comparison of the CoSEB-based datasets at the site scale, we further inter-compared the estimates of land surface radiation and heat fluxes from the CoSEB-based datasets and the mainstream products/datasets, in terms of their global spatial and temporal patterns.

Figures 11–13 show the spatial distributions (excluding Greenland, Antarctic continent, deserts, water bodies and permanent snow) and latitudinal profiles of the global 0.05° mean annual SW_IN, LW_IN and LW_OUT, Rn and LE, as well as H from 2001 to 2018, respectively, as derived from the CoSEB-based datasets and mainstream products/datasets [i.e. GLASS, BESS-Rad, BESSV2.0, FLUXCOM, MOD16A2, PML_V2 and ETMonitor, resampled to 0.05° using arithmetic averaging method or cubic convolutional method if necessary]. Overall, the spatial patterns of the estimates from the CoSEB-based datasets aligned well with those observed in these mainstream products/datasets, though regional discrepancies were present. Specifically, the mean annual LW_IN, LW_OUT, Rn, and LE generally exhibited decreasing trends from the equator towards higher latitudes, peaking in regions such as the Amazon Rainforest, Congo Rainforest, and the Malay Archipelago. In contrast, the higher mean annual SW_IN and H were mainly found in the Tibetan Plateau, southwestern U.S., mid-west Australia, Sahel and Southern Africa, while the lower values were found in high-latitude regions of >50° N. In the region with high values, the mean annual estimates of SW_IN from the CoSEB-based datasets were higher than those from GLASS but lower than those from BESS-Rad, the estimates of LW_IN and LW_OUT from the CoSEB-based datasets were both higher than those from GLASS, the estimates of Rn from the CoSEB-based datasets were significantly higher than those from BESSV2.0, and comparable to or slightly higher than those from FLUXCOM and GLASS, the estimates of LE from the CoSEB-based datasets were close to those from BESSV2.0 and PML_V2, but slightly lower than those from FLUXCOM, MOD16A2 and ETMonitor. Besides, the estimates of H from the CoSEB-based datasets were higher than those from FLUXCOM in regions with high values, while lower than those from FLUXCOM in regions with low values.

Figure 11 Spatial patterns of global mean annual downward shortwave radiation (SW_IN, the first row), downward longwave radiation (LW_IN, the second row) and upward longwave radiation (LW_OUT, the third row) from 2001 to 2018 by CoSEB-based datasets, GLASS and BESS-Rad. The rightmost subfigure of each row represents the latitudinal profiles of mean annual SW_IN, LW_IN and LW_OUT from CoSEB-based datasets, GLASS and BESS-Rad, where the shaded area represents the variation of standard deviation for each product.

Figure 12 Spatial patterns of global mean annual net radiation (Rn, the first row) and latent heat flux (LE, the second and third rows) from 2001 to 2018 by CoSEB-based datasets, FLUXCOM, BESSV2.0, MOD16A2, PML_V2, ETMonitor and GLASS. The last two subfigures of the third row represent the latitudinal profiles of mean annual Rn and LE from CoSEB-based datasets and these mainstream products/datasets, where the shaded area represents the variation of standard deviation for each product.

Figure 13 Spatial patterns of global mean annual sensible heat flux (H) from 2001 to 2018 by CoSEB-based datasets and FLUXCOM. The rightmost subfigure represents the latitudinal profiles of mean annual H from CoSEB-based datasets and FLUXCOM, where the shaded area represents the variation of standard deviation for each product.

The temporal evolutions of the global (excluding Greenland, Antarctic continent, deserts, water bodies and permanent snow) land surface radiation and heat fluxes derived from the CoSEB-based datasets and mainstream products/datasets from 2001 to 2018 were also investigated, as shown in Fig. 14. The results indicated that the temporal variation of each flux from the CoSEB-based datasets generally agreed well with those from mainstream products/datasets. The global annual mean estimates using area weighting average by the CoSEB-based datasets from 2001 to 2018 varied between ∼185.22 and ∼189.50 W m⁻² with the mean of ∼187.23 W m⁻² for SW_IN, between ∼32.67 and ∼33.20 W m⁻² with the mean of ∼32.96 W m⁻² for SW_OUT, between ∼330.24 and ∼334.14 W m⁻² with the mean of ∼331.50 W m⁻² for LW_IN, between ∼387.25 and ∼390.82 W m⁻² with the mean of ∼388.81 W m⁻² for LW_OUT, between ∼95.41 and ∼99.39 W m⁻² with the mean of 97.11 W m⁻² for Rn, between ∼53.24 and ∼56.37 W m⁻² with the mean of ∼54.53 W m⁻² for LE, between ∼40.44 and ∼41.96 W m⁻² with the mean of ∼41.29 W m⁻² for H, and between ∼1.22 and ∼1.52 W m⁻² with the mean of ∼1.33 W m⁻² for G. For each radiation or heat flux, the annual mean estimates from the CoSEB-based datasets were overall higher than those from the mainstream products/datasets. In particular, the annual mean Rn estimates from the CoSEB-based datasets were higher than those from FLUXCOM, GLASS and BESSV2.0 sequentially, and the annual mean LE estimates from the CoSEB-based datasets were marginally higher than those from FLUXCOM, but substantially exceeded those from ETMonitor, PML_V2, MOD16A2 and BESSV2.0 sequentially. The anomaly-based analyses (Fig. S4) reveal clear and coherent temporal trends of these radiation and heat fluxes, which respond well to global climate change, such as increasing atmospheric CO₂ and rising air temperatures.

Figure 14 Temporal variation of annual mean downward shortwave radiation (SW_IN), upward shortwave radiation (SW_OUT), downward longwave radiation (LW_IN), upward longwave radiation (LW_OUT), net radiation (Rn), latent heat flux (LE), sensible heat flux (H) and soil heat flux (G) from 2001 to 2018 from the CoSEB-based datasets, BESS-Rad, GLASS, FLUXCOM, BESSV2.0, PML_V2, MOD16A2 and ETMonitor. The shaded area represents the variation of the standard deviation for each product.

Figures 15–17 show the spatial patterns (excluding Greenland, Antarctic continent, deserts, water bodies and permanent snow) of interannual variability of SW_IN, LW_IN and LW_OUT, Rn and LE, as well as H from 2001 to 2018, respectively, derived from the CoSEB-based datasets and mainstream products/datasets. In general, the estimates from the CoSEB-based datasets displayed similar interannual variability in space with those from the mainstream products/datasets. Specially, the estimates of SW_IN from the CoSEB-based datasets, BESS-Rad, and GLASS exhibited a significant interannual variability mainly in northeastern Australia, eastern South America, Southeast China, and Southwest North America. The interannual variability of LW_IN and LW_OUT by the CoSEB-based datasets and GLASS displayed high values primarily at middle-to-high latitudes of the Northern Hemisphere and parts of Africa and Australia. The interannual variability of Rn observed by the CoSEB-based datasets was generally lower than that of GLASS, but higher than that of BESSV2.0 and FLUXCOM. The CoSEB-based datasets missed the strong interannual variability of LE as observed in MOD16A2, PML_V2 and ETMonitor in parts of Africa, Australia and eastern South America. Furthermore, FLUXCOM exhibited the weakest interannual variability of LE in almost all regions. The interannual variability of H derived from the CoSEB-based datasets was higher than that from FLUXCOM, with stronger interannual variabilities mainly observed in parts of eastern South America, southern Africa, and northeastern Australia.

Figure 15 Spatial distribution of interannual variability (standard deviation) of downward shortwave radiation (SW_IN, the first row), downward longwave radiation (LW_IN, the second row) and upward longwave radiation (LW_OUT, the third row) from 2001 to 2018 by the CoSEB-based datasets, GLASS and BESS-Rad.

Figure 16 Spatial distribution of interannual variability (standard deviation) of net radiation (Rn, the first and second rows) and latent heat flux (LE, the third and fourth rows) from 2001 to 2018 by the CoSEB-based datasets, FLUXCOM, BESSV2.0, MOD16A2, PML_V2, ETMonitor and GLASS.

Figure 17 Spatial distribution of interannual variability (standard deviation) of sensible heat flux (H) from 2001 to 2018 by the CoSEB-based datasets and FLUXCOM.

5 Discussion

Accurately monitoring the spatial and temporal variations of global land surface radiation and heat fluxes is crucial for quantifying the exchange of radiation, heat and water between the land and atmosphere under global climate change (Chen et al., 2020; Du et al., 2024; Kim et al., 2023; Liang et al., 2006; Wang et al., 2020). However, although numerous global RS-based products/datasets of land surface radiation and heat fluxes have been developed using physical and/or statistical methods, they typically provide either merely a single flux or multiple fluxes (see Table 1) that are estimated separately from uncoordinated models (Huang et al., 2024; Jung et al., 2019; Sun et al., 2023; Tang et al., 2019), leading to noticeable radiation imbalance and/or heat imbalance when these products are combined for practical applications. To address these limitations, we generated high-accuracy global datasets of land surface radiation and heat fluxes from 2000 to 2020 that adhere to both radiation and heat conservation laws, using our proposed CoSEB model (Wang et al., 2025).

Our CoSEB model, integrating underlying physical principles of training datasets into machine learning technique to effectively learn the interrelations among multiple targeted outputs, was originally designed for coordinating estimates of global land surface energy balance components (Rn, LE, H and G) to satisfy the energy conservation (Wang et al., 2025). Inspired by the idea of constructing the original CoSEB model, we further incorporated land surface radiation fluxes into our model to simultaneously consider the physical constraints of both surface radiation and heat conservation principles, by renewing the CoSEB using multiple remote sensing and reanalysis datasets, as well as in-situ observations of SW_IN, SW_OUT, LW_IN, LW_OUT, Rn, LE, H and G. In selecting the 19 input variables to accommodate the additional target variables, prior knowledge derived from previous studies was employed to identify factors that exert significant influence on surface radiation and heat flux while maintaining relative inter-independence as much as possible (Jung et al., 2019; Mohan et al., 2020; Wang et al., 2021; Xian et al., 2024). This practice is commonly adopted in data-driven models for estimating land surface water, energy, and carbon fluxes (Bai et al., 2024; Elghawi et al., 2023; Han et al., 2023; O and Orth, 2021). The importance scores of the 19 different feature variables are exhibited in Table S4, and downward solar radiation, the primary source of the energy at the earth surface, is the most important input variable, consistent with the results from our previous study (Wang et al., 2025). Although some of the selected variables may exhibit a certain degree of multi-collinearity, each contributes unique and physically meaningful information, supporting the inclusion of all variables in model construction. Note that the variable importance, derived from the built-in method of the random forests and potentially affected by multicollinearity among the input variables, is presented only as a reference. Retaining all 19 feature variables ensures the model's flexibility and generalization capability, enabling future incorporation of additional representative ground-based observations for further training and improvement. Besides, to investigate the impact of lagged effects of input variables on model performance, experiments were also conducted by adding lagged variables (e.g., the air temperature of the previous day) to the 19 input features. The results (Fig. S5) showed almost no improvement in model accuracy, suggesting that lagged effects on model performance were negligible within the CoSEB framework for estimates of daily surface radiation and heat fluxes. Furthermore, to better illustrate the effect of including additional radiation components (SW_IN, SW_OUT, LW_IN and LW_OUT) in the renewed CoSEB model compared with the original version by Wang et al. (2025), we have tested the performance of a reconstructed model that estimated only Rn, LE, H and G using the same independent variables and samples as those in the renewed CoSEB model. The results (Fig. S6) showed no significant differences in accuracy compared with those of the renewed CoSEB model, indicating the expansion of radiation components did not compromise model performance.

The main advantages of our CoSEB-based datasets of land surface radiation and heat fluxes lie in that [1] they are the first data-driven global datasets that satisfy both surface radiation balance (${\text{SW}}_{\text{IN}}-{\text{SW}}_{\text{OUT}}+{\text{LW}}_{\text{IN}}-{\text{LW}}_{\text{OUT}}=\text{Rn}$) and heat balance ($LE+H+G=\text{Rn}$) among the eight fluxes, as demonstrated by both the RIR and EIR of 0, [2] the radiation and heat fluxes are characterized by high accuracies when validated against in-situ measurements at 44 independent test sites (see the second paragraph in Sect. 2.1), where (1) the RMSEs for daily estimates of SW_IN, SW_OUT, LW_IN, LW_OUT, Rn, LE, H and G from the CoSEB-based datasets were 37.52, 14.20, 22.47, 13.78, 29.66, 30.87, 29.75 and 5.69 W m⁻², respectively, as well as for 8 d estimates were 18.54, 12.19, 18.50, 9.41, 19.12, 22.31, 21.63 and 4.60 W m⁻², respectively, (2) the CoSEB-based datasets, in comparison to the mainstream RS-based products/datasets (i.e. GLASS, BESS-Rad, FLUXCOM, BESSV2.0, MOD16A2, PML_V2 and ETMonitor), better agreed with the in situ observations at the 44 test sites, showing the RMSE reductions ranging from 0.01 to 4.58 W m⁻² for SW_IN, LW_IN, LW_OUT, Rn and LE at daily scale, and 0.24 to 10.48 W m⁻² for SW_IN, LW_IN, LW_OUT, Rn, LE and H at 8 d scale. Furthermore, the CoSEB-based datasets outperformed the ERA5-Land reanalysis datasets in estimating surface energy fluxes (where SW_OUT, LW_OUT, Rn and G for the ERA5-Land were inferred from surface radiation balance and heat balance), particularly for SW_OUT, H and G, with RMSE reductions of 0.13–8.15 W m⁻² when validated against in situ observations at the 44 test sites (Figs. S7 and S8). Preliminary analysis indicates that the CoSEB-based datasets exhibit spatial patterns consistent with those of mainstream RS-based datasets and Earth system model outputs (see Fig. S9), suggesting that the CoSEB-based datasets (or CoSEB framework) more broadly, are capable of reproducing the large-scale spatial features of Earth system models. This capability would be a great benefit to the community given the limitations associated with the high computational cost and long execution time of Earth system models. More detailed analysis about their similarities and differences can be further conducted in future work.

Our developed datasets could be potentially applied in many fields, including but not limited to (1) exploring the spatial-temporal patterns of global land surface radiation and heat flux (es) and their driving mechanisms over the past decades under global change (e.g., rising CO₂ concentration, greening land surface and increasing air temperature), (2) investigating the variability of land surface radiation and heat fluxes caused by extreme events and human activities, e.g. afforestation or deforestation, wildfire, air pollution, weather extremes and urbanization, (3) assessing the resources of solar energy, geothermal energy, surface and ground water at regional and global scales, (4) monitoring natural hazards, e.g. drought in agriculture and forestry.

The uncertainties of our datasets are relevant to (1) the data preprocessing, and (2) the application of the CoSEB model across different spatial scales. Specifically, daily averages of surface radiation and heat fluxes for each day were obtained for analysis from good-quality half-hourly observations when the fraction of these good-quality half-hourly observations was greater than 80 % in a day, due to the lack of consensus on the method for aggregating gapped half-hourly observations to daily data (Tang et al., 2024a; Yao et al., 2017; Zheng et al., 2022). Simple temporal interpolation of half-hourly in situ observations, which could therefore introduce substantial uncertainties, was not applied, because surface radiation and heat fluxes are sensitive to short-term variations in meteorological conditions and their intraday dynamics are often complex. Likewise, since there was no agreement on how to correct for the energy imbalance of turbulent heat fluxes, we adopted the most widely applied Bowen ratio method to enforce energy closure between Rn−G and LE+H (Castelli et al., 2018; Twine et al., 2000; Zhang et al., 2021). Another potential source of uncertainty arises from differences in meteorological reanalysis data caused by spatial downscaling, which, as demonstrated in our previous study (Wang et al., 2025, the last paragraph of Sect. 5.1), has a relatively small impact on model estimates by the machine-learning-based CoSEB model combined with finer-resolution surface-related variables that partially compensate for the spatial heterogeneity and localized variations not captured by the coarse-resolution datasets. These data preprocessing had an effect on the construction of the renewed CoSEB model, which may further affect the global datasets. Moreover, the renewed CoSEB model was constructed at the spatial scale of 500 m to match the footprints of the in situ EC observations, but applied at the spatial resolution of 0.05° to generate global datasets, mainly limited by the computing and storage capabilities of our personal computers. However, the CoSEB-based datasets have also been validated and inter-compared at 44 independent test sites to demonstrate that the difference in spatial scale would not much affect the performance of the datasets. Note that the 302 sites used for training, validation, and testing are predominantly located in the Northern Hemisphere, reflecting the inherent uneven distribution of the global flux networks. Although these sites cover a wide range of land cover types and climate regimes, thereby providing substantial heterogeneity for model development, the limited representation of the Southern Hemisphere may introduce uncertainties in the estimation of surface radiation and heat fluxes for certain ecosystems and soil types. In the future, enhancing the flux observation network coverage in the Southern Hemisphere, particularly in South America and Africa, and incorporating these observations into the CoSEB framework would help further improve the accuracy of surface radiation and heat flux estimates in these regions. Furthermore, the radiation and heat balance in this study refers specifically to the conservation among the eight variables (i.e., SW_IN, LW_IN, SW_OUT, LW_OUT, Rn, LE, G, H), which constitute the major components of the surface energy budget, and does not account for energy introduced by disturbance-related processes such as wildfires and volcanic eruptions. Despite these uncertainties, it is worth emphasizing that our work was the first attempt to innovatively develop data-driven energy-conservation datasets of global land surface radiation and heat fluxes with high accuracies.

6 Data availability

The energy-conservation datasets of global land surface radiation and heat fluxes generated by the CoSEB model with spatial-temporal resolutions of daily and 0.05° from 26 February 2000 to 31 December 2020 are freely available through the National Tibetan Plateau Data Center at https://doi.org/10.11888/Terre.tpdc.302559 (Tang et al., 2025a) and through the Science Data Bank (ScienceDB) at https://doi.org/10.57760/sciencedb.27228 (Tang et al., 2025b).

7 Summary and Conclusion

This study for the first time developed data-driven energy-conservation datasets of global land surface radiation and heat fluxes using our CoSEB model renewed based on GLASS and MODIS products, ERA5-Land reanalysis datasets, topographic data, CO₂ concentration data, and observations at 258 EC sites worldwide.

The CoSEB-based datasets of land surface radiation and heat fluxes are the first data-driven global datasets that satisfy both surface radiation balance (${\text{SW}}_{\text{IN}}-{\text{SW}}_{\text{OUT}}+{\text{LW}}_{\text{IN}}-{\text{LW}}_{\text{OUT}}=\text{Rn}$) and heat balance ($LE+H+G=\text{Rn}$) among the eight fluxes. Meanwhile, the CoSEB-based datasets outperformed the mainstream products/datasets in accuracy. Specifically, at 44 independent test sites, the RMSEs (R²) for daily estimates of SW_IN, SW_OUT, LW_IN, LW_OUT, Rn, LE, H and G from the CoSEB-based datasets were 37.52 W m⁻² (0.81), 14.20 W m⁻² (0.42), 22.47 W m⁻² (0.90), 13.78 W m⁻² (0.95), 29.66 W m⁻² (0.77), 30.87 W m⁻² (0.60), 29.75 W m⁻² (0.44) and 5.69 W m⁻² (0.44), respectively, as well as for 8 d estimates were 18.54 W m⁻² (0.87), 12.19 W m⁻² (0.39), 18.50 W m⁻² (0.92), 9.41 W m⁻² (0.97), 19.12 W m⁻² (0.82), 22.31 W m⁻² (0.67), 21.63 W m⁻² (0.39) and 4.60 W m⁻² (0.47), respectively. Moreover, the estimates from the CoSEB-based datasets in comparison to those from the mainstream products/datasets reduced the RMSE by 0.01 to 4.58 W m⁻² and increased the R² by 0.01 to 0.09 for SW_IN, LW_IN, LW_OUT, Rn and LE at daily scale, and reduced the RMSE by 0.24 to 10.48 W m⁻² and increased the R² by 0.01 to 0.38 for SW_IN, LW_IN, LW_OUT, Rn, LE and H at 8 d scale, when these estimates were validated against in situ observations at 44 independent test sites. Furthermore, the CoSEB-based datasets effectively captured the spatial-temporal variability of global land surface radiation and heat fluxes, aligning well with those from the mainstream products.

Our developed datasets hold significant potential for application across diverse fields such as agriculture, forestry, hydrology, meteorology, ecology, and environmental science. They can facilitate comprehensive studies on the variability, impacts, responses, adaptation strategies, and mitigation measures of global and regional land surface radiation and heat fluxes under the influences of climate change and human activities. These datasets will provide valuable insights and data support for scientific research, policy-making, and environmental management, advancing global solutions to address climate change.