Global hourly, 5 km, all-sky land surface temperature data from 2011 to 2021 based on integrating geostationary and polar-orbiting satellite data

Land surface temperature (LST) plays a dominant role in the surface energy budget (SEB) and hydrological cycling. Thermal infrared (TIR) remote sensing is the primary method of estimating LST globally. However, cloud cover leaves numerous data gaps in satellite LST products, which seriously restricts their applications. Efforts have been made to produce gap-free LST products from polar-orbiting satellites (e.g., Terra and Aqua); however, satellite data from limited overpasses are not suitable for characterizing the diurnal temperature cycle (DTC), which is directly related to heat waves, plant water stress, and soil moisture. Considering the high temporal variability in LST and the importance of the DTC, we refined the SEB-based cloudy-sky LST recovery method by improving its feasibility and efficiency and produced a global hourly, 5 km, all-sky land surface temperature (GHA-LST) dataset from 2011 to 2021. The GHA-LST product was generated using TIR LST products from geostationary and polar-orbiting satellite data from the Copernicus Global Land Service (CGLS) and the Moderate Resolution Imaging Spectroradiometer (MODIS). Based on ground measurements at the 201 global sites from the Surface Radiation Budget (SURFRAD), Baseline Surface Radiation Network (BSRN), Fluxnet, AmeriFlux, Heihe River basin (HRB), and Tibetan Plateau (TP) networks, the overall root-mean-square error (RMSE) of the hourly GHA-LST product was 3.31 K, with a bias of −0.57 K and R² of 0.95. Thus, this product was more accurate than the clear-sky CGLS and MODIS MYD21C1 LST samples. The RMSE value of the daily mean LST was 1.76 K. Validation results at individual sites indicate that the GHA-LST dataset has relatively larger RMSEs for high-elevation regions, which can be attributed to high surface heterogeneity and input data uncertainty. Temporal and spatial analyses suggested that GHA-LST has satisfactory spatiotemporal continuity and reasonable variation and matches the reference data well at hourly and daily scales. Furthermore, the regional comparison of GHA-LST with other gap-free hourly datasets (ERA5 and Global Land Data Assimilation System, GLDAS) demonstrated that GHA-LST can provide more spatial texture information. The monthly anomaly analysis suggests that GHA-LST couples well with global surface air temperature datasets and other LST datasets at daily mean and minimum temperature scales, whereas the maximum temperature and diurnal temperature range of LST and air temperature (AT) have different anomalous magnitudes. The GHA-LST dataset is the first global gap-free LST dataset at an hourly, 5 km scale with high accuracy, and it can be used to estimate global evapotranspiration, monitor extreme weather, and advance meteorological forecasting models. GHA-LST is freely available at https://doi.org/10.5281/zenodo.7487284 (Jia et al., 2022b) and http://glass.umd.edu/allsky_LST/GHA-LST (last access: 10 February 2023; Jia et al., 2022c).

1 Introduction

Land surface temperature (LST) is an essential component of the surface radiation budget and a dominant driving force in atmospheric cycling and hydrological balance (Z.-L. Li et al., 2022, 2013). LST directly reflects the thermal feedback of various land covers towards incoming solar radiation and atmospheric longwave radiation (Liang et al., 2019), and it is employed as an important variable for urban heat island analysis (Liu et al., 2022), permafrost mapping (Zou et al., 2017), and hazard forecasting (Bhardwaj et al., 2017; Mudele et al., 2020; Quintano et al., 2015). Therefore, LST has been extensively utilized as a vital indicator for characterizing regional and global climate change (Zhou et al., 2012; Jin, 2004; Peng et al., 2014). This parameter can be obtained by ground measurements, model simulations, and remote sensing retrievals. However, given the high spatiotemporal heterogeneity caused by various land covers, soil types, topographies, and meteorological conditions (Zhan et al., 2013; Liu et al., 2006; Ma et al., 2021), remote sensing has become the only feasible solution for monitoring LST globally.

LST can be retrieved using thermal infrared (TIR) observations from both polar-orbiting (Wan, 2008; Hulley and Hook, 2009) and geostationary (GEO) satellites (Yu et al., 2008; Freitas et al., 2009). In comparison, GEO satellites provide sub-hourly observations; thus, they can precisely capture the diurnal temperature cycles (DTCs) of the land surface. DTCs characterize the strong temporal variability in LST in a day, which is an important surface thermal property that responds to local environmental changes (Hansen et al., 1995; Sun et al., 2006). Studies have suggested that DTCs are directly related to plant water stress and soil drought (Fensholt et al., 2011; Stisen et al., 2008; Hernandez-Barrera et al., 2017); thus, such a relationship has been utilized for mapping evapotranspiration (Anderson et al., 2011) and soil moisture (Piles et al., 2016). In addition, it helps improve meteorological forecasting through data assimilation (Orth et al., 2017), extreme heat wave assessments (Hrisko et al., 2020; Jiang et al., 2015), crop yield estimations (Anderson et al., 2016), LST spatiotemporal-scale conversions (Hu et al., 2020), orbit drift corrections of Advanced Very High Resolution Radiometer (AVHRR) LST data (Jin and Treadon, 2003), and vegetation phenology analyses (Piao et al., 2015). Considering the great potential of DTCs in scientific applications and the high temporal variability in LSTs, accurate diurnal LST datasets are crucial for the research community and public (Chang et al., 2021; Hrisko et al., 2020; Pinker et al., 2019).

TIR sensors on board GEO satellites, such as the Geostationary Operational Environmental Satellites (GOES)-R Advanced Baseline Imager (ABI) and the Meteosat Second Generation (MSG) Spinning Enhanced Visible and InfraRed Imager (SEVIRI), provide exceptional opportunities to record DTCs. However, two inevitable flaws occur when using GEO satellites to observe diurnal LST variations globally, namely data gaps caused by cloud cover and the limited spatial view fields of individual GEO satellites, which seriously limit the availability of hourly all-sky LST datasets at the global scale. Methods of recovering hourly LSTs have been developed and comprehensively reviewed by Jia et al. (2022a). Currently available gap-free satellite-derived LST products are summarized in Table 1. Some gap-free LST datasets are not listed in the table, such as skin temperature from reanalysis datasets (e.g., ERA5 and the Modern-Era Retrospective Analysis for Research and Applications, version 2 – MERRA-2) (Muñoz-Sabater et al., 2021; Molod et al., 2015) or the results of Coccia et al. (2015), as they assumed that the surface broadband emissivity was equal to 1.

Table 1 Summary of publicly available gap-free LST products.

Table 1 reveals that only a few hourly all-sky LST datasets are currently available; thus, an all-sky hourly LST dataset at the global scale is urgently required. In Table 1, the products were divided into three categories based on their associated methodology: data fusion, mathematical interpolation, and cloudy-sky LST estimation based on the surface energy balance (SEB) theory.

Land surface models and reanalysis datasets release simulated hourly skin temperatures continuously, which have been fused with satellite-retrieved LSTs to generate gap-free LSTs (Dumitrescu et al., 2020; Long et al., 2020; Marullo et al., 2014; Ma et al., 2022; Muñoz-Sabater et al., 2021); however, the accuracy of the recovered LSTs is highly dependent on simulation accuracy, especially for continuous cloudy days. In addition, passive microwave (PMW) observations can penetrate clouds and estimate LSTs under all-sky conditions (Zhang et al., 2019; Wu et al., 2022), and studies have explored fusing such data with TIR LSTs from sensors of polar-orbiting satellites (Zhang et al., 2020; Xu and Cheng, 2021). However, PMW data have limited overpass times in a day; thus, they cannot match well with GEO observations. Mathematical interpolation is a popular method of reconstructing hourly LST because an ideal DTC can be parameterized by a harmonic function in the daytime and an exponential function in the nighttime (Duan et al., 2012). However, such parameterization requires at least four observations per day. Researchers have also tried to improve the feasibility of obtaining gap-free LSTs by combining DTC models with spatial interpolation (Liu et al., 2017) or utilizing a convolutional neural network (CNN) to predict missing values from neighboring clear-sky pixels and texture information (Wu et al., 2019). However, interpolating adjacent clear-sky samples can only obtain theoretical “clear-sky” LSTs because actual LSTs under clouds are impacted by frequent meteorological changes and cloud cooling/warming effects during the daytime/nighttime (Jin, 2000; Jia et al., 2020).

In comparison, cloudy-sky LST estimates based on SEB exhibit advantages in generating all-sky diurnal LST products at large scales. SEB-based methods include two steps: the first step is to reconstruct theoretical clear-sky LST values for cloudy time points, and the second step is to superpose the cloud effect based on the SEB equation (Jin and Dickinson, 2000; Lu et al., 2011). However, traditional SEB-based methods have limited feasibility because of the high input requirements, and they can only be used during the daytime. Therefore, Jia et al. (2021) considerably improved upon these methods by incorporating modeling data into the process, and the improved methods can be applied to larger spatial scales. ERA5 surface longwave radiation data were used to build a spatiotemporally evolving model, and the clear-sky GEO LSTs were assimilated to the evolving model to correct its predictions on cloudy days. Moreover, an optimization method was used to compute the cloud effect during both the daytime and nighttime, and complete DTCs could be recovered for hourly LSTs (Jia et al., 2022a).

The latest SEB scheme has produced all-sky hourly LSTs over the contiguous United States (CONUS) and Mexico from ABI data; however, it has a relatively lower computational efficiency owing to the spatiotemporal assimilation framework, which is not easy to apply globally. Furthermore, single GEO data points have a limited spatial view field, which can be solved by combining data from multiple GEO satellites at middle and low latitudes and from polar-orbiting satellites (Terra and Aqua) at high latitudes. Two polar-orbiting satellites pass high latitudes at a sub-hourly temporal scale and provide observations as frequently as GEO satellites. This strategy has been successfully utilized to generate global hourly Clouds and Earth's Radiant Energy System (CERES) radiation products (Loeb et al., 2018); however, few studies have focused on estimating all-sky LSTs by combining polar-orbiting and GEO satellites.

In this study, we produced a global hourly, all-sky LST dataset (GHA-LST) from 2011 to 2021 at a 5 km scale, and a comprehensive assessment was implemented using 201 ground sites worldwide. Global clear-sky LSTs were obtained by combining GEO LSTs from the Copernicus Global Land Service (CGLS) and Moderate Resolution Imaging Spectroradiometer (MODIS) MxD21 LST swath products, and a more efficient spatiotemporal assimilation scheme was proposed. It represents the first available global all-sky LST scheme on an hourly timescale with satisfactory accuracy based on global site validation; thus, it has great potential for use in analyzing global thermal dynamics, atmospheric cycling, and hydrological budgets.

2 Data and method

2.1 Data

The proposed GHA-LST dataset was recovered from a combination of clear-sky LST products, including the CGLS hourly LSTs, which cover middle and low latitudes, and MOD and MYD21 instantaneous swath LSTs, which cover high latitudes. ERA5 provides dynamic surface temperature signals for building a temperature-time-evolving model, and CERES global hourly surface radiation products were used to compute the cloud effect. In addition, the all-sky LST data were comprehensively assessed based on globally distributed sites collected from the Surface Radiation Budget (SURFRAD), Baseline Surface Radiation Network (BSRN), Fluxnet, AmeriFlux, Heihe River basin (HRB), and Tibetan Plateau (TP) networks.

2.1.1 Input data

CGLS LST provides hourly clear-sky LST retrievals from a constellation of GEO satellites, including multiple generations of Meteosat Second Generation (MSG), Multifunctional Transport Satellite (MTSAT)/Himawari, MSG Indian Ocean Data Coverage (IODC), and GOES. The product is released as a global product that covers land surfaces worldwide within the 60^∘ S to 70^∘ N latitudes. The generalized split-window (GSW) algorithm and dual algorithm (DA) in mono- and dual-channel forms were used to retrieve LST from top-of-atmosphere (TOA) brightness temperatures (BTs) in thermal infrared window channels (Freitas et al., 2013). Based on the ground validation, the accuracy ranges from 1.83 to 3.70 K.

MOD and MYD21 swath instantaneous LST products (Hulley et al., 2016) were used to provide LST over the rest of the land surface, which mainly covered high latitudes. A temperature/emissivity separation (TES) algorithm was used to retrieve the LST in the MOD and MYD21 products. It showed comparable accuracy to that of MOD and MYD11 (Wan, 2008) for most land cover types and performed better in bare land regions (Li et al., 2020; Yao et al., 2020). Level 3 MOD and MYD21 LST products provide gridded LST data in sinusoidal projection such that LST images were available four times a day and pixel locations were fixed. This data format is convenient for users; however, many valid retrieval values are lost at higher latitudes due to reprojection. In fact, dozens of times can be recorded by combining two polar-orbiting satellites, which is comparable to GEO observations at middle and low latitudes. Therefore, to fully incorporate the available clear-sky retrievals, MOD and MYD21 swath instantaneous LST data were used in this study. All swath images were converted to the Climate Modeling Grid (CMG) individually and then aggregated to the same spatial resolution as CGLS LST. The instantaneous observations were used for averaging only when they were within a 30 min window centered on the CGLS recording time (UTC standard time). In addition, to minimize the impact of retrieval uncertainty, records with a view zenith angle greater than 40^∘ were not used in this study (Li et al., 2014; Guillevic et al., 2013). By using this strategy, we can obtain a recording frequency over polar regions that is comparable to that of the CGLS LST data; however, using this strategy does not mean that significantly more clear-sky LST samples will be obtained, as cloud cover persists at high latitudes (King et al., 2013). This process only ensures that clear-sky LSTs at high latitudes are included in as many observations as possible.

To obtain continuous surface thermal variational signals, surface upward longwave radiation (ULW) and downward longwave radiation (DLW) simulated by ERA5 were used to build the LST time-evolving model. Satellite-derived broadband emissivity (BBE) was obtained from the Global LAnd Surface Satellite (GLASS) (Liang et al., 2021). We calculated the LST series using ERA5 DLW and ULW data simulated using clear-sky scenarios, which were generated based on real atmospheric and meteorological conditions, although clouds were assumed to be absent. The ERA5 clear-sky scenario was used because utilizing cloud radiative forcing calculated from global satellite data is more accurate than using the simulated results from the reanalysis (Wang and Dickinson, 2013). The clear-sky LST retrievals were then assimilated into the time-evolving model to obtain continuous LSTs without cloud gaps, and the cloud cooling/warming effect was then estimated and superposed from satellite radiation products. Essentially, in this revised SEB-based recovery method, the temperature change signals under cloud cover were divided into two parts: the evolving model provided LST variations due to real-time meteorological changes under clouds, and satellite radiation products estimated the cloud cooling and warming effects caused by cloud radiative forcing.

Global hourly surface DLW and downward shortwave radiation (DSR) from CERES satellite products were used to estimate the cloud effect. To monitor cloud radiative forcing, the CERES project retrieved global gap-free, hourly DSR and DLW under both all-sky (realistic) and theoretically cloud-free conditions (Doelling et al., 2016). CERES utilized the same strategy as this study to generate global hourly radiation products by combining remote sensing observations from multiple GEO sensors and two MODIS sensors. The CERES surface shortwave radiation and longwave radiation (Doelling et al., 2013) were estimated based on the Langley Fu–Liou radiative transfer theory (Fu et al., 1997), the cloud properties were obtained from microwave cloud products (Minnis et al., 2020), and the aerosol optical depth was based on the MODIS aerosol product (Remer et al., 2006). Surface CERES downward radiation fluxes have an overall bias (standard deviation) of 3.0 W m⁻² (5.7 %) for shortwave and −4.0 W m⁻² (2.9 %) for longwave radiation, which have been validated based on 85 sites (Rutan et al., 2015). CERES has been extensively evaluated and is generally considered a benchmark for satellite radiation products for assessments and intercomparisons (Jia et al., 2018; R. Li et al., 2022; Wang and Dickinson, 2013).

To improve the production efficiency, the complicated downward longwave parameterization schemes in Jia et al. (2022a) were replaced by directly exploiting the CERES dataset and converting its cloud radiative forcing into the corresponding cloud cooling/warming effect. Specifically, the CERES DSR difference between all-sky and clear-sky schemes was considered cloud DSR forcing; combined with the GLASS surface albedo data, the cloud net shortwave forcing was then computed. Cloud DLW forcing represents the difference between CERES all-sky and clear-sky DLW products, and the corresponding net longwave forcing was estimated using an optimization method (Sect. 2.2.4). In addition, according to previous studies (Wang and Dickinson, 2013; Zhang et al., 2015), the impact of the coarse spatial resolution of CERES downward radiation can be ignored because it has less heterogeneity than surface variables. CERES products were bilinearly interpolated to match the spatial scale of the CGLS LST. However, this assumption may introduce a certainty degree of uncertainty in areas with rugged terrain because complicated terrain in a coarse pixel may still affect the downward radiation components and increase the heterogeneity.

The GLASS 0.05^∘ land surface albedo and BBE were used for the net radiation calculation (Liang et al., 2021), and the GLASS leaf area index (LAI) was used for computing ground heat flux from the net radiation. All input data were preprocessed using bilinear resampling to match the CGLS LST. The input metadata are listed in Table 2.

Table 2 Metadata input for production.

2.1.2 Ground measurement

To comprehensively assess the accuracy of the proposed GHA-LST dataset, globally distributed in situ sites must be collected for ground validation. We processed the records from the SURFRAD, BSRN, Fluxnet, AmeriFlux, HRB, and TP networks. SURFRAD was established in 1993 and consistently provides long-term ground measurements of the surface radiation components over the CONUS for climate research and remote sensing retrieval validation (Augustine et al., 2000). The BSRN is a combined network of globally distributed sites from several projects (Driemel et al., 2018), and it provides records with strict data quality maintenance; thus, it is usually used as a reference dataset for radiation product validation at the global scale. Fluxnet includes hundreds of ground sites that have been utilized for global LST validation and analysis (Xing et al., 2021). AmeriFlux measures radiation and carbon fluxes over South and North America (Novick et al., 2018). The HRB network is from the Chinese Heihe Watershed Allied Telemetry Experimental Research (HiWATER) project (X. Li et al., 2013), and the TP network is from the Third Tibetan Plateau Atmospheric Scientific Experiment (TIPEX-III) (Zhao et al., 2018); these two datasets have been used for LST validation at the kilometer scale (Xing et al., 2021).

In addition, only raw observations marked as “good quality” were used for validation. Site-measured ULW and DLW were used to compute LSTs with a GLASS BBE based on the Stefan–Boltzmann law. LST raw records within a 30 min window centered on each UTC time period were then aggregated for hourly LST validation. The daily mean LST was further aggregated as long as 24 hourly LSTs were available in a day, which is sufficient to represent the entire DTC based on Jia et al. (2022a). The validation period was 2011–2020. The following validation metrics are used in this study: N is the sample amount; bias, also called mean bias error (MBE), represents the systematic errors/differences between LST products and ground measurements; root-mean-square error (RMSE) characterizes the actual uncertainty caused by bias and random error; and R² indicates the overall goodness of fit based on a 1:1 line. These metrics are commonly used for LST validation. The standard deviation (SD) of the differences between LST products and site measurements was not used because it provides similar information to the RMSE but cannot reflect errors caused by systematic bias; thus, the SD is generally smaller than the RMSE.

The proposed GHA-LST dataset has a spatial resolution of approximately 5 km, although some sites may not be representative of the corresponding pixels. To remove sites with higher heterogeneity, we utilized two removal strategies. The first is based on 30 m LSTs from the United States Geological Survey (USGS) Landsat 8 Level 2 Collection 2, in which clear-sky Landsat LSTs were extracted from all site locations from 2013 to 2020 and the average 30 m LSTs were extracted within the corresponding 5 km pixel range; then, the RMSE at each site was computed as the site representativeness using the 30 m LSTs paired with the averaged 5 km LSTs. One site was marked if it had a considerably larger RMSE, indicating that there were larger LST differences between the 30 m and 5 km scales. The second strategy considered the MYD21C1 0.05^∘ LST as a benchmark LST product at a 0.05^∘ spatial scale. As MYD21 has been comprehensively validated and produces results with high and stable accuracy at a global scale (Li et al., 2020; Yao et al., 2020; Hulley, 2015), we argue that if one site has a significantly larger RMSE in the validation of MYD21C1 samples with good quality, the site will have low representativeness at the 5 km spatial scale. Preprocessed sites detected by either of these two strategies were excluded from this study, and the analysis results are shown in Fig. 1. The selection threshold of each strategy was equal to the average RMSE + 2SD for all sites.

Figure 1 LST site selection results based on the (a) site representativeness calculated by Landsat samples and (b) site RMSEs of MYD21C1. The thresholds are the black lines, and sites with RMSEs higher than the line were masked out. The abbreviations used in the figure are as follows: SURFRAD – Surface Radiation Budget, BSRN – Baseline Surface Radiation Network, HRB – Heihe River basin (HRB), and TP – Tibetan Plateau.

Based on Fig. 1, 201 global sites were included in this study, including 5 SURFRAD sites, 11 BSRN sites, 91 Fluxnet sites, 89 AmeriFlux sites, 3 HRB sites, and 2 TP sites. The distribution of sites is shown in Fig. 2. Additionally, as site selection influences the final validation statistics, we also validated samples from clear-sky MYD21C1 and CGLS for accuracy comparisons, and the accuracy level of the validation results can be used as the reference for GHA-LST.

Figure 2 Global distribution of the 201 LST sites utilized. The abbreviations used in the figure are as follows: SURFRAD – Surface Radiation Budget, BSRN – Baseline Surface Radiation Network, HRB – Heihe River basin (HRB), and TP – Tibetan Plateau.

2.2 Methods

2.2.1 Production framework

Jia et al. (2022a) developed a three-step framework to generate all-sky hourly LSTs on a regional scale. In this study, we revised this framework to improve its efficiency and feasibility worldwide. A flowchart of the framework is shown in Fig. 3.

Figure 3 Flowchart of the proposed GHA-LST dataset production. The abbreviations used in the figure are as follows: DLW – downward longwave radiation, ULW – upward longwave radiation, DSR – downward shortwave radiation, BBE – broadband emissivity, LST – land surface temperature, and CRE – cloud radiative effect.

In the first step, a time-evolving model of clear-sky LSTs was designed based on the ERA5 LST series at each pixel location (Sect. 2.2.2). The ERA5 LST series was computed from the DLW and ULW from the clear-sky simulation scenarios, and it provides continuous variational information on LST without considering the cloud cooling/warming effect. Such variational information under clouds can be attributed to advective meteorological changes and air movement. The ERA5 skin temperature was not involved because we calculated the cloud effect based on satellite-derived radiation products, which are more reliable at a global scale.

In the second step, the Kalman filter (KF) was used to assimilate available clear-sky LST retrievals into the time-evolving model to correct the predictions for times with cloud cover, and theoretical clear-sky LSTs were then reconstructed. In the original framework of Jia et al. (2022a), three-dimensional data assimilation was utilized to generate a spatiotemporally dynamic model; however, this process is time-consuming, particularly when working at a large spatial scale. We replaced the spatial module in the assimilation by linear regression, which still works well to incorporate spatially adjacent clear-sky retrievals (Sect. 2.2.3). After this step, hypothetical LSTs were reconstructed during times with cloud cover. The cloud effect was further superimposed in the final step.

In the third step, the cloud effect was estimated from the satellite radiation products (Sect. 2.2.4). In the original framework, cloud longwave radiative forcing was computed based on a series of parameterization schemes; however, the scheme was not well assessed at the global scale. Therefore, to simplify the calculation and improve the feasibility of filling gaps over a large spatial scale, we replaced the DLW parameterization with CERES clear-sky and all-sky DLW products that have been assessed globally (Wang and Dickinson, 2013). The cloud effects at daytime and nighttime were determined by searching for the optimal cloud radiative effect (CRE) values to meet the SEB. The final clear-sky LSTs are the assimilated results, and the cloudy-sky LSTs are the reconstructed LSTs from the second step plus the optimal cloud effect.

2.2.2 Time-evolving model

A time-evolving model describes how LSTs change at a certain pixel over time, and it characterizes relative variation based on the ERA5 LST rather than absolute magnitudes. The ERA5 LST series was initially downscaled to match the CGLS LST using elevation (Duan et al., 2017). The evolving model can be mathematically represented by Eqs. (1) and (2):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\text{LST}}_{t,d}=F_{t,d}\times {\text{LST}}_{t,d-\mathrm{1}},\text{(2)}& {\displaystyle }{F}_{t,d}=\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{Z_{t,d}+\mathit{\delta }}}\right){\displaystyle \frac{\mathrm{d}{Z}_{t,d}}{\mathrm{d}t}}.\end{array}$$

Here, LST_t,d is the LST predicted by the model on day d at hour t; F_t,d is the prediction operator, which is generated based on the temperature temporal profile Z_t,d (temperature difference between d and d−1 at hour t, i.e., the difference with the LST 24 h before); and δ=0.01 avoids a null denominator. The model evolves from day to day for each hour of the day (HOD) because the modeling bias is self-correlated at the same HOD on different days (Marullo et al., 2014). The correction of the data assimilation can be better inherited based on the evolving structure of Eqs. (1) and (2). In addition, only the difference information was used in the study, which can minimize the impact of the uncertainty of the ERA5 LST, especially of its systemic bias (Nogueira et al., 2021).

The time-evolving model provides continuous temperature variation; however, the modeling process must be consistently corrected by assimilating available clear-sky retrievals. In addition, partially cloud-covered yet retrieved (“likely cloud-contaminated”) pixels were masked out before assimilation. The detection method follows Jia et al. (2022a). One clear-sky retrieval is excluded if it has a substantially larger absolute difference from the modeled LST (3 SDs) than neighboring days within ±15 d, which assumes that modeled LSTs have fewer anomalies than directly retrieved values. It should be noted that some input data (e.g., CERES and reanalysis data) are not available at near-real time (NRT); moreover, this likely cloud-contaminated detection method also requires a 30 d time window for high-quality clear-sky LST selection, which means that the proposed cloudy-sky LST recovery method cannot be used for NRT all-sky LST production.

2.2.3 Kalman filter

The KF was used to assimilate clear-sky LST retrievals to correct the evolving model prediction because real-time retrievals are discontinuous, while the evolving model prediction is continuous. When a retrieval value is available, a weighted average is implemented between the prediction and the observation based on their individual uncertainties, and the prediction is then corrected. When observations are not available, the prediction will be implemented based on the updated prediction from the former step. Therefore, a continuous LST series can be generated using this iterative process. The KF can be mathematically represented as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\widehat{T}}_{t,d}^{-}=A_{d-\mathrm{1}}^t{\widehat{T}}_{t,d-\mathrm{1}}+{\mathit{\omega }}_{d-\mathrm{1}}^t,\text{(4)}& {\displaystyle }{\widehat{T}}_{t,d}={\widehat{T}}_{t,d}^{-}+K_d^t(T_{t,d}-{\widehat{T}}_{t,d}^{-}),\text{(5)}& {\displaystyle }{K}_{t,d}=P_{t,d}^{-}(P_{t,d}^{-}+R{)}^{-\mathrm{1}},\text{(6)}& {\displaystyle }{P}_{t,d}=(I-K_{t,d})P_{t,d}^{-}.\end{array}$$

Here, ${\widehat{T}}_{t,d}^{-}$ is the temperature prediction at time t on day d from the prediction of d−1, and $A_{d-\mathrm{1}}^t$ is the prediction process (Sect. 2.2.2) with a prediction error of ${\mathit{\omega }}_{d-\mathrm{1}}^t$. The symbol “−” next to a variable indicates that the variable is an initial prediction without assimilation correction. The modeling prediction is propagated to $P_{t,d}^{-}$ after this prediction.

If an observation (T_t,d) is available, ${\widehat{T}}_{t,d}^{-}$ will be corrected using the Kalman gain $K_d^t$ (Eq. 4), which was determined by the relative magnitude of the squared uncertainty of the model prediction $P_d^{-}$ and the satellite retrieval R. R is the squared retrieval uncertainty at each pixel of CGLS and MODIS, and it is calculated based on the “ERRORBAR_LST” or “LST_err” information in the files. The initial modeling uncertainty was calculated based on a comparison between the ERA5 and satellite retrievals at clear-sky time in the data series. The assimilation started in October 2010, and the model prediction reached a stable status before the product release date; thus, the initial value did not affect the output. The prediction error of d will also be corrected to P_t,d (Eq. 5). The next day will then be iteratively predicted. If there are no observations, the LST will be automatically predicted on day d without correction.

Jia et al. (2022a) used a spatial KF module that can assimilate spatially adjacent clear-sky pixels into the evolving model; however, this process is time-consuming and impractical for global production. In this study, if an observation on d was available, the time-evolving model was corrected by KF normally. If an observation was not available, a 30×30 (∼150 km, as per Jia et al., 2022a) spatial window was set for the time-evolving model, and clear-sky pixels and their corresponding ERA5 LST within the spatial window will regress a linear conversion model; the missing LST at the center pixel on d will then be predicted from its corresponding ERA5 LST using such a linear relationship. Essentially, a local linear relationship replaced the spatial KF module, although it still effectively incorporated the adjacent clear-sky retrieval, and the computation efficiency was significantly improved. After the center LST was estimated by linear regression, it was considered the available retrieval for KF correction on d. If the spatial window did not have available clear-sky retrievals, the time-evolving model predicted the LST on d based on the results on d−1.

2.2.4 Cloud effect estimation

After data assimilation, the LST under clouds was initially predicted without considering the cloud effect. The cloud effect is the temperature warming/cooling effect caused by changing the SEB from clouds, which can be estimated using cloud radiative forcing. The SEB equation is as follows:

$$\begin{array}{}\text{(7)}& R_{\mathrm{n}}=R_{\mathrm{s}}^d\left(\mathrm{1}-\mathit{\alpha }\right)+\mathit{\epsilon }{R}_l^d-\mathit{\sigma }\mathit{\epsilon }{T}^{\mathrm{4}}=G+\mathit{LE}+H,\end{array}$$

where R_n is the surface net radiation, $R_{\mathrm{s}}^d$ is the DSR, α is the surface albedo, ε is the surface broadband emissivity (BBE), $R_l^d$ is the DLW, σ is the Stefan–Boltzmann constant, and T is the LST. R_n is partitioned into latent heat (LE), sensible heat (H), and ground heat (G). Cloud cover changes R_n, which is called cloud radiative forcing. By following the land surface analysis (LSA SAF) GEO evapotranspiration product, G can be parameterized as follows:

$$\begin{array}{}\text{(8)}& G=R_{\mathrm{n}}\times \mathrm{0.5}\exp (-\mathrm{2.13}(\mathrm{0.88}-\mathrm{0.78}\exp (-\mathrm{0.6}\text{LAI})\left)\right),\end{array}$$

where the G is set to 0.15 (0.05) R_n for bare land (permanent snow/ice). Based on the conventional “force-restore” method (Jin and Dickinson, 2000), G can be represented as follows:

$$\begin{array}{}\text{(9)}& G=k_{\mathrm{g}}{\displaystyle \frac{\partial T}{\mathrm{\Delta }Z}}=k_{\mathrm{g}}{\displaystyle \frac{T-T_{\mathrm{d}}}{\mathrm{\Delta }Z}},\end{array}$$

where k_g represents surface thermal conductivity (W m⁻¹ K⁻¹) and ΔZ is the responding surface depth, which is set to 0.1 m. The deep-layer temperature (T_d) is assumed to have little response towards SEB; thus, Eq. (9) can be rewritten as follows:

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial G}{\partial T_{\mathrm{s}}}}={\displaystyle \frac{\partial }{\partial T_{\mathrm{s}}}}\left[k_{\mathrm{g}}{\displaystyle \frac{T_{\mathrm{s}}-T_{\mathrm{d}}}{\mathrm{\Delta }Z}}\right]\approx {\displaystyle \frac{k_{\mathrm{g}}}{\mathrm{\Delta }Z}}.\end{array}$$

Accordingly, the change in G (∂G) caused by cloud cover can be directly converted into the variation in LST, and ∂G is determined by partitioned cloud radiative forcing. That is, by knowing any two of the three variables (∂G, ΔT_s, and k_g), the other can be estimated. According to Jia et al. (2022a), k_g was predetermined based on a continuous temperature series from the assimilation results and corresponding radiation data:

$$\begin{array}{}\text{(11)}& k_{\mathrm{g}}=\mathrm{\Delta }h{\displaystyle \frac{\overline{G_{\mathrm{noon}}}-\overline{G_{\mathrm{sr}}}}{\overline{T_{\mathrm{noon}}}-\overline{T_{\mathrm{sr}}}}},\end{array}$$

where $\overline{G_{\mathrm{noon}}}$ ($\overline{T_{\mathrm{noon}}}$) and $\overline{G_{\mathrm{sr}}}$ ($\overline{T_{\mathrm{sr}}}$) are the monthly averaged ground heat (clear-sky LST) within ±15 d at noon and sunrise, respectively, which are considered because morning warming can mainly be attributed to the SEB. The continuous data series from the data assimilation step ensures sufficient sampling for the k_g calculation. Monthly means were used to minimize the disturbance of daily variation. Then, the cloud radiative forcing needs to be determined to estimate ΔT_s:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{R}_{\mathrm{c}}& =\left(\mathrm{1}-\mathit{\alpha }\right)(R_{\mathrm{s},\mathrm{cld}}^d-R_{\mathrm{s},\mathrm{clr}}^d)+\mathit{\epsilon }\left(R_{\mathrm{l},\mathrm{cld}}^d-\mathit{\sigma }{T}_{\mathrm{cld}}^{\mathrm{4}}\right)\\ & -\mathit{\epsilon }\left(R_{\mathrm{l},\mathrm{clr}}^d-\mathit{\sigma }{T}_{\mathrm{clr}}^{\mathrm{4}}\right),\end{array}\end{array}$$

where R_c is the cloud radiative forcing, $R_{\mathrm{s},\mathrm{cld}}^d$ ($R_{\mathrm{s},\mathrm{clr}}^d)$ is the cloudy-sky (clear-sky) DSR, and $R_{\mathrm{l},\mathrm{cld}}^d$ ($R_{\mathrm{l},\mathrm{clr}}^d$) is the cloudy-sky (clear-sky) DLW. The shortwave variables, BBE, and DLW can be obtained from the global radiation products.

The cloudy-sky LST (T_cld) equals the reconstructed LST (T_clr) plus the cloud effect (ΔT_s), and T_clr is reconstructed from step 2. ΔT_s is unknown for the radiative forcing calculation and represents the ultimate target of this step. Therefore, the optimal ΔT_s must be determined based on the optimization method to satisfy the SEB. Following Jia et al. (2022a), ΔT_s was initially assumed to be 0 K, and the initial R_c was obtained based on Eq. (12). After energy partitioning through the LAI, G is computed, and the updated ΔT_s is estimated using k_g; thus, R_c can be recomputed. By iteratively comparing the R_c differences and adjusting ΔT_s (step = 0.05 K), the surface energy budget is balanced (|ΔCRE$|<\mathrm{20}$ W m⁻²; see Fig. 3). The threshold of 20 W m⁻² is the current accuracy level of the longwave radiation products (Wang et al., 2020).

3 Results and discussion

3.1 Overall assessment

Based on all paired samples from the proposed GHA-LST dataset and 201 sites from 2011 to 2020, the overall RMSE of the all-sky GHA-LST is 3.31 K, with a bias of −0.57 K and R² of 0.95. As site selection may influence the accuracy statistics, the accuracy of GHA-LST was also compared with that of the CGLS and MYD21C1 data (Fig. 4). Figure 4a shows the extracted clear-sky samples from the CGLS, whereas the corresponding clear-sky results from GHA-LST are shown in Fig. 4b, which has the same sampling amount as Fig. 4a. The recovered cloudy-sky LST at the corresponding CGLS cloudy time of Fig. 4a is validated in Fig. 4c. Such an accuracy comparison between the GHA-LST and CGLS data is sufficiently fair; thus, we also compared the accuracies between GHA-LST and MYD21C1 (Fig. 4d, e, f). Because MYD21C1 may have a slightly different observation time (<0.5 h) relative to GHA-LST, it was converted to the nearest UTC time based on the diurnal cycle recorded by site observations to match the GHA-LST recording time.

Figure 4 Density scatterplots of hourly LST samples from (a) CGLS clear-sky retrievals, (b) GHA-LST clear-sky samples corresponding to panel (a), (c) GHA-LST cloudy-sky samples corresponding to CGLS cloudy time, (d) MYD21C1 clear-sky retrievals, (e) GHA-LST clear-sky samples corresponding to panel (d), and (f) GHA-LST cloudy-sky samples corresponding to MYD21C1 cloudy time.

The proposed GHA-LST dataset had better accuracy than the CGLS and MYD21C1 data under both clear-sky and cloudy-sky conditions (Fig. 4). Compared to the CGLS (Fig. 4a), the clear-sky GHA-LST had a similar RMSE (3.35 K; Fig. 4b), which is reasonable because most of the clear-sky GHA-LST samples were derived from clear-sky retrievals from the CGLS. The partially cloud-contaminated samples were marked during production (Sect. 2.2.2) and considered cloudy pixels. The results of MYD21C1 (Fig. 4d) only utilized samples that were marked as good quality and passed a cloud contamination test (Sect. 3.6 in Ma et al., 2020). In comparison, GHA-LST produced similar accuracy with high-quality MODIS samples (Fig. 4e) and represented a stable accuracy (RMSE $=\sim \mathrm{3.3}$ K with few biases) under both clear-sky and cloudy-sky conditions. Based on the sampling amount (Fig. 4a, b, c), the number of cloudy pixels was nearly 1.5 times higher than the number of clear-sky pixels, indicating the importance and necessity of the proposed GHA-LST dataset. As MYD21C1 (Fig. 4d, e, f) only includes samples from noon and midnight, it has significantly fewer samples than the CGLS group (Fig. 4a, b, c). Studies on similar topics, such as the recovery of all-sky MSG/SEVIRI LST (Martins et al., 2019), produced RMSEs of 2.1–3.7 K at three sites. However, direct comparisons of the validation statistics are difficult because substantially different sites and sampling amounts were utilized in this study.

Figure 5 illustrates that the GHA-LST process has a considerably higher performance in obtaining daily mean LSTs. GHA-LST has a similar RMSE to the daily mean CGLS, with 24 values in a day retrieved by satellites (Fig. 5a, b); however, GHA-LST has substantially more available samples than the CGLS clear-sky results. In comparison, the daily mean computed from the average of the paired daytime/nighttime MYD21C1 has the largest RMSE, with a clear bias of 2.41 K. However, MODIS LSTs have been more widely used than GEO LSTs, and many studies have obtained daily mean LSTs by simply averaging two instantaneous Aqua retrievals at noon and midnight (Ouyang et al., 2012; Chen et al., 2017; Zou et al., 2017). This study suggests that the proposed GHA-LST dataset can significantly improve the accuracy and data availability of the daily mean LST.

Figure 5 Density scatterplots of the daily mean LST samples from (a) GHA-LST, (b) CGLS, and (c) average daytime and nighttime MYD21C1 pairs.

To demonstrate the accuracy and stability of the GHA-LST under different surface conditions, the hourly samples were further differentiated based on land cover types. Land cover data are obtained from the MCD12Q1 International Geosphere-Biosphere Programme (IGBP) classification. The validation statistics are listed in Table 3.

Table 3 Validation statistics for different land cover types.

Table 3 indicates that the GHA-LST has stable accuracy under both clear-sky and cloudy-sky conditions for various land cover types. However, some clear biases were found for the forest and barren land cover types, which could be caused by split-window retrieval errors under clear-sky conditions due to the large emissivity uncertainty (Z.-L. Li et al., 2022). In comparison, high R² values in these regions reflect the ability of the GHA-LST to capture regional temperature variations.

3.2 Individual site validation

Considering that 201 global sites were utilized, the RMSEs at individual sites can reflect the spatial pattern of GHA-LST accuracy; therefore, site RMSE maps under clear-sky and cloudy-sky conditions are illustrated in Fig. 6.

Figure 6 Global RMSE statistics at individual sites under (a) clear-sky and (b) cloudy-sky conditions.

The GHA-LST had similar accuracy patterns under both clear-sky and cloudy-sky conditions, and it had lower RMSEs (<3.5 K) in eastern America, Europe, and Australia but variable RMSEs from 3.5 to 5 K in western America and the TP. The median site RMSE under clear-sky conditions was 3.18 K, with an SD of 0.81 K, and the median site RMSE under cloudy-sky conditions was 2.97 K, with an SD of 1.01 K, indicating that cloudy-sky results had a slightly larger spatial variance with respect to accuracy. Validation statistics for each site are listed in the Appendix.

3.3 Temporal and spatial analysis

To evaluate the temporal continuity of the proposed GHA-LST dataset, the temporal LST variations from GHA-LST, CGLS, and corresponding ground measurements were compared at hourly and daily mean scales. In Fig. 7, four global sites are shown as the representative sites from SURFRAD (SXF; 43.74^∘ N, 96.62^∘ W), BSRN (CAB; 51.97^∘ N, 4.93^∘ E), Fluxnet (AU-Rig; 36.65^∘ S, 145.58^∘ E), and AmeriFlux (US-Ro1; 44.71^∘ N, 93.09^∘ W) (see the figure caption for full site names). The study period was randomly chosen for different years.

Figure 7 Temporal variation in the (a–d) hourly LST and (e–h) daily mean LST at four sites: (a, e) SURFRAD SXF, (b, f) BSRN CAB, (c, g) Fluxnet AU-Rig, and (d, h) AmeriFlux US-Ro1. The abbreviations used in the figure are as follows: SURFRAD SXF – “Surface Radiation Budget Sioux Falls”; BSRN CAB – “Baseline Surface Radiation Network Cabauw”; AU-Rig – “Australia – Riggs Creek”; US-Ro1 – “United States-Ro1”; GHAT – global hourly, 5 km, all-sky land surface temperature; and CGLS – Copernicus Global Land Service.

The GHA-LST data have satisfactory temporal continuity and correspond to the ground measurements at hourly and daily mean scales. At the hourly scale, the hourly LST exhibits harmonic diurnal variations under clear-sky conditions, especially at Fluxnet AU-Rig (Fig. 7c), where the climate is dry and cloud cover is low. In comparison, a more complicated temporal pattern of LSTs is observed for continuous cloudy time (e.g., BSRN CAB; Fig. 7b), indicating that the harmonic function-based DTC models may not work well in these cases. The GHA-LST data can capture the DTCs under both clear-sky and cloudy-sky conditions and correspond to the ground measurements and clear-sky CGLS. Certain clear-sky CGLS points are scattered and have a clear negative bias (Fig. 7a) because they were detected as partially covered pixels; thus, they were not used in the data assimilation. At the AmeriFlux US-Ro1 site (Fig. 7d), GHA-LST is more consistent with CGLS than the ground measurements on clear-sky days; thus, we infer that US-Ro1 has a larger heterogeneity issue than the other sites.

After temporal aggregation, the daily mean LST variation in different years also demonstrated the continuity and stable accuracy of the GHA-LST. The relatively larger differences between the satellite datasets and ground measurements at noon (hourly scale) and during summer (daily mean scale) can be explained by site representativeness. A temporal variation analysis of accuracy (Jia et al., 2022a) suggests that ground measurements generally have the lowest representativeness at noon and that the RMSE statistics of hourly LST products can increase by more than 1 K from nighttime to noon. As solar radiation increases in the morning, LST has distinct warming responses over different land cover types in a pixel; thus, the spatial heterogeneity of the pixel is enhanced during daytime.

Figure 8 Global LST maps of the (a) GHA-LST annual mean in 2021 and the (b) CGLS, MYD, and MOD21 annual mean in 2021 as well as the corresponding hourly LST maps at 02:00 UTC on 7 September 2021 for (c) GHA-LST and (d) CGLS, MYD, and MOD21. The abbreviations used in the figure are as follows: GHA-LST is the global hourly, 5 km, all-sky land surface temperature, and CGLS is the Copernicus Global Land Service.

A global all-sky LST map was analyzed to demonstrate the spatial continuity of the proposed GHA-LST dataset. The global annual mean maps of all-sky LST from GHA-LST and clear-sky LST from CGLS+MYD/MOD21 were compared (Fig. 8). Overall, the GHA-LST data exhibit high spatial continuity across the globe at different timescales. The annual mean GHA-LST for 2021 (Fig. 8a) illustrates a reasonable spatial pattern. The annual LSTs from CGLS + MYD/MOD21 (Fig. 8b) present a systematic clear-sky bias (Ermida et al., 2019), especially in the connection regions of the CGLS and MODIS data (southwestern Canada and Siberia), where clear artificial lines are shown. In comparison, assimilating the clear-sky results to the time-evolving model can produce more spatially consistent LST maps.

At the hourly scale, the GHA-LST map (Fig. 8c) can also produce the reasonable spatial variation in LST. Compared with the clear-sky pixels in Fig. 8d, the cloudy-sky pixels of the GHA-LST data (Fig. 8c) were well recovered. It should be noted that the clear-sky map (Fig. 8b) had more spatial textures than the all-sky annual mean map (Fig. 8a) because clear-sky LSTs have higher spatial heterogeneity due to solar heating. Furthermore, the various numbers of available clear-sky days in different locations may cause spurious spatial textures (e.g., lines at the connection region of CGLS and MODIS). Additionally, the GHA-LST spatiotemporally filtered clear-sky satellite LSTs using the simulated model series, which may sacrifice spatial textures for data fusion.

To evaluate the ability of GHA-LST to capture spatial textures at regional scales, the GHA-LST annual means in the Alaska and TP regions are shown in Fig. 9. These two areas were selected because they are hot spot regions in terms of their response to climate change (Kuang and Jiao, 2016; Melvin et al., 2017), and GHA-LST was recovered in these two regions from MODIS and CGLS, respectively. The corresponding annual mean skin temperatures of ERA5 and the Global Land Data Assimilation System (GLDAS) were also included for comparison. ERA5 and GLDAS were employed because global hourly all-sky LST is currently only available from reanalysis datasets, and both are widely used in the relevant research (Muñoz-Sabater et al., 2021; Rodell et al., 2004). GHA-LST has a spatial pattern similar to that of the two reanalysis datasets but produces many more spatial details. GHA-LST has a spatial resolution of ∼5 km; therefore, it can provide more spatial texture information than the ERA5 (0.1^∘) and GLDAS (0.25^∘) data. GLDAS data have invalid pixels, mainly because GLDAS ignores all inland lakes.

Figure 9 Temporally averaged all-sky LST maps in 2021 from (a, b) GHA-LST, (c, d) ERA5, and (e, f) GLDAS (a, c, e) over the Tibetan Plateau and (b, d, f) Alaska. The abbreviations used in the figure are as follows: GHA-LST – global hourly, 5 km, all-sky land surface temperature and GLDAS – Global Land Data Assimilation System.

3.4 Global anomaly analysis

To justify the potential use of LST in climate-warming-related issues, the relationship between LST and surface air temperature (AT) has recently been discussed. However, previous studies have either focused on local site scales (Hachem et al., 2012; Good, 2016; Sohrabinia et al., 2015; Mutiibwa et al., 2015) or ignored the clear-sky LST bias (Good et al., 2017). In comparison, the proposed GHA-LST provides an exceptional opportunity to spatiotemporally upscale the all-sky hourly LST to allow for comparisons with AT climate datasets. Monthly anomaly variations in global LST and AT are shown in Fig. 10 by removing the seasonal cycle in the daily mean temperature (T_mean), daily minimum temperature (T_min), daily maximum temperature (T_max), and diurnal temperature range (DTR).

Berkeley Earth Surface Temperatures (BEST; Rohde et al., 2013); NASA Goddard Institute for Space Studies Surface Temperature Analysis, version 4 (GISTEMP v4; Lenssen et al., 2019); and Climatic Research Unit Temperature, version 4 (CRUTEM4; Osborn and Jones, 2014), were used to characterize the AT variation (Fig. 10a). Other all-sky LST datasets were averaged and shown to verify the GHA-LST anomaly, including two MODIS-derived gap-free results, as shown in Table 1 (Hong et al., 2022; Zhang et al., 2022), and the ERA5-Land reanalysis skin temperature. Only BEST and GHA-LST can provide T_max and T_min; thus, they were used in Fig. 10b, c, and d. The reference time period was 2015–2017, and the uncertainty shadow in Fig. 10a is the SD of the averaged LST/AT datasets.

Figure 10 Monthly anomaly variation in the globally geographically weighted LST and land surface air temperature (AT) at different scales: (a) mean temperature (T_mean), (b) maximum temperature (T_max), (c) minimum temperature (T_min), and (d) diurnal temperature range (DTR). The abbreviations used in the figure are as follows: GHA-LST – global hourly, 5 km, all-sky land surface temperature and BEST – Berkeley Earth Surface Temperatures.

The LST anomaly couples well with the global AT anomaly at the T_mean and T_min scales, and LST has a slightly larger amplitude than AT, whereas the T_max and DTR of the two variables can only match the anomalous direction and the magnitude is quite different. At the T_mean scale, anomalies of GHA-LST and other LSTs have very similar variations with the AT datasets, even though they have completely different data sources (Fig. 10a). These findings are consistent with the trend comparison between the ERA5-Land skin temperature and AT (Wang et al., 2022). Both datasets have limited uncertainty (shadowed areas), indicating that they can accurately characterize the land surface thermal dynamics. In comparison, LST had a slightly larger anomaly amplitude than AT, which was mainly caused by the daytime LST. Solar heating increases the temperature difference between LSTs and ATs over different surface land cover types (Good et al., 2017). Accordingly, T_max exhibited the largest difference (Fig. 10b), especially in years with large anomalies, thereby indicating that LST T_max had a stronger response to heat anomalies. T_min showed a higher correlation between the two variables (Fig. 10c). LST has a considerably stronger DTR disturbance than AT owing to the difference in T_max (Fig. 10d). We did not quantify the trend magnitude because the time span was only 11 years and the overall trend was affected by the value in 1 specific year. This analysis demonstrates the potential usefulness of GHA-LST in climate studies and global hourly AT estimates.

4 Discussion

The site RMSE statistics were compared with the corresponding site elevations and latitudes to detect the potential factors that impact the accuracy of the results (Fig. 11).

Figure 11 Scatterplots of the site RMSE against (a, b) elevation and (c, d) latitude under (a, c) clear-sky and (b, d) cloudy-sky conditions. The significant linear relationships (p values <0.01) are drawn (a, b).

The RMSE statistics at each site were mainly affected by the site elevation (Fig. 11). The scatterplots of site RMSE and site elevation under clear-sky and cloudy-sky conditions (Fig. 7a, b) show that the linear relationship was statistically significant (p value <0.01). We suppose that increasing elevation will decrease the spatial representativeness of the sites; therefore, the RMSEs of the clear-sky results showed an increasing trend. In addition, the cloudy-sky results had a larger slope (Fig. 11b), indicating that they were more sensitive to elevation variations. Thus, we inferred that elevation was an essential factor affecting LST recovery. In these regions with high elevation, clear-sky LSTs with larger RMSEs were assimilated in the time-evolving model, thereby affecting the cloudy-sky results. In addition, modeled temperature series include higher uncertainty in these regions, and a relevant “cool bias” issue in highlands was found in regional simulation models and global reanalysis datasets (Jia et al., 2022d; Meng et al., 2018). Although the relationship between the site RMSE and site latitude was not statistically significant (Fig. 11c, d), the GHA-LST data at higher latitudes produced higher RMSEs than the data at lower latitudes, especially under cloudy-sky conditions; thus, we inferred that high latitudes were frequently covered by clouds and fewer clear-sky LSTs could be used in the data assimilation. In addition, sites at higher latitudes are usually located in coastal areas (Fig. 6a), which may limit their spatial representativeness at the 5 km scale.

Furthermore, the spatial continuity at the regional scale was evaluated. A detailed mapping examination suggests that no artificial textures occurred under most conditions at middle and low latitudes; however, at high latitudes where MODIS swath data are the basic input data, swath edges are observed on the map in some cases (Fig. 12).

Figure 12 Case showing the artificial texture (middle of the image) in northwestern Canada at 05:00 UTC on 2 January 2012.

Such spatial discontinuity occurs when the clear-sky LST retrievals within a swath have considerable temperature differences compared with that of spatially neighboring pixels that are not covered by the swath. Regions at high latitude experience longer cloud durations than those at lower latitudes; thus, pixels that are not covered by the swath might accumulate high uncertainties compared with the adjacent clear-sky retrievals. Therefore, an artificial texture remained after the data assimilation. Based on the literature review, spatial continuity is also a major issue for current MODIS LST products due to substantially different view zenith angles and view times of neighboring swaths after reprojection (Figs. 14 and 17 in Z.-L. Li et al., 2022). To address such discontinuity issues, machine learning methods that incorporate additional variables for estimating cloudy-sky LST can be used in the future (Zhao et al., 2019).

5 Data availability

The GHA-LST dataset from 2011 to 2021 is freely available at https://doi.org/10.5281/zenodo.7487284 (Jia et al., 2022b), as well as at the University of Maryland (http://glass.umd.edu/allsky_LST/GHA-LST, last access: 10 February 2023; Jia et al., 2022c). Quality check (QC) flags are also included: the “Bit 0” indicates the sky condition mask (1 for clear sky and 0 for cloudy sky), and the “Bit 1” is the cloud duration that represents the uncertainty level (0 denotes ≤10 d and 1 denotes >10 d). The algorithm produced a stable accuracy within 10 d, as indicated by Jia et al. (2022a).

6 Conclusions

LST is an essential driving factor in the surface radiation budget and hydrological cycling, and TIR-based satellite retrieval is the primary method used to obtain LST globally. However, TIR-derived LST has numerous data gaps, mainly due to cloud cover, which seriously restricts the application of current LST products. Studies have focused on producing all-sky LST products; nevertheless, global all-sky LSTs on an hourly scale are still unavailable. Considering the high temporal variability in the LST and the importance of the DTC in evapotranspiration (ET) estimation, drought detection, and heat wave monitoring, we produced a global hourly, 5 km, all-sky land surface temperature dataset (GHA-LST) from 2011 to 2021. The data are recovered using CGLS LST products from geostationary satellites and the MYD/MOD21 products from polar-orbiting satellites. Moreover, GHA-LST represents the first global gap-free LST product at an hourly scale, and it has been comprehensively validated by in situ measurements at 201 global sites in this study.

Based on the ground measurements from the SURFRAD, BSRN, Fluxnet, AmeriFlux, HRB, and TP networks, the overall RMSE of GHA-LST is 3.31 K, with a bias of −0.57 K and R² of 0.95. The comparisons of individual accuracy suggest that the proposed GHA-LST dataset has better accuracy than the CGLS and MYD21C1 data under both clear-sky and cloudy-sky conditions. In addition, the accuracy is stable under both clear-sky and cloudy-sky conditions (RMSEs of ∼3.3 K with few biases based on different sampling groups). The overall sampling amount was more than 5 million, and we suppose that the overall validation can represent the general accuracy of GHA-LST globally. In addition, after temporal aggregation to the daily mean scale, the GHA-LST dataset produced an RMSE of 1.76 K and significantly improved the accuracy and data availability of the global daily mean LST. The individual site validation indicated that the GHA-LST dataset has similar accuracy in terms of spatial patterns under different sky conditions. In comparison, the cloudy-sky results had a larger spatial variance in accuracy.

Temporal analyses were performed for four representative sites, and the GHA-LST dataset had a high temporal continuity and was consistent with the ground measurements at hourly and daily mean scales. The temporal variation in hourly LST also illustrated that mathematically predictable DTCs cannot be obtained for locations with continuous cloudy days, thereby highlighting the advantage of the time-evolving model-based method used for the GHA-LST product. Spatial analysis suggested that the GHA-LST dataset has satisfactory spatial continuity over clear-sky and cloudy-sky regions, and artificial textures were not observed. Regional mapping analysis of the TP and Alaska regions demonstrated that GHA-LST can capture greater spatial detail than reanalysis datasets, which were the only data source for obtaining hourly gap-free LSTs before this study. The monthly anomaly analysis indicated that the GHA-LST anomalies are consistent with global AT datasets and other LST datasets at the T_mean and T_min scales, whereas the T_max and DTR of the LSTs and ATs are only consistent with the anomalous direction when the magnitudes are different.

In the future, additional clear-sky LST products, such as the Visible Infrared Imaging Radiometer Suite (VIIRS) and Advanced Very High Resolution Radiometer (AVHRR), can also be assimilated to increase the time span and spatial resolution of the proposed dataset. Machine learning can be employed to effectively incorporate information from ground measurements, spatial textures, and related factors (e.g., elevation, soil moisture, land cover, and wind speed). GHA-LST represents the first gap-free LST dataset at an hourly, 5 km scale over the globe, and it has satisfactory accuracy and great potential for use in estimating global plant water stress, monitoring extreme weather, and advancing meteorological forecasting models.

Appendix A

Table A1 Validation statistics of each site.